feat(xc): add SCANL meta-GGA functional with FD Laplacian kernel

source/source_hamilt/module_xc/libxc_abacus.h

@@ -105,11 +105,44 @@ namespace XC_Functional_Libxc
         const double tpiba,
         const Charge* const chr);
 
+    extern void cal_gdr_and_lapl(
+        const int nspin,
+        const std::size_t nrxx,
+        const std::vector<double> &rho,
+        const double tpiba,
+        const Charge* const chr,
+        std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr,
+        std::vector<double> &lapl);
+
     // converting grho (abacus=>libxc)
     extern std::vector<double> convert_sigma(
         const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr);
 
-    // sgn for threshold mask
+    /// Calculate Laplacian of density using spectral method: ∇²ρ = IFFT(-|G|²·FFT(ρ))
+    /// @see cal_lapl() implementation in libxc_tools.cpp for full documentation
+    extern std::vector<double> cal_lapl(
+        const int nspin,
+        const std::size_t nrxx,
+        const std::vector<double> &rho,
+        const double tpiba,
+        const Charge* const chr);
+
+    /// Calculate Laplacian of density using finite-difference kernel (bounded at high G)
+    /// @see cal_lapl_fd() implementation in libxc_tools.cpp for full documentation
+    extern std::vector<double> cal_lapl_fd(
+        const int nspin,
+        const std::size_t nrxx,
+        const std::vector<double> &rho,
+        const Charge* const chr);
+
+    /// Calculate density Hessian: H_ab = ∂²ρ/∂r_a∂r_b (6 independent components per spin)
+    /// @see cal_rho_hessian() implementation in libxc_tools.cpp for full documentation
+    extern std::vector<std::vector<double>> cal_rho_hessian(
+        const int nspin,
+        const std::size_t nrxx,
+        const std::vector<double> &rho,
+        const Charge* const chr);
+
     extern std::vector<double> cal_sgn(
         const double rho_threshold,
         const double grho_threshold,
@@ -209,25 +242,55 @@ namespace XC_Functional_Libxc
 //  libxc_mgga_wrap.cpp
 //-------------------
 
-    // wrapper for the mGGA functionals
+    /// Wrapper for meta-GGA functionals (single spin channel).
+    /// Computes XC energy density and potentials including Laplacian-dependent terms.
+    ///
+    /// @param func_id  LibXC functional IDs (exchange + correlation)
+    /// @param rho      Electron density at grid point
+    /// @param grho     |∇ρ|² at grid point
+    /// @param lapl_rho ∇²ρ at grid point (density Laplacian; pass 0.0 for SCAN)
+    /// @param atau     Kinetic energy density τ at grid point
+    /// @param sxc      [out] XC energy density
+    /// @param v1xc     [out] ∂ε_xc/∂ρ
+    /// @param v2xc     [out] 2·∂ε_xc/∂(|∇ρ|²)
+    /// @param v3xc     [out] ∂ε_xc/∂τ
+    /// @param vlapl    [out] ∂ε_xc/∂(∇²ρ) (Laplacian potential; 0.0 for SCAN)
+    /// @param hybrid_alpha  Exact exchange mixing fraction (0.0 for pure GGA/meta-GGA)
+    /// @param hse_omega     Range-separation parameter for HSE-type functionals
+    ///
+    /// @note For SCAN (non-Laplacian meta-GGA), pass lapl_rho=0.0 and vlapl=0.0.
+    ///       For SCANL, pass the actual ∇²ρ value and vlapl will be populated.
+    /// @note The vlapl output is used by process_vlapl_potential() to compute the
+    ///       FD Laplacian kernel contribution to the XC potential.
     extern void tau_xc(
         const std::vector<int> &func_id,
         const double &rho,
         const double &grho,
+        const double &lapl_rho,
         const double &atau,
         double &sxc,
-        double &v1xc,
-        double &v2xc,
-        double &v3xc,
-        const double &hybrid_alpha,
-        const double &hse_omega);
-
+         double &v1xc,
+         double &v2xc,
+         double &v3xc,
+         double &vlapl,
+         const double &hybrid_alpha,
+         const double &hse_omega);
+
+    /// Wrapper for meta-GGA functionals (spin-polarized, two channels).
+    /// Same as tau_xc but handles spin-up/spin-down channels separately.
+    ///
+    /// @param laplup, lapldw  [in]  ∇²ρ_up, ∇²ρ_dn (Laplacian for each spin channel)
+    /// @param vlaplup, vlapldw [out] ∂ε_xc/∂(∇²ρ_up), ∂ε_xc/∂(∇²ρ_dn)
+    ///
+    /// @note See tau_xc() for other parameter descriptions.
     extern void tau_xc_spin(
         const std::vector<int> &func_id,
         double rhoup,
         double rhodw,
         ModuleBase::Vector3<double> gdr1,
         ModuleBase::Vector3<double> gdr2,
+        double laplup,
+        double lapldw,
         double tauup,
         double taudw,
         double &sxc,
@@ -236,10 +299,12 @@ namespace XC_Functional_Libxc
         double &v2xcup,
         double &v2xcdw,
         double &v2xcud,
-        double &v3xcup,
-        double &v3xcdw,
-        const double &hybrid_alpha,
-        const double &hse_omega);
+         double &v3xcup,
+         double &v3xcdw,
+         double &vlaplup,
+         double &vlapldw,
+         const double &hybrid_alpha,
+         const double &hse_omega);
 
 } // namespace XC_Functional_Libxc
 

