1 //
2 // Author: Hai-Anh Le <anh@u.northwestern.edu>
3 // Date: July 28, 2014
4 //
5
6 #include <iomanip>
7 #include "radialquad.h"
8 #include "so.h"
9
10 using namespace std;
11 using namespace bagel;
12 using namespace test;
13
14 namespace test {
15
16 class Gaussian_Int {
17 protected:
18 double exponent_;
19
20 public:
Gaussian_Int(const double exp)21 Gaussian_Int(const double exp) : exponent_(exp) {}
~Gaussian_Int()22 ~Gaussian_Int() {}
23
compute(const double r)24 double compute(const double r) { return exp(-exponent_ * r * r); }
exponent()25 double exponent() { return exponent_; }
26 };
27
28 }
29
30 using namespace test;
main()31 int main() {
32
33 /* ++++ TEST int r^n exp(-zeta r*r) <phi_A|lm><lm|l|lm'><lm'|phi_C> ++++ */
34 const array<double, 3> centreA = {{0.0000, 0.0000, 0.0000}};
35 const double alphaA = 1.0;
36 const int lA = 1;
37 const array<double, 3> centreC = {{0.0000, 0.0000, 0.0000}};
38 const double alphaC = 1.0;
39 const int lC = 1;
40
41 for (int lza = 0; lza <= lA; ++lza)
42 for (int lya = 0; lya <= lA-lza; ++lya)
43 for (int lzc = 0; lzc <= lC; ++lzc)
44 for (int lyc = 0; lyc <= lC-lzc; ++lyc) {
45 const int lxa = lA-lza-lya;
46 const int lxc = lC-lzc-lyc;
47 const array<int, 3> angA = {{lxa, lya, lza}};
48 auto cargaussA = make_shared<const CartesianGauss>(alphaA, angA, centreA);
49 // cargaussA->print();
50
51 const array<int, 3> angC = {{lxc , lyc, lzc}};
52 auto cargaussC = make_shared<const CartesianGauss>(alphaC, angC, centreC);
53 // cargaussC->print();
54
55 const array<double, 3> centreB = {{0.0000, 0.0000, 0.0000}};
56 const int l = 1;
57 const array<int, 2> lm = {{l, 0}}; // m can be any number st |m| <= l
58 auto rsh = make_shared<const SphHarmonics>(lm, centreB);
59 //rsh->print();
60 //SOIntegral soint(so);
61 //soint.compute(1.461297473945);
62
63 /* ECP Parameters */
64 const double ecp_coef = 1.0;
65 const double ecp_exp = 0.0;
66 const int ecp_r = 0;
67
68 const int max_iter = 100;
69 const double thresh_int = 10e-5;
70 /* Using normalized gaussian... */
71 cout << "lA = (" << lxa << ", " << lya << ", " << lza << ") l = " << l << " ";
72 cout << "lC = (" << lxc << ", " << lyc << ", " << lzc << ") ";
73 cout << "(iaa, rab, iab) = (";
74 for (int ic = 0; ic != 3; ++ic) {
75 tuple<shared_ptr<const CartesianGauss>, shared_ptr<const CartesianGauss>, shared_ptr<const SphHarmonics>, double, int, int>
76 so(cargaussA, cargaussC, rsh, ecp_exp, ecp_r, ic);
77 RadialInt<SOIntegral, tuple<shared_ptr<const CartesianGauss>, shared_ptr<const CartesianGauss>,
78 shared_ptr<const SphHarmonics>, double, int, int>> rad(so, false, max_iter, thresh_int);
79 const double out = ecp_coef*rad.integral();
80 if (ic != 2) {
81 cout << setprecision(12) << out << ", ";
82 } else {
83 cout << setprecision(12) << out << ")" << endl;
84 }
85 }
86 }
87
88 #if 0
89 /* CHECKED */
90 cout << " Test Radial Integration " << endl;
91 const double zeta = 1.0;
92 RadialInt<Gaussian_Int, const double> radial(zeta, true, max_iter, thresh_int);
93 cout << "Analytic = " << std::sqrt(pi__ / zeta) / 2.0 << endl;
94 #endif
95
96 return 0;
97
98 }
99
100