From 9256a78f406d6868f642d9cbc1b225bccff4eee0 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Tue, 26 Sep 2023 16:21:22 -0700 Subject: [PATCH 01/19] Add Lambert-W0 implementation for inverting exponential products --- include/integratorxx/util/lambert_w.hpp | 49 ++++++++++++++ test/1d_quadratures.cxx | 89 +++++++++++++++++++++++++ 2 files changed, 138 insertions(+) create mode 100644 include/integratorxx/util/lambert_w.hpp diff --git a/include/integratorxx/util/lambert_w.hpp b/include/integratorxx/util/lambert_w.hpp new file mode 100644 index 0000000..d2cc4aa --- /dev/null +++ b/include/integratorxx/util/lambert_w.hpp @@ -0,0 +1,49 @@ +#pragma once +#include +#include + +namespace IntegratorXX { + +// Newton-Raphson implementation of Lambert-W function for real arguments +template +RealType lambert_w_newton(RealType x, RealType w0) { + + const auto eps = std::numeric_limits::epsilon(); + const size_t max_iter = 100000; + RealType w; + size_t i; + for(i = 0; i < max_iter; ++i) { + const auto exp_term = std::exp(w0); + const auto w_exp_term = w0 * exp_term; + w = w0 - (w_exp_term - x) / (exp_term + w_exp_term); + if(std::abs(w - w0) < eps) break; + w0 = w; + } + + return w; +} + +template< typename RealType> +RealType lambert_w0(RealType x) { + const RealType one = 1; + const auto e = std::exp(one); + if(x < -one/e) return std::numeric_limits::quiet_NaN(); + if(x == RealType(0)) return 0.0; + if(x == -one/e) return -1.0; + + RealType w0; + if(x > e) { + w0 = std::log(x) - std::log(std::log(x)); + } else if(x > 0.0) { + w0 = x / e; + } else { + const auto xe = x * e; + const auto sqrt_term = std::sqrt(1.0 + xe); + w0 = xe * std::log(1.0 + sqrt_term) / (1.0 + xe + sqrt_term); + } + + return lambert_w_newton(x, w0); + +} + +} diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index aa30f90..75220a2 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -1,5 +1,6 @@ #include "catch2/catch_all.hpp" #include +#include #include #include #include @@ -253,3 +254,91 @@ TEST_CASE( "Delley", "[1d-quad]" ) { test_fn(974); } + +TEST_CASE("LMG", "[1d-quad]") { + using namespace IntegratorXX; + SECTION("Lambert W") { + size_t n = 50; + + SECTION("-1/e") { + REQUIRE(lambert_w0(-1.0 / std::exp(1.0)) == -1.0); + } + + double dx; + double x_min; + std::vector w_ref; + SECTION("(-1/e,0)") { + dx = (1.0/std::exp(1.0)) / (n-1); + x_min = -1.0/std::exp(1.0); + w_ref = { -1.0000000000000000e+00, -8.1043436396882251e-01, -7.3839124833180692e-01, + -6.8534285732547906e-01, -6.4201643964687316e-01, -6.0484100186256395e-01, + -5.7199479483565563e-01, -5.4240102237003140e-01, -5.1536144833012332e-01, + -4.9039331229166899e-01, -4.6714668763808043e-01, -4.4535855202428648e-01, + -4.2482542402809581e-01, -4.0538614978331228e-01, -3.8691058975718701e-01, + -3.6929191010441786e-01, -3.5244116871063264e-01, -3.3628341359547681e-01, + -3.2075480804994166e-01, -3.0580047091961676e-01, -2.9137282629002875e-01, + -2.7743032325972444e-01, -2.6393642934736894e-01, -2.5085882941845555e-01, + -2.3816878116899459e-01, -2.2584059140279072e-01, -2.1385118659837130e-01, + -2.0217975786260850e-01, -1.9080746514340260e-01, -1.7971718907514189e-01, + -1.6889332142980482e-01, -1.5832158709811098e-01, -1.4798889200596840e-01, + -1.3788319250617834e-01, -1.2799338266282390e-01, -1.1830919653009347e-01, + -1.0882112306525034e-01, -9.9520331741533477e-02, -9.0398607266594919e-02, + -8.1448292084908017e-02, -7.2662235562987412e-02, -6.4033748935368034e-02, + -5.5556565235647053e-02, -4.7224803557117209e-02, -3.9032937086775930e-02, + -3.0975764438825409e-02, -2.3048383882382367e-02, -1.5246170115525802e-02, + -7.5647532860519534e-03, 6.4184768611141850e-17 }; + } + SECTION("(0,e)") { + dx = std::exp(1.0) / (n-1); + x_min = 0.0; + w_ref = { 0.0000000000000000e+00, 5.2630934056544072e-02, 1.0035620951623317e-01, + 1.4409241373767120e-01, 1.8451244238478223e-01, 2.2212638256855327e-01, + 2.5733105242307630e-01, 2.9044182654506334e-01, 3.2171389544393414e-01, + 3.5135693697900772e-01, 3.7954552465616309e-01, 4.0642668865609727e-01, + 4.3212552273462468e-01, 4.5674941777919742e-01, 4.8039130984974648e-01, + 5.0313220781385692e-01, 5.2504318560584629e-01, 5.4618697067229260e-01, + 5.6661922372849582e-01, 5.8638957965403482e-01, 6.0554250149726696e-01, + 6.2411798675730124e-01, 6.4215215580659113e-01, 6.5967774546581859e-01, + 6.7672452563807517e-01, 6.9331965306433241e-01, 7.0948797333649716e-01, + 7.2525228005705800e-01, 7.4063353829291223e-01, 7.5565107811039622e-01, + 7.7032276290732105e-01, 7.8466513640828828e-01, 7.9869355151124433e-01, + 8.1242228362814539e-01, 8.2586463072191718e-01, 8.3903300188356378e-01, + 8.5193899600034140e-01, 8.6459347182518476e-01, 8.7700661055879869e-01, + 8.8918797189091481e-01, 9.0114654430979890e-01, 9.1289079037410314e-01, + 9.2442868754455021e-01, 9.3576776509144632e-01, 9.4691513752504819e-01, + 9.5787753493721217e-01, 9.6866133059280313e-01, 9.7927256606663837e-01, + 9.8971697418509541e-01, 9.9999999999999989e-01 }; + } + + SECTION("(e,inf)") { + dx = 1; + x_min = std::exp(1.0); + w_ref = { 1.0000000000000000e+00, 1.1626015113014869e+00, 1.2938344571456706e+00, + 1.4042003708950586e+00, 1.4996204192818812e+00, 1.5837783960587868e+00, + 1.6591292586935935e+00, 1.7273945838205316e+00, 1.7898301426340337e+00, + 1.8473811371144686e+00, 1.9007774734056948e+00, 1.9505949530993991e+00, + 1.9972960780278193e+00, 2.0412581329768691e+00, 2.0827930401857122e+00, + 2.1221617266541868e+00, 2.1595847338181962e+00, 2.1952501935697217e+00, + 2.2293199201838059e+00, 2.2619341295913213e+00, 2.2932151421500375e+00, + 2.3232703215087871e+00, 2.3521944316956782e+00, 2.3800715457358486e+00, + 2.4069766047089463e+00, 2.4329767015585579e+00, 2.4581321461297474e+00, + 2.4824973548124829e+00, 2.5061215984360530e+00, 2.5290496347488083e+00, + 2.5513222462703027e+00, 2.5729767000540669e+00, 2.5940471426162146e+00, + 2.6145649407274134e+00, 2.6345589767577975e+00, 2.6540559056765503e+00, + 2.6730803795436775e+00, 2.6916552443184649e+00, 2.7098017129924283e+00, + 2.7275395183923634e+00, 2.7448870484592360e+00, 2.7618614663662542e+00, + 2.7784788174751269e+00, 2.7947541248280854e+00, 2.8107014746227152e+00, + 2.8263340929076008e+00, 2.8416644145615426e+00, 2.8567041454717463e+00, + 2.8714643187019151e+00, 2.8859553453357312e+00 }; + } + + double x = x_min; + for(auto i = 0; i < w_ref.size(); ++i) { + auto w = lambert_w0(x); + REQUIRE_THAT(w, Catch::Matchers::WithinAbs(w_ref[i], 1e-15)); + //printf("%.16e %.16e %.4e\n", x, lambert_w0(x), + // lambert_w0(x) - boost::math::lambert_w0(x)); + x = x + dx; + } + } +} From 7fe816f389c6cc499ee41d23d8b38a8288025d33 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 27 Sep 2023 13:06:24 -0700 Subject: [PATCH 02/19] Added and tested Lambert Wm1 function --- include/integratorxx/util/lambert_w.hpp | 19 ++++++++ test/1d_quadratures.cxx | 60 ++++++++++++++++++++++--- 2 files changed, 73 insertions(+), 6 deletions(-) diff --git a/include/integratorxx/util/lambert_w.hpp b/include/integratorxx/util/lambert_w.hpp index d2cc4aa..a272514 100644 --- a/include/integratorxx/util/lambert_w.hpp +++ b/include/integratorxx/util/lambert_w.hpp @@ -46,4 +46,23 @@ RealType lambert_w0(RealType x) { } +template< typename RealType> +RealType lambert_wm1(RealType x) { + const RealType one = 1; + const auto e = std::exp(one); + if(x < -one/e) return std::numeric_limits::quiet_NaN(); + if(x >= RealType(0)) return std::numeric_limits::infinity(); + if(x == -one/e) return -1.0; + + RealType wm1; + if(x > -0.25) { + wm1 = std::log(-x) - std::log(-std::log(-x)); + } else { + wm1 = -1.0 - std::sqrt(2.0*(1.0 + x * e)); + } + + return lambert_w_newton(x, wm1); + +} + } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 75220a2..bae823a 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -1,6 +1,7 @@ #include "catch2/catch_all.hpp" #include #include +#include #include #include #include @@ -267,7 +268,8 @@ TEST_CASE("LMG", "[1d-quad]") { double dx; double x_min; std::vector w_ref; - SECTION("(-1/e,0)") { + int branch = 0; + SECTION("W0 (-1/e,0)") { dx = (1.0/std::exp(1.0)) / (n-1); x_min = -1.0/std::exp(1.0); w_ref = { -1.0000000000000000e+00, -8.1043436396882251e-01, -7.3839124833180692e-01, @@ -288,7 +290,7 @@ TEST_CASE("LMG", "[1d-quad]") { -3.0975764438825409e-02, -2.3048383882382367e-02, -1.5246170115525802e-02, -7.5647532860519534e-03, 6.4184768611141850e-17 }; } - SECTION("(0,e)") { + SECTION("W0 (0,e)") { dx = std::exp(1.0) / (n-1); x_min = 0.0; w_ref = { 0.0000000000000000e+00, 5.2630934056544072e-02, 1.0035620951623317e-01, @@ -310,7 +312,7 @@ TEST_CASE("LMG", "[1d-quad]") { 9.8971697418509541e-01, 9.9999999999999989e-01 }; } - SECTION("(e,inf)") { + SECTION("W0 (e,inf)") { dx = 1; x_min = std::exp(1.0); w_ref = { 1.0000000000000000e+00, 1.1626015113014869e+00, 1.2938344571456706e+00, @@ -331,14 +333,60 @@ TEST_CASE("LMG", "[1d-quad]") { 2.8263340929076008e+00, 2.8416644145615426e+00, 2.8567041454717463e+00, 2.8714643187019151e+00, 2.8859553453357312e+00 }; } + SECTION("W1 (-1/e,-0.25)") { + branch = -1; + dx = (1.0/std::exp(1.0) - 0.25) / (n-1); + x_min = -1.0/std::exp(1.0); + w_ref = { -1.0000000000000000e+00, -1.1189650446807282e+00, -1.1711582725750418e+00, + -1.2124855589488057e+00, -1.2482413372948526e+00, -1.2804726429337259e+00, + -1.3102276924123728e+00, -1.3381283510412803e+00, -1.3645793574591196e+00, + -1.3898616085396616e+00, -1.4141795078524806e+00, -1.4376873219364974e+00, + -1.4605049067902260e+00, -1.4827276158806526e+00, -1.5044328217548120e+00, + -1.5256843663192678e+00, -1.5465356909502610e+00, -1.5670320955171273e+00, + -1.5872124053596470e+00, -1.6071102254331719e+00, -1.6267549000718433e+00, + -1.6461722586522056e+00, -1.6653852027928755e+00, -1.6844141744163048e+00, + -1.7032775329676604e+00, -1.7219918624772559e+00, -1.7405722238112493e+00, + -1.7590323636442364e+00, -1.7773848889282362e+00, -1.7956414136079653e+00, + -1.8138126828282584e+00, -1.8319086787493766e+00, -1.8499387112277788e+00, + -1.8679114959619136e+00, -1.8858352221933619e+00, -1.9037176116563133e+00, + -1.9215659701557612e+00, -1.9393872329070991e+00, -1.9571880045721319e+00, + -1.9749745947677455e+00, -1.9927530496951602e+00, -2.0105291804333940e+00, + -2.0283085883553316e+00, -2.0460966880547553e+00, -2.0638987281149239e+00, + -2.0817198100013492e+00, -2.0995649053215297e+00, -2.1174388716610637e+00, + -2.1353464671775746e+00, -2.1532923641103534e+00 }; + } + SECTION("W1 (-0.25,0)") { + branch = -1; + dx = (0.25-1e-10) / (n-1); + x_min = -0.25; + w_ref = { -2.1532923641103499e+00, -2.1914997961914482e+00, -2.2299431488704058e+00, + -2.2686646700353417e+00, -2.3077062240124904e+00, -2.3471096678770036e+00, + -2.3869172040331206e+00, -2.4271717185362274e+00, -2.4679171134091620e+00, + -2.5091986404347137e+00, -2.5510632435016558e+00, -2.5935599164917043e+00, + -2.6367400838934105e+00, -2.6806580118111714e+00, -2.7253712578170455e+00, + -2.7709411692030885e+00, -2.8174334406865955e+00, -2.8649187445794793e+00, + -2.9134734489670580e+00, -2.9631804427033619e+00, -3.0141300902276456e+00, + -3.0664213446240929e+00, -3.1201630543731218e+00, -3.1754755084163060e+00, + -3.2324922762280690e+00, -3.2913624156132215e+00, -3.3522531424379425e+00, + -3.4153530856305858e+00, -3.4808762907481663e+00, -3.5490671909401166e+00, + -3.6202068424308829e+00, -3.6946208337432305e+00, -3.7726894411668441e+00, + -3.8548608453419511e+00, -3.9416685911800489e+00, -4.0337550432658453e+00, + -4.1319034964882038e+00, -4.2370830903073911e+00, -4.3505132002644684e+00, + -4.4737584327019935e+00, -4.6088735602501680e+00, -4.7586337059829287e+00, + -4.9269181859532187e+00, -5.1193905946137361e+00, -5.3448010688929983e+00, + -5.6177518509657292e+00, -5.9654942798986532e+00, -6.4488728297255040e+00, + -7.2605734538864919e+00, -2.6295240572605621e+01 }; + } double x = x_min; for(auto i = 0; i < w_ref.size(); ++i) { - auto w = lambert_w0(x); + auto w = branch == 0 ? lambert_w0(x) : lambert_wm1(x); REQUIRE_THAT(w, Catch::Matchers::WithinAbs(w_ref[i], 1e-15)); - //printf("%.16e %.16e %.4e\n", x, lambert_w0(x), - // lambert_w0(x) - boost::math::lambert_w0(x)); + //printf("%.16e %.16e %.4e\n", x, lambert_wm1(x), + // lambert_wm1(x) - boost::math::lambert_wm1(x)); x = x + dx; } + } + } From 7853fed721220e50092bc69d10be2cec3b599513 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 27 Sep 2023 16:25:31 -0700 Subject: [PATCH 03/19] Added LMG utility functions to compute upper/lower radial bounds --- include/integratorxx/quadratures/radial.hpp | 1 + .../integratorxx/quadratures/radial/lmg.hpp | 54 ++++++++++++++++ test/1d_quadratures.cxx | 62 +++++++++++++++++++ 3 files changed, 117 insertions(+) create mode 100644 include/integratorxx/quadratures/radial/lmg.hpp diff --git a/include/integratorxx/quadratures/radial.hpp b/include/integratorxx/quadratures/radial.hpp index 69a9cca..cebea4d 100644 --- a/include/integratorxx/quadratures/radial.hpp +++ b/include/integratorxx/quadratures/radial.hpp @@ -5,3 +5,4 @@ #include #include #include +#include diff --git a/include/integratorxx/quadratures/radial/lmg.hpp b/include/integratorxx/quadratures/radial/lmg.hpp new file mode 100644 index 0000000..8a37800 --- /dev/null +++ b/include/integratorxx/quadratures/radial/lmg.hpp @@ -0,0 +1,54 @@ +#include +#include + +namespace IntegratorXX { +namespace lmg { + +inline double r_upper_obj(int m, double alpha, double r) { + const double am_term = (m + 1.0) / 2.0; + const double g_term = std::tgamma((m + 3.0) / 2.0); + const double x = alpha * r * r; + return g_term * std::pow(x, am_term) * std::exp(-x); +} + +// Equation 19 in the LMG paper +inline double r_upper(int m, double alpha, double prec) { + // P = G * (ALPHA * R^2)^((M+1)/2) * EXP(-ALPHA * R^2) + // P = G * X^L * EXP(-X): X = ALPHA * R^2, L = (M+1)/2 + // X = -L * LAMBERT_W( - (P/G)^(1/L) / L ) + // R = SQRT(X / ALPHA) + const double am_term = (m + 1.0) / 2.0; + const double g_term = std::tgamma((m + 3.0) / 2.0); + const double arg = std::pow(prec/g_term, 1.0 / am_term) / am_term; + const double wval = lambert_wm1(-arg); // W_(-1) is the larger value here + const double x = -am_term * wval; + const double r = std::sqrt(x / alpha); + return r; +} + +inline double r_lower(int m, double alpha, double prec) { + // Magic numbers from the LMG paper below Eq. 25 + double Dm; + switch(m) { + case -2: + Dm = 9.1; + break; + case 0: + Dm = 1.9; + break; + case 2: + Dm = -1.0; + break; + case 4: + Dm = -2.3; + break; + default: + Dm = NAN; + } + + return std::exp(1.0 / (m + 3.0) * (Dm - std::log(1.0 / prec))) / + std::sqrt(alpha); +} + +} +} diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index bae823a..394e4b1 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -389,4 +389,66 @@ TEST_CASE("LMG", "[1d-quad]") { } + SECTION("LMG Radial Bounds") { + + // aug-cc-pVTZ Hydrogen Data + SECTION("H aug-cc-pVTZ") { + const double alpha_min_0 = 0.02526; + const double alpha_min_1 = 0.10200; + const double alpha_min_2 = 0.24700; + + const double alpha_max_0 = 33.87; + //const double alpha_max_1 = 1.407; + const double alpha_max_2 = 1.057; + + auto r_lower_0 = lmg::r_lower(0, alpha_max_0, 1e-12); + auto r_lower_2 = lmg::r_lower(2, alpha_max_2, 1e-12); + + auto r_upper_0 = lmg::r_upper(0, alpha_min_0, 1e-12); + auto r_upper_1 = lmg::r_upper(1, alpha_min_1, 1e-12); + auto r_upper_2 = lmg::r_upper(2, alpha_min_2, 1e-12); + + + REQUIRE_THAT(r_lower_0, Catch::Matchers::WithinAbs(3.2370214224347012e-05, 1e-15)); + REQUIRE_THAT(r_lower_2, Catch::Matchers::WithinAbs(3.1703237193762383e-03, 1e-15)); + + REQUIRE_THAT(r_upper_0, Catch::Matchers::WithinAbs(3.3998088412766542e+01, 1e-15)); + REQUIRE_THAT(r_upper_1, Catch::Matchers::WithinAbs(1.7452224096146878e+01, 1e-15)); + REQUIRE_THAT(r_upper_2, Catch::Matchers::WithinAbs(1.1588086549342464e+01, 1e-15)); + } + + // aug-cc-pVTZ Carbon Data + SECTION("C aug-cc-pVTZ") { + const double alpha_min_0 = 3.3423e-02; + const double alpha_min_1 = 3.3258e-02; + const double alpha_min_2 = 7.6815e-02; + const double alpha_min_3 = 1.7636e-01;; + const double alpha_min_4 = 5.1329e-01; + + const double alpha_max_0 = 6.177194e+06; + const double alpha_max_2 = 8.857680e+01; + const double alpha_max_4 = 1.445900e+00; + + auto r_lower_0 = lmg::r_lower(0, alpha_max_0, 1e-12); + auto r_lower_2 = lmg::r_lower(2, alpha_max_2, 1e-12); + auto r_lower_4 = lmg::r_lower(4, alpha_max_4, 1e-12); + + auto r_upper_0 = lmg::r_upper(0, alpha_min_0, 1e-12); + auto r_upper_1 = lmg::r_upper(1, alpha_min_1, 1e-12); + auto r_upper_2 = lmg::r_upper(2, alpha_min_2, 1e-12); + auto r_upper_3 = lmg::r_upper(3, alpha_min_3, 1e-12); + auto r_upper_4 = lmg::r_upper(4, alpha_min_4, 1e-12); + + REQUIRE_THAT(r_lower_0, Catch::Matchers::WithinAbs(7.5797965945427875e-08,1e-15)); + REQUIRE_THAT(r_lower_2, Catch::Matchers::WithinAbs(3.4632282095817583e-04,1e-15)); + REQUIRE_THAT(r_lower_4, Catch::Matchers::WithinAbs(1.1559748847781975e-02,1e-15)); + + REQUIRE_THAT(r_upper_0, Catch::Matchers::WithinAbs(2.9556190534662171e+01,1e-15)); + REQUIRE_THAT(r_upper_1, Catch::Matchers::WithinAbs(3.0563480016561439e+01,1e-15)); + REQUIRE_THAT(r_upper_2, Catch::Matchers::WithinAbs(2.0779600292386970e+01,1e-15)); + REQUIRE_THAT(r_upper_3, Catch::Matchers::WithinAbs(1.4179981794582680e+01,1e-15)); + REQUIRE_THAT(r_upper_4, Catch::Matchers::WithinAbs(8.5952200832939525e+00,1e-15)); + } + } + } From 54f13f422283d058327c684b43f338768eadc53d Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 27 Sep 2023 17:06:43 -0700 Subject: [PATCH 04/19] Added LMG utility to compute step size for radial grid of a given precision --- .../integratorxx/quadratures/radial/lmg.hpp | 20 +++++++++++++++++++ test/1d_quadratures.cxx | 10 +++++++++- 2 files changed, 29 insertions(+), 1 deletion(-) diff --git a/include/integratorxx/quadratures/radial/lmg.hpp b/include/integratorxx/quadratures/radial/lmg.hpp index 8a37800..24fbf21 100644 --- a/include/integratorxx/quadratures/radial/lmg.hpp +++ b/include/integratorxx/quadratures/radial/lmg.hpp @@ -50,5 +50,25 @@ inline double r_lower(int m, double alpha, double prec) { std::sqrt(alpha); } +inline double step_size(int m, double prec) { + + // Recast Eqs 17/18 into the form + // R == C * x^L * EXP[-PI/2 * X] with X == PI/h + // X = - 2*L / PI * LAMBERT_W( -PI/(2*L) * (R/C)^(1/L) ) + double C = 4 * M_SQRT2; // This is the constant in Eq 17 + double L; + if( m == 0 ) { + L = 1; + } else { + L = m / 2.0 + 1.0; // Eq 17 -> Eq 18 + C = C * std::tgamma(1.5) / std::tgamma((m + 3.0) / 2.0); + } + + const double L_FAC = M_PI_2 / L; + const double X = -(1.0 / L_FAC) * lambert_wm1( -L_FAC * std::pow(prec/C, 1/L)); + return M_PI / X; + +} + } } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 394e4b1..b30e1c9 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -390,7 +390,6 @@ TEST_CASE("LMG", "[1d-quad]") { } SECTION("LMG Radial Bounds") { - // aug-cc-pVTZ Hydrogen Data SECTION("H aug-cc-pVTZ") { const double alpha_min_0 = 0.02526; @@ -451,4 +450,13 @@ TEST_CASE("LMG", "[1d-quad]") { } } + + SECTION("LMG Step Size") { + REQUIRE_THAT(lmg::step_size(0, 1e-12), Catch::Matchers::WithinAbs(1.5235497564176387e-01,1e-15)); + REQUIRE_THAT(lmg::step_size(1, 1e-12), Catch::Matchers::WithinAbs(1.4579053103691331e-01,1e-15)); + REQUIRE_THAT(lmg::step_size(2, 1e-12), Catch::Matchers::WithinAbs(1.4028895246473957e-01,1e-15)); + REQUIRE_THAT(lmg::step_size(3, 1e-12), Catch::Matchers::WithinAbs(1.3554180272204153e-01,1e-15)); + REQUIRE_THAT(lmg::step_size(4, 1e-12), Catch::Matchers::WithinAbs(1.3136472924710046e-01,1e-15)); + } + } From 26f9c631e165a74624d83c92c24f678445b4874b Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 27 Sep 2023 17:13:42 -0700 Subject: [PATCH 05/19] Remove boost header --- test/1d_quadratures.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index b30e1c9..ec078e1 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -1,7 +1,6 @@ #include "catch2/catch_all.hpp" #include #include -#include #include #include #include From 7b71e32146791f3587eb34bd6dd2fd69e01587bd Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 28 Sep 2023 14:30:43 -0700 Subject: [PATCH 06/19] Added parameter fall through for lmg::r_lower per source provided by Roland --- include/integratorxx/quadratures/radial/lmg.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/integratorxx/quadratures/radial/lmg.hpp b/include/integratorxx/quadratures/radial/lmg.hpp index 24fbf21..809911c 100644 --- a/include/integratorxx/quadratures/radial/lmg.hpp +++ b/include/integratorxx/quadratures/radial/lmg.hpp @@ -43,7 +43,7 @@ inline double r_lower(int m, double alpha, double prec) { Dm = -2.3; break; default: - Dm = NAN; + Dm = 4.0; // Not in the paper, taken from soruce provided by Roland Lindh } return std::exp(1.0 / (m + 3.0) * (Dm - std::log(1.0 / prec))) / From 7e71f9469cd32ebc2af3e749024713e2016959b8 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Mon, 2 Oct 2023 13:42:52 -0700 Subject: [PATCH 07/19] Start LMG RadialTraits + misc dox++ --- .../integratorxx/quadratures/radial/lmg.hpp | 46 ++++++++++++++----- 1 file changed, 35 insertions(+), 11 deletions(-) diff --git a/include/integratorxx/quadratures/radial/lmg.hpp b/include/integratorxx/quadratures/radial/lmg.hpp index 809911c..837a6f4 100644 --- a/include/integratorxx/quadratures/radial/lmg.hpp +++ b/include/integratorxx/quadratures/radial/lmg.hpp @@ -4,6 +4,7 @@ namespace IntegratorXX { namespace lmg { +// Eq 19 of LMG paper (DOI 10.1007/s002140100263) inline double r_upper_obj(int m, double alpha, double r) { const double am_term = (m + 1.0) / 2.0; const double g_term = std::tgamma((m + 3.0) / 2.0); @@ -11,7 +12,7 @@ inline double r_upper_obj(int m, double alpha, double r) { return g_term * std::pow(x, am_term) * std::exp(-x); } -// Equation 19 in the LMG paper +// Solve Eq 19 of LMG paper (DOI 10.1007/s002140100263) inline double r_upper(int m, double alpha, double prec) { // P = G * (ALPHA * R^2)^((M+1)/2) * EXP(-ALPHA * R^2) // P = G * X^L * EXP(-X): X = ALPHA * R^2, L = (M+1)/2 @@ -26,8 +27,9 @@ inline double r_upper(int m, double alpha, double prec) { return r; } +// Solve Eq 25 of LMG paper (DOI 10.1007/s002140100263) inline double r_lower(int m, double alpha, double prec) { - // Magic numbers from the LMG paper below Eq. 25 + // Magic numbers below Eq. 25 double Dm; switch(m) { case -2: @@ -43,26 +45,22 @@ inline double r_lower(int m, double alpha, double prec) { Dm = -2.3; break; default: - Dm = 4.0; // Not in the paper, taken from soruce provided by Roland Lindh + // Not in the paper, taken from source provided by Roland Lindh + Dm = 4.0; } return std::exp(1.0 / (m + 3.0) * (Dm - std::log(1.0 / prec))) / std::sqrt(alpha); } +// Eqs 17 and 18 of the LMG paper (DOI 10.1007/s002140100263) inline double step_size(int m, double prec) { // Recast Eqs 17/18 into the form // R == C * x^L * EXP[-PI/2 * X] with X == PI/h // X = - 2*L / PI * LAMBERT_W( -PI/(2*L) * (R/C)^(1/L) ) - double C = 4 * M_SQRT2; // This is the constant in Eq 17 - double L; - if( m == 0 ) { - L = 1; - } else { - L = m / 2.0 + 1.0; // Eq 17 -> Eq 18 - C = C * std::tgamma(1.5) / std::tgamma((m + 3.0) / 2.0); - } + double C = 4 * M_SQRT2 * std::tgamma(1.5) / std::tgamma(m / 2.0 + 1.5); + double L = m / 2.0 + 1.0; // Eq 17 -> Eq 18 const double L_FAC = M_PI_2 / L; const double X = -(1.0 / L_FAC) * lambert_wm1( -L_FAC * std::pow(prec/C, 1/L)); @@ -71,4 +69,30 @@ inline double step_size(int m, double prec) { } } + + +class LindhMalmqvistGagliardiRadialTraits { + double r_min_; + double r_max_; + double step_size_; + double c_; + +public: + + LindhMalmqvistGagliardiRadialTraits(double r_min, double r_max, int l_max, double prec) : + r_min_(r_min), r_max_(r_max), + // H is determined by L_MAX + step_size_(lmg::step_size(2*l_max, prec)), + c_(r_min_ / (std::exp(step_size_) - 1.0)) { } + + template + inline auto radial_transform(PointType x) const noexcept { + return c_ * (std::exp(x) - 1.0); + } + + template + inline auto radial_jacobian(PointType x) const noexcept { + return c_ * std::exp(x); + } +}; } From 6c41dde7e7dc0cb73ba6efffb2f1d29fd156281b Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Mon, 2 Oct 2023 17:05:11 -0700 Subject: [PATCH 08/19] Added initial implementation of LMG + UTs --- .../integratorxx/quadratures/radial/lmg.hpp | 43 ++++++++++++++++--- .../quadratures/radial/radial_transform.hpp | 29 +++++++++++++ test/1d_quadratures.cxx | 27 ++++++++++++ 3 files changed, 92 insertions(+), 7 deletions(-) diff --git a/include/integratorxx/quadratures/radial/lmg.hpp b/include/integratorxx/quadratures/radial/lmg.hpp index 837a6f4..1a67e51 100644 --- a/include/integratorxx/quadratures/radial/lmg.hpp +++ b/include/integratorxx/quadratures/radial/lmg.hpp @@ -72,18 +72,19 @@ inline double step_size(int m, double prec) { class LindhMalmqvistGagliardiRadialTraits { - double r_min_; - double r_max_; + using self_type = LindhMalmqvistGagliardiRadialTraits; + double step_size_; double c_; public: - LindhMalmqvistGagliardiRadialTraits(double r_min, double r_max, int l_max, double prec) : - r_min_(r_min), r_max_(r_max), - // H is determined by L_MAX - step_size_(lmg::step_size(2*l_max, prec)), - c_(r_min_ / (std::exp(step_size_) - 1.0)) { } + LindhMalmqvistGagliardiRadialTraits(double c, double step_size) : + step_size_(step_size), c_(c) { } + + LindhMalmqvistGagliardiRadialTraits(const self_type&) = default; + self_type& operator=(const self_type&) = default; + template inline auto radial_transform(PointType x) const noexcept { @@ -94,5 +95,33 @@ class LindhMalmqvistGagliardiRadialTraits { inline auto radial_jacobian(PointType x) const noexcept { return c_ * std::exp(x); } + + template + std::enable_if_t< + // LMG really only makes sense with a Trapezoid grid + std::is_same_v< + BaseQuad, + UniformTrapezoid + > + > preprocess_base_quad(typename BaseQuad::point_container& points, + typename BaseQuad::weight_container& weights) const { + + using point_type = typename BaseQuad::point_type; + assert(points.size() > 2); + const auto npts = points.size() - 2; + point_type up = npts * step_size_; + + // Scale points and weights + for(auto& x : points ) x *= step_size_ * (npts + 1); + for(auto& w : weights) w = step_size_; + } + }; + +template +using LindhMalmqvistGagliardi = RadialTransformQuadrature< + UniformTrapezoid, LindhMalmqvistGagliardiRadialTraits +>; + } diff --git a/include/integratorxx/quadratures/radial/radial_transform.hpp b/include/integratorxx/quadratures/radial/radial_transform.hpp index 93f8096..0d44848 100644 --- a/include/integratorxx/quadratures/radial/radial_transform.hpp +++ b/include/integratorxx/quadratures/radial/radial_transform.hpp @@ -4,6 +4,32 @@ namespace IntegratorXX { +namespace detail { + +template +class has_preprocess_base_quad { + +using point_container_ref = typename BaseQuad::point_container&; +using weight_container_ref = typename BaseQuad::weight_container&; + +template (). + template preprocess_base_quad( + std::declval(), + std::declval() + ))> +static std::true_type test(int); + +template static std::false_type test(...); + +public: + +static constexpr bool value = decltype(test(0))::value; + +}; + +} + template struct RadialTransformQuadrature : public Quadrature> { @@ -45,6 +71,9 @@ struct quadrature_traits> { const auto npts_base = base_quad_traits::bound_inclusive ? npts+2 : npts; auto [base_x, base_w] = base_quad_traits::generate(npts_base); + if constexpr (detail::has_preprocess_base_quad::value) + traits.template preprocess_base_quad(base_x, base_w); + const auto ipts_offset = !!base_quad_traits::bound_inclusive; for(size_t i = 0; i < npts; ++i) { const auto xi = base_x[i + ipts_offset]; diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index ec078e1..8ad087d 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -458,4 +458,31 @@ TEST_CASE("LMG", "[1d-quad]") { REQUIRE_THAT(lmg::step_size(4, 1e-12), Catch::Matchers::WithinAbs(1.3136472924710046e-01,1e-15)); } + SECTION("LMG Integrate") { + double alpha_min = 0.3203; + double alpha_max = 80.0; + + // Reference values from OpenMolcas + double r_lower = 2.1062415621095990e-005; + double r_upper = 9.9559233798449043; + double h = 0.15235041341729277; + + double c = r_lower / (std::exp(h) - 1); + int n = std::log(1.0 + r_upper/c) / h + 1.0; + printf("ALPHA_MIN = %.4e ALPHA_MAX = %.4e H = %.4e C = %.6e N = %d\n", + alpha_min, alpha_max, h, c, n); + + + LindhMalmqvistGagliardiRadialTraits traits(c, h); + LindhMalmqvistGagliardi quad(n, traits); + for(int i = 0; i < n; ++i) { + double x = (i+1) * h; + double r = c * (std::exp(x) - 1.0); + double w = h*(r + c); + printf("%d R = %.16e W = %.16e \n", i, r - quad.points()[i], w-quad.weights()[i]); + } + + + } + } From cb6a4f9987d4fbbb318ce3932ff00b1271447fe6 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 12 Oct 2023 11:57:03 -0700 Subject: [PATCH 09/19] Added implementations of integer and half-integer math functions --- include/integratorxx/util/constants.hpp | 7 ++ include/integratorxx/util/factorial.hpp | 72 +++++++++++++++++ include/integratorxx/util/gamma.hpp | 35 +++++++++ include/integratorxx/util/pow_2.hpp | 20 +++++ include/integratorxx/util/type_traits.hpp | 20 +++++ test/CMakeLists.txt | 3 + test/util.cxx | 95 +++++++++++++++++++++++ 7 files changed, 252 insertions(+) create mode 100644 include/integratorxx/util/constants.hpp create mode 100644 include/integratorxx/util/factorial.hpp create mode 100644 include/integratorxx/util/gamma.hpp create mode 100644 include/integratorxx/util/pow_2.hpp create mode 100644 include/integratorxx/util/type_traits.hpp create mode 100644 test/util.cxx diff --git a/include/integratorxx/util/constants.hpp b/include/integratorxx/util/constants.hpp new file mode 100644 index 0000000..2257d1b --- /dev/null +++ b/include/integratorxx/util/constants.hpp @@ -0,0 +1,7 @@ +#pragma once + +namespace IntegratorXX::constants { + template + inline static constexpr T sqrt_pi_v = + 1.7724538509055160272981674833411451827975494561223871282138L; +} diff --git a/include/integratorxx/util/factorial.hpp b/include/integratorxx/util/factorial.hpp new file mode 100644 index 0000000..7c03f04 --- /dev/null +++ b/include/integratorxx/util/factorial.hpp @@ -0,0 +1,72 @@ +#pragma once +#include +#include + +#include + +namespace IntegratorXX { + +namespace detail { + +// Factorials for [0,20] +static constexpr std::array cached_factorials = { + 1ul, 1ul, 2ul, 6ul, 24ul, 120ul, 720ul, 5040ul, 40320ul, 362880ul, 3628800ul, 39916800ul, + 479001600ul, 6227020800ul, 87178291200ul, 1307674368000ul, 20922789888000ul, + 355687428096000ul, 6402373705728000ul, 121645100408832000ul, 2432902008176640000ul +}; + +// Double factorials for [0,32] +static constexpr std::array cached_double_factorials = { + 1ul, 1ul, 2ul, 3ul, 8ul, 15ul, 48ul, 105ul, 384ul, 945ul, 3840ul, 10395ul, 46080ul, 135135ul, + 645120ul, 2027025ul, 10321920ul, 34459425ul, 185794560ul, 654729075ul, 3715891200ul, + 13749310575ul, 81749606400ul, 316234143225ul, 1961990553600ul, 7905853580625ul, + 51011754393600ul, 213458046676875ul, 1428329123020800ul, 6190283353629375ul, + 42849873690624000ul, 191898783962510625ul, 1371195958099968000ul +}; + +} + +template +detail::enable_if_integral_t + factorial(IntegralType n) { + + assert(n >= 0); + if( n < detail::cached_factorials.size() ) { + auto val = detail::cached_factorials[n]; + assert(val <= std::numeric_limits::max()); + return val; + } else { + auto val = n * factorial(n-1); + assert(val <= std::numeric_limits::max()); + return val; + } + +} + +template +detail::enable_if_integral_t + double_factorial(IntegralType n) { + + // If we're signed, check negative arguments + if(std::is_signed_v) { + // Only provide values for which n!! is integral + if(n == -1) return 1; + if(n == -3) return -1; + } + + // This also handles the case when n < 0 is undefied or non-integral + assert(n >= 0); + + if( n < detail::cached_double_factorials.size() ) { + auto val = detail::cached_double_factorials[n]; + assert(val <= std::numeric_limits::max()); + return val; + } else { + auto val = n * double_factorial(n-2); + assert(val <= std::numeric_limits::max()); + return val; + } + +} + +} diff --git a/include/integratorxx/util/gamma.hpp b/include/integratorxx/util/gamma.hpp new file mode 100644 index 0000000..568695c --- /dev/null +++ b/include/integratorxx/util/gamma.hpp @@ -0,0 +1,35 @@ +#pragma once +#include +#include +#include +#include + +namespace IntegratorXX { + +/// GAMMA(n) = (n-1)! +template +FloatType integer_tgamma(IntegralType n) { + if(n <= 0) return std::numeric_limits::quiet_NaN(); + else { + auto val = factorial(n-1); + assert(val <= std::numeric_limits::max()); + return val; + } +} + +/// GAMMA(n/2) = +/// CASE n even -> INTEGER_GAMMA(n/2) +/// CASE n odd -> SQRT_PI / POW_2((n-1)/2) * (n-2)!! +template +FloatType half_integer_tgamma(IntegralType n) { + if(n%2 == 0) return integer_tgamma(n/2); + else { + constexpr auto c = constants::sqrt_pi_v; + const auto df = double_factorial(n-2); + const auto pt = half_integer_pow_2(n-1); + return c * df / pt; + } +} + + +} diff --git a/include/integratorxx/util/pow_2.hpp b/include/integratorxx/util/pow_2.hpp new file mode 100644 index 0000000..e59f101 --- /dev/null +++ b/include/integratorxx/util/pow_2.hpp @@ -0,0 +1,20 @@ +#pragma once +#include + +namespace IntegratorXX { + +template +detail::enable_if_integral_t + integer_pow_2(IntegralType n) { + assert(n >= 0); + return IntegralType(1) << n; +} + +template +FloatType half_integer_pow_2(IntegralType n) { + if(n < 0) return FloatType(1) / half_integer_pow_2(-n); + else if(n%2 == 0) return integer_pow_2(n/2); + else return M_SQRT2 * integer_pow_2(n/2); +} + +} diff --git a/include/integratorxx/util/type_traits.hpp b/include/integratorxx/util/type_traits.hpp new file mode 100644 index 0000000..8b31592 --- /dev/null +++ b/include/integratorxx/util/type_traits.hpp @@ -0,0 +1,20 @@ +#pragma once +#include + +namespace IntegratorXX::detail { + +template +using enable_if_integral = std::enable_if,U>; + +template +using enable_if_integral_t = typename enable_if_integral::type; + +template +using enable_if_float_and_integral = + std::enable_if and std::is_integral_v, V>; + +template +using enable_if_float_and_integral_t = + typename enable_if_float_and_integral::type; + +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 865e52e..424784b 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -21,6 +21,9 @@ endif() add_library(integratorxx_common_ut quad_matcher.cxx) target_link_libraries( integratorxx_common_ut PUBLIC Catch2::Catch2WithMain integratorxx ) +add_executable( util_test util.cxx ) +target_link_libraries( util_test PUBLIC integratorxx_common_ut ) + add_executable( quadrature_manipulation quadrature_manipulation.cxx ) target_link_libraries( quadrature_manipulation PUBLIC integratorxx_common_ut ) diff --git a/test/util.cxx b/test/util.cxx new file mode 100644 index 0000000..7165247 --- /dev/null +++ b/test/util.cxx @@ -0,0 +1,95 @@ +#include "catch2/catch_all.hpp" +#include +#include +#include +#include + +uint64_t dumb_factorial(uint64_t n) { + if( n == 0 or n == 1 ) return 1; + else return n * dumb_factorial(n-1); +} + +uint64_t dumb_double_factorial(uint64_t n) { + if( n == 0 or n == 1 ) return 1; + else return n * dumb_double_factorial(n-2); +} + + +TEST_CASE("Factorial", "[util]") { + using namespace IntegratorXX; + + REQUIRE(std::is_same_v< + decltype(factorial(std::declval())), int>); + REQUIRE(std::is_same_v< + decltype(factorial(std::declval())), unsigned>); + + for(uint64_t i = 0; i <= 20; ++i) { + REQUIRE(factorial(i) == dumb_factorial(i)); + } + + REQUIRE(double_factorial(-1) == 1); + REQUIRE(double_factorial(-3) == -1); + for(uint64_t i = 0; i <= 32; ++i) { + REQUIRE(double_factorial(i) == dumb_double_factorial(i)); + } +} + +TEST_CASE("Pow 2", "[util]") { + using namespace IntegratorXX; + + REQUIRE(std::is_same_v< + decltype(integer_pow_2(std::declval())), int>); + REQUIRE(std::is_same_v< + decltype(integer_pow_2(std::declval())), unsigned>); + + REQUIRE(std::is_same_v< + decltype(half_integer_pow_2(std::declval())), float>); + REQUIRE(std::is_same_v< + decltype(half_integer_pow_2(std::declval())), float>); + REQUIRE(std::is_same_v< + decltype(half_integer_pow_2(std::declval())), double>); + REQUIRE(std::is_same_v< + decltype(half_integer_pow_2(std::declval())), double>); + + for(uint64_t i = 0; i < 64; ++i) { + REQUIRE(integer_pow_2(i) == (1ul << i)); + } + + // 126 to avoid INT64_T overflow + for(int64_t i = 0; i < 126; i += 2) { + const auto v = half_integer_pow_2(i); + REQUIRE( v == (1ul << i/2)); + REQUIRE(half_integer_pow_2(-i) == 1.0/v); + REQUIRE(half_integer_pow_2(i+1) == M_SQRT2*v); + REQUIRE(half_integer_pow_2(-(i+1)) == 1.0/(M_SQRT2*v)); + } + +} + +TEST_CASE("Gamma", "[util]") { + using namespace IntegratorXX; + + REQUIRE(std::is_same_v< + decltype(integer_tgamma(std::declval())), float>); + REQUIRE(std::is_same_v< + decltype(integer_tgamma(std::declval())), float>); + REQUIRE(std::is_same_v< + decltype(integer_tgamma(std::declval())), double>); + REQUIRE(std::is_same_v< + decltype(integer_tgamma(std::declval())), double>); + + REQUIRE(std::isnan(integer_tgamma(-1))); + REQUIRE(std::isnan(integer_tgamma(0))); + for(uint64_t i = 1; i <= 21; ++i) { + REQUIRE(integer_tgamma(i) == factorial(i-1)); + } + + REQUIRE(std::isnan(half_integer_tgamma(0))); + for(int64_t i = 1; i <= 34; ++i) { + const auto v = half_integer_tgamma(i); + const auto r = std::tgamma(i / 2.0l); + REQUIRE_THAT(v, Catch::Matchers::WithinAbs(r, 1e-20)); + } + +} + From bbb2483f39723e6b4d864f5a9c205174b10a28a7 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 12 Oct 2023 12:01:45 -0700 Subject: [PATCH 10/19] Replace tgamma with half_integer_tgamma where appropriate --- include/integratorxx/quadratures/radial/lmg.hpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/include/integratorxx/quadratures/radial/lmg.hpp b/include/integratorxx/quadratures/radial/lmg.hpp index 1a67e51..e36673c 100644 --- a/include/integratorxx/quadratures/radial/lmg.hpp +++ b/include/integratorxx/quadratures/radial/lmg.hpp @@ -1,4 +1,5 @@ #include +#include #include namespace IntegratorXX { @@ -7,7 +8,8 @@ namespace lmg { // Eq 19 of LMG paper (DOI 10.1007/s002140100263) inline double r_upper_obj(int m, double alpha, double r) { const double am_term = (m + 1.0) / 2.0; - const double g_term = std::tgamma((m + 3.0) / 2.0); + //const double g_term = std::tgamma((m + 3.0) / 2.0); + const double g_term = half_integer_tgamma(m + 3); const double x = alpha * r * r; return g_term * std::pow(x, am_term) * std::exp(-x); } @@ -19,7 +21,8 @@ inline double r_upper(int m, double alpha, double prec) { // X = -L * LAMBERT_W( - (P/G)^(1/L) / L ) // R = SQRT(X / ALPHA) const double am_term = (m + 1.0) / 2.0; - const double g_term = std::tgamma((m + 3.0) / 2.0); + //const double g_term = std::tgamma((m + 3.0) / 2.0); + const double g_term = half_integer_tgamma(m + 3); const double arg = std::pow(prec/g_term, 1.0 / am_term) / am_term; const double wval = lambert_wm1(-arg); // W_(-1) is the larger value here const double x = -am_term * wval; @@ -59,7 +62,7 @@ inline double step_size(int m, double prec) { // Recast Eqs 17/18 into the form // R == C * x^L * EXP[-PI/2 * X] with X == PI/h // X = - 2*L / PI * LAMBERT_W( -PI/(2*L) * (R/C)^(1/L) ) - double C = 4 * M_SQRT2 * std::tgamma(1.5) / std::tgamma(m / 2.0 + 1.5); + double C = 4 * M_SQRT2 * half_integer_tgamma(3) / half_integer_tgamma(m + 3); double L = m / 2.0 + 1.0; // Eq 17 -> Eq 18 const double L_FAC = M_PI_2 / L; From 6495dc591551b66ab78648d35fb25a10a9282289 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 12 Oct 2023 12:45:02 -0700 Subject: [PATCH 11/19] Added integer POW routines --- .../integratorxx/quadratures/radial/lmg.hpp | 6 ++--- include/integratorxx/util/pow.hpp | 23 +++++++++++++++++++ test/util.cxx | 12 ++++++++++ 3 files changed, 37 insertions(+), 4 deletions(-) create mode 100644 include/integratorxx/util/pow.hpp diff --git a/include/integratorxx/quadratures/radial/lmg.hpp b/include/integratorxx/quadratures/radial/lmg.hpp index e36673c..d784ba5 100644 --- a/include/integratorxx/quadratures/radial/lmg.hpp +++ b/include/integratorxx/quadratures/radial/lmg.hpp @@ -1,5 +1,6 @@ #include #include +#include #include namespace IntegratorXX { @@ -7,11 +8,9 @@ namespace lmg { // Eq 19 of LMG paper (DOI 10.1007/s002140100263) inline double r_upper_obj(int m, double alpha, double r) { - const double am_term = (m + 1.0) / 2.0; - //const double g_term = std::tgamma((m + 3.0) / 2.0); const double g_term = half_integer_tgamma(m + 3); const double x = alpha * r * r; - return g_term * std::pow(x, am_term) * std::exp(-x); + return g_term * half_integer_pow(x, m+1) * std::exp(-x); } // Solve Eq 19 of LMG paper (DOI 10.1007/s002140100263) @@ -21,7 +20,6 @@ inline double r_upper(int m, double alpha, double prec) { // X = -L * LAMBERT_W( - (P/G)^(1/L) / L ) // R = SQRT(X / ALPHA) const double am_term = (m + 1.0) / 2.0; - //const double g_term = std::tgamma((m + 3.0) / 2.0); const double g_term = half_integer_tgamma(m + 3); const double arg = std::pow(prec/g_term, 1.0 / am_term) / am_term; const double wval = lambert_wm1(-arg); // W_(-1) is the larger value here diff --git a/include/integratorxx/util/pow.hpp b/include/integratorxx/util/pow.hpp new file mode 100644 index 0000000..c6ae500 --- /dev/null +++ b/include/integratorxx/util/pow.hpp @@ -0,0 +1,23 @@ +#pragma once +#include +#include + +namespace IntegratorXX { + +template +FloatType integer_pow(FloatType x, IntegralType n) { + if(n == 0) return 1; + if(n == 1) return x; + FloatType v = x; + for(IntegralType i = 1; i < n; ++i) v *= x; + return v; +} + +template +FloatType half_integer_pow(FloatType x, IntegralType n) { + assert(n >= 0); + if(n%2 == 0) return integer_pow(x, n/2); + else return std::sqrt(x) * integer_pow(x, n/2); +} + +} diff --git a/test/util.cxx b/test/util.cxx index 7165247..b051fea 100644 --- a/test/util.cxx +++ b/test/util.cxx @@ -1,6 +1,7 @@ #include "catch2/catch_all.hpp" #include #include +#include #include #include @@ -34,6 +35,17 @@ TEST_CASE("Factorial", "[util]") { } } +TEST_CASE("Pow", "[util]") { + + using namespace IntegratorXX; + using namespace Catch::Matchers; + for(int i = 0; i < 10; ++i) { + REQUIRE_THAT(integer_pow(2.0, i), WithinAbs(std::pow(2.0,i), 1e-16)); + REQUIRE_THAT(half_integer_pow(2.0, i), WithinAbs(std::pow(2.0,i/2.0), 1e-16)); + } + +} + TEST_CASE("Pow 2", "[util]") { using namespace IntegratorXX; From 2ef9ca36c1f7d6d060a8a9199df1f155cdb9c3e8 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 1 Nov 2023 11:19:29 -0700 Subject: [PATCH 12/19] function wrappers --- .../integratorxx/quadratures/radial/lmg.hpp | 6 +- include/integratorxx/util/factorial.hpp | 63 +++++++++++++++++++ test/1d_quadratures.cxx | 13 +++- 3 files changed, 77 insertions(+), 5 deletions(-) create mode 100644 include/integratorxx/util/factorial.hpp diff --git a/include/integratorxx/quadratures/radial/lmg.hpp b/include/integratorxx/quadratures/radial/lmg.hpp index 1a67e51..ace3c3d 100644 --- a/include/integratorxx/quadratures/radial/lmg.hpp +++ b/include/integratorxx/quadratures/radial/lmg.hpp @@ -1,4 +1,5 @@ #include +#include #include namespace IntegratorXX { @@ -19,8 +20,9 @@ inline double r_upper(int m, double alpha, double prec) { // X = -L * LAMBERT_W( - (P/G)^(1/L) / L ) // R = SQRT(X / ALPHA) const double am_term = (m + 1.0) / 2.0; - const double g_term = std::tgamma((m + 3.0) / 2.0); - const double arg = std::pow(prec/g_term, 1.0 / am_term) / am_term; + //const double g_term = std::tgamma((m + 3.0) / 2.0); + const double g_term = (m%2) ? integral_tgamma((m+3)/2) : half_integral_tgamma(m + 3); + const double arg = std::pow(prec/g_term, 0.0 / am_term) / am_term; const double wval = lambert_wm1(-arg); // W_(-1) is the larger value here const double x = -am_term * wval; const double r = std::sqrt(x / alpha); diff --git a/include/integratorxx/util/factorial.hpp b/include/integratorxx/util/factorial.hpp new file mode 100644 index 0000000..6d47c88 --- /dev/null +++ b/include/integratorxx/util/factorial.hpp @@ -0,0 +1,63 @@ +#pragma once + +namespace IntegratorXX { + +namespace detail { + +constexpr static std::array factorials = { + 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, + 39916800, 479001600, 6227020800, 87178291200, 1307674368000, + 20922789888000, 355687428096000, 6402373705728000, 121645100408832000 +}; + +constexpr static std::array double_factorials = { + 0, 1, 1, 3, 6, 15, 38, 105, 306, 945, 3063, 10395, + 36766, 135135, 514731, 2027025, 8235700, 34459425, + 148242610, 654729075 +}; + + + +} + +template +T factorial(unsigned n) { + if( n < 20 ) { + return detail::factorials[n]; + } else { + return n * factorial(n-1); + } +} + +template +T double_factorial(int n) { + if(n == -3) return -1; + else if(n == -1) return 1; + else if( n < 20 ) { + return detail::double_factorials[size_t(n)]; + } else { + return n * double_factorial(n-2); + } +} + +inline double pow_2(int n) { + return 1ul << n; +} + +// GAMMA(n/2) +template +T half_integral_tgamma(int n) { + assert(n%2); + if(n == 1) return std::sqrt(M_PI); + if(n == 3) return 0.5 * std::sqrt(M_PI); + return double_factorial(n-2) / pow_2(n-3) / M_2_SQRTPI; +} + +template +T integral_tgamma(T n) { + std::cout << n << std::endl; + return factorial(n-1); +} + + +} diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 8ad087d..f8c0001 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -257,6 +257,13 @@ TEST_CASE( "Delley", "[1d-quad]" ) { TEST_CASE("LMG", "[1d-quad]") { using namespace IntegratorXX; + using namespace Catch::Matchers; + + SECTION("HALF INTEGRAL GAMMA") { + REQUIRE_THAT(half_integral_tgamma(1), WithinAbs(std::sqrt(M_PI), 1e-15)); + REQUIRE_THAT(half_integral_tgamma(3), WithinAbs(std::sqrt(M_PI)/2,1e-15)); + } + SECTION("Lambert W") { size_t n = 50; @@ -469,8 +476,8 @@ TEST_CASE("LMG", "[1d-quad]") { double c = r_lower / (std::exp(h) - 1); int n = std::log(1.0 + r_upper/c) / h + 1.0; - printf("ALPHA_MIN = %.4e ALPHA_MAX = %.4e H = %.4e C = %.6e N = %d\n", - alpha_min, alpha_max, h, c, n); + //printf("ALPHA_MIN = %.4e ALPHA_MAX = %.4e H = %.4e C = %.6e N = %d\n", + // alpha_min, alpha_max, h, c, n); LindhMalmqvistGagliardiRadialTraits traits(c, h); @@ -479,7 +486,7 @@ TEST_CASE("LMG", "[1d-quad]") { double x = (i+1) * h; double r = c * (std::exp(x) - 1.0); double w = h*(r + c); - printf("%d R = %.16e W = %.16e \n", i, r - quad.points()[i], w-quad.weights()[i]); + //printf("%d R = %.16e W = %.16e \n", i, r - quad.points()[i], w-quad.weights()[i]); } From bdc342ec27da201e30149fe7464574440ffd84cc Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Mon, 6 Nov 2023 15:24:26 -0800 Subject: [PATCH 13/19] Add util_test to Ctest --- test/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 424784b..e0c92f1 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -47,6 +47,7 @@ target_link_libraries( gausschebyshev2 PUBLIC integratorxx_common_ut ) add_executable( composite_quadratures composite_quadratures.cxx ) target_link_libraries( composite_quadratures PUBLIC integratorxx_common_ut ) +add_test( NAME QUADRATURE_UTIL COMMAND util_test ) add_test( NAME QUADRATURE_MANIP COMMAND quadrature_manipulation ) add_test( NAME QUADRATURES_1D COMMAND 1d_quadratures ) add_test( NAME QUADRATURES_COMPOSITE COMMAND composite_quadratures ) From 80e592c3c71ae1f490d889fe9d9cfd7569c854fd Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 3 Jan 2024 11:03:32 -0800 Subject: [PATCH 14/19] Add INTEGRATORXX_HEADER_ONLY + requisite source refactor to enable header-only builds --- CMakeLists.txt | 89 ++++- include/integratorxx/generators/impl/impl.hpp | 5 + .../generators/impl/pruned_grid.hpp | 110 ++++++ .../generators/impl/radial_factory.hpp | 41 +++ .../generators/impl/radial_types.hpp | 12 + .../generators/impl/robust_pruning.hpp | 76 +++++ .../generators/impl/s2_factory.hpp | 40 +++ .../integratorxx/generators/impl/s2_types.hpp | 11 + .../generators/impl/treutler_pruning.hpp | 69 ++++ .../generators/impl/unpruned_grid.hpp | 59 ++++ .../generators/spherical_factory.hpp | 23 +- src/CMakeLists.txt | 76 ----- src/radial_factory.cxx | 41 +-- src/radial_types.hpp | 23 +- src/s2_factory.cxx | 41 +-- src/s2_types.hpp | 12 +- src/spherical_factory.cxx | 318 +----------------- test/CMakeLists.txt | 3 + test/lib_impl.cxx | 2 + 19 files changed, 550 insertions(+), 501 deletions(-) create mode 100644 include/integratorxx/generators/impl/impl.hpp create mode 100644 include/integratorxx/generators/impl/pruned_grid.hpp create mode 100644 include/integratorxx/generators/impl/radial_factory.hpp create mode 100644 include/integratorxx/generators/impl/radial_types.hpp create mode 100644 include/integratorxx/generators/impl/robust_pruning.hpp create mode 100644 include/integratorxx/generators/impl/s2_factory.hpp create mode 100644 include/integratorxx/generators/impl/s2_types.hpp create mode 100644 include/integratorxx/generators/impl/treutler_pruning.hpp create mode 100644 include/integratorxx/generators/impl/unpruned_grid.hpp create mode 100644 test/lib_impl.cxx diff --git a/CMakeLists.txt b/CMakeLists.txt index d3286f9..c56855e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,14 +3,99 @@ cmake_minimum_required( VERSION 3.17 ) # Require CMake 3.17+ # Set up project definition + version information project( IntegratorXX VERSION 1.2.0 LANGUAGES CXX ) -#add_library( integratorxx INTERFACE ) -add_subdirectory(src) +option(INTEGRATORXX_HEADER_ONLY "Force header-only build" OFF) + +if(INTEGRATORXX_HEADER_ONLY) + add_library( integratorxx INTERFACE ) + set(INTEGRATORXX_TARGET_TYPE INTERFACE) +else() + add_subdirectory(src) + set(INTEGRATORXX_TARGET_TYPE PUBLIC) +endif() + +# Target features +target_compile_features( integratorxx ${INTEGRATORXX_TARGET_TYPE} cxx_std_17 ) +target_include_directories( integratorxx + ${INTEGRATORXX_TARGET_TYPE} + $ + $ +) + +include(CheckCXXCompilerFlag) +check_cxx_compiler_flag("-Wno-missing-braces" INTEGRATORXX_HAS_NO_MISSING_BRACES ) +if( INTEGRATORXX_HAS_NO_MISSING_BRACES ) + target_compile_options( integratorxx INTERFACE $<$: -Wno-missing-braces> ) +endif() + + + + +# INSTALL rules +add_library( IntegratorXX::IntegratorXX ALIAS integratorxx ) +include( GNUInstallDirs ) +set( INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/IntegratorXX ) + +# Targets +install(TARGETS integratorxx + EXPORT integratorxx-targets + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} +) + +set_target_properties( integratorxx PROPERTIES EXPORT_NAME IntegratorXX ) + +# Install Headers +install(DIRECTORY ${PROJECT_SOURCE_DIR}/include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} ) + +# Export target to script +install(EXPORT integratorxx-targets + FILE + IntegratorXXTargets.cmake + NAMESPACE + IntegratorXX:: + DESTINATION + ${INSTALL_CONFIGDIR} +) + +# Export build-tree targets also (e.g. to be usable by FetchContent) +export(EXPORT integratorxx-targets + NAMESPACE IntegratorXX:: + FILE "${PROJECT_BINARY_DIR}/IntegratorXXTargets.cmake") + +# Create a ConfigVersion.cmake file +include(CMakePackageConfigHelpers) +write_basic_package_version_file( + ${CMAKE_CURRENT_BINARY_DIR}/IntegratorXXConfigVersion.cmake + VERSION ${PROJECT_VERSION} + COMPATIBILITY AnyNewerVersion +) + +# Setup IntegratorXXConfig.cmake +configure_package_config_file(${PROJECT_SOURCE_DIR}/cmake/IntegratorXXConfig.cmake.in + ${CMAKE_CURRENT_BINARY_DIR}/IntegratorXXConfig.cmake + INSTALL_DESTINATION ${INSTALL_CONFIGDIR} +) + +#Install the config, configversion and custom find modules +install(DIRECTORY + ${PROJECT_SOURCE_DIR}/cmake/ + DESTINATION ${INSTALL_CONFIGDIR} + FILES_MATCHING PATTERN "*.cmake" +) + + +install(FILES + ${CMAKE_CURRENT_BINARY_DIR}/IntegratorXXConfig.cmake + ${CMAKE_CURRENT_BINARY_DIR}/IntegratorXXConfigVersion.cmake + DESTINATION ${INSTALL_CONFIGDIR} +) # Testing if( NOT DEFINED INTEGRATORXX_ENABLE_TESTS ) set( INTEGRATORXX_ENABLE_TESTS ON ) endif() + if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME) include(CTest) endif() diff --git a/include/integratorxx/generators/impl/impl.hpp b/include/integratorxx/generators/impl/impl.hpp new file mode 100644 index 0000000..f8b061d --- /dev/null +++ b/include/integratorxx/generators/impl/impl.hpp @@ -0,0 +1,5 @@ +#include +#include +#include +#include + diff --git a/include/integratorxx/generators/impl/pruned_grid.hpp b/include/integratorxx/generators/impl/pruned_grid.hpp new file mode 100644 index 0000000..c69e036 --- /dev/null +++ b/include/integratorxx/generators/impl/pruned_grid.hpp @@ -0,0 +1,110 @@ +#pragma once + +#include + +#include +#include + +#include +#include + +namespace IntegratorXX { + +namespace detail { + +template +auto make_pruned_grid_impl(const RadialQuadType& rq, + const std::vector& pruning_regions) { + + RadialGridPartition rgp; + for( auto& region : pruning_regions ) { + rgp.add_quad( rq, region.idx_st, AngularQuadType(region.angular_size) ); + } + rgp.finalize(rq); + + return SphericalGridFactory::generate_pruned_grid(rq, std::move(rgp)); + +} + +template +auto make_pruned_grid(const RadialQuadType& rq, + const std::vector& pruning_regions) { + + if(pruning_regions.size() == 0) + throw std::runtime_error("No Pruning Regions"); + + auto angular_quad = pruning_regions[0].angular_quad; + for(auto r : pruning_regions) { + if(r.angular_quad != angular_quad) + throw std::runtime_error("Mixed Angular Pruning Not Supported"); + } + + switch(angular_quad) { + case AngularQuad::AhrensBeylkin: + return make_pruned_grid_impl(rq, pruning_regions); + case AngularQuad::Delley: + return make_pruned_grid_impl(rq, pruning_regions); + case AngularQuad::LebedevLaikov: + return make_pruned_grid_impl(rq, pruning_regions); + case AngularQuad::Womersley: + return make_pruned_grid_impl(rq, pruning_regions); + default: + throw std::runtime_error("Unsupported Angular Quadrature"); + abort(); + } + + +} + +} // Implementation Details + +SphericalGridFactory::spherical_grid_ptr + SphericalGridFactory::generate_pruned_grid( RadialQuad rq, + const RadialTraits& traits, + const std::vector& pruning_regions) { + + switch( rq ) { + + case RadialQuad::Becke: + return detail::make_pruned_grid( bk_type(traits), pruning_regions ); + case RadialQuad::MuraKnowles: + return detail::make_pruned_grid( mk_type(traits), pruning_regions ); + case RadialQuad::MurrayHandyLaming: + return detail::make_pruned_grid( mhl_type(traits), pruning_regions ); + case RadialQuad::TreutlerAhlrichs: + return detail::make_pruned_grid( ta_type(traits), pruning_regions ); + + default: + throw std::runtime_error("Unsupported Radial Quadrature"); + abort(); + + } + +} + + +PrunedSphericalGridSpecification create_pruned_spec( + PruningScheme scheme, UnprunedSphericalGridSpecification unp +) { + + if(!unp.radial_traits) throw std::runtime_error("RadialTraits Not Set"); + switch(scheme) { + case PruningScheme::Robust: + return robust_psi4_pruning_scheme(unp); + case PruningScheme::Treutler: + return treutler_pruning_scheme(unp); + + // Default to Unpruned Grid + case PruningScheme::Unpruned: + default: + std::vector pruning_regions = { + {0ul, unp.radial_traits->npts(), unp.angular_quad, unp.angular_size} + }; + return PrunedSphericalGridSpecification( + unp.radial_quad, unp.radial_traits->clone(), pruning_regions + ); + } + +} + +} diff --git a/include/integratorxx/generators/impl/radial_factory.hpp b/include/integratorxx/generators/impl/radial_factory.hpp new file mode 100644 index 0000000..ed89d56 --- /dev/null +++ b/include/integratorxx/generators/impl/radial_factory.hpp @@ -0,0 +1,41 @@ +#pragma once +#include + +#include + +#include + +namespace IntegratorXX { + +RadialQuad radial_from_string(std::string name) { + std::transform(name.begin(), name.end(), name.begin(), ::toupper); + if(name == "BECKE") return RadialQuad::Becke; + if(name == "MURAKNOWLES") return RadialQuad::MuraKnowles; + if(name == "MK") return RadialQuad::MuraKnowles; + if(name == "MURRAYHANDYLAMING") return RadialQuad::MurrayHandyLaming; + if(name == "MHL") return RadialQuad::MurrayHandyLaming; + if(name == "TREUTLERAHLRICHS") return RadialQuad::TreutlerAhlrichs; + if(name == "TA") return RadialQuad::TreutlerAhlrichs; + + throw std::runtime_error("Unrecognized Radial Quadrature"); +} + +RadialFactory::radial_grid_ptr RadialFactory::generate(RadialQuad rq, const RadialTraits& traits) { + + switch(rq) { + case RadialQuad::Becke: + return std::make_unique(traits); + case RadialQuad::MuraKnowles: + return std::make_unique(traits); + case RadialQuad::MurrayHandyLaming: + return std::make_unique(traits); + case RadialQuad::TreutlerAhlrichs: + return std::make_unique(traits); + default: + throw std::runtime_error("Unsupported Radial Quadrature"); + abort(); + } +} + + +} diff --git a/include/integratorxx/generators/impl/radial_types.hpp b/include/integratorxx/generators/impl/radial_types.hpp new file mode 100644 index 0000000..4edb1da --- /dev/null +++ b/include/integratorxx/generators/impl/radial_types.hpp @@ -0,0 +1,12 @@ +#pragma once +#include + +namespace IntegratorXX { + +using bk_type = Becke; +using mk_type = MuraKnowles; +using mhl_type = MurrayHandyLaming; +using ta_type = TreutlerAhlrichs; + +} + diff --git a/include/integratorxx/generators/impl/robust_pruning.hpp b/include/integratorxx/generators/impl/robust_pruning.hpp new file mode 100644 index 0000000..f397fea --- /dev/null +++ b/include/integratorxx/generators/impl/robust_pruning.hpp @@ -0,0 +1,76 @@ +#pragma once +#include + +#include + +namespace IntegratorXX { + +namespace detail { + +template +auto get_robust_low_med_sizes(AngularSize asz) { + using traits = quadrature_traits; + const auto base_order = traits::algebraic_order_by_npts(asz); + if( base_order < 0 ) throw std::runtime_error("Invalid Base Grid"); + + const auto med_order = + traits::next_algebraic_order(base_order > 6 ? base_order-6 : base_order); + const auto low_order = traits::next_algebraic_order(7); + + AngularSize med_sz(traits::npts_by_algebraic_order(med_order)); + AngularSize low_sz(traits::npts_by_algebraic_order(low_order)); + + return std::make_pair(low_sz, med_sz); +} + + +PrunedSphericalGridSpecification robust_psi4_pruning_scheme_impl( + size_t low_sz, size_t med_sz, AngularQuad angular_quad, + UnprunedSphericalGridSpecification unp ) { + + const size_t rsz = unp.radial_traits->npts(); + const size_t r_div_4 = rsz / 4ul + 1ul; + const size_t r_div_2 = rsz / 2ul + 1ul; + std::vector pruning_regions = { + {0ul, r_div_4, angular_quad, low_sz}, + {r_div_4, r_div_2, angular_quad, med_sz}, + {r_div_2, rsz, angular_quad, unp.angular_size} + }; + + return PrunedSphericalGridSpecification( + unp.radial_quad, unp.radial_traits->clone(), pruning_regions + ); + +} + +} // Implementation Details + + +PrunedSphericalGridSpecification robust_psi4_pruning_scheme( + UnprunedSphericalGridSpecification unp ) { + + size_t low_sz, med_sz; + const auto angular_quad = unp.angular_quad; + const auto asz = unp.angular_size; + switch(angular_quad) { + case AngularQuad::AhrensBeylkin: + std::tie(low_sz, med_sz) = detail::get_robust_low_med_sizes(asz); + break; + case AngularQuad::Delley: + std::tie(low_sz, med_sz) = detail::get_robust_low_med_sizes(asz); + break; + case AngularQuad::LebedevLaikov: + std::tie(low_sz, med_sz) = detail::get_robust_low_med_sizes(asz); + break; + case AngularQuad::Womersley: + std::tie(low_sz, med_sz) = detail::get_robust_low_med_sizes(asz); + break; + default: + throw std::runtime_error("Unsupported Angular Quadrature"); + abort(); + } + + return detail::robust_psi4_pruning_scheme_impl(low_sz, med_sz, angular_quad, unp); +} + +} diff --git a/include/integratorxx/generators/impl/s2_factory.hpp b/include/integratorxx/generators/impl/s2_factory.hpp new file mode 100644 index 0000000..6dafb3d --- /dev/null +++ b/include/integratorxx/generators/impl/s2_factory.hpp @@ -0,0 +1,40 @@ +#pragma once +#include + +#include + +#include + +namespace IntegratorXX { + +AngularQuad angular_from_string(std::string name) { + std::transform(name.begin(), name.end(), name.begin(), ::toupper); + if(name == "AHRENSBEYLKIN") return AngularQuad::AhrensBeylkin; + if(name == "AB") return AngularQuad::AhrensBeylkin; + if(name == "DELLEY") return AngularQuad::Delley; + if(name == "LEBEDEVLAIKOV") return AngularQuad::LebedevLaikov; + if(name == "LEBEDEV") return AngularQuad::LebedevLaikov; + if(name == "LL") return AngularQuad::LebedevLaikov; + if(name == "WOMERSLEY") return AngularQuad::Womersley; + + throw std::runtime_error("Unrecognized Angular Quadrature"); +} + +S2Factory::s2_grid_ptr S2Factory::generate(AngularQuad aq, size_t npts) { + + switch(aq) { + case AngularQuad::AhrensBeylkin: + return std::make_unique(npts); + case AngularQuad::Delley: + return std::make_unique(npts); + case AngularQuad::LebedevLaikov: + return std::make_unique(npts); + case AngularQuad::Womersley: + return std::make_unique(npts); + default: + throw std::runtime_error("Unsupported Angular Quadrature"); + abort(); + } +} + +} diff --git a/include/integratorxx/generators/impl/s2_types.hpp b/include/integratorxx/generators/impl/s2_types.hpp new file mode 100644 index 0000000..9433881 --- /dev/null +++ b/include/integratorxx/generators/impl/s2_types.hpp @@ -0,0 +1,11 @@ +#pragma once +#include + +namespace IntegratorXX { + +using ah_type = AhrensBeylkin; +using de_type = Delley; +using ll_type = LebedevLaikov; +using wo_type = Womersley; + +} diff --git a/include/integratorxx/generators/impl/treutler_pruning.hpp b/include/integratorxx/generators/impl/treutler_pruning.hpp new file mode 100644 index 0000000..d59a107 --- /dev/null +++ b/include/integratorxx/generators/impl/treutler_pruning.hpp @@ -0,0 +1,69 @@ +#pragma once +#include + +#include + +namespace IntegratorXX { + +namespace detail { + +template +auto get_treutler_low_med_sizes() { + using traits = quadrature_traits; + AngularSize med_sz(traits::npts_by_algebraic_order(traits::next_algebraic_order(11))); + AngularSize low_sz(traits::npts_by_algebraic_order(traits::next_algebraic_order(7 ))); + + return std::make_pair(low_sz, med_sz); +} + +PrunedSphericalGridSpecification treutler_pruning_scheme_impl( + size_t low_sz, size_t med_sz, AngularQuad angular_quad, + UnprunedSphericalGridSpecification unp ) { + + const size_t rsz = unp.radial_traits->npts(); + const size_t r_div_3 = rsz / 3ul + 1ul; + const size_t r_div_2 = rsz / 2ul + 1ul; + std::vector pruning_regions = { + {0ul, r_div_3, angular_quad, low_sz}, + {r_div_3, r_div_2, angular_quad, med_sz}, + {r_div_2, rsz, angular_quad, unp.angular_size} + }; + + return PrunedSphericalGridSpecification( + unp.radial_quad, unp.radial_traits->clone(), pruning_regions + ); + +} + +} // Implementation details + + + +PrunedSphericalGridSpecification treutler_pruning_scheme( + UnprunedSphericalGridSpecification unp ) { + + size_t low_sz, med_sz; + const auto angular_quad = unp.angular_quad; + switch(angular_quad) { + case AngularQuad::AhrensBeylkin: + std::tie(low_sz, med_sz) = detail::get_treutler_low_med_sizes(); + break; + case AngularQuad::Delley: + std::tie(low_sz, med_sz) = detail::get_treutler_low_med_sizes(); + break; + case AngularQuad::LebedevLaikov: + std::tie(low_sz, med_sz) = detail::get_treutler_low_med_sizes(); + break; + case AngularQuad::Womersley: + std::tie(low_sz, med_sz) = detail::get_treutler_low_med_sizes(); + break; + default: + throw std::runtime_error("Unsupported Angular Quadrature"); + abort(); + } + + return detail::treutler_pruning_scheme_impl(low_sz, med_sz, angular_quad, unp); +} + + +} diff --git a/include/integratorxx/generators/impl/unpruned_grid.hpp b/include/integratorxx/generators/impl/unpruned_grid.hpp new file mode 100644 index 0000000..758b1fd --- /dev/null +++ b/include/integratorxx/generators/impl/unpruned_grid.hpp @@ -0,0 +1,59 @@ +#pragma once + +#include + +#include +#include + +namespace IntegratorXX { + +namespace detail { + +template +auto generate_unpruned_grid_impl(RadialQuad rq, const RadialTraits& traits, + AngularQuadType&& ang_quad) { + + switch( rq ) { + + case RadialQuad::Becke: + return SphericalGridFactory::generate_unpruned_grid( bk_type(traits), std::forward(ang_quad) ); + + case RadialQuad::MuraKnowles: + return SphericalGridFactory::generate_unpruned_grid( mk_type(traits), std::forward(ang_quad) ); + + case RadialQuad::MurrayHandyLaming: + return SphericalGridFactory::generate_unpruned_grid( mhl_type(traits), std::forward(ang_quad) ); + + case RadialQuad::TreutlerAhlrichs: + return SphericalGridFactory::generate_unpruned_grid( ta_type(traits), std::forward(ang_quad) ); + + default: + throw std::runtime_error("Unsupported Radial Quadrature"); + abort(); + + } + +} + +} // Implementation details + +SphericalGridFactory::spherical_grid_ptr + SphericalGridFactory::generate_unpruned_grid( RadialQuad rq, + const RadialTraits& traits, AngularQuad aq, AngularSize nang) { + + switch(aq) { + case AngularQuad::AhrensBeylkin: + return detail::generate_unpruned_grid_impl(rq, traits, ah_type(nang)); + case AngularQuad::Delley: + return detail::generate_unpruned_grid_impl(rq, traits, de_type(nang)); + case AngularQuad::LebedevLaikov: + return detail::generate_unpruned_grid_impl(rq, traits, ll_type(nang)); + case AngularQuad::Womersley: + return detail::generate_unpruned_grid_impl(rq, traits, wo_type(nang)); + default: + throw std::runtime_error("Unsupported Angular Quadrature"); + abort(); + } +} + +} diff --git a/include/integratorxx/generators/spherical_factory.hpp b/include/integratorxx/generators/spherical_factory.hpp index 91f391d..43e2844 100644 --- a/include/integratorxx/generators/spherical_factory.hpp +++ b/include/integratorxx/generators/spherical_factory.hpp @@ -30,12 +30,13 @@ struct UnprunedSphericalGridSpecification { AngularQuad angular_quad; /// Angular quadrature specification AngularSize angular_size; /// Number of angular quadrature points - UnprunedSphericalGridSpecification(RadialQuad, const RadialTraits&, AngularQuad, - AngularSize); + UnprunedSphericalGridSpecification( RadialQuad rq, const RadialTraits& traits, + AngularQuad aq, AngularSize as) : radial_quad(rq), radial_traits(traits.clone()), + angular_quad(aq), angular_size(as) {} UnprunedSphericalGridSpecification(const UnprunedSphericalGridSpecification& other) : radial_quad(other.radial_quad), radial_traits(other.radial_traits ? other.radial_traits->clone() : nullptr), - angular_quad(other.angular_quad), angular_size(other.angular_size) {}; + angular_quad(other.angular_quad), angular_size(other.angular_size) {} }; @@ -168,8 +169,20 @@ struct SphericalGridFactory { const std::vector&); - static spherical_grid_ptr generate_grid(UnprunedSphericalGridSpecification gs); - static spherical_grid_ptr generate_grid(PrunedSphericalGridSpecification gs); + static inline spherical_grid_ptr generate_grid( + UnprunedSphericalGridSpecification gs) { + if(!gs.radial_traits) throw std::runtime_error("RadialTraits Not Set"); + return generate_unpruned_grid(gs.radial_quad, *gs.radial_traits, + gs.angular_quad, gs.angular_size); + } + + + static inline spherical_grid_ptr generate_grid( + PrunedSphericalGridSpecification gs) { + if(!gs.radial_traits) throw std::runtime_error("RadialTraits Not Set"); + return generate_pruned_grid( gs.radial_quad, *gs.radial_traits, + gs.pruning_regions ); + } }; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ea07e5c..75733d3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -4,80 +4,4 @@ add_library(integratorxx s2_factory.cxx ) -# Target features -target_compile_features( integratorxx PUBLIC cxx_std_17 ) -target_include_directories( integratorxx - PUBLIC - $ - $ -) - -include(CheckCXXCompilerFlag) -check_cxx_compiler_flag("-Wno-missing-braces" INTEGRATORXX_HAS_NO_MISSING_BRACES ) -if( INTEGRATORXX_HAS_NO_MISSING_BRACES ) - target_compile_options( integratorxx INTERFACE $<$: -Wno-missing-braces> ) -endif() - - - - -# INSTALL rules -add_library( IntegratorXX::IntegratorXX ALIAS integratorxx ) -include( GNUInstallDirs ) -set( INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/IntegratorXX ) - -# Targets -install(TARGETS integratorxx - EXPORT integratorxx-targets - LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} - ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} -) - -set_target_properties( integratorxx PROPERTIES EXPORT_NAME IntegratorXX ) - -# Install Headers -install(DIRECTORY ${PROJECT_SOURCE_DIR}/include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} ) - -# Export target to script -install(EXPORT integratorxx-targets - FILE - IntegratorXXTargets.cmake - NAMESPACE - IntegratorXX:: - DESTINATION - ${INSTALL_CONFIGDIR} -) - -# Export build-tree targets also (e.g. to be usable by FetchContent) -export(EXPORT integratorxx-targets - NAMESPACE IntegratorXX:: - FILE "${PROJECT_BINARY_DIR}/IntegratorXXTargets.cmake") - -# Create a ConfigVersion.cmake file -include(CMakePackageConfigHelpers) -write_basic_package_version_file( - ${CMAKE_CURRENT_BINARY_DIR}/IntegratorXXConfigVersion.cmake - VERSION ${PROJECT_VERSION} - COMPATIBILITY AnyNewerVersion -) - -# Setup IntegratorXXConfig.cmake -configure_package_config_file(${PROJECT_SOURCE_DIR}/cmake/IntegratorXXConfig.cmake.in - ${CMAKE_CURRENT_BINARY_DIR}/IntegratorXXConfig.cmake - INSTALL_DESTINATION ${INSTALL_CONFIGDIR} -) - -#Install the config, configversion and custom find modules -install(DIRECTORY - ${PROJECT_SOURCE_DIR}/cmake/ - DESTINATION ${INSTALL_CONFIGDIR} - FILES_MATCHING PATTERN "*.cmake" -) - - -install(FILES - ${CMAKE_CURRENT_BINARY_DIR}/IntegratorXXConfig.cmake - ${CMAKE_CURRENT_BINARY_DIR}/IntegratorXXConfigVersion.cmake - DESTINATION ${INSTALL_CONFIGDIR} -) diff --git a/src/radial_factory.cxx b/src/radial_factory.cxx index aa0352b..fa8ddb1 100644 --- a/src/radial_factory.cxx +++ b/src/radial_factory.cxx @@ -1,41 +1,2 @@ -#include -#include "radial_types.hpp" +#include -#include - -namespace IntegratorXX { - -RadialQuad radial_from_string(std::string name) { - std::transform(name.begin(), name.end(), name.begin(), ::toupper); - if(name == "BECKE") return RadialQuad::Becke; - if(name == "MURAKNOWLES") return RadialQuad::MuraKnowles; - if(name == "MK") return RadialQuad::MuraKnowles; - if(name == "MURRAYHANDYLAMING") return RadialQuad::MurrayHandyLaming; - if(name == "MHL") return RadialQuad::MurrayHandyLaming; - if(name == "TREUTLERAHLRICHS") return RadialQuad::TreutlerAhlrichs; - if(name == "TA") return RadialQuad::TreutlerAhlrichs; - - throw std::runtime_error("Unrecognized Radial Quadrature"); -} - -using radial_grid_ptr = RadialFactory::radial_grid_ptr; - -radial_grid_ptr RadialFactory::generate(RadialQuad rq, const RadialTraits& traits) { - - switch(rq) { - case RadialQuad::Becke: - return std::make_unique(traits); - case RadialQuad::MuraKnowles: - return std::make_unique(traits); - case RadialQuad::MurrayHandyLaming: - return std::make_unique(traits); - case RadialQuad::TreutlerAhlrichs: - return std::make_unique(traits); - default: - throw std::runtime_error("Unsupported Radial Quadrature"); - abort(); - } -} - - -} diff --git a/src/radial_types.hpp b/src/radial_types.hpp index 4edb1da..8229a03 100644 --- a/src/radial_types.hpp +++ b/src/radial_types.hpp @@ -1,12 +1,13 @@ -#pragma once -#include - -namespace IntegratorXX { - -using bk_type = Becke; -using mk_type = MuraKnowles; -using mhl_type = MurrayHandyLaming; -using ta_type = TreutlerAhlrichs; - -} +#include +//#pragma once +//#include +// +//namespace IntegratorXX { +// +//using bk_type = Becke; +//using mk_type = MuraKnowles; +//using mhl_type = MurrayHandyLaming; +//using ta_type = TreutlerAhlrichs; +// +//} diff --git a/src/s2_factory.cxx b/src/s2_factory.cxx index e3a1d93..ad9844e 100644 --- a/src/s2_factory.cxx +++ b/src/s2_factory.cxx @@ -1,41 +1,2 @@ -#include +#include -#include "s2_types.hpp" - -#include - -namespace IntegratorXX { - -AngularQuad angular_from_string(std::string name) { - std::transform(name.begin(), name.end(), name.begin(), ::toupper); - if(name == "AHRENSBEYLKIN") return AngularQuad::AhrensBeylkin; - if(name == "AB") return AngularQuad::AhrensBeylkin; - if(name == "DELLEY") return AngularQuad::Delley; - if(name == "LEBEDEVLAIKOV") return AngularQuad::LebedevLaikov; - if(name == "LEBEDEV") return AngularQuad::LebedevLaikov; - if(name == "LL") return AngularQuad::LebedevLaikov; - if(name == "WOMERSLEY") return AngularQuad::Womersley; - - throw std::runtime_error("Unrecognized Angular Quadrature"); -} - -using s2_grid_ptr = S2Factory::s2_grid_ptr; - -s2_grid_ptr S2Factory::generate(AngularQuad aq, size_t npts) { - - switch(aq) { - case AngularQuad::AhrensBeylkin: - return std::make_unique(npts); - case AngularQuad::Delley: - return std::make_unique(npts); - case AngularQuad::LebedevLaikov: - return std::make_unique(npts); - case AngularQuad::Womersley: - return std::make_unique(npts); - default: - throw std::runtime_error("Unsupported Angular Quadrature"); - abort(); - } -} - -} diff --git a/src/s2_types.hpp b/src/s2_types.hpp index 801d617..4db0451 100644 --- a/src/s2_types.hpp +++ b/src/s2_types.hpp @@ -1,11 +1 @@ -#include - - -namespace IntegratorXX { - -using ah_type = AhrensBeylkin; -using de_type = Delley; -using ll_type = LebedevLaikov; -using wo_type = Womersley; - -} +#include diff --git a/src/spherical_factory.cxx b/src/spherical_factory.cxx index 513f0e9..fbacc41 100644 --- a/src/spherical_factory.cxx +++ b/src/spherical_factory.cxx @@ -1,317 +1,3 @@ -#include +#include +#include -#include "radial_types.hpp" -#include "s2_types.hpp" - -namespace IntegratorXX { - - -UnprunedSphericalGridSpecification::UnprunedSphericalGridSpecification( - RadialQuad rq, const RadialTraits& traits, AngularQuad aq, AngularSize as) : - radial_quad(rq), radial_traits(traits.clone()), angular_quad(aq), angular_size(as) { } - - - - -using spherical_grid_ptr = SphericalGridFactory::spherical_grid_ptr; - - - -/************************/ -/**** Unpruned Grids ****/ -/************************/ -template -spherical_grid_ptr generate_unpruned_grid_impl(RadialQuad rq, const RadialTraits& traits, - AngularQuadType&& ang_quad) { - - switch( rq ) { - - case RadialQuad::Becke: - return SphericalGridFactory::generate_unpruned_grid( bk_type(traits), std::forward(ang_quad) ); - - case RadialQuad::MuraKnowles: - return SphericalGridFactory::generate_unpruned_grid( mk_type(traits), std::forward(ang_quad) ); - - case RadialQuad::MurrayHandyLaming: - return SphericalGridFactory::generate_unpruned_grid( mhl_type(traits), std::forward(ang_quad) ); - - case RadialQuad::TreutlerAhlrichs: - return SphericalGridFactory::generate_unpruned_grid( ta_type(traits), std::forward(ang_quad) ); - - default: - throw std::runtime_error("Unsupported Radial Quadrature"); - abort(); - - } - - -} - -spherical_grid_ptr SphericalGridFactory::generate_unpruned_grid( RadialQuad rq, - const RadialTraits& traits, AngularQuad aq, AngularSize nang) { - - switch(aq) { - case AngularQuad::AhrensBeylkin: - return generate_unpruned_grid_impl(rq, traits, ah_type(nang)); - case AngularQuad::Delley: - return generate_unpruned_grid_impl(rq, traits, de_type(nang)); - case AngularQuad::LebedevLaikov: - return generate_unpruned_grid_impl(rq, traits, ll_type(nang)); - case AngularQuad::Womersley: - return generate_unpruned_grid_impl(rq, traits, wo_type(nang)); - default: - throw std::runtime_error("Unsupported Angular Quadrature"); - abort(); - } -} - -spherical_grid_ptr SphericalGridFactory::generate_grid( - UnprunedSphericalGridSpecification gs ) { - if(!gs.radial_traits) throw std::runtime_error("RadialTraits Not Set"); - return generate_unpruned_grid(gs.radial_quad, *gs.radial_traits, - gs.angular_quad, gs.angular_size); -} - - - - -/**********************/ -/**** Pruned Grids ****/ -/**********************/ - -template -auto make_pruned_grid_impl(const RadialQuadType& rq, - const std::vector& pruning_regions) { - - RadialGridPartition rgp; - for( auto& region : pruning_regions ) { - rgp.add_quad( rq, region.idx_st, AngularQuadType(region.angular_size) ); - } - rgp.finalize(rq); - - return SphericalGridFactory::generate_pruned_grid(rq, std::move(rgp)); - -} - -template -auto make_pruned_grid(const RadialQuadType& rq, - const std::vector& pruning_regions) { - - if(pruning_regions.size() == 0) - throw std::runtime_error("No Pruning Regions"); - - auto angular_quad = pruning_regions[0].angular_quad; - for(auto r : pruning_regions) { - if(r.angular_quad != angular_quad) - throw std::runtime_error("Mixed Angular Pruning Not Supported"); - } - - switch(angular_quad) { - case AngularQuad::AhrensBeylkin: - return make_pruned_grid_impl(rq, pruning_regions); - case AngularQuad::Delley: - return make_pruned_grid_impl(rq, pruning_regions); - case AngularQuad::LebedevLaikov: - return make_pruned_grid_impl(rq, pruning_regions); - case AngularQuad::Womersley: - return make_pruned_grid_impl(rq, pruning_regions); - default: - throw std::runtime_error("Unsupported Angular Quadrature"); - abort(); - } - - -} - -spherical_grid_ptr SphericalGridFactory::generate_pruned_grid( - RadialQuad rq, const RadialTraits& traits, - const std::vector& pruning_regions) { - - switch( rq ) { - - case RadialQuad::Becke: - return make_pruned_grid( bk_type(traits), pruning_regions ); - case RadialQuad::MuraKnowles: - return make_pruned_grid( mk_type(traits), pruning_regions ); - case RadialQuad::MurrayHandyLaming: - return make_pruned_grid( mhl_type(traits), pruning_regions ); - case RadialQuad::TreutlerAhlrichs: - return make_pruned_grid( ta_type(traits), pruning_regions ); - - default: - throw std::runtime_error("Unsupported Radial Quadrature"); - abort(); - - } - -} - -spherical_grid_ptr SphericalGridFactory::generate_grid( - PrunedSphericalGridSpecification gs ) { - if(!gs.radial_traits) throw std::runtime_error("RadialTraits Not Set"); - return generate_pruned_grid( gs.radial_quad, *gs.radial_traits, - gs.pruning_regions ); -} - - - - -/*************************** - * Default Pruning Schemes * - ***************************/ - - -/*** Psi4 Robust ***/ - - -template -auto get_robust_low_med_sizes(AngularSize asz) { - using traits = quadrature_traits; - const auto base_order = traits::algebraic_order_by_npts(asz); - if( base_order < 0 ) throw std::runtime_error("Invalid Base Grid"); - - const auto med_order = - traits::next_algebraic_order(base_order > 6 ? base_order-6 : base_order); - const auto low_order = traits::next_algebraic_order(7); - - AngularSize med_sz(traits::npts_by_algebraic_order(med_order)); - AngularSize low_sz(traits::npts_by_algebraic_order(low_order)); - - return std::make_pair(low_sz, med_sz); -} - - -PrunedSphericalGridSpecification robust_psi4_pruning_scheme_impl( - size_t low_sz, size_t med_sz, AngularQuad angular_quad, - UnprunedSphericalGridSpecification unp ) { - - const size_t rsz = unp.radial_traits->npts(); - const size_t r_div_4 = rsz / 4ul + 1ul; - const size_t r_div_2 = rsz / 2ul + 1ul; - std::vector pruning_regions = { - {0ul, r_div_4, angular_quad, low_sz}, - {r_div_4, r_div_2, angular_quad, med_sz}, - {r_div_2, rsz, angular_quad, unp.angular_size} - }; - - return PrunedSphericalGridSpecification( - unp.radial_quad, unp.radial_traits->clone(), pruning_regions - ); - -} - -PrunedSphericalGridSpecification robust_psi4_pruning_scheme( - UnprunedSphericalGridSpecification unp ) { - - size_t low_sz, med_sz; - const auto angular_quad = unp.angular_quad; - const auto asz = unp.angular_size; - switch(angular_quad) { - case AngularQuad::AhrensBeylkin: - std::tie(low_sz, med_sz) = get_robust_low_med_sizes(asz); - break; - case AngularQuad::Delley: - std::tie(low_sz, med_sz) = get_robust_low_med_sizes(asz); - break; - case AngularQuad::LebedevLaikov: - std::tie(low_sz, med_sz) = get_robust_low_med_sizes(asz); - break; - case AngularQuad::Womersley: - std::tie(low_sz, med_sz) = get_robust_low_med_sizes(asz); - break; - default: - throw std::runtime_error("Unsupported Angular Quadrature"); - abort(); - } - - return robust_psi4_pruning_scheme_impl(low_sz, med_sz, angular_quad, unp); -} - - - -/*** Treutler-Ahlrichs Pruning ***/ - - -template -auto get_treutler_low_med_sizes() { - using traits = quadrature_traits; - AngularSize med_sz(traits::npts_by_algebraic_order(traits::next_algebraic_order(11))); - AngularSize low_sz(traits::npts_by_algebraic_order(traits::next_algebraic_order(7 ))); - - return std::make_pair(low_sz, med_sz); -} - -PrunedSphericalGridSpecification treutler_pruning_scheme_impl( - size_t low_sz, size_t med_sz, AngularQuad angular_quad, - UnprunedSphericalGridSpecification unp ) { - - const size_t rsz = unp.radial_traits->npts(); - const size_t r_div_3 = rsz / 3ul + 1ul; - const size_t r_div_2 = rsz / 2ul + 1ul; - std::vector pruning_regions = { - {0ul, r_div_3, angular_quad, low_sz}, - {r_div_3, r_div_2, angular_quad, med_sz}, - {r_div_2, rsz, angular_quad, unp.angular_size} - }; - - return PrunedSphericalGridSpecification( - unp.radial_quad, unp.radial_traits->clone(), pruning_regions - ); - -} - - - -PrunedSphericalGridSpecification treutler_pruning_scheme( - UnprunedSphericalGridSpecification unp ) { - - size_t low_sz, med_sz; - const auto angular_quad = unp.angular_quad; - switch(angular_quad) { - case AngularQuad::AhrensBeylkin: - std::tie(low_sz, med_sz) = get_treutler_low_med_sizes(); - break; - case AngularQuad::Delley: - std::tie(low_sz, med_sz) = get_treutler_low_med_sizes(); - break; - case AngularQuad::LebedevLaikov: - std::tie(low_sz, med_sz) = get_treutler_low_med_sizes(); - break; - case AngularQuad::Womersley: - std::tie(low_sz, med_sz) = get_treutler_low_med_sizes(); - break; - default: - throw std::runtime_error("Unsupported Angular Quadrature"); - abort(); - } - - return treutler_pruning_scheme_impl(low_sz, med_sz, angular_quad, unp); -} - - - -PrunedSphericalGridSpecification create_pruned_spec( - PruningScheme scheme, UnprunedSphericalGridSpecification unp -) { - - if(!unp.radial_traits) throw std::runtime_error("RadialTraits Not Set"); - switch(scheme) { - case PruningScheme::Robust: - return robust_psi4_pruning_scheme(unp); - case PruningScheme::Treutler: - return treutler_pruning_scheme(unp); - - // Default to Unpruned Grid - case PruningScheme::Unpruned: - default: - std::vector pruning_regions = { - {0ul, unp.radial_traits->npts(), unp.angular_quad, unp.angular_size} - }; - return PrunedSphericalGridSpecification( - unp.radial_quad, unp.radial_traits->clone(), pruning_regions - ); - } - -} - -} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7adc743..05bd933 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -19,6 +19,9 @@ else() endif() add_library(integratorxx_common_ut quad_matcher.cxx) +if(INTEGRATORXX_HEADER_ONLY) + target_sources(integratorxx_common_ut PRIVATE lib_impl.cxx) +endif() target_link_libraries( integratorxx_common_ut PUBLIC Catch2::Catch2WithMain integratorxx ) add_executable( quadrature_manipulation quadrature_manipulation.cxx ) diff --git a/test/lib_impl.cxx b/test/lib_impl.cxx new file mode 100644 index 0000000..ba83f11 --- /dev/null +++ b/test/lib_impl.cxx @@ -0,0 +1,2 @@ +#include + From 482d8717e9204fd2576b777dcdeaf041a1659954 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 3 Jan 2024 11:15:11 -0800 Subject: [PATCH 15/19] Update README with header-only instructions --- README.md | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/README.md b/README.md index 4f624c4..d416142 100644 --- a/README.md +++ b/README.md @@ -142,6 +142,26 @@ const auto& weights = sph_quad->weights(); // std::vector ``` +### Header-only builds + +With the exception of the runtime grid generator, the entirety +of the grid specification in IntegratorXX is header-only and can operate +without pre-compiled components. By default, the runtime generator is +pre-compiled to improve the efficiency of the compilation process and to +avoid excessive build times in complex projects with aggressive compiler +optimization. **N.B. it is highly recommend that users maintain this default +behavior to avoid excessive compilation sizes and build times**. + +IntegratorXX also allows for header-only use of the runtime generator by +setting `INTEGRATORXX_HEADER_ONLY=ON`. +This feature also allows for circumvention of +the CMake build system by simply including the requisite implementation +header. + +To use the runtime generator header-only, one needs to include +`` **exactly once** per project, +otherwise duplicate / incompatible symbols will occur. + ## Contributing and Bug Reports We welcome any and all contributions and encourage bug reports. Please use the From e92cb07e95df62eb2d557f84cf1db08d76c9965f Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 3 Jan 2024 11:19:29 -0800 Subject: [PATCH 16/19] Add header-only to CI --- .github/workflows/build_and_test_compiler_zoo.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/build_and_test_compiler_zoo.yml b/.github/workflows/build_and_test_compiler_zoo.yml index 4192e5a..be32d71 100644 --- a/.github/workflows/build_and_test_compiler_zoo.yml +++ b/.github/workflows/build_and_test_compiler_zoo.yml @@ -49,6 +49,9 @@ jobs: runs-on: ubuntu-latest container: image: dbwy/chemistry + strategy: + matrix: + header_only: [ON,OFF] steps: - uses: actions/checkout@v3 @@ -64,6 +67,7 @@ jobs: - name: Configure CMake shell: bash run: cmake -S $GITHUB_WORKSPACE -B ${{runner.workspace}}/build + -DINTEGRATORXX_HEADER_ONLY=${{matrix.header_only}} -DCMAKE_INSTALL_PREFIX=${{runner.workspace}}/install -DCMAKE_PREFIX_PATH=${ENV_PREFIX_PATH} -DCMAKE_TOOLCHAIN_FILE=${GITHUB_WORKSPACE}/${GH_ACTIONS_TOOLCHAIN} From bd95ee04531f8230ffc511d1aea8149e26e3745a Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 3 Jan 2024 12:59:05 -0800 Subject: [PATCH 17/19] Testing LMG grids --- test/1d_quadratures.cxx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 2a1b120..699987a 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -261,8 +261,8 @@ TEST_CASE("LMG", "[1d-quad]") { using namespace Catch::Matchers; SECTION("HALF INTEGRAL GAMMA") { - REQUIRE_THAT(half_integral_tgamma(1), WithinAbs(std::sqrt(M_PI), 1e-15)); - REQUIRE_THAT(half_integral_tgamma(3), WithinAbs(std::sqrt(M_PI)/2,1e-15)); + REQUIRE_THAT(half_integer_tgamma(1), WithinAbs(std::sqrt(M_PI), 1e-15)); + REQUIRE_THAT(half_integer_tgamma(3), WithinAbs(std::sqrt(M_PI)/2,1e-15)); } SECTION("Lambert W") { @@ -488,6 +488,8 @@ TEST_CASE("LMG", "[1d-quad]") { double r = c * (std::exp(x) - 1.0); double w = h*(r + c); //printf("%d R = %.16e W = %.16e \n", i, r - quad.points()[i], w-quad.weights()[i]); + REQUIRE_THAT(quad.points()[i], Catch::Matchers::WithinRel(r)); + REQUIRE_THAT(quad.weights()[i], Catch::Matchers::WithinRel(w)); } From dba0abc302bd2ee014b0d44cd8ef6a58ba854ced Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 3 Jan 2024 13:38:53 -0800 Subject: [PATCH 18/19] Make LMG work with new RadialTraits schema --- .../integratorxx/quadratures/radial/lmg.hpp | 33 ++++++++++++++----- test/1d_quadratures.cxx | 6 ++-- 2 files changed, 27 insertions(+), 12 deletions(-) diff --git a/include/integratorxx/quadratures/radial/lmg.hpp b/include/integratorxx/quadratures/radial/lmg.hpp index 75b4a4b..f353092 100644 --- a/include/integratorxx/quadratures/radial/lmg.hpp +++ b/include/integratorxx/quadratures/radial/lmg.hpp @@ -73,19 +73,35 @@ inline double step_size(int m, double prec) { } -class LindhMalmqvistGagliardiRadialTraits { +class LindhMalmqvistGagliardiRadialTraits : public RadialTraits { + using self_type = LindhMalmqvistGagliardiRadialTraits; - double step_size_; + size_t npts_; double c_; + double step_size_; public: - LindhMalmqvistGagliardiRadialTraits(double c, double step_size) : - step_size_(step_size), c_(c) { } + LindhMalmqvistGagliardiRadialTraits(size_t npts, double c, double step_size) : + npts_(npts), step_size_(step_size), c_(c) { } LindhMalmqvistGagliardiRadialTraits(const self_type&) = default; - self_type& operator=(const self_type&) = default; + + size_t npts() const noexcept { return npts_; } + + std::unique_ptr clone() const { + return std::make_unique(*this); + } + + bool compare(const RadialTraits& other) const noexcept { + auto ptr = dynamic_cast(&other); + return ptr ? *this == *ptr : false; + } + + bool operator==(const self_type& other) const noexcept { + return npts_ == other.npts_ and c_ == other.c_ and step_size_ == other.step_size_; + } template @@ -110,12 +126,11 @@ class LindhMalmqvistGagliardiRadialTraits { typename BaseQuad::weight_container& weights) const { using point_type = typename BaseQuad::point_type; - assert(points.size() > 2); - const auto npts = points.size() - 2; - point_type up = npts * step_size_; + assert(points.size() - 2 == npts()); + point_type up = npts() * step_size_; // Scale points and weights - for(auto& x : points ) x *= step_size_ * (npts + 1); + for(auto& x : points ) x *= step_size_ * (npts() + 1); for(auto& w : weights) w = step_size_; } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 699987a..2e651e1 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -481,13 +481,13 @@ TEST_CASE("LMG", "[1d-quad]") { // alpha_min, alpha_max, h, c, n); - LindhMalmqvistGagliardiRadialTraits traits(c, h); - LindhMalmqvistGagliardi quad(n, traits); + LindhMalmqvistGagliardiRadialTraits traits(n, c, h); + LindhMalmqvistGagliardi quad(traits); for(int i = 0; i < n; ++i) { double x = (i+1) * h; double r = c * (std::exp(x) - 1.0); double w = h*(r + c); - //printf("%d R = %.16e W = %.16e \n", i, r - quad.points()[i], w-quad.weights()[i]); + printf("%d R = %.16e W = %.16e \n", i, r - quad.points()[i], w-quad.weights()[i]); REQUIRE_THAT(quad.points()[i], Catch::Matchers::WithinRel(r)); REQUIRE_THAT(quad.weights()[i], Catch::Matchers::WithinRel(w)); } From 15ed81930370779596bfe6aa876d23b6670479c2 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 3 Jan 2024 14:26:29 -0800 Subject: [PATCH 19/19] Add LMG to runtime generator --- .../generators/impl/pruned_grid.hpp | 2 + .../generators/impl/radial_factory.hpp | 4 + .../generators/impl/radial_types.hpp | 1 + .../generators/impl/unpruned_grid.hpp | 3 + .../generators/radial_factory.hpp | 14 ++- .../integratorxx/quadratures/radial/lmg.hpp | 10 ++ test/1d_quadratures.cxx | 4 +- test/spherical_generator.cxx | 93 ++++++++++++++----- 8 files changed, 102 insertions(+), 29 deletions(-) diff --git a/include/integratorxx/generators/impl/pruned_grid.hpp b/include/integratorxx/generators/impl/pruned_grid.hpp index c69e036..8336c31 100644 --- a/include/integratorxx/generators/impl/pruned_grid.hpp +++ b/include/integratorxx/generators/impl/pruned_grid.hpp @@ -73,6 +73,8 @@ SphericalGridFactory::spherical_grid_ptr return detail::make_pruned_grid( mhl_type(traits), pruning_regions ); case RadialQuad::TreutlerAhlrichs: return detail::make_pruned_grid( ta_type(traits), pruning_regions ); + case RadialQuad::LindhMalmqvistGagliardi: + return detail::make_pruned_grid( lmg_type(traits), pruning_regions ); default: throw std::runtime_error("Unsupported Radial Quadrature"); diff --git a/include/integratorxx/generators/impl/radial_factory.hpp b/include/integratorxx/generators/impl/radial_factory.hpp index ed89d56..30f55ec 100644 --- a/include/integratorxx/generators/impl/radial_factory.hpp +++ b/include/integratorxx/generators/impl/radial_factory.hpp @@ -16,6 +16,8 @@ RadialQuad radial_from_string(std::string name) { if(name == "MHL") return RadialQuad::MurrayHandyLaming; if(name == "TREUTLERAHLRICHS") return RadialQuad::TreutlerAhlrichs; if(name == "TA") return RadialQuad::TreutlerAhlrichs; + if(name == "LINDHMALMQVISTGAGLIARDI") return RadialQuad::LindhMalmqvistGagliardi; + if(name == "LMG") return RadialQuad::LindhMalmqvistGagliardi; throw std::runtime_error("Unrecognized Radial Quadrature"); } @@ -31,6 +33,8 @@ RadialFactory::radial_grid_ptr RadialFactory::generate(RadialQuad rq, const Radi return std::make_unique(traits); case RadialQuad::TreutlerAhlrichs: return std::make_unique(traits); + case RadialQuad::LindhMalmqvistGagliardi: + return std::make_unique(traits); default: throw std::runtime_error("Unsupported Radial Quadrature"); abort(); diff --git a/include/integratorxx/generators/impl/radial_types.hpp b/include/integratorxx/generators/impl/radial_types.hpp index 4edb1da..6be87a7 100644 --- a/include/integratorxx/generators/impl/radial_types.hpp +++ b/include/integratorxx/generators/impl/radial_types.hpp @@ -7,6 +7,7 @@ using bk_type = Becke; using mk_type = MuraKnowles; using mhl_type = MurrayHandyLaming; using ta_type = TreutlerAhlrichs; +using lmg_type = LindhMalmqvistGagliardi; } diff --git a/include/integratorxx/generators/impl/unpruned_grid.hpp b/include/integratorxx/generators/impl/unpruned_grid.hpp index 758b1fd..9597608 100644 --- a/include/integratorxx/generators/impl/unpruned_grid.hpp +++ b/include/integratorxx/generators/impl/unpruned_grid.hpp @@ -27,6 +27,9 @@ auto generate_unpruned_grid_impl(RadialQuad rq, const RadialTraits& traits, case RadialQuad::TreutlerAhlrichs: return SphericalGridFactory::generate_unpruned_grid( ta_type(traits), std::forward(ang_quad) ); + case RadialQuad::LindhMalmqvistGagliardi: + return SphericalGridFactory::generate_unpruned_grid( lmg_type(traits), std::forward(ang_quad) ); + default: throw std::runtime_error("Unsupported Radial Quadrature"); abort(); diff --git a/include/integratorxx/generators/radial_factory.hpp b/include/integratorxx/generators/radial_factory.hpp index 2d4f0ab..3004a61 100644 --- a/include/integratorxx/generators/radial_factory.hpp +++ b/include/integratorxx/generators/radial_factory.hpp @@ -5,10 +5,11 @@ namespace IntegratorXX { /// High-level specification of radial quadratures enum class RadialQuad : uint32_t { - Becke = 0x0010, - MurrayHandyLaming = 0x0020, - MuraKnowles = 0x0030, - TreutlerAhlrichs = 0x0040 + Becke = 0x0010, + MurrayHandyLaming = 0x0020, + MuraKnowles = 0x0030, + TreutlerAhlrichs = 0x0040, + LindhMalmqvistGagliardi = 0x0050 }; template @@ -17,6 +18,7 @@ RadialQuad radial_from_type() { if constexpr (detail::is_mk_v ) return RadialQuad::MuraKnowles; if constexpr (detail::is_mhl_v) return RadialQuad::MurrayHandyLaming; if constexpr (detail::is_ta_v) return RadialQuad::TreutlerAhlrichs; + if constexpr (detail::is_lmg_v) return RadialQuad::LindhMalmqvistGagliardi; throw std::runtime_error("Unrecognized Radial Quadrature"); }; @@ -55,6 +57,10 @@ std::unique_ptr make_radial_traits(RadialQuad rq, Args&&... args) ptr = detail::make_radial_traits(std::forward(args)...); break; + case RadialQuad::LindhMalmqvistGagliardi: + ptr = + detail::make_radial_traits(std::forward(args)...); + break; } if(!ptr) throw std::runtime_error("RadialTraits Construction Failed"); diff --git a/include/integratorxx/quadratures/radial/lmg.hpp b/include/integratorxx/quadratures/radial/lmg.hpp index f353092..d4bf445 100644 --- a/include/integratorxx/quadratures/radial/lmg.hpp +++ b/include/integratorxx/quadratures/radial/lmg.hpp @@ -141,4 +141,14 @@ using LindhMalmqvistGagliardi = RadialTransformQuadrature< UniformTrapezoid, LindhMalmqvistGagliardiRadialTraits >; +namespace detail { + +template +static constexpr bool is_lmg_v = std::is_same_v< + QuadType, + LindhMalmqvistGagliardi +>; + +} + } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 2e651e1..a93f991 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -477,7 +477,7 @@ TEST_CASE("LMG", "[1d-quad]") { double c = r_lower / (std::exp(h) - 1); int n = std::log(1.0 + r_upper/c) / h + 1.0; - //printf("ALPHA_MIN = %.4e ALPHA_MAX = %.4e H = %.4e C = %.6e N = %d\n", + //printf("ALPHA_MIN = %.4e ALPHA_MAX = %.4e H = %.4e C = %.16e N = %d\n", // alpha_min, alpha_max, h, c, n); @@ -487,7 +487,7 @@ TEST_CASE("LMG", "[1d-quad]") { double x = (i+1) * h; double r = c * (std::exp(x) - 1.0); double w = h*(r + c); - printf("%d R = %.16e W = %.16e \n", i, r - quad.points()[i], w-quad.weights()[i]); + //printf("%d R = %.16e W = %.16e \n", i, r - quad.points()[i], w-quad.weights()[i]); REQUIRE_THAT(quad.points()[i], Catch::Matchers::WithinRel(r)); REQUIRE_THAT(quad.weights()[i], Catch::Matchers::WithinRel(w)); } diff --git a/test/spherical_generator.cxx b/test/spherical_generator.cxx index 8fbe482..ac87fc6 100644 --- a/test/spherical_generator.cxx +++ b/test/spherical_generator.cxx @@ -12,6 +12,7 @@ using bk_type = IntegratorXX::Becke; using mk_type = IntegratorXX::MuraKnowles; using mhl_type = IntegratorXX::MurrayHandyLaming; using ta_type = IntegratorXX::TreutlerAhlrichs; +using lmg_type = IntegratorXX::LindhMalmqvistGagliardi; using ah_type = IntegratorXX::AhrensBeylkin; using de_type = IntegratorXX::Delley; @@ -39,10 +40,11 @@ std::ostream& operator<<( std::ostream& out, using namespace IntegratorXX; out << "RADIAL QUAD: "; switch(p.radial_quad) { - case RadialQuad::Becke: out << "Becke"; break; - case RadialQuad::MuraKnowles: out << "MK"; break; - case RadialQuad::MurrayHandyLaming: out << "MHL"; break; - case RadialQuad::TreutlerAhlrichs: out << "TA"; break; + case RadialQuad::Becke: out << "Becke"; break; + case RadialQuad::MuraKnowles: out << "MK"; break; + case RadialQuad::MurrayHandyLaming: out << "MHL"; break; + case RadialQuad::TreutlerAhlrichs: out << "TA"; break; + case RadialQuad::LindhMalmqvistGagliardi: out << "LMG"; break; } out << std::endl; @@ -64,6 +66,8 @@ TEST_CASE( "String Getter", "[sph-gen]" ) { REQUIRE(radial_from_type() == radial_from_string("MHL")); REQUIRE(radial_from_type() == radial_from_string("TreutlerAhlrichs")); REQUIRE(radial_from_type() == radial_from_string("TA")); + REQUIRE(radial_from_type() == radial_from_string("LMG")); + REQUIRE(radial_from_type() == radial_from_string("LindhMalmqvistGagliardi")); } SECTION("Angular") { @@ -234,17 +238,29 @@ TEST_CASE( "Pruning Schemes", "[sph-gen]" ) { } using radial_test_types = std::tuple< - bk_type, mk_type, mhl_type, ta_type + bk_type, mk_type, mhl_type, ta_type, lmg_type >; TEMPLATE_LIST_TEST_CASE("Radial Generator", "[sph-gen]", radial_test_types) { using namespace IntegratorXX; using radial_type = TestType; - size_t npts = 10; - radial_type rq(npts, 1.0); + size_t npts; + std::unique_ptr rad_traits = nullptr; + std::unique_ptr rq_ptr = nullptr; auto rad_spec = radial_from_type(); - auto rad_traits = make_radial_traits(rad_spec, npts, 1.0); + if constexpr (std::is_same_v) { + npts = 74; + double h = 0.15235041341729277; + double c = 1.2798590387343324e-04; + rad_traits = make_radial_traits(rad_spec, npts, c, h); + rq_ptr = std::make_unique(*rad_traits); + } else { + npts = 10; + rq_ptr = std::make_unique(npts, 1.0); + rad_traits = make_radial_traits(rad_spec, npts, 1.0); + } + const auto& rq = *rq_ptr; auto rad_grid = RadialFactory::generate(rad_spec, *rad_traits); REQUIRE(rad_grid->npts() == npts); @@ -309,7 +325,12 @@ using sph_test_types = std::tuple< std::tuple, std::tuple, std::tuple, - std::tuple + std::tuple, + + std::tuple, + std::tuple, + std::tuple, + std::tuple >; @@ -322,19 +343,33 @@ TEMPLATE_LIST_TEST_CASE("Unpruned", "[sph-gen]", sph_test_types) { using spherical_type = SphericalQuadrature; - size_t nrad = 10; + + size_t nrad; + std::unique_ptr rad_traits = nullptr; + std::unique_ptr rq_ptr = nullptr; + auto rad_spec = radial_from_type(); + if constexpr (std::is_same_v) { + nrad = 74; + double h = 0.15235041341729277; + double c = 1.2798590387343324e-04; + rad_traits = make_radial_traits(rad_spec, nrad, c, h); + rq_ptr = std::make_unique(*rad_traits); + } else { + nrad = 10; + rq_ptr = std::make_unique(nrad, 1.0); + rad_traits = make_radial_traits(rad_spec, nrad, 1.0); + } + size_t nang = angular_traits::npts_by_algebraic_order( angular_traits::next_algebraic_order(1)); // Smallest possible angular grid + angular_type aq(nang); + // Generate the quadrature manually - radial_type rq(nrad, 1.0); - angular_type aq(nang); - spherical_type sph_ref(rq,aq); + spherical_type sph_ref(*rq_ptr,aq); // Generate via runtime API - auto rad_spec = radial_from_type(); - auto rad_traits = make_radial_traits(rad_spec, nrad, 1.0); UnprunedSphericalGridSpecification unp( rad_spec, *rad_traits, angular_from_type(), nang @@ -368,13 +403,26 @@ TEMPLATE_LIST_TEST_CASE("Pruned", "[sph-gen]", sph_test_types) { using spherical_type = PrunedSphericalQuadrature; - size_t nrad = 99; + size_t nrad; + std::unique_ptr rad_traits = nullptr; + std::unique_ptr rq_ptr = nullptr; + auto rad_spec = radial_from_type(); + if constexpr (std::is_same_v) { + nrad = 74; + double h = 0.15235041341729277; + double c = 1.2798590387343324e-04; + rad_traits = make_radial_traits(rad_spec, nrad, c, h); + rq_ptr = std::make_unique(*rad_traits); + } else { + nrad = 99; + rq_ptr = std::make_unique(nrad, 1.0); + rad_traits = make_radial_traits(rad_spec, nrad, 1.0); + } + size_t nang = angular_traits::npts_by_algebraic_order( angular_traits::next_algebraic_order(29)); // Smallest possible angular grid // Generate pruning scheme - auto rad_spec = radial_from_type(); - auto rad_traits = make_radial_traits(rad_spec, nrad, 1.0); UnprunedSphericalGridSpecification unp( rad_spec, *rad_traits, angular_from_type(), nang @@ -384,15 +432,14 @@ TEMPLATE_LIST_TEST_CASE("Pruned", "[sph-gen]", sph_test_types) { // Generate the quadrature manually - radial_type rq(nrad, 1.0); RadialGridPartition rgp; for(auto pr : pruning_spec.pruning_regions) { angular_type aq(pr.angular_size); - rgp.add_quad(rq, pr.idx_st, aq); + rgp.add_quad(*rq_ptr, pr.idx_st, aq); } - rgp.finalize(rq); + rgp.finalize(*rq_ptr); - spherical_type sph_ref(rq, rgp); + spherical_type sph_ref(*rq_ptr, rgp); // Generate via runtime API @@ -414,5 +461,5 @@ TEMPLATE_LIST_TEST_CASE("Pruned", "[sph-gen]", sph_test_types) { auto w_ref = sph_ref.weights()[i]; REQUIRE_THAT(w, Catch::Matchers::WithinAbs(w_ref,1e-15)); } - + }