diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 33276ebbdd..826e5f79c1 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -228,7 +228,6 @@ OBJS_ELECSTAT=elecstate.o\ pot_xc.o\ OBJS_ELECSTAT_LCAO=elecstate_lcao.o\ - elecstate_lcao_tddft.o\ elecstate_lcao_cal_tau.o\ density_matrix.o\ density_matrix_io.o\ diff --git a/source/module_elecstate/CMakeLists.txt b/source/module_elecstate/CMakeLists.txt index 01d0846ece..78efe4ead0 100644 --- a/source/module_elecstate/CMakeLists.txt +++ b/source/module_elecstate/CMakeLists.txt @@ -36,7 +36,6 @@ list(APPEND objects if(ENABLE_LCAO) list(APPEND objects elecstate_lcao.cpp - elecstate_lcao_tddft.cpp elecstate_lcao_cal_tau.cpp potentials/H_TDDFT_pw.cpp module_dm/density_matrix.cpp diff --git a/source/module_elecstate/elecstate_lcao_tddft.cpp b/source/module_elecstate/elecstate_lcao_tddft.cpp deleted file mode 100644 index 147b790a95..0000000000 --- a/source/module_elecstate/elecstate_lcao_tddft.cpp +++ /dev/null @@ -1,84 +0,0 @@ -#include "elecstate_lcao_tddft.h" - -#include "module_parameter/parameter.h" -#include "cal_dm.h" -#include "module_base/timer.h" -#include "module_elecstate/module_dm/cal_dm_psi.h" -#include "module_hamilt_lcao/module_gint/grid_technique.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -namespace elecstate -{ - -// multi-k case -void ElecStateLCAO_TDDFT::psiToRho_td(const psi::Psi>& psi) -{ - ModuleBase::TITLE("ElecStateLCAO", "psiToRho"); - ModuleBase::timer::tick("ElecStateLCAO", "psiToRho"); - - this->calculate_weights_td(); - this->calEBand(); - - ModuleBase::GlobalFunc::NOTE("Calculate the density matrix."); - - // this part for calculating DMK in 2d-block format, not used for charge now - // psi::Psi> dm_k_2d(); - elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), this->wg, psi, *(this->DM)); - this->DM->cal_DMR(); - - for (int is = 0; is < PARAM.inp.nspin; is++) - { - ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], this->charge->nrxx); // mohan 2009-11-10 - } - - //------------------------------------------------------------ - // calculate the charge density on real space grid. - //------------------------------------------------------------ - - ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!"); - this->gint_k->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint - Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho, PARAM.inp.nspin); // rho calculation - this->gint_k->cal_gint(&inout); - - this->charge->renormalize_rho(); - - ModuleBase::timer::tick("ElecStateLCAO", "psiToRho"); - return; -} - -void ElecStateLCAO_TDDFT::calculate_weights_td() -{ - ModuleBase::TITLE("ElecState", "calculate_weights"); - - if (PARAM.inp.ocp == 1) - { - int num = 0; - num = this->klist->get_nks() * PARAM.inp.nbands; - if (num != PARAM.inp.ocp_kb.size()) - { - ModuleBase::WARNING_QUIT("ElecStateLCAO_TDDFT::calculate_weights_td", - "size of occupation array is wrong , please check ocp_set"); - } - - double num_elec = 0.0; - for (int i = 0; i < PARAM.inp.ocp_kb.size(); i++) - { - num_elec += PARAM.inp.ocp_kb[i]; - } - if (std::abs(num_elec - PARAM.inp.nelec) > 1.0e-5) - { - ModuleBase::WARNING_QUIT("ElecStateLCAO_TDDFT::calculate_weights_td", - "total number of occupations is wrong , please check ocp_set"); - } - - for (int ik = 0; ik < this->klist->get_nks(); ik++) - { - for (int ib = 0; ib < PARAM.inp.nbands; ib++) - { - this->wg(ik, ib) = PARAM.inp.ocp_kb[ik * PARAM.inp.nbands + ib]; - } - } - } - return; -} - -} // namespace elecstate diff --git a/source/module_elecstate/elecstate_lcao_tddft.h b/source/module_elecstate/elecstate_lcao_tddft.h deleted file mode 100644 index 642888ce35..0000000000 --- a/source/module_elecstate/elecstate_lcao_tddft.h +++ /dev/null @@ -1,30 +0,0 @@ -#ifndef W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_ELECSTATE_ELECSTATE_LCAO_TDDFT_H -#define W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_ELECSTATE_ELECSTATE_LCAO_TDDFT_H - -#include "elecstate.h" -#include "elecstate_lcao.h" -#include "module_hamilt_lcao/module_gint/gint_k.h" - -namespace elecstate -{ -class ElecStateLCAO_TDDFT : public ElecStateLCAO> -{ - public: - ElecStateLCAO_TDDFT(Charge* chg_in, - const K_Vectors* klist_in, - int nks_in, - Gint_k* gint_k_in, // mohan add 2024-04-01 - ModulePW::PW_Basis* rhopw_in, - ModulePW::PW_Basis_Big* bigpw_in) - { - init_ks(chg_in, klist_in, nks_in, rhopw_in, bigpw_in); - this->gint_k = gint_k_in; - this->classname = "ElecStateLCAO_TDDFT"; - } - void psiToRho_td(const psi::Psi>& psi); - void calculate_weights_td(); -}; - -} // namespace elecstate - -#endif diff --git a/source/module_esolver/esolver_gets.h b/source/module_esolver/esolver_gets.h index 6565d9d3a1..f675d66051 100644 --- a/source/module_esolver/esolver_gets.h +++ b/source/module_esolver/esolver_gets.h @@ -48,7 +48,7 @@ class ESolver_GetS : public ESolver_KS TwoCenterBundle two_center_bundle_; - // // temporary introduced during removing GlobalC::ORB + // temporary introduced during removing GlobalC::ORB LCAO_Orbitals orb_; }; } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index afd387a1ff..0248fa7e4d 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -222,7 +222,8 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell #endif // 12) set occupations - if (PARAM.inp.ocp) + // tddft does not need to set occupations in the first scf + if (PARAM.inp.ocp && inp.esolver_type != "tddft") { this->pelec->fixed_weights(PARAM.inp.ocp_kb, PARAM.inp.nbands, PARAM.inp.nelec); } diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index fd6767e1e3..942d4b6a15 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -13,7 +13,9 @@ #include "module_base/lapack_connector.h" #include "module_base/scalapack_connector.h" #include "module_elecstate/module_charge/symmetry_rho.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" #include "module_elecstate/module_dm/cal_edm_tddft.h" +#include "module_elecstate/module_dm/density_matrix.h" #include "module_elecstate/occupy.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need divide_HS_in_frag #include "module_hamilt_lcao/module_tddft/evolve_elec.h" @@ -23,7 +25,6 @@ //-----HSolver ElecState Hamilt-------- #include "module_elecstate/elecstate_lcao.h" -#include "module_elecstate/elecstate_lcao_tddft.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" #include "module_hsolver/hsolver_lcao.h" #include "module_parameter/parameter.h" @@ -66,57 +67,11 @@ ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() void ESolver_KS_LCAO_TDDFT::before_all_runners(const Input_para& inp, UnitCell& ucell) { - // 1) run "before_all_runners" in ESolver_KS - ESolver_KS::before_all_runners(inp, ucell); - - // 2) initialize the local pseudopotential with plane wave basis set - GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, pw_rho); - - // 3) initialize the electronic states for TDDFT - if (this->pelec == nullptr) - { - this->pelec = new elecstate::ElecStateLCAO_TDDFT(&this->chr, - &kv, - kv.get_nks(), - &this->GK, // mohan add 2024-04-01 - this->pw_rho, - pw_big); - } - - // 4) read the local orbitals and construct the interpolation tables. - // initialize the pv - LCAO_domain::init_basis_lcao(this->pv, - inp.onsite_radius, - inp.lcao_ecut, - inp.lcao_dk, - inp.lcao_dr, - inp.lcao_rmax, - ucell, - two_center_bundle_, - orb_); - - // 5) allocate H and S matrices according to computational resources - LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, this->pv, kv.get_nks(), orb_); - - // 6) initialize Density Matrix - dynamic_cast>*>(this->pelec) - ->init_DM(&kv, &this->pv, PARAM.inp.nspin); - - // 8) initialize the charge density - this->pelec->charge->allocate(PARAM.inp.nspin); - this->pelec->omega = GlobalC::ucell.omega; // this line is very odd. - - // 9) initializee the potential - this->pelec->pot = new elecstate::Potential(pw_rhod, - pw_rho, - &GlobalC::ucell, - &(GlobalC::ppcell.vloc), - &(sf), - &(pelec->f_en.etxc), - &(pelec->f_en.vtxc)); + // 1) run before_all_runners in ESolver_KS_LCAO + ESolver_KS_LCAO, double>::before_all_runners(inp, ucell); // this line should be optimized - this->pelec_td = dynamic_cast(this->pelec); + // this->pelec = dynamic_cast(this->pelec); } void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int iter, const double ethr) @@ -134,13 +89,13 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int ite this->psi_laststep, this->Hk_laststep, this->Sk_laststep, - this->pelec_td->ekb, + this->pelec->ekb, td_htype, PARAM.inp.propagator, kv.get_nks()); - this->pelec_td->psiToRho_td(this->psi[0]); + this->weight_dm_rho(); } - this->pelec_td->psiToRho_td(this->psi[0]); + this->weight_dm_rho(); } else if (istep >= 2) { @@ -153,11 +108,11 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int ite this->psi_laststep, this->Hk_laststep, this->Sk_laststep, - this->pelec_td->ekb, + this->pelec->ekb, td_htype, PARAM.inp.propagator, kv.get_nks()); - this->pelec_td->psiToRho_td(this->psi[0]); + this->weight_dm_rho(); } else { @@ -168,7 +123,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int ite { bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false; hsolver::HSolverLCAO> hsolver_lcao_obj(&this->pv, PARAM.inp.ks_solver); - hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec_td, skip_charge); + hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, skip_charge); } } @@ -203,8 +158,7 @@ void ESolver_KS_LCAO_TDDFT::iter_finish(const int istep, int& iter) for (int ib = 0; ib < PARAM.inp.nbands; ib++) { std::setprecision(6); - GlobalV::ofs_running << ik + 1 << " " << ib + 1 << " " << this->pelec_td->wg(ik, ib) - << std::endl; + GlobalV::ofs_running << ik + 1 << " " << ib + 1 << " " << this->pelec->wg(ik, ib) << std::endl; } } GlobalV::ofs_running << std::endl; @@ -379,7 +333,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) for (int ib = 0; ib < PARAM.inp.nbands; ib++) { GlobalV::ofs_running << ik + 1 << " " << ib + 1 << " " - << this->pelec_td->ekb(ik, ib) * ModuleBase::Ry_to_eV << std::endl; + << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV << std::endl; } } GlobalV::ofs_running << std::endl; @@ -417,4 +371,21 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) ESolver_KS_LCAO, double>::after_scf(istep); } +void ESolver_KS_LCAO_TDDFT::weight_dm_rho() +{ + if (PARAM.inp.ocp == 1) + { + this->pelec->fixed_weights(PARAM.inp.ocp_kb, PARAM.inp.nbands, PARAM.inp.nelec); + } + this->pelec->calEBand(); + + ModuleBase::GlobalFunc::NOTE("Calculate the density matrix."); + + auto _pes = dynamic_cast>*>(this->pelec); + elecstate::cal_dm_psi(_pes->DM->get_paraV_pointer(), _pes->wg, this->psi[0], *(_pes->DM)); + _pes->DM->cal_DMR(); + + this->pelec->psiToRho(this->psi[0]); +} + } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h index 5dac415a72..c36e41cb9d 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.h +++ b/source/module_esolver/esolver_ks_lcao_tddft.h @@ -2,7 +2,6 @@ #define ESOLVER_KS_LCAO_TDDFT_H #include "esolver_ks.h" #include "esolver_ks_lcao.h" -#include "module_elecstate/elecstate_lcao_tddft.h" #include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h" #include "module_psi/psi.h" @@ -36,10 +35,10 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO, doubl //! Overlap matrix of last time step std::complex** Sk_laststep = nullptr; - //! Electronic states of rt-TDDFT - elecstate::ElecStateLCAO_TDDFT* pelec_td = nullptr; - int td_htype = 1; + + private: + void weight_dm_rho(); }; } // namespace ModuleESolver