Skip to content

Commit 3a842b0

Browse files
1041176461jiyuangQianruipkuESROAMERmaki49
authored
Refactor: the initialization scheme in LRI_CV (#6363)
* Use template to reconstruct parse_expression * Feature: output R matrix at each MD step * Modify'matrix_HS' to 'matrix' for R matrix output * Merge branches 'develop' and 'develop' of https://github.com/1041176461/abacus-develop into develop * Fix: modify index in parse_expression * Fix: add regfree for parse_expression * Doc: update phonopy doc * Doc: update phonopy doc * fix tdos plot for nspin=2 * optimize dosplot for nspin=2 * fix legend for dosplot * Add files via upload * Update cal_edm_tddft.cpp * Refactor: modify exx loop for relax/md * Update result.ref * Fix wrong parameters in integrate test * Update result.ref * Update result.ref * Update result.ref * Update result.ref * Update result.ref * Update result.ref * Update result.ref * Update result.ref * add exx when istep>0 * Update dos.py * Update esolver_sdft_pw.cpp * Update lcao_before_scf.cpp * Update Exx_LRI_interface.h * Update Exx_LRI_interface.hpp * Fix: compile error * Fix: compile error * Fix: change HSE relax/md result.ref for new framework * Fix: compile error * compatible with exx_iter_finish * Add files via upload * Update esolver_ks_lcao_tddft.h * Update esolver_ks_lcao_tddft.cpp * Fix: support negative value in parse_expression * [pre-commit.ci lite] apply automatic fixes * Fix: parse_expression for scientific notation * Update input_conv.h * add complex erf function * refactor ccp framework * fix compile * fix compile * fix compile * fix compile * fix compile * fix compile * fix compile * fix empty error * add revised spencer * add cal_psi_erfc_spencer * add exx_singularity_correction * add exx_singularity_correction test * add lc-corrected hybrid functionals * fix compile * fix compile * fix compile error * fix libxc error * fix libxc error * fix libxc error * add klist.h in RI_Util.h * fix include * fix ptest for hse * fix lc_pbe name for vdw autoset * fix test error * fix cam_pbeh test * fix fock settings * support rsh functionals * fix libxc param * add input * fix input * refactor lri_cv * fix compile * fix compile --------- Co-authored-by: jiyuang <jiyuyang@mail.ustc.com> Co-authored-by: Qianrui <76200646+Qianruipku@users.noreply.github.com> Co-authored-by: HTZhao <104255052+ESROAMER@users.noreply.github.com> Co-authored-by: maki49 <1579492865@qq.com> Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>
1 parent c873197 commit 3a842b0

15 files changed

+281
-110
lines changed

‎source/source_lcao/module_ri/ABFs_Construct-PCA.cpp

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -97,8 +97,12 @@ namespace PCA
9797
}
9898

9999
Matrix_Orbs21 m_abfslcaos_lcaos;
100-
m_abfslcaos_lcaos.init( 1, ucell , orb, kmesh_times, 1 );
101-
m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos );
100+
ORB_gaunt_table MGT;
101+
int Lmax;
102+
m_abfslcaos_lcaos.init( 1, ucell , orb, kmesh_times, orb.get_Rmax(), Lmax );
103+
MGT.init_Gaunt_CH(Lmax);
104+
MGT.init_Gaunt(Lmax);
105+
m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos, MGT );
102106

103107
std::map<std::size_t,std::map<std::size_t,std::set<double>>> delta_R;
104108
for( std::size_t it=0; it!=abfs.size(); ++it ) {

‎source/source_lcao/module_ri/Exx_LRI.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -90,6 +90,7 @@ class Exx_LRI
9090
const Exx_Info::Exx_Info_RI &info;
9191
MPI_Comm mpi_comm;
9292
const K_Vectors *p_kv = nullptr;
93+
ORB_gaunt_table MGT;
9394
std::vector<double> orb_cutoff_;
9495

9596
std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> lcaos;

‎source/source_lcao/module_ri/Exx_LRI.hpp

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -53,12 +53,14 @@ void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in,
5353

5454
this->coulomb_settings = RI_Util::update_coulomb_settings(this->info.coulomb_param, ucell, this->p_kv);
5555

56+
bool init_MGT = true;
5657
for(const auto &settings_list : this->coulomb_settings)
5758
{
5859
this->exx_objs[settings_list.first].abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, settings_list.second.second, this->info.ccp_rmesh_times);
5960
this->exx_objs[settings_list.first].cv.set_orbitals(ucell, orb,
6061
this->lcaos, this->abfs, this->exx_objs[settings_list.first].abfs_ccp,
61-
this->info.kmesh_times, this->info.ccp_rmesh_times );
62+
this->info.kmesh_times, this->MGT, init_MGT, settings_list.second.first );
63+
init_MGT = false; // only init once
6264
}
6365

