Merge branch 'develop' into ZhouXY-PKU-pexsi

Author: ZhouXY-PKU

source/source_estate/elecstate_lcao.cpp

@@ -31,11 +31,11 @@ double ElecStateLCAO<std::complex<double>>::get_spin_constrain_energy()
 }
 
 template <>
-void ElecStateLCAO<double>::dm2Rho(std::vector<double*> pexsi_DM, 
+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)
@@ -72,12 +72,12 @@ void ElecStateLCAO<double>::dm2Rho(std::vector<double*> pexsi_DM,
 
     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,
+void ElecStateLCAO<std::complex<double>>::dm2rho(std::vector<std::complex<double>*> pexsi_DM,
 		std::vector<std::complex<double>*> pexsi_EDM,
 		DensityMatrix<std::complex<double>, double>* dm)
 {

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,