1 //
2 // compute_energy_a.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 #include <stdexcept>
29 #include <scconfig.h>
30 #include <util/misc/math.h>
31 #include <util/misc/formio.h>
32 #include <util/misc/timer.h>
33 #include <math/scmat/abstract.h>
34 #include <chemistry/qc/mbptr12/mbptr12.h>
35 #include <chemistry/qc/mbptr12/mp2r12_energy.h>
36
37 using namespace std;
38 using namespace sc;
39
40
41 void
compute_energy_a_()42 MBPT2_R12::compute_energy_a_()
43 {
44 tim_enter("mp2-r12/a energy");
45
46 if (r12eval_.null()) {
47 Ref<R12IntEvalInfo> r12info = new R12IntEvalInfo(this);
48 r12info->set_dynamic(dynamic_);
49 r12info->set_print_percent(print_percent_);
50 r12info->set_memory(mem_alloc);
51 r12eval_ = new R12IntEval(r12info,gbc_,ebc_,abs_method_,stdapprox_);
52 r12eval_->include_mp1(include_mp1_);
53 r12eval_->set_debug(debug_);
54 }
55 // This will actually compute the intermediates
56 r12eval_->compute();
57
58 double etotal = 0.0;
59
60 // Now we can compute and print pair energies
61 tim_enter("mp2-r12/a pair energies");
62 if (r12a_energy_.null())
63 r12a_energy_ = new MP2R12Energy(r12eval_,LinearR12::StdApprox_A,debug_);
64 r12a_energy_->print_pair_energies(spinadapted_);
65 etotal = r12a_energy_->energy();
66 tim_exit("mp2-r12/a pair energies");
67 if (stdapprox_ == LinearR12::StdApprox_Ap) {
68 tim_enter("mp2-r12/a' pair energies");
69 if (r12ap_energy_.null())
70 r12ap_energy_ = new MP2R12Energy(r12eval_,LinearR12::StdApprox_Ap,debug_);
71 r12ap_energy_->print_pair_energies(spinadapted_);
72
73 /*const double radius = 1.0;
74 SCVector3 r1(0.0,0.0,radius);
75 r12ap_energy_->compute_pair_function_ab(0,r1,r1);
76 ExEnv::out0() << endl<<endl;
77 const int nintervals = 100;
78 const double Phi_start = -(M_PI);
79 const double Phi_end = M_PI;
80 const double dPhi = (Phi_end - Phi_start) / nintervals;
81 const int npts = nintervals + 1;
82 for(int i=-nintervals/2; i<=nintervals/2; i++) {
83 const double Phi = i*dPhi;
84 const double z = radius * cos(Phi);
85 const double x = radius * sin(Phi);
86 SCVector3 r2(x,0.0,z);
87 ExEnv::out0() << indent << Phi;
88 r12ap_energy_->compute_pair_function_ab(0,r1,r2);
89 }*/
90 if (twopdm_grid_aa_.nonnull())
91 r12ap_energy_->compute_pair_function_aa(0,twopdm_grid_aa_);
92 if (twopdm_grid_ab_.nonnull())
93 r12ap_energy_->compute_pair_function_ab(0,twopdm_grid_ab_);
94
95 etotal = r12ap_energy_->energy();
96 tim_exit("mp2-r12/a' pair energies");
97 }
98
99 tim_exit("mp2-r12/a energy");
100
101 etotal += ref_energy();
102 set_energy(etotal);
103 set_actual_value_accuracy(reference_->actual_value_accuracy()
104 *ref_to_mp2_acc);
105
106 return;
107 }
108
109
110 ////////////////////////////////////////////////////////////////////////////
111
112 // Local Variables:
113 // mode: c++
114 // c-file-style: "CLJ-CONDENSED"
115 // End:
116