1 // 2 // Author: Hai-Anh Le <anh@u.northwestern.edu> 3 // Date: January 31, 2014 4 // 5 6 7 #ifndef __SRC_INTEGRAL_ECP_ANG_PROJ_RADIAL_H 8 #define __SRC_INTEGRAL_ECP_ANG_PROJ_RADIAL_H 9 10 #include <cmath> 11 #include <vector> 12 #include <iostream> 13 #include <iomanip> 14 #include <src/util/constants.h> 15 16 using namespace bagel; 17 18 namespace test { 19 20 template<typename T, typename... Value> 21 class RadialInt { 22 23 protected: 24 bool print_intermediate_; 25 int max_iter_; 26 double thresh_int_; 27 std::vector<double> x_, w_, r_; 28 double integral_; 29 30 public: 31 RadialInt(Value... tail, const bool print = false, const int max_iter = 100, const double thresh_int = 1e-10) : print_intermediate_(print)32 print_intermediate_(print), 33 max_iter_(max_iter), 34 thresh_int_(thresh_int) 35 { 36 T function(tail...); 37 integrate(function); 38 } 39 ~RadialInt()40 ~RadialInt() {} 41 integrate(T function)42 void integrate(T function) { 43 double previous = 0.0; 44 int ngrid = 31; 45 for (int iter = 0; iter != max_iter_; ++iter) { 46 transform_Becke(ngrid); 47 // transform_Log(ngrid, 3); //TODO: to be checked 48 // transform_Ahlrichs(ngrid); 49 double ans = 0.0; 50 int cnt = 0; 51 for (auto& it : r_) { 52 ans += function.compute(it) * w_[cnt++]; 53 } 54 const double error = ans - previous; 55 if (print_intermediate_) 56 std::cout << "Iter = " << std::setw(5) << iter << std::setw(10) << "npts = " << std::setw(10) << ngrid 57 << std::setw(10) << "ans = " << std::setw(20) << std::setprecision(10) << ans 58 << std::setw(10) << "err = " << std::setw(20) << std::setprecision(10) << error << std::endl; 59 if (fabs(error) < thresh_int_ && iter != 0) { 60 if (print_intermediate_) { 61 std::cout << "Integration converged..." << std::endl; 62 std::cout << "Radial integral = " << ans << std::endl; 63 } 64 integral_ = ans; 65 break; 66 } else if (iter == max_iter_-1) { 67 std::cout << "Max iteration exceeded..." << std::endl; 68 } 69 previous = ans; 70 x_.clear(); 71 w_.clear(); 72 r_.clear(); 73 ngrid = ngrid*2+1; 74 } 75 } 76 integral()77 double integral() { return integral_; } 78 79 void transform_Log(const int ngrid, const int m = 3) { // Mura and Knowles JCP, 104, 9848. 80 w_.resize(ngrid); 81 r_.resize(ngrid); 82 const double alpha = 5.0; 83 for (int i = 1; i <= ngrid; ++i) { 84 const double x = i / (ngrid + 1.0); 85 const double xm = 1.0 - std::pow(x, m); 86 r_[i-1] = - alpha * std::log(xm); 87 w_[i-1] = std::pow(r_[i-1], 2) * alpha * m * std::pow(x, m-1) / (xm * (ngrid + 1.0)); 88 } 89 } 90 transform_Ahlrichs(const int ngrid)91 void transform_Ahlrichs(const int ngrid) { // Treutler and Ahlrichs JCP, 102, 346. 92 GaussChebyshev2nd(ngrid); 93 r_.resize(ngrid); 94 const double alpha = 1.0; 95 for (int i = 0; i != ngrid; ++i) { 96 const double exp = 0.6; 97 const double prefactor = alpha / std::log(2.0); 98 r_[i] = prefactor * std::pow(1.0 + x_[i], exp) * std::log(2.0 / (1 - x_[i])); 99 w_[i] *= prefactor * (exp * std::pow(1.0 + x_[i], exp - 1.0) * std::log(2.0 / (1.0 - x_[i])) 100 + std::pow(1.0 + x_[i], exp) / (1.0 - x_[i])); 101 } 102 } 103 transform_Becke(const int ngrid)104 void transform_Becke(const int ngrid) { // Becke JCP, 88, 2547. 105 GaussChebyshev2nd(ngrid); 106 r_.resize(ngrid); 107 const double alpha = 1.0; 108 for (int i = 0; i != ngrid; ++i) { 109 r_[i] = alpha * (1.0 + x_[i]) / (1.0 - x_[i]); 110 w_[i] *= 2.0 * alpha / std::pow(1.0 - x_[i], 2); 111 } 112 } 113 GaussChebyshev2nd(const int ngrid)114 void GaussChebyshev2nd(const int ngrid) { 115 x_.resize(ngrid); 116 w_.resize(ngrid); 117 for (int i = 1; i <= ngrid; ++i) { 118 x_[i-1] = std::cos(i * pi__ / (ngrid + 1)); 119 w_[i-1] = pi__ * std::sin(i * pi__ / (ngrid + 1)) / (ngrid + 1); 120 } 121 } 122 123 }; 124 125 } 126 127 #endif 128