Merge pull request #1103 from ouqi0711/develop

Add missed Perdew-Wang correlation part for PBE0 functional

Author: dyzheng

source/input_conv.cpp

@@ -362,6 +362,7 @@ void Input_Conv::Convert(void)
     if (GlobalC::exx_global.info.hybrid_type != Exx_Global::Hybrid_Type::No)
     {
         GlobalC::exx_global.info.hybrid_alpha = INPUT.exx_hybrid_alpha;
+        XC_Functional::get_hybrid_alpha(INPUT.exx_hybrid_alpha);
         GlobalC::exx_global.info.hse_omega = INPUT.exx_hse_omega;
         GlobalC::exx_global.info.separate_loop = INPUT.exx_separate_loop;
         GlobalC::exx_global.info.hybrid_step = INPUT.exx_hybrid_step;

source/module_xc/xc_functional.cpp

@@ -9,6 +9,12 @@ XC_Functional::~XC_Functional(){}
 std::vector<int> XC_Functional::func_id(1);
 int XC_Functional::func_type = 0;
 bool XC_Functional::use_libxc = true;
+double XC_Functional::hybrid_alpha = 0.25;
+
+void XC_Functional::get_hybrid_alpha(const double alpha_in)
+{
+    hybrid_alpha = alpha_in;
+}
 
 int XC_Functional::get_func_type()
 {

source/module_xc/xc_functional.h

@@ -83,6 +83,7 @@ class XC_Functional
 
 	static int get_func_type();
 	static void set_xc_type(const std::string xc_func_in);
+	static void get_hybrid_alpha(const double alpha_in);
 #ifdef USE_LIBXC
 	static void set_xc_type_libxc(const std::string xc_func_in);
 	static std::vector<xc_func_type> init_func(const int xc_polarized);
@@ -94,6 +95,9 @@ class XC_Functional
 	static int func_type; //0:none, 1:lda, 2:gga, 3:mgga, 4:hybrid
 	static bool use_libxc;
 
+	//exx_hybrid_alpha for mixing exx in hybrid functional:
+	static double hybrid_alpha;
+
 //-------------------
 //  xc_functional_wrapper_xc.cpp
 //-------------------

source/module_xc/xc_functional_wrapper_gcxc.cpp

@@ -80,7 +80,9 @@ void XC_Functional::gcxc(const double &rho, const double &grho, double &sxc,
             case XC_HYB_GGA_XC_PBEH: //PBE0
                 double sx, v1x, v2x, sc, v1c, v2c;
                 XC_Functional::pbex(rho, grho, 0, sx, v1x, v2x);
-                sx *= 0.75; v1x *= 0.75; v2x *= 0.75;
+                sx *= (1.0 - XC_Functional::hybrid_alpha); 
+                v1x *= (1.0 - XC_Functional::hybrid_alpha); 
+                v2x *= (1.0 - XC_Functional::hybrid_alpha);
                 XC_Functional::pbec(rho, grho, 0, sc, v1c, v2c);
                 s = sx + sc;
                 v1 = v1x + v1c;
@@ -179,12 +181,16 @@ void XC_Functional::gcx_spin(double rhoup, double rhodw, double grhoup2, double
                 if (rhoup > small && sqrt(fabs(grhoup2)) > small)
                 {
                     XC_Functional::pbex(2.0 * rhoup, 4.0 * grhoup2, 0, sxup, v1xup, v2xup);
-                    sxup *= 0.75; v1xup *= 0.75; v2xup *= 0.75;
+                    sxup *= (1.0 - XC_Functional::hybrid_alpha); 
+                    v1xup *= (1.0 - XC_Functional::hybrid_alpha); 
+                    v2xup *= (1.0 - XC_Functional::hybrid_alpha);
                 }
                 if (rhodw > small && sqrt(fabs(grhodw2)) > small)
                 {
                     XC_Functional::pbex(2.0 * rhodw, 4.0 * grhodw2, 0, sxdw, v1xdw, v2xdw);
-        	    	sxdw *= 0.75; v1xdw *= 0.75; v2xdw *= 0.75;            
+        	    	sxdw *= (1.0 - XC_Functional::hybrid_alpha); 
+                    v1xdw *= (1.0 - XC_Functional::hybrid_alpha); 
+                    v2xdw *= (1.0 - XC_Functional::hybrid_alpha);            
                 }
                 break;
             case XC_GGA_X_PBE_SOL: //PBXsol

source/module_xc/xc_functional_wrapper_xc.cpp

@@ -33,8 +33,13 @@ void XC_Functional::xc(const double &rho, double &exc, double &vxc)
             // Exchange functionals containing attenuated slater exchange
             case XC_HYB_GGA_XC_PBEH:
             //  PBE0
-                XC_Functional::slater(rs, e, v);
-                e *= 0.75; v*= 0.75;
+                double ex, vx, ec, vc;
+                XC_Functional::slater(rs, ex, vx);
+                ex *= (1 - XC_Functional::hybrid_alpha); 
+                vx *= (1 - XC_Functional::hybrid_alpha);
+                XC_Functional::pw(rs, 0, ec, vc);
+                e = ex + ec;
+                v = vx + vc;
                 break;
 
             // Correlation functionals containing PW correlation
@@ -85,8 +90,15 @@ void XC_Functional::xc_spin(const double &rho, const double &zeta,
             // Exchange functionals containing attenuated slater exchange
             case XC_HYB_GGA_XC_PBEH:
             //  PBE0
-                XC_Functional::slater_spin(rho, zeta, e, vup, vdw);
-                e *= 0.75; vup *= 0.75; vdw *=0.75;
+                double ex, vupx, vdwx, ec, vupc, vdwc;
+                XC_Functional::slater_spin(rho, zeta, ex, vupx, vdwx);
+                ex *= (1.0 - XC_Functional::hybrid_alpha); 
+                vupx *= (1.0 - XC_Functional::hybrid_alpha); 
+                vdwx *= (1.0 - XC_Functional::hybrid_alpha);
+                XC_Functional::pw_spin(rs, zeta, ec, vupc, vdwc);
+                e = ex + ec;
+                vup = vupx + vupc;
+                vdw = vdwx + vdwc;
                 break;
 
             // Correlation functionals containing PZ correlation