diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/source/module_hamilt_general/module_xc/xc_functional.cpp
index 0168ecce99..c9730e957b 100644
--- a/source/module_hamilt_general/module_xc/xc_functional.cpp
+++ b/source/module_hamilt_general/module_xc/xc_functional.cpp
@@ -18,11 +18,16 @@ 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)
+void XC_Functional::set_hybrid_alpha(const double alpha_in)
 {
     hybrid_alpha = alpha_in;
 }
 
+double XC_Functional::get_hybrid_alpha()
+{
+    return hybrid_alpha;
+}
+
 int XC_Functional::get_func_type()
 {
     return func_type;
diff --git a/source/module_hamilt_general/module_xc/xc_functional.h b/source/module_hamilt_general/module_xc/xc_functional.h
index 5179c6bf93..ac9c75256a 100644
--- a/source/module_hamilt_general/module_xc/xc_functional.h
+++ b/source/module_hamilt_general/module_xc/xc_functional.h
@@ -66,7 +66,8 @@ class XC_Functional
     static void set_xc_type(const std::string xc_func_in);
 
     // For hybrid functional
-    static void get_hybrid_alpha(const double alpha_in);
+    static void set_hybrid_alpha(const double alpha_in);
+    static double get_hybrid_alpha();
     /// Usually in exx caculation, the first SCF loop should be converged with PBE
     static void set_xc_first_loop(const UnitCell& ucell);
 
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp
index 416b08b65e..8d80f27714 100644
--- a/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp
+++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp
@@ -4,6 +4,10 @@
 #include "module_parameter/parameter.h"
 #include "module_base/tool_quit.h"
 
+#ifdef __EXX
+#include "module_hamilt_pw/hamilt_pwdft/global.h"		// just for GlobalC::exx_info
+#endif
+
 #include <xc.h>
 
 #include <vector>
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp
new file mode 100644
index 0000000000..2ff901e786
--- /dev/null
+++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp
@@ -0,0 +1,292 @@
+#ifdef USE_LIBXC
+
+#include "xc_functional_libxc.h"
+#include "xc_functional.h"
+#include "module_elecstate/module_charge/charge.h"
+
+// converting rho (abacus=>libxc)
+std::vector<double> XC_Functional_Libxc::convert_rho(
+	const int nspin,
+	const std::size_t nrxx,
+	const Charge* const chr)
+{
+	std::vector<double> rho(nrxx*nspin);
+	#ifdef _OPENMP
+	#pragma omp parallel for collapse(2) schedule(static, 1024)
+	#endif
+	for( int is=0; is<nspin; ++is )
+		for( int ir=0; ir<nrxx; ++ir )
+			rho[ir*nspin+is] = chr->rho[is][ir] + 1.0/nspin*chr->rho_core[ir];
+	return rho;
+}
+
+// converting rho (abacus=>libxc)
+std::tuple<std::vector<double>, std::vector<double>>
+XC_Functional_Libxc::convert_rho_amag_nspin4(
+	const int nspin,
+	const std::size_t nrxx,
+	const Charge* const chr)
+{
+	assert(GlobalV::NSPIN==4);
+	std::vector<double> rho(nrxx*nspin);
+	std::vector<double> amag(nrxx);
+	#ifdef _OPENMP
+	#pragma omp parallel for
+	#endif
+	for( int ir=0; ir<nrxx; ++ir )
+	{
+		const double arhox = std::abs( chr->rho[0][ir] + chr->rho_core[ir] );
+		amag[ir] = std::sqrt( std::pow(chr->rho[1][ir],2)
+							+ std::pow(chr->rho[2][ir],2)
+							+ std::pow(chr->rho[3][ir],2) );
+		const double amag_clip = (amag[ir]<arhox) ? amag[ir] : arhox;
+		rho[ir*nspin+0] = (arhox + amag_clip) / 2.0;
+		rho[ir*nspin+1] = (arhox - amag_clip) / 2.0;
+	}
+	return std::make_tuple(std::move(rho), std::move(amag));
+}
+
+// calculating grho
+std::vector<std::vector<ModuleBase::Vector3<double>>>
+XC_Functional_Libxc::cal_gdr(
+	const int nspin,
+	const std::vector<double> &rho,
+	const double tpiba,
+	const Charge* const chr)
+{
+	std::vector<std::vector<ModuleBase::Vector3<double>>> gdr(nspin);
+	const std::size_t nrxx = rho.size();
+	for( int is=0; is!=nspin; ++is )
+	{
+		std::vector<double> rhor(nrxx);
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 1024)
+		#endif
+		for(std::size_t ir=0; ir<nrxx; ++ir)
+			rhor[ir] = rho[ir*nspin+is];
+		//------------------------------------------
+		// initialize the charge density array in reciprocal space
+		// bring electron charge density from real space to reciprocal space
+		//------------------------------------------
+		std::vector<std::complex<double>> rhog(chr->rhopw->npw);
+		chr->rhopw->real2recip(rhor.data(), rhog.data());
+
+		//-------------------------------------------
+		// compute the gradient of charge density and
+		// store the gradient in gdr[is]
+		//-------------------------------------------
+		gdr[is].resize(nrxx);
+		XC_Functional::grad_rho(rhog.data(), gdr[is].data(), chr->rhopw, tpiba);
+	} // end for(is)
+	return gdr;
+}
+
+// converting grho (abacus=>libxc)
+std::vector<double> XC_Functional_Libxc::convert_sigma(
+	const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr)
+{
+	const std::size_t nspin = gdr.size();
+	assert(nspin>0);
+	const std::size_t nrxx = gdr[0].size();
+	for(std::size_t is=1; is<nspin; ++is)
+		assert(nrxx==gdr[is].size());
+
+	std::vector<double> sigma( nrxx * ((1==nspin)?1:3) );
+	if( 1==nspin )
+	{
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 1024)
+		#endif
+		for( std::size_t ir=0; ir<nrxx; ++ir )
+			sigma[ir] = gdr[0][ir]*gdr[0][ir];
+	}
+	else
+	{
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 256)
+		#endif
+		for( std::size_t ir=0; ir<nrxx; ++ir )
+		{
+			sigma[ir*3]   = gdr[0][ir]*gdr[0][ir];
+			sigma[ir*3+1] = gdr[0][ir]*gdr[1][ir];
+			sigma[ir*3+2] = gdr[1][ir]*gdr[1][ir];
+		}
+	}
+	return sigma;
+}
+
+// sgn for threshold mask
+std::vector<double> XC_Functional_Libxc::cal_sgn(
+	const double rho_threshold,
+	const double grho_threshold,
+	const xc_func_type &func,
+	const int nspin,
+	const std::size_t nrxx,
+	const std::vector<double> &rho,
+	const std::vector<double> &sigma)
+{
+	std::vector<double> sgn(nrxx*nspin, 1.0);
+	// 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)
+	{
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 512)
+		#endif
+		for( int ir=0; ir<nrxx; ++ir )
+		{
+			if ( rho[ir*2]<rho_threshold || std::sqrt(std::abs(sigma[ir*3]))<grho_threshold )
+				sgn[ir*2] = 0.0;
+			if ( rho[ir*2+1]<rho_threshold || std::sqrt(std::abs(sigma[ir*3+2]))<grho_threshold )
+				sgn[ir*2+1] = 0.0;
+		}
+	}
+	return sgn;
+}
+
+// converting etxc from exc (libxc=>abacus)
+double XC_Functional_Libxc::convert_etxc(
+	const int nspin,
+	const std::size_t nrxx,
+	const std::vector<double> &sgn,
+	const std::vector<double> &rho,
+	std::vector<double> exc)
+{
+	double etxc = 0.0;
+	#ifdef _OPENMP
+	#pragma omp parallel for collapse(2) reduction(+:etxc) schedule(static, 256)
+	#endif
+	for( int is=0; is<nspin; ++is )
+		for( int ir=0; ir<nrxx; ++ir )
+			etxc += ModuleBase::e2 * exc[ir] * rho[ir*nspin+is] * sgn[ir*nspin+is];
+	return etxc;
+}
+
+// converting etxc from exc (libxc=>abacus)
+double XC_Functional_Libxc::convert_vtxc_v(
+	const xc_func_type &func,
+	const int nspin,
+	const std::size_t nrxx,
+	const std::vector<double> &sgn,
+	const std::vector<double> &rho,
+	const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr,
+	const std::vector<double> &vrho,
+	const std::vector<double> &vsigma,
+	const double tpiba,
+	const Charge* const chr,
+	ModuleBase::matrix &v)
+{
+	assert(v.nr==nspin);
+	assert(v.nc==nrxx);
+
+	double vtxc = 0.0;
+
+	#ifdef _OPENMP
+	#pragma omp parallel for collapse(2) reduction(+:vtxc) schedule(static, 256)
+	#endif
+	for( int is=0; is<nspin; ++is )
+	{
+		for( std::size_t ir=0; ir<nrxx; ++ir )
+		{
+			const double v_tmp = ModuleBase::e2 * vrho[ir*nspin+is] * sgn[ir*nspin+is];
+			v(is,ir) += v_tmp;
+			vtxc += v_tmp * rho[ir*nspin+is];
+		}
+	}
+
+	if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)
+	{
+		const std::vector<std::vector<double>> dh = XC_Functional_Libxc::cal_dh(nspin, nrxx, sgn, gdr, vsigma, tpiba, chr);
+
+		double rvtxc = 0.0;
+		#ifdef _OPENMP
+		#pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256)
+		#endif
+		for( int is=0; is<nspin; ++is )
+		{
+			for( std::size_t ir=0; ir<nrxx; ++ir )
+			{
+				rvtxc += dh[is][ir] * rho[ir*nspin+is];
+				v(is,ir) -= dh[is][ir];
+			}
+		}
+
+		vtxc -= rvtxc;
+	} // end if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA))
+
+	return vtxc;
+}
+
+
+// dh for gga v
+std::vector<std::vector<double>> XC_Functional_Libxc::cal_dh(
+	const int nspin,
+	const std::size_t nrxx,
+	const std::vector<double> &sgn,
+	const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr,
+	const std::vector<double> &vsigma,
+	const double tpiba,
+	const Charge* const chr)
+{
+	std::vector<std::vector<ModuleBase::Vector3<double>>> h(
+		nspin,
+		std::vector<ModuleBase::Vector3<double>>(nrxx) );
+	if( 1==nspin )
+	{
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 1024)
+		#endif
+		for( std::size_t ir=0; ir<nrxx; ++ir )
+			h[0][ir] = 2.0 * gdr[0][ir] * vsigma[ir] * 2.0 * sgn[ir];
+	}
+	else
+	{
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 1024)
+		#endif
+		for( std::size_t ir=0; ir< nrxx; ++ir )
+		{
+			h[0][ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3  ] * sgn[ir*2  ] * 2.0
+							+ gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2]   * sgn[ir*2+1]);
+			h[1][ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0
+							+ gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2]   * sgn[ir*2+1]);
+		}
+	}
+
+	// define two dimensional array dh [ nspin, nrxx ]
+	std::vector<std::vector<double>> dh(nspin, std::vector<double>(nrxx));
+	for( int is=0; is!=nspin; ++is )
+		XC_Functional::grad_dot( h[is].data(), dh[is].data(), chr->rhopw, tpiba);
+
+	return dh;
+}
+
+
+// convert v for NSPIN=4
+ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4(
+	const std::size_t nrxx,
+	const Charge* const chr,
+	const std::vector<double> &amag,
+	const ModuleBase::matrix &v)
+{
+	assert(GlobalV::NSPIN==4);
+	constexpr double vanishing_charge = 1.0e-10;
+	ModuleBase::matrix v_nspin4(GlobalV::NSPIN, nrxx);
+	for( int ir=0; ir<nrxx; ++ir )
+		v_nspin4(0,ir) = 0.5 * (v(0,ir)+v(1,ir));
+	if(GlobalV::DOMAG || GlobalV::DOMAG_Z)
+	{
+		for( int ir=0; ir<nrxx; ++ir )
+		{
+			if ( amag[ir] > vanishing_charge )
+			{
+				const double vs = 0.5 * (v(0,ir)-v(1,ir));
+				for(int ipol=1; ipol<GlobalV::NSPIN; ++ipol)
+					v_nspin4(ipol,ir) = vs * chr->rho[ipol][ir] / amag[ir];
+			}
+		}
+	}
+	return v_nspin4;
+}
+
+#endif
\ No newline at end of file
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp
index 9e8614dfa6..e8b80537f1 100644
--- a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp
+++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp
@@ -260,9 +260,9 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
             for( int ir=0; ir< nrxx; ++ir )
             {
 #ifdef __EXX
-                if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5)
+                if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
                 {
-                    exc[ir] *= (1.0 - XC_Functional::hybrid_alpha);
+                    exc[ir] *= (1.0 - XC_Functional::get_hybrid_alpha());
                 }
 #endif
                 etxc += ModuleBase::e2 * exc[ir] * rho[ir*nspin+is]  * sgn[ir*nspin+is];
