unify gint terms

source/Makefile.Objects

@@ -568,8 +568,7 @@ OBJS_LCAO=evolve_elec.o\
       FORCE_STRESS.o\
 	  FORCE_gamma.o\
 	  FORCE_k.o\
-	  fvl_dphi_gamma.o\
-	  fvl_dphi_k.o\
+      stress_tools.o\
 	  fedm_gamma.o\
 	  fedm_k.o\
 	  ftvnl_dphi_gamma.o\

source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt

@@ -18,8 +18,7 @@ if(ENABLE_LCAO)
         FORCE_STRESS.cpp
         FORCE_gamma.cpp
         FORCE_k.cpp
-        fvl_dphi_gamma.cpp
-        fvl_dphi_k.cpp
+        stress_tools.cpp
         fedm_gamma.cpp
         fedm_k.cpp
         ftvnl_dphi_gamma.cpp

source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h

@@ -138,10 +138,4 @@ class Force_LCAO
                       ModuleBase::matrix& svl_dphi);
 };
 
-// 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/pulay_force_stress.h"
 
 template <>
 void Force_LCAO<double>::allocate(const Parallel_Orbitals& pv,
@@ -231,9 +232,8 @@ 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);
+    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 +314,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/pulay_force_stress.h"
 
 #include <map>
 #include <unordered_map>
@@ -319,7 +320,8 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
     this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi, ra);
 
     // doing on the real space grid.
-    this->cal_fvl_dphi(isforce, isstress, pelec->pot, gint, fvl_dphi, svl_dphi);
+    // this->cal_fvl_dphi(isforce, isstress, pelec->pot, gint, fvl_dphi, svl_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/fedm_gamma.cpp

@@ -7,6 +7,7 @@
 #include "module_psi/psi.h"
 #include "module_elecstate/cal_dm.h"
 #include "module_elecstate/module_dm/cal_dm_psi.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h"
 
 // force due to the overlap matrix.
 // need energy density matrix here.

source/module_hamilt_lcao/hamilt_lcaodft/fedm_k.cpp

@@ -10,6 +10,7 @@
 #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"
 
 #include <map>
 #include <unordered_map>

source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_gamma.cpp

@@ -4,6 +4,7 @@
 #include "module_cell/module_neighbor/sltk_grid_driver.h"
 #include <unordered_map>
 #include "module_base/timer.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h"
 
 template<>
 void Force_LCAO<double>::cal_ftvnl_dphi(

source/module_hamilt_lcao/hamilt_lcaodft/ftvnl_dphi_k.cpp

@@ -14,6 +14,7 @@
 #include "module_elecstate/elecstate_lcao.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/write_HS.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h"
 
 #ifdef __DEEPKS
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h"

source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_gamma.cpp

@@ -3,6 +3,7 @@
 #include "module_parameter/parameter.h"
 #include "module_cell/module_neighbor/sltk_grid_driver.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h"
 
 #include <unordered_map>
 

source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp

@@ -10,6 +10,7 @@
 #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"
 
 #include <map>
 #include <unordered_map>

source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress.h

@@ -0,0 +1,46 @@
+#pragma once
+#include "module_basis/module_nao/two_center_bundle.h"
+#include "module_elecstate/module_dm/density_matrix.h"
+#include "module_hamilt_lcao/module_gint/gint_gamma.h"
+#include "module_hamilt_lcao/module_gint/gint_k.h"
+#include "module_elecstate/potentials/potential_new.h"
+#include "module_cell/unitcell.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h"
+#ifndef TGINT_H
+#define TGINT_H
+template <typename T>
+struct TGint;
+template <> struct TGint<double> { using type = Gint_Gamma; };
+template <> struct TGint<std::complex<double>> { using type = Gint_k; };
+#endif
+
+/// calculate $Tr[D*dH/dx]$ and $1/V Tr[D*(dH/dx_a*x_b)]
+/// where D can be either density matrix or energy density matrix
+namespace PulayForceStress
+{
+    /// for 2-center terms
+    // template<typename TK, typename TR>
+    // void cal_pulay_fs(
+    //     ModuleBase::matrix& f,  ///< [out] force
+    //     ModuleBase::matrix& s,  ///< [out] stress
+    //     const elecstate::DensityMatrix<TK, TR>& dm,  ///< [in] density matrix or energy density matrix
+    //     const UnitCell& ucell,  ///< [in] unit cell
+    //     const ForceStressArrays& fsr,
+    //     const bool& isstress
+    // );
+    /// for grid terms
+    template<typename TK, typename TR>
+    void cal_pulay_fs(
+        ModuleBase::matrix& f,  ///< [out] force
+        ModuleBase::matrix& s,  ///< [out] stress
+        const elecstate::DensityMatrix<TK, TR>& dm,  ///< [in] density matrix or energy density matrix
+        const UnitCell& ucell,  ///< [in] unit cell
+        const elecstate::Potential* pot, ///< [in] potential on grid
+        typename TGint<TK>::type& gint, ///< [in] Gint object
+        const bool& isforce,
+        const bool& isstress,
+        const bool& set_dmr_gint = true
+    );
+}
+
+#include "pulay_force_stress.hpp"