diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp index 6e30414531..5a517b23e0 100644 --- a/source/module_elecstate/elecstate_lcao.cpp +++ b/source/module_elecstate/elecstate_lcao.cpp @@ -62,7 +62,7 @@ void ElecStateLCAO>::psiToRho(const psi::Psigint_k->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); + Gint_inout inout(this->charge->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) @@ -113,7 +113,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi& psi) this->gint_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); + Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho, PARAM.inp.nspin); this->gint_gamma->cal_gint(&inout); @@ -178,7 +178,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); + 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) { diff --git a/source/module_elecstate/elecstate_lcao_cal_tau.cpp b/source/module_elecstate/elecstate_lcao_cal_tau.cpp index 01ee54f03f..e6bd6561a0 100644 --- a/source/module_elecstate/elecstate_lcao_cal_tau.cpp +++ b/source/module_elecstate/elecstate_lcao_cal_tau.cpp @@ -15,7 +15,7 @@ void ElecStateLCAO>::cal_tau(const psi::Psicharge->kin_r[is], this->charge->nrxx); } - Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau); + Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau, PARAM.inp.nspin); this->gint_k->cal_gint(&inout1); ModuleBase::timer::tick("ElecStateLCAO", "cal_tau"); @@ -32,7 +32,7 @@ void ElecStateLCAO::cal_tau(const psi::Psi& psi) { ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx); } - Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau); + Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau, PARAM.inp.nspin); this->gint_gamma->cal_gint(&inout1); ModuleBase::timer::tick("ElecStateLCAO", "cal_tau"); diff --git a/source/module_elecstate/elecstate_lcao_tddft.cpp b/source/module_elecstate/elecstate_lcao_tddft.cpp index 196389404e..147b790a95 100644 --- a/source/module_elecstate/elecstate_lcao_tddft.cpp +++ b/source/module_elecstate/elecstate_lcao_tddft.cpp @@ -36,7 +36,7 @@ void ElecStateLCAO_TDDFT::psiToRho_td(const psi::Psi>& psi) ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!"); this->gint_k->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); // rho calculation + Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho, PARAM.inp.nspin); // rho calculation this->gint_k->cal_gint(&inout); this->charge->renormalize_rho(); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index 650532665d..bff12cb2ca 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -174,13 +174,14 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, pot_in->pot_register(pot_register_in); // effective potential term Operator* veff = new Veff>(GG_in, - this->hsk, - this->kv->kvec_d, - pot_in, - this->hR, // no explicit call yet - &GlobalC::ucell, - orb.cutoffs(), - &GlobalC::GridD + this->hsk, + this->kv->kvec_d, + pot_in, + this->hR, // no explicit call yet + &GlobalC::ucell, + orb.cutoffs(), + &GlobalC::GridD, + PARAM.inp.nspin ); this->getOperator()->add(veff); } @@ -242,13 +243,14 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, pot_in->pot_register(pot_register_in); // Veff term this->getOperator() = new Veff>(GK_in, - this->hsk, - kv->kvec_d, - pot_in, - this->hR, - &GlobalC::ucell, - orb.cutoffs(), - &GlobalC::GridD); + this->hsk, + kv->kvec_d, + pot_in, + this->hR, + &GlobalC::ucell, + orb.cutoffs(), + &GlobalC::GridD, + PARAM.inp.nspin); // reset spin index and real space Hamiltonian matrix int start_spin = -1; GK_in->reset_spin(start_spin); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h index ee92f87777..37ac187224 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h @@ -36,39 +36,41 @@ class Veff> : public OperatorLCAO * @param GK_in: the pointer of Gint_k object, used for grid integration */ Veff>(Gint_k* GK_in, - HS_Matrix_K* hsk_in, - const std::vector>& kvec_d_in, - elecstate::Potential* pot_in, - hamilt::HContainer* hR_in, - const UnitCell* ucell_in, - const std::vector& orb_cutoff, - Grid_Driver* GridD_in) + HS_Matrix_K* hsk_in, + const std::vector>& kvec_d_in, + elecstate::Potential* pot_in, + hamilt::HContainer* hR_in, + const UnitCell* ucell_in, + const std::vector& orb_cutoff, + Grid_Driver* GridD_in, + const int& nspin) : GK(GK_in), orb_cutoff_(orb_cutoff), pot(pot_in), ucell(ucell_in), gd(GridD_in), OperatorLCAO(hsk_in, kvec_d_in, hR_in) { this->cal_type = calculation_type::lcao_gint; this->initialize_HR(ucell_in, GridD_in); - GK_in->initialize_pvpR(*ucell_in, GridD_in); + GK_in->initialize_pvpR(*ucell_in, GridD_in, nspin); } /** * @brief Construct a new Veff object for Gamma-only calculation * @param GG_in: the pointer of Gint_Gamma object, used for grid integration */ Veff>(Gint_Gamma* GG_in, - HS_Matrix_K* hsk_in, - const std::vector>& kvec_d_in, - elecstate::Potential* pot_in, - hamilt::HContainer* hR_in, - const UnitCell* ucell_in, - const std::vector& orb_cutoff, - Grid_Driver* GridD_in) + HS_Matrix_K* hsk_in, + const std::vector>& kvec_d_in, + elecstate::Potential* pot_in, + hamilt::HContainer* hR_in, + const UnitCell* ucell_in, + const std::vector& orb_cutoff, + Grid_Driver* GridD_in, + const int& nspin) : GG(GG_in), orb_cutoff_(orb_cutoff), pot(pot_in), OperatorLCAO(hsk_in, kvec_d_in, hR_in) { this->cal_type = calculation_type::lcao_gint; this->initialize_HR(ucell_in, GridD_in); - GG_in->initialize_pvpR(*ucell_in, GridD_in); + GG_in->initialize_pvpR(*ucell_in, GridD_in, nspin); } ~Veff>(){}; diff --git a/source/module_hamilt_lcao/module_gint/gint.cpp b/source/module_hamilt_lcao/module_gint/gint.cpp index eb49d9f28f..a4dc96b422 100644 --- a/source/module_hamilt_lcao/module_gint/gint.cpp +++ b/source/module_hamilt_lcao/module_gint/gint.cpp @@ -133,15 +133,15 @@ void Gint::prep_grid(const Grid_Technique& gt, return; } -void Gint::initialize_pvpR(const UnitCell& ucell_in, Grid_Driver* gd) { +void Gint::initialize_pvpR(const UnitCell& ucell_in, Grid_Driver* gd, const int& nspin) { ModuleBase::TITLE("Gint", "initialize_pvpR"); int npol = 1; // there is the only resize code of DMRGint if (this->DMRGint.size() == 0) { - this->DMRGint.resize(PARAM.inp.nspin); + this->DMRGint.resize(nspin); } - if (PARAM.inp.nspin != 4) { + if (nspin != 4) { if (this->hRGint != nullptr) { delete this->hRGint; } @@ -153,7 +153,7 @@ void Gint::initialize_pvpR(const UnitCell& ucell_in, Grid_Driver* gd) { } this->hRGintCd = new hamilt::HContainer>(ucell_in.nat); - for (int is = 0; is < PARAM.inp.nspin; is++) { + for (int is = 0; is < nspin; is++) { if (this->DMRGint[is] != nullptr) { delete this->DMRGint[is]; } @@ -186,7 +186,7 @@ void Gint::initialize_pvpR(const UnitCell& ucell_in, Grid_Driver* gd) { } } - if (PARAM.globalv.gamma_only_local && PARAM.inp.nspin != 4) { + if (PARAM.globalv.gamma_only_local && nspin != 4) { this->hRGint->fix_gamma(); } for (int T1 = 0; T1 < ucell_in.ntype; ++T1) { diff --git a/source/module_hamilt_lcao/module_gint/gint.h b/source/module_hamilt_lcao/module_gint/gint.h index dddb11be88..cfd6211443 100644 --- a/source/module_hamilt_lcao/module_gint/gint.h +++ b/source/module_hamilt_lcao/module_gint/gint.h @@ -49,7 +49,7 @@ class Gint { * @brief calculate the neighbor atoms of each atom in this processor * size of BaseMatrix with be the non-parallel version */ - void initialize_pvpR(const UnitCell& unitcell, Grid_Driver* gd); + void initialize_pvpR(const UnitCell& unitcell, Grid_Driver* gd, const int& nspin); /** * @brief transfer DMR (2D para) to DMR (Grid para) in elecstate_lcao.cpp diff --git a/source/module_hamilt_lcao/module_gint/gint_rho_cpu_interface.cpp b/source/module_hamilt_lcao/module_gint/gint_rho_cpu_interface.cpp index 354de74857..160df1377b 100644 --- a/source/module_hamilt_lcao/module_gint/gint_rho_cpu_interface.cpp +++ b/source/module_hamilt_lcao/module_gint/gint_rho_cpu_interface.cpp @@ -54,7 +54,7 @@ void Gint::gint_kernel_rho(Gint_inout* inout) { cal_flag.get_ptr_2D(), psir_ylm.get_ptr_2D()); - for (int is = 0; is < PARAM.inp.nspin; ++is) + for (int is = 0; is < inout->nspin_rho; ++is) { ModuleBase::Array_Pool psir_DM(this->bxyz, LD_pool); ModuleBase::GlobalFunc::ZEROS(psir_DM.get_ptr_1D(), this->bxyz * LD_pool); diff --git a/source/module_hamilt_lcao/module_gint/gint_tools.h b/source/module_hamilt_lcao/module_gint/gint_tools.h index cfc0a2da0a..7466cac2a2 100644 --- a/source/module_hamilt_lcao/module_gint/gint_tools.h +++ b/source/module_hamilt_lcao/module_gint/gint_tools.h @@ -39,6 +39,7 @@ class Gint_inout bool isforce; bool isstress; int ispin; + int nspin_rho; // usually, but not always, equal to global nspin bool if_symm = false; // if true, use dsymv in gint_kernel_rho; if false, use dgemv. // output @@ -48,10 +49,11 @@ class Gint_inout Gint_Tools::job_type job; // electron density and kin_r, multi-k - Gint_inout(double** rho_in, Gint_Tools::job_type job_in, bool if_symm_in = true) + Gint_inout(double** rho_in, Gint_Tools::job_type job_in, const int& nspin_rho_in, bool if_symm_in = true) { rho = rho_in; job = job_in; + nspin_rho = nspin_rho_in; if_symm = if_symm_in; } @@ -110,15 +112,6 @@ class Gint_inout job = job_in; } - // electron density and kin_r, gamma point - Gint_inout(double*** DM_in, double** rho_in, Gint_Tools::job_type job_in, bool if_symm_in = true) - { - DM = DM_in; - rho = rho_in; - job = job_in; - if_symm = if_symm_in; - } - // vlocal, gamma point Gint_inout(const double* vl_in, Gint_Tools::job_type job_in) { diff --git a/source/module_io/get_pchg_lcao.cpp b/source/module_io/get_pchg_lcao.cpp index 7910f22584..33d1c4d645 100644 --- a/source/module_io/get_pchg_lcao.cpp +++ b/source/module_io/get_pchg_lcao.cpp @@ -104,9 +104,9 @@ void IState_Charge::begin(Gint_Gamma& gg, DM.init_DMR(GridD_in, ucell_in); DM.cal_DMR(); - gg.initialize_pvpR(*ucell_in, GridD_in); + gg.initialize_pvpR(*ucell_in, GridD_in, PARAM.inp.nspin); gg.transfer_DM2DtoGrid(DM.get_DMR_vector()); - Gint_inout inout((double***)nullptr, rho, Gint_Tools::job_type::rho); + Gint_inout inout(rho, Gint_Tools::job_type::rho, PARAM.inp.nspin); gg.cal_gint(&inout); // A solution to replace the original implementation of the following code: @@ -243,9 +243,9 @@ void IState_Charge::begin(Gint_k& gk, DM.init_DMR(GridD_in, ucell_in); DM.cal_DMR(ik); - gk.initialize_pvpR(*ucell_in, GridD_in); + gk.initialize_pvpR(*ucell_in, GridD_in, PARAM.inp.nspin); gk.transfer_DM2DtoGrid(DM.get_DMR_vector()); - Gint_inout inout(rho, Gint_Tools::job_type::rho); + Gint_inout inout(rho, Gint_Tools::job_type::rho, PARAM.inp.nspin); gk.cal_gint(&inout); // Using std::vector to replace the original double** rho_save @@ -298,9 +298,9 @@ void IState_Charge::begin(Gint_k& gk, DM.init_DMR(GridD_in, ucell_in); DM.cal_DMR(); - gk.initialize_pvpR(*ucell_in, GridD_in); + gk.initialize_pvpR(*ucell_in, GridD_in, PARAM.inp.nspin); gk.transfer_DM2DtoGrid(DM.get_DMR_vector()); - Gint_inout inout(rho, Gint_Tools::job_type::rho); + Gint_inout inout(rho, Gint_Tools::job_type::rho, PARAM.inp.nspin); gk.cal_gint(&inout); // Using std::vector to replace the original double** rho_save diff --git a/source/module_io/write_eband_terms.hpp b/source/module_io/write_eband_terms.hpp index 97ca13f2f9..41f0a42112 100644 --- a/source/module_io/write_eband_terms.hpp +++ b/source/module_io/write_eband_terms.hpp @@ -87,7 +87,7 @@ namespace ModuleIO hamilt::HContainer v_pp_local_R_ao(pv); if_gamma_fix(v_pp_local_R_ao); std::vector> e_orb_pp_local; - hamilt::Veff> v_pp_local_op(gint, &v_pp_local_k_ao, kv.kvec_d, &pot_local, &v_pp_local_R_ao, &ucell, orb_cutoff, &gd); + hamilt::Veff> v_pp_local_op(gint, &v_pp_local_k_ao, kv.kvec_d, &pot_local, &v_pp_local_R_ao, &ucell, orb_cutoff, &gd, nspin); v_pp_local_op.contributeHR(); for (int ik = 0;ik < kv.get_nks();++ik) { @@ -142,7 +142,7 @@ namespace ModuleIO for (int is = 0; is < nspin0; ++is) { v_hartree_op[is] = new hamilt::Veff>(gint, - &v_hartree_k_ao, kv.kvec_d, &pot_hartree, &v_hartree_R_ao[is], &ucell, orb_cutoff, &gd); + &v_hartree_k_ao, kv.kvec_d, &pot_hartree, &v_hartree_R_ao[is], &ucell, orb_cutoff, &gd, nspin); v_hartree_op[is]->contributeHR(); } std::vector> e_orb_hartree; diff --git a/source/module_io/write_vxc.hpp b/source/module_io/write_vxc.hpp index 2012f77742..d6fad87dad 100644 --- a/source/module_io/write_vxc.hpp +++ b/source/module_io/write_vxc.hpp @@ -238,13 +238,7 @@ void write_Vxc(const int nspin, for (int is = 0; is < nspin0; ++is) { vxcs_op_ao[is] = new hamilt::Veff>(gint, - &vxc_k_ao, - kv.kvec_d, - potxc, - &vxcs_R_ao[is], - &ucell, - orb_cutoff, - &gd); + &vxc_k_ao, kv.kvec_d, potxc, &vxcs_R_ao[is], &ucell, orb_cutoff, &gd, nspin); vxcs_op_ao[is]->contributeHR(); } diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 806c0b628a..9716ce35ac 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -372,7 +372,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu this->pw_rho->startz_current, &ucell, &orb); - this->gint_->initialize_pvpR(ucell, &GlobalC::GridD); + this->gint_->initialize_pvpR(ucell, &GlobalC::GridD, 1); // always use nspin=1 for transition density // if EXX from scratch, init 2-center integral and calculate Cs, Vs #ifdef __EXX diff --git a/source/module_lr/lr_spectrum.cpp b/source/module_lr/lr_spectrum.cpp index dcfefee4ea..78d11d7330 100644 --- a/source/module_lr/lr_spectrum.cpp +++ b/source/module_lr/lr_spectrum.cpp @@ -12,7 +12,7 @@ void LR::LR_Spectrum::cal_gint_rho(double** rho, const int& nspin_solve, cons for (int is = 0;is < nspin_solve;++is) { ModuleBase::GlobalFunc::ZEROS(rho[is], nrxx); } - Gint_inout inout_rho(rho, Gint_Tools::job_type::rho, false); + Gint_inout inout_rho(rho, Gint_Tools::job_type::rho, 1, false); this->gint->cal_gint(&inout_rho); } diff --git a/source/module_lr/operator_casida/operator_lr_hxc.cpp b/source/module_lr/operator_casida/operator_lr_hxc.cpp index d3fffff390..b574ac1f8d 100644 --- a/source/module_lr/operator_casida/operator_lr_hxc.cpp +++ b/source/module_lr/operator_casida/operator_lr_hxc.cpp @@ -101,9 +101,9 @@ namespace LR double** rho_trans; const int& nrxx = this->pot.lock()->nrxx; // LR_Util::new_p2(rho_trans, nspin_solve, nrxx); - LR_Util::new_p2(rho_trans, nspin, nrxx); // currently gint_kernel_rho uses PARAM.inp.nspin, it needs refactor + LR_Util::new_p2(rho_trans, 1, nrxx); // currently gint_kernel_rho uses PARAM.inp.nspin, it needs refactor for (int is = 0;is < nspin_solve;++is)ModuleBase::GlobalFunc::ZEROS(rho_trans[is], nrxx); - Gint_inout inout_rho(rho_trans, Gint_Tools::job_type::rho, false); + Gint_inout inout_rho(rho_trans, Gint_Tools::job_type::rho, 1, false); this->gint->cal_gint(&inout_rho); // 3. v_hxc = f_hxc * rho_trans @@ -150,7 +150,7 @@ namespace LR // LR_Util::new_p2(rho_trans, nspin_solve, nrxx); LR_Util::new_p2(rho_trans, nspin, nrxx); // currently gint_kernel_rho uses PARAM.inp.nspin, it needs refactor for (int is = 0;is < nspin_solve;++is)ModuleBase::GlobalFunc::ZEROS(rho_trans[is], nrxx); - Gint_inout inout_rho(rho_trans, Gint_Tools::job_type::rho, false); + Gint_inout inout_rho(rho_trans, Gint_Tools::job_type::rho, 1, false); this->gint->cal_gint(&inout_rho); // print_grid_nonzero(rho_trans[0], nrxx, 10, "rho_trans");