From d7a2179e76d4dda4f993963106e42322cec74933 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 6 Oct 2026 11:13:54 +0800 Subject: [PATCH 1/2] fix(module_xc): avoid per-grid-point libxc init_func in per-point wrappers issue #6286 reported that init_func was called thousands of times after SCF convergence. The root cause is that the per-point libxc wrappers (xc_spin_libxc, gcxc_libxc, gcxc_spin_libxc, tau_xc, tau_xc_spin) called init_func -> evaluate(1 point) -> finish_func inside every grid point of the OpenMP loops. Since init_func only depends on the functional id, polarization, and runtime parameters (alpha, omega) -- none of which vary across grid points -- it only needs to run once per thread. Changes: - Add funcs-taking overloads to all five per-point wrappers; the old func_id overloads become thin wrappers that init/evaluate/finish. - In gradcorr_xc_kernel, init the functional once per OpenMP thread (inside #pragma omp parallel, before the for) and reuse it for all points that thread handles. - In v_xc's nspin==4 branch, split #pragma omp parallel for into parallel { init; for } so init runs once per thread. - Add a std::atomic init counter with get_and_reset_init_count() for tests to assert per-point paths do not re-initialize. - Add unit tests: funcs-overload vs func_id-overload bit-equality for LDA/GGA/mGGA, and an init-count assertion (old=1, new=0). init_func call count on the per-point paths drops from O(nrxx * SCF) to O(nthreads * SCF). No change to physical results. Verification: make -j 30 (ENABLE_LIBXC=ON, libxc 7.0.0) exit 0; all 10 XC test executables pass (VXC 18, GRADCORR 3, PBE 23, PZ_SPN 14, SCANL_LAPL 5, HSE 9, SCAN 1, LAPL 4, NCGGA_RADIAL 9, NCGGA_DISCRETE_FD 40); code_quality_score 8/8 files pass (avg 80.6, min 61). --- source/source_hamilt/module_xc/libxc_abacus.h | 75 ++++++++++++++ .../module_xc/libxc_gga_wrap.cpp | 77 ++++++++++----- .../module_xc/libxc_lda_wrap.cpp | 32 +++--- .../module_xc/libxc_mgga_wrap.cpp | 97 ++++++++++++++----- .../source_hamilt/module_xc/libxc_setup.cpp | 13 +++ .../source_hamilt/module_xc/test/test_xc2.cpp | 84 ++++++++++++++++ .../source_hamilt/module_xc/test/test_xc6.cpp | 31 ++++++ .../module_xc/xc_grad_kernel.cpp | 32 ++++-- source/source_hamilt/module_xc/xc_pot.cpp | 93 +++++++++++------- 9 files changed, 431 insertions(+), 103 deletions(-) diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index b618bbbe908..4110c414185 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/source/source_hamilt/module_xc/libxc_abacus.h @@ -73,6 +73,12 @@ namespace XC_Functional_Libxc extern void finish_func(std::vector &funcs); + // Returns the number of times init_func has been called since the last + // call to this function, and resets the counter to zero. Intended for + // unit tests that assert per-point evaluation paths do not re-initialize + // the libxc functional on every grid point. + extern int get_and_reset_init_count(); + //------------------- // libxc_pot.cpp @@ -266,6 +272,17 @@ namespace XC_Functional_Libxc const double hybrid_alpha, const double hse_omega); + // Overload accepting an already-initialized functional vector. The caller + // is responsible for init_func/finish_func; this overload does not touch + // the lifetime of funcs. Useful for per-thread reuse inside OpenMP loops. + extern void xc_spin_libxc( + const std::vector& funcs, + const double &rhoup, + const double &rhodw, + double &exc, + double &vxcup, + double &vxcdw); + //------------------- // libxc_gga_wrap.cpp @@ -282,6 +299,15 @@ namespace XC_Functional_Libxc const double hybrid_alpha, const double hse_omega); + // Overload accepting an already-initialized functional vector. + extern void gcxc_libxc( + const std::vector& funcs, + const double &rho, + const double &grho, + double &sxc, + double &v1xc, + double &v2xc); + // the entire GGA functional, for nspin=2 case extern void gcxc_spin_libxc( const std::vector &func_id, @@ -298,6 +324,20 @@ namespace XC_Functional_Libxc const double hybrid_alpha, const double hse_omega); + // Overload accepting an already-initialized functional vector. + extern void gcxc_spin_libxc( + const std::vector& funcs, + const double rhoup, + const double rhodw, + const ModuleBase::Vector3 gdr1, + const ModuleBase::Vector3 gdr2, + double &sxc, + double &v1xcup, + double &v1xcdw, + double &v2xcup, + double &v2xcdw, + double &v2xcud); + //------------------- // libxc_mgga_wrap.cpp @@ -318,6 +358,19 @@ namespace XC_Functional_Libxc const double &hybrid_alpha, const double &hse_omega); + // Overload accepting an already-initialized functional vector. + extern void tau_xc( + const std::vector& funcs, + const double &rho, + const double &grho, + const double &lapl_rho, + const double &atau, + double &sxc, + double &v1xc, + double &v2xc, + double &v3xc, + double &vlaplxc); + extern void tau_xc_spin( const std::vector &func_id, double rhoup, @@ -341,6 +394,28 @@ namespace XC_Functional_Libxc const double &hybrid_alpha, const double &hse_omega); + // Overload accepting an already-initialized functional vector. + extern void tau_xc_spin( + const std::vector& funcs, + double rhoup, + double rhodw, + ModuleBase::Vector3 gdr1, + ModuleBase::Vector3 gdr2, + double laplup, + double lapldw, + double tauup, + double taudw, + double &sxc, + double &v1xcup, + double &v1xcdw, + double &v2xcup, + double &v2xcdw, + double &v2xcud, + double &v3xcup, + double &v3xcdw, + double &vlaplxcup, + double &vlaplxcdw); + } // namespace XC_Functional_Libxc #endif // __LIBXC diff --git a/source/source_hamilt/module_xc/libxc_gga_wrap.cpp b/source/source_hamilt/module_xc/libxc_gga_wrap.cpp index bd86a46657d..e8391fb3f55 100644 --- a/source/source_hamilt/module_xc/libxc_gga_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_gga_wrap.cpp @@ -9,14 +9,12 @@ #include void XC_Functional_Libxc::gcxc_libxc( - const std::vector& func_id, + const std::vector& funcs, const double& rho, const double& grho, double& sxc, double& v1xc, - double& v2xc, - const double hybrid_alpha, - const double hse_omega) + double& v2xc) { sxc = 0.0; v1xc = 0.0; @@ -29,22 +27,35 @@ void XC_Functional_Libxc::gcxc_libxc( return; } - std::vector funcs = XC_Functional_Libxc::init_func( - /* func_id = */ func_id, - /* xc_polarized = */ XC_UNPOLARIZED, - /* hybrid_alpha = */ hybrid_alpha, - /* hse_omega = */ hse_omega); - - for (xc_func_type& func : funcs) + for (const xc_func_type& func : funcs) { double s = 0.0; double v1 = 0.0; double v2 = 0.0; - xc_gga_exc_vxc(&func, 1, &rho, &grho, &s, &v1, &v2); + xc_gga_exc_vxc(const_cast(&func), 1, &rho, &grho, &s, &v1, &v2); sxc += s * rho; v1xc += v1; v2xc += v2 * 2.0; } +} + +void XC_Functional_Libxc::gcxc_libxc( + const std::vector& func_id, + const double& rho, + const double& grho, + double& sxc, + double& v1xc, + double& v2xc, + const double hybrid_alpha, + const double hse_omega) +{ + std::vector funcs = XC_Functional_Libxc::init_func( + /* func_id = */ func_id, + /* xc_polarized = */ XC_UNPOLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); + + gcxc_libxc(funcs, rho, grho, sxc, v1xc, v2xc); XC_Functional_Libxc::finish_func(funcs); } // end subroutine gcxc_libxc @@ -52,7 +63,7 @@ void XC_Functional_Libxc::gcxc_libxc( void XC_Functional_Libxc::gcxc_spin_libxc( - const std::vector& func_id, + const std::vector& funcs, const double rhoup, const double rhodw, const ModuleBase::Vector3 gdr1, @@ -62,9 +73,7 @@ void XC_Functional_Libxc::gcxc_spin_libxc( double& v1xcdw, double& v2xcup, double& v2xcdw, - double& v2xcud, - const double hybrid_alpha, - const double hse_omega) + double& v2xcud) { sxc = 0.0; v1xcup = 0.0; @@ -75,13 +84,7 @@ void XC_Functional_Libxc::gcxc_spin_libxc( const std::array rho = {rhoup, rhodw}; const std::array grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()}; - std::vector funcs = XC_Functional_Libxc::init_func( - /* func_id = */ func_id, - /* xc_polarized = */ XC_POLARIZED, - /* hybrid_alpha = */ hybrid_alpha, - /* hse_omega = */ hse_omega); - - for (xc_func_type& func : funcs) + for (const xc_func_type& func : funcs) { if (func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) { @@ -104,7 +107,7 @@ void XC_Functional_Libxc::gcxc_spin_libxc( std::array v1xc = {0.0, 0.0}; std::array v2xc = {0.0, 0.0, 0.0}; // call Libxc function: xc_gga_exc_vxc - xc_gga_exc_vxc(&func, 1, rho.data(), grho.data(), &s, v1xc.data(), v2xc.data()); + xc_gga_exc_vxc(const_cast(&func), 1, rho.data(), grho.data(), &s, v1xc.data(), v2xc.data()); sxc += s * (rho[0] * sgn[0] + rho[1] * sgn[1]); v1xcup += v1xc[0] * sgn[0]; v1xcdw += v1xc[1] * sgn[1]; @@ -113,6 +116,32 @@ void XC_Functional_Libxc::gcxc_spin_libxc( v2xcdw += 2.0 * v2xc[2] * sgn[1]; } } +} + +void XC_Functional_Libxc::gcxc_spin_libxc( + const std::vector& func_id, + const double rhoup, + const double rhodw, + const ModuleBase::Vector3 gdr1, + const ModuleBase::Vector3 gdr2, + double& sxc, + double& v1xcup, + double& v1xcdw, + double& v2xcup, + double& v2xcdw, + double& v2xcud, + const double hybrid_alpha, + const double hse_omega) +{ + std::vector funcs = XC_Functional_Libxc::init_func( + /* func_id = */ func_id, + /* xc_polarized = */ XC_POLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); + + gcxc_spin_libxc(funcs, rhoup, rhodw, gdr1, gdr2, + sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud); + XC_Functional_Libxc::finish_func(funcs); } diff --git a/source/source_hamilt/module_xc/libxc_lda_wrap.cpp b/source/source_hamilt/module_xc/libxc_lda_wrap.cpp index 4b81e97cf7d..337dbb257fe 100644 --- a/source/source_hamilt/module_xc/libxc_lda_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_lda_wrap.cpp @@ -6,34 +6,42 @@ #endif void XC_Functional_Libxc::xc_spin_libxc( - const std::vector &func_id, + const std::vector& funcs, const double &rhoup, const double &rhodw, - double &exc, double &vxcup, double &vxcdw, - const double hybrid_alpha, - const double hse_omega) + double &exc, double &vxcup, double &vxcdw) { const std::vector rho_ud = {rhoup, rhodw}; exc = vxcup = vxcdw = 0.0; - std::vector funcs = XC_Functional_Libxc::init_func( - /* func_id = */ func_id, - /* xc_polarized = */ XC_POLARIZED, - /* hybrid_alpha = */ hybrid_alpha, - /* hse_omega = */ hse_omega); - - for(xc_func_type &func : funcs) + for(const xc_func_type &func : funcs) { double e = 0.0; std::vector vxc_ud(2); if( func.info->family == XC_FAMILY_LDA) { // call Libxc function: xc_lda_exc_vxc - xc_lda_exc_vxc( &func, 1, rho_ud.data(), &e, vxc_ud.data()); + xc_lda_exc_vxc( const_cast(&func), 1, rho_ud.data(), &e, vxc_ud.data()); } exc += e; vxcup += vxc_ud[0]; vxcdw += vxc_ud[1]; } +} + +void XC_Functional_Libxc::xc_spin_libxc( + const std::vector &func_id, + const double &rhoup, const double &rhodw, + double &exc, double &vxcup, double &vxcdw, + const double hybrid_alpha, + const double hse_omega) +{ + std::vector funcs = XC_Functional_Libxc::init_func( + /* func_id = */ func_id, + /* xc_polarized = */ XC_POLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); + + xc_spin_libxc(funcs, rhoup, rhodw, exc, vxcup, vxcdw); XC_Functional_Libxc::finish_func(funcs); } diff --git a/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp index 3fd7f34b53c..ab52d49b414 100644 --- a/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp @@ -14,7 +14,7 @@ //denote the polarized(nspin=1) or unpolarized(nspin=2) calculations, //definition can be found in xc.h from LIBXC void XC_Functional_Libxc::tau_xc( - const std::vector& func_id, + const std::vector& funcs, const double& rho, const double& grho, const double& lapl_rho, @@ -23,20 +23,13 @@ void XC_Functional_Libxc::tau_xc( double& v1xc, double& v2xc, double& v3xc, - double& vlaplxc, - const double& hybrid_alpha, - const double& hse_omega) + double& vlaplxc) { double s = 0.0; double v1 = 0.0; double v2 = 0.0; double v3 = 0.0; double vlapl_rho = 0.0; - std::vector funcs = XC_Functional_Libxc::init_func( - /* func_id = */ func_id, - /* xc_polarized = */ XC_UNPOLARIZED, - /* hybrid_alpha = */ hybrid_alpha, - /* hse_omega = */ hse_omega); sxc = 0.0; v1xc = 0.0; @@ -44,12 +37,14 @@ void XC_Functional_Libxc::tau_xc( v3xc = 0.0; vlaplxc = 0.0; - for (xc_func_type& func : funcs) + for (const xc_func_type& func : funcs) { - xc_mgga_exc_vxc(&func, 1, &rho, &grho, &lapl_rho, &atau, &s, &v1, &v2, &vlapl_rho, &v3); + xc_mgga_exc_vxc(const_cast(&func), 1, &rho, &grho, + &lapl_rho, &atau, &s, &v1, &v2, &vlapl_rho, &v3); #ifdef __EXX if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) { + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); s *= (1.0 - hybrid_alpha); v1 *= (1.0 - hybrid_alpha); v2 *= (1.0 - hybrid_alpha); @@ -63,6 +58,32 @@ void XC_Functional_Libxc::tau_xc( v3xc += v3; vlaplxc += vlapl_rho; } + + return; +} + +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) +{ + std::vector funcs = XC_Functional_Libxc::init_func( + /* func_id = */ func_id, + /* xc_polarized = */ XC_UNPOLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); + + tau_xc(funcs, rho, grho, lapl_rho, atau, sxc, v1xc, v2xc, v3xc, vlaplxc); + XC_Functional_Libxc::finish_func(funcs); return; @@ -70,7 +91,7 @@ void XC_Functional_Libxc::tau_xc( void XC_Functional_Libxc::tau_xc_spin( - const std::vector& func_id, + const std::vector& funcs, double rhoup, double rhodw, ModuleBase::Vector3 gdr1, @@ -88,9 +109,7 @@ void XC_Functional_Libxc::tau_xc_spin( double& v3xcup, double& v3xcdw, double& vlaplxcup, - double& vlaplxcdw, - const double& hybrid_alpha, - const double& hse_omega) + double& vlaplxcdw) { sxc = 0.0; v1xcup = 0.0; @@ -108,13 +127,7 @@ void XC_Functional_Libxc::tau_xc_spin( const std::array tau = {tauup, taudw}; const std::array lapl = {laplup, lapldw}; - std::vector funcs = XC_Functional_Libxc::init_func( - /* func_id = */ func_id, - /* xc_polarized = */ XC_POLARIZED, - /* hybrid_alpha = */ hybrid_alpha, - /* hse_omega = */ hse_omega); - - for (xc_func_type& func : funcs) + for (const xc_func_type& func : funcs) { if (func.info->family == XC_FAMILY_MGGA || func.info->family == XC_FAMILY_HYB_MGGA) { @@ -139,12 +152,14 @@ void XC_Functional_Libxc::tau_xc_spin( 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_out.data(), v3xc.data()); + xc_mgga_exc_vxc(const_cast(&func), 1, rho.data(), grho.data(), + lapl.data(), tau.data(), &s, 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) { + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); s *= (1.0 - hybrid_alpha); v1xc[0] *= (1.0 - hybrid_alpha); v1xc[1] *= (1.0 - hybrid_alpha); @@ -170,6 +185,40 @@ void XC_Functional_Libxc::tau_xc_spin( vlaplxcdw += vlapl_out[1] * sgn[1]; } } +} + +void XC_Functional_Libxc::tau_xc_spin( + const std::vector& func_id, + double rhoup, + double rhodw, + ModuleBase::Vector3 gdr1, + ModuleBase::Vector3 gdr2, + double laplup, + double lapldw, + double tauup, + double taudw, + double& sxc, + double& v1xcup, + double& v1xcdw, + double& v2xcup, + double& v2xcdw, + double& v2xcud, + double& v3xcup, + double& v3xcdw, + double& vlaplxcup, + double& vlaplxcdw, + const double& hybrid_alpha, + const double& hse_omega) +{ + std::vector funcs = XC_Functional_Libxc::init_func( + /* func_id = */ func_id, + /* xc_polarized = */ XC_POLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); + + tau_xc_spin(funcs, rhoup, rhodw, gdr1, gdr2, laplup, lapldw, tauup, taudw, + sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, + v3xcup, v3xcdw, vlaplxcup, vlaplxcdw); XC_Functional_Libxc::finish_func(funcs); } diff --git a/source/source_hamilt/module_xc/libxc_setup.cpp b/source/source_hamilt/module_xc/libxc_setup.cpp index 9f91c9276e9..25900d32dae 100644 --- a/source/source_hamilt/module_xc/libxc_setup.cpp +++ b/source/source_hamilt/module_xc/libxc_setup.cpp @@ -14,6 +14,7 @@ #include #include #include +#include bool not_supported_xc_with_laplacian(const std::string& xc_func_in) { @@ -282,12 +283,19 @@ const std::vector external_xc_func_ext_params(const int id) return {}; } +// Counter for the number of times init_func is called. Used by tests to +// verify that per-point evaluation paths do not re-initialize the libxc +// functional on every grid point. +static std::atomic g_libxc_init_count{0}; + std::vector XC_Functional_Libxc::init_func(const std::vector &func_id, const int xc_polarized, const double hybrid_alpha, const double hse_omega) { + ++g_libxc_init_count; + std::vector funcs; for (int id : func_id) { @@ -333,4 +341,9 @@ void XC_Functional_Libxc::finish_func(std::vector &funcs) } } +int XC_Functional_Libxc::get_and_reset_init_count() +{ + return g_libxc_init_count.exchange(0); +} + #endif diff --git a/source/source_hamilt/module_xc/test/test_xc2.cpp b/source/source_hamilt/module_xc/test/test_xc2.cpp index 06c1942a4f9..25cbbf59c14 100644 --- a/source/source_hamilt/module_xc/test/test_xc2.cpp +++ b/source/source_hamilt/module_xc/test/test_xc2.cpp @@ -508,4 +508,88 @@ TEST_F(XCTest_PZ_SPN_LibXC, set_xc_type) EXPECT_NEAR(v1_lda[i],v1_lda_ref[i],1.0e-8); EXPECT_NEAR(v2_lda[i],v2_lda_ref[i],1.0e-8); } +} + +// Verify that the already-initialized funcs overload of xc_spin_libxc +// produces the same result as the func_id overload for LDA. +TEST(XC_Libxc_Overload, xc_spin_libxc_funcs_matches_func_id) +{ + XC_Functional::set_xc_type("LDA_X+LDA_C_PZ"); + const std::vector& ids = XC_Functional::get_func_id(); + const double alpha = XC_Functional::get_hybrid_alpha(); + const double omega = XC_Functional::get_hse_omega(); + + const double rhoup = 0.5; + const double rhodw = 0.3; + + double e_old = 0.0, vup_old = 0.0, vdw_old = 0.0; + XC_Functional_Libxc::xc_spin_libxc(ids, rhoup, rhodw, e_old, vup_old, vdw_old, alpha, omega); + + std::vector funcs = XC_Functional_Libxc::init_func(ids, XC_POLARIZED, alpha, omega); + double e_new = 0.0, vup_new = 0.0, vdw_new = 0.0; + XC_Functional_Libxc::xc_spin_libxc(funcs, rhoup, rhodw, e_new, vup_new, vdw_new); + XC_Functional_Libxc::finish_func(funcs); + + EXPECT_NEAR(e_old, e_new, 1.0e-12); + EXPECT_NEAR(vup_old, vup_new, 1.0e-12); + EXPECT_NEAR(vdw_old, vdw_new, 1.0e-12); +} + +// Verify that the already-initialized funcs overload of gcxc_spin_libxc +// produces the same result as the func_id overload for GGA. +TEST(XC_Libxc_Overload, gcxc_spin_libxc_funcs_matches_func_id) +{ + XC_Functional::set_xc_type("GGA_X_PBE+GGA_C_PBE"); + const std::vector& ids = XC_Functional::get_func_id(); + const double alpha = XC_Functional::get_hybrid_alpha(); + const double omega = XC_Functional::get_hse_omega(); + + const double rhoup = 0.5; + const double rhodw = 0.3; + const ModuleBase::Vector3 gdr1(0.1, 0.2, 0.3); + const ModuleBase::Vector3 gdr2(0.05, 0.1, 0.15); + + double s_old = 0.0, v1up_old = 0.0, v1dw_old = 0.0, v2up_old = 0.0, v2dw_old = 0.0, v2ud_old = 0.0; + XC_Functional_Libxc::gcxc_spin_libxc(ids, rhoup, rhodw, gdr1, gdr2, + s_old, v1up_old, v1dw_old, v2up_old, v2dw_old, v2ud_old, alpha, omega); + + std::vector funcs = XC_Functional_Libxc::init_func(ids, XC_POLARIZED, alpha, omega); + double s_new = 0.0, v1up_new = 0.0, v1dw_new = 0.0, v2up_new = 0.0, v2dw_new = 0.0, v2ud_new = 0.0; + XC_Functional_Libxc::gcxc_spin_libxc(funcs, rhoup, rhodw, gdr1, gdr2, + s_new, v1up_new, v1dw_new, v2up_new, v2dw_new, v2ud_new); + XC_Functional_Libxc::finish_func(funcs); + + EXPECT_NEAR(s_old, s_new, 1.0e-12); + EXPECT_NEAR(v1up_old, v1up_new, 1.0e-12); + EXPECT_NEAR(v1dw_old, v1dw_new, 1.0e-12); + EXPECT_NEAR(v2up_old, v2up_new, 1.0e-12); + EXPECT_NEAR(v2dw_old, v2dw_new, 1.0e-12); + EXPECT_NEAR(v2ud_old, v2ud_new, 1.0e-12); +} + +// The func_id overload must initialize the functional exactly once; the +// funcs overload must initialize it zero times (the caller already did). +TEST(XC_Libxc_InitCount, per_point_overloads_do_not_reinit) +{ + XC_Functional::set_xc_type("LDA_X+LDA_C_PZ"); + const std::vector& ids = XC_Functional::get_func_id(); + const double alpha = 0.0; + const double omega = 0.0; + + XC_Functional_Libxc::get_and_reset_init_count(); + + double e = 0.0, vup = 0.0, vdw = 0.0; + XC_Functional_Libxc::xc_spin_libxc(ids, 0.5, 0.3, e, vup, vdw, alpha, omega); + const int count_after_old = XC_Functional_Libxc::get_and_reset_init_count(); + EXPECT_EQ(count_after_old, 1); + + std::vector funcs = XC_Functional_Libxc::init_func(ids, XC_POLARIZED, alpha, omega); + const int count_after_init = XC_Functional_Libxc::get_and_reset_init_count(); + EXPECT_EQ(count_after_init, 1); + + XC_Functional_Libxc::xc_spin_libxc(funcs, 0.5, 0.3, e, vup, vdw); + const int count_after_new = XC_Functional_Libxc::get_and_reset_init_count(); + EXPECT_EQ(count_after_new, 0); + + XC_Functional_Libxc::finish_func(funcs); } \ No newline at end of file diff --git a/source/source_hamilt/module_xc/test/test_xc6.cpp b/source/source_hamilt/module_xc/test/test_xc6.cpp index d4fac723326..fc0af59f175 100644 --- a/source/source_hamilt/module_xc/test/test_xc6.cpp +++ b/source/source_hamilt/module_xc/test/test_xc6.cpp @@ -159,3 +159,34 @@ TEST(XC_ScanL_Reference, tau_xc_wrapper_libxc) EXPECT_NE(v1xc, 0.0); EXPECT_NE(vlaplxc, 0.0); } + +// Verify that the already-initialized funcs overload of tau_xc +// produces the same result as the func_id overload for mGGA. +TEST(XC_Libxc_Overload, tau_xc_funcs_matches_func_id) +{ + XC_Functional::set_xc_type("MGGA_X_SCAN+MGGA_C_SCAN"); + const std::vector& ids = XC_Functional::get_func_id(); + const double alpha = XC_Functional::get_hybrid_alpha(); + const double omega = XC_Functional::get_hse_omega(); + + const double rho = 0.5; + const double grho = 0.01; + const double lapl = 0.0; + const double tau = 0.1; + + double s_old = 0.0, v1_old = 0.0, v2_old = 0.0, v3_old = 0.0, vl_old = 0.0; + XC_Functional_Libxc::tau_xc(ids, rho, grho, lapl, tau, + s_old, v1_old, v2_old, v3_old, vl_old, alpha, omega); + + std::vector funcs = XC_Functional_Libxc::init_func(ids, XC_UNPOLARIZED, alpha, omega); + double s_new = 0.0, v1_new = 0.0, v2_new = 0.0, v3_new = 0.0, vl_new = 0.0; + XC_Functional_Libxc::tau_xc(funcs, rho, grho, lapl, tau, + s_new, v1_new, v2_new, v3_new, vl_new); + XC_Functional_Libxc::finish_func(funcs); + + EXPECT_NEAR(s_old, s_new, 1.0e-12); + EXPECT_NEAR(v1_old, v1_new, 1.0e-12); + EXPECT_NEAR(v2_old, v2_new, 1.0e-12); + EXPECT_NEAR(v3_old, v3_new, 1.0e-12); + EXPECT_NEAR(vl_old, vl_new, 1.0e-12); +} diff --git a/source/source_hamilt/module_xc/xc_grad_kernel.cpp b/source/source_hamilt/module_xc/xc_grad_kernel.cpp index 78d80e4cec9..64939ffd31d 100644 --- a/source/source_hamilt/module_xc/xc_grad_kernel.cpp +++ b/source/source_hamilt/module_xc/xc_grad_kernel.cpp @@ -58,6 +58,20 @@ void gradcorr_xc_kernel(const GradCorrParams& params, double &local_etxcgc = etxcgc; #endif + // Initialize the libxc functional once per thread (not once per grid + // point) and reuse it for all points this thread handles. The + // polarized flag matches the per-point wrappers used below: + // nspin0==1 uses the unpolarized wrappers, otherwise the spin ones. +#ifdef __LIBXC + std::vector tls_funcs; + if (use_libxc) + { + const int xc_polarized = (nspin0 == 1) ? XC_UNPOLARIZED : XC_POLARIZED; + tls_funcs = XC_Functional_Libxc::init_func( + func_id, xc_polarized, hybrid_alpha_in, hse_omega_in); + } +#endif + double grho2a = 0.0; double grho2b = 0.0; double sxc = 0.0; @@ -102,12 +116,12 @@ void gradcorr_xc_kernel(const GradCorrParams& params, double vlaplxc = 0.0; double atau = chr->kin_r[0][ir]/2.0; double lapl_val = (!buf.lapl1.empty()) ? buf.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); + XC_Functional_Libxc::tau_xc( tls_funcs, arho, grho2a, lapl_val, atau, sxc, v1xc, v2xc, v3xc, vlaplxc); if(!buf.vlapl_arr1.empty()) buf.vlapl_arr1[ir] = vlaplxc; } else { - XC_Functional_Libxc::gcxc_libxc( func_id, arho, grho2a, sxc, v1xc, v2xc, hybrid_alpha_in, hse_omega_in); + XC_Functional_Libxc::gcxc_libxc( tls_funcs, arho, grho2a, sxc, v1xc, v2xc); } #endif } @@ -173,19 +187,18 @@ void gradcorr_xc_kernel(const GradCorrParams& params, double laplup_val = (!buf.lapl1.empty()) ? buf.lapl1[ir] : 0.0; double lapldw_val = (!buf.lapl2.empty()) ? buf.lapl2[ir] : 0.0; XC_Functional_Libxc::tau_xc_spin( - func_id, + tls_funcs, buf.rhotmp1[ir], buf.rhotmp2[ir], buf.gdr1[ir], buf.gdr2[ir], - laplup_val, lapldw_val, atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, vlaplxcup, vlaplxcdw, hybrid_alpha_in, hse_omega_in); + laplup_val, lapldw_val, atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, vlaplxcup, vlaplxcdw); if(!buf.vlapl_arr1.empty()) buf.vlapl_arr1[ir] = vlaplxcup; if(!buf.vlapl_arr2.empty()) buf.vlapl_arr2[ir] = vlaplxcdw; } else { XC_Functional_Libxc::gcxc_spin_libxc( - func_id, + tls_funcs, buf.rhotmp1[ir], buf.rhotmp2[ir], buf.gdr1[ir], buf.gdr2[ir], - sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, - hybrid_alpha_in, hse_omega_in); + sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud); } if(is_stress) { @@ -319,6 +332,11 @@ void gradcorr_xc_kernel(const GradCorrParams& params, } } } + + // Release the per-thread functional vector before the reduction. +#ifdef __LIBXC + XC_Functional_Libxc::finish_func(tls_funcs); +#endif #ifdef _OPENMP #pragma omp critical(xc_grad_reduce) { diff --git a/source/source_hamilt/module_xc/xc_pot.cpp b/source/source_hamilt/module_xc/xc_pot.cpp index f8b54bde2dd..82ed6f615f1 100644 --- a/source/source_hamilt/module_xc/xc_pot.cpp +++ b/source/source_hamilt/module_xc/xc_pot.cpp @@ -134,58 +134,79 @@ std::tuple XC_Functional::v_xc( else if(nspin == 4) { #ifdef _OPENMP -#pragma omp parallel for reduction(+:etxc) reduction(+:vtxc) -#endif - for(int ir = 0;irrho[1][ir],2) + pow(chr->rho[2][ir],2) + pow(chr->rho[3][ir],2) ); - double rhox = chr->rho[0][ir] + chr->rho_core[ir]; - double arhox = std::abs( rhox ); - - if ( arhox > vanishing_charge ) +#endif + // Initialize the libxc functional once per thread instead of once + // per grid point; the nspin==4 LDA path uses the spin-polarized + // wrapper. +#ifdef __LIBXC + std::vector tls_funcs; + if (use_libxc) { - double zeta = amag / arhox; - double exc = 0.0; - double vxc[2]; + tls_funcs = XC_Functional_Libxc::init_func( + XC_Functional::get_func_id(), XC_POLARIZED, hybrid_alpha, hse_omega); + } +#endif +#ifdef _OPENMP +#pragma omp for +#endif + for(int ir = 0;irrho[1][ir],2) + pow(chr->rho[2][ir],2) + pow(chr->rho[3][ir],2) ); + double rhox = chr->rho[0][ir] + chr->rho_core[ir]; + double arhox = std::abs( rhox ); - if ( std::abs( zeta ) > 1.0 ) + if ( arhox > vanishing_charge ) { - zeta = (zeta > 0.0) ? 1.0 : (-1.0); - } + double zeta = amag / arhox; + double exc = 0.0; + double vxc[2]; - if(use_libxc) - { + if ( std::abs( zeta ) > 1.0 ) + { + zeta = (zeta > 0.0) ? 1.0 : (-1.0); + } + + if(use_libxc) + { #ifdef __LIBXC - double rhoup = arhox * (1.0+zeta) / 2.0; - double rhodw = arhox * (1.0-zeta) / 2.0; - XC_Functional_Libxc::xc_spin_libxc(XC_Functional::get_func_id(), rhoup, rhodw, exc, vxc[0], vxc[1], hybrid_alpha, hse_omega); + double rhoup = arhox * (1.0+zeta) / 2.0; + double rhodw = arhox * (1.0-zeta) / 2.0; + XC_Functional_Libxc::xc_spin_libxc(tls_funcs, rhoup, rhodw, exc, vxc[0], vxc[1]); #else - ModuleBase::WARNING_QUIT("v_xc", "compile with LIBXC"); + ModuleBase::WARNING_QUIT("v_xc", "compile with LIBXC"); #endif - } - else - { - double rhoup = arhox * (1.0+zeta) / 2.0; - double rhodw = arhox * (1.0-zeta) / 2.0; - XC_Functional::xc_spin(arhox, zeta, exc, vxc[0], vxc[1]); - } + } + else + { + double rhoup = arhox * (1.0+zeta) / 2.0; + double rhodw = arhox * (1.0-zeta) / 2.0; + XC_Functional::xc_spin(arhox, zeta, exc, vxc[0], vxc[1]); + } - etxc += e2 * exc * rhox; + etxc += e2 * exc * rhox; - v(0, ir) = e2*( 0.5 * ( vxc[0] + vxc[1]) ); - vtxc += v(0,ir) * chr->rho[0][ir]; + v(0, ir) = e2*( 0.5 * ( vxc[0] + vxc[1]) ); + vtxc += v(0,ir) * chr->rho[0][ir]; - double vs = 0.5 * ( vxc[0] - vxc[1] ); - if ( amag > vanishing_charge ) - { - for(int ipol = 1;ipol< 4;ipol++) + double vs = 0.5 * ( vxc[0] - vxc[1] ); + if ( amag > vanishing_charge ) { - v(ipol, ir) = e2 * vs * chr->rho[ipol][ir] / amag; - vtxc += v(ipol,ir) * chr->rho[ipol][ir]; + for(int ipol = 1;ipol< 4;ipol++) + { + v(ipol, ir) = e2 * vs * chr->rho[ipol][ir] / amag; + vtxc += v(ipol,ir) * chr->rho[ipol][ir]; + } } } } +#ifdef __LIBXC + XC_Functional_Libxc::finish_func(tls_funcs); +#endif +#ifdef _OPENMP } +#endif } // energy terms, local-density contributions From e0007e8bad84c1420748245eac6198f559ad8bc0 Mon Sep 17 00:00:00 2001 From: Critsium Date: Sat, 10 Oct 2026 17:47:31 +0800 Subject: [PATCH 2/2] fix(module_xc): address review of libxc per-thread init (#8082) - tau_xc/tau_xc_spin funcs overloads take hybrid_alpha explicitly instead of reading the XC_Functional global for SCAN0 scaling - revert the unreachable use_libxc branch change in v_xc (nspin==4) and the now-unused xc_spin_libxc funcs overload - drop the init_func call counter from the production path - only init libxc in gradcorr when a libxc branch runs (nspin0 != 1 || is_stress) - gcxc_libxc func_id overload returns before init_func for negligible rho/grho, as before - remove const_cast: Libxc >= 5.1.7 takes const xc_func_type* - replace self-comparison tests with a libxc gradcorr regression test (references from the per-point implementation) and a SCAN0 test checking the hybrid_alpha argument is used Co-Authored-By: Claude Opus 5.5 (1M context) --- source/source_hamilt/module_xc/libxc_abacus.h | 30 ++--- .../module_xc/libxc_gga_wrap.cpp | 28 ++++- .../module_xc/libxc_lda_wrap.cpp | 32 ++--- .../module_xc/libxc_mgga_wrap.cpp | 16 +-- .../source_hamilt/module_xc/libxc_setup.cpp | 13 -- .../source_hamilt/module_xc/test/test_xc2.cpp | 84 ------------- .../source_hamilt/module_xc/test/test_xc3.cpp | 111 ++++++++++++++++++ .../source_hamilt/module_xc/test/test_xc6.cpp | 88 +++++++++++--- .../module_xc/xc_grad_kernel.cpp | 9 +- source/source_hamilt/module_xc/xc_pot.cpp | 93 ++++++--------- 10 files changed, 275 insertions(+), 229 deletions(-) diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index 4110c414185..d724837db00 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/source/source_hamilt/module_xc/libxc_abacus.h @@ -73,13 +73,6 @@ namespace XC_Functional_Libxc extern void finish_func(std::vector &funcs); - // Returns the number of times init_func has been called since the last - // call to this function, and resets the counter to zero. Intended for - // unit tests that assert per-point evaluation paths do not re-initialize - // the libxc functional on every grid point. - extern int get_and_reset_init_count(); - - //------------------- // libxc_pot.cpp //------------------- @@ -272,17 +265,6 @@ namespace XC_Functional_Libxc const double hybrid_alpha, const double hse_omega); - // Overload accepting an already-initialized functional vector. The caller - // is responsible for init_func/finish_func; this overload does not touch - // the lifetime of funcs. Useful for per-thread reuse inside OpenMP loops. - extern void xc_spin_libxc( - const std::vector& funcs, - const double &rhoup, - const double &rhodw, - double &exc, - double &vxcup, - double &vxcdw); - //------------------- // libxc_gga_wrap.cpp @@ -299,7 +281,9 @@ namespace XC_Functional_Libxc const double hybrid_alpha, const double hse_omega); - // Overload accepting an already-initialized functional vector. + // Overload accepting an already-initialized functional vector. The caller + // is responsible for init_func/finish_func; this overload does not touch + // the lifetime of funcs. Useful for per-thread reuse inside OpenMP loops. extern void gcxc_libxc( const std::vector& funcs, const double &rho, @@ -359,6 +343,7 @@ namespace XC_Functional_Libxc const double &hse_omega); // Overload accepting an already-initialized functional vector. + // hybrid_alpha scales the semilocal SCAN exchange for SCAN0. extern void tau_xc( const std::vector& funcs, const double &rho, @@ -369,7 +354,8 @@ namespace XC_Functional_Libxc double &v1xc, double &v2xc, double &v3xc, - double &vlaplxc); + double &vlaplxc, + const double hybrid_alpha); extern void tau_xc_spin( const std::vector &func_id, @@ -395,6 +381,7 @@ namespace XC_Functional_Libxc const double &hse_omega); // Overload accepting an already-initialized functional vector. + // hybrid_alpha scales the semilocal SCAN exchange for SCAN0. extern void tau_xc_spin( const std::vector& funcs, double rhoup, @@ -414,7 +401,8 @@ namespace XC_Functional_Libxc double &v3xcup, double &v3xcdw, double &vlaplxcup, - double &vlaplxcdw); + double &vlaplxcdw, + const double hybrid_alpha); } // namespace XC_Functional_Libxc diff --git a/source/source_hamilt/module_xc/libxc_gga_wrap.cpp b/source/source_hamilt/module_xc/libxc_gga_wrap.cpp index e8391fb3f55..29a259f7dfe 100644 --- a/source/source_hamilt/module_xc/libxc_gga_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_gga_wrap.cpp @@ -8,6 +8,17 @@ #include #include +namespace +{ +// gcxc_libxc returns zero below these thresholds. +bool gcxc_is_negligible(const double rho, const double grho) +{ + constexpr double small = 1.e-6; + constexpr double smallg = 1.e-10; + return rho <= small || grho < smallg; +} +} // namespace + void XC_Functional_Libxc::gcxc_libxc( const std::vector& funcs, const double& rho, @@ -20,9 +31,7 @@ void XC_Functional_Libxc::gcxc_libxc( v1xc = 0.0; v2xc = 0.0; - constexpr double small = 1.e-6; - constexpr double smallg = 1.e-10; - if (rho <= small || grho < smallg) + if (gcxc_is_negligible(rho, grho)) { return; } @@ -32,7 +41,7 @@ void XC_Functional_Libxc::gcxc_libxc( double s = 0.0; double v1 = 0.0; double v2 = 0.0; - xc_gga_exc_vxc(const_cast(&func), 1, &rho, &grho, &s, &v1, &v2); + xc_gga_exc_vxc(&func, 1, &rho, &grho, &s, &v1, &v2); sxc += s * rho; v1xc += v1; v2xc += v2 * 2.0; @@ -49,6 +58,15 @@ void XC_Functional_Libxc::gcxc_libxc( const double hybrid_alpha, const double hse_omega) { + // Skip the libxc setup entirely for points the evaluation would zero out. + if (gcxc_is_negligible(rho, grho)) + { + sxc = 0.0; + v1xc = 0.0; + v2xc = 0.0; + return; + } + std::vector funcs = XC_Functional_Libxc::init_func( /* func_id = */ func_id, /* xc_polarized = */ XC_UNPOLARIZED, @@ -107,7 +125,7 @@ void XC_Functional_Libxc::gcxc_spin_libxc( std::array v1xc = {0.0, 0.0}; std::array v2xc = {0.0, 0.0, 0.0}; // call Libxc function: xc_gga_exc_vxc - xc_gga_exc_vxc(const_cast(&func), 1, rho.data(), grho.data(), &s, v1xc.data(), v2xc.data()); + xc_gga_exc_vxc(&func, 1, rho.data(), grho.data(), &s, v1xc.data(), v2xc.data()); sxc += s * (rho[0] * sgn[0] + rho[1] * sgn[1]); v1xcup += v1xc[0] * sgn[0]; v1xcdw += v1xc[1] * sgn[1]; diff --git a/source/source_hamilt/module_xc/libxc_lda_wrap.cpp b/source/source_hamilt/module_xc/libxc_lda_wrap.cpp index 337dbb257fe..4b81e97cf7d 100644 --- a/source/source_hamilt/module_xc/libxc_lda_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_lda_wrap.cpp @@ -6,42 +6,34 @@ #endif void XC_Functional_Libxc::xc_spin_libxc( - const std::vector& funcs, + const std::vector &func_id, const double &rhoup, const double &rhodw, - double &exc, double &vxcup, double &vxcdw) + double &exc, double &vxcup, double &vxcdw, + const double hybrid_alpha, + const double hse_omega) { const std::vector rho_ud = {rhoup, rhodw}; exc = vxcup = vxcdw = 0.0; - for(const xc_func_type &func : funcs) + std::vector funcs = XC_Functional_Libxc::init_func( + /* func_id = */ func_id, + /* xc_polarized = */ XC_POLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); + + for(xc_func_type &func : funcs) { double e = 0.0; std::vector vxc_ud(2); if( func.info->family == XC_FAMILY_LDA) { // call Libxc function: xc_lda_exc_vxc - xc_lda_exc_vxc( const_cast(&func), 1, rho_ud.data(), &e, vxc_ud.data()); + xc_lda_exc_vxc( &func, 1, rho_ud.data(), &e, vxc_ud.data()); } exc += e; vxcup += vxc_ud[0]; vxcdw += vxc_ud[1]; } -} - -void XC_Functional_Libxc::xc_spin_libxc( - const std::vector &func_id, - const double &rhoup, const double &rhodw, - double &exc, double &vxcup, double &vxcdw, - const double hybrid_alpha, - const double hse_omega) -{ - std::vector funcs = XC_Functional_Libxc::init_func( - /* func_id = */ func_id, - /* xc_polarized = */ XC_POLARIZED, - /* hybrid_alpha = */ hybrid_alpha, - /* hse_omega = */ hse_omega); - - xc_spin_libxc(funcs, rhoup, rhodw, exc, vxcup, vxcdw); XC_Functional_Libxc::finish_func(funcs); } diff --git a/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp index ab52d49b414..e6031ba5fa6 100644 --- a/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp @@ -23,7 +23,8 @@ void XC_Functional_Libxc::tau_xc( double& v1xc, double& v2xc, double& v3xc, - double& vlaplxc) + double& vlaplxc, + const double hybrid_alpha) { double s = 0.0; double v1 = 0.0; @@ -39,12 +40,11 @@ void XC_Functional_Libxc::tau_xc( for (const xc_func_type& func : funcs) { - xc_mgga_exc_vxc(const_cast(&func), 1, &rho, &grho, + xc_mgga_exc_vxc(&func, 1, &rho, &grho, &lapl_rho, &atau, &s, &v1, &v2, &vlapl_rho, &v3); #ifdef __EXX if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) { - const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); s *= (1.0 - hybrid_alpha); v1 *= (1.0 - hybrid_alpha); v2 *= (1.0 - hybrid_alpha); @@ -82,7 +82,7 @@ void XC_Functional_Libxc::tau_xc( /* hybrid_alpha = */ hybrid_alpha, /* hse_omega = */ hse_omega); - tau_xc(funcs, rho, grho, lapl_rho, atau, sxc, v1xc, v2xc, v3xc, vlaplxc); + tau_xc(funcs, rho, grho, lapl_rho, atau, sxc, v1xc, v2xc, v3xc, vlaplxc, hybrid_alpha); XC_Functional_Libxc::finish_func(funcs); @@ -109,7 +109,8 @@ void XC_Functional_Libxc::tau_xc_spin( double& v3xcup, double& v3xcdw, double& vlaplxcup, - double& vlaplxcdw) + double& vlaplxcdw, + const double hybrid_alpha) { sxc = 0.0; v1xcup = 0.0; @@ -152,14 +153,13 @@ void XC_Functional_Libxc::tau_xc_spin( 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(const_cast(&func), 1, rho.data(), grho.data(), + xc_mgga_exc_vxc(&func, 1, rho.data(), grho.data(), lapl.data(), tau.data(), &s, 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) { - const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); s *= (1.0 - hybrid_alpha); v1xc[0] *= (1.0 - hybrid_alpha); v1xc[1] *= (1.0 - hybrid_alpha); @@ -218,7 +218,7 @@ void XC_Functional_Libxc::tau_xc_spin( tau_xc_spin(funcs, rhoup, rhodw, gdr1, gdr2, laplup, lapldw, tauup, taudw, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, - v3xcup, v3xcdw, vlaplxcup, vlaplxcdw); + v3xcup, v3xcdw, vlaplxcup, vlaplxcdw, hybrid_alpha); XC_Functional_Libxc::finish_func(funcs); } diff --git a/source/source_hamilt/module_xc/libxc_setup.cpp b/source/source_hamilt/module_xc/libxc_setup.cpp index 25900d32dae..9f91c9276e9 100644 --- a/source/source_hamilt/module_xc/libxc_setup.cpp +++ b/source/source_hamilt/module_xc/libxc_setup.cpp @@ -14,7 +14,6 @@ #include #include #include -#include bool not_supported_xc_with_laplacian(const std::string& xc_func_in) { @@ -283,19 +282,12 @@ const std::vector external_xc_func_ext_params(const int id) return {}; } -// Counter for the number of times init_func is called. Used by tests to -// verify that per-point evaluation paths do not re-initialize the libxc -// functional on every grid point. -static std::atomic g_libxc_init_count{0}; - std::vector XC_Functional_Libxc::init_func(const std::vector &func_id, const int xc_polarized, const double hybrid_alpha, const double hse_omega) { - ++g_libxc_init_count; - std::vector funcs; for (int id : func_id) { @@ -341,9 +333,4 @@ void XC_Functional_Libxc::finish_func(std::vector &funcs) } } -int XC_Functional_Libxc::get_and_reset_init_count() -{ - return g_libxc_init_count.exchange(0); -} - #endif diff --git a/source/source_hamilt/module_xc/test/test_xc2.cpp b/source/source_hamilt/module_xc/test/test_xc2.cpp index 25cbbf59c14..06c1942a4f9 100644 --- a/source/source_hamilt/module_xc/test/test_xc2.cpp +++ b/source/source_hamilt/module_xc/test/test_xc2.cpp @@ -508,88 +508,4 @@ TEST_F(XCTest_PZ_SPN_LibXC, set_xc_type) EXPECT_NEAR(v1_lda[i],v1_lda_ref[i],1.0e-8); EXPECT_NEAR(v2_lda[i],v2_lda_ref[i],1.0e-8); } -} - -// Verify that the already-initialized funcs overload of xc_spin_libxc -// produces the same result as the func_id overload for LDA. -TEST(XC_Libxc_Overload, xc_spin_libxc_funcs_matches_func_id) -{ - XC_Functional::set_xc_type("LDA_X+LDA_C_PZ"); - const std::vector& ids = XC_Functional::get_func_id(); - const double alpha = XC_Functional::get_hybrid_alpha(); - const double omega = XC_Functional::get_hse_omega(); - - const double rhoup = 0.5; - const double rhodw = 0.3; - - double e_old = 0.0, vup_old = 0.0, vdw_old = 0.0; - XC_Functional_Libxc::xc_spin_libxc(ids, rhoup, rhodw, e_old, vup_old, vdw_old, alpha, omega); - - std::vector funcs = XC_Functional_Libxc::init_func(ids, XC_POLARIZED, alpha, omega); - double e_new = 0.0, vup_new = 0.0, vdw_new = 0.0; - XC_Functional_Libxc::xc_spin_libxc(funcs, rhoup, rhodw, e_new, vup_new, vdw_new); - XC_Functional_Libxc::finish_func(funcs); - - EXPECT_NEAR(e_old, e_new, 1.0e-12); - EXPECT_NEAR(vup_old, vup_new, 1.0e-12); - EXPECT_NEAR(vdw_old, vdw_new, 1.0e-12); -} - -// Verify that the already-initialized funcs overload of gcxc_spin_libxc -// produces the same result as the func_id overload for GGA. -TEST(XC_Libxc_Overload, gcxc_spin_libxc_funcs_matches_func_id) -{ - XC_Functional::set_xc_type("GGA_X_PBE+GGA_C_PBE"); - const std::vector& ids = XC_Functional::get_func_id(); - const double alpha = XC_Functional::get_hybrid_alpha(); - const double omega = XC_Functional::get_hse_omega(); - - const double rhoup = 0.5; - const double rhodw = 0.3; - const ModuleBase::Vector3 gdr1(0.1, 0.2, 0.3); - const ModuleBase::Vector3 gdr2(0.05, 0.1, 0.15); - - double s_old = 0.0, v1up_old = 0.0, v1dw_old = 0.0, v2up_old = 0.0, v2dw_old = 0.0, v2ud_old = 0.0; - XC_Functional_Libxc::gcxc_spin_libxc(ids, rhoup, rhodw, gdr1, gdr2, - s_old, v1up_old, v1dw_old, v2up_old, v2dw_old, v2ud_old, alpha, omega); - - std::vector funcs = XC_Functional_Libxc::init_func(ids, XC_POLARIZED, alpha, omega); - double s_new = 0.0, v1up_new = 0.0, v1dw_new = 0.0, v2up_new = 0.0, v2dw_new = 0.0, v2ud_new = 0.0; - XC_Functional_Libxc::gcxc_spin_libxc(funcs, rhoup, rhodw, gdr1, gdr2, - s_new, v1up_new, v1dw_new, v2up_new, v2dw_new, v2ud_new); - XC_Functional_Libxc::finish_func(funcs); - - EXPECT_NEAR(s_old, s_new, 1.0e-12); - EXPECT_NEAR(v1up_old, v1up_new, 1.0e-12); - EXPECT_NEAR(v1dw_old, v1dw_new, 1.0e-12); - EXPECT_NEAR(v2up_old, v2up_new, 1.0e-12); - EXPECT_NEAR(v2dw_old, v2dw_new, 1.0e-12); - EXPECT_NEAR(v2ud_old, v2ud_new, 1.0e-12); -} - -// The func_id overload must initialize the functional exactly once; the -// funcs overload must initialize it zero times (the caller already did). -TEST(XC_Libxc_InitCount, per_point_overloads_do_not_reinit) -{ - XC_Functional::set_xc_type("LDA_X+LDA_C_PZ"); - const std::vector& ids = XC_Functional::get_func_id(); - const double alpha = 0.0; - const double omega = 0.0; - - XC_Functional_Libxc::get_and_reset_init_count(); - - double e = 0.0, vup = 0.0, vdw = 0.0; - XC_Functional_Libxc::xc_spin_libxc(ids, 0.5, 0.3, e, vup, vdw, alpha, omega); - const int count_after_old = XC_Functional_Libxc::get_and_reset_init_count(); - EXPECT_EQ(count_after_old, 1); - - std::vector funcs = XC_Functional_Libxc::init_func(ids, XC_POLARIZED, alpha, omega); - const int count_after_init = XC_Functional_Libxc::get_and_reset_init_count(); - EXPECT_EQ(count_after_init, 1); - - XC_Functional_Libxc::xc_spin_libxc(funcs, 0.5, 0.3, e, vup, vdw); - const int count_after_new = XC_Functional_Libxc::get_and_reset_init_count(); - EXPECT_EQ(count_after_new, 0); - - XC_Functional_Libxc::finish_func(funcs); } \ No newline at end of file diff --git a/source/source_hamilt/module_xc/test/test_xc3.cpp b/source/source_hamilt/module_xc/test/test_xc3.cpp index cd7d0e52e34..0171e7635b4 100644 --- a/source/source_hamilt/module_xc/test/test_xc3.cpp +++ b/source/source_hamilt/module_xc/test/test_xc3.cpp @@ -289,6 +289,117 @@ TEST_F(XCTest_GRADCORR_HF, no_local_functional) } } +// gradcorr with a libxc GGA initializes the functional once per OpenMP thread +// and reuses it on every grid point. The references were produced by the +// previous implementation, which initialized the functional per grid point. +// nspin=1 reaches libxc only for stress; nspin=2 uses it for v and stress. +class XCTest_GRADCORR_Libxc : public XCTest +{ + protected: + + double et1 = 0; + double vt1 = 0; + ModuleBase::matrix v1; + std::vector stress1; + + double et2 = 0; + double vt2 = 0; + ModuleBase::matrix v2; + std::vector stress2; + + void SetUp() + { + const int nspin1 = 1; + const int nspin2 = 2; + const bool domag = false; + const bool domag_z = false; + + ModulePW::PW_Basis rhopw; + UnitCell ucell; + Charge chr; + + rhopw.nrxx = 5; + rhopw.npw = 5; + rhopw.nmaxgr = 5; + rhopw.gcar = new ModuleBase::Vector3 [5]; + + ucell.tpiba = 1; + ucell.magnet.lsign_ = true; + unitcell::cal_ux(ucell, 4); + + chr.rho = new double*[2]; + chr.rho[0] = new double[5]; + chr.rho[1] = new double[5]; + chr.rhog = new std::complex*[2]; + chr.rhog[0] = new std::complex[5]; + chr.rhog[1] = new std::complex[5]; + + chr.rho_core = new double[5]; + chr.rhog_core = new std::complex[5]; + + for(int i=0;i<5;i++) + { + chr.rho[0][i] = double(i); + chr.rho[1][i] = 0.1*double(i); + chr.rhog[0][i] = chr.rho[0][i]; + chr.rhog[1][i] = chr.rho[1][i]; + chr.rho_core[i] = 0; + chr.rhog_core[i] = 0; + rhopw.gcar[i]= 1; + } + + v1.create(1,5); + v1.zero_out(); + v2.create(2,5); + v2.zero_out(); + + XC_Functional::set_xc_type("GGA_X_PBE+GGA_C_PBE"); + + const int gga_grad = 0; + const double hybrid_alpha = 0.0; + const double hse_omega = 0.0; + XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,true, nspin1,domag,domag_z,gga_grad,hybrid_alpha,hse_omega); + + XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,false,nspin2,domag,domag_z,gga_grad,hybrid_alpha,hse_omega); + XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,true, nspin2,domag,domag_z,gga_grad,hybrid_alpha,hse_omega); + } +}; + +TEST_F(XCTest_GRADCORR_Libxc, gga_pbe) +{ + EXPECT_EQ(XC_Functional::get_func_type(), 2); + + const double s1 = -0.0253603083703; + std::vector stress1_ref = {s1,0,0,s1,s1,0,s1,s1,s1}; + ASSERT_EQ(stress1.size(), 9); + for(int i=0;i<9;i++) + { + EXPECT_NEAR(stress1[i],stress1_ref[i],1.0e-10); + } + + // the libxc spin wrapper returns the full functional, not only the + // gradient correction, so et2/vt2/v2 differ from the built-in PBE ones + const double et2_ref = -28.9783818855; + const double vt2_ref = -38.154198697; + std::vector v2_ref1 = {0,-2.56088553194,-3.21933929441,-3.67877304186,-4.04360433532}; + std::vector v2_ref2 = {0,-1.39427647251,-1.73902789949,-1.97505893723,-2.16036800258}; + EXPECT_NEAR(et2,et2_ref,1.0e-8); + EXPECT_NEAR(vt2,vt2_ref,1.0e-8); + for(int i=0;i<5;i++) + { + EXPECT_NEAR(v2(0,i),v2_ref1[i],1.0e-8); + EXPECT_NEAR(v2(1,i),v2_ref2[i],1.0e-8); + } + + const double s2 = -0.0280738149619; + std::vector stress2_ref = {s2,0,0,s2,s2,0,s2,s2,s2}; + ASSERT_EQ(stress2.size(), 9); + for(int i=0;i<9;i++) + { + EXPECT_NEAR(stress2[i],stress2_ref[i],1.0e-10); + } +} + class XCTest_GRADWFC : public XCTest { protected: diff --git a/source/source_hamilt/module_xc/test/test_xc6.cpp b/source/source_hamilt/module_xc/test/test_xc6.cpp index fc0af59f175..345e18d6e47 100644 --- a/source/source_hamilt/module_xc/test/test_xc6.cpp +++ b/source/source_hamilt/module_xc/test/test_xc6.cpp @@ -2,6 +2,7 @@ #include "../libxc_abacus.h" #include "gtest/gtest.h" #include "xctest.h" +#include #include #include #include @@ -160,33 +161,84 @@ TEST(XC_ScanL_Reference, tau_xc_wrapper_libxc) EXPECT_NE(vlaplxc, 0.0); } -// Verify that the already-initialized funcs overload of tau_xc -// produces the same result as the func_id overload for mGGA. -TEST(XC_Libxc_Overload, tau_xc_funcs_matches_func_id) +#ifdef __EXX +// SCAN0 scales the semilocal SCAN exchange by (1 - hybrid_alpha). The +// wrappers must use the hybrid_alpha they are given, not the global +// XC_Functional value, so the two are deliberately set to different values. +class XCTest_SCAN0_Alpha : public testing::Test { - XC_Functional::set_xc_type("MGGA_X_SCAN+MGGA_C_SCAN"); - const std::vector& ids = XC_Functional::get_func_id(); - const double alpha = XC_Functional::get_hybrid_alpha(); - const double omega = XC_Functional::get_hse_omega(); + protected: + const std::vector x_only = {XC_MGGA_X_SCAN}; + const double global_alpha = 0.25; + const double hybrid_alpha = 0.4; + const double hse_omega = 0.0; + + void SetUp() override + { + XC_Functional::set_xc_type("SCAN0"); + XC_Functional::set_hybrid_alpha(global_alpha); + } +}; +TEST_F(XCTest_SCAN0_Alpha, tau_xc) +{ + const double no_alpha = 0.0; + const double scale = 1.0 - hybrid_alpha; const double rho = 0.5; const double grho = 0.01; const double lapl = 0.0; const double tau = 0.1; - double s_old = 0.0, v1_old = 0.0, v2_old = 0.0, v3_old = 0.0, vl_old = 0.0; - XC_Functional_Libxc::tau_xc(ids, rho, grho, lapl, tau, - s_old, v1_old, v2_old, v3_old, vl_old, alpha, omega); + std::array full = {0.0, 0.0, 0.0, 0.0, 0.0}; + XC_Functional_Libxc::tau_xc(x_only, rho, grho, lapl, tau, + full[0], full[1], full[2], full[3], full[4], no_alpha, hse_omega); + + std::array scaled = {0.0, 0.0, 0.0, 0.0, 0.0}; + XC_Functional_Libxc::tau_xc(x_only, rho, grho, lapl, tau, + scaled[0], scaled[1], scaled[2], scaled[3], scaled[4], hybrid_alpha, hse_omega); - std::vector funcs = XC_Functional_Libxc::init_func(ids, XC_UNPOLARIZED, alpha, omega); - double s_new = 0.0, v1_new = 0.0, v2_new = 0.0, v3_new = 0.0, vl_new = 0.0; + std::vector funcs = XC_Functional_Libxc::init_func(x_only, XC_UNPOLARIZED, hybrid_alpha, hse_omega); + std::array scaled_funcs = {0.0, 0.0, 0.0, 0.0, 0.0}; XC_Functional_Libxc::tau_xc(funcs, rho, grho, lapl, tau, - s_new, v1_new, v2_new, v3_new, vl_new); + scaled_funcs[0], scaled_funcs[1], scaled_funcs[2], scaled_funcs[3], scaled_funcs[4], + hybrid_alpha); XC_Functional_Libxc::finish_func(funcs); - EXPECT_NEAR(s_old, s_new, 1.0e-12); - EXPECT_NEAR(v1_old, v1_new, 1.0e-12); - EXPECT_NEAR(v2_old, v2_new, 1.0e-12); - EXPECT_NEAR(v3_old, v3_new, 1.0e-12); - EXPECT_NEAR(vl_old, vl_new, 1.0e-12); + EXPECT_NE(full[0], 0.0); + for (int i = 0; i < 5; i++) + { + EXPECT_NEAR(scaled[i], scale * full[i], 1.0e-12); + EXPECT_NEAR(scaled_funcs[i], scale * full[i], 1.0e-12); + } +} + +TEST_F(XCTest_SCAN0_Alpha, tau_xc_spin) +{ + const double no_alpha = 0.0; + const double scale = 1.0 - hybrid_alpha; + const double rhoup = 0.3; + const double rhodw = 0.2; + const ModuleBase::Vector3 gdr1(0.05, 0.02, 0.01); + const ModuleBase::Vector3 gdr2(0.03, 0.01, 0.02); + const double laplup = 0.0; + const double lapldw = 0.0; + const double tauup = 0.06; + const double taudw = 0.04; + + std::array full = {}; + XC_Functional_Libxc::tau_xc_spin(x_only, rhoup, rhodw, gdr1, gdr2, laplup, lapldw, tauup, taudw, + full[0], full[1], full[2], full[3], full[4], full[5], + full[6], full[7], full[8], full[9], no_alpha, hse_omega); + + std::array scaled = {}; + XC_Functional_Libxc::tau_xc_spin(x_only, rhoup, rhodw, gdr1, gdr2, laplup, lapldw, tauup, taudw, + scaled[0], scaled[1], scaled[2], scaled[3], scaled[4], scaled[5], + scaled[6], scaled[7], scaled[8], scaled[9], hybrid_alpha, hse_omega); + + EXPECT_NE(full[0], 0.0); + for (int i = 0; i < 10; i++) + { + EXPECT_NEAR(scaled[i], scale * full[i], 1.0e-12); + } } +#endif diff --git a/source/source_hamilt/module_xc/xc_grad_kernel.cpp b/source/source_hamilt/module_xc/xc_grad_kernel.cpp index 64939ffd31d..09bda321d6b 100644 --- a/source/source_hamilt/module_xc/xc_grad_kernel.cpp +++ b/source/source_hamilt/module_xc/xc_grad_kernel.cpp @@ -62,9 +62,11 @@ void gradcorr_xc_kernel(const GradCorrParams& params, // point) and reuse it for all points this thread handles. The // polarized flag matches the per-point wrappers used below: // nspin0==1 uses the unpolarized wrappers, otherwise the spin ones. + // The nspin0==1 non-stress branch uses XC_Functional::gcxc, not libxc. #ifdef __LIBXC std::vector tls_funcs; - if (use_libxc) + const bool need_libxc_funcs = use_libxc && (nspin0 != 1 || is_stress); + if (need_libxc_funcs) { const int xc_polarized = (nspin0 == 1) ? XC_UNPOLARIZED : XC_POLARIZED; tls_funcs = XC_Functional_Libxc::init_func( @@ -116,7 +118,7 @@ void gradcorr_xc_kernel(const GradCorrParams& params, double vlaplxc = 0.0; double atau = chr->kin_r[0][ir]/2.0; double lapl_val = (!buf.lapl1.empty()) ? buf.lapl1[ir] : 0.0; - XC_Functional_Libxc::tau_xc( tls_funcs, arho, grho2a, lapl_val, atau, sxc, v1xc, v2xc, v3xc, vlaplxc); + XC_Functional_Libxc::tau_xc( tls_funcs, arho, grho2a, lapl_val, atau, sxc, v1xc, v2xc, v3xc, vlaplxc, hybrid_alpha_in); if(!buf.vlapl_arr1.empty()) buf.vlapl_arr1[ir] = vlaplxc; } else @@ -189,7 +191,8 @@ void gradcorr_xc_kernel(const GradCorrParams& params, XC_Functional_Libxc::tau_xc_spin( tls_funcs, buf.rhotmp1[ir], buf.rhotmp2[ir], buf.gdr1[ir], buf.gdr2[ir], - laplup_val, lapldw_val, atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, vlaplxcup, vlaplxcdw); + laplup_val, lapldw_val, atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, vlaplxcup, vlaplxcdw, + hybrid_alpha_in); if(!buf.vlapl_arr1.empty()) buf.vlapl_arr1[ir] = vlaplxcup; if(!buf.vlapl_arr2.empty()) buf.vlapl_arr2[ir] = vlaplxcdw; } diff --git a/source/source_hamilt/module_xc/xc_pot.cpp b/source/source_hamilt/module_xc/xc_pot.cpp index 82ed6f615f1..f8b54bde2dd 100644 --- a/source/source_hamilt/module_xc/xc_pot.cpp +++ b/source/source_hamilt/module_xc/xc_pot.cpp @@ -134,79 +134,58 @@ std::tuple XC_Functional::v_xc( else if(nspin == 4) { #ifdef _OPENMP -#pragma omp parallel reduction(+:etxc) reduction(+:vtxc) - { -#endif - // Initialize the libxc functional once per thread instead of once - // per grid point; the nspin==4 LDA path uses the spin-polarized - // wrapper. -#ifdef __LIBXC - std::vector tls_funcs; - if (use_libxc) - { - tls_funcs = XC_Functional_Libxc::init_func( - XC_Functional::get_func_id(), XC_POLARIZED, hybrid_alpha, hse_omega); - } -#endif -#ifdef _OPENMP -#pragma omp for +#pragma omp parallel for reduction(+:etxc) reduction(+:vtxc) #endif - for(int ir = 0;irrho[1][ir],2) + pow(chr->rho[2][ir],2) + pow(chr->rho[3][ir],2) ); + double rhox = chr->rho[0][ir] + chr->rho_core[ir]; + double arhox = std::abs( rhox ); + + if ( arhox > vanishing_charge ) { - double amag = sqrt( pow(chr->rho[1][ir],2) + pow(chr->rho[2][ir],2) + pow(chr->rho[3][ir],2) ); - double rhox = chr->rho[0][ir] + chr->rho_core[ir]; - double arhox = std::abs( rhox ); + double zeta = amag / arhox; + double exc = 0.0; + double vxc[2]; - if ( arhox > vanishing_charge ) + if ( std::abs( zeta ) > 1.0 ) { - double zeta = amag / arhox; - double exc = 0.0; - double vxc[2]; - - if ( std::abs( zeta ) > 1.0 ) - { - zeta = (zeta > 0.0) ? 1.0 : (-1.0); - } + zeta = (zeta > 0.0) ? 1.0 : (-1.0); + } - if(use_libxc) - { + if(use_libxc) + { #ifdef __LIBXC - double rhoup = arhox * (1.0+zeta) / 2.0; - double rhodw = arhox * (1.0-zeta) / 2.0; - XC_Functional_Libxc::xc_spin_libxc(tls_funcs, rhoup, rhodw, exc, vxc[0], vxc[1]); + double rhoup = arhox * (1.0+zeta) / 2.0; + double rhodw = arhox * (1.0-zeta) / 2.0; + XC_Functional_Libxc::xc_spin_libxc(XC_Functional::get_func_id(), rhoup, rhodw, exc, vxc[0], vxc[1], hybrid_alpha, hse_omega); #else - ModuleBase::WARNING_QUIT("v_xc", "compile with LIBXC"); + ModuleBase::WARNING_QUIT("v_xc", "compile with LIBXC"); #endif - } - else - { - double rhoup = arhox * (1.0+zeta) / 2.0; - double rhodw = arhox * (1.0-zeta) / 2.0; - XC_Functional::xc_spin(arhox, zeta, exc, vxc[0], vxc[1]); - } + } + else + { + double rhoup = arhox * (1.0+zeta) / 2.0; + double rhodw = arhox * (1.0-zeta) / 2.0; + XC_Functional::xc_spin(arhox, zeta, exc, vxc[0], vxc[1]); + } - etxc += e2 * exc * rhox; + etxc += e2 * exc * rhox; - v(0, ir) = e2*( 0.5 * ( vxc[0] + vxc[1]) ); - vtxc += v(0,ir) * chr->rho[0][ir]; + v(0, ir) = e2*( 0.5 * ( vxc[0] + vxc[1]) ); + vtxc += v(0,ir) * chr->rho[0][ir]; - double vs = 0.5 * ( vxc[0] - vxc[1] ); - if ( amag > vanishing_charge ) + double vs = 0.5 * ( vxc[0] - vxc[1] ); + if ( amag > vanishing_charge ) + { + for(int ipol = 1;ipol< 4;ipol++) { - for(int ipol = 1;ipol< 4;ipol++) - { - v(ipol, ir) = e2 * vs * chr->rho[ipol][ir] / amag; - vtxc += v(ipol,ir) * chr->rho[ipol][ir]; - } + v(ipol, ir) = e2 * vs * chr->rho[ipol][ir] / amag; + vtxc += v(ipol,ir) * chr->rho[ipol][ir]; } } } -#ifdef __LIBXC - XC_Functional_Libxc::finish_func(tls_funcs); -#endif -#ifdef _OPENMP } -#endif } // energy terms, local-density contributions