Modularize Radial Quadrature Generation

include/integratorxx/quadratures/gausschebyshev1.hpp

@@ -39,8 +39,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;
@@ -54,8 +53,7 @@ struct quadrature_traits<GaussChebyshev1<PointType, WeightType>> {
   using point_container = std::vector<point_type>;
   using weight_container = std::vector<weight_type>;
 
-  inline static std::tuple<point_container, weight_container> generate(
-      size_t npts, point_type lo, point_type up) {
+  inline static std::tuple<point_container, weight_container> generate(size_t npts) {
     point_container points(npts);
     weight_container weights(npts);
 

include/integratorxx/quadratures/gausschebyshev2.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;
@@ -56,9 +55,10 @@ struct quadrature_traits<GaussChebyshev2<PointType, WeightType>> {
 
   using point_container = std::vector<point_type>;
   using weight_container = std::vector<weight_type>;
+
+  inline static constexpr bool bound_inclusive = false;
 
-  inline static std::tuple<point_container, weight_container> generate(
-      size_t npts, point_type lo, point_type up) {
+  inline static std::tuple<point_container, weight_container> generate(size_t npts) {
     const weight_type pi_ov_npts_p_1 = M_PI / (npts + 1);
 
     weight_container weights(npts);
@@ -96,6 +96,7 @@ struct quadrature_traits<GaussChebyshev2<PointType, WeightType>> {
 
     return std::make_tuple(points, weights);
   }
+
 };
 
 }  // namespace IntegratorXX

include/integratorxx/quadratures/gausschebyshev2modified.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<GaussChebyshev2Modified<PointType, WeightType>> {
   using point_container = std::vector<point_type>;
   using weight_container = std::vector<weight_type>;
 
-  inline static std::tuple<point_container, weight_container> generate(
-      size_t npts, point_type lo, point_type up) {
+  inline static std::tuple<point_container, weight_container> generate(size_t npts) {
     const weight_type oonpp = 1.0 / (npts + 1);
 
     point_container points(npts);

include/integratorxx/quadratures/gausschebyshev3.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<GaussChebyshev3<PointType, WeightType>> {
   using point_container = std::vector<point_type>;
   using weight_container = std::vector<weight_type>;
 
-  inline static std::tuple<point_container, weight_container> generate(
-      size_t npts, point_type lo, point_type up) {
+  inline static std::tuple<point_container, weight_container> 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<GaussChebyshev3<PointType, WeightType>> {
       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;
     }

include/integratorxx/quadratures/mhl.hpp

@@ -2,142 +2,78 @@
 
 #include <integratorxx/quadrature.hpp>
 #include <integratorxx/quadratures/uniform.hpp>
+#include <integratorxx/quadratures/radial_transform.hpp>
 
 namespace IntegratorXX {
 
 /**
- *  @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 
+ *  @tparam M Integer to modulate the MHL transformation. 
+ *            Typically taken to be 2.
  */
-template <typename PointType, typename WeightType>
-class MurrayHandyLaming : 
-  public Quadrature<MurrayHandyLaming<PointType,WeightType>> {
+template <size_t M>
+class MurrayHandyLamingRadialTraits {
 
-  using base_type = Quadrature<MurrayHandyLaming<PointType,WeightType>>;
+  double R_; ///< Radial scaling factor
 
 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;
-
+  MurrayHandyLamingRadialTraits(double R = 1.0) : R_(R) {}
+
   /**
-   *  @brief Construct the Murray-Handy-Laming radial quadrature
+   *  @brief Transformation rule for the MHL radial quadrature
+   *  
+   *  @param[in] x Point in (0,1)
+   *  @return    r = R * (x / (1-x))^M
+   */
+  template <typename PointType>
+  inline auto radial_transform(PointType x) const noexcept {
+    return R_ * std::pow( x / (1.0 - x), M );
+  }
+
+  /**
+   *  @brief Jacobian of the MHL radial transformation
    *
-   *  @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. 
+   *  @param[in] x Point in (0,1)
+   *  @returns   dr/dx (see `radial_transform`)
    */
-  MurrayHandyLaming(size_t npts, weight_type R = 1., int m = 2): 
-    base_type( npts, R, m ) { }
+  template <typename PointType>
+  inline auto radial_jacobian(PointType x) const noexcept {
+    return R_ * M * std::pow(x, M-1) / std::pow(1.0 - x, M+1);
+  }
 
-  MurrayHandyLaming( const MurrayHandyLaming& )     = default;
-  MurrayHandyLaming( MurrayHandyLaming&& ) noexcept = default;
 };
 
 
-
-
-
-
 /**
- *  @brief Quadrature traits for the Murray-Handy-Laming quadrature
+ *  @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 <typename PointType, typename WeightType>
-struct quadrature_traits<
-  MurrayHandyLaming<PointType,WeightType>,
-  std::enable_if_t<
-    std::is_floating_point_v<PointType> and
-    std::is_floating_point_v<WeightType>
-  >
-> {
-
-  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<point_container,weight_container>
-    generate( size_t npts, weight_type R, int m ) {
-
-    point_container  points( npts );
-    weight_container weights( npts );
-
-
-    using base_quad_traits = 
-      quadrature_traits<UniformTrapezoid<PointType,WeightType>>;
-
-    /*
-     * 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 );
-
-  }
-
-};
+using MurrayHandyLaming = RadialTransformQuadrature<
+  UniformTrapezoid<PointType,WeightType>,
+  MurrayHandyLamingRadialTraits<2>
+>;
 
 }
+

include/integratorxx/quadratures/muraknowles.hpp

@@ -2,9 +2,11 @@
 
 #include <integratorxx/quadrature.hpp>
 #include <integratorxx/quadratures/uniform.hpp>
+#include <integratorxx/quadratures/radial_transform.hpp>
 
 namespace IntegratorXX {
 
+#if 0
 /**
  *  @brief Implementation of the Mura-Knowles radial quadrature.
  *
@@ -132,4 +134,35 @@ struct quadrature_traits<
 
 };
 
+#else
+
+class MuraKnowlesRadialTraits {
+
+  double R_; ///< Radial scaling factor
+
+public:
+
+  MuraKnowlesRadialTraits(double R = 1.0) : R_(R) { }
+
+  template <typename PointType>
+  inline auto radial_transform(PointType x) const noexcept {
+    return -R_ * std::log(1.0 - x*x*x);
+  }
+
+  template <typename PointType>
+  inline auto radial_jacobian(PointType x) const noexcept {
+    const auto x2 = x*x;
+    return R_ * 3.0 * x2 / (1.0 - x2 * x);
+  }
+
+}; 
+
+template <typename PointType, typename WeightType>
+using MuraKnowles = RadialTransformQuadrature<
+  UniformTrapezoid<PointType,WeightType>,
+  MuraKnowlesRadialTraits
+>;
+
+#endif
+
 }

include/integratorxx/quadratures/radial_transform.hpp

@@ -0,0 +1,63 @@
+#pragma once
+
+#include <integratorxx/quadrature.hpp>
+
+namespace IntegratorXX {
+
+template <typename BaseQuad, typename RadialTraits>
+struct RadialTransformQuadrature :
+  public Quadrature<RadialTransformQuadrature<BaseQuad, RadialTraits>> {
+
+  using base_type = Quadrature<RadialTransformQuadrature<BaseQuad, RadialTraits>>;
+
+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, const RadialTraits& traits = RadialTraits()) :
+    base_type( npts, traits ) { }
+
+  RadialTransformQuadrature( const RadialTransformQuadrature& )     = default;
+  RadialTransformQuadrature( RadialTransformQuadrature&& ) noexcept = default;
+
+};
+
+template <typename BaseQuad, typename RadialTraits>
+struct quadrature_traits<RadialTransformQuadrature<BaseQuad, RadialTraits>> {
+
+  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<point_container,weight_container>
+    generate( size_t npts, const RadialTraits& traits ) {
+
+    using base_quad_traits = quadrature_traits<BaseQuad>;
+
+    point_container  points( npts );
+    weight_container weights( npts );
+
+
+    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 + ipts_offset];
+      const auto wi = base_w[i + ipts_offset];
+      points[i]  = traits.radial_transform(xi);
+      weights[i] = wi * traits.radial_jacobian(xi);
+    }
+
+    return std::make_tuple(points, weights);
+  }
+
+};
+
+
+
+}