Refactor: the ccp parameters to support additional Coulomb matrix construction methods.

source/module_hamilt_general/module_xc/exx_info.h

@@ -15,7 +15,10 @@ struct Exx_Info
     {
         bool cal_exx = false;
 
-        std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> coulomb_param;
+        std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Method, 
+            std::pair<bool, 
+                std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, 
+                    std::vector<std::map<std::string,std::string>>>>> coulomb_settings;
 		// Fock:
 		//		"alpha":		"0"
 		//		"Rcut_type":	"limits" / "spencer"
@@ -52,7 +55,10 @@ struct Exx_Info
 
     struct Exx_Info_RI
     {
-        const std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> &coulomb_param;
+        const std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Method, 
+            std::pair<bool, 
+                std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, 
+                    std::vector<std::map<std::string,std::string>>>>> &coulomb_settings;
 
         bool real_number = false;
 
@@ -74,7 +80,7 @@ struct Exx_Info
         int abfs_Lmax = 0; // tmp
 
         Exx_Info_RI(const Exx_Info::Exx_Info_Global& info_global)
-            : coulomb_param(info_global.coulomb_param)
+            : coulomb_settings(info_global.coulomb_settings)
         {
         }
     };

source/module_io/input_conv.cpp

@@ -328,35 +328,41 @@ void Input_Conv::Convert()
     {
         GlobalC::exx_info.info_global.cal_exx = true;
         GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf;
-        GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock] = {{
+        std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> coulomb_param;
+        coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock] = {{
             {"alpha", "1"},
             {"Rcut_type", "spencer"},
             {"lambda", std::to_string(PARAM.inp.exx_lambda)} }};
+        GlobalC::exx_info.info_global.coulomb_settings[Conv_Coulomb_Pot_K::Coulomb_Method::Center2] = std::make_pair(true, coulomb_param);
     }
 	// use the error function erf(w|r-r'|), exx just has the short-range part
     else if (dft_functional_lower == "hse"
           || dft_functional_lower == "cwp22")
     {
         GlobalC::exx_info.info_global.cal_exx = true;
         GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Erfc;
-        GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
+        std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> coulomb_param;
+        coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
             {"alpha", "1"},
 			{"omega", std::to_string(PARAM.inp.exx_hse_omega)},
             {"Rcut_type", "limits"} }};
+        GlobalC::exx_info.info_global.coulomb_settings[Conv_Coulomb_Pot_K::Coulomb_Method::Center2] = std::make_pair(true, coulomb_param);
     }
 	// use the error function erf(w|r-r'|), exx just has the long-range part
     else if ( dft_functional_lower == "wp22" )
     {
         GlobalC::exx_info.info_global.cal_exx = true;
         GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Erf;
-        GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock] = {{
+        std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> coulomb_param;
+        coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock] = {{
             {"alpha", "1"},
             {"Rcut_type", "spencer"},
             {"lambda", std::to_string(PARAM.inp.exx_lambda)} }};
-        GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
+        coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
             {"alpha", "-1"},
 			{"omega", std::to_string(PARAM.inp.exx_hse_omega)},
             {"Rcut_type", "limits"} }};
+        GlobalC::exx_info.info_global.coulomb_settings[Conv_Coulomb_Pot_K::Coulomb_Method::Center2] = std::make_pair(true, coulomb_param);
     }
 #ifdef __EXX
     else if (dft_functional_lower == "opt_orb")

source/module_ri/Exx_LRI.h

@@ -15,6 +15,7 @@
 #include <vector>
 #include <array>
 #include <map>
+#include <unordered_map>
 #include <deque>
 #include <mpi.h>
 
@@ -37,6 +38,15 @@
 		class OperatorLREXX;
 	}
 
