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

source/source_estate/module_pot/pot_xc.cpp

@@ -1,8 +1,12 @@
 #include "pot_xc.h"
 
 #include "source_base/timer.h"
+#include "source_base/constants.h"
 #include "source_hamilt/module_xc/xc_functional.h"
 
+#include <vector>
+#include <complex>
+
 #ifdef USE_LIBXC
 #include "source_hamilt/module_xc/libxc_abacus.h"
 #endif
@@ -23,12 +27,65 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module
     if (XC_Functional::get_ked_flag())
     {
 #ifdef USE_LIBXC
-        const std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix> etxc_vtxc_v
+        const std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix, ModuleBase::matrix> etxc_vtxc_v
             = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), nrxx_current, ucell->omega, ucell->tpiba, chg);
         *(this->etxc_) = std::get<0>(etxc_vtxc_v);
         *(this->vtxc_) = std::get<1>(etxc_vtxc_v);
         v_eff += std::get<2>(etxc_vtxc_v);
         *(this->vofk) = std::get<3>(etxc_vtxc_v);
+
+        // Apply Laplacian potential correction using FD kernel
+        // FD kernel matches spectral -|G|² at low G but stays bounded at high G,
+        // preventing |G|² amplification that causes SCF divergence
+        const ModuleBase::matrix& voflapl = std::get<4>(etxc_vtxc_v);
+        const int ng = chg->rhopw->npw;
+        const int nrxx = chg->rhopw->nrxx;
+        const double tpiba2 = ucell->tpiba * ucell->tpiba;
+        const double twopi = 2.0 * M_PI;
+        const int nx = chg->rhopw->nx;
+        const int ny = chg->rhopw->ny;
+        const int nz = chg->rhopw->nz;
+
+        // Pre-compute FD gg for all G vectors
+        std::vector<double> gg_fd(ng);
+        for(int ig = 0; ig < ng; ig++)
+        {
+            const int m0 = chg->rhopw->gdirect[ig].x;
+            const int m1 = chg->rhopw->gdirect[ig].y;
+            const int m2 = chg->rhopw->gdirect[ig].z;
+
+            const double fd00 = 2.0 * nx * nx * (1.0 - std::cos(twopi * m0 / nx));
+            const double fd11 = 2.0 * ny * ny * (1.0 - std::cos(twopi * m1 / ny));
+            const double fd22 = 2.0 * nz * nz * (1.0 - std::cos(twopi * m2 / nz));
+            const double fd01 = nx * ny * std::sin(twopi * m0 / nx) * std::sin(twopi * m1 / ny);
+            const double fd02 = nx * nz * std::sin(twopi * m0 / nx) * std::sin(twopi * m2 / nz);
+            const double fd12 = ny * nz * std::sin(twopi * m1 / ny) * std::sin(twopi * m2 / nz);
+
+            gg_fd[ig] = (chg->rhopw->GGT.e11 * fd00 + chg->rhopw->GGT.e22 * fd11 + chg->rhopw->GGT.e33 * fd22
+                       + 2.0 * chg->rhopw->GGT.e12 * fd01 + 2.0 * chg->rhopw->GGT.e13 * fd02
+                       + 2.0 * chg->rhopw->GGT.e23 * fd12) / (twopi * twopi);
+        }
+
+        std::vector<std::complex<double>> lapl_tmp(chg->rhopw->nmaxgr);
+        for(int is = 0; is < voflapl.nr; is++)
+        {
+            for(int ir = 0; ir < nrxx; ir++)
+                lapl_tmp[ir] = std::complex<double>(voflapl(is, ir), 0.0);
+            for(int ig = ng; ig < chg->rhopw->nmaxgr; ig++)
+                lapl_tmp[ig] = std::complex<double>(0.0, 0.0);
+            chg->rhopw->real2recip(lapl_tmp.data(), lapl_tmp.data());
+            for(int ig = 0; ig < ng; ig++)
+            {
+                lapl_tmp[ig] *= -gg_fd[ig] * tpiba2;
+            }
+            chg->rhopw->recip2real(lapl_tmp.data(), lapl_tmp.data());
+            for(int ir = 0; ir < nrxx; ir++)
+            {
+                double vlapl_corr = ModuleBase::e2 * lapl_tmp[ir].real();
+                v_eff(is, ir) += vlapl_corr;
+                *(this->vtxc_) += vlapl_corr * chg->rho[is][ir] * ucell->omega / chg->rhopw->nxyz;
+            }
+        }
 #else
         ModuleBase::WARNING_QUIT("v_of_rho", "to use mGGA, compile with LIBXC");
 #endif

