Fix: nspin2/4 mismatch with nspin1 with PBE

Author: dyzheng

source/module_hamilt_general/module_xc/xc_funct_corr_gga.cpp

@@ -382,7 +382,7 @@ void XC_Functional::pbec_spin(double rho, double zeta, double grho, const int &i
     // real(kind=DP) :: rho, zeta, grho, sc, v1cup, v1cdw, v2c
     double ga, be[3];//mohan add
     // parameter :
-    ga = 0.0310910;
+    ga = 0.0310906908696548950;
     be[1] = 0.06672455060314922;//zhengdy add 2019-09-12, ensure same parameter with another dft code.
 	be[2] = 0.0460000;//mohan add 2012-05-28
     double third, pi34, xkf, xks;

source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp

@@ -127,6 +127,7 @@ std::vector<double> XC_Functional_Libxc::cal_sgn(
 	const std::vector<double> &sigma)
 {
 	std::vector<double> sgn(nrxx*nspin, 1.0);
+	const double grho_threshold2 = grho_threshold*grho_threshold;
 	// in the case of GGA correlation for polarized case,
 	// a cutoff for grho is required to ensure that libxc gives reasonable results
 	if(nspin==2 && func.info->family != XC_FAMILY_LDA && func.info->kind==XC_CORRELATION)
@@ -136,9 +137,9 @@ std::vector<double> XC_Functional_Libxc::cal_sgn(
 		#endif
 		for( int ir=0; ir<nrxx; ++ir )
 		{
-			if ( rho[ir*2]<rho_threshold || std::sqrt(std::abs(sigma[ir*3]))<grho_threshold )
+			if ( rho[ir*2]<rho_threshold || std::abs(sigma[ir*3])<grho_threshold2 )
 				sgn[ir*2] = 0.0;
-			if ( rho[ir*2+1]<rho_threshold || std::sqrt(std::abs(sigma[ir*3+2]))<grho_threshold )
+			if ( rho[ir*2+1]<rho_threshold || std::abs(sigma[ir*3+2])<grho_threshold2 )
 				sgn[ir*2+1] = 0.0;
 		}
 	}