Refactor: the initialization scheme in LRI_CV

source/source_lcao/module_ri/ABFs_Construct-PCA.cpp

@@ -97,8 +97,12 @@ namespace PCA
 }
 
 		Matrix_Orbs21 m_abfslcaos_lcaos;
-		m_abfslcaos_lcaos.init( 1, ucell , orb, kmesh_times, 1 );
-		m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos );
+		ORB_gaunt_table MGT;
+		int Lmax;
+		m_abfslcaos_lcaos.init( 1, ucell , orb, kmesh_times, orb.get_Rmax(), Lmax );
+		MGT.init_Gaunt_CH(Lmax);
+        MGT.init_Gaunt(Lmax);
+		m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos, MGT );
 
 		std::map<std::size_t,std::map<std::size_t,std::set<double>>> delta_R;
 		for( std::size_t it=0; it!=abfs.size(); ++it ) {

source/source_lcao/module_ri/Exx_LRI.h

@@ -90,6 +90,7 @@ class Exx_LRI
 	const Exx_Info::Exx_Info_RI &info;
 	MPI_Comm mpi_comm;
 	const K_Vectors *p_kv = nullptr;
+	ORB_gaunt_table MGT;
 	std::vector<double> orb_cutoff_;
 
 	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> lcaos;

source/source_lcao/module_ri/Exx_LRI.hpp

@@ -53,12 +53,14 @@ void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in,
 
 	this->coulomb_settings = RI_Util::update_coulomb_settings(this->info.coulomb_param, ucell, this->p_kv);
 
+	bool init_MGT = true;
 	for(const auto &settings_list : this->coulomb_settings)
 	{
 		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);
 		this->exx_objs[settings_list.first].cv.set_orbitals(ucell, orb,
 															this->lcaos, this->abfs, this->exx_objs[settings_list.first].abfs_ccp,
-															this->info.kmesh_times, this->info.ccp_rmesh_times );
+															this->info.kmesh_times, this->MGT, init_MGT, settings_list.second.first );
+		init_MGT = false; // only init once
 	}
 
 	ModuleBase::timer::tick("Exx_LRI", "init");

source/source_lcao/module_ri/LRI_CV.h

@@ -40,7 +40,9 @@ class LRI_CV
 		const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs_in,
 		const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs_ccp_in,
 		const double &kmesh_times,
-		const double &ccp_rmesh_times_in);
+		ORB_gaunt_table& MGT,
+        const bool& init_MGT,
+        const bool& init_C);
 	inline std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
 	cal_Vs(
 		const UnitCell &ucell,
@@ -64,13 +66,13 @@ class LRI_CV
 	size_t get_index_abfs_size(const size_t &iat){return this->index_abfs[iat].count_size; }
 
 private:
-    std::vector<double> orb_cutoff_;
 	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> lcaos;
 	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs;
 	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs_ccp;
 	ModuleBase::Element_Basis_Index::IndexLNM index_lcaos;
 	ModuleBase::Element_Basis_Index::IndexLNM index_abfs;
-	double ccp_rmesh_times;
+	std::vector<double> lcaos_rcut;
+    std::vector<double> abfs_ccp_rcut;
 
 public:
 	std::map<int,std::map<int,std::map<Abfs::Vector3_Order<double>,RI::Tensor<Tdata>>>> Vws;
@@ -92,15 +94,19 @@ class LRI_CV
 		const int it1,
 		const Abfs::Vector3_Order<double> &R,
 		const std::map<std::string,bool> &flags)>;
