From 20683c967d352a6b82b9fa26dfa321073ad96b22 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 6 Jul 2023 17:14:29 -0700 Subject: [PATCH 01/15] Starting modularization of radial quadratures, demonstrated on MHL --- include/integratorxx/quadratures/mhl.hpp | 28 +++++++++ .../quadratures/radial_transform.hpp | 61 +++++++++++++++++++ include/integratorxx/quadratures/uniform.hpp | 4 ++ test/1d_quadratures.cxx | 7 ++- 4 files changed, 98 insertions(+), 2 deletions(-) create mode 100644 include/integratorxx/quadratures/radial_transform.hpp diff --git a/include/integratorxx/quadratures/mhl.hpp b/include/integratorxx/quadratures/mhl.hpp index cec02dc..e469fa0 100644 --- a/include/integratorxx/quadratures/mhl.hpp +++ b/include/integratorxx/quadratures/mhl.hpp @@ -2,6 +2,7 @@ #include #include +#include namespace IntegratorXX { @@ -140,4 +141,31 @@ struct quadrature_traits< }; + + + + + +template +struct MHLRadialTraits { + + template + static auto radial_transform(PointType x) { + return std::pow( x / (1.0 - x), M ); + } + + template + static auto radial_jacobian(PointType x) { + return M * std::pow(x, M-1) / std::pow(1.0 - x, M+1); + } + +}; + + +template +using MHLAuto = RadialTransformQuadrature< + UniformTrapezoid, + MHLRadialTraits<2> +>; + } diff --git a/include/integratorxx/quadratures/radial_transform.hpp b/include/integratorxx/quadratures/radial_transform.hpp new file mode 100644 index 0000000..e43e928 --- /dev/null +++ b/include/integratorxx/quadratures/radial_transform.hpp @@ -0,0 +1,61 @@ +#pragma once + +#include + +namespace IntegratorXX { + +template +struct RadialTransformQuadrature : + public Quadrature> { + + using base_type = Quadrature>; + +public: + + using point_type = typename base_type::point_type; + using weight_type = typename base_type::weight_type; + using point_container = typename base_type::point_container; + using weight_container = typename base_type::weight_container; + + RadialTransformQuadrature(size_t npts, weight_type R = 1.0) : + base_type( npts, R ) { } + + RadialTransformQuadrature( const RadialTransformQuadrature& ) = default; + RadialTransformQuadrature( RadialTransformQuadrature&& ) noexcept = default; + +}; + +template +struct quadrature_traits> { + + using point_type = typename BaseQuad::point_type; + using weight_type = typename BaseQuad::weight_type; + using point_container = typename BaseQuad::point_container; + using weight_container = typename BaseQuad::weight_container; + + inline static std::tuple + generate( size_t npts, weight_type R ) { + + using base_quad_traits = quadrature_traits; + + point_container points( npts ); + weight_container weights( npts ); + + + auto [base_x, base_w] = base_quad_traits::generate(npts+2); + + for(size_t i = 0; i < npts; ++i) { + const auto xi = base_x[i+1]; + const auto wi = base_w[i+1]; + points[i] = R * RadialTraits::radial_transform(xi); + weights[i] = R * wi * RadialTraits::radial_jacobian(xi); + } + + return std::make_tuple(points, weights); + } + +}; + + + +} diff --git a/include/integratorxx/quadratures/uniform.hpp b/include/integratorxx/quadratures/uniform.hpp index e173fe2..27720f1 100644 --- a/include/integratorxx/quadratures/uniform.hpp +++ b/include/integratorxx/quadratures/uniform.hpp @@ -42,6 +42,8 @@ struct quadrature_traits< using point_container = std::vector< point_type >; using weight_container = std::vector< weight_type >; + inline static constexpr bool bound_inclusive = true; + inline static std::tuple generate( size_t npts, point_type lo, point_type up ) { @@ -60,6 +62,8 @@ struct quadrature_traits< } + inline static std::tuple + generate( size_t npts ) { return generate(npts, 0.0, 1.0); } }; } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index c7a28b6..6409cd9 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -166,9 +166,12 @@ TEST_CASE( "Gauss-Legendre Quadratures", "[1d-quad]" ) { TEST_CASE( "Euler-Maclaurin Quadratures", "[1d-quad]" ) { IntegratorXX::MurrayHandyLaming quad( 150 ); + IntegratorXX::MHLAuto quad2(150); - const auto& pts = quad.points(); - const auto& wgt = quad.weights(); + //const auto& pts = quad.points(); + //const auto& wgt = quad.weights(); + const auto& pts = quad2.points(); + const auto& wgt = quad2.weights(); auto f = [=]( double x ){ return gaussian(x); }; From a14393d2d448dd7270d9f4652499f1d5e9da5eb4 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 6 Jul 2023 18:22:30 -0700 Subject: [PATCH 02/15] Moved over MHL entirely over to auto generator, starting working on TA --- .../integratorxx/quadratures/gausscheby2.hpp | 2 +- include/integratorxx/quadratures/mhl.hpp | 11 +++++-- .../quadratures/radial_transform.hpp | 8 +++-- .../quadratures/treutlerahlrichs.hpp | 29 ++++++++++++++++++- test/1d_quadratures.cxx | 9 ++---- 5 files changed, 45 insertions(+), 14 deletions(-) diff --git a/include/integratorxx/quadratures/gausscheby2.hpp b/include/integratorxx/quadratures/gausscheby2.hpp index 8926e04..6b47c20 100644 --- a/include/integratorxx/quadratures/gausscheby2.hpp +++ b/include/integratorxx/quadratures/gausscheby2.hpp @@ -72,6 +72,6 @@ struct quadrature_traits> { return std::make_tuple(points, weights); } -} +}; } // namespace IntegratorXX diff --git a/include/integratorxx/quadratures/mhl.hpp b/include/integratorxx/quadratures/mhl.hpp index e469fa0..880cb8b 100644 --- a/include/integratorxx/quadratures/mhl.hpp +++ b/include/integratorxx/quadratures/mhl.hpp @@ -6,6 +6,7 @@ namespace IntegratorXX { +#if 0 /** * @brief Implementation of the Murray-Handy-Laming radial quadrature. * @@ -140,6 +141,7 @@ struct quadrature_traits< } }; +#else @@ -147,7 +149,7 @@ struct quadrature_traits< template -struct MHLRadialTraits { +struct MurrayHandyLamingRadialTraits { template static auto radial_transform(PointType x) { @@ -163,9 +165,12 @@ struct MHLRadialTraits { template -using MHLAuto = RadialTransformQuadrature< +using MurrayHandyLaming = RadialTransformQuadrature< UniformTrapezoid, - MHLRadialTraits<2> + MurrayHandyLamingRadialTraits<2> >; +#endif + } + diff --git a/include/integratorxx/quadratures/radial_transform.hpp b/include/integratorxx/quadratures/radial_transform.hpp index e43e928..18183e6 100644 --- a/include/integratorxx/quadratures/radial_transform.hpp +++ b/include/integratorxx/quadratures/radial_transform.hpp @@ -42,11 +42,13 @@ struct quadrature_traits> { weight_container weights( npts ); - auto [base_x, base_w] = base_quad_traits::generate(npts+2); + const auto npts_base = base_quad_traits::bound_inclusive ? npts+2 : npts; + auto [base_x, base_w] = base_quad_traits::generate(npts_base); + const auto ipts_offset = !!base_quad_traits::bound_inclusive; for(size_t i = 0; i < npts; ++i) { - const auto xi = base_x[i+1]; - const auto wi = base_w[i+1]; + const auto xi = base_x[i + ipts_offset]; + const auto wi = base_w[i + ipts_offset]; points[i] = R * RadialTraits::radial_transform(xi); weights[i] = R * wi * RadialTraits::radial_jacobian(xi); } diff --git a/include/integratorxx/quadratures/treutlerahlrichs.hpp b/include/integratorxx/quadratures/treutlerahlrichs.hpp index 48ceb91..5610507 100644 --- a/include/integratorxx/quadratures/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/treutlerahlrichs.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include namespace IntegratorXX { @@ -115,12 +116,38 @@ struct quadrature_traits< ( sqrt_term - alpha * log_term / sqrt_term ); } + return std::make_tuple( points, weights ); + + } + +}; + - return std::make_tuple( points, weights ); + +struct TreutlerAhlrichsRadialTraits { + + inline static constexpr double alpha = 0.6; + inline static constexpr double ln_2 = 0.693147180559945309417232; + + template + static auto radial_transform(PointType x) { + const auto pow_term = std::pow(1.0 + x, alpha); + const auto log_term = std::log(0.5 * (1.0 - x)); + return pow_term * log_term / ln_2; + }; + + + template + static auto radial_jacobian(PointType x) { + const auto pow_term = std::pow(1.0 + x, alpha); + const auto log_term = std::log(0.5 * (1.0 - x)); + return pow_term / ln_2 * (-alpha * log_term / (1.0 + x) + 1.0 / (1.0 - x)); } }; + + } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 6409cd9..9864e0b 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -165,13 +165,10 @@ TEST_CASE( "Gauss-Legendre Quadratures", "[1d-quad]" ) { TEST_CASE( "Euler-Maclaurin Quadratures", "[1d-quad]" ) { - IntegratorXX::MurrayHandyLaming quad( 150 ); - IntegratorXX::MHLAuto quad2(150); + IntegratorXX::MurrayHandyLaming quad(150); - //const auto& pts = quad.points(); - //const auto& wgt = quad.weights(); - const auto& pts = quad2.points(); - const auto& wgt = quad2.weights(); + const auto& pts = quad.points(); + const auto& wgt = quad.weights(); auto f = [=]( double x ){ return gaussian(x); }; From 034f8bf6af4039f50f83fe37591efd59d1b01b5b Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Fri, 7 Jul 2023 18:41:37 -0700 Subject: [PATCH 03/15] Fixup of TA radial traits --- .../quadratures/treutlerahlrichs.hpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/include/integratorxx/quadratures/treutlerahlrichs.hpp b/include/integratorxx/quadratures/treutlerahlrichs.hpp index 5610507..ad6d65b 100644 --- a/include/integratorxx/quadratures/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/treutlerahlrichs.hpp @@ -126,28 +126,27 @@ struct quadrature_traits< -struct TreutlerAhlrichsRadialTraits { +struct TreutlerAhlrichsM4RadialTraits { inline static constexpr double alpha = 0.6; + inline static constexpr double a = 1.0; inline static constexpr double ln_2 = 0.693147180559945309417232; template static auto radial_transform(PointType x) { - const auto pow_term = std::pow(1.0 + x, alpha); - const auto log_term = std::log(0.5 * (1.0 - x)); - return pow_term * log_term / ln_2; + const auto pow_term = std::pow(a + x, alpha); + const auto log_term = std::log((a + 1.0) / (1.0 - x)); + return pow_term * log_term / ln_2; }; template static auto radial_jacobian(PointType x) { - const auto pow_term = std::pow(1.0 + x, alpha); - const auto log_term = std::log(0.5 * (1.0 - x)); - return pow_term / ln_2 * (-alpha * log_term / (1.0 + x) + 1.0 / (1.0 - x)); + const auto pow_term = std::pow(a + x, alpha); + const auto log_term = std::log((a + 1.0) / (1.0 - x)); + return pow_term / ln_2 * ( alpha * log_term / (a+x) + (1./(1. - x)) ); } }; - - } From 4e014410a0360dc41bf45c9fdced48a85b77e3ed Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sat, 8 Jul 2023 16:59:39 -0700 Subject: [PATCH 04/15] Added Treutler-Ahlrichs RadialTransformQuadrature, confirmed the same as hand-coded variant --- include/integratorxx/quadratures/gausscheby2.hpp | 8 ++++++++ include/integratorxx/quadratures/treutlerahlrichs.hpp | 7 +++++++ test/1d_quadratures.cxx | 1 + 3 files changed, 16 insertions(+) diff --git a/include/integratorxx/quadratures/gausscheby2.hpp b/include/integratorxx/quadratures/gausscheby2.hpp index 8121112..bfd2b4d 100644 --- a/include/integratorxx/quadratures/gausscheby2.hpp +++ b/include/integratorxx/quadratures/gausscheby2.hpp @@ -56,6 +56,8 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; + + inline static constexpr bool bound_inclusive = false; inline static std::tuple generate(size_t npts, point_type lo, point_type up) { @@ -82,6 +84,12 @@ struct quadrature_traits> { return std::make_tuple(points, weights); } + + inline static std::tuple + generate(size_t npts) { + return generate(npts, -1.0, 1.0); + } + }; } // namespace IntegratorXX diff --git a/include/integratorxx/quadratures/treutlerahlrichs.hpp b/include/integratorxx/quadratures/treutlerahlrichs.hpp index ad6d65b..ba907e5 100644 --- a/include/integratorxx/quadratures/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/treutlerahlrichs.hpp @@ -149,4 +149,11 @@ struct TreutlerAhlrichsM4RadialTraits { }; + +template +using TAAuto = RadialTransformQuadrature< + GaussChebyshev2, + TreutlerAhlrichsM4RadialTraits +>; + } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index e88ab88..5465856 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -223,6 +223,7 @@ TEST_CASE( "Euler-Maclaurin Quadratures", "[1d-quad]" ) { TEST_CASE( "Ahlrichs Quadratures", "[1d-quad]" ) { IntegratorXX::TreutlerAhlrichs quad( 150 ); + IntegratorXX::TAAuto quad2(150); const auto& pts = quad.points(); const auto& wgt = quad.weights(); From 785f6ccd29379993fd8c963d4a3a7f966adc469c Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sat, 8 Jul 2023 17:03:03 -0700 Subject: [PATCH 05/15] Moved TA implementation over to RadialTransformQuadrature --- include/integratorxx/quadratures/treutlerahlrichs.hpp | 8 +++++--- test/1d_quadratures.cxx | 5 ++--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/include/integratorxx/quadratures/treutlerahlrichs.hpp b/include/integratorxx/quadratures/treutlerahlrichs.hpp index ba907e5..026611a 100644 --- a/include/integratorxx/quadratures/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/treutlerahlrichs.hpp @@ -6,6 +6,7 @@ namespace IntegratorXX { +#if 0 /** * @brief Implementation of the Treutler-Ahlrichs radial quadrature. * @@ -123,8 +124,7 @@ struct quadrature_traits< }; - - +#else struct TreutlerAhlrichsM4RadialTraits { @@ -151,9 +151,11 @@ struct TreutlerAhlrichsM4RadialTraits { template -using TAAuto = RadialTransformQuadrature< +using TreutlerAhlrichs = RadialTransformQuadrature< GaussChebyshev2, TreutlerAhlrichsM4RadialTraits >; +#endif + } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 5465856..85b4578 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -220,10 +220,9 @@ TEST_CASE( "Euler-Maclaurin Quadratures", "[1d-quad]" ) { } -TEST_CASE( "Ahlrichs Quadratures", "[1d-quad]" ) { +TEST_CASE( "Treutler-Ahlrichs Quadratures", "[1d-quad]" ) { - IntegratorXX::TreutlerAhlrichs quad( 150 ); - IntegratorXX::TAAuto quad2(150); + IntegratorXX::TreutlerAhlrichs quad(150); const auto& pts = quad.points(); const auto& wgt = quad.weights(); From 81529511a7588246afed976c9927b9ae29b209d1 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sat, 8 Jul 2023 17:13:59 -0700 Subject: [PATCH 06/15] Added MuraKnowles RadialTransformQuadrature, confirmed the same as hand-coded variant --- .../integratorxx/quadratures/muraknowles.hpp | 23 +++++++++++++++++++ test/1d_quadratures.cxx | 4 +++- 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/include/integratorxx/quadratures/muraknowles.hpp b/include/integratorxx/quadratures/muraknowles.hpp index 5847ce0..bb9b3d6 100644 --- a/include/integratorxx/quadratures/muraknowles.hpp +++ b/include/integratorxx/quadratures/muraknowles.hpp @@ -2,6 +2,7 @@ #include #include +#include namespace IntegratorXX { @@ -132,4 +133,26 @@ struct quadrature_traits< }; + +struct MuraKnowlesRadialTraits { + + template + static auto radial_transform(PointType x) { + return -std::log(1.0 - x*x*x); + } + + template + static auto radial_jacobian(PointType x) { + const auto x2 = x*x; + return 3.0 * x2 / (1.0 - x2 * x); + } + +}; + +template +using MKAuto = RadialTransformQuadrature< + UniformTrapezoid, + MuraKnowlesRadialTraits +>; + } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 85b4578..62202bb 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -241,6 +241,7 @@ TEST_CASE( "Treutler-Ahlrichs Quadratures", "[1d-quad]" ) { TEST_CASE( "Knowles Quadratures", "[1d-quad]" ) { IntegratorXX::MuraKnowles quad( 150 ); + IntegratorXX::MKAuto quad2( 150 ); const auto& pts = quad.points(); const auto& wgt = quad.weights(); @@ -248,8 +249,9 @@ TEST_CASE( "Knowles Quadratures", "[1d-quad]" ) { auto f = [=]( double x ){ return gaussian(x); }; double res = 0.; - for( auto i = 0; i < quad.npts(); ++i ) + for( auto i = 0; i < quad.npts(); ++i ) { res += wgt[i] * f(pts[i]); + CHECK( res == Catch::Approx(ref_gaussian_int(0.,inf)) ); From 55f2f88ac22a860e320674572e83b86b47047531 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sat, 8 Jul 2023 17:17:52 -0700 Subject: [PATCH 07/15] Moved MK implementation over to RadialTransformQuadrature --- include/integratorxx/quadratures/muraknowles.hpp | 6 +++++- test/1d_quadratures.cxx | 5 ++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/include/integratorxx/quadratures/muraknowles.hpp b/include/integratorxx/quadratures/muraknowles.hpp index bb9b3d6..06b9588 100644 --- a/include/integratorxx/quadratures/muraknowles.hpp +++ b/include/integratorxx/quadratures/muraknowles.hpp @@ -6,6 +6,7 @@ namespace IntegratorXX { +#if 0 /** * @brief Implementation of the Mura-Knowles radial quadrature. * @@ -133,6 +134,7 @@ struct quadrature_traits< }; +#else struct MuraKnowlesRadialTraits { @@ -150,9 +152,11 @@ struct MuraKnowlesRadialTraits { }; template -using MKAuto = RadialTransformQuadrature< +using MuraKnowles = RadialTransformQuadrature< UniformTrapezoid, MuraKnowlesRadialTraits >; +#endif + } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 62202bb..12dbd15 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -240,8 +240,7 @@ TEST_CASE( "Treutler-Ahlrichs Quadratures", "[1d-quad]" ) { TEST_CASE( "Knowles Quadratures", "[1d-quad]" ) { - IntegratorXX::MuraKnowles quad( 150 ); - IntegratorXX::MKAuto quad2( 150 ); + IntegratorXX::MuraKnowles quad(150); const auto& pts = quad.points(); const auto& wgt = quad.weights(); @@ -249,7 +248,7 @@ TEST_CASE( "Knowles Quadratures", "[1d-quad]" ) { auto f = [=]( double x ){ return gaussian(x); }; double res = 0.; - for( auto i = 0; i < quad.npts(); ++i ) { + for( auto i = 0; i < quad.npts(); ++i ) res += wgt[i] * f(pts[i]); From c9936ea8a85a83e942135b4378299bd77dfe5522 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sat, 8 Jul 2023 17:35:45 -0700 Subject: [PATCH 08/15] Add dox to new MHL implemenatation --- include/integratorxx/quadratures/mhl.hpp | 171 +++++------------------ 1 file changed, 34 insertions(+), 137 deletions(-) diff --git a/include/integratorxx/quadratures/mhl.hpp b/include/integratorxx/quadratures/mhl.hpp index 880cb8b..0e684cc 100644 --- a/include/integratorxx/quadratures/mhl.hpp +++ b/include/integratorxx/quadratures/mhl.hpp @@ -6,156 +6,37 @@ namespace IntegratorXX { -#if 0 /** - * @brief Implementation of the Murray-Handy-Laming radial quadrature. - * - * Generates a quadrature on the bounds (0, inf). Suitable for integrands - * which tend to zero as their argument tends to 0 and inf. Tailored for - * radial integrands, i.e. r^2 * f(r), with lim_{r->inf} f(r) = 0. + * @brief Implementation of the Murray-Handy-Laming radial quadrature + * transformation rules. * * Reference: * Molecular Physics, 78:4, 997-1014, * DOI: https://doi.org/10.1080/00268979300100651 * - * @tparam PointType Type describing the quadrature points - * @tparam WeightType Type describing the quadrature weights - */ -template -class MurrayHandyLaming : - public Quadrature> { - - using base_type = Quadrature>; - -public: - - using point_type = typename base_type::point_type; - using weight_type = typename base_type::weight_type; - using point_container = typename base_type::point_container; - using weight_container = typename base_type::weight_container; - - /** - * @brief Construct the Murray-Handy-Laming radial quadrature - * - * @param[in] npts Number of quadrature points to generate - * @param[in] R Radial scaling factor. Suggested to be set to - * correspond to the Bragg-Slater radius of the - * appropriate nucleus - * @param[in] m Exponent factor for the quadrature. 2 was the - * default suggested in the reference. - */ - MurrayHandyLaming(size_t npts, weight_type R = 1., int m = 2): - base_type( npts, R, m ) { } - - MurrayHandyLaming( const MurrayHandyLaming& ) = default; - MurrayHandyLaming( MurrayHandyLaming&& ) noexcept = default; -}; - - - - - - -/** - * @brief Quadrature traits for the Murray-Handy-Laming quadrature - * - * @tparam PointType Type describing the quadrature points - * @tparam WeightType Type describing the quadrature weights + * @tparam M Integer to modulate the MHL transformation. + * Typically taken to be 2. */ -template -struct quadrature_traits< - MurrayHandyLaming, - std::enable_if_t< - std::is_floating_point_v and - std::is_floating_point_v - > -> { - - using point_type = PointType; - using weight_type = WeightType; - - using point_container = std::vector< point_type >; - using weight_container = std::vector< weight_type >; - - - /** - * @brief Generator for the Murray-Handy-Laming quadrature - * - * @param[in] npts Number of quadrature points to generate - * @param[in] R Radial scaling factor. Suggested to be set to - * correspond to the Bragg-Slater radius of the - * appropriate nucleus - * @param[in] m Exponent fector for the quadrature. 2 was the - * default suggested in the original reference. - * - * @returns Tuple of quadrature points and weights - */ - inline static std::tuple - generate( size_t npts, weight_type R, int m ) { - - point_container points( npts ); - weight_container weights( npts ); - - - using base_quad_traits = - quadrature_traits>; - - /* - * Generate uniform trapezoid points on [0,1] - * ux(j) = j/(m-1) - * uw(j) = 1/(m-1) - * m = npts + 2, j in [0, npts+2) - */ - auto [ux, uw] = base_quad_traits::generate( npts+2, 0., 1. ); - - /* - * Perform Murray, Handy and Laming transformation - * - * Original reference: - * Molecular Physics, 78:4, 997-1014, - * DOI: https://doi.org/10.1080/00268979300100651 - * - * Closed form for points / weights obtained from: - * Journal of Computational Chemistry, 24: 732–740, 2003 - * DOI: https://doi.org/10.1002/jcc.10211 - * - * x(i) = ux(i+1) - * r(i) = [x(i) / (1 - x(i))]^m - * w(i) = uw(i+1) * m * x(i)^(m-1) * [1 - x(i)]^(-m-1) - * i in [0, npts) - * - * XXX: i+1 offset on trapezoid points ignores enpoints of trapezoid - * quadrature (f(r) = 0 with r in {0,inf}) - */ - for( size_t i = 0; i < npts; ++i ) { - const auto xi = ux[i+1]; - const auto one_m_xi = 1. - xi; - points[i] = R * std::pow( xi / one_m_xi, m ); - weights[i] = R * uw[i+1] * m * std::pow( xi, m-1 ) / std::pow( one_m_xi, m+1); - } - - - - return std::make_tuple( points, weights ); - - } - -}; -#else - - - - - - template struct MurrayHandyLamingRadialTraits { + /** + * @brief Transformation rule for the MHL radial quadrature + * + * @param[in] x Point in (0,1) + * @return r = (x / (1-x))^M, r \in (0,inf) + */ template static auto radial_transform(PointType x) { return std::pow( x / (1.0 - x), M ); } + /** + * @brief Jacobian of the MHL radial transformation + * + * @param[in] x Point in (0,1) + * @returns dr = M * x^(M-1) / (1-x)^(M+1) + */ template static auto radial_jacobian(PointType x) { return M * std::pow(x, M-1) / std::pow(1.0 - x, M+1); @@ -164,13 +45,29 @@ struct MurrayHandyLamingRadialTraits { }; +/** + * @brief Implementation of the Murray-Handy-Laming radial quadrature. + * + * Taken as the convolution of the Uniform (Trapezoid) quadrature with + * the MHL radial transformation. See MurrayHandyLamingRadialTraits for + * details. + * + * Suitable for integrands which tend to zero as their argument tends + * to 0 and inf. Tailored for radial integrands, i.e. r^2 * f(r), with + * lim_{r->inf} f(r) = 0. + * + * Reference: + * Molecular Physics, 78:4, 997-1014, + * DOI: https://doi.org/10.1080/00268979300100651 + * + * @tparam PointType Type describing the quadrature points + * @tparam WeightType Type describing the quadrature weights + */ template using MurrayHandyLaming = RadialTransformQuadrature< UniformTrapezoid, MurrayHandyLamingRadialTraits<2> >; -#endif - } From 01a1ce471c5367f4582737df25683e8f00b5b78e Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sat, 8 Jul 2023 17:50:00 -0700 Subject: [PATCH 09/15] Add dox to new TA implemenatation --- include/integratorxx/quadratures/mhl.hpp | 4 +- .../quadratures/treutlerahlrichs.hpp | 153 +++++------------- 2 files changed, 39 insertions(+), 118 deletions(-) diff --git a/include/integratorxx/quadratures/mhl.hpp b/include/integratorxx/quadratures/mhl.hpp index 0e684cc..02473b1 100644 --- a/include/integratorxx/quadratures/mhl.hpp +++ b/include/integratorxx/quadratures/mhl.hpp @@ -24,7 +24,7 @@ struct MurrayHandyLamingRadialTraits { * @brief Transformation rule for the MHL radial quadrature * * @param[in] x Point in (0,1) - * @return r = (x / (1-x))^M, r \in (0,inf) + * @return r = (x / (1-x))^M */ template static auto radial_transform(PointType x) { @@ -35,7 +35,7 @@ struct MurrayHandyLamingRadialTraits { * @brief Jacobian of the MHL radial transformation * * @param[in] x Point in (0,1) - * @returns dr = M * x^(M-1) / (1-x)^(M+1) + * @returns dr/dx (see `radial_transform`) */ template static auto radial_jacobian(PointType x) { diff --git a/include/integratorxx/quadratures/treutlerahlrichs.hpp b/include/integratorxx/quadratures/treutlerahlrichs.hpp index 026611a..fe621df 100644 --- a/include/integratorxx/quadratures/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/treutlerahlrichs.hpp @@ -2,136 +2,35 @@ #include #include +#include namespace IntegratorXX { - -#if 0 /** - * @brief Implementation of the Treutler-Ahlrichs radial quadrature. - * - * Generates a quadrature on the bounds (0, inf). Suitable for integrands - * which tend to zero as their argument tends to 0 and inf. Tailored for - * radial integrands, i.e. r^2 * f(r), with lim_{r->inf} f(r) = 0. + * @brief Implementation of the Treutler-Ahlrichs M4 radial + * quadrature transformation rules. * * Reference: * J. Chem. Phys. 102, 346 (1995) * DOI: https://doi.org/10.1063/1.469408 * - * @tparam PointType Type describing the quadrature points - * @tparam WeightType Type describing the quadrature weights - */ -template -class TreutlerAhlrichs : - public Quadrature> { - - using base_type = Quadrature>; - -public: - - using point_type = typename base_type::point_type; - using weight_type = typename base_type::weight_type; - using point_container = typename base_type::point_container; - using weight_container = typename base_type::weight_container; - - /** - * @brief Construct the Treutler-Ahlrichs radial quadrature - * - * @param[in] npts Number of quadrature points to generate - * @param[in] R Radial scaling factor. Table for suggested - * values is given in the original reference. - * @param[in] alpha Exponent factor for the quadrature. 0.6 was the - * default suggested in the reference. - */ - TreutlerAhlrichs(size_t npts, weight_type R = 1., weight_type alpha = 0.6): - base_type( npts, R, alpha ) { } - - TreutlerAhlrichs( const TreutlerAhlrichs& ) = default; - TreutlerAhlrichs( TreutlerAhlrichs&& ) noexcept = default; -}; - - - - - - -/** - * @brief Quadrature traits for the Treutler-Ahlrichs quadrature - * - * @tparam PointType Type describing the quadrature points - * @tparam WeightType Type describing the quadrature weights + * @tparam M Integer to modulate the MHL transformation. + * Typically taken to be 2. */ - -template -struct quadrature_traits< - TreutlerAhlrichs -> { - - using point_type = PointType; - using weight_type = WeightType; - - using point_container = std::vector< point_type >; - using weight_container = std::vector< weight_type >; - - /** - * @brief Generator for the Treutler-Ahlrichs quadrature - * - * @param[in] npts Number of quadrature points to generate - * @param[in] R Radial scaling factor. Table for suggested - * values is given in the original reference. - * @param[in] alpha Exponent factor for the quadrature. 0.6 was the - * - * @returns Tuple of quadrature points and weights - */ - inline static std::tuple - generate( size_t npts, weight_type R, weight_type alpha ) { - - - const point_type ln_2 = std::log(2.); - const point_type pi_ov_npts_p1 = M_PI / (npts + 1); - - point_container points( npts ); - weight_container weights( npts ); - - /* - * Treutler-Ahlrichs quadrature - * - * Original reference: - * J. Chem. Phys. 102, 346 (1995) - * DOI: https://doi.org/10.1063/1.469408 - * - * Closed form for points / weights obtained from: - * Journal of Computational Chemistry, 24: 732–740, 2003 - * DOI: https://doi.org/10.1002/jcc.10211 - * - */ - for( size_t i = 0; i < npts; ++i ) { - const auto xi = std::cos( (i+1) * pi_ov_npts_p1 ); - - const auto pow_term = std::pow( 1. + xi, alpha ) / ln_2; - const auto log_term = std::log( (1. - xi)/2. ); - const auto sqrt_term = std::sqrt( (1.+xi)/(1.-xi) ); - - points[i] = -R * pow_term * log_term; - weights[i] = R * pi_ov_npts_p1 * pow_term * - ( sqrt_term - alpha * log_term / sqrt_term ); - } - - return std::make_tuple( points, weights ); - - } - -}; - - -#else - struct TreutlerAhlrichsM4RadialTraits { inline static constexpr double alpha = 0.6; inline static constexpr double a = 1.0; inline static constexpr double ln_2 = 0.693147180559945309417232; + /** + * @brief Transformation rule for the TA M4 radial quadrature + * + * Equation (19) of J. Chem. Phys. 102, 346 (1995) + * + * @param[in] x Point in (-1,1) + * @return r = (a+x)^alpha * log((a+1)/(1-x)) / ln(2) + */ template static auto radial_transform(PointType x) { const auto pow_term = std::pow(a + x, alpha); @@ -140,6 +39,12 @@ struct TreutlerAhlrichsM4RadialTraits { }; + /** + * @brief Jacobian of the TA M4 radial transformation + * + * @param[in] x Point in (-1,1) + * @returns dr/dx (see `radial_transform`) + */ template static auto radial_jacobian(PointType x) { const auto pow_term = std::pow(a + x, alpha); @@ -150,12 +55,28 @@ struct TreutlerAhlrichsM4RadialTraits { }; +/** + * @brief Implementation of the Treutler-Ahlrichs M4 radial quadrature. + * + * Taken as the convolution of the Gauss-Chebyshev (second kind) quadrature + * with the TA M4 radial transformation. See TreutlerAhlrichsM4RadialTraits for + * details. + * + * Suitable for integrands which tend to zero as their argument tends + * to 0 and inf. Tailored for radial integrands, i.e. r^2 * f(r), with + * lim_{r->inf} f(r) = 0. + * + * Reference: + * J. Chem. Phys. 102, 346 (1995) + * DOI: https://doi.org/10.1063/1.469408 + * + * @tparam PointType Type describing the quadrature points + * @tparam WeightType Type describing the quadrature weights + */ template using TreutlerAhlrichs = RadialTransformQuadrature< GaussChebyshev2, TreutlerAhlrichsM4RadialTraits >; -#endif - } From 7331e2c9436f31b0c24e40ccce0bb9effdba7bef Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 9 Jul 2023 10:57:35 -0700 Subject: [PATCH 10/15] Dox fixup --- include/integratorxx/quadratures/treutlerahlrichs.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/include/integratorxx/quadratures/treutlerahlrichs.hpp b/include/integratorxx/quadratures/treutlerahlrichs.hpp index fe621df..e81d5f5 100644 --- a/include/integratorxx/quadratures/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/treutlerahlrichs.hpp @@ -13,9 +13,6 @@ namespace IntegratorXX { * Reference: * J. Chem. Phys. 102, 346 (1995) * DOI: https://doi.org/10.1063/1.469408 - * - * @tparam M Integer to modulate the MHL transformation. - * Typically taken to be 2. */ struct TreutlerAhlrichsM4RadialTraits { From ffad94c7b88dba57ebf5f552938cf2c69fabfd6b Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 9 Jul 2023 11:27:05 -0700 Subject: [PATCH 11/15] RadialTraits are now runtime configurable --- include/integratorxx/quadratures/mhl.hpp | 18 ++++++++++----- .../integratorxx/quadratures/muraknowles.hpp | 16 +++++++++---- .../quadratures/radial_transform.hpp | 10 ++++---- .../quadratures/treutlerahlrichs.hpp | 23 ++++++++++++------- 4 files changed, 43 insertions(+), 24 deletions(-) diff --git a/include/integratorxx/quadratures/mhl.hpp b/include/integratorxx/quadratures/mhl.hpp index 02473b1..1359eb5 100644 --- a/include/integratorxx/quadratures/mhl.hpp +++ b/include/integratorxx/quadratures/mhl.hpp @@ -18,17 +18,23 @@ namespace IntegratorXX { * Typically taken to be 2. */ template -struct MurrayHandyLamingRadialTraits { +class MurrayHandyLamingRadialTraits { + + double R_; ///< Radial scaling factor + +public: + + MurrayHandyLamingRadialTraits(double R = 1.0) : R_(R) {} /** * @brief Transformation rule for the MHL radial quadrature * * @param[in] x Point in (0,1) - * @return r = (x / (1-x))^M + * @return r = R * (x / (1-x))^M */ template - static auto radial_transform(PointType x) { - return std::pow( x / (1.0 - x), M ); + inline auto radial_transform(PointType x) const noexcept { + return R_ * std::pow( x / (1.0 - x), M ); } /** @@ -38,8 +44,8 @@ struct MurrayHandyLamingRadialTraits { * @returns dr/dx (see `radial_transform`) */ template - static auto radial_jacobian(PointType x) { - return M * std::pow(x, M-1) / std::pow(1.0 - x, M+1); + inline auto radial_jacobian(PointType x) const noexcept { + return R_ * M * std::pow(x, M-1) / std::pow(1.0 - x, M+1); } }; diff --git a/include/integratorxx/quadratures/muraknowles.hpp b/include/integratorxx/quadratures/muraknowles.hpp index 06b9588..534070a 100644 --- a/include/integratorxx/quadratures/muraknowles.hpp +++ b/include/integratorxx/quadratures/muraknowles.hpp @@ -136,17 +136,23 @@ struct quadrature_traits< #else -struct MuraKnowlesRadialTraits { +class MuraKnowlesRadialTraits { + + double R_; ///< Radial scaling factor + +public: + + MuraKnowlesRadialTraits(double R = 1.0) : R_(R) { } template - static auto radial_transform(PointType x) { - return -std::log(1.0 - x*x*x); + inline auto radial_transform(PointType x) const noexcept { + return -R_ * std::log(1.0 - x*x*x); } template - static auto radial_jacobian(PointType x) { + inline auto radial_jacobian(PointType x) const noexcept { const auto x2 = x*x; - return 3.0 * x2 / (1.0 - x2 * x); + return R_ * 3.0 * x2 / (1.0 - x2 * x); } }; diff --git a/include/integratorxx/quadratures/radial_transform.hpp b/include/integratorxx/quadratures/radial_transform.hpp index 18183e6..93f8096 100644 --- a/include/integratorxx/quadratures/radial_transform.hpp +++ b/include/integratorxx/quadratures/radial_transform.hpp @@ -17,8 +17,8 @@ struct RadialTransformQuadrature : using point_container = typename base_type::point_container; using weight_container = typename base_type::weight_container; - RadialTransformQuadrature(size_t npts, weight_type R = 1.0) : - base_type( npts, R ) { } + RadialTransformQuadrature(size_t npts, const RadialTraits& traits = RadialTraits()) : + base_type( npts, traits ) { } RadialTransformQuadrature( const RadialTransformQuadrature& ) = default; RadialTransformQuadrature( RadialTransformQuadrature&& ) noexcept = default; @@ -34,7 +34,7 @@ struct quadrature_traits> { using weight_container = typename BaseQuad::weight_container; inline static std::tuple - generate( size_t npts, weight_type R ) { + generate( size_t npts, const RadialTraits& traits ) { using base_quad_traits = quadrature_traits; @@ -49,8 +49,8 @@ struct quadrature_traits> { for(size_t i = 0; i < npts; ++i) { const auto xi = base_x[i + ipts_offset]; const auto wi = base_w[i + ipts_offset]; - points[i] = R * RadialTraits::radial_transform(xi); - weights[i] = R * wi * RadialTraits::radial_jacobian(xi); + points[i] = traits.radial_transform(xi); + weights[i] = wi * traits.radial_jacobian(xi); } return std::make_tuple(points, weights); diff --git a/include/integratorxx/quadratures/treutlerahlrichs.hpp b/include/integratorxx/quadratures/treutlerahlrichs.hpp index e81d5f5..b10c47d 100644 --- a/include/integratorxx/quadratures/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/treutlerahlrichs.hpp @@ -14,12 +14,19 @@ namespace IntegratorXX { * J. Chem. Phys. 102, 346 (1995) * DOI: https://doi.org/10.1063/1.469408 */ -struct TreutlerAhlrichsM4RadialTraits { +class TreutlerAhlrichsM4RadialTraits { + + double R_; + double alpha_; + +public: - inline static constexpr double alpha = 0.6; inline static constexpr double a = 1.0; inline static constexpr double ln_2 = 0.693147180559945309417232; + TreutlerAhlrichsM4RadialTraits(double R = 1.0, double alpha = 0.6) : + R_(R), alpha_(alpha) { } + /** * @brief Transformation rule for the TA M4 radial quadrature * @@ -29,10 +36,10 @@ struct TreutlerAhlrichsM4RadialTraits { * @return r = (a+x)^alpha * log((a+1)/(1-x)) / ln(2) */ template - static auto radial_transform(PointType x) { - const auto pow_term = std::pow(a + x, alpha); + inline auto radial_transform(PointType x) const noexcept { + const auto pow_term = std::pow(a + x, alpha_); const auto log_term = std::log((a + 1.0) / (1.0 - x)); - return pow_term * log_term / ln_2; + return R_ * pow_term * log_term / ln_2; }; @@ -43,10 +50,10 @@ struct TreutlerAhlrichsM4RadialTraits { * @returns dr/dx (see `radial_transform`) */ template - static auto radial_jacobian(PointType x) { - const auto pow_term = std::pow(a + x, alpha); + inline auto radial_jacobian(PointType x) const noexcept { + const auto pow_term = std::pow(a + x, alpha_); const auto log_term = std::log((a + 1.0) / (1.0 - x)); - return pow_term / ln_2 * ( alpha * log_term / (a+x) + (1./(1. - x)) ); + return R_ * pow_term / ln_2 * ( alpha_ * log_term / (a+x) + (1./(1.-x)) ); } }; From 44e6c97c3ab24a169b534408643f68cf45a6c0f2 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 9 Jul 2023 12:26:47 -0700 Subject: [PATCH 12/15] Remove internal bounds transformations from GC quadratures --- include/integratorxx/quadratures/gausscheby1.hpp | 6 ++---- include/integratorxx/quadratures/gausscheby2.hpp | 11 ++--------- include/integratorxx/quadratures/gausscheby2_mod.hpp | 6 ++---- include/integratorxx/quadratures/gausscheby3.hpp | 7 +++---- test/1d_quadratures.cxx | 8 ++++---- 5 files changed, 13 insertions(+), 25 deletions(-) diff --git a/include/integratorxx/quadratures/gausscheby1.hpp b/include/integratorxx/quadratures/gausscheby1.hpp index 5912325..c57cee4 100644 --- a/include/integratorxx/quadratures/gausscheby1.hpp +++ b/include/integratorxx/quadratures/gausscheby1.hpp @@ -40,8 +40,7 @@ class GaussChebyshev1 using point_container = typename base_type::point_container; using weight_container = typename base_type::weight_container; - GaussChebyshev1(size_t npts, point_type lo, point_type up) - : base_type(npts, lo, up) {} + GaussChebyshev1(size_t npts) : base_type(npts) {} GaussChebyshev1(const GaussChebyshev1 &) = default; GaussChebyshev1(GaussChebyshev1 &&) noexcept = default; @@ -55,8 +54,7 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; - inline static std::tuple generate( - size_t npts, point_type lo, point_type up) { + inline static std::tuple generate(size_t npts) { point_container points(npts); weight_container weights(npts); diff --git a/include/integratorxx/quadratures/gausscheby2.hpp b/include/integratorxx/quadratures/gausscheby2.hpp index 7c414a4..5336eb5 100644 --- a/include/integratorxx/quadratures/gausscheby2.hpp +++ b/include/integratorxx/quadratures/gausscheby2.hpp @@ -42,8 +42,7 @@ class GaussChebyshev2 using point_container = typename base_type::point_container; using weight_container = typename base_type::weight_container; - GaussChebyshev2(size_t npts, point_type lo, point_type up) - : base_type(npts, lo, up) {} + GaussChebyshev2(size_t npts) : base_type(npts) {} GaussChebyshev2(const GaussChebyshev2 &) = default; GaussChebyshev2(GaussChebyshev2 &&) noexcept = default; @@ -59,8 +58,7 @@ struct quadrature_traits> { inline static constexpr bool bound_inclusive = false; - inline static std::tuple generate( - size_t npts, point_type lo, point_type up) { + inline static std::tuple generate(size_t npts) { const weight_type pi_ov_npts_p_1 = M_PI / (npts + 1); weight_container weights(npts); @@ -87,11 +85,6 @@ struct quadrature_traits> { return std::make_tuple(points, weights); } - inline static std::tuple - generate(size_t npts) { - return generate(npts, -1.0, 1.0); - } - }; } // namespace IntegratorXX diff --git a/include/integratorxx/quadratures/gausscheby2_mod.hpp b/include/integratorxx/quadratures/gausscheby2_mod.hpp index 6ebee05..818dc07 100644 --- a/include/integratorxx/quadratures/gausscheby2_mod.hpp +++ b/include/integratorxx/quadratures/gausscheby2_mod.hpp @@ -43,8 +43,7 @@ class GaussChebyshev2Modified using point_container = typename base_type::point_container; using weight_container = typename base_type::weight_container; - GaussChebyshev2Modified(size_t npts, point_type lo, point_type up) - : base_type(npts, lo, up) {} + GaussChebyshev2Modified(size_t npts) : base_type(npts) {} GaussChebyshev2Modified(const GaussChebyshev2Modified &) = default; GaussChebyshev2Modified(GaussChebyshev2Modified &&) noexcept = default; @@ -58,8 +57,7 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; - inline static std::tuple generate( - size_t npts, point_type lo, point_type up) { + inline static std::tuple generate(size_t npts) { const weight_type oonpp = 1.0 / (npts + 1); point_container points(npts); diff --git a/include/integratorxx/quadratures/gausscheby3.hpp b/include/integratorxx/quadratures/gausscheby3.hpp index 9b8f799..2e4684d 100644 --- a/include/integratorxx/quadratures/gausscheby3.hpp +++ b/include/integratorxx/quadratures/gausscheby3.hpp @@ -39,8 +39,7 @@ class GaussChebyshev3 using point_container = typename base_type::point_container; using weight_container = typename base_type::weight_container; - GaussChebyshev3(size_t npts, point_type lo, point_type up) - : base_type(npts, lo, up) {} + GaussChebyshev3(size_t npts) : base_type(npts) {} GaussChebyshev3(const GaussChebyshev3 &) = default; GaussChebyshev3(GaussChebyshev3 &&) noexcept = default; @@ -54,8 +53,7 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; - inline static std::tuple generate( - size_t npts, point_type lo, point_type up) { + inline static std::tuple generate(size_t npts) { const weight_type pi_ov_2n_p_1 = M_PI / (2 * npts + 1); weight_container weights(npts); @@ -77,6 +75,7 @@ struct quadrature_traits> { wi *= std::sqrt((1.0 - xi) / xi); // Finally, convert from [0,1] to [-1,1] + // TODO: This should be done externally points[idx] = 2 * xi - 1.0; weights[idx] = 2.0 * wi; } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 16c9a88..01d1fb9 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -180,19 +180,19 @@ TEST_CASE( "Gauss-Chebyshev Quadratures", "[1d-quad]") { }; SECTION("First Kind") { - IntegratorXX::GaussChebyshev1 quad(order, -1., 1.); + IntegratorXX::GaussChebyshev1 quad(order); integrate(quad); } SECTION("Second Kind") { - IntegratorXX::GaussChebyshev2 quad(order, -1., 1.); + IntegratorXX::GaussChebyshev2 quad(order); integrate(quad); } SECTION("Second Kind (Modified)") { - IntegratorXX::GaussChebyshev2Modified quad(order, -1., 1.); + IntegratorXX::GaussChebyshev2Modified quad(order); integrate(quad); } SECTION("Third Kind") { - IntegratorXX::GaussChebyshev3 quad(order, -1., 1.); + IntegratorXX::GaussChebyshev3 quad(order); integrate(quad); } } From ec2283a25437dc069bfcad9c10eacf460c8d3a67 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Mon, 10 Jul 2023 13:21:54 -0700 Subject: [PATCH 13/15] Update GC header conventions --- include/integratorxx/quadratures/treutlerahlrichs.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/integratorxx/quadratures/treutlerahlrichs.hpp b/include/integratorxx/quadratures/treutlerahlrichs.hpp index b10c47d..d2f3abd 100644 --- a/include/integratorxx/quadratures/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/treutlerahlrichs.hpp @@ -1,7 +1,7 @@ #pragma once #include -#include +#include #include namespace IntegratorXX { From 8caf841f685dbd3630333a28d34d5bad2a659799 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Tue, 11 Jul 2023 16:12:45 -0700 Subject: [PATCH 14/15] Add reference MK data from Molpro (HT @pjknowles) --- test/CMakeLists.txt | 2 ++ test/muraknowles.cxx | 65 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 67 insertions(+) create mode 100644 test/muraknowles.cxx diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 543719a..fdb1566 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -30,6 +30,8 @@ target_link_libraries( gausslegendre PUBLIC Catch2::Catch2WithMain integratorxx add_executable( gausslobatto gausslobatto.cxx ) target_link_libraries( gausslobatto PUBLIC Catch2::Catch2WithMain integratorxx ) +add_executable( muraknowles muraknowles.cxx ) +target_link_libraries( muraknowles PUBLIC Catch2::Catch2WithMain integratorxx ) add_executable( composite_quadratures composite_quadratures.cxx ) target_link_libraries( composite_quadratures PUBLIC Catch2::Catch2WithMain integratorxx ) diff --git a/test/muraknowles.cxx b/test/muraknowles.cxx new file mode 100644 index 0000000..198eb4a --- /dev/null +++ b/test/muraknowles.cxx @@ -0,0 +1,65 @@ +// clang-format off +#include "catch2/catch_all.hpp" +#include +#include +#include +// clang-format on + +const double x_tolerance = std::numeric_limits::epsilon(); +const double w_tolerance = std::numeric_limits::epsilon(); + +// Reference data provided by Peter Knowles +// https://github.com/wavefunction91/IntegratorXX/pull/27#discussion_r1260345076 +// ALPHA = 7.0 +// M = 3 +// NPTS = 48 +TEST_CASE("48 point MuraKnowles", "[1d-quad]") { + std::array ref_pts = { + 5.9499271133939913E-005, 4.7600833032953915E-004, 1.6066578611874361E-003, + 3.8089732866334377E-003, 7.4413311256047504E-003, 1.2863600152286464E-002, + 2.0437970723343253E-002, 3.0529977777593031E-002, 4.3509725746882791E-002, + 5.9753327005997357E-002, 7.9644569620570760E-002, 0.10357683514824530, + 0.13195529327830460, 0.16519940738611444, 0.20374579393428979, + 0.24805148947543415, 0.29859769233640943, 0.35589406260236178, + 0.42048368472012887, 0.49294882318131050, 0.57391763505319449, + 0.66407204595981006, 0.76415705172765447, 0.87499178083047879, + 0.99748274937600856, 1.1326398697743505, 1.2815959495630309, + 1.4456306574120197, 1.6262002677869758, 1.8249749675453268, + 2.0438861838399740, 2.2851873784253796, 2.5515332181111825, + 2.8460842535535971, 3.1726476906379575, 3.5358703430273275, + 3.9415088946346315, 4.3968179569776069, 4.9111235220893743, + 5.4966995004301378, 6.1701626127089160, 6.9548035249514397, + 7.8847270471719106, 9.0128012057250508, 10.427589422901853, + 12.295111781246881, 14.988083164786680, 19.695799083317340}; + + std::array ref_wgt = { + 6.3191408757371128E-013, 1.6178925855460011E-010, 4.1478214302682156E-009, + 4.1457509371301411E-008, 2.4736257856667049E-007, 1.0652626596586621E-006, + 3.6641284849845252E-006, 1.0694442455644551E-005, 2.7541538616163128E-005, + 6.4278015117145217E-005, 1.3857027486190222E-004, 2.7986264953446042E-004, + 5.3525184281221503E-004, 9.7758247883188328E-004, 1.7164520965350302E-003, + 2.9130257653514745E-003, 4.7998471771280651E-003, 7.7072235182964691E-003, + 1.2098296752658244E-002, 1.8615653176199027E-002, 2.8143351896444287E-002, + 4.1889696901313278E-002, 6.1498123098954999E-002, 8.9196494268593482E-002, + 0.12799934660813353, 0.18198381348671958, 0.25666916599665018, + 0.35954374153109814, 0.50080416940484207, 0.69440463378479145, + 0.95956585458830923, 1.3229773146815327, 1.8220646710507786, + 2.5099284939617785, 3.4629679021480624, 4.7929337979063913, + 6.6665154642672082, 9.3381941641618038, 13.207427859666060, + 18.922636048645305, 27.580468647622130, 41.132941288449700, + 63.288684869795212, 101.72526373304323, 174.33040888875871, + 330.69364765060254, 753.72105349957747, 2659.6767864598341}; + + IntegratorXX::MuraKnowlesRadialTraits traits{7.0}; + IntegratorXX::MuraKnowles quad(48, traits); + + const auto& pts = quad.points(); + const auto& wgt = quad.weights(); + for(auto i = 0ul; i < 2; i++) { + REQUIRE_THAT(pts[i], Catch::Matchers::WithinAbs(ref_pts[i], x_tolerance)); + // Reference data contains the spherical Jacobian (r^2) + const auto rsq = pts[i] * pts[i]; + REQUIRE_THAT(rsq * wgt[i], + Catch::Matchers::WithinAbs(ref_wgt[i], w_tolerance)); + } +} From 5842dbb0ef8ab2ad08bd74d756b4633c65e64bb4 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Tue, 11 Jul 2023 16:36:32 -0700 Subject: [PATCH 15/15] Rename TreutlerAhlrichs{M4,}RadialTraits, add dox for M3 code path --- .../quadratures/treutlerahlrichs.hpp | 29 ++++++++++++------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/include/integratorxx/quadratures/treutlerahlrichs.hpp b/include/integratorxx/quadratures/treutlerahlrichs.hpp index d2f3abd..5759d36 100644 --- a/include/integratorxx/quadratures/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/treutlerahlrichs.hpp @@ -7,14 +7,14 @@ namespace IntegratorXX { /** - * @brief Implementation of the Treutler-Ahlrichs M4 radial + * @brief Implementation of the Treutler-Ahlrichs M3+M4 radial * quadrature transformation rules. * * Reference: * J. Chem. Phys. 102, 346 (1995) * DOI: https://doi.org/10.1063/1.469408 */ -class TreutlerAhlrichsM4RadialTraits { +class TreutlerAhlrichsRadialTraits { double R_; double alpha_; @@ -24,13 +24,22 @@ class TreutlerAhlrichsM4RadialTraits { inline static constexpr double a = 1.0; inline static constexpr double ln_2 = 0.693147180559945309417232; - TreutlerAhlrichsM4RadialTraits(double R = 1.0, double alpha = 0.6) : + /** + * Specify Treutler-Ahlrichs quadrature parameters + * + * M3: Equation (18) of J. Chem. Phys. 102, 346 (1995) + * M4: Equation (19) of J. Chem. Phys. 102, 346 (1995) + * + * Default to M4 (alpha = 0.6). M3 resolved with alpha = 0.0 + * + * @param[in] R Radial scaling factor + * @param[in] alpha TA exponential factor + */ + TreutlerAhlrichsRadialTraits(double R = 1.0, double alpha = 0.6) : R_(R), alpha_(alpha) { } /** - * @brief Transformation rule for the TA M4 radial quadrature - * - * Equation (19) of J. Chem. Phys. 102, 346 (1995) + * @brief Transformation rule for the TA M3+M4 radial quadratures * * @param[in] x Point in (-1,1) * @return r = (a+x)^alpha * log((a+1)/(1-x)) / ln(2) @@ -44,7 +53,7 @@ class TreutlerAhlrichsM4RadialTraits { /** - * @brief Jacobian of the TA M4 radial transformation + * @brief Jacobian of the TA M3+M4 radial transformations * * @param[in] x Point in (-1,1) * @returns dr/dx (see `radial_transform`) @@ -63,11 +72,11 @@ class TreutlerAhlrichsM4RadialTraits { * @brief Implementation of the Treutler-Ahlrichs M4 radial quadrature. * * Taken as the convolution of the Gauss-Chebyshev (second kind) quadrature - * with the TA M4 radial transformation. See TreutlerAhlrichsM4RadialTraits for + * with the TA M4 radial transformation. See TreutlerAhlrichsRadialTraits for * details. * * Suitable for integrands which tend to zero as their argument tends - * to 0 and inf. Tailored for radial integrands, i.e. r^2 * f(r), with + f(r), with * lim_{r->inf} f(r) = 0. * * Reference: @@ -80,7 +89,7 @@ class TreutlerAhlrichsM4RadialTraits { template using TreutlerAhlrichs = RadialTransformQuadrature< GaussChebyshev2, - TreutlerAhlrichsM4RadialTraits + TreutlerAhlrichsRadialTraits >; }