Refactor: remove elecstate_lcao_tddft

YuLiu98 · YuLiu98 · commit dcb3efd7d232 · 2024-11-19T16:48:09.000+08:00
diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -228,7 +228,6 @@ OBJS_ELECSTAT=elecstate.o\
     pot_xc.o\
 
 OBJS_ELECSTAT_LCAO=elecstate_lcao.o\
-      elecstate_lcao_tddft.o\
       elecstate_lcao_cal_tau.o\
       density_matrix.o\
       density_matrix_io.o\
diff --git a/source/module_elecstate/CMakeLists.txt b/source/module_elecstate/CMakeLists.txt
@@ -36,7 +36,6 @@ list(APPEND objects
 if(ENABLE_LCAO)
   list(APPEND objects
       elecstate_lcao.cpp
-      elecstate_lcao_tddft.cpp
       elecstate_lcao_cal_tau.cpp
       potentials/H_TDDFT_pw.cpp
       module_dm/density_matrix.cpp
diff --git a/source/module_elecstate/elecstate_lcao_tddft.cpp b/source/module_elecstate/elecstate_lcao_tddft.cpp
diff --git a/source/module_elecstate/elecstate_lcao_tddft.h b/source/module_elecstate/elecstate_lcao_tddft.h
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
@@ -56,7 +56,6 @@
 
 //-----HSolver ElecState Hamilt--------
 #include "module_elecstate/elecstate_lcao.h"
-#include "module_elecstate/elecstate_lcao_tddft.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
 #include "module_hsolver/hsolver_lcao.h"
 // function used by deepks
@@ -129,26 +128,14 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(const Input_para& inp, UnitCell
     // autoset nbands in ElecState, it should before basis_init (for Psi 2d division)
     if (this->pelec == nullptr)
     {
-        if (inp.esolver_type == "tddft")
-        {
-            this->pelec = new elecstate::ElecStateLCAO_TDDFT(&(this->chr),
-                                                             &(this->kv),
-                                                             this->kv.get_nks(),
-                                                             &(this->GK), // mohan add 2024-04-01
-                                                             this->pw_rho,
-                                                             this->pw_big);
-        }
-        else
-        {
-            // TK stands for double and complex<double>?
-            this->pelec = new elecstate::ElecStateLCAO<TK>(&(this->chr), // use which parameter?
-                                                           &(this->kv),
-                                                           this->kv.get_nks(),
-                                                           &(this->GG), // mohan add 2024-04-01
-                                                           &(this->GK), // mohan add 2024-04-01
-                                                           this->pw_rho,
-                                                           this->pw_big);
-        }
+        // TK stands for double and complex<double>?
+        this->pelec = new elecstate::ElecStateLCAO<TK>(&(this->chr), // use which parameter?
+                                                       &(this->kv),
+                                                       this->kv.get_nks(),
+                                                       &(this->GG), // mohan add 2024-04-01
+                                                       &(this->GK), // mohan add 2024-04-01
+                                                       this->pw_rho,
+                                                       this->pw_big);
     }
 
     // 3) init LCAO basis
diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp
@@ -13,7 +13,9 @@
 #include "module_base/lapack_connector.h"
 #include "module_base/scalapack_connector.h"
 #include "module_elecstate/module_charge/symmetry_rho.h"
+#include "module_elecstate/module_dm/cal_dm_psi.h"
 #include "module_elecstate/module_dm/cal_edm_tddft.h"
+#include "module_elecstate/module_dm/density_matrix.h"
 #include "module_elecstate/occupy.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need divide_HS_in_frag
 #include "module_hamilt_lcao/module_tddft/evolve_elec.h"
@@ -23,7 +25,6 @@
 
 //-----HSolver ElecState Hamilt--------
 #include "module_elecstate/elecstate_lcao.h"
-#include "module_elecstate/elecstate_lcao_tddft.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
 #include "module_hsolver/hsolver_lcao.h"
 #include "module_parameter/parameter.h"
@@ -70,7 +71,7 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(const Input_para& inp, UnitCell&
     ESolver_KS_LCAO<std::complex<double>, double>::before_all_runners(inp, ucell);
 
     // this line should be optimized
-    this->pelec_td = dynamic_cast<elecstate::ElecStateLCAO_TDDFT*>(this->pelec);
+    // this->pelec = dynamic_cast<elecstate::ElecStateLCAO_TDDFT*>(this->pelec);
 }
 
 void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int iter, const double ethr)
@@ -88,13 +89,13 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int ite
                                                  this->psi_laststep,
                                                  this->Hk_laststep,
                                                  this->Sk_laststep,
-                                                 this->pelec_td->ekb,
+                                                 this->pelec->ekb,
                                                  td_htype,
                                                  PARAM.inp.propagator,
                                                  kv.get_nks());
-            this->pelec_td->psiToRho_td(this->psi[0]);
+            this->weight_dm_rho();
         }
