fix(xc): implement proper Laplacian calculation for SCAN-L functional

- Add laplacian_rho() function to calculate Laplacian in reciprocal space
- Update tau_xc() and tau_xc_spin() to accept laplacian parameter
- Fix v_xc_meta() batch path which used sigma as laplacian
- Add laplacian computation in gradcorr() for mGGA functionals
- Add unit test verifying laplacian parameter affects XC energy

Co-authored-by: monkeycode-ai <monkeycode-ai@chaitin.com>
Co-authored-by: monkeycode-ai <monkeycode-ai@chaitin.com>

Author: mintleaf84

source/source_hamilt/module_xc/libxc_abacus.h

@@ -198,6 +198,7 @@ namespace XC_Functional_Libxc
         const std::vector<int> &func_id,
         const double &rho,
         const double &grho,
+        const double &lapl_rho,
         const double &atau,
         double &sxc,
         double &v1xc,
@@ -211,6 +212,8 @@ namespace XC_Functional_Libxc
         double rhodw,
         ModuleBase::Vector3<double> gdr1,
         ModuleBase::Vector3<double> gdr2,
+        double laplup,
+        double lapldw,
         double tauup,
         double taudw,
         double &sxc,

source/source_hamilt/module_xc/libxc_mgga_wrap.cpp

@@ -17,6 +17,7 @@ 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,
@@ -28,7 +29,6 @@ void XC_Functional_Libxc::tau_xc(
     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,
@@ -68,6 +68,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,
@@ -119,7 +121,7 @@ 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> lapl = {laplup, lapldw};
             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

source/source_hamilt/module_xc/libxc_pot.cpp

@@ -12,6 +12,7 @@
 #include <xc.h>
 
 #include <vector>
+#include <complex>
 
 std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		// Peize Lin update for nspin==4 at 2023.01.14
         const std::vector<int> &func_id,
@@ -237,6 +238,27 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
         = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr);
     const std::vector<double> sigma = XC_Functional_Libxc::convert_sigma(gdr);
 
+    // compute laplacian for mGGA functionals
+    std::vector<double> lapl(nrxx * nspin, 0.0);
+    {
+        std::vector<std::complex<double>> rhog_tmp(chr->rhopw->npw);
+        std::vector<double> rhor(nrxx);
+        for( int is=0; is<nspin; ++is )
+        {
+            for( int ir=0; ir<nrxx; ++ir )
+            {
+                rhor[ir] = rho[ir*nspin+is];
+            }
+            chr->rhopw->real2recip(rhor.data(), rhog_tmp.data());
+            std::vector<double> lapl_spin(nrxx);
+            XC_Functional::laplacian_rho(rhog_tmp.data(), lapl_spin.data(), chr->rhopw, tpiba);
+            for( int ir=0; ir<nrxx; ++ir )
+            {
+                lapl[ir*nspin+is] = lapl_spin[ir];
+            }
+        }
+    }
+
     //converting kin_r
     std::vector<double> kin_r;
     kin_r.resize(nrxx*nspin);
@@ -315,7 +337,7 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
                 nrxx_thread,
                 rho.data() + ir_start * nspin,
                 sigma.data() + ir_start * ((1==nspin)?1:3),
-                sigma.data() + ir_start * ((1==nspin)?1:3),
+                lapl.data() + ir_start * nspin,
                 kin_r.data() + ir_start * nspin,
                 exc.data() + ir_start,
                 vrho.data() + ir_start * nspin,

source/source_hamilt/module_xc/test/CMakeLists.txt

@@ -82,4 +82,17 @@ AddTest(
     ../../../source_base/module_fft/fft_bundle.cpp
     ../../../source_base/module_fft/fft_cpu.cpp
     ${FFT_SRC}
+)
+
+AddTest(
+  TARGET MODULE_HAMILT_XCTest_SCANL_LAPL
+  LIBS parameter MPI::MPI_CXX Libxc::xc
+  SOURCES test_xc6.cpp ../xc_functional.cpp ../xc_lda_wrap.cpp
+    ../xc_gga_wrap.cpp
+    ../libxc_setup.cpp
+    ../libxc_lda_wrap.cpp
+    ../libxc_gga_wrap.cpp
+    ../libxc_mgga_wrap.cpp
+    ../xc_gga_corr.cpp ../xc_lda_corr.cpp
+    ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp
 )

source/source_hamilt/module_xc/test/test_xc4.cpp

@@ -47,7 +47,8 @@ class XCTest_SCAN : public XCTest
             {
                 double e,v,v1,v2,v3;
                 double hybrid_alpha = 0.0;
-                XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3,hybrid_alpha);
+                double lapl_rho = grho[i];
+                XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],lapl_rho,tau[i],e,v1,v2,v3,hybrid_alpha);
                 e_.push_back(e);
                 v1_.push_back(v1);
                 v2_.push_back(v2);

