Refactor: simplify snap_psibeta_half_tddft.cpp (#6846)

* simplify snap_psibeta_half_tddft.cpp

* fix: use dk_1 instead of inv_dk_1 in Polynomial_Interpolation call

The Polynomial_Interpolation function expects table_interval (dk) as the parameter,
not its inverse. The function computes position = x / table_interval internally.
This fixes the 05_rtTDDFT test failures.

Author: dzzz2001

source/source_lcao/module_rt/snap_psibeta_half_tddft.cpp

@@ -5,26 +5,27 @@
 #include "source_base/math_polyint.h"
 #include "source_base/timer.h"
 #include "source_base/ylm.h"
-#include <vector>
-#include <complex>
+
 #include <cmath>
+#include <complex>
+#include <vector>
 
 namespace module_rt
 {
 
 /**
  * @brief Initialize Gauss-Legendre grid points and weights.
  *        Thread-safe initialization using static local variable.
- * 
+ *
  * @param grid_size Number of grid points (140)
  * @param gl_x Output: Grid points in [-1, 1]
  * @param gl_w Output: Weights
  */
 static void init_gauss_legendre_grid(int grid_size, std::vector<double>& gl_x, std::vector<double>& gl_w)
 {
     static bool init = false;
-    // Thread-safe initialization
-    #pragma omp critical(init_gauss_legendre)
+// Thread-safe initialization
+#pragma omp critical(init_gauss_legendre)
     {
         if (!init)
         {
@@ -34,48 +35,17 @@ static void init_gauss_legendre_grid(int grid_size, std::vector<double>& gl_x, s
     }
 }
 
-/**
- * @brief Interpolate radial function value at a given distance using cubic interpolation.
- * 
- * @param psi Pointer to radial function array
- * @param mesh Number of mesh points
- * @param inv_dk Inverse of grid spacing (1/dk)
- * @param distance Distance to interpolate at
- * @return double Interpolated value
- */
-inline double interpolate_radial(
-    const double* psi, 
-    int mesh, 
-    double inv_dk, 
-    double distance)
-{
-    double position = distance * inv_dk;
-    int iq = static_cast<int>(position);
-    
-    // Boundary check safe-guard
-    if (iq > mesh - 4) return 0.0;
-
-    const double x0 = position - static_cast<double>(iq);
-    const double x1 = 1.0 - x0;
-    const double x2 = 2.0 - x0;
-    const double x3 = 3.0 - x0;
-    
-    // Cubic interpolation formula
-    return x1*x2*(psi[iq]*x3+psi[iq+3]*x0)/6.0
-         + x0*x3*(psi[iq+1]*x2-psi[iq+2]*x1)/2.0;
-}
-
 /**
  * @brief Main function to calculate overlap integrals <phi|exp^{-iAr}|beta>
  *        and its derivatives (if calc_r is true).
- * 
+ *
  *        This function integrates the overlap between a local orbital phi (at R1)
  *        and a non-local projector beta (at R0), modulated by a plane-wave-like phase factor
  *        exp^{-iAr}, where A is a vector potential.
- * 
+ *
  * @param orb LCAO Orbitals information
  * @param infoNL_ Non-local pseudopotential information
- * @param nlm Output: 
+ * @param nlm Output:
  *            nlm[0] : <phi|exp^{-iAr}|beta>
  *            nlm[1, 2, 3] : <phi|r_a * exp^{-iAr}|beta>, a = x, y, z (if calc_r=true)
  * @param R1 Position of atom 1 (orbital phi)
@@ -96,7 +66,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
                              const int& L1,
                              const int& m1,
                              const int& N1,
-                             const ModuleBase::Vector3<double>& R0, 
+                             const ModuleBase::Vector3<double>& R0,
                              const int& T0,
                              const ModuleBase::Vector3<double>& A,
                              const bool& calc_r)
@@ -105,27 +75,29 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
 
     // 1. Initialization and Early Exits
     const int nproj = infoNL_.nproj[T0];
-    
+
     // Resize output vector based on whether position operator matrix elements are needed
     int required_size = calc_r ? 4 : 1;
-    if (nlm.size() != required_size) nlm.resize(required_size);
-    
-    if (nproj == 0) return;
+    if (nlm.size() != required_size)
+        nlm.resize(required_size);
+
+    if (nproj == 0)
+        return;
 
     // 2. Determine total number of projectors and identify active ones based on cutoff
     int natomwfc = 0;
     std::vector<bool> calproj(nproj, false);
-    
+
     const double Rcut1 = orb.Phi[T1].getRcut();
     const ModuleBase::Vector3<double> dRa = R0 - R1;
     const double distance10 = dRa.norm();
-    
+
     bool any_active = false;
     for (int ip = 0; ip < nproj; ip++)
     {
         const int L0 = infoNL_.Beta[T0].Proj[ip].getL();
         natomwfc += 2 * L0 + 1;
-        
+
         const double Rcut0 = infoNL_.Beta[T0].Proj[ip].getRcut();
         if (distance10 <= (Rcut1 + Rcut0))
         {
@@ -135,7 +107,8 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
     }
 
     // Initialize output values to zero and resize inner vectors
-    for (auto& x : nlm) {
+    for (auto& x: nlm)
+    {
         x.assign(natomwfc, 0.0);
     }
 
@@ -150,12 +123,11 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
     const int mesh_r1 = phi_ln.getNr();
     const double* psi_1 = phi_ln.getPsi();
     const double dk_1 = phi_ln.getDk();
-    const double inv_dk_1 = 1.0 / dk_1;
 
     // 4. Prepare Integration Grids
     const int radial_grid_num = 140;
     const int angular_grid_num = 110;
-    
+
     // Cached standard Gauss-Legendre grid
     static std::vector<double> gl_x(radial_grid_num);
     static std::vector<double> gl_w(radial_grid_num);
@@ -169,17 +141,17 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
     std::vector<double> A_dot_lebedev(angular_grid_num);
     for (int ian = 0; ian < angular_grid_num; ++ian)
     {
-        A_dot_lebedev[ian] = A.x * ModuleBase::Integral::Lebedev_Laikov_grid110_x[ian] +
-                             A.y * ModuleBase::Integral::Lebedev_Laikov_grid110_y[ian] +
-                             A.z * ModuleBase::Integral::Lebedev_Laikov_grid110_z[ian];
+        A_dot_lebedev[ian] = A.x * ModuleBase::Integral::Lebedev_Laikov_grid110_x[ian]
+                             + A.y * ModuleBase::Integral::Lebedev_Laikov_grid110_y[ian]
+                             + A.z * ModuleBase::Integral::Lebedev_Laikov_grid110_z[ian];
     }
 
     // Reuseable buffers for inner loops to avoid allocation
     std::vector<std::complex<double>> result_angular; // Accumulator for angular integration
     // Accumulators for position operator components
     std::vector<std::complex<double>> res_ang_x, res_ang_y, res_ang_z;
-    
-    std::vector<double> rly1((L1 + 1) * (L1 + 1)); // Spherical harmonics buffer for L1
+
+    std::vector<double> rly1((L1 + 1) * (L1 + 1));                 // Spherical harmonics buffer for L1
     std::vector<std::vector<double>> rly0_cache(angular_grid_num); // Cache for L0 Ylm
 
     // 5. Loop over Projectors (Beta)
@@ -207,8 +179,8 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
         double r_max = radial0[mesh_r0 - 1];
         double xl = (r_max - r_min) * 0.5;
         double xmean = (r_max + r_min) * 0.5;
-        
-        for(int i=0; i<radial_grid_num; ++i)
+
+        for (int i = 0; i < radial_grid_num; ++i)
         {
             r_radial[i] = xmean + xl * gl_x[i];
             w_radial[i] = xl * gl_w[i];
@@ -219,12 +191,13 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
 
         // 5.2 Precompute Spherical Harmonics (Ylm) for L0 on angular grid
         // Since L0 changes with projector, we compute this per projector loop.
-        for(int ian = 0; ian < angular_grid_num; ++ian) {
-            ModuleBase::Ylm::rl_sph_harm(L0, 
-                ModuleBase::Integral::Lebedev_Laikov_grid110_x[ian], 
-                ModuleBase::Integral::Lebedev_Laikov_grid110_y[ian], 
-                ModuleBase::Integral::Lebedev_Laikov_grid110_z[ian], 
-                rly0_cache[ian]);
+        for (int ian = 0; ian < angular_grid_num; ++ian)
+        {
+            ModuleBase::Ylm::rl_sph_harm(L0,
+                                         ModuleBase::Integral::Lebedev_Laikov_grid110_x[ian],
+                                         ModuleBase::Integral::Lebedev_Laikov_grid110_y[ian],
+                                         ModuleBase::Integral::Lebedev_Laikov_grid110_z[ian],
+                                         rly0_cache[ian]);
         }
 
         // Resize accumulators if needed
@@ -260,21 +233,22 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
                 const double y = ModuleBase::Integral::Lebedev_Laikov_grid110_y[ian];
                 const double z = ModuleBase::Integral::Lebedev_Laikov_grid110_z[ian];
                 const double w_ang = ModuleBase::Integral::Lebedev_Laikov_grid110_w[ian];
-                
+
                 // Vector r = r_val * u_angle
                 double rx = r_val * x;
                 double ry = r_val * y;
                 double rz = r_val * z;
-                
+
                 // Vector r' = r + R0 - R1 = r + dRa
                 double tx = rx + dRa.x;
                 double ty = ry + dRa.y;
                 double tz = rz + dRa.z;
-                
-                double tnorm = std::sqrt(tx*tx + ty*ty + tz*tz);
-                
+
+                double tnorm = std::sqrt(tx * tx + ty * ty + tz * tz);
+
                 // If r' is outside the cutoff of Phi(r'), skip
-                if (tnorm > Rcut1) continue;
+                if (tnorm > Rcut1)
+                    continue;
 
                 // Compute Ylm for L1 at direction r'
                 if (tnorm > 1e-10)
@@ -285,7 +259,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
                 else
                 {
                     // At origin, only Y_00 is non-zero (if using real spherical harmonics convention)
-                    ModuleBase::Ylm::rl_sph_harm(L1, 0.0, 0.0, 1.0, rly1); 
+                    ModuleBase::Ylm::rl_sph_harm(L1, 0.0, 0.0, 1.0, rly1);
                 }
 
                 // Calculate common phase and weight factor
@@ -294,11 +268,11 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
                 const std::complex<double> exp_iAr = std::exp(ModuleBase::IMAG_UNIT * phase);
 
                 // Interpolate Psi at |r'|
-                double interp_psi = interpolate_radial(psi_1, mesh_r1, inv_dk_1, tnorm);
-                
+                double interp_psi = ModuleBase::PolyInt::Polynomial_Interpolation(psi_1, mesh_r1, dk_1, tnorm);
+
                 const int offset_L1 = L1 * L1 + m1;
                 const double ylm_L1_val = rly1[offset_L1];
-                
+
                 // Combined factor: exp(iAr) * Y_L1m1(r') * Psi(|r'|) * weight_angle
                 const std::complex<double> common_factor = exp_iAr * ylm_L1_val * interp_psi * w_ang;
 
@@ -319,7 +293,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
                         double r_op_x = rx + R0.x;
                         double r_op_y = ry + R0.y;
                         double r_op_z = rz + R0.z;
-                        
+
                         res_ang_x[m0] += term * r_op_x;
                         res_ang_y[m0] += term * r_op_y;
                         res_ang_z[m0] += term * r_op_z;
@@ -329,16 +303,14 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
 
             // 5.5 Combine Radial and Angular parts
             // Interpolate Beta(|r|)
-            // Note: The original code implies beta_r stores values that might need scaling or are just the function values.
-            // Typically radial integration is \int f(r) r^2 dr.
-            // Here we have factor: beta_val * r_radial[ir] * w_radial[ir]
-            // w_radial includes the Jacobian for the change of variable from [-1,1] to [r_min, r_max].
-            // The extra r_radial[ir] suggests either beta is stored as r*beta, or we are doing \int ... r dr (2D?), 
-            // or Jacobian r^2 is split. Assuming original logic is correct.
-            
-            double beta_val = ModuleBase::PolyInt::Polynomial_Interpolation(
-                                beta_r, mesh_r0, dk_0, r_radial[ir]);
-            
+            // Note: The original code implies beta_r stores values that might need scaling or are just the function
+            // values. Typically radial integration is \int f(r) r^2 dr. Here we have factor: beta_val * r_radial[ir] *
+            // w_radial[ir] w_radial includes the Jacobian for the change of variable from [-1,1] to [r_min, r_max]. The
+            // extra r_radial[ir] suggests either beta is stored as r*beta, or we are doing \int ... r dr (2D?), or
+            // Jacobian r^2 is split. Assuming original logic is correct.
+
+            double beta_val = ModuleBase::PolyInt::Polynomial_Interpolation(beta_r, mesh_r0, dk_0, r_radial[ir]);
+
             double radial_factor = beta_val * r_radial[ir] * w_radial[ir];
 
             int current_idx = index_offset;
@@ -347,7 +319,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
                 // Final accumulation into global nlm array
                 // Add phase exp(i A * R0)
                 nlm[0][current_idx] += radial_factor * result_angular[m0] * exp_iAR0;
-                
+
                 if (calc_r)
                 {
                     nlm[1][current_idx] += radial_factor * res_ang_x[m0] * exp_iAR0;
@@ -364,11 +336,11 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
 
     // 6. Final Conjugation
     // Apply conjugation to all elements as per convention <phi|beta> = <beta|phi>*
-    for(int dim = 0; dim < nlm.size(); dim++)
+    for (int dim = 0; dim < nlm.size(); dim++)
     {
-        for (auto &x : nlm[dim])
+        for (auto& x: nlm[dim])
         {
-            x = std::conj(x); 
+            x = std::conj(x);
         }
     }