From 0f9438740f520f21e84705ef5b763eb8b3ea4c9d Mon Sep 17 00:00:00 2001 From: jiebin chen Date: Tue, 9 Jun 2026 10:30:47 +0000 Subject: [PATCH 1/2] 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 Co-authored-by: monkeycode-ai --- source/source_hamilt/module_xc/libxc_abacus.h | 3 + .../module_xc/libxc_mgga_wrap.cpp | 6 +- source/source_hamilt/module_xc/libxc_pot.cpp | 24 +++++- .../module_xc/test/CMakeLists.txt | 13 +++ .../source_hamilt/module_xc/test/test_xc4.cpp | 3 +- .../source_hamilt/module_xc/test/test_xc6.cpp | 84 +++++++++++++++++++ .../source_hamilt/module_xc/xc_functional.h | 6 ++ source/source_hamilt/module_xc/xc_grad.cpp | 78 ++++++++++++++++- 8 files changed, 209 insertions(+), 8 deletions(-) create mode 100644 source/source_hamilt/module_xc/test/test_xc6.cpp diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index 652293cd22..3ea823d797 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/source/source_hamilt/module_xc/libxc_abacus.h @@ -198,6 +198,7 @@ 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, @@ -211,6 +212,8 @@ namespace XC_Functional_Libxc double rhodw, ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, + double laplup, + double lapldw, double tauup, double taudw, double &sxc, diff --git a/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp index 038f965d15..c40f4327b0 100644 --- a/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp @@ -17,6 +17,7 @@ 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, @@ -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 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 gdr1, ModuleBase::Vector3 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 v1xc = {0.0, 0.0}; std::array v3xc = {0.0, 0.0}; - std::array lapl = {0.0, 0.0}; + std::array lapl = {laplup, lapldw}; std::array vlapl = {0.0, 0.0}; std::array v2xc = {0.0, 0.0, 0.0}; // call Libxc function: xc_mgga_exc_vxc diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index 9ff7190764..9ef5555595 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -12,6 +12,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, @@ -237,6 +238,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); @@ -315,7 +337,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, diff --git a/source/source_hamilt/module_xc/test/CMakeLists.txt b/source/source_hamilt/module_xc/test/CMakeLists.txt index 644f99dde9..ca205f45e4 100644 --- a/source/source_hamilt/module_xc/test/CMakeLists.txt +++ b/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 ) \ No newline at end of file diff --git a/source/source_hamilt/module_xc/test/test_xc4.cpp b/source/source_hamilt/module_xc/test/test_xc4.cpp index 467cc2667c..7d17256fab 100644 --- a/source/source_hamilt/module_xc/test/test_xc4.cpp +++ b/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); 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..76a054de8f --- /dev/null +++ b/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 +#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; + 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; +} \ No newline at end of file diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index 37695900a9..c4e63c8bce 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/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* rhog, + double* lapl, + const ModulePW::PW_Basis* rho_basis, + const double tpiba); + 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 0129662543..c2f98d68fc 100644 --- a/source/source_hamilt/module_xc/xc_grad.cpp +++ b/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* rhog, + double* lapl, + const ModulePW::PW_Basis* rho_basis, + const double tpiba) +{ + std::complex* lapl_tmp = new std::complex[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, lapl_tmp); + for(int ir=0; irnrxx; ir++) + { + lapl[ir] += lapl_tmp[ir].real() * tpiba * tpiba; + } + } + + delete[] lapl_tmp; +} + + void XC_Functional::noncolin_rho( double *rhoout1, double *rhoout2, From 4a546e64324dbbe2b4c4dad47b399b521de92294 Mon Sep 17 00:00:00 2001 From: jiebin chen Date: Tue, 9 Jun 2026 11:54:35 +0000 Subject: [PATCH 2/2] fix: use SCAN-L instead of SCAN for laplacian sensitivity test Co-authored-by: monkeycode-ai --- source/source_hamilt/module_xc/test/test_xc6.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/source_hamilt/module_xc/test/test_xc6.cpp b/source/source_hamilt/module_xc/test/test_xc6.cpp index 76a054de8f..16f1dfc1bc 100644 --- a/source/source_hamilt/module_xc/test/test_xc6.cpp +++ b/source/source_hamilt/module_xc/test/test_xc6.cpp @@ -36,7 +36,7 @@ class XCTest_SCANL_Laplacian : public XCTest void SetUp() { - XC_Functional::set_xc_type("SCAN"); + XC_Functional::set_xc_type("MGGA_X_SCANL+MGGA_C_SCANL"); const double rho = 0.17E+01; const double grho = 0.81E-11;