1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: RelCASPT2.cc
4 // Copyright (C) 2014 Toru Shiozaki
5 //
6 // Author: Toru Shiozaki <shiozaki@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 #ifdef COMPILE_SMITH
27 
28 
29 #include <cstdio>
30 #include <src/smith/relcaspt2/RelCASPT2.h>
31 #include <src/util/math/linearRM.h>
32 #include <src/prop/pseudospin/pseudospin.h>
33 
34 using namespace std;
35 using namespace bagel;
36 using namespace bagel::SMITH;
37 
RelCASPT2(shared_ptr<const SMITH_Info<std::complex<double>>> ref)38 RelCASPT2::RelCASPT2::RelCASPT2(shared_ptr<const SMITH_Info<std::complex<double>>> ref) : SpinFreeMethod(ref) {
39   nstates_ = info_->nact() ? ref->ciwfn()->nstates() : 1;
40   auto eig = f1_->diag();
41   eig_.resize(eig.size());
42   for (int i = 0; i != eig.size(); ++i)
43     eig_[i] = real(eig[i]);
44 
45   // MS-CASPT2: t2 and s as MultiTensor (t2all, sall)
46   for (int i = 0; i != nstates_; ++i) {
47     auto tmp = make_shared<MultiTensor>(nstates_);
48     auto tmp2 = make_shared<MultiTensor>(nstates_);
49 
50     for (int j = 0; j != nstates_; ++j)
51       if (!info_->sssr() || i == j) {
52         (*tmp)[j] = init_amplitude();
53         (*tmp2)[j] = init_residual();
54       }
55 
56     t2all_.push_back(tmp);
57     sall_.push_back(tmp2);
58     rall_.push_back(tmp2->clone());
59   }
60 
61   energy_.resize(nstates_);
62   pt2energy_.resize(nstates_);
63 }
64 
65 
solve()66 void RelCASPT2::RelCASPT2::solve() {
67   Timer timer;
68   print_iteration();
69 
70   // <proj_jst|H|0_K> set to sall in ms-caspt2
71   for (int istate = 0; istate != nstates_; ++istate) { //K states
72     if (istate > info_->state_begin() || (istate == info_->state_begin() && info_->restart_iter() == 0))
73       t2all_[istate]->fac(istate) = 0.0;
74     sall_[istate]->fac(istate)  = 0.0;
75 
76     for (int jst=0; jst != nstates_; ++jst) { // <jst|
77       if (info_->sssr() && jst != istate)
78         continue;
79       set_rdm(jst, istate);
80       s = sall_[istate]->at(jst);
81       shared_ptr<Queue> sourceq = make_sourceq(false, jst == istate);
82       while(!sourceq->done())
83         sourceq->next_compute();
84     }
85   }
86 
87   // solve linear equation for t amplitudes
88   t2all_ = solve_linear(sall_, t2all_);
89 
90   // print out energies
91   vector<complex<double>> shift_correction(nstates_*nstates_);
92   cout << endl;
93   for (int istate = 0; istate != nstates_; ++istate) {
94     if (info_->shift() == 0.0) {
95       pt2energy_[istate] = energy_[istate] + std::real((*eref_)(istate,istate));
96       assert(std::abs(std::imag((*eref_)(istate,istate))) < 1.0e-8);
97       cout << "    * RelCASPT2 energy : state " << setw(2) << istate << fixed << setw(20) << setprecision(10) << pt2energy_[istate] <<endl;
98     } else {
99       // will be used in normq
100       n = init_residual();
101       complex<double> norm = 0.0;
102       for (int jstate = 0; jstate != nstates_; ++jstate) {
103         if (istate != jstate && info_->shift_diag())
104           continue;
105         complex<double> nn = 0.0;
106         for (int jst = 0; jst != nstates_; ++jst) { // bra
107           for (int ist = 0; ist != nstates_; ++ist) { // ket
108             if (info_->sssr() && (jst != istate || ist != jstate))
109               continue;
110             set_rdm(jst, ist);
111             t2 = t2all_[jstate]->at(ist);
112             shared_ptr<Queue> normq = make_normq(true, jst == ist);
113             while (!normq->done())
114               normq->next_compute();
115             nn += conj(dot_product_transpose(n, t2all_[istate]->at(jst)));
116           }
117         }
118         shift_correction[istate + nstates_*jstate] = nn;
119         if (jstate == istate) norm = nn;
120       }
121 
122       pt2energy_[istate] = energy_[istate] + std::real((*eref_)(istate,istate)) - info_->shift()*std::real(norm);
123       assert(std::abs(std::imag((*eref_)(istate,istate))) < 1.0e-8);
124       assert(std::abs(std::imag(norm)) < 1.0e-8);
125       cout << "    * RelCASPT2 energy : state " << setw(2) << istate << fixed << setw(20) << setprecision(10) << pt2energy_[istate] << endl;
126       cout << "        w/o shift correction  " << fixed << setw(20) << setprecision(10) << energy_[istate]+(*eref_)(istate,istate) <<endl;
127       cout <<endl;
128     }
129   }
130 
131   // MS-CASPT2
132   if (info_->do_ms() && nstates_ > 1) {
133     heff_ = make_shared<ZMatrix>(nstates_, nstates_);
134 
135     if (info_->shift())
136       cout << "    MS-RelCASPT2:  Applying levelshift correction to " << (info_->shift_diag() ? "diagonal" : "all" ) <<  " elements of the effective Hamiltonian."  << endl << endl;
137 
138     for (int ist = 0; ist != nstates_; ++ist) {
139       auto sist = make_shared<MultiTensor>(nstates_);
140       for (int jst=0; jst != nstates_; ++jst) {
141         if (sall_[ist]->at(jst)) {
142           sist->at(jst) = sall_[ist]->at(jst);
143         } else {
144           set_rdm(jst, ist);
145           s = init_residual();
146           shared_ptr<Queue> sourceq = make_sourceq(false, jst == ist);
147           while(!sourceq->done())
148             sourceq->next_compute();
149           sist->at(jst) = s;
150         }
151       }
152 
153       for (int jst = 0; jst != nstates_; ++jst) {
154         if (ist == jst) {
155           // set diagonal elements
156           (*heff_)(ist, ist) = pt2energy_[ist];
157         } else {
158           // set off-diag elements
159           // 1/2 [ <1g | H | Oe> + <0g |H | 1e > ]
160           (*heff_)(jst, ist) = std::conj(dot_product_transpose(sist, t2all_[jst])) + (*eref_)(jst, ist);
161           if (!info_->shift_diag())
162             (*heff_)(jst, ist) -= shift_correction[jst+nstates_*ist]*info_->shift();
163         }
164       }
165     }
166     heff_->hermite();
167 
168     // print out the effective Hamiltonian
169     cout << endl;
170     cout << "    * MS-RelCASPT2 Heff";
171     for (int ist = 0; ist != nstates_; ++ist) {
172       cout << endl << "      ";
173       for (int jst = 0; jst != nstates_; ++jst)
174         cout << setw(20) << setprecision(10) << (*heff_)(ist, jst);
175     }
176     cout << endl << endl;
177 
178     VectorB eig(nstates_);
179     heff_->diagonalize(eig);
180     copy_n(eig.data(), nstates_, pt2energy_.data());
181 
182     // print out the eigen vector
183     cout << endl;
184     cout << "    * MS-RelCASPT2 rotation matrix";
185     for (int ist = 0; ist != nstates_; ++ist) {
186       cout << endl << "      ";
187       for (int jst = 0; jst != nstates_; ++jst)
188         cout << setw(20) << setprecision(10) << (*heff_)(ist, jst);
189     }
190     cout << endl << endl;
191 
192     // energy printout
193     for (int istate = 0; istate != nstates_; ++istate)
194       cout << "    * MS-RelCASPT2 energy : state " << setw(2) << istate << fixed << setw(20) << setprecision(10) << pt2energy_[istate] << endl;
195     cout << endl << endl;
196   } else {
197     heff_ = make_shared<ZMatrix>(1,1);
198     heff_->element(0,0) = 1.0;
199   }
200   energy_ = pt2energy_;
201 
202   // TODO When the Property class is implemented, this should be one
203   if (info_->aniso_data()) {
204     if (info_->geom()->magnetism()) {
205       cout << "  ** Magnetic anisotropy analysis is currently only available for zero-field calculations; sorry." << endl;
206     } else {
207       shared_ptr<const RelCIWfn> new_ciwfn = info_->do_ms() ? rotate_ciwfn(info_->ciwfn(), *heff_) : info_->ciwfn();
208       const int nspin = info_->aniso_data()->get<int>("nspin", nstates_-1);
209       Pseudospin ps(nspin, info_->geom(), new_ciwfn, info_->aniso_data());
210       ps.compute(energy_, info_->relref()->relcoeff()->block_format()->active_part());
211     }
212   }
213 }
214 
215 
216 // function to solve linear equation
solve_linear(vector<shared_ptr<MultiTensor_<complex<double>>>> s,vector<shared_ptr<MultiTensor_<complex<double>>>> t)217 vector<shared_ptr<MultiTensor_<complex<double>>>> RelCASPT2::RelCASPT2::solve_linear(vector<shared_ptr<MultiTensor_<complex<double>>>> s,
218                                                                                      vector<shared_ptr<MultiTensor_<complex<double>>>> t) {
219   Timer mtimer;
220 #ifndef DISABLE_SERIALIZATION
221   if (info_->restart()) {
222     OArchive archive("RelSMITH_info");
223     archive << info_;
224     mtimer.tick_print("Save RelSMITH info Archive");
225   }
226 #endif
227 
228   // ms-caspt2: R_K = <proj_jst| H0 - E0_K |1_ist> + <proj_jst| H |0_K> is set to rall
229   // loop over state of interest
230   bool converged = true;
231   for (int i = 0; i != nstates_; ++i) {  // K states
232     bool conv = false;
233     double error = 0.0;
234     e0_ = e0all_[i] - info_->shift();
235     energy_[i] = 0.0;
236     // set guess vector
237     if (i > info_->state_begin() || (i == info_->state_begin() && info_->restart_iter() == 0))
238       t[i]->zero();
239     if (s[i]->rms() < 1.0e-15) {
240       print_iteration(0, 0.0, 0.0, mtimer.tick());
241       if (i+1 != nstates_) cout << endl;
242       continue;
243     } else {
244       if (i > info_->state_begin() || (i == info_->state_begin() && info_->restart_iter() == 0))
245         update_amplitude(t[i], s[i]);
246     }
247 
248     auto solver = make_shared<LinearRM<MultiTensor, ZMatrix>>(info_->davidson_subspace(), s[i]);
249     for (int iter = 0; iter != info_->maxiter(); ++iter) {
250       rall_[i]->zero();
251       const double norm = t[i]->norm();
252       t[i]->scale(1.0/norm);
253 
254       // compute residuals named r for each K
255       for (int jst = 0; jst != nstates_; ++jst) { // jst bra vector
256         for (int ist = 0; ist != nstates_; ++ist) { // ist ket vector
257           if (info_->sssr() && (jst != i || ist != i))
258             continue;
259           // first term <proj_jst| H0 - E0_K |1_ist>
260           set_rdm(jst, ist);
261           t2 = t[i]->at(ist);
262           r = rall_[i]->at(jst);
263 
264           // compute residuals named r for each K
265           e0_ = e0all_[i] - info_->shift();
266           shared_ptr<Queue> queue = make_residualq(false, jst == ist);
267           while (!queue->done())
268             queue->next_compute();
269           diagonal(r, t2, jst == ist);
270         }
271       }
272       // solve using subspace updates
273       rall_[i] = solver->compute_residual(t[i], rall_[i]);
274       t[i] = solver->civec();
275       t2all_[i] = t[i];
276 
277       // energy is now the Hylleraas energy
278       energy_[i] = detail::real(dot_product_transpose(s[i], t[i]));
279       energy_[i] += detail::real(dot_product_transpose(rall_[i], t[i]));
280 
281       // compute rms for state i
282       error = rall_[i]->norm() / pow(rall_[i]->size(), 0.25);
283       print_iteration(iter, energy_[i], error, mtimer.tick());
284       conv = error < info_->thresh();
285 
286 #ifndef DISABLE_SERIALIZATION
287       if (info_->restart() && (conv || (info_->restart_each_iter() && iter > 0))) {
288         string arch = "RelCASPT2_";
289         if (mpi__->rank() == 0)
290           arch += "t2_" + to_string(i) + (conv ? "_converged" : "_iter_" + to_string(iter));
291         else
292           arch += "temp_trash";
293         {
294           OArchive archive(arch);
295           archive << t2all_[i];
296         }
297         mpi__->barrier();
298         if (conv)
299           mtimer.tick_print("Save T-amplitude Archive (RelSMITH)");
300         remove("RelCASPT2_temp_trash.archive");
301       }
302 #endif
303 
304       // compute delta t2 and update amplitude
305       if (!conv) {
306         t[i]->zero();
307         update_amplitude(t[i], rall_[i]);
308       }
309       if (conv) break;
310     }
311     if (i+1 != nstates_) cout << endl;
312     converged &= conv;
313   }
314   print_iteration(!converged);
315   return t;
316 }
317 
318 
load_t2all(shared_ptr<MultiTensor> t2in,const int ist)319 void RelCASPT2::RelCASPT2::load_t2all(shared_ptr<MultiTensor> t2in, const int ist) {
320   assert(ist <= info_->state_begin());
321   t2all_[ist] = t2in;
322 }
323 
324 
solve_gradient(const int targetJ,const int targetI,shared_ptr<const NacmType> nacmtype,const bool nocider)325 void RelCASPT2::RelCASPT2::solve_gradient(const int targetJ, const int targetI, shared_ptr<const NacmType> nacmtype, const bool nocider) {
326   throw std::logic_error("Nuclear gradients not implemented for RelCASPT2");
327 }
328 
329 
330 #endif
331