Merge branch 'deepmodeling:develop' into pchg_wf_param

python/pyabacus/src/py_diago_dav_subspace.hpp

@@ -110,11 +110,14 @@ class PyDiagoDavSubspace
         bool scf_type,
         hsolver::diag_comm_info comm_info
     ) {
-        auto hpsi_func = [mm_op] (std::complex<double> *hpsi_out,
-                    std::complex<double> *psi_in, const int nband_in,
-                    const int nbasis_in, const int band_index1,
-                    const int band_index2) 
-        {
+        auto hpsi_func = [mm_op] (
+            std::complex<double> *psi_in,
+            std::complex<double> *hpsi_out, 
+            const int nband_in,
+            const int nbasis_in, 
+            const int band_index1,
+            const int band_index2
+        ) {
             // Note: numpy's py::array_t is row-major, but
             //       our raw pointer-array is column-major
             py::array_t<std::complex<double>, py::array::f_style> psi({nbasis_in, band_index2 - band_index1 + 1});

python/pyabacus/src/py_diago_david.hpp

@@ -109,8 +109,8 @@ class PyDiagoDavid
         hsolver::diag_comm_info comm_info
     ) {
         auto hpsi_func = [mm_op] (
-            std::complex<double> *hpsi_out,
-            std::complex<double> *psi_in, 
+            std::complex<double> *psi_in,
+            std::complex<double> *hpsi_out, 
             const int nband_in, 
             const int nbasis_in, 
             const int band_index1, 

python/pyabacus/src/pyabacus/hsolver/_hsolver.py

@@ -16,7 +16,7 @@ def rank(self) -> int: ...
     def nproc(self) -> int: ...
 
 def dav_subspace(
-    mm_op: Callable[[NDArray[np.complex128]], NDArray[np.complex128]],
+    mvv_op: Callable[[NDArray[np.complex128]], NDArray[np.complex128]],
     init_v: NDArray[np.complex128],
     dim: int,
     num_eigs: int,
@@ -32,9 +32,10 @@ def dav_subspace(
 
     Parameters
     ----------
-    mm_op : Callable[[NDArray[np.complex128]], NDArray[np.complex128]],
-        The operator to be diagonalized, which is a function that takes a matrix as input
-        and returns a matrix mv_op(X) = H * X as output.
+    mvv_op : Callable[[NDArray[np.complex128]], NDArray[np.complex128]],
+        The operator to be diagonalized, which is a function that takes a set of 
+        vectors X = [x1, ..., xN] as input and returns a matrix(vector block)
+        mvv_op(X) = H * X ([Hx1, ..., HxN]) as output.
     init_v : NDArray[np.complex128]
         The initial guess for the eigenvectors.
     dim : int
@@ -68,8 +69,8 @@ def dav_subspace(
     v : NDArray[np.complex128]
         The eigenvectors corresponding to the eigenvalues.
     """
-    if not callable(mm_op):
-        raise TypeError("mm_op must be a callable object.")
+    if not callable(mvv_op):
+        raise TypeError("mvv_op must be a callable object.")
 
     if is_occupied is None:
         is_occupied = [True] * num_eigs
@@ -86,7 +87,7 @@ def dav_subspace(
     assert dav_ndim * num_eigs < dim * comm_info.nproc, "dav_ndim * num_eigs must be less than dim * comm_info.nproc."
 
     _ = _diago_obj_dav_subspace.diag(
-        mm_op,
+        mvv_op,
         pre_condition,
         dav_ndim,
         tol,
@@ -103,7 +104,7 @@ def dav_subspace(
     return e, v
 
 def davidson(
-    mm_op: Callable[[NDArray[np.complex128]], NDArray[np.complex128]],
+    mvv_op: Callable[[NDArray[np.complex128]], NDArray[np.complex128]],
     init_v: NDArray[np.complex128],
     dim: int,
     num_eigs: int,
@@ -119,9 +120,10 @@ def davidson(
 
     Parameters
     ----------
-    mm_op : Callable[[NDArray[np.complex128]], NDArray[np.complex128]],
-        The operator to be diagonalized, which is a function that takes a matrix as input
-        and returns a matrix mv_op(X) = H * X as output.
+    mvv_op : Callable[[NDArray[np.complex128]], NDArray[np.complex128]],
+        The operator to be diagonalized, which is a function that takes a set of 
+        vectors X = [x1, ..., xN] as input and returns a matrix(vector block)
+        mvv_op(X) = H * X ([Hx1, ..., HxN]) as output.
     init_v : NDArray[np.complex128]
         The initial guess for the eigenvectors.
     dim : int
@@ -146,8 +148,8 @@ def davidson(
     v : NDArray[np.complex128]
         The eigenvectors corresponding to the eigenvalues.
     """
-    if not callable(mm_op):
-        raise TypeError("mm_op must be a callable object.")
+    if not callable(mvv_op):
+        raise TypeError("mvv_op must be a callable object.")
 
     if init_v.ndim != 1 or init_v.dtype != np.complex128:
         init_v = init_v.flatten().astype(np.complex128, order='C')
@@ -159,7 +161,7 @@ def davidson(
     comm_info = hsolver.diag_comm_info(0, 1)
 
     _ = _diago_obj_dav_subspace.diag(
-        mm_op,
+        mvv_op,
         pre_condition,
         dav_ndim,
         tol,

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