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