diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index b618bbbe90..d724837db0 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/source/source_hamilt/module_xc/libxc_abacus.h @@ -73,7 +73,6 @@ namespace XC_Functional_Libxc extern void finish_func(std::vector &funcs); - //------------------- // libxc_pot.cpp //------------------- @@ -282,6 +281,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 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 +308,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 +342,21 @@ namespace XC_Functional_Libxc const double &hybrid_alpha, 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, + const double &grho, + const double &lapl_rho, + const double &atau, + double &sxc, + double &v1xc, + double &v2xc, + double &v3xc, + double &vlaplxc, + const double hybrid_alpha); + extern void tau_xc_spin( const std::vector &func_id, double rhoup, @@ -341,6 +380,30 @@ namespace XC_Functional_Libxc const double &hybrid_alpha, 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, + 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); + } // 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 bd86a46657..29a259f7df 100644 --- a/source/source_hamilt/module_xc/libxc_gga_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_gga_wrap.cpp @@ -8,34 +8,35 @@ #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& 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; 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; } - 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; @@ -45,6 +46,34 @@ void XC_Functional_Libxc::gcxc_libxc( 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) +{ + // 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, + /* 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 +81,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 +91,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 +102,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) { @@ -113,6 +134,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_mgga_wrap.cpp b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp index 3fd7f34b53..e6031ba5fa 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, @@ -24,19 +24,13 @@ void XC_Functional_Libxc::tau_xc( double& v2xc, double& v3xc, double& vlaplxc, - const double& hybrid_alpha, - const double& hse_omega) + const double hybrid_alpha) { 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,9 +38,10 @@ 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(&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) { @@ -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, hybrid_alpha); + 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, @@ -89,8 +110,7 @@ void XC_Functional_Libxc::tau_xc_spin( double& v3xcdw, double& vlaplxcup, double& vlaplxcdw, - const double& hybrid_alpha, - const double& hse_omega) + const double hybrid_alpha) { sxc = 0.0; v1xcup = 0.0; @@ -108,13 +128,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,8 +153,9 @@ 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(&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) @@ -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, hybrid_alpha); XC_Functional_Libxc::finish_func(funcs); } diff --git a/source/source_hamilt/module_xc/test/test_xc3.cpp b/source/source_hamilt/module_xc/test/test_xc3.cpp index cd7d0e52e3..0171e7635b 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 d4fac72332..345e18d6e4 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 @@ -159,3 +160,85 @@ TEST(XC_ScanL_Reference, tau_xc_wrapper_libxc) EXPECT_NE(v1xc, 0.0); EXPECT_NE(vlaplxc, 0.0); } + +#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 +{ + 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; + + 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(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, + scaled_funcs[0], scaled_funcs[1], scaled_funcs[2], scaled_funcs[3], scaled_funcs[4], + hybrid_alpha); + XC_Functional_Libxc::finish_func(funcs); + + 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 78d80e4cec..09bda321d6 100644 --- a/source/source_hamilt/module_xc/xc_grad_kernel.cpp +++ b/source/source_hamilt/module_xc/xc_grad_kernel.cpp @@ -58,6 +58,22 @@ 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. + // The nspin0==1 non-stress branch uses XC_Functional::gcxc, not libxc. +#ifdef __LIBXC + std::vector tls_funcs; + 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( + 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 +118,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, hybrid_alpha_in); 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 +189,19 @@ 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, + hybrid_alpha_in); 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 +335,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) {