Fix DM in PEXSI (#6689)

* 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

* fix dm bug in PEXSI

* fix PEXSI

* Update edm.cpp

* fix edm

* Update hsolver_lcao.cpp

* Update hsolver_lcao.cpp

---------

Co-authored-by: Levi Zhou <[email protected]>

Author: mohanchen

source/source_estate/elecstate_lcao.cpp

@@ -1,6 +1,5 @@
-#include "elecstate_lcao.h"
-
-#include "cal_dm.h"
+#include "source_estate/elecstate_lcao.h"
+#include "source_estate/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"
@@ -31,25 +30,28 @@ double ElecStateLCAO<std::complex<double>>::get_spin_constrain_energy()
     return sc.cal_escon();
 }
 
-#ifdef __PEXSI
 template <>
-void ElecStateLCAO<double>::dm2Rho(std::vector<double*> pexsi_DM, std::vector<double*> pexsi_EDM)
+void ElecStateLCAO<double>::dm2rho(std::vector<double*> pexsi_DM, 
+		std::vector<double*> pexsi_EDM,
+		DensityMatrix<double, double>* dm)
 {
-    ModuleBase::timer::tick("ElecStateLCAO", "dm2Rho");
+    ModuleBase::timer::tick("ElecStateLCAO", "dm2rho");
 
     int nspin = PARAM.inp.nspin;
     if (PARAM.inp.nspin == 4)
     {
         nspin = 1;
     }
 
-    this->get_DM()->pexsi_EDM = pexsi_EDM;
+#ifdef __PEXSI
+    dm->pexsi_EDM = pexsi_EDM;
+#endif
 
     for (int is = 0; is < nspin; is++)
     {
-        this->DM->set_DMK_pointer(is, pexsi_DM[is]);
+        dm->set_DMK_pointer(is, pexsi_DM[is]);
     }
-    DM->cal_DMR();
+    dm->cal_DMR();
 
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {
@@ -58,30 +60,30 @@ void ElecStateLCAO<double>::dm2Rho(std::vector<double*> pexsi_DM, std::vector<do
     }
 
     ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
-    ModuleGint::cal_gint_rho(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho);
+    ModuleGint::cal_gint_rho(dm->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho);
     if (XC_Functional::get_ked_flag())
     {
         for (int is = 0; is < PARAM.inp.nspin; is++)
         {
             ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[0], this->charge->nrxx);
         }
-        ModuleGint::cal_gint_tau(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->kin_r);
+        ModuleGint::cal_gint_tau(dm->get_DMR_vector(), PARAM.inp.nspin, this->charge->kin_r);
     }
 
     this->charge->renormalize_rho();
 
-    ModuleBase::timer::tick("ElecStateLCAO", "dm2Rho");
+    ModuleBase::timer::tick("ElecStateLCAO", "dm2rho");
     return;
 }
 
 template <>
 void ElecStateLCAO<std::complex<double>>::dm2rho(std::vector<std::complex<double>*> pexsi_DM,
-                                                  std::vector<std::complex<double>*> pexsi_EDM)
+		std::vector<std::complex<double>*> pexsi_EDM,
+		DensityMatrix<std::complex<double>, double>* dm)
 {
     ModuleBase::WARNING_QUIT("ElecStateLCAO", "pexsi is not completed for multi-k case");
 }
 
-#endif
 
 template class ElecStateLCAO<double>;               // Gamma_only case
 template class ElecStateLCAO<std::complex<double>>; // multi-k case

source/source_estate/elecstate_lcao.h

@@ -37,7 +37,6 @@ class ElecStateLCAO : public ElecState
 
     double get_spin_constrain_energy() override;
 
-#ifdef __PEXSI
     // use for pexsi
 
     /**
@@ -46,8 +45,9 @@ 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);
-#endif
+	void dm2rho(std::vector<TK*> pexsi_DM, 
+			std::vector<TK*> pexsi_EDM, 
+			DensityMatrix<TK, double>* dm);
 
 };
 

source/source_hsolver/hsolver_lcao.cpp

@@ -112,7 +112,7 @@ void HSolverLCAO<TK, Device>::solve(hamilt::Hamilt<TK>* pHamilt,
     else if (this->method == "pexsi")
     {
 #ifdef __PEXSI // other purification methods should follow this routine
-        DiagoPexsi<T> pe(ParaV);
+        DiagoPexsi<TK> pe(ParaV);
         for (int ik = 0; ik < psi.get_nk(); ++ik)
         {
             /// update H(k) for each k point
@@ -121,10 +121,10 @@ void HSolverLCAO<TK, Device>::solve(hamilt::Hamilt<TK>* pHamilt,
             // solve eigenvector and eigenvalue for H(k)
             pe.diag(pHamilt, psi, nullptr);
         }
-        auto _pes = dynamic_cast<elecstate::ElecStateLCAO<T>*>(pes);
+        auto _pes = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(pes);
         pes->f_en.eband = pe.totalFreeEnergy;
         // maybe eferm could be dealt with in the future
-        _pes->dm2rho(dm, pe.EDM);
+        _pes->dm2rho(pe.DM, pe.EDM, &dm);
 #endif
     }
 

source/source_lcao/edm.cpp

@@ -31,10 +31,10 @@ elecstate::DensityMatrix<double, double> Force_LCAO<double>::cal_edm(const elecs
 #ifdef __PEXSI
     if (PARAM.inp.ks_solver == "pexsi")
     {
-        auto pes = dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec);
+        // auto pes = dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec);
         for (int ik = 0; ik < nspin; ik++)
         {
-            edm.set_DMK_pointer(ik, pes->get_DM()->pexsi_EDM[ik]);
+            edm.set_DMK_pointer(ik, dm.pexsi_EDM[ik]);
         }
 
     }
@@ -47,7 +47,8 @@ elecstate::DensityMatrix<double, double> Force_LCAO<double>::cal_edm(const elecs
 }
 
 template<>
-elecstate::DensityMatrix<std::complex<double>, double> Force_LCAO<std::complex<double>>::cal_edm(const elecstate::ElecState* pelec,
+elecstate::DensityMatrix<std::complex<double>, double> Force_LCAO<std::complex<double>>::cal_edm(
+    const elecstate::ElecState* pelec,
     const psi::Psi<std::complex<double>>& psi,
     const elecstate::DensityMatrix<std::complex<double>, double>& dm,
     const K_Vectors& kv,