From 01f017922eea18f9a3196a43f2039be464bff453 Mon Sep 17 00:00:00 2001 From: linpz Date: Sun, 14 Sep 2025 18:18:44 +0800 Subject: [PATCH 1/2] Refactor: change cal_proj() to cal_mul() in Exx_Opt_Orb --- source/source_lcao/module_ri/exx_opt_orb.cpp | 96 ++++++++++---------- source/source_lcao/module_ri/exx_opt_orb.h | 12 +-- 2 files changed, 51 insertions(+), 57 deletions(-) diff --git a/source/source_lcao/module_ri/exx_opt_orb.cpp b/source/source_lcao/module_ri/exx_opt_orb.cpp index c22c2f8939..65ac08b9ed 100644 --- a/source/source_lcao/module_ri/exx_opt_orb.cpp +++ b/source/source_lcao/module_ri/exx_opt_orb.cpp @@ -174,22 +174,19 @@ void Exx_Opt_Orb::generate_matrix( const std::vector>> ms_abfs_abfs_I = cal_I( ms_abfs_abfs, T,I,T,I ); // < lcaos lcaos | lcaos lcaos > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | lcaos lcaos > const RI::Tensor m_lcaoslcaos_lcaoslcaos_proj = - cal_proj_22( - ms_lcaoslcaos_lcaoslcaos.at(T).at(I).at(T).at(I), + ms_lcaoslcaos_lcaoslcaos.at(T).at(I).at(T).at(I) - cal_mul_22( ms_lcaoslcaos_abfs.at(T).at(I).at(T).at(I), ms_abfs_abfs_I, ms_lcaoslcaos_abfs.at(T).at(I).at(T).at(I)); // < lcaos lcaos | jys > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector> m_lcaoslcaos_jys_proj = - {cal_proj_21( - ms_lcaoslcaos_jys.at(T).at(I).at(T).at(I)[0], + {ms_lcaoslcaos_jys.at(T).at(I).at(T).at(I)[0] - cal_mul_21( ms_lcaoslcaos_abfs.at(T).at(I).at(T).at(I), ms_abfs_abfs_I, {ms_jys_abfs.at(T).at(I).at(T).at(I)})}; // < jys | jys > - < jys | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector>> m_jys_jys_proj = - {{cal_proj_11( - ms_jys_jys.at(T).at(I).at(T).at(I), + {{ms_jys_jys.at(T).at(I).at(T).at(I) - cal_mul_11( {ms_jys_abfs.at(T).at(I).at(T).at(I)}, ms_abfs_abfs_I, {ms_jys_abfs.at(T).at(I).at(T).at(I)})}}; @@ -228,42 +225,35 @@ void Exx_Opt_Orb::generate_matrix( const std::vector>> ms_abfs_abfs_I = cal_I( ms_abfs_abfs, TA,IA,TB,IB ); // < lcaos lcaos | lcaos lcaos > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | lcaos lcaos > const RI::Tensor m_lcaoslcaos_lcaoslcaos_proj = - cal_proj_22( - ms_lcaoslcaos_lcaoslcaos.at(TA).at(IA).at(TB).at(IB), + ms_lcaoslcaos_lcaoslcaos.at(TA).at(IA).at(TB).at(IB) - cal_mul_22( ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB), ms_abfs_abfs_I, ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB)); // < lcaos lcaos | jys > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector> m_lcaoslcaos_jys_proj = - {cal_proj_21( - ms_lcaoslcaos_jys.at(TA).at(IA).at(TB).at(IB)[0], + {ms_lcaoslcaos_jys.at(TA).at(IA).at(TB).at(IB)[0] - cal_mul_21( ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB), ms_abfs_abfs_I, { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }), - cal_proj_21( - ms_lcaoslcaos_jys.at(TA).at(IA).at(TB).at(IB)[1], + ms_lcaoslcaos_jys.at(TA).at(IA).at(TB).at(IB)[1] - cal_mul_21( ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB), ms_abfs_abfs_I, { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) })}; // < jys | jys > - < jys | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector>> m_jys_jys_proj = - {{cal_proj_11( - ms_jys_jys.at(TA).at(IA).at(TA).at(IA), + {{ms_jys_jys.at(TA).at(IA).at(TA).at(IA) - cal_mul_11( { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }), - cal_proj_11( - ms_jys_jys.at(TA).at(IA).at(TB).at(IB), + ms_jys_jys.at(TA).at(IA).at(TB).at(IB) - cal_mul_11( { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }) }, - {cal_proj_11( - ms_jys_jys.at(TB).at(IB).at(TA).at(IA), + {ms_jys_jys.at(TB).at(IB).at(TA).at(IA) - cal_mul_11( { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }), - cal_proj_11( - ms_jys_jys.at(TB).at(IB).at(TB).at(IB), + ms_jys_jys.at(TB).at(IB).at(TB).at(IB) - cal_mul_11( { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }) }}; @@ -301,86 +291,94 @@ void Exx_Opt_Orb::generate_matrix( } } -// m_big - m_left * m_middle * m_right.T -RI::Tensor Exx_Opt_Orb::cal_proj_22( - const RI::Tensor & m_big, +// m_left * m_middle * m_right.T +RI::Tensor Exx_Opt_Orb::cal_mul_22( const std::vector> & m_left, const std::vector>> & m_middle, const std::vector> & m_right ) const { - ModuleBase::TITLE("Exx_Opt_Orb::cal_proj_22"); - RI::Tensor m_proj = m_big.copy(); + ModuleBase::TITLE("Exx_Opt_Orb::cal_mul_22"); + RI::Tensor m_mul; for( size_t il=0; il!=m_left.size(); ++il ) { for( size_t ir=0; ir!=m_right.size(); ++ir ) { - // m_proj = m_proj - m_left[il] * m_middle[il][ir] * m_right[ir].T; + // m_mul += m_left[il] * m_middle[il][ir] * m_right[ir].T; const RI::Tensor m_lm = RI::Tensor_Multiply::x0x1y1_x0x1a_ay1(m_left[il], m_middle[il][ir]); const RI::Tensor m_lmr = RI::Tensor_Multiply::x0x1y0y1_x0x1a_y0y1a(m_lm, m_right[ir]); - m_proj -= m_lmr; + if(m_mul.empty()) + { m_mul = std::move(m_lmr); } + else + { m_mul += m_lmr; } } } - return m_proj; + return m_mul; } -RI::Tensor Exx_Opt_Orb::cal_proj_21( - const RI::Tensor & m_big, +RI::Tensor Exx_Opt_Orb::cal_mul_21( const std::vector> & m_left, const std::vector>> & m_middle, const std::vector> & m_right ) const { - ModuleBase::TITLE("Exx_Opt_Orb::cal_proj_21"); - RI::Tensor m_proj = m_big.copy(); + ModuleBase::TITLE("Exx_Opt_Orb::cal_mul_21"); + RI::Tensor m_mul; for( size_t il=0; il!=m_left.size(); ++il ) { for( size_t ir=0; ir!=m_right.size(); ++ir ) { - // m_proj = m_proj - m_left[il] * m_middle[il][ir] * m_right[ir].T; + // m_mul += m_left[il] * m_middle[il][ir] * m_right[ir].T; const RI::Tensor m_lm = RI::Tensor_Multiply::x0x1y1_x0x1a_ay1(m_left[il], m_middle[il][ir]); const RI::Tensor m_lmr = RI::Tensor_Multiply::x0x1y0_x0x1a_y0a(m_lm, m_right[ir]); - m_proj -= m_lmr; + if(m_mul.empty()) + { m_mul = std::move(m_lmr); } + else + { m_mul += m_lmr; } } } - return m_proj; + return m_mul; } -RI::Tensor Exx_Opt_Orb::cal_proj_12( - const RI::Tensor & m_big, +RI::Tensor Exx_Opt_Orb::cal_mul_12( const std::vector> & m_left, const std::vector>> & m_middle, const std::vector> & m_right ) const { - ModuleBase::TITLE("Exx_Opt_Orb::cal_proj_12"); - RI::Tensor m_proj = m_big.copy(); + ModuleBase::TITLE("Exx_Opt_Orb::cal_mul_12"); + RI::Tensor m_mul; for( size_t il=0; il!=m_left.size(); ++il ) { for( size_t ir=0; ir!=m_right.size(); ++ir ) { - // m_proj = m_proj - m_left[il] * m_middle[il][ir] * m_right[ir].T; + // m_mul += m_left[il] * m_middle[il][ir] * m_right[ir].T; const RI::Tensor m_lm = RI::Tensor_Multiply::x0y1_x0a_ay1(m_left[il], m_middle[il][ir]); const RI::Tensor m_lmr = RI::Tensor_Multiply::x0y0y1_x0a_y0y1a(m_lm, m_right[ir]); - m_proj -= m_lmr; + if(m_mul.empty()) + { m_mul = std::move(m_lmr); } + else + { m_mul += m_lmr; } } } - return m_proj; + return m_mul; } -RI::Tensor Exx_Opt_Orb::cal_proj_11( - const RI::Tensor & m_big, +RI::Tensor Exx_Opt_Orb::cal_mul_11( const std::vector> & m_left, const std::vector>> & m_middle, const std::vector> & m_right ) const { - ModuleBase::TITLE("Exx_Opt_Orb::cal_proj_11"); - RI::Tensor m_proj = m_big.copy(); + ModuleBase::TITLE("Exx_Opt_Orb::cal_mul_11"); + RI::Tensor m_mul; for( size_t il=0; il!=m_left.size(); ++il ) { for( size_t ir=0; ir!=m_right.size(); ++ir ) { - // m_proj = m_proj - m_left[il] * m_middle[il][ir] * m_right[ir].T; + // m_mul += m_left[il] * m_middle[il][ir] * m_right[ir].T; const RI::Tensor m_lm = RI::Tensor_Multiply::x0y1_x0a_ay1(m_left[il], m_middle[il][ir]); const RI::Tensor m_lmr = RI::Tensor_Multiply::x0y0_x0a_y0a(m_lm, m_right[ir]); - m_proj -= m_lmr; + if(m_mul.empty()) + { m_mul = std::move(m_lmr); } + else + { m_mul += m_lmr; } } } - return m_proj; + return m_mul; } std::vector>> Exx_Opt_Orb::cal_I( diff --git a/source/source_lcao/module_ri/exx_opt_orb.h b/source/source_lcao/module_ri/exx_opt_orb.h index 2520a584ba..7ccaf0854b 100644 --- a/source/source_lcao/module_ri/exx_opt_orb.h +++ b/source/source_lcao/module_ri/exx_opt_orb.h @@ -23,23 +23,19 @@ class Exx_Opt_Orb std::vector>> cal_I( const std::map>>>> &ms, const size_t TA, const size_t IA, const size_t TB, const size_t IB ) const; - RI::Tensor cal_proj_22( - const RI::Tensor & m_big, + RI::Tensor cal_mul_22( const std::vector> & m_left, const std::vector>> & m_middle, const std::vector> & m_right ) const; - RI::Tensor cal_proj_21( - const RI::Tensor & m_big, + RI::Tensor cal_mul_21( const std::vector> & m_left, const std::vector>> & m_middle, const std::vector> & m_right ) const; - RI::Tensor cal_proj_12( - const RI::Tensor & m_big, + RI::Tensor cal_mul_12( const std::vector> & m_left, const std::vector>> & m_middle, const std::vector> & m_right ) const; - RI::Tensor cal_proj_11( - const RI::Tensor & m_big, + RI::Tensor cal_mul_11( const std::vector> & m_left, const std::vector>> & m_middle, const std::vector> & m_right ) const; From 71bafa1e643166d3cb2b1b2532b050212839d0db Mon Sep 17 00:00:00 2001 From: linpz Date: Mon, 15 Sep 2025 20:46:03 +0800 Subject: [PATCH 2/2] 1. Improve output accuracy of Exx_Opt_Orb 2. Change default value of exx_opt_orb_tolerence --- docs/advanced/input_files/input-main.md | 2 +- source/source_hamilt/module_xc/exx_info.h | 2 +- source/source_io/module_parameter/input_parameter.h | 2 +- source/source_io/test/read_input_ptest.cpp | 2 +- source/source_io/test/support/INPUT | 2 +- source/source_lcao/module_ri/exx_opt_orb-print.cpp | 2 ++ 6 files changed, 7 insertions(+), 5 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index a2024dec66..c407df99d1 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -2977,7 +2977,7 @@ These variables are relevant when using hybrid functionals with *[basis_type](#b - **Type**: Real - **Availability**: *[calculation](#calculation)==gen_opt_abfs* - **Description**: The threshold when solving for the zeros of spherical Bessel functions. A reasonable choice is 1e-12. -- **Default**: 0 +- **Default**: 1E-12 ### exx_real_number diff --git a/source/source_hamilt/module_xc/exx_info.h b/source/source_hamilt/module_xc/exx_info.h index 0b427b438d..91c4652707 100644 --- a/source/source_hamilt/module_xc/exx_info.h +++ b/source/source_hamilt/module_xc/exx_info.h @@ -81,7 +81,7 @@ struct Exx_Info { int abfs_Lmax = 0; // tmp double ecut_exx = 60; - double tolerence = 1E-2; + double tolerence = 1E-12; double kmesh_times = 4; }; Exx_Info_Opt_ABFs info_opt_abfs; diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index db2b5d5a72..b540bfb895 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -546,7 +546,7 @@ struct Input_para ///< calculating Columb potential is to that of atomic orbitals int exx_opt_orb_lmax = 0; ///< the maximum l of the spherical Bessel functions for opt ABFs double exx_opt_orb_ecut = 0.0; ///< the cut-off of plane wave expansion for opt ABFs - double exx_opt_orb_tolerence = 0.0; ///< the threshold when solving for the zeros of spherical Bessel + double exx_opt_orb_tolerence = 1E-12; ///< the threshold when solving for the zeros of spherical Bessel ///< functions for opt ABFs bool exx_symmetry_realspace = true; ///< whether to reduce the real-space sector in when using symmetry=1 in EXX calculation diff --git a/source/source_io/test/read_input_ptest.cpp b/source/source_io/test/read_input_ptest.cpp index 14b0ef5a93..ebb276265d 100644 --- a/source/source_io/test/read_input_ptest.cpp +++ b/source/source_io/test/read_input_ptest.cpp @@ -296,7 +296,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_DOUBLE_EQ(param.inp.rpa_ccp_rmesh_times, 10.0); EXPECT_EQ(param.inp.exx_opt_orb_lmax, 0); EXPECT_DOUBLE_EQ(param.inp.exx_opt_orb_ecut, 0.0); - EXPECT_DOUBLE_EQ(param.inp.exx_opt_orb_tolerence, 0.0); + EXPECT_DOUBLE_EQ(param.inp.exx_opt_orb_tolerence, 1E-12); EXPECT_FALSE(param.inp.noncolin); EXPECT_FALSE(param.inp.lspinorb); EXPECT_DOUBLE_EQ(param.inp.soc_lambda, 1.0); diff --git a/source/source_io/test/support/INPUT b/source/source_io/test/support/INPUT index 33526863a8..90357fb14f 100644 --- a/source/source_io/test/support/INPUT +++ b/source/source_io/test/support/INPUT @@ -291,7 +291,7 @@ exx_v_grad_r_threshold 0 # exx_ccp_rmesh_times default # exx_opt_orb_lmax 0 # exx_opt_orb_ecut 0 # -exx_opt_orb_tolerence 0 # +exx_opt_orb_tolerence 1E-12 # #Parameters (16.tddft) td_vext 0 #add extern potential or not diff --git a/source/source_lcao/module_ri/exx_opt_orb-print.cpp b/source/source_lcao/module_ri/exx_opt_orb-print.cpp index a851d2c579..d3d5444ff4 100644 --- a/source/source_lcao/module_ri/exx_opt_orb-print.cpp +++ b/source/source_lcao/module_ri/exx_opt_orb-print.cpp @@ -1,6 +1,7 @@ #include "exx_opt_orb.h" #include "../../source_pw/module_pwdft/global.h" #include "exx_abfs-jle.h" +#include void Exx_Opt_Orb::print_matrix( const Exx_Info::Exx_Info_Opt_ABFs &info, @@ -204,6 +205,7 @@ void Exx_Opt_Orb::print_matrix( std::ofstream ofs(file_name+"_"+std::to_string(TA)+"_"+std::to_string(IA)+"_"+std::to_string(TB)+"_"+std::to_string(IB)); print_header(ofs); + ofs<