source/source_hamilt/module_xc/test/test_xc6.cpp

@@ -0,0 +1,84 @@
+#include "../xc_functional.h"
+#include "../libxc_abacus.h"
+#include "gtest/gtest.h"
+#include "xctest.h"
+#include "../exx_info.h"
+#include <cmath>
+#include <iostream>
+#include <iomanip>
+
+namespace ModuleBase
+{
+    void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);}
+    void TITLE(const std::string &class_function_name,bool disable){};
+    void TITLE(const std::string &class_name,const std::string &function_name,bool disable){};
+}
+
+namespace GlobalV
+{
+    std::string BASIS_TYPE = "";
+    bool CAL_STRESS = false;
+    int CAL_FORCE = 0;
+    int NSPIN = 1;
+}
+
+namespace GlobalC
+{
+    Exx_Info exx_info;
+}
+
+class XCTest_SCANL_Laplacian : public XCTest
+{
+    protected:
+        double e_base, v1_base, v2_base, v3_base;
+        double e_modified, v1_modified, v2_modified, v3_modified;
+        double e_scaled, v1_scaled, v2_scaled, v3_scaled;
+
+        void SetUp()
+        {
+            XC_Functional::set_xc_type("SCAN");
+
+            const double rho = 0.17E+01;
+            const double grho = 0.81E-11;
+            const double tau = 0.02403590412;
+            const double lapl_base = 0.15E+01;
+            double hybrid_alpha = 0.0;
+
+            XC_Functional_Libxc::tau_xc(
+                XC_Functional::get_func_id(),
+                rho, grho, lapl_base, tau,
+                e_base, v1_base, v2_base, v3_base, hybrid_alpha
+            );
+
+            XC_Functional_Libxc::tau_xc(
+                XC_Functional::get_func_id(),
+                rho, grho, lapl_base + 1.0, tau,
+                e_modified, v1_modified, v2_modified, v3_modified, hybrid_alpha
+            );
+
+            XC_Functional_Libxc::tau_xc(
+                XC_Functional::get_func_id(),
+                rho, grho, 2.0 * lapl_base, tau,
+                e_scaled, v1_scaled, v2_scaled, v3_scaled, hybrid_alpha
+            );
+        }
+};
+
+TEST_F(XCTest_SCANL_Laplacian, laplacian_affects_energy)
+{
+    EXPECT_NE(e_base, e_modified);
+    EXPECT_NE(e_base, e_scaled);
+
+    std::cout << std::scientific << std::setprecision(15);
+    std::cout << "\n=== SCAN Laplacian Sensitivity Test ===" << std::endl;
+    std::cout << "Base Laplacian:    " << 0.15E+01 << std::endl;
+    std::cout << "  E_xc =           " << e_base << std::endl;
+    std::cout << "  dE/dlapl ≈       " << (e_modified - e_base) / 1.0 << std::endl;
+    std::cout << "Modified Laplacian:" << 0.15E+01 + 1.0 << std::endl;
+    std::cout << "  E_xc =           " << e_modified << std::endl;
+    std::cout << "  Delta E =        " << e_modified - e_base << std::endl;
+    std::cout << "Scaled Laplacian:  " << 2.0 * 0.15E+01 << std::endl;
+    std::cout << "  E_xc =           " << e_scaled << std::endl;
+    std::cout << "  Delta E =        " << e_scaled - e_base << std::endl;
+    std::cout << "========================================" << std::endl;
+}

source/source_hamilt/module_xc/xc_functional.h

@@ -229,6 +229,12 @@ class XC_Functional
         const ModulePW::PW_Basis* rho_basis,
         const double tpiba);
 