+template<typename Tdata>
+class Exx_Obj
+{
+	// match with Conv_Coulomb_Pot_K::Coulomb_Method
+	public:
+		LRI_CV<Tdata> cv;
+		std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs_ccp;
+};
+
 template<typename Tdata>
 class Exx_LRI
 {
@@ -85,9 +95,9 @@ class Exx_LRI
 
 	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;
-
-	LRI_CV<Tdata> cv;
+	//std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs_ccp;
+	std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Method, Exx_Obj<Tdata>> exx_objs;
+	//LRI_CV<Tdata> cv;
 	RI::Exx<TA,Tcell,Ndim,Tdata> exx_lri;
 
 	void post_process_Hexx( std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> &Hexxs_io ) const;

source/module_ri/Exx_LRI.hpp

@@ -48,19 +48,18 @@ void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in,
 		{ this->abfs = Exx_Abfs::IO::construct_abfs( abfs_same_atom, orb, this->info.files_abfs, this->info.kmesh_times ); 	}
 	Exx_Abfs::Construct_Orbs::print_orbs_size(ucell, this->abfs, GlobalV::ofs_running);
 
-	const std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>>
-		coulomb_param_updated = RI_Util::update_coulomb_param(this->info.coulomb_param, ucell.omega, this->p_kv->get_nkstot_full());
-	this->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, coulomb_param_updated, this->info.ccp_rmesh_times);
-
 	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->cv.set_orbitals(
-		ucell,
-		orb,
-		this->lcaos, this->abfs, this->abfs_ccp,
-		this->info.kmesh_times, this->info.ccp_rmesh_times );
-
+	for(const auto &settings_list : this->info.coulomb_settings)
+	{
+		const std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>>
+			coulomb_param_updated = RI_Util::update_coulomb_param(settings_list.second.second, ucell.omega, this->p_kv->get_nkstot_full());
+		this->exx_objs[settings_list.first].abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, coulomb_param_updated, 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 );
+	}
 	ModuleBase::timer::tick("Exx_LRI", "init");
 }
 
@@ -95,26 +94,40 @@ void Exx_LRI<Tdata>::cal_exx_ions(const UnitCell& ucell,
 	const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA,std::array<Tcell,Ndim>>>>>
 		list_As_Vs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false);
 
-	std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
-		Vs = this->cv.cal_Vs(ucell,
-			list_As_Vs.first, list_As_Vs.second[0],
-			{{"writable_Vws",true}});
-	this->cv.Vws = LRI_CV_Tools::get_CVws(ucell,Vs);
+	std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> Vs;
+	std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, Ndim>>> dVs;
+	for(const auto &settings_list : this->info.coulomb_settings)
+	{
+		std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
+			Vs_temp = this->exx_objs[settings_list.first].cv.cal_Vs(ucell,
+				list_As_Vs.first, list_As_Vs.second[0],
+				{{"writable_Vws",true}});
+		this->exx_objs[settings_list.first].cv.Vws = LRI_CV_Tools::get_CVws(ucell,Vs_temp);
+		Vs = Vs.empty() ? Vs_temp : LRI_CV_Tools::add(Vs, Vs_temp);
+
+		if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
+		{
+			std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, Ndim>>>
+				dVs_temp = this->exx_objs[settings_list.first].cv.cal_dVs(ucell,
+					list_As_Vs.first, list_As_Vs.second[0],
+					{{"writable_dVws",true}});
+			this->exx_objs[settings_list.first].cv.dVws = LRI_CV_Tools::get_dCVws(ucell,dVs_temp);
+			dVs = dVs.empty() ? dVs_temp : LRI_CV_Tools::add(dVs, dVs_temp);
+		}
+
+	}
 	if (write_cv && GlobalV::MY_RANK == 0)
 		{ LRI_CV_Tools::write_Vs_abf(Vs, PARAM.globalv.global_out_dir + "Vs"); }
 	this->exx_lri.set_Vs(std::move(Vs), this->info.V_threshold);
 
 	if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
 	{
-		std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>
-			dVs = this->cv.cal_dVs(ucell,
-				list_As_Vs.first, list_As_Vs.second[0],
-				{{"writable_dVws",true}});
-		this->cv.dVws = LRI_CV_Tools::get_dCVws(ucell,dVs);
-		this->exx_lri.set_dVs(std::move(dVs), this->info.V_grad_threshold);
+		std::array<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>, Ndim>
+			dVs_order = LRI_CV_Tools::change_order(std::move(dVs));
+		this->exx_lri.set_dVs(std::move(dVs_order), this->info.V_grad_threshold);
 		if(PARAM.inp.cal_stress)
 		{
-			std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dVRs = LRI_CV_Tools::cal_dMRs(ucell,dVs);
+			std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dVRs = LRI_CV_Tools::cal_dMRs(ucell,dVs_order);
 			this->exx_lri.set_dVRs(std::move(dVRs), this->info.V_grad_R_threshold);
 		}
 	}
