Remove DM from elecstate (#6675)

* add update_pot in source_estate

* split update_pot in rt-TDDFT into two functions, one is the original update_pot, which is replaced by update_pot function in estate, the other is save2 function used in after_iter

* fix bugs

* update

* update esolver_ks

* small updates

* small bug fixed

* update E_bandgap to E_gap(k)

* delete rho_restart file

* change setup_parameters to print_parameters

* move setup_pw.cpp to setup_pwwfc.cpp in module_pwdft

* update esolver, small things

* delete some old_gint codes in esolver

* fix bugs, update outputs

* fix bug in tests

* change efermi to Efermi in fp_energy.h

* fix bugs

* fix exd and exc in ctrl_runner_lcao

* fix cmake

* fix bug

* delete lcao_after_scf and lcao_before_scf

* modify codes including read_wfc_nao to delete pelec, add setup_psi_lcao codes

* fix bugs

* delete grid integral in FORCE_STRESS

* move read psi to esolver, not in psi

* fix system bug

* fix bug in read_wfc_nao.h

* change PARAM.inp to inp

* fix bugs

* delete omega in pelec, remove some deepks codes to setup_deepks.cpp

* remove exx redundant parameters in Hamilt_LCAO

* Fix circular dependencies of header files.

* fix bugs in MLALGO

* refactor exx lcao

* add rho_tau_lcao files

* fix bug in setup_deepks

* move rho_tau_lcao to source_lcao

* remove old psiToRho function in elecstate_lcao and replace with new rho_tau_lcao

* update estate

* fix bug in ELF

* fix bug

* move DM outside of pelec

* deal with init_dm functions

* remove pelec in dftu, replace with dmk_d(ouble) and dmk_c(omplex)

* fix cal_ldos with DM

* update delta_spin

* update delta spin run_lambda_loop

* the big issue to remove dm may be caused by cal_mw_from_lambda

* remove some DM

* keep updating DM

* refactor some DM in double_xc

* now all functions can be compiled

* fix it

* update again

* add and compile setup_dm

* update fix bugs

* simplify dm

* one more fix

* now I can compile the code

* fix some bugs in EXX version

* fix FORCE_STRESS

* fix 2

* fix 3

* fix 4

* fix 5

* improve esolver_ks_lcao

* improve again

* comment out psi2rho in init_dm

* clean esolver_ks_lcao

* update structure_factor

* update sf

* update

* fix sf

* fix

* fix structure factor

Author: mohanchen

source/Makefile.Objects

@@ -643,6 +643,7 @@ OBJS_LCAO=evolve_elec.o\
       LCAO_init_basis.o\
       setup_exx.o\
       setup_deepks.o\
+      setup_dm.o\
       rho_tau_lcao.o\
       center2_orb.o\
       center2_orb-orb11.o\

source/source_esolver/esolver_dm2rho.cpp

@@ -46,30 +46,22 @@ void ESolver_DM2rho<TK, TR>::runner(UnitCell& ucell, const int istep)
 
     ESolver_KS_LCAO<TK, TR>::before_scf(ucell, istep);
 
-    auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
-
-    if(!estate)
-    {
-        ModuleBase::WARNING_QUIT("ESolver_DM2rho::after_scf","pelec does not exist");
-    }
-
     // file name of DM
     std::string zipname = "output_DM0.npz";
-    elecstate::DensityMatrix<TK, double>* dm = estate->get_DM();
 
     // read DM from file
-    ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(1)));
+    ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(this->dmat.dm->get_DMR_pointer(1)));
 
     // if nspin=2, need extra reading
     if (PARAM.inp.nspin == 2)
     {
         zipname = "output_DM1.npz";
-        ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(2)));
+        ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(this->dmat.dm->get_DMR_pointer(2)));
     }
 
     // it's dangerous to design psiToRho function like this, mohan note 20251024
     // this->pelec->psiToRho(*this->psi);
-    LCAO_domain::dm2rho(estate->DM->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
+    LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
 
     int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1;
 

source/source_esolver/esolver_double_xc.cpp

@@ -84,7 +84,7 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(UnitCell& ucell, const Input_p
     }
 
     // 6) initialize the density matrix
-    dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->init_DM(&this->kv, &(this->pv), PARAM.inp.nspin);
+    this->dmat_base.allocate_dm(&this->kv, &this->pv, PARAM.inp.nspin);
 
     // 10) inititlize the charge density
     this->chr_base.allocate(PARAM.inp.nspin);
@@ -138,8 +138,6 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
     if (this->p_hamilt_base == nullptr)
     {
-        elecstate::DensityMatrix<TK, double>* DM = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->get_DM();
-
         this->p_hamilt_base = new hamilt::HamiltLCAO<TK, TR>(
             ucell,
             this->gd,
@@ -148,7 +146,7 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
             this->kv,
             this->two_center_bundle_,
             this->orb_,
-            DM,
+            this->dmat_base.dm,
 			this->deepks,
 			istep,
 			this->exx_nao);
@@ -159,13 +157,11 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); 
 
     // DMR should be same size with Hamiltonian(R)
-    dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)
-    ->get_DM()
-    ->init_DMR(*(dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt_base)->getHR()));
+    this->dmat_base.dm->init_DMR(*(dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt_base)->getHR()));
 
     if (istep > 0)
     {
-        dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->get_DM()->cal_DMR();
+        this->dmat_base.dm->cal_DMR();
     }
 
     ModuleBase::timer::tick("ESolver_DoubleXC", "before_scf");
