fix(xc): implement proper Laplacian calculation for SCAN-L functional (Useful Information for discussing implementation of SCAN-L functional)

source/source_estate/module_pot/pot_xc.cpp

@@ -1,6 +1,7 @@
 #include "pot_xc.h"
 
 #include "source_base/timer.h"
+#include "source_base/constants.h"
 #include "source_hamilt/module_xc/xc_functional.h"
 #include "source_io/module_parameter/parameter.h"
 

source/source_hamilt/module_xc/libxc_abacus.h

@@ -105,6 +105,16 @@ 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,
+        const bool need_laplacian = true);
+
     // converting grho (abacus=>libxc)
     extern std::vector<double> convert_sigma(
         const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr);
@@ -214,11 +224,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,
         const double &hse_omega);
 
@@ -228,6 +240,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,
@@ -238,6 +252,8 @@ namespace XC_Functional_Libxc
         double &v2xcud,
         double &v3xcup,
         double &v3xcdw,
+        double &vlaplxcup,
+        double &vlaplxcdw,
         const double &hybrid_alpha,
         const double &hse_omega);
 

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& vlaplxc,
     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;
+    vlaplxc = 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;
+        vlaplxc += 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& vlaplxcup,
+    double& vlaplxcdw,
     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;
+    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()};
     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,11 @@ 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, 2> vlapl_out = {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());
+			    v1xc.data(), v2xc.data(), vlapl_out.data(), v3xc.data());
 
 #ifdef __EXX
             if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
@@ -143,6 +153,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_out[0] *= (1.0 - hybrid_alpha);
+                vlapl_out[1] *= (1.0 - hybrid_alpha);
             }
 #endif
 
@@ -154,6 +166,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_out[0] * sgn[0];
+            vlaplxcdw += vlapl_out[1] * sgn[1];
         }
     }
 

source/source_hamilt/module_xc/libxc_pot.cpp

@@ -15,6 +15,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,
@@ -225,13 +226,12 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
     ModuleBase::TITLE("XC_Functional_Libxc","v_xc_meta");
     ModuleBase::timer::start("XC_Functional_Libxc","v_xc_meta");
 
-    double e2 = 2.0;
-
     //output of the subroutine
     double etxc = 0.0;
     double vtxc = 0.0;
     ModuleBase::matrix v(nspin,nrxx);
     ModuleBase::matrix vofk(nspin,nrxx);
+    ModuleBase::matrix voflapl(nspin,nrxx);
 
     //----------------------------------------------------------
     // xc_func_type is defined in Libxc package
@@ -246,8 +246,10 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
         /* hse_omega = */ hse_omega);
 
     const std::vector<double> rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr);
-    const std::vector<std::vector<ModuleBase::Vector3<double>>> gdr
-        = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr);
+    const bool need_laplacian = XC_Functional::get_need_laplacian();
+    std::vector<std::vector<ModuleBase::Vector3<double>>> gdr;
+    std::vector<double> lapl;
+    XC_Functional_Libxc::cal_gdr_and_lapl(nspin, nrxx, rho, tpiba, chr, gdr, lapl, need_laplacian);
     const std::vector<double> sigma = XC_Functional_Libxc::convert_sigma(gdr);
 
     //converting kin_r
@@ -328,7 +330,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,
@@ -441,7 +443,7 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
         }
         vtxc -= rvtxc;
 
-        //process vtau
+        //process vtau and vlapl
 #ifdef _OPENMP
 #pragma omp parallel for collapse(2) schedule(static, 1024)
 #endif
@@ -453,9 +455,41 @@ 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];
+            }
+        }
+}
+
+    // v_xc += nabla^2(vlapl) where vlapl = d(rho*eps_xc)/d(nabla^2 rho)
+    if (need_laplacian)
+    {
+        const int ng = chr->rhopw->npw;
+        const double tpiba2 = tpiba * tpiba;
+        std::vector<std::complex<double>> lapl_tmp(chr->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 < chr->rhopw->nmaxgr; ig++)
+                lapl_tmp[ig] = std::complex<double>(0.0, 0.0);
+            chr->rhopw->real2recip(lapl_tmp.data(), lapl_tmp.data());
+            for(int ig = 0; ig < ng; ig++)
+            {
+                double g2 = 0.0;
+                for(int i = 0; i < 3; i++)
+                    g2 += chr->rhopw->gcar[ig][i] * chr->rhopw->gcar[ig][i];
+                lapl_tmp[ig] *= -g2 * tpiba2;
+            }
+            chr->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(is, ir) += vlapl_corr;
+                vtxc += vlapl_corr * chr->rho[is][ir];
             }
         }
     }

