center2 terms

source/Makefile.Objects

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

source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt

@@ -15,6 +15,7 @@ if(ENABLE_LCAO)
         operator_lcao/td_nonlocal_lcao.cpp
         operator_lcao/sc_lambda_lcao.cpp
         operator_lcao/dftu_lcao.cpp
+        pulay_force_stress/pulay_force_stress_center2.cpp
         FORCE_STRESS.cpp
         FORCE_gamma.cpp
         FORCE_k.cpp

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp

@@ -214,12 +214,20 @@ 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,
+    // 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);
+    PulayForceStress::cal_pulay_fs(foverlap, soverlap,
         this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra), 
-        psi, pv, pelec, foverlap, soverlap);
+        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 };
+    // this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi);
+    PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress);
 
     this->cal_fvnl_dbeta(dm,
                          pv,

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp

@@ -310,13 +310,20 @@ 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, 
+    // 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);
+    PulayForceStress::cal_pulay_fs(foverlap, soverlap,
         this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra), 
-        psi, pv, pelec, foverlap, soverlap, kv, 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
+    // this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi, ra);
+    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);

source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp

@@ -97,6 +97,6 @@ elecstate::DensityMatrix<std::complex<double>, double> Force_LCAO<std::complex<d
     // cal_dm_2d
     edm.init_DMR(ra, &ucell);
     edm.cal_DMR();
-    edm.sum_DMR_spin();
+    // edm.sum_DMR_spin();
     return edm;
 }

source/module_hamilt_lcao/hamilt_lcaodft/force_stress_arrays.h

@@ -60,6 +60,10 @@ class ForceStressArrays
     double* DSloc_23;
     double* DSloc_33;
 
+    //缺少DSlocR_11, DSlocR_12, DSlocR_13, DSlocR_22, DSlocR_23, DSlocR_33
+    //DH_r更省内存
+    // 数组接口需要重构
+
     double* DHloc_fixed_11;
     double* DHloc_fixed_12;
     double* DHloc_fixed_13;

source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress/pulay_force_stress.h

@@ -14,21 +14,45 @@ 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
+/// calculate the abstract formulas: 
+/// $Tr[D*dH/dx]$ (force) and $1/V Tr[D*(dH/dx_a*x_b)]$ (stress)
+/// where D can be any (energy) density matrix
+/// and H can be any operator
 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
+    /// for 2-center-integration terms, provided force and stress derivatives
+    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 Parallel_Orbitals& pv,  ///< [in] parallel orbitals
+        const double* (&dHSx)[3],  ///< [in] dHSx x, y, z, for force
+        const double* (&dHSxy)[6],  ///< [in] dHSxy 11, 12, 13, 22, 23, 33, for stress
+        const bool& isforce,
+        const bool& isstress,
+        Record_adj* ra = nullptr,
+        const double& factor_force = 1.0,
+        const double& factor_stress = 1.0);
+
+    /// for 2-center-integration terms, provided force derivatives and coordinate difference
+    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 Parallel_Orbitals& pv,  ///< [in] parallel orbitals
+        const double* (&dHSx)[3],  ///< [in] dHSx x, y, z, for force and stress
+        const double* dtau,  ///< [in] dr x, y, z, for stress
+        const bool& isforce,
+        const bool& isstress,
+        Record_adj* ra = nullptr,
+        const double& factor_force = 1.0,
+        const double& factor_stress = 1.0);
+
+    /// for grid-integration terms
     template<typename TK, typename TR>
     void cal_pulay_fs(
         ModuleBase::matrix& f,  ///< [out] force
@@ -39,8 +63,7 @@ namespace PulayForceStress
         typename TGint<TK>::type& gint, ///< [in] Gint object
         const bool& isforce,
         const bool& isstress,
-        const bool& set_dmr_gint = true
-    );
+        const bool& set_dmr_gint = true);
 }
-
-#include "pulay_force_stress.hpp"
+#include "pulay_force_stress_center2_template.hpp"
+#include "pulay_force_stress_gint.hpp"