refactor LCAO Force

source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h

@@ -0,0 +1,142 @@
+#ifndef FORCE_LCAO_GAMMA_H
+#define FORCE_LCAO_GAMMA_H
+
+#include "module_base/global_function.h"
+#include "module_base/global_variable.h"
+#include "module_base/matrix.h"
+#include "module_elecstate/module_dm/density_matrix.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h"
+#include "module_psi/psi.h"
+#include "module_hamilt_lcao/module_gint/gint_gamma.h"
+#include "module_hamilt_lcao/module_gint/gint_k.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
+
+template<typename T> class Force_Stress_LCAO;
+
+template<typename T>
+class Force_LCAO
+{
+public:
+    friend class Force_Stress_LCAO<T>;
+
+    Force_LCAO() {};
+    ~Force_LCAO() {};
+
+private:
+
+    const Parallel_Orbitals* ParaV;
+
+    elecstate::Potential* pot;
+
+    // orthonormal force + contribution from T and VNL
+    void ftable(const bool isforce,
+        const bool isstress,
+        const psi::Psi<T>* psi,
+        const elecstate::ElecState* pelec,
+        ModuleBase::matrix& foverlap,
+        ModuleBase::matrix& ftvnl_dphi,
+        ModuleBase::matrix& fvnl_dbeta,
+        ModuleBase::matrix& fvl_dphi,
+        ModuleBase::matrix& soverlap,
+        ModuleBase::matrix& stvnl_dphi,
+        ModuleBase::matrix& svnl_dbeta,
+        ModuleBase::matrix& svl_dphi,
+#ifdef __DEEPKS
+        ModuleBase::matrix& svnl_dalpha,
+#endif
+        LCAO_gen_fixedH& gen_h, // mohan add 2024-04-02
+        typename TGint<T>::type& gint,
+        const Parallel_Orbitals& pv,
+        LCAO_Matrix& lm,
+        const K_Vectors* kv = nullptr,
+        Record_adj* ra = nullptr);
+
+
+    // get the ds, dt, dvnl.
+    void allocate(const Parallel_Orbitals& pv,
+        LCAO_Matrix& lm,
+        LCAO_gen_fixedH& gen_h,
+        const int& nks = 0,
+        const std::vector<ModuleBase::Vector3<double>>& kvec_d = {});
+
+
+    void finish_ftable(LCAO_Matrix& lm);
+
+    void average_force(double* fm);
+
+    void test(Parallel_Orbitals& pv, double* mm, const std::string& name);
+
+    //-------------------------------------------------------------
+    // forces reated to overlap matrix
+    // forces related to energy density matrix
+    //-------------------------------------------------------------
+
+    void cal_foverlap(const bool isforce,
+        const bool isstress,
+        const psi::Psi<T>* psi,
+        const Parallel_Orbitals& pv,
+        const elecstate::ElecState* pelec,
+        LCAO_Matrix& lm,
+        ModuleBase::matrix& foverlap,
+        ModuleBase::matrix& soverlap,
+        const K_Vectors* kv = nullptr,
+        Record_adj* ra = nullptr,
+        const elecstate::DensityMatrix<T, double>* DM = nullptr);
+
+    //-------------------------------------------------------------
+    // forces related to kinetic and non-local pseudopotentials
+    //--------------------------------------------------------------
+    void cal_ftvnl_dphi(const elecstate::DensityMatrix<T, double>* DM,
+        const Parallel_Orbitals& pv,
+        const UnitCell& ucell,
+        LCAO_Matrix& lm,
+        const bool isforce,
+        const bool isstress,
+        ModuleBase::matrix& ftvnl_dphi,
+        ModuleBase::matrix& stvnl_dphi,
+        Record_adj* ra = nullptr);
+
+    void cal_fvnl_dbeta(const elecstate::DensityMatrix<T, double>* DM,
+        const Parallel_Orbitals& pv,
+        const UnitCell& ucell,
+        const LCAO_Orbitals& orb,
+        const ORB_gen_tables& uot,
+        Grid_Driver& gd,
+        const bool isforce,
+        const bool isstress,
+        ModuleBase::matrix& fvnl_dbeta,
+        ModuleBase::matrix& svnl_dbeta);
+
+    //-------------------------------------------
+    // forces related to local pseudopotentials
+    //-------------------------------------------
+    void cal_fvl_dphi(const bool isforce,
+        const bool isstress,
+        const elecstate::Potential* pot_in,
+        typename TGint<T>::type& gint,
+        ModuleBase::matrix& fvl_dphi,
+        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_STRESS.cpp

@@ -255,7 +255,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
 
     if (!GlobalV::GAMMA_ONLY_LOCAL)
     {
-        this->flk.finish_k(lm);
+        this->flk.finish_ftable(lm);
     }
 
 #ifdef __EXX
@@ -751,15 +751,14 @@ void Force_Stress_LCAO<double>::integral_part(
     LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 	Gint_Gamma &gint_gamma, // mohan add 2024-04-01
 	Gint_k &gint_k, // mohan add 2024-04-01
-	Parallel_Orbitals &pv,
+    const Parallel_Orbitals& pv,
     LCAO_Matrix &lm,
     const K_Vectors& kv)
 {
 
-    flk.ftable_gamma(isforce,
+    flk.ftable(isforce,
         isstress,
         psi,
-        loc,
         pelec,
         foverlap,
         ftvnl_dphi,
@@ -768,14 +767,13 @@ void Force_Stress_LCAO<double>::integral_part(
         soverlap,
         stvnl_dphi,
         svnl_dbeta,
-#if __DEEPKS
         svl_dphi,
+#if __DEEPKS
         svnl_dalpha,
-#else
-        svl_dphi,
 #endif
         gen_h,
         gint_gamma,
+        pv,
         lm);
     return;
 }
@@ -805,34 +803,31 @@ void Force_Stress_LCAO<std::complex<double>>::integral_part(
     LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 	Gint_Gamma &gint_gamma,
 	Gint_k &gint_k,
-	Parallel_Orbitals &pv,
+    const Parallel_Orbitals& pv,
 	LCAO_Matrix &lm,
 	const K_Vectors& kv)
 {
-        flk.ftable_k(isforce,
-                     isstress,
-                     *this->RA,
-                     psi,
-                     loc,
-                     pelec,
-                     foverlap,
-                     ftvnl_dphi,
-                     fvnl_dbeta,
-                     fvl_dphi,
-                     soverlap,
-                     stvnl_dphi,
-                     svnl_dbeta,
+    flk.ftable(isforce,
+        isstress,
+        psi,
+        pelec,
+        foverlap,
+        ftvnl_dphi,
+        fvnl_dbeta,
+        fvl_dphi,
+        soverlap,
+        stvnl_dphi,
+        svnl_dbeta,
+        svl_dphi,
 #if __DEEPKS
-                     svl_dphi,
-                     svnl_dalpha,
-#else
-                     svl_dphi,
+        svnl_dalpha,
 #endif
-					 gen_h,
-                     gint_k,
-					 pv,
-					 lm,
-                     kv);
+        gen_h,
+        gint_k,
+        pv,
+        lm,
+        & kv,
+        this->RA);
     return;
 }
 

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h

@@ -1,7 +1,7 @@
 #ifndef FORCE_STRESS_LCAO_H
 #define FORCE_STRESS_LCAO_H
 
-#include "FORCE_k.h"
+#include "FORCE.h"
 #include "module_base/global_function.h"
 #include "module_base/global_variable.h"
 #include "module_base/matrix.h"
@@ -54,8 +54,7 @@ class Force_Stress_LCAO
   private:
     int nat;
     Record_adj* RA;
-    Force_LCAO_k flk;
-    //	Force_LCAO_gamma flg;
+    Force_LCAO<T> flk;
     Stress_Func<double> sc_pw;
     Forces<double> f_pw;
 
@@ -95,7 +94,7 @@ class Force_Stress_LCAO
         LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 		Gint_Gamma &gint_gamma,
 		Gint_k &gint_k,
-	    Parallel_Orbitals &pv,
+	    const Parallel_Orbitals &pv,
 		LCAO_Matrix &lm,
 		const K_Vectors& kv);