1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: finite.cc
4 // Copyright (C) 2012 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 <src/grad/gradeval.h>
26 #include <src/grad/finite.h>
27 #include <src/util/timer.h>
28 #include <src/wfn/get_energy.h>
29
30 using namespace std;
31 using namespace bagel;
32
compute()33 shared_ptr<GradFile> FiniteGrad::compute() {
34 for (auto& m : *idata_) {
35 const string title = to_lower(m->get<string>("title", ""));
36 tie(energy_, ref_) = get_energy(title, m, geom_, ref_, target_state_);
37 }
38
39 const int natom = geom_->natom();
40 cout << " Gradient evaluation with respect to " << natom * 3 << " DOFs" << endl;
41 cout << " Finite difference size (dx) is " << setprecision(8) << dx_ << " Bohr" << endl;
42
43 Timer timer;
44 muffle_ = make_shared<Muffle>("finite.log");
45
46 auto grad = make_shared<GradFile>(natom);
47
48 const int ncomm = mpi__->world_size() / nproc_;
49 const int icomm = mpi__->world_rank() / nproc_;
50
51 mpi__->split(nproc_);
52
53 for (int i = 0, counter = 0; i != natom; ++i) { // atom i
54 for (int j = 0; j != 3; ++j, ++counter) { // xyz
55 if (counter % ncomm == icomm && ncomm != icomm) {
56 muffle_->mute();
57
58 double energy_plus = 0.0;
59 {
60 auto displ = make_shared<XYZFile>(natom);
61 displ->element(j,i) = dx_;
62 auto geom_plus = make_shared<Geometry>(*geom_, displ, make_shared<PTree>(), false, false);
63 geom_plus->print_atoms();
64
65 shared_ptr<const Reference> ref_plus;
66 if (ref_)
67 ref_plus = ref_->project_coeff(geom_plus);
68
69 for (auto& m : *idata_) {
70 const string title = to_lower(m->get<string>("title", ""));
71 tie(energy_plus, ref_plus) = get_energy(title, m, geom_plus, ref_plus, target_state_);
72 }
73 }
74
75 double energy_minus = 0.0;
76 {
77 auto displ = make_shared<XYZFile>(natom);
78 displ->element(j,i) = -dx_;
79 auto geom_minus = make_shared<Geometry>(*geom_, displ, make_shared<PTree>(), false, false);
80 geom_minus->print_atoms();
81
82 shared_ptr<const Reference> ref_minus;
83 if (ref_)
84 ref_minus = ref_->project_coeff(geom_minus);
85
86 for (auto& m : *idata_) {
87 const string title = to_lower(m->get<string>("title", ""));
88 tie(energy_minus, ref_minus) = get_energy(title, m, geom_minus, ref_minus, target_state_);
89 }
90 }
91
92 if (mpi__->rank() == 0)
93 grad->element(j,i) = (energy_plus - energy_minus) / (2.0 * dx_);
94 muffle_->unmute();
95 stringstream ss; ss << "Finite difference evaluation (" << setw(2) << i*3+j+1 << " / " << geom_->natom() * 3 << ")";
96 timer.tick_print(ss.str());
97 }
98 }
99 }
100 mpi__->merge();
101 grad->allreduce();
102
103 grad->print(": Calculated with finite difference", 0);
104 return grad;
105 }
106
107 template<>
compute()108 shared_ptr<GradFile> FiniteNacm<CASSCF>::compute() {
109 const int natom = geom_->natom();
110 cout << " Derivative coupling evaluation with respect to " << natom * 3 << " DOFs" << endl;
111 cout << " Finite difference size (dx) is " << setprecision(8) << dx_ << " Bohr" << endl;
112
113 shared_ptr<Dvec> civ_ref = ref_->civectors()->copy();
114 civ_ref->print (/*sort=*/false);
115
116 const int nclosed = ref_->nclosed();
117 const int nocc = ref_->nocc();
118
119 Timer timer;
120 muffle_ = make_shared<Muffle>("finite.log");
121
122 auto acoeff_ref = make_shared<Matrix>(ref_->coeff()->slice(nclosed, nocc));
123 auto grad = make_shared<GradFile>(natom);
124 const int norb = civ_ref->det()->norb();
125 auto gmo = make_shared<Matrix>(norb, norb);
126 gmo->zero();
127
128 const int ncomm = mpi__->world_size() / nproc_;
129 const int icomm = mpi__->world_rank() / nproc_;
130
131 mpi__->split(nproc_);
132
133 for (int i = 0, counter = 0; i != natom; ++i) { // atom i
134 for (int j = 0; j != 3; ++j, ++counter) { // xyz
135 if (counter % ncomm == icomm && ncomm != icomm) {
136 muffle_->mute();
137
138 shared_ptr<Matrix> acoeff_plus;
139 shared_ptr<Dvec> civ_plus;
140 {
141 auto displ = make_shared<XYZFile>(natom);
142 displ->element(j,i) = dx_;
143 auto geom_plus = make_shared<Geometry>(*geom_, displ, make_shared<PTree>(), false, false);
144 geom_plus->print_atoms();
145
146 shared_ptr<const Reference> ref_plus;
147 if (ref_)
148 ref_plus = ref_->project_coeff(geom_plus);
149
150 double energy_plus;
151 tie(energy_plus, ref_plus) = get_energy("casscf", idata_, geom_plus, ref_plus);
152 acoeff_plus = make_shared<Matrix>(ref_plus->coeff()->slice(nclosed, nocc));
153 for (int im = 0; im != acoeff_ref->mdim(); ++im) {
154 double dmatch = blas::dot_product(acoeff_ref->element_ptr(0, im), acoeff_ref->ndim(), acoeff_plus->element_ptr(0,im));
155 if (dmatch < 0.0)
156 blas::scale_n(-1.0, acoeff_plus->element_ptr(0, im), acoeff_ref->ndim());
157 }
158 civ_plus = ref_plus->civectors()->copy();
159 civ_plus->match(civ_ref);
160 }
161
162 shared_ptr<Matrix> acoeff_minus;
163 shared_ptr<Dvec> civ_minus;
164 {
165 auto displ = make_shared<XYZFile>(natom);
166 displ->element(j,i) = -dx_;
167 auto geom_minus = make_shared<Geometry>(*geom_, displ, make_shared<PTree>(), false, false);
168 geom_minus->print_atoms();
169
170 shared_ptr<const Reference> ref_minus;
171 if (ref_)
172 ref_minus = ref_->project_coeff(geom_minus);
173
174 double energy_minus;
175 tie(energy_minus, ref_minus) = get_energy("casscf", idata_, geom_minus, ref_minus);
176 acoeff_minus = make_shared<Matrix>(ref_minus->coeff()->slice(nclosed, nocc));
177 for (int im = 0; im != acoeff_ref->mdim(); ++im) {
178 double dmatch = blas::dot_product(acoeff_ref->element_ptr(0, im), acoeff_ref->ndim(), acoeff_minus->element_ptr(0,im));
179 if (dmatch < 0.0)
180 blas::scale_n(-1.0, acoeff_minus->element_ptr(0, im), acoeff_ref->ndim());
181 }
182 civ_minus = ref_minus->civectors()->copy();
183 civ_minus->match(civ_ref);
184 }
185
186 auto civ_diff = civ_plus->copy();
187 *civ_diff -= *civ_minus;
188 civ_diff->scale(1.0 / (2.0 * dx_));
189 auto acoeff_diff = make_shared<Matrix>(*acoeff_plus - *acoeff_minus);
190 acoeff_diff->scale(1.0 / (2.0 * dx_));
191
192 auto Smn = make_shared<Overlap>(geom_);
193 auto Uij = make_shared<Matrix>(*acoeff_ref % *Smn * *acoeff_diff);
194 if (mpi__->rank() == 0) {
195 grad->element(j,i) = civ_ref->data(target_state1_)->dot_product(civ_diff->data(target_state2_));
196
197 for (int ii = 0; ii != norb; ++ii) {
198 for (int ij = 0; ij != norb; ++ij) {
199 const int lena = civ_ref->det()->lena();
200 const int lenb = civ_ref->det()->lenb();
201 if (ii != ij) {
202 for (auto& iter : civ_ref->det()->phia(ii, ij)) {
203 size_t iaA = iter.source;
204 size_t iaB = iter.target;
205 double sign = static_cast<double>(iter.sign);
206
207 for (size_t ib = 0; ib != lenb; ++ib) {
208 double factor = civ_ref->data(target_state1_)->data(ib+iaB*lenb) * civ_ref->data(target_state2_)->data(ib+iaA*lenb) * sign;
209 grad->element(j,i) += factor * (Uij->element(ij, ii) - Uij->element(ii, ij)) * .5;
210 if ((i + j * 3) == 0) {
211 gmo->element(ij, ii) += factor * .5;
212 gmo->element(ii, ij) -= factor * .5;
213 }
214 }
215 }
216 for (size_t ia = 0; ia != lena; ++ia) {
217 for (auto& iter : civ_ref->det()->phib(ii, ij)) {
218 size_t ibA = iter.source;
219 size_t ibB = iter.target;
220 double sign = static_cast<double>(iter.sign);
221 double factor = civ_ref->data(target_state1_)->data(ibB+ia*lenb) * civ_ref->data(target_state2_)->data(ibA+ia*lenb) * sign;
222 grad->element(j,i) += factor * (Uij->element(ij, ii) - Uij->element(ii, ij)) * .5;
223 if ((i + j * 3) == 0) {
224 gmo->element(ij, ii) += factor * .5;
225 gmo->element(ii, ij) -= factor * .5;
226 }
227 }
228 }
229 }
230 }
231 }
232 }
233 }
234 muffle_->unmute();
235 stringstream ss; ss << "Finite difference evaluation (" << setw(2) << i*3+j+1 << " / " << geom_->natom() * 3 << ")";
236 timer.tick_print(ss.str());
237 }
238 }
239 mpi__->merge();
240 grad->allreduce();
241 gmo->allreduce();
242 auto gfin = make_shared<Matrix>(*acoeff_ref * *gmo ^ *acoeff_ref);
243 auto grad_basis = make_shared<GradFile>(natom);
244 grad_basis = contract_gradient(nullptr, nullptr, nullptr, nullptr, gfin, /*numerical=*/true);
245 *grad += *grad_basis;
246
247 grad->print(": NACME calculated with finite difference", 0);
248
249 return grad;
250 }
251