Merge branch 'LibRI_loop3' of https://github.com/PeizeLin/abacus-develop

into symexxnew
maki49 · Jul 30, 2024 · 0bf4896 · 0bf4896
2 parents fe4a6b5 + ab01840
commit 0bf4896
Show file tree

Hide file tree

Showing 9 changed files with 78 additions and 10 deletions.
diff --git a/deps/LibRI b/deps/LibRI
diff --git a/source/module_hamilt_general/module_xc/exx_info.h b/source/module_hamilt_general/module_xc/exx_info.h
@@ -48,6 +48,8 @@ struct Exx_Info
         double cauchy_threshold = 0;
         double C_grad_threshold = 0;
         double V_grad_threshold = 0;
+        double C_grad_R_threshold = 0;
+        double V_grad_R_threshold = 0;
         double cauchy_force_threshold = 0;
         double cauchy_stress_threshold = 0;
         double ccp_rmesh_times = 10;

diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp
@@ -507,6 +507,8 @@ void Input_Conv::Convert()
         GlobalC::exx_info.info_ri.cauchy_threshold = PARAM.inp.exx_cauchy_threshold;
         GlobalC::exx_info.info_ri.C_grad_threshold = PARAM.inp.exx_c_grad_threshold;
         GlobalC::exx_info.info_ri.V_grad_threshold = PARAM.inp.exx_v_grad_threshold;
+        GlobalC::exx_info.info_ri.C_grad_R_threshold = PARAM.inp.exx_c_grad_r_threshold;
+        GlobalC::exx_info.info_ri.V_grad_R_threshold = PARAM.inp.exx_v_grad_r_threshold;
         GlobalC::exx_info.info_ri.cauchy_force_threshold = PARAM.inp.exx_cauchy_force_threshold;
         GlobalC::exx_info.info_ri.cauchy_stress_threshold = PARAM.inp.exx_cauchy_stress_threshold;
         GlobalC::exx_info.info_ri.ccp_rmesh_times = std::stod(PARAM.inp.exx_ccp_rmesh_times);

diff --git a/source/module_io/read_input_item_exx_dftu.cpp b/source/module_io/read_input_item_exx_dftu.cpp
@@ -149,6 +149,18 @@ void ReadInput::item_exx()
         read_sync_double(input.exx_v_grad_threshold);
         this->add_item(item);
     }
+    {
+        Input_Item item("exx_c_grad_r_threshold");
+        item.annotation = "threshold to screen nabla C matrix in exx";
+        read_sync_double(input.exx_c_grad_r_threshold);
+        this->add_item(item);
+    }
+    {
+        Input_Item item("exx_v_grad_r_threshold");
+        item.annotation = "threshold to screen nabla V matrix in exx";
+        read_sync_double(input.exx_v_grad_r_threshold);
+        this->add_item(item);
+    }
     {
         Input_Item item("exx_cauchy_force_threshold");
         item.annotation = "threshold to screen exx force using Cauchy-Schwartz inequality";

diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h
@@ -470,6 +470,8 @@ struct Input_para
     double exx_cauchy_threshold = 1e-07;        ///< threshold to screen exx using Cauchy-Schwartz inequality
     double exx_c_grad_threshold = 0.0001;       ///< threshold to screen nabla C matrix in exx
     double exx_v_grad_threshold = 0.1;          ///< threshold to screen nabla V matrix in exx
+    double exx_c_grad_r_threshold = 0.0001;     ///< threshold to screen nabla C matrix in exx
+    double exx_v_grad_r_threshold = 0.1;        ///< threshold to screen nabla V matrix in exx
     double exx_cauchy_force_threshold = 1e-07;  ///< threshold to screen exx force using Cauchy-Schwartz
                                                 ///< inequality
     double exx_cauchy_stress_threshold = 1e-07; ///< threshold to screen exx stress using Cauchy-Schwartz

diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp
@@ -140,6 +140,11 @@ void Exx_LRI<Tdata>::cal_exx_ions()
 				{{"writable_dVws",true}});
 		this->cv.dVws = LRI_CV_Tools::get_dCVws(dVs);
 		this->exx_lri.set_dVs(std::move(dVs), this->info.V_grad_threshold);
