1 // 2 // BAGEL - Brilliantly Advanced General Electronic Structure Library 3 // Filename: civec.h 4 // Copyright (C) 2011 Toru Shiozaki 5 // 6 // Author: Toru Shiozaki <shiozaki@northwestern.edu> 7 // Maintainer: Shiozaki group 8 // 9 // This file is part of the BAGEL package. 10 // 11 // This program is free software: you can redistribute it and/or modify 12 // it under the terms of the GNU General Public License as published by 13 // the Free Software Foundation, either version 3 of the License, or 14 // (at your option) any later version. 15 // 16 // This program is distributed in the hope that it will be useful, 17 // but WITHOUT ANY WARRANTY; without even the implied warranty of 18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 // GNU General Public License for more details. 20 // 21 // You should have received a copy of the GNU General Public License 22 // along with this program. If not, see <http://www.gnu.org/licenses/>. 23 // 24 25 26 #ifndef BAGEL_FCI_CIVEC_H 27 #define BAGEL_FCI_CIVEC_H 28 29 #include <list> 30 #include <src/util/math/algo.h> 31 #include <src/util/f77.h> 32 #include <src/util/parallel/staticdist.h> 33 #include <src/ci/fci/determinants.h> 34 #include <src/ci/fci/dvector_base.h> 35 36 namespace bagel { 37 38 template<typename DataType> 39 class Civector { 40 public: using DetType = Determinants; // used to automatically determine type for Determinants object in templates 41 public: using LocalizedType = std::true_type; 42 43 protected: 44 // The determinant space in which this Civec object is defined 45 mutable std::shared_ptr<const Determinants> det_; 46 47 size_t lena_; 48 size_t lenb_; 49 50 // !!CAUTION!! 51 // cc is formated so that B runs first. 52 // Also, cc_ can be null if this is constructed by Dvec. 53 std::unique_ptr<DataType[]> cc_; 54 55 DataType* cc_ptr_; 56 cc(const size_t & i)57 DataType& cc(const size_t& i) { return *(cc_ptr_+i); } cc(const size_t & i)58 const DataType& cc(const size_t& i) const { return *(cc_ptr_+i); } cc()59 DataType* cc() { return cc_ptr_; } cc()60 const DataType* cc() const { return cc_ptr_; } 61 62 private: 63 friend class boost::serialization::access; 64 template<class Archive> serialize(Archive & ar,const unsigned int version)65 void serialize(Archive& ar, const unsigned int version) { 66 boost::serialization::split_member(ar, *this, version); 67 } 68 template<class Archive> save(Archive & ar,const unsigned int)69 void save(Archive& ar, const unsigned int) const { 70 if (!cc_.get()) 71 throw std::logic_error("illegal call of Civector<T>::save"); 72 ar << det_ << lena_ << lenb_ << make_array(cc(), size()); 73 } 74 template<class Archive> load(Archive & ar,const unsigned int)75 void load(Archive& ar, const unsigned int) { 76 ar >> det_ >> lena_ >> lenb_; 77 cc_ = std::unique_ptr<DataType[]>(new DataType[size()]); 78 cc_ptr_ = cc_.get(); 79 ar >> make_array(cc(), size()); 80 } 81 82 public: Civector()83 Civector() { } 84 Civector(std::shared_ptr<const Determinants> det); 85 // constructor that is called by Dvec. 86 Civector(std::shared_ptr<const Determinants> det, DataType* din_); 87 // copy constructor 88 Civector(const Civector<DataType>& o); 89 // this is not a copy constructor. 90 Civector(std::shared_ptr<Civector<DataType>> o, std::shared_ptr<const Determinants> det); 91 clone()92 std::shared_ptr<Civector<DataType>> clone() const { return std::make_shared<Civector<DataType>>(det_); } copy()93 std::shared_ptr<Civector<DataType>> copy() const { return std::make_shared<Civector<DataType>>(*this); } 94 data()95 DataType* data() { return cc(); } element(size_t i,size_t j)96 DataType& element(size_t i, size_t j) { return cc(i+j*lenb_); } element_ptr(size_t i,size_t j)97 DataType* element_ptr(size_t i, size_t j) { return cc()+i+j*lenb_; } 98 data(const size_t & i)99 DataType& data(const size_t& i) { return cc(i); } data(const size_t & i)100 const DataType& data(const size_t& i) const { return cc(i); } 101 data()102 const DataType* data() const { return cc(); } element_ptr(size_t i,size_t j)103 const DataType* element_ptr(size_t i, size_t j) const { return cc()+i+j*lenb_; } 104 det()105 std::shared_ptr<const Determinants> det() const { return det_; } set_det(std::shared_ptr<const Determinants> o)106 void set_det(std::shared_ptr<const Determinants> o) const { det_ = o; } 107 zero()108 void zero() { std::fill(cc(), cc()+lena_*lenb_, DataType(0.0)); } 109 size()110 size_t size() const { return lena_*lenb_; } 111 112 std::shared_ptr<Civector<DataType>> transpose(std::shared_ptr<const Determinants> det = nullptr) const; 113 lena()114 size_t lena() const { return lena_; } lenb()115 size_t lenb() const { return lenb_; } asize()116 size_t asize() const { return lena_; } 117 118 // some functions for convenience ax_plus_y(DataType a,const Civector<DataType> & other)119 void ax_plus_y(DataType a, const Civector<DataType>& other) { 120 assert((lena_ == other.lena_) && (lenb_ == other.lenb_)); 121 blas::ax_plus_y_n(a, other.data(), size(), cc()); 122 } ax_plus_y(DataType a,std::shared_ptr<const Civector> other)123 void ax_plus_y(DataType a, std::shared_ptr<const Civector> other) { ax_plus_y(a, *other); } 124 dot_product(const Civector<DataType> & other)125 DataType dot_product(const Civector<DataType>& other) const { 126 assert((lena_ == other.lena_) && (lenb_ == other.lenb_)); 127 return blas::dot_product(cc(), size(), other.data()); 128 } dot_product(std::shared_ptr<const Civector> other)129 DataType dot_product(std::shared_ptr<const Civector> other) const { return dot_product(*other); } 130 norm()131 double norm() const { return std::sqrt(detail::real(dot_product(*this))); } variance()132 double variance() const { return detail::real(dot_product(*this)) / size(); } rms()133 double rms() const { return std::sqrt(variance()); } 134 scale(const DataType a)135 void scale(const DataType a) { blas::scale_n(a, cc(), size()); } 136 137 // Spin functions are only implememted as specialized functions for double (see civec.cc) spin_expectation()138 DataType spin_expectation() const { // returns < S^2 > 139 std::shared_ptr<Civector<DataType>> S2 = spin(); 140 return dot_product(*S2); 141 } 142 std::shared_ptr<Civector<DataType>> spin() const; 143 144 // S_- = \sum_i i_beta^\dagger i_alpha 145 std::shared_ptr<Civector<DataType>> spin_lower(std::shared_ptr<const Determinants> target_det = nullptr) const; 146 147 // S_+ = \sum_i i_alpha^\dagger i_beta 148 std::shared_ptr<Civector<DataType>> spin_raise(std::shared_ptr<const Determinants> target_det = nullptr) const; 149 150 void spin_decontaminate(const double thresh = 1.0e-12); 151 152 std::shared_ptr<Civector<DataType>> apply(const int orbital, const bool action, const bool spin) const; 153 154 Civector<DataType>& operator*=(const double& a) { scale(a); return *this; } 155 Civector<DataType>& operator+=(const double& a) { std::for_each(cc(), cc()+size(), [&a](DataType& p){ p += a; }); return *this; } 156 Civector<DataType>& operator-=(const double& a) { std::for_each(cc(), cc()+size(), [&a](DataType& p){ p -= a; }); return *this; } 157 158 Civector<DataType>& operator=(const Civector<DataType>& o) { 159 assert(det()->lena() == o.det()->lena() && det()->lenb() == o.det()->lenb()); 160 std::copy_n(o.cc(), size(), cc()); 161 return *this; 162 } 163 Civector<DataType>& operator+=(const Civector<DataType>& o) { ax_plus_y( 1.0, o); return *this; } 164 Civector<DataType>& operator-=(const Civector<DataType>& o) { ax_plus_y(-1.0, o); return *this; } 165 Civector<DataType>& operator/=(const Civector<DataType>& o) { 166 for (size_t i = 0; i != size(); ++i) 167 data(i) /= o.data(i); 168 return *this; 169 } 170 Civector<DataType> operator/(const Civector<DataType>& o) const { 171 Civector<DataType> out(*this); 172 out /= o; 173 return out; 174 } 175 176 // assumes that Civec's in c are already orthogonal with each other. 177 // returns scaling factor (see implementation) 178 orthog(std::list<std::shared_ptr<const Civector<DataType>>> c)179 double orthog(std::list<std::shared_ptr<const Civector<DataType>>> c) { 180 for (auto& iter : c) 181 project_out(iter); 182 return normalize(); 183 } 184 orthog(std::shared_ptr<const Civector<DataType>> o)185 double orthog(std::shared_ptr<const Civector<DataType>> o) { 186 return orthog(std::list<std::shared_ptr<const Civector<DataType>>>{o}); 187 } 188 normalize()189 double normalize() { 190 const double norm = this->norm(); 191 const double scal = (norm*norm<1.0e-60 ? 0.0 : 1.0/norm); 192 scale(scal); 193 return 1.0/scal; 194 } 195 project_out(std::shared_ptr<const Civector<DataType>> o)196 void project_out(std::shared_ptr<const Civector<DataType>> o) { ax_plus_y(-detail::conj(dot_product(*o)), *o); } 197 198 void print(const double thr, const bool sort = true) const; 199 200 void synchronize(const int root = 0) { 201 mpi__->broadcast(cc_ptr_, size(), root); 202 } 203 }; 204 205 template<> void Civector<double>::spin_decontaminate(const double thresh); 206 template<> void Civector<std::complex<double>>::spin_decontaminate(const double thresh); 207 208 using Civec = Civector<double>; 209 using ZCivec = Civector<std::complex<double>>; 210 211 template <> 212 void Dvector_base<Civec>::apply_and_fill(std::shared_ptr<const Dvector_base<Civec>> source_dvec, const int orbital, const bool action, const bool spin); 213 214 using CASDvec = Dvector_base<Civec>; 215 216 } 217 218 extern template class bagel::Civector<double>; 219 extern template class bagel::Civector<std::complex<double>>; 220 221 #endif 222