diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 1899cb9dfd..4276cb5fb3 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -41,6 +41,9 @@ ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() { // this->orb_con.clear_after_ions(GlobalC::UOT, GlobalC::ORB, GlobalV::deepks_setorb, GlobalC::ucell.infoNL.nproj); delete psi_laststep; + for(int ik = 0; ik < GlobalC::kv.nks; ++ik){ + delete Hk_laststep[ik];} + delete Hk_laststep; } void ESolver_KS_LCAO_TDDFT::Init(Input& inp, UnitCell& ucell) @@ -188,7 +191,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr) if (GlobalV::ESOLVER_TYPE == "tddft" && istep >= 2 && !GlobalV::GAMMA_ONLY_LOCAL) { - ELEC_evolve::evolve_psi(istep, this->p_hamilt, this->LOWF, this->psi, this->psi_laststep, this->pelec_td->ekb); + ELEC_evolve::evolve_psi(istep, this->p_hamilt, this->LOWF, this->psi, this->psi_laststep, this->Hk_laststep, this->pelec_td->ekb); this->pelec_td->psiToRho_td(this->psi[0]); // this->pelec_td->psiToRho(this->psi[0]); } @@ -351,7 +354,7 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) GlobalC::en.cal_converged(this->pelec); } - // store wfc + // store wfc and Hk laststep if (istep >= 1 && this->conv_elec) { if (this->psi_laststep == nullptr) @@ -364,6 +367,13 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) this->psi_laststep = new psi::Psi>(GlobalC::kv.nks, GlobalV::NBANDS, GlobalV::NLOCAL, nullptr); #endif + if (this->Hk_laststep == nullptr){ + this->Hk_laststep = new std::complex*[GlobalC::kv.nks]; + for(int ik = 0; ik < GlobalC::kv.nks; ++ik){ + this->Hk_laststep[ik] = new std::complex[this->LOC.ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(Hk_laststep[ik], this->LOC.ParaV->nloc); + } + } for (int ik = 0; ik < GlobalC::kv.nks; ++ik) { @@ -372,8 +382,15 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) int size0 = psi->get_nbands() * psi->get_nbasis(); for (int index = 0; index < size0; ++index) psi_laststep[0].get_pointer()[index] = psi[0].get_pointer()[index]; + //use for new propagator + this->p_hamilt->updateHk(ik); + //store Hk laststep + hamilt::MatrixBlock> h_mat, s_mat; + this->p_hamilt->matrix(h_mat, s_mat); + BlasConnector::copy(this->LOC.ParaV->nloc, h_mat.p, 1, Hk_laststep[ik], 1); } + if (istep > 1 && ELEC_evolve::td_edm == 0) this->cal_edm_tddft(); } diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h index a943cebe5a..6f9da4a185 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.h +++ b/source/module_esolver/esolver_ks_lcao_tddft.h @@ -22,6 +22,7 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO void Init(Input& inp, UnitCell& cell) override; psi::Psi>* psi_laststep = nullptr; + std::complex** Hk_laststep = nullptr; //same as pelec elecstate::ElecStateLCAO_TDDFT* pelec_td = nullptr; diff --git a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp index 66a8fcb1a7..3cfdd30b1d 100644 --- a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp +++ b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp @@ -25,6 +25,7 @@ void ELEC_evolve::evolve_psi(const int& istep, Local_Orbital_wfc& lowf, psi::Psi>* psi, psi::Psi>* psi_laststep, + std::complex** Hk_laststep, ModuleBase::matrix& ekb) { ModuleBase::TITLE("ELEC_evolve", "eveolve_psi"); @@ -39,7 +40,7 @@ void ELEC_evolve::evolve_psi(const int& istep, Evolve_LCAO_Matrix ELM(lowf.ParaV); psi->fix_k(ik); psi_laststep->fix_k(ik); - ELM.evolve_complex_matrix(ik, phm, psi, psi_laststep, &(ekb(ik, 0))); + ELM.evolve_complex_matrix(ik, phm, psi, psi_laststep, Hk_laststep[ik], &(ekb(ik, 0))); ModuleBase::timer::tick("Efficience", "evolve_k"); } // end k diff --git a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h index 6571e72033..eb647b83c3 100644 --- a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h +++ b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h @@ -44,6 +44,7 @@ class ELEC_evolve Local_Orbital_wfc& lowf, psi::Psi>* psi, psi::Psi>* psi_laststep, + std::complex** Hk_laststep, ModuleBase::matrix& ekb); }; diff --git a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp index 66325ee7cf..49b190f6af 100644 --- a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp +++ b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp @@ -27,6 +27,7 @@ void Evolve_LCAO_Matrix::evolve_complex_matrix(const int& ik, hamilt::Hamilt* p_hamilt, psi::Psi>* psi_k, psi::Psi>* psi_k_laststep, + std::complex* Hk_laststep, double* ekb) const { ModuleBase::TITLE("Evolve_LCAO_Matrix", "evolve_complex_matrix"); @@ -36,7 +37,7 @@ void Evolve_LCAO_Matrix::evolve_complex_matrix(const int& ik, if (GlobalV::ESOLVER_TYPE == "tddft") { #ifdef __MPI - this->using_ScaLAPACK_complex(GlobalV::NBANDS, GlobalV::NLOCAL, psi_k_laststep[0].get_pointer(), p_hamilt, psi_k[0].get_pointer(), ekb); + this->using_ScaLAPACK_complex(GlobalV::NBANDS, GlobalV::NLOCAL, psi_k_laststep[0].get_pointer(), Hk_laststep, p_hamilt, psi_k[0].get_pointer(), ekb); #else this->using_LAPACK_complex(ik, p_hamilt, psi_k[0].get_pointer(), psi_k_laststep[0].get_pointer(), ekb); #endif @@ -344,6 +345,7 @@ void Evolve_LCAO_Matrix::using_ScaLAPACK_complex( const int nband, const int nlocal, const std::complex* psi_k_laststep, + const std::complex* Hk_laststep, hamilt::Hamilt* p_hamilt, std::complex* psi_k, double* ekb) const @@ -361,7 +363,51 @@ void Evolve_LCAO_Matrix::using_ScaLAPACK_complex( complex* Htmp = new complex[this->ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(Htmp, this->ParaV->nloc); BlasConnector::copy(this->ParaV->nloc, h_mat.p, 1, Htmp, 1); - + + //print test + /*GlobalV::ofs_running << " H(t) matrix :" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << Hk_laststep[i * this->ParaV->ncol + j].real() << "+" << Hk_laststep[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl;} + GlobalV::ofs_running << " H(t+dt) matrix :" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << Htmp[i * this->ParaV->ncol + j].real() << "+" << Htmp[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl;} + // +*/ + + complex mixing = {0.5, 0.0}; + ScalapackConnector::geadd( + 'N', + nlocal, + nlocal, + mixing, + Hk_laststep, + 1, + 1, + this->ParaV->desc, + mixing, + Htmp, + 1, + 1, + this->ParaV->desc + ); + /*GlobalV::ofs_running << " H(t+dt/2) matrix :" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << Htmp[i * this->ParaV->ncol + j].real() << "+" << Htmp[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl;}*/ complex* U_operator = new complex[this->ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(U_operator, this->ParaV->nloc); diff --git a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h index dd13e45e2f..1f1f4e1fc0 100644 --- a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h +++ b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h @@ -22,6 +22,7 @@ class Evolve_LCAO_Matrix hamilt::Hamilt* p_hamilt, psi::Psi>* psi_k, psi::Psi>* psi_k_laststep, + std::complex* Hk_laststep, double* ekb) const; private: @@ -38,6 +39,7 @@ class Evolve_LCAO_Matrix const int nband, const int nlocal, const std::complex* psi_k_laststep, + const std::complex* Hk_laststep, hamilt::Hamilt* p_hamilt, std::complex* psi_k, double* ekb) const;