Refactor: namespace Conv_Coulomb_Pot_K

source/module_ri/Exx_LRI.hpp

@@ -12,7 +12,6 @@
 #include "module_ri/exx_abfs-construct_orbs.h"
 #include "module_ri/exx_abfs-io.h"
 #include "module_ri/conv_coulomb_pot_k.h"
-#include "module_ri/conv_coulomb_pot_k-template.h"
 #include "module_base/tool_title.h"
 #include "module_base/timer.h"
 #include "module_ri/serialization_cereal.h"
@@ -71,14 +70,19 @@ void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in)
 			case Conv_Coulomb_Pot_K::Ccp_Type::Ccp:
 				return {};
 			case Conv_Coulomb_Pot_K::Ccp_Type::Hf:
-				return {};
+			{
+				// 4/3 * pi * Rcut^3 = V_{supercell} = V_{unitcell} * Nk
+				const int nspin0 = (GlobalV::NSPIN==2) ? 2 : 1;
+				const double hf_Rcut = std::pow(0.75 * this->p_kv->nkstot_full/nspin0 * GlobalC::ucell.omega / (ModuleBase::PI), 1.0/3.0);
+				return {{"hf_Rcut", hf_Rcut}};
+			}
 			case Conv_Coulomb_Pot_K::Ccp_Type::Hse:
 				return {{"hse_omega", this->info.hse_omega}};
 			default:
 				throw std::domain_error(std::string(__FILE__)+" line "+std::to_string(__LINE__));	break;
 		}
 	};
-    this->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, this->info.ccp_type, get_ccp_parameter(), this->info.ccp_rmesh_times, this->p_kv->nkstot_full);
+    this->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, this->info.ccp_type, get_ccp_parameter(), this->info.ccp_rmesh_times);
 
 
 	for( size_t T=0; T!=this->abfs.size(); ++T )

source/module_ri/conv_coulomb_pot_k.cpp

@@ -2,104 +2,109 @@
 #include "../module_base/constants.h"
 #include "../module_basis/module_ao/ORB_atomic_lm.h"
 #include "../module_hamilt_pw/hamilt_pwdft/global.h"
-std::vector<double> Conv_Coulomb_Pot_K::cal_psi_ccp( const std::vector<double> & psif )
+
+namespace Conv_Coulomb_Pot_K
 {
-	std::vector<double> psik2_ccp(psif.size());
-	for( size_t ik=0; ik<psif.size(); ++ik )
-		psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik];
-	return psik2_ccp;
-}
 
-// rongshi add 2022-07-27
-// Sphere truction -- Spencer
-std::vector<double> Conv_Coulomb_Pot_K::cal_psi_hf(const int& nks, const std::vector<double> &psif,
-                                                   const std::vector<double> &k_radial,
-                                                   const double omega = 0)
-{	
-    const int nspin0 = (GlobalV::NSPIN==2) ? 2 : 1;
-    const double Rc = std::pow(0.75 * nks/nspin0 * GlobalC::ucell.omega / (ModuleBase::PI), 1.0/3.0);
-    std::vector<double> psik2_ccp(psif.size());
-    for (size_t ik = 0; ik < psif.size(); ++ik)
-        psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik] * (1 - std::cos(k_radial[ik] * Rc));
-    return psik2_ccp;
-}
+	std::vector<double> cal_psi_ccp(
+		const std::vector<double> & psif)
+	{
+		std::vector<double> psik2_ccp(psif.size());
+		for( size_t ik=0; ik<psif.size(); ++ik )
+			psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik];
+		return psik2_ccp;
+	}
 
+	// rongshi add 2022-07-27
+	// Sphere truction -- Spencer
+	std::vector<double> cal_psi_hf(
+		const std::vector<double> &psif,
+		const std::vector<double> &k_radial,
+		const double hf_Rcut)
+	{
+		std::vector<double> psik2_ccp(psif.size());
+		for (size_t ik = 0; ik < psif.size(); ++ik)
+			psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik] * (1 - std::cos(k_radial[ik] * hf_Rcut));
+		return psik2_ccp;
+	}
 
-std::vector<double> Conv_Coulomb_Pot_K::cal_psi_hse( 
-	const std::vector<double> & psif,
-	const std::vector<double> & k_radial,
-	const double omega)
-{
-	std::vector<double> psik2_ccp(psif.size());
-	for( size_t ik=0; ik<psif.size(); ++ik )
-		psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik] * (1-std::exp(-(k_radial[ik]*k_radial[ik])/(4*omega*omega)));
-	return psik2_ccp;
-}
 