+    static void laplacian_rho(
+        const std::complex<double>* rhog,
+        double* lapl,
+        const ModulePW::PW_Basis* rho_basis,
+        const double tpiba);
+
     static void noncolin_rho(
         double* rhoout1,
         double* rhoout2,

source/source_hamilt/module_xc/xc_grad.cpp

@@ -90,6 +90,8 @@ void XC_Functional::gradcorr(
     double* neg = nullptr;
     double** vsave = nullptr;
     double** vgg = nullptr;
+    double* lapl1 = nullptr;
+    double* lapl2 = nullptr;
 
     // for spin unpolarized case,
     // calculate the gradient of (rho_core+rho) in reciprocal space.
@@ -118,6 +120,12 @@ void XC_Functional::gradcorr(
 
     XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
 
+    if(func_type == 3 || func_type == 5)
+    {
+        lapl1 = new double[rhopw->nrxx];
+        XC_Functional::laplacian_rho(rhogsum1, lapl1, rhopw, ucell->tpiba);
+    }
+
     // for spin polarized case;
     // calculate the gradient of (rho_core+rho) in reciprocal space.
     if(PARAM.inp.nspin==2)
@@ -146,6 +154,12 @@ void XC_Functional::gradcorr(
         }
 
         XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba);
+
+        if(func_type == 3 || func_type == 5)
+        {
+            lapl2 = new double[rhopw->nrxx];
+            XC_Functional::laplacian_rho(rhogsum2, lapl2, rhopw, ucell->tpiba);
+        }
     }
 
     if(PARAM.inp.nspin == 4&&(PARAM.globalv.domag||PARAM.globalv.domag_z))
@@ -222,6 +236,20 @@ void XC_Functional::gradcorr(
 
         XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
         XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba);
+
+        if(func_type == 3 || func_type == 5)
+        {
+            if(lapl1 == nullptr)
+            {
+                lapl1 = new double[rhopw->nrxx];
+            }
+            XC_Functional::laplacian_rho(rhogsum1, lapl1, rhopw, ucell->tpiba);
+            if(lapl2 == nullptr)
+            {
+                lapl2 = new double[rhopw->nrxx];
+            }
+            XC_Functional::laplacian_rho(rhogsum2, lapl2, rhopw, ucell->tpiba);
+        }
     }
 
     const double epsr = 1.0e-6;
@@ -297,7 +325,8 @@ void XC_Functional::gradcorr(
 #ifdef __EXX
                             hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha;
 #endif
-                            XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc, hybrid_alpha);
+                            double lapl_val = (lapl1 != nullptr) ? lapl1[ir] : 0.0;
+                            XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, lapl_val, atau, sxc, v1xc, v2xc, v3xc, hybrid_alpha);
                         }
                         else
                         {
@@ -366,10 +395,12 @@ void XC_Functional::gradcorr(
 #ifdef __EXX
                         hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha;
 #endif
+                        double laplup_val = (lapl1 != nullptr) ? lapl1[ir] : 0.0;
+                        double lapldw_val = (lapl2 != nullptr) ? lapl2[ir] : 0.0;
                         XC_Functional_Libxc::tau_xc_spin(
                             func_id,
                             rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir],
-                            atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, hybrid_alpha);
+                            laplup_val, lapldw_val, atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, hybrid_alpha);
                     }
                     else
                     {
@@ -653,14 +684,18 @@ void XC_Functional::gradcorr(
         delete[] rhotmp2;
         delete[] rhogsum2;
         delete[] gdr2;
+        delete[] lapl2;
         if(!is_stress)
         {
             delete[] h2;
         }
+        delete[] lapl1;
     }
-    if(PARAM.inp.nspin == 4 && (PARAM.globalv.domag||PARAM.globalv.domag_z))
+    else if(PARAM.inp.nspin == 4 && (PARAM.globalv.domag||PARAM.globalv.domag_z))
     {
         delete[] neg;
+        delete[] lapl1;
+        delete[] lapl2;
         if(!is_stress)
         {
             for(int i=0; i<nspin0; i++)
@@ -679,6 +714,10 @@ void XC_Functional::gradcorr(
         delete[] rhogsum2;
         delete[] gdr2;
     }
+    else
+    {
+        delete[] lapl1;
+    }
 
     return;
 }
@@ -816,11 +855,42 @@ void XC_Functional::grad_dot(
         dh[ir] = aux[ir].real() * tpiba;
     }
 
-    delete[] aux;
+delete[] aux;
     delete[] gaux;
     return;
 }
 
+
+void XC_Functional::laplacian_rho(
+    const std::complex<double>* rhog,
+    double* lapl,
+    const ModulePW::PW_Basis* rho_basis,
+    const double tpiba)
+{
+    std::complex<double>* lapl_tmp = new std::complex<double>[rho_basis->nmaxgr];
+
+    for(int ir=0; ir<rho_basis->nrxx; ir++)
+    {
+        lapl[ir] = 0.0;
+    }
+
+    for(int i=0; i<3; i++)
+    {
+        for(int ig=0; ig<rho_basis->npw; ig++)
+        {
+            lapl_tmp[ig] = -rhog[ig] * rho_basis->gcar[ig][i] * rho_basis->gcar[ig][i];
+        }
+        rho_basis->recip2real(lapl_tmp, lapl_tmp);
+        for(int ir=0; ir<rho_basis->nrxx; ir++)
+        {
+            lapl[ir] += lapl_tmp[ir].real() * tpiba * tpiba;
+        }
+    }
+
+    delete[] lapl_tmp;
+}
+
+
 void XC_Functional::noncolin_rho(
     double *rhoout1,
     double *rhoout2,