1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: caspt2energy.cc
4 // Copyright (C) 2016 Toru Shiozaki
5 //
6 // Author: Jae Woo Park <jwpk1201@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 #include <bagel_config.h>
26 
27 #include <src/scf/hf/fock.h>
28 #include <src/grad/cpcasscf.h>
29 #include <src/grad/gradeval.h>
30 #include <src/grad/finite.h>
31 #include <src/multi/casscf/cassecond.h>
32 #include <src/multi/casscf/casnoopt.h>
33 #include <src/multi/casscf/qvec.h>
34 #include <src/smith/smith.h>
35 #include <src/smith/caspt2grad.h>
36 #include <src/prop/multipole.h>
37 #include <src/prop/hyperfine.h>
38 
39 
40 using namespace std;
41 using namespace bagel;
42 
CASPT2Energy(shared_ptr<const PTree> inp,shared_ptr<const Geometry> geom,shared_ptr<const Reference> ref)43 CASPT2Energy::CASPT2Energy(shared_ptr<const PTree> inp, shared_ptr<const Geometry> geom, shared_ptr<const Reference> ref)
44   : Method(inp, geom, ref) {
45 #ifdef COMPILE_SMITH
46   Timer timer;
47 
48   // compute CASSCF first
49   if (inp->get<string>("algorithm", "") != "noopt") {
50     auto cas = make_shared<CASSecond>(inp, geom, ref);
51     cas->compute();
52     ref_ = cas->conv_to_ref();
53     fci_ = cas->fci();
54     thresh_ = cas->thresh();
55   } else {
56     auto cas = make_shared<CASNoopt>(inp, geom, ref);
57     cas->compute();
58     ref_ = cas->conv_to_ref();
59     fci_ = cas->fci();
60     thresh_ = cas->thresh();
61   }
62 
63   // gradient/property calculation
64   do_hyperfine_ = inp->get<bool>("hyperfine", false);
65   target_state1_ = inp->get<int>("_target", -1);
66   target_state2_ = inp->get<int>("_target2", -1);
67 
68   timer.tick_print("Reference calculation");
69 
70   cout << endl << "  === DF-CASPT2 calculation ===" << endl << endl;
71 #else
72   throw logic_error("single point CASPT2 require SMITH-generated code. Please compile BAGEL with --enable-smith");
73 #endif
74 }
75 
76 
77 // compute smith and set rdms and ci deriv to a member
compute()78 void CASPT2Energy::compute() {
79 #ifdef COMPILE_SMITH
80   shared_ptr<PTree> smithinput = idata_->get_child("smith");
81   smithinput->put("_grad", false);
82   smithinput->put("_hyperfine", do_hyperfine_);
83 
84   if (target_state1_!=-1) {
85     smithinput->put("_target", target_state1_);
86     smithinput->put("_target2", target_state2_);
87   }
88   auto smith = make_shared<Smith>(smithinput, ref_->geom(), ref_);
89   smith->compute();
90 
91   energy_  = smith->algo()->energyvec();
92   ref_->energy() = energy_;
93 
94   if (target_state1_!=-1) {
95     smith->compute_gradient(target_state1_, target_state2_, make_shared<const NacmType>());
96     ncore_   = smith->algo()->info()->ncore();
97     coeff_   = smith->coeff();
98     msrot_   = smith->msrot();
99 
100     auto d1set = [this](shared_ptr<const Matrix> d1t) {
101       if (!ncore_) {
102         return d1t->copy();
103       } else {
104         auto out = make_shared<Matrix>(coeff_->mdim(), coeff_->mdim());
105         out->copy_block(ncore_, ncore_, coeff_->mdim()-ncore_, coeff_->mdim()-ncore_, d1t);
106         return out;
107       }
108     };
109     auto vd1tmp = make_shared<Matrix>(*smith->vd1());
110     vd1_ = d1set(vd1tmp);
111   }
112 #endif
113 }
114 
115 template<>
compute()116 shared_ptr<GradFile> FiniteNacm<CASPT2Energy>::compute() {
117 #ifdef COMPILE_SMITH
118   const int natom = geom_->natom();
119   cout << "  Derivative coupling evaluation with respect to " << natom * 3 << " DOFs" << endl;
120   cout << "  Finite difference size (dx) is " << setprecision(8) << dx_ << " Bohr" << endl;
121 
122   shared_ptr<Dvec> civ_ref = ref_->civectors()->copy();
123   civ_ref->rotate (task_->msrot());
124   civ_ref->print (/*sort=*/false);
125 
126   const int nclosed = ref_->nclosed();
127   const int nocc = ref_->nocc();
128 
129   Timer timer;
130   muffle_ = make_shared<Muffle>("finite.log");
131 
132   auto grad = make_shared<GradFile>(natom);
133 
134   auto acoeff_ref = make_shared<Matrix>(task_->coeff()->slice(nclosed, nocc));
135   auto coeff_ref  = make_shared<Matrix>(*task_->coeff());
136 
137   const int norb = civ_ref->det()->norb();
138   auto gmo = make_shared<Matrix>(norb, norb);
139   auto vd1a = make_shared<Matrix>(*task_->vd1());
140 
141   auto idata_out = std::make_shared<PTree>(*idata_);
142   idata_out->put("_target", target_state1_);
143   idata_out->put("_target2", target_state2_);
144 
145   const int ncomm = mpi__->world_size() / nproc_;
146   const int icomm = mpi__->world_rank() / nproc_;
147 
148   mpi__->split(nproc_);
149 
150   for (int i = 0, counter = 0; i != natom; ++i) {    // atom i
151     for (int j = 0; j != 3; ++j, ++counter) {      // xyz
152       if (counter % ncomm == icomm && ncomm != icomm) {
153         muffle_->mute();
154         shared_ptr<Dvec> civ_plus;
155         shared_ptr<Matrix> coeff_plus;
156         {
157           auto displ = make_shared<XYZFile>(natom);
158           displ->element(j,i) = dx_;
159           auto geom_plus = make_shared<Geometry>(*geom_, displ, idata_, false, false);
160           geom_plus->print_atoms();
161 
162           shared_ptr<const Reference> ref_plus;
163           if (ref_)
164             ref_plus = ref_->project_coeff(geom_plus);
165 
166           task_ = std::make_shared<CASPT2Energy>(idata_out, geom_plus, ref_plus);
167           task_->compute();
168           ref_plus  = task_->conv_to_ref();
169 
170           coeff_plus = make_shared<Matrix>(*task_->coeff());
171 
172           for (int im = 0; im != coeff_ref->mdim(); ++im) {
173             double dmatch = blas::dot_product(coeff_ref->element_ptr(0, im), coeff_ref->ndim(), coeff_plus->element_ptr(0,im));
174             if (dmatch < 0.0) {
175               blas::scale_n(-1.0, coeff_plus->element_ptr(0, im), coeff_ref->ndim());
176             }
177           }
178 
179           civ_plus = ref_plus->civectors()->copy();
180           civ_plus->rotate(task_->msrot());
181           civ_plus->match(civ_ref);
182         }
183 
184         shared_ptr<Dvec> civ_minus;
185         shared_ptr<Matrix> coeff_minus;
186         {
187           auto displ = make_shared<XYZFile>(natom);
188           displ->element(j,i) = -dx_;
189           auto geom_minus = make_shared<Geometry>(*geom_, displ, idata_, false, false);
190           geom_minus->print_atoms();
191 
192           shared_ptr<const Reference> ref_minus;
193           if (ref_)
194             ref_minus = ref_->project_coeff(geom_minus);
195 
196           task_ = std::make_shared<CASPT2Energy>(idata_out, geom_minus, ref_minus);
197           task_->compute();
198           ref_minus  = task_->conv_to_ref();
199 
200           coeff_minus = make_shared<Matrix>(*task_->coeff());
201 
202           for (int im = 0; im != coeff_ref->mdim(); ++im) {
203             double dmatch = blas::dot_product(coeff_ref->element_ptr(0, im), coeff_ref->ndim(), coeff_minus->element_ptr(0,im));
204             if (dmatch < 0.0) {
205               blas::scale_n(-1.0, coeff_minus->element_ptr(0, im), coeff_ref->ndim());
206             }
207           }
208 
209           civ_minus = ref_minus->civectors()->copy();
210           civ_minus->rotate(task_->msrot());
211           civ_minus->match(civ_ref);
212         }
213 
214         auto civ_diff = civ_plus->copy();
215         *civ_diff -= *civ_minus;
216         civ_diff->scale(1.0 / (2.0 * dx_));
217         civ_diff->print(/*sort=*/false);
218         auto coeff_diff = make_shared<Matrix>(*coeff_plus - *coeff_minus);
219         coeff_diff->scale(1.0 / (2.0 * dx_));
220         auto acoeff_diff = make_shared<Matrix>(coeff_diff->slice(nclosed, nocc));
221 
222         auto Smn = make_shared<Overlap>(geom_);
223         auto Uij = make_shared<Matrix>(*acoeff_ref % *Smn * *acoeff_diff);
224         if (mpi__->rank() == 0) {
225           grad->element(j,i) = civ_ref->data(target_state1_)->dot_product(civ_diff->data(target_state2_));
226           for (int ii = 0; ii != norb; ++ii) {
227             for (int ij = 0; ij != norb; ++ij) {
228               const int lena = civ_ref->det()->lena();
229               const int lenb = civ_ref->det()->lenb();
230               if (ii != ij) {
231                 for (auto& iter : civ_ref->det()->phia(ii, ij)) {
232                   size_t iaA = iter.source;
233                   size_t iaB = iter.target;
234                   double sign = static_cast<double>(iter.sign);
235                   for (size_t ib = 0; ib != lenb; ++ib) {
236                     double factor = civ_ref->data(target_state1_)->data(ib+iaB*lenb) * civ_ref->data(target_state2_)->data(ib+iaA*lenb) * sign;
237                     grad->element(j,i) += factor * (Uij->element(ij, ii) - Uij->element(ii, ij)) * .5;
238                     if ((i + j * 3) == 0) {
239                       gmo->element(ij, ii) += factor * .5;
240                       gmo->element(ii, ij) -= factor * .5;
241                     }
242                   }
243                 }
244                 for (size_t ia = 0; ia != lena; ++ia) {
245                   for (auto& iter : civ_ref->det()->phib(ii, ij)) {
246                     size_t ibA = iter.source;
247                     size_t ibB = iter.target;
248                     double sign = static_cast<double>(iter.sign);
249                     double factor = civ_ref->data(target_state1_)->data(ibB+ia*lenb) * civ_ref->data(target_state2_)->data(ibA+ia*lenb) * sign;
250                     grad->element(j,i) += factor * (Uij->element(ij, ii) - Uij->element(ii, ij)) * .5;
251                     if ((i + j * 3) == 0) {
252                       gmo->element(ij, ii) += factor * .5;
253                       gmo->element(ii, ij) -= factor * .5;
254                     }
255                   }
256                 }
257               }
258             }
259           }
260           const int nmobasis = task_->coeff()->ndim();
261           auto Ifactor = make_shared<Matrix>(*coeff_ref % *Smn * *coeff_diff);
262           for (int ii = 0; ii != nmobasis; ++ii) {
263             for (int ij = 0; ij != nmobasis; ++ij) {
264               grad->element(j,i) += vd1a->element(ij, ii) * Ifactor->element(ij, ii);
265             }
266           }
267         }
268       }
269       muffle_->unmute();
270       stringstream ss; ss << "Finite difference evaluation (" << setw(2) << i*3+j+1 << " / " << geom_->natom() * 3 << ")";
271       timer.tick_print(ss.str());
272     }
273   }
274   mpi__->merge();
275   grad->allreduce();
276   gmo->allreduce();
277 
278   auto gfin = make_shared<Matrix>((*acoeff_ref * (*gmo) ^ *acoeff_ref) + (*coeff_ref * (*vd1a) ^ *coeff_ref));
279   auto grad_basis = make_shared<GradFile>(geom_->natom());
280   grad_basis = contract_gradient(nullptr, nullptr, nullptr, nullptr, gfin, /*numerical=*/true);
281 
282   *grad += *grad_basis;
283   grad->print(": NACME calculated with finite difference", 0);
284   return grad;
285 #else
286   return nullptr;
287 #endif
288 }
289 
conv_to_ref() const290 shared_ptr<const Reference> CASPT2Energy::conv_to_ref() const {
291  return std::make_shared<Reference>(ref_->geom(), ref_->coeff(), ref_->nclosed(), ref_->nact(), ref_->nvirt(), energy_,
292                                fci_->rdm1(), fci_->rdm2(), fci_->rdm1_av(), fci_->rdm2_av(), fci_->conv_to_ciwfn());
293 }
294