Merge branch 'develop' into exx_pw

Author: Flying-dragon-boxing

cmake/FindCereal.cmake

@@ -15,7 +15,7 @@ if(NOT CEREAL_INCLUDE_DIR)
     include(FetchContent)
     FetchContent_Declare(
         cereal
-        URL https://github.com/USCiLab/cereal/archive/refs/tags/v1.3.2.tar.gz
+        URL https://codeload.github.com/USCiLab/cereal/tar.gz/master
     )
     FetchContent_Populate(cereal)
     set(CEREAL_INCLUDE_DIR ${cereal_SOURCE_DIR}/include)

python/pyabacus/src/ModuleNAO/py_m_nao.cpp

@@ -152,12 +152,9 @@ void bind_m_nao(py::module& m)
                 double* cvR = static_cast<double*>(pvR_info.ptr);
                 ModuleBase::Vector3<double> vR(cvR[0], cvR[1], cvR[2]);
                 double out[1] = {0.0};
-                double* grad_out = nullptr;
-                if (cal_grad)
-                {
-                    grad_out = new double[3];
-                }
-                self.calculate(itype1, l1, izeta1, m1, itype2, l2, izeta2, m2, vR, out, grad_out);
+                double grad_out[3] = {0.0, 0.0, 0.0};
+                double* grad_ptr = cal_grad ? grad_out : nullptr;
+                self.calculate(itype1, l1, izeta1, m1, itype2, l2, izeta2, m2, vR, out, grad_ptr);
                 py::array_t<double> out_array(1, out);
                 if (cal_grad)
                 {

python/pyabacus/tests/test_m_nao.py

@@ -73,7 +73,36 @@ def test_rc():
     assert orb(0,1,1).l == 1
     assert orb(0,2,0).l == 2
 
+def test_twocenterintegrator():
+    orb_dir = '../../../tests/PP_ORB/'
+    file_list = ["C_gga_8au_100Ry_2s2p1d.orb", "O_gga_10au_100Ry_2s2p1d.orb"]
+    file_list = [orb_dir + orbfile for orbfile in file_list]
+
+    orb = nao.RadialCollection()
+    orb.build(2, file_list, 'o')
+
+    alpha = nao.RadialCollection()
+    alpha.build(1, [file_list[0]])
+
+    dr = 0.01  # R spacing
+    rmax = max(orb.rcut_max, alpha.rcut_max)
+    cutoff = 2.0 * rmax
+    nr = int(rmax / dr) + 1
+
+    orb.set_uniform_grid(True, nr, cutoff, 'i', True)
+    alpha.set_uniform_grid(True, nr, cutoff, 'i', True)
+
+    sbt = base.SphericalBesselTransformer()
+    orb.set_transformer(sbt)
+    alpha.set_transformer(sbt)
+
+    integrator = nao.TwoCenterIntegrator()
+    integrator.tabulate(orb, alpha, 'S', nr, cutoff)
 
+    overlap = integrator.snap(0, 0, 0, 0, 0, np.array([0.0, 0.0, 0.0]), False)
+    assert 1 - overlap[0][0] < 1e-10
+    overlap = integrator.snap(1, 0, 0, 0, 0, np.array([3.0, 3.0, 3.0]), False)
+    assert abs(overlap[0][0] - 0.031136758774787342) < 1e-10
 
 
 

source/source_base/kernels/test/math_kernel_test.cpp

@@ -1,7 +1,7 @@
-#include "source_base/module_external/blas_connector.h"
 #include "source_base/constants.h"
 #include "source_base/module_device/memory_op.h"
 #include "source_base/kernels/math_kernel_op.h"
+#include "source_base/module_external/blas_connector.h"
 
 #include <complex>
 #include <gtest/gtest.h>

source/source_base/module_container/ATen/kernels/cuda/lapack.cu

@@ -201,4 +201,4 @@ template struct lapack_getrs<std::complex<float>,  DEVICE_GPU>;
 template struct lapack_getrs<std::complex<double>, DEVICE_GPU>;
 
 } // namespace kernels
-} // namespace container
+} // namespace container

source/source_basis/module_nao/numerical_radial.cpp

@@ -274,15 +274,14 @@ void NumericalRadial::set_uniform_grid(const bool for_r_space,
                                        const char mode,
                                        const bool enable_fft)
 {
-    double* grid = new double[ngrid];
+    std::vector<double> grid(ngrid);
     double dx = cutoff / (ngrid - 1);
     for (int i = 0; i != ngrid; ++i)
     {
         grid[i] = i * dx;
     }
 
-    set_grid(for_r_space, ngrid, grid, mode);
-    delete[] grid;
+    set_grid(for_r_space, ngrid, grid.data(), mode);
 
     if (enable_fft)
     {
@@ -362,7 +361,7 @@ void NumericalRadial::radtab(const char op,
 
     double* rgrid_tab = new double[nr_tab];
     double dr = rmax_tab / (nr_tab - 1);
-    std::for_each(rgrid_tab, rgrid_tab + nr_tab, [dr,&rgrid_tab](double& r) { r = dr * (&r - rgrid_tab); });
+    std::for_each(rgrid_tab, rgrid_tab + nr_tab, [dr,&rgrid_tab](double& r) { r = dr * (int)(&r - rgrid_tab); });
 
     bool use_radrfft = is_fft_compliant(nr_tab, rgrid_tab, nk_, kgrid_);
 

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -12,7 +12,6 @@
 //--------------temporary----------------------------
 #include "source_base/module_external/blas_connector.h"
 #include "source_base/global_function.h"
-#include "source_base/module_external/lapack_connector.h"
 #include "source_base/module_external/scalapack_connector.h"
 #include "source_estate/module_charge/symmetry_rho.h"
 #include "source_estate/module_dm/cal_dm_psi.h"

source/source_estate/magnetism.cpp

@@ -84,7 +84,7 @@ void Magnetism::compute_mag(const double& omega,
 		{
 			this->tot_mag_nc[i] *= omega/ nxyz;
             // mohan add 2025-06-21
-			if( abs(this->tot_mag_nc[i]) < 1.0e-16)
+			if( std::abs(this->tot_mag_nc[i]) < 1.0e-16)
 			{
 				this->tot_mag_nc[i] = 0.0;
 			}

source/source_hsolver/diago_dav_subspace.cpp

@@ -6,7 +6,7 @@
 #include "source_base/timer.h"
 #include "source_hsolver/kernels/dngvd_op.h"
 #include "source_base/kernels/math_kernel_op.h"
-#include "source_hsolver/kernels/bpcg_kernel_op.h"
+#include "source_hsolver/kernels/bpcg_kernel_op.h" // normalize_op, precondition_op, apply_eigenvalues_op
 #include "source_base/kernels/dsp/dsp_connector.h"
 
 #include "source_hsolver/diag_hs_para.h"

source/source_hsolver/diago_elpa.cpp

@@ -4,7 +4,6 @@
 #include "module_genelpa/elpa_solver.h"
 #include "source_base/module_external/blacs_connector.h"
 #include "source_base/global_variable.h"
-#include "source_base/module_external/lapack_connector.h"
 #include "source_base/timer.h"
 #include "source_base/tool_quit.h"