Merge pull request #170 from Srlive1201/develop

add Sphere truction of coulomb in HF/PBE0
abacusmodeling · Jul 27, 2022 · 44c02d5 · 44c02d5
2 parents cb47661 + 34ffbbe
commit 44c02d5
Show file tree

Hide file tree

Showing 3 changed files with 26 additions and 3 deletions.
diff --git a/source/src_ri/conv_coulomb_pot_k.cpp b/source/src_ri/conv_coulomb_pot_k.cpp
@@ -1,7 +1,7 @@
 #include "conv_coulomb_pot_k.h"
 #include "../module_base/constants.h"
 #include "../module_orbital/ORB_atomic_lm.h"
-
+#include "../src_pw/global.h"
 std::vector<double> Conv_Coulomb_Pot_K::cal_psi_ccp( const std::vector<double> & psif )
 {
 	std::vector<double> psik2_ccp(psif.size());
@@ -10,6 +10,21 @@ std::vector<double> Conv_Coulomb_Pot_K::cal_psi_ccp( const std::vector<double> &
 	return psik2_ccp;
 }
 
+// rongshi add 2022-07-27
+// Sphere truction -- Spencer
+std::vector<double> Conv_Coulomb_Pot_K::cal_psi_hf(const std::vector<double> &psif,
+                                                   const std::vector<double> &k_radial,
+                                                   const double omega = 0)
+{	
+    double Rc = pow(0.75 * GlobalC::kv.nks * GlobalC::ucell.omega / (ModuleBase::PI), 0.3333334);
+   // cout << "hf_Rc:  " << Rc << endl;
+    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> Conv_Coulomb_Pot_K::cal_psi_hse( 
 	const std::vector<double> & psif,
 	const std::vector<double> & k_radial,
@@ -35,6 +50,8 @@ Numerical_Orbital_Lm Conv_Coulomb_Pot_K::cal_orbs_ccp<Numerical_Orbital_Lm>(
 	{
 		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());      break;
 		case Ccp_Type::Hse:
 			psik2_ccp = cal_psi_hse( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hse_omega") );		break;
 		default:

diff --git a/source/src_ri/conv_coulomb_pot_k.h b/source/src_ri/conv_coulomb_pot_k.h
@@ -10,7 +10,8 @@ class Conv_Coulomb_Pot_K
 public:
 
 	enum class Ccp_Type{		//  parameter:
-		Ccp,					//
+		Ccp,                 // 
+		Hf,					//
 		Hse};					//  	"hse_omega"
 
 	template<typename T> static T cal_orbs_ccp(
@@ -35,6 +36,10 @@ class Conv_Coulomb_Pot_K
 private:
 
 	static std::vector<double> cal_psi_ccp( const std::vector<double> & psif );
+
+	static std::vector<double> cal_psi_hf(const std::vector<double> &psif,
+                                          const std::vector<double> &k_radial,
+                                          const double omega);
 
 	static std::vector<double> cal_psi_hse( 
 		const std::vector<double> & psif,

diff --git a/source/src_ri/exx_lcao.cpp b/source/src_ri/exx_lcao.cpp
@@ -643,8 +643,9 @@ gettimeofday( &t_start, NULL);
 	switch(info.hybrid_type)
 	{
 		case Exx_Global::Hybrid_Type::HF:
+			abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, Conv_Coulomb_Pot_K::Ccp_Type::Hf, {}, info.ccp_rmesh_times);   break;
 		case Exx_Global::Hybrid_Type::PBE0:
-			abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp( this->abfs, Conv_Coulomb_Pot_K::Ccp_Type::Ccp, {}, info.ccp_rmesh_times );		break;
+			abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp( this->abfs, Conv_Coulomb_Pot_K::Ccp_Type::Hf, {}, info.ccp_rmesh_times );		break;
 		case Exx_Global::Hybrid_Type::HSE:
 			abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp( this->abfs, Conv_Coulomb_Pot_K::Ccp_Type::Hse, {{"hse_omega",info.hse_omega}}, info.ccp_rmesh_times );	break;
 		default: