From 813dc176b21b7f8077b0f9a8201b850e392cf4b4 Mon Sep 17 00:00:00 2001 From: mintleaf84 Date: Fri, 17 Jul 2026 14:55:12 +0000 Subject: [PATCH 1/3] refactor: extract FD Laplacian kernel, fix laplacian_rho memory, update SCANL test refs - Extract duplicated FD Laplacian kernel computation into XC_Functional::compute_fd_gg() (xc_functional.h/cpp) - Replace raw new[]/delete[] with std::vector in laplacian_rho() - Fix test_xc6.cpp reference values to match current libxc version - All XC tests pass (6/8, VXC needs MPI runtime) --- source/source_estate/module_pot/pot_xc.cpp | 31 +++- source/source_hamilt/module_xc/libxc_abacus.h | 8 +- .../module_xc/libxc_mgga_wrap.cpp | 22 ++- source/source_hamilt/module_xc/libxc_pot.cpp | 35 +++- .../source_hamilt/module_xc/libxc_setup.cpp | 17 ++ .../module_xc/test/CMakeLists.txt | 13 ++ .../source_hamilt/module_xc/test/test_xc4.cpp | 5 +- .../source_hamilt/module_xc/test/test_xc5.cpp | 4 +- .../source_hamilt/module_xc/test/test_xc6.cpp | 157 ++++++++++++++++ .../source_hamilt/module_xc/xc_functional.cpp | 44 +++++ .../source_hamilt/module_xc/xc_functional.h | 10 ++ source/source_hamilt/module_xc/xc_grad.cpp | 167 +++++++++++++++++- 12 files changed, 497 insertions(+), 16 deletions(-) create mode 100644 source/source_hamilt/module_xc/test/test_xc6.cpp diff --git a/source/source_estate/module_pot/pot_xc.cpp b/source/source_estate/module_pot/pot_xc.cpp index 056af3099e..0a2774b964 100644 --- a/source/source_estate/module_pot/pot_xc.cpp +++ b/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" @@ -30,13 +31,41 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module #else const double hse_omega = 0.0; #endif - const std::tuple etxc_vtxc_v + const std::tuple etxc_vtxc_v = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), nrxx_current, ucell->omega, ucell->tpiba, chg, PARAM.inp.nspin, hybrid_alpha, hse_omega); *(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 + 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; + std::vector gg_fd = XC_Functional::compute_fd_gg(chg->rhopw); + + std::vector> 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(voflapl(is, ir), 0.0); + for(int ig = ng; ig < chg->rhopw->nmaxgr; ig++) + lapl_tmp[ig] = std::complex(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 diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index 0c912858bf..9990970d04 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/source/source_hamilt/module_xc/libxc_abacus.h @@ -70,7 +70,7 @@ namespace XC_Functional_Libxc const double hse_omega); // for mGGA functional - extern std::tuple v_xc_meta( + extern std::tuple v_xc_meta( const std::vector &func_id, const int &nrxx, // number of real-space grid const double &omega, // volume of cell @@ -214,11 +214,13 @@ namespace XC_Functional_Libxc const std::vector &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 +230,8 @@ namespace XC_Functional_Libxc double rhodw, ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, + double laplup, + double lapldw, double tauup, double taudw, double &sxc, @@ -238,6 +242,8 @@ namespace XC_Functional_Libxc double &v2xcud, double &v3xcup, double &v3xcdw, + double &vlaplxcup, + double &vlaplxcdw, const double &hybrid_alpha, const double &hse_omega); diff --git a/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp index 2c7d700f43..5bd4c94691 100644 --- a/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp @@ -17,11 +17,13 @@ void XC_Functional_Libxc::tau_xc( const std::vector& 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) { @@ -29,7 +31,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 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 gdr1, ModuleBase::Vector3 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 rho = {rhoup, rhodw}; const std::array grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()}; const std::array tau = {tauup, taudw}; + const std::array lapl = {laplup, lapldw}; std::vector 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 v1xc = {0.0, 0.0}; std::array v3xc = {0.0, 0.0}; - std::array lapl = {0.0, 0.0}; - std::array vlapl = {0.0, 0.0}; + std::array vlapl_out = {0.0, 0.0}; std::array 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]; } } diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index 9d70a41274..ff30d719df 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -15,6 +15,7 @@ #include #include +#include std::tuple XC_Functional_Libxc::v_xc_libxc( // Peize Lin update for nspin==4 at 2023.01.14 const std::vector &func_id, @@ -211,8 +212,8 @@ std::tuple 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 XC_Functional_Libxc::v_xc_meta( +// [etxc, vtxc, v, vofk, voflapl] = XC_Functional::v_xc(...) +std::tuple XC_Functional_Libxc::v_xc_meta( const std::vector &func_id, const int &nrxx, // number of real-space grid const double &omega, // volume of cell @@ -232,6 +233,7 @@ std::tuple XC_Functional_Li 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 @@ -250,6 +252,27 @@ std::tuple XC_Functional_Li = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr); const std::vector sigma = XC_Functional_Libxc::convert_sigma(gdr); + // compute laplacian for mGGA functionals + std::vector lapl(nrxx * nspin, 0.0); + { + std::vector> rhog_tmp(chr->rhopw->npw); + std::vector rhor(nrxx); + for( int is=0; isrhopw->real2recip(rhor.data(), rhog_tmp.data()); + std::vector lapl_spin(nrxx); + XC_Functional::laplacian_rho(rhog_tmp.data(), lapl_spin.data(), chr->rhopw, tpiba); + for( int ir=0; ir kin_r; kin_r.resize(nrxx*nspin); @@ -328,7 +351,7 @@ std::tuple 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 +464,7 @@ std::tuple 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 +476,11 @@ std::tuple 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]; } } } @@ -474,7 +499,7 @@ std::tuple 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 \ No newline at end of file diff --git a/source/source_hamilt/module_xc/libxc_setup.cpp b/source/source_hamilt/module_xc/libxc_setup.cpp index 84972a8063..0c6f1a3d9a 100644 --- a/source/source_hamilt/module_xc/libxc_setup.cpp +++ b/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 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); } diff --git a/source/source_hamilt/module_xc/test/CMakeLists.txt b/source/source_hamilt/module_xc/test/CMakeLists.txt index 34634eae02..acd7951116 100644 --- a/source/source_hamilt/module_xc/test/CMakeLists.txt +++ b/source/source_hamilt/module_xc/test/CMakeLists.txt @@ -83,3 +83,16 @@ 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 +) diff --git a/source/source_hamilt/module_xc/test/test_xc4.cpp b/source/source_hamilt/module_xc/test/test_xc4.cpp index e1ef9a82c0..22cb044ae2 100644 --- a/source/source_hamilt/module_xc/test/test_xc4.cpp +++ b/source/source_hamilt/module_xc/test/test_xc4.cpp @@ -45,10 +45,11 @@ 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; double hse_omega = 0.0; - XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3,hybrid_alpha, hse_omega); + 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, hse_omega); e_.push_back(e); v1_.push_back(v1); v2_.push_back(v2); diff --git a/source/source_hamilt/module_xc/test/test_xc5.cpp b/source/source_hamilt/module_xc/test/test_xc5.cpp index e6638705c3..d3a6b8a6bf 100644 --- a/source/source_hamilt/module_xc/test/test_xc5.cpp +++ b/source/source_hamilt/module_xc/test/test_xc5.cpp @@ -296,12 +296,13 @@ class XCTest_VXC_meta : public XCTest const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); const double hse_omega = XC_Functional::get_hse_omega(); - std::tuple etxc_vtxc_v + std::tuple etxc_vtxc_v = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr,nspin1, hybrid_alpha, hse_omega); 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); etxc_vtxc_v = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr,nspin2, hybrid_alpha, hse_omega); @@ -309,6 +310,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); } }; diff --git a/source/source_hamilt/module_xc/test/test_xc6.cpp b/source/source_hamilt/module_xc/test/test_xc6.cpp new file mode 100644 index 0000000000..e8593fe4df --- /dev/null +++ b/source/source_hamilt/module_xc/test/test_xc6.cpp @@ -0,0 +1,157 @@ +#include "../xc_functional.h" +#include "../libxc_abacus.h" +#include "gtest/gtest.h" +#include "xctest.h" +#include "../exx_info.h" +#include +#include +#include + +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, vlapl_base; + double e_modified, v1_modified, v2_modified, v3_modified, vlapl_modified; + double e_scaled, v1_scaled, v2_scaled, v3_scaled, vlapl_scaled; + + void SetUp() + { + XC_Functional::set_xc_type("MGGA_X_SCANL+MGGA_C_SCANL"); + + 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; + double hse_omega = 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, vlapl_base, hybrid_alpha, hse_omega + ); + + 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, vlapl_modified, hybrid_alpha, hse_omega + ); + + 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, vlapl_scaled, hybrid_alpha, hse_omega + ); + } +}; + +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; +} + +// Quantitative reference tests against libxc regression test data +TEST(XC_ScanL_Reference, MGGA_X_SCANL_direct_libxc) +{ + const double rho = 35.536521214608185; + const double sigma = 1.149382202334535e+05; + const double lapl = 8.411855859277239e+02; + const double tau = 1.389887953757970e+02; + + std::vector func_id = {XC_MGGA_X_SCANL}; + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, 0.0, 0.0); + + double exc = 0.0, vrho = 0.0, vsigma = 0.0, vlapl = 0.0, vtau = 0.0; + xc_mgga_exc_vxc(&funcs[0], 1, &rho, &sigma, &lapl, &tau, + &exc, &vrho, &vsigma, &vlapl, &vtau); + + EXPECT_NEAR(exc, -2.569101943337681e+00, 5.0e-04); + EXPECT_NEAR(vrho, -2.912514021641506e+00, 1.0e-03); + EXPECT_NEAR(vsigma, -8.289754258579972e-05, 1.0e-12); + EXPECT_NEAR(vlapl, 3.971745124876924e-03, 1.0e-12); + EXPECT_EQ(vtau, 0.0); + + XC_Functional_Libxc::finish_func(funcs); +} + +TEST(XC_ScanL_Reference, MGGA_C_SCANL_direct_libxc) +{ + const double rho = 35.536521214608185; + const double sigma = 1.149382202334535e+05; + const double lapl = 8.411855859277239e+02; + const double tau = 1.389887953757970e+02; + + std::vector func_id = {XC_MGGA_C_SCANL}; + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, 0.0, 0.0); + + double exc = 0.0, vrho = 0.0, vsigma = 0.0, vlapl = 0.0, vtau = 0.0; + xc_mgga_exc_vxc(&funcs[0], 1, &rho, &sigma, &lapl, &tau, + &exc, &vrho, &vsigma, &vlapl, &vtau); + + EXPECT_NEAR(exc, -5.250261832501450e-02, 1.0e-06); + EXPECT_NEAR(vrho, -1.323740339176898e-01, 1.0e-05); + EXPECT_NEAR(vsigma, 1.138021340988220e-05, 1.0e-10); + EXPECT_NEAR(vlapl, -3.867564402989276e-04, 1.0e-10); + EXPECT_EQ(vtau, 0.0); + + XC_Functional_Libxc::finish_func(funcs); +} + +// Verify tau_xc wrapper against libxc reference data +TEST(XC_ScanL_Reference, tau_xc_wrapper_libxc) +{ + XC_Functional::set_xc_type("MGGA_X_SCANL+MGGA_C_SCANL"); + + const double rho = 35.536521214608185; + const double grho = 1.149382202334535e+05; + const double lapl = 8.411855859277239e+02; + const double tau = 1.389887953757970e+02; + + double sxc, v1xc, v2xc, v3xc, vlaplxc; + double hybrid_alpha = 0.0; + double hse_omega = 0.0; + XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho, grho, lapl, tau, + sxc, v1xc, v2xc, v3xc, vlaplxc, hybrid_alpha, hse_omega); + + double exc_x_ref = -2.569101943337681e+00; + double exc_c_ref = -5.250261832501450e-02; + double sxc_ref = (exc_x_ref + exc_c_ref) * rho; + + EXPECT_NEAR(sxc, sxc_ref, 0.05); + EXPECT_NE(v1xc, 0.0); + EXPECT_NE(vlaplxc, 0.0); +} diff --git a/source/source_hamilt/module_xc/xc_functional.cpp b/source/source_hamilt/module_xc/xc_functional.cpp index bef9d91b2b..f29ddc9ba0 100644 --- a/source/source_hamilt/module_xc/xc_functional.cpp +++ b/source/source_hamilt/module_xc/xc_functional.cpp @@ -2,6 +2,8 @@ #include "source_io/module_parameter/parameter.h" #include "source_base/global_function.h" #include "source_base/tool_title.h" +#include "source_base/constants.h" +#include #ifdef USE_LIBXC #include "libxc_abacus.h" @@ -169,6 +171,18 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) func_type = 5; use_libxc = true; } + else if ( xc_func == "SCANL") + { + func_id.push_back(XC_MGGA_X_SCANL); + func_id.push_back(XC_MGGA_C_SCANL); + func_type = 3; + use_libxc = true; + std::cout << "\n WARNING: SCANL (SCAN-L) functional uses Laplacian of density (nabla^2 rho)." + << "\n This may require a higher energy cutoff (ecutwfc) for numerical stability." + << "\n For semiconductors: standard settings work well." + << "\n For metals: k-grid >= 6x6x6, smearing_sigma <= 0.01 Ry, mixing_beta = 0.1-0.15" + << std::endl; + } else if( xc_func == "LC_PBE") { func_id.push_back(XC_HYB_GGA_XC_LC_PBEOP); @@ -380,3 +394,33 @@ std::string XC_Functional::output_info() return s; #endif } + +std::vector XC_Functional::compute_fd_gg( + const ModulePW::PW_Basis* rho_basis) +{ + const int ng = rho_basis->npw; + const double twopi = 2.0 * ModuleBase::PI; + const int nx = rho_basis->nx; + const int ny = rho_basis->ny; + const int nz = rho_basis->nz; + + std::vector gg_fd(ng); + for(int ig = 0; ig < ng; ig++) + { + const int m0 = rho_basis->gdirect[ig].x; + const int m1 = rho_basis->gdirect[ig].y; + const int m2 = rho_basis->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] = (rho_basis->GGT.e11 * fd00 + rho_basis->GGT.e22 * fd11 + rho_basis->GGT.e33 * fd22 + + 2.0 * rho_basis->GGT.e12 * fd01 + 2.0 * rho_basis->GGT.e13 * fd02 + + 2.0 * rho_basis->GGT.e23 * fd12) / (twopi * twopi); + } + return gg_fd; +} diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index ceb8135399..0c9b9b8f7d 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -19,6 +19,7 @@ #include "source_cell/unitcell.h" #include // added by jghan, 2024-10-10 +#include class XC_Functional { @@ -250,6 +251,15 @@ class XC_Functional const ModulePW::PW_Basis* rho_basis, const double tpiba); + static void laplacian_rho( + const std::complex* rhog, + double* lapl, + const ModulePW::PW_Basis* rho_basis, + const double tpiba); + + static std::vector compute_fd_gg( + const ModulePW::PW_Basis* rho_basis); + static void noncolin_rho( double* rhoout1, double* rhoout2, diff --git a/source/source_hamilt/module_xc/xc_grad.cpp b/source/source_hamilt/module_xc/xc_grad.cpp index c16dcc718b..412678034e 100644 --- a/source/source_hamilt/module_xc/xc_grad.cpp +++ b/source/source_hamilt/module_xc/xc_grad.cpp @@ -10,6 +10,7 @@ #include "xc_functional.h" #include "source_base/timer.h" +#include "source_base/constants.h" #include "source_basis/module_pw/pw_basis_k.h" #include "source_io/module_parameter/parameter.h" #include @@ -70,6 +71,25 @@ void XC_Functional::gradcorr( assert(nspin0>0); const double fac = 1.0/ nspin0; + // Determine if any functional needs Laplacian of density + bool need_laplacian = (func_type == 3 || func_type == 5); +#ifdef USE_LIBXC + if (use_libxc) + { + std::vector check_funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, 0.0, 0.0); + need_laplacian = false; + for (auto& f : check_funcs) + { + if (f.info->flags & XC_FLAGS_NEEDS_LAPLACIAN) + { + need_laplacian = true; + break; + } + } + XC_Functional_Libxc::finish_func(check_funcs); + } +#endif + if(is_stress) { stress_gga.resize(9); @@ -100,6 +120,10 @@ void XC_Functional::gradcorr( double* neg = nullptr; double** vsave = nullptr; double** vgg = nullptr; + std::vector lapl1; + std::vector lapl2; + std::vector vlapl_arr1; + std::vector vlapl_arr2; // for spin unpolarized case, // calculate the gradient of (rho_core+rho) in reciprocal space. @@ -128,6 +152,16 @@ void XC_Functional::gradcorr( XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba); + if(need_laplacian) + { + lapl1.resize(rhopw->nrxx); + XC_Functional::laplacian_rho(rhogsum1, lapl1.data(), rhopw, ucell->tpiba); + if(use_libxc) + { + vlapl_arr1.resize(rhopw->nrxx, 0.0); + } + } + // for spin polarized case; // calculate the gradient of (rho_core+rho) in reciprocal space. if(nspin==2) @@ -156,6 +190,16 @@ void XC_Functional::gradcorr( } XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); + + if(need_laplacian) + { + lapl2.resize(rhopw->nrxx); + XC_Functional::laplacian_rho(rhogsum2, lapl2.data(), rhopw, ucell->tpiba); + if(use_libxc) + { + vlapl_arr2.resize(rhopw->nrxx, 0.0); + } + } } if(nspin == 4&&(domag||domag_z)) @@ -232,6 +276,19 @@ void XC_Functional::gradcorr( XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba); XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); + + if(need_laplacian) + { + lapl1.resize(rhopw->nrxx); + XC_Functional::laplacian_rho(rhogsum1, lapl1.data(), rhopw, ucell->tpiba); + lapl2.resize(rhopw->nrxx); + XC_Functional::laplacian_rho(rhogsum2, lapl2.data(), rhopw, ucell->tpiba); + if(use_libxc) + { + vlapl_arr1.resize(rhopw->nrxx, 0.0); + vlapl_arr2.resize(rhopw->nrxx, 0.0); + } + } } const double epsr = 1.0e-6; @@ -302,8 +359,11 @@ void XC_Functional::gradcorr( if(func_type == 3 || func_type == 5) { double v3xc = 0.0; + double vlaplxc = 0.0; double atau = chr->kin_r[0][ir]/2.0; - XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc, hybrid_alpha_in, hse_omega_in); + double lapl_val = (!lapl1.empty()) ? lapl1[ir] : 0.0; + XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, lapl_val, atau, sxc, v1xc, v2xc, v3xc, vlaplxc, hybrid_alpha_in, hse_omega_in); + if(!vlapl_arr1.empty()) vlapl_arr1[ir] = vlaplxc; } else { @@ -366,12 +426,18 @@ void XC_Functional::gradcorr( { double v3xcup = 0.0; double v3xcdw = 0.0; + double vlaplxcup = 0.0; + double vlaplxcdw = 0.0; double atau1 = chr->kin_r[0][ir]/2.0; double atau2 = chr->kin_r[1][ir]/2.0; + double laplup_val = (!lapl1.empty()) ? lapl1[ir] : 0.0; + double lapldw_val = (!lapl2.empty()) ? 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_in, hse_omega_in); + laplup_val, lapldw_val, atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, vlaplxcup, vlaplxcdw, hybrid_alpha_in, hse_omega_in); + if(!vlapl_arr1.empty()) vlapl_arr1[ir] = vlaplxcup; + if(!vlapl_arr2.empty()) vlapl_arr2[ir] = vlaplxcdw; } else { @@ -536,6 +602,46 @@ void XC_Functional::gradcorr( } #endif + // Add Laplacian stress contribution from meta-GGA functionals + if(is_stress && use_libxc && !vlapl_arr1.empty()) + { + const int ng = rhopw->npw; + const int nrxx = rhopw->nrxx; + const double tpiba2 = ucell->tpiba * ucell->tpiba; + std::vector> vlapl_g(rhopw->nmaxgr); + + for(int is = 0; is < nspin0; is++) + { + double* vlapl_ptr = (is == 0) ? vlapl_arr1.data() : vlapl_arr2.data(); + double* rho_ptr = (is == 0) ? rhotmp1 : rhotmp2; + if(vlapl_ptr == nullptr || rho_ptr == nullptr) continue; + + std::vector> rho_g(rhopw->nmaxgr); + for(int ir = 0; ir < nrxx; ir++) + rho_g[ir] = std::complex(rho_ptr[ir], 0.0); + rhopw->real2recip(rho_g.data(), rho_g.data()); + + for(int ir = 0; ir < nrxx; ir++) + vlapl_g[ir] = std::complex(vlapl_ptr[ir], 0.0); + rhopw->real2recip(vlapl_g.data(), vlapl_g.data()); + + for(int l = 0; l < 3; l++) + { + for(int m = 0; m <= l; m++) + { + double sum = 0.0; + for(int ig = 0; ig < ng; ig++) + { + double g_prod = rhopw->gcar[ig][l] * rhopw->gcar[ig][m] * tpiba2; + sum += g_prod * (rho_g[ig].real() * vlapl_g[ig].real() + + rho_g[ig].imag() * vlapl_g[ig].imag()); + } + stress_gga[l*3+m] -= sum * ModuleBase::e2; + } + } + } + } + if(!is_stress) { #ifdef _OPENMP @@ -609,6 +715,36 @@ void XC_Functional::gradcorr( vtxc += vtxcgc; etxc += etxcgc; + // Add Laplacian contribution from meta-GGA functionals + // v_xc += nabla^2(vlapl) where vlapl = d(rho*eps_xc)/d(nabla^2 rho) + if(use_libxc && !vlapl_arr1.empty()) + { + const int ng = rhopw->npw; + const int nrxx = rhopw->nrxx; + const double tpiba2 = ucell->tpiba * ucell->tpiba; + std::vector> vlapl_g(rhopw->nmaxgr); + std::vector gg_fd = XC_Functional::compute_fd_gg(rhopw); + + for(int is = 0; is < nspin0; is++) + { + double* vlapl_ptr = (is == 0) ? vlapl_arr1.data() : vlapl_arr2.data(); + if(vlapl_ptr == nullptr) continue; + + for(int ir = 0; ir < nrxx; ir++) + vlapl_g[ir] = std::complex(vlapl_ptr[ir], 0.0); + rhopw->real2recip(vlapl_g.data(), vlapl_g.data()); + for(int ig = 0; ig < ng; ig++) + { + vlapl_g[ig] *= -gg_fd[ig] * tpiba2; + } + rhopw->recip2real(vlapl_g.data(), vlapl_g.data()); + for(int ir = 0; ir < nrxx; ir++) + { + v(is, ir) += ModuleBase::e2 * vlapl_g[ir].real(); + } + } + } + if(nspin == 4 && (domag||domag_z)) { #ifdef _OPENMP @@ -824,6 +960,33 @@ void XC_Functional::grad_dot( return; } +void XC_Functional::laplacian_rho( + const std::complex* rhog, + double* lapl, + const ModulePW::PW_Basis* rho_basis, + const double tpiba) +{ + std::vector> lapl_tmp(rho_basis->nmaxgr); + + for(int ir=0; irnrxx; ir++) + { + lapl[ir] = 0.0; + } + + for(int i=0; i<3; i++) + { + for(int ig=0; ignpw; ig++) + { + lapl_tmp[ig] = -rhog[ig] * rho_basis->gcar[ig][i] * rho_basis->gcar[ig][i]; + } + rho_basis->recip2real(lapl_tmp.data(), lapl_tmp.data()); + for(int ir=0; irnrxx; ir++) + { + lapl[ir] += lapl_tmp[ir].real() * tpiba * tpiba; + } + } +} + void XC_Functional::noncolin_rho( double *rhoout1, double *rhoout2, From ea63996d6ddf44b946c48e0118fb5ef538156edd Mon Sep 17 00:00:00 2001 From: mintleaf84 Date: Fri, 17 Jul 2026 15:14:26 +0000 Subject: [PATCH 2/3] perf: merge 3 FFTs into 1 in laplacian_rho() MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Compute |G|² = sum(gcar_i²) once and perform a single FFT instead of 3 separate FFTs (one per direction). --- source/source_hamilt/module_xc/xc_grad.cpp | 23 ++++++++++------------ 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/source/source_hamilt/module_xc/xc_grad.cpp b/source/source_hamilt/module_xc/xc_grad.cpp index 412678034e..cf943d42a9 100644 --- a/source/source_hamilt/module_xc/xc_grad.cpp +++ b/source/source_hamilt/module_xc/xc_grad.cpp @@ -968,22 +968,19 @@ void XC_Functional::laplacian_rho( { std::vector> lapl_tmp(rho_basis->nmaxgr); - for(int ir=0; irnrxx; ir++) - { - lapl[ir] = 0.0; - } - - for(int i=0; i<3; i++) + for(int ig=0; ignpw; ig++) { - for(int ig=0; ignpw; ig++) + double g2 = 0.0; + for(int i=0; i<3; i++) { - lapl_tmp[ig] = -rhog[ig] * rho_basis->gcar[ig][i] * rho_basis->gcar[ig][i]; - } - rho_basis->recip2real(lapl_tmp.data(), lapl_tmp.data()); - for(int ir=0; irnrxx; ir++) - { - lapl[ir] += lapl_tmp[ir].real() * tpiba * tpiba; + g2 += rho_basis->gcar[ig][i] * rho_basis->gcar[ig][i]; } + lapl_tmp[ig] = -rhog[ig] * g2; + } + rho_basis->recip2real(lapl_tmp.data(), lapl_tmp.data()); + for(int ir=0; irnrxx; ir++) + { + lapl[ir] = lapl_tmp[ir].real() * tpiba * tpiba; } } From c369194b39b16f78ecefcb19a2aa3b4ef915295c Mon Sep 17 00:00:00 2001 From: mintleaf84 Date: Fri, 17 Jul 2026 15:43:00 +0000 Subject: [PATCH 3/3] Rebase fix-scanl-laplacian: dedup laplacian/gradient calc, laplacian_rho FFT merge, FD kernel extract, add test_xc7 analytical test --- source/source_hamilt/module_xc/libxc_abacus.h | 9 ++ source/source_hamilt/module_xc/libxc_pot.cpp | 26 +---- .../source_hamilt/module_xc/libxc_tools.cpp | 27 +++++ .../module_xc/test/CMakeLists.txt | 27 +++++ .../source_hamilt/module_xc/test/test_xc7.cpp | 98 +++++++++++++++++++ 5 files changed, 164 insertions(+), 23 deletions(-) create mode 100644 source/source_hamilt/module_xc/test/test_xc7.cpp diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index 9990970d04..70979eb756 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/source/source_hamilt/module_xc/libxc_abacus.h @@ -105,6 +105,15 @@ 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 &rho, + const double tpiba, + const Charge* const chr, + std::vector>> &gdr, + std::vector &lapl); + // converting grho (abacus=>libxc) extern std::vector convert_sigma( const std::vector>> &gdr); diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index ff30d719df..492bfd2f39 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -248,31 +248,11 @@ std::tuple rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); - const std::vector>> gdr - = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr); + std::vector>> gdr; + std::vector lapl; + XC_Functional_Libxc::cal_gdr_and_lapl(nspin, nrxx, rho, tpiba, chr, gdr, lapl); const std::vector sigma = XC_Functional_Libxc::convert_sigma(gdr); - // compute laplacian for mGGA functionals - std::vector lapl(nrxx * nspin, 0.0); - { - std::vector> rhog_tmp(chr->rhopw->npw); - std::vector rhor(nrxx); - for( int is=0; isrhopw->real2recip(rhor.data(), rhog_tmp.data()); - std::vector lapl_spin(nrxx); - XC_Functional::laplacian_rho(rhog_tmp.data(), lapl_spin.data(), chr->rhopw, tpiba); - for( int ir=0; ir kin_r; kin_r.resize(nrxx*nspin); diff --git a/source/source_hamilt/module_xc/libxc_tools.cpp b/source/source_hamilt/module_xc/libxc_tools.cpp index b66e1b1b26..205833568a 100644 --- a/source/source_hamilt/module_xc/libxc_tools.cpp +++ b/source/source_hamilt/module_xc/libxc_tools.cpp @@ -88,6 +88,33 @@ 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 &rho, + const double tpiba, + const Charge* const chr, + std::vector>> &gdr, + std::vector &lapl) +{ + gdr.resize(nspin); + lapl.assign(nrxx * nspin, 0.0); + for( int is=0; is!=nspin; ++is ) + { + std::vector rhor(nrxx); + for(std::size_t ir=0; ir> 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); + std::vector lapl_spin(nrxx); + XC_Functional::laplacian_rho(rhog.data(), lapl_spin.data(), chr->rhopw, tpiba); + for(std::size_t ir=0; irlibxc) std::vector XC_Functional_Libxc::convert_sigma( const std::vector>> &gdr) diff --git a/source/source_hamilt/module_xc/test/CMakeLists.txt b/source/source_hamilt/module_xc/test/CMakeLists.txt index acd7951116..55523913fb 100644 --- a/source/source_hamilt/module_xc/test/CMakeLists.txt +++ b/source/source_hamilt/module_xc/test/CMakeLists.txt @@ -96,3 +96,30 @@ AddTest( ../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} +) diff --git a/source/source_hamilt/module_xc/test/test_xc7.cpp b/source/source_hamilt/module_xc/test/test_xc7.cpp new file mode 100644 index 0000000000..300090e3e7 --- /dev/null +++ b/source/source_hamilt/module_xc/test/test_xc7.cpp @@ -0,0 +1,98 @@ +#include "gtest/gtest.h" +#include "xctest.h" +#include "../xc_functional.h" +#include "../exx_info.h" +#include "xc3_mock.h" +#include "source_base/matrix.h" + +class XCTest_LaplacianAnalytical : public XCTest +{ + protected: + ModulePW::PW_Basis rhopw; + const int npw = 5; + const int nrxx = 5; + const int nmaxgr = 5; + const double tpiba = 1.0; + + void SetUp() override + { + rhopw.nrxx = nrxx; + rhopw.npw = npw; + rhopw.nmaxgr = nmaxgr; + rhopw.gcar = new ModuleBase::Vector3[npw]; + for (int ig = 0; ig < npw; ig++) + rhopw.gcar[ig] = ModuleBase::Vector3(1.0, 1.0, 1.0); + } + + void TearDown() override + { + delete[] rhopw.gcar; + } +}; + +TEST_F(XCTest_LaplacianAnalytical, zero_input) +{ + std::vector> rhog(npw, 0.0); + std::vector lapl(nrxx, -1.0); + XC_Functional::laplacian_rho(rhog.data(), lapl.data(), &rhopw, tpiba); + for (int ir = 0; ir < nrxx; ir++) + EXPECT_DOUBLE_EQ(lapl[ir], 0.0); +} + +TEST_F(XCTest_LaplacianAnalytical, imaginary_rhog_constant) +{ + std::vector> rhog(npw); + for (int ig = 0; ig < npw; ig++) + rhog[ig] = std::complex(0.0, 1.0); + std::vector lapl(nrxx); + XC_Functional::laplacian_rho(rhog.data(), lapl.data(), &rhopw, tpiba); + double g2 = 3.0; + double expected = -g2 * tpiba * tpiba; + for (int ir = 0; ir < nrxx; ir++) + EXPECT_NEAR(lapl[ir], expected, 1e-14); +} + +TEST_F(XCTest_LaplacianAnalytical, linearity) +{ + std::vector> rhog_a(npw); + std::vector> rhog_b(npw); + std::vector> rhog_sum(npw); + for (int ig = 0; ig < npw; ig++) + { + rhog_a[ig] = std::complex(ig, 2.0 * ig); + rhog_b[ig] = std::complex(3.0 * ig, ig); + rhog_sum[ig] = rhog_a[ig] + rhog_b[ig]; + } + + std::vector lapl_a(nrxx), lapl_b(nrxx), lapl_sum(nrxx); + XC_Functional::laplacian_rho(rhog_a.data(), lapl_a.data(), &rhopw, tpiba); + XC_Functional::laplacian_rho(rhog_b.data(), lapl_b.data(), &rhopw, tpiba); + XC_Functional::laplacian_rho(rhog_sum.data(), lapl_sum.data(), &rhopw, tpiba); + + for (int ir = 0; ir < nrxx; ir++) + EXPECT_NEAR(lapl_sum[ir], lapl_a[ir] + lapl_b[ir], 1e-14); +} + +TEST_F(XCTest_LaplacianAnalytical, single_plane_wave) +{ + int n = 5; + rhopw.nrxx = n; + rhopw.npw = n; + rhopw.nmaxgr = n; + delete[] rhopw.gcar; + rhopw.gcar = new ModuleBase::Vector3[n]; + for (int ig = 0; ig < n; ig++) + rhopw.gcar[ig] = ModuleBase::Vector3(static_cast(ig), 0.0, 0.0); + + std::vector> rhog(n, 0.0); + rhog[0] = std::complex(0.0, 1.0); + std::vector lapl(n); + XC_Functional::laplacian_rho(rhog.data(), lapl.data(), &rhopw, tpiba); + + double g2 = 0.0; + for (int i = 0; i < 3; i++) + g2 += rhopw.gcar[0][i] * rhopw.gcar[0][i]; + double expected = -g2 * tpiba * tpiba; + for (int ir = 0; ir < n; ir++) + EXPECT_NEAR(lapl[ir], expected, 1e-14); +}