@@ -226,23 +222,23 @@ void ESolver_DoubleXC<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int
         std::shared_ptr<LCAO_Deepks<TK>> ld_shared_ptr(&this->deepks.ld, [](LCAO_Deepks<TK>*) {});
         LCAO_Deepks_Interface<TK, TR> deepks_interface(ld_shared_ptr);
 
-        deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
-                                            this->kv.get_nks(),
-                                            ucell.nat,
-                                            PARAM.globalv.nlocal,
-                                            this->pelec->ekb,
-                                            this->kv.kvec_d,
-                                            ucell,
-                                            this->orb_,
-                                            this->gd,
-                                            &(this->pv),
-                                            *(this->psi),
-                                            dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
-                                            p_ham_deepks,
-                                            iter,
-                                            conv_esolver,
-                                            GlobalV::MY_RANK,
-                                            GlobalV::ofs_running);
+		deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
+				this->kv.get_nks(),
+				ucell.nat,
+				PARAM.globalv.nlocal,
+				this->pelec->ekb,
+				this->kv.kvec_d,
+				ucell,
+				this->orb_,
+				this->gd,
+				&(this->pv),
+				*(this->psi),
+				this->dmat.dm,
+				p_ham_deepks,
+				iter,
+				conv_esolver,
+				GlobalV::MY_RANK,
+				GlobalV::ofs_running);
 #endif
 
         // restore to density after charge mixing
@@ -352,9 +348,12 @@ void ESolver_DoubleXC<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int
             auto _pes_lcao = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
             for (int ik = 0; ik < nks; ik++)
             {
-                _pes_lcao_base->get_DM()->set_DMK_pointer(ik, _pes_lcao->get_DM()->get_DMK_pointer(ik));
+// mohan update 2025-11-03
+                this->dmat_base.dm->set_DMK_pointer(ik, this->dmat.dm->get_DMK_pointer(ik));
+//                _pes_lcao_base->get_DM()->set_DMK_pointer(ik, _pes_lcao->get_DM()->get_DMK_pointer(ik));
             }
-            _pes_lcao_base->get_DM()->cal_DMR();
+            this->dmat_base.dm->cal_DMR();
+//            _pes_lcao_base->get_DM()->cal_DMR();
             _pes_lcao_base->ekb = _pes_lcao->ekb;
             _pes_lcao_base->wg = _pes_lcao->wg;          
         }        
@@ -386,6 +385,7 @@ void ESolver_DoubleXC<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
                        this->gd,
                        this->pv,
                        this->pelec_base,
+                       this->dmat_base, // mohan add 2025-11-03
                        this->psi,
                        this->two_center_bundle_,
                        this->orb_,

source/source_esolver/esolver_double_xc.h

@@ -32,6 +32,9 @@ class ESolver_DoubleXC : public ESolver_KS_LCAO<TK, TR>
     //! Electronic states
     elecstate::ElecState* pelec_base = nullptr;
 
+    //! Density Matrix, mohan add 2025-11-03
+    LCAO_domain::Setup_DM<TK> dmat_base;
+
     //! Electorn charge density
     Charge chr_base;
 };

source/source_esolver/esolver_fp.cpp

@@ -54,9 +54,7 @@ void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp)
 
     // write geometry file
     ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif",
-                               ucell,
-                               "# Generated by ABACUS ModuleIO::CifParser",
-                               "data_?");
+      ucell, "# Generated by ABACUS ModuleIO::CifParser", "data_?");
 
     // init charge extrapolation
     this->CE.Init_CE(inp.nspin, ucell.nat, this->pw_rhod->nrxx, inp.chg_extrap);

source/source_esolver/esolver_ks.cpp

@@ -45,7 +45,7 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
 {
     ModuleBase::TITLE("ESolver_KS", "before_all_runners");
 
-    //! 1) init "before_all_runniers" in ESolver_FP
+    //! 1) setup "before_all_runniers" in ESolver_FP
     ESolver_FP::before_all_runners(ucell, inp);
 
     //! 2) setup some parameters
@@ -67,7 +67,7 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
 
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
 
-    //! 4) setup Exc for the first element '0', because all elements have same exc 
+    //! 4) setup Exc for the first element '0' (all elements have same exc) 
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
 
     //! 5) setup the charge mixing parameters
@@ -84,7 +84,7 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
     }
 
-    //! 7) Setup the k points according to symmetry.
+    //! 7) setup k points in the Brillouin zone according to symmetry.
     this->kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
 
@@ -99,19 +99,12 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
 			this->pw_rhod->nplane, this->pw_rhod->nrxx, pw_big->nbz, pw_big->bz);
 
     //! 11) calculate the structure factor
-    this->sf.setup_structure_factor(&ucell, Pgrid, this->pw_rhod);
+    this->sf.setup(&ucell, Pgrid, this->pw_rhod);
 }
 
 template <typename T, typename Device>
 void ESolver_KS<T, Device>::hamilt2rho_single(UnitCell& ucell, const int istep, const int iter, const double ethr)
-{
-    ModuleBase::timer::tick(this->classname, "hamilt2rho_single");
-    // Temporarily, before HSolver is constructed, it should be overrided by
-    // LCAO, PW, SDFT and TDDFT.
-    // After HSolver is constructed, LCAO, PW, SDFT should delete their own
-    // hamilt2rho_single() and use:
-    ModuleBase::timer::tick(this->classname, "hamilt2rho_single");
-}
+{}
 
 template <typename T, typename Device>
 void ESolver_KS<T, Device>::hamilt2rho(UnitCell& ucell, const int istep, const int iter, const double ethr)