Feature: update LibRI v0.2.0 (deepmodeling#4830)

* Fix C++ compiler warning

* Refactor Exx_LRI::cal_exx_stress() to LibRI loop3

* Feature: add input and info CV_grad_R_threshold for exx

* Delete exx_lri.set_csm_threshold()

* Update LibRI in toolchain to v0.2.0

* Update submodule deps/LibRI to v0.2.0

* fix include

* change ref in integrate_test

* change ref in integrate_test

* change ref results for no cauchy screening (dyzheng#3)

---------

Co-authored-by: LUNASEA <33978601+maki49@users.noreply.github.com>

Author: PeizeLin

deps/LibRI

@@ -47,6 +47,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;

source/module_hamilt_general/module_xc/exx_info.h

@@ -497,6 +497,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);

source/module_io/input_conv.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";

source/module_io/read_input_item_exx_dftu.cpp

@@ -471,6 +471,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

source/module_parameter/input_parameter.h

@@ -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");
 }
@@ -172,8 +182,6 @@ void Exx_LRI<Tdata>::cal_exx_elec(const std::vector<std::map<TA,std::map<TAC,RI:
 
 	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)
@@ -209,7 +217,7 @@ void Exx_LRI<Tdata>::post_process_Hexx( std::map<TA, std::map<TAC, RI::Tensor<Td
 }
 
 template<typename Tdata>
-double Exx_LRI<Tdata>::post_process_Eexx(const double& Eexx_in) const
+double Exx_LRI<Tdata>::post_process_Eexx(const double& Eexx_in) const
 {
 	ModuleBase::TITLE("Exx_LRI","post_process_Eexx");
     const double SPIN_multiple = std::map<int, double>{ {1,2}, {2,1}, {4,1} }.at(GlobalV::NSPIN);				// why?
@@ -239,8 +247,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)
@@ -263,8 +269,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)

source/module_ri/Exx_LRI.hpp

@@ -8,10 +8,11 @@
 
 #include "LRI_CV.h"
 #include "LRI_CV_Tools.h"
-#include "module_ri/exx_abfs-abfs_index.h"
+#include "exx_abfs-abfs_index.h"
 #include "RI_Util.h"
-#include "module_base/tool_title.h"
-#include "module_base/timer.h"
+#include "../module_base/tool_title.h"
+#include "../module_base/timer.h"
+#include "../module_hamilt_pw/hamilt_pwdft/global.h"
 #include <RI/global/Global_Func-1.h>
 #include <omp.h>
 

source/module_ri/LRI_CV.hpp

@@ -6,11 +6,14 @@
 #ifndef LRI_CV_TOOLS_H
 #define LRI_CV_TOOLS_H
 
+#include "module_base/abfs-vector3_order.h"
+
 #include <RI/global/Tensor.h>
 
 #include <cstddef>
 #include <array>
 #include <vector>
+#include <map>
 
 namespace LRI_CV_Tools
 {
@@ -84,7 +87,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"

source/module_ri/LRI_CV_Tools.h

@@ -9,6 +9,7 @@
 #include "LRI_CV_Tools.h"
 #include "Inverse_Matrix.h"
 #include "../module_base/mathzone.h"
+#include "../module_hamilt_pw/hamilt_pwdft/global.h"
 
 template<typename Tdata>
 RI::Tensor<Tdata>
@@ -54,8 +55,10 @@ template<typename T, std::size_t N>
 bool LRI_CV_Tools::exist(const std::array<T,N> &dV)
 {
 	for(size_t i=0; i<3; ++i)
+	{
 		if(!dV[i].empty())
 			return true;
+	}
 	return false;
 }
 
@@ -312,4 +315,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

source/module_ri/LRI_CV_Tools.hpp

@@ -1,6 +1,6 @@
-6 (-1.039229,-3.6984385e-16) (-1.4337629e-08,8.1948535e-16) (-1.6549262e-15,-2.4275964e-17) (-7.3477682e-14,-2.8317688e-17) (-2.2743308e-08,6.0642134e-16) (-4.7365417e-10,2.7036238e-14)
- (-0.96491479,1.2723875e-15) (-6.2469742e-09,-1.8654043e-12) (-8.1614632e-11,8.9132756e-13) (0.12279288,7.8499298e-15) (9.3436543e-14,2.364213e-16)
- (-1.064337,2.3693393e-16) (1.0547119e-15,1.5451914e-16) (-6.6444545e-09,2.0880435e-12) (-5.7610897e-10,-3.0200736e-12)
- (-1.064337,2.3700845e-16) (6.3717573e-10,-9.9770901e-13) (4.9778198e-10,-3.4757565e-13)
- (-0.65449032,1.2140463e-16) (1.8561136e-14,4.810951e-17)
- (-0.65138796,-7.4288639e-18)
+6 (-1.0394425,-3.7363449e-16) (3.9703113e-10,8.3455762e-16) (4.4417161e-14,5.5568749e-17) (-2.4341466e-13,4.3883506e-17) (-2.0039369e-09,5.8451215e-16) (1.9596449e-11,6.8294051e-12)
+ (-0.96326621,1.2640312e-15) (-1.4831543e-09,-4.0086346e-12) (8.370459e-10,-5.1125035e-13) (0.12202161,-5.3763309e-15) (1.7561473e-14,3.428681e-15)
+ (-1.0644129,2.3743998e-16) (-7.1678774e-15,-2.6869918e-16) (-1.4937099e-09,4.0465037e-12) (2.5345497e-09,8.0997135e-10)
+ (-1.0644129,2.3735151e-16) (8.1785897e-10,5.1608006e-13) (-7.4156745e-10,-2.3462101e-10)
+ (-0.65468695,1.1556487e-16) (4.4158931e-15,1.0408529e-15)
+ (-0.65149008,-4.4004079e-17)