@@ -278,9 +278,9 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
             for( int ir=0; ir< nrxx; ++ir )
             {
 #ifdef __EXX
-                if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5)
+                if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
                 {
-                    vrho[ir*nspin+is] *= (1.0 - XC_Functional::hybrid_alpha);
+                    vrho[ir*nspin+is] *= (1.0 - XC_Functional::get_hybrid_alpha());
                 }
 #endif
                 const double v_tmp = ModuleBase::e2 * vrho[ir*nspin+is]  * sgn[ir*nspin+is];
@@ -301,9 +301,9 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
             for( int ir=0; ir< nrxx; ++ir )
             {
 #ifdef __EXX
-                if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5)
+                if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
                 {
-                    vsigma[ir] *= (1.0 - XC_Functional::hybrid_alpha);
+                    vsigma[ir] *= (1.0 - XC_Functional::get_hybrid_alpha());
                 }
 #endif
                 h[0][ir] = 2.0 * gdr[0][ir] * vsigma[ir] * 2.0 * sgn[ir];
@@ -317,11 +317,11 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
             for( int ir=0; ir< nrxx; ++ir )
             {
 #ifdef __EXX
-                if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5)
+                if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
                 {
-                    vsigma[ir*3]   *= (1.0 - XC_Functional::hybrid_alpha);
-                    vsigma[ir*3+1] *= (1.0 - XC_Functional::hybrid_alpha);
-                    vsigma[ir*3+2] *= (1.0 - XC_Functional::hybrid_alpha);
+                    vsigma[ir*3]   *= (1.0 - XC_Functional::get_hybrid_alpha());
+                    vsigma[ir*3+1] *= (1.0 - XC_Functional::get_hybrid_alpha());
+                    vsigma[ir*3+2] *= (1.0 - XC_Functional::get_hybrid_alpha());
                 }
 #endif
                 h[0][ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3  ] * sgn[ir*2  ] * 2.0
