From ccb8a5d2ac62c1e94fd334d6a473c58ae972a395 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Fri, 26 Sep 2025 18:06:46 +0800 Subject: [PATCH 1/8] add esolver_of_tddft --- source/Makefile.Objects | 1 + source/source_esolver/esolver.cpp | 9 + source/source_esolver/esolver_of.h | 4 +- source/source_esolver/esolver_of_tddft.cpp | 404 ++++++++++++++++++++ source/source_esolver/esolver_of_tddft.h | 30 ++ source/source_io/read_input_item_system.cpp | 4 +- 6 files changed, 448 insertions(+), 4 deletions(-) create mode 100644 source/source_esolver/esolver_of_tddft.cpp create mode 100644 source/source_esolver/esolver_of_tddft.h diff --git a/source/Makefile.Objects b/source/Makefile.Objects index db5585dd3a..35b585c508 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -266,6 +266,7 @@ OBJS_ESOLVER=esolver.o\ esolver_lj.o\ esolver_dp.o\ esolver_of.o\ + esolver_of_tddft.o\ esolver_of_tool.o\ esolver_of_interface.o\ pw_others.o\ diff --git a/source/source_esolver/esolver.cpp b/source/source_esolver/esolver.cpp index 5d513b398c..d5b09bee12 100644 --- a/source/source_esolver/esolver.cpp +++ b/source/source_esolver/esolver.cpp @@ -20,6 +20,7 @@ extern "C" #include "esolver_dp.h" #include "esolver_lj.h" #include "esolver_of.h" +#include "esolver_of_tddft.h" #include "source_io/module_parameter/md_parameter.h" #include @@ -40,6 +41,10 @@ std::string determine_type() { esolver_type = "ofdft"; } + else if (PARAM.inp.esolver_type == "tdofdft") + { + esolver_type = "tdofdft"; + } else if (PARAM.inp.esolver_type == "ksdft") { esolver_type = "ksdft_pw"; @@ -321,6 +326,10 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell) { return new ESolver_OF(); } + else if (esolver_type == "tdofdft") + { + return new ESolver_OF_TDDFT(); + } else if (esolver_type == "lj_pot") { return new ESolver_LJ(); diff --git a/source/source_esolver/esolver_of.h b/source/source_esolver/esolver_of.h index 508a18b3db..f35808cf1b 100644 --- a/source/source_esolver/esolver_of.h +++ b/source/source_esolver/esolver_of.h @@ -27,7 +27,7 @@ class ESolver_OF : public ESolver_FP virtual void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) override; - private: + protected: // ======================= variables ========================== // ---------- the kinetic energy density functionals ---------- KEDF_Manager* kedf_manager_ = nullptr; // KEDF manager, which will be initialized in before_all_runners @@ -83,7 +83,7 @@ class ESolver_OF : public ESolver_FP // ============================ tools =============================== // --------------------- initialize --------------------------------- - void init_elecstate(UnitCell& ucell); + void init_elecstate(UnitCell& ucell); void allocate_array(); // --------------------- calculate physical qualities --------------- diff --git a/source/source_esolver/esolver_of_tddft.cpp b/source/source_esolver/esolver_of_tddft.cpp new file mode 100644 index 0000000000..3b2ecdb42d --- /dev/null +++ b/source/source_esolver/esolver_of_tddft.cpp @@ -0,0 +1,404 @@ +#include "esolver_of_tddft.h" + +#include "source_io/module_parameter/parameter.h" +#include "source_io/cube_io.h" +#include "source_io/output_log.h" +#include "source_io/write_elecstat_pot.h" +//-----------temporary------------------------- +#include "source_base/global_function.h" +#include "source_estate/module_charge/symmetry_rho.h" +#include "source_hamilt/module_ewald/H_Ewald_pw.h" +#include "source_pw/module_pwdft/global.h" +#include "source_io/print_info.h" +#include "source_estate/cal_ux.h" +//-----force------------------- +#include "source_pw/module_pwdft/forces.h" +//-----stress------------------ +#include "source_pw/module_ofdft/of_stress_pw.h" + +namespace ModuleESolver +{ + +ESolver_OF_TDDFT::ESolver_OF_TDDFT() +{ + this->classname = "ESolver_OF_TDDFT"; + /*this->pphi_td= new std::complex*[PARAM.inp.nspin]; + + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + this->pphi_td[is]= new std::complex[pw_rho->nrxx]; + }*/ +} + +ESolver_OF_TDDFT::~ESolver_OF_TDDFT() +{ + //delete psi_td; + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + delete[] this->pphi_td[is]; + } + delete[] this->pphi_td; +} + + +void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) +{ + ModuleBase::timer::tick("ESolver_OF_TDDFT", "runner"); + // get Ewald energy, initial rho and phi if necessary + this->before_opt(istep, ucell); + this->iter_ = 0; + + bool conv_esolver = false; // this conv_esolver is added by mohan 20250302 +#ifdef __MPI + this->iter_time = MPI_Wtime(); +#else + this->iter_time = std::chrono::system_clock::now(); +#endif + + if (istep==0) + { + this->pphi_td= new std::complex*[PARAM.inp.nspin]; + + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + this->pphi_td[is]= new std::complex[pw_rho->nrxx]; + } + } + + if ((istep<1) && PARAM.inp.init_chg != "file") + { + while (true) + { + // once we get a new rho and phi, update potential + this->update_potential(ucell); + + // calculate the energy of new rho and phi + this->energy_llast_ = this->energy_last_; + this->energy_last_ = this->energy_current_; + this->energy_current_ = this->cal_energy(); + + + // check if the job is done + if (this->check_exit(conv_esolver)) + { + break; + } + + // find the optimization direction and step lenghth theta according to the potential + this->optimize(ucell); + + std::cout<<"optimize------"<update_rho(); + + this->iter_++; + + ESolver_FP::iter_finish(ucell, istep, this->iter_, conv_esolver); + } + + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + pphi_td[is][ir]=pphi_[is][ir]; + } + } + } + else + { + std::cout<<"propagate_psi -------"<propagate_psi(ucell,this->pphi_td,this->pw_rho); + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + pphi_[is][ir]=abs(pphi_td[is][ir]); + } + } + conv_esolver=true; + } + + this->after_opt(istep, ucell, conv_esolver); + + ModuleBase::timer::tick("ESolver_OF_TDDFT", "runner"); +} + +/** + * @brief Prepare to optimize the charge density, + * update elecstate, kedf, and opts if needed + * calculate ewald energy, initialize the rho, phi, theta + * + * @param istep + * @param ucell + */ +void ESolver_OF_TDDFT::before_opt(const int istep, UnitCell& ucell) +{ + ModuleBase::TITLE("ESolver_OF", "before_opt"); + ModuleBase::timer::tick("ESolver_OF", "before_opt"); + + //! 1) call before_scf() of ESolver_FP + ESolver_FP::before_scf(ucell, istep); + + if (ucell.cell_parameter_updated) + { + this->dV_ = ucell.omega / this->pw_rho->nxyz; + + // initialize elecstate, including potential + this->init_elecstate(ucell); + + // Initialize KEDF + this->kedf_manager_->init(PARAM.inp, this->pw_rho, this->dV_, this->nelec_[0]); + + // Initialize optimization methods + this->init_opt(); + + // Refresh the arrays + delete this->psi_; + delete this->psi_td; + this->psi_ = new psi::Psi(1, + PARAM.inp.nspin, + this->pw_rho->nrxx, + this->pw_rho->nrxx, + true); + /*this->psi_td = new psi::Psi>(1, + PARAM.inp.nspin, + this->pw_rho->nrxx, + this->pw_rho->nrxx, + true);*/ + //this->pphi_td= new std::complex*[PARAM.inp.nspin]; + + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + this->pphi_[is] = this->psi_->get_pointer(is); + //this->pphi_td[is] = this->psi_td->get_pointer(is); + //this->pphi_td[is]= new std::complex[pw_rho->nrxx]; + } + + delete this->ptemp_rho_; + this->ptemp_rho_ = new Charge(); + this->ptemp_rho_->set_rhopw(this->pw_rho); + this->ptemp_rho_->allocate(PARAM.inp.nspin); + + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + delete[] this->pdLdphi_[is]; + delete[] this->pdEdphi_[is]; + delete[] this->pdirect_[is]; + delete[] this->precip_dir_[is]; + this->pdLdphi_[is] = new double[this->pw_rho->nrxx]; + this->pdEdphi_[is] = new double[this->pw_rho->nrxx]; + this->pdirect_[is] = new double[this->pw_rho->nrxx]; + this->precip_dir_[is] = new std::complex[pw_rho->npw]; + } + } + + this->pelec->init_scf(istep, ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm); + + Symmetry_rho srho; + for (int is = 0; is < PARAM.inp.nspin; is++) + { + srho.begin(is, this->chr, this->pw_rho, ucell.symm); + } + + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + if (PARAM.inp.init_chg != "file") + { + for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs) + { + // Here we initialize rho to be uniform, + // because the rho got by pot.init_pot -> Charge::atomic_rho may contain minus elements. + this->chr.rho[is][ibs] = this->nelec_[is] / this->pelec->omega; + this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]); + } + } + else + { + for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs) + { + this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]); + } + } + } + + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + this->pelec->eferm.set_efval(is, 0); + this->theta_[is] = 0.; + ModuleBase::GlobalFunc::ZEROS(this->pdLdphi_[is], this->pw_rho->nrxx); + ModuleBase::GlobalFunc::ZEROS(this->pdEdphi_[is], this->pw_rho->nrxx); + ModuleBase::GlobalFunc::ZEROS(this->pdirect_[is], this->pw_rho->nrxx); + } + if (PARAM.inp.nspin == 1) + { + this->theta_[0] = 0.2; + } + + ModuleBase::timer::tick("ESolver_OF", "before_opt"); +} + +void ESolver_OF_TDDFT::get_Hpsi(UnitCell& ucell, const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi) +{ + // update rho + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + this->chr.rho[is][ir] = abs(psi_[is][ir])*abs(psi_[is][ir]); + } + } + + this->pelec->pot->update_from_charge(&this->chr, &ucell); // Hartree + XC + external + this->get_tf_potential(this->chr.rho,pw_rho ,this->pelec->pot->get_effective_v()); // TF potential + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + const double* vr_eff = this->pelec->pot->get_effective_v(is); + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + Hpsi[is][ir] = vr_eff[ir]*psi_[is][ir]; + } + } + this->get_vw_potential_phi(psi_, pw_rho, Hpsi); +} + +void ESolver_OF_TDDFT::get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) +{ + const double c_tf_ + = 3.0 / 10.0 * std::pow(3 * std::pow(M_PI, 2.0), 2.0 / 3.0) + * 2; + if (PARAM.inp.nspin == 1) + { + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + rpot(0, ir) += 5.0 / 3.0 * c_tf_ * std::pow(prho[0][ir], 2. / 3.); + } + } + else if (PARAM.inp.nspin == 2) + { + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + rpot(is, ir) += 5.0 / 3.0 * c_tf_ * std::pow(2. * prho[is][ir], 2. / 3.); + } + } + } +} + +void ESolver_OF_TDDFT::get_vw_potential_phi(const std::complex* const* pphi, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi) +{ + std::complex** rLapPhi = new std::complex*[PARAM.inp.nspin]; + for (int is = 0; is < PARAM.inp.nspin; ++is) { + rLapPhi[is] = new std::complex[pw_rho->nrxx]; + } + std::complex** recipPhi = new std::complex*[PARAM.inp.nspin]; + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + recipPhi[is] = new std::complex[pw_rho->npw]; + + pw_rho->real2recip(pphi[is], recipPhi[is]); + for (int ik = 0; ik < pw_rho->npw; ++ik) + { + recipPhi[is][ik] *= pw_rho->gg[ik] * pw_rho->tpiba2; + } + pw_rho->recip2real(recipPhi[is], rLapPhi[is]); + for (int ik = 0; ik < pw_rho->npw; ++ik) + { + Hpsi[is][ik]+=rLapPhi[is][ik]; + } + } + + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + delete[] recipPhi[is]; + delete[] rLapPhi[is]; + } + delete[] recipPhi; + delete[] rLapPhi; +} + +void ESolver_OF_TDDFT::get_CD_potential(const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) +{ + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + //recipCurrent = new std::complex[pw_rho->npw]; + //delete[] recipCurrent; + } +} + +void ESolver_OF_TDDFT::propagate_psi(UnitCell& ucell, std::complex** pphi_, ModulePW::PW_Basis* pw_rho) +{ + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); + + std::complex imag(0.0,1.0); + double dt=PARAM.inp.mdp.md_dt; + std::complex** K1 = new std::complex*[PARAM.inp.nspin]; + std::complex** K2 = new std::complex*[PARAM.inp.nspin]; + std::complex** K3 = new std::complex*[PARAM.inp.nspin]; + std::complex** K4 = new std::complex*[PARAM.inp.nspin]; + std::complex** psi1 = new std::complex*[PARAM.inp.nspin]; + std::complex** psi2 = new std::complex*[PARAM.inp.nspin]; + std::complex** psi3 = new std::complex*[PARAM.inp.nspin]; + for (int is = 0; is < PARAM.inp.nspin; ++is) { + K1[is] = new std::complex[pw_rho->nrxx]; + K2[is] = new std::complex[pw_rho->nrxx]; + K3[is] = new std::complex[pw_rho->nrxx]; + K4[is] = new std::complex[pw_rho->nrxx]; + psi1[is] = new std::complex[pw_rho->nrxx]; + psi2[is] = new std::complex[pw_rho->nrxx]; + psi3[is] = new std::complex[pw_rho->nrxx]; + } + get_Hpsi(ucell,pphi_,pw_rho,K1); + for (int is = 0; is < PARAM.inp.nspin; ++is){ + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + K1[is][ir]=-1.0*K1[is][ir]*dt*imag; + psi1[is][ir]=pphi_[is][ir]+0.5*K1[is][ir]; + } + } + get_Hpsi(ucell,psi1,pw_rho,K2); + for (int is = 0; is < PARAM.inp.nspin; ++is){ + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + K2[is][ir]=-1.0*K2[is][ir]*dt*imag; + psi2[is][ir]=pphi_[is][ir]+0.5*K2[is][ir]; + } + } + get_Hpsi(ucell,psi2,pw_rho,K3); + for (int is = 0; is < PARAM.inp.nspin; ++is){ + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + K3[is][ir]=-1.0*K3[is][ir]*dt*imag; + psi3[is][ir]=pphi_[is][ir]+K3[is][ir]; + } + } + get_Hpsi(ucell,psi3,pw_rho,K4); + for (int is = 0; is < PARAM.inp.nspin; ++is){ + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + K4[is][ir]=-1.0*K4[is][ir]*dt*imag; + pphi_[is][ir]+=1.0/6.0*(K1[is][ir]+2.0*K2[is][ir]+2.0*K3[is][ir]+K4[is][ir]); + } + } + + for (int is = 0; is < PARAM.inp.nspin; ++is) { + delete[] K1[is]; + delete[] K2[is]; + delete[] K3[is]; + delete[] K4[is]; + delete[] psi1[is]; + delete[] psi2[is]; + delete[] psi3[is]; + } + delete[] K1; + delete[] K2; + delete[] K3; + delete[] K4; + delete[] psi1; + delete[] psi2; + delete[] psi3; + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); +} + +} // namespace ModuleESolver diff --git a/source/source_esolver/esolver_of_tddft.h b/source/source_esolver/esolver_of_tddft.h new file mode 100644 index 0000000000..04c15eb9c3 --- /dev/null +++ b/source/source_esolver/esolver_of_tddft.h @@ -0,0 +1,30 @@ +#ifndef ESOLVER_OF_TDDFT_H +#define ESOLVER_OF_TDDFT_H + +#include "esolver_of.h" + +namespace ModuleESolver +{ +class ESolver_OF_TDDFT : public ESolver_OF +{ + public: + ESolver_OF_TDDFT(); + ~ESolver_OF_TDDFT(); + + virtual void runner(UnitCell& ucell, const int istep) override; + + protected: + std::complex** pphi_td = nullptr; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). + psi::Psi>* psi_td = nullptr; + //std::complex** pdEdphi_phi_ = nullptr; // (dE/dphi)*phi + // ==================== main process of TDOFDFT ====================== + void before_opt(const int istep, UnitCell& ucell); + void get_Hpsi(UnitCell& ucell, const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi); + void get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); + void get_vw_potential_phi(const std::complex* const* pphi, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi); // -1/2 \nabla^2 \phi + void get_CD_potential(const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); + void propagate_psi(UnitCell& ucell, std::complex** pphi_, ModulePW::PW_Basis* pw_rho); +}; +} // namespace ModuleESolver + +#endif diff --git a/source/source_io/read_input_item_system.cpp b/source/source_io/read_input_item_system.cpp index c9854c9ece..bd123b4d51 100644 --- a/source/source_io/read_input_item_system.cpp +++ b/source/source_io/read_input_item_system.cpp @@ -108,10 +108,10 @@ void ReadInput::item_system() } { Input_Item item("esolver_type"); - item.annotation = "the energy solver: ksdft, sdft, ofdft, tddft, lj, dp, ks-lr, lr"; + item.annotation = "the energy solver: ksdft, sdft, ofdft, tdofdft, tddft, lj, dp, ks-lr, lr"; read_sync_string(input.esolver_type); item.check_value = [](const Input_Item& item, const Parameter& para) { - const std::vector esolver_types = { "ksdft", "sdft", "ofdft", "tddft", "lj", "dp", "lr", "ks-lr" }; + const std::vector esolver_types = { "ksdft", "sdft", "ofdft", "tdofdft", "tddft", "lj", "dp", "lr", "ks-lr" }; if (std::find(esolver_types.begin(), esolver_types.end(), para.input.esolver_type) == esolver_types.end()) { const std::string warningstr = nofound_str(esolver_types, "esolver_type"); From 96dd8cc69e27d76f5794f1b15f6490538d2fc1df Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Fri, 26 Sep 2025 19:13:03 +0800 Subject: [PATCH 2/8] esolver_of_tddft cmakelist --- source/source_esolver/CMakeLists.txt | 1 + source/source_esolver/esolver_of_tddft.cpp | 1 - source/source_esolver/esolver_of_tddft.h | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) diff --git a/source/source_esolver/CMakeLists.txt b/source/source_esolver/CMakeLists.txt index f58e47b3aa..acfd72e009 100644 --- a/source/source_esolver/CMakeLists.txt +++ b/source/source_esolver/CMakeLists.txt @@ -8,6 +8,7 @@ list(APPEND objects esolver_lj.cpp esolver_dp.cpp esolver_of.cpp + esolver_of_tddft.cpp esolver_of_interface.cpp esolver_of_tool.cpp pw_others.cpp diff --git a/source/source_esolver/esolver_of_tddft.cpp b/source/source_esolver/esolver_of_tddft.cpp index 3b2ecdb42d..5aee3d53d9 100644 --- a/source/source_esolver/esolver_of_tddft.cpp +++ b/source/source_esolver/esolver_of_tddft.cpp @@ -107,7 +107,6 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) } else { - std::cout<<"propagate_psi -------"<propagate_psi(ucell,this->pphi_td,this->pw_rho); for (int is = 0; is < PARAM.inp.nspin; ++is) { diff --git a/source/source_esolver/esolver_of_tddft.h b/source/source_esolver/esolver_of_tddft.h index 04c15eb9c3..2cf30a9c66 100644 --- a/source/source_esolver/esolver_of_tddft.h +++ b/source/source_esolver/esolver_of_tddft.h @@ -18,7 +18,7 @@ class ESolver_OF_TDDFT : public ESolver_OF psi::Psi>* psi_td = nullptr; //std::complex** pdEdphi_phi_ = nullptr; // (dE/dphi)*phi // ==================== main process of TDOFDFT ====================== - void before_opt(const int istep, UnitCell& ucell); + void before_opt(const int istep, UnitCell& ucell); void get_Hpsi(UnitCell& ucell, const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi); void get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); void get_vw_potential_phi(const std::complex* const* pphi, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi); // -1/2 \nabla^2 \phi From c87639c185b1139c799aad23110bbfaee4c8e05f Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sat, 27 Sep 2025 04:35:33 +0800 Subject: [PATCH 3/8] move memeber function of esolver_of_tddft to module_ofdft --- source/Makefile.Objects | 1 + source/source_esolver/esolver_of_tddft.cpp | 288 +------------------ source/source_esolver/esolver_of_tddft.h | 11 +- source/source_pw/module_ofdft/CMakeLists.txt | 1 + source/source_pw/module_ofdft/evolve_phi.cpp | 166 +++++++++++ source/source_pw/module_ofdft/evolve_phi.h | 40 +++ 6 files changed, 213 insertions(+), 294 deletions(-) create mode 100644 source/source_pw/module_ofdft/evolve_phi.cpp create mode 100644 source/source_pw/module_ofdft/evolve_phi.h diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 35b585c508..460a474ea5 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -361,6 +361,7 @@ OBJS_HAMILT_OF=kedf_tf.o\ kedf_xwm.o\ kedf_lkt.o\ kedf_manager.o\ + evolve_phi.o\ OBJS_HAMILT_LCAO=hamilt_lcao.o\ operator_lcao.o\ diff --git a/source/source_esolver/esolver_of_tddft.cpp b/source/source_esolver/esolver_of_tddft.cpp index 5aee3d53d9..7585d42ed2 100644 --- a/source/source_esolver/esolver_of_tddft.cpp +++ b/source/source_esolver/esolver_of_tddft.cpp @@ -22,17 +22,12 @@ namespace ModuleESolver ESolver_OF_TDDFT::ESolver_OF_TDDFT() { this->classname = "ESolver_OF_TDDFT"; - /*this->pphi_td= new std::complex*[PARAM.inp.nspin]; - - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - this->pphi_td[is]= new std::complex[pw_rho->nrxx]; - }*/ + this->evolve_psi=new EVOLVE_PSI(); } ESolver_OF_TDDFT::~ESolver_OF_TDDFT() { - //delete psi_td; + delete this->evolve_psi; for (int is = 0; is < PARAM.inp.nspin; ++is) { delete[] this->pphi_td[is]; @@ -107,7 +102,7 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) } else { - this->propagate_psi(ucell,this->pphi_td,this->pw_rho); + this->evolve_psi->propagate_psi(this->pelec, this->chr, ucell, this->pphi_td, this->pw_rho); for (int is = 0; is < PARAM.inp.nspin; ++is) { for (int ir = 0; ir < pw_rho->nrxx; ++ir) @@ -123,281 +118,4 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) ModuleBase::timer::tick("ESolver_OF_TDDFT", "runner"); } -/** - * @brief Prepare to optimize the charge density, - * update elecstate, kedf, and opts if needed - * calculate ewald energy, initialize the rho, phi, theta - * - * @param istep - * @param ucell - */ -void ESolver_OF_TDDFT::before_opt(const int istep, UnitCell& ucell) -{ - ModuleBase::TITLE("ESolver_OF", "before_opt"); - ModuleBase::timer::tick("ESolver_OF", "before_opt"); - - //! 1) call before_scf() of ESolver_FP - ESolver_FP::before_scf(ucell, istep); - - if (ucell.cell_parameter_updated) - { - this->dV_ = ucell.omega / this->pw_rho->nxyz; - - // initialize elecstate, including potential - this->init_elecstate(ucell); - - // Initialize KEDF - this->kedf_manager_->init(PARAM.inp, this->pw_rho, this->dV_, this->nelec_[0]); - - // Initialize optimization methods - this->init_opt(); - - // Refresh the arrays - delete this->psi_; - delete this->psi_td; - this->psi_ = new psi::Psi(1, - PARAM.inp.nspin, - this->pw_rho->nrxx, - this->pw_rho->nrxx, - true); - /*this->psi_td = new psi::Psi>(1, - PARAM.inp.nspin, - this->pw_rho->nrxx, - this->pw_rho->nrxx, - true);*/ - //this->pphi_td= new std::complex*[PARAM.inp.nspin]; - - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - this->pphi_[is] = this->psi_->get_pointer(is); - //this->pphi_td[is] = this->psi_td->get_pointer(is); - //this->pphi_td[is]= new std::complex[pw_rho->nrxx]; - } - - delete this->ptemp_rho_; - this->ptemp_rho_ = new Charge(); - this->ptemp_rho_->set_rhopw(this->pw_rho); - this->ptemp_rho_->allocate(PARAM.inp.nspin); - - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - delete[] this->pdLdphi_[is]; - delete[] this->pdEdphi_[is]; - delete[] this->pdirect_[is]; - delete[] this->precip_dir_[is]; - this->pdLdphi_[is] = new double[this->pw_rho->nrxx]; - this->pdEdphi_[is] = new double[this->pw_rho->nrxx]; - this->pdirect_[is] = new double[this->pw_rho->nrxx]; - this->precip_dir_[is] = new std::complex[pw_rho->npw]; - } - } - - this->pelec->init_scf(istep, ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm); - - Symmetry_rho srho; - for (int is = 0; is < PARAM.inp.nspin; is++) - { - srho.begin(is, this->chr, this->pw_rho, ucell.symm); - } - - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - if (PARAM.inp.init_chg != "file") - { - for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs) - { - // Here we initialize rho to be uniform, - // because the rho got by pot.init_pot -> Charge::atomic_rho may contain minus elements. - this->chr.rho[is][ibs] = this->nelec_[is] / this->pelec->omega; - this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]); - } - } - else - { - for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs) - { - this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]); - } - } - } - - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - this->pelec->eferm.set_efval(is, 0); - this->theta_[is] = 0.; - ModuleBase::GlobalFunc::ZEROS(this->pdLdphi_[is], this->pw_rho->nrxx); - ModuleBase::GlobalFunc::ZEROS(this->pdEdphi_[is], this->pw_rho->nrxx); - ModuleBase::GlobalFunc::ZEROS(this->pdirect_[is], this->pw_rho->nrxx); - } - if (PARAM.inp.nspin == 1) - { - this->theta_[0] = 0.2; - } - - ModuleBase::timer::tick("ESolver_OF", "before_opt"); -} - -void ESolver_OF_TDDFT::get_Hpsi(UnitCell& ucell, const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi) -{ - // update rho - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - for (int ir = 0; ir < pw_rho->nrxx; ++ir) - { - this->chr.rho[is][ir] = abs(psi_[is][ir])*abs(psi_[is][ir]); - } - } - - this->pelec->pot->update_from_charge(&this->chr, &ucell); // Hartree + XC + external - this->get_tf_potential(this->chr.rho,pw_rho ,this->pelec->pot->get_effective_v()); // TF potential - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - const double* vr_eff = this->pelec->pot->get_effective_v(is); - for (int ir = 0; ir < pw_rho->nrxx; ++ir) - { - Hpsi[is][ir] = vr_eff[ir]*psi_[is][ir]; - } - } - this->get_vw_potential_phi(psi_, pw_rho, Hpsi); -} - -void ESolver_OF_TDDFT::get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) -{ - const double c_tf_ - = 3.0 / 10.0 * std::pow(3 * std::pow(M_PI, 2.0), 2.0 / 3.0) - * 2; - if (PARAM.inp.nspin == 1) - { - for (int ir = 0; ir < pw_rho->nrxx; ++ir) - { - rpot(0, ir) += 5.0 / 3.0 * c_tf_ * std::pow(prho[0][ir], 2. / 3.); - } - } - else if (PARAM.inp.nspin == 2) - { - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - for (int ir = 0; ir < pw_rho->nrxx; ++ir) - { - rpot(is, ir) += 5.0 / 3.0 * c_tf_ * std::pow(2. * prho[is][ir], 2. / 3.); - } - } - } -} - -void ESolver_OF_TDDFT::get_vw_potential_phi(const std::complex* const* pphi, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi) -{ - std::complex** rLapPhi = new std::complex*[PARAM.inp.nspin]; - for (int is = 0; is < PARAM.inp.nspin; ++is) { - rLapPhi[is] = new std::complex[pw_rho->nrxx]; - } - std::complex** recipPhi = new std::complex*[PARAM.inp.nspin]; - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - recipPhi[is] = new std::complex[pw_rho->npw]; - - pw_rho->real2recip(pphi[is], recipPhi[is]); - for (int ik = 0; ik < pw_rho->npw; ++ik) - { - recipPhi[is][ik] *= pw_rho->gg[ik] * pw_rho->tpiba2; - } - pw_rho->recip2real(recipPhi[is], rLapPhi[is]); - for (int ik = 0; ik < pw_rho->npw; ++ik) - { - Hpsi[is][ik]+=rLapPhi[is][ik]; - } - } - - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - delete[] recipPhi[is]; - delete[] rLapPhi[is]; - } - delete[] recipPhi; - delete[] rLapPhi; -} - -void ESolver_OF_TDDFT::get_CD_potential(const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) -{ - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - //recipCurrent = new std::complex[pw_rho->npw]; - //delete[] recipCurrent; - } -} - -void ESolver_OF_TDDFT::propagate_psi(UnitCell& ucell, std::complex** pphi_, ModulePW::PW_Basis* pw_rho) -{ - ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); - - std::complex imag(0.0,1.0); - double dt=PARAM.inp.mdp.md_dt; - std::complex** K1 = new std::complex*[PARAM.inp.nspin]; - std::complex** K2 = new std::complex*[PARAM.inp.nspin]; - std::complex** K3 = new std::complex*[PARAM.inp.nspin]; - std::complex** K4 = new std::complex*[PARAM.inp.nspin]; - std::complex** psi1 = new std::complex*[PARAM.inp.nspin]; - std::complex** psi2 = new std::complex*[PARAM.inp.nspin]; - std::complex** psi3 = new std::complex*[PARAM.inp.nspin]; - for (int is = 0; is < PARAM.inp.nspin; ++is) { - K1[is] = new std::complex[pw_rho->nrxx]; - K2[is] = new std::complex[pw_rho->nrxx]; - K3[is] = new std::complex[pw_rho->nrxx]; - K4[is] = new std::complex[pw_rho->nrxx]; - psi1[is] = new std::complex[pw_rho->nrxx]; - psi2[is] = new std::complex[pw_rho->nrxx]; - psi3[is] = new std::complex[pw_rho->nrxx]; - } - get_Hpsi(ucell,pphi_,pw_rho,K1); - for (int is = 0; is < PARAM.inp.nspin; ++is){ - for (int ir = 0; ir < pw_rho->nrxx; ++ir) - { - K1[is][ir]=-1.0*K1[is][ir]*dt*imag; - psi1[is][ir]=pphi_[is][ir]+0.5*K1[is][ir]; - } - } - get_Hpsi(ucell,psi1,pw_rho,K2); - for (int is = 0; is < PARAM.inp.nspin; ++is){ - for (int ir = 0; ir < pw_rho->nrxx; ++ir) - { - K2[is][ir]=-1.0*K2[is][ir]*dt*imag; - psi2[is][ir]=pphi_[is][ir]+0.5*K2[is][ir]; - } - } - get_Hpsi(ucell,psi2,pw_rho,K3); - for (int is = 0; is < PARAM.inp.nspin; ++is){ - for (int ir = 0; ir < pw_rho->nrxx; ++ir) - { - K3[is][ir]=-1.0*K3[is][ir]*dt*imag; - psi3[is][ir]=pphi_[is][ir]+K3[is][ir]; - } - } - get_Hpsi(ucell,psi3,pw_rho,K4); - for (int is = 0; is < PARAM.inp.nspin; ++is){ - for (int ir = 0; ir < pw_rho->nrxx; ++ir) - { - K4[is][ir]=-1.0*K4[is][ir]*dt*imag; - pphi_[is][ir]+=1.0/6.0*(K1[is][ir]+2.0*K2[is][ir]+2.0*K3[is][ir]+K4[is][ir]); - } - } - - for (int is = 0; is < PARAM.inp.nspin; ++is) { - delete[] K1[is]; - delete[] K2[is]; - delete[] K3[is]; - delete[] K4[is]; - delete[] psi1[is]; - delete[] psi2[is]; - delete[] psi3[is]; - } - delete[] K1; - delete[] K2; - delete[] K3; - delete[] K4; - delete[] psi1; - delete[] psi2; - delete[] psi3; - ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); -} - } // namespace ModuleESolver diff --git a/source/source_esolver/esolver_of_tddft.h b/source/source_esolver/esolver_of_tddft.h index 2cf30a9c66..597a9d7f17 100644 --- a/source/source_esolver/esolver_of_tddft.h +++ b/source/source_esolver/esolver_of_tddft.h @@ -2,6 +2,7 @@ #define ESOLVER_OF_TDDFT_H #include "esolver_of.h" +#include "source_pw/module_ofdft/evolve_psi.h" namespace ModuleESolver { @@ -15,15 +16,7 @@ class ESolver_OF_TDDFT : public ESolver_OF protected: std::complex** pphi_td = nullptr; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). - psi::Psi>* psi_td = nullptr; - //std::complex** pdEdphi_phi_ = nullptr; // (dE/dphi)*phi - // ==================== main process of TDOFDFT ====================== - void before_opt(const int istep, UnitCell& ucell); - void get_Hpsi(UnitCell& ucell, const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi); - void get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); - void get_vw_potential_phi(const std::complex* const* pphi, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi); // -1/2 \nabla^2 \phi - void get_CD_potential(const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); - void propagate_psi(UnitCell& ucell, std::complex** pphi_, ModulePW::PW_Basis* pw_rho); + EVOLVE_PSI* evolve_psi=nullptr; }; } // namespace ModuleESolver diff --git a/source/source_pw/module_ofdft/CMakeLists.txt b/source/source_pw/module_ofdft/CMakeLists.txt index d34ac86492..967a673713 100644 --- a/source/source_pw/module_ofdft/CMakeLists.txt +++ b/source/source_pw/module_ofdft/CMakeLists.txt @@ -6,6 +6,7 @@ list(APPEND hamilt_ofdft_srcs kedf_lkt.cpp kedf_manager.cpp of_stress_pw.cpp + evolve_phi.cpp ) add_library( diff --git a/source/source_pw/module_ofdft/evolve_phi.cpp b/source/source_pw/module_ofdft/evolve_phi.cpp new file mode 100644 index 0000000000..646a82d04f --- /dev/null +++ b/source/source_pw/module_ofdft/evolve_phi.cpp @@ -0,0 +1,166 @@ +#include "evolve_phi.h" + +#include "source_io/module_parameter/parameter.h" +#include + +#include "source_base/parallel_reduce.h" + +void EVOLVE_PHI::get_Hpsi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi) +{ + // update rho + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + chr.rho[is][ir] = abs(psi_[is][ir])*abs(psi_[is][ir]); + } + } + + pelec->pot->update_from_charge(&chr, &ucell); // Hartree + XC + external + this->get_tf_potential(chr.rho,pw_rho ,pelec->pot->get_effective_v()); // TF potential + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + const double* vr_eff = pelec->pot->get_effective_v(is); + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + Hpsi[is][ir] = vr_eff[ir]*psi_[is][ir]; + } + } + this->get_vw_potential_phi(psi_, pw_rho, Hpsi); +} + +void EVOLVE_PHI::get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) +{ + if (PARAM.inp.nspin == 1) + { + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + rpot(0, ir) += 5.0 / 3.0 * this->c_tf_ * std::pow(prho[0][ir], 2. / 3.); + } + } + else if (PARAM.inp.nspin == 2) + { + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + rpot(is, ir) += 5.0 / 3.0 * this->c_tf_ * std::pow(2. * prho[is][ir], 2. / 3.); + } + } + } +} + +void EVOLVE_PHI::get_vw_potential_phi(const std::complex* const* pphi, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi) +{ + std::complex** rLapPhi = new std::complex*[PARAM.inp.nspin]; + for (int is = 0; is < PARAM.inp.nspin; ++is) { + rLapPhi[is] = new std::complex[pw_rho->nrxx]; + } + std::complex** recipPhi = new std::complex*[PARAM.inp.nspin]; + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + recipPhi[is] = new std::complex[pw_rho->npw]; + + pw_rho->real2recip(pphi[is], recipPhi[is]); + for (int ik = 0; ik < pw_rho->npw; ++ik) + { + recipPhi[is][ik] *= pw_rho->gg[ik] * pw_rho->tpiba2; + } + pw_rho->recip2real(recipPhi[is], rLapPhi[is]); + for (int ik = 0; ik < pw_rho->npw; ++ik) + { + Hpsi[is][ik]+=rLapPhi[is][ik]; + } + } + + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + delete[] recipPhi[is]; + delete[] rLapPhi[is]; + } + delete[] recipPhi; + delete[] rLapPhi; +} + +void EVOLVE_PHI::get_CD_potential(const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) +{ + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + //recipCurrent = new std::complex[pw_rho->npw]; + //delete[] recipCurrent; + } +} + +void EVOLVE_PHI::propagate_psi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::complex** pphi_, ModulePW::PW_Basis* pw_rho) +{ + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); + + std::complex imag(0.0,1.0); + double dt=PARAM.inp.mdp.md_dt; + std::complex** K1 = new std::complex*[PARAM.inp.nspin]; + std::complex** K2 = new std::complex*[PARAM.inp.nspin]; + std::complex** K3 = new std::complex*[PARAM.inp.nspin]; + std::complex** K4 = new std::complex*[PARAM.inp.nspin]; + std::complex** psi1 = new std::complex*[PARAM.inp.nspin]; + std::complex** psi2 = new std::complex*[PARAM.inp.nspin]; + std::complex** psi3 = new std::complex*[PARAM.inp.nspin]; + for (int is = 0; is < PARAM.inp.nspin; ++is) { + K1[is] = new std::complex[pw_rho->nrxx]; + K2[is] = new std::complex[pw_rho->nrxx]; + K3[is] = new std::complex[pw_rho->nrxx]; + K4[is] = new std::complex[pw_rho->nrxx]; + psi1[is] = new std::complex[pw_rho->nrxx]; + psi2[is] = new std::complex[pw_rho->nrxx]; + psi3[is] = new std::complex[pw_rho->nrxx]; + } + get_Hpsi(pelec,chr,ucell,pphi_,pw_rho,K1); + for (int is = 0; is < PARAM.inp.nspin; ++is){ + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + K1[is][ir]=-1.0*K1[is][ir]*dt*imag; + psi1[is][ir]=pphi_[is][ir]+0.5*K1[is][ir]; + } + } + get_Hpsi(pelec,chr,ucell,psi1,pw_rho,K2); + for (int is = 0; is < PARAM.inp.nspin; ++is){ + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + K2[is][ir]=-1.0*K2[is][ir]*dt*imag; + psi2[is][ir]=pphi_[is][ir]+0.5*K2[is][ir]; + } + } + get_Hpsi(pelec,chr,ucell,psi2,pw_rho,K3); + for (int is = 0; is < PARAM.inp.nspin; ++is){ + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + K3[is][ir]=-1.0*K3[is][ir]*dt*imag; + psi3[is][ir]=pphi_[is][ir]+K3[is][ir]; + } + } + get_Hpsi(pelec,chr,ucell,psi3,pw_rho,K4); + for (int is = 0; is < PARAM.inp.nspin; ++is){ + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + K4[is][ir]=-1.0*K4[is][ir]*dt*imag; + pphi_[is][ir]+=1.0/6.0*(K1[is][ir]+2.0*K2[is][ir]+2.0*K3[is][ir]+K4[is][ir]); + } + } + + for (int is = 0; is < PARAM.inp.nspin; ++is) { + delete[] K1[is]; + delete[] K2[is]; + delete[] K3[is]; + delete[] K4[is]; + delete[] psi1[is]; + delete[] psi2[is]; + delete[] psi3[is]; + } + delete[] K1; + delete[] K2; + delete[] K3; + delete[] K4; + delete[] psi1; + delete[] psi2; + delete[] psi3; + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); +} \ No newline at end of file diff --git a/source/source_pw/module_ofdft/evolve_phi.h b/source/source_pw/module_ofdft/evolve_phi.h new file mode 100644 index 0000000000..fb6737cd44 --- /dev/null +++ b/source/source_pw/module_ofdft/evolve_phi.h @@ -0,0 +1,40 @@ +#ifndef EVOLVE_PHI_H +#define EVOLVE_PHI_H +#include +#include + +#include "source_base/global_function.h" +#include "source_base/global_variable.h" +#include "source_base/matrix.h" +#include "source_base/timer.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_estate/elecstate.h" // electronic states +#include "source_estate/module_charge/charge.h" + +/** + * @brief TDOFDFT + * @author liyuanbo on 2025-09 + */ +class EVOLVE_PHI +{ + public: + EVOLVE_PHI() + { + } + ~EVOLVE_PHI() + { + } + void propagate_psi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::complex** pphi_, ModulePW::PW_Basis* pw_rho); + + private: + const double c_tf_ + = 3.0 / 10.0 * std::pow(3 * std::pow(M_PI, 2.0), 2.0 / 3.0) + * 2; // 10/3*(3*pi^2)^{2/3}, multiply by 2 to convert unit from Hartree to Ry, finally in Ry*Bohr^(-2) + + void get_Hpsi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi); + void get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); + void get_vw_potential_phi(const std::complex* const* pphi, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi); // -1/2 \nabla^2 \phi + void get_CD_potential(const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); + +}; +#endif \ No newline at end of file From d311eb37b2632a9a630e7294c4b6e7715a62a959 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sat, 27 Sep 2025 04:45:21 +0800 Subject: [PATCH 4/8] fix bug (confuse evolve_psi with evolve_phi) --- source/source_esolver/esolver_of_tddft.cpp | 6 +++--- source/source_esolver/esolver_of_tddft.h | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/source/source_esolver/esolver_of_tddft.cpp b/source/source_esolver/esolver_of_tddft.cpp index 7585d42ed2..11c72c178b 100644 --- a/source/source_esolver/esolver_of_tddft.cpp +++ b/source/source_esolver/esolver_of_tddft.cpp @@ -22,12 +22,12 @@ namespace ModuleESolver ESolver_OF_TDDFT::ESolver_OF_TDDFT() { this->classname = "ESolver_OF_TDDFT"; - this->evolve_psi=new EVOLVE_PSI(); + this->evolve_phi=new EVOLVE_PHI(); } ESolver_OF_TDDFT::~ESolver_OF_TDDFT() { - delete this->evolve_psi; + delete this->evolve_phi; for (int is = 0; is < PARAM.inp.nspin; ++is) { delete[] this->pphi_td[is]; @@ -102,7 +102,7 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) } else { - this->evolve_psi->propagate_psi(this->pelec, this->chr, ucell, this->pphi_td, this->pw_rho); + this->evolve_phi->propagate_psi(this->pelec, this->chr, ucell, this->pphi_td, this->pw_rho); for (int is = 0; is < PARAM.inp.nspin; ++is) { for (int ir = 0; ir < pw_rho->nrxx; ++ir) diff --git a/source/source_esolver/esolver_of_tddft.h b/source/source_esolver/esolver_of_tddft.h index 597a9d7f17..6794f80cad 100644 --- a/source/source_esolver/esolver_of_tddft.h +++ b/source/source_esolver/esolver_of_tddft.h @@ -2,7 +2,7 @@ #define ESOLVER_OF_TDDFT_H #include "esolver_of.h" -#include "source_pw/module_ofdft/evolve_psi.h" +#include "source_pw/module_ofdft/evolve_phi.h" namespace ModuleESolver { @@ -16,7 +16,7 @@ class ESolver_OF_TDDFT : public ESolver_OF protected: std::complex** pphi_td = nullptr; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). - EVOLVE_PSI* evolve_psi=nullptr; + EVOLVE_PHI* evolve_phi=nullptr; }; } // namespace ModuleESolver From 262c4e5cb9c2b8898b1f5f9a6bf4d5165b1c00f0 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sat, 27 Sep 2025 13:10:04 +0800 Subject: [PATCH 5/8] change evolve_phi with evolve_ofdft; change to vector --- source/Makefile.Objects | 2 +- source/source_esolver/esolver_of_tddft.cpp | 19 +++--- source/source_esolver/esolver_of_tddft.h | 6 +- source/source_pw/module_ofdft/CMakeLists.txt | 2 +- .../{evolve_phi.cpp => evolve_ofdft.cpp} | 58 ++++++------------- .../{evolve_phi.h => evolve_ofdft.h} | 18 +++--- 6 files changed, 40 insertions(+), 65 deletions(-) rename source/source_pw/module_ofdft/{evolve_phi.cpp => evolve_ofdft.cpp} (63%) rename source/source_pw/module_ofdft/{evolve_phi.h => evolve_ofdft.h} (56%) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 460a474ea5..f51c21d625 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -361,7 +361,7 @@ OBJS_HAMILT_OF=kedf_tf.o\ kedf_xwm.o\ kedf_lkt.o\ kedf_manager.o\ - evolve_phi.o\ + evolve_ofdft.o\ OBJS_HAMILT_LCAO=hamilt_lcao.o\ operator_lcao.o\ diff --git a/source/source_esolver/esolver_of_tddft.cpp b/source/source_esolver/esolver_of_tddft.cpp index 11c72c178b..ed85b67b1e 100644 --- a/source/source_esolver/esolver_of_tddft.cpp +++ b/source/source_esolver/esolver_of_tddft.cpp @@ -22,17 +22,12 @@ namespace ModuleESolver ESolver_OF_TDDFT::ESolver_OF_TDDFT() { this->classname = "ESolver_OF_TDDFT"; - this->evolve_phi=new EVOLVE_PHI(); + this->evolve_ofdft=new Evolve_OFDFT(); } ESolver_OF_TDDFT::~ESolver_OF_TDDFT() { - delete this->evolve_phi; - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - delete[] this->pphi_td[is]; - } - delete[] this->pphi_td; + delete this->evolve_ofdft; } @@ -52,11 +47,11 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) if (istep==0) { - this->pphi_td= new std::complex*[PARAM.inp.nspin]; + this->pphi_td.resize(PARAM.inp.nspin); for (int is = 0; is < PARAM.inp.nspin; ++is) { - this->pphi_td[is]= new std::complex[pw_rho->nrxx]; + this->pphi_td[is].resize(this->pw_rho->nrxx); } } @@ -94,7 +89,7 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) for (int is = 0; is < PARAM.inp.nspin; ++is) { - for (int ir = 0; ir < pw_rho->nrxx; ++ir) + for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) { pphi_td[is][ir]=pphi_[is][ir]; } @@ -102,10 +97,10 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) } else { - this->evolve_phi->propagate_psi(this->pelec, this->chr, ucell, this->pphi_td, this->pw_rho); + this->evolve_ofdft->propagate_psi(this->pelec, this->chr, ucell, this->pphi_td, this->pw_rho); for (int is = 0; is < PARAM.inp.nspin; ++is) { - for (int ir = 0; ir < pw_rho->nrxx; ++ir) + for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) { pphi_[is][ir]=abs(pphi_td[is][ir]); } diff --git a/source/source_esolver/esolver_of_tddft.h b/source/source_esolver/esolver_of_tddft.h index 6794f80cad..85570573f2 100644 --- a/source/source_esolver/esolver_of_tddft.h +++ b/source/source_esolver/esolver_of_tddft.h @@ -2,7 +2,7 @@ #define ESOLVER_OF_TDDFT_H #include "esolver_of.h" -#include "source_pw/module_ofdft/evolve_phi.h" +#include "source_pw/module_ofdft/evolve_ofdft.h" namespace ModuleESolver { @@ -15,8 +15,8 @@ class ESolver_OF_TDDFT : public ESolver_OF virtual void runner(UnitCell& ucell, const int istep) override; protected: - std::complex** pphi_td = nullptr; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). - EVOLVE_PHI* evolve_phi=nullptr; + std::vector>> pphi_td; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). + Evolve_OFDFT* evolve_ofdft=nullptr; }; } // namespace ModuleESolver diff --git a/source/source_pw/module_ofdft/CMakeLists.txt b/source/source_pw/module_ofdft/CMakeLists.txt index 967a673713..7465b47a98 100644 --- a/source/source_pw/module_ofdft/CMakeLists.txt +++ b/source/source_pw/module_ofdft/CMakeLists.txt @@ -6,7 +6,7 @@ list(APPEND hamilt_ofdft_srcs kedf_lkt.cpp kedf_manager.cpp of_stress_pw.cpp - evolve_phi.cpp + evolve_ofdft.cpp ) add_library( diff --git a/source/source_pw/module_ofdft/evolve_phi.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp similarity index 63% rename from source/source_pw/module_ofdft/evolve_phi.cpp rename to source/source_pw/module_ofdft/evolve_ofdft.cpp index 646a82d04f..1f4d07658f 100644 --- a/source/source_pw/module_ofdft/evolve_phi.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -1,11 +1,11 @@ -#include "evolve_phi.h" +#include "evolve_ofdft.h" #include "source_io/module_parameter/parameter.h" #include #include "source_base/parallel_reduce.h" -void EVOLVE_PHI::get_Hpsi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi) +void Evolve_OFDFT::get_Hpsi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::vector>> psi_, ModulePW::PW_Basis* pw_rho, std::vector>> Hpsi) { // update rho for (int is = 0; is < PARAM.inp.nspin; ++is) @@ -29,7 +29,7 @@ void EVOLVE_PHI::get_Hpsi(elecstate::ElecState* pelec, const Charge& chr, UnitCe this->get_vw_potential_phi(psi_, pw_rho, Hpsi); } -void EVOLVE_PHI::get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) +void Evolve_OFDFT::get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) { if (PARAM.inp.nspin == 1) { @@ -50,18 +50,22 @@ void EVOLVE_PHI::get_tf_potential(const double* const* prho, ModulePW::PW_Basis* } } -void EVOLVE_PHI::get_vw_potential_phi(const std::complex* const* pphi, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi) +void Evolve_OFDFT::get_vw_potential_phi(std::vector>> pphi, ModulePW::PW_Basis* pw_rho, std::vector>> Hpsi) { std::complex** rLapPhi = new std::complex*[PARAM.inp.nspin]; for (int is = 0; is < PARAM.inp.nspin; ++is) { rLapPhi[is] = new std::complex[pw_rho->nrxx]; + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + rLapPhi[is][ir]=pphi[is][ir]; + } } std::complex** recipPhi = new std::complex*[PARAM.inp.nspin]; for (int is = 0; is < PARAM.inp.nspin; ++is) { recipPhi[is] = new std::complex[pw_rho->npw]; - pw_rho->real2recip(pphi[is], recipPhi[is]); + pw_rho->real2recip(rLapPhi[is], recipPhi[is]); for (int ik = 0; ik < pw_rho->npw; ++ik) { recipPhi[is][ik] *= pw_rho->gg[ik] * pw_rho->tpiba2; @@ -82,7 +86,7 @@ void EVOLVE_PHI::get_vw_potential_phi(const std::complex* const* pphi, M delete[] rLapPhi; } -void EVOLVE_PHI::get_CD_potential(const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) +void Evolve_OFDFT::get_CD_potential(std::vector>> psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) { for (int is = 0; is < PARAM.inp.nspin; ++is) { @@ -91,28 +95,20 @@ void EVOLVE_PHI::get_CD_potential(const std::complex* const* psi_, Modul } } -void EVOLVE_PHI::propagate_psi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::complex** pphi_, ModulePW::PW_Basis* pw_rho) +void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::vector>> pphi_, ModulePW::PW_Basis* pw_rho) { ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); std::complex imag(0.0,1.0); double dt=PARAM.inp.mdp.md_dt; - std::complex** K1 = new std::complex*[PARAM.inp.nspin]; - std::complex** K2 = new std::complex*[PARAM.inp.nspin]; - std::complex** K3 = new std::complex*[PARAM.inp.nspin]; - std::complex** K4 = new std::complex*[PARAM.inp.nspin]; - std::complex** psi1 = new std::complex*[PARAM.inp.nspin]; - std::complex** psi2 = new std::complex*[PARAM.inp.nspin]; - std::complex** psi3 = new std::complex*[PARAM.inp.nspin]; - for (int is = 0; is < PARAM.inp.nspin; ++is) { - K1[is] = new std::complex[pw_rho->nrxx]; - K2[is] = new std::complex[pw_rho->nrxx]; - K3[is] = new std::complex[pw_rho->nrxx]; - K4[is] = new std::complex[pw_rho->nrxx]; - psi1[is] = new std::complex[pw_rho->nrxx]; - psi2[is] = new std::complex[pw_rho->nrxx]; - psi3[is] = new std::complex[pw_rho->nrxx]; - } + std::vector>> K1(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); + std::vector>> K2(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); + std::vector>> K3(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); + std::vector>> K4(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); + std::vector>> psi1(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); + std::vector>> psi2(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); + std::vector>> psi3(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); + get_Hpsi(pelec,chr,ucell,pphi_,pw_rho,K1); for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) @@ -146,21 +142,5 @@ void EVOLVE_PHI::propagate_psi(elecstate::ElecState* pelec, const Charge& chr, U } } - for (int is = 0; is < PARAM.inp.nspin; ++is) { - delete[] K1[is]; - delete[] K2[is]; - delete[] K3[is]; - delete[] K4[is]; - delete[] psi1[is]; - delete[] psi2[is]; - delete[] psi3[is]; - } - delete[] K1; - delete[] K2; - delete[] K3; - delete[] K4; - delete[] psi1; - delete[] psi2; - delete[] psi3; ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); } \ No newline at end of file diff --git a/source/source_pw/module_ofdft/evolve_phi.h b/source/source_pw/module_ofdft/evolve_ofdft.h similarity index 56% rename from source/source_pw/module_ofdft/evolve_phi.h rename to source/source_pw/module_ofdft/evolve_ofdft.h index fb6737cd44..3ba9427643 100644 --- a/source/source_pw/module_ofdft/evolve_phi.h +++ b/source/source_pw/module_ofdft/evolve_ofdft.h @@ -1,5 +1,5 @@ -#ifndef EVOLVE_PHI_H -#define EVOLVE_PHI_H +#ifndef Evolve_OFDFT_H +#define Evolve_OFDFT_H #include #include @@ -15,26 +15,26 @@ * @brief TDOFDFT * @author liyuanbo on 2025-09 */ -class EVOLVE_PHI +class Evolve_OFDFT { public: - EVOLVE_PHI() + Evolve_OFDFT() { } - ~EVOLVE_PHI() + ~Evolve_OFDFT() { } - void propagate_psi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::complex** pphi_, ModulePW::PW_Basis* pw_rho); + void propagate_psi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::vector>> pphi_, ModulePW::PW_Basis* pw_rho); private: const double c_tf_ = 3.0 / 10.0 * std::pow(3 * std::pow(M_PI, 2.0), 2.0 / 3.0) * 2; // 10/3*(3*pi^2)^{2/3}, multiply by 2 to convert unit from Hartree to Ry, finally in Ry*Bohr^(-2) - void get_Hpsi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi); + void get_Hpsi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::vector>> psi_, ModulePW::PW_Basis* pw_rho, std::vector>> Hpsi); void get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); - void get_vw_potential_phi(const std::complex* const* pphi, ModulePW::PW_Basis* pw_rho, std::complex** Hpsi); // -1/2 \nabla^2 \phi - void get_CD_potential(const std::complex* const* psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); + void get_vw_potential_phi(std::vector>> pphi, ModulePW::PW_Basis* pw_rho, std::vector>> Hpsi); // -1/2 \nabla^2 \phi + void get_CD_potential(std::vector>> psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); }; #endif \ No newline at end of file From d38224d10d17423e021b43d8659ae2d9270d2bcf Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sun, 28 Sep 2025 00:47:04 +0800 Subject: [PATCH 6/8] details --- source/source_esolver/esolver_of_tddft.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/source/source_esolver/esolver_of_tddft.cpp b/source/source_esolver/esolver_of_tddft.cpp index ed85b67b1e..15886de7b8 100644 --- a/source/source_esolver/esolver_of_tddft.cpp +++ b/source/source_esolver/esolver_of_tddft.cpp @@ -77,8 +77,6 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) // find the optimization direction and step lenghth theta according to the potential this->optimize(ucell); - std::cout<<"optimize------"<update_rho(); @@ -102,7 +100,7 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) { for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) { - pphi_[is][ir]=abs(pphi_td[is][ir]); + pphi_[is][ir]=std::abs(pphi_td[is][ir]); } } conv_esolver=true; From cde47d38bad32a909af675f44070e1980ff884ed Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Mon, 29 Sep 2025 15:04:39 +0800 Subject: [PATCH 7/8] code optimization --- source/source_esolver/esolver_of_tddft.cpp | 19 +-- source/source_esolver/esolver_of_tddft.h | 2 +- .../source_pw/module_ofdft/evolve_ofdft.cpp | 112 +++++++++++++----- source/source_pw/module_ofdft/evolve_ofdft.h | 24 +++- 4 files changed, 111 insertions(+), 46 deletions(-) diff --git a/source/source_esolver/esolver_of_tddft.cpp b/source/source_esolver/esolver_of_tddft.cpp index 15886de7b8..116550e7fa 100644 --- a/source/source_esolver/esolver_of_tddft.cpp +++ b/source/source_esolver/esolver_of_tddft.cpp @@ -47,12 +47,7 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) if (istep==0) { - this->pphi_td.resize(PARAM.inp.nspin); - - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - this->pphi_td[is].resize(this->pw_rho->nrxx); - } + this->phi_td.resize(PARAM.inp.nspin*this->pw_rho->nrxx); } if ((istep<1) && PARAM.inp.init_chg != "file") @@ -85,22 +80,28 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) ESolver_FP::iter_finish(ucell, istep, this->iter_, conv_esolver); } +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif for (int is = 0; is < PARAM.inp.nspin; ++is) { for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) { - pphi_td[is][ir]=pphi_[is][ir]; + phi_td[is*this->pw_rho->nrxx+ir]=pphi_[is][ir]; } } } else { - this->evolve_ofdft->propagate_psi(this->pelec, this->chr, ucell, this->pphi_td, this->pw_rho); + this->evolve_ofdft->propagate_psi(this->pelec, this->chr, ucell, this->phi_td, this->pw_rho); +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif for (int is = 0; is < PARAM.inp.nspin; ++is) { for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) { - pphi_[is][ir]=std::abs(pphi_td[is][ir]); + pphi_[is][ir]=std::abs(phi_td[is*this->pw_rho->nrxx+ir]); } } conv_esolver=true; diff --git a/source/source_esolver/esolver_of_tddft.h b/source/source_esolver/esolver_of_tddft.h index 85570573f2..99a4238c0f 100644 --- a/source/source_esolver/esolver_of_tddft.h +++ b/source/source_esolver/esolver_of_tddft.h @@ -15,7 +15,7 @@ class ESolver_OF_TDDFT : public ESolver_OF virtual void runner(UnitCell& ucell, const int istep) override; protected: - std::vector>> pphi_td; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). + std::vector> phi_td; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). Evolve_OFDFT* evolve_ofdft=nullptr; }; } // namespace ModuleESolver diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index 1f4d07658f..4c158998cc 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -5,34 +5,49 @@ #include "source_base/parallel_reduce.h" -void Evolve_OFDFT::get_Hpsi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::vector>> psi_, ModulePW::PW_Basis* pw_rho, std::vector>> Hpsi) +void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec, + const Charge& chr, + UnitCell& ucell, + std::vector> psi_, + ModulePW::PW_Basis* pw_rho, + std::vector> Hpsi) { // update rho +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif for (int is = 0; is < PARAM.inp.nspin; ++is) { for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - chr.rho[is][ir] = abs(psi_[is][ir])*abs(psi_[is][ir]); + chr.rho[is][ir] = abs(psi_[is * pw_rho->nrxx + ir])*abs(psi_[is * pw_rho->nrxx + ir]); } } pelec->pot->update_from_charge(&chr, &ucell); // Hartree + XC + external - this->get_tf_potential(chr.rho,pw_rho ,pelec->pot->get_effective_v()); // TF potential + this->cal_tf_potential(chr.rho,pw_rho ,pelec->pot->get_effective_v()); // TF potential + +#ifdef _OPENMP +#pragma omp parallel for +#endif for (int is = 0; is < PARAM.inp.nspin; ++is) { const double* vr_eff = pelec->pot->get_effective_v(is); for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - Hpsi[is][ir] = vr_eff[ir]*psi_[is][ir]; + Hpsi[is * pw_rho->nrxx + ir] = vr_eff[ir]*psi_[is * pw_rho->nrxx + ir]; } } - this->get_vw_potential_phi(psi_, pw_rho, Hpsi); + this->cal_vw_potential_phi(psi_, pw_rho, Hpsi); } -void Evolve_OFDFT::get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) +void Evolve_OFDFT::cal_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) { if (PARAM.inp.nspin == 1) { +#ifdef _OPENMP +#pragma omp parallel for +#endif for (int ir = 0; ir < pw_rho->nrxx; ++ir) { rpot(0, ir) += 5.0 / 3.0 * this->c_tf_ * std::pow(prho[0][ir], 2. / 3.); @@ -40,6 +55,9 @@ void Evolve_OFDFT::get_tf_potential(const double* const* prho, ModulePW::PW_Basi } else if (PARAM.inp.nspin == 2) { +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif for (int is = 0; is < PARAM.inp.nspin; ++is) { for (int ir = 0; ir < pw_rho->nrxx; ++ir) @@ -50,17 +68,26 @@ void Evolve_OFDFT::get_tf_potential(const double* const* prho, ModulePW::PW_Basi } } -void Evolve_OFDFT::get_vw_potential_phi(std::vector>> pphi, ModulePW::PW_Basis* pw_rho, std::vector>> Hpsi) +void Evolve_OFDFT::cal_vw_potential_phi(std::vector> pphi, + ModulePW::PW_Basis* pw_rho, + std::vector> Hpsi) { std::complex** rLapPhi = new std::complex*[PARAM.inp.nspin]; +#ifdef _OPENMP +#pragma omp parallel for +#endif for (int is = 0; is < PARAM.inp.nspin; ++is) { rLapPhi[is] = new std::complex[pw_rho->nrxx]; for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - rLapPhi[is][ir]=pphi[is][ir]; + rLapPhi[is][ir]=pphi[is * pw_rho->nrxx + ir]; } } std::complex** recipPhi = new std::complex*[PARAM.inp.nspin]; + +#ifdef _OPENMP +#pragma omp parallel for +#endif for (int is = 0; is < PARAM.inp.nspin; ++is) { recipPhi[is] = new std::complex[pw_rho->npw]; @@ -71,12 +98,15 @@ void Evolve_OFDFT::get_vw_potential_phi(std::vectorgg[ik] * pw_rho->tpiba2; } pw_rho->recip2real(recipPhi[is], rLapPhi[is]); - for (int ik = 0; ik < pw_rho->npw; ++ik) + for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - Hpsi[is][ik]+=rLapPhi[is][ik]; + Hpsi[is * pw_rho->nrxx + ir]+=rLapPhi[is][ir]; } } +#ifdef _OPENMP +#pragma omp parallel for +#endif for (int is = 0; is < PARAM.inp.nspin; ++is) { delete[] recipPhi[is]; @@ -86,7 +116,9 @@ void Evolve_OFDFT::get_vw_potential_phi(std::vector>> psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) +void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, + ModulePW::PW_Basis* pw_rho, + ModuleBase::matrix& rpot) { for (int is = 0; is < PARAM.inp.nspin; ++is) { @@ -95,50 +127,68 @@ void Evolve_OFDFT::get_CD_potential(std::vector } } -void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::vector>> pphi_, ModulePW::PW_Basis* pw_rho) +void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, + const Charge& chr, UnitCell& ucell, + std::vector> pphi_, + ModulePW::PW_Basis* pw_rho) { ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); std::complex imag(0.0,1.0); double dt=PARAM.inp.mdp.md_dt; - std::vector>> K1(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); - std::vector>> K2(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); - std::vector>> K3(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); - std::vector>> K4(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); - std::vector>> psi1(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); - std::vector>> psi2(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); - std::vector>> psi3(PARAM.inp.nspin,std::vector>(pw_rho->nrxx)); + const int nspin = PARAM.inp.nspin; + const int nrxx = pw_rho->nrxx; + const int total_size = nspin * nrxx; + std::vector> K1(total_size); + std::vector> K2(total_size); + std::vector> K3(total_size); + std::vector> K4(total_size); + std::vector> psi1(total_size); + std::vector> psi2(total_size); + std::vector> psi3(total_size); - get_Hpsi(pelec,chr,ucell,pphi_,pw_rho,K1); + cal_Hpsi(pelec,chr,ucell,pphi_,pw_rho,K1); +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K1[is][ir]=-1.0*K1[is][ir]*dt*imag; - psi1[is][ir]=pphi_[is][ir]+0.5*K1[is][ir]; + K1[is * nrxx + ir]=-1.0*K1[is * nrxx + ir]*dt*imag; + psi1[is * nrxx + ir]=pphi_[is * nrxx + ir]+0.5*K1[is * nrxx + ir]; } } - get_Hpsi(pelec,chr,ucell,psi1,pw_rho,K2); + cal_Hpsi(pelec,chr,ucell,psi1,pw_rho,K2); +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K2[is][ir]=-1.0*K2[is][ir]*dt*imag; - psi2[is][ir]=pphi_[is][ir]+0.5*K2[is][ir]; + K2[is * nrxx + ir]=-1.0*K2[is * nrxx + ir]*dt*imag; + psi2[is * nrxx + ir]=pphi_[is * nrxx + ir]+0.5*K2[is * nrxx + ir]; } } - get_Hpsi(pelec,chr,ucell,psi2,pw_rho,K3); + cal_Hpsi(pelec,chr,ucell,psi2,pw_rho,K3); +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K3[is][ir]=-1.0*K3[is][ir]*dt*imag; - psi3[is][ir]=pphi_[is][ir]+K3[is][ir]; + K3[is * nrxx + ir]=-1.0*K3[is * nrxx + ir]*dt*imag; + psi3[is * nrxx + ir]=pphi_[is * nrxx + ir]+K3[is * nrxx + ir]; } } - get_Hpsi(pelec,chr,ucell,psi3,pw_rho,K4); + cal_Hpsi(pelec,chr,ucell,psi3,pw_rho,K4); +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K4[is][ir]=-1.0*K4[is][ir]*dt*imag; - pphi_[is][ir]+=1.0/6.0*(K1[is][ir]+2.0*K2[is][ir]+2.0*K3[is][ir]+K4[is][ir]); + K4[is * nrxx + ir]=-1.0*K4[is * nrxx + ir]*dt*imag; + pphi_[is * nrxx + ir]+=1.0/6.0*(K1[is * nrxx + ir]+2.0*K2[is * nrxx + ir]+2.0*K3[is * nrxx + ir]+K4[is * nrxx + ir]); } } diff --git a/source/source_pw/module_ofdft/evolve_ofdft.h b/source/source_pw/module_ofdft/evolve_ofdft.h index 3ba9427643..d3cd6d358b 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.h +++ b/source/source_pw/module_ofdft/evolve_ofdft.h @@ -24,17 +24,31 @@ class Evolve_OFDFT ~Evolve_OFDFT() { } - void propagate_psi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::vector>> pphi_, ModulePW::PW_Basis* pw_rho); + void propagate_psi(elecstate::ElecState* pelec, + const Charge& chr, UnitCell& ucell, + std::vector> pphi_, + ModulePW::PW_Basis* pw_rho); private: const double c_tf_ = 3.0 / 10.0 * std::pow(3 * std::pow(M_PI, 2.0), 2.0 / 3.0) * 2; // 10/3*(3*pi^2)^{2/3}, multiply by 2 to convert unit from Hartree to Ry, finally in Ry*Bohr^(-2) - void get_Hpsi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::vector>> psi_, ModulePW::PW_Basis* pw_rho, std::vector>> Hpsi); - void get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); - void get_vw_potential_phi(std::vector>> pphi, ModulePW::PW_Basis* pw_rho, std::vector>> Hpsi); // -1/2 \nabla^2 \phi - void get_CD_potential(std::vector>> psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); + void cal_Hpsi(elecstate::ElecState* pelec, + const Charge& chr, + UnitCell& ucell, + std::vector> psi_, + ModulePW::PW_Basis* pw_rho, + std::vector> Hpsi); + void cal_tf_potential(const double* const* prho, + ModulePW::PW_Basis* pw_rho, + ModuleBase::matrix& rpot); + void cal_vw_potential_phi(std::vector> pphi, + ModulePW::PW_Basis* pw_rho, + std::vector> Hpsi); // -1/2 \nabla^2 \phi + void cal_CD_potential(std::vector> psi_, + ModulePW::PW_Basis* pw_rho, + ModuleBase::matrix& rpot); }; #endif \ No newline at end of file From a35a7ae5d7478d5b0b572ade669d82a776fc77bf Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Tue, 30 Sep 2025 03:47:54 +0800 Subject: [PATCH 8/8] add situation of reading charge file for tdofdft --- source/source_esolver/esolver_of_tddft.cpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/source/source_esolver/esolver_of_tddft.cpp b/source/source_esolver/esolver_of_tddft.cpp index 116550e7fa..a978a1267c 100644 --- a/source/source_esolver/esolver_of_tddft.cpp +++ b/source/source_esolver/esolver_of_tddft.cpp @@ -91,6 +91,20 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) } } } + else if ((istep<1) && PARAM.inp.init_chg == "file") + { +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) + { + phi_td[is*this->pw_rho->nrxx+ir]=pphi_[is][ir]; + } + } + conv_esolver=true; + } else { this->evolve_ofdft->propagate_psi(this->pelec, this->chr, ucell, this->phi_td, this->pw_rho);