@@ -123,29 +136,47 @@ void Exx_LRI<Tdata>::cal_exx_ions(const UnitCell& ucell,
 	const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA,std::array<Tcell,Ndim>>>>>
 		list_As_Cs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false);
 
-	std::pair<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>, std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>>
-		Cs_dCs = this->cv.cal_Cs_dCs(
-			ucell,
-			list_As_Cs.first, list_As_Cs.second[0],
-			{{"cal_dC",PARAM.inp.cal_force||PARAM.inp.cal_stress},
-			 {"writable_Cws",true}, {"writable_dCws",true}, {"writable_Vws",false}, {"writable_dVws",false}});
-	std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> &Cs = std::get<0>(Cs_dCs);
-	this->cv.Cws = LRI_CV_Tools::get_CVws(ucell,Cs);
+	std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> Cs;
+	std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>> dCs;
+	for(const auto &settings_list : this->info.coulomb_settings)
+	{
+		if(settings_list.second.first)
+		{
+			std::pair<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>,
+				std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>>
+					Cs_dCs = this->exx_objs[settings_list.first].cv.cal_Cs_dCs(
+						ucell,
+						list_As_Cs.first, list_As_Cs.second[0],
+						{{"cal_dC",PARAM.inp.cal_force||PARAM.inp.cal_stress},
+						{"writable_Cws",true}, {"writable_dCws",true}, {"writable_Vws",false}, {"writable_dVws",false}});
+			std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> &Cs_temp = std::get<0>(Cs_dCs);
+			this->exx_objs[settings_list.first].cv.Cws = LRI_CV_Tools::get_CVws(ucell,Cs_temp);
+			Cs = Cs.empty() ? Cs_temp : LRI_CV_Tools::add(Cs, Cs_temp);
+
+			if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
+			{
+				std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>> &dCs_temp = std::get<1>(Cs_dCs);
+				this->exx_objs[settings_list.first].cv.dCws = LRI_CV_Tools::get_dCVws(ucell,dCs_temp);
+				dCs = dCs.empty() ? dCs_temp : LRI_CV_Tools::add(dCs, dCs_temp);
+			}
+		}
+	}
 	if (write_cv && GlobalV::MY_RANK == 0)
 		{ LRI_CV_Tools::write_Cs_ao(Cs, PARAM.globalv.global_out_dir + "Cs"); }
 	this->exx_lri.set_Cs(std::move(Cs), this->info.C_threshold);
 
 	if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
 	{
-		std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3> &dCs = std::get<1>(Cs_dCs);
-		this->cv.dCws = LRI_CV_Tools::get_dCVws(ucell,dCs);
-		this->exx_lri.set_dCs(std::move(dCs), this->info.C_grad_threshold);
+		std::array<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>, Ndim>
+			dCs_order = LRI_CV_Tools::change_order(std::move(dCs));
+		this->exx_lri.set_dCs(std::move(dCs_order), this->info.C_grad_threshold);
 		if(PARAM.inp.cal_stress)
 		{
-			std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dCRs = LRI_CV_Tools::cal_dMRs(ucell,dCs);
+			std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dCRs = LRI_CV_Tools::cal_dMRs(ucell,dCs_order);
 			this->exx_lri.set_dCRs(std::move(dCRs), this->info.C_grad_R_threshold);
 		}
 	}
+
 	ModuleBase::timer::tick("Exx_LRI", "cal_exx_ions");
 }
 

source/module_ri/LRI_CV.h

@@ -47,14 +47,14 @@ class LRI_CV
 		const std::vector<TA> &list_A0,
 		const std::vector<TAC> &list_A1,
 		const std::map<std::string,bool> &flags);						// "writable_Vws"
-	inline std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>
+	inline std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>
 	cal_dVs(
 		const UnitCell &ucell,
 		const std::vector<TA> &list_A0,
 		const std::vector<TAC> &list_A1,
 		const std::map<std::string,bool> &flags);						// "writable_dVws"
-	std::pair<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,
-	          std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>>
+    std::pair<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>,
+              std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>>
 	cal_Cs_dCs(
 		const UnitCell &ucell,
 		const std::vector<TA> &list_A0,

source/module_ri/LRI_CV.hpp

@@ -146,15 +146,14 @@ auto LRI_CV<Tdata>::cal_dVs(
 	const std::vector<TA> &list_A0,
 	const std::vector<TAC> &list_A1,
 	const std::map<std::string,bool> &flags)					// + "writable_dVws"
--> std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>
+-> std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>
 {
 	ModuleBase::TITLE("LRI_CV","cal_dVs");
 	const T_func_DPcal_data<std::array<RI::Tensor<Tdata>,3>>
 		func_DPcal_dV = std::bind(
 			&LRI_CV<Tdata>::DPcal_dV, this,
 			std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
-	return LRI_CV_Tools::change_order(
-		this->cal_datas(ucell,list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_dV));
+	return this->cal_datas(ucell,list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_dV);
 }
 
 template<typename Tdata>
@@ -163,7 +162,9 @@ auto LRI_CV<Tdata>::cal_Cs_dCs(
 	const std::vector<TA> &list_A0,
 	const std::vector<TAC> &list_A1,
 	const std::map<std::string,bool> &flags)					// "cal_dC" + "writable_Cws", "writable_dCws", "writable_Vws", "writable_dVws"
--> std::pair<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>, std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>>
+-> std::pair<
+        std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>,
+        std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>>
 {
 	ModuleBase::TITLE("LRI_CV","cal_Cs_dCs");
 	const T_func_DPcal_data<std::pair<RI::Tensor<Tdata>, std::array<RI::Tensor<Tdata>,3>>>
@@ -173,16 +174,16 @@ auto LRI_CV<Tdata>::cal_Cs_dCs(
 	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);
 
-	std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> Cs;
-	std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3> dCs;
-	for(auto &Cs_dCs_A : Cs_dCs_tmp)
-		for(auto &Cs_dCs_B : Cs_dCs_A.second)
-		{
-			Cs[Cs_dCs_A.first][Cs_dCs_B.first] = std::move(std::get<0>(Cs_dCs_B.second));
-			if(flags.at("cal_dC"))
-				for(int ix=0; ix<3; ++ix)
-					dCs[ix][Cs_dCs_A.first][Cs_dCs_B.first] = std::move(std::get<1>(Cs_dCs_B.second)[ix]);
-		}
+	std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Cs;
+    std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>> dCs;
+    for (auto& Cs_dCs_A: Cs_dCs_tmp)
+        for (auto& Cs_dCs_B: Cs_dCs_A.second) {
+            Cs[Cs_dCs_A.first][Cs_dCs_B.first]
+                = std::move(std::get<0>(Cs_dCs_B.second));
+            if (flags.at("cal_dC"))
+                dCs[Cs_dCs_A.first][Cs_dCs_B.first]
+                    = std::move(std::get<1>(Cs_dCs_B.second));
+        }
 	return std::make_pair(Cs, dCs);
 }