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