+	using T_func_cal_Rcut = std::function<double(const int it0, const int it1)>;
 	template<typename Tresult>
 	std::map<TA,std::map<TAC,Tresult>>
 	cal_datas(
 		const UnitCell &ucell,
-		const std::vector<TA> &list_A0,
-		const std::vector<TAC> &list_A1,
-		const std::map<std::string,bool> &flags,
-		const double &rmesh_times,
-		const T_func_DPcal_data<Tresult> &func_DPcal_data);
+        const std::vector<TA>& list_A0,
+        const std::vector<TAC>& list_A1,
+        const std::map<std::string, bool>& flags,
+        const T_func_cal_Rcut& func_cal_Rcut,
+        const T_func_DPcal_data<Tresult>& func_DPcal_data);
+
+    inline double cal_V_Rcut(const int it0, const int it1);
+    inline double cal_C_Rcut(const int it0, const int it1);
 
 	inline RI::Tensor<Tdata>
 	DPcal_V(

source/source_lcao/module_ri/LRI_CV.hpp

@@ -9,6 +9,7 @@
 #include "LRI_CV.h"
 #include "LRI_CV_Tools.h"
 #include "exx_abfs-abfs_index.h"
+#include "exx_abfs-construct_orbs.h"
 #include "RI_Util.h"
 #include "../../source_base/tool_title.h"
 #include "../../source_base/timer.h"
@@ -43,16 +44,22 @@ void LRI_CV<Tdata>::set_orbitals(
 	const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs_in,
 	const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs_ccp_in,
 	const double &kmesh_times,
-	const double &ccp_rmesh_times_in)
+	ORB_gaunt_table& MGT,
+    const bool& init_MGT,
+    const bool& init_C)
 {
 	ModuleBase::TITLE("LRI_CV", "set_orbitals");
 	ModuleBase::timer::tick("LRI_CV", "set_orbitals");
 
-	this->orb_cutoff_ = orb.cutoffs();
 	this->lcaos = lcaos_in;
 	this->abfs = abfs_in;
 	this->abfs_ccp = abfs_ccp_in;
-	this->ccp_rmesh_times = ccp_rmesh_times_in;
+
+	this->lcaos_rcut = Exx_Abfs::Construct_Orbs::get_Rcut(this->lcaos);
+    this->abfs_ccp_rcut = Exx_Abfs::Construct_Orbs::get_Rcut(this->abfs_ccp);
+    const double lcaos_rmax = Exx_Abfs::Construct_Orbs::get_Rmax(this->lcaos);
+    const double abfs_ccp_rmax
+        = Exx_Abfs::Construct_Orbs::get_Rmax(this->abfs_ccp);
 
 	const ModuleBase::Element_Basis_Index::Range
 		range_lcaos = Exx_Abfs::Abfs_Index::construct_range( lcaos );
@@ -62,27 +69,50 @@ void LRI_CV<Tdata>::set_orbitals(
 		range_abfs = Exx_Abfs::Abfs_Index::construct_range( abfs );
 	this->index_abfs = ModuleBase::Element_Basis_Index::construct_index( range_abfs );
 
-	this->m_abfs_abfs.init( 2, ucell,orb, kmesh_times, (1+this->ccp_rmesh_times)/2.0 );
-	this->m_abfs_abfs.init_radial( this->abfs_ccp, this->abfs );
-	this->m_abfs_abfs.init_radial_table();
-
-	this->m_abfslcaos_lcaos.init( 1, ucell , orb, kmesh_times, 1 );
-	this->m_abfslcaos_lcaos.init_radial( this->abfs_ccp, this->lcaos, this->lcaos );
-	this->m_abfslcaos_lcaos.init_radial_table();
+    int Lmax_v = std::numeric_limits<double>::min();
+    this->m_abfs_abfs.init(2, ucell, orb, kmesh_times, lcaos_rmax + abfs_ccp_rmax, Lmax_v);
+    int Lmax_c = std::numeric_limits<double>::min();
+    if (init_C)
+        this->m_abfslcaos_lcaos.init(1, ucell, orb, kmesh_times, lcaos_rmax, Lmax_c);
+    int Lmax = std::max(Lmax_v, Lmax_c);
+
+    if (init_MGT) {
+        MGT.init_Gaunt_CH(Lmax);
+        MGT.init_Gaunt(Lmax);
+    }
+
+    this->m_abfs_abfs.init_radial(this->abfs_ccp, this->abfs, MGT);
+    this->m_abfs_abfs.init_radial_table();
+    if (init_C) {
+        this->m_abfslcaos_lcaos.init_radial(this->abfs_ccp,
+                                            this->lcaos,
+                                            this->lcaos,
+                                            MGT);
+        this->m_abfslcaos_lcaos.init_radial_table();
+    }
 
 	ModuleBase::timer::tick("LRI_CV", "set_orbitals");
 }
 
+template <typename Tdata>
+double LRI_CV<Tdata>::cal_V_Rcut(const int it0, const int it1) {
+    return this->abfs_ccp_rcut[it0] + this->lcaos_rcut[it1];
+}
 
+template <typename Tdata>
+double LRI_CV<Tdata>::cal_C_Rcut(const int it0, const int it1) {
+    return std::min(this->abfs_ccp_rcut[it0], this->lcaos_rcut[it0])
+           + this->lcaos_rcut[it1];
+}
 
 template<typename Tdata> template<typename Tresult>
 auto LRI_CV<Tdata>::cal_datas(
 	const UnitCell &ucell,
-	const std::vector<TA> &list_A0,
-	const std::vector<TAC> &list_A1,
-	const std::map<std::string,bool> &flags,
-	const double &rmesh_times,
-	const T_func_DPcal_data<Tresult> &func_DPcal_data)
+    const std::vector<TA>& list_A0,
+    const std::vector<TAC>& list_A1,
+    const std::map<std::string, bool>& flags,
+    const T_func_cal_Rcut& func_cal_Rcut,
+    const T_func_DPcal_data<Tresult>& func_DPcal_data)
 -> std::map<TA,std::map<TAC,Tresult>>
 {
 	ModuleBase::TITLE("LRI_CV","cal_datas");
@@ -104,9 +134,8 @@ auto LRI_CV<Tdata>::cal_datas(
 			const int ia1 = ucell.iat2ia[iat1];
 			const ModuleBase::Vector3<double> tau0 = ucell.atoms[it0].tau[ia0];
 			const ModuleBase::Vector3<double> tau1 = ucell.atoms[it1].tau[ia1];
-			const double Rcut = std::min(
-				orb_cutoff_[it0] * rmesh_times + orb_cutoff_[it1],
-				orb_cutoff_[it1] * rmesh_times + orb_cutoff_[it0]);
+			const double Rcut
+                = std::min(func_cal_Rcut(it0, it1), func_cal_Rcut(it1, it0));
 			const Abfs::Vector3_Order<double> R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*ucell.latvec);
 			if( R_delta.norm()*ucell.lat0 < Rcut )
 			{
@@ -137,7 +166,12 @@ auto LRI_CV<Tdata>::cal_Vs(
 		func_DPcal_V = std::bind(
 			&LRI_CV<Tdata>::DPcal_V, this,
 			std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
-	return this->cal_datas(ucell,list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_V);
+	const T_func_cal_Rcut func_cal_Rcut = std::bind(&LRI_CV<Tdata>::cal_V_Rcut,
+                                                    this,
+                                                    std::placeholders::_1,
+                                                    std::placeholders::_2);
+
+	return this->cal_datas(ucell,list_A0, list_A1, flags, func_cal_Rcut, func_DPcal_V);
 }
 
 template<typename Tdata>
@@ -153,7 +187,13 @@ auto LRI_CV<Tdata>::cal_dVs(
 		func_DPcal_dV = std::bind(
 			&LRI_CV<Tdata>::DPcal_dV, this,
 			std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
-		return this->cal_datas(ucell,list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_dV);
+
+	const T_func_cal_Rcut func_cal_Rcut = std::bind(&LRI_CV<Tdata>::cal_V_Rcut,
+                                                    this,
+                                                    std::placeholders::_1,
+                                                    std::placeholders::_2);
+
+	return this->cal_datas(ucell,list_A0, list_A1, flags, func_cal_Rcut, func_DPcal_dV);
 }
 
 template<typename Tdata>
@@ -171,8 +211,13 @@ auto LRI_CV<Tdata>::cal_Cs_dCs(
 		func_DPcal_C_dC = std::bind(
 			&LRI_CV<Tdata>::DPcal_C_dC, this,
 			std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
+    const T_func_cal_Rcut func_cal_Rcut = std::bind(&LRI_CV<Tdata>::cal_C_Rcut,
+                                                    this,
+                                                    std::placeholders::_1,
+                                                    std::placeholders::_2);
+
 	std::map<TA,std::map<TAC, std::pair<RI::Tensor<Tdata>, std::array<RI::Tensor<Tdata>,3>>>>
-		Cs_dCs_tmp = this->cal_datas(ucell,list_A0, list_A1, flags, std::min(1.0,this->ccp_rmesh_times), func_DPcal_C_dC);
+		Cs_dCs_tmp = this->cal_datas(ucell,list_A0, list_A1, flags, func_cal_Rcut, func_DPcal_C_dC);
 
 	std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Cs;
     std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>> dCs;

source/source_lcao/module_ri/Matrix_Orbs11.cpp

@@ -13,12 +13,13 @@ void Matrix_Orbs11::init(const int mode,
                          const UnitCell& ucell,
                          const LCAO_Orbitals& orb, 
                          const double kmesh_times, 
-                         const double rmesh_times)
+                         const double rmax, 
+                         int& Lmax)
 {
     ModuleBase::TITLE("Matrix_Orbs11", "init");
     ModuleBase::timer::tick("Matrix_Orbs11", "init");
 
-    int Lmax_used, Lmax;
+    int Lmax_used;
     this->lat0 = &ucell.lat0;
     const int ntype = orb.get_ntype();
     int lmax_orb = -1, lmax_beta = -1;
@@ -30,7 +31,7 @@ void Matrix_Orbs11::init(const int mode,
     const double dr = orb.get_dR();
     const double dk = orb.get_dk();
     const int kmesh = orb.get_kmesh() * kmesh_times + 1;
-    int Rmesh = static_cast<int>(orb.get_Rmax() * rmesh_times / dr) + 4;
+    int Rmesh = static_cast<int>(rmax / dr) + 4;
     Rmesh += 1 - Rmesh % 2;
 
     Center2_Orb::init_Table_Spherical_Bessel(2,
@@ -49,14 +50,15 @@ void Matrix_Orbs11::init(const int mode,
     //=========================================
     // (3) make Gaunt coefficients table
     //=========================================
-    this->MGT.init_Gaunt_CH(Lmax);
-    this->MGT.init_Gaunt(Lmax);
+    // this->MGT.init_Gaunt_CH(Lmax);
+    // this->MGT.init_Gaunt(Lmax);
 
     ModuleBase::timer::tick("Matrix_Orbs11", "init");
 }
 
 void Matrix_Orbs11::init_radial(const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_A,
-                                const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B)
+                                const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B,
+                                const ORB_gaunt_table& MGT)
 {
     ModuleBase::TITLE("Matrix_Orbs11", "init_radial");
     ModuleBase::timer::tick("Matrix_Orbs11", "init_radial");
@@ -68,7 +70,7 @@ void Matrix_Orbs11::init_radial(const std::vector<std::vector<std::vector<Numeri
                         for (size_t NB = 0; NB != orb_B[TB][LB].size(); ++NB) {
                             center2_orb11_s[TA][TB][LA][NA][LB].insert(std::make_pair(
                                 NB,
-                                Center2_Orb::Orb11(orb_A[TA][LA][NA], orb_B[TB][LB][NB], psb_, this->MGT)));
+                                Center2_Orb::Orb11(orb_A[TA][LA][NA], orb_B[TB][LB][NB], psb_, MGT)));
                         }
                     }
                 }
@@ -78,7 +80,7 @@ void Matrix_Orbs11::init_radial(const std::vector<std::vector<std::vector<Numeri
     ModuleBase::timer::tick("Matrix_Orbs11", "init_radial");
 }
 
-void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B)
+void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B, const ORB_gaunt_table& MGT)
 {
     ModuleBase::TITLE("Matrix_Orbs11", "init_radial");
     ModuleBase::timer::tick("Matrix_Orbs11", "init_radial");
@@ -93,7 +95,7 @@ void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals&
                                                Center2_Orb::Orb11(orb_A.Phi[TA].PhiLN(LA, NA),
                                                                   orb_B.Phi[TB].PhiLN(LB, NB),
                                                                   psb_,
-                                                                  this->MGT)));
+                                                                  MGT)));
                         }
                     }
                 }

source/source_lcao/module_ri/Matrix_Orbs11.h

@@ -28,11 +28,13 @@ class Matrix_Orbs11
               const UnitCell& ucell,
               const LCAO_Orbitals& orb,
               const double kmesh_times,  // extend Kcut, keep dK
-              const double rmesh_times); // extend Rcut, keep dR
+              const double rmax,
+              int& Lmax); 
 
     void init_radial(const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_A,
-                     const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B);
-    void init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B);
+                     const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B,
+                     const ORB_gaunt_table& MGT);
+    void init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B, const ORB_gaunt_table& MGT);
 
     void init_radial_table();
     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
 
   private:
     ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr;
-    ORB_gaunt_table MGT;
     const double lcao_dr_ = 0.01;
     double* lat0=nullptr;                                         // restore ucell.lat0
     std::map<size_t,                                              // TA