+	std::vector<double> cal_psi_hse(
+		const std::vector<double> & psif,
+		const std::vector<double> & k_radial,
+		const double hse_omega)
+	{
+		std::vector<double> psik2_ccp(psif.size());
+		for( size_t ik=0; ik<psif.size(); ++ik )
+			psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik] * (1-std::exp(-(k_radial[ik]*k_radial[ik])/(4*hse_omega*hse_omega)));
+		return psik2_ccp;
+	}
 
 
-template<>
-Numerical_Orbital_Lm Conv_Coulomb_Pot_K::cal_orbs_ccp<Numerical_Orbital_Lm>(
-	const Numerical_Orbital_Lm &orbs,
-	const Ccp_Type &ccp_type,
-	const std::map<std::string,double> &parameter,
-    const double rmesh_times,
-    const int& nks)
-{
-	std::vector<double> psik2_ccp;
-	switch(ccp_type)
+
+	template<>
+	Numerical_Orbital_Lm cal_orbs_ccp<Numerical_Orbital_Lm>(
+		const Numerical_Orbital_Lm &orbs,
+		const Ccp_Type &ccp_type,
+		const std::map<std::string,double> &parameter,
+		const double rmesh_times)
 	{
-		case Ccp_Type::Ccp:
-			psik2_ccp = cal_psi_ccp( orbs.get_psif() );		break;
-		case Ccp_Type::Hf:
-        	psik2_ccp = cal_psi_hf(nks, orbs.get_psif(), orbs.get_k_radial());      break;
-		case Ccp_Type::Hse:
-			psik2_ccp = cal_psi_hse( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hse_omega") );		break;
-		default:
-			throw( ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__) );		break;
-	}
+		std::vector<double> psik2_ccp;
+		switch(ccp_type)
+		{
+			case Ccp_Type::Ccp:
+				psik2_ccp = cal_psi_ccp( orbs.get_psif() );		break;
+			case Ccp_Type::Hf:
+				psik2_ccp = cal_psi_hf( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hf_Rcut"));      break;
+			case Ccp_Type::Hse:
+				psik2_ccp = cal_psi_hse( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hse_omega") );		break;
+			default:
+				throw( ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__) );		break;
+		}
 
-	const double dr = orbs.get_rab().back();
-	const int Nr = (static_cast<int>(orbs.getNr()*rmesh_times)) | 1;
-	std::vector<double> rab(Nr);
-	for( size_t ir=0; ir<std::min(orbs.getNr(),Nr); ++ir )
-		rab[ir] = orbs.getRab(ir);
-	for( size_t ir=orbs.getNr(); ir<Nr; ++ir )
-		rab[ir] = dr;
-	std::vector<double> r_radial(Nr);
-	for( size_t ir=0; ir<std::min(orbs.getNr(),Nr); ++ir )
-		r_radial[ir] = orbs.getRadial(ir);
-	for( size_t ir=orbs.getNr(); ir<Nr; ++ir )
-        r_radial[ir] = orbs.get_r_radial().back() + (ir - orbs.getNr() + 1) * dr;
-
-	Numerical_Orbital_Lm orbs_ccp;
-	orbs_ccp.set_orbital_info(
- 		orbs.getLabel(),
-	 	orbs.getType(),
-		orbs.getL(),
-		orbs.getChi(),
-	    Nr,
-		ModuleBase::GlobalFunc::VECTOR_TO_PTR(rab),
-		ModuleBase::GlobalFunc::VECTOR_TO_PTR(r_radial),
-		Numerical_Orbital_Lm::Psi_Type::Psik2,
-		ModuleBase::GlobalFunc::VECTOR_TO_PTR(psik2_ccp),
-		orbs.getNk(),
-		orbs.getDk(),
-		orbs.getDruniform(),
-		false,
-		true, GlobalV::CAL_FORCE);
-	return orbs_ccp;
-}
+		const double dr = orbs.get_rab().back();
+		const int Nr = (static_cast<int>(orbs.getNr()*rmesh_times)) | 1;
+		std::vector<double> rab(Nr);
+		for( size_t ir=0; ir<std::min(orbs.getNr(),Nr); ++ir )
+			rab[ir] = orbs.getRab(ir);
+		for( size_t ir=orbs.getNr(); ir<Nr; ++ir )
+			rab[ir] = dr;
+		std::vector<double> r_radial(Nr);
+		for( size_t ir=0; ir<std::min(orbs.getNr(),Nr); ++ir )
+			r_radial[ir] = orbs.getRadial(ir);
+		for( size_t ir=orbs.getNr(); ir<Nr; ++ir )
+			r_radial[ir] = orbs.get_r_radial().back() + (ir - orbs.getNr() + 1) * dr;
 
-template<>
-double Conv_Coulomb_Pot_K::get_rmesh_proportion(
-	const Numerical_Orbital_Lm &orbs,
-	const double psi_threshold)
-{
-	for(int ir=orbs.getNr()-1; ir>=0; --ir)
+		Numerical_Orbital_Lm orbs_ccp;
+		orbs_ccp.set_orbital_info(
+			orbs.getLabel(),
+			orbs.getType(),
+			orbs.getL(),
+			orbs.getChi(),
+			Nr,
+			ModuleBase::GlobalFunc::VECTOR_TO_PTR(rab),
+			ModuleBase::GlobalFunc::VECTOR_TO_PTR(r_radial),
+			Numerical_Orbital_Lm::Psi_Type::Psik2,
+			ModuleBase::GlobalFunc::VECTOR_TO_PTR(psik2_ccp),
+			orbs.getNk(),
+			orbs.getDk(),
+			orbs.getDruniform(),
+			false,
+			true, GlobalV::CAL_FORCE);
+		return orbs_ccp;
+	}
+
+	template<>
+	double get_rmesh_proportion(
+		const Numerical_Orbital_Lm &orbs,
+		const double psi_threshold)
 	{
-		if(std::abs(orbs.getPsi(ir))>=psi_threshold)
-			return static_cast<double>(ir)/orbs.getNr();
+		for(int ir=orbs.getNr()-1; ir>=0; --ir)
+		{
+			if(std::abs(orbs.getPsi(ir))>=psi_threshold)
+				return static_cast<double>(ir)/orbs.getNr();
+		}
+		return 0.0;
 	}
-	return 0.0;
+
 }

