1 //
2 // mp2r12_energy.cc
3 //
4 // Copyright (C) 2003 Edward Valeev
5 //
6 // Author: Edward Valeev <edward.valeev@chemistry.gatech.edu>
7 // Maintainer: EV
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit 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 Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #ifdef __GNUG__
29 #pragma implementation
30 #endif
31
32 #include <ostream>
33 #include <fstream>
34 #include <sstream>
35 #include <stdexcept>
36 #include <util/misc/string.h>
37 #include <util/misc/formio.h>
38 #include <util/misc/timer.h>
39 #include <util/ref/ref.h>
40 #include <math/scmat/local.h>
41 #include <chemistry/qc/mbpt/bzerofast.h>
42 #include <chemistry/qc/mbptr12/mp2r12_energy.h>
43 #include <chemistry/qc/mbptr12/pairiter.h>
44 #include <chemistry/qc/mbptr12/vxb_eval_info.h>
45 #include <chemistry/qc/mbptr12/svd.h>
46 #include <chemistry/qc/mbptr12/print_scmat_norms.h>
47
48 using namespace std;
49 using namespace sc;
50 using namespace sc::exp;
51
max(int a,int b)52 inline int max(int a,int b) { return (a > b) ? a : b;}
53
54 #define USE_INVERT 0
55
56 /*-------------
57 MP2R12Energy
58 -------------*/
59 static ClassDesc MP2R12Energy_cd(
60 typeid(MP2R12Energy),"MP2R12Energy",1,"virtual public SavableState",
61 0, 0, create<MP2R12Energy>);
62
MP2R12Energy(Ref<R12IntEval> & r12eval,LinearR12::StandardApproximation stdapp,int debug)63 MP2R12Energy::MP2R12Energy(Ref<R12IntEval>& r12eval, LinearR12::StandardApproximation stdapp, int debug)
64 {
65 r12eval_ = r12eval;
66 stdapprox_ = stdapp;
67 if (debug >= 0)
68 debug_ = debug;
69 else
70 debug_ = 0;
71 evaluated_ = false;
72
73 init_();
74 }
75
76 void
init_()77 MP2R12Energy::init_()
78 {
79
80 RefSCDimension dim_oo_aa = r12eval_->dim_oo_aa();
81 RefSCDimension dim_oo_ab = r12eval_->dim_oo_ab();
82 Ref<SCMatrixKit> kit = r12eval_->r12info()->matrixkit();
83 er12_aa_ = kit->vector(dim_oo_aa);
84 er12_ab_ = kit->vector(dim_oo_ab);
85 emp2r12_aa_ = kit->vector(dim_oo_aa);
86 emp2r12_ab_ = kit->vector(dim_oo_ab);
87
88 RefSCDimension dim_vv_aa = r12eval_->dim_vv_aa();
89 RefSCDimension dim_vv_ab = r12eval_->dim_vv_ab();
90 Caa_ = kit->matrix(dim_oo_aa, dim_oo_aa);
91 Cab_ = kit->matrix(dim_oo_ab, dim_oo_ab);
92
93 }
94
MP2R12Energy(StateIn & si)95 MP2R12Energy::MP2R12Energy(StateIn& si) : SavableState(si)
96 {
97 r12eval_ << SavableState::restore_state(si);
98
99 init_();
100
101 er12_aa_.restore(si);
102 er12_ab_.restore(si);
103 emp2r12_aa_.restore(si);
104 emp2r12_ab_.restore(si);
105
106 Caa_.restore(si);
107 Cab_.restore(si);
108
109 int stdapprox;
110 si.get(stdapprox);
111 stdapprox_ = (LinearR12::StandardApproximation) stdapprox;
112 si.get(debug_);
113 int evaluated;
114 si.get(evaluated);
115 evaluated_ = (bool) evaluated;
116 }
117
~MP2R12Energy()118 MP2R12Energy::~MP2R12Energy()
119 {
120 r12eval_ = 0;
121 }
122
save_data_state(StateOut & so)123 void MP2R12Energy::save_data_state(StateOut& so)
124 {
125 SavableState::save_state(r12eval_.pointer(),so);
126
127 er12_aa_.save(so);
128 er12_ab_.save(so);
129 emp2r12_aa_.save(so);
130 emp2r12_ab_.save(so);
131
132 Caa_.save(so);
133 Cab_.save(so);
134
135 so.put((int)stdapprox_);
136 so.put(debug_);
137 so.put((int)evaluated_);
138 }
139
obsolete()140 void MP2R12Energy::obsolete()
141 {
142 evaluated_ = false;
143 }
144
r12eval() const145 Ref<R12IntEval> MP2R12Energy::r12eval() const { return r12eval_; };
ebc() const146 bool MP2R12Energy::ebc() const { return ebc_; };
gbc() const147 bool MP2R12Energy::gbc() const { return r12eval_->gbc(); };
stdapp() const148 LinearR12::StandardApproximation MP2R12Energy::stdapp() const { return stdapprox_; };
set_debug(int debug)149 void MP2R12Energy::set_debug(int debug) { debug_ = debug; };
get_debug() const150 int MP2R12Energy::get_debug() const { return debug_; };
151
energy()152 double MP2R12Energy::energy()
153 {
154 double value = emp2tot_aa_() + emp2tot_ab_() + er12tot_aa_() + er12tot_ab_();
155 return value;
156 }
157
emp2tot_aa_() const158 double MP2R12Energy::emp2tot_aa_() const
159 {
160 RefSCVector emp2_aa = r12eval_->emp2_aa();
161 int nij = emp2_aa.dim().n();
162 double value = 0;
163 for(int ij=0; ij<nij; ij++)
164 value += emp2_aa.get_element(ij);
165
166 return value;
167 }
168
emp2tot_ab_() const169 double MP2R12Energy::emp2tot_ab_() const
170 {
171 RefSCVector emp2_ab = r12eval_->emp2_ab();
172 int nij = emp2_ab.dim().n();
173 double value = 0;
174 for(int ij=0; ij<nij; ij++)
175 value += emp2_ab.get_element(ij);
176
177 return value;
178 }
179
er12tot_aa_()180 double MP2R12Energy::er12tot_aa_()
181 {
182 compute();
183 int nij = er12_aa_.dim().n();
184 double value = 0;
185 for(int ij=0; ij<nij; ij++)
186 value += er12_aa_.get_element(ij);
187
188 return value;
189 }
190
er12tot_ab_()191 double MP2R12Energy::er12tot_ab_()
192 {
193 compute();
194 int nij = er12_ab_.dim().n();
195 double value = 0;
196 for(int ij=0; ij<nij; ij++)
197 value += er12_ab_.get_element(ij);
198
199 return value;
200 }
201
compute()202 void MP2R12Energy::compute()
203 {
204 if (evaluated_)
205 return;
206
207 Ref<R12IntEvalInfo> r12info = r12eval_->r12info();
208 Ref<MessageGrp> msg = r12info->msg();
209 const int me = msg->me();
210 const int ntasks = msg->n();
211
212 const bool ebc = r12eval_->ebc();
213 const bool follow_ks_ebcfree = r12eval_->follow_ks_ebcfree();
214
215 //
216 // Evaluate pair energies:
217 // distribute workload among nodes by pair index
218 //
219
220 // Need eigenvalues
221 const int nocc = r12info->nocc();
222 const int nfzc = r12info->nfzc();
223 const int nocc_act = r12info->nocc_act();
224 const int nvir_act = r12info->nvir_act();
225 RefDiagSCMatrix evalmat = r12eval_->evals();
226 vector<double> evals_act_occ(nocc_act);
227 vector<double> evals_act_vir(nvir_act);
228 for(int i=nfzc; i<nocc; i++)
229 evals_act_occ[i-nfzc] = evalmat(i);
230 for(int i=0; i<nvir_act; i++)
231 evals_act_vir[i] = evalmat(i+nocc);
232 evalmat = 0;
233
234 // Get the intermediates
235 RefSCMatrix Vaa = r12eval_->V_aa();
236 RefSCMatrix Xaa = r12eval_->X_aa();
237 RefSymmSCMatrix Baa = r12eval_->B_aa();
238 RefSCMatrix Aaa = r12eval_->A_aa();
239 RefSCMatrix Vab = r12eval_->V_ab();
240 RefSCMatrix Xab = r12eval_->X_ab();
241 RefSymmSCMatrix Bab = r12eval_->B_ab();
242 RefSCMatrix Aab = r12eval_->A_ab();
243 RefSCMatrix Ac_aa, Ac_ab;
244 if (follow_ks_ebcfree) {
245 Ac_aa = r12eval_->Ac_aa();
246 Ac_ab = r12eval_->Ac_ab();
247 }
248 RefSCVector emp2_aa = r12eval_->emp2_aa();
249 RefSCVector emp2_ab = r12eval_->emp2_ab();
250
251 // Prepare total and R12 pairs
252 Ref<SCMatrixKit> localkit = Vaa.kit();
253 RefSCDimension dim_oo_aa = r12eval_->dim_oo_aa();
254 RefSCDimension dim_oo_ab = r12eval_->dim_oo_ab();
255 int naa = dim_oo_aa.n();
256 int nab = dim_oo_ab.n();
257 emp2r12_aa_ = localkit->vector(dim_oo_aa);
258 emp2r12_ab_ = localkit->vector(dim_oo_ab);
259 er12_aa_ = localkit->vector(dim_oo_aa);
260 er12_ab_ = localkit->vector(dim_oo_ab);
261 double* er12_aa_vec = new double[naa];
262 double* er12_ab_vec = new double[nab];
263 bzerofast(er12_aa_vec,naa);
264 bzerofast(er12_ab_vec,nab);
265
266 //
267 // Alpha-alpha pairs
268 //
269 if (naa > 0) {
270 if (debug_ > 0) {
271 print_scmat_norms(Vaa,"Alpha-alpha V matrix");
272 print_scmat_norms(Baa,"Alpha-alpha MP2-R12/A B matrix");
273 if (ebc == false)
274 print_scmat_norms(Aaa,"Alpha-alpha A matrix");
275 }
276 if (debug_ > 1) {
277 Vaa.print("Alpha-alpha V matrix");
278 Baa.print("Alpha-alpha MP2-R12/A B matrix");
279 if (ebc == false)
280 Aaa.print("Alpha-alpha A matrix");
281 }
282
283 // Allocate the B matrix:
284 // 1) in MP2-R12/A the B matrix is the same for all pairs
285 // 2) int MP2-R12/A' the B matrix is pair-specific
286 RefSymmSCMatrix Baa_ij = Baa.clone();
287 if (stdapprox_ == LinearR12::StdApprox_A) {
288 #if USE_INVERT
289 Baa_ij->assign(Baa);
290 Baa_ij->gen_invert_this();
291 if (debug_ > 0)
292 print_scmat_norms(Baa_ij,"Inverse alpha-alpha MP2-R12/A B matrix");
293 if (debug_ > 1)
294 Baa_ij.print("Inverse alpha-alpha MP2-R12/A B matrix");
295 #else
296 // solve B * C = V
297 RefSCMatrix Caa_kl_by_ij = Caa_.clone();
298 sc::exp::lapack_linsolv_symmnondef(Baa, Caa_kl_by_ij, Vaa);
299 Caa_kl_by_ij = Caa_kl_by_ij.t();
300 Caa_.assign(Caa_kl_by_ij); Caa_kl_by_ij = 0;
301 Caa_.scale(-1.0);
302 #endif
303 }
304
305 int ij=0;
306 for(int i=0; i<nocc_act; i++)
307 for(int j=0; j<i; j++, ij++) {
308
309 if (ij%ntasks != me)
310 continue;
311
312 RefSCVector Vaa_ij = Vaa.get_column(ij);
313
314 // In MP2-R12/A' matrices B are pair-specific:
315 // Form B(ij)kl,ow = Bkl,ow + 1/2(ek + el + eo + ew - 2ei - 2ej)Xkl,ow
316 if (stdapprox_ == LinearR12::StdApprox_Ap) {
317 Baa_ij.assign(Baa);
318 int kl=0;
319 for(int k=0; k<nocc_act; k++)
320 for(int l=0; l<k; l++, kl++) {
321 int ow=0;
322 for(int o=0; o<nocc_act; o++)
323 for(int w=0; w<o; w++, ow++) {
324
325 if (ow > kl)
326 continue;
327
328 double fx = 0.5 * (evals_act_occ[k] + evals_act_occ[l] + evals_act_occ[o] + evals_act_occ[w]
329 - 2.0*evals_act_occ[i] - 2.0*evals_act_occ[j]) *
330 Xaa.get_element(kl,ow);
331
332 Baa_ij.accumulate_element(kl,ow,fx);
333
334 // If EBC is not assumed add Akl,cd*Acd,ow/(ec+ed-ei-ej)
335 if (ebc == false) {
336 double fy = 0.0;
337 int cd=0;
338 if (follow_ks_ebcfree) {
339 for(int c=0; c<nvir_act; c++)
340 for(int d=0; d<c; d++, cd++) {
341 fy -= 0.5 * (Aaa.get_element(kl,cd)*Ac_aa.get_element(ow,cd) + Ac_aa.get_element(kl,cd)*Aaa.get_element(ow,cd))/(evals_act_vir[c] + evals_act_vir[d]
342 - evals_act_occ[i] - evals_act_occ[j]);
343 }
344 }
345 else {
346 for(int c=0; c<nvir_act; c++)
347 for(int d=0; d<c; d++, cd++) {
348 fy -= Aaa.get_element(kl,cd)*Aaa.get_element(ow,cd)/(evals_act_vir[c] + evals_act_vir[d]
349 - evals_act_occ[i] - evals_act_occ[j]);
350 }
351 }
352
353 Baa_ij.accumulate_element(kl,ow,fy);
354 }
355
356 }
357 }
358
359 std::ostringstream oss;
360 oss << "Alpha-alpha MP2-R12/A' B(ij=" << i << "," << j << ") matrix";
361 std::string label(oss.str());
362
363 if (debug_ > 0)
364 print_scmat_norms(Baa_ij,label);
365 if (debug_ > 1)
366 Baa_ij.print(label.c_str());
367
368 #if USE_INVERT
369 Baa_ij->gen_invert_this();
370 std::string invlabel("Inverse "); invlabel += label;
371 if (debug_ > 0)
372 print_scmat_norms(Baa_ij,invlavel);
373 if (debug_ > 1)
374 Baa_ij.print(invlabel.c_str());
375 #endif
376
377 }
378
379 #if USE_INVERT
380 // The r12 amplitudes B^-1 * V
381 RefSCVector Cij = -1.0*(Baa_ij * Vaa_ij);
382 const int nkl = Cij.dim().n();
383 for(int kl=0; kl<nkl; kl++)
384 Caa_.set_element(ij,kl,Cij.get_element(kl));
385 #else
386 RefSCVector Cij = Vaa_ij.clone();
387 if (stdapprox_ == LinearR12::StdApprox_A) {
388 double* v = new double[Cij.n()];
389 Caa_.get_row(ij).convert(v);
390 Cij.assign(v);
391 delete[] v;
392 }
393 else {
394 // solve B * C = V
395 Cij = Vaa_ij.clone();
396 sc::exp::lapack_linsolv_symmnondef(Baa_ij, Cij, Vaa_ij);
397 Cij.scale(-1.0);
398 const int nkl = Cij.dim().n();
399 for(int kl=0; kl<nkl; kl++)
400 Caa_.set_element(ij,kl,Cij.get_element(kl));
401 }
402 #endif
403 double eaa_ij = 2.0*Vaa_ij.dot(Cij);
404 er12_aa_vec[ij] = eaa_ij;
405 }
406 Baa_ij = 0;
407 msg->sum(er12_aa_vec,naa,0,-1);
408 er12_aa_->assign(er12_aa_vec);
409 emp2r12_aa_->assign(emp2_aa);
410 emp2r12_aa_->accumulate(er12_aa_);
411 delete[] er12_aa_vec;
412 }
413 if (debug_ > 0)
414 print_scmat_norms(Caa_,"Alpha-alpha R12 amplitudes");
415
416 //
417 // Alpha-beta pairs
418 //
419 if (nab > 0) {
420 if (debug_ > 0) {
421 print_scmat_norms(Vab,"Alpha-beta V matrix");
422 print_scmat_norms(Bab,"Alpha-beta MP2-R12/A B matrix");
423 if (ebc == false)
424 print_scmat_norms(Aab,"Alpha-beta A matrix");
425 }
426 if (debug_ > 1) {
427 Vab.print("Alpha-beta V matrix");
428 Bab.print("Alpha-beta MP2-R12/A B matrix");
429 if (ebc == false)
430 Aab.print("Alpha-beta A matrix");
431 }
432
433 RefSymmSCMatrix Bab_ij = Bab.clone();
434 // In MP2-R12/A the B matrix is the same for all pairs
435 if (stdapprox_ == LinearR12::StdApprox_A) {
436 #if USE_INVERT
437 Bab_ij.assign(Bab);
438 Bab_ij->gen_invert_this();
439 if (debug_ > 0)
440 print_scmat_norms(Bab_ij,"Inverse alpha-beta MP2-R12/A B matrix");
441 if (debug_ > 1)
442 Bab_ij.print("Inverse alpha-beta MP2-R12/A B matrix");
443 #else
444 // solve B * C = V
445 RefSCMatrix Cab_kl_by_ij = Cab_.clone();
446 sc::exp::lapack_linsolv_symmnondef(Bab, Cab_kl_by_ij, Vab);
447 Cab_kl_by_ij = Cab_kl_by_ij.t();
448 Cab_.assign(Cab_kl_by_ij); Cab_kl_by_ij = 0;
449 Cab_.scale(-1.0);
450 #endif
451 }
452
453 int ij=0;
454 for(int i=0; i<nocc_act; i++)
455 for(int j=0; j<nocc_act; j++, ij++) {
456
457 if (ij%ntasks != me)
458 continue;
459
460 RefSCVector Vab_ij = Vab.get_column(ij);
461
462 // In MP2-R12/A' matrices B are pair-specific:
463 // Form B(ij)kl,ow = Bkl,ow + 1/2(ek + el + eo + ew - 2ei - 2ej)Xkl,ow
464 if (stdapprox_ == LinearR12::StdApprox_Ap) {
465 Bab_ij.assign(Bab);
466 int kl=0;
467 for(int k=0; k<nocc_act; k++)
468 for(int l=0; l<nocc_act; l++, kl++) {
469 int ow=0;
470 for(int o=0; o<nocc_act; o++)
471 for(int w=0; w<nocc_act; w++, ow++) {
472
473 if (ow > kl)
474 continue;
475
476 double fx = 0.5 * (evals_act_occ[k] + evals_act_occ[l] + evals_act_occ[o] + evals_act_occ[w]
477 - 2.0*evals_act_occ[i] - 2.0*evals_act_occ[j]) *
478 Xab.get_element(kl,ow);
479 Bab_ij.accumulate_element(kl,ow,fx);
480
481 // If EBC is not assumed add Akl,cd*Acd,ow/(ec+ed-ei-ej)
482 if (ebc == false) {
483 double fy = 0.0;
484 int cd=0;
485 if (follow_ks_ebcfree) {
486 for(int c=0; c<nvir_act; c++)
487 for(int d=0; d<nvir_act; d++, cd++) {
488 fy -= 0.5 * (Aab.get_element(kl,cd)*Ac_ab.get_element(ow,cd) + Ac_ab.get_element(kl,cd)*Aab.get_element(ow,cd))/(evals_act_vir[c] + evals_act_vir[d]
489 - evals_act_occ[i] - evals_act_occ[j]);
490 }
491 }
492 else {
493 for(int c=0; c<nvir_act; c++)
494 for(int d=0; d<nvir_act; d++, cd++) {
495 fy -= Aab.get_element(kl,cd)*Aab.get_element(ow,cd)/(evals_act_vir[c] + evals_act_vir[d]
496 - evals_act_occ[i] - evals_act_occ[j]);
497 }
498 }
499
500 Bab_ij.accumulate_element(kl,ow,fy);
501 }
502
503 }
504 }
505
506 std::ostringstream oss;
507 oss << "Alpha-beta MP2-R12/A' B(ij=" << i << "," << j << ") matrix";
508 std::string label(oss.str());
509
510 if (debug_ > 0)
511 print_scmat_norms(Bab_ij,label);
512 if (debug_ > 1)
513 Bab_ij.print(label.c_str());
514
515 #if USE_INVERT
516 Bab_ij->gen_invert_this();
517 std::string invlabel("Inverse "); invlabel += label;
518 if (debug_ > 0)
519 print_scmat_norms(Bab_ij,invlavel);
520 if (debug_ > 1)
521 Bab_ij.print(invlabel.c_str());
522 #endif
523
524 }
525
526 #if USE_INVERT
527 // the r12 amplitudes B^-1 * V
528 RefSCVector Cij = -1.0*(Bab_ij * Vab_ij);
529 const int nkl = Cij.dim().n();
530 for(int kl=0; kl<nkl; kl++)
531 Cab_.set_element(ij,kl,Cij.get_element(kl));
532 #else
533 RefSCVector Cij = Vab_ij.clone();
534 if (stdapprox_ == LinearR12::StdApprox_A) {
535 double* v = new double[Cij.n()];
536 Cab_.get_row(ij).convert(v);
537 Cij.assign(v);
538 delete[] v;
539 }
540 else {
541 // solve B * C = V
542 Cij = Vab_ij.clone();
543 sc::exp::lapack_linsolv_symmnondef(Bab_ij, Cij, Vab_ij);
544 Cij.scale(-1.0);
545 const int nkl = Cij.dim().n();
546 for(int kl=0; kl<nkl; kl++)
547 Cab_.set_element(ij,kl,Cij.get_element(kl));
548 }
549 #endif
550 double eab_ij = 1.0*Vab_ij.dot(Cij);
551 er12_ab_vec[ij] = eab_ij;
552 }
553 Bab_ij=0;
554 msg->sum(er12_ab_vec,nab,0,-1);
555 er12_ab_->assign(er12_ab_vec);
556 emp2r12_ab_->assign(emp2_ab);
557 emp2r12_ab_->accumulate(er12_ab_);
558 delete[] er12_ab_vec;
559 }
560 if (debug_ > 0)
561 print_scmat_norms(Cab_,"Alpha-beta R12 amplitudes");
562
563 evaluated_ = true;
564
565 return;
566 }
567
print_psi_values(std::ostream & fout,const SCVector3 & r1,const SCVector3 & r2,double phi_0,double phi_1_mp2,double phi_1_r12)568 static void print_psi_values(std::ostream& fout, const SCVector3& r1, const SCVector3& r2, double phi_0, double phi_1_mp2, double phi_1_r12)
569 {
570 fout << scprintf("%9.5lf %9.5lf %9.5lf %9.5lf %9.5lf %9.5lf %12.8lf %12.8lf %12.8lf",
571 r1.x(),r1.y(),r1.z(),r2.x(),r2.y(),r2.z(),phi_0,phi_1_mp2,phi_1_r12) << endl;
572 }
573
574 double
compute_pair_function_aa(int ij,const SCVector3 & r1,const SCVector3 & r2)575 MP2R12Energy::compute_pair_function_aa(int ij, const SCVector3& r1, const SCVector3& r2)
576 {
577 Ref<R12Amplitudes> Amps = r12eval_->amps();
578 RefSCMatrix T2aa = Amps->T2_aa();
579 RefSCMatrix Rvv_aa = Amps->Rvv_aa();
580 RefSCMatrix Roo_aa = Amps->Roo_aa();
581 RefSCMatrix Rvo_aa = Amps->Rvo_aa();
582 RefSCMatrix Rxo_aa = Amps->Rxo_aa();
583
584 Ref<SCMatrixKit> localkit = new LocalSCMatrixKit;
585 RefSCMatrix Caa = localkit->matrix(Caa_.rowdim(),Caa_.coldim());
586 double* caa = new double[Caa_.rowdim().n()*Caa_.coldim().n()];
587 Caa_.convert(caa);
588 Caa.assign(caa);
589 delete[] caa;
590 RefSCMatrix Cvv = Caa * Rvv_aa;
591 RefSCMatrix Coo = Caa * Roo_aa;
592 RefSCMatrix Cov = Caa * Rvo_aa;
593 RefSCMatrix Cox = Caa * Rxo_aa;
594
595 Ref<R12IntEvalInfo> r12info = r12eval_->r12info();
596 Ref<MOIndexSpace> act_vir_space = r12info->act_vir_space();
597 Ref<MOIndexSpace> act_occ_space = r12info->act_occ_space();
598 Ref<MOIndexSpace> occ_space = r12info->occ_space();
599 Ref<MOIndexSpace> ribs_space = r12info->ribs_space();
600
601 RefSCVector phi_aa = compute_2body_values_(true,act_occ_space,act_occ_space,r1,r2);
602 RefSCVector phi_vv = compute_2body_values_(true,act_vir_space,act_vir_space,r1,r2);
603 RefSCVector phi_oo = compute_2body_values_(true,occ_space,occ_space,r1,r2);
604 RefSCVector phi_ov = compute_2body_values_(true,occ_space,act_vir_space,r1,r2);
605 RefSCVector phi_ox = compute_2body_values_(true,occ_space,ribs_space,r1,r2);
606
607 double phi_t2 = T2aa.get_row(ij).dot(phi_vv);
608
609 SCVector3 r12 = r1 - r2;
610 const double dist12 = r12.norm();
611 double phi_r12;
612 phi_r12 = 0.5 * Caa.get_row(ij).dot(phi_aa) * dist12;
613 phi_r12 -= 0.5 * Cvv.get_row(ij).dot(phi_vv);
614 phi_r12 -= 0.5 * Coo.get_row(ij).dot(phi_oo);
615 phi_r12 -= 1.0 * Cov.get_row(ij).dot(phi_ov);
616 phi_r12 -= 1.0 * Cox.get_row(ij).dot(phi_ox);
617
618 print_psi_values(ExEnv::out0(),r1,r2,phi_aa.get_element(ij),phi_t2,phi_r12);
619
620 return phi_t2 + phi_r12;
621 }
622
623 void
compute_pair_function_aa(int ij,const Ref<TwoBodyGrid> & tbgrid)624 MP2R12Energy::compute_pair_function_aa(int ij, const Ref<TwoBodyGrid>& tbgrid)
625 {
626 Ref<R12Amplitudes> Amps = r12eval_->amps();
627 RefSCMatrix T2aa = Amps->T2_aa();
628 RefSCMatrix Rvv_aa = Amps->Rvv_aa();
629 RefSCMatrix Roo_aa = Amps->Roo_aa();
630 RefSCMatrix Rvo_aa = Amps->Rvo_aa();
631 RefSCMatrix Rxo_aa = Amps->Rxo_aa();
632
633 Ref<SCMatrixKit> localkit = new LocalSCMatrixKit;
634 RefSCMatrix Caa = localkit->matrix(Caa_.rowdim(),Caa_.coldim());
635 double* caa = new double[Caa_.rowdim().n()*Caa_.coldim().n()];
636 Caa_.convert(caa);
637 Caa.assign(caa);
638 delete[] caa;
639 RefSCMatrix Cvv = Caa * Rvv_aa;
640 RefSCMatrix Coo = Caa * Roo_aa;
641 RefSCMatrix Cov = Caa * Rvo_aa;
642 RefSCMatrix Cox = Caa * Rxo_aa;
643
644 Ref<R12IntEvalInfo> r12info = r12eval_->r12info();
645 Ref<MOIndexSpace> act_vir_space = r12info->act_vir_space();
646 Ref<MOIndexSpace> act_occ_space = r12info->act_occ_space();
647 Ref<MOIndexSpace> occ_space = r12info->occ_space();
648 Ref<MOIndexSpace> ribs_space = r12info->ribs_space();
649
650 const int nelem = tbgrid->nelem();
651 std::stringstream output_file_name;
652 output_file_name << tbgrid->name() << ".ab.pair"
653 << ij << ".txt";
654 ofstream ofile(output_file_name.str().c_str());
655
656 for(int i=0; i<nelem; i++) {
657 RefSCVector phi_aa = compute_2body_values_(true,act_occ_space,act_occ_space,tbgrid->xyz1(i),tbgrid->xyz2(i));
658 RefSCVector phi_vv = compute_2body_values_(true,act_vir_space,act_vir_space,tbgrid->xyz1(i),tbgrid->xyz2(i));
659 RefSCVector phi_oo = compute_2body_values_(true,occ_space,occ_space,tbgrid->xyz1(i),tbgrid->xyz2(i));
660 RefSCVector phi_ov = compute_2body_values_(true,occ_space,act_vir_space,tbgrid->xyz1(i),tbgrid->xyz2(i));
661 RefSCVector phi_ox = compute_2body_values_(true,occ_space,ribs_space,tbgrid->xyz1(i),tbgrid->xyz2(i));
662
663 double phi_t2 = T2aa.get_row(ij).dot(phi_vv);
664
665 SCVector3 r12 = tbgrid->xyz1(i) - tbgrid->xyz2(i);
666 const double dist12 = r12.norm();
667 double phi_r12;
668 phi_r12 = 0.5 * Caa.get_row(ij).dot(phi_aa) * dist12;
669 phi_r12 -= 0.5 * Cvv.get_row(ij).dot(phi_vv);
670 phi_r12 -= 0.5 * Coo.get_row(ij).dot(phi_oo);
671 phi_r12 -= 1.0 * Cov.get_row(ij).dot(phi_ov);
672 phi_r12 -= 1.0 * Cox.get_row(ij).dot(phi_ox);
673
674 print_psi_values(ofile,tbgrid->xyz1(i),tbgrid->xyz2(i),phi_aa.get_element(ij),phi_t2,phi_r12);
675 }
676 }
677
678 double
compute_pair_function_ab(int ij,const SCVector3 & r1,const SCVector3 & r2)679 MP2R12Energy::compute_pair_function_ab(int ij, const SCVector3& r1, const SCVector3& r2)
680 {
681 Ref<R12Amplitudes> Amps = r12eval_->amps();
682 RefSCMatrix T2ab = Amps->T2_ab();
683 RefSCMatrix Rvv_ab = Amps->Rvv_ab();
684 RefSCMatrix Roo_ab = Amps->Roo_ab();
685 RefSCMatrix Rvo_ab = Amps->Rvo_ab();
686 RefSCMatrix Rxo_ab = Amps->Rxo_ab();
687
688 Ref<SCMatrixKit> localkit = new LocalSCMatrixKit;
689 RefSCMatrix Cab = localkit->matrix(Cab_.rowdim(),Cab_.coldim());
690 double* cab = new double[Cab_.rowdim().n()*Cab_.coldim().n()];
691 Cab_.convert(cab);
692 Cab.assign(cab);
693 delete[] cab;
694 RefSCMatrix Cvv = Cab * Rvv_ab;
695 RefSCMatrix Coo = Cab * Roo_ab;
696 RefSCMatrix Cov = Cab * Rvo_ab;
697 RefSCMatrix Cox = Cab * Rxo_ab;
698
699 Ref<R12IntEvalInfo> r12info = r12eval_->r12info();
700 Ref<MOIndexSpace> act_vir_space = r12info->act_vir_space();
701 Ref<MOIndexSpace> act_occ_space = r12info->act_occ_space();
702 Ref<MOIndexSpace> occ_space = r12info->occ_space();
703 Ref<MOIndexSpace> ribs_space = r12info->ribs_space();
704
705 RefSCVector phi_aa = compute_2body_values_(false,act_occ_space,act_occ_space,r1,r2);
706 RefSCVector phi_vv = compute_2body_values_(false,act_vir_space,act_vir_space,r1,r2);
707 RefSCVector phi_oo = compute_2body_values_(false,occ_space,occ_space,r1,r2);
708 RefSCVector phi_ov = compute_2body_values_(false,occ_space,act_vir_space,r1,r2);
709 RefSCVector phi_ox = compute_2body_values_(false,occ_space,ribs_space,r1,r2);
710
711 double phi_t2 = T2ab.get_row(ij).dot(phi_vv);
712
713 SCVector3 r12 = r1 - r2;
714 const double dist12 = r12.norm();
715 double phi_r12;
716 phi_r12 = 0.5*Cab.get_row(ij).dot(phi_aa) * dist12;
717 phi_r12 -= 0.5 * Cvv.get_row(ij).dot(phi_vv);
718 phi_r12 -= 0.5 * Coo.get_row(ij).dot(phi_oo);
719 phi_r12 -= 1.0 * Cov.get_row(ij).dot(phi_ov);
720 phi_r12 -= 1.0 * Cox.get_row(ij).dot(phi_ox);
721
722 print_psi_values(ExEnv::out0(),r1,r2,phi_aa.get_element(ij),phi_t2,phi_r12);
723
724 return phi_t2 + phi_r12;
725 }
726
727 void
compute_pair_function_ab(int ij,const Ref<TwoBodyGrid> & tbgrid)728 MP2R12Energy::compute_pair_function_ab(int ij, const Ref<TwoBodyGrid>& tbgrid)
729 {
730 Ref<R12Amplitudes> Amps = r12eval_->amps();
731 RefSCMatrix T2ab = Amps->T2_ab();
732 RefSCMatrix Rvv_ab = Amps->Rvv_ab();
733 RefSCMatrix Roo_ab = Amps->Roo_ab();
734 RefSCMatrix Rvo_ab = Amps->Rvo_ab();
735 RefSCMatrix Rxo_ab = Amps->Rxo_ab();
736
737 Ref<SCMatrixKit> localkit = new LocalSCMatrixKit;
738 RefSCMatrix Cab = localkit->matrix(Cab_.rowdim(),Cab_.coldim());
739 double* cab = new double[Cab_.rowdim().n()*Cab_.coldim().n()];
740 Cab_.convert(cab);
741 Cab.assign(cab);
742 delete[] cab;
743 RefSCMatrix Cvv = Cab * Rvv_ab;
744 RefSCMatrix Coo = Cab * Roo_ab;
745 RefSCMatrix Cov = Cab * Rvo_ab;
746 RefSCMatrix Cox = Cab * Rxo_ab;
747
748 Ref<R12IntEvalInfo> r12info = r12eval_->r12info();
749 Ref<MOIndexSpace> act_vir_space = r12info->act_vir_space();
750 Ref<MOIndexSpace> act_occ_space = r12info->act_occ_space();
751 Ref<MOIndexSpace> occ_space = r12info->occ_space();
752 Ref<MOIndexSpace> ribs_space = r12info->ribs_space();
753
754 const int nelem = tbgrid->nelem();
755 std::stringstream output_file_name;
756 output_file_name << tbgrid->name() << ".ab.pair"
757 << ij << ".txt";
758 ofstream ofile(output_file_name.str().c_str());
759
760 for(int i=0; i<nelem; i++) {
761 RefSCVector phi_aa = compute_2body_values_(false,act_occ_space,act_occ_space,tbgrid->xyz1(i),tbgrid->xyz2(i));
762 RefSCVector phi_vv = compute_2body_values_(false,act_vir_space,act_vir_space,tbgrid->xyz1(i),tbgrid->xyz2(i));
763 RefSCVector phi_oo = compute_2body_values_(false,occ_space,occ_space,tbgrid->xyz1(i),tbgrid->xyz2(i));
764 RefSCVector phi_ov = compute_2body_values_(false,occ_space,act_vir_space,tbgrid->xyz1(i),tbgrid->xyz2(i));
765 RefSCVector phi_ox = compute_2body_values_(false,occ_space,ribs_space,tbgrid->xyz1(i),tbgrid->xyz2(i));
766
767 double phi_t2 = T2ab.get_row(ij).dot(phi_vv);
768
769 SCVector3 r12 = tbgrid->xyz1(i) - tbgrid->xyz2(i);
770 const double dist12 = r12.norm();
771 double phi_r12;
772 phi_r12 = 0.5*Cab.get_row(ij).dot(phi_aa) * dist12;
773 phi_r12 -= 0.5 * Cvv.get_row(ij).dot(phi_vv);
774 phi_r12 -= 0.5 * Coo.get_row(ij).dot(phi_oo);
775 phi_r12 -= 1.0 * Cov.get_row(ij).dot(phi_ov);
776 phi_r12 -= 1.0 * Cox.get_row(ij).dot(phi_ox);
777
778 print_psi_values(ofile,tbgrid->xyz1(i),tbgrid->xyz2(i),phi_aa.get_element(ij),phi_t2,phi_r12);
779 }
780 }
781
782 RefSCVector
compute_2body_values_(bool equiv,const Ref<MOIndexSpace> & space1,const Ref<MOIndexSpace> & space2,const SCVector3 & r1,const SCVector3 & r2) const783 MP2R12Energy::compute_2body_values_(bool equiv, const Ref<MOIndexSpace>& space1, const Ref<MOIndexSpace>& space2,
784 const SCVector3& r1, const SCVector3& r2) const
785 {
786 const Ref<Integral> ints = r12eval_->r12info()->integral();
787 const Ref<GaussianBasisSet> bs1 = space1->basis();
788 const Ref<GaussianBasisSet> bs2 = space2->basis();
789 ints->set_basis(bs1,bs2);
790 GaussianBasisSet::ValueData* vdata1 = new GaussianBasisSet::ValueData(bs1,ints);
791 GaussianBasisSet::ValueData* vdata2 = new GaussianBasisSet::ValueData(bs2,ints);
792
793 const bool space1_eq_space2 = (space1 == space2);
794 const int nbasis1 = bs1->nbasis();
795 const int nbasis2 = bs2->nbasis();
796 const int rank1 = space1->rank();
797 const int rank2 = space2->rank();
798
799 const int npair = (space1_eq_space2 && equiv) ? rank1*(rank1-1)/2 : rank1*rank2;
800 RefSCDimension pairdim = new SCDimension(npair);
801
802 double* values11 = new double[nbasis1];
803 double* values12 = new double[nbasis1];
804 double* values21 = new double[nbasis2];
805 double* values22 = new double[nbasis2];
806
807 bs1->values(r1,vdata1,values11);
808 bs1->values(r2,vdata1,values12);
809 bs2->values(r1,vdata2,values21);
810 bs2->values(r2,vdata2,values22);
811
812 RefSCMatrix ao2mo_1 = space1->coefs().t();
813 RefSCMatrix ao2mo_2 = space2->coefs().t();
814
815 Ref<SCMatrixKit> kit = ao2mo_1.kit();
816 RefSCVector vals11 = kit->vector(ao2mo_1.coldim());
817 RefSCVector vals12 = kit->vector(ao2mo_1.coldim());
818 RefSCVector vals21 = kit->vector(ao2mo_2.coldim());
819 RefSCVector vals22 = kit->vector(ao2mo_2.coldim());
820 vals11.assign(values11);
821 vals12.assign(values12);
822 vals21.assign(values21);
823 vals22.assign(values22);
824 delete[] values11;
825 delete[] values12;
826 delete[] values21;
827 delete[] values22;
828
829 RefSCVector movals11 = ao2mo_1 * vals11;
830 RefSCVector movals12 = ao2mo_1 * vals12;
831 RefSCVector movals21 = ao2mo_2 * vals21;
832 RefSCVector movals22 = ao2mo_2 * vals22;
833
834 kit = new LocalSCMatrixKit;
835 RefSCVector vals = kit->vector(pairdim);
836
837 MOPairIterFactory PIFactory;
838 Ref<SpatialMOPairIter> ij_iter = PIFactory.mopairiter(space1,space2);
839 for(ij_iter->start();int(*ij_iter.pointer());ij_iter->next()) {
840 const int i = ij_iter->i();
841 const int j = ij_iter->j();
842 const int ij_aa = ij_iter->ij_aa();
843 const int ij_ab = ij_iter->ij_ab();
844 const int ij_ba = ij_iter->ij_ba();
845
846 if (equiv) {
847 if (ij_aa != -1) {
848 const double value = movals11.get_element(i) * movals22.get_element(j) -
849 movals12.get_element(i) * movals21.get_element(j);
850 vals.set_element(ij_aa,value);
851 }
852 }
853 else {
854 const double value = movals11.get_element(i) * movals22.get_element(j);
855 vals.set_element(ij_ab,value);
856 if (space1_eq_space2 && ij_ab != ij_ba) {
857 const double value = movals11.get_element(j) * movals22.get_element(i);
858 vals.set_element(ij_ba,value);
859 }
860 }
861
862 }
863
864 vdata1->~ValueData();
865 vdata2->~ValueData();
866
867 return vals;
868 }
869
print(std::ostream & so) const870 void MP2R12Energy::print(std::ostream& so) const
871 {
872 }
873
print_pair_energies(bool spinadapted,std::ostream & so)874 void MP2R12Energy::print_pair_energies(bool spinadapted, std::ostream& so)
875 {
876 compute();
877
878 char* SA_str;
879 switch (stdapprox_) {
880 case LinearR12::StdApprox_A:
881 SA_str = strdup("A");
882 break;
883
884 case LinearR12::StdApprox_Ap:
885 SA_str = strdup("A'");
886 break;
887
888 case LinearR12::StdApprox_B:
889 SA_str = strdup("B");
890 break;
891
892 default:
893 throw std::runtime_error("MP2R12Energy::print_pair_energies -- stdapprox_ is not valid");
894 }
895
896 Ref<R12IntEvalInfo> r12info = r12eval_->r12info();
897 int nocc_act = r12info->nocc_act();
898 double escf = r12info->ref()->energy();
899
900 double emp2tot_aa = 0.0;
901 double emp2tot_ab = 0.0;
902 double er12tot_aa = 0.0;
903 double er12tot_ab = 0.0;
904 double emp2tot_0 = 0.0;
905 double emp2tot_1 = 0.0;
906 double er12tot_0 = 0.0;
907 double er12tot_1 = 0.0;
908
909 RefSCVector emp2_aa = r12eval_->emp2_aa();
910 RefSCVector emp2_ab = r12eval_->emp2_ab();
911
912 /*---------------------------------------
913 Spin-adapt pair energies, if necessary
914 ---------------------------------------*/
915 if (!spinadapted) {
916
917 so << endl << indent << "Alpha-alpha MBPT2-R12/" << SA_str << " pair energies:" << endl;
918 so << indent << scprintf(" i j mp2(ij) r12(ij) mp2-r12(ij)") << endl;
919 so << indent << scprintf(" ----- ----- ------------ ------------ ------------") << endl;
920 for(int i=0,ij=0;i<nocc_act;i++)
921 for(int j=0;j<i;j++,ij++) {
922 so << indent << scprintf(" %3d %3d %12.9lf %12.9lf %12.9lf",i+1,j+1,
923 emp2_aa->get_element(ij),
924 er12_aa_->get_element(ij),
925 emp2r12_aa_->get_element(ij)) << endl;
926 }
927
928 so << endl << indent << "Alpha-beta MBPT2-R12/" << SA_str << " pair energies:" << endl;
929 so << indent << scprintf(" i j mp2(ij) r12(ij) mp2-r12(ij)") << endl;
930 so << indent << scprintf(" ----- ----- ------------ ------------ ------------") << endl;
931 for(int i=0,ij=0;i<nocc_act;i++)
932 for(int j=0;j<nocc_act;j++,ij++) {
933 so << indent << scprintf(" %3d %3d %12.9lf %12.9lf %12.9lf",i+1,j+1,
934 emp2_ab->get_element(ij),
935 er12_ab_->get_element(ij),
936 emp2r12_ab_->get_element(ij)) << endl;
937 }
938
939 }
940 else {
941
942 Ref<SCMatrixKit> localkit = er12_aa_.kit();
943 RefSCVector emp2r12_0 = localkit->vector(r12eval_->dim_oo_s());
944 RefSCVector emp2r12_1 = localkit->vector(r12eval_->dim_oo_t());
945 RefSCVector emp2_0 = localkit->vector(r12eval_->dim_oo_s());
946 RefSCVector emp2_1 = localkit->vector(r12eval_->dim_oo_t());
947 RefSCVector er12_0 = localkit->vector(r12eval_->dim_oo_s());
948 RefSCVector er12_1 = localkit->vector(r12eval_->dim_oo_t());
949
950 // Triplet pairs are easy
951 emp2r12_1->assign(emp2r12_aa_);
952 emp2r12_1->scale(1.5);
953 emp2_1->assign(emp2_aa);
954 emp2_1->scale(1.5);
955 er12_1->assign(er12_aa_);
956 er12_1->scale(1.5);
957
958 // Singlet pairs are a bit trickier
959 int ij_s = 0;
960 for(int i=0; i<nocc_act; i++)
961 for(int j=0; j<=i; j++, ij_s++) {
962 int ij_ab = i*nocc_act + j;
963 int ij_aa = i*(i-1)/2 + j;
964 double eab, eaa, e_s;
965
966 eab = emp2r12_ab_->get_element(ij_ab);
967 if (i != j)
968 eaa = emp2r12_aa_->get_element(ij_aa);
969 else
970 eaa = 0.0;
971 e_s = (i != j ? 2.0 : 1.0) * eab - 0.5 * eaa;
972 emp2r12_0->set_element(ij_s,e_s);
973
974 eab = emp2_ab->get_element(ij_ab);
975 if (i != j)
976 eaa = emp2_aa->get_element(ij_aa);
977 else
978 eaa = 0.0;
979 e_s = (i != j ? 2.0 : 1.0) * eab - 0.5 * eaa;
980 emp2_0->set_element(ij_s,e_s);
981
982 eab = er12_ab_->get_element(ij_ab);
983 if (i != j)
984 eaa = er12_aa_->get_element(ij_aa);
985 else
986 eaa = 0.0;
987 e_s = (i != j ? 2.0 : 1.0) * eab - 0.5 * eaa;
988 er12_0->set_element(ij_s,e_s);
989 }
990
991 // compute total singlet and triplet energies
992 RefSCVector unit_0 = localkit->vector(r12eval_->dim_oo_s());
993 RefSCVector unit_1 = localkit->vector(r12eval_->dim_oo_t());
994 unit_0->assign(1.0);
995 unit_1->assign(1.0);
996 emp2tot_0 = emp2_0.dot(unit_0);
997 emp2tot_1 = emp2_1.dot(unit_1);
998 er12tot_0 = er12_0.dot(unit_0);
999 er12tot_1 = er12_1.dot(unit_1);
1000
1001 so << endl << indent << "Singlet MBPT2-R12/" << SA_str << " pair energies:" << endl;
1002 so << indent << scprintf(" i j mp2(ij) r12(ij) mp2-r12(ij)") << endl;
1003 so << indent << scprintf(" ----- ----- ------------ ------------ ------------") << endl;
1004 for(int i=0,ij=0;i<nocc_act;i++)
1005 for(int j=0;j<=i;j++,ij++) {
1006 so << indent << scprintf(" %3d %3d %12.9lf %12.9lf %12.9lf",i+1,j+1,
1007 emp2_0->get_element(ij),
1008 er12_0->get_element(ij),
1009 emp2r12_0->get_element(ij)) << endl;
1010 }
1011
1012 so << endl << indent << "Triplet MBPT2-R12/" << SA_str << " pair energies:" << endl;
1013 so << indent << scprintf(" i j mp2(ij) r12(ij) mp2-r12(ij)") << endl;
1014 so << indent << scprintf(" ----- ----- ------------ ------------ ------------") << endl;
1015 for(int i=0,ij=0;i<nocc_act;i++)
1016 for(int j=0;j<i;j++,ij++) {
1017 so << indent << scprintf(" %3d %3d %12.9lf %12.9lf %12.9lf",i+1,j+1,
1018 emp2_1->get_element(ij),
1019 er12_1->get_element(ij),
1020 emp2r12_1->get_element(ij)) << endl;
1021 }
1022
1023 }
1024
1025
1026 double mp2_corr_energy_ = emp2tot_aa_() + emp2tot_ab_();
1027 double r12_corr_energy_ = er12tot_aa_() + er12tot_ab_();
1028
1029 ///////////////////////////////////////////////////////////////
1030 // The computation of the MP2 energy is now complete on each
1031 // node;
1032 ///////////////////////////////////////////////////////////////
1033
1034 if (spinadapted) {
1035 so <<endl<<indent
1036 <<scprintf("Singlet MP2 correlation energy [au]: %17.12lf\n", emp2tot_0);
1037 so <<indent
1038 <<scprintf("Triplet MP2 correlation energy [au]: %17.12lf\n", emp2tot_1);
1039 so <<indent
1040 <<scprintf("Singlet (MP2)-R12/%2s correlation energy [au]: %17.12lf\n", SA_str, er12tot_0);
1041 so <<indent
1042 <<scprintf("Triplet (MP2)-R12/%2s correlation energy [au]: %17.12lf\n", SA_str, er12tot_1);
1043 so <<indent
1044 <<scprintf("Singlet MP2-R12/%2s correlation energy [au]: %17.12lf\n", SA_str,
1045 emp2tot_0 + er12tot_0);
1046 so <<indent
1047 <<scprintf("Triplet MP2-R12/%2s correlation energy [au]: %17.12lf\n", SA_str,
1048 emp2tot_1 + er12tot_1);
1049 }
1050
1051 double etotal = escf + mp2_corr_energy_ + r12_corr_energy_;
1052 so <<endl<<indent
1053 <<scprintf("RHF energy [au]: %17.12lf\n", escf);
1054 so <<indent
1055 <<scprintf("MP2 correlation energy [au]: %17.12lf\n", mp2_corr_energy_);
1056 so <<indent
1057 <<scprintf("(MBPT2)-R12/%2s correlation energy [au]: %17.12lf\n", SA_str, r12_corr_energy_);
1058 so <<indent
1059 <<scprintf("MBPT2-R12/%2s correlation energy [au]: %17.12lf\n", SA_str,
1060 mp2_corr_energy_ + r12_corr_energy_);
1061 so <<indent
1062 <<scprintf("MBPT2-R12/%2s energy [au]: %17.12lf\n", SA_str, etotal) << endl;
1063
1064 so.flush();
1065
1066 free(SA_str);
1067
1068 return;
1069 }
1070
1071
1072 /////////////////////////////////////////////////////////////////////////////
1073
1074 // Local Variables:
1075 // mode: c++
1076 // c-file-style: "CLJ-CONDENSED"
1077 // End:
1078