6466
ModuleBase::timer::tick("Exx_LRI", "init");

‎source/source_lcao/module_ri/LRI_CV.h

Lines changed: 14 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,9 @@ class LRI_CV
4040
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs_in,
4141
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs_ccp_in,
4242
const double &kmesh_times,
43-
const double &ccp_rmesh_times_in);
43+
ORB_gaunt_table& MGT,
44+
const bool& init_MGT,
45+
const bool& init_C);
4446
inline std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
4547
cal_Vs(
4648
const UnitCell &ucell,
@@ -64,13 +66,13 @@ class LRI_CV
6466
size_t get_index_abfs_size(const size_t &iat){return this->index_abfs[iat].count_size; }
6567

6668
private:
67-
std::vector<double> orb_cutoff_;
6869
std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> lcaos;
6970
std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs;
7071
std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs_ccp;
7172
ModuleBase::Element_Basis_Index::IndexLNM index_lcaos;
7273
ModuleBase::Element_Basis_Index::IndexLNM index_abfs;
73-
double ccp_rmesh_times;
74+
std::vector<double> lcaos_rcut;
75+
std::vector<double> abfs_ccp_rcut;
7476

7577
public:
7678
std::map<int,std::map<int,std::map<Abfs::Vector3_Order<double>,RI::Tensor<Tdata>>>> Vws;
@@ -92,15 +94,19 @@ class LRI_CV
9294
const int it1,
9395
const Abfs::Vector3_Order<double> &R,
9496
const std::map<std::string,bool> &flags)>;
97+
using T_func_cal_Rcut = std::function<double(const int it0, const int it1)>;
9598
template<typename Tresult>
9699
std::map<TA,std::map<TAC,Tresult>>
97100
cal_datas(
98101
const UnitCell &ucell,
99-
const std::vector<TA> &list_A0,
100-
const std::vector<TAC> &list_A1,
101-
const std::map<std::string,bool> &flags,
102-
const double &rmesh_times,
103-
const T_func_DPcal_data<Tresult> &func_DPcal_data);
102+
const std::vector<TA>& list_A0,
103+
const std::vector<TAC>& list_A1,
104+
const std::map<std::string, bool>& flags,
105+
const T_func_cal_Rcut& func_cal_Rcut,
106+
const T_func_DPcal_data<Tresult>& func_DPcal_data);
107+
108+
inline double cal_V_Rcut(const int it0, const int it1);
109+
inline double cal_C_Rcut(const int it0, const int it1);
104110