source/module_ri/conv_coulomb_pot_k.h

@@ -5,40 +5,37 @@
 #include <map>
 #include <string>
 
-class Conv_Coulomb_Pot_K
+namespace Conv_Coulomb_Pot_K
 {
-public:
+	enum class Ccp_Type{		//	parameter:
+		Ccp,					//
+		Hf,						//		"hf_Rcut"
+		Hse};					//		"hse_omega"
 
-	enum class Ccp_Type{		//  parameter:
-		Ccp,                 // 
-		Hf,					//
-		Hse};					//  	"hse_omega"
-
-	template<typename T> static T cal_orbs_ccp(
+	template<typename T> T cal_orbs_ccp(
 		const T &orbs,
 		const Ccp_Type &ccp_type,
 		const std::map<std::string,double> &parameter,
-        const double rmesh_times,
-        const int& nks);
-
-private:
-
-	template< typename T > static double get_rmesh_proportion(
+        const double rmesh_times);
+
+  //private:
+	template< typename T > double get_rmesh_proportion(
 		const T &orbs,
 		const double psi_threshold);
-
-private:
 
-	static std::vector<double> cal_psi_ccp( const std::vector<double> & psif );
-
-	static std::vector<double> cal_psi_hf(const int& nks, const std::vector<double> &psif,
-                                          const std::vector<double> &k_radial,
-                                          const double omega);
-
-	static std::vector<double> cal_psi_hse( 
+  //private:
+	std::vector<double> cal_psi_ccp(
+		const std::vector<double> & psif);
+	std::vector<double> cal_psi_hf(
+		const std::vector<double> &psif,
+		const std::vector<double> &k_radial,
+		const double hf_Rcut);
+	std::vector<double> cal_psi_hse(
 		const std::vector<double> & psif,
 		const std::vector<double> & k_radial,
-		const double omega);
-};
+		const double hse_omega);
+}
+
+#include "conv_coulomb_pot_k.hpp"
 
 #endif

source/module_ri/conv_coulomb_pot_k.hpp

@@ -0,0 +1,37 @@
+#ifndef CONV_COULOMB_POT_K_HPP
+#define CONV_COULOMB_POT_K_HPP
+
+#include "conv_coulomb_pot_k.h"
+#include <vector>
+#include <cmath>
+
+namespace Conv_Coulomb_Pot_K
+{
+
+	template< typename T >
+	std::vector<T> cal_orbs_ccp(
+		const std::vector<T> & orbs,
+		const Ccp_Type &ccp_type,
+		const std::map<std::string,double> &parameter,
+		const double rmesh_times)
+	{
+		std::vector<T> orbs_ccp(orbs.size());
+		for( size_t i=0; i!=orbs.size(); ++i )
+			orbs_ccp[i] = cal_orbs_ccp(orbs[i], ccp_type, parameter, rmesh_times);
+		return orbs_ccp;
+	}
+
+	template< typename T >
+	double get_rmesh_proportion(
+		const std::vector<T> & orbs,
+		const double psi_threshold)
+	{
+		double rmesh_proportion=0;
+		for( const auto &orb : orbs )
+			rmesh_proportion = std::max(rmesh_proportion, get_rmesh_proportion(orb,psi_threshold));
+		return rmesh_proportion;
+	}
+
+}
+
+#endif