Fix: Numerical error in pyabacus. (#6414)

* Fix numerical error in pyabacus.

* Add test for TwoCenterIntegrator for pyabacus.

* Update numerica_radial.cpp.

Author: ErjieWu

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_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_);