1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: smith.cc
4 // Copyright (C) 2013 Toru Shiozaki
5 //
6 // Author: Matthew K. MacLeod <matthew.macleod@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/smith/smith.h>
28
29 #ifdef COMPILE_SMITH
30 #include <src/smith/mrci/MRCI.h>
31 #include <src/smith/relmrci/RelMRCI.h>
32 #include <src/smith/caspt2/CASPT2.h>
33 #include <src/smith/caspt2/SPCASPT2.h>
34 #include <src/smith/relcaspt2/RelCASPT2.h>
35 #include <src/smith/casa/CASA.h>
36 #include <src/smith/relcasa/RelCASA.h>
37 using namespace bagel::SMITH;
38 #endif
39 using namespace std;
40 using namespace bagel;
41
Smith(const shared_ptr<const PTree> idata,shared_ptr<const Geometry> g,shared_ptr<const Reference> r)42 Smith::Smith(const shared_ptr<const PTree> idata, shared_ptr<const Geometry> g, shared_ptr<const Reference> r) : Method(idata, g, r) {
43 const string method = to_lower(idata_->get<string>("method", "caspt2"));
44
45 #ifdef COMPILE_SMITH
46 // print a header
47 if (!idata->get<bool>("_grad", false))
48 cout << " === SMITH program ===" << endl << endl;
49
50 // make a smith_info class
51 auto info = make_shared<const SMITH_Info<double>>(r, idata);
52
53 {
54 // print memory requirements
55 const size_t nact = r->nact();
56 const size_t nact2 = nact * nact;
57 const size_t nact4 = nact2 * nact2;
58 const size_t nact6 = nact4 * nact2;
59 const size_t nstate = r->nstate();
60 const size_t nclosed = r->nclosed() - info->ncore();
61 const size_t nvirtual = r->nvirt() - info->nfrozenvirt();
62 const size_t nocc = nclosed + nact;
63 const size_t norb = nocc + nvirtual;
64 const size_t rsize = nclosed*nocc*nvirtual*nvirtual + nclosed*nclosed*nact*(nvirtual+nact) + nact2*(norb*nvirtual+nact*nclosed);
65 cout << " * Approximate memory requirement for SMITH calulation per MPI process:" << endl;
66 cout << " o Storage requirement for T-amplitude, lambda, and residual is ";
67 if (info->orthogonal_basis()) {
68 // T-amp (r, o), Lambda (r, o), Residual (r, o)
69 cout << setprecision(2) << rsize*(info->sssr() ? nstate : nstate*nstate) * (info->grad() ? 6 : 4) * 8.e-9 / mpi__->size() << " GB" << endl;
70 } else {
71 // T-amp (r), Lambda (r), Residual (r)
72 cout << setprecision(2) << rsize*(info->sssr() ? nstate : nstate*nstate) * (info->grad() ? 3 : 2) * 8.e-9 / mpi__->size() << " GB" << endl;
73 }
74 cout << " o Storage requirement for MO integrals is ";
75 cout << setprecision(2) << (norb*norb*2 + nocc*nocc*(nact+nvirtual)*(nact+nvirtual)) * 8.e-9 / mpi__->size() << " GB" << endl;
76 if (info->grad()) {
77 cout << " o Storage requirement for SMITH-computed gradient tensors is ";
78 cout << setprecision(2) << nstate*nstate*(nact6*2 + nact4 + nact2 + 1) * 8.e-9 << " GB" << endl;
79 }
80 }
81
82 if (info->restart())
83 cout << " ** Restarting calculations is currently unavailable for non-relativistic SMITH methods. Serialized archives will not be generated." << endl << endl;
84
85 if (method == "caspt2") {
86 algo_ = make_shared<CASPT2::CASPT2>(info);
87 } else if (method == "casa") {
88 algo_ = make_shared<CASA::CASA>(info);
89 } else if (method == "mrci") {
90 algo_ = make_shared<MRCI::MRCI>(info);
91 } else {
92 #else
93 {
94 #endif
95 stringstream ss; ss << method << " method is not implemented in SMITH";
96 throw logic_error(ss.str());
97 }
98 }
99
100
101 void Smith::compute() {
102 #ifdef COMPILE_SMITH
103 algo_->solve();
104 #else
105 throw logic_error("You must enable SMITH during compilation for this method to be available.");
106 #endif
107 }
108
109
110 void Smith::compute_gradient(const int istate, const int jstate, shared_ptr<const NacmType> nacmtype, const bool nocider) {
111 #ifdef COMPILE_SMITH
112 if (!(algo_->info()->grad())) {
113 auto algop = dynamic_pointer_cast<CASPT2::CASPT2>(algo_);
114 algop->solve_dm(istate, jstate);
115 msrot_ = algop->msrot();
116 coeff_ = algop->coeff();
117 vd1_ = algop->vden1();
118 } else {
119 auto algop = make_shared<CASPT2::CASPT2>(*(dynamic_pointer_cast<CASPT2::CASPT2>(algo_)));
120 assert(algop);
121
122 algop->solve_gradient(istate, jstate, nacmtype, nocider);
123 dm1_ = algop->rdm12();
124 dm11_ = algop->rdm11();
125 dm2_ = algop->rdm21();
126 dcheck_ = algop->dcheck();
127 if (istate != jstate)
128 vd1_ = algop->vden1();
129
130 // compute <1|1>
131 wf1norm_ = algop->correlated_norm();
132 // convert ci derivative tensor to civec
133 cider_ = algop->ci_deriv();
134 msrot_ = algop->msrot();
135 coeff_ = algop->coeff();
136
137 // if spin-density is requested...
138 if (idata_->get<bool>("_hyperfine")) {
139 auto sp = make_shared<SPCASPT2::SPCASPT2>(*algop);
140 sp->solve();
141 sdm1_ = make_shared<Matrix>(*sp->rdm12() * 2.0 - *dm1_); // CAUTION! dm1 includes <1|1>D0 where as sp->rdm12() does not
142 sdm11_ = make_shared<Matrix>(*sp->rdm11() * 2.0 - *dm11_);
143 }
144 }
145 #else
146 throw logic_error("You must enable SMITH during compilation for this method to be available.");
147 #endif
148 }
149
150
151 RelSmith::RelSmith(const shared_ptr<const PTree> idata, shared_ptr<const Geometry> g, shared_ptr<const Reference> r) : Method(idata, g, r) {
152 #ifdef COMPILE_SMITH
153 // print a header
154 if (!idata->get<bool>("_grad", false))
155 cout << " === SMITH program ===" << endl << endl;
156
157 const string method = to_lower(idata_->get<string>("method", "caspt2"));
158 if (!dynamic_pointer_cast<const RelReference>(r) && method != "continue")
159 throw runtime_error("Relativistic correlation methods require a fully relativistic reference wavefunction.");
160
161 if (method != "continue") {
162 // make a smith_info class
163 auto info = make_shared<const SMITH_Info<complex<double>>>(r, idata);
164 if (method == "caspt2") {
165 algo_ = make_shared<RelCASPT2::RelCASPT2>(info);
166 } else if (method == "casa") {
167 algo_ = make_shared<RelCASA::RelCASA>(info);
168 } else if (method == "mrci") {
169 algo_ = make_shared<RelMRCI::RelMRCI>(info);
170 } else {
171 stringstream ss; ss << method << " method is not implemented in RelSMITH";
172 throw logic_error(ss.str());
173 }
174 } else {
175 #ifndef DISABLE_SERIALIZATION
176 // method == "continue" - so load SMITH_Info and T2 amplitudes from Archives
177 Timer mtimer;
178
179 shared_ptr<const SMITH_Info<complex<double>>> info;
180 string arch_prefix = idata_->get<string>("archive_location", "");
181 if (arch_prefix.size() > 2 && arch_prefix.back() != '/')
182 arch_prefix += "/";
183 {
184 IArchive archive(arch_prefix + "RelSMITH_info");
185 shared_ptr<SMITH_Info<complex<double>>> ptr;
186 archive >> ptr;
187 const int state_begin = idata_->get<int>("state_begin", 0);
188 const int restart_iter = idata_->get<int>("restart_iter", 0);
189 ptr->set_restart_params(state_begin, restart_iter);
190 info = shared_ptr<SMITH_Info<complex<double>>>(ptr);
191 }
192 mtimer.tick_print("Load RelSMITH info Archive");
193
194 ref_ = info->ref();
195 geom_ = ref_->geom();
196 const string method = to_lower(info->method());
197 assert(method == "caspt2" || method == "casa");
198 if (method == "caspt2") {
199 arch_prefix += "RelCASPT2";
200 } else if (method == "casa") {
201 arch_prefix += "RelCASA";
202 }
203
204 if (method == "caspt2")
205 algo_ = make_shared<RelCASPT2::RelCASPT2>(info);
206 else if (method == "casa")
207 algo_ = make_shared<RelCASA::RelCASA>(info);
208 mtimer.tick_print("Construct " + method + " architecture");
209
210 for (int ist = 0; ist <= info->state_begin(); ++ist) {
211 if (ist < info->state_begin() || info->restart_iter() > 0) {
212 string arch = arch_prefix + "_t2_" + to_string(ist);
213 if (ist < info->state_begin())
214 arch += "_converged";
215 else
216 arch += "_iter_" + to_string(info->restart_iter());
217 IArchive archive(arch);
218 shared_ptr<MultiTensor_<complex<double>>> t2in;
219 archive >> t2in;
220 if (method == "caspt2")
221 (dynamic_pointer_cast<RelCASPT2::RelCASPT2>(algo_))->load_t2all(t2in, ist);
222 else
223 (dynamic_pointer_cast<RelCASA::RelCASA>(algo_))->load_t2all(t2in, ist);
224 }
225 }
226 mtimer.tick_print("Load T-amplitude Archive (RelSMITH)");
227 #else
228 throw logic_error("You must enable serialization to make \"continue\" to be available.");
229 #endif
230 }
231 #else
232 throw logic_error("You must enable SMITH during compilation for this method to be available.");
233 #endif
234 }
235