Simplify Lmax and ORB_gaunt_table in Center2_Orb and Matrix_Orbs

source/source_basis/module_ao/ORB_gaunt_table.cpp

@@ -20,6 +20,7 @@ void ORB_gaunt_table::init_Gaunt(const int &lmax)
     ModuleBase::TITLE("ORB_gaunt_table", "init_Gaunt");
     ModuleBase::timer::tick("ORB_gaunt_table", "init_Gaunt");
 
+	this->Lmax_Gaunt_Coefficients = lmax;
 	const int nlm = (lmax * 2 + 1) * (lmax * 2 + 1);
 	this->Gaunt_Coefficients.create(nlm, nlm, nlm);
 
@@ -184,6 +185,7 @@ void ORB_gaunt_table::init_Gaunt_CH(const int& Lmax)
 	ModuleBase::TITLE("ORB_gaunt_table","init_Gaunt_CH");
 	ModuleBase::timer::tick("ORB_gaunt_table","init_Gaunt_CH");
 
+	this->Lmax_Gaunt_CH = Lmax;
 	int L = 2*Lmax + 1;
 	int Eff_Np = this->EP_EL(L);
 

source/source_basis/module_ao/ORB_gaunt_table.h

@@ -93,6 +93,9 @@ class ORB_gaunt_table
 
 	static int Index_M(const int& m);
 
+	int get_Lmax_Gaunt_Coefficients() const { return Lmax_Gaunt_Coefficients; }
+	int get_Lmax_Gaunt_CH()           const { return Lmax_Gaunt_CH; }
+
 	private:
 
 	// Index Function
@@ -122,5 +125,8 @@ class ORB_gaunt_table
 
 	//direct integral
 	ModuleBase::matrix Ylm_Gaunt;
+
+	int Lmax_Gaunt_Coefficients = -1;
+	int Lmax_Gaunt_CH = -1;
 };
 #endif

source/source_hamilt/module_xc/exx_info.h

@@ -69,8 +69,6 @@ struct Exx_Info
         double kmesh_times = 4;
         double Cs_inv_thr = -1;
 
