deepmodeling · mohanchen · Oct 26, 2025 · Oct 16, 2025 · Oct 17, 2025 · Oct 17, 2025
diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -228,11 +228,11 @@ OBJS_ELECSTAT=elecstate.o\
     elecstate_energy.o\
     elecstate_exx.o\
     elecstate_print.o\
+    elecstate_tools.o\
     elecstate_pw.o\
     elecstate_pw_sdft.o\
     elecstate_pw_cal_tau.o\
     elecstate_op.o\
-    elecstate_tools.o\
     efield.o\
     gatefield.o\
     potential_new.o\
@@ -252,8 +252,7 @@ OBJS_ELECSTAT=elecstate.o\
     update_pot.o\
 
 OBJS_ELECSTAT_LCAO=elecstate_lcao.o\
-      elecstate_lcao_cal_tau.o\
-      setup_dm.o\
+      init_dm.o\
       density_matrix.o\
       density_matrix_io.o\
       cal_dm_psi.o\
@@ -643,6 +642,7 @@ OBJS_LCAO=evolve_elec.o\
       LCAO_init_basis.o\
       setup_exx.o\
       setup_deepks.o\
+      rho_tau_lcao.o\
       center2_orb.o\
       center2_orb-orb11.o\
       center2_orb-orb21.o\

diff --git a/source/source_esolver/esolver_dm2rho.cpp b/source/source_esolver/esolver_dm2rho.cpp
@@ -10,6 +10,7 @@
 #include "source_io/cube_io.h"
 #include "source_io/io_npz.h"
 #include "source_io/print_info.h"
+#include "source_lcao/rho_tau_lcao.h" // mohan add 2025-10-24
 
 namespace ModuleESolver
 {
@@ -45,9 +46,16 @@ 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 = dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
+    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)));
@@ -59,7 +67,9 @@ void ESolver_DM2rho<TK, TR>::runner(UnitCell& ucell, const int istep)
         ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(2)));
     }
 
-    this->pelec->psiToRho(*this->psi);
+    // it's dangrous 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);
 
     int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1;
 

diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp
@@ -88,7 +88,6 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(UnitCell& ucell, const Input_p
 
     // 10) inititlize the charge density
     this->chr_base.allocate(PARAM.inp.nspin);
-    this->pelec_base->omega = ucell.omega;
 
     // 11) initialize the potential
     if (this->pelec_base->pot == nullptr)
@@ -114,7 +113,6 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
 
     ESolver_KS_LCAO<TK,TR>::before_scf(ucell, istep);
 
-    this->pelec_base->omega = ucell.omega;
     //----------------------------------------------------------
     //! calculate D2 or D3 vdW
     //----------------------------------------------------------
@@ -151,16 +149,10 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
             this->two_center_bundle_,
             this->orb_,
             DM,
-            this->deepks
-#ifdef __EXX
-            ,
-            istep,
-            GlobalC::exx_info.info_ri.real_number ? &this->exx_nao.exd->two_level_step : &this->exx_nao.exc->two_level_step,
-            GlobalC::exx_info.info_ri.real_number ? &this->exx_nao.exd->get_Hexxs() : nullptr,
-            GlobalC::exx_info.info_ri.real_number ? nullptr : &this->exx_nao.exc->get_Hexxs()
-#endif
-        );
-    }
+			this->deepks,
+			istep,
+			this->exx_nao);
+	}
 
     XC_Functional::set_xc_type(PARAM.inp.deepks_out_base);
     this->pelec_base->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);

diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp
@@ -108,8 +108,6 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
         this->locpp.init_vloc(ucell, this->pw_rhod);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