source/source_hamilt/module_xc/libxc_abacus.h

@@ -63,7 +63,7 @@ namespace XC_Functional_Libxc
         const std::map<int, double>* scaling_factor = nullptr); // added by jghan, 2024-10-10
 
     // for mGGA functional
-    extern std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix> v_xc_meta(
+    extern std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix, ModuleBase::matrix> v_xc_meta(
         const std::vector<int> &func_id,
         const int &nrxx, // number of real-space grid
         const double &omega, // volume of cell
@@ -198,11 +198,13 @@ 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,
         double &v2xc,
         double &v3xc,
+        double &vlaplxc,
         const double &hybrid_alpha);
 
     extern void tau_xc_spin(
@@ -211,6 +213,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,
@@ -221,6 +225,8 @@ namespace XC_Functional_Libxc
         double &v2xcud,
         double &v3xcup,
         double &v3xcdw,
+        double &vlaplxcup,
+        double &vlaplxcdw,
         const double &hybrid_alpha);
 
 } // namespace XC_Functional_Libxc

source/source_hamilt/module_xc/libxc_mgga_wrap.cpp

@@ -17,18 +17,19 @@ 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& vlaplxc,
     const double& hybrid_alpha)
 {
     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,
@@ -38,6 +39,7 @@ void XC_Functional_Libxc::tau_xc(
     v1xc = 0.0;
     v2xc = 0.0;
     v3xc = 0.0;
+    vlaplxc = 0.0;
 
     for (xc_func_type& func : funcs)
     {
@@ -49,12 +51,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;
+        vlaplxc += vlapl_rho;
     }
     XC_Functional_Libxc::finish_func(funcs);
 
@@ -68,6 +72,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,
@@ -78,6 +84,8 @@ void XC_Functional_Libxc::tau_xc_spin(
     double& v2xcud,
     double& v3xcup,
     double& v3xcdw,
+    double& vlaplxcup,
+    double& vlaplxcdw,
     const double& hybrid_alpha)
 {
     sxc = 0.0;
@@ -88,6 +96,8 @@ void XC_Functional_Libxc::tau_xc_spin(
     v2xcud = 0.0;
     v3xcup = 0.0;
     v3xcdw = 0.0;
+    vlaplxcup = 0.0;
+    vlaplxcdw = 0.0;
 
     const std::array<double, 2> rho = {rhoup, rhodw};
     const std::array<double, 3> grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()};
@@ -119,7 +129,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
@@ -137,6 +147,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
 
@@ -148,6 +160,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];
+            vlaplxcup += vlapl[0] * sgn[0];
+            vlaplxcdw += vlapl[1] * sgn[1];
         }
     }
 

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,
@@ -202,7 +203,7 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
 //XC_POLARIZED, XC_UNPOLARIZED: internal flags used in LIBXC, denote the polarized(nspin=1) or unpolarized(nspin=2) calculations, definition can be found in xc.h from LIBXC
 
 // [etxc, vtxc, v, vofk] = XC_Functional::v_xc(...)