-        int abfs_Lmax = 0; // tmp
-
         Exx_Info_RI(const Exx_Info::Exx_Info_Global& info_global)
             : coulomb_param(info_global.coulomb_param)
         {
@@ -80,7 +78,7 @@ struct Exx_Info
 
     struct Exx_Info_Opt_ABFs
     {
-        int abfs_Lmax = 0; // tmp
+        int abfs_Lmax = 0;
         double ecut_exx = 60;
         double tolerence = 1E-12;
         std::vector<std::string> files_jles;

source/source_lcao/LCAO_init_basis.cpp

@@ -62,11 +62,6 @@ void init_basis_lcao(Parallel_Orbitals& pv,
         two_center_bundle.build_beta(ucell.ntype, ucell.infoNL.Beta);
     }
 
-    int Lmax = 0;
-#ifdef __EXX
-    Lmax = GlobalC::exx_info.info_ri.abfs_Lmax;
-#endif
-
 #ifdef USE_NEW_TWO_CENTER
     two_center_bundle.tabulate();
 #else

source/source_lcao/center2_orb.cpp

@@ -49,8 +49,7 @@ std::pair<int,int> Center2_Orb::init_Lmax_2_2(const int& lmax_exx)
 // used in berryphase by jingan
 std::pair<int,int> Center2_Orb::init_Lmax_2_3(const int lmax_orb)
 {
-    int Lmax = std::max(-1, lmax_orb);
-    Lmax++;
+    const int Lmax = std::max(-1, lmax_orb) + 1;
     const int Lmax_used = 2 * Lmax + 1;
     assert(Lmax_used >= 1);
     return {Lmax_used, Lmax};
@@ -172,15 +171,14 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
 
     // double* integrated_func = new double[kmesh];
 
-    int ll = 0;
-    if (l != 0)
-    {
-        ll = l - 1;
-    }
-
-    const std::vector<std::vector<double>>& jlm1 = psb->get_jlx()[ll];
-    const std::vector<std::vector<double>>& jl = psb->get_jlx()[l];
-    const std::vector<std::vector<double>>& jlp1 = psb->get_jlx()[l + 1];
+    assert(psb->get_jlx().size()>=l+2);
+    const int lml = (l>0) ? (l-1) : 0;
+    const std::vector<std::vector<double>>& jlm1 = psb->get_jlx().at(lml);
+    const std::vector<std::vector<double>>& jl = psb->get_jlx().at(l);
+    const std::vector<std::vector<double>>& jlp1 = psb->get_jlx().at(l+1);
+    assert(jlm1.size()>=rmesh);
+    assert(jl.size()>=rmesh);
+    assert(jlp1.size()>=rmesh);
 
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static)
@@ -189,6 +187,7 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
     {
         std::vector<double> integrated_func(kmesh);
         const std::vector<double>& jl_r = jl[ir];
+        assert(jl_r.size()>=kmesh);
         for (int ik = 0; ik < kmesh; ++ik)
         {
             integrated_func[ik] = jl_r[ik] * k1_dot_k2[ik];
@@ -202,6 +201,8 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
         // Peize Lin accelerate 2017-10-02
         const std::vector<double>& jlm1_r = jlm1[ir];
         const std::vector<double>& jlp1_r = jlp1[ir];
+        assert(jlm1_r.size()>=kmesh);
+        assert(jlp1_r.size()>=kmesh);
         const double fac = l / (l + 1.0);
         if (l == 0)
         {
@@ -305,10 +306,11 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
 
     std::vector<double> integrated_func(kmesh);
 
-    const int lm1 = (l > 0 ? l - 1 : 0);
-    const std::vector<std::vector<double>>& jlm1 = psb->get_jlx()[lm1];
-    const std::vector<std::vector<double>>& jl = psb->get_jlx()[l];
-    const std::vector<std::vector<double>>& jlp1 = psb->get_jlx()[l + 1];
+    assert(psb->get_jlx().size()>=l+2);
+    const int lm1 = (l>0) ? (l-1) : 0;
+    const std::vector<std::vector<double>>& jlm1 = psb->get_jlx().at(lm1);
+    const std::vector<std::vector<double>>& jl = psb->get_jlx().at(l);
+    const std::vector<std::vector<double>>& jlp1 = psb->get_jlx().at(l+1);
 
     for (const size_t& ir: radials)
     {
@@ -320,6 +322,7 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
         }
 
         const std::vector<double>& jl_r = jl[ir];
+        assert(jl_r.size()>=kmesh);
         for (int ik = 0; ik < kmesh; ++ik)
         {
             integrated_func[ik] = jl_r[ik] * k1_dot_k2[ik];
@@ -329,8 +332,10 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
         ModuleBase::Integral::Simpson_Integral(kmesh, ModuleBase::GlobalFunc::VECTOR_TO_PTR(integrated_func), dk, temp);
         rs[ir] = temp * ModuleBase::FOUR_PI;
 
-        const std::vector<double>& jlm1_r = jlm1[ir];
-        const std::vector<double>& jlp1_r = jlp1[ir];
+        const std::vector<double>& jlm1_r = jlm1.at(ir);
+        const std::vector<double>& jlp1_r = jlp1.at(ir);
+        assert(jlm1_r.size()>=kmesh);
+        assert(jlp1_r.size()>=kmesh);
         const double fac = l / (l + 1.0);
         if (l == 0)
         {

source/source_lcao/module_ri/ABFs_Construct-PCA.cpp

@@ -172,29 +172,13 @@ RI::Tensor<double> get_column_mean0_matrix(const RI::Tensor<double>& m)
     const ModuleBase::Element_Basis_Index::IndexLNM index_abfs
         = ModuleBase::Element_Basis_Index::construct_index(range_abfs);
 
-		const int Lmax_bak = GlobalC::exx_info.info_ri.abfs_Lmax;
-		GlobalC::exx_info.info_ri.abfs_Lmax = std::numeric_limits<int>::min();
-    for (std::size_t T = 0; T != abfs.size(); ++T)
-    {
-        GlobalC::exx_info.info_ri.abfs_Lmax
-            = std::max(GlobalC::exx_info.info_ri.abfs_Lmax, static_cast<int>(abfs[T].size()) - 1);
-}
-
-		Matrix_Orbs21 m_abfslcaos_lcaos;
-		ORB_gaunt_table MGT;
-    const int Lmax = m_abfslcaos_lcaos.init(1, ucell, orb, kmesh_times, orb.get_Rmax());
-		MGT.init_Gaunt_CH(Lmax);
-    MGT.init_Gaunt(Lmax);
-    m_abfslcaos_lcaos.init_radial(abfs, lcaos, lcaos, MGT);
+	Matrix_Orbs21 m_abfslcaos_lcaos;
+    m_abfslcaos_lcaos.init(abfs, lcaos, lcaos, ucell, orb, kmesh_times, orb.get_Rmax());
 
     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)
-    {
-			delta_R[it][it] = {0.0};
-}
+        { delta_R[it][it] = {0.0}; }
 		m_abfslcaos_lcaos.init_radial_table(delta_R);
-
-		GlobalC::exx_info.info_ri.abfs_Lmax = Lmax_bak;
 
     std::vector<std::vector<std::pair<std::vector<double>, RI::Tensor<double>>>> eig(abfs.size());
     for (std::size_t T = 0; T != abfs.size(); ++T)

source/source_lcao/module_ri/Exx_LRI.h

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

@@ -39,28 +39,26 @@ void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in,
 	this->orb_cutoff_ = orb.cutoffs();
 
 	this->lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->info.kmesh_times );
+	Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->lcaos);
 
 	const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>
 		abfs_same_atom = Exx_Abfs::Construct_Orbs::abfs_same_atom(ucell, orb, this->lcaos, this->info.kmesh_times, this->info.pca_threshold );
 	if(this->info.files_abfs.empty())
 		{ this->abfs = abfs_same_atom;}
 	else
 		{ this->abfs = Exx_Abfs::IO::construct_abfs( abfs_same_atom, orb, this->info.files_abfs, this->info.kmesh_times ); 	}
+	Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->abfs);
 	Exx_Abfs::Construct_Orbs::print_orbs_size(ucell, this->abfs, GlobalV::ofs_running);
 
-	for( size_t T=0; T!=this->abfs.size(); ++T )
-		{ GlobalC::exx_info.info_ri.abfs_Lmax = std::max( GlobalC::exx_info.info_ri.abfs_Lmax, static_cast<int>(this->abfs[T].size())-1 ); }
-
 	this->coulomb_settings = RI_Util::update_coulomb_settings(this->info.coulomb_param, ucell, this->p_kv);
 
-	bool init_MGT = true;
+	std::shared_ptr<ORB_gaunt_table> MGT = std::make_shared<ORB_gaunt_table>();
 	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->MGT, init_MGT, settings_list.second.first );
-		init_MGT = false; // only init once
+															this->info.kmesh_times, MGT, settings_list.second.first );
 	}
 
 	ModuleBase::timer::tick("Exx_LRI", "init");