-        this->pelec_td->psiToRho_td(this->psi[0]);
+        this->weight_dm_rho();
     }
     else if (istep >= 2)
     {
@@ -107,11 +108,11 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int ite
                                              this->psi_laststep,
                                              this->Hk_laststep,
                                              this->Sk_laststep,
-                                             this->pelec_td->ekb,
+                                             this->pelec->ekb,
                                              td_htype,
                                              PARAM.inp.propagator,
                                              kv.get_nks());
-        this->pelec_td->psiToRho_td(this->psi[0]);
+        this->weight_dm_rho();
     }
     else
     {
@@ -122,7 +123,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int ite
         {
             bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
             hsolver::HSolverLCAO<std::complex<double>> hsolver_lcao_obj(&this->pv, PARAM.inp.ks_solver);
-            hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec_td, skip_charge);
+            hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, skip_charge);
         }
     }
 
@@ -157,8 +158,7 @@ void ESolver_KS_LCAO_TDDFT::iter_finish(const int istep, int& iter)
             for (int ib = 0; ib < PARAM.inp.nbands; ib++)
             {
                 std::setprecision(6);
-                GlobalV::ofs_running << ik + 1 << "     " << ib + 1 << "      " << this->pelec_td->wg(ik, ib)
-                                     << std::endl;
+                GlobalV::ofs_running << ik + 1 << "     " << ib + 1 << "      " << this->pelec->wg(ik, ib) << std::endl;
             }
         }
         GlobalV::ofs_running << std::endl;
@@ -333,7 +333,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
             for (int ib = 0; ib < PARAM.inp.nbands; ib++)
             {
                 GlobalV::ofs_running << ik + 1 << "     " << ib + 1 << "      "
-                                     << this->pelec_td->ekb(ik, ib) * ModuleBase::Ry_to_eV << std::endl;
+                                     << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV << std::endl;
             }
         }
         GlobalV::ofs_running << std::endl;
@@ -371,4 +371,21 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep)
     ESolver_KS_LCAO<std::complex<double>, double>::after_scf(istep);
 }
 
+void ESolver_KS_LCAO_TDDFT::weight_dm_rho()
+{
+    if (PARAM.inp.ocp == 1)
+    {
+        this->pelec->fixed_weights(PARAM.inp.ocp_kb, PARAM.inp.nbands, PARAM.inp.nelec);
+    }
+    this->pelec->calEBand();
+
+    ModuleBase::GlobalFunc::NOTE("Calculate the density matrix.");
+
+    auto _pes = dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec);
+    elecstate::cal_dm_psi(_pes->DM->get_paraV_pointer(), _pes->wg, this->psi[0], *(_pes->DM));
+    _pes->DM->cal_DMR();
+
+    this->pelec->psiToRho(this->psi[0]);
+}
+
 } // namespace ModuleESolver
diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h
@@ -2,7 +2,6 @@
 #define ESOLVER_KS_LCAO_TDDFT_H
 #include "esolver_ks.h"
 #include "esolver_ks_lcao.h"
-#include "module_elecstate/elecstate_lcao_tddft.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h"
 #include "module_psi/psi.h"
 
@@ -36,10 +35,10 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO<std::complex<double>, doubl
     //! Overlap matrix of last time step
     std::complex<double>** Sk_laststep = nullptr;
 
-    //! Electronic states of rt-TDDFT
-    elecstate::ElecStateLCAO_TDDFT* pelec_td = nullptr;
-
     int td_htype = 1;
+
+  private:
+    void weight_dm_rho();
 };
 
 } // namespace ModuleESolver