-std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Libxc::v_xc_meta(
+std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Libxc::v_xc_meta(
     const std::vector<int> &func_id,
     const int &nrxx, // number of real-space grid
     const double &omega, // volume of cell
@@ -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);
@@ -255,6 +277,7 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
     std::vector<double> vrho   ( nrxx * nspin            );
     std::vector<double> vsigma ( nrxx * ((1==nspin)?1:3) );
     std::vector<double> vtau   ( nrxx * nspin            );
+    ModuleBase::matrix voflapl(PARAM.inp.nspin,nrxx);
     std::vector<double> vlapl  ( nrxx * nspin            );
 
     constexpr double rho_th  = 1e-8;
@@ -315,7 +338,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,
@@ -440,9 +463,11 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
                 if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
                 {
                     vtau[ir*nspin+is] *= (1.0 - XC_Functional::get_hybrid_alpha());
+                    vlapl[ir*nspin+is] *= (1.0 - XC_Functional::get_hybrid_alpha());
                 }
 #endif
                 vofk(is,ir) += vtau[ir*nspin+is]  * sgn[ir*nspin+is];
+                voflapl(is,ir) += vlapl[ir*nspin+is] * sgn[ir*nspin+is];
             }
         }
     }
@@ -461,7 +486,7 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
     XC_Functional_Libxc::finish_func(funcs);
 
     ModuleBase::timer::end("XC_Functional_Libxc","v_xc_meta");
-    return std::make_tuple( etxc, vtxc, std::move(v), std::move(vofk) );
+    return std::make_tuple( etxc, vtxc, std::move(v), std::move(vofk), std::move(voflapl) );
 }
 
 #endif

source/source_hamilt/module_xc/libxc_setup.cpp

@@ -180,6 +180,23 @@ XC_Functional_Libxc::set_xc_type_libxc(const std::string& xc_func_in)
         ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc", message);
     }
 
+    // warn if any functional needs Laplacian of density
+    {
+        std::vector<xc_func_type> tmp_funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED);
+        for (auto& f : tmp_funcs)
+        {
+            if (f.info->flags & XC_FLAGS_NEEDS_LAPLACIAN)
+            {
+                std::cout << " WARNING: XC functional \"" << f.info->name
+                          << "\" requires Laplacian of density (nabla^2 rho)."
+                          << " This may require a higher energy cutoff for numerical stability."
+                          << std::endl;
+                break;
+            }
+        }
+        XC_Functional_Libxc::finish_func(tmp_funcs);
+    }
+
     // return
     return std::make_pair(func_type, func_id);
 }

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

@@ -45,9 +45,10 @@ class XCTest_SCAN : public XCTest
 
             for(int i=0;i<5;i++)
             {
-                double e,v,v1,v2,v3;
+                double e,v,v1,v2,v3,vlapl;
                 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 = 0.0; // SCAN does not use Laplacian; value has no effect on result
+                XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],lapl_rho,tau[i],e,v1,v2,v3,vlapl,hybrid_alpha);
                 e_.push_back(e);
                 v1_.push_back(v1);
                 v2_.push_back(v2);

source/source_hamilt/module_xc/test/test_xc5.cpp

@@ -282,12 +282,13 @@ class XCTest_VXC_meta : public XCTest
             XC_Functional::set_xc_type("SCAN");
 
             PARAM.input.nspin = 1;
-            std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix> etxc_vtxc_v
+            std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix, ModuleBase::matrix> etxc_vtxc_v
                 = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr);
             et1 = std::get<0>(etxc_vtxc_v);
             vt1 = std::get<1>(etxc_vtxc_v);
             v1  = std::get<2>(etxc_vtxc_v);
             vtau1 = std::get<3>(etxc_vtxc_v);
+            ModuleBase::matrix vlapl1 = std::get<4>(etxc_vtxc_v);
 
             PARAM.input.nspin = 2;
             etxc_vtxc_v
@@ -296,6 +297,7 @@ class XCTest_VXC_meta : public XCTest
             vt2 = std::get<1>(etxc_vtxc_v);
             v2  = std::get<2>(etxc_vtxc_v);
             vtau2 = std::get<3>(etxc_vtxc_v);
+            ModuleBase::matrix vlapl2 = std::get<4>(etxc_vtxc_v);
         }
 };