21 #ifndef _libint2_src_lib_libint_boys_h_ 22 #define _libint2_src_lib_libint_boys_h_ 24 #if defined(__cplusplus) 30 #include <libint2/util/vector.h> 36 #include <type_traits> 40 #include <libint2/util/cxxstd.h> 41 #if LIBINT2_CPLUSPLUS_STD < 2011 42 # error "Libint2 C++ API requires C++11 support" 45 #include <libint2/boys_fwd.h> 46 #include <libint2/numeric.h> 47 #include <libint2/initialize.h> 49 #if HAVE_LAPACK // use F77-type interface for now, switch to LAPACKE later 50 extern "C" void dgesv_(
const int* n,
51 const int* nrhs,
double* A,
const int* lda,
52 int* ipiv,
double* b,
const int* ldb,
59 template<
typename Real>
66 for (
int i = 1; i <= ifac; i++) {
67 fac[i] = i * fac[i - 1];
79 for (
int i = 3; i <= idf; i++) {
80 df[i] = (i - 1) * df[i - 2];
85 bc_.resize((ibc+1)*(ibc+1));
86 std::fill(bc_.begin(), bc_.end(), Real(0));
89 for(
int i=1; i<=ibc; ++i)
90 bc[i] = bc[i-1] + (ibc+1);
92 for(
int i=0; i<=ibc; i++)
94 for(
int i=0; i<=ibc; i++)
95 for(
int j=1; j<=i; ++j)
96 bc[i][j] = bc[i][j-1] * Real(i-j+1) / Real(j);
99 for (
int i = 0; i < 128; i++) {
100 twoi1[i] = 1.0 / (Real(2.0) * i + Real(1.0));
101 ihalf[i] = Real(i) - Real(0.5);
109 std::vector<Real> fac;
110 std::vector<Real> df;
111 std::vector<Real*> bc;
120 std::vector<Real> bc_;
123 #define _local_min_macro(a,b) ((a) > (b) ? (a) : (b)) 143 template<
typename Real>
146 static std::shared_ptr<const FmEval_Reference> instance(
int , Real ) {
149 static auto instance_ = std::make_shared<const FmEval_Reference>();
155 static Real
eval(Real T,
size_t m) {
157 static const Real T_crit = std::numeric_limits<Real>::is_bounded ==
true ? -log( std::numeric_limits<Real>::min() * 100.5 / 2. ) : Real(0) ;
158 if (std::numeric_limits<Real>::is_bounded && T > T_crit)
159 throw std::overflow_error(
"FmEval_Reference<Real>::eval: Real lacks precision for the given value of argument T");
160 static const Real half = Real(1)/2;
161 Real denom = (m + half);
163 Real term = exp(-T) / (2 * denom);
166 const Real epsilon = get_epsilon(T);
167 const Real epsilon_divided_10 = epsilon / 10;
171 term = old_term * T / denom;
175 }
while (term > sum * epsilon_divided_10 || old_term < term);
184 static void eval(Real* Fm, Real T,
size_t mmax) {
189 for(
size_t m=0; m<=mmax; ++m)
206 template<
typename Real>
209 static std::shared_ptr<const FmEval_Reference2> instance(
int , Real ) {
212 static auto instance_ = std::make_shared<const FmEval_Reference2>();
221 static void eval(Real* Fm, Real t,
size_t mmax) {
227 const Real two_over_sqrt_pi{1.128379167095512573896158903121545171688101258657997713688171443421284936882986828973487320404214727};
228 const Real K = 1/two_over_sqrt_pi;
232 auto sqrt_t = sqrt(t);
233 Fm[0] = K*erf(sqrt_t)/sqrt_t;
235 for(
size_t m=0; m<=mmax-1; m++) {
236 Fm[m+1] = ((2*m + 1)*Fm[m] - et)/(t2);
246 template <
typename Real =
double>
249 #include <libint2/boys_cheb7.h> 251 static_assert(std::is_same<Real,double>::value,
"FmEval_Chebyshev7 only supports double as the real type");
253 static constexpr
const int ORDER = interpolation_order;
254 static constexpr
const int ORDERp1 = ORDER+1;
256 static constexpr
const Real T_crit = cheb_table_tmax;
257 static constexpr Real delta = cheb_table_delta;
258 static constexpr Real one_over_delta = 1/delta;
270 mmax(m_max), numbers_(14) {
271 #if defined(__x86_64__) || defined(_M_X64) || defined(__i386) || defined(_M_IX86) 272 # if !defined(__AVX__) && defined(NDEBUG) 274 static bool printed_performance_warning =
false;
275 if (!printed_performance_warning) {
277 <<
"libint2::FmEval_Chebyshev7 on x86(-64) platforms needs AVX support for best performance" 279 printed_performance_warning =
true;
284 if (precision < std::numeric_limits<double>::epsilon())
285 throw std::invalid_argument(std::string(
"FmEval_Chebyshev7 does not support precision smaller than ") + std::to_string(std::numeric_limits<double>::epsilon()));
286 if (mmax > cheb_table_mmax)
287 throw std::invalid_argument(
288 "FmEval_Chebyshev7::init() : requested mmax exceeds the " 300 static std::shared_ptr<const FmEval_Chebyshev7>
instance(
int m_max,
double = 0.0) {
304 static auto instance_ = std::make_shared<const FmEval_Chebyshev7>(m_max);
306 while (instance_->max_m() < m_max) {
307 static std::mutex mtx;
308 std::lock_guard<std::mutex> lck(mtx);
309 if (instance_->max_m() < m_max) {
310 auto new_instance = std::make_shared<const FmEval_Chebyshev7>(m_max);
311 instance_ = new_instance;
325 inline void eval(Real* Fm, Real x,
int m_max)
const {
330 const double one_over_x = 1/x;
331 Fm[0] = 0.88622692545275801365 * sqrt(one_over_x);
335 for (
int i = 1; i <= m_max; i++)
336 Fm[i] = Fm[i - 1] * numbers_.ihalf[i] * one_over_x;
344 const Real x_over_delta = x * one_over_delta;
345 const int iv = int(x_over_delta);
346 const Real xd = x_over_delta - (Real)iv - 0.5;
350 const auto x2 = xd*xd;
351 const auto x3 = x2*xd;
352 const auto x4 = x2*x2;
353 const auto x5 = x2*x3;
354 const auto x6 = x3*x3;
355 const auto x7 = x3*x4;
360 const Real *d = c + (ORDERp1) * (iv * (mmax+1) + m_min);
364 const int unroll_size = 4;
365 const int m_fence = (m_max + 2 - unroll_size);
366 for(; m<m_fence; m+=unroll_size, d+=ORDERp1*unroll_size) {
368 d20v, d21v, d30v, d31v;
382 const int unroll_size = 2;
383 const int m_fence = (m_max + 2 - unroll_size);
384 for(; m<m_fence; m+=unroll_size, d+=ORDERp1*unroll_size) {
397 for(; m<=m_max; ++m, d+=ORDERp1) {
401 Fm[m] = horizontal_add(d0v * x0vec + d1v * x1vec);
404 #else // AVX not available 405 for(m=m_min; m<=m_max; ++m, d+=ORDERp1) {
435 int status = posix_memalign(&result, ORDERp1*
sizeof(Real), (mmax + 1) * cheb_table_nintervals * ORDERp1 *
sizeof(Real));
437 if (status == EINVAL)
438 throw std::logic_error(
439 "FmEval_Chebyshev7::init() : posix_memalign failed, alignment must be a power of 2 at least as large as sizeof(void *)");
440 if (status == ENOMEM)
441 throw std::bad_alloc();
444 c = static_cast<Real*>(result);
448 for (std::size_t iv = 0; iv < cheb_table_nintervals; ++iv) {
450 std::copy(cheb_table[iv], cheb_table[iv]+(mmax+1)*ORDERp1, c+(iv*(mmax+1))*ORDERp1);
456 #if LIBINT2_CONSTEXPR_STATICS 457 template <
typename Real>
459 double FmEval_Chebyshev7<Real>::cheb_table[FmEval_Chebyshev7<Real>::cheb_table_nintervals][(FmEval_Chebyshev7<Real>::cheb_table_mmax+1)*(FmEval_Chebyshev7<Real>::interpolation_order+1)];
463 template <
typename Real>
464 double FmEval_Chebyshev7<Real>::cheb_table[FmEval_Chebyshev7<Real>::cheb_table_nintervals][(FmEval_Chebyshev7<Real>::cheb_table_mmax+1)*(FmEval_Chebyshev7<Real>::interpolation_order+1)];
471 const double oon[] = {0.0, 1.0, 1.0/2.0, 1.0/3.0, 1.0/4.0, 1.0/5.0, 1.0/6.0, 1.0/7.0, 1.0/8.0, 1.0/9.0, 1.0/10.0, 1.0/11.0};
480 template<
typename Real =
double,
int INTERPOLATION_ORDER = 7>
483 static_assert(std::is_same<Real,double>::value,
"FmEval_Taylor only supports double as the real type");
485 static const int max_interp_order = 8;
486 static const bool INTERPOLATION_AND_RECURSION =
false;
487 const double soft_zero_;
491 soft_zero_(1e-6), cutoff_(
precision), numbers_(
492 INTERPOLATION_ORDER + 1, 2 * (mmax + INTERPOLATION_ORDER - 1)) {
494 #if defined(__x86_64__) || defined(_M_X64) || defined(__i386) || defined(_M_IX86) 495 # if !defined(__AVX__) && defined(NDEBUG) 497 static bool printed_performance_warning =
false;
498 if (!printed_performance_warning) {
500 <<
"libint2::FmEval_Taylor on x86(-64) platforms needs AVX support for best performance" 502 printed_performance_warning =
true;
510 const double sqrt_pi = std::sqrt(M_PI);
518 * std::pow(cutoff_ * numbers_.fac[INTERPOLATION_ORDER + 1],
519 1.0 / INTERPOLATION_ORDER);
520 oodelT_ = 1.0 / delT_;
521 max_m_ = mmax + INTERPOLATION_ORDER - 1;
523 T_crit_ =
new Real[max_m_ + 1];
526 for (
int m = max_m_; m >= 0; --m) {
533 double T = -log(cutoff_);
534 const double egamma = cutoff_ * sqrt_pi * numbers_.df[2 * m]
539 const double damping_factor = 0.2;
542 func = std::pow(T, m - 0.5) * std::exp(-T) - egamma;
543 const double dfuncdT = ((m - 0.5) * std::pow(T, m - 1.5)
544 - std::pow(T, m - 0.5)) * std::exp(-T);
550 double deltaT = -func / dfuncdT;
551 const double sign_deltaT = (deltaT > 0.0) ? 1.0 : -1.0;
552 const double max_deltaT = damping_factor * T;
553 if (std::fabs(deltaT) > max_deltaT)
554 deltaT = sign_deltaT * max_deltaT;
560 }
while (std::fabs(func / egamma) >= soft_zero_);
562 const int T_idx = (int) std::floor(T_new / delT_);
563 max_T_ = std::max(max_T_, T_idx);
568 const int nrow = max_T_ + 1;
569 const int ncol = max_m_ + 1;
570 grid_ =
new Real*[nrow];
571 grid_[0] =
new Real[nrow * ncol];
573 for (
int r = 1; r < nrow; ++r)
574 grid_[r] = grid_[r - 1] + ncol;
583 for (
int T_idx = max_T_; T_idx >= 0; --T_idx) {
584 const double T = T_idx * delT_;
597 static std::shared_ptr<const FmEval_Taylor>
instance(
unsigned int mmax, Real
precision = std::numeric_limits<Real>::epsilon()) {
601 static auto instance_ = std::make_shared<const FmEval_Taylor>(mmax,
precision);
603 while (instance_->max_m() < mmax || instance_->precision() >
precision) {
604 static std::mutex mtx;
605 std::lock_guard<std::mutex> lck(mtx);
606 if (instance_->max_m() < mmax || instance_->precision() >
precision) {
607 auto new_instance = std::make_shared<const FmEval_Taylor>(mmax,
precision);
608 instance_ = new_instance;
616 int max_m()
const {
return max_m_ - INTERPOLATION_ORDER + 1; }
625 void eval(Real* Fm, Real T,
int mmax)
const {
626 const double sqrt_pio2 = 1.2533141373155002512;
627 const double two_T = 2.0 * T;
630 const int mmin = INTERPOLATION_AND_RECURSION ? mmax : 0;
634 const bool use_upward_recursion =
true;
635 if (use_upward_recursion) {
637 if (T > T_crit_[0]) {
638 const double one_over_x = 1.0/T;
639 Fm[0] = 0.88622692545275801365 * sqrt(one_over_x);
643 for (
int i = 1; i <= mmax; i++)
644 Fm[i] = Fm[i - 1] * numbers_.ihalf[i] * one_over_x;
650 if (T > T_crit_[mmax]) {
651 double pow_two_T_to_minusjp05 = std::pow(two_T, -mmax - 0.5);
652 for (
int m = mmax; m >= mmin; --m) {
654 Fm[m] = numbers_.df[2 * m] * sqrt_pio2 * pow_two_T_to_minusjp05;
655 pow_two_T_to_minusjp05 *= two_T;
660 const int T_ind = (int) (0.5 + T * oodelT_);
661 const Real h = T_ind * delT_ - T;
662 const Real* F_row = grid_[T_ind] + mmin;
664 #if defined (__AVX__) 666 if (INTERPOLATION_ORDER == 3 || INTERPOLATION_ORDER == 7) {
667 const double h2 = h*h*oon[2];
668 const double h3 = h2*h*oon[3];
670 if (INTERPOLATION_ORDER == 7) {
671 const double h4 = h3*h*oon[4];
672 const double h5 = h4*h*oon[5];
673 const double h6 = h5*h*oon[6];
674 const double h7 = h6*h*oon[7];
684 const int unroll_size = 2;
685 const int m_fence = (mmax + 2 - unroll_size);
686 for(; m<m_fence; m+=unroll_size, F_row+=unroll_size) {
689 if (INTERPOLATION_ORDER == 3 || INTERPOLATION_ORDER == 7) {
693 if (INTERPOLATION_ORDER == 7) {
696 fm01 += horizontal_add(fr0_4567*h4567, fr1_4567*h4567);
702 Real total0 = 0.0, total1 = 0.0;
703 for(
int i=INTERPOLATION_ORDER; i>=1; --i) {
704 total0 = oon[i]*h*(F_row[i] + total0);
705 total1 = oon[i]*h*(F_row[i+1] + total1);
707 Fm[m] = F_row[0] + total0;
708 Fm[m+1] = F_row[1] + total1;
716 if (INTERPOLATION_ORDER == 3 || INTERPOLATION_ORDER == 7) {
718 if (INTERPOLATION_ORDER == 7) {
722 Fm[m] = horizontal_add(fr0123*h0123 + fr4567*h4567);
725 Fm[m] = horizontal_add(fr0123*h0123);
731 for(
int i=INTERPOLATION_ORDER; i>=1; --i) {
732 total = oon[i]*h*(F_row[i] + total);
734 Fm[m] = F_row[0] + total;
754 if (INTERPOLATION_AND_RECURSION && mmin > 0) {
755 const Real exp_mT = std::exp(-T);
756 for (
int m = mmin - 1; m >= 0; --m) {
757 Fm[m] = (exp_mT + two_T * Fm[m+1]) * numbers_.twoi1[m];
787 static double truncation_error(
unsigned int m,
double T) {
788 const double m2= m * m;
789 const double m3= m2 * m;
790 const double m4= m2 * m2;
791 const double T2= T * T;
792 const double T3= T2 * T;
793 const double T4= T2 * T2;
794 const double T5= T2 * T3;
796 const double result = exp(-T) * (105 + 16*m4 + 16*m3*(T - 8) - 30*T + 12*T2
797 - 8*T3 + 16*T4 + 8*m2*(43 - 9*T + 2*T2) +
798 4*m*(-88 + 23*T - 8*T2 + 4*T3))/(32*T5);
809 static double truncation_error(
double T) {
810 const double result = exp(-T) /(2*T);
817 constexpr
double pow10(
long exp) {
818 return (exp == 0) ? 1. : ((exp > 0) ? 10. * pow10(exp-1) : 0.1 * pow10(exp+1));
833 template<
typename Real =
double>
838 #include <libint2/tenno_cheb.h> 840 static_assert(std::is_same<Real,double>::value,
"TennoGmEval only supports double as the real type");
842 static const int mmin_ = -1;
843 static constexpr
const Real Tmax = (1 << cheb_table_tmaxlog2);
844 static constexpr
const Real Umax = detail::pow10(cheb_table_umaxlog10);
845 static constexpr
const Real Umin = detail::pow10(cheb_table_uminlog10);
846 static constexpr
const std::size_t ORDERp1 = interpolation_order + 1;
847 static constexpr
const Real maxabsprecision = 1.4e-14;
857 #if defined(__x86_64__) || defined(_M_X64) || defined(__i386) || defined(_M_IX86) 858 # if !defined(__AVX__) && defined(NDEBUG) 860 static bool printed_performance_warning =
false;
861 if (!printed_performance_warning) {
863 <<
"libint2::TennoGmEval on x86(-64) platforms needs AVX support for best performance" 865 printed_performance_warning =
true;
877 if (mmax > cheb_table_mmax)
878 throw std::invalid_argument(
879 "TennoGmEval::init() : requested mmax exceeds the " 890 static std::shared_ptr<const TennoGmEval>
instance(
int m_max,
double = 0) {
894 static auto instance_ = std::make_shared<const TennoGmEval>(m_max);
896 while (instance_->max_m() < m_max) {
897 static std::mutex mtx;
898 std::lock_guard<std::mutex> lck(mtx);
899 if (instance_->max_m() < m_max) {
900 auto new_instance = std::make_shared<const TennoGmEval>(m_max);
901 instance_ = new_instance;
908 unsigned int max_m()
const {
return mmax_; }
922 void eval_yukawa(Real* Gm, Real one_over_rho, Real T,
size_t mmax, Real zeta)
const {
923 assert(mmax <= mmax_);
925 const auto U = 0.25 * zeta * zeta * one_over_rho;
927 if (T > Tmax || U < Umin) {
928 eval_Gm_urr(Gm, T, U, mmax, 0);
930 interpolate_Gm<false>(Gm, T, U, 0, mmax);
944 void eval_slater(Real* Gm, Real one_over_rho, Real T,
size_t mmax, Real zeta)
const {
945 assert(mmax <= mmax_);
947 const auto U = 0.25 * zeta * zeta * one_over_rho;
949 const auto zeta_over_two_rho = 0.5 * zeta * one_over_rho;
951 eval_Gm_urr(Gm, T, U, mmax, -1);
953 interpolate_Gm<true>(Gm, T, U, zeta_over_two_rho, mmax);
964 template <
bool compute_Gm1>
965 static inline std::tuple<Real,Real> eval_G0_and_maybe_Gm1(Real T, Real U) {
968 const Real sqrtPi(1.7724538509055160272981674833411451827975494561224);
973 assert(compute_Gm1 ==
false);
974 if (T < std::numeric_limits<Real>::epsilon()) {
978 const Real sqrt_T = sqrt(T);
979 const Real sqrtPi_over_2 = sqrtPi * 0.5;
980 G_0 = sqrtPi_over_2 * erf(sqrt_T) / sqrt_T;
983 else if (T > std::numeric_limits<Real>::epsilon()) {
984 const Real sqrt_U = sqrt(U);
985 const Real sqrt_T = sqrt(T);
986 const Real oo_sqrt_T = 1 / sqrt_T;
987 const Real oo_sqrt_U = compute_Gm1 ? (1 / sqrt_U) : 0;
988 const Real kappa = sqrt_U - sqrt_T;
989 const Real lambda = sqrt_U + sqrt_T;
990 const Real sqrtPi_over_4 = sqrtPi * 0.25;
991 const Real pfac = sqrtPi_over_4;
992 const Real erfc_k = exp(kappa * kappa - T) * erfc(kappa);
993 const Real erfc_l = exp(lambda * lambda - T) * erfc(lambda);
995 G_m1 = compute_Gm1 ? pfac * (erfc_k + erfc_l) * oo_sqrt_U : 0.0;
996 G_0 = pfac * (erfc_k - erfc_l) * oo_sqrt_T;
999 const Real exp_U = exp(U);
1000 const Real sqrt_U = sqrt(U);
1002 const Real sqrtPi_over_2 = sqrtPi * 0.5;
1003 const Real oo_sqrt_U = compute_Gm1 ? (1 / sqrt_U) : 0;
1005 G_m1 = exp_U * sqrtPi_over_2 * erfc(sqrt_U) * oo_sqrt_U;
1007 G_0 = 1 - exp_U * sqrtPi * erfc(sqrt_U) * sqrt_U;
1010 return std::make_tuple(G_0, G_m1);
1020 static void eval_Gm_urr(Real* Gm, Real T, Real U,
size_t mmax,
long mmin) {
1021 assert(mmin == 0 || mmin == -1);
1025 const Real sqrt_U = sqrt(U);
1026 const Real sqrt_T = sqrt(T);
1027 const Real oo_sqrt_T = 1 / sqrt_T;
1028 const Real oo_sqrt_U = 1 / sqrt_U;
1029 const Real kappa = sqrt_U - sqrt_T;
1030 const Real lambda = sqrt_U + sqrt_T;
1031 const Real sqrtPi_over_4(0.44311346272637900682454187083528629569938736403060);
1032 const Real pfac = sqrtPi_over_4;
1033 const Real erfc_k = exp(kappa*kappa - T) * erfc(kappa);
1034 const Real erfc_l = exp(lambda*lambda - T) * erfc(lambda);
1037 Real& G_m1 = (mmin == -1) ? Gm[0] : G_m1_value;
1038 G_m1 = pfac * (erfc_k + erfc_l) * oo_sqrt_U;
1039 Real& G_0 = (mmin == -1) ? Gm[1] : Gm[0];
1040 G_0 = pfac * (erfc_k - erfc_l) * oo_sqrt_T;
1045 const Real oo_two_T = 0.5 / T;
1046 const Real two_U = 2.0 * U;
1047 const Real exp_mT = exp(-T);
1051 Real* Gmp1 = Gm0 + 1;
1052 for(
unsigned int m=1, two_m_minus_1=1; m<=mmax; ++m, two_m_minus_1+=2, ++Gmp1) {
1053 *Gmp1 = oo_two_T * ( two_m_minus_1 * *Gm0 + two_U * *Gmm1 - exp_mT);
1070 void interpolate_Gm(Real* Gm_vec, Real T, Real U, Real zeta_over_two_rho,
long mmax)
const {
1072 assert(U >= Umin && U <= Umax);
1075 auto linear_map_02 = [](Real x) {
1076 assert(x >= 0 && x <= 2);
1077 return (x - 1) * 0.5;
1080 auto log2_map = [](Real x, Real one_over_a) {
1081 return std::log2(x * one_over_a) - 0.5;
1084 auto log10_map = [](Real x, Real one_over_a) {
1085 return std::log10(x * one_over_a) - 0.5;
1089 const int Tint = T < 2 ? 0 : int(std::floor(std::log2(T))); assert(Tint >= 0 && Tint < 10);
1091 const constexpr Real one_over_2K[] = {1, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.0078125, 0.00390625, .001953125};
1093 const int Uint = int(std::floor(std::log10(U / Umin))); assert(Uint >= 0 && Uint < 10);
1095 const constexpr Real one_over_10K[] = {1. / detail::pow10(cheb_table_uminlog10),
1096 1. / detail::pow10(cheb_table_uminlog10+1),
1097 1. / detail::pow10(cheb_table_uminlog10+2),
1098 1. / detail::pow10(cheb_table_uminlog10+3),
1099 1. / detail::pow10(cheb_table_uminlog10+4),
1100 1. / detail::pow10(cheb_table_uminlog10+5),
1101 1. / detail::pow10(cheb_table_uminlog10+6),
1102 1. / detail::pow10(cheb_table_uminlog10+7),
1103 1. / detail::pow10(cheb_table_uminlog10+8),
1104 1. / detail::pow10(cheb_table_uminlog10+9)};
1105 const Real t = Tint == 0 ? linear_map_02(T) : log2_map(T, one_over_2K[Tint]);
1106 const Real u = log10_map(U, one_over_10K[Uint]);
1108 const int interval = Tint * 10 + Uint;
1110 #if defined(__AVX__) 1111 assert(ORDERp1 == 16);
1112 const auto t2 = t*t;
1113 const auto t3 = t2*t;
1114 const auto t4 = t2*t2;
1115 const auto t5 = t2*t3;
1116 const auto t6 = t3*t3;
1117 const auto t7 = t3*t4;
1118 const auto t8 = t4*t4;
1119 const auto t9 = t4*t5;
1120 const auto t10 = t5*t5;
1121 const auto t11 = t5*t6;
1122 const auto t12 = t6*t6;
1123 const auto t13 = t6*t7;
1124 const auto t14 = t7*t7;
1125 const auto t15 = t7*t8;
1126 const auto u2 = u*u;
1127 const auto u3 = u2*u;
1128 const auto u4 = u2*u2;
1129 const auto u5 = u2*u3;
1130 const auto u6 = u3*u3;
1131 const auto u7 = u3*u4;
1132 const auto u8 = u4*u4;
1133 const auto u9 = u4*u5;
1134 const auto u10 = u5*u5;
1135 const auto u11 = u5*u6;
1136 const auto u12 = u6*u6;
1137 const auto u13 = u6*u7;
1138 const auto u14 = u7*u7;
1139 const auto u15 = u7*u8;
1148 constexpr
const bool compute_Gmm10 =
true;
1150 if (compute_Gmm10) {
1153 std::tie(exp ? G0 : Gm_vec[0], Gmm1) = eval_G0_and_maybe_Gm1<exp>(T, U);
1155 Gm_vec[0] = (Gmm1 - G0) * zeta_over_two_rho;
1161 mmin_interp = (exp ==
true) ? -1 : 0;
1164 for(
long m=mmin_interp; m<=mmax; ++m){
1165 const Real *c_tuint = c_ + (ORDERp1) * (ORDERp1) * (interval * (mmax_ - mmin_ + 1) + (m - mmin_));
1166 #if defined(__AVX__) 1168 c20v, c21v, c22v, c23v, c30v, c31v, c32v, c33v,
1169 c40v, c41v, c42v, c43v, c50v, c51v, c52v, c53v,
1170 c60v, c61v, c62v, c63v, c70v, c71v, c72v, c73v;
1172 ca0v, ca1v, ca2v, ca3v, cb0v, cb1v, cb2v, cb3v,
1173 cc0v, cc1v, cc2v, cc3v, cd0v, cd1v, cd2v, cd3v,
1174 ce0v, ce1v, ce2v, ce3v, cf0v, cf1v, cf2v, cf3v;
1255 for(
int i=1; i!=16; ++i) {
1256 uvec[i] = uvec[i-1] * u;
1257 tvec[i] = tvec[i-1] * t;
1260 for(
int i=0, ij=0; i!=16; ++i) {
1261 for (
int j = 0; j != 16; ++j, ++ij) {
1262 Gm += c_tuint[ij] * tvec[i] * uvec[j];
1272 Gm_vec[m] = (Gmm1 - Gm) * zeta_over_two_rho;
1284 ExpensiveNumbers<Real> numbers_;
1291 int status = posix_memalign(&result, std::max(
sizeof(Real), 32ul), (mmax_ - mmin_ + 1) * cheb_table_nintervals * ORDERp1 * ORDERp1 *
sizeof(Real));
1293 if (status == EINVAL)
1294 throw std::logic_error(
1295 "TennoGmEval::init() : posix_memalign failed, alignment must be a power of 2 at least as large as sizeof(void *)");
1296 if (status == ENOMEM)
1297 throw std::bad_alloc();
1300 c_ = static_cast<Real*>(result);
1304 for (std::size_t iv = 0; iv < cheb_table_nintervals; ++iv) {
1306 std::copy(cheb_table[iv], cheb_table[iv]+((mmax_-mmin_)+1)*ORDERp1*ORDERp1, c_+(iv*((mmax_-mmin_)+1))*ORDERp1*ORDERp1);
1313 #if LIBINT2_CONSTEXPR_STATICS 1314 template <
typename Real>
1316 double TennoGmEval<Real>::cheb_table[TennoGmEval<Real>::cheb_table_nintervals][(TennoGmEval<Real>::cheb_table_mmax+2)*(TennoGmEval<Real>::interpolation_order+1)*(TennoGmEval<Real>::interpolation_order+1)];
1320 template <
typename Real>
1321 double TennoGmEval<Real>::cheb_table[TennoGmEval<Real>::cheb_table_nintervals][(TennoGmEval<Real>::cheb_table_mmax+2)*(TennoGmEval<Real>::interpolation_order+1)*(TennoGmEval<Real>::interpolation_order+1)];
1325 template<
typename Real,
int k>
1337 template <
typename Real>
1339 std::vector<Real> Fm_;
1340 std::vector<Real> g_i;
1341 std::vector<Real> r_i;
1342 std::vector<Real> oorhog_i;
1349 void init(
int mmax) {
1353 oorhog_i.resize(mmax+1);
1385 template<
typename Real,
int k>
1388 #ifndef LIBINT_USER_DEFINED_REAL 1398 detail::CoreEvalScratch<
GaussianGmEval<Real, k>>(mmax), mmax_(mmax),
1399 precision_(precision), fm_eval_(),
1400 numbers_(-1,-1,mmax) {
1401 assert(k == -1 || k == 0 || k == 2);
1404 fm_eval_ = FmEvalType::instance(mmax_, precision_);
1413 static std::shared_ptr<GaussianGmEval> instance(
unsigned int mmax, Real precision = std::numeric_limits<Real>::epsilon()) {
1415 assert(precision >= 0);
1417 static auto instance_ = std::make_shared<GaussianGmEval>(mmax, precision);
1419 while (instance_->max_m() < mmax || instance_->precision() > precision) {
1420 static std::mutex mtx;
1421 std::lock_guard<std::mutex> lck(mtx);
1422 if (instance_->max_m() < mmax || instance_->precision() > precision) {
1423 auto new_instance = std::make_shared<GaussianGmEval>(mmax, precision);
1424 instance_ = new_instance;
1432 int max_m()
const {
return mmax_; }
1434 Real precision()
const {
return precision_; }
1450 template <
typename AnyReal>
1451 void eval(Real* Gm, Real rho, Real T,
size_t mmax,
1452 const std::vector<std::pair<AnyReal, AnyReal> >& geminal,
1455 std::fill(Gm, Gm+mmax+1, Real(0));
1457 const auto sqrt_rho = sqrt(rho);
1458 const auto oo_sqrt_rho = 1/sqrt_rho;
1460 void* _scr = (scr == 0) ?
this : scr;
1461 auto& scratch = *(
reinterpret_cast<detail::CoreEvalScratch<GaussianGmEval<Real, -1
>>*>(_scr));
1462 for(
int i=1; i<=mmax; i++) {
1463 scratch.r_i[i] = scratch.r_i[i-1] * rho;
1467 typedef typename std::vector<std::pair<AnyReal, AnyReal> >::const_iterator citer;
1468 const citer gend = geminal.end();
1469 for(citer i=geminal.begin(); i!= gend; ++i) {
1471 const auto gamma = i->first;
1472 const auto gcoef = i->second;
1473 const auto rhog = rho + gamma;
1474 const auto oorhog = 1/rhog;
1476 const auto gorg = gamma * oorhog;
1477 const auto rorg = rho * oorhog;
1478 const auto sqrt_rho_org = sqrt_rho * oorhog;
1479 const auto sqrt_rhog = sqrt(rhog);
1480 const auto sqrt_rorg = sqrt_rho_org * sqrt_rhog;
1483 constexpr Real const_SQRTPI_2(0.88622692545275801364908374167057259139877472806119);
1484 const auto SS_K0G12_SS = gcoef * oo_sqrt_rho * const_SQRTPI_2 * rorg * sqrt_rorg * exp(-gorg*T);
1487 void* _scr = (scr == 0) ?
this : scr;
1488 auto& scratch = *(
reinterpret_cast<detail::CoreEvalScratch<GaussianGmEval<Real, -1
>>*>(_scr));
1490 const auto rorgT = rorg * T;
1491 fm_eval_->eval(&scratch.Fm_[0], rorgT, mmax);
1494 constexpr Real const_2_SQRTPI(1.12837916709551257389615890312154517);
1495 const auto pfac = const_2_SQRTPI * sqrt_rhog * SS_K0G12_SS;
1496 scratch.oorhog_i[0] = pfac;
1497 for(
int i=1; i<=mmax; i++) {
1498 scratch.g_i[i] = scratch.g_i[i-1] * gamma;
1499 scratch.oorhog_i[i] = scratch.oorhog_i[i-1] * oorhog;
1501 for(
int m=0; m<=mmax; m++) {
1503 Real* bcm = numbers_.bc[m];
1504 for(
int n=0; n<=m; n++) {
1505 ssss += bcm[n] * scratch.r_i[n] * scratch.g_i[m-n] * scratch.Fm_[n];
1507 Gm[m] += ssss * scratch.oorhog_i[m];
1513 auto ss_oper_ss_m = SS_K0G12_SS;
1514 Gm[0] += ss_oper_ss_m;
1515 for(
int m=1; m<=mmax; ++m) {
1516 ss_oper_ss_m *= gorg;
1517 Gm[m] += ss_oper_ss_m;
1524 const auto rorgT = rorg * T;
1525 const auto SS_K2G12_SS_0 = (1.5 + rorgT) * (SS_K0G12_SS * oorhog);
1526 const auto SS_K2G12_SS_m1 = rorg * (SS_K0G12_SS * oorhog);
1528 auto SS_K2G12_SS_gorg_m = SS_K2G12_SS_0 ;
1529 auto SS_K2G12_SS_gorg_m1 = SS_K2G12_SS_m1;
1530 Gm[0] += SS_K2G12_SS_gorg_m;
1531 for(
int m=1; m<=mmax; ++m) {
1532 SS_K2G12_SS_gorg_m *= gorg;
1533 Gm[m] += SS_K2G12_SS_gorg_m - m * SS_K2G12_SS_gorg_m1;
1534 SS_K2G12_SS_gorg_m1 *= gorg;
1545 std::shared_ptr<const FmEvalType> fm_eval_;
1547 ExpensiveNumbers<Real> numbers_;
1550 template <
typename GmEvalFunction>
1551 struct GenericGmEval :
private GmEvalFunction {
1553 typedef typename GmEvalFunction::value_type Real;
1555 GenericGmEval(
int mmax, Real precision) : GmEvalFunction(mmax, precision),
1556 mmax_(mmax), precision_(precision) {}
1558 static std::shared_ptr<const GenericGmEval> instance(
int mmax, Real precision = std::numeric_limits<Real>::epsilon()) {
1559 return std::make_shared<const GenericGmEval>(mmax, precision);
1562 template <
typename Real,
typename... ExtraArgs>
1563 void eval(Real* Gm, Real rho, Real T,
int mmax, ExtraArgs... args)
const {
1564 assert(mmax <= mmax_);
1565 (GmEvalFunction(*
this))(Gm, rho, T, mmax, std::forward<ExtraArgs>(args)...);
1569 int max_m()
const {
return mmax_; }
1571 Real precision()
const {
return precision_; }
1579 namespace os_core_ints {
1580 template <
typename Real,
int K>
struct r12_xx_K_gm_eval;
1581 template <
typename Real>
struct erfc_coulomb_gm_eval;
1586 template <
typename Real>
1587 struct CoreEvalScratch<os_core_ints::r12_xx_K_gm_eval<Real, 1>> {
1588 std::vector<Real> Fm_;
1589 CoreEvalScratch(
const CoreEvalScratch&) =
default;
1590 CoreEvalScratch(CoreEvalScratch&&) =
default;
1592 explicit CoreEvalScratch(
int mmax) { Fm_.resize(mmax+2); }
1595 template <
typename Real>
1596 struct CoreEvalScratch<os_core_ints::erfc_coulomb_gm_eval<Real>> {
1597 std::vector<Real> Fm_;
1598 CoreEvalScratch(
const CoreEvalScratch&) =
default;
1599 CoreEvalScratch(CoreEvalScratch&&) =
default;
1601 explicit CoreEvalScratch(
int mmax) { Fm_.resize(mmax+1); }
1606 namespace os_core_ints {
1609 template <
typename Real>
1610 struct delta_gm_eval {
1611 typedef Real value_type;
1613 delta_gm_eval(
unsigned int, Real) {}
1614 void operator()(Real* Gm, Real rho, Real T,
int mmax)
const {
1615 constexpr
static auto one_over_two_pi = 1.0 / (2.0 * M_PI);
1616 const auto G0 = exp(-T) * rho * one_over_two_pi;
1617 std::fill(Gm, Gm + mmax + 1, G0);
1626 template <
typename Real,
int K>
1627 struct r12_xx_K_gm_eval;
1629 template <
typename Real>
1630 struct r12_xx_K_gm_eval<Real, 1>
1631 :
private detail::CoreEvalScratch<r12_xx_K_gm_eval<Real, 1>> {
1632 typedef detail::CoreEvalScratch<r12_xx_K_gm_eval<Real, 1>> base_type;
1633 typedef Real value_type;
1635 #ifndef LIBINT_USER_DEFINED_REAL 1641 r12_xx_K_gm_eval(
unsigned int mmax, Real precision)
1643 fm_eval_ = FmEvalType::instance(mmax + 1, precision);
1645 void operator()(Real* Gm, Real rho, Real T,
int mmax) {
1646 fm_eval_->eval(&base_type::Fm_[0], T, mmax + 1);
1647 auto T_plus_m_plus_one = T + 1.0;
1648 Gm[0] = T_plus_m_plus_one * base_type::Fm_[0] - T * base_type::Fm_[1];
1649 auto minus_m = -1.0;
1650 T_plus_m_plus_one += 1.0;
1651 for (
auto m = 1; m <= mmax;
1652 ++m, minus_m -= 1.0, T_plus_m_plus_one += 1.0) {
1654 minus_m * base_type::Fm_[m - 1] + T_plus_m_plus_one * base_type::Fm_[m] - T * base_type::Fm_[m + 1];
1659 std::shared_ptr<const FmEvalType> fm_eval_;
1663 template <
typename Real>
1664 struct erf_coulomb_gm_eval {
1665 typedef Real value_type;
1667 #ifndef LIBINT_USER_DEFINED_REAL 1673 erf_coulomb_gm_eval(
unsigned int mmax, Real precision) {
1674 fm_eval_ = FmEvalType::instance(mmax, precision);
1676 void operator()(Real* Gm, Real rho, Real T,
int mmax, Real omega)
const {
1678 auto omega2 = omega * omega;
1679 auto omega2_over_omega2_plus_rho = omega2 / (omega2 + rho);
1680 fm_eval_->eval(Gm, T * omega2_over_omega2_plus_rho,
1684 auto ooversqrto2prho_exp_2mplus1 =
1685 sqrt(omega2_over_omega2_plus_rho);
1686 for (
auto m = 0; m <= mmax;
1687 ++m, ooversqrto2prho_exp_2mplus1 *= omega2_over_omega2_plus_rho) {
1688 Gm[m] *= ooversqrto2prho_exp_2mplus1;
1692 std::fill(Gm, Gm+mmax+1, Real{0});
1697 std::shared_ptr<const FmEvalType> fm_eval_;
1703 template <
typename Real>
1704 struct erfc_coulomb_gm_eval :
private 1705 detail::CoreEvalScratch<erfc_coulomb_gm_eval<Real>> {
1706 typedef detail::CoreEvalScratch<erfc_coulomb_gm_eval<Real>> base_type;
1707 typedef Real value_type;
1709 #ifndef LIBINT_USER_DEFINED_REAL 1715 erfc_coulomb_gm_eval(
unsigned int mmax, Real precision)
1717 fm_eval_ = FmEvalType::instance(mmax, precision);
1719 void operator()(Real* Gm, Real rho, Real T,
int mmax, Real omega) {
1720 fm_eval_->eval(&base_type::Fm_[0], T, mmax);
1721 std::copy(base_type::Fm_.cbegin(), base_type::Fm_.cbegin() + mmax + 1, Gm);
1723 auto omega2 = omega * omega;
1724 auto omega2_over_omega2_plus_rho = omega2 / (omega2 + rho);
1725 fm_eval_->eval(&base_type::Fm_[0], T * omega2_over_omega2_plus_rho,
1729 auto ooversqrto2prho_exp_2mplus1 =
1730 sqrt(omega2_over_omega2_plus_rho);
1731 for (
auto m = 0; m <= mmax;
1732 ++m, ooversqrto2prho_exp_2mplus1 *= omega2_over_omega2_plus_rho) {
1733 Gm[m] -= ooversqrto2prho_exp_2mplus1 * base_type::Fm_[m];
1739 std::shared_ptr<const FmEvalType> fm_eval_;
1753 template <
typename Real>
1757 return -std::exp(-zeta*x)/zeta;
1760 template <
typename Real>
1762 fngtg(
const std::vector<Real>& cc,
1763 const std::vector<Real>& aa,
1766 const Real x2 = x * x;
1767 const unsigned int n = cc.size();
1768 for(
unsigned int i=0; i<n; ++i)
1769 value += cc[i] * std::exp(- aa[i] * x2);
1776 template <
typename Real>
1778 wwtewklopper(Real x) {
1779 const Real x2 = x * x;
1780 return x2 * std::exp(-2 * x2);
1782 template <
typename Real>
1785 const Real x2 = x * x;
1786 const Real x6 = x2 * x2 * x2;
1787 return std::exp(-0.005 * x6);
1790 template <
typename Real>
1797 template <
typename Real>
1799 norm(
const std::vector<Real>& vec) {
1801 const unsigned int n = vec.size();
1802 for(
unsigned int i=0; i<n; ++i)
1803 value += vec[i] * vec[i];
1807 template <
typename Real>
1808 void LinearSolveDamped(
const std::vector<Real>& A,
1809 const std::vector<Real>& b,
1811 std::vector<Real>& x) {
1812 const size_t n = b.size();
1813 std::vector<Real> Acopy(A);
1814 for(
size_t m=0; m<n; ++m) Acopy[m*n + m] *= (1 + lambda);
1815 std::vector<Real> e(b);
1819 std::vector<int> ipiv(n);
1823 dgesv_(&n, &one, &Acopy[0], &n, &ipiv[0], &e[0], &n, &info);
1840 template <
typename Real>
1841 void stg_ng_fit(
unsigned int n,
1843 std::vector< std::pair<Real, Real> >& geminal,
1846 unsigned int npts = 1001) {
1849 std::vector<Real> cc(n, 1.0);
1850 std::vector<Real> aa(n);
1851 for(
unsigned int i=0; i<n; ++i)
1852 aa[i] = std::pow(3.0, (i + 2 - (n + 1)/2.0));
1855 Real ffnormfac = 0.0;
1857 for(
unsigned int i=0; i<n; ++i)
1858 for(
unsigned int j=0; j<n; ++j)
1859 ffnormfac += cc[i] * cc[j]/sqrt(aa[i] + aa[j]);
1860 const Real Nf = sqrt(2.0 * zeta) * zeta;
1861 const Real Nff = sqrt(2.0) / (sqrt(ffnormfac) *
1863 for(
unsigned int i=0; i<n; ++i) cc[i] *= -Nff/Nf;
1865 Real lambda0 = 1000;
1866 const Real nu = 3.0;
1867 const Real epsilon = 1e-15;
1868 const unsigned int maxniter = 200;
1871 std::vector<Real> xi(npts);
1872 for(
unsigned int i=0; i<npts; ++i) xi[i] = xmin + (xmax - xmin)*i/(npts - 1);
1874 std::vector<Real> err(npts);
1876 const size_t nparams = 2*n;
1877 std::vector<Real> J( npts * nparams );
1878 std::vector<Real> delta(nparams);
1885 Real errnormIm1 = 1e3;
1886 bool converged =
false;
1887 unsigned int iter = 0;
1888 while (!converged && iter < maxniter) {
1891 for(
unsigned int i=0; i<npts; ++i) {
1892 const Real x = xi[i];
1893 err[i] = (fstg(zeta, x) - fngtg(cc, aa, x)) * sqrt(ww(x));
1895 errnormI = norm(err)/sqrt((Real)npts);
1898 converged = std::abs((errnormI - errnormIm1)/errnormIm1) <= epsilon;
1899 if (converged)
break;
1900 errnormIm1 = errnormI;
1902 for(
unsigned int i=0; i<npts; ++i) {
1903 const Real x2 = xi[i] * xi[i];
1904 const Real sqrt_ww_x = sqrt(ww(xi[i]));
1905 const unsigned int ioffset = i * nparams;
1906 for(
unsigned int j=0; j<n; ++j)
1907 J[ioffset+j] = (std::exp(-aa[j] * x2)) * sqrt_ww_x;
1908 const unsigned int ioffsetn = ioffset+n;
1909 for(
unsigned int j=0; j<n; ++j)
1910 J[ioffsetn+j] = - sqrt_ww_x * x2 * cc[j] * std::exp(-aa[j] * x2);
1913 std::vector<Real> A( nparams * nparams);
1914 for(
size_t r=0, rc=0; r<nparams; ++r) {
1915 for(
size_t c=0; c<nparams; ++c, ++rc) {
1917 for(
size_t i=0, ir=r, ic=c; i<npts; ++i, ir+=nparams, ic+=nparams)
1918 Arc += J[ir] * J[ic];
1923 std::vector<Real> b( nparams );
1924 for(
size_t r=0; r<nparams; ++r) {
1926 for(
size_t i=0, ir=r; i<npts; ++i, ir+=nparams)
1927 br += J[ir] * err[i];
1934 for(
int l=-1; l<1000; ++l) {
1936 LinearSolveDamped(A, b, lambda0, delta );
1938 std::vector<double> cc_0(cc);
for(
unsigned int i=0; i<n; ++i) cc_0[i] += delta[i];
1939 std::vector<double> aa_0(aa);
for(
unsigned int i=0; i<n; ++i) aa_0[i] += delta[i+n];
1942 bool step_too_large =
false;
1943 for(
unsigned int i=0; i<n; ++i)
1944 if (aa_0[i] < 0.0) {
1945 step_too_large =
true;
1948 if (!step_too_large) {
1949 std::vector<double> err_0(npts);
1950 for(
unsigned int i=0; i<npts; ++i) {
1951 const double x = xi[i];
1952 err_0[i] = (fstg(zeta, x) - fngtg(cc_0, aa_0, x)) * sqrt(ww(x));
1954 const double errnorm_0 = norm(err_0)/sqrt((
double)npts);
1955 if (errnorm_0 < errnormI) {
1970 assert(not (iter == maxniter && errnormI > 1e-10));
1972 for(
unsigned int i=0; i<n; ++i)
1973 geminal[i] = std::make_pair(aa[i], cc[i]);
1980 #endif // header guard holds tables of expensive quantities
Definition: boys.h:60
double horizontal_add(VectorSSEDouble const &a)
Horizontal add.
Definition: vector_x86.h:232
static std::shared_ptr< const FmEval_Chebyshev7 > instance(int m_max, double=0.0)
Singleton interface allows to manage the lone instance; adjusts max m values as needed in thread-safe...
Definition: boys.h:300
int max_m() const
Definition: boys.h:616
SIMD vector of 4 double-precision floating-point real numbers, operations on which use AVX instructio...
Definition: vector_x86.h:639
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
std::ostream & verbose_stream()
Accessor for the disgnostics stream.
Definition: initialize.h:98
void convert(T *a) const
writes this to a
Definition: vector_x86.h:729
bool verbose()
Accessor for the verbose flag.
Definition: initialize.h:104
some evaluators need thread-local scratch, but most don't
Definition: boys.h:1331
Computes the Boys function, $ F_m (T) = \int_0^1 u^{2m} \exp(-T u^2) \, {\rm d}u $,...
Definition: boys.h:207
void eval_slater(Real *Gm, Real one_over_rho, Real T, size_t mmax, Real zeta) const
and
Definition: boys.h:944
void eval(Real *Fm, Real T, int mmax) const
computes Boys function values with m index in range [0,mmax]
Definition: boys.h:625
void load(T const *a)
loads a to this
Definition: vector_x86.h:720
static std::shared_ptr< const FmEval_Taylor > instance(unsigned int mmax, Real precision=std::numeric_limits< Real >::epsilon())
Singleton interface allows to manage the lone instance; adjusts max m and precision values as needed ...
Definition: boys.h:597
static std::shared_ptr< const TennoGmEval > instance(int m_max, double=0)
Singleton interface allows to manage the lone instance; adjusts max m values as needed in thread-safe...
Definition: boys.h:890
Computes the Boys function, $ F_m (T) = \int_0^1 u^{2m} \exp(-T u^2) \, {\rm d}u $,...
Definition: boys.h:247
static void eval(Real *Fm, Real t, size_t mmax)
fills up an array of Fm(T) for m in [0,mmax]
Definition: boys.h:221
void load_aligned(T const *a)
loads a to this
Definition: vector_x86.h:725
static void eval(Real *Fm, Real T, size_t mmax)
fills up an array of Fm(T) for m in [0,mmax]
Definition: boys.h:184
Real precision() const
Definition: boys.h:910
void eval(Real *Fm, Real x, int m_max) const
fills in Fm with computed Boys function values for m in [0,mmax]
Definition: boys.h:325
Computes the Boys function, , using single algorithm (asymptotic expansion).
Definition: boys.h:144
core integral for Yukawa and exponential interactions
Definition: boys.h:834
FmEval_Chebyshev7(int m_max, double precision=std::numeric_limits< double >::epsilon())
Definition: boys.h:269
TennoGmEval(unsigned int mmax, Real precision=-1)
Definition: boys.h:854
FmEval_Taylor(unsigned int mmax, Real precision=std::numeric_limits< Real >::epsilon())
Constructs the object to be able to compute Boys funcion for m in [0,mmax], with relative precision.
Definition: boys.h:490
void eval_yukawa(Real *Gm, Real one_over_rho, Real T, size_t mmax, Real zeta) const
Definition: boys.h:922
static Real eval(Real T, size_t m)
computes a single value of using MacLaurin series to full precision of Real
Definition: boys.h:155
Real precision() const
Definition: boys.h:618
Computes the Boys function, $ F_m (T) = \int_0^1 u^{2m} \exp(-T u^2) \, {\rm d}u $,...
Definition: boys.h:481
int max_m() const
Definition: boys.h:319