source/source_hamilt/module_xc/libxc_setup.cpp

@@ -17,21 +17,10 @@
 
 bool not_supported_xc_with_laplacian(const std::string& xc_func_in)
 {
-    // see Pyscf: https://github.com/pyscf/pyscf/blob/master/pyscf/dft/libxc.py#L1062
-    // ABACUS issue: https://github.com/deepmodeling/abacus-develop/issues/5372
-    const std::vector<std::string> not_supported = {
-        "MGGA_XC_CC06",
-        "MGGA_C_CS",
-        "MGGA_X_BR89",
-        "MGGA_X_MK00"
-    };
-    for (const std::string& s : not_supported)
-    {
-        if (xc_func_in.find(s) != std::string::npos)
-        {
-            return true;
-        }
-    }
+    // Laplacian of density is now supported for meta-GGA functionals.
+    // The following functionals were previously blocked but are now supported:
+    //   MGGA_XC_CC06, MGGA_C_CS, MGGA_X_BR89, MGGA_X_MK00
+    // Ensure PW stress path handles vlapl contribution correctly if using PW basis.
     return false;
 }
 
@@ -179,6 +168,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, 0.0, 0.0);
+        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/libxc_tools.cpp

@@ -88,6 +88,37 @@ XC_Functional_Libxc::cal_gdr(
 	return gdr;
 }
 
+void XC_Functional_Libxc::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,
+	const bool need_laplacian)
+{
+	gdr.resize(nspin);
+	lapl.assign(nrxx * nspin, 0.0);
+	for( int is=0; is!=nspin; ++is )
+	{
+		std::vector<double> rhor(nrxx);
+		for(std::size_t ir=0; ir<nrxx; ++ir)
+			rhor[ir] = rho[ir*nspin+is];
+		std::vector<std::complex<double>> rhog(chr->rhopw->npw);
+		chr->rhopw->real2recip(rhor.data(), rhog.data());
+		gdr[is].resize(nrxx);
+		XC_Functional::grad_rho(rhog.data(), gdr[is].data(), chr->rhopw, tpiba);
+		if (need_laplacian)
+		{
+			std::vector<double> lapl_spin(nrxx);
+			XC_Functional::laplacian_rho(rhog.data(), lapl_spin.data(), chr->rhopw, tpiba);
+			for(std::size_t ir=0; ir<nrxx; ++ir)
+				lapl[ir*nspin+is] = lapl_spin[ir];
+		}
+	}
+}
+
 // converting grho (abacus=>libxc)
 std::vector<double> XC_Functional_Libxc::convert_sigma(
 	const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr)

source/source_hamilt/module_xc/test/CMakeLists.txt

@@ -83,3 +83,43 @@ AddTest(
     ../../../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
+)
+
+if (USE_CUDA)
+list(APPEND FFT_SRC ../../../source_base/module_fft/fft_cuda.cpp)
+endif()
+if (USE_ROCM)
+list(APPEND FFT_SRC ../../../source_base/module_fft/fft_rocm.cpp)
+endif()
+AddTest(
+  TARGET MODULE_HAMILT_XCTest_LAPL
+  LIBS parameter MPI::MPI_CXX Libxc::xc psi device container
+  SOURCES test_xc7.cpp ../xc_grad.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_base/matrix.cpp
+  ../../../source_base/memory_recorder.cpp
+  ../../../source_base/libm/branred.cpp
+  ../../../source_base/libm/sincos.cpp
+  ../../../source_base/module_external/blas_connector_base.cpp ../../../source_base/module_external/blas_connector_vector.cpp ../../../source_base/module_external/blas_connector_matrix.cpp
+  ../../../source_base/module_fft/fft_bundle.cpp
+  ../../../source_base/module_fft/fft_cpu.cpp
+  ${FFT_SRC}
+)