-        this->pelec->omega = ucell.omega;
-
         // perform symmetry analysis
         if (ModuleSymmetry::Symmetry::symm_flag == 1)
         {

diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp
@@ -347,12 +347,14 @@ void ESolver_KS<T, Device>::after_scf(UnitCell& ucell, const int istep, const bo
 {
     ModuleBase::TITLE("ESolver_KS", "after_scf");
 
+/*
     // 1) calculate the kinetic energy density tau
     if (PARAM.inp.out_elf[0] > 0)
     {
         assert(this->psi != nullptr);
         this->pelec->cal_tau(*(this->psi));
     }
+*/
 
     // 2) call after_scf() of ESolver_FP
     ESolver_FP::after_scf(ucell, istep, conv_esolver);

diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp
@@ -1,9 +1,7 @@
 #include "esolver_ks_lcao.h"
 #include "source_estate/elecstate_tools.h"
 #include "source_lcao/module_deltaspin/spin_constrain.h"
-#include "source_io/read_wfc_nao.h"
 #include "source_lcao/hs_matrix_k.hpp" // there may be multiple definitions if using hpp
-#include "source_estate/cal_ux.h"
 #include "source_estate/module_charge/symmetry_rho.h"
 #include "source_lcao/LCAO_domain.h" // need DeePKS_init
 #include "source_lcao/module_dftu/dftu.h"
@@ -16,13 +14,14 @@
 #endif
 #include "source_lcao/module_rdmft/rdmft.h"
 #include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006 
-#include "source_estate/module_dm/setup_dm.h" // setup dm from electronic wave functions
+#include "source_estate/module_dm/init_dm.h" // init dm from electronic wave functions
 #include "source_io/ctrl_runner_lcao.h" // use ctrl_runner_lcao() 
 #include "source_io/ctrl_iter_lcao.h" // use ctrl_iter_lcao() 
 #include "source_io/ctrl_scf_lcao.h" // use ctrl_scf_lcao()
 #include "source_psi/setup_psi.h" // mohan add 20251019
 #include "source_io/read_wfc_nao.h" 
 #include "source_io/print_info.h"
+#include "source_lcao/rho_tau_lcao.h" // mohan add 20251024
 
 namespace ModuleESolver
 {
@@ -110,7 +109,6 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
 
     // 11) init charge density
     this->chr.allocate(inp.nspin);
-    this->pelec->omega = ucell.omega;
 
     // 12) init potentials
     if (this->pelec->pot == nullptr)
@@ -199,15 +197,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
 
         this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
             ucell, this->gd, &this->pv, this->pelec->pot, this->kv,
-            two_center_bundle_, orb_, DM, this->deepks
-#ifdef __EXX
-            ,
-            istep,
-            GlobalC::exx_info.info_ri.real_number ? &this->exx_nao.exd->two_level_step : &this->exx_nao.exc->two_level_step,
-            GlobalC::exx_info.info_ri.real_number ? &this->exx_nao.exd->get_Hexxs() : nullptr,
-            GlobalC::exx_info.info_ri.real_number ? nullptr : &this->exx_nao.exc->get_Hexxs()
-#endif
-        );
+            two_center_bundle_, orb_, DM, this->deepks, istep, exx_nao);
     }
 
     // 9) for each ionic step, the overlap <phi|alpha> must be rebuilt
@@ -224,19 +214,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
 
     // 11) set xc type before the first cal of xc in pelec->init_scf, Peize Lin add 2016-12-03
-#ifdef __EXX
-    if (PARAM.inp.calculation != "nscf")
-    {
-        if (GlobalC::exx_info.info_ri.real_number)
-        {
-            this->exx_nao.exd->exx_beforescf(istep, this->kv, *this->p_chgmix, ucell, orb_);
-        }
-        else
-        {
-            this->exx_nao.exc->exx_beforescf(istep, this->kv, *this->p_chgmix, ucell, orb_);
-        }
-    }
-#endif
+    this->exx_nao.before_scf(ucell, this->kv, orb_, this->p_chgmix, istep, PARAM.inp);
 
     // 12) init_scf, should be before_scf? mohan add 2025-03-10
     this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
@@ -318,10 +296,6 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
     ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_force");
 }
 
