Lindh-Malmqvist-Gagliardi Radial Grid

include/integratorxx/quadratures/radial.hpp

@@ -5,3 +5,4 @@
 #include <integratorxx/quadratures/radial/mhl.hpp>
 #include <integratorxx/quadratures/radial/muraknowles.hpp>
 #include <integratorxx/quadratures/radial/treutlerahlrichs.hpp>
+#include <integratorxx/quadratures/radial/lmg.hpp>

include/integratorxx/quadratures/radial/lmg.hpp

@@ -0,0 +1,127 @@
+#include <integratorxx/util/lambert_w.hpp>
+#include <integratorxx/quadratures/radial/radial_transform.hpp>
+
+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);
+  const double x = alpha * r * r;
+  return g_term * std::pow(x, am_term) * std::exp(-x);
+}
+
+// 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
+  // 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;
+}
+
+// Solve Eq 25 of LMG paper (DOI 10.1007/s002140100263)
+inline double r_lower(int m, double alpha, double prec) {
+  // Magic numbers  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:
+      // 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 * 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));
+  return M_PI / X;
+
+}
+
+}
+
+
+class LindhMalmqvistGagliardiRadialTraits {
+  using self_type = LindhMalmqvistGagliardiRadialTraits;
+
+  double step_size_;
+  double c_;
+
+public:
+
+  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 <typename PointType>
+  inline auto radial_transform(PointType x) const noexcept {
+    return c_ * (std::exp(x) - 1.0);
+  }
+
+  template <typename PointType>
+  inline auto radial_jacobian(PointType x) const noexcept {
+    return c_ * std::exp(x);
+  }
+
+  template <typename BaseQuad>
+  std::enable_if_t<
+    // LMG really only makes sense with a Trapezoid grid
+    std::is_same_v<
+      BaseQuad,
+      UniformTrapezoid<typename BaseQuad::point_type,
+                       typename BaseQuad::weight_type>
+    >
+  > 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 <typename PointType, typename WeightType>
+using LindhMalmqvistGagliardi = RadialTransformQuadrature<
+  UniformTrapezoid<PointType, WeightType>, LindhMalmqvistGagliardiRadialTraits
+>;
+
+}

include/integratorxx/quadratures/radial/radial_transform.hpp

@@ -4,6 +4,32 @@
 
 namespace IntegratorXX {
 
+namespace detail {
+
+template <typename TraitsType, typename BaseQuad>
+class has_preprocess_base_quad {
+
+using point_container_ref  = typename BaseQuad::point_container&;
+using weight_container_ref = typename BaseQuad::weight_container&;
+
+template <typename C, typename = 
+  decltype(std::declval<C>().
+    template preprocess_base_quad<BaseQuad>(
+      std::declval<point_container_ref>(), 
+      std::declval<weight_container_ref>()
+    ))>
+static std::true_type test(int);
+
+template<typename C> static std::false_type test(...);
+
+public:
+
+static constexpr bool value = decltype(test<TraitsType>(0))::value;
+
+};
+
+}
+
 template <typename BaseQuad, typename RadialTraits>
 struct RadialTransformQuadrature :
   public Quadrature<RadialTransformQuadrature<BaseQuad, RadialTraits>> {
@@ -45,6 +71,9 @@ struct quadrature_traits<RadialTransformQuadrature<BaseQuad, RadialTraits>> {
     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<RadialTraits,BaseQuad>::value)
+      traits.template preprocess_base_quad<BaseQuad>(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];

include/integratorxx/util/lambert_w.hpp

@@ -0,0 +1,68 @@
+#pragma once
+#include <cmath>
+#include <iostream>
+
+namespace IntegratorXX {
+
+// Newton-Raphson implementation of Lambert-W function for real arguments
+template <typename RealType>
+RealType lambert_w_newton(RealType x, RealType w0) {
+
+  const auto eps = std::numeric_limits<RealType>::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<RealType>::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);
+
+}
+
+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<RealType>::quiet_NaN();
+  if(x >= RealType(0)) return std::numeric_limits<RealType>::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);
+
+}
+
+}