Fix(module_xc): compute real density laplacian for meta-GGA functionals via reciprocal space

SCANL, R2SCANL, and REVSCAN require the density laplacian (∇²ρ) but
v_xc_meta was passing |∇ρ|² (sigma) instead. This caused ~150 kbar
stress error for Na with SCANL vs SCAN.

Laplacian computation method:
  ∇²ρ(r) = -∑_G |G|² ρ(G) e^{iGr}
  - ρ(r) is Fourier-transformed to ρ(G) via real2recip()
  - Each G-component is multiplied by -|G|² (i.e. -gg[ig] * tpiba²)
  - Inverse FFT via recip2real() yields ∇²ρ(r) on the real-space grid
  - Only the real part is retained (imaginary part should be ~0 for real ρ)
  This is exact within the plane-wave basis, not numerical differentiation.

Spin handling:
  - cal_lapl(nspin, nrxx, rho, tpiba2, chr) computes ∇²ρ for each spin
    channel: rho is interleaved [rho_up[0], rho_dw[0], rho_up[1], ...]
    and output lapl is [lapl_up[0], lapl_dw[0], lapl_up[1], ...]
  - nspin=1: single channel (spin-unpolarized), loop runs once for is=0
  - nspin=2: separate laplacian for spin-up and spin-down channels,
    each computed independently from its own real-space density
  - nspin=4 (noncollinear): effective up/down densities are constructed
    via noncolin_rho() diagonalization of the spin density matrix,
    then ∇²ρ is computed for each effective channel

Changes:
  - Add cal_lapl() in xc_functional_libxc_tools.cpp as reusable function
    for reciprocal-space Laplacian computation per spin channel
  - Fix v_xc_meta to call cal_lapl() and pass real laplacian to LibXC;
    handle vlapl output by computing ∇²(vlapl) via same G-space method
  - Fix tau_xc wrapper — added lapl_rho parameter, removed grho hack
  - Fix tau_xc_spin wrapper — initialize lapl array to zero instead of
    leaving uninitialized
  - Fix gradcorr stress/force paths — replace ~50 lines of duplicated
    inline Laplacian code with cal_lapl() calls; handles nspin=1/2/4
    spin channel selection consistently with existing grad_rho pattern
  - Update test_xc4.cpp and test CMakeLists for new tau_xc signature

Author: linpeize

source/module_esolver/esolver_ks.cpp

@@ -185,6 +185,7 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
     {
         XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
     }
+    GlobalV::ofs_running<<XC_Functional::output_info()<<std::endl;
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
 
     //! 4) setup the charge mixing parameters

source/module_hamilt_general/module_xc/test/CMakeLists.txt

@@ -11,6 +11,7 @@ AddTest(
   TARGET XCTest_HSE
   LIBS parameter MPI::MPI_CXX Libxc::xc  # required by global.h; for details, `remove_definitions(-D__MPI)`.
   SOURCES test_xc1.cpp ../xc_functional.cpp ../xc_functional_libxc.cpp
+    ../xc_functional_libxc_tools.cpp
 )
 
 
@@ -36,6 +37,7 @@ AddTest(
     ../xc_functional_libxc_wrapper_xc.cpp
     ../xc_functional_libxc_wrapper_gcxc.cpp
     ../xc_functional_libxc_wrapper_tauxc.cpp
+    ../xc_functional_libxc_tools.cpp
     ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp
     ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp
     ../../../module_base/matrix.cpp
@@ -57,6 +59,7 @@ AddTest(
     ../xc_functional_libxc_wrapper_xc.cpp
     ../xc_functional_libxc_wrapper_gcxc.cpp
     ../xc_functional_libxc_wrapper_tauxc.cpp
+    ../xc_functional_libxc_tools.cpp
     ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp 
     ../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp 
 )

source/module_hamilt_general/module_xc/test/test_xc4.cpp

@@ -44,7 +44,8 @@ class XCTest_SCAN : public XCTest
             for(int i=0;i<5;i++)
             {
                 double e,v,v1,v2,v3;
-                XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3);
+                // SCAN does not depend on laplacian, pass 0.0
+                XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],0.0,tau[i],e,v1,v2,v3);
                 e_.push_back(e);
                 v1_.push_back(v1);
                 v2_.push_back(v2);

source/module_hamilt_general/module_xc/xc_functional.cpp

@@ -34,24 +34,37 @@ int XC_Functional::get_func_type()
 {
     return func_type;
 }
