update pulay_fs files

Author: mohanchen

source/Makefile.Objects

@@ -655,7 +655,7 @@ OBJS_LCAO=evolve_elec.o\
 	  FORCE_k.o\
       stress_tools.o\
       edm.o\
-      pulay_force_stress_center2.o\
+      pulay_fs_center2.o\
 	  grid_init.o\
       spar_dh.o\
       spar_exx.o\

source/source_lcao/CMakeLists.txt

@@ -22,7 +22,7 @@ if(ENABLE_LCAO)
         module_operator_lcao/td_pot_hybrid.cpp
         module_operator_lcao/dspin_lcao.cpp
         module_operator_lcao/dftu_lcao.cpp
-        pulay_force_stress_center2.cpp
+        pulay_fs_center2.cpp
         FORCE_STRESS.cpp
         FORCE_gamma.cpp
         FORCE_k.cpp

source/source_lcao/FORCE_gamma.cpp

@@ -12,7 +12,7 @@
 #include "source_cell/module_neighbor/sltk_grid_driver.h" //GridD
 #include "source_estate/elecstate_lcao.h"
 #include "source_lcao/LCAO_domain.h"
-#include "source_lcao/pulay_force_stress.h"
+#include "source_lcao/pulay_fs.h"
 #include "source_io/write_HS.h"
 
 template <>

source/source_lcao/FORCE_k.cpp

@@ -9,7 +9,7 @@
 #include "source_estate/elecstate_lcao.h"
 #include "source_estate/module_dm/cal_dm_psi.h"
 #include "source_lcao/LCAO_domain.h"
-#include "source_lcao/pulay_force_stress.h"
+#include "source_lcao/pulay_fs.h"
 #include "source_pw/module_pwdft/global.h"
 #include "source_io/write_HS.h"
 #include "source_io/module_parameter/parameter.h"

source/source_lcao/pulay_force_stress_center2.cpp

@@ -65,5 +65,5 @@ namespace PulayForceStress
         const bool& isstress,
         const bool& set_dmr_gint = true);
 }
-#include "pulay_force_stress_center2_template.hpp"
-#include "pulay_force_stress_gint.hpp"
+#include "pulay_fs_temp.hpp"
+#include "pulay_fs_gint.hpp"

source/source_lcao/pulay_force_stress.h renamed to source/source_lcao/pulay_fs.h