@@ -363,9 +363,9 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
             for( int ir=0; ir< nrxx; ++ir )
             {
 #ifdef __EXX
-                if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5)
+                if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
                 {
-                    vtau[ir*nspin+is] *= (1.0 - XC_Functional::hybrid_alpha);
+                    vtau[ir*nspin+is] *= (1.0 - XC_Functional::get_hybrid_alpha());
                 }
 #endif
                 vofk(is,ir) += vtau[ir*nspin+is]  * sgn[ir*nspin+is];
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp
index 0826f10c5f..48bd301ae6 100644
--- a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp
+++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp
@@ -25,7 +25,7 @@ void XC_Functional_Libxc::tau_xc(
     {
         xc_mgga_exc_vxc(&func,1,&rho,&grho,&lapl_rho,&atau,&s,&v1,&v2,&vlapl_rho,&v3);
 #ifdef __EXX
-        if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5)
+        if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
         {
             s *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha);
             v1 *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha);
diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp
index dd6063d5d0..63e4539e3f 100644
--- a/source/module_io/input_conv.cpp
+++ b/source/module_io/input_conv.cpp
@@ -469,7 +469,7 @@ void Input_Conv::Convert()
         // EXX case, convert all EXX related variables
         // GlobalC::exx_info.info_global.cal_exx = true;
         GlobalC::exx_info.info_global.hybrid_alpha = std::stod(PARAM.inp.exx_hybrid_alpha);
-        XC_Functional::get_hybrid_alpha(std::stod(PARAM.inp.exx_hybrid_alpha));
+        XC_Functional::set_hybrid_alpha(std::stod(PARAM.inp.exx_hybrid_alpha));
         GlobalC::exx_info.info_global.hse_omega = PARAM.inp.exx_hse_omega;
         GlobalC::exx_info.info_global.separate_loop = PARAM.inp.exx_separate_loop;
         GlobalC::exx_info.info_global.hybrid_step = PARAM.inp.exx_hybrid_step;