fix bug in ELF

mohanchen · mohanchen · commit 8e5677f79bc7 · 2025-10-24T10:28:21.000+08:00
diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -252,7 +252,7 @@ OBJS_ELECSTAT=elecstate.o\
     update_pot.o\
 
 OBJS_ELECSTAT_LCAO=elecstate_lcao.o\
-      setup_dm.o\
+      init_dm.o\
       density_matrix.o\
       density_matrix_io.o\
       cal_dm_psi.o\
diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp
@@ -14,7 +14,7 @@
 #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()
@@ -383,7 +383,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
                   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
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_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt
@@ -46,7 +46,7 @@ if(ENABLE_LCAO)
   list(APPEND objects
       elecstate_lcao.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_lcao.cpp b/source/source_estate/elecstate_lcao.cpp
@@ -81,7 +81,7 @@ void ElecStateLCAO<double>::dm2Rho(std::vector<double*> pexsi_DM, std::vector<do
 }
 
 template <>
-void ElecStateLCAO<std::complex<double>>::dm2Rho(std::vector<std::complex<double>*> pexsi_DM,
+void ElecStateLCAO<std::complex<double>>::dm2rho(std::vector<std::complex<double>*> pexsi_DM,
                                                   std::vector<std::complex<double>*> pexsi_EDM)
 {
     ModuleBase::WARNING_QUIT("ElecStateLCAO", "pexsi is not completed for multi-k case");
diff --git a/source/source_estate/elecstate_lcao.h b/source/source_estate/elecstate_lcao.h
@@ -56,7 +56,7 @@ class ElecStateLCAO : public ElecState
      * @param pexsi_EDM: pointers of energy-weighed density matrix (EDMK) calculated by pexsi, needed by MD, will be
      * stored in DensityMatrix::pexsi_EDM
      */
-    void dm2Rho(std::vector<TK*> pexsi_DM, std::vector<TK*> pexsi_EDM);
+    void dm2rho(std::vector<TK*> pexsi_DM, std::vector<TK*> pexsi_EDM);
 #endif
 
     DensityMatrix<TK, double>* DM = nullptr;
diff --git a/source/source_estate/module_dm/init_dm.cpp b/source/source_estate/module_dm/init_dm.cpp
@@ -1,17 +1,17 @@
-#include "source_estate/module_dm/setup_dm.h"
+#include "source_estate/module_dm/init_dm.h"
 #include "source_estate/module_dm/cal_dm_psi.h"
 #include "source_estate/elecstate_tools.h"
 #include "source_estate/cal_ux.h"
 
 template <typename TK>
-void elecstate::setup_dm(UnitCell& ucell,
+void elecstate::init_dm(UnitCell& ucell,
 		elecstate::ElecStateLCAO<TK>* pelec,
         psi::Psi<TK>* psi,
 		Charge &chr,
         const int iter,
         const int exx_two_level_step)
 {
-    ModuleBase::TITLE("elecstate", "setup_dm");
+    ModuleBase::TITLE("elecstate", "init_dm");
 
 	if (iter == 1 && exx_two_level_step == 0)
 	{
@@ -49,14 +49,14 @@ void elecstate::setup_dm(UnitCell& ucell,
 }
 
 
-template void elecstate::setup_dm<double>(UnitCell& ucell,
+template void elecstate::init_dm<double>(UnitCell& ucell,
 		elecstate::ElecStateLCAO<double>* pelec,
         psi::Psi<double>* psi,
 		Charge &chr,
         const int iter,
         const int exx_two_level_step);
 
-template void elecstate::setup_dm<std::complex<double>>(UnitCell& ucell,
+template void elecstate::init_dm<std::complex<double>>(UnitCell& ucell,
 		elecstate::ElecStateLCAO<std::complex<double>>* pelec,
         psi::Psi<std::complex<double>>* psi,
 		Charge &chr,
diff --git a/source/source_estate/module_dm/init_dm.h b/source/source_estate/module_dm/init_dm.h
@@ -1,5 +1,5 @@
-#ifndef SETUP_DM_H
-#define SETUP_DM_H
+#ifndef INIT_DM_H
+#define INIT_DM_H
 
 #include "source_cell/unitcell.h" // use unitcell
 #include "source_estate/elecstate_lcao.h"// use ElecStateLCAO
@@ -10,7 +10,7 @@ namespace elecstate
 {
 
 template <typename TK> 
-void setup_dm(UnitCell& ucell,
+void init_dm(UnitCell& ucell,
 		ElecStateLCAO<TK>* pelec,
         psi::Psi<TK>* psi,
 		Charge &chr,
diff --git a/source/source_hsolver/hsolver_lcao.cpp b/source/source_hsolver/hsolver_lcao.cpp
@@ -124,7 +124,7 @@ void HSolverLCAO<T, Device>::solve(hamilt::Hamilt<T>* pHamilt,
         auto _pes = dynamic_cast<elecstate::ElecStateLCAO<T>*>(pes);
         pes->f_en.eband = pe.totalFreeEnergy;
         // maybe eferm could be dealt with in the future
-        _pes->dmToRho(pe.DM, pe.EDM);
+        _pes->dm2rho(pe.DM, pe.EDM);
 #endif
     }
 
diff --git a/source/source_io/ctrl_scf_lcao.cpp b/source/source_io/ctrl_scf_lcao.cpp
@@ -29,6 +29,7 @@
 #endif
 #include "source_lcao/module_rdmft/rdmft.h" // use RDMFT codes
 #include "source_io/to_qo.h" // use toQO
+#include "source_lcao/rho_tau_lcao.h" // mohan add 2025-10-24
 
 
 template <typename TK, typename TR>
@@ -368,10 +369,16 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
 		std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wave function to Wannier90");
 	}
 
+    // 14) calculate the kinetic energy density tau
+    // mohan add 2025-10-24
+    if (inp.out_elf[0] > 0)
+	{
+		LCAO_domain::dm2tau(pelec->DM->get_DMR_vector(), inp.nspin, pelec->charge);
+	}
 
 #ifdef __EXX
     //------------------------------------------------------------------
-    //! 14) Output Hexx matrix in LCAO basis
+    //! 15) Output Hexx matrix in LCAO basis
     // (see `out_chg` in docs/advanced/input_files/input-main.md)
     //------------------------------------------------------------------
     if (inp.out_chg[0])
@@ -392,7 +399,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
     }
 
     //------------------------------------------------------------------
-    //! 15) Write RPA information in LCAO basis
+    //! 16) Write RPA information in LCAO basis
     //------------------------------------------------------------------
     if (inp.rpa)
     {
@@ -408,7 +415,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
 #endif
 
     //------------------------------------------------------------------
-    //! 16) Perform RDMFT calculations, added by jghan, 2024-10-17
+    //! 17) Perform RDMFT calculations, added by jghan, 2024-10-17
     //------------------------------------------------------------------
     if (inp.rdmft == true)
     {
diff --git a/source/source_lcao/rho_tau_lcao.cpp b/source/source_lcao/rho_tau_lcao.cpp
@@ -19,7 +19,10 @@ void LCAO_domain::dm2rho(std::vector<hamilt::HContainer<double>*> &dmr,
     chr->renormalize_rho();
 
     // should be moved somewhere else, mohan 20251024
-    dm2tau(dmr, nspin, chr);
+	if (XC_Functional::get_ked_flag())
+	{
+		dm2tau(dmr, nspin, chr);
+	}
 
     // symmetrize of charge density should be here, mohan 20251023
 
@@ -35,14 +38,11 @@ void LCAO_domain::dm2tau(std::vector<hamilt::HContainer<double>*> &dmr,
     ModuleBase::TITLE("LCAO_domain", "dm2tau");
     ModuleBase::timer::tick("LCAO_domain", "dm2tau");
 
-    if (XC_Functional::get_ked_flag())
-    {
-        for (int is = 0; is < nspin; is++)
-        {
-            ModuleBase::GlobalFunc::ZEROS(chr->kin_r[is], chr->nrxx);
-        }
-        ModuleGint::cal_gint_tau(dmr, nspin, chr->kin_r);
-    }
+	for (int is = 0; is < nspin; is++)
+	{
+		ModuleBase::GlobalFunc::ZEROS(chr->kin_r[is], chr->nrxx);
+	}
+	ModuleGint::cal_gint_tau(dmr, nspin, chr->kin_r);
 
     ModuleBase::timer::tick("LCAO_domain", "dm2tau");
 }
diff --git a/source/source_lcao/rho_tau_lcao.h b/source/source_lcao/rho_tau_lcao.h
@@ -14,7 +14,6 @@ namespace LCAO_domain
 	void dm2tau(std::vector<hamilt::HContainer<double>*> &dmr,
 			const int nspin,
 			Charge* chr);
-
 }
 
 #endif
diff --git a/source/source_lcao/setup_dm.cpp b/source/source_lcao/setup_dm.cpp
@@ -0,0 +1,28 @@
+#include "source_lcao/setup_dm.h"
+
+#include "cal_dm.h"
+#include "source_base/timer.h"
+#include "source_estate/module_dm/cal_dm_psi.h"
+#include "source_hamilt/module_xc/xc_functional.h"
+#include "source_lcao/module_deltaspin/spin_constrain.h"
+#include "source_pw/module_pwdft/global.h"
+#include "source_io/module_parameter/parameter.h"
+
+#include "source_lcao/module_gint/gint_interface.h"
+
+#include <vector>
+
+namespace LCAO_domain
+{
+
+template <typename TK>
+void Setup_DM<TK>::init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin)
+{
+    const int nspin_dm = nspin == 2 ? 2 : 1;
+    this->dm = new DensityMatrix<TK, double>(paraV, nspin_dm, kv->kvec_d, kv->get_nks() / nspin_dm);
+}
+
+template class Setup_DM<double>;               // Gamma_only case
+template class Setup_DM<std::complex<double>>; // multi-k case
+
+} // namespace elecstate
diff --git a/source/source_lcao/setup_dm.h b/source/source_lcao/setup_dm.h
@@ -0,0 +1,44 @@
+#ifndef SETUP_DM_H
+#define SETUP_DM_H 
+
+#include "source_estate/module_dm/density_matrix.h"
+
+#include <vector>
+
+namespace LCAO_DOMAIN
+{
+template <typename TK>
+class Setup_DM
+{
+	public:
+
+    Setup_DM()
+    {
+    } // will be called by ElecStateLCAO_TDDFT
+
+    ~Setup_DM()
+    {
+        if (this->DM != nullptr)
+        {
+            delete this->dm;
+        }
+    }
+
+    // initial density matrix
+    void init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin);
+
+    DensityMatrix<TK, double>* get_dm() const
+    {
+        return const_cast<DensityMatrix<TK, double>*>(this->dm);
+    }
+
+    private:
+
+    DensityMatrix<TK, double>* dm = nullptr;
+
+};
+
+
+} // namespace elecstate
+
+#endif

Original file line number	Diff line number	Diff line change
`@@ -359,6 +359,7 @@ void ESolver_KS_PW<T, Device>::after_scf(UnitCell& ucell, const int istep, const`
`359`	`359`	`if (PARAM.inp.out_elf[0] > 0)`
`360`	`360`	`{`
`361`	`361`	`this->ESolver_KS<T, Device>::psi = new psi::Psi<T>(this->stp.psi_cpu[0]);`
	`362`	`+ this->pelec->cal_tau(*(this->psi));`
`362`	`363`	`}`
`363`	`364`
`364`	`365`	`// Call 'after_scf' of ESolver_KS`
Original file line number	Diff line number	Diff line change
`@@ -81,7 +81,7 @@ void ElecStateLCAO<double>::dm2Rho(std::vector<double*> pexsi_DM, std::vector<do`
`81`	`81`	`}`
`82`	`82`
`83`	`83`	`template <>`
`84`		`-void ElecStateLCAO<std::complex<double>>::dm2Rho(std::vector<std::complex<double>*> pexsi_DM,`
	`84`	`+void ElecStateLCAO<std::complex<double>>::dm2rho(std::vector<std::complex<double>*> pexsi_DM,`
`85`	`85`	`std::vector<std::complex<double>*> pexsi_EDM)`
`86`	`86`	`{`
`87`	`87`	`ModuleBase::WARNING_QUIT("ElecStateLCAO", "pexsi is not completed for multi-k case");`
Original file line number	Diff line number	Diff line change
`@@ -124,7 +124,7 @@ void HSolverLCAO<T, Device>::solve(hamilt::Hamilt<T>* pHamilt,`
`124`	`124`	`auto _pes = dynamic_cast<elecstate::ElecStateLCAO<T>*>(pes);`
`125`	`125`	`pes->f_en.eband = pe.totalFreeEnergy;`
`126`	`126`	`// maybe eferm could be dealt with in the future`
`127`		`- _pes->dmToRho(pe.DM, pe.EDM);`
	`127`	`+ _pes->dm2rho(pe.DM, pe.EDM);`
`128`	`128`	`#endif`
`129`	`129`	`}`
`130`	`130`