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