extract edm

source/Makefile.Objects

@@ -569,6 +569,7 @@ OBJS_LCAO=evolve_elec.o\
 	  FORCE_gamma.o\
 	  FORCE_k.o\
       stress_tools.o\
+      edm.o\
 	  fedm_gamma.o\
 	  fedm_k.o\
 	  ftvnl_dphi_gamma.o\

source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt

@@ -19,6 +19,7 @@ if(ENABLE_LCAO)
         FORCE_gamma.cpp
         FORCE_k.cpp
         stress_tools.cpp
+        edm.cpp
         fedm_gamma.cpp
         fedm_k.cpp
         ftvnl_dphi_gamma.cpp

source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h

@@ -94,7 +94,7 @@ class Force_LCAO
                   const bool isstress,
                   ForceStressArrays& fsr,
                   const UnitCell& ucell,
-                  const elecstate::DensityMatrix<T, double>* dm,
+                  const elecstate::DensityMatrix<T, double>& dm,
                   const psi::Psi<T>* psi,
                   const Parallel_Orbitals& pv,
                   const elecstate::ElecState* pelec,
@@ -136,6 +136,16 @@ class Force_LCAO
                       typename TGint<T>::type& gint,
                       ModuleBase::matrix& fvl_dphi,
                       ModuleBase::matrix& svl_dphi);
+
+    elecstate::DensityMatrix<T, double> cal_edm(const elecstate::ElecState* pelec,
+        const psi::Psi<T>& psi,
+        const elecstate::DensityMatrix<T, double>& dm,
+        const K_Vectors& kv,
+        const Parallel_Orbitals& pv,
+        const int& nspin, 
+        const int& nbands,
+        const UnitCell& ucell,
+        Record_adj& ra) const;
 };
 
 #endif

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp

@@ -217,7 +217,9 @@ void Force_LCAO<double>::ftable(const bool isforce,
     this->allocate(pv, fsr, two_center_bundle, orb);
 
     // calculate the force related to 'energy density matrix'.
-    this->cal_fedm(isforce, isstress, fsr, ucell, dm, psi, pv, pelec, foverlap, soverlap);
+    this->cal_fedm(isforce, isstress, fsr, ucell,
+        this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra), 
+        psi, pv, pelec, foverlap, soverlap);
 
     this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi);
 

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp

@@ -315,7 +315,9 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
 
     // calculate the energy density matrix
     // and the force related to overlap matrix and energy density matrix.
-    this->cal_fedm(isforce, isstress, fsr, ucell, dm, psi, pv, pelec, foverlap, soverlap, kv, ra);
+    this->cal_fedm(isforce, isstress, fsr, ucell, 
+        this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra), 
+        psi, pv, pelec, foverlap, soverlap, kv, ra);
 
     this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi, ra);
 

source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp

@@ -0,0 +1,102 @@
+#include "FORCE.h"
+#include "module_elecstate/module_dm/cal_dm_psi.h"
+#include "module_base/memory.h"
+template<>
+elecstate::DensityMatrix<double, double> Force_LCAO<double>::cal_edm(const elecstate::ElecState* pelec,
+    const psi::Psi<double>& psi,
+    const elecstate::DensityMatrix<double, double>& dm,
+    const K_Vectors& kv,
+    const Parallel_Orbitals& pv,
+    const int& nspin, 
+    const int& nbands,
+    const UnitCell& ucell,
+    Record_adj& ra) const
+{
+    ModuleBase::matrix wg_ekb;
+    wg_ekb.create(nspin, nbands);
+
+    for(int is=0; is<nspin; is++)
+    {
+        for(int ib=0; ib<nbands; ib++)
+        {
+            wg_ekb(is,ib) = pelec->wg(is,ib) * pelec->ekb(is, ib);
+        }
+    }
+
+    // construct a DensityMatrix for Gamma-Only
+    elecstate::DensityMatrix<double, double> edm(&pv, nspin);
+
+#ifdef __PEXSI
+    if (PARAM.inp.ks_solver == "pexsi")
+    {
+        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]);
+        }
+
+    }
+    else
+#endif
+    {
+        elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi, edm);
+    }
+    return edm;
+}
+
+template<>
+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,
+    const Parallel_Orbitals& pv,
+    const int& nspin, 
+    const int& nbands,
+    const UnitCell& ucell,
+    Record_adj& ra) const
+{
+
+    // construct a DensityMatrix object
+    elecstate::DensityMatrix<std::complex<double>, double> edm(&kv, &pv, nspin);
+
+    //--------------------------------------------
+    // calculate the energy density matrix here.
+    //--------------------------------------------
+
+    ModuleBase::matrix wg_ekb;
+    wg_ekb.create(kv.get_nks(), nbands);
+    ModuleBase::Memory::record("Force::wg_ekb", sizeof(double) * kv.get_nks() * nbands);
+#ifdef _OPENMP
+#pragma omp parallel for collapse(2) schedule(static, 1024)
+#endif
+    for (int ik = 0; ik < kv.get_nks(); ik++)
+    {
+        for (int ib = 0; ib < nbands; ib++)
+        {
+            wg_ekb(ik, ib) = pelec->wg(ik, ib) * pelec->ekb(ik, ib);
+        }
+    }
+
+    // use the original formula (Hamiltonian matrix) to calculate energy density matrix
+    if (dm.EDMK.size())
+    {
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static)
+#endif
+        for (int ik = 0; ik < kv.get_nks(); ++ik)
+        {
+            edm.set_DMK_pointer(ik, dm.EDMK[ik].c);
+        }
+    }
+    else
+    {
+        // cal_dm_psi
+        elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi, edm);
+    }
+
+    // cal_dm_2d
+    edm.init_DMR(ra, &ucell);
+    edm.cal_DMR();
+    edm.sum_DMR_spin();
+    return edm;
+}

source/module_hamilt_lcao/hamilt_lcaodft/fedm_k.cpp

@@ -7,7 +7,6 @@
 #include "module_cell/module_neighbor/sltk_grid_driver.h"
 #include "module_elecstate/cal_dm.h"
 #include "module_elecstate/elecstate_lcao.h"
-#include "module_elecstate/module_dm/cal_dm_psi.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/write_HS.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h"
@@ -30,7 +29,7 @@ void Force_LCAO<std::complex<double>>::cal_fedm(const bool isforce,
                                                 const bool isstress,
                                                 ForceStressArrays& fsr,
                                                 const UnitCell& ucell,
-                                                const elecstate::DensityMatrix<std::complex<double>, double>* dm,
+                                                const elecstate::DensityMatrix<std::complex<double>, double>& dm,
                                                 const psi::Psi<std::complex<double>>* psi,
                                                 const Parallel_Orbitals& pv,
                                                 const elecstate::ElecState* pelec,
@@ -43,6 +42,7 @@ void Force_LCAO<std::complex<double>>::cal_fedm(const bool isforce,
     ModuleBase::timer::tick("Force_LCAO", "cal_fedm");
 
     const int nspin = PARAM.inp.nspin;
+<<<<<<< HEAD
     const int nbands = PARAM.inp.nbands;
 
     // construct a DensityMatrix object
@@ -88,6 +88,8 @@ void Force_LCAO<std::complex<double>>::cal_fedm(const bool isforce,
     edm.init_DMR(*ra, &ucell);
     edm.cal_DMR();
     edm.sum_DMR_spin();
+=======
+>>>>>>> 27869bac3 (extract edm)
     //--------------------------------------------
     // summation \sum_{i,j} E(i,j)*dS(i,j)
     // BEGIN CALCULATION OF FORCE OF EACH ATOM
@@ -147,7 +149,7 @@ void Force_LCAO<std::complex<double>>::cal_fedm(const bool isforce,
                 double Ry = ra->info[iat][cb][1];
                 double Rz = ra->info[iat][cb][2];
                 // get BaseMatrix
-                hamilt::BaseMatrix<double>* tmp_matrix = edm.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz);
+                hamilt::BaseMatrix<double>* tmp_matrix = dm.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz);
                 if (tmp_matrix == nullptr)
                 {
                     continue;