-//------------------------------------------------------------------------------
-//! the 7th function of ESolver_KS_LCAO: cal_stress
-//! mohan add 2024-05-11
-//------------------------------------------------------------------------------
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress)
 {
@@ -406,10 +380,11 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
 		{
 			// the following steps are only needed in the first outer exx loop
 			exx_two_level_step
-				= GlobalC::exx_info.info_ri.real_number ? this->exx_nao.exd->two_level_step : this->exx_nao.exc->two_level_step;
+				= GlobalC::exx_info.info_ri.real_number ? 
+                  this->exx_nao.exd->two_level_step : this->exx_nao.exc->two_level_step;
 		}
 #endif
-		elecstate::setup_dm<TK>(ucell, estate, this->psi, this->chr, iter, exx_two_level_step);
+		elecstate::init_dm<TK>(ucell, estate, this->psi, this->chr, iter, exx_two_level_step);
 	}
 
 #ifdef __EXX
@@ -495,7 +470,7 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int
     if (!skip_solve)
     {
         hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
-        hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, skip_charge);
+        hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, this->chr, PARAM.inp.nspin, skip_charge);
     }
 
     // 4) EXX
@@ -563,15 +538,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
     }
 
     // 2) for deepks, calculate delta_e, output labels during electronic steps
-#ifdef __MLALGO
-    if (PARAM.inp.deepks_scf)
-    {
-        this->deepks.ld.dpks_cal_e_delta_band(dm_vec, this->kv.get_nks());
-        DeePKS_domain::update_dmr(this->kv.kvec_d, dm_vec, ucell, orb_, this->pv, this->gd, this->deepks.ld.dm_r);
-        estate->f_en.edeepks_scf = this->deepks.ld.E_delta - this->deepks.ld.e_delta_band;
-        estate->f_en.edeepks_delta = this->deepks.ld.E_delta;
-    }
-#endif
+    this->deepks.delta_e(ucell, this->kv, this->orb_, this->pv, this->gd, dm_vec, this->pelec->f_en, PARAM.inp);
 
     // 3) for delta spin
     if (PARAM.inp.sc_mag_switch)
@@ -617,9 +584,6 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
     ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf");
     ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf");
 
-    //! 1) call after_scf() of ESolver_KS
-    ESolver_KS<TK>::after_scf(ucell, istep, conv_esolver);
-
     auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
     auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
 
@@ -633,6 +597,15 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
         ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_scf","p_hamilt does not exist");
     }
 
+    if (PARAM.inp.out_elf[0] > 0)
+	{
+		LCAO_domain::dm2tau(estate->DM->get_DMR_vector(), PARAM.inp.nspin, estate->charge);
+	}
+
+    //! 1) call after_scf() of ESolver_KS
+    ESolver_KS<TK>::after_scf(ucell, istep, conv_esolver);
+
+
     //! 2) output of lcao every few ionic steps
     ModuleIO::ctrl_scf_lcao<TK, TR>(ucell,
             PARAM.inp, this->kv, estate, this->pv,
@@ -654,7 +627,6 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
     ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf");
 }
 
-
 template class ESolver_KS_LCAO<double, double>;
 template class ESolver_KS_LCAO<std::complex<double>, double>;
 template class ESolver_KS_LCAO<std::complex<double>, std::complex<double>>;

diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp
@@ -31,10 +31,9 @@
 #include "source_lcao/hamilt_lcao.h"
 #include "source_psi/psi.h"
 
-//-----force& stress-------------------
 #include "source_lcao/FORCE_STRESS.h"
+#include "source_lcao/rho_tau_lcao.h" // mohan add 2025-10-24
 
-//---------------------------------------------------
 
 namespace ModuleESolver
 {
@@ -293,7 +292,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell,
         {
             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, skip_charge);
+            hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, this->chr, PARAM.inp.nspin, skip_charge);
         }
     }
 
@@ -574,8 +573,8 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::weight_dm_rho(const UnitCell& ucell)
          _pes->DM->cal_DMR();
     }
 
-    // get the real-space charge density
-    this->pelec->psiToRho(this->psi[0]);
+    // get the real-space charge density, mohan add 2025-10-24
+    LCAO_domain::dm2rho(_pes->DM->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
 }
 
 template class ESolver_KS_LCAO_TDDFT<double, base_device::DEVICE_CPU>;

diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp
@@ -359,6 +359,7 @@ void ESolver_KS_PW<T, Device>::after_scf(UnitCell& ucell, const int istep, const
     if (PARAM.inp.out_elf[0] > 0)
     {
         this->ESolver_KS<T, Device>::psi = new psi::Psi<T>(this->stp.psi_cpu[0]);
+        this->pelec->cal_tau(*(this->psi));
     }
 
     // Call 'after_scf' of ESolver_KS

diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp
@@ -273,7 +273,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
             {
                 // Here we initialize rho to be uniform,
                 // because the rho got by pot.init_pot -> Charge::atomic_rho may contain minus elements.
-                this->chr.rho[is][ibs] = this->nelec_[is] / this->pelec->omega;
+                this->chr.rho[is][ibs] = this->nelec_[is] / ucell.omega;
                 this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]);
             }
         }
@@ -587,7 +587,7 @@ void ESolver_OF::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress)
 {
     ModuleBase::matrix kinetic_stress_;
     kinetic_stress_.create(3, 3);
-    this->kedf_manager_->get_stress(this->pelec->omega,
+    this->kedf_manager_->get_stress(ucell.omega,
                                     this->chr.rho,
                                     this->pphi_,
                                     this->pw_rho,

diff --git a/source/source_esolver/esolver_of_tool.cpp b/source/source_esolver/esolver_of_tool.cpp
@@ -22,7 +22,6 @@ void ESolver_OF::init_elecstate(UnitCell& ucell)
         this->pelec = new elecstate::ElecState((Charge*)(&chr), this->pw_rho, pw_big);
         this->chr.allocate(PARAM.inp.nspin);
     }
-    this->pelec->omega = ucell.omega;
 
     delete this->pelec->pot;
     this->pelec->pot = new elecstate::Potential(this->pw_rhod,

diff --git a/source/source_esolver/lcao_others.cpp b/source/source_esolver/lcao_others.cpp
@@ -183,16 +183,10 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
             two_center_bundle_,
             orb_,
             DM,
-            this->deepks
-#ifdef __EXX
-            ,
-            istep,
-            GlobalC::exx_info.info_ri.real_number ? &this->exx_nao.exd->two_level_step : &this->exx_nao.exc->two_level_step,
-            GlobalC::exx_info.info_ri.real_number ? &this->exx_nao.exd->get_Hexxs() : nullptr,
-            GlobalC::exx_info.info_ri.real_number ? nullptr : &this->exx_nao.exc->get_Hexxs()
-#endif
-        );
-    }
+			this->deepks,
+			istep,
+			this->exx_nao);
+	}
 
 
     // for each ionic step, the overlap <phi|alpha> must be rebuilt

diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt
@@ -45,9 +45,8 @@ list(APPEND objects
 if(ENABLE_LCAO)
   list(APPEND objects
       elecstate_lcao.cpp
-      elecstate_lcao_cal_tau.cpp
       module_pot/H_TDDFT_pw.cpp
-      module_dm/setup_dm.cpp
+      module_dm/init_dm.cpp
       module_dm/density_matrix.cpp
       module_dm/density_matrix_io.cpp
       module_dm/cal_dm_psi.cpp

diff --git a/source/source_estate/elecstate.cpp b/source/source_estate/elecstate.cpp
@@ -54,13 +54,13 @@ void ElecState::init_scf(const int istep,
 }
 
 
-void ElecState::init_ks(Charge* chg_in, // pointer for class Charge
+void ElecState::init_ks(Charge* chr_in, // pointer for class Charge
                         const K_Vectors* klist_in,
                         int nk_in,
                         ModulePW::PW_Basis* rhopw_in,
                         const ModulePW::PW_Basis_Big* bigpw_in)
 {
-    this->charge = chg_in;
+    this->charge = chr_in;
     this->klist = klist_in;
     this->charge->set_rhopw(rhopw_in);
     this->bigpw = bigpw_in;