+		if(GlobalV::CAL_STRESS)
+		{
+			std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dVRs = LRI_CV_Tools::cal_dMRs(dVs);
+			this->exx_lri.set_dVRs(std::move(dVRs), this->info.V_grad_R_threshold);
+		}
 	}
 
 	const std::array<Tcell,Ndim> period_Cs = LRI_CV_Tools::cal_latvec_range<Tcell>(2);
@@ -160,6 +165,11 @@ void Exx_LRI<Tdata>::cal_exx_ions()
 		std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3> &dCs = std::get<1>(Cs_dCs);
 		this->cv.dCws = LRI_CV_Tools::get_dCVws(dCs);
 		this->exx_lri.set_dCs(std::move(dCs), this->info.C_grad_threshold);
+		if(GlobalV::CAL_STRESS)
+		{
+			std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dCRs = LRI_CV_Tools::cal_dMRs(dCs);
+			this->exx_lri.set_dCRs(std::move(dCRs), this->info.C_grad_R_threshold);
+		}
 	}
 	ModuleBase::timer::tick("Exx_LRI", "cal_exx_ions");
 }
@@ -174,8 +184,6 @@ void Exx_LRI<Tdata>::cal_exx_elec(const std::vector<std::map<TA, std::map<TAC, R
 
 	const std::vector<std::tuple<std::set<TA>, std::set<TA>>> judge = RI_2D_Comm::get_2D_judge(pv);
 
-	this->exx_lri.set_csm_threshold(this->info.cauchy_threshold);
-
 	this->Hexxs.resize(GlobalV::NSPIN);
 	this->Eexx = 0;
 	for(int is=0; is<GlobalV::NSPIN; ++is)
@@ -247,8 +255,6 @@ void Exx_LRI<Tdata>::cal_exx_force()
 {
 	ModuleBase::TITLE("Exx_LRI","cal_exx_force");
 	ModuleBase::timer::tick("Exx_LRI", "cal_exx_force");
-
-	this->exx_lri.set_csm_threshold(this->info.cauchy_force_threshold);
 
 	this->force_exx.create(GlobalC::ucell.nat, Ndim);
 	for(int is=0; is<GlobalV::NSPIN; ++is)
@@ -271,8 +277,6 @@ void Exx_LRI<Tdata>::cal_exx_stress()
 {
 	ModuleBase::TITLE("Exx_LRI","cal_exx_stress");
 	ModuleBase::timer::tick("Exx_LRI", "cal_exx_stress");
-
-	this->exx_lri.set_csm_threshold(this->info.cauchy_stress_threshold);
 
 	this->stress_exx.create(Ndim, Ndim);
 	for(int is=0; is<GlobalV::NSPIN; ++is)

diff --git a/source/module_ri/LRI_CV_Tools.h b/source/module_ri/LRI_CV_Tools.h
@@ -84,7 +84,12 @@ namespace LRI_CV_Tools
 	template<typename TA, typename Tcell, typename Tdata>
 	extern std::map<int,std::map<int,std::map<Abfs::Vector3_Order<double>,std::array<RI::Tensor<Tdata>,3>>>>
 	get_dCVws(
-		const std::array<std::map<TA,std::map<std::pair<TA,std::array<Tcell,3>>,RI::Tensor<Tdata>>>,3> &dCVs);	
+		const std::array<std::map<TA,std::map<std::pair<TA,std::array<Tcell,3>>,RI::Tensor<Tdata>>>,3> &dCVs);
+
+	template<typename TA, typename TC, typename Tdata>
+	extern std::array<std::array<std::map<TA,std::map<std::pair<TA,TC>,RI::Tensor<Tdata>>>,3>,3>
+	cal_dMRs(
+		const std::array<std::map<TA,std::map<std::pair<TA,TC>,RI::Tensor<Tdata>>>,3> &dMs);
 }
 
 #include "LRI_CV_Tools.hpp"

diff --git a/source/module_ri/LRI_CV_Tools.hpp b/source/module_ri/LRI_CV_Tools.hpp
@@ -312,4 +312,45 @@ LRI_CV_Tools::get_dCVws(
 	return dCVws;
 }
 
+
+// dMRs[ipos0][ipos1] = \nabla_{ipos0} M R_{ipos1}
+template<typename TA, typename TC, typename Tdata>
+std::array<std::array<std::map<TA,std::map<std::pair<TA,TC>,RI::Tensor<Tdata>>>,3>,3>
+LRI_CV_Tools::cal_dMRs(
+	const std::array<std::map<TA,std::map<std::pair<TA,TC>,RI::Tensor<Tdata>>>,3> &dMs)
+{
+	auto get_R_delta = [&](const TA &iat0, const std::pair<TA,TC> &A1) -> std::array<Tdata,3>
+	{
+		const TA iat1 = A1.first;
+		const TC &cell1 = A1.second;
+		const int it0 = GlobalC::ucell.iat2it[iat0];
+		const int ia0 = GlobalC::ucell.iat2ia[iat0];
+		const int it1 = GlobalC::ucell.iat2it[iat1];
+		const int ia1 = GlobalC::ucell.iat2ia[iat1];
+		const ModuleBase::Vector3<double> tau0 = GlobalC::ucell.atoms[it0].tau[ia0];
+		const ModuleBase::Vector3<double> tau1 = GlobalC::ucell.atoms[it1].tau[ia1];
+		const Abfs::Vector3_Order<double> R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*GlobalC::ucell.latvec);
+		return std::array<Tdata,3>{R_delta.x, R_delta.y, R_delta.z};
+	};
+	constexpr int Npos = 3;
+	std::array<std::array<std::map<TA,std::map<std::pair<TA,TC>,RI::Tensor<Tdata>>>,Npos>,Npos> dMRs;
+	for(int ipos0=0; ipos0<Npos; ++ipos0)
+	{
+		for(int ipos1=0; ipos1<Npos; ++ipos1)
+		{
+			for(const auto &dMs_A : dMs[ipos0])
+			{
+				const TA iat0 = dMs_A.first;
+				for(const auto &dMs_B : dMs_A.second)
+				{
+					const std::pair<TA,TC> A1 = dMs_B.first;
+					const RI::Tensor<Tdata> &dM = dMs_B.second;
+					const std::array<Tdata,3> R_delta = get_R_delta(iat0, A1);
+					dMRs[ipos0][ipos1][iat0][A1] = dM * R_delta[ipos1];
+				}
+			}
+		}
+	}
+	return dMRs;
+}
 #endif
diff --git a/toolchain/scripts/stage4/install_libri.sh b/toolchain/scripts/stage4/install_libri.sh
@@ -11,8 +11,8 @@
 [ "${BASH_SOURCE[0]}" ] && SCRIPT_NAME="${BASH_SOURCE[0]}" || SCRIPT_NAME=$0
 SCRIPT_DIR="$(cd "$(dirname "$SCRIPT_NAME")/.." && pwd -P)"
 
-libri_ver="0.1.1"
-libri_sha256="51deb08aa373e54d2c123b57bfd4b3507accac0d496a94b766eaeadccd9e4bd0"
+libri_ver="0.2.0"
+libri_sha256="ad79dfbc3ed8ff066c85549a2737d29205dbf755b38ea216ab2ab42754f80389"
 source "${SCRIPT_DIR}"/common_vars.sh
 source "${SCRIPT_DIR}"/tool_kit.sh
 source "${SCRIPT_DIR}"/signal_trap.sh