Refactor: Unify the interfaces of Pulay terms of force and stress

source/Makefile.Objects

@@ -568,12 +568,9 @@ OBJS_LCAO=evolve_elec.o\
       FORCE_STRESS.o\
 	  FORCE_gamma.o\
 	  FORCE_k.o\
-	  fvl_dphi_gamma.o\
-	  fvl_dphi_k.o\
-	  fedm_gamma.o\
-	  fedm_k.o\
-	  ftvnl_dphi_gamma.o\
-	  ftvnl_dphi_k.o\
+      stress_tools.o\
+      edm.o\
+      pulay_force_stress_center2.o\
 	  fvnl_dbeta_gamma.o\
 	  fvnl_dbeta_k.o\
 	  grid_init.o\

source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt

@@ -15,15 +15,12 @@ if(ENABLE_LCAO)
         operator_lcao/td_nonlocal_lcao.cpp
         operator_lcao/sc_lambda_lcao.cpp
         operator_lcao/dftu_lcao.cpp
+        pulay_force_stress_center2.cpp
         FORCE_STRESS.cpp
         FORCE_gamma.cpp
         FORCE_k.cpp
-        fvl_dphi_gamma.cpp
-        fvl_dphi_k.cpp
-        fedm_gamma.cpp
-        fedm_k.cpp
-        ftvnl_dphi_gamma.cpp
-        ftvnl_dphi_k.cpp
+        stress_tools.cpp
+        edm.cpp
         fvnl_dbeta_gamma.cpp
         fvnl_dbeta_k.cpp
         grid_init.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,12 +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;
 };
 
-// this namespace used to store global function for some stress operation
-namespace StressTools
-{
-// set upper matrix to whole matrix
-void stress_fill(const double& lat0_, const double& omega_, ModuleBase::matrix& stress_matrix);
-} // namespace StressTools
 #endif

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp

@@ -13,6 +13,7 @@
 #include "module_elecstate/elecstate_lcao.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
 #include "module_io/write_HS.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h"
 
 template <>
 void Force_LCAO<double>::allocate(const Parallel_Orbitals& pv,
@@ -215,10 +216,17 @@ void Force_LCAO<double>::ftable(const bool isforce,
     // allocate DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z
     this->allocate(pv, fsr, two_center_bundle, orb);
 
+    const double* dSx[3] = { fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z };
+    const double* dSxy[6] = { fsr.DSloc_11, fsr.DSloc_12, fsr.DSloc_13, fsr.DSloc_22, fsr.DSloc_23, fsr.DSloc_33 };
     // calculate the force related to 'energy density matrix'.
-    this->cal_fedm(isforce, isstress, fsr, ucell, dm, psi, pv, pelec, foverlap, soverlap);
+    PulayForceStress::cal_pulay_fs(foverlap, soverlap,
+        this->cal_edm(pelec, *psi, *dm, *kv, pv, PARAM.inp.nspin, PARAM.inp.nbands, ucell, *ra),
+        ucell, pv, dSx, dSxy, isforce, isstress);
 
-    this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi);
+    const double* dHx[3] = { fsr.DHloc_fixed_x, fsr.DHloc_fixed_y, fsr.DHloc_fixed_z };
+    const double* dHxy[6] = { fsr.DHloc_fixed_11, fsr.DHloc_fixed_12, fsr.DHloc_fixed_13, fsr.DHloc_fixed_22, fsr.DHloc_fixed_23, fsr.DHloc_fixed_33 };
+    //tvnl_dphi
+    PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress);
 
     this->cal_fvnl_dbeta(dm,
                          pv,
@@ -231,9 +239,9 @@ void Force_LCAO<double>::ftable(const bool isforce,
                          fvnl_dbeta,
                          svnl_dbeta);
 
-    this->cal_fvl_dphi(isforce, isstress, pelec->pot, gint, fvl_dphi, svl_dphi);
+    // vl_dphi
+    PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, pelec->pot, gint, isforce, isstress, false/*reset dm to gint*/);
 
-    // caoyu add for DeePKS
 #ifdef __DEEPKS
     if (PARAM.inp.deepks_scf)
     {
@@ -314,23 +322,3 @@ void Force_LCAO<double>::ftable(const bool isforce,
     ModuleBase::timer::tick("Force_LCAO", "ftable");
     return;
 }
-
-namespace StressTools
-{
-void stress_fill(const double& lat0_, const double& omega_, ModuleBase::matrix& stress_matrix)
-{
-    assert(omega_ > 0.0);
-    double weight = lat0_ / omega_;
-    for (int i = 0; i < 3; ++i)
-    {
-        for (int j = 0; j < 3; ++j)
-        {
-            if (j > i)
-            {
-                stress_matrix(j, i) = stress_matrix(i, j);
-            }
-            stress_matrix(i, j) *= weight;
-        }
-    }
-}
-} // namespace StressTools

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp

@@ -12,6 +12,7 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/write_HS.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h"
 
 #include <map>
 #include <unordered_map>
@@ -312,14 +313,21 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
                    kv->get_nks(),
                    kv->kvec_d);
 
+    const double* dSx[3] = { fsr.DSloc_Rx, fsr.DSloc_Ry, fsr.DSloc_Rz };
     // 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);
+    PulayForceStress::cal_pulay_fs(foverlap, soverlap,
+        this->cal_edm(pelec, *psi, *dm, *kv, pv, PARAM.inp.nspin, PARAM.inp.nbands, ucell, *ra),
+        ucell, pv, dSx, fsr.DH_r, isforce, isstress, ra, -1.0, 1.0);
 
-    this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi, ra);
+    const double* dHx[3] = { fsr.DHloc_fixedR_x, fsr.DHloc_fixedR_y, fsr.DHloc_fixedR_z }; // T+Vnl
+    const double* dHxy[6] = { fsr.stvnl11, fsr.stvnl12, fsr.stvnl13, fsr.stvnl22, fsr.stvnl23, fsr.stvnl33 };   //T
+    // tvnl_dphi
+    PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress, ra, 1.0, -1.0);
 
     // doing on the real space grid.
-    this->cal_fvl_dphi(isforce, isstress, pelec->pot, gint, fvl_dphi, svl_dphi);
+    // vl_dphi
+    PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, pelec->pot, gint, isforce, isstress, false/*reset dm to gint*/);
 
     this->cal_fvnl_dbeta(dm,
                          pv,

source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp

@@ -0,0 +1,103 @@
+#include "FORCE.h"
+#include "module_elecstate/module_dm/cal_dm_psi.h"
+#include "module_base/memory.h"
+#include "module_parameter/parameter.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;
+}