diff --git a/source/driver.cpp b/source/driver.cpp index 40b314b543..67476ddd5c 100644 --- a/source/driver.cpp +++ b/source/driver.cpp @@ -19,9 +19,6 @@ Driver::Driver() Driver::~Driver() { - // Release the device memory within singleton object GlobalC::ppcell - // before the main function exits. - GlobalC::ppcell.release_memory(); } void Driver::init() diff --git a/source/module_elecstate/elecstate.cpp b/source/module_elecstate/elecstate.cpp index 231d1c7728..c84ce504d2 100644 --- a/source/module_elecstate/elecstate.cpp +++ b/source/module_elecstate/elecstate.cpp @@ -208,13 +208,14 @@ void ElecState::calEBand() void ElecState::init_scf(const int istep, const ModuleBase::ComplexMatrix& strucfac, + const bool* numeric, ModuleSymmetry::Symmetry& symm, const void* wfcpw) { //! core correction potential. if (!PARAM.inp.use_paw) { - this->charge->set_rho_core(strucfac); + this->charge->set_rho_core(strucfac, numeric); } else { diff --git a/source/module_elecstate/elecstate.h b/source/module_elecstate/elecstate.h index 624ef355f5..a90555a249 100644 --- a/source/module_elecstate/elecstate.h +++ b/source/module_elecstate/elecstate.h @@ -110,6 +110,7 @@ class ElecState */ void init_scf(const int istep, const ModuleBase::ComplexMatrix& strucfac, + const bool* numeric, ModuleSymmetry::Symmetry& symm, const void* wfcpw = nullptr); std::string classname = "elecstate"; diff --git a/source/module_elecstate/module_charge/charge.h b/source/module_elecstate/module_charge/charge.h index 1593c9fa47..bc8a1aa4ec 100644 --- a/source/module_elecstate/module_charge/charge.h +++ b/source/module_elecstate/module_charge/charge.h @@ -84,7 +84,7 @@ class Charge const ModuleBase::ComplexMatrix& strucFac, const UnitCell& ucell) const; - void set_rho_core(const ModuleBase::ComplexMatrix &structure_factor); + void set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, const bool* numeric); void set_rho_core_paw(); void renormalize_rho(); diff --git a/source/module_elecstate/module_charge/charge_init.cpp b/source/module_elecstate/module_charge/charge_init.cpp index b41db4047f..b38c5cdd5c 100644 --- a/source/module_elecstate/module_charge/charge_init.cpp +++ b/source/module_elecstate/module_charge/charge_init.cpp @@ -199,9 +199,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, //========================================================== // computes the core charge on the real space 3D mesh. //========================================================== -void Charge::set_rho_core( - const ModuleBase::ComplexMatrix &structure_factor -) +void Charge::set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, const bool* numeric) { ModuleBase::TITLE("Charge","set_rho_core"); ModuleBase::timer::tick("Charge","set_rho_core"); @@ -249,7 +247,7 @@ void Charge::set_rho_core( // each shell of g vec //---------------------------------------------------------- this->non_linear_core_correction( - GlobalC::ppcell.numeric, + numeric, GlobalC::ucell.atoms[it].ncpp.msh, GlobalC::ucell.atoms[it].ncpp.r.data(), GlobalC::ucell.atoms[it].ncpp.rab.data(), diff --git a/source/module_elecstate/test/elecstate_base_test.cpp b/source/module_elecstate/test/elecstate_base_test.cpp index 19402a173f..050a50b74a 100644 --- a/source/module_elecstate/test/elecstate_base_test.cpp +++ b/source/module_elecstate/test/elecstate_base_test.cpp @@ -60,7 +60,7 @@ void ModulePW::PW_Basis::initgrids(double, ModuleBase::Matrix3, int, int, int) void ModulePW::PW_Basis::distribute_r() { } -void Charge::set_rho_core(ModuleBase::ComplexMatrix const&) +void Charge::set_rho_core(ModuleBase::ComplexMatrix const&, const bool*) { } void Charge::set_rho_core_paw() @@ -249,7 +249,7 @@ TEST_F(ElecStateTest, InitSCF) ModuleBase::ComplexMatrix strucfac; elecstate->eferm = efermi; ModuleSymmetry::Symmetry symm; - EXPECT_NO_THROW(elecstate->init_scf(istep, strucfac, symm)); + EXPECT_NO_THROW(elecstate->init_scf(istep, strucfac, nullptr, symm)); // delete elecstate->pot is done in the destructor of elecstate delete charge; } diff --git a/source/module_elecstate/test/elecstate_pw_test.cpp b/source/module_elecstate/test/elecstate_pw_test.cpp index b491da01a0..c8692a01e0 100644 --- a/source/module_elecstate/test/elecstate_pw_test.cpp +++ b/source/module_elecstate/test/elecstate_pw_test.cpp @@ -139,7 +139,7 @@ K_Vectors::~K_Vectors() { } -void Charge::set_rho_core(ModuleBase::ComplexMatrix const&) +void Charge::set_rho_core(ModuleBase::ComplexMatrix const&, const bool*) { } void Charge::set_rho_core_paw() diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index a336f0b864..776e837d9a 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -281,7 +281,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) this->pw_rhod->collect_uniqgg(); } - GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, pw_rhod); + this->p_locpp->init_vloc(this->pw_rhod); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); this->pelec->omega = ucell.omega; diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h index 708fda3c35..fc66cb42a4 100644 --- a/source/module_esolver/esolver_fp.h +++ b/source/module_esolver/esolver_fp.h @@ -53,6 +53,9 @@ namespace ModuleESolver ModulePW::PW_Basis* pw_rho; + //! pointer to pseudopotential + pseudopot_cell_vl* p_locpp = nullptr; + /** * @brief same as pw_rho for ncpp. Here 'd' stands for 'dense' * dense grid for for uspp, used for ultrasoft augmented charge density. diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 53136de819..d8fa460ba1 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -68,6 +68,8 @@ ESolver_KS::ESolver_KS() ///---------------------------------------------------------- p_chgmix = new Charge_Mixing(); p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod); + this->ppcell.cell_factor = PARAM.inp.cell_factor; + this->p_locpp = &this->ppcell; } //------------------------------------------------------------------------------ @@ -81,6 +83,7 @@ ESolver_KS::~ESolver_KS() delete this->pw_wfc; delete this->p_hamilt; delete this->p_chgmix; + this->ppcell.release_memory(); } //------------------------------------------------------------------------------ diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 892b948276..9c9b53c4af 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -66,6 +66,9 @@ class ESolver_KS : public ESolver_FP //! Charge mixing method, only used in KDSFT, not in OFDFT Charge_Mixing* p_chgmix = nullptr; + //! nonlocal pseudo potential + pseudopot_cell_vnl ppcell; + //! Electronic wavefunctions psi::Psi* psi = nullptr; diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 437e634d21..a6fd6dc0ab 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -197,7 +197,8 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa } // 8) initialize ppcell - GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rho); + this->ppcell.init_vloc(this->pw_rho); + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); // 9) inititlize the charge density this->pelec->charge->allocate(PARAM.inp.nspin); @@ -209,7 +210,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho, &ucell, - &(GlobalC::ppcell.vloc), + &(this->ppcell.vloc), &(this->sf), &(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc)); @@ -302,6 +303,7 @@ void ESolver_KS_LCAO::cal_force(UnitCell& ucell, ModuleBase::matrix& for orb_, force, this->scs, + this->ppcell, this->sf, this->kv, this->pw_rho, @@ -447,7 +449,7 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) this->sf, *this->pw_rho, *this->pw_rhod, - GlobalC::ppcell.vloc, + this->ppcell.vloc, *this->pelec->charge, this->GG, this->GK, @@ -474,7 +476,7 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) this->sf, *this->pw_rho, *this->pw_rhod, - GlobalC::ppcell.vloc, + this->ppcell.vloc, *this->pelec->charge, this->GG, this->GK, diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp index eb54bd9b66..45bcac2785 100644 --- a/source/module_esolver/esolver_ks_lcaopw.cpp +++ b/source/module_esolver/esolver_ks_lcaopw.cpp @@ -64,7 +64,7 @@ namespace ModuleESolver template void ESolver_KS_LIP::allocate_hamilt() { - this->p_hamilt = new hamilt::HamiltLIP(this->pelec->pot, this->pw_wfc, &this->kv + this->p_hamilt = new hamilt::HamiltLIP(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell #ifdef __EXX , *this->exx_lip #endif @@ -250,7 +250,7 @@ namespace ModuleESolver *this->pw_wfc, *this->pw_rho, *this->pw_rhod, - GlobalC::ppcell.vloc, + this->ppcell.vloc, *this->pelec->charge, this->kv, this->pelec->wg diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index f081f5596d..3d6a993bb2 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -117,7 +117,7 @@ ESolver_KS_PW::~ESolver_KS_PW() template void ESolver_KS_PW::allocate_hamilt() { - this->p_hamilt = new hamilt::HamiltPW(this->pelec->pot, this->pw_wfc, &this->kv); + this->p_hamilt = new hamilt::HamiltPW(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell); } template void ESolver_KS_PW::deallocate_hamilt() @@ -145,7 +145,7 @@ void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_p &(this->chr), &(this->kv), &ucell, - &(GlobalC::ppcell), + &(this->ppcell), this->pw_rhod, this->pw_rho, this->pw_big); @@ -156,7 +156,7 @@ void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_p &(this->chr), &(this->kv), &ucell, - &GlobalC::ppcell, + &this->ppcell, this->pw_rhod, this->pw_rho, this->pw_big); @@ -175,7 +175,7 @@ void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_p this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho, &ucell, - &GlobalC::ppcell.vloc, + &this->ppcell.vloc, &(this->sf), &(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc)); @@ -194,7 +194,7 @@ void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_p &GlobalC::Pkpoints, GlobalV::MY_RANK, #endif - &GlobalC::ppcell); + &this->ppcell); if (this->psi != nullptr) { @@ -202,14 +202,14 @@ void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_p } //! init pseudopotential - GlobalC::ppcell.init(ucell.ntype, &this->sf, this->pw_wfc); + this->ppcell.init(ucell.ntype, &this->sf, this->pw_wfc); //! initalize local pseudopotential - GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rhod); + this->ppcell.init_vloc(this->pw_rhod); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); //! Initalize non-local pseudopotential - GlobalC::ppcell.init_vnl(ucell, this->pw_rhod); + this->ppcell.init_vnl(ucell, this->pw_rhod); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL"); //! Allocate psi @@ -218,7 +218,8 @@ void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_p this->kv.get_nks(), this->kv.ngk.data(), this->pw_wfc->npwk_max, - &(this->sf)); + &this->sf, + &this->ppcell); this->kspw_psi = PARAM.inp.device == "gpu" || PARAM.inp.precision == "single" ? new psi::Psi(this->psi[0]) @@ -247,7 +248,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) if (ucell.cell_parameter_updated) { - GlobalC::ppcell.init_vnl(ucell, this->pw_rhod); + this->ppcell.init_vnl(ucell, this->pw_rhod); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL"); this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, this->pw_wfc->nx, this->pw_wfc->ny, this->pw_wfc->nz); @@ -256,7 +257,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) this->pw_wfc->collect_local_pw(PARAM.inp.erf_ecut, PARAM.inp.erf_height, PARAM.inp.erf_sigma); - this->p_wf_init->make_table(this->kv.get_nks(), &this->sf); + this->p_wf_init->make_table(this->kv.get_nks(), &this->sf, &this->ppcell); } if (ucell.ionic_position_updated) { @@ -303,7 +304,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) } //! calculate the total local pseudopotential in real space - this->pelec->init_scf(istep, this->sf.strucFac, ucell.symm, (void*)this->pw_wfc); + this->pelec->init_scf(istep, this->sf.strucFac, this->ppcell.numeric, ucell.symm, (void*)this->pw_wfc); //! output the initial charge density if (PARAM.inp.out_chg[0] == 2) @@ -358,7 +359,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) // projectors ModuleBase::matrix veff = this->pelec->pot->get_effective_v(); - GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, ucell); + this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell); // after init_rho (in pelec->init_scf), we have rho now. // before hamilt2density, we update Hk and initialize psi @@ -374,6 +375,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) this->p_wf_init->initialize_psi(this->psi, this->kspw_psi, this->p_hamilt, + this->ppcell, GlobalV::ofs_running, this->already_initpsi); @@ -500,7 +502,7 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int if (PARAM.globalv.use_uspp) { ModuleBase::matrix veff = this->pelec->pot->get_effective_v(); - GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, ucell); + this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell); } if (this->out_freq_elec && iter % this->out_freq_elec == 0) @@ -623,7 +625,16 @@ void ESolver_KS_PW::cal_force(UnitCell& ucell, ModuleBase::matrix& fo : reinterpret_cast, Device>*>(this->kspw_psi); // Calculate forces - ff.cal_force(force, *this->pelec, this->pw_rhod, &ucell.symm, &this->sf, &this->kv, this->pw_wfc, this->__kspw_psi); + ff.cal_force(force, + *this->pelec, + this->pw_rhod, + &ucell.symm, + &this->sf, + &this->ppcell, + &this->ppcell, + &this->kv, + this->pw_wfc, + this->__kspw_psi); } template @@ -641,7 +652,7 @@ void ESolver_KS_PW::cal_stress(UnitCell& ucell, ModuleBase::matrix& s : reinterpret_cast, Device>*>(this->kspw_psi); ss.cal_stress(stress, ucell, - &GlobalC::ppcell, + this->ppcell, this->pw_rhod, &ucell.symm, &this->sf, @@ -781,7 +792,7 @@ void ESolver_KS_PW::after_all_runners(UnitCell& ucell) //! 7) Use Kubo-Greenwood method to compute conductivities if (PARAM.inp.cal_cond) { - EleCond elec_cond(&ucell, &this->kv, this->pelec, this->pw_wfc, this->psi, &GlobalC::ppcell); + EleCond elec_cond(&ucell, &this->kv, this->pelec, this->pw_wfc, this->psi, &this->ppcell); elec_cond.KG(PARAM.inp.cond_smear, PARAM.inp.cond_fwhm, PARAM.inp.cond_wcut, diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index a5a68ba14a..210b5e5926 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -23,6 +23,7 @@ ESolver_OF::ESolver_OF() { this->classname = "ESolver_OF"; this->task_ = new char[60]; + this->p_locpp = &this->locpp; } ESolver_OF::~ESolver_OF() @@ -111,24 +112,20 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp) ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS"); // initialize local pseudopotential - GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, pw_rho); + this->locpp.init_vloc(pw_rho); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); - // initialize non local pseudopotential - GlobalC::ppcell.init_vnl(ucell, pw_rho); - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL"); // initialize elecstate, including potential this->init_elecstate(ucell); // calculate the total local pseudopotential in real space - this->pelec->init_scf(0, sf.strucFac, ucell.symm); // atomic_rho, v_of_rho, set_vrs + this->pelec->init_scf(0, sf.strucFac, locpp.numeric, ucell.symm); // atomic_rho, v_of_rho, set_vrs // liuyu move here 2023-10-09 // D in uspp need vloc, thus behind init_scf() // calculate the effective coefficient matrix for non-local pseudopotential projectors ModuleBase::matrix veff = this->pelec->pot->get_effective_v(); - GlobalC::ppcell.cal_effective_D(veff, this->pw_rho, ucell); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT POTENTIAL"); @@ -212,7 +209,6 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) // initialize elecstate, including potential this->init_elecstate(ucell); - GlobalC::ppcell.init_vnl(ucell, pw_rho); // Initialize KEDF this->init_kedf(PARAM.inp); @@ -259,7 +255,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) GlobalV::ofs_warning); } - this->pelec->init_scf(istep, sf.strucFac, ucell.symm); + this->pelec->init_scf(istep, sf.strucFac, locpp.numeric, ucell.symm); // calculate ewald energy this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rho, sf.strucFac); @@ -545,7 +541,7 @@ double ESolver_OF::cal_energy() void ESolver_OF::cal_force(UnitCell& ucell, ModuleBase::matrix& force) { Forces ff(ucell.nat); - ff.cal_force(force, *pelec, this->pw_rho, &ucell.symm, &sf); + ff.cal_force(force, *pelec, this->pw_rho, &ucell.symm, &sf, &this->locpp); } /** @@ -560,6 +556,6 @@ void ESolver_OF::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) this->kinetic_stress(kinetic_stress_); OF_Stress_PW ss(this->pelec, this->pw_rho); - ss.cal_stress(stress, kinetic_stress_, ucell, &ucell.symm, &sf, &kv); + ss.cal_stress(stress, kinetic_stress_, ucell, &ucell.symm, this->locpp, &sf, &kv); } } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_of.h b/source/module_esolver/esolver_of.h index 081a077606..51e014886f 100644 --- a/source/module_esolver/esolver_of.h +++ b/source/module_esolver/esolver_of.h @@ -30,6 +30,8 @@ class ESolver_OF : public ESolver_FP virtual void cal_force(UnitCell& ucell, ModuleBase::matrix& force) override; virtual void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) override; + protected: + pseudopot_cell_vl locpp; private: // ======================= variables ========================== diff --git a/source/module_esolver/esolver_of_tool.cpp b/source/module_esolver/esolver_of_tool.cpp index 85c6b640fe..4fc93e005d 100644 --- a/source/module_esolver/esolver_of_tool.cpp +++ b/source/module_esolver/esolver_of_tool.cpp @@ -28,7 +28,7 @@ void ESolver_OF::init_elecstate(UnitCell& ucell) this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho, &ucell, - &(GlobalC::ppcell.vloc), + &(this->locpp.vloc), &(this->sf), &(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc)); diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index 9cb35acbc3..76b510a36a 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -99,6 +99,7 @@ void ESolver_SDFT_PW::before_scf(UnitCell& ucell, const int istep) this->p_hamilt = new hamilt::HamiltSdftPW(this->pelec->pot, this->pw_wfc, &this->kv, + &this->ppcell, PARAM.globalv.npol, &this->stoche.emin_sto, &this->stoche.emax_sto); @@ -220,7 +221,7 @@ void ESolver_SDFT_PW::cal_force(UnitCell& ucell, ModuleBase::matrix& &this->sf, &this->kv, this->pw_wfc, - GlobalC::ppcell, + this->ppcell, ucell, *this->kspw_psi, this->stowf); @@ -240,7 +241,7 @@ void ESolver_SDFT_PW::cal_stress(UnitCell& ucell, ModuleBase::matrix& *this->kspw_psi, this->stowf, this->pelec->charge, - &GlobalC::ppcell, + &this->ppcell, ucell); } @@ -290,7 +291,7 @@ void ESolver_SDFT_PW, base_device::DEVICE_CPU>::after_all_r this->pelec, this->pw_wfc, this->psi, - &GlobalC::ppcell, + &this->ppcell, this->p_hamilt, this->stoche, &stowf); diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index aa5d623dbc..b7ec3ce705 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -265,7 +265,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) } #endif // __EXX - this->pelec->init_scf(istep, this->sf.strucFac, ucell.symm); + this->pelec->init_scf(istep, this->sf.strucFac, this->ppcell.numeric, ucell.symm); //! output the initial charge density if (PARAM.inp.out_chg[0] == 2) @@ -375,7 +375,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) if( PARAM.inp.rdmft == true ) { // necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp - rdmft_solver.update_ion(ucell, *(this->pw_rho), GlobalC::ppcell.vloc, this->sf.strucFac); // add by jghan, 2024-03-16/2024-10-08 + rdmft_solver.update_ion(ucell, *(this->pw_rho), this->ppcell.vloc, this->sf.strucFac); // add by jghan, 2024-03-16/2024-10-08 } return; diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 75f4164e2f..10d688dc00 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -252,7 +252,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) } // pelec should be initialized before these calculations - this->pelec->init_scf(istep, this->sf.strucFac, ucell.symm); + this->pelec->init_scf(istep, this->sf.strucFac, this->ppcell.numeric, ucell.symm); // self consistent calculations for electronic ground state if (cal_type == "get_pchg") { diff --git a/source/module_hamilt_general/module_surchem/sol_force.cpp b/source/module_hamilt_general/module_surchem/sol_force.cpp index 68ba6613c5..80a30040dd 100644 --- a/source/module_hamilt_general/module_surchem/sol_force.cpp +++ b/source/module_hamilt_general/module_surchem/sol_force.cpp @@ -2,7 +2,10 @@ #include "module_base/timer.h" #include "module_parameter/parameter.h" -void force_cor_one(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, ModuleBase::matrix& forcesol) +void force_cor_one(const UnitCell& cell, + const ModulePW::PW_Basis* rho_basis, + const ModuleBase::matrix& vloc, + ModuleBase::matrix& forcesol) { @@ -26,7 +29,7 @@ void force_cor_one(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, Mo { complex phase = exp( ModuleBase::NEG_IMAG_UNIT *ModuleBase::TWO_PI * ( rho_basis->gcar[ig] * cell.atoms[it].tau[ia])); //vloc for each atom - vloc_at[ig] = GlobalC::ppcell.vloc(it, rho_basis->ig2igg[ig]) * phase; + vloc_at[ig] = vloc(it, rho_basis->ig2igg[ig]) * phase; if(rho_basis->ig_gge0 == ig) { N[ig] = GlobalC::ucell.atoms[it].ncpp.zv / GlobalC::ucell.omega; @@ -146,7 +149,10 @@ void force_cor_two(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, Mo } -void surchem::cal_force_sol(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, ModuleBase::matrix& forcesol) +void surchem::cal_force_sol(const UnitCell& cell, + const ModulePW::PW_Basis* rho_basis, + const ModuleBase::matrix& vloc, + ModuleBase::matrix& forcesol) { ModuleBase::TITLE("surchem", "cal_force_sol"); ModuleBase::timer::tick("surchem", "cal_force_sol"); @@ -155,7 +161,7 @@ void surchem::cal_force_sol(const UnitCell& cell, const ModulePW::PW_Basis* rho_ ModuleBase::matrix force1(nat, 3); ModuleBase::matrix force2(nat, 3); - force_cor_one(cell, rho_basis,force1); + force_cor_one(cell, rho_basis, vloc, force1); force_cor_two(cell, rho_basis,force2); int iat = 0; diff --git a/source/module_hamilt_general/module_surchem/surchem.h b/source/module_hamilt_general/module_surchem/surchem.h index 5696157345..1a6add9e8c 100644 --- a/source/module_hamilt_general/module_surchem/surchem.h +++ b/source/module_hamilt_general/module_surchem/surchem.h @@ -106,7 +106,10 @@ class surchem const ModulePW::PW_Basis* rho_basis, const double* const* const rho); - void cal_force_sol(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, ModuleBase::matrix& forcesol); + void cal_force_sol(const UnitCell& cell, + const ModulePW::PW_Basis* rho_basis, + const ModuleBase::matrix& vloc, + ModuleBase::matrix& forcesol); void get_totn_reci(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, complex* totn_reci); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index 97e34a78b3..7423ade290 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -44,6 +44,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, const LCAO_Orbitals& orb, ModuleBase::matrix& fcs, ModuleBase::matrix& scs, + const pseudopot_cell_vnl& nlpp, const Structure_Factor& sf, const K_Vectors& kv, ModulePW::PW_Basis* rhopw, @@ -101,6 +102,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, pelec->vnew_exist, pelec->charge, rhopw, + nlpp, sf); } @@ -144,6 +146,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, pelec->f_en.etxc, pelec->charge, rhopw, + nlpp, sf); } @@ -274,7 +277,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, if (PARAM.inp.imp_sol && isforce) { fsol.create(nat, 3); - GlobalC::solvent_model.cal_force_sol(GlobalC::ucell, rhopw, fsol); + GlobalC::solvent_model.cal_force_sol(GlobalC::ucell, rhopw, nlpp.vloc, fsol); } //! atomic forces from DFT+U (Quxin version) @@ -852,6 +855,7 @@ void Force_Stress_LCAO::calForcePwPart(ModuleBase::matrix& fvl_dvl, const bool vnew_exist, const Charge* const chr, ModulePW::PW_Basis* rhopw, + const pseudopot_cell_vnl& nlpp, const Structure_Factor& sf) { ModuleBase::TITLE("Force_Stress_LCAO", "calForcePwPart"); @@ -859,7 +863,7 @@ void Force_Stress_LCAO::calForcePwPart(ModuleBase::matrix& fvl_dvl, // local pseudopotential force: // use charge density; plane wave; local pseudopotential; //-------------------------------------------------------- - f_pw.cal_force_loc(fvl_dvl, rhopw, chr); + f_pw.cal_force_loc(fvl_dvl, rhopw, nlpp.vloc, chr); //-------------------------------------------------------- // ewald force: use plane wave only. //-------------------------------------------------------- @@ -868,11 +872,11 @@ void Force_Stress_LCAO::calForcePwPart(ModuleBase::matrix& fvl_dvl, //-------------------------------------------------------- // force due to core correlation. //-------------------------------------------------------- - f_pw.cal_force_cc(fcc, rhopw, chr,GlobalC::ucell); + f_pw.cal_force_cc(fcc, rhopw, chr, nlpp.numeric, GlobalC::ucell); //-------------------------------------------------------- // force due to self-consistent charge. //-------------------------------------------------------- - f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist,GlobalC::ucell); + f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist, nlpp.numeric, GlobalC::ucell); return; } @@ -988,6 +992,7 @@ void Force_Stress_LCAO::calStressPwPart(ModuleBase::matrix& sigmadvl, const double& etxc, const Charge* const chr, ModulePW::PW_Basis* rhopw, + const pseudopot_cell_vnl& nlpp, const Structure_Factor& sf) { ModuleBase::TITLE("Force_Stress_LCAO", "calStressPwPart"); @@ -995,7 +1000,7 @@ void Force_Stress_LCAO::calStressPwPart(ModuleBase::matrix& sigmadvl, // local pseudopotential stress: // use charge density; plane wave; local pseudopotential; //-------------------------------------------------------- - sc_pw.stress_loc(sigmadvl, rhopw, &sf, 0, chr); + sc_pw.stress_loc(sigmadvl, rhopw, nlpp.vloc, &sf, 0, chr); //-------------------------------------------------------- // hartree term @@ -1010,7 +1015,7 @@ void Force_Stress_LCAO::calStressPwPart(ModuleBase::matrix& sigmadvl, //-------------------------------------------------------- // stress due to core correlation. //-------------------------------------------------------- - sc_pw.stress_cc(sigmacc, rhopw, &sf, 0, chr); + sc_pw.stress_cc(sigmacc, rhopw, &sf, 0, nlpp.numeric, chr); //-------------------------------------------------------- // stress due to self-consistent charge. diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h index cfb1dd87fa..35a596bdcb 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h @@ -42,6 +42,7 @@ class Force_Stress_LCAO const LCAO_Orbitals& orb, ModuleBase::matrix& fcs, ModuleBase::matrix& scs, + const pseudopot_cell_vnl& nlpp, const Structure_Factor& sf, const K_Vectors& kv, ModulePW::PW_Basis* rhopw, @@ -69,6 +70,7 @@ class Force_Stress_LCAO const bool vnew_exist, const Charge* const chr, ModulePW::PW_Basis* rhopw, + const pseudopot_cell_vnl& nlpp, const Structure_Factor& sf); void integral_part(const bool isGammaOnly, @@ -103,6 +105,7 @@ class Force_Stress_LCAO const double& etxc, const Charge* const chr, ModulePW::PW_Basis* rhopw, + const pseudopot_cell_vnl& nlpp, const Structure_Factor& sf); static double force_invalid_threshold_ev; diff --git a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp index 7d196e8f05..1f28939c09 100644 --- a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp @@ -11,6 +11,7 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, ModuleBase::matrix& kinetic_stress, UnitCell& ucell, ModuleSymmetry::Symmetry* p_symm, + const pseudopot_cell_vl& locpp, Structure_Factor* p_sf, K_Vectors* p_kv) { @@ -75,10 +76,10 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, stress_gga(sigmaxc, this->rhopw, pelec->charge); // local contribution - stress_loc(sigmaloc, this->rhopw, p_sf, true, pelec->charge); + stress_loc(sigmaloc, this->rhopw, locpp.vloc, p_sf, true, pelec->charge); // nlcc - stress_cc(sigmaxcc, this->rhopw, p_sf, true, pelec->charge); + stress_cc(sigmaxcc, this->rhopw, p_sf, true, locpp.numeric, pelec->charge); // vdw term stress_vdw(sigmavdw, ucell); diff --git a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.h b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.h index 81402561cd..05265915c8 100644 --- a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.h +++ b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.h @@ -15,6 +15,7 @@ class OF_Stress_PW : public Stress_Func ModuleBase::matrix& kinetic_stress, UnitCell& ucell, ModuleSymmetry::Symmetry* p_symm, + const pseudopot_cell_vl& locpp, Structure_Factor* p_sf, K_Vectors* p_kv); diff --git a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp index 2516dadfe4..98f4198036 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp @@ -17,9 +17,10 @@ pseudopot_cell_vl::~pseudopot_cell_vl() delete[] zp; } -void pseudopot_cell_vl::init_vloc(ModuleBase::matrix& vloc_in, const ModulePW::PW_Basis* rho_basis) +void pseudopot_cell_vl::init_vloc(const ModulePW::PW_Basis* rho_basis) { - if(PARAM.inp.use_paw) return; + if(PARAM.inp.use_paw) { return; +} ModuleBase::TITLE("pseudopot_cell_vl","init_vloc"); // This routine computes the fourier coefficient of the local @@ -59,9 +60,9 @@ void pseudopot_cell_vl::init_vloc(ModuleBase::matrix& vloc_in, const ModulePW::P ModuleBase::WARNING_QUIT("init_vloc","not available now."); } - if(it>=0 && it=0) + if(it>=0 && itvloc.nr && this->vloc.nc>=0) { - ModuleBase::GlobalFunc::COPYARRAY(vloc1d, &vloc_in(it, 0), rho_basis->ngg); + ModuleBase::GlobalFunc::COPYARRAY(vloc1d, &this->vloc(it, 0), rho_basis->ngg); } } @@ -77,8 +78,10 @@ void pseudopot_cell_vl::init_vloc(ModuleBase::matrix& vloc_in, const ModulePW::P void pseudopot_cell_vl::allocate(const int ngg) { - if(PARAM.inp.test_pp>0) ModuleBase::TITLE("pseudopot_cell_vl","allocate"); - if(PARAM.inp.use_paw) return; + if(PARAM.inp.test_pp>0) { ModuleBase::TITLE("pseudopot_cell_vl","allocate"); +} + if(PARAM.inp.use_paw) { return; +} this->vloc.create(GlobalC::ucell.ntype, ngg); delete[] numeric; @@ -252,7 +255,8 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh, void pseudopot_cell_vl::print_vloc(const ModulePW::PW_Basis* rho_basis) const { - if(GlobalV::MY_RANK!=0) return; //mohan fix bug 2011-10-13 + if(GlobalV::MY_RANK!=0) { return; //mohan fix bug 2011-10-13 +} bool check_vl = PARAM.inp.out_element_info; if(check_vl) { diff --git a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.h b/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.h index b9042f3223..a5f91b0f81 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.h +++ b/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.h @@ -13,7 +13,13 @@ class pseudopot_cell_vl pseudopot_cell_vl(); ~pseudopot_cell_vl(); - void init_vloc(ModuleBase::matrix& vloc_in, const ModulePW::PW_Basis* rho_basis); + /** + * @brief init local potential + * + * @param rho_basis pw basis + * @return this->vloc + */ + void init_vloc(const ModulePW::PW_Basis* rho_basis); ModuleBase::matrix vloc; //(ntype,ngl),the local potential for each atom type(ntype,ngl) bool *numeric; //[ntype], =true diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp index 41fb91da6f..0ffce5313a 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp @@ -29,6 +29,8 @@ void Forces::cal_force(ModuleBase::matrix& force, ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry* p_symm, Structure_Factor* p_sf, + const pseudopot_cell_vl* locpp, + const pseudopot_cell_vnl* p_nlpp, K_Vectors* pkv, ModulePW::PW_Basis_K* wfc_basis, const psi::Psi, Device>* psi_in) @@ -52,7 +54,7 @@ void Forces::cal_force(ModuleBase::matrix& force, // For PAW, calculated together in paw_cell.calculate_force if (!PARAM.inp.use_paw) { - this->cal_force_loc(forcelc, rho_basis, chr); + this->cal_force_loc(forcelc, rho_basis, locpp->vloc, chr); } else { @@ -68,11 +70,11 @@ void Forces::cal_force(ModuleBase::matrix& force, if (!PARAM.inp.use_paw) { this->npwx = wfc_basis->npwk_max; - Forces::cal_force_nl(forcenl, wg, ekb, pkv, wfc_basis, p_sf, &GlobalC::ppcell, GlobalC::ucell, psi_in); + Forces::cal_force_nl(forcenl, wg, ekb, pkv, wfc_basis, p_sf, *p_nlpp, GlobalC::ucell, psi_in); if (PARAM.globalv.use_uspp) { - this->cal_force_us(forcenl, rho_basis, &GlobalC::ppcell, elec, GlobalC::ucell); + this->cal_force_us(forcenl, rho_basis, *p_nlpp, elec, GlobalC::ucell); } } else @@ -159,7 +161,7 @@ void Forces::cal_force(ModuleBase::matrix& force, // not relevant for PAW if (!PARAM.inp.use_paw) { - Forces::cal_force_cc(forcecc, rho_basis, chr,GlobalC::ucell); + Forces::cal_force_cc(forcecc, rho_basis, chr, locpp->numeric, GlobalC::ucell); } else { @@ -170,7 +172,7 @@ void Forces::cal_force(ModuleBase::matrix& force, // For PAW, calculated together in paw_cell.calculate_force if (!PARAM.inp.use_paw) { - this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist,GlobalC::ucell); + this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, locpp->numeric, GlobalC::ucell); } else { @@ -222,7 +224,7 @@ void Forces::cal_force(ModuleBase::matrix& force, if (PARAM.inp.imp_sol) { forcesol.create(this->nat, 3); - GlobalC::solvent_model.cal_force_sol(GlobalC::ucell, rho_basis, forcesol); + GlobalC::solvent_model.cal_force_sol(GlobalC::ucell, rho_basis, locpp->vloc, forcesol); if (PARAM.inp.test_force) { ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "IMP_SOL FORCE (Ry/Bohr)", forcesol); @@ -463,6 +465,7 @@ void Forces::cal_force(ModuleBase::matrix& force, template void Forces::cal_force_loc(ModuleBase::matrix& forcelc, ModulePW::PW_Basis* rho_basis, + const ModuleBase::matrix& vloc, const Charge* const chr) { ModuleBase::TITLE("Forces", "cal_force_loc"); @@ -520,7 +523,7 @@ void Forces::cal_force_loc(ModuleBase::matrix& forcelc, double sinp, cosp; ModuleBase::libm::sincos(phase, &sinp, &cosp); const double factor - = GlobalC::ppcell.vloc(it, rho_basis->ig2igg[ig]) * (cosp * aux[ig].imag() + sinp * aux[ig].real()); + = vloc(it, rho_basis->ig2igg[ig]) * (cosp * aux[ig].imag() + sinp * aux[ig].real()); forcelc(iat, 0) += rho_basis->gcar[ig][0] * factor; forcelc(iat, 1) += rho_basis->gcar[ig][1] * factor; forcelc(iat, 2) += rho_basis->gcar[ig][2] * factor; diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.h b/source/module_hamilt_pw/hamilt_pwdft/forces.h index 712cceac66..eca52c95cf 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces.h +++ b/source/module_hamilt_pw/hamilt_pwdft/forces.h @@ -38,6 +38,8 @@ class Forces ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry* p_symm, Structure_Factor* p_sf, + const pseudopot_cell_vl* locpp, + const pseudopot_cell_vnl* nlpp = nullptr, K_Vectors* pkv = nullptr, ModulePW::PW_Basis_K* psi_basis = nullptr, const psi::Psi, Device>* psi_in = nullptr); @@ -46,9 +48,16 @@ class Forces int nat = 0; int npwx = 0; - void cal_force_loc(ModuleBase::matrix& forcelc, ModulePW::PW_Basis* rho_basis, const Charge* const chr); + void cal_force_loc(ModuleBase::matrix& forcelc, + ModulePW::PW_Basis* rho_basis, + const ModuleBase::matrix& vloc, + const Charge* const chr); void cal_force_ew(ModuleBase::matrix& forceion, ModulePW::PW_Basis* rho_basis, const Structure_Factor* p_sf); - void cal_force_cc(ModuleBase::matrix& forcecc, ModulePW::PW_Basis* rho_basis, const Charge* const chr, UnitCell& ucell_in); + void cal_force_cc(ModuleBase::matrix& forcecc, + ModulePW::PW_Basis* rho_basis, + const Charge* const chr, + const bool* numeric, + UnitCell& ucell_in); /** * @brief This routine computes the atomic force of non-local pseudopotential * F^{NL}_i = \sum_{n,k}f_{nk}\sum_I \sum_{lm,l'm'}D_{l,l'}^{I} [ @@ -65,17 +74,18 @@ class Forces const K_Vectors* p_kv, const ModulePW::PW_Basis_K* psi_basis, const Structure_Factor* p_sf, - pseudopot_cell_vnl* nlpp_in, + const pseudopot_cell_vnl& nlpp_in, const UnitCell& ucell_in, const psi::Psi, Device>* psi_in = nullptr); void cal_force_scc(ModuleBase::matrix& forcescc, ModulePW::PW_Basis* rho_basis, const ModuleBase::matrix& v_current, const bool vnew_exist, + const bool* numeric, const UnitCell& ucell_in); void cal_force_us(ModuleBase::matrix& forcenl, ModulePW::PW_Basis* rho_basis, - pseudopot_cell_vnl* ppcell_in, + const pseudopot_cell_vnl& ppcell_in, const elecstate::ElecState& elec, const UnitCell& ucell); void cal_ylm(int lmax, int npw, const FPTYPE* gk_in, FPTYPE* ylm); diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp index 417c812059..541f8f316c 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp @@ -1,7 +1,6 @@ #include "forces.h" #include "stress_func.h" #include "module_parameter/parameter.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/output_log.h" // new #include "module_base/complexmatrix.h" @@ -33,6 +32,7 @@ template void Forces::cal_force_cc(ModuleBase::matrix& forcecc, ModulePW::PW_Basis* rho_basis, const Charge* const chr, + const bool* numeric, UnitCell& ucell_in) { ModuleBase::TITLE("Forces", "cal_force_cc"); @@ -155,13 +155,13 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, if (ucell_in.atoms[it].ncpp.nlcc) { - // chr->non_linear_core_correction(GlobalC::ppcell.numeric, + // chr->non_linear_core_correction(numeric.numeric, // ucell_in.atoms[it].ncpp.msh, // ucell_in.atoms[it].ncpp.r, // ucell_in.atoms[it].ncpp.rab, // ucell_in.atoms[it].ncpp.rho_atc, // rhocg); - this->deriv_drhoc(GlobalC::ppcell.numeric, + this->deriv_drhoc(numeric, ucell_in.atoms[it].ncpp.msh, ucell_in.atoms[it].ncpp.r.data(), ucell_in.atoms[it].ncpp.rab.data(), diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces_nl.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_nl.cpp index a8fc6976c6..8ecba030f3 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces_nl.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces_nl.cpp @@ -14,12 +14,12 @@ void Forces::cal_force_nl(ModuleBase::matrix& forcenl, const K_Vectors* p_kv, const ModulePW::PW_Basis_K* wfc_basis, const Structure_Factor* p_sf, - pseudopot_cell_vnl* nlpp_in, + const pseudopot_cell_vnl& nlpp, const UnitCell& ucell_in, const psi::Psi, Device>* psi_in) { ModuleBase::TITLE("Forces", "cal_force_nl"); - if (nlpp_in->nkb == 0 || psi_in == nullptr || wfc_basis == nullptr) + if (nlpp.nkb == 0 || psi_in == nullptr || wfc_basis == nullptr) { return; } @@ -30,7 +30,7 @@ void Forces::cal_force_nl(ModuleBase::matrix& forcenl, resmem_var_op()(this->ctx, force, ucell_in.nat * 3); base_device::memory::set_memory_op()(this->ctx, force, 0.0, ucell_in.nat * 3); - hamilt::FS_Nonlocal_tools nl_tools(nlpp_in, &ucell_in, p_kv, wfc_basis, p_sf, wg, &ekb); + hamilt::FS_Nonlocal_tools nl_tools(&nlpp, &ucell_in, p_kv, wfc_basis, p_sf, wg, &ekb); const int nks = wfc_basis->nks; const int max_nbands = wg.nc; diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces_scc.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_scc.cpp index aef02eaa2d..f670ad9b27 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces_scc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces_scc.cpp @@ -1,5 +1,4 @@ #include "forces.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/output_log.h" #include "stress_func.h" // new @@ -28,6 +27,7 @@ void Forces::cal_force_scc(ModuleBase::matrix& forcescc, ModulePW::PW_Basis* rho_basis, const ModuleBase::matrix& vnew, const bool vnew_exist, + const bool* numeric, const UnitCell& ucell_in) { ModuleBase::TITLE("Forces", "cal_force_scc"); ModuleBase::timer::tick("Forces", "cal_force_scc"); @@ -85,7 +85,7 @@ void Forces::cal_force_scc(ModuleBase::matrix& forcescc, for (int nt = 0; nt < ucell_in.ntype; nt++) { // Here we compute the G.ne.0 term const int mesh = ucell_in.atoms[nt].ncpp.msh; - this->deriv_drhoc_scc(GlobalC::ppcell.numeric, + this->deriv_drhoc_scc(numeric, mesh, ucell_in.atoms[nt].ncpp.r.data(), ucell_in.atoms[nt].ncpp.rab.data(), diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces_us.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_us.cpp index aca975431f..4323ee6c4f 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces_us.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces_us.cpp @@ -16,7 +16,7 @@ template void Forces::cal_force_us(ModuleBase::matrix& forcenl, ModulePW::PW_Basis* rho_basis, - pseudopot_cell_vnl* ppcell_in, + const pseudopot_cell_vnl& nlpp, const elecstate::ElecState& elec, const UnitCell& ucell) { @@ -24,7 +24,7 @@ void Forces::cal_force_us(ModuleBase::matrix& forcenl, ModuleBase::timer::tick("Forces", "cal_force_us"); const int npw = rho_basis->npw; - const int nh_tot = ppcell_in->nhm * (ppcell_in->nhm + 1) / 2; + const int nh_tot = nlpp.nhm * (nlpp.nhm + 1) / 2; const std::complex fac = ModuleBase::NEG_IMAG_UNIT * ucell.tpiba; const std::complex ci_tpi = ModuleBase::IMAG_UNIT * ModuleBase::TWO_PI; double* becsum = static_cast, Device>&>(elec).becsum; @@ -39,8 +39,8 @@ void Forces::cal_force_us(ModuleBase::matrix& forcenl, rho_basis->real2recip(&veff.c[is * veff.nc], &vg(is, 0)); } - ModuleBase::matrix ylmk0(ppcell_in->lmaxq * ppcell_in->lmaxq, npw); - ModuleBase::YlmReal::Ylm_Real(ppcell_in->lmaxq * ppcell_in->lmaxq, npw, rho_basis->gcar, ylmk0); + ModuleBase::matrix ylmk0(nlpp.lmaxq * nlpp.lmaxq, npw); + ModuleBase::YlmReal::Ylm_Real(nlpp.lmaxq * nlpp.lmaxq, npw, rho_basis->gcar, ylmk0); double* qnorm = new double[npw]; for (int ig = 0; ig < npw; ig++) @@ -65,7 +65,7 @@ void Forces::cal_force_us(ModuleBase::matrix& forcenl, { for (int jh = ih; jh < atom->ncpp.nh; jh++) { - ppcell_in->radial_fft_q(npw, ih, jh, it, qnorm, ylmk0, &qgm(ijh, 0)); + nlpp.radial_fft_q(npw, ih, jh, it, qnorm, ylmk0, &qgm(ijh, 0)); ijh++; } } diff --git a/source/module_hamilt_pw/hamilt_pwdft/global.cpp b/source/module_hamilt_pw/hamilt_pwdft/global.cpp index 4d569f45ea..cd2568792d 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/global.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/global.cpp @@ -7,7 +7,6 @@ namespace GlobalC #ifdef __EXX Exx_Info exx_info; #endif -pseudopot_cell_vnl ppcell; UnitCell ucell; Parallel_Grid Pgrid; //mohan add 2010-06-06 Parallel_Kpoints Pkpoints; // mohan add 2010-06-07 diff --git a/source/module_hamilt_pw/hamilt_pwdft/global.h b/source/module_hamilt_pw/hamilt_pwdft/global.h index 457bb4a762..bed170b32b 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/global.h +++ b/source/module_hamilt_pw/hamilt_pwdft/global.h @@ -258,7 +258,6 @@ namespace GlobalC #ifdef __EXX extern Exx_Info exx_info; #endif -extern pseudopot_cell_vnl ppcell; } // namespace GlobalC #include "module_cell/parallel_kpoints.h" diff --git a/source/module_hamilt_pw/hamilt_pwdft/hamilt_lcaopw.h b/source/module_hamilt_pw/hamilt_pwdft/hamilt_lcaopw.h index ca0cd823de..eff63d8c0b 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/hamilt_lcaopw.h +++ b/source/module_hamilt_pw/hamilt_pwdft/hamilt_lcaopw.h @@ -13,12 +13,19 @@ namespace hamilt class HamiltLIP : public HamiltPW { public: - HamiltLIP(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv) - : HamiltPW(pot_in, wfc_basis, p_kv) {}; + HamiltLIP(elecstate::Potential* pot_in, + ModulePW::PW_Basis_K* wfc_basis, + K_Vectors* p_kv, + pseudopot_cell_vnl* nlpp) + : HamiltPW(pot_in, wfc_basis, p_kv, nlpp){}; #ifdef __EXX - HamiltLIP(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, Exx_Lip& exx_lip_in) - : HamiltPW(pot_in, wfc_basis, p_kv), exx_lip(exx_lip_in) {}; - Exx_Lip& exx_lip; + HamiltLIP(elecstate::Potential* pot_in, + ModulePW::PW_Basis_K* wfc_basis, + K_Vectors* p_kv, + pseudopot_cell_vnl* nlpp, + Exx_Lip& exx_lip_in) + : HamiltPW(pot_in, wfc_basis, p_kv, nlpp), exx_lip(exx_lip_in){}; + Exx_Lip& exx_lip; #endif }; diff --git a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp index db877bac07..abcc87ffc0 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp @@ -18,11 +18,14 @@ namespace hamilt { -template -HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv) +template +HamiltPW::HamiltPW(elecstate::Potential* pot_in, + ModulePW::PW_Basis_K* wfc_basis, + K_Vectors* pkv, + pseudopot_cell_vnl* nlpp) { this->classname = "HamiltPW"; - this->ppcell = &GlobalC::ppcell; + this->ppcell = nlpp; this->qq_nt = this->ppcell->template get_qq_nt_data(); this->qq_so = this->ppcell->template get_qq_so_data(); this->vkb = this->ppcell->template get_vkb_data(); @@ -100,7 +103,7 @@ HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K if (PARAM.inp.vnl_in_h) { Operator* nonlocal - = new Nonlocal>(isk, &GlobalC::ppcell, &GlobalC::ucell, wfc_basis); + = new Nonlocal>(isk, this->ppcell, &GlobalC::ucell, wfc_basis); if(this->ops == nullptr) { this->ops = nonlocal; diff --git a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h index 07e526161b..7843c0d3f1 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h +++ b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h @@ -20,7 +20,7 @@ class HamiltPW : public Hamilt // otherwise return the real type of T(complex, complex) using Real = typename GetTypeReal::type; public: - HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv); + HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl* nlpp); template explicit HamiltPW(const HamiltPW* hamilt); ~HamiltPW(); diff --git a/source/module_hamilt_pw/hamilt_pwdft/kernels/force_op.h b/source/module_hamilt_pw/hamilt_pwdft/kernels/force_op.h index ce6ae7afb3..539bdebe88 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/kernels/force_op.h +++ b/source/module_hamilt_pw/hamilt_pwdft/kernels/force_op.h @@ -62,8 +62,8 @@ struct cal_force_nl_op /// @param d_wg - input parameter wg /// @param occ - if use the occupation of the bands /// @param d_ekb - input parameter ekb - /// @param qq_nt - GlobalC::ppcell.qq_nt - /// @param deeq - GlobalC::ppcell.deeq + /// @param qq_nt - ppcell.qq_nt + /// @param deeq - ppcell.deeq /// @param becp - intermediate matrix with PARAM.inp.nbands * nkb /// @param dbecp - intermediate matrix with 3 * PARAM.inp.nbands * nkb /// diff --git a/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.h b/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.h index 800bb3eba3..64cc01514a 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.h +++ b/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.h @@ -74,8 +74,8 @@ struct cal_stress_nl_op /// @param d_wg - input parameter wg /// @param occ - if use the occupation of the bands /// @param d_ekb - input parameter ekb - /// @param qq_nt - GlobalC::ppcell.qq_nt - /// @param deeq - GlobalC::ppcell.deeq + /// @param qq_nt - ppcell.qq_nt + /// @param deeq - ppcell.deeq /// @param becp - intermediate matrix with PARAM.inp.nbands * nkb /// @param dbecp - intermediate matrix with 3 * PARAM.inp.nbands * nkb /// diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func.h b/source/module_hamilt_pw/hamilt_pwdft/stress_func.h index ddd1c6ca44..fbd51930e3 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func.h +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func.h @@ -84,6 +84,7 @@ class Stress_Func // 4) the stress from the local pseudopotentials void stress_loc(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, + const ModuleBase::matrix& vloc, const Structure_Factor* p_sf, const bool is_pw, const Charge* const chr); // local pseudopotential part in PW or LCAO @@ -111,6 +112,7 @@ class Stress_Func ModulePW::PW_Basis* rho_basis, const Structure_Factor* p_sf, const bool is_pw, + const bool *numeric, const Charge* const chr); // nonlinear core correction stress in PW or LCAO basis void deriv_drhoc(const bool& numeric, @@ -152,7 +154,7 @@ class Stress_Func ModuleSymmetry::Symmetry* p_symm, ModulePW::PW_Basis_K* wfc_basis, const psi::Psi, Device>* psi_in, - pseudopot_cell_vnl* nlpp_in, + const pseudopot_cell_vnl& nlpp_in, const UnitCell& ucell_in); // nonlocal part in PW basis void get_dvnl1(ModuleBase::ComplexMatrix& vkb, @@ -195,7 +197,7 @@ class Stress_Func * @param dylmk0 [in] derivetives of spherical harmonics * @param dqg [out] the Fourier transform of interest */ - void dqvan2(const pseudopot_cell_vnl* ppcell_in, + void dqvan2(const pseudopot_cell_vnl& ppcell_in, const int ih, const int jh, const int itype, 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 898657febc..91eb74a287 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp @@ -17,6 +17,7 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const Structure_Factor* p_sf, const bool is_pw, + const bool *numeric, const Charge* const chr) { ModuleBase::TITLE("Stress_Func","stress_cc"); @@ -111,7 +112,7 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, { //drhoc(); this->deriv_drhoc( - GlobalC::ppcell.numeric, + numeric, GlobalC::ucell.atoms[nt].ncpp.msh, GlobalC::ucell.atoms[nt].ncpp.r.data(), GlobalC::ucell.atoms[nt].ncpp.rab.data(), @@ -136,7 +137,7 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, sigmadiag += local_sigmadiag.real(); } this->deriv_drhoc ( - GlobalC::ppcell.numeric, + numeric, GlobalC::ucell.atoms[nt].ncpp.msh, GlobalC::ucell.atoms[nt].ncpp.r.data(), GlobalC::ucell.atoms[nt].ncpp.rab.data(), diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_loc.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_loc.cpp index d214c99ba1..8eca06d879 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_loc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_loc.cpp @@ -10,6 +10,7 @@ template void Stress_Func::stress_loc(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, + const ModuleBase::matrix& vloc, const Structure_Factor* p_sf, const bool is_pw, const Charge* const chr) @@ -73,10 +74,10 @@ void Stress_Func::stress_loc(ModuleBase::matrix& sigma, for (int ig=0; ignpw; ig++) { if (rho_basis->ig_gge0 == ig) { - evloc += GlobalC::ppcell.vloc(it, rho_basis->ig2igg[ig]) + evloc += vloc(it, rho_basis->ig2igg[ig]) * (p_sf->strucFac(it, ig) * conj(aux[ig])).real(); } else { - evloc += GlobalC::ppcell.vloc(it, rho_basis->ig2igg[ig]) + evloc += vloc(it, rho_basis->ig2igg[ig]) * (p_sf->strucFac(it, ig) * conj(aux[ig]) * fact).real(); } } diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_nl.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_nl.cpp index 2a761d8118..73b9e08a82 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_nl.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_nl.cpp @@ -18,12 +18,12 @@ void Stress_Func::stress_nl(ModuleBase::matrix& sigma, ModuleSymmetry::Symmetry* p_symm, ModulePW::PW_Basis_K* wfc_basis, const psi::Psi, Device>* psi_in, - pseudopot_cell_vnl* nlpp_in, + const pseudopot_cell_vnl& nlpp_in, const UnitCell& ucell_in) { ModuleBase::TITLE("Stress_Func", "stress_nl"); // skip nkb==0 - if (nlpp_in->nkb == 0 || psi_in == nullptr || wfc_basis == nullptr) + if (nlpp_in.nkb == 0 || psi_in == nullptr || wfc_basis == nullptr) { return; } @@ -34,7 +34,7 @@ void Stress_Func::stress_nl(ModuleBase::matrix& sigma, setmem_var_op()(this->ctx, stress_device, 0, 9); std::vector sigmanlc(9, 0.0); - hamilt::FS_Nonlocal_tools nl_tools(nlpp_in, &ucell_in, p_kv, wfc_basis, p_sf, wg, &ekb); + hamilt::FS_Nonlocal_tools nl_tools(&nlpp_in, &ucell_in, p_kv, wfc_basis, p_sf, wg, &ekb); const int nks = p_kv->get_nks(); const int max_nbands = wg.nc; diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_us.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_us.cpp index 1af5e857c4..464308f317 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_us.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_us.cpp @@ -12,14 +12,14 @@ template void Stress_PW::stress_us(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, - pseudopot_cell_vnl* ppcell_in, + const pseudopot_cell_vnl& nlpp, const UnitCell& ucell) { ModuleBase::TITLE("Stress_Func", "stress_us"); ModuleBase::timer::tick("Stress_Func", "stress_us"); const int npw = rho_basis->npw; - const int nh_tot = ppcell_in->nhm * (ppcell_in->nhm + 1) / 2; + const int nh_tot = nlpp.nhm * (nlpp.nhm + 1) / 2; const std::complex fac = ModuleBase::NEG_IMAG_UNIT * ucell.tpiba; const std::complex ci_tpi = ModuleBase::IMAG_UNIT * ModuleBase::TWO_PI; double* becsum = static_cast, Device>*>(this->pelec)->becsum; @@ -34,8 +34,8 @@ void Stress_PW::stress_us(ModuleBase::matrix& sigma, rho_basis->real2recip(&veff.c[is * veff.nc], &vg(is, 0)); } - ModuleBase::matrix ylmk0(ppcell_in->lmaxq * ppcell_in->lmaxq, npw); - ModuleBase::YlmReal::Ylm_Real(ppcell_in->lmaxq * ppcell_in->lmaxq, npw, rho_basis->gcar, ylmk0); + ModuleBase::matrix ylmk0(nlpp.lmaxq * nlpp.lmaxq, npw); + ModuleBase::YlmReal::Ylm_Real(nlpp.lmaxq * nlpp.lmaxq, npw, rho_basis->gcar, ylmk0); // double* qnorm = new double[npw]; std::vector qnorm_vec(npw); @@ -48,11 +48,11 @@ void Stress_PW::stress_us(ModuleBase::matrix& sigma, // here we compute the integral Q*V for each atom, // I = sum_G G_a exp(-iR.G) Q_nm v^* // (no contribution from G=0) - ModuleBase::matrix dylmk0(ppcell_in->lmaxq * ppcell_in->lmaxq, npw); + ModuleBase::matrix dylmk0(nlpp.lmaxq * nlpp.lmaxq, npw); for (int ipol = 0; ipol < 3; ipol++) { double* gcar_ptr = reinterpret_cast(rho_basis->gcar); - hamilt::Nonlocal_maths::dylmr2(ppcell_in->lmaxq * ppcell_in->lmaxq, + hamilt::Nonlocal_maths::dylmr2(nlpp.lmaxq * nlpp.lmaxq, npw, gcar_ptr, dylmk0.c, @@ -75,7 +75,7 @@ void Stress_PW::stress_us(ModuleBase::matrix& sigma, { for (int jh = ih; jh < atom->ncpp.nh; jh++) { - this->dqvan2(ppcell_in, + this->dqvan2(nlpp, ih, jh, it, @@ -189,7 +189,7 @@ void Stress_PW::stress_us(ModuleBase::matrix& sigma, } template -void Stress_Func::dqvan2(const pseudopot_cell_vnl* ppcell_in, +void Stress_Func::dqvan2(const pseudopot_cell_vnl& nlpp, const int ih, const int jh, const int itype, @@ -207,10 +207,10 @@ void Stress_Func::dqvan2(const pseudopot_cell_vnl* ppcell_in, } // computes the indices which correspond to ih,jh - const int nb = ppcell_in->indv(itype, ih); - const int mb = ppcell_in->indv(itype, jh); - assert(nb < ppcell_in->nbetam); - assert(mb < ppcell_in->nbetam); + const int nb = nlpp.indv(itype, ih); + const int mb = nlpp.indv(itype, jh); + assert(nb < nlpp.nbetam); + assert(mb < nlpp.nbetam); int ijv = 0; if (nb >= mb) { @@ -220,8 +220,8 @@ void Stress_Func::dqvan2(const pseudopot_cell_vnl* ppcell_in, { ijv = mb * (mb + 1) / 2 + nb; } - const int ivl = ppcell_in->nhtolm(itype, ih); - const int jvl = ppcell_in->nhtolm(itype, jh); + const int ivl = nlpp.nhtolm(itype, ih); + const int jvl = nlpp.nhtolm(itype, jh); for (int ig = 0; ig < ng; ig++) { @@ -231,9 +231,9 @@ void Stress_Func::dqvan2(const pseudopot_cell_vnl* ppcell_in, // make the sum over the non zero LM int l = -1; std::complex pref(0.0, 0.0); - for (int lm = 0; lm < ppcell_in->lpx(ivl, jvl); lm++) + for (int lm = 0; lm < nlpp.lpx(ivl, jvl); lm++) { - int lp = ppcell_in->lpl(ivl, jvl, lm); + int lp = nlpp.lpl(ivl, jvl, lm); assert(lp >= 0); assert(lp < 49); if (lp == 0) @@ -264,7 +264,7 @@ void Stress_Func::dqvan2(const pseudopot_cell_vnl* ppcell_in, { l = 6; } - pref = pow(ModuleBase::NEG_IMAG_UNIT, l) * ppcell_in->ap(lp, ivl, jvl); + pref = pow(ModuleBase::NEG_IMAG_UNIT, l) * nlpp.ap(lp, ivl, jvl); double qm1 = -1.0; // any number smaller than qnorm double work = 0.0, work1 = 0.0; @@ -272,14 +272,14 @@ void Stress_Func::dqvan2(const pseudopot_cell_vnl* ppcell_in, { if (std::abs(qnorm[ig] - qm1) > 1e-6) { - work = ModuleBase::PolyInt::Polynomial_Interpolation(ppcell_in->qrad, + work = ModuleBase::PolyInt::Polynomial_Interpolation(nlpp.qrad, itype, l, ijv, PARAM.globalv.nqxq, PARAM.globalv.dq, qnorm[ig]); - work1 = this->Polynomial_Interpolation_nl(ppcell_in->qrad, itype, l, ijv, PARAM.globalv.dq, qnorm[ig]); + work1 = this->Polynomial_Interpolation_nl(nlpp.qrad, itype, l, ijv, PARAM.globalv.dq, qnorm[ig]); qm1 = qnorm[ig]; } dqg[ig] += pref * work * dylmk0(lp, ig) / tpiba; diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp index f32f0067dc..713a233990 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp @@ -8,7 +8,7 @@ template void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, UnitCell& ucell, - pseudopot_cell_vnl* nlpp, + const pseudopot_cell_vnl& nlpp, ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry* p_symm, Structure_Factor* p_sf, @@ -89,10 +89,10 @@ void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, } // local contribution - this->stress_loc(sigmaloc, rho_basis, p_sf, 1, pelec->charge); + this->stress_loc(sigmaloc, rho_basis, nlpp.vloc, p_sf, 1, pelec->charge); // nlcc - this->stress_cc(sigmaxcc, rho_basis, p_sf, 1, pelec->charge); + this->stress_cc(sigmaxcc, rho_basis, p_sf, 1, nlpp.numeric, pelec->charge); // nonlocal this->stress_nl(sigmanl, this->pelec->wg, this->pelec->ekb, p_sf, p_kv, p_symm, wfc_basis, d_psi_in, nlpp, ucell); @@ -100,7 +100,7 @@ void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, // add US term from augmentation charge derivatives if (PARAM.globalv.use_uspp) { - this->stress_us(sigmanl, rho_basis, &GlobalC::ppcell, ucell); + this->stress_us(sigmanl, rho_basis, nlpp, ucell); } // vdw term diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.h b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.h index 4671d25567..70f3ff8d1e 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.h +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.h @@ -13,7 +13,7 @@ class Stress_PW : public Stress_Func // calculate the stress in PW basis void cal_stress(ModuleBase::matrix& smearing_sigmatot, UnitCell& ucell, - pseudopot_cell_vnl* nlpp, + const pseudopot_cell_vnl& nlpp, ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry* p_symm, Structure_Factor* p_sf, @@ -30,7 +30,7 @@ class Stress_PW : public Stress_Func // which is due to the dependence of the Q function on the atomic position void stress_us(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, - pseudopot_cell_vnl* ppcell_in, + const pseudopot_cell_vnl& nlpp, const UnitCell& ucell); // nonlocal part of uspp in PW basis const elecstate::ElecState* pelec = nullptr; diff --git a/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.cpp b/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.cpp index 7f0077871f..6dd2d845f9 100644 --- a/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.cpp @@ -9,10 +9,11 @@ template HamiltSdftPW::HamiltSdftPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, + pseudopot_cell_vnl* nlpp, const int& npol, double* emin_in, double* emax_in) - : HamiltPW(pot_in, wfc_basis, p_kv), ngk(p_kv->ngk) + : HamiltPW(pot_in, wfc_basis, p_kv, nlpp), ngk(p_kv->ngk) { this->classname = "HamiltSdftPW"; this->npwk_max = wfc_basis->npwk_max; diff --git a/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.h b/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.h index ce07ea3222..599060dca4 100644 --- a/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.h +++ b/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.h @@ -23,6 +23,7 @@ class HamiltSdftPW : public HamiltPW HamiltSdftPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, + pseudopot_cell_vnl* nlpp, const int& npol, double* emin_in, double* emax_in); diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp index c652b38635..e3289273a6 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp @@ -41,11 +41,11 @@ void Sto_Forces::cal_stoforce(ModuleBase::matrix& force, ModuleBase::matrix forcecc(this->nat, 3); ModuleBase::matrix forcenl(this->nat, 3); ModuleBase::matrix forcescc(this->nat, 3); - this->cal_force_loc(forcelc, rho_basis, chr); + this->cal_force_loc(forcelc, rho_basis, nlpp.vloc, chr); this->cal_force_ew(forceion, rho_basis, p_sf); this->cal_sto_force_nl(forcenl, wg, pkv, wfc_basis, p_sf, nlpp, ucell, psi, stowf); - this->cal_force_cc(forcecc, rho_basis, chr, ucell); - this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, ucell); + this->cal_force_cc(forcecc, rho_basis, chr, nlpp.numeric, ucell); + this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, nlpp.numeric, ucell); // impose total force = 0 ModuleBase::matrix force_e; diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp index 1ea3b51321..0ca1d4dce2 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp @@ -52,10 +52,10 @@ void Sto_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, this->stress_gga(sigmaxc, rho_basis, chr); // local contribution - this->stress_loc(sigmaloc, rho_basis, p_sf, true, chr); + this->stress_loc(sigmaloc, rho_basis, nlpp->vloc, p_sf, true, chr); // nlcc - this->stress_cc(sigmaxcc, rho_basis, p_sf, true, chr); + this->stress_cc(sigmaxcc, rho_basis, p_sf, true, nlpp->numeric, chr); // nonlocal this->sto_stress_nl(sigmanl, wg, p_sf, p_symm, p_kv, wfc_basis, *nlpp, ucell_in, psi_in, stowf); diff --git a/source/module_hamilt_pw/hamilt_stodft/test/test_hamilt_sto.cpp b/source/module_hamilt_pw/hamilt_stodft/test/test_hamilt_sto.cpp index ed40e17d80..bdbad81516 100644 --- a/source/module_hamilt_pw/hamilt_stodft/test/test_hamilt_sto.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/test/test_hamilt_sto.cpp @@ -11,7 +11,7 @@ void elecstate::Potential::cal_v_eff(Charge const*, UnitCell const*, ModuleBase: void elecstate::Potential::cal_fixed_v(double*){} template -hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv){} +hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl*){} template hamilt::HamiltPW::~HamiltPW(){ delete this->ops; @@ -70,7 +70,7 @@ class TestHamiltSto : public ::testing::Test p_kv = new K_Vectors(); std::vector ngk = {2}; p_kv->ngk = ngk; - hamilt_sto = new hamilt::HamiltSdftPW, base_device::DEVICE_CPU>(pot, wfc_basis, p_kv, npol, &emin, &emax); + hamilt_sto = new hamilt::HamiltSdftPW, base_device::DEVICE_CPU>(pot, wfc_basis, p_kv, nullptr, npol, &emin, &emax); hamilt_sto->ops = new TestOp, base_device::DEVICE_CPU>(); } diff --git a/source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp b/source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp index 817cdb7b52..d947801700 100644 --- a/source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp @@ -8,7 +8,7 @@ ***********************************************/ template -hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv){} +hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv,pseudopot_cell_vnl*){} template hamilt::HamiltPW::~HamiltPW(){}; template @@ -20,9 +20,10 @@ template hamilt::HamiltSdftPW::HamiltSdftPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, + pseudopot_cell_vnl* nlpp, const int& npol, double* emin_in, - double* emax_in): HamiltPW(pot_in, wfc_basis, p_kv), ngk(p_kv->ngk){} + double* emax_in): HamiltPW(pot_in, wfc_basis, p_kv, nlpp), ngk(p_kv->ngk){} template void hamilt::HamiltSdftPW::hPsi_norm(const T* psi_in, T* hpsi, const int& nbands){} diff --git a/source/module_hsolver/test/diago_bpcg_test.cpp b/source/module_hsolver/test/diago_bpcg_test.cpp index 47b51c26dd..39eacf5db6 100644 --- a/source/module_hsolver/test/diago_bpcg_test.cpp +++ b/source/module_hsolver/test/diago_bpcg_test.cpp @@ -97,7 +97,7 @@ class DiagoBPCGPrepare double *en = new double[npw]; int ik = 1; hamilt::Hamilt>* ha; - ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr); + ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr); int* ngk = new int [1]; //psi::Psi> psi(ngk,ik,nband,npw); psi::Psi> psi; diff --git a/source/module_hsolver/test/diago_cg_float_test.cpp b/source/module_hsolver/test/diago_cg_float_test.cpp index fcb1599b1d..bab85fb179 100644 --- a/source/module_hsolver/test/diago_cg_float_test.cpp +++ b/source/module_hsolver/test/diago_cg_float_test.cpp @@ -108,7 +108,7 @@ class DiagoCGPrepare float *en = new float[npw]; int ik = 1; hamilt::Hamilt>* ha; - ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr); + ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr); int* ngk = new int [1]; //psi::Psi> psi(ngk,ik,nband,npw); psi::Psi> psi; diff --git a/source/module_hsolver/test/diago_cg_real_test.cpp b/source/module_hsolver/test/diago_cg_real_test.cpp index 66037cb27f..b6fc2e0401 100644 --- a/source/module_hsolver/test/diago_cg_real_test.cpp +++ b/source/module_hsolver/test/diago_cg_real_test.cpp @@ -108,7 +108,7 @@ class DiagoCGPrepare double* en = new double[npw]; int ik = 1; hamilt::Hamilt* ha; - ha = new hamilt::HamiltPW(nullptr, nullptr, nullptr); + ha = new hamilt::HamiltPW(nullptr, nullptr, nullptr, nullptr); int* ngk = new int[1]; psi::Psi psi; psi.resize(ik, nband, npw); diff --git a/source/module_hsolver/test/diago_cg_test.cpp b/source/module_hsolver/test/diago_cg_test.cpp index 0061350d1e..ada6b90210 100644 --- a/source/module_hsolver/test/diago_cg_test.cpp +++ b/source/module_hsolver/test/diago_cg_test.cpp @@ -104,7 +104,7 @@ class DiagoCGPrepare double *en = new double[npw]; int ik = 1; hamilt::Hamilt>* ha; - ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr); + ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr); int* ngk = new int [1]; //psi::Psi> psi(ngk,ik,nband,npw); psi::Psi> psi; diff --git a/source/module_hsolver/test/diago_david_float_test.cpp b/source/module_hsolver/test/diago_david_float_test.cpp index da52b7272b..98a5d7c775 100644 --- a/source/module_hsolver/test/diago_david_float_test.cpp +++ b/source/module_hsolver/test/diago_david_float_test.cpp @@ -82,7 +82,7 @@ class DiagoDavPrepare //do Diago_David::diag() float* en = new float[npw]; hamilt::Hamilt> *phm; - phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr); + phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr); #ifdef __MPI const hsolver::diag_comm_info comm_info = {MPI_COMM_WORLD, mypnum, nprocs}; diff --git a/source/module_hsolver/test/diago_david_real_test.cpp b/source/module_hsolver/test/diago_david_real_test.cpp index e451d45fbe..67ca3fde8a 100644 --- a/source/module_hsolver/test/diago_david_real_test.cpp +++ b/source/module_hsolver/test/diago_david_real_test.cpp @@ -81,7 +81,7 @@ class DiagoDavPrepare //do Diago_David::diag() double* en = new double[npw]; hamilt::Hamilt* phm; - phm = new hamilt::HamiltPW(nullptr, nullptr, nullptr); + phm = new hamilt::HamiltPW(nullptr, nullptr, nullptr, nullptr); #ifdef __MPI const hsolver::diag_comm_info comm_info = {MPI_COMM_WORLD, mypnum, nprocs}; diff --git a/source/module_hsolver/test/diago_david_test.cpp b/source/module_hsolver/test/diago_david_test.cpp index 4911d40a4f..237233a941 100644 --- a/source/module_hsolver/test/diago_david_test.cpp +++ b/source/module_hsolver/test/diago_david_test.cpp @@ -81,7 +81,7 @@ class DiagoDavPrepare //do Diago_David::diag() double* en = new double[npw]; hamilt::Hamilt> *phm; - phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr); + phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr); #ifdef __MPI const hsolver::diag_comm_info comm_info = {MPI_COMM_WORLD, mypnum, nprocs}; diff --git a/source/module_hsolver/test/diago_mock.h b/source/module_hsolver/test/diago_mock.h index 07d500fdf0..973d01593b 100644 --- a/source/module_hsolver/test/diago_mock.h +++ b/source/module_hsolver/test/diago_mock.h @@ -574,7 +574,7 @@ template<> void hamilt::HamiltPW::updateHk(const int ik) return; } -template<> hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv) +template<> hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*) { this->ops = new OperatorMock; } @@ -589,7 +589,7 @@ template<> void hamilt::HamiltPW>::updateHk(const int ik) return; } -template<> hamilt::HamiltPW>::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv) +template<> hamilt::HamiltPW>::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*) { this->ops = new OperatorMock>; } @@ -604,7 +604,7 @@ template<> void hamilt::HamiltPW>::updateHk(const int ik) return; } -template<> hamilt::HamiltPW>::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv) +template<> hamilt::HamiltPW>::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*) { this->ops = new OperatorMock>; } diff --git a/source/module_hsolver/test/hsolver_pw_sup.h b/source/module_hsolver/test/hsolver_pw_sup.h index bc32b777e4..89ee061c64 100644 --- a/source/module_hsolver/test/hsolver_pw_sup.h +++ b/source/module_hsolver/test/hsolver_pw_sup.h @@ -190,6 +190,8 @@ void diago_PAO_in_pw_k2( psi::Psi, base_device::DEVICE_CPU>& wvf, ModulePW::PW_Basis_K* wfc_basis, wavefunc* p_wf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, hamilt::Hamilt, base_device::DEVICE_CPU>* phm_in) { for (int i = 0; i < wvf.size(); i++) { wvf.get_pointer()[i] = std::complex((float)i + 1, 0); @@ -203,6 +205,8 @@ void diago_PAO_in_pw_k2( psi::Psi, base_device::DEVICE_CPU>& wvf, ModulePW::PW_Basis_K* wfc_basis, wavefunc* p_wf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, hamilt::Hamilt, base_device::DEVICE_CPU>* phm_in) { for (int i = 0; i < wvf.size(); i++) { wvf.get_pointer()[i] = std::complex((double)i + 1, 0); diff --git a/source/module_hsolver/test/hsolver_supplementary_mock.h b/source/module_hsolver/test/hsolver_supplementary_mock.h index 6f7b0df338..2ee381ad48 100644 --- a/source/module_hsolver/test/hsolver_supplementary_mock.h +++ b/source/module_hsolver/test/hsolver_supplementary_mock.h @@ -44,6 +44,7 @@ void ElecState::print_eigenvalue(std::ofstream& ofs) void ElecState::init_scf(const int istep, const ModuleBase::ComplexMatrix& strucfac, + const bool*, ModuleSymmetry::Symmetry&, const void*) { diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 9c853f9040..ac93d06bd9 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -441,7 +441,6 @@ void Input_Conv::Convert() } #endif // __LCAO #endif // __EXX - GlobalC::ppcell.cell_factor = PARAM.inp.cell_factor; // LiuXh add 20180619 //---------------------------------------------------------- // reset symmetry flag to avoid error diff --git a/source/module_psi/psi_init.cpp b/source/module_psi/psi_init.cpp index 419a396cdf..fdd08958ba 100644 --- a/source/module_psi/psi_init.cpp +++ b/source/module_psi/psi_init.cpp @@ -92,7 +92,8 @@ void PSIInit::allocate_psi(Psi>*& psi, const int nks, const int* ngk, const int npwx, - Structure_Factor* p_sf) + Structure_Factor* p_sf, + pseudopot_cell_vnl* p_ppcell) { // allocate memory for std::complex datatype psi // New psi initializer in ABACUS, Developer's note: @@ -125,7 +126,7 @@ void PSIInit::allocate_psi(Psi>*& psi, // however, init_at_1 does not actually initialize the psi, instead, it is a // function to calculate a interpolate table saving overlap intergral or say // Spherical Bessel Transform of atomic orbitals. - this->wf_old.init_at_1(p_sf); + this->wf_old.init_at_1(p_sf, &p_ppcell->tab_at); // similarly, wfcinit not really initialize any wavefunction, instead, it initialize // the mapping from ixy, the 1d flattened index of point on fft grid (x, y) plane, // to the index of "stick", composed of grid points. @@ -134,7 +135,7 @@ void PSIInit::allocate_psi(Psi>*& psi, } template -void PSIInit::make_table(const int nks, Structure_Factor* p_sf) +void PSIInit::make_table(const int nks, Structure_Factor* p_sf, pseudopot_cell_vnl* p_ppcell) { if (this->use_psiinitializer) { @@ -142,7 +143,7 @@ void PSIInit::make_table(const int nks, Structure_Factor* p_sf) else // old initialization method, used in EXX calculation { this->wf_old.init_after_vc(nks); // reallocate wanf2, the planewave expansion of lcao - this->wf_old.init_at_1(p_sf); // re-calculate tab_at, the overlap matrix between atomic pswfc and jlq + this->wf_old.init_at_1(p_sf, &p_ppcell->tab_at); // re-calculate tab_at, the overlap matrix between atomic pswfc and jlq } } @@ -150,6 +151,7 @@ template void PSIInit::initialize_psi(Psi>* psi, psi::Psi* kspw_psi, hamilt::Hamilt* p_hamilt, + const pseudopot_cell_vnl& nlpp, std::ofstream& ofs_running, const bool is_already_initpsi) { @@ -269,7 +271,14 @@ void PSIInit::initialize_psi(Psi>* psi, kspw_psi->fix_k(ik); /// for psi init guess!!!! - hamilt::diago_PAO_in_pw_k2(this->ctx, ik, *(kspw_psi), this->pw_wfc, &this->wf_old, p_hamilt); + hamilt::diago_PAO_in_pw_k2(this->ctx, + ik, + *(kspw_psi), + this->pw_wfc, + &this->wf_old, + nlpp.tab_at, + nlpp.lmaxkb, + p_hamilt); } } else @@ -278,7 +287,14 @@ void PSIInit::initialize_psi(Psi>* psi, kspw_psi->fix_k(ik); /// for psi init guess!!!! - hamilt::diago_PAO_in_pw_k2(this->ctx, ik, *(kspw_psi), this->pw_wfc, &this->wf_old, p_hamilt); + hamilt::diago_PAO_in_pw_k2(this->ctx, + ik, + *(kspw_psi), + this->pw_wfc, + &this->wf_old, + nlpp.tab_at, + nlpp.lmaxkb, + p_hamilt); } } } diff --git a/source/module_psi/psi_init.h b/source/module_psi/psi_init.h index a631939dfe..df0ec83446 100644 --- a/source/module_psi/psi_init.h +++ b/source/module_psi/psi_init.h @@ -35,10 +35,11 @@ class PSIInit const int nks, //< number of k-points in the current pool const int* ngk, //< number of G-vectors in the current pool const int npwx, //< max number of plane waves of all pools - Structure_Factor* p_sf); //< structure factor + Structure_Factor* p_sf, //< structure factor + pseudopot_cell_vnl* p_ppcell); //< nonlocal pseudopotential // make interpolate table - void make_table(const int nks, Structure_Factor* p_sf); + void make_table(const int nks, Structure_Factor* p_sf, pseudopot_cell_vnl* p_ppcell); //------------------------ only for psi_initializer -------------------- /** @@ -52,6 +53,7 @@ class PSIInit void initialize_psi(Psi>* psi, psi::Psi* kspw_psi, hamilt::Hamilt* p_hamilt, + const pseudopot_cell_vnl& nlpp, std::ofstream& ofs_running, const bool is_already_initpsi); diff --git a/source/module_psi/wavefunc.cpp b/source/module_psi/wavefunc.cpp index 87686f8edc..c51bcb3048 100644 --- a/source/module_psi/wavefunc.cpp +++ b/source/module_psi/wavefunc.cpp @@ -188,6 +188,8 @@ void diago_PAO_in_pw_k2(const int& ik, psi::Psi>& wvf, ModulePW::PW_Basis_K* wfc_basis, wavefunc* p_wf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, hamilt::Hamilt>* phm_in) { ModuleBase::TITLE("wavefunc", "diago_PAO_in_pw_k2"); @@ -290,9 +292,10 @@ void diago_PAO_in_pw_k2(const int& ik, p_wf->atomic_wfc(ik, current_nbasis, GlobalC::ucell.lmax_ppwf, + lmaxkb, wfc_basis, wfcatom, - GlobalC::ppcell.tab_at, + tab_at, PARAM.globalv.nqx, PARAM.globalv.dq); @@ -347,6 +350,8 @@ void diago_PAO_in_pw_k2(const int& ik, psi::Psi>& wvf, ModulePW::PW_Basis_K* wfc_basis, wavefunc* p_wf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, hamilt::Hamilt>* phm_in) { ModuleBase::TITLE("wavefunc", "diago_PAO_in_pw_k2"); @@ -447,9 +452,10 @@ void diago_PAO_in_pw_k2(const int& ik, p_wf->atomic_wfc(ik, current_nbasis, GlobalC::ucell.lmax_ppwf, + lmaxkb, wfc_basis, wfcatom, - GlobalC::ppcell.tab_at, + tab_at, PARAM.globalv.nqx, PARAM.globalv.dq); @@ -501,9 +507,11 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_CPU* ctx, psi::Psi, base_device::DEVICE_CPU>& wvf, ModulePW::PW_Basis_K* wfc_basis, wavefunc* p_wf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, hamilt::Hamilt, base_device::DEVICE_CPU>* phm_in) { - diago_PAO_in_pw_k2(ik, wvf, wfc_basis, p_wf, phm_in); + diago_PAO_in_pw_k2(ik, wvf, wfc_basis, p_wf, tab_at, lmaxkb, phm_in); } template <> @@ -512,9 +520,11 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_CPU* ctx, psi::Psi, base_device::DEVICE_CPU>& wvf, ModulePW::PW_Basis_K* wfc_basis, wavefunc* p_wf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, hamilt::Hamilt, base_device::DEVICE_CPU>* phm_in) { - diago_PAO_in_pw_k2(ik, wvf, wfc_basis, p_wf, phm_in); + diago_PAO_in_pw_k2(ik, wvf, wfc_basis, p_wf, tab_at, lmaxkb, phm_in); } #if ((defined __CUDA) || (defined __ROCM)) @@ -524,6 +534,8 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_GPU* ctx, psi::Psi, base_device::DEVICE_GPU>& wvf, ModulePW::PW_Basis_K* wfc_basis, wavefunc* p_wf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, hamilt::Hamilt, base_device::DEVICE_GPU>* phm_in) { ModuleBase::TITLE("wavefunc", "diago_PAO_in_pw_k2"); @@ -555,9 +567,10 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_GPU* ctx, p_wf->atomic_wfc(ik, current_nbasis, GlobalC::ucell.lmax_ppwf, + lmaxkb, wfc_basis, wfcatom, - GlobalC::ppcell.tab_at, + tab_at, PARAM.globalv.nqx, PARAM.globalv.dq); if (PARAM.inp.init_wfc == "atomic+random" && starting_nw == GlobalC::ucell.natomwfc) // added by qianrui 2021-5-16 @@ -628,6 +641,8 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_GPU* ctx, psi::Psi, base_device::DEVICE_GPU>& wvf, ModulePW::PW_Basis_K* wfc_basis, wavefunc* p_wf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, hamilt::Hamilt, base_device::DEVICE_GPU>* phm_in) { ModuleBase::TITLE("wavefunc", "diago_PAO_in_pw_k2"); @@ -658,9 +673,10 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_GPU* ctx, p_wf->atomic_wfc(ik, current_nbasis, GlobalC::ucell.lmax_ppwf, + lmaxkb, wfc_basis, wfcatom, - GlobalC::ppcell.tab_at, + tab_at, PARAM.globalv.nqx, PARAM.globalv.dq); if (PARAM.inp.init_wfc == "atomic+random" && starting_nw == GlobalC::ucell.natomwfc) // added by qianrui 2021-5-16 diff --git a/source/module_psi/wavefunc.h b/source/module_psi/wavefunc.h index 44ca10a849..adc9146bc6 100644 --- a/source/module_psi/wavefunc.h +++ b/source/module_psi/wavefunc.h @@ -33,13 +33,21 @@ void diago_PAO_in_pw_k2(const int& ik, psi::Psi>& wvf, ModulePW::PW_Basis_K* wfc_basis, wavefunc* p_wf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, hamilt::Hamilt>* phm_in = nullptr); void diago_PAO_in_pw_k2(const int& ik, psi::Psi>& wvf, ModulePW::PW_Basis_K* wfc_basis, wavefunc* p_wf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, hamilt::Hamilt>* phm_in = nullptr); -void diago_PAO_in_pw_k2(const int& ik, ModuleBase::ComplexMatrix& wvf, wavefunc* p_wf); +void diago_PAO_in_pw_k2(const int& ik, + ModuleBase::ComplexMatrix& wvf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, + wavefunc* p_wf); template void diago_PAO_in_pw_k2(const Device* ctx, @@ -47,6 +55,8 @@ void diago_PAO_in_pw_k2(const Device* ctx, psi::Psi, Device>& wvf, ModulePW::PW_Basis_K* wfc_basis, wavefunc* p_wf, + const ModuleBase::realArray& tab_at, + const int& lmaxkb, hamilt::Hamilt, Device>* phm_in = nullptr); } // namespace hamilt diff --git a/source/module_psi/wf_atomic.cpp b/source/module_psi/wf_atomic.cpp index 484fc3b847..c5d2acf21d 100644 --- a/source/module_psi/wf_atomic.cpp +++ b/source/module_psi/wf_atomic.cpp @@ -31,12 +31,7 @@ WF_atomic::~WF_atomic() } } -//========================================================== -// MEMBER FUNCTION : -// NAME : init_at_1(init a table with the radial Fourier -// transform of the atomic WF_atomictions) -//========================================================== -void WF_atomic::init_at_1(Structure_Factor *sf_in) +void WF_atomic::init_at_1(Structure_Factor *sf_in, ModuleBase::realArray* tab_at) { if(PARAM.inp.use_paw) { return; } @@ -60,7 +55,7 @@ void WF_atomic::init_at_1(Structure_Factor *sf_in) // needed to normalize atomic wfcs (not a bad idea in general and // necessary to compute correctly lda+U projections) - GlobalC::ppcell.tab_at.zero_out(); + tab_at->zero_out(); //---------------------------------------------------------- // EXPLAIN : If use gauss orbitals to represent aotmic // orbitals (controlled by parameters) @@ -199,7 +194,7 @@ void WF_atomic::init_at_1(Structure_Factor *sf_in) double vqint = 0.0; ModuleBase::Integral::Simpson_Integral(atom->ncpp.msh, vchi, atom->ncpp.rab.data(), vqint); - GlobalC::ppcell.tab_at(it, ic, iq) = vqint * pref; + tab_at->operator()(it, ic, iq) = vqint * pref; // if( it == 0 && ic == 0 ) // { // @@ -252,9 +247,10 @@ void WF_atomic::print_PAOs()const return; } -void WF_atomic::atomic_wfc(const int ik, - const int np, - const int lmax_wfc, +void WF_atomic::atomic_wfc(const int& ik, + const int& np, + const int& lmax_wfc, + const int& lmaxkb, const ModulePW::PW_Basis_K* wfc_basis, ModuleBase::ComplexMatrix& wfcatom, const ModuleBase::realArray& table_q, @@ -335,7 +331,7 @@ void WF_atomic::atomic_wfc(const int ik, { if(fabs(fact[is])>1e-8) { - const int ind = GlobalC::ppcell.lmaxkb + soc.sph_ind(l,j,m,is); + const int ind = lmaxkb + soc.sph_ind(l,j,m,is); ModuleBase::GlobalFunc::ZEROS(aux, np); for(int n1=0;n1<2*l+1;n1++){ const int lm = l*l +n1; diff --git a/source/module_psi/wf_atomic.h b/source/module_psi/wf_atomic.h index babd02ffa8..de2bc27efb 100644 --- a/source/module_psi/wf_atomic.h +++ b/source/module_psi/wf_atomic.h @@ -32,16 +32,22 @@ class WF_atomic ModuleBase::ComplexMatrix *wanf2 = nullptr; // wannier functions in the PW basis - void init_at_1(Structure_Factor *sf_in); // from init_at_1.f90 + /** + * @brief init a table with the radial Fourier transform of the atomic WF_atomictions + * @param sf_in [out] the structure factor + * @param tab_at [out] atomic table + */ + void init_at_1(Structure_Factor *sf_in, ModuleBase::realArray* tab_at); void print_PAOs()const; public: //template change to public, will be refactor later. added by zhengdy 20230302 int *irindex = nullptr; - void atomic_wfc(const int ik, - const int np, - const int lmax_wfc, + void atomic_wfc(const int& ik, + const int& np, + const int& lmax_wfc, + const int& lmaxkb, const ModulePW::PW_Basis_K* wfc_basis, ModuleBase::ComplexMatrix& wfcatom, const ModuleBase::realArray& table_q,