@@ -73,11 +71,6 @@ void Exx_LRI<Tdata>::cal_exx_ions(const UnitCell& ucell,
 	ModuleBase::TITLE("Exx_LRI","cal_exx_ions");
 	ModuleBase::timer::tick("Exx_LRI", "cal_exx_ions");
 
-	// init_radial_table_ions( cal_atom_centres_core(atom_pairs_core_origin), atom_pairs_core_origin );
-
-	// this->m_abfsabfs.init_radial_table(Rradial);
-	// this->m_abfslcaos_lcaos.init_radial_table(Rradial);
-
 	std::vector<TA> atoms(ucell.nat);
 	for(int iat=0; iat<ucell.nat; ++iat)
 		{ atoms[iat] = iat; }

source/source_lcao/module_ri/Exx_LRI_interface.hpp

@@ -111,7 +111,7 @@ void Exx_LRI_Interface<T, Tdata>::exx_before_all_runners(
         this->symrot_.find_irreducible_sector(
             ucell.symm, ucell.atoms, ucell.st,
             RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat);
-        this->symrot_.set_abfs_Lmax(GlobalC::exx_info.info_ri.abfs_Lmax);
+        this->symrot_.set_abfs_Lmax(Exx_Abfs::get_Lmax(this->exx_ptr->abfs));
         this->symrot_.cal_Ms(kv, ucell, pv);
     }
 }

source/source_lcao/module_ri/LRI_CV.h

@@ -40,8 +40,7 @@ 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,
-		ORB_gaunt_table& MGT,
-        const bool& init_MGT,
+		std::shared_ptr<ORB_gaunt_table> MGT,
         const bool& init_C);
 	inline std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
 	cal_Vs(

source/source_lcao/module_ri/LRI_CV.hpp

@@ -44,8 +44,7 @@ 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,
-	ORB_gaunt_table& MGT,
-    const bool& init_MGT,
+	std::shared_ptr<ORB_gaunt_table> MGT,
     const bool& init_C)
 {
 	ModuleBase::TITLE("LRI_CV", "set_orbitals");
@@ -69,24 +68,17 @@ void LRI_CV<Tdata>::set_orbitals(
 		range_abfs = ModuleBase::Element_Basis_Index::construct_range( abfs );
 	this->index_abfs = ModuleBase::Element_Basis_Index::construct_index( range_abfs );
 
-    const int Lmax_v = this->m_abfs_abfs.init(2, ucell, orb, kmesh_times, lcaos_rmax + abfs_ccp_rmax);
-    int Lmax_c = std::numeric_limits<double>::min();
+	this->m_abfs_abfs.MGT = this->m_abfslcaos_lcaos.MGT = MGT;
+    this->m_abfs_abfs.init(
+		this->abfs_ccp, this->abfs,
+		ucell, orb, kmesh_times, lcaos_rmax + abfs_ccp_rmax);
     if (init_C)
-        Lmax_c = this->m_abfslcaos_lcaos.init(1, ucell, orb, kmesh_times, lcaos_rmax);
-    const int Lmax = std::max(Lmax_v, Lmax_c);
+        this->m_abfslcaos_lcaos.init(
+			this->abfs_ccp, this->lcaos, this->lcaos,
+			ucell, orb, kmesh_times, lcaos_rmax);
 
-    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();
     }