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