diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 9e3c736b89..bf7321d5bc 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -651,7 +651,6 @@ OBJS_SRCPW=H_Ewald_pw.o\ vnl_op.o\ global.o\ magnetism.o\ - elecstate_getters.o\ occupy.o\ structure_factor.o\ structure_factor_k.o\ diff --git a/source/module_elecstate/CMakeLists.txt b/source/module_elecstate/CMakeLists.txt index 2d0148a0d1..2964d3bc44 100644 --- a/source/module_elecstate/CMakeLists.txt +++ b/source/module_elecstate/CMakeLists.txt @@ -1,6 +1,5 @@ list(APPEND objects elecstate.cpp - elecstate_getters.cpp elecstate_energy_terms.cpp elecstate_energy.cpp elecstate_exx.cpp diff --git a/source/module_elecstate/elecstate_energy.cpp b/source/module_elecstate/elecstate_energy.cpp index 4f6e3e1756..3248b14917 100644 --- a/source/module_elecstate/elecstate_energy.cpp +++ b/source/module_elecstate/elecstate_energy.cpp @@ -1,7 +1,7 @@ #include "elecstate.h" -#include "elecstate_getters.h" #include "module_base/global_variable.h" #include "module_base/parallel_reduce.h" +#include "module_hamilt_general/module_xc/xc_functional.h" #include "module_parameter/parameter.h" #include @@ -103,7 +103,7 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const const double* v_eff = this->pot->get_effective_v(0); const double* v_fixed = this->pot->get_fixed_v(); const double* v_ofk = nullptr; - const bool v_ofk_flag = (get_xc_func_type() == 3 || get_xc_func_type() == 5); + const bool v_ofk_flag = (XC_Functional::get_ked_flag()); #ifdef USE_PAW if (PARAM.inp.use_paw) { @@ -208,14 +208,14 @@ double ElecState::cal_delta_escf() const const double* v_fixed = this->pot->get_fixed_v(); const double* v_ofk = nullptr; - if (get_xc_func_type() == 3 || get_xc_func_type() == 5) + if (XC_Functional::get_ked_flag()) { v_ofk = this->pot->get_effective_vofk(0); } for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) { descf -= (this->charge->rho[0][ir] - this->charge->rho_save[0][ir]) * (v_eff[ir] - v_fixed[ir]); - if (get_xc_func_type() == 3 || get_xc_func_type() == 5) + if (XC_Functional::get_ked_flag()) { // cause in the get_effective_vofk, the func will return nullptr assert(v_ofk != nullptr); @@ -226,14 +226,14 @@ double ElecState::cal_delta_escf() const if (PARAM.inp.nspin == 2) { v_eff = this->pot->get_effective_v(1); - if (get_xc_func_type() == 3 || get_xc_func_type() == 5) + if (XC_Functional::get_ked_flag()) { v_ofk = this->pot->get_effective_vofk(1); } for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) { descf -= (this->charge->rho[1][ir] - this->charge->rho_save[1][ir]) * (v_eff[ir] - v_fixed[ir]); - if (get_xc_func_type() == 3 || get_xc_func_type() == 5) + if (XC_Functional::get_ked_flag()) { descf -= (this->charge->kin_r[1][ir] - this->charge->kin_r_save[1][ir]) * v_ofk[ir]; } diff --git a/source/module_elecstate/elecstate_getters.cpp b/source/module_elecstate/elecstate_getters.cpp deleted file mode 100644 index bf39413498..0000000000 --- a/source/module_elecstate/elecstate_getters.cpp +++ /dev/null @@ -1,17 +0,0 @@ -#include "module_elecstate/elecstate_getters.h" - -#include "module_cell/unitcell.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_parameter/parameter.h" -#include "module_hamilt_general/module_xc/xc_functional.h" - -namespace elecstate -{ - -int get_xc_func_type() -{ - return XC_Functional::get_func_type(); -} - - -} // namespace elecstate diff --git a/source/module_elecstate/elecstate_getters.h b/source/module_elecstate/elecstate_getters.h deleted file mode 100644 index e9a7cdf6c1..0000000000 --- a/source/module_elecstate/elecstate_getters.h +++ /dev/null @@ -1,16 +0,0 @@ -#ifndef ELECSTATE_GETTERS_H -#define ELECSTATE_GETTERS_H - -#include - -// Description: Getters for elecstate module -namespace elecstate -{ - -/// @brief get the value of XC_Functional::func_type -int get_xc_func_type(); - - -} // namespace elecstate - -#endif // ELECSTATE_GETTERS_H diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp index 748ef7a9b8..498c1735a6 100644 --- a/source/module_elecstate/elecstate_lcao.cpp +++ b/source/module_elecstate/elecstate_lcao.cpp @@ -63,7 +63,7 @@ void ElecStateLCAO>::psiToRho(const psi::Psicharge->rho, Gint_Tools::job_type::rho, PARAM.inp.nspin); this->gint_k->cal_gint(&inout); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { elecstate::lcao_cal_tau_k(gint_k, this->charge); } @@ -98,7 +98,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi& psi) this->gint_gamma->cal_gint(&inout); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { elecstate::lcao_cal_tau_gamma(gint_gamma, this->charge); } @@ -161,7 +161,7 @@ void ElecStateLCAO::dmToRho(std::vector pexsi_DM, std::vectorgint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho, PARAM.inp.nspin); this->gint_gamma->cal_gint(&inout); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { for (int is = 0; is < PARAM.inp.nspin; is++) { diff --git a/source/module_elecstate/elecstate_print.cpp b/source/module_elecstate/elecstate_print.cpp index fc08baea4e..a223b4eac6 100644 --- a/source/module_elecstate/elecstate_print.cpp +++ b/source/module_elecstate/elecstate_print.cpp @@ -1,5 +1,4 @@ #include "elecstate.h" -#include "elecstate_getters.h" #include "module_base/formatter.h" #include "module_base/global_variable.h" #include "module_base/parallel_common.h" @@ -461,7 +460,7 @@ void ElecState::print_etot(const Magnetism& magnet, break; } std::vector drho = {scf_thr}; - if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5) + if (XC_Functional::get_ked_flag()) { drho.push_back(scf_thr_kin); } diff --git a/source/module_elecstate/elecstate_pw.cpp b/source/module_elecstate/elecstate_pw.cpp index 944bcddef4..f01c70e4e1 100644 --- a/source/module_elecstate/elecstate_pw.cpp +++ b/source/module_elecstate/elecstate_pw.cpp @@ -1,12 +1,13 @@ #include "elecstate_pw.h" -#include "module_parameter/parameter.h" -#include "elecstate_getters.h" + #include "module_base/constants.h" #include "module_base/libm/libm.h" #include "module_base/math_ylmreal.h" +#include "module_base/module_device/device.h" #include "module_base/parallel_reduce.h" #include "module_base/timer.h" -#include "module_base/module_device/device.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_parameter/parameter.h" namespace elecstate { @@ -41,7 +42,7 @@ ElecStatePW::~ElecStatePW() delmem_complex_op()(this->rhog_data); delete[] this->rhog; } - if (get_xc_func_type() == 3 || PARAM.inp.out_elf[0] > 0) + if (XC_Functional::get_func_type() == 3 || PARAM.inp.out_elf[0] > 0) { delmem_var_op()(this->kin_r_data); delete[] this->kin_r; @@ -80,7 +81,7 @@ void ElecStatePW::init_rho_data() this->rhog[ii] = this->rhog_data + ii * this->charge->rhopw->npw; } } - if (get_xc_func_type() == 3 || PARAM.inp.out_elf[0] > 0) + if (XC_Functional::get_func_type() == 3 || PARAM.inp.out_elf[0] > 0) { this->kin_r = new Real*[this->charge->nspin]; resmem_var_op()(this->kin_r_data, this->charge->nspin * this->charge->nrxx); @@ -96,7 +97,7 @@ void ElecStatePW::init_rho_data() { this->rhog = reinterpret_cast(this->charge->rhog); } - if (get_xc_func_type() == 3 || PARAM.inp.out_elf[0] > 0) + if (XC_Functional::get_func_type() == 3 || PARAM.inp.out_elf[0] > 0) { this->kin_r = reinterpret_cast(this->charge->kin_r); } @@ -119,7 +120,7 @@ void ElecStatePW::psiToRho(const psi::Psi& psi) // denghui replaced at 20221110 // ModuleBase::GlobalFunc::ZEROS(this->rho[is], this->charge->nrxx); setmem_var_op()(this->rho[is], 0, this->charge->nrxx); - if (get_xc_func_type() == 3) + if (XC_Functional::get_func_type() == 3) { // ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx); setmem_var_op()(this->kin_r[is], 0, this->charge->nrxx); @@ -143,7 +144,7 @@ void ElecStatePW::psiToRho(const psi::Psi& psi) for (int ii = 0; ii < PARAM.inp.nspin; ii++) { castmem_var_d2h_op()(this->charge->rho[ii], this->rho[ii], this->charge->nrxx); - if (get_xc_func_type() == 3) + if (XC_Functional::get_func_type() == 3) { castmem_var_d2h_op()(this->charge->kin_r[ii], this->kin_r[ii], this->charge->nrxx); } @@ -240,7 +241,7 @@ void ElecStatePW::rhoBandK(const psi::Psi& psi) } // kinetic energy density - if (get_xc_func_type() == 3) + if (XC_Functional::get_func_type() == 3) { for (int j = 0; j < 3; j++) { diff --git a/source/module_elecstate/elecstate_pw_cal_tau.cpp b/source/module_elecstate/elecstate_pw_cal_tau.cpp index 5c225c3d62..628dd25aef 100644 --- a/source/module_elecstate/elecstate_pw_cal_tau.cpp +++ b/source/module_elecstate/elecstate_pw_cal_tau.cpp @@ -1,5 +1,4 @@ #include "elecstate_pw.h" -#include "elecstate_getters.h" namespace elecstate { diff --git a/source/module_elecstate/magnetism.cpp b/source/module_elecstate/magnetism.cpp index 80d4856fc1..b41f5261ad 100644 --- a/source/module_elecstate/magnetism.cpp +++ b/source/module_elecstate/magnetism.cpp @@ -1,7 +1,7 @@ #include "magnetism.h" -#include "elecstate_getters.h" -#include "module_parameter/parameter.h" + #include "module_base/parallel_reduce.h" +#include "module_parameter/parameter.h" Magnetism::Magnetism() { diff --git a/source/module_elecstate/module_charge/charge.cpp b/source/module_elecstate/module_charge/charge.cpp index f152ca850e..dfd9b4d633 100644 --- a/source/module_elecstate/module_charge/charge.cpp +++ b/source/module_elecstate/module_charge/charge.cpp @@ -17,11 +17,9 @@ // even in a LSDA calculation. //---------------------------------------------------------- #include "charge.h" -#include -#include "module_parameter/parameter.h" + #include "module_base/global_function.h" #include "module_base/global_variable.h" -#include "module_parameter/parameter.h" #include "module_base/libm/libm.h" #include "module_base/math_integral.h" #include "module_base/memory.h" @@ -29,8 +27,11 @@ #include "module_base/timer.h" #include "module_base/tool_threading.h" #include "module_cell/unitcell.h" -#include "module_elecstate/elecstate_getters.h" #include "module_elecstate/magnetism.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_parameter/parameter.h" + +#include #ifdef USE_PAW #include "module_cell/module_paw/paw_cell.h" @@ -80,7 +81,7 @@ void Charge::destroy() delete[] _space_rhog_save; delete[] _space_kin_r; delete[] _space_kin_r_save; - if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0) + if (XC_Functional::get_ked_flag() || PARAM.inp.out_elf[0] > 0) { delete[] kin_r; delete[] kin_r_save; @@ -121,7 +122,7 @@ void Charge::allocate(const int& nspin_in) _space_rho_save = new double[nspin * nrxx]; _space_rhog = new std::complex[nspin * ngmc]; _space_rhog_save = new std::complex[nspin * ngmc]; - if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0) + if (XC_Functional::get_ked_flag() || PARAM.inp.out_elf[0] > 0) { _space_kin_r = new double[nspin * nrxx]; _space_kin_r_save = new double[nspin * nrxx]; @@ -130,7 +131,7 @@ void Charge::allocate(const int& nspin_in) rhog = new std::complex*[nspin]; rho_save = new double*[nspin]; rhog_save = new std::complex*[nspin]; - if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0) + if (XC_Functional::get_ked_flag() || PARAM.inp.out_elf[0] > 0) { kin_r = new double*[nspin]; kin_r_save = new double*[nspin]; @@ -151,7 +152,7 @@ void Charge::allocate(const int& nspin_in) ModuleBase::GlobalFunc::ZEROS(rhog[is], ngmc); ModuleBase::GlobalFunc::ZEROS(rho_save[is], nrxx); ModuleBase::GlobalFunc::ZEROS(rhog_save[is], ngmc); - if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0) + if (XC_Functional::get_ked_flag() || PARAM.inp.out_elf[0] > 0) { kin_r[is] = _space_kin_r + is * nrxx; ModuleBase::GlobalFunc::ZEROS(kin_r[is], nrxx); @@ -171,7 +172,7 @@ void Charge::allocate(const int& nspin_in) ModuleBase::Memory::record("Chg::rho_save", sizeof(double) * nspin * nrxx); ModuleBase::Memory::record("Chg::rhog", sizeof(double) * nspin * ngmc); ModuleBase::Memory::record("Chg::rhog_save", sizeof(double) * nspin * ngmc); - if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0) + if (XC_Functional::get_ked_flag() || PARAM.inp.out_elf[0] > 0) { ModuleBase::Memory::record("Chg::kin_r", sizeof(double) * nspin * ngmc); ModuleBase::Memory::record("Chg::kin_r_save", sizeof(double) * nspin * ngmc); @@ -701,9 +702,10 @@ void Charge::save_rho_before_sum_band() for (int is = 0; is < PARAM.inp.nspin; is++) { ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is], this->rhopw->nrxx); - if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5) { + if (XC_Functional::get_ked_flag()) + { ModuleBase::GlobalFunc::DCOPY(kin_r[is], kin_r_save[is], this->rhopw->nrxx); -} + } #ifdef USE_PAW if(PARAM.inp.use_paw) { ModuleBase::GlobalFunc::DCOPY(nhat[is], nhat_save[is], this->rhopw->nrxx); diff --git a/source/module_elecstate/module_charge/charge_init.cpp b/source/module_elecstate/module_charge/charge_init.cpp index f6a0b400d4..d162a458ff 100644 --- a/source/module_elecstate/module_charge/charge_init.cpp +++ b/source/module_elecstate/module_charge/charge_init.cpp @@ -117,7 +117,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, } } - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { // If the charge density is not read in, then the kinetic energy density is not read in either if (!read_error) @@ -188,7 +188,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, } // wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { if (PARAM.inp.init_chg == "atomic" || read_kin_error) { diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 5bc3188190..78fe8ea70c 100644 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -144,7 +144,7 @@ void Charge_Mixing::init_mixing() } // initailize tau_mdata - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) + if ((XC_Functional::get_ked_flag()) && mixing_tau) { if (PARAM.inp.scf_thr_type == 1) { @@ -180,7 +180,7 @@ void Charge_Mixing::mix_reset() this->mixing->reset(); this->rho_mdata.reset(); // initailize tau_mdata - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) + if ((XC_Functional::get_ked_flag()) && mixing_tau) { this->tau_mdata.reset(); } diff --git a/source/module_elecstate/module_charge/charge_mixing_residual.cpp b/source/module_elecstate/module_charge/charge_mixing_residual.cpp index b2d1425597..5d97dc0fce 100644 --- a/source/module_elecstate/module_charge/charge_mixing_residual.cpp +++ b/source/module_elecstate/module_charge/charge_mixing_residual.cpp @@ -72,7 +72,7 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec) double Charge_Mixing::get_dkin(Charge* chr, const double nelec) { - if (!(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)) + if (!(XC_Functional::get_ked_flag())) { return 0.0; }; diff --git a/source/module_elecstate/module_charge/charge_mixing_rho.cpp b/source/module_elecstate/module_charge/charge_mixing_rho.cpp index a3d632c211..772cd69666 100644 --- a/source/module_elecstate/module_charge/charge_mixing_rho.cpp +++ b/source/module_elecstate/module_charge/charge_mixing_rho.cpp @@ -249,7 +249,7 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) } // For kinetic energy density - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) + if ((XC_Functional::get_ked_flag()) && mixing_tau) { std::vector> kin_g(PARAM.inp.nspin * rhodpw->npw); std::vector> kin_g_save(PARAM.inp.nspin * rhodpw->npw); @@ -485,7 +485,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr) } double *taur_out, *taur_in; - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) + if ((XC_Functional::get_ked_flag()) && mixing_tau) { taur_in = chr->kin_r_save[0]; taur_out = chr->kin_r[0]; @@ -521,7 +521,7 @@ void Charge_Mixing::mix_rho(Charge* chr) } } std::vector kin_r123; - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) + if ((XC_Functional::get_ked_flag()) && mixing_tau) { kin_r123.resize(PARAM.inp.nspin * nrxx); for (int is = 0; is < PARAM.inp.nspin; ++is) @@ -581,7 +581,7 @@ void Charge_Mixing::mix_rho(Charge* chr) } } - if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau) + if ((XC_Functional::get_ked_flag()) && mixing_tau) { for (int is = 0; is < PARAM.inp.nspin; ++is) { diff --git a/source/module_elecstate/module_charge/charge_mpi.cpp b/source/module_elecstate/module_charge/charge_mpi.cpp index 61f01ea58b..f55841be6e 100644 --- a/source/module_elecstate/module_charge/charge_mpi.cpp +++ b/source/module_elecstate/module_charge/charge_mpi.cpp @@ -3,7 +3,7 @@ #include "module_base/global_variable.h" #include "module_base/parallel_reduce.h" #include "module_base/timer.h" -#include "module_elecstate/elecstate_getters.h" +#include "module_hamilt_general/module_xc/xc_functional.h" #include "module_parameter/parameter.h" #ifdef __MPI void Charge::init_chgmpi() @@ -132,7 +132,7 @@ void Charge::rho_mpi() for (int is = 0; is < PARAM.inp.nspin; ++is) { reduce_diff_pools(this->rho[is]); - if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0) + if (XC_Functional::get_ked_flag() || PARAM.inp.out_elf[0] > 0) { reduce_diff_pools(this->kin_r[is]); } diff --git a/source/module_elecstate/module_charge/symmetry_rho.cpp b/source/module_elecstate/module_charge/symmetry_rho.cpp index 3455402da5..9bf0c6b3d5 100644 --- a/source/module_elecstate/module_charge/symmetry_rho.cpp +++ b/source/module_elecstate/module_charge/symmetry_rho.cpp @@ -20,27 +20,27 @@ void Symmetry_rho::begin(const int& spin_now, if (ModuleSymmetry::Symmetry::symm_flag != 1) { return; } - // both parallel and serial - // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space - // { - // psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm); - // if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now], - // rho_basis,Pgrid,symm); - // } - // else //space group, do rho_symm in reciprocal space - { - rho_basis->real2recip(CHR.rho[spin_now], CHR.rhog[spin_now]); - psymmg(CHR.rhog[spin_now], rho_basis, symm); // need to modify - rho_basis->recip2real(CHR.rhog[spin_now], CHR.rho[spin_now]); +// both parallel and serial +// if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space +// { +// psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm); +// if(XC_Functional::get_ked_flag()) psymm(CHR.kin_r[spin_now], +// rho_basis,Pgrid,symm); +// } +// else //space group, do rho_symm in reciprocal space +{ + rho_basis->real2recip(CHR.rho[spin_now], CHR.rhog[spin_now]); + psymmg(CHR.rhog[spin_now], rho_basis, symm); // need to modify + rho_basis->recip2real(CHR.rhog[spin_now], CHR.rho[spin_now]); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5 || CHR.cal_elf) - { - // Use std::vector to manage kin_g instead of raw pointer - std::vector> kin_g(CHR.ngmc); - rho_basis->real2recip(CHR.kin_r[spin_now], kin_g.data()); - psymmg(kin_g.data(), rho_basis, symm); - rho_basis->recip2real(kin_g.data(), CHR.kin_r[spin_now]); - } + if (XC_Functional::get_ked_flag() || CHR.cal_elf) + { + // Use std::vector to manage kin_g instead of raw pointer + std::vector> kin_g(CHR.ngmc); + rho_basis->real2recip(CHR.kin_r[spin_now], kin_g.data()); + psymmg(kin_g.data(), rho_basis, symm); + rho_basis->recip2real(kin_g.data(), CHR.kin_r[spin_now]); + } } return; } @@ -63,7 +63,7 @@ void Symmetry_rho::begin(const int& spin_now, // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space // { // psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm); - // if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now], + // if(XC_Functional::get_ked_flag()) psymm(CHR.kin_r[spin_now], // rho_basis,Pgrid,symm); // } // else //space group, do rho_symm in reciprocal space @@ -72,7 +72,7 @@ void Symmetry_rho::begin(const int& spin_now, psymmg(rhog[spin_now], rho_basis, symm); rho_basis->recip2real(rhog[spin_now], rho[spin_now]); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { // Use std::vector to manage kin_g instead of raw pointer std::vector> kin_g(ngmc); diff --git a/source/module_elecstate/occupy.cpp b/source/module_elecstate/occupy.cpp index c4cc61213d..eedf1466ce 100644 --- a/source/module_elecstate/occupy.cpp +++ b/source/module_elecstate/occupy.cpp @@ -1,10 +1,9 @@ #include "occupy.h" -#include "module_parameter/parameter.h" #include "module_base/constants.h" #include "module_base/mymath.h" #include "module_base/parallel_reduce.h" -#include "module_elecstate/elecstate_getters.h" +#include "module_parameter/parameter.h" Occupy::Occupy() { @@ -142,10 +141,12 @@ void Occupy::iweights( { assert(is < 2); double degspin = 2.0; - if (PARAM.inp.nspin == 4) + if (PARAM.inp.nspin == 4) { degspin = 1.0; - if (is != -1) +} + if (is != -1) { degspin = 1.0; +} double ib_mind = nelec / degspin; int ib_min = std::ceil(ib_mind); @@ -225,8 +226,9 @@ void Occupy::gweights(const int nks, for (int ik = 0; ik < nks; ik++) { // mohan add 2011-04-03 - if (is != -1 && is != isk[ik]) + if (is != -1 && is != isk[ik]) { continue; +} for (int ib = 0; ib < PARAM.inp.nbands; ib++) { @@ -400,8 +402,9 @@ double Occupy::sumkg(const ModuleBase::matrix& ekb, double sum2 = 0.0; for (int ik = 0; ik < nks; ik++) { - if (is != -1 && is != isk[ik]) + if (is != -1 && is != isk[ik]) { continue; +} double sum1 = 0.0; for (int ib = 0; ib < nband; ib++) diff --git a/source/module_elecstate/potentials/pot_xc.cpp b/source/module_elecstate/potentials/pot_xc.cpp index b04370b2f5..f6b4607f2d 100644 --- a/source/module_elecstate/potentials/pot_xc.cpp +++ b/source/module_elecstate/potentials/pot_xc.cpp @@ -20,7 +20,7 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module // calculate the exchange-correlation potential //---------------------------------------------------------- - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { #ifdef USE_LIBXC const std::tuple etxc_vtxc_v diff --git a/source/module_elecstate/potentials/potential_new.cpp b/source/module_elecstate/potentials/potential_new.cpp index f3d68df05a..a5e992beee 100644 --- a/source/module_elecstate/potentials/potential_new.cpp +++ b/source/module_elecstate/potentials/potential_new.cpp @@ -1,18 +1,18 @@ #include "potential_new.h" -#include "module_parameter/parameter.h" #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "module_base/memory.h" #include "module_base/timer.h" #include "module_base/tool_quit.h" #include "module_base/tool_title.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_parameter/parameter.h" #ifdef USE_PAW #include "module_hamilt_general/module_xc/xc_functional.h" #include "module_cell/module_paw/paw_cell.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #endif -#include "module_elecstate/elecstate_getters.h" #include @@ -123,7 +123,7 @@ void Potential::allocate() ModuleBase::Memory::record("Pot::vxc", sizeof(double) * PARAM.inp.nspin * nrxx); } - if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5) + if (XC_Functional::get_ked_flag()) { this->vofk_effective.create(PARAM.inp.nspin, nrxx); ModuleBase::Memory::record("Pot::vofk", sizeof(double) * PARAM.inp.nspin * nrxx); @@ -320,7 +320,7 @@ void Potential::interpolate_vrs() rho_basis_smooth_->recip2real(&vrs(is, 0), &veff_smooth(is, 0)); } - if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5) + if (XC_Functional::get_ked_flag()) { ModuleBase::ComplexMatrix vrs_ofk(PARAM.inp.nspin, rho_basis_->npw); for (int is = 0; is < PARAM.inp.nspin; is++) diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index 049ef74fc2..2580730cf8 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -1,21 +1,18 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" #define private public -#include "module_parameter/parameter.h" -#undef private -#include "module_base/module_mixing/broyden_mixing.h" -#define private public -#include "module_parameter/parameter.h" #include "../module_charge/charge_mixing.h" +#include "module_base/module_mixing/broyden_mixing.h" #include "module_basis/module_pw/pw_basis.h" #include "module_hamilt_general/module_xc/xc_functional.h" -#undef private +#include "module_parameter/parameter.h" #ifdef _OPENMP #include #endif -int FUNC_TYPE = 1; +int XC_Functional::func_type = 1; +bool XC_Functional::ked_flag = false; // mock function Magnetism::~Magnetism() @@ -35,10 +32,6 @@ void Charge::set_rhopw(ModulePW::PW_Basis* rhopw_in) { this->rhopw = rhopw_in; } -int XC_Functional::get_func_type() -{ - return FUNC_TYPE; -} #ifdef __LCAO InfoNonlocal::InfoNonlocal() { @@ -232,7 +225,8 @@ TEST_F(ChargeMixingTest, InitMixingTest) omp_set_num_threads(1); #endif PARAM.input.nspin = 1; - FUNC_TYPE = 1; + XC_Functional::func_type = 1; + XC_Functional::ked_flag = false; Charge_Mixing CMtest; CMtest.set_rhopw(&pw_basis, &pw_basis); @@ -275,7 +269,8 @@ TEST_F(ChargeMixingTest, InitMixingTest) PARAM.input.mixing_dmr, ucell.omega, ucell.tpiba); - FUNC_TYPE = 3; + XC_Functional::func_type = 3; + XC_Functional::ked_flag = true; CMtest.init_mixing(); EXPECT_EQ(CMtest.tau_mdata.length, pw_basis.nrxx); @@ -826,7 +821,8 @@ TEST_F(ChargeMixingTest, MixRhoTest) charge.set_rhopw(&pw_basis); const int nspin = PARAM.input.nspin = 1; PARAM.sys.domag_z = false; - FUNC_TYPE = 3; + XC_Functional::func_type = 3; + XC_Functional::ked_flag = true; PARAM.input.mixing_beta = 0.7; PARAM.input.mixing_ndim = 1; PARAM.input.mixing_gg0 = 0.0; @@ -964,7 +960,8 @@ TEST_F(ChargeMixingTest, MixDoubleGridRhoTest) charge.set_rhopw(&pw_dbasis); const int nspin = PARAM.input.nspin = 1; PARAM.sys.domag_z = false; - FUNC_TYPE = 3; + XC_Functional::func_type = 3; + XC_Functional::ked_flag = true; PARAM.input.mixing_beta = 0.7; PARAM.input.mixing_ndim = 1; PARAM.input.mixing_gg0 = 0.0; diff --git a/source/module_elecstate/test/charge_test.cpp b/source/module_elecstate/test/charge_test.cpp index b6c98f85a2..85a0d965c6 100644 --- a/source/module_elecstate/test/charge_test.cpp +++ b/source/module_elecstate/test/charge_test.cpp @@ -3,13 +3,11 @@ #define private public #define protected public -#include "module_parameter/parameter.h" #include "module_cell/unitcell.h" #include "module_elecstate/module_charge/charge.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_parameter/parameter.h" #include "prepare_unitcell.h" -#undef protected -#undef private - // mock functions for UnitCell #ifdef __LCAO InfoNonlocal::InfoNonlocal() @@ -31,13 +29,10 @@ Magnetism::~Magnetism() } // mock functions for Charge +int XC_Functional::func_type = 1; +bool XC_Functional::ked_flag = false; namespace elecstate { -int tmp_xc_func_type = 1; -int get_xc_func_type() -{ - return tmp_xc_func_type; -} double tmp_ucell_omega = 500.0; double tmp_gridecut = 80.0; void Set_GlobalV_Default() @@ -118,7 +113,8 @@ TEST_F(ChargeTest, Allocate) EXPECT_EQ(rhopw->npwtot, 3143); // call Charge::allocate() PARAM.input.test_charge = 2; - elecstate::tmp_xc_func_type = 3; + XC_Functional::func_type = 3; + XC_Functional::ked_flag = true; charge->set_rhopw(rhopw); EXPECT_FALSE(charge->allocate_rho); charge->allocate(PARAM.input.nspin); @@ -202,7 +198,8 @@ TEST_F(ChargeTest, SaveRhoBeforeSumBand) } } EXPECT_EQ(PARAM.input.nelec, 8); - elecstate::tmp_xc_func_type = 3; + XC_Functional::func_type = 3; + XC_Functional::ked_flag = true; charge->set_omega(&ucell->omega);; charge->renormalize_rho(); charge->save_rho_before_sum_band(); @@ -212,7 +209,8 @@ TEST_F(ChargeTest, SaveRhoBeforeSumBand) TEST_F(ChargeTest, InitFinalScf) { charge->set_rhopw(rhopw); - elecstate::tmp_xc_func_type = 1; + XC_Functional::func_type = 1; + XC_Functional::ked_flag = false; PARAM.input.test_charge = 2; charge->init_final_scf(); EXPECT_TRUE(charge->allocate_rho_final_scf); diff --git a/source/module_elecstate/test/elecstate_energy_test.cpp b/source/module_elecstate/test/elecstate_energy_test.cpp index bc2f5eb4be..b2040447ac 100644 --- a/source/module_elecstate/test/elecstate_energy_test.cpp +++ b/source/module_elecstate/test/elecstate_energy_test.cpp @@ -2,21 +2,18 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" #define private public -#include "module_parameter/parameter.h" -#undef private #include "module_elecstate/elecstate.h" -#include "module_elecstate/elecstate_getters.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_parameter/parameter.h" + #include Parameter PARMA; // mock functions +int XC_Functional::func_type = 1; +bool XC_Functional::ked_flag = false; namespace elecstate { -int tmp_xc_func_type = 1; -int get_xc_func_type() -{ - return tmp_xc_func_type; -} void Potential::get_vnew(Charge const*, ModuleBase::matrix&) { return; diff --git a/source/module_elecstate/test/elecstate_magnetism_test.cpp b/source/module_elecstate/test/elecstate_magnetism_test.cpp index a527de29d5..5416b1dff2 100644 --- a/source/module_elecstate/test/elecstate_magnetism_test.cpp +++ b/source/module_elecstate/test/elecstate_magnetism_test.cpp @@ -5,7 +5,6 @@ #define private public #include "module_parameter/parameter.h" #undef private -#include "module_elecstate/elecstate_getters.h" /************************************************ * unit test of magnetism.cpp diff --git a/source/module_elecstate/test/elecstate_occupy_test.cpp b/source/module_elecstate/test/elecstate_occupy_test.cpp index dd9342987e..9f3f1546a9 100644 --- a/source/module_elecstate/test/elecstate_occupy_test.cpp +++ b/source/module_elecstate/test/elecstate_occupy_test.cpp @@ -4,7 +4,6 @@ #define private public #include "module_parameter/parameter.h" #undef private -#include "module_elecstate/elecstate_getters.h" /*************************************************************** * unit test of class Occupy diff --git a/source/module_elecstate/test/elecstate_print_test.cpp b/source/module_elecstate/test/elecstate_print_test.cpp index e05759f6ca..dae3b5904f 100644 --- a/source/module_elecstate/test/elecstate_print_test.cpp +++ b/source/module_elecstate/test/elecstate_print_test.cpp @@ -3,14 +3,13 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" #define private public -#include "module_parameter/parameter.h" -#undef private +#include "module_cell/klist.h" #include "module_elecstate/elecstate.h" -#include "module_elecstate/elecstate_getters.h" +#include "module_elecstate/module_charge/charge.h" #include "module_elecstate/potentials/efield.h" #include "module_elecstate/potentials/gatefield.h" -#include "module_elecstate/module_charge/charge.h" -#include "module_cell/klist.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_parameter/parameter.h" /*************************************************************** * mock functions @@ -41,10 +40,9 @@ Charge::Charge() Charge::~Charge() { } -int elecstate::get_xc_func_type() -{ - return 0; -} + +int XC_Functional::func_type = 0; +bool XC_Functional::ked_flag = false; /*************************************************************** * unit test of functions in elecstate_print.cpp diff --git a/source/module_elecstate/test/elecstate_pw_test.cpp b/source/module_elecstate/test/elecstate_pw_test.cpp index a2a3533ade..d45417794d 100644 --- a/source/module_elecstate/test/elecstate_pw_test.cpp +++ b/source/module_elecstate/test/elecstate_pw_test.cpp @@ -3,20 +3,15 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" #define private public -#include "module_parameter/parameter.h" -#undef private #define protected public #include "module_elecstate/elecstate_pw.h" +#include "module_hamilt_general/module_xc/xc_functional.h" #include "module_hamilt_pw/hamilt_pwdft/VL_in_pw.h" -#undef protected +#include "module_parameter/parameter.h" // mock functions for testing +int XC_Functional::func_type = 1; namespace elecstate { -int tmp_xc_func_type = 1; -int get_xc_func_type() -{ - return tmp_xc_func_type; -} void Potential::init_pot(int, Charge const*) { } @@ -272,7 +267,7 @@ TEST_F(ElecStatePWTest, ConstructorSingle) TEST_F(ElecStatePWTest, InitRhoDataDouble) { - elecstate::tmp_xc_func_type = 3; + XC_Functional::func_type = 3; chg->nrxx = 1000; elecstate_pw_d = new elecstate::ElecStatePW, base_device::DEVICE_CPU>(wfcpw, chg, @@ -291,7 +286,7 @@ TEST_F(ElecStatePWTest, InitRhoDataDouble) TEST_F(ElecStatePWTest, InitRhoDataSingle) { PARAM.input.precision = "single"; - elecstate::tmp_xc_func_type = 3; + XC_Functional::func_type = 3; chg->nspin = PARAM.input.nspin; chg->nrxx = 1000; elecstate_pw_s = new elecstate::ElecStatePW, base_device::DEVICE_CPU>(wfcpw, diff --git a/source/module_elecstate/test/potential_new_test.cpp b/source/module_elecstate/test/potential_new_test.cpp index 7115886bff..4e40756daa 100644 --- a/source/module_elecstate/test/potential_new_test.cpp +++ b/source/module_elecstate/test/potential_new_test.cpp @@ -2,9 +2,9 @@ #include #define private public -#include "module_parameter/parameter.h" #include "module_elecstate/potentials/potential_new.h" -#undef private +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_parameter/parameter.h" // mock functions Structure_Factor::Structure_Factor() { @@ -44,14 +44,10 @@ surchem::surchem() surchem::~surchem() { } - +int XC_Functional::func_type = 1; +bool XC_Functional::ked_flag = false; namespace elecstate { -int tmp_xc_func_type = 1; -int get_xc_func_type() -{ - return tmp_xc_func_type; -} PotBase* Potential::get_pot_type(const std::string& pot_type) { @@ -197,7 +193,8 @@ TEST_F(PotentialNewTest, ConstructorNRXX0) TEST_F(PotentialNewTest, ConstructorXC3) { - elecstate::tmp_xc_func_type = 3; + XC_Functional::func_type = 3; + XC_Functional::ked_flag = true; rhopw->nrxx = 100; pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc); EXPECT_TRUE(pot->fixed_mode); @@ -463,7 +460,8 @@ TEST_F(PotentialNewTest, GetEffectiveVarrayNullptr) TEST_F(PotentialNewTest, GetEffectiveVofkmatrix) { // construct potential - elecstate::tmp_xc_func_type = 3; + XC_Functional::func_type = 3; + XC_Functional::ked_flag = true; rhopw->nrxx = 100; pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc); // @@ -532,7 +530,8 @@ TEST_F(PotentialNewTest, GetVeffSmooth) { // construct potential rhopw->nrxx = 100; - elecstate::tmp_xc_func_type = 3; + XC_Functional::func_type = 3; + XC_Functional::ked_flag = true; pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc); // ModuleBase::matrix veff_smooth_tmp = pot->get_veff_smooth(); @@ -576,28 +575,29 @@ TEST_F(PotentialNewTest, GetVofkSmooth) TEST_F(PotentialNewTest, InterpolateVrsDoubleGrids) { PARAM.sys.double_grid = true; - elecstate::tmp_xc_func_type = 3; - // Init pw_basis - rhopw->initgrids(4, ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1), 4); - rhopw->initparameters(false, 4); - rhopw->setuptransform(); - rhopw->collect_local_pw(); - - rhodpw->initgrids(4, ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1), 6); - rhodpw->initparameters(false, 6); - static_cast(rhodpw)->setuptransform(rhopw); - rhodpw->collect_local_pw(); - - pot = new elecstate::Potential(rhodpw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc); - - for (int ir = 0; ir < pot->v_effective.nr; ir++) - { - for (int ic = 0; ic < pot->v_effective.nc; ic++) - { - pot->v_effective(ir,ic) = ir+ic; - pot->vofk_effective(ir,ic) = ir+2*ic; - } - } + XC_Functional::func_type = 3; + XC_Functional::ked_flag = true; + // Init pw_basis + rhopw->initgrids(4, ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1), 4); + rhopw->initparameters(false, 4); + rhopw->setuptransform(); + rhopw->collect_local_pw(); + + rhodpw->initgrids(4, ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1), 6); + rhodpw->initparameters(false, 6); + static_cast(rhodpw)->setuptransform(rhopw); + rhodpw->collect_local_pw(); + + pot = new elecstate::Potential(rhodpw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc); + + for (int ir = 0; ir < pot->v_effective.nr; ir++) + { + for (int ic = 0; ic < pot->v_effective.nc; ic++) + { + pot->v_effective(ir, ic) = ir + ic; + pot->vofk_effective(ir, ic) = ir + 2 * ic; + } + } pot->interpolate_vrs(); @@ -641,23 +641,24 @@ TEST_F(PotentialNewTest, InterpolateVrsWarningQuit) TEST_F(PotentialNewTest, InterpolateVrsSingleGrids) { PARAM.sys.double_grid = false; - elecstate::tmp_xc_func_type = 3; - // Init pw_basis - rhopw->initgrids(4, ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1), 4); - rhopw->initparameters(false, 4); - rhopw->setuptransform(); - rhopw->collect_local_pw(); - - pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc); - - for (int ir = 0; ir < pot->v_effective.nr; ir++) - { - for (int ic = 0; ic < pot->v_effective.nc; ic++) - { - pot->v_effective(ir,ic) = ir+ic; - pot->vofk_effective(ir,ic) = ir+2*ic; - } - } + XC_Functional::func_type = 3; + XC_Functional::ked_flag = true; + // Init pw_basis + rhopw->initgrids(4, ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1), 4); + rhopw->initparameters(false, 4); + rhopw->setuptransform(); + rhopw->collect_local_pw(); + + pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, solvent, etxc, vtxc); + + for (int ir = 0; ir < pot->v_effective.nr; ir++) + { + for (int ic = 0; ic < pot->v_effective.nc; ic++) + { + pot->v_effective(ir, ic) = ir + ic; + pot->vofk_effective(ir, ic) = ir + 2 * ic; + } + } pot->interpolate_vrs(); diff --git a/source/module_elecstate/test_mpi/charge_mpi_test.cpp b/source/module_elecstate/test_mpi/charge_mpi_test.cpp index 040f5953d6..fda923cb86 100644 --- a/source/module_elecstate/test_mpi/charge_mpi_test.cpp +++ b/source/module_elecstate/test_mpi/charge_mpi_test.cpp @@ -1,12 +1,13 @@ +#include "gmock/gmock.h" +#include "gtest/gtest.h" +#define private public #include "module_base/matrix3.h" #include "module_base/parallel_global.h" -#define private public -#include "module_parameter/parameter.h" -#undef private #include "module_elecstate/module_charge/charge.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_parameter/parameter.h" -#include "gmock/gmock.h" -#include "gtest/gtest.h" +bool XC_Functional::ked_flag = false; Charge::Charge() { } @@ -15,14 +16,6 @@ Charge::~Charge() delete[] rec; delete[] dis; } -namespace elecstate -{ -int tmp_xc_func_type = 3; -int get_xc_func_type() -{ - return tmp_xc_func_type; -} -} // namespace elecstate auto sum_array = [](const double* v, const int& nv) { double sum = 0; diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 80ad1ad349..989e4cbfba 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -169,7 +169,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep) &(ucell), PARAM.inp.out_chg[1], 1); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_TAU.cube"; ModuleIO::write_vdata_palgrid(Pgrid, @@ -310,7 +310,7 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter) GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { std::vector> kin_g_space(PARAM.inp.nspin * this->pelec->charge->ngmc, {0.0, 0.0}); std::vector*> kin_g; diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 8de7f56a05..53b55cc67d 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -666,7 +666,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i // get mtaGGA related parameters double dkin = 0.0; // for meta-GGA - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { dkin = p_chgmix->get_dkin(pelec->charge, PARAM.inp.nelec); } diff --git a/source/module_hamilt_general/module_xc/test/test_xc1.cpp b/source/module_hamilt_general/module_xc/test/test_xc1.cpp index 9a5c45eebd..917a6e4a2c 100644 --- a/source/module_hamilt_general/module_xc/test/test_xc1.cpp +++ b/source/module_hamilt_general/module_xc/test/test_xc1.cpp @@ -60,6 +60,7 @@ class XCTest_SCAN0 : public XCTest TEST_F(XCTest_SCAN0, set_xc_type) { EXPECT_EQ(XC_Functional::get_func_type(),5); + EXPECT_TRUE(XC_Functional::get_ked_flag()); } class XCTest_KSDT : public XCTest @@ -102,6 +103,7 @@ class XCTest_R2SCAN : public XCTest TEST_F(XCTest_R2SCAN, set_xc_type) { EXPECT_EQ(XC_Functional::get_func_type(),3); + EXPECT_TRUE(XC_Functional::get_ked_flag()); } class XCTest_LB07 : public XCTest @@ -130,6 +132,7 @@ class XCTest_BMK : public XCTest TEST_F(XCTest_BMK, set_xc_type) { EXPECT_EQ(XC_Functional::get_func_type(),5); + EXPECT_TRUE(XC_Functional::get_ked_flag()); } class XCTest_HF : public XCTest diff --git a/source/module_hamilt_general/module_xc/test/test_xc4.cpp b/source/module_hamilt_general/module_xc/test/test_xc4.cpp index 41fd60f425..946af8a1cc 100644 --- a/source/module_hamilt_general/module_xc/test/test_xc4.cpp +++ b/source/module_hamilt_general/module_xc/test/test_xc4.cpp @@ -56,6 +56,7 @@ class XCTest_SCAN : public XCTest TEST_F(XCTest_SCAN, set_xc_type) { EXPECT_EQ(XC_Functional::get_func_type(),3); + EXPECT_TRUE(XC_Functional::get_ked_flag()); std::vector e_ref = {-1.802228724,-1.802210367,-1.526620682,-0.03583190143,-19035.36827}; std::vector v1_ref = {-1.40580492,-1.405589907,-1.349177581,-0.5343722373,-14.09114168}; std::vector v2_ref = {-0.002302953645,-0.002306239202,-0.002688619206,-0.07677599434,-3.0633237e-07}; diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/source/module_hamilt_general/module_xc/xc_functional.cpp index 4b9b95a2ce..5fc030ec13 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional.cpp @@ -16,6 +16,7 @@ XC_Functional::~XC_Functional(){} std::vector XC_Functional::func_id(1); 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; std::map XC_Functional::scaling_factor_xc = { {1, 1.0} }; // added by jghan, 2024-10-10 @@ -25,15 +26,6 @@ void XC_Functional::set_hybrid_alpha(const double alpha_in) hybrid_alpha = alpha_in; } -double XC_Functional::get_hybrid_alpha() -{ - return hybrid_alpha; -} - -int XC_Functional::get_func_type() -{ - return func_type; -} void XC_Functional::set_xc_first_loop(const UnitCell& ucell) { /** In the special "two-level" calculation case, @@ -266,10 +258,15 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) #endif } - if (func_id[0] == XC_GGA_X_OPTX) - { - std::cerr << "\n OPTX untested please test,"; - } + if (func_type == 3 || func_type == 5) + { + ked_flag = true; + } + + if (func_id[0] == XC_GGA_X_OPTX) + { + std::cerr << "\n OPTX untested please test,"; + } if((func_type == 4 || func_type == 5) && PARAM.inp.basis_type == "pw") { @@ -279,10 +276,6 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) { ModuleBase::WARNING_QUIT("set_xc_type","meta-GGA has not been implemented for nspin = 4 yet"); } - //if((func_type == 3 || func_type == 5) && PARAM.inp.cal_stress == 1 && PARAM.inp.nspin!=1) - //{ - // ModuleBase::WARNING_QUIT("set_xc_type","mgga stress not implemented for polarized case yet"); - //} #ifndef __EXX if(func_type == 4 || func_type == 5) diff --git a/source/module_hamilt_general/module_xc/xc_functional.h b/source/module_hamilt_general/module_xc/xc_functional.h index 04ca90216f..c0fbc91744 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.h +++ b/source/module_hamilt_general/module_xc/xc_functional.h @@ -64,12 +64,22 @@ class XC_Functional // func_type, which is as specified in get_func_type // use_libxc, whether to use LIBXC. The rule is to NOT use it for functionals that we already have. - static int get_func_type(); + static int get_func_type() + { + return func_type; + }; static void set_xc_type(const std::string xc_func_in); // For hybrid functional static void set_hybrid_alpha(const double alpha_in); - static double get_hybrid_alpha(); + static double get_hybrid_alpha() + { + return hybrid_alpha; + }; + static bool get_ked_flag() + { + return ked_flag; + }; /// Usually in exx caculation, the first SCF loop should be converged with PBE static void set_xc_first_loop(const UnitCell& ucell); @@ -77,10 +87,11 @@ class XC_Functional static std::vector func_id; // libxc id of functional static int func_type; //0:none, 1:lda, 2:gga, 3:mgga, 4:hybrid lda/gga, 5:hybrid mgga - static bool use_libxc; + static bool ked_flag; // whether the functional has kinetic energy density + static bool use_libxc; - //exx_hybrid_alpha for mixing exx in hybrid functional: - static double hybrid_alpha; + // exx_hybrid_alpha for mixing exx in hybrid functional: + static double hybrid_alpha; // added by jghan, 2024-07-07 // as a scaling factor for different xc-functionals diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp index aeb5d55c01..e4e04a6a20 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp @@ -77,7 +77,7 @@ void Veff>::contributeHR() // if you change the place of the following code, // rememeber to delete the #include - if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) + if (XC_Functional::get_ked_flag()) { Gint_inout inout(vr_eff1, vofk_eff1, 0, Gint_Tools::job_type::vlocal_meta); this->GK->cal_gint(&inout); @@ -96,12 +96,12 @@ void Veff>::contributeHR() for (int is = 1; is < 4; is++) { vr_eff1 = this->pot->get_effective_v(is); - if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) + if (XC_Functional::get_ked_flag()) { vofk_eff1 = this->pot->get_effective_vofk(is); } - - if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) + + if (XC_Functional::get_ked_flag()) { Gint_inout inout(vr_eff1, vofk_eff1, is, Gint_Tools::job_type::vlocal_meta); this->GK->cal_gint(&inout); @@ -141,7 +141,7 @@ void Veff>::contributeHR(void) // and diagonalize the H matrix (T+Vl+Vnl). //-------------------------------------------- - if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) + if (XC_Functional::get_ked_flag()) { Gint_inout inout(vr_eff1, vofk_eff1, Gint_Tools::job_type::vlocal_meta); this->GG->cal_vlocal(&inout, this->new_e_iteration); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp index feb567eca2..56f7e561e2 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp @@ -23,7 +23,7 @@ namespace PulayForceStress { const double* vr_eff1 = pot->get_effective_v(is); const double* vofk_eff1 = nullptr; - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { vofk_eff1 = pot->get_effective_vofk(is); Gint_inout inout(is, vr_eff1, vofk_eff1, isforce, isstress, &f, &s, Gint_Tools::job_type::force_meta); diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp index 41184b11d0..9c83581f5e 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp @@ -57,7 +57,7 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, ModuleBase::matrix v(PARAM.inp.nspin, rho_basis->nrxx); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { #ifdef USE_LIBXC const auto etxc_vtxc_v diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp index bbdefb737a..af297fba60 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp @@ -51,8 +51,8 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, //recalculate the exchange-correlation potential ModuleBase::matrix vxc; - if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { + 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); diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp index 36ea369d60..f0ca34ea3a 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp @@ -82,7 +82,7 @@ void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, sigmaxc(i, i) = -(pelec->f_en.etxc - pelec->f_en.vtxc) / ucell.omega; } this->stress_gga(ucell, sigmaxc, rho_basis, pelec->charge); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { this->stress_mgga(ucell, sigmaxc, diff --git a/source/module_io/test/read_wfc_to_rho_test.cpp b/source/module_io/test/read_wfc_to_rho_test.cpp index 7717b1bbc0..9272a7c0a9 100644 --- a/source/module_io/test/read_wfc_to_rho_test.cpp +++ b/source/module_io/test/read_wfc_to_rho_test.cpp @@ -4,25 +4,23 @@ #undef __LCAO #define private public -#include "module_elecstate/module_charge/charge.h" -#include "module_parameter/parameter.h" -#undef private - -#ifdef __MPI -#include "module_base/parallel_global.h" -#include "module_basis/module_pw/test/test_tool.h" -#include "mpi.h" -#endif #include "module_cell/klist.h" #include "module_cell/unitcell.h" -#include "module_elecstate/elecstate_getters.h" +#include "module_elecstate/module_charge/charge.h" #include "module_elecstate/module_charge/symmetry_rho.h" #include "module_hamilt_general/module_xc/xc_functional.h" #include "module_hamilt_pw/hamilt_pwdft/parallel_grid.h" #include "module_io/read_wfc_to_rho.h" #include "module_io/write_wfc_pw.h" +#include "module_parameter/parameter.h" #include "module_psi/psi.h" +#ifdef __MPI +#include "module_base/parallel_global.h" +#include "module_basis/module_pw/test/test_tool.h" +#include "mpi.h" +#endif + Parallel_Grid::Parallel_Grid() { } @@ -47,14 +45,9 @@ Magnetism::Magnetism() Magnetism::~Magnetism() { } -int elecstate::get_xc_func_type() -{ - return 0; -} -int XC_Functional::get_func_type() -{ - return 0; -} +int XC_Functional::func_type = 0; +bool XC_Functional::ked_flag = false; + Symmetry_rho::Symmetry_rho() { } diff --git a/source/module_rdmft/update_state_rdmft.cpp b/source/module_rdmft/update_state_rdmft.cpp index dc0398e8c9..6d9f829f06 100644 --- a/source/module_rdmft/update_state_rdmft.cpp +++ b/source/module_rdmft/update_state_rdmft.cpp @@ -110,7 +110,7 @@ void RDMFT::update_charge(UnitCell& ucell) Gint_inout inout(charge->rho, Gint_Tools::job_type::rho, nspin); GG->cal_gint(&inout); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { // for (int is = 0; is < nspin; is++) // { @@ -140,7 +140,7 @@ void RDMFT::update_charge(UnitCell& ucell) Gint_inout inout(charge->rho, Gint_Tools::job_type::rho, nspin); GK->cal_gint(&inout); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + if (XC_Functional::get_ked_flag()) { // for (int is = 0; is < nspin; is++) // {