@@ -0,0 +1,144 @@
+#include "pulay_fs.h"
+
+template<>  // gamma-only, provided xy
+void PulayForceStress::cal_pulay_fs(
+		ModuleBase::matrix& force,
+		ModuleBase::matrix& stress,
+		const elecstate::DensityMatrix<double, double>& dm,
+		const UnitCell& ucell,
+		const Parallel_Orbitals& pv,
+		const double* (&dHSx)[3],
+		const double* (&dHSxy)[6],
+		const bool& isforce,
+		const bool& isstress,
+		Record_adj* ra,
+		const double& factor_force,
+		const double& factor_stress)
+{
+	ModuleBase::TITLE("Forces", "cal_pulay_fs");
+	ModuleBase::timer::tick("Forces", "cal_pulay_fs");
+
+	const int nspin = PARAM.inp.nspin;
+	const int nlocal = PARAM.globalv.nlocal;
+
+	for (int i = 0; i < nlocal; ++i)
+	{
+		const int iat = ucell.iwt2iat[i];
+		for (int j = 0; j < nlocal; ++j)
+		{
+			const int mu = pv.global2local_row(j);
+			const int nu = pv.global2local_col(i);
+
+			if (mu >= 0 && nu >= 0)
+			{
+				const int index = mu * pv.ncol + nu;
+				double sum = 0.0;
+				for (int is = 0; is < nspin; ++is) 
+				{ 
+					sum += dm.get_DMK(is + 1, 0, nu, mu); 
+				}
+				if (isforce)
+				{
+					const double sumf = sum * factor_force;
+					for (int i = 0; i < 3; ++i) 
+					{ 
+						force(iat, i) += sumf * 2.0 * dHSx[i][index]; 
+					}
+				}
+				if (isstress)
+				{
+					const double sums = sum * factor_stress;
+					int ij = 0;
+					for (int i = 0; i < 3;++i) 
+					{ 
+						for (int j = i; j < 3; ++j) 
+						{ 
+							stress(i, j) += sums * dHSxy[ij++][index]; 
+						} 
+					}
+				}
+			}
+		}
+	}
+
+	if (isstress) 
+	{ 
+		StressTools::stress_fill(ucell.lat0, ucell.omega, stress);
+	}
+
+	ModuleBase::timer::tick("Forces", "cal_pulay_fs");
+}
+
+
+template<>  //multi-k, provided xy
+void PulayForceStress::cal_pulay_fs(
+		ModuleBase::matrix& force,
+		ModuleBase::matrix& stress,
+		const elecstate::DensityMatrix<std::complex<double>, double>& dm,
+		const UnitCell& ucell,
+		const Parallel_Orbitals& pv,
+		const double* (&dHSx)[3],
+		const double* (&dHSxy)[6],
+		const bool& isforce,
+		const bool& isstress,
+		Record_adj* ra,
+		const double& factor_force,
+		const double& factor_stress)
+{
+	auto stress_func = [](ModuleBase::matrix& local_s, 
+			const double& dm2d1_s, 
+			const double** dHSx, 
+			const double** dHSxy, 
+			const double* dtau, 
+			const int& irr)
+	{
+		int ij = 0;
+		for (int i = 0; i < 3; ++i) 
+		{ 
+			for (int j = i; j < 3; ++j) 
+			{ 
+				local_s(i, j) += dm2d1_s * dHSxy[ij++][irr]; 
+			} 
+		}
+	};
+	cal_pulay_fs(force, stress, dm, ucell, pv, dHSx, dHSxy, 
+			nullptr, isforce, isstress, ra, 
+			factor_force, factor_stress, stress_func);
+}
+
+
+template<>  // multi-k, provided x
+void PulayForceStress::cal_pulay_fs(
+        ModuleBase::matrix& force,
+        ModuleBase::matrix& stress,
+        const elecstate::DensityMatrix<std::complex<double>, double>& dm,
+        const UnitCell& ucell,
+        const Parallel_Orbitals& pv,
+        const double* (&dHSx)[3],
+        const double* dtau,
+        const bool& isforce,
+        const bool& isstress,
+        Record_adj* ra,
+        const double& factor_force,
+        const double& factor_stress)
+{
+	auto stress_func = [](ModuleBase::matrix& local_s, 
+			const double& dm2d1_s, 
+			const double** dHSx, 
+			const double** dHSxy, 
+			const double* dtau, 
+			const int& irr)
+	{
+		for (int i = 0; i < 3; ++i) 
+		{ 
+			for (int j = i; j < 3; ++j) 
+			{ 
+				local_s(i, j) += dm2d1_s * dHSx[i][irr] * dtau[irr * 3 + j]; 
+			} 
+		}
+	};
+	cal_pulay_fs(force, stress, dm, ucell, pv, dHSx, 
+			nullptr, dtau, isforce, isstress, ra, 
+			factor_force, factor_stress, stress_func);
+}
+

source/source_lcao/pulay_fs_center2.cpp

@@ -1,5 +1,5 @@
 #pragma once
-#include "pulay_force_stress.h"
+#include "pulay_fs.h"
 #include "source_lcao/stress_tools.h"
 #include "source_hamilt/module_xc/xc_functional.h"
 #include "source_io/module_parameter/parameter.h"
@@ -62,4 +62,4 @@ namespace PulayForceStress
 
         if (isstress) { StressTools::stress_fill(-1.0, ucell.omega, s); }
     }
-}
+}

source/source_lcao/pulay_force_stress_gint.hpp renamed to source/source_lcao/pulay_fs_gint.hpp

@@ -1,5 +1,5 @@
 #pragma once
-#include "pulay_force_stress.h"
+#include "pulay_fs.h"
 #include "source_base/timer.h"
 #include "source_io/module_parameter/parameter.h"
 namespace PulayForceStress