+
 void XC_Functional::set_xc_first_loop(const UnitCell& ucell)
 {
+    ModuleBase::TITLE("XC_Functional", "set_xc_first_loop");
     /** In the special "two-level" calculation case,
 the first scf iteration only calculate the functional without exact
 exchange. but in "nscf" calculation, there is no need of "two-level"
 method. */
-    if (ucell.atoms[0].ncpp.xc_func == "HF" || ucell.atoms[0].ncpp.xc_func == "HSE" 
-        || ucell.atoms[0].ncpp.xc_func == "PBE0"|| ucell.atoms[0].ncpp.xc_func == "LC_PBE" 
-        || ucell.atoms[0].ncpp.xc_func == "LC_WPBE" || ucell.atoms[0].ncpp.xc_func == "LRC_WPBEH" 
-        || ucell.atoms[0].ncpp.xc_func == "CAM_PBEH") {
+    if (ucell.atoms[0].ncpp.xc_func == "HF" || ucell.atoms[0].ncpp.xc_func == "PBE0" || ucell.atoms[0].ncpp.xc_func == "HSE")
+    {
         XC_Functional::set_xc_type("pbe");
     }
-    else if (ucell.atoms[0].ncpp.xc_func == "SCAN0") {
-        XC_Functional::set_xc_type("scan");
+    else if ( ucell.atoms[0].ncpp.xc_func == "LC_PBE" || ucell.atoms[0].ncpp.xc_func == "LC_WPBE"
+        || ucell.atoms[0].ncpp.xc_func == "LRC_WPBEH" || ucell.atoms[0].ncpp.xc_func == "CAM_PBEH" )
+    {
+        XC_Functional::set_xc_type("pbe");
     }
-    else if (ucell.atoms[0].ncpp.xc_func == "B3LYP") {
+    // added by jghan, 2024-07-07
+    else if ( ucell.atoms[0].ncpp.xc_func == "MULLER" || ucell.atoms[0].ncpp.xc_func == "POWER"
+        || ucell.atoms[0].ncpp.xc_func == "WP22" || ucell.atoms[0].ncpp.xc_func == "CWP22" )
+    {
+        XC_Functional::set_xc_type("pbe");
+    }
+    else if (ucell.atoms[0].ncpp.xc_func == "B3LYP")
+    {
         XC_Functional::set_xc_type("blyp");
     }
+    else if (ucell.atoms[0].ncpp.xc_func == "SCAN0")
+    {
+        XC_Functional::set_xc_type("scan");
+    }
 }
 
 // The setting values of functional id according to the index in LIBXC
@@ -346,4 +359,44 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
     use_libxc = false;
 #endif
 
+}
+
+
+std::string XC_Functional::output_info()
+{
+  #ifdef USE_LIBXC
+    if(use_libxc)
+    {
+        std::stringstream ss;
+        ss<<" Libxc v"<<xc_version_string()<<std::endl;
+        ss<<"\t"<<xc_reference()<<std::endl;
+
+        std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED);
+        for(const auto &func : funcs)
+        {
+            const xc_func_info_type *info = xc_func_get_info(&func);
+            ss<<" XC: "<<xc_func_info_get_name(info)<<std::endl;
+            for(int i=0; i<XC_MAX_REFERENCES; ++i)
+            {
+                const func_reference_type *ref = xc_func_info_get_references(func.info, i);
+                if(ref)
+                    ss<<"\t"<<xc_func_reference_get_ref(ref)<<std::endl;
+            }
+        }
+        XC_Functional_Libxc::finish_func(funcs);
+        return ss.str();
+    }
+    else
+    {
+        std::string s = " XC:\t";
+        for(const auto &id: func_id)
+            s += std::string(xc_functional_get_name(id))+"\t";
+        return s;
+    }
+  #else
+    std::string s = " XC:\t";
+    for(const auto &id: func_id)
+        s += std::to_string(id)+"\t";
+    return s;
+  #endif
 }

source/module_hamilt_general/module_xc/xc_functional.h

@@ -88,6 +88,7 @@ class XC_Functional
 
 	public:
 	static std::vector<int> get_func_id() { return func_id; }
+	static std::string output_info();
 
 //-------------------
 //  xc_functional_wrapper_xc.cpp

source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp

@@ -67,6 +67,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 	ModuleBase::Vector3<double>* gdr2 = nullptr;
 	ModuleBase::Vector3<double>* h1 = nullptr;
 	ModuleBase::Vector3<double>* h2 = nullptr;
+	double* lapl1 = nullptr;
+	double* lapl2 = nullptr;
 	double* neg = nullptr;
 	double** vsave = nullptr;
 	double** vgg = nullptr;
@@ -95,6 +97,27 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 
 	XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
 
+	XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
+
+	// compute laplacian of total density (rho[0] + rho_core) for nspin=1,2
+	// and as the "up" channel equivalent for nspin=4 noncollinear case
+	{
+		std::vector<double> rho_total(rhopw->nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+		for(int ir=0; ir<rhopw->nrxx; ++ir)
+			rho_total[ir] = rhotmp1[ir];
+		std::vector<double> lapl_all = XC_Functional_Libxc::cal_lapl(
+			1, rhopw->nrxx, rho_total, ucell->tpiba2, chr);
+		lapl1 = new double[rhopw->nrxx];
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+		for(int ir=0; ir<rhopw->nrxx; ++ir)
+			lapl1[ir] = lapl_all[ir];
+	}
+
 	// for spin polarized case;
 	// calculate the gradient of (rho_core+rho) in reciprocal space.
 	if(PARAM.inp.nspin==2)
@@ -120,6 +143,24 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 		if(!is_stress) { h2 = new ModuleBase::Vector3<double>[rhopw->nrxx]; }
 
 		XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba);
+
+		// compute laplacian of spin-down density (rho[1] + rho_core)
+		{
+			std::vector<double> rho_dw(rhopw->nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				rho_dw[ir] = rhotmp2[ir];
+			std::vector<double> lapl_all = XC_Functional_Libxc::cal_lapl(
+				1, rhopw->nrxx, rho_dw, ucell->tpiba2, chr);
+			lapl2 = new double[rhopw->nrxx];
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				lapl2[ir] = lapl_all[ir];
+		}
 	}
 
 	if(PARAM.inp.nspin == 4&&(PARAM.globalv.domag||PARAM.globalv.domag_z))
@@ -189,6 +230,40 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 		XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
 		XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba);
 
