From ccb8a5d2ac62c1e94fd334d6a473c58ae972a395 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Fri, 26 Sep 2025 18:06:46 +0800 Subject: [PATCH 01/25] 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 02/25] 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 03/25] 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 04/25] 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 05/25] 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 06/25] 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 07/25] 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 08/25] 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); From ef291598ee5863f20a77f692d0047b8efdd2cdc1 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sat, 11 Oct 2025 21:18:35 +0800 Subject: [PATCH 09/25] add CD Potential --- .../source_pw/module_ofdft/evolve_ofdft.cpp | 76 ++++++++++++++++++- 1 file changed, 74 insertions(+), 2 deletions(-) diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index 4c158998cc..a31768e9e8 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -120,11 +120,83 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot) { + std::complex imag(0.0,1.0); + std::complex** recipPhi = new std::complex*[PARAM.inp.nspin]; + std::complex** rPhi = new std::complex*[PARAM.inp.nspin]; +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int is = 0; is < PARAM.inp.nspin; ++is) { + rPhi[is] = new std::complex[pw_rho->nrxx]; + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + rPhi[is][ir]=psi_[is * pw_rho->nrxx + ir]; + } + } + + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + std::complex *recipCurrent_x=new std::complex[pw_rho->npw]; + std::complex *recipCurrent_y=new std::complex[pw_rho->npw]; + std::complex *recipCurrent_z=new std::complex[pw_rho->npw]; + std::complex *recipCDPotential=new std::complex[pw_rho->npw]; + std::complex *rCurrent_x=new std::complex[pw_rho->nrxx]; + std::complex *rCurrent_y=new std::complex[pw_rho->nrxx]; + std::complex *rCurrent_z=new std::complex[pw_rho->nrxx]; + std::complex *kF_r=new std::complex[pw_rho->nrxx]; + std::complex *rCDPotential=new std::complex[pw_rho->nrxx]; + recipPhi[is] = new std::complex[pw_rho->npw]; + + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + kF_r[ir]=std::pow(3*std::pow(ModuleBase::PI*std::abs(rPhi[is][ir]),2),1/3); + } + + pw_rho->real2recip(rPhi[is], recipPhi[is]); + for (int ik = 0; ik < pw_rho->npw; ++ik) + { + recipCurrent_x[ik]=imag*pw_rho->gcar[ik].x*recipPhi[is][ik]* pw_rho->tpiba; + recipCurrent_y[ik]=imag*pw_rho->gcar[ik].y*recipPhi[is][ik]* pw_rho->tpiba; + recipCurrent_z[ik]=imag*pw_rho->gcar[ik].z*recipPhi[is][ik]* pw_rho->tpiba; + } + pw_rho->recip2real(recipCurrent_x,rCurrent_x); + pw_rho->recip2real(recipCurrent_y,rCurrent_y); + pw_rho->recip2real(recipCurrent_z,rCurrent_z); + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + rCurrent_x[ir]=std::imag(rCurrent_x[ir]*std::conj(rPhi[is][ir])); + rCurrent_y[ir]=std::imag(rCurrent_y[ir]*std::conj(rPhi[is][ir])); + rCurrent_z[ir]=std::imag(rCurrent_z[ir]*std::conj(rPhi[is][ir])); + } + pw_rho->real2recip(rCurrent_x,recipCurrent_x); + pw_rho->real2recip(rCurrent_y,recipCurrent_y); + pw_rho->real2recip(rCurrent_z,recipCurrent_z); + for (int ik = 0; ik < pw_rho->npw; ++ik) + { + recipCDPotential[ik]=recipCurrent_x[ik]*pw_rho->gcar[ik].x+recipCurrent_y[ik]*pw_rho->gcar[ik].y+recipCurrent_z[ik]*pw_rho->gcar[ik].z; + recipCDPotential[ik]*=imag/pw_rho->gg[ik]; + } + pw_rho->recip2real(recipCDPotential,rCDPotential); + + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + rpot(0, ir) -= std::real(2*rCDPotential[ir]*std::pow(ModuleBase::PI,3) / (2*std::pow(kF_r[ir],2))); + } + delete[] recipCurrent_x; + delete[] recipCurrent_y; + delete[] recipCurrent_z; + delete[] rCurrent_x; + delete[] rCurrent_y; + delete[] rCurrent_z; + } + for (int is = 0; is < PARAM.inp.nspin; ++is) { - //recipCurrent = new std::complex[pw_rho->npw]; - //delete[] recipCurrent; + delete[] recipPhi[is]; + delete[] rPhi[is]; } + delete[] recipPhi; + delete[] rPhi; } void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, From 2eeb5e9dce3775e42aa040a44aa328bb7a31a260 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sun, 12 Oct 2025 07:10:08 +0800 Subject: [PATCH 10/25] add INPUT parameters for tdofdft --- source/source_esolver/esolver_of.cpp | 16 +++++++++------- .../source_io/module_parameter/input_parameter.h | 5 +++++ source/source_io/read_input.cpp | 1 + source/source_io/read_input.h | 2 ++ source/source_io/read_input_item_tddft.cpp | 16 ++++++++++++++++ source/source_io/test/read_input_ptest.cpp | 2 ++ source/source_io/test/support/INPUT | 4 ++++ source/source_pw/module_ofdft/evolve_ofdft.cpp | 11 ++++++++--- source/source_pw/module_ofdft/evolve_ofdft.h | 3 ++- 9 files changed, 49 insertions(+), 11 deletions(-) diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index 45f26babcb..1b0fed366f 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -25,7 +25,10 @@ ESolver_OF::ESolver_OF() ESolver_OF::~ESolver_OF() { - delete psi_; + //**************************************************** + // do not add any codes in this deconstructor funcion + //**************************************************** + delete psi_; delete[] this->pphi_; for (int i = 0; i < PARAM.inp.nspin; ++i) @@ -492,18 +495,17 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_eso this->kedf_manager_->get_energy_density(this->chr.rho, this->pphi_, this->pw_rho, this->chr.kin_r); } - //------------------------------------------------------------------ - // 2) call after_scf() of ESolver_FP - //------------------------------------------------------------------ - ESolver_FP::after_scf(ucell, istep, conv_esolver); - - // should not be here? mohan note 2025-03-03 for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) { this->chr.rho_save[0][ir] = this->chr.rho[0][ir]; } + //------------------------------------------------------------------ + // 2) call after_scf() of ESolver_FP + //------------------------------------------------------------------ + ESolver_FP::after_scf(ucell, istep, conv_esolver); + #ifdef __MLALGO //------------------------------------------------------------------ // Generate data if needed diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 15ab2d3a3c..c917e97df1 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -686,5 +686,10 @@ struct Input_para // src/gga_c_pbe.c std::vector xc_corr_ext = { 130, 0.06672455060314922, 0.031090690869654895034, 1.00000}; + + // ============== #Parameters (24.td-ofdft) =========================== + bool of_cd = false; ///< add CD potential or not https://doi.org/10.1103/PhysRevB.98.144302 + double of_mCD_alpha = 1.0; /// parameter of modified CD Potential + }; #endif diff --git a/source/source_io/read_input.cpp b/source/source_io/read_input.cpp index c59c31d8c6..14f25d115d 100644 --- a/source/source_io/read_input.cpp +++ b/source/source_io/read_input.cpp @@ -101,6 +101,7 @@ ReadInput::ReadInput(const int& rank) this->item_sdft(); this->item_deepks(); this->item_rt_tddft(); + this->item_tdofdft(); this->item_lr_tddft(); this->item_output(); this->item_postprocess(); diff --git a/source/source_io/read_input.h b/source/source_io/read_input.h index 8f1f2e0ec0..92bda04bdb 100644 --- a/source/source_io/read_input.h +++ b/source/source_io/read_input.h @@ -108,6 +108,8 @@ class ReadInput // items for real time tddft void item_rt_tddft(); // items for linear response tddft + void item_tdofdft(); + // items for td-ofdft void item_lr_tddft(); // items for output void item_output(); diff --git a/source/source_io/read_input_item_tddft.cpp b/source/source_io/read_input_item_tddft.cpp index 89751fd224..66aa47164e 100644 --- a/source/source_io/read_input_item_tddft.cpp +++ b/source/source_io/read_input_item_tddft.cpp @@ -305,6 +305,22 @@ void ReadInput::item_rt_tddft() } +} +void ReadInput::item_tdofdft() +{ + // TD-OFDFT + { + Input_Item item("of_cd"); + item.annotation = "add CD Potential or not"; + read_sync_bool(input.of_cd); + this->add_item(item); + } + { + Input_Item item("of_mcd_alpha"); + item.annotation = "parameter of modified CD Potential"; + read_sync_double(input.of_mCD_alpha); + this->add_item(item); + } } void ReadInput::item_lr_tddft() { diff --git a/source/source_io/test/read_input_ptest.cpp b/source/source_io/test/read_input_ptest.cpp index 7e0287d90b..529aa88656 100644 --- a/source/source_io/test/read_input_ptest.cpp +++ b/source/source_io/test/read_input_ptest.cpp @@ -364,6 +364,8 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.of_full_pw_dim, 0); EXPECT_FALSE(param.inp.of_read_kernel); EXPECT_EQ(param.inp.of_kernel_file, "WTkernel.txt"); + EXPECT_FALSE(param.inp.of_cd); + EXPECT_DOUBLE_EQ(param.inp.of_mCD_alpha,1.0); EXPECT_DOUBLE_EQ(param.inp.of_xwm_kappa, 1.); EXPECT_DOUBLE_EQ(param.inp.of_xwm_rho_ref, 1.); EXPECT_EQ(param.inp.device, "cpu"); diff --git a/source/source_io/test/support/INPUT b/source/source_io/test/support/INPUT index d69b67a6b3..46840bd9da 100644 --- a/source/source_io/test/support/INPUT +++ b/source/source_io/test/support/INPUT @@ -387,3 +387,7 @@ nsc_min 4 #Minimum number of spin-constrained iteration sc_scf_nmin 4 #Minimum number of outer scf loop before initializing lambda loop alpha_trial 0.02 #Initial trial step size for lambda in eV/uB^2 sccut 4 #Maximal step size for lambda in eV/uB + +#Parameters (23. Time-dependent orbital-free DFT) +of_cd 0 #0: no CD potential; 1: add CD potential +of_mCD_alpha 1.0 # parameter of modified CD potential \ No newline at end of file diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index a31768e9e8..bf5fbb5e78 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -25,7 +25,11 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec, } pelec->pot->update_from_charge(&chr, &ucell); // Hartree + XC + external - this->cal_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 + if (PARAM.inp.of_cd) + { + this->cal_CD_potential(psi_, pw_rho, pelec->pot->get_effective_v(), PARAM.inp.of_mCD_alpha); // CD potential + } #ifdef _OPENMP #pragma omp parallel for @@ -118,7 +122,8 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector> pphi, void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, ModulePW::PW_Basis* pw_rho, - ModuleBase::matrix& rpot) + ModuleBase::matrix& rpot, + double mCD_para) { std::complex imag(0.0,1.0); std::complex** recipPhi = new std::complex*[PARAM.inp.nspin]; @@ -180,7 +185,7 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - rpot(0, ir) -= std::real(2*rCDPotential[ir]*std::pow(ModuleBase::PI,3) / (2*std::pow(kF_r[ir],2))); + rpot(0, ir) -= mCD_para*std::real(2*rCDPotential[ir]*std::pow(ModuleBase::PI,3) / (2*std::pow(kF_r[ir],2))); } delete[] recipCurrent_x; delete[] recipCurrent_y; diff --git a/source/source_pw/module_ofdft/evolve_ofdft.h b/source/source_pw/module_ofdft/evolve_ofdft.h index d3cd6d358b..45ffffe7ad 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.h +++ b/source/source_pw/module_ofdft/evolve_ofdft.h @@ -48,7 +48,8 @@ class Evolve_OFDFT std::vector> Hpsi); // -1/2 \nabla^2 \phi void cal_CD_potential(std::vector> psi_, ModulePW::PW_Basis* pw_rho, - ModuleBase::matrix& rpot); + ModuleBase::matrix& rpot, + double mCD_para); }; #endif \ No newline at end of file From 37d3bf5f90499e803c57290df244d245448b7e3e Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sun, 12 Oct 2025 08:07:30 +0800 Subject: [PATCH 11/25] add integrate test for TDOFDFT --- tests/17_TDOFDFT/01_TDOFDFT_Al/INPUT | 36 ++++++++++++++++++ tests/17_TDOFDFT/01_TDOFDFT_Al/KPT | 4 ++ tests/17_TDOFDFT/01_TDOFDFT_Al/README | 1 + tests/17_TDOFDFT/01_TDOFDFT_Al/STRU | 19 ++++++++++ tests/17_TDOFDFT/01_TDOFDFT_Al/result.ref | 5 +++ tests/17_TDOFDFT/02_TDOFDFT_Al_CD/INPUT | 37 +++++++++++++++++++ tests/17_TDOFDFT/02_TDOFDFT_Al_CD/KPT | 4 ++ tests/17_TDOFDFT/02_TDOFDFT_Al_CD/README | 1 + tests/17_TDOFDFT/02_TDOFDFT_Al_CD/STRU | 19 ++++++++++ tests/17_TDOFDFT/02_TDOFDFT_Al_CD/result.ref | 5 +++ tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/INPUT | 37 +++++++++++++++++++ tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/KPT | 4 ++ tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/README | 1 + tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/STRU | 19 ++++++++++ tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/result.ref | 5 +++ tests/17_TDOFDFT/CASES_CPU.txt | 3 ++ tests/17_TDOFDFT/CMakeLists.txt | 16 ++++++++ tests/CMakeLists.txt | 1 + 18 files changed, 217 insertions(+) create mode 100644 tests/17_TDOFDFT/01_TDOFDFT_Al/INPUT create mode 100644 tests/17_TDOFDFT/01_TDOFDFT_Al/KPT create mode 100644 tests/17_TDOFDFT/01_TDOFDFT_Al/README create mode 100644 tests/17_TDOFDFT/01_TDOFDFT_Al/STRU create mode 100644 tests/17_TDOFDFT/01_TDOFDFT_Al/result.ref create mode 100644 tests/17_TDOFDFT/02_TDOFDFT_Al_CD/INPUT create mode 100644 tests/17_TDOFDFT/02_TDOFDFT_Al_CD/KPT create mode 100644 tests/17_TDOFDFT/02_TDOFDFT_Al_CD/README create mode 100644 tests/17_TDOFDFT/02_TDOFDFT_Al_CD/STRU create mode 100644 tests/17_TDOFDFT/02_TDOFDFT_Al_CD/result.ref create mode 100644 tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/INPUT create mode 100644 tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/KPT create mode 100644 tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/README create mode 100644 tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/STRU create mode 100644 tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/result.ref create mode 100644 tests/17_TDOFDFT/CASES_CPU.txt create mode 100644 tests/17_TDOFDFT/CMakeLists.txt diff --git a/tests/17_TDOFDFT/01_TDOFDFT_Al/INPUT b/tests/17_TDOFDFT/01_TDOFDFT_Al/INPUT new file mode 100644 index 0000000000..17fa540e80 --- /dev/null +++ b/tests/17_TDOFDFT/01_TDOFDFT_Al/INPUT @@ -0,0 +1,36 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation md +esolver_type tdofdft + +pseudo_dir ../../PP_ORB +pseudo_rcut 16 +cal_force 1 +cal_stress 1 + +#Parameters (2.Iteration) +ecutwfc 17 +scf_nmax 100 + +#OFDFT +of_kinetic tf+ +of_method tn +of_full_pw 1 +of_full_pw_dim 1 +of_cd 0 + +#Parameters (3.Basis) +basis_type pw + +init_vel 1 + +md_restart 0 +md_type nve +md_nstep 2 +md_dt 0.25 +md_tfirst 58022.52706 +md_dumpfreq 10 +md_tfreq 1.08 +md_tchain 1 +nbspline 10 # be sure fft dimension is odd diff --git a/tests/17_TDOFDFT/01_TDOFDFT_Al/KPT b/tests/17_TDOFDFT/01_TDOFDFT_Al/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/17_TDOFDFT/01_TDOFDFT_Al/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/17_TDOFDFT/01_TDOFDFT_Al/README b/tests/17_TDOFDFT/01_TDOFDFT_Al/README new file mode 100644 index 0000000000..be31e1c21f --- /dev/null +++ b/tests/17_TDOFDFT/01_TDOFDFT_Al/README @@ -0,0 +1 @@ +Test case without current dependent potential in TDOFDFT diff --git a/tests/17_TDOFDFT/01_TDOFDFT_Al/STRU b/tests/17_TDOFDFT/01_TDOFDFT_Al/STRU new file mode 100644 index 0000000000..1498bc2cc1 --- /dev/null +++ b/tests/17_TDOFDFT/01_TDOFDFT_Al/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Al 26.98 al.lda.lps blps + +LATTICE_CONSTANT +7.50241114482312 // add lattice constant + +LATTICE_VECTORS +0.000000000000 0.500000000000 0.500000000000 +0.500000000000 0.000000000000 0.500000000000 +0.500000000000 0.500000000000 0.000000000000 + +ATOMIC_POSITIONS +Direct + +Al +0 +2 +0.000000000000 0.000000000000 0.000000000000 1 1 1 v 0 0 0.5 +0.500000000000 0.500000000000 0.500000000000 1 1 1 v 0 0 -0.5 diff --git a/tests/17_TDOFDFT/01_TDOFDFT_Al/result.ref b/tests/17_TDOFDFT/01_TDOFDFT_Al/result.ref new file mode 100644 index 0000000000..404887187c --- /dev/null +++ b/tests/17_TDOFDFT/01_TDOFDFT_Al/result.ref @@ -0,0 +1,5 @@ +etotref -100.6027239728818614 +etotperatomref -50.3013619864 +totalforceref 5.786358 +totalstressref 11367.508779 +totaltimeref 0.10 diff --git a/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/INPUT b/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/INPUT new file mode 100644 index 0000000000..b14e49962e --- /dev/null +++ b/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/INPUT @@ -0,0 +1,37 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation md +esolver_type tdofdft + +pseudo_dir ../../PP_ORB +pseudo_rcut 16 +cal_force 1 +cal_stress 1 + +#Parameters (2.Iteration) +ecutwfc 17 +scf_nmax 100 + +#OFDFT +of_kinetic tf+ +of_method tn +of_full_pw 1 +of_full_pw_dim 1 +of_cd 1 +of_mcd_alpha 1.0 + +#Parameters (3.Basis) +basis_type pw + +init_vel 1 + +md_restart 0 +md_type nve +md_nstep 2 +md_dt 0.25 +md_tfirst 58022.52706 +md_dumpfreq 10 +md_tfreq 1.08 +md_tchain 1 +nbspline 10 # be sure fft dimension is odd diff --git a/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/KPT b/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/README b/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/README new file mode 100644 index 0000000000..d09f6a4ab9 --- /dev/null +++ b/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/README @@ -0,0 +1 @@ +Test case with current dependent potential in TDOFDFT diff --git a/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/STRU b/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/STRU new file mode 100644 index 0000000000..1498bc2cc1 --- /dev/null +++ b/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Al 26.98 al.lda.lps blps + +LATTICE_CONSTANT +7.50241114482312 // add lattice constant + +LATTICE_VECTORS +0.000000000000 0.500000000000 0.500000000000 +0.500000000000 0.000000000000 0.500000000000 +0.500000000000 0.500000000000 0.000000000000 + +ATOMIC_POSITIONS +Direct + +Al +0 +2 +0.000000000000 0.000000000000 0.000000000000 1 1 1 v 0 0 0.5 +0.500000000000 0.500000000000 0.500000000000 1 1 1 v 0 0 -0.5 diff --git a/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/result.ref b/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/result.ref new file mode 100644 index 0000000000..b94d694de2 --- /dev/null +++ b/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/result.ref @@ -0,0 +1,5 @@ +etotref -100.6027239728818614 +etotperatomref -50.3013619864 +totalforceref 5.786358 +totalstressref 11367.508779 +totaltimeref 0.13 diff --git a/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/INPUT b/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/INPUT new file mode 100644 index 0000000000..b4c1593e9f --- /dev/null +++ b/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/INPUT @@ -0,0 +1,37 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation md +esolver_type tdofdft + +pseudo_dir ../../PP_ORB +pseudo_rcut 16 +cal_force 1 +cal_stress 1 + +#Parameters (2.Iteration) +ecutwfc 17 +scf_nmax 100 + +#OFDFT +of_kinetic tf+ +of_method tn +of_full_pw 1 +of_full_pw_dim 1 +of_cd 1 +of_mcd_alpha 2.0 + +#Parameters (3.Basis) +basis_type pw + +init_vel 1 + +md_restart 0 +md_type nve +md_nstep 2 +md_dt 0.25 +md_tfirst 58022.52706 +md_dumpfreq 10 +md_tfreq 1.08 +md_tchain 1 +nbspline 10 # be sure fft dimension is odd diff --git a/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/KPT b/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/README b/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/README new file mode 100644 index 0000000000..cedbc26afd --- /dev/null +++ b/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/README @@ -0,0 +1 @@ +Test case with modified current dependent potential (alpha=2.0) in TDOFDFT diff --git a/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/STRU b/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/STRU new file mode 100644 index 0000000000..1498bc2cc1 --- /dev/null +++ b/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Al 26.98 al.lda.lps blps + +LATTICE_CONSTANT +7.50241114482312 // add lattice constant + +LATTICE_VECTORS +0.000000000000 0.500000000000 0.500000000000 +0.500000000000 0.000000000000 0.500000000000 +0.500000000000 0.500000000000 0.000000000000 + +ATOMIC_POSITIONS +Direct + +Al +0 +2 +0.000000000000 0.000000000000 0.000000000000 1 1 1 v 0 0 0.5 +0.500000000000 0.500000000000 0.500000000000 1 1 1 v 0 0 -0.5 diff --git a/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/result.ref b/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/result.ref new file mode 100644 index 0000000000..3b61275a93 --- /dev/null +++ b/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/result.ref @@ -0,0 +1,5 @@ +etotref -100.6027239728818614 +etotperatomref -50.3013619864 +totalforceref 5.786358 +totalstressref 11367.508779 +totaltimeref 0.11 diff --git a/tests/17_TDOFDFT/CASES_CPU.txt b/tests/17_TDOFDFT/CASES_CPU.txt new file mode 100644 index 0000000000..5f7ca560b1 --- /dev/null +++ b/tests/17_TDOFDFT/CASES_CPU.txt @@ -0,0 +1,3 @@ +01_TDOFDFT_Al +02_TDOFDFT_Al_CD +03_TDOFDFT_Al_mCD diff --git a/tests/17_TDOFDFT/CMakeLists.txt b/tests/17_TDOFDFT/CMakeLists.txt new file mode 100644 index 0000000000..b784d88eea --- /dev/null +++ b/tests/17_TDOFDFT/CMakeLists.txt @@ -0,0 +1,16 @@ +enable_testing() + +find_program(BASH bash) +if(ENABLE_ASAN) + add_test( + NAME 17_TDOFDFT_test_with_asan + COMMAND ${BASH} ../integrate/Autotest.sh -a ${ABACUS_BIN_PATH} -n 2 -s true + WORKING_DIRECTORY ${ABACUS_TEST_DIR}/17_TDOFDFT + ) +else() + add_test( + NAME 17_TDOFDFT + COMMAND ${BASH} ../integrate/Autotest.sh -a ${ABACUS_BIN_PATH} -n 4 + WORKING_DIRECTORY ${ABACUS_TEST_DIR}/17_TDOFDFT + ) +endif() diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5b1022bddc..e48f2315b9 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -9,6 +9,7 @@ add_subdirectory(07_OFDFT) add_subdirectory(08_EXX) add_subdirectory(10_others) add_subdirectory(11_PW_GPU) +add_subdirectory(17_TDOFDFT) if(ENABLE_MLALGO) add_subdirectory(09_DeePKS) From d8dead8c1a70023adc1ae657f986fc20d70f0b19 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sun, 12 Oct 2025 08:08:27 +0800 Subject: [PATCH 12/25] update doc file for tdofdft --- docs/advanced/input_files/input-main.md | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 665c9f9da6..6968bc99d7 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -256,6 +256,9 @@ - [of\_ml\_chi\_pnl](#of_ml_chi_pnl) - [of\_ml\_chi\_qnl](#of_ml_chi_qnl) - [of\_ml\_local\_test](#of_ml_local_test) + - [TD-OFDFT: time dependent orbital free density functional theory](#tdofdft-time-dependent-orbital-free-density-functional-theory) + - [of\_cd](#of_cd) + - [of\_mcd\_alpha](#of_mcd_alpha) - [Electric Field and Dipole Correction](#electric-field-and-dipole-correction) - [efield\_flag](#efield_flag) - [dip\_cor\_flag](#dip_cor_flag) @@ -528,6 +531,7 @@ These variables are used to control general system parameters. - **Description**: choose the energy solver. - ksdft: Kohn-Sham density functional theory - ofdft: orbital-free density functional theory + - tdofdft: time-dependent orbital-free density functional theory - sdft: [stochastic density functional theory](#electronic-structure-sdft) - tddft: real-time time-dependent density functional theory (TDDFT) - lj: Leonard Jones potential @@ -2729,6 +2733,25 @@ Warning: this function is not robust enough for the current version. Please try [back to top](#full-list-of-input-keywords) +## TDOFDFT: time dependent orbital free density functional theory + +### of_cd + +- **Type**: Boolean +- **Availability**: TDOFDFT +- **Type**: Boolean +- **Description**: Added the current dependent(CD) potential. (https://doi.org/10.1103/PhysRevB.98.144302) + - True: Added the CD potential. + - False: Not added the CD potential. +- **Default**: False + +### of_mcd_alpha + +- **Type**: Real +- **Availability**: TDOFDFT +- **Description**: The value of the parameter alpha in modified CD potential method. mCDPotenial=alpha*CDPotenial(proposed in paper PhysRevB.98.144302) +- **Default**: 1.0 + ## Electric field and dipole correction These variables are relevant to electric field and dipole correction From 51644c2e75f17aae44d30ec51a23c305f5a2859d Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sun, 12 Oct 2025 08:24:07 +0800 Subject: [PATCH 13/25] solve conflict --- source/source_pw/module_ofdft/evolve_ofdft.cpp | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index 4023ecfbdf..bf5fbb5e78 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -122,7 +122,6 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector> pphi, void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, ModulePW::PW_Basis* pw_rho, -<<<<<<< HEAD ModuleBase::matrix& rpot, double mCD_para) { @@ -203,15 +202,6 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, } delete[] recipPhi; delete[] rPhi; -======= - ModuleBase::matrix& rpot) -{ - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - //recipCurrent = new std::complex[pw_rho->npw]; - //delete[] recipCurrent; - } ->>>>>>> c66790bd1da1068b928178ce993d68a8fa4965b3 } void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, From 22b88fcd7bbe350c72190ea861cea081f7ae5eb1 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sun, 12 Oct 2025 08:36:10 +0800 Subject: [PATCH 14/25] type transformation --- source/source_pw/module_ofdft/evolve_ofdft.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index bf5fbb5e78..89065ad59c 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -185,7 +185,7 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - rpot(0, ir) -= mCD_para*std::real(2*rCDPotential[ir]*std::pow(ModuleBase::PI,3) / (2*std::pow(kF_r[ir],2))); + rpot(0, ir) -= mCD_para*std::real(2.0*rCDPotential[ir]*std::pow(ModuleBase::PI,3) / (2*std::pow(kF_r[ir],2))); } delete[] recipCurrent_x; delete[] recipCurrent_y; From a917e0b58ea0eb5d6f403437499dd66f208f085a Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sun, 12 Oct 2025 08:41:16 +0800 Subject: [PATCH 15/25] type transformation --- source/source_pw/module_ofdft/evolve_ofdft.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index 89065ad59c..bfc61cd5cf 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -185,7 +185,7 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - rpot(0, ir) -= mCD_para*std::real(2.0*rCDPotential[ir]*std::pow(ModuleBase::PI,3) / (2*std::pow(kF_r[ir],2))); + rpot(0, ir) -= mCD_para*2.0*std::real(rCDPotential[ir])*std::pow(ModuleBase::PI,3) / (2.0*std::pow(std::real(kF_r[ir]),2)); } delete[] recipCurrent_x; delete[] recipCurrent_y; From 4c0a683aa9e3ce3c299fb154f668d9652bf7739b Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sun, 12 Oct 2025 17:36:37 +0800 Subject: [PATCH 16/25] move integration test for TDOFDFT to 07_OFDFT --- .../01_TDOFDFT_Al => 07_OFDFT/30_TDDFT_Al}/INPUT | 0 .../01_TDOFDFT_Al => 07_OFDFT/30_TDDFT_Al}/KPT | 0 .../30_TDDFT_Al}/README | 0 .../01_TDOFDFT_Al => 07_OFDFT/30_TDDFT_Al}/STRU | 0 .../30_TDDFT_Al}/result.ref | 0 .../31_TDDFT_Al_CD}/INPUT | 0 .../31_TDDFT_Al_CD}/KPT | 0 .../31_TDDFT_Al_CD}/README | 0 .../31_TDDFT_Al_CD}/STRU | 0 .../31_TDDFT_Al_CD}/result.ref | 0 .../32_TDDFT_Al_mCD}/INPUT | 0 .../32_TDDFT_Al_mCD}/KPT | 0 .../32_TDDFT_Al_mCD}/README | 0 .../32_TDDFT_Al_mCD}/STRU | 0 .../32_TDDFT_Al_mCD}/result.ref | 0 tests/07_OFDFT/CASES_CPU.txt | 5 ++++- tests/17_TDOFDFT/CASES_CPU.txt | 3 --- tests/17_TDOFDFT/CMakeLists.txt | 16 ---------------- tests/CMakeLists.txt | 1 - 19 files changed, 4 insertions(+), 21 deletions(-) rename tests/{17_TDOFDFT/01_TDOFDFT_Al => 07_OFDFT/30_TDDFT_Al}/INPUT (100%) rename tests/{17_TDOFDFT/01_TDOFDFT_Al => 07_OFDFT/30_TDDFT_Al}/KPT (100%) rename tests/{17_TDOFDFT/01_TDOFDFT_Al => 07_OFDFT/30_TDDFT_Al}/README (100%) rename tests/{17_TDOFDFT/01_TDOFDFT_Al => 07_OFDFT/30_TDDFT_Al}/STRU (100%) rename tests/{17_TDOFDFT/01_TDOFDFT_Al => 07_OFDFT/30_TDDFT_Al}/result.ref (100%) rename tests/{17_TDOFDFT/02_TDOFDFT_Al_CD => 07_OFDFT/31_TDDFT_Al_CD}/INPUT (100%) rename tests/{17_TDOFDFT/02_TDOFDFT_Al_CD => 07_OFDFT/31_TDDFT_Al_CD}/KPT (100%) rename tests/{17_TDOFDFT/02_TDOFDFT_Al_CD => 07_OFDFT/31_TDDFT_Al_CD}/README (100%) rename tests/{17_TDOFDFT/02_TDOFDFT_Al_CD => 07_OFDFT/31_TDDFT_Al_CD}/STRU (100%) rename tests/{17_TDOFDFT/02_TDOFDFT_Al_CD => 07_OFDFT/31_TDDFT_Al_CD}/result.ref (100%) rename tests/{17_TDOFDFT/03_TDOFDFT_Al_mCD => 07_OFDFT/32_TDDFT_Al_mCD}/INPUT (100%) rename tests/{17_TDOFDFT/03_TDOFDFT_Al_mCD => 07_OFDFT/32_TDDFT_Al_mCD}/KPT (100%) rename tests/{17_TDOFDFT/03_TDOFDFT_Al_mCD => 07_OFDFT/32_TDDFT_Al_mCD}/README (100%) rename tests/{17_TDOFDFT/03_TDOFDFT_Al_mCD => 07_OFDFT/32_TDDFT_Al_mCD}/STRU (100%) rename tests/{17_TDOFDFT/03_TDOFDFT_Al_mCD => 07_OFDFT/32_TDDFT_Al_mCD}/result.ref (100%) delete mode 100644 tests/17_TDOFDFT/CASES_CPU.txt delete mode 100644 tests/17_TDOFDFT/CMakeLists.txt diff --git a/tests/17_TDOFDFT/01_TDOFDFT_Al/INPUT b/tests/07_OFDFT/30_TDDFT_Al/INPUT similarity index 100% rename from tests/17_TDOFDFT/01_TDOFDFT_Al/INPUT rename to tests/07_OFDFT/30_TDDFT_Al/INPUT diff --git a/tests/17_TDOFDFT/01_TDOFDFT_Al/KPT b/tests/07_OFDFT/30_TDDFT_Al/KPT similarity index 100% rename from tests/17_TDOFDFT/01_TDOFDFT_Al/KPT rename to tests/07_OFDFT/30_TDDFT_Al/KPT diff --git a/tests/17_TDOFDFT/01_TDOFDFT_Al/README b/tests/07_OFDFT/30_TDDFT_Al/README similarity index 100% rename from tests/17_TDOFDFT/01_TDOFDFT_Al/README rename to tests/07_OFDFT/30_TDDFT_Al/README diff --git a/tests/17_TDOFDFT/01_TDOFDFT_Al/STRU b/tests/07_OFDFT/30_TDDFT_Al/STRU similarity index 100% rename from tests/17_TDOFDFT/01_TDOFDFT_Al/STRU rename to tests/07_OFDFT/30_TDDFT_Al/STRU diff --git a/tests/17_TDOFDFT/01_TDOFDFT_Al/result.ref b/tests/07_OFDFT/30_TDDFT_Al/result.ref similarity index 100% rename from tests/17_TDOFDFT/01_TDOFDFT_Al/result.ref rename to tests/07_OFDFT/30_TDDFT_Al/result.ref diff --git a/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/INPUT b/tests/07_OFDFT/31_TDDFT_Al_CD/INPUT similarity index 100% rename from tests/17_TDOFDFT/02_TDOFDFT_Al_CD/INPUT rename to tests/07_OFDFT/31_TDDFT_Al_CD/INPUT diff --git a/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/KPT b/tests/07_OFDFT/31_TDDFT_Al_CD/KPT similarity index 100% rename from tests/17_TDOFDFT/02_TDOFDFT_Al_CD/KPT rename to tests/07_OFDFT/31_TDDFT_Al_CD/KPT diff --git a/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/README b/tests/07_OFDFT/31_TDDFT_Al_CD/README similarity index 100% rename from tests/17_TDOFDFT/02_TDOFDFT_Al_CD/README rename to tests/07_OFDFT/31_TDDFT_Al_CD/README diff --git a/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/STRU b/tests/07_OFDFT/31_TDDFT_Al_CD/STRU similarity index 100% rename from tests/17_TDOFDFT/02_TDOFDFT_Al_CD/STRU rename to tests/07_OFDFT/31_TDDFT_Al_CD/STRU diff --git a/tests/17_TDOFDFT/02_TDOFDFT_Al_CD/result.ref b/tests/07_OFDFT/31_TDDFT_Al_CD/result.ref similarity index 100% rename from tests/17_TDOFDFT/02_TDOFDFT_Al_CD/result.ref rename to tests/07_OFDFT/31_TDDFT_Al_CD/result.ref diff --git a/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/INPUT b/tests/07_OFDFT/32_TDDFT_Al_mCD/INPUT similarity index 100% rename from tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/INPUT rename to tests/07_OFDFT/32_TDDFT_Al_mCD/INPUT diff --git a/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/KPT b/tests/07_OFDFT/32_TDDFT_Al_mCD/KPT similarity index 100% rename from tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/KPT rename to tests/07_OFDFT/32_TDDFT_Al_mCD/KPT diff --git a/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/README b/tests/07_OFDFT/32_TDDFT_Al_mCD/README similarity index 100% rename from tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/README rename to tests/07_OFDFT/32_TDDFT_Al_mCD/README diff --git a/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/STRU b/tests/07_OFDFT/32_TDDFT_Al_mCD/STRU similarity index 100% rename from tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/STRU rename to tests/07_OFDFT/32_TDDFT_Al_mCD/STRU diff --git a/tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/result.ref b/tests/07_OFDFT/32_TDDFT_Al_mCD/result.ref similarity index 100% rename from tests/17_TDOFDFT/03_TDOFDFT_Al_mCD/result.ref rename to tests/07_OFDFT/32_TDDFT_Al_mCD/result.ref diff --git a/tests/07_OFDFT/CASES_CPU.txt b/tests/07_OFDFT/CASES_CPU.txt index 67e095d710..e94290c382 100644 --- a/tests/07_OFDFT/CASES_CPU.txt +++ b/tests/07_OFDFT/CASES_CPU.txt @@ -26,4 +26,7 @@ 26_OF_MD_LibxcPBE 27_OF_CR_LDA 28_OF_KE_XWM -29_OF_XWM_para \ No newline at end of file +29_OF_XWM_para +30_TDOFDFT_Al +31_TDOFDFT_Al_CD +32_TDOFDFT_Al_mCD diff --git a/tests/17_TDOFDFT/CASES_CPU.txt b/tests/17_TDOFDFT/CASES_CPU.txt deleted file mode 100644 index 5f7ca560b1..0000000000 --- a/tests/17_TDOFDFT/CASES_CPU.txt +++ /dev/null @@ -1,3 +0,0 @@ -01_TDOFDFT_Al -02_TDOFDFT_Al_CD -03_TDOFDFT_Al_mCD diff --git a/tests/17_TDOFDFT/CMakeLists.txt b/tests/17_TDOFDFT/CMakeLists.txt deleted file mode 100644 index b784d88eea..0000000000 --- a/tests/17_TDOFDFT/CMakeLists.txt +++ /dev/null @@ -1,16 +0,0 @@ -enable_testing() - -find_program(BASH bash) -if(ENABLE_ASAN) - add_test( - NAME 17_TDOFDFT_test_with_asan - COMMAND ${BASH} ../integrate/Autotest.sh -a ${ABACUS_BIN_PATH} -n 2 -s true - WORKING_DIRECTORY ${ABACUS_TEST_DIR}/17_TDOFDFT - ) -else() - add_test( - NAME 17_TDOFDFT - COMMAND ${BASH} ../integrate/Autotest.sh -a ${ABACUS_BIN_PATH} -n 4 - WORKING_DIRECTORY ${ABACUS_TEST_DIR}/17_TDOFDFT - ) -endif() diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index e48f2315b9..5b1022bddc 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -9,7 +9,6 @@ add_subdirectory(07_OFDFT) add_subdirectory(08_EXX) add_subdirectory(10_others) add_subdirectory(11_PW_GPU) -add_subdirectory(17_TDOFDFT) if(ENABLE_MLALGO) add_subdirectory(09_DeePKS) From 824e014f074e973658787b7598422243146656cb Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sun, 12 Oct 2025 17:45:34 +0800 Subject: [PATCH 17/25] add judgement for nspin in TDOFDFT --- source/source_pw/module_ofdft/evolve_ofdft.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index bfc61cd5cf..d6424348d1 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -126,6 +126,10 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, double mCD_para) { std::complex imag(0.0,1.0); + + if (nspin <= 0) { + throw std::invalid_argument("nspin must be positive"); + } std::complex** recipPhi = new std::complex*[PARAM.inp.nspin]; std::complex** rPhi = new std::complex*[PARAM.inp.nspin]; #ifdef _OPENMP @@ -139,6 +143,9 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, } } +#ifdef _OPENMP +#pragma omp parallel for +#endif for (int is = 0; is < PARAM.inp.nspin; ++is) { std::complex *recipCurrent_x=new std::complex[pw_rho->npw]; @@ -195,6 +202,9 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, delete[] rCurrent_z; } +#ifdef _OPENMP +#pragma omp parallel for +#endif for (int is = 0; is < PARAM.inp.nspin; ++is) { delete[] recipPhi[is]; From 040af900644cf456c04d28b95bb2518841a961b7 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sun, 12 Oct 2025 17:47:35 +0800 Subject: [PATCH 18/25] add judgement for nspin in TDOFDFT --- source/source_pw/module_ofdft/evolve_ofdft.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index d6424348d1..8b4448655b 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -76,6 +76,9 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector> pphi, ModulePW::PW_Basis* pw_rho, std::vector> Hpsi) { + if (PARAM.inp.nspin <= 0) { + throw std::invalid_argument("nspin must be positive"); + } std::complex** rLapPhi = new std::complex*[PARAM.inp.nspin]; #ifdef _OPENMP #pragma omp parallel for @@ -127,7 +130,7 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, { std::complex imag(0.0,1.0); - if (nspin <= 0) { + if (PARAM.inp.nspin <= 0) { throw std::invalid_argument("nspin must be positive"); } std::complex** recipPhi = new std::complex*[PARAM.inp.nspin]; From e3fdd3e51920ec058f88b930cfa510c201de9777 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Sun, 12 Oct 2025 18:35:10 +0800 Subject: [PATCH 19/25] TDOFDFT integration test --- tests/07_OFDFT/{30_TDDFT_Al => 30_TDOFDFT_Al}/INPUT | 0 tests/07_OFDFT/{30_TDDFT_Al => 30_TDOFDFT_Al}/KPT | 0 tests/07_OFDFT/{30_TDDFT_Al => 30_TDOFDFT_Al}/README | 0 tests/07_OFDFT/{30_TDDFT_Al => 30_TDOFDFT_Al}/STRU | 0 tests/07_OFDFT/{30_TDDFT_Al => 30_TDOFDFT_Al}/result.ref | 0 tests/07_OFDFT/{31_TDDFT_Al_CD => 31_TDOFDFT_Al_CD}/INPUT | 0 tests/07_OFDFT/{31_TDDFT_Al_CD => 31_TDOFDFT_Al_CD}/KPT | 0 tests/07_OFDFT/{31_TDDFT_Al_CD => 31_TDOFDFT_Al_CD}/README | 0 tests/07_OFDFT/{31_TDDFT_Al_CD => 31_TDOFDFT_Al_CD}/STRU | 0 tests/07_OFDFT/{31_TDDFT_Al_CD => 31_TDOFDFT_Al_CD}/result.ref | 0 tests/07_OFDFT/{32_TDDFT_Al_mCD => 32_TDOFDFT_Al_mCD}/INPUT | 0 tests/07_OFDFT/{32_TDDFT_Al_mCD => 32_TDOFDFT_Al_mCD}/KPT | 0 tests/07_OFDFT/{32_TDDFT_Al_mCD => 32_TDOFDFT_Al_mCD}/README | 0 tests/07_OFDFT/{32_TDDFT_Al_mCD => 32_TDOFDFT_Al_mCD}/STRU | 0 tests/07_OFDFT/{32_TDDFT_Al_mCD => 32_TDOFDFT_Al_mCD}/result.ref | 0 15 files changed, 0 insertions(+), 0 deletions(-) rename tests/07_OFDFT/{30_TDDFT_Al => 30_TDOFDFT_Al}/INPUT (100%) rename tests/07_OFDFT/{30_TDDFT_Al => 30_TDOFDFT_Al}/KPT (100%) rename tests/07_OFDFT/{30_TDDFT_Al => 30_TDOFDFT_Al}/README (100%) rename tests/07_OFDFT/{30_TDDFT_Al => 30_TDOFDFT_Al}/STRU (100%) rename tests/07_OFDFT/{30_TDDFT_Al => 30_TDOFDFT_Al}/result.ref (100%) rename tests/07_OFDFT/{31_TDDFT_Al_CD => 31_TDOFDFT_Al_CD}/INPUT (100%) rename tests/07_OFDFT/{31_TDDFT_Al_CD => 31_TDOFDFT_Al_CD}/KPT (100%) rename tests/07_OFDFT/{31_TDDFT_Al_CD => 31_TDOFDFT_Al_CD}/README (100%) rename tests/07_OFDFT/{31_TDDFT_Al_CD => 31_TDOFDFT_Al_CD}/STRU (100%) rename tests/07_OFDFT/{31_TDDFT_Al_CD => 31_TDOFDFT_Al_CD}/result.ref (100%) rename tests/07_OFDFT/{32_TDDFT_Al_mCD => 32_TDOFDFT_Al_mCD}/INPUT (100%) rename tests/07_OFDFT/{32_TDDFT_Al_mCD => 32_TDOFDFT_Al_mCD}/KPT (100%) rename tests/07_OFDFT/{32_TDDFT_Al_mCD => 32_TDOFDFT_Al_mCD}/README (100%) rename tests/07_OFDFT/{32_TDDFT_Al_mCD => 32_TDOFDFT_Al_mCD}/STRU (100%) rename tests/07_OFDFT/{32_TDDFT_Al_mCD => 32_TDOFDFT_Al_mCD}/result.ref (100%) diff --git a/tests/07_OFDFT/30_TDDFT_Al/INPUT b/tests/07_OFDFT/30_TDOFDFT_Al/INPUT similarity index 100% rename from tests/07_OFDFT/30_TDDFT_Al/INPUT rename to tests/07_OFDFT/30_TDOFDFT_Al/INPUT diff --git a/tests/07_OFDFT/30_TDDFT_Al/KPT b/tests/07_OFDFT/30_TDOFDFT_Al/KPT similarity index 100% rename from tests/07_OFDFT/30_TDDFT_Al/KPT rename to tests/07_OFDFT/30_TDOFDFT_Al/KPT diff --git a/tests/07_OFDFT/30_TDDFT_Al/README b/tests/07_OFDFT/30_TDOFDFT_Al/README similarity index 100% rename from tests/07_OFDFT/30_TDDFT_Al/README rename to tests/07_OFDFT/30_TDOFDFT_Al/README diff --git a/tests/07_OFDFT/30_TDDFT_Al/STRU b/tests/07_OFDFT/30_TDOFDFT_Al/STRU similarity index 100% rename from tests/07_OFDFT/30_TDDFT_Al/STRU rename to tests/07_OFDFT/30_TDOFDFT_Al/STRU diff --git a/tests/07_OFDFT/30_TDDFT_Al/result.ref b/tests/07_OFDFT/30_TDOFDFT_Al/result.ref similarity index 100% rename from tests/07_OFDFT/30_TDDFT_Al/result.ref rename to tests/07_OFDFT/30_TDOFDFT_Al/result.ref diff --git a/tests/07_OFDFT/31_TDDFT_Al_CD/INPUT b/tests/07_OFDFT/31_TDOFDFT_Al_CD/INPUT similarity index 100% rename from tests/07_OFDFT/31_TDDFT_Al_CD/INPUT rename to tests/07_OFDFT/31_TDOFDFT_Al_CD/INPUT diff --git a/tests/07_OFDFT/31_TDDFT_Al_CD/KPT b/tests/07_OFDFT/31_TDOFDFT_Al_CD/KPT similarity index 100% rename from tests/07_OFDFT/31_TDDFT_Al_CD/KPT rename to tests/07_OFDFT/31_TDOFDFT_Al_CD/KPT diff --git a/tests/07_OFDFT/31_TDDFT_Al_CD/README b/tests/07_OFDFT/31_TDOFDFT_Al_CD/README similarity index 100% rename from tests/07_OFDFT/31_TDDFT_Al_CD/README rename to tests/07_OFDFT/31_TDOFDFT_Al_CD/README diff --git a/tests/07_OFDFT/31_TDDFT_Al_CD/STRU b/tests/07_OFDFT/31_TDOFDFT_Al_CD/STRU similarity index 100% rename from tests/07_OFDFT/31_TDDFT_Al_CD/STRU rename to tests/07_OFDFT/31_TDOFDFT_Al_CD/STRU diff --git a/tests/07_OFDFT/31_TDDFT_Al_CD/result.ref b/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref similarity index 100% rename from tests/07_OFDFT/31_TDDFT_Al_CD/result.ref rename to tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref diff --git a/tests/07_OFDFT/32_TDDFT_Al_mCD/INPUT b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/INPUT similarity index 100% rename from tests/07_OFDFT/32_TDDFT_Al_mCD/INPUT rename to tests/07_OFDFT/32_TDOFDFT_Al_mCD/INPUT diff --git a/tests/07_OFDFT/32_TDDFT_Al_mCD/KPT b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/KPT similarity index 100% rename from tests/07_OFDFT/32_TDDFT_Al_mCD/KPT rename to tests/07_OFDFT/32_TDOFDFT_Al_mCD/KPT diff --git a/tests/07_OFDFT/32_TDDFT_Al_mCD/README b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/README similarity index 100% rename from tests/07_OFDFT/32_TDDFT_Al_mCD/README rename to tests/07_OFDFT/32_TDOFDFT_Al_mCD/README diff --git a/tests/07_OFDFT/32_TDDFT_Al_mCD/STRU b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/STRU similarity index 100% rename from tests/07_OFDFT/32_TDDFT_Al_mCD/STRU rename to tests/07_OFDFT/32_TDOFDFT_Al_mCD/STRU diff --git a/tests/07_OFDFT/32_TDDFT_Al_mCD/result.ref b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref similarity index 100% rename from tests/07_OFDFT/32_TDDFT_Al_mCD/result.ref rename to tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref From b278f264102be38d9fc86ebb52f832f8d806b03c Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Mon, 13 Oct 2025 08:13:26 +0800 Subject: [PATCH 20/25] change type of warning --- source/source_pw/module_ofdft/evolve_ofdft.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index 8b4448655b..f4aaea2cec 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -77,7 +77,7 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector> pphi, std::vector> Hpsi) { if (PARAM.inp.nspin <= 0) { - throw std::invalid_argument("nspin must be positive"); + ModuleBase::WARNING_QUIT("Evolve_OFDFT","nspin must be positive"); } std::complex** rLapPhi = new std::complex*[PARAM.inp.nspin]; #ifdef _OPENMP @@ -131,7 +131,7 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, std::complex imag(0.0,1.0); if (PARAM.inp.nspin <= 0) { - throw std::invalid_argument("nspin must be positive"); + ModuleBase::WARNING_QUIT("Evolve_OFDFT","nspin must be positive"); } std::complex** recipPhi = new std::complex*[PARAM.inp.nspin]; std::complex** rPhi = new std::complex*[PARAM.inp.nspin]; From 1e56b712bf562f644b2c59b1522c4aca01369f9b Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Mon, 12 Jan 2026 13:35:54 +0800 Subject: [PATCH 21/25] fix bug of TDOFDFT --- source/source_esolver/esolver_of_tddft.cpp | 13 +- source/source_esolver/esolver_of_tddft.h | 2 +- .../source_pw/module_ofdft/evolve_ofdft.cpp | 181 ++++++++++++++---- source/source_pw/module_ofdft/evolve_ofdft.h | 27 ++- 4 files changed, 166 insertions(+), 57 deletions(-) diff --git a/source/source_esolver/esolver_of_tddft.cpp b/source/source_esolver/esolver_of_tddft.cpp index a978a1267c..daeda628cb 100644 --- a/source/source_esolver/esolver_of_tddft.cpp +++ b/source/source_esolver/esolver_of_tddft.cpp @@ -16,6 +16,8 @@ //-----stress------------------ #include "source_pw/module_ofdft/of_stress_pw.h" +#include + namespace ModuleESolver { @@ -45,12 +47,13 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) this->iter_time = std::chrono::system_clock::now(); #endif - if (istep==0) + if (this->phi_td.empty()) { - this->phi_td.resize(PARAM.inp.nspin*this->pw_rho->nrxx); + const int size = PARAM.inp.nspin * this->pw_rho->nrxx; + this->phi_td.resize(size, std::complex(0.0, 0.0)); } - if ((istep<1) && PARAM.inp.init_chg != "file") + if ((istep==0) && PARAM.inp.init_chg != "file") { while (true) { @@ -91,7 +94,7 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) } } } - else if ((istep<1) && PARAM.inp.init_chg == "file") + else if ((istep==0) && PARAM.inp.init_chg == "file") { #ifdef _OPENMP #pragma omp parallel for collapse(2) @@ -107,7 +110,7 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) } else { - this->evolve_ofdft->propagate_psi(this->pelec, this->chr, ucell, this->phi_td, this->pw_rho); + this->evolve_ofdft->propagate_psi_RK4(this->pelec, this->chr, ucell, this->phi_td, this->pw_rho); #ifdef _OPENMP #pragma omp parallel for collapse(2) #endif diff --git a/source/source_esolver/esolver_of_tddft.h b/source/source_esolver/esolver_of_tddft.h index 99a4238c0f..85293b1761 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> phi_td; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). + std::vector> phi_td; // time dependent wavefunction 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 f4aaea2cec..e849a69868 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -6,11 +6,11 @@ #include "source_base/parallel_reduce.h" void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec, - const Charge& chr, + Charge& chr, UnitCell& ucell, - std::vector> psi_, + std::vector>& psi_, ModulePW::PW_Basis* pw_rho, - std::vector> Hpsi) + std::vector>& Hpsi) { // update rho #ifdef _OPENMP @@ -23,6 +23,7 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec, chr.rho[is][ir] = abs(psi_[is * pw_rho->nrxx + ir])*abs(psi_[is * pw_rho->nrxx + ir]); } } + this->renormalize_psi(chr, pw_rho, psi_); pelec->pot->update_from_charge(&chr, &ucell); // Hartree + XC + external this->cal_tf_potential(chr.rho, pw_rho, pelec->pot->get_effective_v()); // TF potential @@ -45,6 +46,26 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec, this->cal_vw_potential_phi(psi_, pw_rho, Hpsi); } +void Evolve_OFDFT::renormalize_psi(Charge& chr, ModulePW::PW_Basis* pw_rho, std::vector>& pphi_) +{ + const double sr = chr.sum_rho(); + const double normalize_factor = PARAM.inp.nelec / sr; + + std::cout<<"sr="<> pphi, pw_rho->recip2real(recipPhi[is], rLapPhi[is]); for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - Hpsi[is * pw_rho->nrxx + ir]+=rLapPhi[is][ir]; + Hpsi[is * pw_rho->nrxx + ir] += rLapPhi[is][ir]; } } @@ -123,7 +144,7 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector> pphi, delete[] rLapPhi; } -void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, +void Evolve_OFDFT::cal_CD_potential(std::vector>& psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot, double mCD_para) @@ -151,20 +172,20 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, #endif for (int is = 0; is < PARAM.inp.nspin; ++is) { - std::complex *recipCurrent_x=new std::complex[pw_rho->npw]; - std::complex *recipCurrent_y=new std::complex[pw_rho->npw]; - std::complex *recipCurrent_z=new std::complex[pw_rho->npw]; - std::complex *recipCDPotential=new std::complex[pw_rho->npw]; - std::complex *rCurrent_x=new std::complex[pw_rho->nrxx]; - std::complex *rCurrent_y=new std::complex[pw_rho->nrxx]; - std::complex *rCurrent_z=new std::complex[pw_rho->nrxx]; - std::complex *kF_r=new std::complex[pw_rho->nrxx]; - std::complex *rCDPotential=new std::complex[pw_rho->nrxx]; + std::vector> recipCurrent_x(pw_rho->npw); + std::vector> recipCurrent_y(pw_rho->npw); + std::vector> recipCurrent_z(pw_rho->npw); + std::vector> recipCDPotential(pw_rho->npw); + std::vector> rCurrent_x(pw_rho->nrxx); + std::vector> rCurrent_y(pw_rho->nrxx); + std::vector> rCurrent_z(pw_rho->nrxx); + std::vector> kF_r(pw_rho->nrxx); + std::vector> rCDPotential(pw_rho->nrxx); recipPhi[is] = new std::complex[pw_rho->npw]; for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - kF_r[ir]=std::pow(3*std::pow(ModuleBase::PI*std::abs(rPhi[is][ir]),2),1/3); + kF_r[ir]=std::pow(3*std::pow(ModuleBase::PI*std::abs(rPhi[is][ir]),2),1.0/3.0); } pw_rho->real2recip(rPhi[is], recipPhi[is]); @@ -174,35 +195,40 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, recipCurrent_y[ik]=imag*pw_rho->gcar[ik].y*recipPhi[is][ik]* pw_rho->tpiba; recipCurrent_z[ik]=imag*pw_rho->gcar[ik].z*recipPhi[is][ik]* pw_rho->tpiba; } - pw_rho->recip2real(recipCurrent_x,rCurrent_x); - pw_rho->recip2real(recipCurrent_y,rCurrent_y); - pw_rho->recip2real(recipCurrent_z,rCurrent_z); + pw_rho->recip2real(recipCurrent_x.data(),rCurrent_x.data()); + pw_rho->recip2real(recipCurrent_y.data(),rCurrent_y.data()); + pw_rho->recip2real(recipCurrent_z.data(),rCurrent_z.data()); for (int ir = 0; ir < pw_rho->nrxx; ++ir) { rCurrent_x[ir]=std::imag(rCurrent_x[ir]*std::conj(rPhi[is][ir])); rCurrent_y[ir]=std::imag(rCurrent_y[ir]*std::conj(rPhi[is][ir])); rCurrent_z[ir]=std::imag(rCurrent_z[ir]*std::conj(rPhi[is][ir])); } - pw_rho->real2recip(rCurrent_x,recipCurrent_x); - pw_rho->real2recip(rCurrent_y,recipCurrent_y); - pw_rho->real2recip(rCurrent_z,recipCurrent_z); + pw_rho->real2recip(rCurrent_x.data(),recipCurrent_x.data()); + pw_rho->real2recip(rCurrent_y.data(),recipCurrent_y.data()); + pw_rho->real2recip(rCurrent_z.data(),recipCurrent_z.data()); for (int ik = 0; ik < pw_rho->npw; ++ik) { recipCDPotential[ik]=recipCurrent_x[ik]*pw_rho->gcar[ik].x+recipCurrent_y[ik]*pw_rho->gcar[ik].y+recipCurrent_z[ik]*pw_rho->gcar[ik].z; - recipCDPotential[ik]*=imag/pw_rho->gg[ik]; + if (pw_rho->gg[ik]==0) + { + recipCDPotential[ik]=0.0; + } + else + { + recipCDPotential[ik]*=imag/sqrt(pw_rho->gg[ik]); + } } - pw_rho->recip2real(recipCDPotential,rCDPotential); + pw_rho->recip2real(recipCDPotential.data(),rCDPotential.data()); for (int ir = 0; ir < pw_rho->nrxx; ++ir) { rpot(0, ir) -= mCD_para*2.0*std::real(rCDPotential[ir])*std::pow(ModuleBase::PI,3) / (2.0*std::pow(std::real(kF_r[ir]),2)); + if (isnan(rpot(0, ir))) + { + rpot(0, ir)=0.0; + } } - delete[] recipCurrent_x; - delete[] recipCurrent_y; - delete[] recipCurrent_z; - delete[] rCurrent_x; - delete[] rCurrent_y; - delete[] rCurrent_z; } #ifdef _OPENMP @@ -217,15 +243,16 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, delete[] rPhi; } -void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, - const Charge& chr, UnitCell& ucell, - std::vector> pphi_, +void Evolve_OFDFT::propagate_psi_RK4(elecstate::ElecState* pelec, + Charge& chr, + UnitCell& ucell, + std::vector>& pphi_, ModulePW::PW_Basis* pw_rho) { - ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagate_psi_RK4"); std::complex imag(0.0,1.0); - double dt=PARAM.inp.mdp.md_dt; + double dt=PARAM.inp.mdp.md_dt / ModuleBase::AU_to_FS; const int nspin = PARAM.inp.nspin; const int nrxx = pw_rho->nrxx; const int total_size = nspin * nrxx; @@ -244,7 +271,7 @@ void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K1[is * nrxx + ir]=-1.0*K1[is * nrxx + ir]*dt*imag; + K1[is * nrxx + ir]=-0.5*K1[is * nrxx + ir]*dt*imag; // 0.5 convert Ry to Hartree psi1[is * nrxx + ir]=pphi_[is * nrxx + ir]+0.5*K1[is * nrxx + ir]; } } @@ -255,7 +282,7 @@ void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K2[is * nrxx + ir]=-1.0*K2[is * nrxx + ir]*dt*imag; + K2[is * nrxx + ir]=-0.5*K2[is * nrxx + ir]*dt*imag; psi2[is * nrxx + ir]=pphi_[is * nrxx + ir]+0.5*K2[is * nrxx + ir]; } } @@ -266,7 +293,7 @@ void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K3[is * nrxx + ir]=-1.0*K3[is * nrxx + ir]*dt*imag; + K3[is * nrxx + ir]=-0.5*K3[is * nrxx + ir]*dt*imag; psi3[is * nrxx + ir]=pphi_[is * nrxx + ir]+K3[is * nrxx + ir]; } } @@ -277,10 +304,80 @@ void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K4[is * nrxx + ir]=-1.0*K4[is * nrxx + ir]*dt*imag; + K4[is * nrxx + ir]=-0.5*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]); } } - ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); -} \ No newline at end of file +#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(pphi_[is * pw_rho->nrxx + ir])*abs(pphi_[is * pw_rho->nrxx + ir]); + } + } + this->renormalize_psi(chr, pw_rho, pphi_); + + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagate_psi_RK4"); +} + +void Evolve_OFDFT::propagate_psi_RK2(elecstate::ElecState* pelec, + Charge& chr, + UnitCell& ucell, + std::vector>& pphi_, + ModulePW::PW_Basis* pw_rho) +{ + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagate_psi_RK2"); + + const std::complex imag(0.0, 1.0); + double dt=PARAM.inp.mdp.md_dt / ModuleBase::AU_to_FS; + 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> psi_mid(total_size); + + cal_Hpsi(pelec, chr, ucell, pphi_, pw_rho, K1); + +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif + for (int is = 0; is < nspin; ++is) { + for (int ir = 0; ir < nrxx; ++ir) { + const int idx = is * nrxx + ir; + K1[idx] = -0.5 * K1[idx] * dt * imag; + psi_mid[idx] = pphi_[idx] + 0.5 * K1[idx]; + } + } + + cal_Hpsi(pelec, chr, ucell, psi_mid, pw_rho, K2); + +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif + for (int is = 0; is < nspin; ++is) { + for (int ir = 0; ir < nrxx; ++ir) { + const int idx = is * nrxx + ir; + K2[idx] = -0.5 * K2[idx] * dt * imag; + pphi_[idx] += K2[idx]; + } + } + +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif + for (int is = 0; is < nspin; ++is) { + for (int ir = 0; ir < nrxx; ++ir) { + chr.rho[is][ir] = std::norm(pphi_[is * nrxx + ir]); + } + } + + this->renormalize_psi(chr, pw_rho, pphi_); + + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagate_psi_RK2"); +} diff --git a/source/source_pw/module_ofdft/evolve_ofdft.h b/source/source_pw/module_ofdft/evolve_ofdft.h index 45ffffe7ad..8d592a56a4 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.h +++ b/source/source_pw/module_ofdft/evolve_ofdft.h @@ -24,10 +24,19 @@ class Evolve_OFDFT ~Evolve_OFDFT() { } - void propagate_psi(elecstate::ElecState* pelec, - const Charge& chr, UnitCell& ucell, - std::vector> pphi_, + void propagate_psi_RK4(elecstate::ElecState* pelec, + Charge& chr, + UnitCell& ucell, + std::vector>& pphi_, ModulePW::PW_Basis* pw_rho); + + void propagate_psi_RK2(elecstate::ElecState* pelec, + Charge& chr, + UnitCell& ucell, + std::vector>& pphi_, + ModulePW::PW_Basis* pw_rho); + + void renormalize_psi(Charge& chr, ModulePW::PW_Basis* pw_rho, std::vector>& pphi_); private: const double c_tf_ @@ -35,18 +44,18 @@ class Evolve_OFDFT * 2; // 10/3*(3*pi^2)^{2/3}, multiply by 2 to convert unit from Hartree to Ry, finally in Ry*Bohr^(-2) void cal_Hpsi(elecstate::ElecState* pelec, - const Charge& chr, + Charge& chr, UnitCell& ucell, - std::vector> psi_, + std::vector>& psi_, ModulePW::PW_Basis* pw_rho, - std::vector> Hpsi); + 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, + 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_, + std::vector>& Hpsi); // -1/2 \nabla^2 \phi + void cal_CD_potential(std::vector>& psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot, double mCD_para); From 9238edbb9833967aef5fb68395e17e54604d8ac6 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Mon, 12 Jan 2026 18:52:13 +0800 Subject: [PATCH 22/25] solve conflicts --- .../source_pw/module_ofdft/evolve_ofdft.cpp | 40 ------------------- 1 file changed, 40 deletions(-) diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index ca27d61ff2..6b88825124 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -26,17 +26,10 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec, this->renormalize_psi(chr, pw_rho, psi_); pelec->pot->update_from_charge(&chr, &ucell); // Hartree + XC + external -<<<<<<< HEAD this->cal_tf_potential(chr.rho, pw_rho, pelec->pot->get_effective_v()); // TF potential if (PARAM.inp.of_cd) { this->cal_CD_potential(psi_, pw_rho, pelec->pot->get_effective_v(), PARAM.inp.of_mCD_alpha); // CD potential -======= - this->cal_tf_potential(chr.rho, pw_rho, pelec->pot->get_eff_v()); // TF potential - if (PARAM.inp.of_cd) - { - this->cal_CD_potential(psi_, pw_rho, pelec->pot->get_eff_v(), PARAM.inp.of_mCD_alpha); // CD potential ->>>>>>> cd4d234bc2ce329cfc85091fcdad56a4acafae94 } #ifdef _OPENMP @@ -179,7 +172,6 @@ void Evolve_OFDFT::cal_CD_potential(std::vector>& psi_, #endif for (int is = 0; is < PARAM.inp.nspin; ++is) { -<<<<<<< HEAD std::vector> recipCurrent_x(pw_rho->npw); std::vector> recipCurrent_y(pw_rho->npw); std::vector> recipCurrent_z(pw_rho->npw); @@ -189,26 +181,11 @@ void Evolve_OFDFT::cal_CD_potential(std::vector>& psi_, std::vector> rCurrent_z(pw_rho->nrxx); std::vector> kF_r(pw_rho->nrxx); std::vector> rCDPotential(pw_rho->nrxx); -======= - std::complex *recipCurrent_x=new std::complex[pw_rho->npw]; - std::complex *recipCurrent_y=new std::complex[pw_rho->npw]; - std::complex *recipCurrent_z=new std::complex[pw_rho->npw]; - std::complex *recipCDPotential=new std::complex[pw_rho->npw]; - std::complex *rCurrent_x=new std::complex[pw_rho->nrxx]; - std::complex *rCurrent_y=new std::complex[pw_rho->nrxx]; - std::complex *rCurrent_z=new std::complex[pw_rho->nrxx]; - std::complex *kF_r=new std::complex[pw_rho->nrxx]; - std::complex *rCDPotential=new std::complex[pw_rho->nrxx]; ->>>>>>> cd4d234bc2ce329cfc85091fcdad56a4acafae94 recipPhi[is] = new std::complex[pw_rho->npw]; for (int ir = 0; ir < pw_rho->nrxx; ++ir) { -<<<<<<< HEAD kF_r[ir]=std::pow(3*std::pow(ModuleBase::PI*std::abs(rPhi[is][ir]),2),1.0/3.0); -======= - kF_r[ir]=std::pow(3*std::pow(ModuleBase::PI*std::abs(rPhi[is][ir]),2),1/3); ->>>>>>> cd4d234bc2ce329cfc85091fcdad56a4acafae94 } pw_rho->real2recip(rPhi[is], recipPhi[is]); @@ -218,15 +195,9 @@ void Evolve_OFDFT::cal_CD_potential(std::vector>& psi_, recipCurrent_y[ik]=imag*pw_rho->gcar[ik].y*recipPhi[is][ik]* pw_rho->tpiba; recipCurrent_z[ik]=imag*pw_rho->gcar[ik].z*recipPhi[is][ik]* pw_rho->tpiba; } -<<<<<<< HEAD pw_rho->recip2real(recipCurrent_x.data(),rCurrent_x.data()); pw_rho->recip2real(recipCurrent_y.data(),rCurrent_y.data()); pw_rho->recip2real(recipCurrent_z.data(),rCurrent_z.data()); -======= - pw_rho->recip2real(recipCurrent_x,rCurrent_x); - pw_rho->recip2real(recipCurrent_y,rCurrent_y); - pw_rho->recip2real(recipCurrent_z,rCurrent_z); ->>>>>>> cd4d234bc2ce329cfc85091fcdad56a4acafae94 for (int ir = 0; ir < pw_rho->nrxx; ++ir) { rCurrent_x[ir]=std::imag(rCurrent_x[ir]*std::conj(rPhi[is][ir])); @@ -333,7 +304,6 @@ void Evolve_OFDFT::propagate_psi_RK4(elecstate::ElecState* pelec, for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { -<<<<<<< HEAD K4[is * nrxx + ir]=-0.5*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]); } @@ -410,14 +380,4 @@ void Evolve_OFDFT::propagate_psi_RK2(elecstate::ElecState* pelec, this->renormalize_psi(chr, pw_rho, pphi_); ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagate_psi_RK2"); -======= - 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]); - } - } - - ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); ->>>>>>> cd4d234bc2ce329cfc85091fcdad56a4acafae94 } From 707811ed6526a68b9c0bde5d581e6d3407264ac3 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Mon, 12 Jan 2026 23:01:58 +0800 Subject: [PATCH 23/25] solve conflicts --- source/source_pw/module_ofdft/evolve_ofdft.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index 6b88825124..44c0f8b5c6 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -26,10 +26,10 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec, this->renormalize_psi(chr, pw_rho, psi_); pelec->pot->update_from_charge(&chr, &ucell); // Hartree + XC + external - this->cal_tf_potential(chr.rho, pw_rho, pelec->pot->get_effective_v()); // TF potential + this->cal_tf_potential(chr.rho, pw_rho, pelec->pot->get_eff_v()); // TF potential if (PARAM.inp.of_cd) { - this->cal_CD_potential(psi_, pw_rho, pelec->pot->get_effective_v(), PARAM.inp.of_mCD_alpha); // CD potential + this->cal_CD_potential(psi_, pw_rho, pelec->pot->get_eff_v(), PARAM.inp.of_mCD_alpha); // CD potential } #ifdef _OPENMP From ac336aecbbc19912a7bb03973b69e3f29dc47613 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Tue, 13 Jan 2026 01:00:03 +0800 Subject: [PATCH 24/25] update integration of tdofdft --- tests/07_OFDFT/30_TDOFDFT_Al/result.ref | 10 +++++----- tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref | 10 +++++----- tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref | 10 +++++----- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/tests/07_OFDFT/30_TDOFDFT_Al/result.ref b/tests/07_OFDFT/30_TDOFDFT_Al/result.ref index 404887187c..b53afa1a10 100644 --- a/tests/07_OFDFT/30_TDOFDFT_Al/result.ref +++ b/tests/07_OFDFT/30_TDOFDFT_Al/result.ref @@ -1,5 +1,5 @@ -etotref -100.6027239728818614 -etotperatomref -50.3013619864 -totalforceref 5.786358 -totalstressref 11367.508779 -totaltimeref 0.10 +etotref 226.9979568201541724 +etotperatomref 113.4989784101 +totalforceref 562.524006 +totalstressref 194931.772028 +totaltimeref 0.45 diff --git a/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref b/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref index b94d694de2..286a081395 100644 --- a/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref +++ b/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref @@ -1,5 +1,5 @@ -etotref -100.6027239728818614 -etotperatomref -50.3013619864 -totalforceref 5.786358 -totalstressref 11367.508779 -totaltimeref 0.13 +etotref 225.6175831066692297 +etotperatomref 112.8087915533 +totalforceref 488.530072 +totalstressref 185999.479726 +totaltimeref 0.09 diff --git a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref index 3b61275a93..a6e3cab7f3 100644 --- a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref +++ b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref @@ -1,5 +1,5 @@ -etotref -100.6027239728818614 -etotperatomref -50.3013619864 -totalforceref 5.786358 -totalstressref 11367.508779 -totaltimeref 0.11 +etotref 223.1193797466232809 +etotperatomref 111.5596898733 +totalforceref 386.437466 +totalstressref 181323.311310 +totaltimeref 0.08 From 7ef3bc2d190d4a27514189a7d5e62ec12dd89c5a Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Tue, 13 Jan 2026 09:48:20 +0800 Subject: [PATCH 25/25] add threshold for TDOFDFT integration test --- tests/07_OFDFT/30_TDOFDFT_Al/threshold | 4 ++++ tests/07_OFDFT/31_TDOFDFT_Al_CD/threshold | 4 ++++ tests/07_OFDFT/32_TDOFDFT_Al_mCD/threshold | 4 ++++ 3 files changed, 12 insertions(+) create mode 100644 tests/07_OFDFT/30_TDOFDFT_Al/threshold create mode 100644 tests/07_OFDFT/31_TDOFDFT_Al_CD/threshold create mode 100644 tests/07_OFDFT/32_TDOFDFT_Al_mCD/threshold diff --git a/tests/07_OFDFT/30_TDOFDFT_Al/threshold b/tests/07_OFDFT/30_TDOFDFT_Al/threshold new file mode 100644 index 0000000000..a137b5cae7 --- /dev/null +++ b/tests/07_OFDFT/30_TDOFDFT_Al/threshold @@ -0,0 +1,4 @@ +threshold 0.0001 +force_threshold 0.001 +stress_threshold 0.01 +fatal_threshold 1 diff --git a/tests/07_OFDFT/31_TDOFDFT_Al_CD/threshold b/tests/07_OFDFT/31_TDOFDFT_Al_CD/threshold new file mode 100644 index 0000000000..a137b5cae7 --- /dev/null +++ b/tests/07_OFDFT/31_TDOFDFT_Al_CD/threshold @@ -0,0 +1,4 @@ +threshold 0.0001 +force_threshold 0.001 +stress_threshold 0.01 +fatal_threshold 1 diff --git a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/threshold b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/threshold new file mode 100644 index 0000000000..a137b5cae7 --- /dev/null +++ b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/threshold @@ -0,0 +1,4 @@ +threshold 0.0001 +force_threshold 0.001 +stress_threshold 0.01 +fatal_threshold 1