105111
inline RI::Tensor<Tdata>
106112
DPcal_V(

‎source/source_lcao/module_ri/LRI_CV.hpp

Lines changed: 66 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
#include "LRI_CV.h"
1010
#include "LRI_CV_Tools.h"
1111
#include "exx_abfs-abfs_index.h"
12+
#include "exx_abfs-construct_orbs.h"
1213
#include "RI_Util.h"
1314
#include "../../source_base/tool_title.h"
1415
#include "../../source_base/timer.h"
@@ -43,16 +44,22 @@ void LRI_CV<Tdata>::set_orbitals(
4344
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs_in,
4445
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs_ccp_in,
4546
const double &kmesh_times,
46-
const double &ccp_rmesh_times_in)
47+
ORB_gaunt_table& MGT,
48+
const bool& init_MGT,
49+
const bool& init_C)
4750
{
4851
ModuleBase::TITLE("LRI_CV", "set_orbitals");
4952
ModuleBase::timer::tick("LRI_CV", "set_orbitals");
5053

51-
this->orb_cutoff_ = orb.cutoffs();
5254
this->lcaos = lcaos_in;
5355
this->abfs = abfs_in;
5456
this->abfs_ccp = abfs_ccp_in;
55-
this->ccp_rmesh_times = ccp_rmesh_times_in;
57+
58+
this->lcaos_rcut = Exx_Abfs::Construct_Orbs::get_Rcut(this->lcaos);
59+
this->abfs_ccp_rcut = Exx_Abfs::Construct_Orbs::get_Rcut(this->abfs_ccp);
60+
const double lcaos_rmax = Exx_Abfs::Construct_Orbs::get_Rmax(this->lcaos);
61+
const double abfs_ccp_rmax
62+
= Exx_Abfs::Construct_Orbs::get_Rmax(this->abfs_ccp);
5663

5764
const ModuleBase::Element_Basis_Index::Range
5865
range_lcaos = Exx_Abfs::Abfs_Index::construct_range( lcaos );
@@ -62,27 +69,50 @@ void LRI_CV<Tdata>::set_orbitals(
6269
range_abfs = Exx_Abfs::Abfs_Index::construct_range( abfs );
6370
this->index_abfs = ModuleBase::Element_Basis_Index::construct_index( range_abfs );
6471

65-
this->m_abfs_abfs.init( 2, ucell,orb, kmesh_times, (1+this->ccp_rmesh_times)/2.0 );
66-
this->m_abfs_abfs.init_radial( this->abfs_ccp, this->abfs );
67-
this->m_abfs_abfs.init_radial_table();
68-
69-
this->m_abfslcaos_lcaos.init( 1, ucell , orb, kmesh_times, 1 );
70-
this->m_abfslcaos_lcaos.init_radial( this->abfs_ccp, this->lcaos, this->lcaos );
71-
this->m_abfslcaos_lcaos.init_radial_table();
72+
int Lmax_v = std::numeric_limits<double>::min();
73+
this->m_abfs_abfs.init(2, ucell, orb, kmesh_times, lcaos_rmax + abfs_ccp_rmax, Lmax_v);
74+
int Lmax_c = std::numeric_limits<double>::min();
75+
if (init_C)
76+
this->m_abfslcaos_lcaos.init(1, ucell, orb, kmesh_times, lcaos_rmax, Lmax_c);
77+
int Lmax = std::max(Lmax_v, Lmax_c);
78+
79+
if (init_MGT) {
80+
MGT.init_Gaunt_CH(Lmax);
81+
MGT.init_Gaunt(Lmax);
82+
}
83+
84+
this->m_abfs_abfs.init_radial(this->abfs_ccp, this->abfs, MGT);
85+
this->m_abfs_abfs.init_radial_table();
86+
if (init_C) {
87+
this->m_abfslcaos_lcaos.init_radial(this->abfs_ccp,
88+
this->lcaos,
89+
this->lcaos,
90+
MGT);
91+
this->m_abfslcaos_lcaos.init_radial_table();
92+
}
7293

7394
ModuleBase::timer::tick("LRI_CV", "set_orbitals");
7495
}
7596

97+
template <typename Tdata>
98+
double LRI_CV<Tdata>::cal_V_Rcut(const int it0, const int it1) {
99+
return this->abfs_ccp_rcut[it0] + this->lcaos_rcut[it1];
100+
}
76101

102+
template <typename Tdata>
103+
double LRI_CV<Tdata>::cal_C_Rcut(const int it0, const int it1) {
104+
return std::min(this->abfs_ccp_rcut[it0], this->lcaos_rcut[it0])
105+
+ this->lcaos_rcut[it1];
106+
}
77107

78108
template<typename Tdata> template<typename Tresult>
79109
auto LRI_CV<Tdata>::cal_datas(
80110
const UnitCell &ucell,
81-
const std::vector<TA> &list_A0,
82-
const std::vector<TAC> &list_A1,
83-
const std::map<std::string,bool> &flags,
84-
const double &rmesh_times,
85-
const T_func_DPcal_data<Tresult> &func_DPcal_data)
111+
const std::vector<TA>& list_A0,
112+
const std::vector<TAC>& list_A1,
113+
const std::map<std::string, bool>& flags,
114+
const T_func_cal_Rcut& func_cal_Rcut,
115+
const T_func_DPcal_data<Tresult>& func_DPcal_data)
86116
-> std::map<TA,std::map<TAC,Tresult>>
87117
{
88118
ModuleBase::TITLE("LRI_CV","cal_datas");
@@ -104,9 +134,8 @@ auto LRI_CV<Tdata>::cal_datas(
104134
const int ia1 = ucell.iat2ia[iat1];
105135
const ModuleBase::Vector3<double> tau0 = ucell.atoms[it0].tau[ia0];
106136
const ModuleBase::Vector3<double> tau1 = ucell.atoms[it1].tau[ia1];
107-
const double Rcut = std::min(
108-
orb_cutoff_[it0] * rmesh_times + orb_cutoff_[it1],
109-
orb_cutoff_[it1] * rmesh_times + orb_cutoff_[it0]);
137+
const double Rcut
138+
= std::min(func_cal_Rcut(it0, it1), func_cal_Rcut(it1, it0));
110139
const Abfs::Vector3_Order<double> R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*ucell.latvec);
111140
if( R_delta.norm()*ucell.lat0 < Rcut )
112141
{
@@ -137,7 +166,12 @@ auto LRI_CV<Tdata>::cal_Vs(
137166
func_DPcal_V = std::bind(
138167
&LRI_CV<Tdata>::DPcal_V, this,
139168
std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
140-
return this->cal_datas(ucell,list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_V);
169+
const T_func_cal_Rcut func_cal_Rcut = std::bind(&LRI_CV<Tdata>::cal_V_Rcut,
170+
this,
171+
std::placeholders::_1,
172+
std::placeholders::_2);
173+
174+
return this->cal_datas(ucell,list_A0, list_A1, flags, func_cal_Rcut, func_DPcal_V);
141175
}
142176

143177
template<typename Tdata>
@@ -153,7 +187,13 @@ auto LRI_CV<Tdata>::cal_dVs(
153187
func_DPcal_dV = std::bind(
154188
&LRI_CV<Tdata>::DPcal_dV, this,
155189
std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
156-
return this->cal_datas(ucell,list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_dV);
190+
191+
const T_func_cal_Rcut func_cal_Rcut = std::bind(&LRI_CV<Tdata>::cal_V_Rcut,
192+
this,
193+
std::placeholders::_1,
194+
std::placeholders::_2);
195+
196+
return this->cal_datas(ucell,list_A0, list_A1, flags, func_cal_Rcut, func_DPcal_dV);
157197
}
158198

159199
template<typename Tdata>
@@ -171,8 +211,13 @@ auto LRI_CV<Tdata>::cal_Cs_dCs(
171211
func_DPcal_C_dC = std::bind(
172212
&LRI_CV<Tdata>::DPcal_C_dC, this,
173213
std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
214+
const T_func_cal_Rcut func_cal_Rcut = std::bind(&LRI_CV<Tdata>::cal_C_Rcut,
215+
this,
216+
std::placeholders::_1,
217+
std::placeholders::_2);
218+
174219
std::map<TA,std::map<TAC, std::pair<RI::Tensor<Tdata>, std::array<RI::Tensor<Tdata>,3>>>>
175-
Cs_dCs_tmp = this->cal_datas(ucell,list_A0, list_A1, flags, std::min(1.0,this->ccp_rmesh_times), func_DPcal_C_dC);
220+
Cs_dCs_tmp = this->cal_datas(ucell,list_A0, list_A1, flags, func_cal_Rcut, func_DPcal_C_dC);
176221

177222
std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Cs;
178223
std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>> dCs;

‎source/source_lcao/module_ri/Matrix_Orbs11.cpp

Lines changed: 11 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -13,12 +13,13 @@ void Matrix_Orbs11::init(const int mode,
1313
const UnitCell& ucell,
1414
const LCAO_Orbitals& orb,
1515
const double kmesh_times,
16-
const double rmesh_times)
16+
const double rmax,
17+
int& Lmax)
1718
{
1819
ModuleBase::TITLE("Matrix_Orbs11", "init");
1920
ModuleBase::timer::tick("Matrix_Orbs11", "init");
2021

21-
int Lmax_used, Lmax;
22+
int Lmax_used;
2223
this->lat0 = &ucell.lat0;
2324
const int ntype = orb.get_ntype();
2425
int lmax_orb = -1, lmax_beta = -1;
@@ -30,7 +31,7 @@ void Matrix_Orbs11::init(const int mode,
3031
const double dr = orb.get_dR();
3132
const double dk = orb.get_dk();
3233
const int kmesh = orb.get_kmesh() * kmesh_times + 1;
33-
int Rmesh = static_cast<int>(orb.get_Rmax() * rmesh_times / dr) + 4;
34+
int Rmesh = static_cast<int>(rmax / dr) + 4;
3435
Rmesh += 1 - Rmesh % 2;
3536

3637
Center2_Orb::init_Table_Spherical_Bessel(2,
@@ -49,14 +50,15 @@ void Matrix_Orbs11::init(const int mode,
4950
//=========================================
5051
// (3) make Gaunt coefficients table
5152
//=========================================
52-
this->MGT.init_Gaunt_CH(Lmax);
53-
this->MGT.init_Gaunt(Lmax);
53+
// this->MGT.init_Gaunt_CH(Lmax);
54+
// this->MGT.init_Gaunt(Lmax);
5455

5556
ModuleBase::timer::tick("Matrix_Orbs11", "init");
5657
}
5758

5859
void Matrix_Orbs11::init_radial(const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_A,
59-
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B)
60+
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B,
61+
const ORB_gaunt_table& MGT)
6062
{
6163
ModuleBase::TITLE("Matrix_Orbs11", "init_radial");
6264
ModuleBase::timer::tick("Matrix_Orbs11", "init_radial");
@@ -68,7 +70,7 @@ void Matrix_Orbs11::init_radial(const std::vector<std::vector<std::vector<Numeri
6870
for (size_t NB = 0; NB != orb_B[TB][LB].size(); ++NB) {
6971
center2_orb11_s[TA][TB][LA][NA][LB].insert(std::make_pair(
7072
NB,
71-
Center2_Orb::Orb11(orb_A[TA][LA][NA], orb_B[TB][LB][NB], psb_, this->MGT)));
73+
Center2_Orb::Orb11(orb_A[TA][LA][NA], orb_B[TB][LB][NB], psb_, MGT)));
7274
}
7375
}
7476
}
@@ -78,7 +80,7 @@ void Matrix_Orbs11::init_radial(const std::vector<std::vector<std::vector<Numeri
7880
ModuleBase::timer::tick("Matrix_Orbs11", "init_radial");
7981
}
8082

81-
void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B)
83+
void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B, const ORB_gaunt_table& MGT)
8284
{
8385
ModuleBase::TITLE("Matrix_Orbs11", "init_radial");
8486
ModuleBase::timer::tick("Matrix_Orbs11", "init_radial");
@@ -93,7 +95,7 @@ void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals&
9395
Center2_Orb::Orb11(orb_A.Phi[TA].PhiLN(LA, NA),
9496
orb_B.Phi[TB].PhiLN(LB, NB),
9597
psb_,
96-
this->MGT)));
98+
MGT)));
9799
}
98100
}
99101
}

‎source/source_lcao/module_ri/Matrix_Orbs11.h

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -28,11 +28,13 @@ class Matrix_Orbs11
2828
const UnitCell& ucell,
2929
const LCAO_Orbitals& orb,
3030
const double kmesh_times, // extend Kcut, keep dK
31-
const double rmesh_times); // extend Rcut, keep dR
31+
const double rmax,
32+
int& Lmax);
3233

3334
void init_radial(const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_A,
34-
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B);
35-
void init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B);
35+
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B,
36+
const ORB_gaunt_table& MGT);
37+
void init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B, const ORB_gaunt_table& MGT);
3638

3739
void init_radial_table();
3840
void init_radial_table(const std::map<size_t, std::map<size_t, std::set<double>>>& Rs); // unit: ucell.lat0
@@ -69,7 +71,6 @@ class Matrix_Orbs11
6971

7072
private:
7173
ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr;
72-
ORB_gaunt_table MGT;
7374
const double lcao_dr_ = 0.01;
7475
double* lat0=nullptr; // restore ucell.lat0
7576
std::map<size_t, // TA

0 commit comments

Comments
 (0)