diff --git a/source/Makefile.Objects b/source/Makefile.Objects index ba4370a9f9..504349de48 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -233,6 +233,7 @@ OBJS_ELECSTAT_LCAO=elecstate_lcao.o\ density_matrix.o\ density_matrix_io.o\ cal_dm_psi.o\ + cal_edm_tddft.o\ OBJS_ESOLVER=esolver.o\ esolver_ks.o\ @@ -259,7 +260,6 @@ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\ lcao_others.o\ lcao_init_after_vc.o\ lcao_fun.o\ - cal_edm_tddft.o\ OBJS_GINT=gint.o\ gint_gamma_env.o\ diff --git a/source/module_elecstate/CMakeLists.txt b/source/module_elecstate/CMakeLists.txt index 37b5da39ff..01d0846ece 100644 --- a/source/module_elecstate/CMakeLists.txt +++ b/source/module_elecstate/CMakeLists.txt @@ -42,6 +42,7 @@ if(ENABLE_LCAO) module_dm/density_matrix.cpp module_dm/density_matrix_io.cpp module_dm/cal_dm_psi.cpp + module_dm/cal_edm_tddft.cpp ) endif() diff --git a/source/module_esolver/cal_edm_tddft.cpp b/source/module_elecstate/module_dm/cal_edm_tddft.cpp similarity index 61% rename from source/module_esolver/cal_edm_tddft.cpp rename to source/module_elecstate/module_dm/cal_edm_tddft.cpp index a9cdd2c051..f3dd47ffe5 100644 --- a/source/module_esolver/cal_edm_tddft.cpp +++ b/source/module_elecstate/module_dm/cal_edm_tddft.cpp @@ -1,81 +1,46 @@ -#include "esolver_ks_lcao_tddft.h" - -#include "module_io/cal_r_overlap_R.h" -#include "module_io/dipole_io.h" -#include "module_io/td_current_io.h" -#include "module_io/write_HS.h" -#include "module_io/write_HS_R.h" -#include "module_io/write_wfc_nao.h" - -//--------------temporary---------------------------- -#include "module_base/blas_connector.h" -#include "module_base/global_function.h" -#include "module_base/scalapack_connector.h" +#include "cal_edm_tddft.h" + #include "module_base/lapack_connector.h" -#include "module_elecstate/module_charge/symmetry_rho.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" -#include "module_hamilt_lcao/module_tddft/td_velocity.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_io/print_info.h" - -//-----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" -#include "module_psi/psi.h" - -//-----force& stress------------------- -#include "module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h" - -//--------------------------------------------------- - -namespace ModuleESolver +#include "module_base/scalapack_connector.h" +namespace elecstate { -// use the original formula (Hamiltonian matrix) to calculate energy density -// matrix -void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() +// use the original formula (Hamiltonian matrix) to calculate energy density matrix +void cal_edm_tddft(Parallel_Orbitals& pv, + elecstate::ElecState* pelec, + K_Vectors& kv, + hamilt::Hamilt>* p_hamilt) { // mohan add 2024-03-27 const int nlocal = PARAM.globalv.nlocal; assert(nlocal >= 0); - dynamic_cast>*>(this->pelec) - ->get_DM() - ->EDMK.resize(kv.get_nks()); - for (int ik = 0; ik < kv.get_nks(); ++ik) { + auto _pelec = dynamic_cast>*>(pelec); - p_hamilt->updateHk(ik); + _pelec->get_DM()->EDMK.resize(kv.get_nks()); - std::complex* tmp_dmk - = dynamic_cast>*>(this->pelec)->get_DM()->get_DMK_pointer(ik); - - ModuleBase::ComplexMatrix& tmp_edmk - = dynamic_cast>*>(this->pelec)->get_DM()->EDMK[ik]; - - const Parallel_Orbitals* tmp_pv - = dynamic_cast>*>(this->pelec)->get_DM()->get_paraV_pointer(); + for (int ik = 0; ik < kv.get_nks(); ++ik) + { + p_hamilt->updateHk(ik); + std::complex* tmp_dmk = _pelec->get_DM()->get_DMK_pointer(ik); + ModuleBase::ComplexMatrix& tmp_edmk = _pelec->get_DM()->EDMK[ik]; #ifdef __MPI // mohan add 2024-03-27 //! be careful, the type of nloc is 'long' //! whether the long type is safe, needs more discussion - const long nloc = this->pv.nloc; - const int ncol = this->pv.ncol; - const int nrow = this->pv.nrow; + const long nloc = pv.nloc; + const int ncol = pv.ncol; + const int nrow = pv.nrow; tmp_edmk.create(ncol, nrow); - complex* Htmp = new complex[nloc]; - complex* Sinv = new complex[nloc]; - complex* tmp1 = new complex[nloc]; - complex* tmp2 = new complex[nloc]; - complex* tmp3 = new complex[nloc]; - complex* tmp4 = new complex[nloc]; + std::complex* Htmp = new std::complex[nloc]; + std::complex* Sinv = new std::complex[nloc]; + std::complex* tmp1 = new std::complex[nloc]; + std::complex* tmp2 = new std::complex[nloc]; + std::complex* tmp3 = new std::complex[nloc]; + std::complex* tmp4 = new std::complex[nloc]; ModuleBase::GlobalFunc::ZEROS(Htmp, nloc); ModuleBase::GlobalFunc::ZEROS(Sinv, nloc); @@ -86,8 +51,8 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() const int inc = 1; - hamilt::MatrixBlock> h_mat; - hamilt::MatrixBlock> s_mat; + hamilt::MatrixBlock> h_mat; + hamilt::MatrixBlock> s_mat; p_hamilt->matrix(h_mat, s_mat); zcopy_(&nloc, h_mat.p, &inc, Htmp, &inc); @@ -97,7 +62,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() int info = 0; const int one_int = 1; - pzgetrf_(&nlocal, &nlocal, Sinv, &one_int, &one_int, this->pv.desc, ipiv.data(), &info); + pzgetrf_(&nlocal, &nlocal, Sinv, &one_int, &one_int, pv.desc, ipiv.data(), &info); int lwork = -1; int liwork = -1; @@ -112,7 +77,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() Sinv, &one_int, &one_int, - this->pv.desc, + pv.desc, ipiv.data(), work.data(), &lwork, @@ -129,7 +94,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() Sinv, &one_int, &one_int, - this->pv.desc, + pv.desc, ipiv.data(), work.data(), &lwork, @@ -139,9 +104,9 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() const char N_char = 'N'; const char T_char = 'T'; - const complex one_float = {1.0, 0.0}; - const complex zero_float = {0.0, 0.0}; - const complex half_float = {0.5, 0.0}; + const std::complex one_float = {1.0, 0.0}; + const std::complex zero_float = {0.0, 0.0}; + const std::complex half_float = {0.5, 0.0}; pzgemm_(&N_char, &N_char, @@ -152,16 +117,16 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() Htmp, &one_int, &one_int, - this->pv.desc, + pv.desc, Sinv, &one_int, &one_int, - this->pv.desc, + pv.desc, &zero_float, tmp1, &one_int, &one_int, - this->pv.desc); + pv.desc); pzgemm_(&T_char, &N_char, @@ -172,16 +137,16 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() tmp1, &one_int, &one_int, - this->pv.desc, + pv.desc, tmp_dmk, &one_int, &one_int, - this->pv.desc, + pv.desc, &zero_float, tmp2, &one_int, &one_int, - this->pv.desc); + pv.desc); pzgemm_(&N_char, &N_char, @@ -192,16 +157,16 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() Sinv, &one_int, &one_int, - this->pv.desc, + pv.desc, Htmp, &one_int, &one_int, - this->pv.desc, + pv.desc, &zero_float, tmp3, &one_int, &one_int, - this->pv.desc); + pv.desc); pzgemm_(&N_char, &T_char, @@ -212,16 +177,16 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() tmp_dmk, &one_int, &one_int, - this->pv.desc, + pv.desc, tmp3, &one_int, &one_int, - this->pv.desc, + pv.desc, &zero_float, tmp4, &one_int, &one_int, - this->pv.desc); + pv.desc); pzgeadd_(&N_char, &nlocal, @@ -230,12 +195,12 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() tmp2, &one_int, &one_int, - this->pv.desc, + pv.desc, &half_float, tmp4, &one_int, &one_int, - this->pv.desc); + pv.desc); zcopy_(&nloc, tmp4, &inc, tmp_edmk.c, &inc); @@ -247,12 +212,12 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() delete[] tmp4; #else // for serial version - tmp_edmk.create(this->pv.ncol, this->pv.nrow); + tmp_edmk.create(pv.ncol, pv.nrow); ModuleBase::ComplexMatrix Sinv(nlocal, nlocal); ModuleBase::ComplexMatrix Htmp(nlocal, nlocal); - hamilt::MatrixBlock> h_mat; - hamilt::MatrixBlock> s_mat; + hamilt::MatrixBlock> h_mat; + hamilt::MatrixBlock> s_mat; p_hamilt->matrix(h_mat, s_mat); // cout<<"hmat "<>* p_hamilt); +} // namespace elecstate +#endif // CAL_EDM_TDDFT_H diff --git a/source/module_esolver/CMakeLists.txt b/source/module_esolver/CMakeLists.txt index 61968c7734..0a818e1464 100644 --- a/source/module_esolver/CMakeLists.txt +++ b/source/module_esolver/CMakeLists.txt @@ -26,7 +26,6 @@ if(ENABLE_LCAO) lcao_others.cpp lcao_init_after_vc.cpp lcao_fun.cpp - cal_edm_tddft.cpp ) endif() diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 2221e291ef..fd6767e1e3 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -13,6 +13,7 @@ #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_edm_tddft.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" @@ -358,7 +359,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) // calculate energy density matrix for tddft if (istep >= (wf.init_wfc == "file" ? 0 : 2) && module_tddft::Evolve_elec::td_edm == 0) { - this->cal_edm_tddft(); + elecstate::cal_edm_tddft(this->pv, this->pelec, this->kv, this->p_hamilt); } } diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h index 7589d77587..e672ccc25f 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.h +++ b/source/module_esolver/esolver_ks_lcao_tddft.h @@ -37,8 +37,6 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO, doubl virtual void iter_finish(const int istep, int& iter) override; virtual void after_scf(const int istep) override; - - void cal_edm_tddft(); }; } // namespace ModuleESolver