diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 8d4d9933a1..d2b302b84b 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\ @@ -360,6 +361,7 @@ OBJS_HAMILT_OF=kedf_tf.o\ kedf_xwm.o\ kedf_lkt.o\ kedf_manager.o\ + evolve_ofdft.o\ OBJS_HAMILT_LCAO=hamilt_lcao.o\ operator_lcao.o\ 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.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..a978a1267c --- /dev/null +++ b/source/source_esolver/esolver_of_tddft.cpp @@ -0,0 +1,129 @@ +#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->evolve_ofdft=new Evolve_OFDFT(); +} + +ESolver_OF_TDDFT::~ESolver_OF_TDDFT() +{ + delete this->evolve_ofdft; +} + + +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->phi_td.resize(PARAM.inp.nspin*this->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); + + // update the rho and phi based on the direction and theta + this->update_rho(); + + this->iter_++; + + 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) + { + phi_td[is*this->pw_rho->nrxx+ir]=pphi_[is][ir]; + } + } + } + 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); +#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(phi_td[is*this->pw_rho->nrxx+ir]); + } + } + conv_esolver=true; + } + + this->after_opt(istep, ucell, conv_esolver); + + ModuleBase::timer::tick("ESolver_OF_TDDFT", "runner"); +} + +} // 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..99a4238c0f --- /dev/null +++ b/source/source_esolver/esolver_of_tddft.h @@ -0,0 +1,23 @@ +#ifndef ESOLVER_OF_TDDFT_H +#define ESOLVER_OF_TDDFT_H + +#include "esolver_of.h" +#include "source_pw/module_ofdft/evolve_ofdft.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::vector> phi_td; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). + Evolve_OFDFT* evolve_ofdft=nullptr; +}; +} // 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"); diff --git a/source/source_pw/module_ofdft/CMakeLists.txt b/source/source_pw/module_ofdft/CMakeLists.txt index d34ac86492..7465b47a98 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_ofdft.cpp ) add_library( diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp new file mode 100644 index 0000000000..4c158998cc --- /dev/null +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -0,0 +1,196 @@ +#include "evolve_ofdft.h" + +#include "source_io/module_parameter/parameter.h" +#include + +#include "source_base/parallel_reduce.h" + +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 * pw_rho->nrxx + ir])*abs(psi_[is * pw_rho->nrxx + ir]); + } + } + + 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 + +#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 * pw_rho->nrxx + ir] = vr_eff[ir]*psi_[is * pw_rho->nrxx + ir]; + } + } + this->cal_vw_potential_phi(psi_, pw_rho, Hpsi); +} + +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.); + } + } + 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) + { + rpot(is, ir) += 5.0 / 3.0 * this->c_tf_ * std::pow(2. * prho[is][ir], 2. / 3.); + } + } + } +} + +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 * 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]; + + 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; + } + 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]; + } + } + +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + delete[] recipPhi[is]; + delete[] rLapPhi[is]; + } + delete[] recipPhi; + delete[] rLapPhi; +} + +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) + { + //recipCurrent = new std::complex[pw_rho->npw]; + //delete[] recipCurrent; + } +} + +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; + 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); + + 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 * nrxx + ir]=-1.0*K1[is * nrxx + ir]*dt*imag; + psi1[is * nrxx + ir]=pphi_[is * nrxx + ir]+0.5*K1[is * nrxx + ir]; + } + } + 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 * nrxx + ir]=-1.0*K2[is * nrxx + ir]*dt*imag; + psi2[is * nrxx + ir]=pphi_[is * nrxx + ir]+0.5*K2[is * nrxx + ir]; + } + } + 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 * nrxx + ir]=-1.0*K3[is * nrxx + ir]*dt*imag; + psi3[is * nrxx + ir]=pphi_[is * nrxx + ir]+K3[is * nrxx + ir]; + } + } + 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 * 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"); +} \ No newline at end of file diff --git a/source/source_pw/module_ofdft/evolve_ofdft.h b/source/source_pw/module_ofdft/evolve_ofdft.h new file mode 100644 index 0000000000..d3cd6d358b --- /dev/null +++ b/source/source_pw/module_ofdft/evolve_ofdft.h @@ -0,0 +1,54 @@ +#ifndef Evolve_OFDFT_H +#define Evolve_OFDFT_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_OFDFT +{ + public: + Evolve_OFDFT() + { + } + ~Evolve_OFDFT() + { + } + 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 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