diff --git a/source/source_estate/module_pot/pot_xc.cpp b/source/source_estate/module_pot/pot_xc.cpp index 4888380f23c..c456b0b2065 100644 --- a/source/source_estate/module_pot/pot_xc.cpp +++ b/source/source_estate/module_pot/pot_xc.cpp @@ -2,6 +2,7 @@ #include "source_base/timer.h" #include "source_hamilt/module_xc/xc_functional.h" +#include "source_io/module_parameter/parameter.h" #ifdef USE_LIBXC #include "source_hamilt/module_xc/libxc_abacus.h" @@ -24,7 +25,8 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module { #ifdef USE_LIBXC const std::tuple etxc_vtxc_v - = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), nrxx_current, ucell->omega, ucell->tpiba, chg); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), nrxx_current, ucell->omega, ucell->tpiba, chg, + PARAM.inp.nspin); *(this->etxc_) = std::get<0>(etxc_vtxc_v); *(this->vtxc_) = std::get<1>(etxc_vtxc_v); v_eff += std::get<2>(etxc_vtxc_v); @@ -36,7 +38,10 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module else { const std::tuple etxc_vtxc_v - = XC_Functional::v_xc(nrxx_current, chg, ucell); + = XC_Functional::v_xc(nrxx_current, chg, ucell, + PARAM.inp.nspin, + PARAM.globalv.domag, + PARAM.globalv.domag_z); *(this->etxc_) = std::get<0>(etxc_vtxc_v); *(this->vtxc_) = std::get<1>(etxc_vtxc_v); v_eff += std::get<2>(etxc_vtxc_v); diff --git a/source/source_estate/module_pot/pot_xc_fdm.cpp b/source/source_estate/module_pot/pot_xc_fdm.cpp index d0a49cc1012..aa950669d5b 100644 --- a/source/source_estate/module_pot/pot_xc_fdm.cpp +++ b/source/source_estate/module_pot/pot_xc_fdm.cpp @@ -5,6 +5,7 @@ #include "pot_xc_fdm.h" #include "source_hamilt/module_xc/xc_functional.h" +#include "source_io/module_parameter/parameter.h" namespace elecstate { @@ -20,7 +21,10 @@ PotXC_FDM::PotXC_FDM( this->fixed_mode = false; const std::tuple etxc_vtxc_v_0 - = XC_Functional::v_xc(this->chg_0->nrxx, this->chg_0, ucell); + = XC_Functional::v_xc(this->chg_0->nrxx, this->chg_0, ucell, + PARAM.inp.nspin, + PARAM.globalv.domag, + PARAM.globalv.domag_z); this->v_xc_0 = std::get<2>(etxc_vtxc_v_0); } @@ -47,7 +51,10 @@ void PotXC_FDM::cal_v_eff( } const std::tuple etxc_vtxc_v_01 - = XC_Functional::v_xc(chg_01.nrxx, &chg_01, ucell); + = XC_Functional::v_xc(chg_01.nrxx, &chg_01, ucell, + PARAM.inp.nspin, + PARAM.globalv.domag, + PARAM.globalv.domag_z); const ModuleBase::matrix &v_xc_01 = std::get<2>(etxc_vtxc_v_01); v_eff += v_xc_01 - this->v_xc_0; diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index 652293cd22f..4e63737cf44 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/source/source_hamilt/module_xc/libxc_abacus.h @@ -60,7 +60,10 @@ namespace XC_Functional_Libxc const double &omega, // volume of cell const double tpiba, const Charge* const chr, // charge density - const std::map* scaling_factor = nullptr); // added by jghan, 2024-10-10 + const int nspin, + const bool domag, + const bool domag_z, + const std::map* scaling_factor); // for mGGA functional extern std::tuple v_xc_meta( @@ -68,7 +71,8 @@ namespace XC_Functional_Libxc const int &nrxx, // number of real-space grid const double &omega, // volume of cell const double tpiba, - const Charge* const chr); + const Charge* const chr, + const int nspin); //------------------- diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index 9ff7190764f..1380346142b 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -19,13 +19,16 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / const double &omega, // volume of cell const double tpiba, const Charge* const chr, + const int nspin_in, + const bool domag, + const bool domag_z, const std::map* scaling_factor) { ModuleBase::TITLE("XC_Functional_Libxc","v_xc_libxc"); ModuleBase::timer::start("XC_Functional_Libxc","v_xc_libxc"); const int nspin = - (PARAM.inp.nspin == 1 || ( PARAM.inp.nspin ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) + (nspin_in == 1 || ( nspin_in ==4 && !domag && !domag_z)) ? 1 : 2; //---------------------------------------------------------- @@ -56,7 +59,7 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / // converting rho std::vector rho; std::vector amag; - if(1==nspin || 2==PARAM.inp.nspin) + if(1==nspin || 2==nspin_in) { rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); } @@ -167,7 +170,7 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / v += std::get<1>(vtxc_v) * factor; } // end for( xc_func_type &func : funcs ) - if(4==PARAM.inp.nspin) + if(4==nspin_in) { v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v); } @@ -207,7 +210,8 @@ std::tuple XC_Functional_Li const int &nrxx, // number of real-space grid const double &omega, // volume of cell const double tpiba, - const Charge* const chr) + const Charge* const chr, + const int nspin) { ModuleBase::TITLE("XC_Functional_Libxc","v_xc_meta"); ModuleBase::timer::start("XC_Functional_Libxc","v_xc_meta"); @@ -217,8 +221,8 @@ std::tuple XC_Functional_Li //output of the subroutine double etxc = 0.0; double vtxc = 0.0; - ModuleBase::matrix v(PARAM.inp.nspin,nrxx); - ModuleBase::matrix vofk(PARAM.inp.nspin,nrxx); + ModuleBase::matrix v(nspin,nrxx); + ModuleBase::matrix vofk(nspin,nrxx); //---------------------------------------------------------- // xc_func_type is defined in Libxc package @@ -226,8 +230,6 @@ std::tuple XC_Functional_Li // use can check on website, for example: // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ //---------------------------------------------------------- - - const int nspin = PARAM.inp.nspin; std::vector funcs = XC_Functional_Libxc::init_func( /* func_id = */ func_id, /* xc_polarized = */ (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED); diff --git a/source/source_hamilt/module_xc/test/test_xc3.cpp b/source/source_hamilt/module_xc/test/test_xc3.cpp index 96a67e785cd..fd17aaa8d71 100644 --- a/source/source_hamilt/module_xc/test/test_xc3.cpp +++ b/source/source_hamilt/module_xc/test/test_xc3.cpp @@ -1,8 +1,5 @@ #include "gtest/gtest.h" #include "xctest.h" -#define private public -#include "source_io/module_parameter/parameter.h" -#undef private #include "../xc_functional.h" #include "../exx_info.h" #include "xc3_mock.h" @@ -37,6 +34,14 @@ class XCTest_GRADCORR : public XCTest void SetUp() { + // Define variables for parameters + int nspin1 = 1; + int nspin2 = 2; + int nspin4 = 4; + bool domag = false; + bool domag_z = false; + bool domag_true = true; + ModulePW::PW_Basis rhopw; UnitCell ucell; Charge chr; @@ -84,17 +89,13 @@ class XCTest_GRADCORR : public XCTest XC_Functional::set_xc_type("PBE"); - PARAM.input.nspin = 1; - XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,false); - XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,true); + XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,false,nspin1,domag,domag_z); + XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,true,nspin1,domag,domag_z); - PARAM.input.nspin = 2; - XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,false); - XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,true); + XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,false,nspin2,domag,domag_z); + XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,true,nspin2,domag,domag_z); - PARAM.input.nspin = 4; - PARAM.sys.domag = true; - XC_Functional::gradcorr(et4,vt4,v4,&chr,&rhopw,&ucell,stress4,false); + XC_Functional::gradcorr(et4,vt4,v4,&chr,&rhopw,&ucell,stress4,false,nspin4,domag_true,domag_z); } }; diff --git a/source/source_hamilt/module_xc/test/test_xc5.cpp b/source/source_hamilt/module_xc/test/test_xc5.cpp index 469b3dcc207..555f9be8702 100644 --- a/source/source_hamilt/module_xc/test/test_xc5.cpp +++ b/source/source_hamilt/module_xc/test/test_xc5.cpp @@ -1,9 +1,6 @@ #include "../xc_functional.h" #include "../libxc_abacus.h" #include "gtest/gtest.h" -#define private public -#include "source_io/module_parameter/parameter.h" -#undef private #include "xctest.h" #include "../exx_info.h" #include "xc3_mock.h" @@ -32,6 +29,12 @@ class XCTest_VXC : public XCTest void SetUp() { + // Define variables for parameters + int nspin1 = 1; + int nspin2 = 2; + bool domag = false; + bool domag_z = false; + ModulePW::PW_Basis rhopw; UnitCell ucell; Charge chr; @@ -75,16 +78,14 @@ class XCTest_VXC : public XCTest XC_Functional::set_xc_type("PBE"); - PARAM.input.nspin = 1; std::tuple etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin1,domag,domag_z); et1 = std::get<0>(etxc_vtxc_v); vt1 = std::get<1>(etxc_vtxc_v); v1 = std::get<2>(etxc_vtxc_v); - PARAM.input.nspin = 2; etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin2,domag,domag_z); et2 = std::get<0>(etxc_vtxc_v); vt2 = std::get<1>(etxc_vtxc_v); v2 = std::get<2>(etxc_vtxc_v); @@ -130,6 +131,12 @@ class XCTest_VXC_Libxc : public XCTest void SetUp() { + // Define variables for parameters + int nspin1 = 1; + int nspin2 = 2; + bool domag = false; + bool domag_z = false; + ModulePW::PW_Basis rhopw; UnitCell ucell; Charge chr; @@ -173,16 +180,14 @@ class XCTest_VXC_Libxc : public XCTest XC_Functional::set_xc_type("GGA_X_PBE+GGA_C_PBE"); - PARAM.input.nspin = 1; std::tuple etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin1,domag,domag_z); et1 = std::get<0>(etxc_vtxc_v); vt1 = std::get<1>(etxc_vtxc_v); v1 = std::get<2>(etxc_vtxc_v); - PARAM.input.nspin = 2; etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin2,domag,domag_z); et2 = std::get<0>(etxc_vtxc_v); vt2 = std::get<1>(etxc_vtxc_v); v2 = std::get<2>(etxc_vtxc_v); @@ -228,6 +233,10 @@ class XCTest_VXC_meta : public XCTest void SetUp() { + // Define variables for parameters + int nspin1 = 1; + int nspin2 = 2; + ModulePW::PW_Basis rhopw; UnitCell ucell; Charge chr; @@ -281,17 +290,15 @@ class XCTest_VXC_meta : public XCTest XC_Functional::set_xc_type("SCAN"); - PARAM.input.nspin = 1; std::tuple etxc_vtxc_v - = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr,nspin1); et1 = std::get<0>(etxc_vtxc_v); vt1 = std::get<1>(etxc_vtxc_v); v1 = std::get<2>(etxc_vtxc_v); vtau1 = std::get<3>(etxc_vtxc_v); - PARAM.input.nspin = 2; etxc_vtxc_v - = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr,nspin2); et2 = std::get<0>(etxc_vtxc_v); vt2 = std::get<1>(etxc_vtxc_v); v2 = std::get<2>(etxc_vtxc_v); diff --git a/source/source_hamilt/module_xc/test/xctest.h b/source/source_hamilt/module_xc/test/xctest.h index 07ae4512c3f..608edbe72b4 100644 --- a/source/source_hamilt/module_xc/test/xctest.h +++ b/source/source_hamilt/module_xc/test/xctest.h @@ -1,17 +1,9 @@ #ifndef XCTEST_H #define XCTEST_H #include "gtest/gtest.h" -#define private public -#include "source_io/module_parameter/parameter.h" -#undef private class XCTest: public testing::Test { public: - XCTest() - { - PARAM.input.basis_type = ""; - PARAM.input.cal_force = 0; - PARAM.input.cal_stress = 0; - } + XCTest() {} }; #endif \ No newline at end of file diff --git a/source/source_hamilt/module_xc/xc_functional.cpp b/source/source_hamilt/module_xc/xc_functional.cpp index a2042b69333..4290b4ff080 100644 --- a/source/source_hamilt/module_xc/xc_functional.cpp +++ b/source/source_hamilt/module_xc/xc_functional.cpp @@ -301,21 +301,18 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) std::cerr << "\n OPTX untested please test,"; } - // if((func_type == 4 || func_type == 5) && PARAM.inp.basis_type == "pw") + // if((func_type == 4 || func_type == 5) && basis_type == "pw") // { // ModuleBase::WARNING_QUIT("set_xc_type","hybrid functional not realized for planewave yet"); // } - if((func_type == 3 || func_type == 5) && PARAM.inp.nspin==4) - { - ModuleBase::WARNING_QUIT("set_xc_type","meta-GGA has not been implemented for nspin = 4 yet"); - } -#ifndef __EXX - if((func_type == 4 || func_type == 5) && PARAM.inp.basis_type == "lcao") - { - ModuleBase::WARNING_QUIT("set_xc_type","compile with libri to use hybrid functional in lcao basis"); - } -#endif + // Hybrid functional is now supported for both PW and LCAO basis + // #ifndef __EXX + // if((func_type == 4 || func_type == 5) && basis_type == "lcao") + // { + // ModuleBase::WARNING_QUIT("set_xc_type","compile with libri to use hybrid functional in lcao basis"); + // } + // #endif #ifndef USE_LIBXC if(xc_func == "SCAN" || xc_func == "HSE" || xc_func == "SCAN0" diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index 37695900a9b..aaba3e945b7 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -47,7 +47,10 @@ class XC_Functional static std::tuple v_xc( const int &nrxx, // number of real-space grid const Charge* const chr, - const UnitCell *ucell); // charge density + const UnitCell *ucell, // charge density + const int nspin, + const bool domag, + const bool domag_z); //------------------- // xc_functional.cpp @@ -206,7 +209,11 @@ class XC_Functional ModulePW::PW_Basis* rhopw, const UnitCell* ucell, std::vector& stress_gga, - const bool is_stress = false); + const bool is_stress, + const int nspin, + const bool domag, + const bool domag_z); + template ::type> diff --git a/source/source_hamilt/module_xc/xc_grad.cpp b/source/source_hamilt/module_xc/xc_grad.cpp index 0129662543c..8104c03a61c 100644 --- a/source/source_hamilt/module_xc/xc_grad.cpp +++ b/source/source_hamilt/module_xc/xc_grad.cpp @@ -32,10 +32,18 @@ void XC_Functional::gradcorr( ModulePW::PW_Basis* rhopw, const UnitCell *ucell, std::vector &stress_gga, - const bool is_stress) + const bool is_stress, + const int nspin, + const bool domag, + const bool domag_z) { ModuleBase::TITLE("XC_Functional","gradcorr"); + if((func_type == 3 || func_type == 5) && nspin==4) + { + ModuleBase::WARNING_QUIT("gradcorr","meta-GGA has not been implemented for nspin = 4 yet"); + } + if(func_type == 0 || func_type == 1) { return; @@ -47,12 +55,12 @@ void XC_Functional::gradcorr( igcc_is_lyp = true; } - int nspin0 = PARAM.inp.nspin; - if(PARAM.inp.nspin==4) + int nspin0 = nspin; + if(nspin==4) { nspin0 =1; } - if(PARAM.inp.nspin==4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) + if(nspin==4&&(domag||domag_z)) { nspin0 = 2; } @@ -71,7 +79,7 @@ void XC_Functional::gradcorr( // doing FFT to get rho in G space: rhog1 rhopw->real2recip(chr->rho[0], chr->rhog[0]); - if(PARAM.inp.nspin==2) + if(nspin==2) { rhopw->real2recip(chr->rho[1], chr->rhog[1]); } @@ -120,7 +128,7 @@ void XC_Functional::gradcorr( // for spin polarized case; // calculate the gradient of (rho_core+rho) in reciprocal space. - if(PARAM.inp.nspin==2) + if(nspin==2) { rhotmp2 = new double[rhopw->nrxx]; rhogsum2 = new std::complex[rhopw->npw]; @@ -148,7 +156,7 @@ void XC_Functional::gradcorr( XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); } - if(PARAM.inp.nspin == 4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) + if(nspin == 4&&(domag||domag_z)) { rhotmp2 = new double[rhopw->nrxx]; rhogsum2 = new std::complex[rhopw->npw]; @@ -172,15 +180,15 @@ void XC_Functional::gradcorr( } if(!is_stress) { - vsave = new double* [PARAM.inp.nspin]; - for(int is = 0;isnrxx]; } #ifdef _OPENMP #pragma omp parallel for collapse(2) schedule(static, 1024) #endif - for(int is = 0;isnrxx;ir++) { @@ -445,7 +453,7 @@ void XC_Functional::gradcorr( else { double zeta = ( rhotmp1[ir] - rhotmp2[ir] ) / rh; - if(PARAM.inp.nspin==4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) + if(nspin==4&&(domag||domag_z)) { zeta = fabs(zeta) * neg[ir]; } @@ -606,12 +614,12 @@ void XC_Functional::gradcorr( vtxc += vtxcgc; etxc += etxcgc; - if(PARAM.inp.nspin == 4 && (PARAM.globalv.domag||PARAM.globalv.domag_z)) + if(nspin == 4 && (domag||domag_z)) { #ifdef _OPENMP #pragma omp parallel for collapse(2) schedule(static, 1024) #endif - for(int is=0;isnrxx;ir++) { @@ -648,7 +656,7 @@ void XC_Functional::gradcorr( delete[] h1; } - if(PARAM.inp.nspin==2) + if(nspin==2) { delete[] rhotmp2; delete[] rhogsum2; @@ -658,7 +666,7 @@ void XC_Functional::gradcorr( delete[] h2; } } - if(PARAM.inp.nspin == 4 && (PARAM.globalv.domag||PARAM.globalv.domag_z)) + if(nspin == 4 && (domag||domag_z)) { delete[] neg; if(!is_stress) @@ -668,7 +676,7 @@ void XC_Functional::gradcorr( delete[] vgg[i]; } delete[] vgg; - for(int i=0; i XC_Functional::v_xc( const int& nrxx, const Charge* const chr, - const UnitCell* ucell) + const UnitCell* ucell, + const int nspin, + const bool domag, + const bool domag_z) { ModuleBase::TITLE("XC_Functional", "v_xc"); @@ -29,6 +32,9 @@ std::tuple XC_Functional::v_xc( ucell->omega, ucell->tpiba, chr, + nspin, + domag, + domag_z, &(scaling_factor_xc)); #else ModuleBase::WARNING_QUIT("v_xc", "compile with LIBXC"); @@ -40,14 +46,14 @@ std::tuple XC_Functional::v_xc( //Exchange-Correlation potential Vxc(r) from n(r) double etxc = 0.0; double vtxc = 0.0; - ModuleBase::matrix v(PARAM.inp.nspin, nrxx); + ModuleBase::matrix v(nspin, nrxx); // the square of the e charge // in Rydeberg unit, so * 2.0. double e2 = 2.0; double vanishing_charge = 1.0e-10; - if (PARAM.inp.nspin == 1 || ( PARAM.inp.nspin ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) + if (nspin == 1 || ( nspin ==4 && !domag && !domag_z)) { // spin-unpolarized case #ifdef _OPENMP @@ -71,7 +77,7 @@ std::tuple XC_Functional::v_xc( } // endif } //enddo } - else if(PARAM.inp.nspin ==2) + else if(nspin ==2) { // spin-polarized case #ifdef _OPENMP @@ -97,7 +103,7 @@ std::tuple XC_Functional::v_xc( double vxc[2]; XC_Functional::xc_spin(arhox, zeta, exc, vxc[0], vxc[1]); - for (int is = 0;is < PARAM.inp.nspin;is++) + for (int is = 0;is < nspin;is++) { v(is, ir) = e2 * vxc[is]; } @@ -107,7 +113,7 @@ std::tuple XC_Functional::v_xc( } } } - else if(PARAM.inp.nspin == 4) + else if(nspin == 4) { #ifdef _OPENMP #pragma omp parallel for reduction(+:etxc) reduction(+:vtxc) @@ -171,7 +177,7 @@ std::tuple XC_Functional::v_xc( // the dummy variable dum contains gradient correction to stress // which is not used here std::vector dum; - gradcorr(etxc, vtxc, v, chr, chr->rhopw, ucell, dum); + gradcorr(etxc, vtxc, v, chr, chr->rhopw, ucell, dum, false, nspin, domag, domag_z); // parallel code : collect vtxc,etxc // mohan add 2008-06-01 diff --git a/source/source_pw/module_pwdft/forces_cc.cpp b/source/source_pw/module_pwdft/forces_cc.cpp index 45f7ea2fb16..41322e62e55 100644 --- a/source/source_pw/module_pwdft/forces_cc.cpp +++ b/source/source_pw/module_pwdft/forces_cc.cpp @@ -60,7 +60,8 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, { #ifdef USE_LIBXC const auto etxc_vtxc_v - = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rho_basis->nrxx, ucell_in.omega, ucell_in.tpiba, chr); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rho_basis->nrxx, ucell_in.omega, ucell_in.tpiba, chr, + PARAM.inp.nspin); // etxc = std::get<0>(etxc_vtxc_v); // vtxc = std::get<1>(etxc_vtxc_v); @@ -72,7 +73,10 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, else { elecstate::cal_ux(ucell_in); - const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell_in); + const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell_in, + PARAM.inp.nspin, + PARAM.globalv.domag, + PARAM.globalv.domag_z); // etxc = std::get<0>(etxc_vtxc_v); // vtxc = std::get<1>(etxc_vtxc_v); diff --git a/source/source_pw/module_pwdft/stress_cc.cpp b/source/source_pw/module_pwdft/stress_cc.cpp index bc986cdda76..3b0f291e889 100644 --- a/source/source_pw/module_pwdft/stress_cc.cpp +++ b/source/source_pw/module_pwdft/stress_cc.cpp @@ -55,7 +55,8 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, { #ifdef USE_LIBXC const auto etxc_vtxc_v - = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rho_basis->nrxx, ucell.omega, ucell.tpiba, chr); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rho_basis->nrxx, ucell.omega, ucell.tpiba, chr, + PARAM.inp.nspin); // etxc = std::get<0>(etxc_vtxc_v); // vtxc = std::get<1>(etxc_vtxc_v); @@ -67,7 +68,10 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, else { elecstate::cal_ux(ucell); - const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell); + const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell, + PARAM.inp.nspin, + PARAM.globalv.domag, + PARAM.globalv.domag_z); // etxc = std::get<0>(etxc_vtxc_v); // may delete? // vtxc = std::get<1>(etxc_vtxc_v); // may delete? vxc = std::get<2>(etxc_vtxc_v); diff --git a/source/source_pw/module_pwdft/stress_gga.cpp b/source/source_pw/module_pwdft/stress_gga.cpp index 0c23f58ef05..198e4cc843a 100644 --- a/source/source_pw/module_pwdft/stress_gga.cpp +++ b/source/source_pw/module_pwdft/stress_gga.cpp @@ -2,6 +2,7 @@ #include "source_base/parallel_reduce.h" #include "source_hamilt/module_xc/xc_functional.h" #include "source_base/timer.h" +#include "source_io/module_parameter/parameter.h" //calculate the GGA stress correction in PW and LCAO template @@ -25,10 +26,14 @@ void Stress_Func::stress_gga(const UnitCell& ucell, FPTYPE dum1=0.0; FPTYPE dum2=0.0; ModuleBase::matrix dum3; + const bool is_stress = true; // call gradcorr to evaluate gradient correction to stress // the first three terms are etxc, vtxc and v, which // is not used here, so dummy variables are used. - XC_Functional::gradcorr(dum1, dum2, dum3, chr, rho_basis, &ucell, stress_gga, 1); + XC_Functional::gradcorr( + dum1, dum2, dum3, chr, rho_basis, &ucell, + stress_gga, is_stress, + PARAM.inp.nspin, PARAM.globalv.domag, PARAM.globalv.domag_z); for(int l = 0;l< 3;l++) {