diff --git a/deps/LibRI b/deps/LibRI index 553c91c0be..3883b699f5 160000 --- a/deps/LibRI +++ b/deps/LibRI @@ -1 +1 @@ -Subproject commit 553c91c0be1d60a86e7666f0502ef866c366c600 +Subproject commit 3883b699f50eccd84993e849843452e22d035015 diff --git a/source/module_hamilt_general/module_xc/exx_info.h b/source/module_hamilt_general/module_xc/exx_info.h index 161e4518e3..d94f5e6375 100644 --- 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 index db0c3adab1..4f3ee9c6f4 100644 --- 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 index 79edc45632..057e14298f 100644 --- 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 index 6b9b454daa..4add50fa0a 100644 --- 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 index f1b4a7934b..d60c91117a 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -140,6 +140,11 @@ void Exx_LRI::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>>,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 period_Cs = LRI_CV_Tools::cal_latvec_range(2); @@ -160,6 +165,11 @@ void Exx_LRI::cal_exx_ions() std::array>>,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>>,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::cal_exx_elec(const std::vector, std::set>> 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::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::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 extern std::map,std::array,3>>>> get_dCVws( - const std::array>,RI::Tensor>>,3> &dCVs); + const std::array>,RI::Tensor>>,3> &dCVs); + + template + extern std::array,RI::Tensor>>,3>,3> + cal_dMRs( + const std::array,RI::Tensor>>,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 index 8ad95c3715..6333c208b8 100644 --- 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 +std::array,RI::Tensor>>,3>,3> +LRI_CV_Tools::cal_dMRs( + const std::array,RI::Tensor>>,3> &dMs) +{ + auto get_R_delta = [&](const TA &iat0, const std::pair &A1) -> std::array + { + 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 tau0 = GlobalC::ucell.atoms[it0].tau[ia0]; + const ModuleBase::Vector3 tau1 = GlobalC::ucell.atoms[it1].tau[ia1]; + const Abfs::Vector3_Order R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*GlobalC::ucell.latvec); + return std::array{R_delta.x, R_delta.y, R_delta.z}; + }; + constexpr int Npos = 3; + std::array,RI::Tensor>>,Npos>,Npos> dMRs; + for(int ipos0=0; ipos0 A1 = dMs_B.first; + const RI::Tensor &dM = dMs_B.second; + const std::array R_delta = get_R_delta(iat0, A1); + dMRs[ipos0][ipos1][iat0][A1] = dM * R_delta[ipos1]; + } + } + } + } + return dMRs; +} #endif \ No newline at end of file diff --git a/toolchain/scripts/stage4/install_libri.sh b/toolchain/scripts/stage4/install_libri.sh index d33fe3068f..78af9ab2ab 100755 --- 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