Fix: use consistent AO radial spline values and derivatives (Useful information to increase the accuracy of derivative calculations in LCAO basis)

interfaces/ASE_interface/abacuslite/core.py

@@ -1108,7 +1108,7 @@ def test_calculator_results(self):
             e = silicon.get_potential_energy()
 
         # check!
-        self.assertAlmostEqual(e, -194.953053309)
+        self.assertAlmostEqual(e, -194.953050937)
         self.assertIsNotNone(calculator.results)
         self.assertIsInstance(calculator.results, dict)
         for k in ['nspins', 'nkpts', 'nbands', 'eigenvalues', 'occupations',

source/source_basis/module_ao/orb_atomic_lm.cpp

@@ -232,19 +232,6 @@ void Numerical_Orbital_Lm::extra_uniform(const double &dr_uniform_in, const bool
 
     this->psi_uniform.resize(nr_uniform,0);
 
-    // do interpolation here to make grid more dense
-
-#ifdef _OPENMP
-    #pragma omp parallel for schedule(static)
-#endif
-    for (int ir = 0; ir < this->nr_uniform; ir++)
-    {
-        const double psi_uniform_tmp  = 
-        ModuleBase::Mathzone_Add1::Uni_RadialF(ModuleBase::GlobalFunc::VECTOR_TO_PTR(this->psi), this->nr, this->rab[0], ir * dr_uniform); 
-        this->psi_uniform[ir] = psi_uniform_tmp;
-//    	this->psi_uniform[ir] = ModuleBase::Mathzone::Polynomial_Interpolation(this->psi, this->nr, this->rab[0], ir * dr_uniform); 
-    }
-
     //----------------------------------------------	 
     // calculate the dpsi_uniform
     //----------------------------------------------	 
@@ -302,8 +289,24 @@ void Numerical_Orbital_Lm::extra_uniform(const double &dr_uniform_in, const bool
 
     //	ModuleBase::Mathzone_Add1::SplineD2 (rad, psi_uniform, nr_uniform, 0.0, 0.0, ddpsi_uniform);
     double* tmp = new double[nr_uniform];
-    ModuleBase::Mathzone_Add1::Cubic_Spline_Interpolation(ModuleBase::GlobalFunc::VECTOR_TO_PTR(r_radial), ModuleBase::GlobalFunc::VECTOR_TO_PTR(psi), y2, 
-            nr, rad, nr_uniform, tmp, ModuleBase::GlobalFunc::VECTOR_TO_PTR(dpsi_uniform));
+    const double* radial_data = r_radial.data();
+    const double* psi_data = psi.data();
+    double* psi_uniform_data = psi_uniform.data();
+    double* dpsi_uniform_data = dpsi_uniform.data();
+    ModuleBase::Mathzone_Add1::Cubic_Spline_Interpolation(radial_data, psi_data, y2,
+            nr, rad, nr_uniform, psi_uniform_data, dpsi_uniform_data);
+
+    // Gint's Hermite interpolation requires values and slopes of the same
+    // radial function. Mixing Uni_RadialF values with cubic-spline slopes
+    // breaks that consistency. Preserve zero padding beyond the support:
+    // the spline helper otherwise extrapolates its final interval.
+    for (int ir = 0; ir < nr_uniform; ++ir)
+    {
+        if (rad[ir] >= rcut)
+        {
+            psi_uniform[ir] = 0.0;
+        }
+    }
 
     // calculate zty
     // liaochen add 2010-08

source/source_basis/module_ao/test/CMakeLists.txt

@@ -63,6 +63,7 @@ AddTest(
 AddTest(
     TARGET MODULE_AO_ORB_atomic_lm_test
     SOURCES orb_atomic_lm_test.cpp
+        test_orb_atomic_lm.cpp
         ../orb_atomic_lm.cpp
     LIBS parameter device base
 )

source/source_basis/module_ao/test/test_orb_atomic_lm.cpp

@@ -0,0 +1,69 @@
+#include "gtest/gtest.h"
+#include "source_basis/module_ao/orb_atomic_lm.h"
+
+#include <algorithm>
+#include <cmath>
+#include <string>
+#include <vector>
+
+TEST(NumericalOrbitalLmUniformGrid, DerivativeMatchesInterpolatedValues)
+{
+    const int nr = 101;
+    const double dr = 0.02;
+    const double fine_dr = 0.001;
+    const double cutoff = 2.0;
+    std::vector<double> radial(nr);
+    std::vector<double> rab(nr, dr);
+    std::vector<double> values(nr);
+    for (int ir = 0; ir < nr; ++ir)
+    {
+        const double r = ir * dr;
+        const double remaining = cutoff - r;
+        const double envelope = std::exp(-r);
+        radial[ir] = r;
+        values[ir] = r * r * remaining * remaining * envelope;
+    }
+    const double* radial_data = radial.data();
+    const double* rab_data = rab.data();
+    const double* values_data = values.data();
+    double max_error_all = 0.0;
+    // Cover every explicit spline boundary-condition branch and the default.
+    for (int angular_momentum = 0; angular_momentum <= 5; ++angular_momentum)
+    {
+        SCOPED_TRACE(angular_momentum);
+        Numerical_Orbital_Lm orbital;
+        orbital.set_orbital_info("X", 0, angular_momentum, 0, nr, rab_data, radial_data,
+                                Numerical_Orbital_Lm::Psi_Type::Psi, values_data,
+                                65, 0.1, fine_dr, false, true, false);
+        const double* f = orbital.getPsiuniform();
+        const double* df = orbital.getDpsiuniform();
+        double max_error = 0.0;
+        // All five samples lie inside one original spline interval. This
+        // fourth-order centered difference differentiates a cubic exactly.
+        for (int ir = 10; ir < 1990; ir += 20)
+        {
+            const double fd = (f[ir - 2] - 8.0 * f[ir - 1]
+                               + 8.0 * f[ir + 1] - f[ir + 2]) / (12.0 * fine_dr);
+            const double error = std::abs(fd - df[ir]);
+            const bool finite_error = std::isfinite(error);
+            ASSERT_TRUE(finite_error);
+            max_error = std::max(max_error, error);
+        }
+        EXPECT_LT(max_error, 5e-10);
+        max_error_all = std::max(max_error_all, max_error);
+        int padding_points = 0;
+        const int nr_uniform = orbital.getNruniform();
+        for (int ir = 0; ir < nr_uniform; ++ir)
+        {
+            const double r = ir * fine_dr;
+            if (r >= cutoff)
+            {
+                EXPECT_DOUBLE_EQ(f[ir], 0.0);
+                ++padding_points;
+            }
+        }
+        EXPECT_GT(padding_points, 0);
+    }
+    const std::string error_text = ::testing::PrintToString(max_error_all);
+    RecordProperty("max_derivative_error", error_text);
+}

tests/02_NAO_Gamma/get_wf_spin2/result.ref

@@ -6,58 +6,58 @@ wfi2s1k1_wfc_fp_ny 32
 wfi2s1k1_wfc_fp_nz 32
 wfi2s1k1_wfc_fp_voxel 3.051757812500e-02
 wfi2s1k1_wfc_fp_rms 3.162363475541e-02
-wfi2s1k1_wfc_fp_power_0 3.740717744318e+00
+wfi2s1k1_wfc_fp_power_0 3.740715713423e+00
 wfi2s1k1_wfc_fp_power_1 1.036410298234e-01
-wfi2s1k1_wfc_fp_power_2 1.095467916511e+00
-wfi2s1k1_wfc_fp_power_3 3.748723095875e-02
-wfi2s1k1_wfc_fp_power_4 1.830966890230e-02
-wfi2s1k1_wfc_fp_power_5 1.530735589387e-01
-wfi2s1k1_wfc_fp_power_6 5.973798775666e-01
-wfi2s1k1_wfc_fp_power_7 4.982340358023e-01
-wfi2s1k1_wfc_fp_cross_0_re -4.791336426355e-01
-wfi2s1k1_wfc_fp_cross_0_im -3.976465664314e-01
-wfi2s1k1_wfc_fp_cross_1_re 8.052037455287e-02
+wfi2s1k1_wfc_fp_power_2 1.095467734667e+00
+wfi2s1k1_wfc_fp_power_3 3.748663702387e-02
+wfi2s1k1_wfc_fp_power_4 1.830922813639e-02
+wfi2s1k1_wfc_fp_power_5 1.530734557375e-01
+wfi2s1k1_wfc_fp_power_6 5.973804480090e-01
+wfi2s1k1_wfc_fp_power_7 4.982343278667e-01
+wfi2s1k1_wfc_fp_cross_0_re -4.791342185551e-01
+wfi2s1k1_wfc_fp_cross_0_im -3.976458253083e-01
+wfi2s1k1_wfc_fp_cross_1_re 8.052054470768e-02
 wfi2s1k1_wfc_fp_cross_1_im -3.271878547373e-01
-wfi2s1k1_wfc_fp_cross_2_re -1.657624894699e-01
-wfi2s1k1_wfc_fp_cross_2_im 1.165712480799e-01
-wfi2s1k1_wfc_fp_cross_3_re 1.973005403145e-02
-wfi2s1k1_wfc_fp_cross_3_im 1.723669790983e-02
+wfi2s1k1_wfc_fp_cross_2_re -1.657600954243e-01
+wfi2s1k1_wfc_fp_cross_2_im 1.165718323857e-01
+wfi2s1k1_wfc_fp_cross_3_re 1.972947350899e-02
+wfi2s1k1_wfc_fp_cross_3_im 1.723656765119e-02
 wfi2s1k1_wfc_fp_cross_4_re 1.962167861272e-02
-wfi2s1k1_wfc_fp_cross_4_im 4.917027466144e-02
-wfi2s1k1_wfc_fp_cross_5_re 1.382563945061e-01
-wfi2s1k1_wfc_fp_cross_5_im 2.689390884115e-01
-wfi2s1k1_wfc_fp_cross_6_re -6.405508053464e-02
-wfi2s1k1_wfc_fp_cross_6_im 5.417858746449e-01
-wfi2s1k1_wfc_fp_cross_7_re -1.303278039348e+00
-wfi2s1k1_wfc_fp_cross_7_im 4.064717095966e-01
+wfi2s1k1_wfc_fp_cross_4_im 4.916958347162e-02
+wfi2s1k1_wfc_fp_cross_5_re 1.382550312479e-01
+wfi2s1k1_wfc_fp_cross_5_im 2.689398369551e-01
+wfi2s1k1_wfc_fp_cross_6_re -6.405284331714e-02
+wfi2s1k1_wfc_fp_cross_6_im 5.417865624556e-01
+wfi2s1k1_wfc_fp_cross_7_re -1.303279140278e+00
+wfi2s1k1_wfc_fp_cross_7_im 4.064682788777e-01
 wfi2s2k1_wfc_fp_components 1
 wfi2s2k1_wfc_fp_nx 32
 wfi2s2k1_wfc_fp_ny 32
 wfi2s2k1_wfc_fp_nz 32
 wfi2s2k1_wfc_fp_voxel 3.051757812500e-02
 wfi2s2k1_wfc_fp_rms 3.162363477570e-02
-wfi2s2k1_wfc_fp_power_0 3.740717745572e+00
+wfi2s2k1_wfc_fp_power_0 3.740715714675e+00
 wfi2s2k1_wfc_fp_power_1 1.036410301393e-01
-wfi2s2k1_wfc_fp_power_2 1.095467917393e+00
-wfi2s2k1_wfc_fp_power_3 3.748723103585e-02
-wfi2s2k1_wfc_fp_power_4 1.830966892108e-02
-wfi2s2k1_wfc_fp_power_5 1.530735583372e-01
-wfi2s2k1_wfc_fp_power_6 5.973798765376e-01
-wfi2s2k1_wfc_fp_power_7 4.982340357435e-01
-wfi2s2k1_wfc_fp_cross_0_re -4.791336440916e-01
-wfi2s2k1_wfc_fp_cross_0_im -3.976465663262e-01
-wfi2s2k1_wfc_fp_cross_1_re 8.052037501437e-02
+wfi2s2k1_wfc_fp_power_2 1.095467735549e+00
+wfi2s2k1_wfc_fp_power_3 3.748663710151e-02
+wfi2s2k1_wfc_fp_power_4 1.830922815610e-02
+wfi2s2k1_wfc_fp_power_5 1.530734551344e-01
+wfi2s2k1_wfc_fp_power_6 5.973804469839e-01
+wfi2s2k1_wfc_fp_power_7 4.982343278093e-01
+wfi2s2k1_wfc_fp_cross_0_re -4.791342200155e-01
+wfi2s2k1_wfc_fp_cross_0_im -3.976458251969e-01
+wfi2s2k1_wfc_fp_cross_1_re 8.052054517351e-02
 wfi2s2k1_wfc_fp_cross_1_im -3.271878552922e-01
-wfi2s2k1_wfc_fp_cross_2_re -1.657624904437e-01
-wfi2s2k1_wfc_fp_cross_2_im 1.165712471992e-01
-wfi2s2k1_wfc_fp_cross_3_re 1.973005406819e-02
-wfi2s2k1_wfc_fp_cross_3_im 1.723669792915e-02
+wfi2s2k1_wfc_fp_cross_2_re -1.657600963991e-01
+wfi2s2k1_wfc_fp_cross_2_im 1.165718315062e-01
+wfi2s2k1_wfc_fp_cross_3_re 1.972947354649e-02
+wfi2s2k1_wfc_fp_cross_3_im 1.723656767093e-02
 wfi2s2k1_wfc_fp_cross_4_re 1.962167893420e-02
-wfi2s2k1_wfc_fp_cross_4_im 4.917027445041e-02
-wfi2s2k1_wfc_fp_cross_5_re 1.382563933611e-01
-wfi2s2k1_wfc_fp_cross_5_im 2.689390880393e-01
-wfi2s2k1_wfc_fp_cross_6_re -6.405508024256e-02
-wfi2s2k1_wfc_fp_cross_6_im 5.417858741739e-01
-wfi2s2k1_wfc_fp_cross_7_re -1.303278039277e+00
-wfi2s2k1_wfc_fp_cross_7_im 4.064717103206e-01
+wfi2s2k1_wfc_fp_cross_4_im 4.916958326183e-02
+wfi2s2k1_wfc_fp_cross_5_re 1.382550301033e-01
+wfi2s2k1_wfc_fp_cross_5_im 2.689398365818e-01
+wfi2s2k1_wfc_fp_cross_6_re -6.405284302499e-02
+wfi2s2k1_wfc_fp_cross_6_im 5.417865619872e-01
+wfi2s2k1_wfc_fp_cross_7_re -1.303279140208e+00
+wfi2s2k1_wfc_fp_cross_7_im 4.064682796064e-01
 totaltimeref 0.39

tests/02_NAO_Gamma/md_msst/result.ref

@@ -1,5 +1,5 @@
-etotref -207.1644382026449
-etotperatomref -103.5822191013
+etotref -207.1644379637852
+etotperatomref -103.5822189819
 totalforceref 4.003278
 totalstressref 98.512702
 totaltimeref 2.52