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