+		// re-compute laplacians after noncollinear rho update
+		{
+			std::vector<double> rho_up(rhopw->nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				rho_up[ir] = rhotmp1[ir];
+			std::vector<double> lapl_up = XC_Functional_Libxc::cal_lapl(
+				1, rhopw->nrxx, rho_up, ucell->tpiba2, chr);
+			lapl1 = new double[rhopw->nrxx];
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				lapl1[ir] = lapl_up[ir];
+		}
+		{
+			std::vector<double> rho_dw(rhopw->nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				rho_dw[ir] = rhotmp2[ir];
+			std::vector<double> lapl_dw = XC_Functional_Libxc::cal_lapl(
+				1, rhopw->nrxx, rho_dw, ucell->tpiba2, chr);
+			lapl2 = new double[rhopw->nrxx];
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				lapl2[ir] = lapl_dw[ir];
+		}
+
 	}
 
 	const double epsr = 1.0e-6;
@@ -251,7 +326,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 					{
 						double v3xc;
 						double atau = chr->kin_r[0][ir]/2.0;
-						XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc);
+						XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, lapl1[ir], atau, sxc, v1xc, v2xc, v3xc);
 					}
 					else
 					{
@@ -582,6 +657,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 	delete[] rhotmp1;
 	delete[] rhogsum1;
 	delete[] gdr1;
+	delete[] lapl1;
 	if(!is_stress) { delete[] h1;
 }
 
@@ -590,6 +666,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 		delete[] rhotmp2;
 		delete[] rhogsum2;
 		delete[] gdr2;
+		delete[] lapl2;
 		if(!is_stress) { delete[] h2;
 }
 	}
@@ -609,6 +686,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 		delete[] rhotmp2;
 		delete[] rhogsum2;
 		delete[] gdr2;
+		delete[] lapl2;
 	}
 
 	return;

source/module_hamilt_general/module_xc/xc_functional_libxc.h

@@ -82,6 +82,14 @@ namespace XC_Functional_Libxc
     extern std::vector<double> convert_sigma(
         const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr);
 
+    // calculating laplacian of density: ∇²ρ(r)
+    extern std::vector<double> cal_lapl(
+        const int nspin,
+        const std::size_t nrxx,
+        const std::vector<double> &rho,
+        const double tpiba2,
+        const Charge* const chr);
+
     // sgn for threshold mask
     extern std::vector<double> cal_sgn(
         const double rho_threshold,
@@ -166,7 +174,7 @@ namespace XC_Functional_Libxc
     // wrapper for the mGGA functionals
     extern void tau_xc(
         const std::vector<int> &func_id,
-        const double &rho, const double &grho, const double &atau, double &sxc,
+        const double &rho, const double &grho, const double &lapl_rho, const double &atau, double &sxc,
         double &v1xc, double &v2xc, double &v3xc);
 
     extern void tau_xc_spin(

source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp

@@ -82,6 +82,46 @@ XC_Functional_Libxc::cal_gdr(
 	return gdr;
 }
 
+// calculating laplacian of density: ∇²ρ(r)
+std::vector<double> XC_Functional_Libxc::cal_lapl(
+	const int nspin,
+	const std::size_t nrxx,
+	const std::vector<double> &rho,
+	const double tpiba2,
+	const Charge* const chr)
+{
+	std::vector<double> lapl(nspin * nrxx, 0.0);
+	for( int is=0; is!=nspin; ++is )
+	{
+		std::vector<double> rhor(nrxx);
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 1024)
+		#endif
+		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());
+
+		// ∇²ρ(r) = -∑_G |G|² ρ(G) e^{iGr}
+		std::vector<std::complex<double>> rhog_lapl(chr->rhopw->npw);
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 1024)
+		#endif
+		for(int ig=0; ig<chr->rhopw->npw; ++ig)
+			rhog_lapl[ig] = -rhog[ig] * chr->rhopw->gg[ig] * tpiba2;
+
+		std::vector<std::complex<double>> aux(chr->rhopw->nmaxgr);
+		chr->rhopw->recip2real(rhog_lapl.data(), aux.data());
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 1024)
+		#endif
+		for(std::size_t ir=0; ir<nrxx; ++ir)
+			lapl[is * nrxx + ir] = aux[ir].real();
+	} // end for(is)
+	return lapl;
+}
+
 // converting grho (abacus=>libxc)
 std::vector<double> XC_Functional_Libxc::convert_sigma(
 	const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr)