diff --git a/source/source_cell/sep_cell.cpp b/source/source_cell/sep_cell.cpp index e7ec5f1bafa..35d0c8f2ab8 100644 --- a/source/source_cell/sep_cell.cpp +++ b/source/source_cell/sep_cell.cpp @@ -9,11 +9,6 @@ #include #include -// namespace GlobalC -// { -// Sep_Cell sep_cell; -// } - Sep_Cell::Sep_Cell() noexcept : ntype(0), omega(0.0), tpiba2(0.0) { } diff --git a/source/source_cell/sep_cell.h b/source/source_cell/sep_cell.h index 03cddca2ea0..253ae4905e8 100644 --- a/source/source_cell/sep_cell.h +++ b/source/source_cell/sep_cell.h @@ -66,9 +66,4 @@ class Sep_Cell double tpiba2; // tpiba ^ 2 }; -// namespace GlobalC -// { -// extern Sep_Cell sep_cell; -// } - #endif // SEP_CEll diff --git a/source/source_cell/unitcell.cpp b/source/source_cell/unitcell.cpp index f92b37b1f8d..0af6fbf6fb8 100644 --- a/source/source_cell/unitcell.cpp +++ b/source/source_cell/unitcell.cpp @@ -253,9 +253,6 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) // readl sep potential, currently using the pseudopotential folder (pseudo_dir in INPUT) //========================== if (PARAM.inp.dfthalf_type > 0) { - // GlobalC::sep_cell.init(this->ntype); - // ok3 = GlobalC::sep_cell.read_sep_potentials(ifa, PARAM.inp.pseudo_dir, GlobalV::ofs_warning, this->atom_label); - sep_cell.init(this->ntype); ok3 = sep_cell.read_sep_potentials(ifa, PARAM.inp.pseudo_dir, GlobalV::ofs_warning, this->atom_label); } @@ -285,7 +282,6 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) #ifdef __MPI unitcell::bcast_unitcell(*this); - // GlobalC::sep_cell.bcast_sep_cell(); sep_cell.bcast_sep_cell(); #endif @@ -350,7 +346,6 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) //=================================== this->set_iat2itia(); - // GlobalC::sep_cell.set_omega(this->omega, this->tpiba2); sep_cell.set_omega(this->omega, this->tpiba2); return; diff --git a/source/source_esolver/esolver_ks_lcaopw.cpp b/source/source_esolver/esolver_ks_lcaopw.cpp index b77a0345cf4..4bbc9878aee 100644 --- a/source/source_esolver/esolver_ks_lcaopw.cpp +++ b/source/source_esolver/esolver_ks_lcaopw.cpp @@ -2,6 +2,8 @@ #include "source_pw/module_pwdft/elecond.h" #include "source_io/module_parameter/input_conv.h" +#include "source_io/module_output/output_log.h" +#include "source_hamilt/module_xc/exx_info.h" #include @@ -141,9 +143,11 @@ namespace ModuleESolver // add exx #ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) + bool cal_exx = GlobalC::exx_info.info_global.cal_exx; + double hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; + if (cal_exx) { - this->pelec->set_exx(this->exx_lip->get_exx_energy()); // Peize Lin add 2019-03-09 + this->pelec->set_exx(this->exx_lip->get_exx_energy(), cal_exx, hybrid_alpha); // Peize Lin add 2019-03-09 } #endif @@ -227,6 +231,8 @@ namespace ModuleESolver #ifdef __LCAO if (PARAM.inp.out_mat_xc) { + bool cal_exx = GlobalC::exx_info.info_global.cal_exx; + double hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; ModuleIO::write_Vxc(PARAM.inp.nspin, PARAM.globalv.nlocal, GlobalV::DRANK, @@ -240,7 +246,9 @@ namespace ModuleESolver this->locpp.vloc, this->chr, this->kv, - this->pelec->wg + this->pelec->wg, + cal_exx, + hybrid_alpha #ifdef __EXX , *this->exx_lip diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index 79448f713ad..f08d2aa99f3 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -239,9 +239,11 @@ template void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) { // Related to EXX - if (GlobalC::exx_info.info_global.cal_exx && !exx_helper->get_op_first_iter()) + bool cal_exx = GlobalC::exx_info.info_global.cal_exx; + double hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; + if (cal_exx && !exx_helper->get_op_first_iter()) { - this->pelec->set_exx(exx_helper->cal_exx_energy(this->stp.template get_psi_t())); + this->pelec->set_exx(exx_helper->cal_exx_energy(this->stp.template get_psi_t()), cal_exx, hybrid_alpha); } // deband is calculated from "output" charge density diff --git a/source/source_esolver/esolver_of_interface.cpp b/source/source_esolver/esolver_of_interface.cpp index d428b795069..1c507e7a95f 100644 --- a/source/source_esolver/esolver_of_interface.cpp +++ b/source/source_esolver/esolver_of_interface.cpp @@ -207,19 +207,6 @@ void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& // while(true) // { // this->pelec->f_en.calculate_etot(this->pw_rho->nrxx, - // this->pw_rho->nxyz); temp_energy = - // this->pelec->f_en.etot; kinetic_energy = - // this->kinetic_energy(); pseudopot_energy = 0.; for (int - // is = 0; is < PARAM.inp.nspin; ++is) { - // pseudopot_energy += - // this->inner_product(GlobalC::pot.vltot, - // ptemp_rho_[is], this->pw_rho->nrxx, this->dV_); - // } - // Parallel_Reduce::reduce_all(pseudopot_energy); - // temp_energy += kinetic_energy + pseudopot_energy; - // this->opt_dcsrch_->dcSrch(temp_energy, dEdalpha, - // thetaAlpha, this->task_); numDC++; - // if (strncmp(this->task_, "FG", 2) == 0) // { // for (int is = 0; is < PARAM.inp.nspin; ++is) diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index a4191df38d8..228d5248368 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -115,8 +115,8 @@ class ElecState bool vnew_exist = false; void cal_converged(); void cal_energies(const int type); - void set_exx(const double& Eexx); - void set_exx(const std::complex& Eexx); + void set_exx(const double& Eexx, const bool cal_exx, const double hybrid_alpha); + void set_exx(const std::complex& Eexx, const bool cal_exx, const double hybrid_alpha); double get_hartree_energy(); double get_etot_efield(); diff --git a/source/source_estate/elecstate_exx.cpp b/source/source_estate/elecstate_exx.cpp index addc7a03da0..ce538ddebd4 100644 --- a/source/source_estate/elecstate_exx.cpp +++ b/source/source_estate/elecstate_exx.cpp @@ -1,20 +1,27 @@ #include "source_estate/elecstate.h" -#include "source_hamilt/module_xc/exx_info.h" // use GlobalC::exx_info +#include "source_base/tool_quit.h" +#include // use std::complex namespace elecstate { /// @brief calculation if converged /// @date Peize Lin add 2016-12-03 -void ElecState::set_exx(const double& Eexx) +void ElecState::set_exx(const double& Eexx, const bool cal_exx, const double hybrid_alpha) { ModuleBase::TITLE("energy", "set_exx"); - if (GlobalC::exx_info.info_global.cal_exx) + if (cal_exx) { - this->f_en.exx = GlobalC::exx_info.info_global.hybrid_alpha * Eexx; + this->f_en.exx = hybrid_alpha * Eexx; } return; } +void ElecState::set_exx(const std::complex& Eexx, const bool cal_exx, const double hybrid_alpha) +{ + ModuleBase::WARNING_QUIT("ElecState::set_exx", + "std::complex version is not implemented yet"); +} + } diff --git a/source/source_estate/module_pot/pot_sep.cpp b/source/source_estate/module_pot/pot_sep.cpp index 2d5e0ee0bb0..118675e1bef 100644 --- a/source/source_estate/module_pot/pot_sep.cpp +++ b/source/source_estate/module_pot/pot_sep.cpp @@ -10,10 +10,6 @@ void PotSep::cal_fixed_v(double* vl_pseudo) ModuleBase::TITLE("PotSep", "cal_fixed_v"); ModuleBase::timer::start("PotSep", "cal_fixed_v"); - // GlobalC::vsep_cell.generate_vsep_r(this->rho_basis_[0], this->sf_[0]); - - // const_cast(this->vsep_)->generate_vsep_r(this->rho_basis_[0], this->sf_[0]); - if (vsep_cell != nullptr) { for (int ir = 0; ir < this->rho_basis_->nrxx; ++ir) diff --git a/source/source_estate/module_pot/pot_xc.cpp b/source/source_estate/module_pot/pot_xc.cpp index c456b0b2065..056af3099e9 100644 --- a/source/source_estate/module_pot/pot_xc.cpp +++ b/source/source_estate/module_pot/pot_xc.cpp @@ -24,9 +24,15 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module if (XC_Functional::get_ked_flag()) { #ifdef USE_LIBXC + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); +#ifdef __EXX + const double hse_omega = XC_Functional::get_hse_omega(); +#else + const double hse_omega = 0.0; +#endif const std::tuple etxc_vtxc_v = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), nrxx_current, ucell->omega, ucell->tpiba, chg, - PARAM.inp.nspin); + PARAM.inp.nspin, hybrid_alpha, hse_omega); *(this->etxc_) = std::get<0>(etxc_vtxc_v); *(this->vtxc_) = std::get<1>(etxc_vtxc_v); v_eff += std::get<2>(etxc_vtxc_v); @@ -37,11 +43,19 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module } else { + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); +#ifdef __EXX + const double hse_omega = XC_Functional::get_hse_omega(); +#else + const double hse_omega = 0.0; +#endif const std::tuple etxc_vtxc_v = XC_Functional::v_xc(nrxx_current, chg, ucell, PARAM.inp.nspin, PARAM.globalv.domag, - PARAM.globalv.domag_z); + PARAM.globalv.domag_z, + hybrid_alpha, + hse_omega); *(this->etxc_) = std::get<0>(etxc_vtxc_v); *(this->vtxc_) = std::get<1>(etxc_vtxc_v); v_eff += std::get<2>(etxc_vtxc_v); diff --git a/source/source_estate/module_pot/pot_xc_fdm.cpp b/source/source_estate/module_pot/pot_xc_fdm.cpp index bec10b8fa19..b0103db5f8b 100644 --- a/source/source_estate/module_pot/pot_xc_fdm.cpp +++ b/source/source_estate/module_pot/pot_xc_fdm.cpp @@ -1,65 +1,82 @@ -//======================= -// AUTHOR : Peize Lin -// DATE : 2025-10-01 -//======================= - -#include "pot_xc_fdm.h" -#include "source_hamilt/module_xc/xc_functional.h" -#include "source_io/module_parameter/parameter.h" - -namespace elecstate -{ - -PotXC_FDM::PotXC_FDM( - const ModulePW::PW_Basis* rho_basis_in, - const Charge*const chg_0_in, - const UnitCell*const ucell) - : chg_0(chg_0_in) -{ - this->rho_basis_ = rho_basis_in; - this->dynamic_mode = true; - this->fixed_mode = false; - - const std::tuple etxc_vtxc_v_0 - = 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); -} - -void PotXC_FDM::cal_v_eff( - const Charge*const chg_1, - const UnitCell*const ucell, - ModuleBase::matrix& v_eff) -{ - ModuleBase::TITLE("PotXC_FDM", "cal_veff"); - ModuleBase::timer::start("PotXC_FDM", "cal_veff"); - - assert(this->chg_0->nrxx == chg_1->nrxx); - assert(this->chg_0->nspin == chg_1->nspin); - - Charge chg_01; - chg_01.set_rhopw(chg_1->rhopw); - chg_01.allocate(chg_1->nspin, chg_01.kin_density()); - - for(int ir=0; irrho[is][ir] + chg_1->rho[is][ir]; } - chg_01.rho_core[ir] = chg_0->rho_core[ir] + chg_1->rho_core[ir]; - } - - const std::tuple etxc_vtxc_v_01 - = 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; - - ModuleBase::timer::end("PotXC_FDM", "cal_veff"); -} - -} // namespace elecstate +//======================= +// AUTHOR : Peize Lin +// DATE : 2025-10-01 +//======================= + +#include "pot_xc_fdm.h" +#include "source_hamilt/module_xc/xc_functional.h" +#include "source_io/module_parameter/parameter.h" + +namespace elecstate +{ + +PotXC_FDM::PotXC_FDM( + const ModulePW::PW_Basis* rho_basis_in, + const Charge*const chg_0_in, + const UnitCell*const ucell) + : chg_0(chg_0_in) +{ + this->rho_basis_ = rho_basis_in; + this->dynamic_mode = true; + this->fixed_mode = false; + + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); +#ifdef __EXX + const double hse_omega = XC_Functional::get_hse_omega(); +#else + const double hse_omega = 0.0; +#endif + const std::tuple etxc_vtxc_v_0 + = XC_Functional::v_xc(this->chg_0->nrxx, this->chg_0, ucell, + PARAM.inp.nspin, + PARAM.globalv.domag, + PARAM.globalv.domag_z, + hybrid_alpha, + hse_omega); + this->v_xc_0 = std::get<2>(etxc_vtxc_v_0); +} + +void PotXC_FDM::cal_v_eff( + const Charge*const chg_1, + const UnitCell*const ucell, + ModuleBase::matrix& v_eff) +{ + ModuleBase::TITLE("PotXC_FDM", "cal_veff"); + ModuleBase::timer::start("PotXC_FDM", "cal_veff"); + + assert(this->chg_0->nrxx == chg_1->nrxx); + assert(this->chg_0->nspin == chg_1->nspin); + + Charge chg_01; + chg_01.set_rhopw(chg_1->rhopw); + chg_01.allocate(chg_1->nspin, chg_01.kin_density()); + + for(int ir=0; irrho[is][ir] + chg_1->rho[is][ir]; } + chg_01.rho_core[ir] = chg_0->rho_core[ir] + chg_1->rho_core[ir]; + } + + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); +#ifdef __EXX + const double hse_omega = XC_Functional::get_hse_omega(); +#else + const double hse_omega = 0.0; +#endif + const std::tuple etxc_vtxc_v_01 + = XC_Functional::v_xc(chg_01.nrxx, &chg_01, ucell, + PARAM.inp.nspin, + PARAM.globalv.domag, + PARAM.globalv.domag_z, + hybrid_alpha, + hse_omega); + const ModuleBase::matrix &v_xc_01 = std::get<2>(etxc_vtxc_v_01); + + v_eff += v_xc_01 - this->v_xc_0; + + ModuleBase::timer::end("PotXC_FDM", "cal_veff"); +} + +} // namespace elecstate + diff --git a/source/source_hamilt/module_surchem/sol_force.cpp b/source/source_hamilt/module_surchem/sol_force.cpp index 8c5e6b6096b..c998243531d 100644 --- a/source/source_hamilt/module_surchem/sol_force.cpp +++ b/source/source_hamilt/module_surchem/sol_force.cpp @@ -17,7 +17,6 @@ void surchem::force_cor_one(const UnitCell& cell, //ModuleBase::GlobalFunc::ZEROS(delta_phi_g, rho_basis->npw); rho_basis->real2recip(this->delta_phi, delta_phi_g); - // GlobalC::UFFT.ToReciSpace(this->delta_phi, delta_phi_g,rho_basis); // double Ael = 0; // double Ael1 = 0; // ModuleBase::GlobalFunc::ZEROS(vg, ngmc); diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index 4e63737cf44..0c912858bf4 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/source/source_hamilt/module_xc/libxc_abacus.h @@ -45,7 +45,9 @@ namespace XC_Functional_Libxc */ extern std::vector init_func( const std::vector &func_id, - const int xc_polarized); + const int xc_polarized, + const double hybrid_alpha, + const double hse_omega); extern void finish_func(std::vector &funcs); @@ -63,7 +65,9 @@ namespace XC_Functional_Libxc const int nspin, const bool domag, const bool domag_z, - const std::map* scaling_factor); + const std::map* scaling_factor, + const double hybrid_alpha, + const double hse_omega); // for mGGA functional extern std::tuple v_xc_meta( @@ -72,7 +76,9 @@ namespace XC_Functional_Libxc const double &omega, // volume of cell const double tpiba, const Charge* const chr, - const int nspin); + const int nspin, + const double hybrid_alpha, + const double hse_omega); //------------------- @@ -162,7 +168,9 @@ namespace XC_Functional_Libxc const double &rhodw, double &exc, double &vxcup, - double &vxcdw); + double &vxcdw, + const double hybrid_alpha, + const double hse_omega); //------------------- @@ -176,7 +184,9 @@ namespace XC_Functional_Libxc const double &grho, double &sxc, double &v1xc, - double &v2xc); + double &v2xc, + const double hybrid_alpha, + const double hse_omega); // the entire GGA functional, for nspin=2 case extern void gcxc_spin_libxc( @@ -190,7 +200,9 @@ namespace XC_Functional_Libxc double &v1xcdw, double &v2xcup, double &v2xcdw, - double &v2xcud); + double &v2xcud, + const double hybrid_alpha, + const double hse_omega); //------------------- @@ -207,7 +219,8 @@ namespace XC_Functional_Libxc double &v1xc, double &v2xc, double &v3xc, - const double &hybrid_alpha); + const double &hybrid_alpha, + const double &hse_omega); extern void tau_xc_spin( const std::vector &func_id, @@ -225,7 +238,8 @@ namespace XC_Functional_Libxc double &v2xcud, double &v3xcup, double &v3xcdw, - const double &hybrid_alpha); + const double &hybrid_alpha, + const double &hse_omega); } // 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 8077469166b..cd5aaeb858b 100644 --- a/source/source_hamilt/module_xc/libxc_gga_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_gga_wrap.cpp @@ -1,6 +1,9 @@ #ifdef USE_LIBXC #include "libxc_abacus.h" +#ifdef __EXX +#include "source_hamilt/module_xc/exx_info.h" +#endif #include #include @@ -11,7 +14,9 @@ void XC_Functional_Libxc::gcxc_libxc( const double& grho, double& sxc, double& v1xc, - double& v2xc) + double& v2xc, + const double hybrid_alpha, + const double hse_omega) { sxc = 0.0; v1xc = 0.0; @@ -26,7 +31,9 @@ void XC_Functional_Libxc::gcxc_libxc( std::vector funcs = XC_Functional_Libxc::init_func( /* func_id = */ func_id, - /* xc_polarized = */ XC_UNPOLARIZED); + /* xc_polarized = */ XC_UNPOLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); for (xc_func_type& func : funcs) { @@ -55,7 +62,9 @@ void XC_Functional_Libxc::gcxc_spin_libxc( double& v1xcdw, double& v2xcup, double& v2xcdw, - double& v2xcud) + double& v2xcud, + const double hybrid_alpha, + const double hse_omega) { sxc = 0.0; v1xcup = 0.0; @@ -68,7 +77,9 @@ void XC_Functional_Libxc::gcxc_spin_libxc( std::vector funcs = XC_Functional_Libxc::init_func( /* func_id = */ func_id, - /* xc_polarized = */ XC_POLARIZED); + /* xc_polarized = */ XC_POLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); for (xc_func_type& 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 57509c4e753..646d65edbf4 100644 --- a/source/source_hamilt/module_xc/libxc_lda_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_lda_wrap.cpp @@ -1,18 +1,25 @@ #ifdef USE_LIBXC #include "libxc_abacus.h" +#ifdef __EXX +#include "source_hamilt/module_xc/exx_info.h" +#endif void XC_Functional_Libxc::xc_spin_libxc( 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; std::vector funcs = XC_Functional_Libxc::init_func( /* func_id = */ func_id, - /* xc_polarized = */ XC_POLARIZED); + /* xc_polarized = */ XC_POLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); for(xc_func_type &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 038f965d155..2c7d700f43a 100644 --- a/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp @@ -22,7 +22,8 @@ void XC_Functional_Libxc::tau_xc( double& v1xc, double& v2xc, double& v3xc, - const double& hybrid_alpha) + const double& hybrid_alpha, + const double& hse_omega) { double s = 0.0; double v1 = 0.0; @@ -32,7 +33,9 @@ void XC_Functional_Libxc::tau_xc( double vlapl_rho = 0.0; std::vector funcs = XC_Functional_Libxc::init_func( /* func_id = */ func_id, - /* xc_polarized = */ XC_UNPOLARIZED); + /* xc_polarized = */ XC_UNPOLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); sxc = 0.0; v1xc = 0.0; @@ -78,7 +81,8 @@ void XC_Functional_Libxc::tau_xc_spin( double& v2xcud, double& v3xcup, double& v3xcdw, - const double& hybrid_alpha) + const double& hybrid_alpha, + const double& hse_omega) { sxc = 0.0; v1xcup = 0.0; @@ -95,7 +99,9 @@ void XC_Functional_Libxc::tau_xc_spin( std::vector funcs = XC_Functional_Libxc::init_func( /* func_id = */ func_id, - /* xc_polarized = */ XC_POLARIZED); + /* xc_polarized = */ XC_POLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); for (xc_func_type& func : funcs) { diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index 1380346142b..9d70a41274f 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -8,6 +8,9 @@ #include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_base/tool_title.h" +#ifdef __EXX +#include "source_hamilt/module_xc/exx_info.h" +#endif #include @@ -22,7 +25,9 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / const int nspin_in, const bool domag, const bool domag_z, - const std::map* scaling_factor) + const std::map* scaling_factor, + const double hybrid_alpha, + const double hse_omega) { ModuleBase::TITLE("XC_Functional_Libxc","v_xc_libxc"); ModuleBase::timer::start("XC_Functional_Libxc","v_xc_libxc"); @@ -40,7 +45,9 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / std::vector funcs = XC_Functional_Libxc::init_func( /* func_id = */ func_id, - /* xc_polarized = */ (1==nspin) ? XC_UNPOLARIZED : XC_POLARIZED); + /* xc_polarized = */ (1==nspin) ? XC_UNPOLARIZED : XC_POLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); const bool is_gga = [&funcs]() { @@ -211,7 +218,9 @@ std::tuple XC_Functional_Li const double &omega, // volume of cell const double tpiba, const Charge* const chr, - const int nspin) + const int nspin, + const double hybrid_alpha, + const double hse_omega) { ModuleBase::TITLE("XC_Functional_Libxc","v_xc_meta"); ModuleBase::timer::start("XC_Functional_Libxc","v_xc_meta"); @@ -232,7 +241,9 @@ std::tuple XC_Functional_Li //---------------------------------------------------------- std::vector funcs = XC_Functional_Libxc::init_func( /* func_id = */ func_id, - /* xc_polarized = */ (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED); + /* xc_polarized = */ (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED, + /* hybrid_alpha = */ hybrid_alpha, + /* hse_omega = */ hse_omega); const std::vector rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); const std::vector>> gdr diff --git a/source/source_hamilt/module_xc/libxc_setup.cpp b/source/source_hamilt/module_xc/libxc_setup.cpp index 892aab5ba71..84972a80635 100644 --- a/source/source_hamilt/module_xc/libxc_setup.cpp +++ b/source/source_hamilt/module_xc/libxc_setup.cpp @@ -184,7 +184,9 @@ XC_Functional_Libxc::set_xc_type_libxc(const std::string& xc_func_in) return std::make_pair(func_type, func_id); } -const std::vector in_built_xc_func_ext_params(const int id) +const std::vector in_built_xc_func_ext_params(const int id, + const double hybrid_alpha, + const double hse_omega) { switch(id) { @@ -198,27 +200,23 @@ const std::vector in_built_xc_func_ext_params(const int id) #ifdef __EXX // hybrid functionals case XC_HYB_GGA_XC_PBEH: - return {GlobalC::exx_info.info_global.hybrid_alpha, - GlobalC::exx_info.info_global.hse_omega, - GlobalC::exx_info.info_global.hse_omega}; + return {hybrid_alpha, hse_omega, hse_omega}; case XC_HYB_GGA_XC_HSE06: - return {GlobalC::exx_info.info_global.hybrid_alpha, - GlobalC::exx_info.info_global.hse_omega, - GlobalC::exx_info.info_global.hse_omega}; + return {hybrid_alpha, hse_omega, hse_omega}; // short-range of B88_X case XC_GGA_X_ITYH: - return {GlobalC::exx_info.info_global.hse_omega}; + return {hse_omega}; // short-range of LYP_C case XC_GGA_C_LYPR: return {0.04918, 0.132, 0.2533, 0.349, - 0.35/2.29, 2.0/2.29, GlobalC::exx_info.info_global.hse_omega}; + 0.35/2.29, 2.0/2.29, hse_omega}; // Long-range corrected functionals: case XC_HYB_GGA_XC_LC_PBEOP: // LC version of PBE { // This is a range-separated hybrid functional with range-separation constant 0.330, // and 0.0% short-range and 100.0% long-range exact exchange, // using the error function kernel. - return { GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.33 + return { hse_omega }; //Range separation constant: 0.33 } case XC_HYB_GGA_XC_LC_WPBE: // Long-range corrected PBE (LC-wPBE) by Vydrov and Scuseria { @@ -227,7 +225,7 @@ const std::vector in_built_xc_func_ext_params(const int id) // using the error function kernel. return { std::stod(PARAM.inp.exx_fock_alpha[0]), //Fraction of Hartree-Fock exchange: 1.0 std::stod(PARAM.inp.exx_erfc_alpha[0]), //Fraction of short-range exact exchange: -1.0 - GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.4 + hse_omega }; //Range separation constant: 0.4 } case XC_HYB_GGA_XC_LRC_WPBE: // Long-range corrected PBE (LRC-wPBE) by by Rohrdanz, Martins and Herbert { @@ -236,7 +234,7 @@ const std::vector in_built_xc_func_ext_params(const int id) // using the error function kernel. return { std::stod(PARAM.inp.exx_fock_alpha[0]), //Fraction of Hartree-Fock exchange: 1.0 std::stod(PARAM.inp.exx_erfc_alpha[0]), //Fraction of short-range exact exchange: -1.0 - GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.3 + hse_omega }; //Range separation constant: 0.3 } case XC_HYB_GGA_XC_LRC_WPBEH: // Long-range corrected short-range hybrid PBE (LRC-wPBEh) by Rohrdanz, Martins and Herbert { @@ -245,7 +243,7 @@ const std::vector in_built_xc_func_ext_params(const int id) // using the error function kernel. return { std::stod(PARAM.inp.exx_fock_alpha[0]), //Fraction of Hartree-Fock exchange: 1.0 std::stod(PARAM.inp.exx_erfc_alpha[0]), //Fraction of short-range exact exchange: -0.8 - GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.2 + hse_omega }; //Range separation constant: 0.2 } case XC_HYB_GGA_XC_CAM_PBEH: // CAM hybrid screened exchange PBE version { @@ -254,7 +252,7 @@ const std::vector in_built_xc_func_ext_params(const int id) // using the error function kernel. return { std::stod(PARAM.inp.exx_fock_alpha[0]), //Fraction of Hartree-Fock exchange: 0.2 std::stod(PARAM.inp.exx_erfc_alpha[0]), //Fraction of short-range exact exchange: 0.8 - GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.7 + hse_omega }; //Range separation constant: 0.7 } #endif default: @@ -282,7 +280,9 @@ const std::vector external_xc_func_ext_params(const int id) std::vector XC_Functional_Libxc::init_func(const std::vector &func_id, - const int xc_polarized) + const int xc_polarized, + const double hybrid_alpha, + const double hse_omega) { std::vector funcs; for (int id : func_id) @@ -291,7 +291,7 @@ XC_Functional_Libxc::init_func(const std::vector &func_id, xc_func_init(&funcs.back(), id, xc_polarized); // instantiate the XC term // search for external parameters - const std::vector in_built_ext_params = in_built_xc_func_ext_params(id); + const std::vector in_built_ext_params = in_built_xc_func_ext_params(id, hybrid_alpha, hse_omega); const std::vector external_ext_params = external_xc_func_ext_params(id); // for temporary use, I name their size as n1 and n2 const int n1 = in_built_ext_params.size(); diff --git a/source/source_hamilt/module_xc/test/test_xc.cpp b/source/source_hamilt/module_xc/test/test_xc.cpp index 0efd997bfa9..72b2f547bc8 100644 --- a/source/source_hamilt/module_xc/test/test_xc.cpp +++ b/source/source_hamilt/module_xc/test/test_xc.cpp @@ -881,13 +881,15 @@ class XCTest_PBE_LibXC : public XCTest std::vector rho = {0.17E+01, 0.17E+01, 0.15E+01, 0.88E-01, 0.18E+04}; std::vector grho = {0.81E-11, 0.17E+01, 0.36E+02, 0.87E-01, 0.55E+00}; + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); + const double hse_omega = XC_Functional::get_hse_omega(); for(int i=0;i<5;i++) { double e,v,v1,v2; XC_Functional::xc(rho[i],e,v); e_lda.push_back(e); v_lda.push_back(v); - XC_Functional_Libxc::gcxc_libxc(XC_Functional::get_func_id(), rho[i],grho[i],e,v1,v2); + XC_Functional_Libxc::gcxc_libxc(XC_Functional::get_func_id(), rho[i],grho[i],e,v1,v2, hybrid_alpha, hse_omega); e_gga.push_back(e); v1_gga.push_back(v1); v2_gga.push_back(v2); diff --git a/source/source_hamilt/module_xc/test/test_xc2.cpp b/source/source_hamilt/module_xc/test/test_xc2.cpp index 3925df7d92d..5a6401e1145 100644 --- a/source/source_hamilt/module_xc/test/test_xc2.cpp +++ b/source/source_hamilt/module_xc/test/test_xc2.cpp @@ -458,7 +458,9 @@ class XCTest_PBE_SPN_LibXC : public XCTest double e,v1,v2,v3,v4,v5; double r1 = rho[i] * (1+zeta[i]) / 2.0; double r2 = rho[i] * (1-zeta[i]) / 2.0; - XC_Functional_Libxc::gcxc_spin_libxc(XC_Functional::get_func_id(), r1,r2,gdr[i],gdr[i],e,v1,v2,v3,v4,v5); + double hybrid_alpha = 0.0; + double hse_omega = 0.0; + XC_Functional_Libxc::gcxc_spin_libxc(XC_Functional::get_func_id(), r1,r2,gdr[i],gdr[i],e,v1,v2,v3,v4,v5, hybrid_alpha, hse_omega); e_gga.push_back(e); v1_gga.push_back(v1+v3); v2_gga.push_back(v2+v4); @@ -492,12 +494,14 @@ class XCTest_PZ_SPN_LibXC : public XCTest std::vector rho = {-1, 0.17E+01, 0.17E+01, 0.15E+01, 0.88E-01, 0.18E+04}; std::vector zeta = {0.0, 0.0, 0.2, 0.5, 0.8, 1.0}; + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); + const double hse_omega = XC_Functional::get_hse_omega(); for(int i=0;i<5;i++) { double e,v1,v2; double r1 = rho[i] * (1+zeta[i]) / 2.0; double r2 = rho[i] * (1-zeta[i]) / 2.0; - XC_Functional_Libxc::xc_spin_libxc(XC_Functional::get_func_id(), r1,r2,e,v1,v2); + XC_Functional_Libxc::xc_spin_libxc(XC_Functional::get_func_id(), r1,r2,e,v1,v2, hybrid_alpha, hse_omega); e_lda.push_back(e); v1_lda.push_back(v1); v2_lda.push_back(v2); diff --git a/source/source_hamilt/module_xc/test/test_xc3.cpp b/source/source_hamilt/module_xc/test/test_xc3.cpp index 71a4aaf27a6..43e1e895355 100644 --- a/source/source_hamilt/module_xc/test/test_xc3.cpp +++ b/source/source_hamilt/module_xc/test/test_xc3.cpp @@ -89,13 +89,15 @@ class XCTest_GRADCORR : public XCTest XC_Functional::set_xc_type("PBE"); - 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); + double hybrid_alpha = 0.0; + double hse_omega = 0.0; + XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,false,nspin1,domag,domag_z, hybrid_alpha, hse_omega); + XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,true,nspin1,domag,domag_z, hybrid_alpha, hse_omega); - 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); + XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,false,nspin2,domag,domag_z, hybrid_alpha, hse_omega); + XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,true,nspin2,domag,domag_z, hybrid_alpha, hse_omega); - XC_Functional::gradcorr(et4,vt4,v4,&chr,&rhopw,&ucell,stress4,false,nspin4,domag_true,domag_z); + XC_Functional::gradcorr(et4,vt4,v4,&chr,&rhopw,&ucell,stress4,false,nspin4,domag_true,domag_z, hybrid_alpha, hse_omega); } }; diff --git a/source/source_hamilt/module_xc/test/test_xc4.cpp b/source/source_hamilt/module_xc/test/test_xc4.cpp index 467cc2667cd..e1ef9a82c08 100644 --- a/source/source_hamilt/module_xc/test/test_xc4.cpp +++ b/source/source_hamilt/module_xc/test/test_xc4.cpp @@ -47,7 +47,8 @@ class XCTest_SCAN : public XCTest { double e,v,v1,v2,v3; double hybrid_alpha = 0.0; - XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3,hybrid_alpha); + double hse_omega = 0.0; + XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3,hybrid_alpha, hse_omega); e_.push_back(e); v1_.push_back(v1); v2_.push_back(v2); diff --git a/source/source_hamilt/module_xc/test/test_xc5.cpp b/source/source_hamilt/module_xc/test/test_xc5.cpp index 60811051364..e6638705c31 100644 --- a/source/source_hamilt/module_xc/test/test_xc5.cpp +++ b/source/source_hamilt/module_xc/test/test_xc5.cpp @@ -78,14 +78,16 @@ class XCTest_VXC : public XCTest XC_Functional::set_xc_type("PBE"); + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); + const double hse_omega = XC_Functional::get_hse_omega(); std::tuple etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin1,domag,domag_z); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin1,domag,domag_z, hybrid_alpha, hse_omega); et1 = std::get<0>(etxc_vtxc_v); vt1 = std::get<1>(etxc_vtxc_v); v1 = std::get<2>(etxc_vtxc_v); etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin2,domag,domag_z); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin2,domag,domag_z, hybrid_alpha, hse_omega); et2 = std::get<0>(etxc_vtxc_v); vt2 = std::get<1>(etxc_vtxc_v); v2 = std::get<2>(etxc_vtxc_v); @@ -180,14 +182,16 @@ class XCTest_VXC_Libxc : public XCTest XC_Functional::set_xc_type("GGA_X_PBE+GGA_C_PBE"); + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); + const double hse_omega = XC_Functional::get_hse_omega(); std::tuple etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin1,domag,domag_z); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin1,domag,domag_z, hybrid_alpha, hse_omega); et1 = std::get<0>(etxc_vtxc_v); vt1 = std::get<1>(etxc_vtxc_v); v1 = std::get<2>(etxc_vtxc_v); etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin2,domag,domag_z); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin2,domag,domag_z, hybrid_alpha, hse_omega); et2 = std::get<0>(etxc_vtxc_v); vt2 = std::get<1>(etxc_vtxc_v); v2 = std::get<2>(etxc_vtxc_v); @@ -290,15 +294,17 @@ class XCTest_VXC_meta : public XCTest XC_Functional::set_xc_type("SCAN"); + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); + const double hse_omega = XC_Functional::get_hse_omega(); std::tuple etxc_vtxc_v - = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr,nspin1); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr,nspin1, hybrid_alpha, hse_omega); et1 = std::get<0>(etxc_vtxc_v); vt1 = std::get<1>(etxc_vtxc_v); v1 = std::get<2>(etxc_vtxc_v); vtau1 = std::get<3>(etxc_vtxc_v); etxc_vtxc_v - = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr,nspin2); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr,nspin2, hybrid_alpha, hse_omega); 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/xc_functional.cpp b/source/source_hamilt/module_xc/xc_functional.cpp index 4290b4ff080..bef9d91b2b8 100644 --- a/source/source_hamilt/module_xc/xc_functional.cpp +++ b/source/source_hamilt/module_xc/xc_functional.cpp @@ -16,6 +16,7 @@ int XC_Functional::func_type = 0; bool XC_Functional::ked_flag = false; bool XC_Functional::use_libxc = true; double XC_Functional::hybrid_alpha = 0.25; +double XC_Functional::hse_omega = 0.0; std::map XC_Functional::scaling_factor_xc = { {1, 1.0} }; // added by jghan, 2024-10-10 void XC_Functional::set_hybrid_alpha(const double alpha_in) @@ -23,6 +24,11 @@ void XC_Functional::set_hybrid_alpha(const double alpha_in) hybrid_alpha = alpha_in; } +void XC_Functional::set_hse_omega(const double omega_in) +{ + hse_omega = omega_in; +} + void XC_Functional::set_xc_first_loop(const UnitCell& ucell) { ModuleBase::TITLE("XC_Functional", "set_xc_first_loop"); @@ -337,7 +343,9 @@ std::string XC_Functional::output_info() ss<<" Libxc v"< funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED); + double hybrid_alpha = 0.0; + double hse_omega = 0.0; + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, hybrid_alpha, hse_omega); for(const auto &func : funcs) { const xc_func_info_type *info = xc_func_get_info(&func); diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index aaba3e945b7..ceb81353992 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -50,7 +50,9 @@ class XC_Functional const UnitCell *ucell, // charge density const int nspin, const bool domag, - const bool domag_z); + const bool domag_z, + const double hybrid_alpha, + const double hse_omega); //------------------- // xc_functional.cpp @@ -80,6 +82,13 @@ class XC_Functional return hybrid_alpha; }; + static void set_hse_omega(const double omega_in); + + static double get_hse_omega() + { + return hse_omega; + }; + static bool get_ked_flag() { return ked_flag; @@ -100,6 +109,9 @@ class XC_Functional // exx_hybrid_alpha for mixing exx in hybrid functional: static double hybrid_alpha; + // hse_omega for HSE functional: + static double hse_omega; + // added by jghan, 2024-07-07 // as a scaling factor for different xc-functionals static std::map scaling_factor_xc; @@ -212,7 +224,9 @@ class XC_Functional const bool is_stress, const int nspin, const bool domag, - const bool domag_z); + const bool domag_z, + const double hybrid_alpha, + const double hse_omega); template ::type> diff --git a/source/source_hamilt/module_xc/xc_grad.cpp b/source/source_hamilt/module_xc/xc_grad.cpp index 8104c03a61c..c16dcc718bd 100644 --- a/source/source_hamilt/module_xc/xc_grad.cpp +++ b/source/source_hamilt/module_xc/xc_grad.cpp @@ -35,7 +35,9 @@ void XC_Functional::gradcorr( const bool is_stress, const int nspin, const bool domag, - const bool domag_z) + const bool domag_z, + const double hybrid_alpha_in, + const double hse_omega_in) { ModuleBase::TITLE("XC_Functional","gradcorr"); @@ -301,15 +303,11 @@ void XC_Functional::gradcorr( { double v3xc = 0.0; double atau = chr->kin_r[0][ir]/2.0; - double hybrid_alpha = 0.0; -#ifdef __EXX - hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; -#endif - XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc, hybrid_alpha); + XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc, hybrid_alpha_in, hse_omega_in); } else { - XC_Functional_Libxc::gcxc_libxc( func_id, arho, grho2a, sxc, v1xc, v2xc); + XC_Functional_Libxc::gcxc_libxc( func_id, arho, grho2a, sxc, v1xc, v2xc, hybrid_alpha_in, hse_omega_in); } #endif } @@ -370,21 +368,18 @@ void XC_Functional::gradcorr( double v3xcdw = 0.0; double atau1 = chr->kin_r[0][ir]/2.0; double atau2 = chr->kin_r[1][ir]/2.0; - double hybrid_alpha = 0.0; -#ifdef __EXX - hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; -#endif XC_Functional_Libxc::tau_xc_spin( func_id, rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], - atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, hybrid_alpha); + atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, hybrid_alpha_in, hse_omega_in); } else { XC_Functional_Libxc::gcxc_spin_libxc( func_id, rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], - sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud); + sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, + hybrid_alpha_in, hse_omega_in); } if(is_stress) { diff --git a/source/source_hamilt/module_xc/xc_pot.cpp b/source/source_hamilt/module_xc/xc_pot.cpp index 6e56f642a16..a0e905e5915 100644 --- a/source/source_hamilt/module_xc/xc_pot.cpp +++ b/source/source_hamilt/module_xc/xc_pot.cpp @@ -11,6 +11,9 @@ #ifdef USE_LIBXC #include "libxc_abacus.h" +#ifdef __EXX +#include "source_hamilt/module_xc/exx_info.h" +#endif #endif // [etxc, vtxc, v] = XC_Functional::v_xc(...) @@ -20,7 +23,9 @@ std::tuple XC_Functional::v_xc( const UnitCell* ucell, const int nspin, const bool domag, - const bool domag_z) + const bool domag_z, + const double hybrid_alpha, + const double hse_omega) { ModuleBase::TITLE("XC_Functional", "v_xc"); @@ -35,7 +40,9 @@ std::tuple XC_Functional::v_xc( nspin, domag, domag_z, - &(scaling_factor_xc)); + &(scaling_factor_xc), + hybrid_alpha, + hse_omega); #else ModuleBase::WARNING_QUIT("v_xc", "compile with LIBXC"); #endif @@ -140,7 +147,7 @@ std::tuple XC_Functional::v_xc( #ifdef USE_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]); + 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"); #endif @@ -177,7 +184,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, false, nspin, domag, domag_z); + gradcorr(etxc, vtxc, v, chr, chr->rhopw, ucell, dum, false, nspin, domag, domag_z, hybrid_alpha, hse_omega); // parallel code : collect vtxc,etxc // mohan add 2008-06-01 diff --git a/source/source_hsolver/hsolver_lcaopw.cpp b/source/source_hsolver/hsolver_lcaopw.cpp index cb0f31a4b3e..f6ce57ea2d9 100644 --- a/source/source_hsolver/hsolver_lcaopw.cpp +++ b/source/source_hsolver/hsolver_lcaopw.cpp @@ -44,21 +44,23 @@ void HSolverLIP::solve(hamilt::Hamilt* pHamilt, // ESolver_KS_PW::p_hamilt #ifdef __EXX auto& exx_lip = dynamic_cast*>(pHamilt)->exx_lip; - auto add_exx_to_subspace_hamilt = [&ik, &exx_lip](T* hcc, const int naos) -> void { - if (GlobalC::exx_info.info_global.cal_exx) + bool cal_exx = GlobalC::exx_info.info_global.cal_exx; + double hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; + auto add_exx_to_subspace_hamilt = [&ik, &exx_lip, cal_exx, hybrid_alpha](T* hcc, const int naos) -> void { + if (cal_exx) { for (int n = 0; n < naos; ++n) { for (int m = 0; m < naos; ++m) { hcc[n * naos + m] - += (T)GlobalC::exx_info.info_global.hybrid_alpha * exx_lip.get_exx_matrix()[ik][m][n]; + += (T)hybrid_alpha * exx_lip.get_exx_matrix()[ik][m][n]; } } } }; - auto set_exxlip_lcaowfc = [&ik, &exx_lip](const T* const vcc, const int naos, const int nbands) -> void { - if (GlobalC::exx_info.info_global.cal_exx) + auto set_exxlip_lcaowfc = [&ik, &exx_lip, cal_exx](const T* const vcc, const int naos, const int nbands) -> void { + if (cal_exx) { exx_lip.set_hvec(ik, vcc, naos, nbands); } diff --git a/source/source_io/module_chgpot/write_libxc_r.cpp b/source/source_io/module_chgpot/write_libxc_r.cpp index 6373e283814..0579650aa61 100644 --- a/source/source_io/module_chgpot/write_libxc_r.cpp +++ b/source/source_io/module_chgpot/write_libxc_r.cpp @@ -50,10 +50,14 @@ void ModuleIO::write_libxc_r( // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ //---------------------------------------------------------- + double hybrid_alpha = 0.0; + double hse_omega = 0.0; std::vector funcs = XC_Functional_Libxc::init_func( func_id, - (1==nspin) ? XC_UNPOLARIZED : XC_POLARIZED + (1==nspin) ? XC_UNPOLARIZED : XC_POLARIZED, + hybrid_alpha, + hse_omega ); const bool is_gga = [&funcs]() diff --git a/source/source_io/module_ctrl/ctrl_iter_lcao.cpp b/source/source_io/module_ctrl/ctrl_iter_lcao.cpp index a4de81385f2..7fd133d414b 100644 --- a/source/source_io/module_ctrl/ctrl_iter_lcao.cpp +++ b/source/source_io/module_ctrl/ctrl_iter_lcao.cpp @@ -47,12 +47,14 @@ void ctrl_iter_lcao(UnitCell& ucell, // unit cell * } #ifdef __EXX - // save exx matrix + bool cal_exx = GlobalC::exx_info.info_global.cal_exx; + bool real_number = GlobalC::exx_info.info_ri.real_number; + if (inp.calculation != "nscf") { - if (GlobalC::exx_info.info_global.cal_exx) + if (cal_exx) { - GlobalC::exx_info.info_ri.real_number ? + real_number ? exx_nao.exd->exx_iter_finish(kv, ucell, *p_hamilt, *pelec, &dm, *p_chgmix, scf_ene_thr, iter, istep, conv_esolver) : exx_nao.exc->exx_iter_finish(kv, ucell, *p_hamilt, *pelec, &dm, diff --git a/source/source_io/module_ctrl/ctrl_runner_lcao.cpp b/source/source_io/module_ctrl/ctrl_runner_lcao.cpp index 849b2925555..0fa74574c64 100644 --- a/source/source_io/module_ctrl/ctrl_runner_lcao.cpp +++ b/source/source_io/module_ctrl/ctrl_runner_lcao.cpp @@ -2,6 +2,7 @@ #include "source_estate/elecstate_lcao.h" // use elecstate::ElecState #include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO +#include "source_hamilt/module_xc/exx_info.h" #include "../module_energy/write_proj_band_lcao.h" // projcted band structure #include "../module_dos/cal_ldos.h" // cal LDOS @@ -55,6 +56,7 @@ void ctrl_runner_lcao(UnitCell& ucell, // unitcell // 3) print out exchange-correlation potential if (inp.out_mat_xc) { + bool cal_exx = GlobalC::exx_info.info_global.cal_exx; ModuleIO::write_Vxc(inp.nspin, PARAM.globalv.nlocal, GlobalV::DRANK, @@ -70,7 +72,8 @@ void ctrl_runner_lcao(UnitCell& ucell, // unitcell kv, orb.cutoffs(), pelec->wg, - gd + gd, + cal_exx #ifdef __EXX , exx_nao.exd ? &exx_nao.exd->get_Hexxs() : nullptr, @@ -81,6 +84,9 @@ void ctrl_runner_lcao(UnitCell& ucell, // unitcell if (inp.out_mat_xc2[0]) { + bool cal_exx = GlobalC::exx_info.info_global.cal_exx; + double hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; + bool real_number = GlobalC::exx_info.info_ri.real_number; ModuleIO::write_Vxc_R(inp.nspin, &pv, ucell, @@ -92,7 +98,10 @@ void ctrl_runner_lcao(UnitCell& ucell, // unitcell chr, kv, orb.cutoffs(), - gd + gd, + cal_exx, + hybrid_alpha, + real_number #ifdef __EXX , exx_nao.exd ? &exx_nao.exd->get_Hexxs() : nullptr, diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp index d109b667aaf..35a280139d4 100644 --- a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp +++ b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp @@ -441,12 +441,15 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, //! 15) Output Hexx matrix in LCAO basis // (see `out_chg` in docs/advanced/input_files/input-main.md) //------------------------------------------------------------------ + bool cal_exx = GlobalC::exx_info.info_global.cal_exx; + bool real_number = GlobalC::exx_info.info_ri.real_number; + if (inp.out_chg[0]) { - if (GlobalC::exx_info.info_global.cal_exx && inp.calculation != "nscf") // Peize Lin add if 2022.11.14 + if (cal_exx && inp.calculation != "nscf") // Peize Lin add if 2022.11.14 { const std::string file_name_exx = global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK); - if (GlobalC::exx_info.info_ri.real_number) + if (real_number) { ModuleIO::write_Hexxs_csr(file_name_exx, ucell, exx_nao.exd->get_Hexxs()); } diff --git a/source/source_io/module_current/td_current_io_comm.cpp b/source/source_io/module_current/td_current_io_comm.cpp index b987643d49a..68b05012397 100644 --- a/source/source_io/module_current/td_current_io_comm.cpp +++ b/source/source_io/module_current/td_current_io_comm.cpp @@ -148,8 +148,6 @@ void ModuleIO::set_rR_from_hR(const UnitCell& ucell, auto col_indexes = pv->get_indexes_col(iat2); const ModuleBase::Vector3& tau1 = ucell.get_tau(iat1); - // std::cout << "tau1: " << tau1 << " tau2: " << GlobalC::ucell.get_tau(iat2) << " r_index: " << r_index - // << std::endl; const ModuleBase::Vector3 tau2 = tau1 + dtau; for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) { @@ -389,19 +387,6 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell, { hamilt::folding_HR(sR, sk, kv.kvec_d[ik], nrow, 1); } - // for (int ir = 0; ir < pv->nrow; ir++) - // { - // const int iwt1 = pv->local2global_row(ir); - // const int iat1 = GlobalC::ucell.iwt2iat[iwt1]; - // for (int ic = 0; ic < pv->ncol; ic++) - // { - // const int iwt2 = pv->local2global_col(ic); - // const int iat2 = GlobalC::ucell.iwt2iat[iwt2]; - // const int irc = ic * pv->nrow + ir; - // std::cout << "ik: " << ik << " iat1:" << iat1 << " iat2:" << iat2 << " iwt1: " << iwt1 - // << " iwt2: " << iwt2 << " hk: " << hk[irc] << std::endl; - // } - // } // 2. set inverse S(k) -> sk will be changed to sk_inv int* ipiv = new int[pv->nloc]; int info = 0; @@ -470,22 +455,6 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell, { module_rt::folding_partial_HR(ucell, sR, partial_sk, kv.kvec_d[ik], i_alpha, nrow, 1); } - // if(i_alpha == 2) - // { - // for(int ir=0;ir< pv->nrow; ir++) - // { - // const int iwt1 = pv->local2global_row(ir); - // const int iat1 = GlobalC::ucell.iwt2iat[iwt1]; - // for(int ic=0;ic< pv->ncol; ic++) - // { - // const int iwt2 = pv->local2global_col(ic); - // const int iat2 = GlobalC::ucell.iwt2iat[iwt2]; - // const int irc=ic*pv->nrow + ir; - // std::cout<<"ik: "<nloc); @@ -498,24 +467,6 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell, { hamilt::folding_HR(*rR[i_alpha], rk, kv.kvec_d[ik], nrow, 1); // set r(k) } - // if (i_alpha == 2) - // { - // std::cout << "ik: " << ik << " i_alpha: " << i_alpha << std::endl; - // for (int ir = 0; ir < pv->nrow; ir++) - // { - // const int iwt1 = pv->local2global_row(ir); - // const int iat1 = GlobalC::ucell.iwt2iat[iwt1]; - // for (int ic = 0; ic < pv->ncol; ic++) - // { - // const int iwt2 = pv->local2global_col(ic); - // const int iat2 = GlobalC::ucell.iwt2iat[iwt2]; - // const int irc = ic * pv->nrow + ir; - // std::cout << " iat1: " << iat1 << " iat2: " << iat2 << " iw1: " << - // GlobalC::ucell.iwt2iw[iwt1] - // << " iw2: " << GlobalC::ucell.iwt2iw[iwt2] << " rk: " << rk[irc] << std::endl; - // } - // } - // } // 4. calculate <\vu,k|v_a|\mu,k> = partial_Hk + IMAG_UNIT * (Hk * Sk_inv * rk) - IMAG_UNIT * (rk * Sk_inv * // Hk) - Hk * Sk_inv * partial_Sk // 4.1.1 Hk * Sk_inv (note 2.) @@ -667,23 +618,7 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell, pv->desc); // 5. copy h_is_ps to velocity_basis_k[ik][i_alpha] BlasConnector::copy(pv->nloc, h_is_ps, 1, velocity_basis_k[ik][i_alpha], 1); - // if(i_alpha == 2) - // { - // for(int ir=0;ir< pv->nrow; ir++) - // { - // const int iwt1 = pv->local2global_row(ir); - // const int iat1 = GlobalC::ucell.iwt2iat[iwt1]; - // for(int ic=0;ic< pv->ncol; ic++) - // { - // const int iwt2 = pv->local2global_col(ic); - // const int iat2 = GlobalC::ucell.iwt2iat[iwt2]; - // const int irc=ic*pv->nrow + ir; - // std::cout<<"ik: "<>* psi, for (int ir = 0; ir < PARAM.inp.nbands; ++ir) { - // const int iwt1 = pv->local2global_row(ir); - // const int iat1 = GlobalC::ucell.iwt2iat[iwt1]; for (int ic = 0; ic < PARAM.inp.nbands; ++ic) { const int irc = ic * pv->nrow + ir; if (pv->in_this_processor(ir, ic)) { - // const int iwt2 = pv->local2global_col(ic); - // const int iat2 = GlobalC::ucell.iwt2iat[iwt2]; velocity_k[ik][i_alpha](ir, ic) = vk_c[irc]; - // if (i_alpha == 0) - // { - // std::cout<<"ik: "<(nspin, nbasis, drank, @@ -196,7 +198,8 @@ void write_eband_terms(const int nspin, kv, orb_cutoff, wg, - gd + gd, + cal_exx #ifdef __EXX , Hexxd, diff --git a/source/source_io/module_hs/write_vxc.hpp b/source/source_io/module_hs/write_vxc.hpp index 881f267e43e..d4719a17874 100644 --- a/source/source_io/module_hs/write_vxc.hpp +++ b/source/source_io/module_hs/write_vxc.hpp @@ -152,7 +152,8 @@ void write_Vxc(const int nspin, const K_Vectors& kv, const std::vector& orb_cutoff, const ModuleBase::matrix& wg, - Grid_Driver& gd + Grid_Driver& gd, + bool cal_exx #ifdef __EXX , std::vector>>>* Hexxd = nullptr, @@ -224,7 +225,7 @@ void write_Vxc(const int nspin, const std::vector& vlocxc_k_mo = cVc(vxc_k_ao.get_hk(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); #ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) + if (cal_exx) { e_orb_locxc.emplace_back(orbital_energy(ik, nbands, vlocxc_k_mo, p2d)); ModuleBase::GlobalFunc::ZEROS(vexxonly_k_ao.get_hk(), pv->nloc); @@ -232,9 +233,6 @@ void write_Vxc(const int nspin, vexxonly_op_ao.contributeHk(ik); std::vector vexx_k_mo = cVc(vexxonly_k_ao.get_hk(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); e_orb_exx.emplace_back(orbital_energy(ik, nbands, vexx_k_mo, p2d)); - // ======test======= - // exx_energy += all_band_energy(ik, vexx_k_mo, p2d, wg); - // ======test======= } #endif if (PARAM.inp.dft_plus_u) @@ -289,7 +287,7 @@ void write_Vxc(const int nspin, { write_orb_energy(kv, nspin0, nbands, e_orb_tot, "vxc", ""); #ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) + if (cal_exx) { write_orb_energy(kv, nspin0, nbands, e_orb_locxc, "vxc", "local"); write_orb_energy(kv, nspin0, nbands, e_orb_exx, "vxc", "exx"); diff --git a/source/source_io/module_hs/write_vxc_lip.hpp b/source/source_io/module_hs/write_vxc_lip.hpp index 30a7e043d01..5909094a2f7 100644 --- a/source/source_io/module_hs/write_vxc_lip.hpp +++ b/source/source_io/module_hs/write_vxc_lip.hpp @@ -105,7 +105,6 @@ namespace ModuleIO int naos, int drank, const psi::Psi>& psi_pw, - // const psi::Psi& psi_lcao, const UnitCell& ucell, Structure_Factor& sf, surchem& solvent, @@ -115,7 +114,9 @@ namespace ModuleIO const ModuleBase::matrix& vloc, const Charge& chg, const K_Vectors& kv, - const ModuleBase::matrix& wg + const ModuleBase::matrix& wg, + bool cal_exx, + double hybrid_alpha #ifdef __EXX , const Exx_Lip>& exx_lip @@ -186,24 +187,19 @@ namespace ModuleIO std::vector vxc_tot_k_mo(std::move(vxc_local_k_mo)); std::vector vexx_k_ao(naos * naos); #if((defined __LCAO)&&(defined __EXX) && !(defined __CUDA)&& !(defined __ROCM)) - if (GlobalC::exx_info.info_global.cal_exx) + if (cal_exx) { for (int n = 0; n < naos; ++n) { for (int m = 0; m < naos; ++m) { - vexx_k_ao[n * naos + m] += (T)GlobalC::exx_info.info_global.hybrid_alpha + vexx_k_ao[n * naos + m] += (T)hybrid_alpha * exx_lip.get_exx_matrix()[ik][m][n]; } } std::vector vexx_k_mo = cVc(vexx_k_ao.data(), &(exx_lip.get_hvec()(ik, 0, 0)), naos, nbands); Parallel_Reduce::reduce_pool(vexx_k_mo.data(), nbands * nbands); e_orb_exx.emplace_back(orbital_energy(ik, nbands, vexx_k_mo)); - // ======test======= - // std::cout << "exx_energy from matrix:" << all_band_energy(ik, nbands, vexx_k_mo, wg) << std::endl; - // std::cout << "exx_energy from orbitals: " << all_band_energy(ik, e_orb_exx.at(ik), wg) << std::endl; - // std::cout << "exx_energy from exx_lip: " << GlobalC::exx_info.info_global.hybrid_alpha * exx_lip.get_exx_energy() << std::endl; - // ======test======= container::BlasConnector::axpy(nbands * nbands, 1.0, vexx_k_mo.data(), 1, vxc_tot_k_mo.data(), 1); } #endif @@ -243,18 +239,6 @@ namespace ModuleIO // std::cout << "xc all-bands energy by rho =" << exc_by_rho << std::endl; //===== test total xc energy ======= //===== test total exx energy ======= -// #if((defined __LCAO)&&(defined __EXX) && !(defined __CUDA)&& !(defined __ROCM)) -// if (GlobalC::exx_info.info_global.cal_exx) -// { -// FPTYPE exx_by_orb = 0.0; -// for (int ik = 0;ik < e_orb_exx.size();++ik) -// exx_by_orb += all_band_energy(ik, e_orb_exx[ik], wg); -// exx_by_orb /= 2; -// std::cout << "exx all-bands energy by orbital =" << exx_by_orb << std::endl; -// FPTYPE exx_from_lip = GlobalC::exx_info.info_global.hybrid_alpha * exx_lip.get_exx_energy(); -// std::cout << "exx all-bands energy from exx_lip =" << exx_from_lip << std::endl; -// } -// #endif //===== test total exx energy ======= // write the orbital energy for xc and exx in LibRPA format const int nspin0 = (nspin == 2) ? 2 : 1; @@ -284,7 +268,7 @@ namespace ModuleIO { write_orb_energy(e_orb_tot, ""); #if((defined __LCAO)&&(defined __EXX) && !(defined __CUDA)&& !(defined __ROCM)) - if (GlobalC::exx_info.info_global.cal_exx) + if (cal_exx) { write_orb_energy(e_orb_locxc, "local"); write_orb_energy(e_orb_exx, "exx"); diff --git a/source/source_io/module_hs/write_vxc_r.hpp b/source/source_io/module_hs/write_vxc_r.hpp index a06ec6e8062..81d4b38a102 100644 --- a/source/source_io/module_hs/write_vxc_r.hpp +++ b/source/source_io/module_hs/write_vxc_r.hpp @@ -36,10 +36,15 @@ void write_Vxc_R(const int nspin, const K_Vectors& kv, const std::vector& orb_cutoff, Grid_Driver& gd, + bool cal_exx, + double hybrid_alpha, + bool real_number #ifdef __EXX + , const std::vector>>>* const Hexxd, - const std::vector>>>>* const Hexxc, + const std::vector>>>>* const Hexxc #endif + , const double sparse_thr = 1e-10) { ModuleBase::TITLE("ModuleIO", "write_Vxc_R"); @@ -69,9 +74,9 @@ void write_Vxc_R(const int nspin, vxcs_R_ao[is].fix_gamma(); } #ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) + if (cal_exx) { - GlobalC::exx_info.info_ri.real_number + real_number ? hamilt::reallocate_hcontainer(*Hexxd, &vxcs_R_ao[is], &cell_nearest) : hamilt::reallocate_hcontainer(*Hexxc, &vxcs_R_ao[is], &cell_nearest); } @@ -93,22 +98,22 @@ void write_Vxc_R(const int nspin, vxcs_op_ao.set_current_spin(is); vxcs_op_ao.contributeHR(); #ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) + if (cal_exx) { - GlobalC::exx_info.info_ri.real_number ? RI_2D_Comm::add_HexxR(is, - GlobalC::exx_info.info_global.hybrid_alpha, - *Hexxd, - *pv, - ucell.get_npol(), - vxcs_R_ao[is], - &cell_nearest) - : RI_2D_Comm::add_HexxR(is, - GlobalC::exx_info.info_global.hybrid_alpha, - *Hexxc, - *pv, - ucell.get_npol(), - vxcs_R_ao[is], - &cell_nearest); + real_number ? RI_2D_Comm::add_HexxR(is, + hybrid_alpha, + *Hexxd, + *pv, + ucell.get_npol(), + vxcs_R_ao[is], + &cell_nearest) + : RI_2D_Comm::add_HexxR(is, + hybrid_alpha, + *Hexxc, + *pv, + ucell.get_npol(), + vxcs_R_ao[is], + &cell_nearest); } #endif } diff --git a/source/source_io/module_parameter/input_conv.cpp b/source/source_io/module_parameter/input_conv.cpp index 6a1d06c0d46..703485e1d0b 100644 --- a/source/source_io/module_parameter/input_conv.cpp +++ b/source/source_io/module_parameter/input_conv.cpp @@ -479,6 +479,7 @@ void Input_Conv::Convert() XC_Functional::set_hybrid_alpha(GlobalC::exx_info.info_global.hybrid_alpha); if(!PARAM.inp.exx_erfc_omega.empty()) { GlobalC::exx_info.info_global.hse_omega = std::stod(PARAM.inp.exx_erfc_omega[0]); } + XC_Functional::set_hse_omega(GlobalC::exx_info.info_global.hse_omega); if(!PARAM.inp.exx_fock_lambda.empty()) { GlobalC::exx_info.info_lip.lambda = std::stod(PARAM.inp.exx_fock_lambda[0]); } GlobalC::exx_info.info_global.separate_loop = PARAM.inp.exx_separate_loop; diff --git a/source/source_lcao/FORCE_STRESS.cpp b/source/source_lcao/FORCE_STRESS.cpp index c5a0ca6e6b2..d24c563e568 100644 --- a/source/source_lcao/FORCE_STRESS.cpp +++ b/source/source_lcao/FORCE_STRESS.cpp @@ -478,35 +478,38 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, // } #ifdef __EXX - // Force and Stress contribution from exx + bool cal_exx = GlobalC::exx_info.info_global.cal_exx; + bool real_number = GlobalC::exx_info.info_ri.real_number; + double hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; + ModuleBase::matrix force_exx; ModuleBase::matrix stress_exx; - if (GlobalC::exx_info.info_global.cal_exx) + if (cal_exx) { if (isforce) { - if (GlobalC::exx_info.info_ri.real_number) + if (real_number) { exx_nao.exd->cal_exx_force(ucell.nat); - force_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_nao.exd->get_force(); + force_exx = hybrid_alpha * exx_nao.exd->get_force(); } else { exx_nao.exc->cal_exx_force(ucell.nat); - force_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_nao.exc->get_force(); + force_exx = hybrid_alpha * exx_nao.exc->get_force(); } } if (isstress) { - if (GlobalC::exx_info.info_ri.real_number) + if (real_number) { exx_nao.exd->cal_exx_stress(ucell.omega, ucell.lat0); - stress_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_nao.exd->get_stress(); + stress_exx = hybrid_alpha * exx_nao.exd->get_stress(); } else { exx_nao.exc->cal_exx_stress(ucell.omega, ucell.lat0); - stress_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_nao.exc->get_stress(); + stress_exx = hybrid_alpha * exx_nao.exc->get_stress(); } } } diff --git a/source/source_lcao/module_lr/potentials/xc_kernel.cpp b/source/source_lcao/module_lr/potentials/xc_kernel.cpp index c5d5f9bb667..8a6d36cf765 100644 --- a/source/source_lcao/module_lr/potentials/xc_kernel.cpp +++ b/source/source_lcao/module_lr/potentials/xc_kernel.cpp @@ -120,9 +120,13 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl assert(nspin == 1 || nspin == 2); + double hybrid_alpha = 0.0; + double hse_omega = 0.0; std::vector funcs = XC_Functional_Libxc::init_func( XC_Functional::get_func_id(), - (1 == nspin) ? XC_UNPOLARIZED : XC_POLARIZED); + (1 == nspin) ? XC_UNPOLARIZED : XC_POLARIZED, + hybrid_alpha, + hse_omega); const int& nrxx = rho_basis_.nrxx; // converting rho (extract it as a subfuntion in the future) diff --git a/source/source_lcao/module_rdmft/rdmft.cpp b/source/source_lcao/module_rdmft/rdmft.cpp index d62a3ee3716..d0ecde8b5f0 100644 --- a/source/source_lcao/module_rdmft/rdmft.cpp +++ b/source/source_lcao/module_rdmft/rdmft.cpp @@ -346,24 +346,10 @@ void RDMFT::cal_Energy(const int cal_type) this->pelec->cal_energies(2); Etotal = this->pelec->f_en.etot; - // if( GlobalC::exx_info.info_global.cal_exx ) - // { - // ModuleBase::matrix Exc_n_k(wg.nr, wg.nc, true); - // // because we have got wk_fun_occNum, we can use symbol=1 realize it - // occNum_Mul_wfcHwfc(wk_fun_occNum, wfcHwfc_XC, Exc_n_k, 1); - // E_RDMFT[2] = getEnergy(Exc_n_k); - // Parallel_Reduce::reduce_all(E_RDMFT[2]); - - // // test - // Etotal -= E_RDMFT[2]; - // } - } + } // // print results // std::cout << "\n\nfrom class RDMFT: \nXC_fun: " << XC_func_rdmft << std::endl; -// #ifdef __EXX -// if( GlobalC::exx_info.info_global.cal_exx ) std::cout << "alpha_power: " << alpha_power << std::endl; -// #endif // std::cout << std::fixed << std::setprecision(10) // << "******\nE(TV + Hartree + XC) by RDMFT: " << E_RDMFT[3] // << "\n\nE_TV_RDMFT: " << E_RDMFT[0] diff --git a/source/source_lcao/module_ri/Exx_LRI_interface.hpp b/source/source_lcao/module_ri/Exx_LRI_interface.hpp index 6a0c74aa091..d9155e6197d 100644 --- a/source/source_lcao/module_ri/Exx_LRI_interface.hpp +++ b/source/source_lcao/module_ri/Exx_LRI_interface.hpp @@ -215,7 +215,9 @@ void Exx_LRI_Interface::exx_hamilt2rho(elecstate::ElecState& elec, con Parallel_Common::bcast_double(this->exx_ptr->Eexx); this->exx_ptr->Eexx /= GlobalC::exx_info.info_global.hybrid_alpha; } - elec.set_exx(this->get_Eexx()); + bool cal_exx = GlobalC::exx_info.info_global.cal_exx; + double hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; + elec.set_exx(this->get_Eexx(), cal_exx, hybrid_alpha); } else { diff --git a/source/source_lcao/module_ri/RPA_LRI.hpp b/source/source_lcao/module_ri/RPA_LRI.hpp index 0323df970d0..31b6181850b 100644 --- a/source/source_lcao/module_ri/RPA_LRI.hpp +++ b/source/source_lcao/module_ri/RPA_LRI.hpp @@ -152,8 +152,6 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix GlobalC::exx_info.sync_from_global(); // reserve exx_ccp_rmesh_times to calculate full Coulomb this->ccp_rmesh_times_ewald = GlobalC::exx_info.info_ri.ccp_rmesh_times; - // rpa=1 set - // GlobalC::exx_info.info_ri.ccp_rmesh_times=rpa_ccp_rmesh_times // Using this->info.ccp_rmesh_times to calculate cut Coulomb this->Vs_period GlobalC::exx_info.info_ri.ccp_rmesh_times = PARAM.inp.rpa_ccp_rmesh_times; if (!exx_cut_coulomb) @@ -1249,13 +1247,10 @@ void RPA_LRI::out_coulomb_k(const UnitCell& ucell, // this->info.kmesh_times, this->info.ccp_rmesh_times ); // } -// // for( size_t T=0; T!=this->abfs.size(); ++T ) -// // GlobalC::exx_info.info_ri.abfs_Lmax = std::max( -// GlobalC::exx_info.info_ri.abfs_Lmax, static_cast(this->abfs[T].size())-1 -// ); - // } + + // template // void RPA_LRI::cal_rpa_ions() // { diff --git a/source/source_lcao/module_ri/conv_coulomb_pot_k.h b/source/source_lcao/module_ri/conv_coulomb_pot_k.h index 8793f573ca5..ca3ce4ce40f 100644 --- a/source/source_lcao/module_ri/conv_coulomb_pot_k.h +++ b/source/source_lcao/module_ri/conv_coulomb_pot_k.h @@ -48,6 +48,8 @@ namespace Conv_Coulomb_Pot_K const double rcut); } +using CoulombParam = std::map>>; + #include "conv_coulomb_pot_k.hpp" #endif \ No newline at end of file diff --git a/source/source_lcao/module_ri/exx_lip.h b/source/source_lcao/module_ri/exx_lip.h index ce4a27bfba3..d3c537b9200 100644 --- a/source/source_lcao/module_ri/exx_lip.h +++ b/source/source_lcao/module_ri/exx_lip.h @@ -110,9 +110,6 @@ class Exx_Lip void b_sum(const int iq, const int ib); void sum_all(const int ik); void exx_energy_cal(); - // void read_q_pack(const ModuleSymmetry::Symmetry& symm, - // const ModulePW::PW_Basis_K* wfc_basis, - // const Structure_Factor& sf); //2*pi*i const T two_pi_i = Treal(ModuleBase::TWO_PI) * T(0.0, 1.0); diff --git a/source/source_lcao/module_ri/exx_lip.hpp b/source/source_lcao/module_ri/exx_lip.hpp index 2dfd54f96af..fb84383825d 100644 --- a/source/source_lcao/module_ri/exx_lip.hpp +++ b/source/source_lcao/module_ri/exx_lip.hpp @@ -127,11 +127,6 @@ Exx_Lip::Exx_Lip(const Exx_Info_Lip& info_in, { this->q_pack = this->k_pack; } - // else if(PARAM.inp.init_chg=="file") - // { - // read_q_pack(symm, this->wfc_basis, sf); - // } - this->phi.resize(PARAM.globalv.nlocal); for (int iw = 0; iw < PARAM.globalv.nlocal; ++iw) { this->phi[iw].resize(this->rho_basis->nrxx); } @@ -538,88 +533,4 @@ void Exx_Lip::write_q_pack() const ModuleBase::timer::end("Exx_Lip", "write_q_pack"); } -/* -void Exx_Lip::read_q_pack(const ModuleSymmetry::Symmetry& symm, - const ModulePW::PW_Basis_K* this->wfc_basis, - const Structure_Factor& sf) -{ - const std::string exx_q_pack = "exx_q_pack/"; - this->q_pack = new k_package(); - this->q_pack->kv_ptr = new K_Vectors(); - const std::string exx_kpoint_card = PARAM.globalv.global_out_dir + exx_q_pack + PARAM.inp.kpoint_file; - this->q_pack->kv_ptr->set( symm, exx_kpoint_card, PARAM.inp.nspin, this->ucell_ptr->G, this->ucell_ptr->latvec, GlobalV::ofs_running ); - this->q_pack->wf_ptr = new wavefunc(); - this->q_pack->wf_ptr->allocate(this->q_pack->kv_ptr->get_nkstot(), - this->q_pack->kv_ptr->get_nks(), - this->q_pack->kv_ptr->ngk.data(), - this->wfc_basis->npwk_max); // mohan update 2021-02-25 - // this->q_pack->wf_ptr->init(this->q_pack->kv_ptr->get_nks(),this->q_pack->kv_ptr,this->ucell_ptr,old_pwptr,&ppcell,&GlobalC::ORB,&hm,&Pkpoints); - this->q_pack->wf_ptr->table_local.create(ucell.ntype, ucell.nmax_total, PARAM.globalv.nqx); - // this->q_pack->wf_ptr->table_local.create(this->q_pack->wf_ptr->this->ucell_ptr->ntype, this->q_pack->wf_ptr->this->ucell_ptr->nmax_total, PARAM.globalv.nqx); - #ifdef __LCAO - Wavefunc_in_pw::make_table_q(GlobalC::ORB.orbital_file, this->q_pack->wf_ptr->table_local); - // Wavefunc_in_pw::make_table_q(this->q_pack->wf_ptr->ORB_ptr->orbital_file, this->q_pack->wf_ptr->table_local, this->q_pack->wf_ptr); - for(int iq=0; iqq_pack->kv_ptr->get_nks(); ++iq) - { - Wavefunc_in_pw::produce_local_basis_in_pw(iq, - this->wfc_basis, - sf, - this->q_pack->wf_ptr->wanf2[iq], - this->q_pack->wf_ptr->table_local); - // Wavefunc_in_pw::produce_local_basis_in_pw(iq, this->q_pack->wf_ptr->wanf2[iq], this->q_pack->wf_ptr->table_local, - // this->q_pack->wf_ptr); - } - #endif - this->q_pack->wf_wg.create(this->q_pack->kv_ptr->get_nks(),PARAM.inp.nbands); - if(!GlobalV::RANK_IN_POOL) - { - std::stringstream ss_wf_wg; - ss_wf_wg << PARAM.globalv.global_out_dir << exx_q_pack << "wf_wg_" << GlobalV::MY_POOL; - std::ifstream ifs_wf_wg(ss_wf_wg.str().c_str()); - for( int iq = 0; iq < this->q_pack->kv_ptr->get_nks(); ++iq) - { - for( int ib=0; ib>this->q_pack->wf_wg(iq,ib); - } - } - ifs_wf_wg.close(); - } - #ifdef __MPI - MPI_Bcast( this->q_pack->wf_wg.c, this->q_pack->kv_ptr->get_nks()*PARAM.inp.nbands, MPI_DOUBLE, 0, POOL_WORLD); - #endif - this->q_pack->hvec_array = new ModuleBase::ComplexMatrix [this->q_pack->kv_ptr->get_nks()]; - for( int iq=0; iqq_pack->kv_ptr->get_nks(); ++iq) - { - this->q_pack->hvec_array[iq].create(PARAM.globalv.nlocal,PARAM.inp.nbands); - } - if(!GlobalV::RANK_IN_POOL) - { - std::stringstream ss_hvec; - ss_hvec << PARAM.globalv.global_out_dir << exx_q_pack << "hvec_" << GlobalV::MY_POOL; - std::ifstream ifs_hvec(ss_hvec.str().c_str()); - for( int iq=0; iqq_pack->kv_ptr->get_nks(); ++iq) - { - for( int iw=0; iwb; - ifs_hvec>>a>>this->b; - this->q_pack->hvec_array[iq](iw,ib) = {a,this->b}; - } - } - } - ifs_hvec.close(); - } - #ifdef __MPI - for( int iq=0; iqq_pack->kv_ptr->get_nks(); ++iq) - { - MPI_Bcast( this->q_pack->hvec_array[iq].c, PARAM.globalv.nlocal*PARAM.inp.nbands, MPI_DOUBLE_COMPLEX, 0, POOL_WORLD); - } - #endif - return; -} -*/ - #endif diff --git a/source/source_lcao/spar_exx.cpp b/source/source_lcao/spar_exx.cpp index 65efc194814..57715b25ccc 100644 --- a/source/source_lcao/spar_exx.cpp +++ b/source/source_lcao/spar_exx.cpp @@ -36,12 +36,13 @@ void cal_HR_exx( const int& current_spin, const double& sparse_threshold, const int (&nmp)[3], - const std::vector>, RI::Tensor>>>& Hexxs) + const std::vector>, RI::Tensor>>>& Hexxs, + const double hybrid_alpha) { ModuleBase::TITLE("sparse_format", "cal_HR_exx"); ModuleBase::timer::start("sparse_format", "cal_HR_exx"); - const Tdata frac = GlobalC::exx_info.info_global.hybrid_alpha; + const Tdata frac = hybrid_alpha; std::map> atoms_pos; for (int iat = 0; iat < ucell.nat; ++iat) @@ -157,7 +158,8 @@ template void cal_HR_exx( const int& current_spin, const double& sparse_thr, const int (&nmp)[3], - const std::vector>, RI::Tensor>>>& Hexxs); + const std::vector>, RI::Tensor>>>& Hexxs, + const double hybrid_alpha); template void cal_HR_exx>( const UnitCell& ucell, @@ -166,7 +168,8 @@ template void cal_HR_exx>( const int& current_spin, const double& sparse_thr, const int (&nmp)[3], - const std::vector>, RI::Tensor>>>>& Hexxs); + const std::vector>, RI::Tensor>>>>& Hexxs, + const double hybrid_alpha); } // namespace sparse_format diff --git a/source/source_lcao/spar_exx.h b/source/source_lcao/spar_exx.h index 4aca75aa5ae..ffbc29fc2b9 100644 --- a/source/source_lcao/spar_exx.h +++ b/source/source_lcao/spar_exx.h @@ -46,7 +46,8 @@ void cal_HR_exx( const int& current_spin, const double& sparse_thr, const int (&nmp)[3], - const std::vector>, RI::Tensor>>>& Hexxs); + const std::vector>, RI::Tensor>>>& Hexxs, + const double hybrid_alpha); // Explicit instantiations for double and complex types extern template void cal_HR_exx( @@ -56,7 +57,8 @@ extern template void cal_HR_exx( const int& current_spin, const double& sparse_thr, const int (&nmp)[3], - const std::vector>, RI::Tensor>>>& Hexxs); + const std::vector>, RI::Tensor>>>& Hexxs, + const double hybrid_alpha); extern template void cal_HR_exx>( const UnitCell& ucell, @@ -65,7 +67,8 @@ extern template void cal_HR_exx>( const int& current_spin, const double& sparse_thr, const int (&nmp)[3], - const std::vector>, RI::Tensor>>>>& Hexxs); + const std::vector>, RI::Tensor>>>>& Hexxs, + const double hybrid_alpha); } diff --git a/source/source_pw/module_pwdft/exx_helper.cpp b/source/source_pw/module_pwdft/exx_helper.cpp index fa80a7ff102..68ab9d8f9da 100644 --- a/source/source_pw/module_pwdft/exx_helper.cpp +++ b/source/source_pw/module_pwdft/exx_helper.cpp @@ -1,5 +1,5 @@ #include "exx_helper.h" -#include "source_io/module_parameter/parameter.h" // use PARAM +#include "source_io/module_parameter/parameter.h" #include "source_hamilt/module_xc/exx_info.h" // use GlobalC::exx_info #include "source_hamilt/module_xc/xc_functional.h" // use XC_Functional #include "source_pw/module_pwdft/hamilt_pw.h" // use HamiltPW @@ -61,13 +61,13 @@ bool Exx_Helper::iter_finish(void* p_elec, Charge* p_charge, void* ps bool& conv_esolver, int& iter) { /// Return if EXX is not enabled - if (!GlobalC::exx_info.info_global.cal_exx) + if (op_exx == nullptr) { return false; } /// Handle separate_loop mode - if (GlobalC::exx_info.info_global.separate_loop) + if (op_exx->separate_loop) { if (conv_esolver) { @@ -131,7 +131,7 @@ bool Exx_Helper::exx_after_converge(int &iter, bool ene_conv) { op_exx->first_iter = false; } - else if (!GlobalC::exx_info.info_global.separate_loop) + else if (!op_exx->separate_loop) { return true; } @@ -162,8 +162,10 @@ void Exx_Helper::set_psi(void* psi_) if (psi == nullptr) return; this->psi = psi; + if (op_exx == nullptr) + return; op_exx->set_psi(*psi); - if (PARAM.inp.exxace && GlobalC::exx_info.info_global.separate_loop) + if (PARAM.inp.exxace && op_exx->separate_loop) { op_exx->construct_ace(); } diff --git a/source/source_pw/module_pwdft/forces_cc.cpp b/source/source_pw/module_pwdft/forces_cc.cpp index ed962be3c1c..9b743c5700f 100644 --- a/source/source_pw/module_pwdft/forces_cc.cpp +++ b/source/source_pw/module_pwdft/forces_cc.cpp @@ -55,12 +55,18 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, ModuleBase::matrix v(PARAM.inp.nspin, rho_basis->nrxx); + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); +#ifdef __EXX + const double hse_omega = XC_Functional::get_hse_omega(); +#else + const double hse_omega = 0.0; +#endif if (XC_Functional::get_ked_flag()) { #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, - PARAM.inp.nspin); + PARAM.inp.nspin, hybrid_alpha, hse_omega); // etxc = std::get<0>(etxc_vtxc_v); // vtxc = std::get<1>(etxc_vtxc_v); @@ -75,7 +81,9 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, 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); + PARAM.globalv.domag_z, + hybrid_alpha, + hse_omega); // etxc = std::get<0>(etxc_vtxc_v); // vtxc = std::get<1>(etxc_vtxc_v); diff --git a/source/source_pw/module_pwdft/hamilt_pw.cpp b/source/source_pw/module_pwdft/hamilt_pw.cpp index 47fae9bcb96..7c06f971690 100644 --- a/source/source_pw/module_pwdft/hamilt_pw.cpp +++ b/source/source_pw/module_pwdft/hamilt_pw.cpp @@ -132,7 +132,10 @@ HamiltPW::HamiltPW(elecstate::Potential* pot_in, } if (GlobalC::exx_info.info_global.cal_exx) { - auto exx = new OperatorEXXPW(isk, wfc_basis, pot_in->get_rho_basis(), pkv, ucell); + bool separate_loop = GlobalC::exx_info.info_global.separate_loop; + double hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; + auto coulomb_param = GlobalC::exx_info.info_global.coulomb_param; + auto exx = new OperatorEXXPW(isk, wfc_basis, pot_in->get_rho_basis(), pkv, ucell, separate_loop, hybrid_alpha, coulomb_param); if (this->ops == nullptr) { this->ops = exx; diff --git a/source/source_pw/module_pwdft/kernels/cuda/force_op.cu b/source/source_pw/module_pwdft/kernels/cuda/force_op.cu index 48b0f10b5ed..cdcd477e82e 100644 --- a/source/source_pw/module_pwdft/kernels/cuda/force_op.cu +++ b/source/source_pw/module_pwdft/kernels/cuda/force_op.cu @@ -98,7 +98,6 @@ __global__ void cal_force_nl( { ps_qq = - ekb_now * qq_nt[it * deeq_3 * deeq_4 + ip * deeq_4 + ip]; } - // FPTYPE ps = GlobalC::ppcell.deeq[spin, iat, ip, ip]; FPTYPE ps = deeq[((spin * deeq_2 + iat) * deeq_3 + ip) * deeq_4 + ip] + ps_qq; const int inkb = sum + ip; //out<<"\n ps = "< #include @@ -37,8 +36,13 @@ OperatorEXXPW::OperatorEXXPW(const int* isk_in, const ModulePW::PW_Basis_K* wfcpw_in, const ModulePW::PW_Basis* rhopw_in, K_Vectors *kv_in, - const UnitCell *ucell) - : isk(isk_in), wfcpw(wfcpw_in), rhopw(rhopw_in), kv(kv_in), ucell(ucell) + const UnitCell *ucell, + const bool separate_loop_in, + const Real hybrid_alpha_in, + const CoulombParam& coulomb_param_in) + : isk(isk_in), wfcpw(wfcpw_in), rhopw(rhopw_in), kv(kv_in), ucell(ucell), + separate_loop(separate_loop_in), hybrid_alpha(hybrid_alpha_in), + coulomb_param(coulomb_param_in) { if (GlobalV::KPAR != 1 && PARAM.inp.exxace == false) { @@ -98,7 +102,7 @@ OperatorEXXPW::OperatorEXXPW(const int* isk_in, rhopw_dev->setuptransform(); rhopw_dev->collect_local_pw(); - auto param_fock = GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock]; + auto param_fock = this->coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock]; for (auto param: param_fock) { fock_div.push_back(exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Fock, @@ -110,7 +114,7 @@ OperatorEXXPW::OperatorEXXPW(const int* isk_in, gamma_extrapolation, ucell->omega)); } - auto param_erfc = GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc]; + auto param_erfc = this->coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc]; for (auto param: param_erfc) { erfc_div.push_back(exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc, @@ -186,7 +190,7 @@ void OperatorEXXPW::act(const int nbands, setmem_complex_op()(tmhpsi, 0, nbasis*nbands/npol); } - if (PARAM.inp.exxace && GlobalC::exx_info.info_global.separate_loop) + if (PARAM.inp.exxace && this->separate_loop) { act_op_ace(nbands, nbasis, npol, tmpsi_in, tmhpsi, ngk_ik, is_first_node); } @@ -234,10 +238,9 @@ void OperatorEXXPW::act_op(const int nbands, Real nqs = q_points.size(); for (int iq: q_points) { - get_exx_potential(kv, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, this->ik, iq % nk); + get_exx_potential(kv, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, this->ik, iq % nk, false, this->coulomb_param); for (int m_iband = 0; m_iband < psi.get_nbands(); m_iband++) { - // double wg_mqb_real = GlobalC::exx_helper.wg(iq, m_iband); double wg_mqb_real = (*wg)(this->ik, m_iband); T wg_mqb = wg_mqb_real; if (wg_mqb_real < 1e-12) @@ -283,8 +286,7 @@ void OperatorEXXPW::act_op(const int nbands, } // end of iq T* h_psi_nk = tmhpsi + n_iband * nbasis; - Real hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; - wfcpw->real_to_recip(ctx, h_psi_real, h_psi_nk, this->ik, true, hybrid_alpha); + wfcpw->real_to_recip(ctx, h_psi_real, h_psi_nk, this->ik, true, this->hybrid_alpha); setmem_complex_op()(h_psi_real, 0, rhopw_dev->nrxx); } @@ -320,7 +322,7 @@ void OperatorEXXPW::act_op_kpar(const int nbands, for (int iq = 0; iq < nqs; iq++) { // for \psi_nk, get the pw of iq and band m - get_exx_potential(kv, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, this->ik, iq); + get_exx_potential(kv, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, this->ik, iq, false, this->coulomb_param); // decide which pool does the iq belong to int iq_pool = kv->para_k.whichpool[iq]; @@ -384,8 +386,7 @@ void OperatorEXXPW::act_op_kpar(const int nbands, Real tmp_scalar = wg_mqb / wk_ik / nqs; // wk_ik works for now, but wrong for symmetry. T* h_psi_nk = tmhpsi + n_iband * nbasis; - Real hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; - wfcpw->real_to_recip(ctx, density_real, h_psi_nk, this->ik, true, hybrid_alpha * tmp_scalar); + wfcpw->real_to_recip(ctx, density_real, h_psi_nk, this->ik, true, this->hybrid_alpha * tmp_scalar); } // end of m_iband @@ -496,7 +497,7 @@ OperatorEXXPW::OperatorEXXPW(const OperatorEXXPW *op template double OperatorEXXPW::cal_exx_energy(psi::Psi *psi_) const { - if (PARAM.inp.exxace && GlobalC::exx_info.info_global.separate_loop) + if (PARAM.inp.exxace && this->separate_loop) { return cal_exx_energy_ace(psi_); } @@ -534,7 +535,6 @@ double OperatorEXXPW::cal_exx_energy_op(psi::Psi *ppsi_) c setmem_complex_op()(density_real, 0, rhopw_dev->nrxx); setmem_complex_op()(density_recip, 0, rhopw_dev->npw); - // double wg_ikb_real = GlobalC::exx_helper.wg(this->ik, n_iband); double wg_ikb_real = (*wg)(ik, n_iband); T wg_ikb = wg_ikb_real; if (wg_ikb_real < 1e-12) @@ -580,10 +580,9 @@ double OperatorEXXPW::cal_exx_energy_op(psi::Psi *ppsi_) c for (int iq: q_points) { int nk = wfcpw->nks / nk_fac; - get_exx_potential(kv, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, ik, iq % nk); + get_exx_potential(kv, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, ik, iq % nk, false, this->coulomb_param); for (int m_iband = 0; m_iband < psi.get_nbands(); m_iband++) { - // double wg_f = GlobalC::exx_helper.wg(iq, m_iband); double wg_iqb_real = (*wg)(iq, m_iband); T wg_iqb = wg_iqb_real; if (wg_iqb_real < 1e-12) diff --git a/source/source_pw/module_pwdft/op_pw_exx.h b/source/source_pw/module_pwdft/op_pw_exx.h index 0329de84779..f3a108d5957 100644 --- a/source/source_pw/module_pwdft/op_pw_exx.h +++ b/source/source_pw/module_pwdft/op_pw_exx.h @@ -30,7 +30,10 @@ class OperatorEXXPW : public OperatorPW const ModulePW::PW_Basis_K* wfcpw_in, const ModulePW::PW_Basis* rhopw_in, K_Vectors* kv_in, - const UnitCell* ucell); + const UnitCell* ucell, + const bool separate_loop_in, + const Real hybrid_alpha_in, + const CoulombParam& coulomb_param_in); template explicit OperatorEXXPW(const OperatorEXXPW *op_exx); @@ -54,6 +57,9 @@ class OperatorEXXPW : public OperatorPW void construct_ace() const; bool first_iter = true; + bool separate_loop = false; + Real hybrid_alpha = 0.0; + CoulombParam coulomb_param; static std::vector fock_div, erfc_div; @@ -179,7 +185,8 @@ void get_exx_potential(const K_Vectors* kv, double ucell_omega, int ik, int iq, - bool is_stress = false); + bool is_stress, + const CoulombParam& coulomb_param_in); template void get_exx_stress_potential(const K_Vectors* kv, @@ -190,7 +197,8 @@ void get_exx_stress_potential(const K_Vectors* kv, bool gamma_extrapolation, double ucell_omega, int ik, - int iq); + int iq, + const CoulombParam& coulomb_param_in); double exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type coulomb_type, double erfc_omega, diff --git a/source/source_pw/module_pwdft/op_pw_exx_ace.cpp b/source/source_pw/module_pwdft/op_pw_exx_ace.cpp index 8685b355f08..071c4e4440d 100644 --- a/source/source_pw/module_pwdft/op_pw_exx_ace.cpp +++ b/source/source_pw/module_pwdft/op_pw_exx_ace.cpp @@ -159,7 +159,7 @@ void OperatorEXXPW::construct_ace() const iq = iq0 + ispin * nk; // iq in the same spin channel // for \psi_nk, get the pw of iq and band m - get_exx_potential(kv, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, ik, iq); + get_exx_potential(kv, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, ik, iq, false, this->coulomb_param); // decide which pool does the iq belong to int iq_pool = kv->para_k.whichpool[iq0]; @@ -295,7 +295,7 @@ double OperatorEXXPW::cal_exx_energy_ace(psi::Psi* ppsi_) psi::Psi psi_ = *ppsi_; int* ik_ = const_cast(&this->ik); int ik_save = this->ik; - Real hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; + Real hybrid_alpha = this->hybrid_alpha; for (int i = 0; i < wfcpw->nks; i++) { setmem_complex_op()(h_psi_ace, 0, psi_.get_nbands() * psi_.get_nbasis()); diff --git a/source/source_pw/module_pwdft/op_pw_exx_pot.cpp b/source/source_pw/module_pwdft/op_pw_exx_pot.cpp index 9f31bf642e4..672b7d48b78 100644 --- a/source/source_pw/module_pwdft/op_pw_exx_pot.cpp +++ b/source/source_pw/module_pwdft/op_pw_exx_pot.cpp @@ -1,7 +1,6 @@ #include "op_pw_exx.h" #include "source_base/parallel_reduce.h" #include "source_io/module_parameter/parameter.h" -#include "source_hamilt/module_xc/exx_info.h" // use GlobalC::exx_info namespace hamilt { @@ -15,7 +14,8 @@ void get_exx_potential(const K_Vectors* kv, double ucell_omega, int ik, int iq, - bool is_stress) + bool is_stress, + const CoulombParam& coulomb_param_in) { using setmem_real_cpu_op = base_device::memory::set_memory_op; using syncmem_real_c2d_op = base_device::memory::synchronize_memory_op; @@ -46,148 +46,154 @@ void get_exx_potential(const K_Vectors* kv, } // calculate Fock pot - auto param_fock = GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock]; - for (int i = 0; i < param_fock.size(); i++) + auto it_fock = coulomb_param_in.find(Conv_Coulomb_Pot_K::Coulomb_Type::Fock); + if (it_fock != coulomb_param_in.end()) { - auto param = param_fock[i]; - double exx_div = OperatorEXXPW, Device>::fock_div[i]; - double alpha = std::stod(param["alpha"]); - const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; - const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; - const ModuleBase::Vector3 q_c = qvec_c[iq]; - const ModuleBase::Vector3 q_d = qvec_d[iq]; + for (int i = 0; i < it_fock->second.size(); i++) + { + auto param = it_fock->second[i]; + double exx_div = OperatorEXXPW, Device>::fock_div[i]; + double alpha = std::stod(param["alpha"]); + const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; + const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; + const ModuleBase::Vector3 q_c = qvec_c[iq]; + const ModuleBase::Vector3 q_d = qvec_d[iq]; #ifdef _OPENMP #pragma omp parallel for schedule(static) #endif - for (int ig = 0; ig < rhopw_dev->npw; ig++) - { - const ModuleBase::Vector3 g_d = rhopw_dev->gdirect[ig]; - const ModuleBase::Vector3 kqg_d = k_d - q_d + g_d; - // For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114) - // 7/8 of the points in the grid are "activated" and 1/8 are disabled. - // grid_factor is designed for the 7/8 of the grid to function like all of the points - Real grid_factor = 1; - double extrapolate_grid = 8.0 / 7.0; - if (gamma_extrapolation) + for (int ig = 0; ig < rhopw_dev->npw; ig++) { - // if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3) - auto isint = [](double x) { - double epsilon = 1e-6; // this follows the isint judgement in q-e - return std::abs(x - std::round(x)) < epsilon; - }; - if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)) + const ModuleBase::Vector3 g_d = rhopw_dev->gdirect[ig]; + const ModuleBase::Vector3 kqg_d = k_d - q_d + g_d; + // For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114) + // 7/8 of the points in the grid are "activated" and 1/8 are disabled. + // grid_factor is designed for the 7/8 of the grid to function like all of the points + Real grid_factor = 1; + double extrapolate_grid = 8.0 / 7.0; + if (gamma_extrapolation) { - grid_factor = 0; + // if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3) + auto isint = [](double x) { + double epsilon = 1e-6; // this follows the isint judgement in q-e + return std::abs(x - std::round(x)) < epsilon; + }; + if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)) + { + grid_factor = 0; + } + else + { + grid_factor = extrapolate_grid; + } + } + + const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + const int nk = nks / nk_fac; + + Real gg = (k_c - q_c + rhopw_dev->gcar[ig]).norm2() * tpiba2; + // if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2 + if (gg >= 1e-8) + { + Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg; + pot_cpu[ig] += fac * grid_factor * alpha; } + // } else { - grid_factor = extrapolate_grid; + pot_cpu[ig] += exx_div * alpha; } + // assert(is_finite(density_recip[ig])); } - - const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; - const int nk = nks / nk_fac; - - Real gg = (k_c - q_c + rhopw_dev->gcar[ig]).norm2() * tpiba2; - // if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2 - if (gg >= 1e-8) - { - Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg; - pot_cpu[ig] += fac * grid_factor * alpha; - } - // } - else - { - pot_cpu[ig] += exx_div * alpha; - } - // assert(is_finite(density_recip[ig])); } } // calculate erfc pot - auto param_erfc = GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc]; - for (int i = 0; i < param_erfc.size(); i++) + auto it_erfc = coulomb_param_in.find(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc); + if (it_erfc != coulomb_param_in.end()) { - auto param = param_erfc[i]; - double erfc_omega = std::stod(param["omega"]); - double erfc_omega2 = erfc_omega * erfc_omega; - double alpha = std::stod(param["alpha"]); - // double exx_div = OperatorEXXPW, Device>::erfc_div[i]; - double exx_div = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc, - erfc_omega, - kv, - wfcpw, - rhopw_dev, - tpiba, - gamma_extrapolation, - ucell_omega); - const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; - const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; - const ModuleBase::Vector3 q_c = qvec_c[iq]; - const ModuleBase::Vector3 q_d = qvec_d[iq]; + for (int i = 0; i < it_erfc->second.size(); i++) + { + auto param = it_erfc->second[i]; + double erfc_omega = std::stod(param["omega"]); + double erfc_omega2 = erfc_omega * erfc_omega; + double alpha = std::stod(param["alpha"]); + // double exx_div = OperatorEXXPW, Device>::erfc_div[i]; + double exx_div = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc, + erfc_omega, + kv, + wfcpw, + rhopw_dev, + tpiba, + gamma_extrapolation, + ucell_omega); + const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; + const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; + const ModuleBase::Vector3 q_c = qvec_c[iq]; + const ModuleBase::Vector3 q_d = qvec_d[iq]; #ifdef _OPENMP #pragma omp parallel for schedule(static) #endif - for (int ig = 0; ig < rhopw_dev->npw; ig++) - { - const ModuleBase::Vector3 g_d = rhopw_dev->gdirect[ig]; - const ModuleBase::Vector3 kqg_d = k_d - q_d + g_d; - // For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114) - // 7/8 of the points in the grid are "activated" and 1/8 are disabled. - // grid_factor is designed for the 7/8 of the grid to function like all of the points - Real grid_factor = 1; - double extrapolate_grid = 8.0 / 7.0; - if (gamma_extrapolation) + for (int ig = 0; ig < rhopw_dev->npw; ig++) { - // if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3) - auto isint = [](double x) { - double epsilon = 1e-6; // this follows the isint judgement in q-e - return std::abs(x - std::round(x)) < epsilon; - }; - if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)) - { - grid_factor = 0; - } - else + const ModuleBase::Vector3 g_d = rhopw_dev->gdirect[ig]; + const ModuleBase::Vector3 kqg_d = k_d - q_d + g_d; + // For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114) + // 7/8 of the points in the grid are "activated" and 1/8 are disabled. + // grid_factor is designed for the 7/8 of the grid to function like all of the points + Real grid_factor = 1; + double extrapolate_grid = 8.0 / 7.0; + if (gamma_extrapolation) { - grid_factor = extrapolate_grid; + // if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3) + auto isint = [](double x) { + double epsilon = 1e-6; // this follows the isint judgement in q-e + return std::abs(x - std::round(x)) < epsilon; + }; + if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)) + { + grid_factor = 0; + } + else + { + grid_factor = extrapolate_grid; + } } - } - const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; - const int nk = nks / nk_fac; - // const int ig_kq = ik * nks * npw + iq * npw + ig; - - Real gg = (k_c - q_c + rhopw_dev->gcar[ig]).norm2() * tpiba2; - // if (ig == 0 && GlobalV::MY_RANK==1) - // { - // printf("k-q+G: %f %f %f\n", (k_c - q_c + rhopw_dev->gcar[ig])[0], (k_c - q_c + rhopw_dev->gcar[ig])[1], (k_c - q_c + rhopw_dev->gcar[ig])[2]); - // } - // if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2 - if (gg >= 1e-8) - { - Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg; - pot_cpu[ig] += fac * (1.0 - std::exp(-gg / 4.0 / erfc_omega2)) * grid_factor * alpha; - } - // } - else - { - // if (PARAM.inp.dft_functional == "hse") - if (!gamma_extrapolation) + const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + const int nk = nks / nk_fac; + // const int ig_kq = ik * nks * npw + iq * npw + ig; + + Real gg = (k_c - q_c + rhopw_dev->gcar[ig]).norm2() * tpiba2; + // if (ig == 0 && GlobalV::MY_RANK==1) + // { + // printf("k-q+G: %f %f %f\n", (k_c - q_c + rhopw_dev->gcar[ig])[0], (k_c - q_c + rhopw_dev->gcar[ig])[1], (k_c - q_c + rhopw_dev->gcar[ig])[2]); + // } + // if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2 + if (gg >= 1e-8) { - if (is_stress) - pot_cpu[ig] += (- ModuleBase::PI * ModuleBase::e2 / erfc_omega2) * alpha; - else - pot_cpu[ig] += (exx_div - ModuleBase::PI * ModuleBase::e2 / erfc_omega2) * alpha; + Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg; + pot_cpu[ig] += fac * (1.0 - std::exp(-gg / 4.0 / erfc_omega2)) * grid_factor * alpha; } + // } else { - pot_cpu[ig] += exx_div * alpha; + // if (PARAM.inp.dft_functional == "hse") + if (!gamma_extrapolation) + { + if (is_stress) + pot_cpu[ig] += (- ModuleBase::PI * ModuleBase::e2 / erfc_omega2) * alpha; + else + pot_cpu[ig] += (exx_div - ModuleBase::PI * ModuleBase::e2 / erfc_omega2) * alpha; + } + else + { + pot_cpu[ig] += exx_div * alpha; + } } + // assert(is_finite(density_recip[ig])); } - // assert(is_finite(density_recip[ig])); } } @@ -221,7 +227,8 @@ void get_exx_stress_potential(const K_Vectors* kv, bool gamma_extrapolation, double ucell_omega, int ik, - int iq) + int iq, + const CoulombParam& coulomb_param_in) { using setmem_real_cpu_op = base_device::memory::set_memory_op; using syncmem_real_c2d_op = base_device::memory::synchronize_memory_op; @@ -238,139 +245,145 @@ void get_exx_stress_potential(const K_Vectors* kv, setmem_real_cpu_op()(pot_cpu, 0, npw); // calculate Fock pot - auto param_fock = GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock]; - for (auto param: param_fock) + auto it_fock = coulomb_param_in.find(Conv_Coulomb_Pot_K::Coulomb_Type::Fock); + if (it_fock != coulomb_param_in.end()) { - // double exx_div = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Fock, - // 0.0, - // kv, - // wfcpw, - // rhopw_dev, - // tpiba, - // gamma_extrapolation, - // ucell_omega); - double alpha = std::stod(param["alpha"]); - - const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; - const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; - const ModuleBase::Vector3 q_c = wfcpw->kvec_c[iq]; - const ModuleBase::Vector3 q_d = wfcpw->kvec_d[iq]; + for (auto param: it_fock->second) + { + // double exx_div = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Fock, + // 0.0, + // kv, + // wfcpw, + // rhopw_dev, + // tpiba, + // gamma_extrapolation, + // ucell_omega); + double alpha = std::stod(param["alpha"]); + + const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; + const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; + const ModuleBase::Vector3 q_c = wfcpw->kvec_c[iq]; + const ModuleBase::Vector3 q_d = wfcpw->kvec_d[iq]; #ifdef _OPENMP #pragma omp parallel for schedule(static) #endif - for (int ig = 0; ig < rhopw_dev->npw; ig++) - { - const ModuleBase::Vector3 g_d = rhopw_dev->gdirect[ig]; - const ModuleBase::Vector3 kqg_d = k_d - q_d + g_d; - // For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114) - // 7/8 of the points in the grid are "activated" and 1/8 are disabled. - // grid_factor is designed for the 7/8 of the grid to function like all of the points - Real grid_factor = 1; - double extrapolate_grid = 8.0 / 7.0; - if (gamma_extrapolation) + for (int ig = 0; ig < rhopw_dev->npw; ig++) { - // if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3) - auto isint = [](double x) { - double epsilon = 1e-6; // this follows the isint judgement in q-e - return std::abs(x - std::round(x)) < epsilon; - }; - if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)) - { - grid_factor = 0; - } - else + const ModuleBase::Vector3 g_d = rhopw_dev->gdirect[ig]; + const ModuleBase::Vector3 kqg_d = k_d - q_d + g_d; + // For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114) + // 7/8 of the points in the grid are "activated" and 1/8 are disabled. + // grid_factor is designed for the 7/8 of the grid to function like all of the points + Real grid_factor = 1; + double extrapolate_grid = 8.0 / 7.0; + if (gamma_extrapolation) { - grid_factor = extrapolate_grid; + // if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3) + auto isint = [](double x) { + double epsilon = 1e-6; // this follows the isint judgement in q-e + return std::abs(x - std::round(x)) < epsilon; + }; + if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)) + { + grid_factor = 0; + } + else + { + grid_factor = extrapolate_grid; + } } - } - const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; - const int nk = nks / nk_fac; - // const int ig_kq = ik * nks * npw + iq * npw + ig; + const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + const int nk = nks / nk_fac; + // const int ig_kq = ik * nks * npw + iq * npw + ig; - Real gg = (k_c - q_c + rhopw_dev->gcar[ig]).norm2() * tpiba2; - // if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2 - if (gg >= 1e-8) - { - Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg; - pot_cpu[ig] += 1.0 / gg * grid_factor * alpha; + Real gg = (k_c - q_c + rhopw_dev->gcar[ig]).norm2() * tpiba2; + // if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2 + if (gg >= 1e-8) + { + Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg; + pot_cpu[ig] += 1.0 / gg * grid_factor * alpha; + } } } } // calculate erfc pot - auto param_erfc = GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc]; - for (auto param: param_erfc) + auto it_erfc = coulomb_param_in.find(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc); + if (it_erfc != coulomb_param_in.end()) { - double erfc_omega = std::stod(param["omega"]); - double erfc_omega2 = erfc_omega * erfc_omega; - double alpha = std::stod(param["alpha"]); - // double exx_div = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc, - // erfc_omega, - // kv, - // wfcpw, - // rhopw_dev, - // tpiba, - // gamma_extrapolation, - // ucell_omega); - - const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; - const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; - const ModuleBase::Vector3 q_c = wfcpw->kvec_c[iq]; - const ModuleBase::Vector3 q_d = wfcpw->kvec_d[iq]; + for (auto param: it_erfc->second) + { + double erfc_omega = std::stod(param["omega"]); + double erfc_omega2 = erfc_omega * erfc_omega; + double alpha = std::stod(param["alpha"]); + // double exx_div = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc, + // erfc_omega, + // kv, + // wfcpw, + // rhopw_dev, + // tpiba, + // gamma_extrapolation, + // ucell_omega); + + const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; + const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; + const ModuleBase::Vector3 q_c = wfcpw->kvec_c[iq]; + const ModuleBase::Vector3 q_d = wfcpw->kvec_d[iq]; #ifdef _OPENMP #pragma omp parallel for schedule(static) #endif - for (int ig = 0; ig < rhopw_dev->npw; ig++) - { - const ModuleBase::Vector3 g_d = rhopw_dev->gdirect[ig]; - const ModuleBase::Vector3 kqg_d = k_d - q_d + g_d; - // For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114) - // 7/8 of the points in the grid are "activated" and 1/8 are disabled. - // grid_factor is designed for the 7/8 of the grid to function like all of the points - Real grid_factor = 1; - double extrapolate_grid = 8.0 / 7.0; - if (gamma_extrapolation) + for (int ig = 0; ig < rhopw_dev->npw; ig++) { - // if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3) - auto isint = [](double x) { - double epsilon = 1e-6; // this follows the isint judgement in q-e - return std::abs(x - std::round(x)) < epsilon; - }; - if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)) + const ModuleBase::Vector3 g_d = rhopw_dev->gdirect[ig]; + const ModuleBase::Vector3 kqg_d = k_d - q_d + g_d; + // For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114) + // 7/8 of the points in the grid are "activated" and 1/8 are disabled. + // grid_factor is designed for the 7/8 of the grid to function like all of the points + Real grid_factor = 1; + double extrapolate_grid = 8.0 / 7.0; + if (gamma_extrapolation) { - grid_factor = 0; - } - else - { - grid_factor = extrapolate_grid; + // if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3) + auto isint = [](double x) { + double epsilon = 1e-6; // this follows the isint judgement in q-e + return std::abs(x - std::round(x)) < epsilon; + }; + if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)) + { + grid_factor = 0; + } + else + { + grid_factor = extrapolate_grid; + } } - } - const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; - const int nk = nks / nk_fac; - // const int ig_kq = ik * nks * npw + iq * npw + ig; + const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + const int nk = nks / nk_fac; + // const int ig_kq = ik * nks * npw + iq * npw + ig; - Real gg = (k_c - q_c + rhopw_dev->gcar[ig]).norm2() * tpiba2; - // if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2 - if (gg >= 1e-8) - { - Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg; - pot_cpu[ig] += (1.0 - (1.0 + gg / 4.0 / erfc_omega2) * std::exp(-gg / 4.0 / erfc_omega2)) - / (1.0 - std::exp(-gg / 4.0 / erfc_omega2)) / gg * grid_factor * alpha; - } - // } - else - { - // if (PARAM.inp.dft_functional == "hse") - if (!gamma_extrapolation) + Real gg = (k_c - q_c + rhopw_dev->gcar[ig]).norm2() * tpiba2; + // if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2 + if (gg >= 1e-8) { - pot_cpu[ig] += 1.0 / 4.0 / erfc_omega2 * alpha; + Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg; + pot_cpu[ig] += (1.0 - (1.0 + gg / 4.0 / erfc_omega2) * std::exp(-gg / 4.0 / erfc_omega2)) + / (1.0 - std::exp(-gg / 4.0 / erfc_omega2)) / gg * grid_factor * alpha; } + // } + else + { + // if (PARAM.inp.dft_functional == "hse") + if (!gamma_extrapolation) + { + pot_cpu[ig] += 1.0 / 4.0 / erfc_omega2 * alpha; + } + } + // assert(is_finite(density_recip[ig])); } - // assert(is_finite(density_recip[ig])); } } @@ -529,7 +542,8 @@ template void get_exx_potential(const K_Vectors* double, int, int, - bool); + bool, + const CoulombParam&); template void get_exx_potential(const K_Vectors*, const ModulePW::PW_Basis_K*, ModulePW::PW_Basis*, @@ -539,7 +553,8 @@ template void get_exx_potential(const K_Vectors double, int, int, - bool); + bool, + const CoulombParam&); template void get_exx_stress_potential(const K_Vectors*, const ModulePW::PW_Basis_K*, ModulePW::PW_Basis*, @@ -548,7 +563,8 @@ template void get_exx_stress_potential(const K_V bool, double, int, - int); + int, + const CoulombParam&); template void get_exx_stress_potential(const K_Vectors*, const ModulePW::PW_Basis_K*, ModulePW::PW_Basis*, @@ -557,7 +573,8 @@ template void get_exx_stress_potential(const K_ bool, double, int, - int); + int, + const CoulombParam&); #if ((defined __CUDA) || (defined __ROCM)) template class OperatorEXXPW, base_device::DEVICE_GPU>; template class OperatorEXXPW, base_device::DEVICE_GPU>; @@ -570,7 +587,8 @@ template void get_exx_potential(const K_Vectors* double, int, int, - bool); + bool, + const CoulombParam&); template void get_exx_potential(const K_Vectors*, const ModulePW::PW_Basis_K*, ModulePW::PW_Basis*, @@ -580,7 +598,8 @@ template void get_exx_potential(const K_Vectors double, int, int, - bool); + bool, + const CoulombParam&); template void get_exx_stress_potential(const K_Vectors*, const ModulePW::PW_Basis_K*, ModulePW::PW_Basis*, @@ -589,7 +608,8 @@ template void get_exx_stress_potential(const K_V bool, double, int, - int); + int, + const CoulombParam&); template void get_exx_stress_potential(const K_Vectors*, const ModulePW::PW_Basis_K*, ModulePW::PW_Basis*, @@ -598,6 +618,7 @@ template void get_exx_stress_potential(const K_ bool, double, int, - int); + int, + const CoulombParam&); #endif } // namespace hamilt diff --git a/source/source_pw/module_pwdft/stress_cc.cpp b/source/source_pw/module_pwdft/stress_cc.cpp index 3b0f291e889..bdba3fcae22 100644 --- a/source/source_pw/module_pwdft/stress_cc.cpp +++ b/source/source_pw/module_pwdft/stress_cc.cpp @@ -51,12 +51,18 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, //recalculate the exchange-correlation potential ModuleBase::matrix vxc; + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); +#ifdef __EXX + const double hse_omega = XC_Functional::get_hse_omega(); +#else + const double hse_omega = 0.0; +#endif if (XC_Functional::get_ked_flag()) { #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, - PARAM.inp.nspin); + PARAM.inp.nspin, hybrid_alpha, hse_omega); // etxc = std::get<0>(etxc_vtxc_v); // vtxc = std::get<1>(etxc_vtxc_v); @@ -71,7 +77,9 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell, PARAM.inp.nspin, PARAM.globalv.domag, - PARAM.globalv.domag_z); + PARAM.globalv.domag_z, + hybrid_alpha, + hse_omega); // 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_exx.cpp b/source/source_pw/module_pwdft/stress_exx.cpp index 29900f92f67..9e7e99e2d01 100644 --- a/source/source_pw/module_pwdft/stress_exx.cpp +++ b/source/source_pw/module_pwdft/stress_exx.cpp @@ -1,4 +1,4 @@ -#include "source_hamilt/module_xc/exx_info.h" + #include "op_pw_exx.h" #include "source_base/parallel_common.h" #include "source_base/parallel_reduce.h" @@ -10,7 +10,9 @@ void Stress_PW::stress_exx(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rhopw, ModulePW::PW_Basis_K* wfcpw, const K_Vectors *p_kv, - const psi::Psi , Device>* d_psi_in, const UnitCell& ucell) + const psi::Psi , Device>* d_psi_in, const UnitCell& ucell, + const double hybrid_alpha, + const CoulombParam& coulomb_param) { bool gamma_extrapolation = PARAM.inp.exx_gamma_extrapolation; bool is_mp = p_kv->get_is_mp(); @@ -73,8 +75,8 @@ void Stress_PW::stress_exx(ModuleBase::matrix& sigma, for (int iq = 0; iq < nqs; iq++) { - hamilt::get_exx_potential(p_kv, wfcpw, rhopw, pot, tpiba, gamma_extrapolation, omega, ik, iq, true); - hamilt::get_exx_stress_potential(p_kv, wfcpw, rhopw, pot_stress, tpiba, gamma_extrapolation, omega, ik, iq); + hamilt::get_exx_potential(p_kv, wfcpw, rhopw, pot, tpiba, gamma_extrapolation, omega, ik, iq, true, coulomb_param); + hamilt::get_exx_stress_potential(p_kv, wfcpw, rhopw, pot_stress, tpiba, gamma_extrapolation, omega, ik, iq, coulomb_param); for (int mband = 0; mband < d_psi_in->get_nbands(); mband++) { // psi_mq in real space @@ -118,8 +120,7 @@ void Stress_PW::stress_exx(ModuleBase::matrix& sigma, } - // 0.5 in the following line is caused by 2x in the pot - sigma(alpha, beta) -= GlobalC::exx_info.info_global.hybrid_alpha + sigma(alpha, beta) -= hybrid_alpha * 0.25 * sigma_ab_loc * wg(ik, nband) * wg(iq, mband) / nqs / p_kv->wk[ik]; } diff --git a/source/source_pw/module_pwdft/stress_gga.cpp b/source/source_pw/module_pwdft/stress_gga.cpp index 198e4cc843a..bb175514803 100644 --- a/source/source_pw/module_pwdft/stress_gga.cpp +++ b/source/source_pw/module_pwdft/stress_gga.cpp @@ -1,68 +1,67 @@ #include "stress_func.h" #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 void Stress_Func::stress_gga(const UnitCell& ucell, - ModuleBase::matrix& sigma, + ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const Charge* const chr) { ModuleBase::TITLE("Stress","stress_gga"); - ModuleBase::timer::start("Stress","stress_gga"); - - int func_type = XC_Functional::get_func_type(); - if (func_type == 0 || func_type == 1) - { - ModuleBase::timer::end("Stress","stress_gga"); - return; - } + ModuleBase::timer::start("Stress","stress_gga"); - FPTYPE sigma_gradcorr[3][3]; - std::vector stress_gga; - FPTYPE dum1=0.0; + int func_type = XC_Functional::get_func_type(); + if (func_type == 0 || func_type == 1) + { + ModuleBase::timer::end("Stress","stress_gga"); + return; + } + + FPTYPE sigma_gradcorr[3][3]; + std::vector stress_gga; + 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. + ModuleBase::matrix dum3; + const bool is_stress = true; + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); + const double hse_omega = XC_Functional::get_hse_omega(); XC_Functional::gradcorr( dum1, dum2, dum3, chr, rho_basis, &ucell, stress_gga, is_stress, - PARAM.inp.nspin, PARAM.globalv.domag, PARAM.globalv.domag_z); + PARAM.inp.nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + hybrid_alpha, hse_omega); for(int l = 0;l< 3;l++) - { - for(int m = 0;m< l+1;m++) - { - int ind = l*3 + m; - sigma_gradcorr[l][m] = stress_gga [ind]; - sigma_gradcorr[m][l] = sigma_gradcorr[l][m]; - } - } + { + for(int m = 0;m< l+1;m++) + { + int ind = l*3 + m; + sigma_gradcorr[l][m] = stress_gga [ind]; + sigma_gradcorr[m][l] = sigma_gradcorr[l][m]; + } + } - for(int l = 0;l<3;l++) - { - for(int m = 0;m<3;m++) - { + for(int l = 0;l<3;l++) + { + for(int m = 0;m<3;m++) + { Parallel_Reduce::reduce_pool(sigma_gradcorr[l][m]); - } - } - - for(int i=0;i<3;i++) - { - for(int j=0;j<3;j++) - { + } + } + + for(int i=0;i<3;i++) + { + for(int j=0;j<3;j++) + { sigma(i, j) += sigma_gradcorr[i][j] / rho_basis->nxyz; } - } + } - ModuleBase::timer::end("Stress","stress_gga"); - return; + ModuleBase::timer::end("Stress","stress_gga"); + return; } template class Stress_Func; diff --git a/source/source_pw/module_pwdft/stress_pw.cpp b/source/source_pw/module_pwdft/stress_pw.cpp index d511361a147..4b26b7ec9fc 100644 --- a/source/source_pw/module_pwdft/stress_pw.cpp +++ b/source/source_pw/module_pwdft/stress_pw.cpp @@ -126,9 +126,12 @@ void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, } // EXX PW stress - if (GlobalC::exx_info.info_global.cal_exx) + bool cal_exx = GlobalC::exx_info.info_global.cal_exx; + double hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha; + auto coulomb_param = GlobalC::exx_info.info_global.coulomb_param; + if (cal_exx) { - this->stress_exx(sigmaexx, this->pelec->wg, rho_basis, wfc_basis, p_kv, d_psi_in, ucell); + this->stress_exx(sigmaexx, this->pelec->wg, rho_basis, wfc_basis, p_kv, d_psi_in, ucell, hybrid_alpha, coulomb_param); } @@ -169,7 +172,7 @@ void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, { ModuleIO::print_stress("ONSITE STRESS", sigmaonsite, screen, ry, GlobalV::ofs_running); } - if (GlobalC::exx_info.info_global.cal_exx) + if (cal_exx) { ModuleIO::print_stress("EXX STRESS", sigmaexx, screen, ry, GlobalV::ofs_running); } diff --git a/source/source_pw/module_pwdft/stress_pw.h b/source/source_pw/module_pwdft/stress_pw.h index dde4688681c..3defa091286 100644 --- a/source/source_pw/module_pwdft/stress_pw.h +++ b/source/source_pw/module_pwdft/stress_pw.h @@ -5,6 +5,7 @@ #include "source_pw/module_pwdft/vl_pw.h" #include "stress_func.h" #include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-07 +#include "source_lcao/module_ri/conv_coulomb_pot_k.h" template class Stress_PW : public Stress_Func @@ -45,7 +46,9 @@ class Stress_PW : public Stress_Func ModulePW::PW_Basis_K* wfc_basis, const K_Vectors* p_kv, const psi::Psi , Device>* d_psi_in, - const UnitCell& ucell); // exx stress in PW basis + const UnitCell& ucell, + const double hybrid_alpha, + const CoulombParam& coulomb_param); // exx stress in PW basis const elecstate::ElecState* pelec = nullptr; }; diff --git a/source/source_pw/module_pwdft/vsep_pw.cpp b/source/source_pw/module_pwdft/vsep_pw.cpp index 9e7a792bc9f..d5433e944ef 100644 --- a/source/source_pw/module_pwdft/vsep_pw.cpp +++ b/source/source_pw/module_pwdft/vsep_pw.cpp @@ -13,10 +13,6 @@ #include #include -// namespace GlobalC -// { -// VSep vsep_cell; -// } namespace { double sphere_cut(double r, double r_out, double r_power) diff --git a/source/source_pw/module_pwdft/vsep_pw.h b/source/source_pw/module_pwdft/vsep_pw.h index 49d24e019f1..9bec4fe0495 100644 --- a/source/source_pw/module_pwdft/vsep_pw.h +++ b/source/source_pw/module_pwdft/vsep_pw.h @@ -23,10 +23,4 @@ class VSep private: int nrxx = 0; }; -// -// namespace GlobalC -// { -// extern VSep vsep_cell; -// } - #endif /* ifndef VSEP_IN_PW */