From 62a8191d26e6b3488b91612b90385d01a086f8d9 Mon Sep 17 00:00:00 2001
From: Susi Lehtola <susi.lehtola@gmail.com>
Date: Tue, 3 Oct 2023 11:39:08 +0300
Subject: [PATCH 1/2] First steps towards the implementation of the Baker et al
 radial grid.

---
 .../integratorxx/quadratures/radial/baker.hpp | 79 +++++++++++++++++++
 1 file changed, 79 insertions(+)
 create mode 100644 include/integratorxx/quadratures/radial/baker.hpp

diff --git a/include/integratorxx/quadratures/radial/baker.hpp b/include/integratorxx/quadratures/radial/baker.hpp
new file mode 100644
index 0000000..af062ec
--- /dev/null
+++ b/include/integratorxx/quadratures/radial/baker.hpp
@@ -0,0 +1,79 @@
+#pragma once
+
+#include <integratorxx/quadratures/primitive/uniform.hpp>
+#include <integratorxx/quadratures/radial/radial_transform.hpp>
+
+namespace IntegratorXX {
+
+/**
+ *  @brief Implementation of the Baker-Andzelm-Scheiner-Delley radial
+ *  quadrature transformation rules.
+ *
+ *  Reference:
+ *  J. Chem. Phys. 101, 8894 (1994)
+ *  DOI: https://doi.org/10.1063/1.468081
+ */
+class BakerRadialTraits {
+
+  double Rmax_;
+
+public:
+
+
+/**
+ *  Specify Baker-Andzelm-Scheiner-Delley quadrature parameters.
+ *
+ *  @param[in] Rmax  Maximum radial extent, constant 12 in Equation (12b) of J. Chem. Phys. 101, 8894 (1994).
+ */
+  BakerRadialTraits(double Rmax = 12.0) :
+    Rmax_(Rmax) { }
+
+  /**
+   *  @brief Transformation rule for the Baker-Andzelm-Scheiner-Delley radial quadratures
+   *  
+   *  @param[in] x Point in (-1,1)
+   *  @return    r = Rfac * log [1 - x_i^2 ]
+   */
+  template <typename PointType>
+  inline auto radial_transform(PointType x) const noexcept {
+    const auto log_term = std::log(1.0 - x*x);
+    return Rfac_ * log_term * log_term;
+  };
+
+
+  /**
+   *  @brief Jacobian of the Baker-Andzelm-Scheiner-Delley radial transformations
+   *
+   *  @param[in] x Point in (-1,1)
+   *  @returns   dr/dx (see `radial_transform`)
+   */
+  template <typename PointType>
+  inline auto radial_jacobian(PointType x) const noexcept {
+    const auto log_term = std::log(1.0 - x*x);
+    return - 2.0 * Rfac_ * x / (1.0 - x*x);
+  }
+
+};
+
+
+/**
+ *  @brief Implementation of the Baker-Andzelm-Scheiner-Delley radial quadrature.
+ *
+ *  Taken as the convolution of the uniform trapezoidal quadrature
+ *  with the Baker-Andzelm-Scheiner-Delley radial transformation. See
+ *  BakerRadialTraits for details.
+ *
+ *  Reference:
+ *  J. Chem. Phys. 101, 8894 (1994)
+ *  DOI: https://doi.org/10.1063/1.468081
+ *
+ *  @tparam PointType  Type describing the quadrature points  
+ *  @tparam WeightType Type describing the quadrature weights 
+ */
+template <typename PointType, typename WeightType>
+using Baker = RadialTransformQuadrature<
+  UniformTrapezoid<PointType, WeightType>,
+  BakerRadialTraits 
+>;
+
+}

From c420c93292112d92abf80d29d3b1f792675cab85 Mon Sep 17 00:00:00 2001
From: Susi Lehtola <susi.lehtola@gmail.com>
Date: Tue, 3 Oct 2023 12:05:28 +0300
Subject: [PATCH 2/2] Erroneous minus sign

---
 include/integratorxx/quadratures/radial/baker.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/integratorxx/quadratures/radial/baker.hpp b/include/integratorxx/quadratures/radial/baker.hpp
index af062ec..dce0ae8 100644
--- a/include/integratorxx/quadratures/radial/baker.hpp
+++ b/include/integratorxx/quadratures/radial/baker.hpp
@@ -50,7 +50,7 @@ class BakerRadialTraits {
   template <typename PointType>
   inline auto radial_jacobian(PointType x) const noexcept {
     const auto log_term = std::log(1.0 - x*x);
-    return - 2.0 * Rfac_ * x / (1.0 - x*x);
+    return 2.0 * Rfac_ * x / (1.0 - x*x);
   }
 
 };