Merge pull request #1120 from dyzheng/develop

Merge: abacusmodeling/develop PR166, PR167, PR170

Author: caic99

source/module_orbital/ORB_atomic_lm.cpp

@@ -178,7 +178,7 @@ void Numerical_Orbital_Lm::copy_parameter(
     // assert(r_radial_in[nr-1]<50);						// Peize Lin delete 2017-08-18
     this->rcut = r_radial_in[nr-1];
     assert(nk_in>1);
-    assert(nk_in<10000);
+    //assert(nk_in<10000);				// Jiyy delete 2022-07-18
     this->nk = nk_in;
     assert(nk%2!=0);
     assert(dk_in>0);

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:

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,

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::SCAN0:
 			abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp( this->abfs, Conv_Coulomb_Pot_K::Ccp_Type::Ccp, {}, info.ccp_rmesh_times );		break;
 		case Exx_Global::Hybrid_Type::HSE:

tools/plot-tools/abacus_plot/dos.py

@@ -391,7 +391,7 @@ def _parial_plot(self,
         Returns:
             DOSPlot object: for manually plotting picture with dosplot.ax 
         """
-        if not isinstance(ax, list):
+        if isinstance(ax, axes.Axes):
             ax = [ax]
 
         dos, totnum = parse_projected_data(self.orbitals, species, keyname)