diff --git a/source/Makefile.Objects b/source/Makefile.Objects index e03f39948e4..45903f3a733 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -496,6 +496,10 @@ OBJS_LCAO=DM_gamma.o\ FORCE_gamma_tvnl.o\ FORCE_gamma_vl.o\ FORCE_k.o\ + foverlap_k.o\ + ftvnl_dphi_k.o\ + fvl_dphi_k.o\ + fvnl_dbeta_k.o\ LCAO_gen_fixedH.o\ LCAO_hamilt.o\ LCAO_matrix.o\ diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp index 6dd00a7f4f7..a70153fe687 100644 --- a/source/module_elecstate/elecstate_lcao.cpp +++ b/source/module_elecstate/elecstate_lcao.cpp @@ -142,9 +142,9 @@ void ElecStateLCAO>::psiToRho(const psi::Psiuhm->GK.transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint + this->gint_k->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint Gint_inout inout(this->loc->DM_R, this->charge->rho, Gint_Tools::job_type::rho); - this->uhm->GK.cal_gint(&inout); + this->gint_k->cal_gint(&inout); if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { @@ -153,7 +153,7 @@ void ElecStateLCAO>::psiToRho(const psi::Psicharge->kin_r[is], this->charge->nrxx); } Gint_inout inout1(this->loc->DM_R, this->charge->kin_r, Gint_Tools::job_type::tau); - this->uhm->GK.cal_gint(&inout1); + this->gint_k->cal_gint(&inout1); } this->charge->renormalize_rho(); @@ -216,9 +216,13 @@ void ElecStateLCAO::psiToRho(const psi::Psi& psi) // calculate the charge density on real space grid. //------------------------------------------------------------ ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!"); - this->uhm->GG.transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint + + this->gint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint + Gint_inout inout(this->loc->DM, this->charge->rho, Gint_Tools::job_type::rho); - this->uhm->GG.cal_gint(&inout); + + this->gint_gamma->cal_gint(&inout); + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { for (int is = 0; is < GlobalV::NSPIN; is++) @@ -226,7 +230,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi& psi) ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx); } Gint_inout inout1(this->loc->DM, this->charge->kin_r, Gint_Tools::job_type::tau); - this->uhm->GG.cal_gint(&inout1); + this->gint_gamma->cal_gint(&inout1); } this->charge->renormalize_rho(); diff --git a/source/module_elecstate/elecstate_lcao.h b/source/module_elecstate/elecstate_lcao.h index c47e297b027..977e0a2fdb3 100644 --- a/source/module_elecstate/elecstate_lcao.h +++ b/source/module_elecstate/elecstate_lcao.h @@ -2,7 +2,8 @@ #define ELECSTATELCAO_H #include "elecstate.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h" +#include "module_hamilt_lcao/module_gint/gint_gamma.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h" #include "module_elecstate/module_dm/density_matrix.h" @@ -18,6 +19,8 @@ class ElecStateLCAO : public ElecState const K_Vectors* klist_in , int nks_in, Local_Orbital_Charge* loc_in , + Gint_Gamma* gint_gamma_in, //mohan add 2024-04-01 + Gint_k* gint_k_in, //mohan add 2024-04-01 LCAO_Hamilt* uhm_in , Local_Orbital_wfc* lowf_in , ModulePW::PW_Basis* rhopw_in , @@ -25,15 +28,18 @@ class ElecStateLCAO : public ElecState { init_ks(chg_in, klist_in, nks_in, rhopw_in, bigpw_in); this->loc = loc_in; - this->uhm = uhm_in; + this->gint_gamma = gint_gamma_in; // mohan add 2024-04-01 + this->gint_k = gint_k_in; // mohan add 2024-04-01 this->lowf = lowf_in; this->classname = "ElecStateLCAO"; } virtual ~ElecStateLCAO() { - if (this->DM != nullptr) - delete this->DM; + if (this->DM != nullptr) + { + delete this->DM; + } } // void init(Charge* chg_in):charge(chg_in){} override; @@ -69,7 +75,9 @@ class ElecStateLCAO : public ElecState // void rhoBandK(const psi::Psi>& psi); Local_Orbital_Charge* loc = nullptr; - LCAO_Hamilt* uhm = nullptr; + Gint_Gamma* gint_gamma = nullptr; // mohan add 2024-04-01 + Gint_k* gint_k = nullptr; // mohan add 2024-04-01 + //LCAO_Hamilt* uhm = nullptr; // mohan modify 2024-04-01 Local_Orbital_wfc* lowf = nullptr; DensityMatrix* DM = nullptr; @@ -86,4 +94,4 @@ bool ElecStateLCAO::need_psi_grid = 1; } // namespace elecstate -#endif \ No newline at end of file +#endif diff --git a/source/module_elecstate/elecstate_lcao_tddft.cpp b/source/module_elecstate/elecstate_lcao_tddft.cpp index 3d3539fd1a9..cc0e7133e4e 100644 --- a/source/module_elecstate/elecstate_lcao_tddft.cpp +++ b/source/module_elecstate/elecstate_lcao_tddft.cpp @@ -63,9 +63,9 @@ void ElecStateLCAO_TDDFT::psiToRho_td(const psi::Psi>& psi) //------------------------------------------------------------ ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!"); - this->uhm->GK.transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint + this->gint_k->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint Gint_inout inout(this->loc->DM_R, this->charge->rho, Gint_Tools::job_type::rho); // rho calculation - this->uhm->GK.cal_gint(&inout); + this->gint_k->cal_gint(&inout); this->charge->renormalize_rho(); @@ -110,5 +110,4 @@ void ElecStateLCAO_TDDFT::calculate_weights_td() } - -} // namespace elecstate \ No newline at end of file +} // namespace elecstate diff --git a/source/module_elecstate/elecstate_lcao_tddft.h b/source/module_elecstate/elecstate_lcao_tddft.h index f84a2c5ec21..140759e9d47 100644 --- a/source/module_elecstate/elecstate_lcao_tddft.h +++ b/source/module_elecstate/elecstate_lcao_tddft.h @@ -3,7 +3,7 @@ #include "elecstate.h" #include "elecstate_lcao.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h" @@ -16,14 +16,14 @@ class ElecStateLCAO_TDDFT : public ElecStateLCAO> const K_Vectors* klist_in , int nks_in, Local_Orbital_Charge* loc_in , - LCAO_Hamilt* uhm_in , + Gint_k* gint_k_in, //mohan add 2024-04-01 Local_Orbital_wfc* lowf_in , ModulePW::PW_Basis* rhopw_in , ModulePW::PW_Basis_Big* bigpw_in ) { init_ks(chg_in, klist_in, nks_in, rhopw_in, bigpw_in); this->loc = loc_in; - this->uhm = uhm_in; + this->gint_k = gint_k_in; this->lowf = lowf_in; this->classname = "ElecStateLCAO_TDDFT"; } @@ -33,4 +33,4 @@ class ElecStateLCAO_TDDFT : public ElecStateLCAO> } // namespace elecstate -#endif \ No newline at end of file +#endif diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 3379777f905..d5c7a9ae0e1 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -111,10 +111,13 @@ void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) // autoset nbands in ElecState, it should before basis_init (for Psi 2d divid) if (this->pelec == nullptr) { - this->pelec = new elecstate::ElecStateLCAO(&(this->chr), + this->pelec = new elecstate::ElecStateLCAO( + &(this->chr), &(this->kv), this->kv.nks, &(this->LOC), + &(this->GG), // mohan add 2024-04-01 + &(this->GK), // mohan add 2024-04-01 &(this->UHM), &(this->LOWF), this->pw_rho, @@ -255,6 +258,8 @@ void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) &(this->kv), this->kv.nks, &(this->LOC), + &(this->GG), // mohan add 2024-04-01 + &(this->GK), // mohan add 2024-04-01 &(this->UHM), &(this->LOWF), this->pw_rho, @@ -313,6 +318,8 @@ void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) this->pelec, this->psi, this->UHM, + this->GG, // mohan add 2024-04-01 + this->GK, // mohan add 2024-04-01 force, this->scs, this->sf, @@ -665,7 +672,7 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) // update Gint_K if (!GlobalV::GAMMA_ONLY_LOCAL) { - this->UHM.GK.renew(); + this->GK.renew(); } // update real space Hamiltonian this->p_hamilt->refresh(); @@ -803,7 +810,7 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) { if (!GlobalV::GAMMA_ONLY_LOCAL && hsolver::HSolverLCAO::out_mat_hs[0]) { - this->UHM.GK.renew(true); + this->GK.renew(true); } for (int ik = 0; ik < this->kv.nks; ++ik) { @@ -1016,9 +1023,24 @@ void ESolver_KS_LCAO::after_scf(const int istep) bool out_exc = true; // tmp, add parameter! if (GlobalV::out_mat_xc) { - ModuleIO::write_Vxc(GlobalV::NSPIN, GlobalV::NLOCAL, GlobalV::DRANK, - *this->psi, GlobalC::ucell, this->sf, *this->pw_rho, *this->pw_rhod, GlobalC::ppcell.vloc, - *this->pelec->charge, this->UHM, this->LM, this->LOC, this->kv, this->pelec->wg, GlobalC::GridD); + ModuleIO::write_Vxc( + GlobalV::NSPIN, + GlobalV::NLOCAL, + GlobalV::DRANK, + *this->psi, + GlobalC::ucell, + this->sf, + *this->pw_rho, + *this->pw_rhod, + GlobalC::ppcell.vloc, + *this->pelec->charge, + this->GG, + this->GK, + this->LM, + this->LOC, + this->kv, + this->pelec->wg, + GlobalC::GridD); } #ifdef __EXX @@ -1181,7 +1203,8 @@ ModuleIO::Output_DM1 ESolver_KS_LCAO::create_Output_DM1(int istep) template ModuleIO::Output_Mat_Sparse ESolver_KS_LCAO::create_Output_Mat_Sparse(int istep) { - return ModuleIO::Output_Mat_Sparse(hsolver::HSolverLCAO::out_mat_hsR, + return ModuleIO::Output_Mat_Sparse( + hsolver::HSolverLCAO::out_mat_hsR, hsolver::HSolverLCAO::out_mat_dh, hsolver::HSolverLCAO::out_mat_t, INPUT.out_mat_r, @@ -1189,6 +1212,7 @@ ModuleIO::Output_Mat_Sparse ESolver_KS_LCAO::create_Output_Mat_Spars this->pelec->pot->get_effective_v(), *this->LOWF.ParaV, this->UHM, + this->GK, // mohan add 2024-04-01 this->LM, this->kv, this->p_hamilt); diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 49d95cdfe55..c8c81ab203d 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -6,6 +6,10 @@ #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h" +// for grid integration +#include "module_hamilt_lcao/module_gint/gint_gamma.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" + #include "module_basis/module_ao/ORB_control.h" #ifdef __EXX #include "module_ri/Mix_DMk_2D.h" @@ -75,6 +79,12 @@ namespace ModuleESolver // we will get rid of this class soon, don't use it, mohan 2024-03-28 LCAO_Hamilt UHM; + // used for k-dependent grid integration. + Gint_k GK; + + // used for gamma only algorithms. + Gint_Gamma GG; + // we will get rid of this class soon, don't use it, mohan 2024-03-28 LCAO_Matrix LM; diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 57a0bcd3eac..290feed8b30 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -126,7 +126,12 @@ void ESolver_KS_LCAO::beforesolver(const int istep) } // prepare grid in Gint - this->UHM.grid_prepare(this->GridT, *this->pw_rho, *this->pw_big); + this->UHM.grid_prepare( + this->GridT, + this->GG, + this->GK, + *this->pw_rho, + *this->pw_big); // init Hamiltonian if (this->p_hamilt != nullptr) @@ -137,8 +142,8 @@ void ESolver_KS_LCAO::beforesolver(const int istep) if (this->p_hamilt == nullptr) { elecstate::DensityMatrix* DM = dynamic_cast*>(this->pelec)->get_DM(); - this->p_hamilt = new hamilt::HamiltLCAO(GlobalV::GAMMA_ONLY_LOCAL ? &(this->UHM.GG) : nullptr, - GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->UHM.GK), + this->p_hamilt = new hamilt::HamiltLCAO(GlobalV::GAMMA_ONLY_LOCAL ? &(this->GG) : nullptr, + GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->GK), &(this->UHM.genH), &(this->LM), &(this->LOC), @@ -194,7 +199,7 @@ void ESolver_KS_LCAO::beforesolver(const int istep) if (GlobalV::GAMMA_ONLY_LOCAL) { Gint_inout inout(this->LOC.DM, this->pelec->charge->rho, Gint_Tools::job_type::rho); - this->UHM.GG.cal_gint(&inout); + this->GG.cal_gint(&inout); if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { for (int is = 0; is < GlobalV::NSPIN; is++) @@ -202,13 +207,13 @@ void ESolver_KS_LCAO::beforesolver(const int istep) ModuleBase::GlobalFunc::ZEROS(this->pelec->charge->kin_r[0], this->pw_rho->nrxx); } Gint_inout inout1(this->LOC.DM, this->pelec->charge->kin_r, Gint_Tools::job_type::tau); - this->UHM.GG.cal_gint(&inout1); + this->GG.cal_gint(&inout1); } } else { Gint_inout inout(this->LOC.DM_R, this->pelec->charge->rho, Gint_Tools::job_type::rho); - this->UHM.GK.cal_gint(&inout); + this->GK.cal_gint(&inout); if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { for (int is = 0; is < GlobalV::NSPIN; is++) @@ -216,7 +221,7 @@ void ESolver_KS_LCAO::beforesolver(const int istep) ModuleBase::GlobalFunc::ZEROS(this->pelec->charge->kin_r[0], this->pw_rho->nrxx); } Gint_inout inout1(this->LOC.DM_R, this->pelec->charge->kin_r, Gint_Tools::job_type::tau); - this->UHM.GK.cal_gint(&inout1); + this->GK.cal_gint(&inout1); } } @@ -389,7 +394,7 @@ void ESolver_KS_LCAO::others(const int istep) else if (GlobalV::CALCULATION == "get_pchg") { IState_Charge ISC(this->psi, this->LOC); - ISC.begin(this->UHM.GG, + ISC.begin(this->GG, this->pelec, this->pw_rho, this->pw_big, @@ -411,7 +416,7 @@ void ESolver_KS_LCAO::others(const int istep) this->pw_wfc, this->pw_big, this->LOWF, - this->UHM.GG, + this->GG, INPUT.out_wfc_pw, this->wf.out_wfc_r, this->kv, @@ -427,7 +432,7 @@ void ESolver_KS_LCAO::others(const int istep) this->pw_wfc, this->pw_big, this->LOWF, - this->UHM.GK, + this->GK, INPUT.out_wfc_pw, this->wf.out_wfc_r, this->kv, diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index a9bd2e44313..182242eb6a3 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -79,8 +79,8 @@ void ESolver_KS_LCAO_TDDFT::init(Input& inp, UnitCell& ucell) this->pelec = new elecstate::ElecStateLCAO_TDDFT(&(this->chr), &(kv), kv.nks, - &(this->LOC), - &(this->UHM), + &(this->LOC), + &(this->GK), // mohan add 2024-04-01 &(this->LOWF), this->pw_rho, pw_big); @@ -264,7 +264,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) { if (!GlobalV::GAMMA_ONLY_LOCAL) { - this->UHM.GK.renew(true); + this->GK.renew(true); } for (int ik = 0; ik < kv.nks; ++ik) { diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt index 10b1df34bde..3a6bcdd3fd6 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt +++ b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt @@ -21,6 +21,10 @@ if(ENABLE_LCAO) FORCE_gamma_tvnl.cpp FORCE_gamma_vl.cpp FORCE_k.cpp + foverlap_k.cpp + ftvnl_dphi_k.cpp + fvl_dphi_k.cpp + fvnl_dbeta_k.cpp LCAO_gen_fixedH.cpp LCAO_hamilt.cpp LCAO_matrix.cpp @@ -48,4 +52,4 @@ if(ENABLE_LCAO) add_coverage(hamilt_lcao) endif() -endif() \ No newline at end of file +endif() diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index 4fa805e78de..e6944bd2cc7 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -34,8 +34,10 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, LCAO_Matrix &lm, const elecstate::ElecState* pelec, const psi::Psi* psi, - LCAO_Hamilt& uhm, - ModuleBase::matrix& fcs, + LCAO_Hamilt& uhm, + Gint_Gamma &gint_gamma, // mohan add 2024-04-01 + Gint_k &gint_k, // mohan add 2024-04-01 + ModuleBase::matrix& fcs, ModuleBase::matrix& scs, const Structure_Factor& sf, const K_Vectors& kv, @@ -138,7 +140,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, //-------------------------------------------------------- // implement four terms which needs integration //-------------------------------------------------------- - this->calForceStressIntegralPart(GlobalV::GAMMA_ONLY_LOCAL, + this->integral_part(GlobalV::GAMMA_ONLY_LOCAL, isforce, isstress, loc, @@ -158,9 +160,12 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, svl_dphi, #endif uhm, + gint_gamma, + gint_k, pv, lm, kv); + // implement vdw force or stress here // Peize Lin add 2014-04-04, update 2021-03-09 // jiyy add 2019-05-18, update 2021-05-02 @@ -718,7 +723,7 @@ void Force_Stress_LCAO::calForcePwPart(ModuleBase::matrix& fvl_dvl, // overlap, kinetic, nonlocal pseudopotential, Local potential terms in force and stress template<> -void Force_Stress_LCAO::calForceStressIntegralPart( +void Force_Stress_LCAO::integral_part( const bool isGammaOnly, const bool isforce, const bool isstress, @@ -738,7 +743,9 @@ void Force_Stress_LCAO::calForceStressIntegralPart( #else ModuleBase::matrix& svl_dphi, #endif - LCAO_Hamilt &uhm, + LCAO_Hamilt &uhm, + Gint_Gamma &gint_gamma, // mohan add 2024-04-01 + Gint_k &gint_k, // mohan add 2024-04-01 Parallel_Orbitals &pv, LCAO_Matrix &lm, const K_Vectors& kv) @@ -763,13 +770,14 @@ void Force_Stress_LCAO::calForceStressIntegralPart( svl_dphi, #endif uhm, + gint_gamma, lm); return; } template<> -void Force_Stress_LCAO>::calForceStressIntegralPart( +void Force_Stress_LCAO>::integral_part( const bool isGammaOnly, const bool isforce, const bool isstress, @@ -790,9 +798,11 @@ void Force_Stress_LCAO>::calForceStressIntegralPart( ModuleBase::matrix& svl_dphi, #endif LCAO_Hamilt &uhm, + Gint_Gamma &gint_gamma, + Gint_k &gint_k, Parallel_Orbitals &pv, LCAO_Matrix &lm, - const K_Vectors& kv) + const K_Vectors& kv) { flk.ftable_k(isforce, isstress, @@ -814,6 +824,7 @@ void Force_Stress_LCAO>::calForceStressIntegralPart( svl_dphi, #endif uhm, + gint_k, pv, lm, kv); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h index 8aaf0d4a387..1604fd76cc8 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h @@ -13,6 +13,8 @@ #ifdef __EXX #include "module_ri/Exx_LRI.h" #endif +#include "module_hamilt_lcao/module_gint/gint_gamma.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" template class Force_Stress_LCAO @@ -35,7 +37,9 @@ class Force_Stress_LCAO LCAO_Matrix &lm, const elecstate::ElecState* pelec, const psi::Psi* psi, - LCAO_Hamilt& uhm, + LCAO_Hamilt& uhm, + Gint_Gamma &gint_gamma, // mohan add 2024-04-01 + Gint_k &gint_k, // mohan add 2024-04-01 ModuleBase::matrix& fcs, ModuleBase::matrix& scs, const Structure_Factor& sf, @@ -68,7 +72,8 @@ class Force_Stress_LCAO ModulePW::PW_Basis* rhopw, const Structure_Factor& sf); - void calForceStressIntegralPart(const bool isGammaOnly, + void integral_part( + const bool isGammaOnly, const bool isforce, const bool isstress, Local_Orbital_Charge& loc, @@ -88,6 +93,8 @@ class Force_Stress_LCAO ModuleBase::matrix& svl_dphi, #endif LCAO_Hamilt &uhm, + Gint_Gamma &gint_gamma, + Gint_k &gint_k, Parallel_Orbitals &pv, LCAO_Matrix &lm, const K_Vectors& kv); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp index 2bc04ebd51b..0b5835446b1 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp @@ -39,6 +39,7 @@ void Force_LCAO_gamma::ftable_gamma(const bool isforce, ModuleBase::matrix& svl_dphi, #endif LCAO_Hamilt &uhm, + Gint_Gamma &gint_gamma, LCAO_Matrix &lm) { ModuleBase::TITLE("Force_LCAO_gamma", "ftable"); @@ -50,11 +51,10 @@ void Force_LCAO_gamma::ftable_gamma(const bool isforce, this->ParaV = DM->get_paraV_pointer(); //const Parallel_Orbitals* pv = loc.ParaV; - this->UHM = &uhm; // allocate DSloc_x, DSloc_y, DSloc_z // allocate DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z - this->allocate_gamma(*this->ParaV, lm); + this->allocate_gamma(*this->ParaV, uhm.genH, lm); // calculate the 'energy density matrix' here. this->cal_foverlap(isforce, isstress, psid, pelec, lm, foverlap, soverlap); @@ -65,7 +65,7 @@ void Force_LCAO_gamma::ftable_gamma(const bool isforce, this->cal_ftvnl_dphi(DM, lm, isforce, isstress, ftvnl_dphi, stvnl_dphi); this->cal_fvnl_dbeta(DM, isforce, isstress, fvnl_dbeta, svnl_dbeta); - this->cal_fvl_dphi(loc.DM, isforce, isstress, pelec->pot, fvl_dphi, svl_dphi); + this->cal_fvl_dphi(loc.DM, isforce, isstress, pelec->pot, gint_gamma, fvl_dphi, svl_dphi); // caoyu add for DeePKS #ifdef __DEEPKS @@ -128,6 +128,7 @@ void Force_LCAO_gamma::ftable_gamma(const bool isforce, void Force_LCAO_gamma::allocate_gamma( const Parallel_Orbitals &pv, + LCAO_gen_fixedH &gen_h, LCAO_Matrix &lm) { ModuleBase::TITLE("Force_LCAO_gamma", "allocate_gamma"); @@ -179,7 +180,7 @@ void Force_LCAO_gamma::allocate_gamma( } // calculate dS in LCAO basis // ModuleBase::timer::tick("Force_LCAO_gamma","build_S_new"); - this->UHM->genH.build_ST_new('S', cal_deri, GlobalC::ucell, lm.Sloc.data()); + gen_h.build_ST_new('S', cal_deri, GlobalC::ucell, lm.Sloc.data()); // ModuleBase::timer::tick("Force_LCAO_gamma","build_S_new"); // calculate dT in LCAP @@ -196,12 +197,12 @@ void Force_LCAO_gamma::allocate_gamma( // calculate dT // calculate T + VNL(P1) in LCAO basis // ModuleBase::timer::tick("Force_LCAO_gamma","build_T_new"); - this->UHM->genH.build_ST_new('T', cal_deri, GlobalC::ucell, lm.Hloc_fixed.data()); + gen_h.build_ST_new('T', cal_deri, GlobalC::ucell, lm.Hloc_fixed.data()); // ModuleBase::timer::tick("Force_LCAO_gamma","build_T_new"); // test_gamma(lm.DHloc_fixed_x, "dHloc_fixed_x T part"); // ModuleBase::timer::tick("Force_LCAO_gamma","build_Nonlocal_mu"); - this->UHM->genH.build_Nonlocal_mu_new(lm.Hloc_fixed.data(), cal_deri); + gen_h.build_Nonlocal_mu_new(lm.Hloc_fixed.data(), cal_deri); // ModuleBase::timer::tick("Force_LCAO_gamma","build_Nonlocal_mu"); // test_gamma(lm.DHloc_fixed_x, "dHloc_fixed_x Vnl part"); @@ -212,7 +213,7 @@ void Force_LCAO_gamma::allocate_gamma( lm.zeros_HSgamma('S'); - this->UHM->genH.build_ST_new('S', cal_deri, GlobalC::ucell, + gen_h.build_ST_new('S', cal_deri, GlobalC::ucell, lm.Sloc.data(), INPUT.cal_syns, INPUT.dmax); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h index a3358382ee7..c8fe00c040c 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h @@ -9,6 +9,7 @@ #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h" #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" #include "module_psi/psi.h" +#include "module_hamilt_lcao/module_gint/gint_gamma.h" class Force_LCAO_gamma { @@ -21,8 +22,6 @@ class Force_LCAO_gamma private: - LCAO_Hamilt* UHM; - const Parallel_Orbitals* ParaV; elecstate::Potential* pot; @@ -46,14 +45,16 @@ class Force_LCAO_gamma #else ModuleBase::matrix& svl_dphi, #endif - LCAO_Hamilt& uhm, - LCAO_Matrix& lm); + LCAO_Hamilt& uhm, + Gint_Gamma &gint_gamma, + LCAO_Matrix& lm); // get the ds, dt, dvnl. void allocate_gamma( const Parallel_Orbitals& pv, + LCAO_gen_fixedH &gen_h, LCAO_Matrix &lm); void finish_ftable_gamma(LCAO_Matrix &lm); @@ -98,6 +99,7 @@ class Force_LCAO_gamma const bool isforce, const bool isstress, const elecstate::Potential* pot_in, + Gint_Gamma &gint_gamma, // mohan add 2024-04-01 ModuleBase::matrix& fvl_dphi, ModuleBase::matrix& svl_dphi); }; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_vl.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_vl.cpp index 6e97123cd07..5ce367c5682 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_vl.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_vl.cpp @@ -7,6 +7,7 @@ void Force_LCAO_gamma::cal_fvl_dphi( const bool isforce, const bool isstress, const elecstate::Potential* pot_in, + Gint_Gamma &gint_gamma, ModuleBase::matrix& fvl_dphi, ModuleBase::matrix& svl_dphi) { @@ -27,13 +28,22 @@ void Force_LCAO_gamma::cal_fvl_dphi( if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) { - Gint_inout inout(DM_in, is, vr_eff1, vofk_eff1, isforce, isstress, &fvl_dphi, &svl_dphi, Gint_Tools::job_type::force_meta); - this->UHM->GG.cal_gint(&inout); + Gint_inout inout(DM_in, + is, + vr_eff1, + vofk_eff1, + isforce, + isstress, + &fvl_dphi, + &svl_dphi, + Gint_Tools::job_type::force_meta); + + gint_gamma.cal_gint(&inout); } else { Gint_inout inout(DM_in, is, vr_eff1, isforce, isstress, &fvl_dphi, &svl_dphi, Gint_Tools::job_type::force); - this->UHM->GG.cal_gint(&inout); + gint_gamma.cal_gint(&inout); } } @@ -50,4 +60,4 @@ void Force_LCAO_gamma::cal_fvl_dphi( } } ModuleBase::timer::tick("Force_LCAO_gamma","cal_fvl_dphi"); -} \ No newline at end of file +} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp index 947d5bc5be8..98fd9379793 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp @@ -51,6 +51,7 @@ void Force_LCAO_k::ftable_k(const bool isforce, ModuleBase::matrix& svl_dphi, #endif LCAO_Hamilt &uhm, + Gint_k &gint_k, Parallel_Orbitals &pv, LCAO_Matrix &lm, const K_Vectors& kv) @@ -58,27 +59,61 @@ void Force_LCAO_k::ftable_k(const bool isforce, ModuleBase::TITLE("Force_LCAO_k", "ftable_k"); ModuleBase::timer::tick("Force_LCAO_k", "ftable_k"); - this->UHM = &uhm; - elecstate::DensityMatrix,double>* DM = dynamic_cast>*>(pelec)->get_DM(); - this->allocate_k(pv, lm, kv.nks, kv.kvec_d); + this->allocate_k( + pv, + lm, + uhm.genH, + kv.nks, + kv.kvec_d); // calculate the energy density matrix // and the force related to overlap matrix and energy density matrix. - this->cal_foverlap_k(isforce, isstress, ra, psi, loc, pv, lm, DM, foverlap, soverlap, pelec, kv.nks, kv); - - //DM->sum_DMR_spin(); - - this->cal_ftvnl_dphi_k(DM, pv, lm, isforce, isstress, ra, ftvnl_dphi, stvnl_dphi); - - // --------------------------------------- - // doing on the real space grid. - // --------------------------------------- - this->cal_fvl_dphi_k(isforce, isstress, lm, pelec->pot, fvl_dphi, svl_dphi, loc.DM_R); - - this->cal_fvnl_dbeta_k(DM, isforce, isstress, pv, fvnl_dbeta, svnl_dbeta); + this->cal_foverlap_k( + isforce, + isstress, + ra, + psi, + loc, + pv, + lm, + DM, + foverlap, + soverlap, + pelec, + kv.nks, + kv); + + this->cal_ftvnl_dphi_k( + DM, + pv, + lm, + isforce, + isstress, + ra, + ftvnl_dphi, + stvnl_dphi); + + // doing on the real space grid. + this->cal_fvl_dphi_k( + isforce, + isstress, + lm, + gint_k, + pelec->pot, + fvl_dphi, + svl_dphi, + loc.DM_R); + + this->cal_fvnl_dbeta_k( + DM, + isforce, + isstress, + pv, + fvnl_dbeta, + svnl_dbeta); #ifdef __DEEPKS if (GlobalV::deepks_scf) @@ -149,6 +184,7 @@ void Force_LCAO_k::ftable_k(const bool isforce, void Force_LCAO_k::allocate_k(const Parallel_Orbitals& pv, LCAO_Matrix &lm, + LCAO_gen_fixedH &gen_h, const int& nks, const std::vector>& kvec_d) { @@ -211,7 +247,7 @@ void Force_LCAO_k::allocate_k(const Parallel_Orbitals& pv, // calculate dS = //----------------------------- bool cal_deri = true; - this->UHM->genH.build_ST_new('S', cal_deri, GlobalC::ucell, this->UHM->genH.LM->SlocR.data()); + gen_h.build_ST_new('S', cal_deri, GlobalC::ucell, gen_h.LM->SlocR.data()); //----------------------------------------- // (2) allocate for @@ -222,7 +258,8 @@ void Force_LCAO_k::allocate_k(const Parallel_Orbitals& pv, // mohan add lm on 2024-03-31 const auto init_DHloc_fixedR_xyz = [this, nnr, &lm](int num_threads, int thread_id) { - int beg, len; + int beg=0; + int len=0; ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, nnr, 1024, beg, len); ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixedR_x + beg, len); ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixedR_y + beg, len); @@ -233,18 +270,18 @@ void Force_LCAO_k::allocate_k(const Parallel_Orbitals& pv, // calculate dT= in LCAO // calculate T + VNL(P1) in LCAO basis - this->UHM->genH.build_ST_new('T', cal_deri, GlobalC::ucell, this->UHM->genH.LM->Hloc_fixedR.data()); + gen_h.build_ST_new('T', cal_deri, GlobalC::ucell, gen_h.LM->Hloc_fixedR.data()); // calculate dVnl= in LCAO - this->UHM->genH.build_Nonlocal_mu_new(this->UHM->genH.LM->Hloc_fixed.data(), cal_deri); + gen_h.build_Nonlocal_mu_new(gen_h.LM->Hloc_fixed.data(), cal_deri); // calculate asynchronous S matrix to output for Hefei-NAMD if (INPUT.cal_syns) { cal_deri = false; - // this->UHM->genH.build_ST_new('S', cal_deri, GlobalC::ucell, this->UHM->genH.LM->SlocR.data(), + // gen_h.build_ST_new('S', cal_deri, GlobalC::ucell, gen_h.LM->SlocR.data(), // INPUT.cal_syns); - this->UHM->genH.build_ST_new('S', + gen_h.build_ST_new('S', cal_deri, GlobalC::ucell, lm.SlocR.data(), @@ -309,371 +346,6 @@ void Force_LCAO_k::finish_k(LCAO_Matrix &lm) return; } -#include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h" -void Force_LCAO_k::cal_foverlap_k(const bool isforce, - const bool isstress, - Record_adj &ra, - const psi::Psi> *psi, - Local_Orbital_Charge &loc, - Parallel_Orbitals &pv, - LCAO_Matrix &lm, - const elecstate::DensityMatrix, double> *DM, - ModuleBase::matrix &foverlap, - ModuleBase::matrix &soverlap, - const elecstate::ElecState *pelec, - const int &nks, - const K_Vectors &kv) -{ - ModuleBase::TITLE("Force_LCAO_k", "cal_foverlap_k"); - ModuleBase::timer::tick("Force_LCAO_k", "cal_foverlap_k"); - - // construct a DensityMatrix object - elecstate::DensityMatrix, double> EDM(&kv,&pv,GlobalV::NSPIN); - - //-------------------------------------------- - // calculate the energy density matrix here. - //-------------------------------------------- - ModuleBase::timer::tick("Force_LCAO_k", "cal_edm_2d"); - - ModuleBase::matrix wgEkb; - wgEkb.create(nks, GlobalV::NBANDS); - ModuleBase::Memory::record("Force::wgEkb", sizeof(double) * nks * GlobalV::NBANDS); -#ifdef _OPENMP -#pragma omp parallel for collapse(2) schedule(static, 1024) -#endif - for (int ik = 0; ik < nks; ik++) - { - for (int ib = 0; ib < GlobalV::NBANDS; ib++) - { - wgEkb(ik, ib) = pelec->wg(ik, ib) * pelec->ekb(ik, ib); - } - } - std::vector edm_k(nks); - - // use the original formula (Hamiltonian matrix) to calculate energy density matrix - if (DM->EDMK.size()) - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 1024) -#endif - for (int ik = 0; ik < nks; ++ik) - { - //edm_k[ik] = loc.edm_k_tddft[ik]; - EDM.set_DMK_pointer(ik,DM->EDMK[ik].c); - } - } - else - { - //elecstate::cal_dm(pv, wgEkb, psi[0], edm_k); - // cal_dm_psi - elecstate::cal_dm_psi(EDM.get_paraV_pointer(), wgEkb, psi[0], EDM); - } - - //loc.cal_dm_R(edm_k, ra, edm2d, kv); - - // cal_dm_2d - EDM.init_DMR(ra,&GlobalC::ucell); - EDM.cal_DMR(); - EDM.sum_DMR_spin(); - // - ModuleBase::timer::tick("Force_LCAO_k", "cal_edm_2d"); - //-------------------------------------------- - // summation \sum_{i,j} E(i,j)*dS(i,j) - // BEGIN CALCULATION OF FORCE OF EACH ATOM - //-------------------------------------------- - int total_irr = 0; -#ifdef _OPENMP -#pragma omp parallel - { - int num_threads = omp_get_num_threads(); - ModuleBase::matrix local_soverlap(3, 3); - int local_total_irr = 0; -#else - ModuleBase::matrix& local_soverlap = soverlap; - int& local_total_irr = total_irr; -#endif - - ModuleBase::Vector3 tau1, dtau, tau2; - -#ifdef _OPENMP -#pragma omp for schedule(dynamic) -#endif - for (int iat = 0; iat < GlobalC::ucell.nat; iat++) - { - const int T1 = GlobalC::ucell.iat2it[iat]; - Atom* atom1 = &GlobalC::ucell.atoms[T1]; - const int I1 = GlobalC::ucell.iat2ia[iat]; - // get iat1 - int iat1 = GlobalC::ucell.itia2iat(T1, I1); - double* foverlap_iat; - if (isforce) - { - foverlap_iat = &foverlap(iat, 0); - } - -#ifdef _OPENMP - // using local stack to avoid false sharing in multi-threaded case - double foverlap_temp[3] = {0.0, 0.0, 0.0}; - if (num_threads > 1) - { - foverlap_iat = foverlap_temp; - } -#endif - int irr = pv.nlocstart[iat]; - const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0); - for (int cb = 0; cb < ra.na_each[iat]; ++cb) - { - const int T2 = ra.info[iat][cb][3]; - const int I2 = ra.info[iat][cb][4]; - - const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0); - - Atom* atom2 = &GlobalC::ucell.atoms[T2]; - - // get iat2 - int iat2 = GlobalC::ucell.itia2iat(T2, I2); - double Rx = ra.info[iat][cb][0]; - double Ry = ra.info[iat][cb][1]; - double Rz = ra.info[iat][cb][2]; - // get BaseMatrix - hamilt::BaseMatrix* tmp_matrix = EDM.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz); - if(tmp_matrix == nullptr) - { - continue; - } - int row_ap = pv.atom_begin_row[iat1]; - int col_ap = pv.atom_begin_col[iat2]; - // get DMR - for (int mu = 0; mu < pv.get_row_size(iat1); ++mu) - { - for (int nu = 0; nu < pv.get_col_size(iat2); ++nu) - { - // here do not sum over spin due to EDM.sum_DMR_spin(); - double edm2d1 = tmp_matrix->get_value(mu,nu); - double edm2d2 = 2.0 * edm2d1; - - if (isforce) - { - foverlap_iat[0] -= edm2d2 * lm.DSloc_Rx[irr]; - foverlap_iat[1] -= edm2d2 * lm.DSloc_Ry[irr]; - foverlap_iat[2] -= edm2d2 * lm.DSloc_Rz[irr]; - } - if (isstress) - { - for (int ipol = 0; ipol < 3; ipol++) - { - local_soverlap(0, ipol) += edm2d1 * lm.DSloc_Rx[irr] - * lm.DH_r[irr * 3 + ipol]; - if (ipol < 1) - { - continue; - } - local_soverlap(1, ipol) += edm2d1 * lm.DSloc_Ry[irr] - * lm.DH_r[irr * 3 + ipol]; - if (ipol < 2) - { - continue; - } - local_soverlap(2, ipol) += edm2d1 * lm.DSloc_Rz[irr] - * lm.DH_r[irr * 3 + ipol]; - } - } - //} - ++local_total_irr; - ++irr; - } // end kk - } // end jj - } // end cb -#ifdef _OPENMP - if (isforce && num_threads > 1) - { - foverlap(iat, 0) += foverlap_iat[0]; - foverlap(iat, 1) += foverlap_iat[1]; - foverlap(iat, 2) += foverlap_iat[2]; - } -#endif - } // end iat -#ifdef _OPENMP -#pragma omp critical(cal_foverlap_k_reduce) - { - total_irr += local_total_irr; - if (isstress) - { - for (int ipol = 0; ipol < 3; ipol++) - { - soverlap(0, ipol) += local_soverlap(0, ipol); - if (ipol < 1) - { - continue; - } - soverlap(1, ipol) += local_soverlap(1, ipol); - if (ipol < 2) - { - continue; - } - soverlap(2, ipol) += local_soverlap(2, ipol); - } - } - } - } -#endif - - if (isstress) - { - StressTools::stress_fill(GlobalC::ucell.lat0, GlobalC::ucell.omega, soverlap); - } - - if (total_irr != pv.nnr) - { - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wrong irr", total_irr); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wrong LNNR.nnr", pv.nnr); - ModuleBase::WARNING_QUIT("Force_LCAO_k::cal_foverlap_k", "irr!=LNNR.nnr"); - } - - ModuleBase::timer::tick("Force_LCAO_k", "cal_foverlap_k"); - return; -} - -void Force_LCAO_k::cal_ftvnl_dphi_k(const elecstate::DensityMatrix, double>* DM, - const Parallel_Orbitals &pv, - LCAO_Matrix &lm, - const bool isforce, - const bool isstress, - Record_adj& ra, - ModuleBase::matrix& ftvnl_dphi, - ModuleBase::matrix& stvnl_dphi) -{ - ModuleBase::TITLE("Force_LCAO_k", "cal_ftvnl_dphi_k"); - ModuleBase::timer::tick("Force_LCAO_k", "cal_ftvnl_dphi_k"); - - int total_irr = 0; - // get the adjacent atom's information. - - // GlobalV::ofs_running << " calculate the ftvnl_dphi_k force" << std::endl; -#ifdef _OPENMP -#pragma omp parallel - { - int num_threads = omp_get_num_threads(); - ModuleBase::matrix local_stvnl_dphi(3, 3); - int local_total_irr = 0; -#pragma omp for schedule(dynamic) -#else - ModuleBase::matrix& local_stvnl_dphi = stvnl_dphi; - int& local_total_irr = total_irr; -#endif - for (int iat = 0; iat < GlobalC::ucell.nat; iat++) - { - const int T1 = GlobalC::ucell.iat2it[iat]; - Atom* atom1 = &GlobalC::ucell.atoms[T1]; - const int I1 = GlobalC::ucell.iat2ia[iat]; - // get iat1 - int iat1 = GlobalC::ucell.itia2iat(T1, I1); - // - int irr = pv.nlocstart[iat]; - const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0); - double* ftvnl_dphi_iat; - if (isforce) - { - ftvnl_dphi_iat = &ftvnl_dphi(iat, 0); - } -#ifdef _OPENMP - // using local stack to avoid false sharing in multi-threaded case - double ftvnl_dphi_temp[3] = {0.0, 0.0, 0.0}; - if (num_threads > 1) - { - ftvnl_dphi_iat = ftvnl_dphi_temp; - } -#endif - for (int cb = 0; cb < ra.na_each[iat]; ++cb) - { - const int T2 = ra.info[iat][cb][3]; - const int I2 = ra.info[iat][cb][4]; - const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0); - Atom* atom2 = &GlobalC::ucell.atoms[T2]; - // get iat2 - int iat2 = GlobalC::ucell.itia2iat(T2, I2); - double Rx = ra.info[iat][cb][0]; - double Ry = ra.info[iat][cb][1]; - double Rz = ra.info[iat][cb][2]; - // get BaseMatrix - if (pv.get_row_size(iat1) <= 0 || pv.get_col_size(iat2) <= 0) - { - continue; - } - std::vector*> tmp_matrix; - for (int is = 0; is < GlobalV::NSPIN; ++is) - { - tmp_matrix.push_back(DM->get_DMR_pointer(is+1)->find_matrix(iat1, iat2, Rx, Ry, Rz)); - } - //hamilt::BaseMatrix* tmp_matrix = DM->get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz); - for (int mu = 0; mu < pv.get_row_size(iat1); ++mu) - { - for (int nu = 0; nu < pv.get_col_size(iat2); ++nu) - { - // get value from DM - double dm2d1 = 0.0; - for (int is = 0; is < GlobalV::NSPIN; ++is) - { - dm2d1 += tmp_matrix[is]->get_value(mu, nu); - } - double dm2d2 = 2.0 * dm2d1; - // - if (isforce) - { - ftvnl_dphi_iat[0] += dm2d2 * lm.DHloc_fixedR_x[irr]; - ftvnl_dphi_iat[1] += dm2d2 * lm.DHloc_fixedR_y[irr]; - ftvnl_dphi_iat[2] += dm2d2 * lm.DHloc_fixedR_z[irr]; - } - if (isstress) - { - local_stvnl_dphi(0, 0) -= dm2d1 * lm.stvnl11[irr]; - local_stvnl_dphi(0, 1) -= dm2d1 * lm.stvnl12[irr]; - local_stvnl_dphi(0, 2) -= dm2d1 * lm.stvnl13[irr]; - local_stvnl_dphi(1, 1) -= dm2d1 * lm.stvnl22[irr]; - local_stvnl_dphi(1, 2) -= dm2d1 * lm.stvnl23[irr]; - local_stvnl_dphi(2, 2) -= dm2d1 * lm.stvnl33[irr]; - } - //} - ++local_total_irr; - ++irr; - } // end kk - } // end jj - } // end cb -#ifdef _OPENMP - if (isforce && num_threads > 1) - { - ftvnl_dphi(iat, 0) += ftvnl_dphi_iat[0]; - ftvnl_dphi(iat, 1) += ftvnl_dphi_iat[1]; - ftvnl_dphi(iat, 2) += ftvnl_dphi_iat[2]; - } -#endif - } // end iat -#ifdef _OPENMP -#pragma omp critical(cal_ftvnl_dphi_k_reduce) - { - total_irr += local_total_irr; - if (isstress) - { - stvnl_dphi(0, 0) += local_stvnl_dphi(0, 0); - stvnl_dphi(0, 1) += local_stvnl_dphi(0, 1); - stvnl_dphi(0, 2) += local_stvnl_dphi(0, 2); - stvnl_dphi(1, 1) += local_stvnl_dphi(1, 1); - stvnl_dphi(1, 2) += local_stvnl_dphi(1, 2); - stvnl_dphi(2, 2) += local_stvnl_dphi(2, 2); - } - } - } -#endif - assert(total_irr == pv.nnr); - - if (isstress) - { - StressTools::stress_fill(GlobalC::ucell.lat0, GlobalC::ucell.omega, stvnl_dphi); - } - - ModuleBase::timer::tick("Force_LCAO_k", "cal_ftvnl_dphi_k"); - return; -} void Force_LCAO_k::test( Parallel_Orbitals &pv, @@ -737,11 +409,15 @@ void Force_LCAO_k::test( { for (int j = 0; j < GlobalV::NLOCAL; j++) { - if (std::abs(test[i * GlobalV::NLOCAL + j]) > 1.0e-5) - std::cout << std::setw(12) << test[i * GlobalV::NLOCAL + j]; - else - std::cout << std::setw(12) << "0"; - } + if (std::abs(test[i * GlobalV::NLOCAL + j]) > 1.0e-5) + { + std::cout << std::setw(12) << test[i * GlobalV::NLOCAL + j]; + } + else + { + std::cout << std::setw(12) << "0"; + } + } std::cout << std::endl; } delete[] test; @@ -750,563 +426,5 @@ void Force_LCAO_k::test( return; } -typedef std::tuple key_tuple; - -// must consider three-center H matrix. -void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix, double>* DM, - const bool isforce, - const bool isstress, - const Parallel_Orbitals &pv, - ModuleBase::matrix& fvnl_dbeta, - ModuleBase::matrix& svnl_dbeta) -{ - ModuleBase::TITLE("Force_LCAO_k", "cal_fvnl_dbeta_k_new"); - ModuleBase::timer::tick("Force_LCAO_k", "cal_fvnl_dbeta_k_new"); - - // Data structure for storing , for a detailed description - // check out the same data structure in build_Nonlocal_mu_new - std::vector>>>> nlm_tot; - nlm_tot.resize(GlobalC::ucell.nat); -#ifdef _OPENMP -// use schedule(dynamic) for load balancing because adj_num is various -#pragma omp parallel for schedule(dynamic) -#endif - for (int iat = 0; iat < GlobalC::ucell.nat; iat++) - { - - const int it = GlobalC::ucell.iat2it[iat]; - const int ia = GlobalC::ucell.iat2ia[iat]; - - // Step 1 : generate - // type of atom; distance; atomic basis; projectors - - const double Rcut_Beta = GlobalC::ucell.infoNL.Beta[it].get_rcut_max(); - const ModuleBase::Vector3 tau = GlobalC::ucell.atoms[it].tau[ia]; - AdjacentAtomInfo adjs; - GlobalC::GridD.Find_atom(GlobalC::ucell, tau, it, ia, &adjs); - - nlm_tot[iat].clear(); - - for (int ad = 0; ad < adjs.adj_num + 1; ++ad) - { - const int T1 = adjs.ntype[ad]; - const int I1 = adjs.natom[ad]; - const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0); - const double Rcut_AO1 = GlobalC::ORB.Phi[T1].getRcut(); - - const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad]; - const Atom* atom1 = &GlobalC::ucell.atoms[T1]; - const int nw1_tot = atom1->nw * GlobalV::NPOL; - - const ModuleBase::Vector3 dtau = tau1 - tau; - const double dist1 = dtau.norm2() * pow(GlobalC::ucell.lat0, 2); - if (dist1 > pow(Rcut_Beta + Rcut_AO1, 2)) - { - continue; - } - - std::unordered_map>> nlm_cur; - nlm_cur.clear(); - - for (int iw1 = 0; iw1 < nw1_tot; ++iw1) - { - const int iw1_all = start1 + iw1; - const int iw1_local = pv.global2local_row(iw1_all); - const int iw2_local = pv.global2local_col(iw1_all); - if (iw1_local < 0 && iw2_local < 0) - { - continue; - } - const int iw1_0 = iw1 / GlobalV::NPOL; - std::vector> nlm; -#ifdef USE_NEW_TWO_CENTER - //================================================================= - // new two-center integral (temporary) - //================================================================= - int L1 = atom1->iw2l[ iw1_0 ]; - int N1 = atom1->iw2n[ iw1_0 ]; - int m1 = atom1->iw2m[ iw1_0 ]; - - // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l) - int M1 = (m1 % 2 == 0) ? -m1/2 : (m1+1)/2; - - ModuleBase::Vector3 dtau = tau - tau1; - GlobalC::UOT.two_center_bundle->overlap_orb_beta->snap( - T1, L1, N1, M1, it, dtau * GlobalC::ucell.lat0, true, nlm); -#else - GlobalC::UOT.snap_psibeta_half(GlobalC::ORB, - GlobalC::ucell.infoNL, - nlm, - tau1, - T1, - atom1->iw2l[iw1_0], // L1 - atom1->iw2m[iw1_0], // m1 - atom1->iw2n[iw1_0], // N1 - tau, - it, - 1); // R0,T0 -#endif - nlm_cur.insert({iw1_all, nlm}); - } // end iw - const int iat1 = GlobalC::ucell.itia2iat(T1, I1); - const int rx1 = adjs.box[ad].x; - const int ry1 = adjs.box[ad].y; - const int rz1 = adjs.box[ad].z; - key_tuple key_1(iat1, rx1, ry1, rz1); - nlm_tot[iat][key_1] = nlm_cur; - } // end ad - } - - //======================================================= - // Step2: - // calculate sum_(L0,M0) beta - // and accumulate the value to Hloc_fixedR(i,j) - //======================================================= - int total_nnr = 0; -#ifdef _OPENMP -#pragma omp parallel reduction(+ : total_nnr) - { - ModuleBase::matrix local_svnl_dbeta(3, 3); - const int num_threads = omp_get_num_threads(); -#else - ModuleBase::matrix& local_svnl_dbeta = svnl_dbeta; -#endif - - ModuleBase::Vector3 tau1; - ModuleBase::Vector3 tau2; - ModuleBase::Vector3 dtau; - ModuleBase::Vector3 tau0; - ModuleBase::Vector3 dtau1; - ModuleBase::Vector3 dtau2; - - double rcut; - double distance; - - double rcut1; - double rcut2; - double distance1; - double distance2; - -#ifdef _OPENMP -// use schedule(dynamic) for load balancing because adj_num is various -#pragma omp for schedule(dynamic) -#endif - for (int iat1 = 0; iat1 < GlobalC::ucell.nat; iat1++) - { - const int T1 = GlobalC::ucell.iat2it[iat1]; - const Atom* atom1 = &GlobalC::ucell.atoms[T1]; - - { - const int I1 = GlobalC::ucell.iat2ia[iat1]; - tau1 = atom1->tau[I1]; - AdjacentAtomInfo adjs; - GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1, &adjs); - const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0); - int nnr = pv.nlocstart[iat1]; - - /* - !!!!!!!!!!!!!!!! - This optimization is also improving the performance of single thread. - Making memory access more linearly in the core loop - */ - bool iat_recorded = false; - bool force_updated = false; - // record iat of adjs - std::vector adj_iat; - // record fvnl_dbeta diff of adjs - std::vector adj_fvnl_dbeta; - if (isforce) - { - adj_iat.resize(adjs.adj_num + 1); - adj_fvnl_dbeta.resize((adjs.adj_num + 1) * 3, 0.0); - } - - for (int ad2 = 0; ad2 < adjs.adj_num + 1; ++ad2) - { - const int T2 = adjs.ntype[ad2]; - const Atom* atom2 = &GlobalC::ucell.atoms[T2]; - - const int I2 = adjs.natom[ad2]; - const int iat2 = GlobalC::ucell.itia2iat(T2, I2); - const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0); - tau2 = adjs.adjacent_tau[ad2]; - - const int rx2 = adjs.box[ad2].x; - const int ry2 = adjs.box[ad2].y; - const int rz2 = adjs.box[ad2].z; - - dtau = tau2 - tau1; - distance = dtau.norm2() * pow(GlobalC::ucell.lat0, 2); - // this rcut is in order to make nnr consistent - // with other matrix. - rcut = pow(GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ORB.Phi[T2].getRcut(), 2); - - // check if this a adjacent atoms. - bool is_adj = false; - if (distance < rcut) - is_adj = true; - else if (distance >= rcut) - { - for (int ad0 = 0; ad0 < adjs.adj_num + 1; ++ad0) - { - const int T0 = adjs.ntype[ad0]; - if (GlobalC::ucell.infoNL.nproj[T0] == 0) - { - continue; - } - const int I0 = adjs.natom[ad0]; - // const int iat0 = GlobalC::ucell.itia2iat(T0, I0); - // const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); - - tau0 = adjs.adjacent_tau[ad0]; - dtau1 = tau0 - tau1; - distance1 = dtau1.norm() * GlobalC::ucell.lat0; - rcut1 = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); - - dtau2 = tau0 - tau2; - distance2 = dtau2.norm() * GlobalC::ucell.lat0; - rcut2 = GlobalC::ORB.Phi[T2].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); - - if (distance1 < rcut1 && distance2 < rcut2) - { - is_adj = true; - break; - } - } - } - - if (is_adj) - { - // basematrix and its data pointer - if (pv.get_row_size(iat1) <= 0 || pv.get_col_size(iat2) <= 0) - { - continue; - } - std::vector tmp_matrix_ptr; - for (int is = 0; is < GlobalV::NSPIN; ++is) - { - auto* tmp_base_matrix = DM->get_DMR_pointer(is+1)->find_matrix(iat1, iat2, rx2, ry2, rz2); - tmp_matrix_ptr.push_back(tmp_base_matrix->get_pointer()); - } - //hamilt::BaseMatrix* tmp_matrix = DM->get_DMR_pointer(1)->find_matrix(iat1, iat2, rx2, ry2, rz2); - //double* tmp_matrix_ptr = tmp_matrix->get_pointer(); - for (int ad0 = 0; ad0 < adjs.adj_num + 1; ++ad0) - { - const int T0 = adjs.ntype[ad0]; - const int I0 = adjs.natom[ad0]; - const int iat = GlobalC::ucell.itia2iat(T0, I0); - if (!iat_recorded && isforce) - adj_iat[ad0] = iat; - - // mohan add 2010-12-19 - if (GlobalC::ucell.infoNL.nproj[T0] == 0) - continue; - - // const int I0 = GlobalC::GridD.getNatom(ad0); - // const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); - tau0 = adjs.adjacent_tau[ad0]; - - dtau1 = tau0 - tau1; - dtau2 = tau0 - tau2; - const double distance1 = dtau1.norm2() * pow(GlobalC::ucell.lat0, 2); - const double distance2 = dtau2.norm2() * pow(GlobalC::ucell.lat0, 2); - - // seems a bug here!! mohan 2011-06-17 - rcut1 = pow(GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(), - 2); - rcut2 = pow(GlobalC::ORB.Phi[T2].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(), - 2); - - double r0[3]; - double r1[3]; - r1[0] = (tau1.x - tau0.x); - r1[1] = (tau1.y - tau0.y); - r1[2] = (tau1.z - tau0.z); - r0[0] = (tau2.x - tau0.x); - r0[1] = (tau2.y - tau0.y); - r0[2] = (tau2.z - tau0.z); - - if (distance1 >= rcut1 || distance2 >= rcut2) - { - continue; - } - - const int rx0 = adjs.box[ad0].x; - const int ry0 = adjs.box[ad0].y; - const int rz0 = adjs.box[ad0].z; - key_tuple key1(iat1, -rx0, -ry0, -rz0); - key_tuple key2(iat2, rx2 - rx0, ry2 - ry0, rz2 - rz0); - - int nnr_inner = 0; - for (int j = 0; j < atom1->nw * GlobalV::NPOL; j++) - { - const int j0 = j / GlobalV::NPOL; // added by zhengdy-soc - const int iw1_all = start1 + j; - const int mu = pv.global2local_row(iw1_all); - if (mu < 0) - { - continue; - } - - for (int k = 0; k < atom2->nw * GlobalV::NPOL; k++) - { - const int k0 = k / GlobalV::NPOL; - const int iw2_all = start2 + k; - const int nu = pv.global2local_col(iw2_all); - if (nu < 0) - { - continue; - } - - // const Atom* atom0 = &GlobalC::ucell.atoms[T0]; - double nlm[3] = {0, 0, 0}; - std::vector nlm_1 = nlm_tot[iat][key2][iw2_all][0]; - std::vector> nlm_2; - nlm_2.resize(3); - for (int i = 0; i < 3; i++) - { - nlm_2[i] = nlm_tot[iat][key1][iw1_all][i + 1]; - } - - assert(nlm_1.size() == nlm_2[0].size()); - - const int nproj = GlobalC::ucell.infoNL.nproj[T0]; - int ib = 0; - for (int nb = 0; nb < nproj; nb++) - { - const int L0 = GlobalC::ucell.infoNL.Beta[T0].Proj[nb].getL(); - for (int m = 0; m < 2 * L0 + 1; m++) - { - for (int ir = 0; ir < 3; ir++) - { - nlm[ir] += nlm_2[ir][ib] * nlm_1[ib] - * GlobalC::ucell.atoms[T0].ncpp.dion(nb, nb); - } - ib += 1; - } - } - assert(ib == nlm_1.size()); - - double nlm1[3] = {0, 0, 0}; - if (isstress) - { - std::vector nlm_1 = nlm_tot[iat][key1][iw1_all][0]; - std::vector> nlm_2; - nlm_2.resize(3); - for (int i = 0; i < 3; i++) - { - nlm_2[i] = nlm_tot[iat][key2][iw2_all][i + 1]; - } - - assert(nlm_1.size() == nlm_2[0].size()); - - const int nproj = GlobalC::ucell.infoNL.nproj[T0]; - int ib = 0; - for (int nb = 0; nb < nproj; nb++) - { - const int L0 = GlobalC::ucell.infoNL.Beta[T0].Proj[nb].getL(); - for (int m = 0; m < 2 * L0 + 1; m++) - { - for (int ir = 0; ir < 3; ir++) - { - nlm1[ir] += nlm_2[ir][ib] * nlm_1[ib] - * GlobalC::ucell.atoms[T0].ncpp.dion(nb, nb); - } - ib += 1; - } - } - assert(ib == nlm_1.size()); - } - - /// only one projector for each atom force, but another projector for stress - force_updated = true; - // get DMR - double dm2d1 = 0.0; - for (int is = 0; is < GlobalV::NSPIN; ++is) - { - dm2d1 += tmp_matrix_ptr[is][nnr_inner]; - } - double dm2d2 = 2.0 * dm2d1; - // - for (int jpol = 0; jpol < 3; jpol++) - { - if (isforce) - { - adj_fvnl_dbeta[ad0 * 3 + jpol] -= dm2d2 * nlm[jpol]; - } - if (isstress) - { - for (int ipol = jpol; ipol < 3; ipol++) - { - local_svnl_dbeta(jpol, ipol) - += dm2d1 - * (nlm[jpol] * r1[ipol] + nlm1[jpol] * r0[ipol]); - } - } - } - //} - nnr_inner++; - } // k - } // j - } // ad0 - - // outer circle : accumulate nnr - for (int j = 0; j < atom1->nw * GlobalV::NPOL; j++) - { - const int j0 = j / GlobalV::NPOL; // added by zhengdy-soc - const int iw1_all = start1 + j; - const int mu = pv.global2local_row(iw1_all); - if (mu < 0) - { - continue; - } - - // fix a serious bug: atom2[T2] -> atom2 - // mohan 2010-12-20 - for (int k = 0; k < atom2->nw * GlobalV::NPOL; k++) - { - const int k0 = k / GlobalV::NPOL; - const int iw2_all = start2 + k; - const int nu = pv.global2local_col(iw2_all); - if (nu < 0) - { - continue; - } - total_nnr++; - nnr++; - } - } - iat_recorded = true; - } // is_adj - } // ad2 - - // sum the diff to fvnl_dbeta - if (force_updated && isforce) - { -#ifdef _OPENMP - if (num_threads > 1) - { - for (int ad0 = 0; ad0 < adjs.adj_num + 1; ++ad0) - { -#pragma omp atomic - fvnl_dbeta(adj_iat[ad0], 0) += adj_fvnl_dbeta[ad0 * 3 + 0]; -#pragma omp atomic - fvnl_dbeta(adj_iat[ad0], 1) += adj_fvnl_dbeta[ad0 * 3 + 1]; -#pragma omp atomic - fvnl_dbeta(adj_iat[ad0], 2) += adj_fvnl_dbeta[ad0 * 3 + 2]; - } - } - else -#endif - { - for (int ad0 = 0; ad0 < adjs.adj_num + 1; ++ad0) - { - fvnl_dbeta(adj_iat[ad0], 0) += adj_fvnl_dbeta[ad0 * 3 + 0]; - fvnl_dbeta(adj_iat[ad0], 1) += adj_fvnl_dbeta[ad0 * 3 + 1]; - fvnl_dbeta(adj_iat[ad0], 2) += adj_fvnl_dbeta[ad0 * 3 + 2]; - } - } - } - } // I1 - } // T1 - -#ifdef _OPENMP - if (isstress) - { -#pragma omp critical(cal_fvnl_dbeta_k_new_reduce) - { - for (int l = 0; l < 3; l++) - { - for (int m = 0; m < 3; m++) - { - svnl_dbeta(l, m) += local_svnl_dbeta(l, m); - } - } - } - } - } -#endif - - assert(total_nnr == pv.nnr); - - if (isstress) - { - StressTools::stress_fill(GlobalC::ucell.lat0, GlobalC::ucell.omega, svnl_dbeta); - } - - ModuleBase::timer::tick("Force_LCAO_k", "cal_fvnl_dbeta_k_new"); - return; -} - -// calculate the force due to < phi | Vlocal | dphi > -void Force_LCAO_k::cal_fvl_dphi_k(const bool isforce, - const bool isstress, - LCAO_Matrix &lm, - const elecstate::Potential* pot_in, - ModuleBase::matrix& fvl_dphi, - ModuleBase::matrix& svl_dphi, - double** DM_R) -{ - ModuleBase::TITLE("Force_LCAO_k", "cal_fvl_dphi_k"); - ModuleBase::timer::tick("Force_LCAO_k", "cal_fvl_dphi_k"); - - if (!isforce && !isstress) - { - ModuleBase::timer::tick("Force_LCAO_k", "cal_fvl_dphi_k"); - return; - } - - assert(lm.DHloc_fixedR_x != NULL); - assert(lm.DHloc_fixedR_y != NULL); - assert(lm.DHloc_fixedR_z != NULL); - - int istep = 1; - - for (int is = 0; is < GlobalV::NSPIN; ++is) - { - GlobalV::CURRENT_SPIN = is; - - const double* vr_eff1 = pot_in->get_effective_v(GlobalV::CURRENT_SPIN); - const double* vofk_eff1 = nullptr; - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { - vofk_eff1 = pot_in->get_effective_vofk(GlobalV::CURRENT_SPIN); - } - - //-------------------------------- - // Grid integration here. - //-------------------------------- - // fvl_dphi can not be set to zero here if Vna is used - if (isstress || isforce) - { - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { - Gint_inout inout(DM_R, - is, - vr_eff1, - vofk_eff1, - isforce, - isstress, - &fvl_dphi, - &svl_dphi, - Gint_Tools::job_type::force_meta); - this->UHM->GK.cal_gint(&inout); - } - else - { - Gint_inout inout(DM_R, is, vr_eff1, isforce, isstress, &fvl_dphi, &svl_dphi, Gint_Tools::job_type::force); - this->UHM->GK.cal_gint(&inout); - } - } - } - - if (isstress) - { - StressTools::stress_fill(-1.0, GlobalC::ucell.omega, svl_dphi); - } - - ModuleBase::timer::tick("Force_LCAO_k", "cal_fvl_dphi_k"); - return; -} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h index 2d9b3bc2c30..9fa60e2ad5e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h @@ -9,6 +9,7 @@ #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h" #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" class Force_LCAO_k : public Force_LCAO_gamma { @@ -20,7 +21,6 @@ class Force_LCAO_k : public Force_LCAO_gamma ~Force_LCAO_k(); private: - LCAO_Hamilt* UHM; // orthonormal force + contribution from T and VNL void ftable_k(const bool isforce, @@ -43,6 +43,7 @@ class Force_LCAO_k : public Force_LCAO_gamma ModuleBase::matrix& svl_dphi, #endif LCAO_Hamilt &uhm, + Gint_k &gint_k, Parallel_Orbitals &pv, LCAO_Matrix &lm, const K_Vectors &kv); @@ -50,6 +51,7 @@ class Force_LCAO_k : public Force_LCAO_gamma // get the ds, dt, dvnl. void allocate_k(const Parallel_Orbitals& pv, LCAO_Matrix &lm, + LCAO_gen_fixedH &gen_h, const int& nks, const std::vector>& kvec_d); @@ -84,6 +86,7 @@ class Force_LCAO_k : public Force_LCAO_gamma void cal_fvl_dphi_k(const bool isforce, const bool isstress, LCAO_Matrix &lm, + Gint_k &gint_k, const elecstate::Potential* pot_in, ModuleBase::matrix& fvl_dphi, ModuleBase::matrix& svl_dphi, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp index 419f1726e71..b1881849f48 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp @@ -33,6 +33,8 @@ LCAO_Hamilt::~LCAO_Hamilt() //-------------------------------------------- void LCAO_Hamilt::grid_prepare( const Grid_Technique& gt, + Gint_Gamma &gint_gamma, + Gint_k &gint_k, const ModulePW::PW_Basis& rhopw, const ModulePW::PW_Basis_Big& bigpw) { @@ -41,7 +43,7 @@ void LCAO_Hamilt::grid_prepare( if(GlobalV::GAMMA_ONLY_LOCAL) { - this->GG.prep_grid(gt, bigpw.nbx, bigpw.nby, bigpw.nbzp, bigpw.nbzp_start, + gint_gamma.prep_grid(gt, bigpw.nbx, bigpw.nby, bigpw.nbzp, bigpw.nbzp_start, rhopw.nxyz, bigpw.bx, bigpw.by, bigpw.bz, bigpw.bxyz, bigpw.nbxx, rhopw.ny, rhopw.nplane, rhopw.startz_current); @@ -49,7 +51,7 @@ void LCAO_Hamilt::grid_prepare( else // multiple k-points { // cal the grid integration of 'Vl' matrix for l-points algorithms. - this->GK.prep_grid(gt, bigpw.nbx, bigpw.nby, bigpw.nbzp, bigpw.nbzp_start, + gint_k.prep_grid(gt, bigpw.nbx, bigpw.nby, bigpw.nbzp, bigpw.nbzp_start, rhopw.nxyz, bigpw.bx, bigpw.by, bigpw.bz, bigpw.bxyz, bigpw.nbxx, rhopw.ny, rhopw.nplane, rhopw.startz_current); } @@ -58,7 +60,7 @@ void LCAO_Hamilt::grid_prepare( return; } -void LCAO_Hamilt::set_R_range_sparse(void) +void LCAO_Hamilt::set_R_range_sparse(LCAO_Matrix &lm) { int R_minX = int(GlobalC::GridD.getD_minX()); int R_minY = int(GlobalC::GridD.getD_minY()); @@ -75,7 +77,7 @@ void LCAO_Hamilt::set_R_range_sparse(void) for(int iz = 0; iz < R_z; iz++) { Abfs::Vector3_Order temp_R(ix+R_minX, iy+R_minY, iz+R_minZ); - this->LM->all_R_coor.insert(temp_R); + lm.all_R_coor.insert(temp_R); } } } @@ -83,7 +85,13 @@ void LCAO_Hamilt::set_R_range_sparse(void) return; } -void LCAO_Hamilt::cal_STN_R_sparse_for_S(const double &sparse_threshold) +void LCAO_Hamilt::cal_STN_R_sparse_for_S( + const Parallel_Orbitals pv, + std::vector &slocR, + std::vector> &slocR_soc, + std::map, std::map>> &SR_sparse, + std::map, std::map>>> &SR_soc_sparse, + const double &sparse_threshold) { ModuleBase::TITLE("LCAO_Hamilt","cal_STN_R_sparse_for_S"); @@ -91,8 +99,8 @@ void LCAO_Hamilt::cal_STN_R_sparse_for_S(const double &sparse_threshold) ModuleBase::Vector3 dtau, tau1, tau2; ModuleBase::Vector3 dtau1, dtau2, tau0; - double temp_value_double = 0.0; - std::complex temp_value_complex; + double tmp=0.0; + std::complex tmpc=complex(0.0,0.0); for(int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1) { @@ -117,7 +125,10 @@ void LCAO_Hamilt::cal_STN_R_sparse_for_S(const double &sparse_threshold) bool adj = false; - if(distance < rcut) adj = true; + if(distance < rcut) + { + adj = true; + } else if(distance >= rcut) { for(int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0) @@ -146,36 +157,45 @@ void LCAO_Hamilt::cal_STN_R_sparse_for_S(const double &sparse_threshold) { const int start2 = GlobalC::ucell.itiaiw2iwt(T2,I2,0); - Abfs::Vector3_Order dR(GlobalC::GridD.getBox(ad).x, GlobalC::GridD.getBox(ad).y, GlobalC::GridD.getBox(ad).z); + Abfs::Vector3_Order dR( + GlobalC::GridD.getBox(ad).x, + GlobalC::GridD.getBox(ad).y, + GlobalC::GridD.getBox(ad).z); for(int ii=0; iinw*GlobalV::NPOL; ii++) { const int iw1_all = start + ii; - const int mu = this->LM->ParaV->global2local_row(iw1_all); + const int mu = pv.global2local_row(iw1_all); - if(mu<0)continue; + if(mu<0) + { + continue; + } for(int jj=0; jjnw*GlobalV::NPOL; jj++) { int iw2_all = start2 + jj; - const int nu = this->LM->ParaV->global2local_col(iw2_all); + const int nu = pv.global2local_col(iw2_all); - if(nu<0)continue; + if(nu<0) + { + continue; + } if(GlobalV::NSPIN!=4) { - temp_value_double = this->LM->SlocR[index]; - if (std::abs(temp_value_double) > sparse_threshold) + tmp = slocR[index]; + if (std::abs(tmp) > sparse_threshold) { - this->LM->SR_sparse[dR][iw1_all][iw2_all] = temp_value_double; + SR_sparse[dR][iw1_all][iw2_all] = tmp; } } else { - temp_value_complex = this->LM->SlocR_soc[index]; - if(std::abs(temp_value_complex) > sparse_threshold) + tmpc = slocR_soc[index]; + if(std::abs(tmpc) > sparse_threshold) { - this->LM->SR_soc_sparse[dR][iw1_all][iw2_all] = temp_value_complex; + SR_soc_sparse[dR][iw1_all][iw2_all] = tmpc; } } @@ -286,7 +306,7 @@ void LCAO_Hamilt::cal_HSR_sparse( { ModuleBase::TITLE("LCAO_Hamilt","cal_HSR_sparse"); - set_R_range_sparse(); + set_R_range_sparse(*this->LM); //cal_STN_R_sparse(current_spin, sparse_threshold); if(GlobalV::NSPIN!=4) @@ -352,11 +372,14 @@ void LCAO_Hamilt::cal_HSR_sparse( clear_zero_elements(current_spin, sparse_threshold); } -void LCAO_Hamilt::cal_dH_sparse(const int ¤t_spin, const double &sparse_threshold) +void LCAO_Hamilt::cal_dH_sparse( + const int ¤t_spin, + const double &sparse_threshold, + Gint_k &gint_k) { ModuleBase::TITLE("LCAO_Hamilt","cal_dH_sparse"); - set_R_range_sparse(); + set_R_range_sparse(*this->LM); const int nnr = this->LM->ParaV->nnr; this->LM->DHloc_fixedR_x = new double[nnr]; @@ -386,9 +409,10 @@ void LCAO_Hamilt::cal_dH_sparse(const int ¤t_spin, const double &sparse_th delete[] this->LM->DHloc_fixedR_y; delete[] this->LM->DHloc_fixedR_z; - GK.cal_dvlocal_R_sparseMatrix(current_spin, sparse_threshold, this->LM); + gint_k.cal_dvlocal_R_sparseMatrix(current_spin, sparse_threshold, this->LM); } + void LCAO_Hamilt::cal_STN_R_sparse_for_T(const double &sparse_threshold) { ModuleBase::TITLE("LCAO_Hamilt","cal_STN_R_sparse_for_T"); @@ -397,8 +421,8 @@ void LCAO_Hamilt::cal_STN_R_sparse_for_T(const double &sparse_threshold) ModuleBase::Vector3 dtau, tau1, tau2; ModuleBase::Vector3 dtau1, dtau2, tau0; - double temp_value_double; - std::complex temp_value_complex; + double tmp=0.0; + std::complex tmpc=complex(0.0,0.0); for(int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1) { @@ -473,18 +497,18 @@ void LCAO_Hamilt::cal_STN_R_sparse_for_T(const double &sparse_threshold) if(GlobalV::NSPIN!=4) { - temp_value_double = this->LM->Hloc_fixedR[index]; - if (std::abs(temp_value_double) > sparse_threshold) + tmp = this->LM->Hloc_fixedR[index]; + if (std::abs(tmp) > sparse_threshold) { - this->LM->TR_sparse[dR][iw1_all][iw2_all] = temp_value_double; + this->LM->TR_sparse[dR][iw1_all][iw2_all] = tmp; } } else { - temp_value_complex = this->LM->Hloc_fixedR_soc[index]; - if(std::abs(temp_value_complex) > sparse_threshold) + tmpc = this->LM->Hloc_fixedR_soc[index]; + if(std::abs(tmpc) > sparse_threshold) { - this->LM->TR_soc_sparse[dR][iw1_all][iw2_all] = temp_value_complex; + this->LM->TR_soc_sparse[dR][iw1_all][iw2_all] = tmpc; } } @@ -502,7 +526,7 @@ void LCAO_Hamilt::cal_STN_R_sparse_for_T(const double &sparse_threshold) void LCAO_Hamilt::cal_SR_sparse(const double &sparse_threshold, hamilt::Hamilt>* p_ham) { ModuleBase::TITLE("LCAO_Hamilt","cal_SR_sparse"); - set_R_range_sparse(); + set_R_range_sparse(*this->LM); //cal_STN_R_sparse(current_spin, sparse_threshold); if(GlobalV::NSPIN!=4) { @@ -526,7 +550,7 @@ void LCAO_Hamilt::cal_TR_sparse(const double &sparse_threshold) this->LM->Hloc_fixedR.resize(this->LM->ParaV->nnr); this->LM->zeros_HSR('T'); this->genH.build_ST_new('T', 0, GlobalC::ucell, this->LM->Hloc_fixedR.data()); - this->set_R_range_sparse(); + this->set_R_range_sparse(*this->LM); this->cal_STN_R_sparse_for_T(sparse_threshold); return; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h index 1422752c15d..b57ebf43e7b 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h @@ -4,10 +4,13 @@ #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "LCAO_gen_fixedH.h" -#include "module_hamilt_lcao/module_gint/gint_gamma.h" -#include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" #include "module_hamilt_general/hamilt.h" +#include "module_hamilt_lcao/module_gint/grid_technique.h" + +#include "module_hamilt_lcao/module_gint/gint_gamma.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" + #ifdef __EXX #include @@ -21,10 +24,15 @@ class LCAO_Hamilt ~LCAO_Hamilt(); - void grid_prepare(const Grid_Technique& gt, const ModulePW::PW_Basis& rhopw, const ModulePW::PW_Basis_Big& bigpw); + void grid_prepare( + const Grid_Technique& gt, + Gint_Gamma &gint_gamma, + Gint_k &gint_k, + const ModulePW::PW_Basis& rhopw, + const ModulePW::PW_Basis_Big& bigpw); // jingan add 2021-6-4 - void set_R_range_sparse(); + void set_R_range_sparse(LCAO_Matrix &lm); void cal_HContainer_sparse_d(const int ¤t_spin, const double &sparse_threshold, @@ -38,7 +46,14 @@ class LCAO_Hamilt void cal_dSTN_R_sparse(const int ¤t_spin, const double &sparse_threshold); - void cal_STN_R_sparse_for_S(const double &sparse_threshold); + // mohan update 2024-04-01 + void cal_STN_R_sparse_for_S( + const Parallel_Orbitals pv, + std::vector &slocR, + std::vector> &slocR_soc, + std::map, std::map>> &SR_sparse, + std::map, std::map>>> &SR_soc_sparse, + const double &sparse_threshold); void cal_STN_R_sparse_for_T(const double &sparse_threshold); @@ -70,16 +85,13 @@ class LCAO_Hamilt void destroy_TR_sparse(void); - void cal_dH_sparse(const int ¤t_spin, const double &sparse_threshold); + void cal_dH_sparse( + const int ¤t_spin, + const double &sparse_threshold, + Gint_k &gint_k); // mohan add 2024-04-01 void destroy_dH_R_sparse(void); - // used for gamma only algorithms. - Gint_Gamma GG; - - // used for k-dependent grid integration. - Gint_k GK; - // use overlap matrix to generate fixed Hamiltonian LCAO_gen_fixedH genH; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.cpp index b394ccd032f..38bf7b138f3 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.cpp @@ -45,11 +45,13 @@ void LCAO_Matrix::divide_HS_in_frag(const bool isGamma, Parallel_Orbitals &pv, c { na[it] = GlobalC::ucell.atoms[it].na; } + GlobalC::ld.init(GlobalC::ORB, GlobalC::ucell.nat, GlobalC::ucell.ntype, pv, na); + if(GlobalV::deepks_scf) { if(isGamma) @@ -72,7 +74,11 @@ void LCAO_Matrix::allocate_HS_gamma(const long &nloc) ModuleBase::TITLE("LCAO_Matrix","allocate_HS_gamma"); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nloc",nloc); - if(nloc==0) return; //mohan fix bug 2012-05-25 + + if(nloc==0) + { + return; + } // because we initilize in the constructor function // with dimension '1', so here we reconstruct these @@ -95,7 +101,11 @@ void LCAO_Matrix::allocate_HS_k(const long &nloc) ModuleBase::TITLE("LCAO_Matrix","allocate_HS_k"); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nloc",nloc); - if(nloc==0) return; //mohan fix bug 2012-05-25 + + if(nloc==0) + { + return; //mohan fix bug 2012-05-25 + } // because we initilize in the constructor function // with dimension '1', so here we reconstruct these @@ -116,19 +126,33 @@ void LCAO_Matrix::allocate_HS_R(const int &nnR) if(GlobalV::NSPIN!=4) { this->SlocR.resize(nnR); - if(GlobalV::CALCULATION!="get_S") this->Hloc_fixedR.resize(nnR); + if(GlobalV::CALCULATION!="get_S") + { + this->Hloc_fixedR.resize(nnR); + } ModuleBase::GlobalFunc::ZEROS(SlocR.data(), nnR); - if(GlobalV::CALCULATION!="get_S") ModuleBase::GlobalFunc::ZEROS(Hloc_fixedR.data(), nnR); + + if(GlobalV::CALCULATION!="get_S") + { + ModuleBase::GlobalFunc::ZEROS(Hloc_fixedR.data(), nnR); + } } else { this->SlocR_soc.resize(nnR); - if(GlobalV::CALCULATION!="get_S") this->Hloc_fixedR_soc.resize(nnR); + if(GlobalV::CALCULATION!="get_S") + { + this->Hloc_fixedR_soc.resize(nnR); + } ModuleBase::GlobalFunc::ZEROS(SlocR_soc.data(), nnR); - if(GlobalV::CALCULATION!="get_S") ModuleBase::GlobalFunc::ZEROS(Hloc_fixedR_soc.data(), nnR); - + + if(GlobalV::CALCULATION!="get_S") + { + ModuleBase::GlobalFunc::ZEROS(Hloc_fixedR_soc.data(), nnR); + } + } return; @@ -139,16 +163,26 @@ void LCAO_Matrix::set_HSgamma(const int& iw1_all, const int& iw2_all, const doub LCAO_Matrix::set_mat2d(iw1_all, iw2_all, v, *this->ParaV, HSloc); return; } + + void LCAO_Matrix::set_HSk(const int& iw1_all, const int& iw2_all, const std::complex& v, const char& dtype, const int spin) { - if (dtype=='S')//overlap Hamiltonian. - LCAO_Matrix::set_mat2d>(iw1_all, iw2_all, v, *this->ParaV, this->Sloc2.data()); - else if (dtype=='T' || dtype=='N')// kinetic and nonlocal Hamiltonian. - LCAO_Matrix::set_mat2d>(iw1_all, iw2_all, v, *this->ParaV, this->Hloc_fixed2.data()); - else if (dtype=='L') // Local potential Hamiltonian. - LCAO_Matrix::set_mat2d>(iw1_all, iw2_all, v, *this->ParaV, this->Hloc2.data()); - else - ModuleBase::WARNING_QUIT("LCAO_Matrix", "set_HSk"); + if (dtype=='S')//overlap Hamiltonian. + { + LCAO_Matrix::set_mat2d>(iw1_all, iw2_all, v, *this->ParaV, this->Sloc2.data()); + } + else if (dtype=='T' || dtype=='N')// kinetic and nonlocal Hamiltonian. + { + LCAO_Matrix::set_mat2d>(iw1_all, iw2_all, v, *this->ParaV, this->Hloc_fixed2.data()); + } + else if (dtype=='L') // Local potential Hamiltonian. + { + LCAO_Matrix::set_mat2d>(iw1_all, iw2_all, v, *this->ParaV, this->Hloc2.data()); + } + else + { + ModuleBase::WARNING_QUIT("LCAO_Matrix", "set_HSk"); + } return; } @@ -267,17 +301,38 @@ void LCAO_Matrix::zeros_HSgamma(const char &mtype) long long beg, len; if (mtype=='S') { - ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, (long long)this->Sloc.size(), (long long)512, beg, len); - ModuleBase::GlobalFunc::ZEROS(this->Sloc.data() + beg, len); + ModuleBase::BLOCK_TASK_DIST_1D( + num_threads, + thread_id, + (long long)this->Sloc.size(), + (long long)512, + beg, + len); + + ModuleBase::GlobalFunc::ZEROS(this->Sloc.data() + beg, len); } else if (mtype=='T') { - ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, (long long)this->Hloc_fixed.size(), (long long)512, beg, len); + ModuleBase::BLOCK_TASK_DIST_1D( + num_threads, + thread_id, + (long long)this->Hloc_fixed.size(), + (long long)512, + beg, + len); + ModuleBase::GlobalFunc::ZEROS(this->Hloc_fixed.data() + beg, len); } else if (mtype=='H') { - ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, (long long)this->Hloc.size(), (long long)512, beg, len); + ModuleBase::BLOCK_TASK_DIST_1D( + num_threads, + thread_id, + (long long)this->Hloc.size(), + (long long)512, + beg, + len); + ModuleBase::GlobalFunc::ZEROS(this->Hloc.data() + beg, len); } }; @@ -404,12 +459,10 @@ void LCAO_Matrix::print_HSk(const char &mtype, const char &vtype, const double & if( std::abs(v) > accuracy ) { - // os << std::setw(15) << v; os << v << "\t"; } else { - // os << std::setw(15) << "0"; os << "0" << "\t"; } } @@ -556,9 +609,6 @@ void LCAO_Matrix::allocate_Hloc_fixedR_tr(void) { ModuleBase::TITLE("LCAO_Matrix","allocate_Hloc_fixedR_tr"); - //int R_x = 10; - //int R_y = 10; - //int R_z = 10; int R_x = GlobalC::GridD.getCellX(); int R_y = GlobalC::GridD.getCellY(); int R_z = GlobalC::GridD.getCellZ(); @@ -566,26 +616,16 @@ void LCAO_Matrix::allocate_Hloc_fixedR_tr(void) if(GlobalV::NSPIN!=4) { Hloc_fixedR_tr = new double***[R_x]; - //HR_tr = new double***[R_x]; - //SlocR_tr = new double***[R_x]; for(int ix=0; ixParaV->nloc]; - //HR_tr[ix][iy][iz] = new double[this->ParaV->nloc]; - //SlocR_tr[ix][iy][iz] = new double[this->ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(Hloc_fixedR_tr[ix][iy][iz], this->ParaV->nloc); - //ModuleBase::GlobalFunc::ZEROS(HR_tr[ix][iy][iz], this->ParaV->nloc); - //ModuleBase::GlobalFunc::ZEROS(SlocR_tr[ix][iy][iz], this->ParaV->nloc); } } } @@ -593,36 +633,20 @@ void LCAO_Matrix::allocate_Hloc_fixedR_tr(void) else { Hloc_fixedR_tr_soc = new std::complex***[R_x]; - //HR_tr = new double***[R_x]; - //SlocR_tr = new double***[R_x]; for(int ix=0; ix**[R_y]; - //HR_tr[ix] = new double**[R_y]; - //SlocR_tr[ix] = new double**[R_y]; for(int iy=0; iy*[R_z]; - //HR_tr[ix][iy] = new double*[R_z]; - //SlocR_tr[ix][iy] = new double*[R_z]; for(int iz=0; iz[this->ParaV->nloc]; - //HR_tr[ix][iy][iz] = new double[this->ParaV->nloc]; - //SlocR_tr[ix][iy][iz] = new double[this->ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(Hloc_fixedR_tr_soc[ix][iy][iz], this->ParaV->nloc); - //ModuleBase::GlobalFunc::ZEROS(HR_tr[ix][iy][iz], this->ParaV->nloc); - //ModuleBase::GlobalFunc::ZEROS(SlocR_tr[ix][iy][iz], this->ParaV->nloc); } } } } -//std::cout<<"R_x: "<ParaV->nloc: "<ParaV->nloc< &v) +void LCAO_Matrix::set_HR_tr_soc( + const int &Rx, + const int &Ry, + const int &Rz, + const int &iw1_all, + const int &iw2_all, + const std::complex &v) { const int ir = this->ParaV->global2local_row(iw1_all); const int ic = this->ParaV->global2local_col(iw2_all); -//std::cout<<"ir: "<ParaV->nrow+ir; -//std::cout<<"index: "<ParaV->ncol+ic; -//std::cout<<"index: "<ParaV->nloc: "<ParaV->nloc<ParaV->nloc); -//std::cout<<"Rx: "< +#include + +#include "module_base/memory.h" +#include "module_base/parallel_reduce.h" +#include "module_base/timer.h" +#include "module_base/tool_threading.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_elecstate/cal_dm.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" +#include "module_elecstate/elecstate_lcao.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_io/write_HS.h" + +#ifdef __DEEPKS +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#endif + +#ifdef _OPENMP +#include +#endif + +#include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h" + + +void Force_LCAO_k::cal_foverlap_k(const bool isforce, + const bool isstress, + Record_adj &ra, + const psi::Psi> *psi, + Local_Orbital_Charge &loc, + Parallel_Orbitals &pv, + LCAO_Matrix &lm, + const elecstate::DensityMatrix, double> *DM, + ModuleBase::matrix &foverlap, + ModuleBase::matrix &soverlap, + const elecstate::ElecState *pelec, + const int &nks, + const K_Vectors &kv) +{ + ModuleBase::TITLE("Force_LCAO_k", "cal_foverlap_k"); + ModuleBase::timer::tick("Force_LCAO_k", "cal_foverlap_k"); + + // construct a DensityMatrix object + elecstate::DensityMatrix, double> EDM(&kv,&pv,GlobalV::NSPIN); + + //-------------------------------------------- + // calculate the energy density matrix here. + //-------------------------------------------- + ModuleBase::timer::tick("Force_LCAO_k", "cal_edm_2d"); + + ModuleBase::matrix wgEkb; + wgEkb.create(nks, GlobalV::NBANDS); + ModuleBase::Memory::record("Force::wgEkb", sizeof(double) * nks * GlobalV::NBANDS); +#ifdef _OPENMP +#pragma omp parallel for collapse(2) schedule(static, 1024) +#endif + for (int ik = 0; ik < nks; ik++) + { + for (int ib = 0; ib < GlobalV::NBANDS; ib++) + { + wgEkb(ik, ib) = pelec->wg(ik, ib) * pelec->ekb(ik, ib); + } + } + std::vector edm_k(nks); + + // use the original formula (Hamiltonian matrix) to calculate energy density matrix + if (DM->EDMK.size()) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for (int ik = 0; ik < nks; ++ik) + { + //edm_k[ik] = loc.edm_k_tddft[ik]; + EDM.set_DMK_pointer(ik,DM->EDMK[ik].c); + } + } + else + { + //elecstate::cal_dm(pv, wgEkb, psi[0], edm_k); + // cal_dm_psi + elecstate::cal_dm_psi(EDM.get_paraV_pointer(), wgEkb, psi[0], EDM); + } + + //loc.cal_dm_R(edm_k, ra, edm2d, kv); + + // cal_dm_2d + EDM.init_DMR(ra,&GlobalC::ucell); + EDM.cal_DMR(); + EDM.sum_DMR_spin(); + // + ModuleBase::timer::tick("Force_LCAO_k", "cal_edm_2d"); + //-------------------------------------------- + // summation \sum_{i,j} E(i,j)*dS(i,j) + // BEGIN CALCULATION OF FORCE OF EACH ATOM + //-------------------------------------------- + int total_irr = 0; +#ifdef _OPENMP +#pragma omp parallel + { + int num_threads = omp_get_num_threads(); + ModuleBase::matrix local_soverlap(3, 3); + int local_total_irr = 0; +#else + ModuleBase::matrix& local_soverlap = soverlap; + int& local_total_irr = total_irr; +#endif + + ModuleBase::Vector3 tau1, dtau, tau2; + +#ifdef _OPENMP +#pragma omp for schedule(dynamic) +#endif + for (int iat = 0; iat < GlobalC::ucell.nat; iat++) + { + const int T1 = GlobalC::ucell.iat2it[iat]; + Atom* atom1 = &GlobalC::ucell.atoms[T1]; + const int I1 = GlobalC::ucell.iat2ia[iat]; + // get iat1 + int iat1 = GlobalC::ucell.itia2iat(T1, I1); + double* foverlap_iat; + if (isforce) + { + foverlap_iat = &foverlap(iat, 0); + } + +#ifdef _OPENMP + // using local stack to avoid false sharing in multi-threaded case + double foverlap_temp[3] = {0.0, 0.0, 0.0}; + if (num_threads > 1) + { + foverlap_iat = foverlap_temp; + } +#endif + int irr = pv.nlocstart[iat]; + const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0); + for (int cb = 0; cb < ra.na_each[iat]; ++cb) + { + const int T2 = ra.info[iat][cb][3]; + const int I2 = ra.info[iat][cb][4]; + + const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0); + + Atom* atom2 = &GlobalC::ucell.atoms[T2]; + + // get iat2 + int iat2 = GlobalC::ucell.itia2iat(T2, I2); + double Rx = ra.info[iat][cb][0]; + double Ry = ra.info[iat][cb][1]; + double Rz = ra.info[iat][cb][2]; + // get BaseMatrix + hamilt::BaseMatrix* tmp_matrix = EDM.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz); + if(tmp_matrix == nullptr) + { + continue; + } + int row_ap = pv.atom_begin_row[iat1]; + int col_ap = pv.atom_begin_col[iat2]; + // get DMR + for (int mu = 0; mu < pv.get_row_size(iat1); ++mu) + { + for (int nu = 0; nu < pv.get_col_size(iat2); ++nu) + { + // here do not sum over spin due to EDM.sum_DMR_spin(); + double edm2d1 = tmp_matrix->get_value(mu,nu); + double edm2d2 = 2.0 * edm2d1; + + if (isforce) + { + foverlap_iat[0] -= edm2d2 * lm.DSloc_Rx[irr]; + foverlap_iat[1] -= edm2d2 * lm.DSloc_Ry[irr]; + foverlap_iat[2] -= edm2d2 * lm.DSloc_Rz[irr]; + } + if (isstress) + { + for (int ipol = 0; ipol < 3; ipol++) + { + local_soverlap(0, ipol) += edm2d1 * lm.DSloc_Rx[irr] + * lm.DH_r[irr * 3 + ipol]; + if (ipol < 1) + { + continue; + } + local_soverlap(1, ipol) += edm2d1 * lm.DSloc_Ry[irr] + * lm.DH_r[irr * 3 + ipol]; + if (ipol < 2) + { + continue; + } + local_soverlap(2, ipol) += edm2d1 * lm.DSloc_Rz[irr] + * lm.DH_r[irr * 3 + ipol]; + } + } + //} + ++local_total_irr; + ++irr; + } // end kk + } // end jj + } // end cb +#ifdef _OPENMP + if (isforce && num_threads > 1) + { + foverlap(iat, 0) += foverlap_iat[0]; + foverlap(iat, 1) += foverlap_iat[1]; + foverlap(iat, 2) += foverlap_iat[2]; + } +#endif + } // end iat +#ifdef _OPENMP +#pragma omp critical(cal_foverlap_k_reduce) + { + total_irr += local_total_irr; + if (isstress) + { + for (int ipol = 0; ipol < 3; ipol++) + { + soverlap(0, ipol) += local_soverlap(0, ipol); + if (ipol < 1) + { + continue; + } + soverlap(1, ipol) += local_soverlap(1, ipol); + if (ipol < 2) + { + continue; + } + soverlap(2, ipol) += local_soverlap(2, ipol); + } + } + } + } +#endif + + if (isstress) + { + StressTools::stress_fill(GlobalC::ucell.lat0, GlobalC::ucell.omega, soverlap); + } + + if (total_irr != pv.nnr) + { + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wrong irr", total_irr); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wrong LNNR.nnr", pv.nnr); + ModuleBase::WARNING_QUIT("Force_LCAO_k::cal_foverlap_k", "irr!=LNNR.nnr"); + } + + ModuleBase::timer::tick("Force_LCAO_k", "cal_foverlap_k"); + return; +} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_k.cpp new file mode 100644 index 00000000000..ab95392b117 --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_k.cpp @@ -0,0 +1,165 @@ +#include "FORCE_k.h" + +#include +#include + +#include "module_base/memory.h" +#include "module_base/parallel_reduce.h" +#include "module_base/timer.h" +#include "module_base/tool_threading.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_elecstate/cal_dm.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" +#include "module_elecstate/elecstate_lcao.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_io/write_HS.h" + +#ifdef __DEEPKS +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#endif + +#ifdef _OPENMP +#include +#endif + + +void Force_LCAO_k::cal_ftvnl_dphi_k(const elecstate::DensityMatrix, double>* DM, + const Parallel_Orbitals &pv, + LCAO_Matrix &lm, + const bool isforce, + const bool isstress, + Record_adj& ra, + ModuleBase::matrix& ftvnl_dphi, + ModuleBase::matrix& stvnl_dphi) +{ + ModuleBase::TITLE("Force_LCAO_k", "cal_ftvnl_dphi_k"); + ModuleBase::timer::tick("Force_LCAO_k", "cal_ftvnl_dphi_k"); + + int total_irr = 0; + // get the adjacent atom's information. + + // GlobalV::ofs_running << " calculate the ftvnl_dphi_k force" << std::endl; +#ifdef _OPENMP +#pragma omp parallel + { + int num_threads = omp_get_num_threads(); + ModuleBase::matrix local_stvnl_dphi(3, 3); + int local_total_irr = 0; +#pragma omp for schedule(dynamic) +#else + ModuleBase::matrix& local_stvnl_dphi = stvnl_dphi; + int& local_total_irr = total_irr; +#endif + for (int iat = 0; iat < GlobalC::ucell.nat; iat++) + { + const int T1 = GlobalC::ucell.iat2it[iat]; + Atom* atom1 = &GlobalC::ucell.atoms[T1]; + const int I1 = GlobalC::ucell.iat2ia[iat]; + // get iat1 + int iat1 = GlobalC::ucell.itia2iat(T1, I1); + // + int irr = pv.nlocstart[iat]; + const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0); + double* ftvnl_dphi_iat; + if (isforce) + { + ftvnl_dphi_iat = &ftvnl_dphi(iat, 0); + } +#ifdef _OPENMP + // using local stack to avoid false sharing in multi-threaded case + double ftvnl_dphi_temp[3] = {0.0, 0.0, 0.0}; + if (num_threads > 1) + { + ftvnl_dphi_iat = ftvnl_dphi_temp; + } +#endif + for (int cb = 0; cb < ra.na_each[iat]; ++cb) + { + const int T2 = ra.info[iat][cb][3]; + const int I2 = ra.info[iat][cb][4]; + const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0); + Atom* atom2 = &GlobalC::ucell.atoms[T2]; + // get iat2 + int iat2 = GlobalC::ucell.itia2iat(T2, I2); + double Rx = ra.info[iat][cb][0]; + double Ry = ra.info[iat][cb][1]; + double Rz = ra.info[iat][cb][2]; + // get BaseMatrix + if (pv.get_row_size(iat1) <= 0 || pv.get_col_size(iat2) <= 0) + { + continue; + } + std::vector*> tmp_matrix; + for (int is = 0; is < GlobalV::NSPIN; ++is) + { + tmp_matrix.push_back(DM->get_DMR_pointer(is+1)->find_matrix(iat1, iat2, Rx, Ry, Rz)); + } + //hamilt::BaseMatrix* tmp_matrix = DM->get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz); + for (int mu = 0; mu < pv.get_row_size(iat1); ++mu) + { + for (int nu = 0; nu < pv.get_col_size(iat2); ++nu) + { + // get value from DM + double dm2d1 = 0.0; + for (int is = 0; is < GlobalV::NSPIN; ++is) + { + dm2d1 += tmp_matrix[is]->get_value(mu, nu); + } + double dm2d2 = 2.0 * dm2d1; + // + if (isforce) + { + ftvnl_dphi_iat[0] += dm2d2 * lm.DHloc_fixedR_x[irr]; + ftvnl_dphi_iat[1] += dm2d2 * lm.DHloc_fixedR_y[irr]; + ftvnl_dphi_iat[2] += dm2d2 * lm.DHloc_fixedR_z[irr]; + } + if (isstress) + { + local_stvnl_dphi(0, 0) -= dm2d1 * lm.stvnl11[irr]; + local_stvnl_dphi(0, 1) -= dm2d1 * lm.stvnl12[irr]; + local_stvnl_dphi(0, 2) -= dm2d1 * lm.stvnl13[irr]; + local_stvnl_dphi(1, 1) -= dm2d1 * lm.stvnl22[irr]; + local_stvnl_dphi(1, 2) -= dm2d1 * lm.stvnl23[irr]; + local_stvnl_dphi(2, 2) -= dm2d1 * lm.stvnl33[irr]; + } + //} + ++local_total_irr; + ++irr; + } // end kk + } // end jj + } // end cb +#ifdef _OPENMP + if (isforce && num_threads > 1) + { + ftvnl_dphi(iat, 0) += ftvnl_dphi_iat[0]; + ftvnl_dphi(iat, 1) += ftvnl_dphi_iat[1]; + ftvnl_dphi(iat, 2) += ftvnl_dphi_iat[2]; + } +#endif + } // end iat +#ifdef _OPENMP +#pragma omp critical(cal_ftvnl_dphi_k_reduce) + { + total_irr += local_total_irr; + if (isstress) + { + stvnl_dphi(0, 0) += local_stvnl_dphi(0, 0); + stvnl_dphi(0, 1) += local_stvnl_dphi(0, 1); + stvnl_dphi(0, 2) += local_stvnl_dphi(0, 2); + stvnl_dphi(1, 1) += local_stvnl_dphi(1, 1); + stvnl_dphi(1, 2) += local_stvnl_dphi(1, 2); + stvnl_dphi(2, 2) += local_stvnl_dphi(2, 2); + } + } + } +#endif + assert(total_irr == pv.nnr); + + if (isstress) + { + StressTools::stress_fill(GlobalC::ucell.lat0, GlobalC::ucell.omega, stvnl_dphi); + } + + ModuleBase::timer::tick("Force_LCAO_k", "cal_ftvnl_dphi_k"); + return; +} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/fvl_dphi_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/fvl_dphi_k.cpp new file mode 100644 index 00000000000..f6d267c511b --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/fvl_dphi_k.cpp @@ -0,0 +1,96 @@ +#include "FORCE_k.h" + +#include +#include + +#include "module_base/memory.h" +#include "module_base/parallel_reduce.h" +#include "module_base/timer.h" +#include "module_base/tool_threading.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_elecstate/cal_dm.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" +#include "module_elecstate/elecstate_lcao.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_io/write_HS.h" + +#ifdef __DEEPKS +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#endif + +#ifdef _OPENMP +#include +#endif + + +// calculate the force due to < phi | Vlocal | dphi > +void Force_LCAO_k::cal_fvl_dphi_k(const bool isforce, + const bool isstress, + LCAO_Matrix &lm, + Gint_k &gint_k, + const elecstate::Potential* pot_in, + ModuleBase::matrix& fvl_dphi, + ModuleBase::matrix& svl_dphi, + double** DM_R) +{ + ModuleBase::TITLE("Force_LCAO_k", "cal_fvl_dphi_k"); + ModuleBase::timer::tick("Force_LCAO_k", "cal_fvl_dphi_k"); + + if (!isforce && !isstress) + { + ModuleBase::timer::tick("Force_LCAO_k", "cal_fvl_dphi_k"); + return; + } + + assert(lm.DHloc_fixedR_x != NULL); + assert(lm.DHloc_fixedR_y != NULL); + assert(lm.DHloc_fixedR_z != NULL); + + int istep = 1; + + for (int is = 0; is < GlobalV::NSPIN; ++is) + { + GlobalV::CURRENT_SPIN = is; + + const double* vr_eff1 = pot_in->get_effective_v(GlobalV::CURRENT_SPIN); + const double* vofk_eff1 = nullptr; + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + { + vofk_eff1 = pot_in->get_effective_vofk(GlobalV::CURRENT_SPIN); + } + + //-------------------------------- + // Grid integration here. + //-------------------------------- + // fvl_dphi can not be set to zero here if Vna is used + if (isstress || isforce) + { + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + { + Gint_inout inout(DM_R, + is, + vr_eff1, + vofk_eff1, + isforce, + isstress, + &fvl_dphi, + &svl_dphi, + Gint_Tools::job_type::force_meta); + gint_k.cal_gint(&inout); + } + else + { + Gint_inout inout(DM_R, is, vr_eff1, isforce, isstress, &fvl_dphi, &svl_dphi, Gint_Tools::job_type::force); + gint_k.cal_gint(&inout); + } + } + } + + if (isstress) + { + StressTools::stress_fill(-1.0, GlobalC::ucell.omega, svl_dphi); + } + + ModuleBase::timer::tick("Force_LCAO_k", "cal_fvl_dphi_k"); + return; +} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp index 6cfac40f2d7..65e3ebb290e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp @@ -1,3 +1,29 @@ +#include "FORCE_k.h" + +#include +#include + +#include "module_base/memory.h" +#include "module_base/parallel_reduce.h" +#include "module_base/timer.h" +#include "module_base/tool_threading.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_elecstate/cal_dm.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" +#include "module_elecstate/elecstate_lcao.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_io/write_HS.h" + +#ifdef __DEEPKS +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#endif + +#ifdef _OPENMP +#include +#endif + + +typedef std::tuple key_tuple; // must consider three-center H matrix. void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix, double>* DM, diff --git a/source/module_io/output_mat_sparse.cpp b/source/module_io/output_mat_sparse.cpp index 693528a43ee..0c82f71a350 100644 --- a/source/module_io/output_mat_sparse.cpp +++ b/source/module_io/output_mat_sparse.cpp @@ -7,7 +7,8 @@ namespace ModuleIO { template - Output_Mat_Sparse::Output_Mat_Sparse(int out_mat_hsR, + Output_Mat_Sparse::Output_Mat_Sparse( + int out_mat_hsR, int out_mat_dh, int out_mat_t, int out_mat_r, @@ -15,6 +16,7 @@ namespace ModuleIO const ModuleBase::matrix& v_eff, const Parallel_Orbitals& pv, LCAO_Hamilt& UHM, + Gint_k& gint_k, // mohan add 2024-04-01 LCAO_Matrix& LM, const K_Vectors& kv, hamilt::Hamilt* p_ham) @@ -26,15 +28,21 @@ namespace ModuleIO _v_eff(v_eff), _pv(pv), _UHM(UHM), + _gint_k(gint_k), // mohan add 2024-04-01 _LM(LM), _kv(kv), _p_ham(p_ham) { - } - template<> - void Output_Mat_Sparse::write() {} - template<> - void Output_Mat_Sparse>::write() +} + +template<> +void Output_Mat_Sparse::write(void) +{ +} + + +template<> +void Output_Mat_Sparse>::write(void) { if (_out_mat_hsR) { @@ -48,8 +56,14 @@ namespace ModuleIO if (_out_mat_dh) { - output_dH_R(_istep, this->_v_eff, this->_UHM, _kv); // LiuXh add 2019-07-15 - } + output_dH_R( + _istep, + this->_v_eff, + this->_UHM, + this->_gint_k, // mohan add 2024-04-01 + this->_LM, + _kv); // LiuXh add 2019-07-15 + } // add by jingan for out r_R matrix 2019.8.14 if (_out_mat_r) @@ -67,6 +81,7 @@ namespace ModuleIO } } - template class Output_Mat_Sparse; - template class Output_Mat_Sparse>; -} // namespace ModuleIO \ No newline at end of file +template class Output_Mat_Sparse; +template class Output_Mat_Sparse>; + +} // namespace ModuleIO diff --git a/source/module_io/output_mat_sparse.h b/source/module_io/output_mat_sparse.h index 4d5bcde5c26..8012b108bca 100644 --- a/source/module_io/output_mat_sparse.h +++ b/source/module_io/output_mat_sparse.h @@ -4,6 +4,7 @@ #include "module_basis/module_ao/parallel_orbitals.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hsolver/hsolver_lcao.h" #include "output_interface.h" @@ -22,9 +23,11 @@ namespace ModuleIO const ModuleBase::matrix& v_eff, const Parallel_Orbitals& pv, LCAO_Hamilt& UHM, + Gint_k &gint_k, // mohan add 2024-04-01 LCAO_Matrix& LM, const K_Vectors& kv, hamilt::Hamilt* p_ham); + void write() override; private: @@ -36,10 +39,11 @@ namespace ModuleIO const ModuleBase::matrix& _v_eff; const Parallel_Orbitals& _pv; LCAO_Hamilt& _UHM; + Gint_k& _gint_k; // mohan add 2024-04-01 LCAO_Matrix& _LM; const K_Vectors& _kv; hamilt::Hamilt* _p_ham; }; } // namespace ModuleIO -#endif // OUTPUT_MAT_SPARSE_H \ No newline at end of file +#endif // OUTPUT_MAT_SPARSE_H diff --git a/source/module_io/write_HS_R.cpp b/source/module_io/write_HS_R.cpp index 81346621f78..76bf11c7cae 100644 --- a/source/module_io/write_HS_R.cpp +++ b/source/module_io/write_HS_R.cpp @@ -59,6 +59,8 @@ void ModuleIO::output_HS_R(const int& istep, void ModuleIO::output_dH_R(const int& istep, const ModuleBase::matrix& v_eff, LCAO_Hamilt& UHM, + Gint_k& gint_k, // mohan add 2024-04-01 + LCAO_Matrix &lm, // mohan add 2024-04-01 const K_Vectors& kv, const bool& binary, const double& sparse_threshold) @@ -66,11 +68,13 @@ void ModuleIO::output_dH_R(const int& istep, ModuleBase::TITLE("ModuleIO","output_dH_R"); ModuleBase::timer::tick("ModuleIO","output_dH_R"); - UHM.LM->Hloc_fixedR.resize(UHM.LM->ParaV->nnr); - UHM.GK.allocate_pvdpR(); + lm.Hloc_fixedR.resize(lm.ParaV->nnr); + gint_k.allocate_pvdpR(); if(GlobalV::NSPIN==1||GlobalV::NSPIN==4) { - UHM.cal_dH_sparse(0, sparse_threshold); + // mohan add 2024-04-01 + assert(GlobalV::CURRENT_SPIN==0); + UHM.cal_dH_sparse(GlobalV::CURRENT_SPIN, sparse_threshold, gint_k); } else if(GlobalV::NSPIN==2) { @@ -91,19 +95,20 @@ void ModuleIO::output_dH_R(const int& istep, if(GlobalV::VL_IN_H) { Gint_inout inout(vr_eff1, GlobalV::CURRENT_SPIN, Gint_Tools::job_type::dvlocal); - UHM.GK.cal_gint(&inout); + gint_k.cal_gint(&inout); } } - UHM.cal_dH_sparse(GlobalV::CURRENT_SPIN, sparse_threshold); + UHM.cal_dH_sparse(GlobalV::CURRENT_SPIN, sparse_threshold, gint_k); } } } - ModuleIO::save_dH_sparse(istep, *UHM.LM, sparse_threshold, binary); + // mohan update 2024-04-01 + ModuleIO::save_dH_sparse(istep, lm, sparse_threshold, binary); UHM.destroy_dH_R_sparse(); - UHM.GK.destroy_pvdpR(); + gint_k.destroy_pvdpR(); ModuleBase::timer::tick("ModuleIO","output_HS_R"); return; diff --git a/source/module_io/write_HS_R.h b/source/module_io/write_HS_R.h index 5389af2f912..738af1d3a7d 100644 --- a/source/module_io/write_HS_R.h +++ b/source/module_io/write_HS_R.h @@ -23,7 +23,9 @@ namespace ModuleIO void output_dH_R( const int &istep, const ModuleBase::matrix& v_eff, - LCAO_Hamilt &UHM, + LCAO_Hamilt &UHM, + Gint_k& gint_k, // mohan add 2024-04-01 + LCAO_Matrix &lm, // mohan add 2024-04-01 const K_Vectors& kv, const bool& binary = false, const double& sparse_threshold = 1e-10); diff --git a/source/module_io/write_Vxc.hpp b/source/module_io/write_Vxc.hpp index baa24ceef3e..b9f5c435b96 100644 --- a/source/module_io/write_Vxc.hpp +++ b/source/module_io/write_Vxc.hpp @@ -3,25 +3,45 @@ #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h" #include "module_base/scalapack_connector.h" #include "module_base/parallel_reduce.h" + template struct TGint; + template <> struct TGint { using type = Gint_Gamma; }; + template <> struct TGint> { using type = Gint_k; }; + + namespace ModuleIO { - inline void gint_vl(Gint_Gamma& gg, Gint_inout& io, LCAO_Matrix& lm) { gg.cal_vlocal(&io, &lm, false); }; - inline void gint_vl(Gint_k& gk, Gint_inout& io, LCAO_Matrix& lm, ModuleBase::matrix& wg) { gk.cal_gint(&io); }; + + inline void gint_vl( + Gint_Gamma& gg, + Gint_inout& io, + LCAO_Matrix& lm) + { + gg.cal_vlocal(&io, &lm, false); + }; + + inline void gint_vl( + Gint_k& gk, + Gint_inout& io, + LCAO_Matrix& lm, + ModuleBase::matrix& wg) + { + gk.cal_gint(&io); + }; void set_para2d_MO(const Parallel_Orbitals& pv, const int nbands, Parallel_2D& p2d) { std::ofstream ofs; #ifdef __MPI - int dsize; + int dsize=0; MPI_Comm_size(MPI_COMM_WORLD, &dsize); p2d.set_block_size(pv.nb); p2d.set_proc_dim(dsize); @@ -37,8 +57,14 @@ namespace ModuleIO #endif } - std::vector> cVc(std::complex* V, std::complex* c, int nbasis, int nbands, const Parallel_Orbitals& pv, const Parallel_2D& p2d) - { + std::vector> cVc( + std::complex* V, + std::complex* c, + int nbasis, + int nbands, + const Parallel_Orbitals& pv, + const Parallel_2D& p2d) + { std::vector> Vc(pv.nloc_wfc, 0.0); char transa = 'N'; char transb = 'N'; @@ -46,21 +72,31 @@ namespace ModuleIO const std::complex beta(0.0, 0.0); #ifdef __MPI const int i1 = 1; - pzgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &i1, &i1, pv.desc, c, &i1, &i1, pv.desc_wfc, &beta, Vc.data(), &i1, &i1, pv.desc_wfc); + pzgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &i1, + &i1, pv.desc, c, &i1, &i1, pv.desc_wfc, &beta, Vc.data(), &i1, &i1, pv.desc_wfc); #else - zgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &nbasis, c, &nbasis, &beta, Vc.data(), &nbasis); + zgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &nbasis, + c, &nbasis, &beta, Vc.data(), &nbasis); #endif std::vector> cVc(p2d.nloc, 0.0); transa = 'C'; #ifdef __MPI - pzgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &i1, &i1, pv.desc_wfc, Vc.data(), &i1, &i1, pv.desc_wfc, &beta, cVc.data(), &i1, &i1, p2d.desc); + pzgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &i1, &i1, + pv.desc_wfc, Vc.data(), &i1, &i1, pv.desc_wfc, &beta, cVc.data(), &i1, &i1, p2d.desc); #else - zgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &nbasis, Vc.data(), &nbasis, &beta, cVc.data(), &nbasis); + zgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &nbasis, + Vc.data(), &nbasis, &beta, cVc.data(), &nbasis); #endif return cVc; } - std::vector cVc(double* V, double* c, int nbasis, int nbands, const Parallel_Orbitals& pv, const Parallel_2D& p2d) + std::vector cVc( + double* V, + double* c, + int nbasis, + int nbands, + const Parallel_Orbitals& pv, + const Parallel_2D& p2d) { std::vector Vc(pv.nloc_wfc, 0.0); char transa = 'N'; @@ -69,57 +105,98 @@ namespace ModuleIO const double beta = 0.0; #ifdef __MPI const int i1 = 1; - pdgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &i1, &i1, pv.desc, c, &i1, &i1, pv.desc_wfc, &beta, Vc.data(), &i1, &i1, pv.desc_wfc); + pdgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &i1, + &i1, pv.desc, c, &i1, &i1, pv.desc_wfc, &beta, Vc.data(), &i1, &i1, pv.desc_wfc); #else - dgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &nbasis, c, &nbasis, &beta, Vc.data(), &nbasis); + dgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &nbasis, + c, &nbasis, &beta, Vc.data(), &nbasis); #endif std::vector cVc(p2d.nloc, 0.0); transa = 'T'; #ifdef __MPI - pdgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &i1, &i1, pv.desc_wfc, Vc.data(), &i1, &i1, pv.desc_wfc, &beta, cVc.data(), &i1, &i1, p2d.desc); + pdgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &i1, &i1, + pv.desc_wfc, Vc.data(), &i1, &i1, pv.desc_wfc, &beta, cVc.data(), &i1, &i1, p2d.desc); #else - dgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &nbasis, Vc.data(), &nbasis, &beta, cVc.data(), &nbasis); + dgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, + &nbasis, Vc.data(), &nbasis, &beta, cVc.data(), &nbasis); #endif return cVc; } inline double get_real(const std::complex& c) { return c.real(); } + inline double get_real(const double& d) { return d; } template - double all_band_energy(const int ik, const std::vector& mat_mo, const Parallel_2D& p2d, const ModuleBase::matrix& wg) + double all_band_energy( + const int ik, + const std::vector& mat_mo, + const Parallel_2D& p2d, + const ModuleBase::matrix& wg) { double e = 0.0; for (int i = 0;i < p2d.get_row_size();++i) - for (int j = 0;j < p2d.get_col_size();++j) - if (p2d.local2global_row(i) == p2d.local2global_col(j)) - e += get_real(mat_mo[j * p2d.get_row_size() + i]) * wg(ik, p2d.local2global_row(i)); + { + for (int j = 0;j < p2d.get_col_size();++j) + { + if (p2d.local2global_row(i) == p2d.local2global_col(j)) + { + e += get_real(mat_mo[j * p2d.get_row_size() + i]) * wg(ik, p2d.local2global_row(i)); + } + } + } Parallel_Reduce::reduce_all(e); return e; } + // mohan update 2024-04-01 template - void set_gint_pointer(LCAO_Hamilt& uhm, typename TGint::type*& gint); + void set_gint_pointer( + Gint_Gamma &gint_gamma, + Gint_k &gint_k, + typename TGint::type*& gint); - template <> - void set_gint_pointer(LCAO_Hamilt& uhm, typename TGint::type*& gint) - { - gint = &uhm.GG; + // mohan update 2024-04-01 + template <> + void set_gint_pointer( + Gint_Gamma &gint_gamma, + Gint_k &gint_k, + typename TGint::type*& gint) + { + gint = &gint_gamma; } + + // mohan update 2024-04-01 template <> - void set_gint_pointer>(LCAO_Hamilt& uhm, typename TGint>::type*& gint) + void set_gint_pointer>( + Gint_Gamma &gint_gamma, + Gint_k &gint_k, + typename TGint>::type*& gint) { - gint = &uhm.GK; + gint = &gint_k; } - /// @brief write the Vxc matrix in KS orbital representation, usefull for GW calculation /// including terms: local/semi-local XC, EXX, DFTU template - void write_Vxc(int nspin, int nbasis, int drank, const psi::Psi& psi, const UnitCell& ucell, Structure_Factor& sf, - const ModulePW::PW_Basis& rho_basis, const ModulePW::PW_Basis& rhod_basis, const ModuleBase::matrix& vloc, - const Charge& chg, LCAO_Hamilt& uhm, LCAO_Matrix& lm, Local_Orbital_Charge& loc, - const K_Vectors& kv, const ModuleBase::matrix& wg, Grid_Driver& gd) + void write_Vxc( + int nspin, + int nbasis, + int drank, + const psi::Psi& psi, + const UnitCell& ucell, + Structure_Factor& sf, + const ModulePW::PW_Basis& rho_basis, + const ModulePW::PW_Basis& rhod_basis, + const ModuleBase::matrix& vloc, + const Charge& chg, + Gint_Gamma &gint_gamma, // mohan add 2024-04-01 + Gint_k &gint_k, // mohan add 2024-04-01 + LCAO_Matrix& lm, + Local_Orbital_Charge& loc, + const K_Vectors& kv, + const ModuleBase::matrix& wg, + Grid_Driver& gd) { ModuleBase::TITLE("ModuleIO", "write_Vxc"); const Parallel_Orbitals* pv = lm.ParaV; @@ -146,11 +223,24 @@ namespace ModuleIO // 3. allocate operators and contribute HR // op (corresponding to hR) typename TGint::type* gint = nullptr; - set_gint_pointer(uhm, gint); + + set_gint_pointer(gint_gamma, gint_k, gint); + std::vector>*> vxcs_op_ao(nspin0); for (int is = 0;is < nspin0;++is) { - vxcs_op_ao[is] = new hamilt::Veff>(gint, &loc, &lm, kv.kvec_d, potxc, &vxcs_R_ao[is], &vxc_k_ao, &ucell, &gd, pv); + vxcs_op_ao[is] = new hamilt::Veff>( + gint, + &loc, + &lm, + kv.kvec_d, + potxc, + &vxcs_R_ao[is], + &vxc_k_ao, + &ucell, + &gd, + pv); + GlobalV::CURRENT_SPIN = is; //caution: Veff::contributeHR depends on GlobalV::CURRENT_SPIN vxcs_op_ao[is]->contributeHR(); } @@ -177,7 +267,10 @@ namespace ModuleIO int is = GlobalV::CURRENT_SPIN = kv.isk[ik]; dynamic_cast*>(vxcs_op_ao[is])->contributeHk(ik); #ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) vexx_op_ao.contributeHk(ik); + if (GlobalC::exx_info.info_global.cal_exx) + { + vexx_op_ao.contributeHk(ik); + } // ======test======= // ModuleBase::GlobalFunc::ZEROS(test_vexxonly_k_ao.data(), pv->nloc); // if (GlobalC::exx_info.info_global.cal_exx) test_vexxonly_op_ao.contributeHk(ik); @@ -185,8 +278,11 @@ namespace ModuleIO // exx_energy += all_band_energy(ik, test_vexxonly_k_mo, p2d, wg); // ======test======= #endif - if (GlobalV::dft_plus_u) vdftu_op_ao.contributeHk(ik); - std::vector vxc_k_mo = cVc(vxc_k_ao.data(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); + if (GlobalV::dft_plus_u) + { + vdftu_op_ao.contributeHk(ik); + } + std::vector vxc_k_mo = cVc(vxc_k_ao.data(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); // write ModuleIO::save_mat(-1, vxc_k_mo.data(), nbands, false/*binary*/, GlobalV::out_ndigits, true/*triangle*/, false/*append*/, @@ -204,6 +300,9 @@ namespace ModuleIO // std::cout << "exx_energy: " << 0.5 * exx_energy << std::endl; // ======test======= delete potxc; - for (int is = 0;is < nspin0;++is) delete vxcs_op_ao[is]; + for (int is = 0;is < nspin0;++is) + { + delete vxcs_op_ao[is]; + } } -} \ No newline at end of file +}