[Feature] TD-OFDFT

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\

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

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 <stdexcept>
@@ -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();

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 ---------------

source/source_esolver/esolver_of_tddft.cpp

@@ -0,0 +1,116 @@
+#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->pphi_td.resize(PARAM.inp.nspin);
+
+        for (int is = 0; is < PARAM.inp.nspin; ++is)
+        {
+            this->pphi_td[is].resize(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);
+
+            std::cout<<"optimize------"<<std::endl;
+
+            // 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);
+        }
+
+        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];
+            }
+        }
+    }
+    else
+    {
+        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 < this->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");
+}
+
+} // namespace ModuleESolver

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<std::vector<std::complex<double>>> pphi_td;                     // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi().
+    Evolve_OFDFT* evolve_ofdft=nullptr;
+};
+} // namespace ModuleESolver
+
+#endif

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<std::string> esolver_types = { "ksdft", "sdft", "ofdft", "tddft", "lj", "dp", "lr", "ks-lr" };
+            const std::vector<std::string> 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");

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(

source/source_pw/module_ofdft/evolve_ofdft.cpp

@@ -0,0 +1,146 @@
+#include "evolve_ofdft.h"
+
+#include "source_io/module_parameter/parameter.h"
+#include <iostream>
+
+#include "source_base/parallel_reduce.h"
+
+void Evolve_OFDFT::get_Hpsi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::vector<std::vector<std::complex<double>>> psi_, ModulePW::PW_Basis* pw_rho, std::vector<std::vector<std::complex<double>>> 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_OFDFT::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_OFDFT::get_vw_potential_phi(std::vector<std::vector<std::complex<double>>> pphi, ModulePW::PW_Basis* pw_rho, std::vector<std::vector<std::complex<double>>> Hpsi)
+{
+    std::complex<double>** rLapPhi = new std::complex<double>*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is) {
+        rLapPhi[is] = new std::complex<double>[pw_rho->nrxx];
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            rLapPhi[is][ir]=pphi[is][ir];
+        }
+    }
+    std::complex<double>** recipPhi = new std::complex<double>*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        recipPhi[is] = new std::complex<double>[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 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_OFDFT::get_CD_potential(std::vector<std::vector<std::complex<double>>> psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot)
+{
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        //recipCurrent = new std::complex<double>[pw_rho->npw];
+        //delete[] recipCurrent;
+    }
+}
+
+void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, const Charge& chr, UnitCell& ucell, std::vector<std::vector<std::complex<double>>> pphi_, ModulePW::PW_Basis* pw_rho)
+{
+    ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi");
+
+    std::complex<double> imag(0.0,1.0);
+    double dt=PARAM.inp.mdp.md_dt;
+    std::vector<std::vector<std::complex<double>>> K1(PARAM.inp.nspin,std::vector<std::complex<double>>(pw_rho->nrxx));
+    std::vector<std::vector<std::complex<double>>> K2(PARAM.inp.nspin,std::vector<std::complex<double>>(pw_rho->nrxx));
+    std::vector<std::vector<std::complex<double>>> K3(PARAM.inp.nspin,std::vector<std::complex<double>>(pw_rho->nrxx));
+    std::vector<std::vector<std::complex<double>>> K4(PARAM.inp.nspin,std::vector<std::complex<double>>(pw_rho->nrxx));
+    std::vector<std::vector<std::complex<double>>> psi1(PARAM.inp.nspin,std::vector<std::complex<double>>(pw_rho->nrxx));
+    std::vector<std::vector<std::complex<double>>> psi2(PARAM.inp.nspin,std::vector<std::complex<double>>(pw_rho->nrxx));
+    std::vector<std::vector<std::complex<double>>> psi3(PARAM.inp.nspin,std::vector<std::complex<double>>(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]);
+        }
+    }
+
+    ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi");
+}