source/source_hamilt/module_xc/libxc_mgga_wrap.cpp

@@ -17,19 +17,20 @@ void XC_Functional_Libxc::tau_xc(
     const std::vector<int>& func_id,
     const double& rho,
     const double& grho,
+    const double& lapl_rho,
     const double& atau,
     double& sxc,
     double& v1xc,
     double& v2xc,
     double& v3xc,
+    double& vlapl,
     const double& hybrid_alpha,
     const double& hse_omega)
 {
     double s = 0.0;
     double v1 = 0.0;
     double v2 = 0.0;
     double v3 = 0.0;
-    double lapl_rho = grho;
     double vlapl_rho = 0.0;
     std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(
         /* func_id = */ func_id,
@@ -41,6 +42,7 @@ void XC_Functional_Libxc::tau_xc(
     v1xc = 0.0;
     v2xc = 0.0;
     v3xc = 0.0;
+    vlapl = 0.0;
 
     for (xc_func_type& func : funcs)
     {
@@ -52,12 +54,14 @@ void XC_Functional_Libxc::tau_xc(
             v1 *= (1.0 - hybrid_alpha);
             v2 *= (1.0 - hybrid_alpha);
             v3 *= (1.0 - hybrid_alpha);
+            vlapl_rho *= (1.0 - hybrid_alpha);
         }
 #endif
         sxc += s * rho;
         v2xc += v2 * 2.0;
         v1xc += v1;
         v3xc += v3;
+        vlapl += vlapl_rho;
     }
     XC_Functional_Libxc::finish_func(funcs);
 
@@ -71,6 +75,8 @@ void XC_Functional_Libxc::tau_xc_spin(
     double rhodw,
     ModuleBase::Vector3<double> gdr1,
     ModuleBase::Vector3<double> gdr2,
+    double laplup,
+    double lapldw,
     double tauup,
     double taudw,
     double& sxc,
@@ -81,6 +87,8 @@ void XC_Functional_Libxc::tau_xc_spin(
     double& v2xcud,
     double& v3xcup,
     double& v3xcdw,
+    double& vlaplup,
+    double& vlapldw,
     const double& hybrid_alpha,
     const double& hse_omega)
 {
@@ -92,10 +100,13 @@ void XC_Functional_Libxc::tau_xc_spin(
     v2xcud = 0.0;
     v3xcup = 0.0;
     v3xcdw = 0.0;
+    vlaplup = 0.0;
+    vlapldw = 0.0;
 
     const std::array<double, 2> rho = {rhoup, rhodw};
     const std::array<double, 3> grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()};
     const std::array<double, 2> tau = {tauup, taudw};
+    const std::array<double, 2> lapl = {laplup, lapldw};
 
     std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(
         /* func_id = */ func_id,
@@ -125,12 +136,10 @@ void XC_Functional_Libxc::tau_xc_spin(
             double s = 0.0;
             std::array<double, 2> v1xc = {0.0, 0.0};
             std::array<double, 2> v3xc = {0.0, 0.0};
-            std::array<double, 2> lapl = {0.0, 0.0};
             std::array<double, 2> vlapl = {0.0, 0.0};
             std::array<double, 3> v2xc = {0.0, 0.0, 0.0};
-            // call Libxc function: xc_mgga_exc_vxc
-            xc_mgga_exc_vxc(&func, 1, rho.data(), grho.data(), lapl.data(), tau.data(), &s, 
-			    v1xc.data(), v2xc.data(), vlapl.data(), v3xc.data());
+            xc_mgga_exc_vxc(&func, 1, rho.data(), grho.data(), lapl.data(), tau.data(), &s,
+                v1xc.data(), v2xc.data(), vlapl.data(), v3xc.data());
 
 #ifdef __EXX
             if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
@@ -143,6 +152,8 @@ void XC_Functional_Libxc::tau_xc_spin(
                 v2xc[2] *= (1.0 - hybrid_alpha);
                 v3xc[0] *= (1.0 - hybrid_alpha);
                 v3xc[1] *= (1.0 - hybrid_alpha);
+                vlapl[0] *= (1.0 - hybrid_alpha);
+                vlapl[1] *= (1.0 - hybrid_alpha);
             }
 #endif
 
@@ -154,6 +165,8 @@ void XC_Functional_Libxc::tau_xc_spin(
             v2xcdw += 2.0 * v2xc[2] * sgn[1];
             v3xcup += v3xc[0] * sgn[0];
             v3xcdw += v3xc[1] * sgn[1];
+            vlaplup += vlapl[0] * sgn[0];
+            vlapldw += vlapl[1] * sgn[1];
         }
     }