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