[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_phi.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,121 @@
+#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_phi=new EVOLVE_PHI();
+}
+
+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;
+}
+
+
+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<double>*[PARAM.inp.nspin];
+
+        for (int is = 0; is < PARAM.inp.nspin; ++is)
+        {
+            this->pphi_td[is]= new std::complex<double>[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 < pw_rho->nrxx; ++ir)
+            {
+                pphi_td[is][ir]=pphi_[is][ir];
+            }
+        }
+    }
+    else
+    {
+        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)
+            {
+                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_phi.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<double>** pphi_td = nullptr;                     // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi().
+    EVOLVE_PHI* evolve_phi=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_phi.cpp
 )
 
 add_library(