Refactor: merge before_all_runners of tddft and lcao

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\

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

source/module_esolver/esolver_gets.h

@@ -48,7 +48,7 @@ class ESolver_GetS : public ESolver_KS<TK>
 
     TwoCenterBundle two_center_bundle_;
 
-    //   // temporary introduced during removing GlobalC::ORB
+    // temporary introduced during removing GlobalC::ORB
     LCAO_Orbitals orb_;
 };
 } // namespace ModuleESolver

source/module_esolver/esolver_ks_lcao.cpp

@@ -222,7 +222,8 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(const Input_para& inp, UnitCell
 #endif
 
     // 12) set occupations
-    if (PARAM.inp.ocp)
+    // tddft does not need to set occupations in the first scf
+    if (PARAM.inp.ocp && inp.esolver_type != "tddft")
     {
         this->pelec->fixed_weights(PARAM.inp.ocp_kb, PARAM.inp.nbands, PARAM.inp.nelec);
     }

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"
@@ -66,57 +67,11 @@ ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT()
 
 void ESolver_KS_LCAO_TDDFT::before_all_runners(const Input_para& inp, UnitCell& ucell)
 {
-    // 1) run "before_all_runners" in ESolver_KS
-    ESolver_KS::before_all_runners(inp, ucell);
-
-    // 2) initialize the local pseudopotential with plane wave basis set
-    GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, pw_rho);
-
-    // 3) initialize the electronic states for TDDFT
-    if (this->pelec == nullptr)
-    {
-        this->pelec = new elecstate::ElecStateLCAO_TDDFT(&this->chr,
-                                                         &kv,
-                                                         kv.get_nks(),
-                                                         &this->GK, // mohan add 2024-04-01
-                                                         this->pw_rho,
-                                                         pw_big);
-    }
-
-    // 4) read the local orbitals and construct the interpolation tables.
-    // initialize the pv
-    LCAO_domain::init_basis_lcao(this->pv,
-                                 inp.onsite_radius,
-                                 inp.lcao_ecut,
-                                 inp.lcao_dk,
-                                 inp.lcao_dr,
-                                 inp.lcao_rmax,
-                                 ucell,
-                                 two_center_bundle_,
-                                 orb_);
-
-    // 5) allocate H and S matrices according to computational resources
-    LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, this->pv, kv.get_nks(), orb_);
-
-    // 6) initialize Density Matrix
-    dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)
-        ->init_DM(&kv, &this->pv, PARAM.inp.nspin);
-
-    // 8) initialize the charge density
-    this->pelec->charge->allocate(PARAM.inp.nspin);
-    this->pelec->omega = GlobalC::ucell.omega; // this line is very odd.
-
-    // 9) initializee the potential
-    this->pelec->pot = new elecstate::Potential(pw_rhod,
-                                                pw_rho,
-                                                &GlobalC::ucell,
-                                                &(GlobalC::ppcell.vloc),
-                                                &(sf),
-                                                &(pelec->f_en.etxc),
-                                                &(pelec->f_en.vtxc));
+    // 1) run before_all_runners in ESolver_KS_LCAO
+    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)
@@ -134,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)
     {
@@ -153,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
     {
@@ -168,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);
         }
     }
 
@@ -203,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;
@@ -379,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;
@@ -417,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

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