Refactor: Two center integral implemented in module_nao

source/module_basis/module_nao/numerical_radial.cpp

@@ -365,17 +365,21 @@ void NumericalRadial::radtab(const char op,
                              const int l,
                              double* const table,
                              const int nr_tab,
-                             const double* const rgrid_tab,
+                             const double rmax_tab,
                              const bool deriv) const
 {
     assert(op == 'S' || op == 'I' || op == 'T' || op == 'U');
     assert(l >= 0);
-    assert(rgrid_tab && nr_tab > 0);
+    assert(rmax_tab > 0 && nr_tab > 0);
 
     // radtab requires that two NumericalRadial objects have exactly the same (non-null) kgrid_
     assert(nk_ > 0 && nk_ == ket.nk_);
     assert(std::equal(kgrid_, kgrid_ + nk_, ket.kgrid_));
 
+    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); });
+
     bool use_radrfft = is_fft_compliant(nr_tab, rgrid_tab, nk_, kgrid_);
 
     // function to undergo a spherical Bessel transform:
@@ -407,6 +411,7 @@ void NumericalRadial::radtab(const char op,
     }
 
     delete[] fk;
+    delete[] rgrid_tab;
 
     // spherical Bessel transform has a prefactor of sqrt(2/pi) while the prefactor for the radial table
     // of two-center integrals is 4*pi

source/module_basis/module_nao/numerical_radial.h

@@ -191,17 +191,17 @@ class NumericalRadial
      *                          /  0                  l
      *
      *                                                                                  */
-    void radtab(const char op,                 //!< [in] operator, could be:
-                                               //!<        - 'S' or 'I': overlap
-                                               //!<        - 'T': kinetic
-                                               //!<        - 'U': Coulomb
-                const NumericalRadial& ket,    //!< [in] the other NumericalRadial object with which
-                                               //!       the two-center integral is computed
-                const int l,                   //!< [in] angular momentum of the table
-                double* const table,           //!< [out] on finish, contain the computed table
-                const int nr_tab,              //!< [in] size of table grid
-                const double* const rgrid_tab, //!< [in] grid on which the table is calculated.
-                const bool deriv = false       //!< [in] if true, calculates the derivative of the table
+    void radtab(const char op,              //!< [in] operator, could be:
+                                            //!<        - 'S' or 'I': overlap
+                                            //!<        - 'T': kinetic
+                                            //!<        - 'U': Coulomb
+                const NumericalRadial& ket, //!< [in] the other NumericalRadial object with which
+                                            //!       the two-center integral is computed
+                const int l,                //!< [in] angular momentum of the table
+                double* const table,        //!< [out] on finish, contain the computed table
+                const int nr_tab,           //!< [in] size of table grid
+                const double rmax_tab,      //!< [in] cutoff radius of table grid
+                const bool deriv = false    //!< [in] if true, calculates the derivative of the table
     ) const;
 
     //! Normalizes the radial function.

source/module_basis/module_nao/test/CMakeLists.txt

@@ -60,3 +60,17 @@ AddTest(
   LIBS ${math_libs} device base container kernels
 )
 
+AddTest(
+  TARGET two_center_integrator
+  SOURCES
+    two_center_integrator_test.cpp
+    ../two_center_integrator.cpp
+    ../two_center_table.cpp
+    ../real_gaunt_table.cpp
+    ../radial_collection.cpp
+    ../atomic_radials.cpp
+    ../beta_radials.cpp
+    ../radial_set.cpp
+    ../numerical_radial.cpp
+  LIBS ${math_libs} device base container kernels orb
+)

source/module_basis/module_nao/test/numerical_radial_test.cpp

@@ -472,58 +472,59 @@ TEST_F(NumericalRadialTest, RadialTable)
 
     double* table = new double[sz];
     double table_pref = ModuleBase::FOUR_PI * std::sqrt(ModuleBase::PI / 2.0);
+    double rmax_tab = chi1.rcut();
 
-    chi1.radtab('S', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid());
+    chi1.radtab('S', chi2, 0, table, chi1.nr(), rmax_tab);
     for (int i = 0; i != sz; ++i)
     {
         double R = i * dr;
         EXPECT_NEAR(table[i], table_pref * (3 - R * R) / 32 * std::exp(-R * R / 2), tol);
     }
 
-    chi1.radtab('S', chi2, 2, table, chi1.nr(), chi1.ptr_rgrid());
+    chi1.radtab('S', chi2, 2, table, chi1.nr(), rmax_tab);
     for (int i = 0; i != sz; ++i)
     {
         double R = i * dr;
         EXPECT_NEAR(table[i], table_pref * R * R / 32 * std::exp(-R * R / 2), tol);
     }
 
-    chi1.radtab('T', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid());
+    chi1.radtab('T', chi2, 0, table, chi1.nr(), rmax_tab);
     for (int i = 0; i != sz; ++i)
     {
         double R = i * dr;
         EXPECT_NEAR(table[i], table_pref * (std::pow(R, 4) - 10 * R * R + 15) / 32 * std::exp(-R * R / 2), tol);
     }
 
-    chi1.radtab('U', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid());
+    chi1.radtab('U', chi2, 0, table, chi1.nr(), rmax_tab);
     for (int i = 0; i != sz; ++i)
     {
         double R = i * dr;
         EXPECT_NEAR(table[i], table_pref * 1. / 32 * std::exp(-R * R / 2), tol);
     }
 
     // derivative of radial tables
-    chi1.radtab('S', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid(), true);
+    chi1.radtab('S', chi2, 0, table, chi1.nr(), rmax_tab, true);
     for (int i = 0; i != sz; ++i)
     {
         double R = i * dr;
         EXPECT_NEAR(table[i], table_pref * (std::pow(R, 3) - 5 * R) / 32 * std::exp(-R * R / 2), tol);
     }
 
-    chi1.radtab('S', chi2, 2, table, chi1.nr(), chi1.ptr_rgrid(), true);
+    chi1.radtab('S', chi2, 2, table, chi1.nr(), rmax_tab, true);
     for (int i = 0; i != sz; ++i)
     {
         double R = i * dr;
         EXPECT_NEAR(table[i], table_pref * (2 * R - std::pow(R, 3)) / 32 * std::exp(-R * R / 2), tol);
     }
 
-    chi1.radtab('T', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid(), true);
+    chi1.radtab('T', chi2, 0, table, chi1.nr(), rmax_tab, true);
     for (int i = 0; i != sz; ++i)
     {
         double R = i * dr;
         EXPECT_NEAR(table[i], table_pref * (-std::pow(R, 5) + 14 * std::pow(R, 3) - 35 * R) / 32 * std::exp(-R * R / 2), tol);
     }
 
-    chi1.radtab('U', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid(), true);
+    chi1.radtab('U', chi2, 0, table, chi1.nr(), rmax_tab, true);
     for (int i = 0; i != sz; ++i)
     {
         double R = i * dr;