Radial quadrature grid (#5173)

* initial commit

* ta3

* add unit test

* more tests

* minor change

Author: jinzx10

source/module_base/grid/delley.cpp

@@ -1,20 +1,21 @@
 #include "module_base/grid/delley.h"
 
+#include <functional>
 #include <algorithm>
 #include <cmath>
 #include <cassert>
 
 namespace {
 
-struct Table {
+struct DelleyTable {
     const int lmax_;
     const int ngrid_;
     const int ntype_[6];
     const std::vector<double> data_;
 };
 
 // Delley's table from the original article
-const std::vector<Table> table = {
+const std::vector<DelleyTable> delley_table = {
     {
         17, 110, {1, 1, 0, 3, 1, 0},
         {
@@ -209,7 +210,7 @@ const std::vector<Table> table = {
             0.63942796347491023, 0.06424549224220589, 0.0009158016174693465,
         }
     }
-}; // end of the definition of "table"
+}; // end of the definition of "delley_table"
 
 // size of each group of points with octahedral symmetry
 // 6: (1, 0, 0) x sign x permutation (vertices)
@@ -324,15 +325,15 @@ const std::vector<Fill_t> fill = {
     },
 }; // end of the definition of "fill"
 
-const Table* _find_table(int& lmax) {
-    // NOTE: this function assumes elements in "Delley::table_" are
+const DelleyTable* _find_delley(int& lmax) {
+    // NOTE: this function assumes elements in "delley_table" are
     // arranged such that their members "lmax_" are in ascending order.
-    auto tab = std::find_if(table.begin(), table.end(),
-            [lmax](const Table& t) { return t.lmax_ >= lmax; });
-    return tab == table.end() ? nullptr : (lmax = tab->lmax_, &(*tab));
+    auto tab = std::find_if(delley_table.begin(), delley_table.end(),
+            [lmax](const DelleyTable& t) { return t.lmax_ >= lmax; });
+    return tab == delley_table.end() ? nullptr : (lmax = tab->lmax_, &(*tab));
 }
 
-void _get(const Table* tab, double* grid, double* weight) {
+void _delley(const DelleyTable* tab, double* grid, double* weight) {
     assert(tab);
     const double* ptr = &tab->data_[0];
     for (int itype = 0; itype < 6; ++itype) {
@@ -348,29 +349,29 @@ void _get(const Table* tab, double* grid, double* weight) {
 } // end of anonymous namespace
 
 
-int Grid::Delley::ngrid(int& lmax) {
-    auto tab = _find_table(lmax);
+int Grid::Angular::ngrid_delley(int& lmax) {
+    auto tab = _find_delley(lmax);
     return tab ? tab->ngrid_ : -1;
 }
 
 
-int Grid::Delley::get(int& lmax, double* grid, double* weight) {
-    auto tab = _find_table(lmax);
-    return tab ? _get(tab, grid, weight), 0 : -1;
+int Grid::Angular::delley(int& lmax, double* grid, double* weight) {
+    auto tab = _find_delley(lmax);
+    return tab ? _delley(tab, grid, weight), 0 : -1;
 }
 
 
-int Grid::Delley::get(
+int Grid::Angular::delley(
     int& lmax,
     std::vector<double>& grid,
     std::vector<double>& weight
 ) {
-    auto tab = _find_table(lmax);
+    auto tab = _find_delley(lmax);
     if (!tab) {
         return -1;
     }
     grid.resize(3 * tab->ngrid_);
     weight.resize(tab->ngrid_);
-    _get(tab, grid.data(), weight.data());
+    _delley(tab, grid.data(), weight.data());
     return 0;
 }

source/module_base/grid/delley.h

@@ -1,19 +1,10 @@
-#ifndef GRID_DELLEY_H
-#define GRID_DELLEY_H
+#ifndef GRID_ANGULAR_DELLEY_H
+#define GRID_ANGULAR_DELLEY_H
 
 #include <vector>
-#include <functional>
 
-/**
- * @brief Delley's grid for quadrature on the unit sphere.
- *
- * Reference:
- * Delley, B. (1996). High order integration schemes on the unit sphere.
- * Journal of computational chemistry, 17(9), 1152-1155.
- *
- */
 namespace Grid {
-namespace Delley {
+namespace Angular {
 
 /**
  * @brief Number of Delley's grid points for a certain order of accuracy.
@@ -27,7 +18,7 @@ namespace Delley {
  * lmax will be set to 23.
  *
  */
-int ngrid(int& lmax);
+int ngrid_delley(int& lmax);
 
 
 /**
@@ -45,13 +36,16 @@ int ngrid(int& lmax);
  * ngrid(lmax) elements, respectively. The function will return 0 if
  * successful, or -1 if the requested order cannot be fulfilled.
  *
+ * Reference:
+ * Delley, B. (1996). High order integration schemes on the unit sphere.
+ * Journal of computational chemistry, 17(9), 1152-1155.
  */
-int get(int& lmax, double* grid, double* weight);
+int delley(int& lmax, double* grid, double* weight);
 
 // a handy wrapper doing the same as above
-int get(int& lmax, std::vector<double>& grid, std::vector<double>& weight);
+int delley(int& lmax, std::vector<double>& grid, std::vector<double>& weight);
 
-} // end of namespace Delley
+} // end of namespace Angular
 } // end of namespace Grid
 
 #endif

source/module_base/grid/radial.cpp

@@ -0,0 +1,69 @@
+#include "module_base/grid/radial.h"
+
+#include <cmath>
+
+namespace {
+
+const double pi = std::acos(-1.0);
+const double inv_ln2 = 1.0 / std::log(2.0);
+
+} // end of anonymous namespace
+
+
+namespace Grid {
+namespace Radial {
+
+void baker(int nbase, double R, double* r, double* w, int mult) {
+    int n = (nbase+1) * mult - 1;
+    double r0 = -R / std::log((2.0 * nbase + 1.0) / ((nbase+1)*(nbase+1)));
+    for (int i = 1; i <= n; ++i) {
+        r[i-1] = -r0 * std::log(1.0 - static_cast<double>(i)*i/((n+1)*(n+1)));
+        w[i-1] = 2.0 * i * r0 * r[i-1] * r[i-1] / ((n+1+i)*(n+1-i));
+    }
+}
+
+
+void baker(int nbase, double R, std::vector<double>& r,
+           std::vector<double>& w, int mult) {
+    int n = (nbase+1) * mult - 1;
+    r.resize(n);
+    w.resize(n);
+    baker(nbase, R, r.data(), w.data(), mult);
+}
+
+
+void murray(int n, double R, double* r, double* w) {
+    for (int i = 1; i <= n; ++i) {
+        double x = static_cast<double>(i) / (n + 1);
+        r[i-1] = std::pow(x / (1.0 - x), 2) * R;
+        w[i-1] = 2.0 / (n + 1) * std::pow(R, 3) * std::pow(x, 5)
+                 / std::pow(1.0 - x, 7);
+    }
+}
+
+
+void treutler_m4(int n, double R, double* r, double* w, double alpha) {
+    for (int i = 1; i <= n; ++i) {
+        double x = std::cos(i * pi / (n + 1));
+        double beta = std::sqrt((1.0 + x) / (1.0 - x));
+        double gamma = std::log((1.0 - x) / 2.0);
+        double delta = std::pow(1.0 + x, alpha);
+        r[i-1] = -R * inv_ln2 * delta * gamma;
+        w[i-1] = pi / (n + 1) * std::pow(delta * R * inv_ln2, 3)
+                 * gamma * gamma * (beta - alpha / beta * gamma);
+    }
+}
+
+
+void mura(int n, double R, double* r, double* w) {
+    for (int i = 1; i <= n; ++i) {
+        double x = static_cast<double>(i) / (n + 1);
+        double alpha = 1.0 - x * x * x;
+        r[i-1] = -R * std::log(alpha);
+        w[i-1] = 3.0 * R * std::pow(x * r[i-1], 2) / ((n+1) * alpha);
+    }
+}
+
+
+} // end of namespace Radial
+} // end of namespace Grid

source/module_base/grid/radial.h

@@ -0,0 +1,70 @@
+#ifndef GRID_RADIAL_H
+#define GRID_RADIAL_H
+
+#include <vector>
+
+namespace Grid {
+namespace Radial {
+
+/**
+ * @brief Radial quadratures.
+ *
+ * This namespace contains functions that generate grids and weights
+ * for numerical integration
+ *
+ *      / inf     2
+ *      |     dr r  g(r) ~ \sum_i w[i] g(r[i])
+ *      /  0
+ *
+ */
+
+/**
+ * Baker, J., Andzelm, J., Scheiner, A., & Delley, B. (1994).
+ * The effect of grid quality and weight derivatives in
+ * density functional calculations.
+ * The Journal of chemical physics, 101(10), 8894-8902.
+ *
+ * Zhang, I. Y., Ren, X., Rinke, P., Blum, V., & Scheffler, M. (2013).
+ * Numeric atom-centered-orbital basis sets with valence-correlation
+ * consistency from H to Ar.
+ * New Journal of Physics, 15(12), 123033.
+ *
+ * @note nbase is the number of points of the "base" grid, i.e.,
+ * before applying the "radial multiplier" introduced by Zhang et al.
+ * The true number of grid points is (nbase+1) * mult - 1.
+ */
+void baker(int nbase, double R, double* r, double* w, int mult = 1);
+void baker(int nbase, double R, std::vector<double>& r,
+           std::vector<double>& w, int mult = 1);
+
+
+/**
+ * Murray, C. W., Handy, N. C., & Laming, G. J. (1993).
+ * Quadrature schemes for integrals of density functional theory.
+ * Molecular Physics, 78(4), 997-1014.
+ */
+void murray(int n, double R, double* r, double* w);
+
+
+/**
+ * Treutler, O., & Ahlrichs, R. (1995).
+ * Efficient molecular numerical integration schemes.
+ * The Journal of Chemical Physics, 102(1), 346-354.
+ *
+ * @note M4 reduces to M3 at alpha = 0.
+ */
+void treutler_m4(int n, double R, double* r, double* w, double alpha = 0.6);
+
+
+/**
+ * Mura, M. E., & Knowles, P. J. (1996).
+ * Improved radial grids for quadrature in molecular
+ * density‐functional calculations.
+ * The Journal of chemical physics, 104(24), 9848-9858.
+ */
+void mura(int n, double R, double* r, double* w);
+
+} // end of namespace Radial
+} // end of namespace Grid
+
+#endif

source/module_base/grid/test/CMakeLists.txt

@@ -7,3 +7,8 @@ AddTest(
   ../../ylm.cpp
 )
 
+AddTest(
+  TARGET test_radial
+  SOURCES test_radial.cpp
+  ../radial.cpp
+)

source/module_base/grid/test/test_delley.cpp

@@ -7,7 +7,7 @@
 #include <mpi.h>
 #endif
 
-using namespace Grid;
+using namespace Grid::Angular;
 
 // mock the function to prevent unnecessary dependency
 namespace ModuleBase {
@@ -47,27 +47,27 @@ void DelleyTest::randgen(int lmax, std::vector<double>& coef) {
 
 TEST_F(DelleyTest, NumGrid) {
     int lmax = 5;
-    int ngrid = Delley::ngrid(lmax);
+    int ngrid = ngrid_delley(lmax);
     EXPECT_EQ(lmax, 17);
     EXPECT_EQ(ngrid, 110);
 
     lmax = 17;
-    ngrid = Delley::ngrid(lmax);
+    ngrid = ngrid_delley(lmax);
     EXPECT_EQ(lmax, 17);
     EXPECT_EQ(ngrid, 110);
 
     lmax = 20;
-    ngrid = Delley::ngrid(lmax);
+    ngrid = ngrid_delley(lmax);
     EXPECT_EQ(lmax, 23);
     EXPECT_EQ(ngrid, 194);
 
     lmax = 59;
-    ngrid = Delley::ngrid(lmax);
+    ngrid = ngrid_delley(lmax);
     EXPECT_EQ(lmax, 59);
     EXPECT_EQ(ngrid, 1202);
 
     lmax = 60;
-    ngrid = Delley::ngrid(lmax);
+    ngrid = ngrid_delley(lmax);
     EXPECT_EQ(lmax, 60);
     EXPECT_EQ(ngrid, -1);
 }
@@ -91,7 +91,7 @@ TEST_F(DelleyTest, Accuracy) {
     std::vector<double> grid, weight, coef;
 
     for (int grid_lmax = 17; grid_lmax < 60; grid_lmax +=6) {
-        Delley::get(grid_lmax, grid, weight);
+        delley(grid_lmax, grid, weight);
         int func_lmax = grid_lmax / 2;
         randgen(func_lmax, coef);