Refactor little codes in XC_Functional_Libxc

source/module_hamilt_general/module_xc/test/test_xc.cpp

@@ -1,5 +1,6 @@
 #include "gtest/gtest.h"
 #include "../xc_functional.h"
+#include "../xc_functional_libxc.h"
 #include "../exx_info.h"
 #include "xctest.h"
 
@@ -829,7 +830,7 @@ class XCTest_PBE0 : public XCTest
         void SetUp()
         {
             XC_Functional::set_xc_type("PBE0");
-            XC_Functional::get_hybrid_alpha(0.5);
+            XC_Functional::set_hybrid_alpha(0.5);
             std::vector<double> rho  = {0.17E+01, 0.17E+01, 0.15E+01, 0.88E-01, 0.18E+04};
             std::vector<double> grho = {0.81E-11, 0.17E+01, 0.36E+02, 0.87E-01, 0.55E+00};
 
@@ -884,7 +885,7 @@ class XCTest_PBE_LibXC : public XCTest
                 XC_Functional::xc(rho[i],e,v);
                 e_lda.push_back(e);
                 v_lda.push_back(v);
-                XC_Functional::gcxc_libxc(rho[i],grho[i],e,v1,v2);
+                XC_Functional_Libxc::gcxc_libxc(XC_Functional::get_func_id(), rho[i],grho[i],e,v1,v2);
                 e_gga.push_back(e);
                 v1_gga.push_back(v1);
                 v2_gga.push_back(v2);

source/module_hamilt_general/module_xc/test/test_xc1.cpp

@@ -38,7 +38,7 @@ class XCTest_HSE : public XCTest
         void SetUp()
         {
             XC_Functional::set_xc_type("HSE");
-            XC_Functional::get_hybrid_alpha(0.5);
+            XC_Functional::set_hybrid_alpha(0.5);
         }
 };
 
@@ -53,7 +53,7 @@ class XCTest_SCAN0 : public XCTest
         void SetUp()
         {
             XC_Functional::set_xc_type("SCAN0");
-            XC_Functional::get_hybrid_alpha(0.5);
+            XC_Functional::set_hybrid_alpha(0.5);
         }
 };
 

source/module_hamilt_general/module_xc/test/test_xc2.cpp

@@ -1,6 +1,7 @@
 #include "gtest/gtest.h"
 #include "xctest.h"
 #include "../xc_functional.h"
+#include "../xc_functional_libxc.h"
 #include "../exx_info.h"
 /************************************************
 *  unit test of functionals
@@ -455,7 +456,7 @@ class XCTest_PBE_SPN_LibXC : public XCTest
                 double e,v1,v2,v3,v4,v5;
                 double r1 = rho[i] * (1+zeta[i]) / 2.0;
                 double r2 = rho[i] * (1-zeta[i]) / 2.0;
-                XC_Functional::gcxc_spin_libxc(r1,r2,gdr[i],gdr[i],e,v1,v2,v3,v4,v5);
+                XC_Functional_Libxc::gcxc_spin_libxc(XC_Functional::get_func_id(), r1,r2,gdr[i],gdr[i],e,v1,v2,v3,v4,v5);
                 e_gga.push_back(e);
                 v1_gga.push_back(v1+v3);
                 v2_gga.push_back(v2+v4);
@@ -494,7 +495,7 @@ class XCTest_PZ_SPN_LibXC : public XCTest
                 double e,v1,v2;
                 double r1 = rho[i] * (1+zeta[i]) / 2.0;
                 double r2 = rho[i] * (1-zeta[i]) / 2.0;
-                XC_Functional::xc_spin_libxc(r1,r2,e,v1,v2);
+                XC_Functional_Libxc::xc_spin_libxc(XC_Functional::get_func_id(), r1,r2,e,v1,v2);
                 e_lda.push_back(e);
                 v1_lda.push_back(v1);
                 v2_lda.push_back(v2);

source/module_hamilt_general/module_xc/test/test_xc4.cpp

@@ -1,4 +1,5 @@
 #include "../xc_functional.h"
+#include "../xc_functional_libxc.h"
 #include "gtest/gtest.h"
 #include "xctest.h"
 #include "../exx_info.h"
@@ -43,7 +44,7 @@ class XCTest_SCAN : public XCTest
             for(int i=0;i<5;i++)
             {
                 double e,v,v1,v2,v3;
-                XC_Functional::tau_xc(rho[i],grho[i],tau[i],e,v1,v2,v3);
+                XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3);
                 e_.push_back(e);
                 v1_.push_back(v1);
                 v2_.push_back(v2);

source/module_hamilt_general/module_xc/test/test_xc5.cpp

@@ -1,4 +1,5 @@
 #include "../xc_functional.h"
+#include "../xc_functional_libxc.h"
 #include "gtest/gtest.h"
 #include "xctest.h"
 #include "../exx_info.h"
@@ -279,15 +280,15 @@ class XCTest_VXC_meta : public XCTest
 
             GlobalV::NSPIN = 1;
             std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix> etxc_vtxc_v
-                = XC_Functional::v_xc_meta(rhopw.nrxx,ucell.omega,ucell.tpiba,&chr);
+                = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr);
             et1 = std::get<0>(etxc_vtxc_v);
             vt1 = std::get<1>(etxc_vtxc_v);
             v1  = std::get<2>(etxc_vtxc_v);
             vtau1 = std::get<3>(etxc_vtxc_v);
 
             GlobalV::NSPIN = 2;
             etxc_vtxc_v
-                = XC_Functional::v_xc_meta(rhopw.nrxx,ucell.omega,ucell.tpiba,&chr);
+                = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr);
             et2 = std::get<0>(etxc_vtxc_v);
             vt2 = std::get<1>(etxc_vtxc_v);
             v2  = std::get<2>(etxc_vtxc_v);

source/module_hamilt_general/module_xc/xc_functional.cpp

@@ -58,7 +58,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
     //such as for PBE we have:
     //        func_id.push_back(XC_GGA_X_PBE);
     //        func_id.push_back(XC_GGA_C_PBE);
-    
+
     func_id.clear();
     std::string xc_func = xc_func_in;
     std::transform(xc_func.begin(), xc_func.end(), xc_func.begin(), (::toupper));
@@ -128,7 +128,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
         func_id.push_back(XC_GGA_C_PBE);
         func_type = 2;
         use_libxc = false;
-	}	
+	}
 	else if ( xc_func == "BLYP") //B88+LYP
 	{
         func_id.push_back(XC_GGA_X_B88);
@@ -142,14 +142,14 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
         func_id.push_back(XC_GGA_C_P86);
         func_type = 2;
         use_libxc = false;
-	} 
+	}
 	else if ( xc_func == "PW91") //PW91_X+PW91_C
 	{
         func_id.push_back(XC_GGA_X_PW91);
         func_id.push_back(XC_GGA_C_PW91);
         func_type = 2;
         use_libxc = false;
-	} 
+	}
 	else if ( xc_func == "HCTH") //HCTH_X+HCTH_C
 	{
         func_id.push_back(XC_GGA_X_HCTH_A);

source/module_hamilt_general/module_xc/xc_functional_libxc.cpp

@@ -50,7 +50,7 @@ std::vector<xc_func_type> XC_Functional_Libxc::init_func(const std::vector<int>
 	std::vector<xc_func_type> funcs;
 
 	//-------------------------------------------
-	// define a function named 'add_func', which 
+	// define a function named 'add_func', which
 	// initialize a functional according to its ID
 	//-------------------------------------------
 	auto add_func = [&]( const int func_id )
@@ -71,17 +71,17 @@ std::vector<xc_func_type> XC_Functional_Libxc::init_func(const std::vector<int>
 #ifdef __EXX
 		else if( id == XC_HYB_GGA_XC_PBEH ) // PBE0
 		{
-			add_func( XC_HYB_GGA_XC_PBEH );		
-			double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, 
-				GlobalC::exx_info.info_global.hse_omega, 
+			add_func( XC_HYB_GGA_XC_PBEH );
+			double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha,
+				GlobalC::exx_info.info_global.hse_omega,
 				GlobalC::exx_info.info_global.hse_omega };
-			xc_func_set_ext_params(&funcs.back(), parameter_hse);	
+			xc_func_set_ext_params(&funcs.back(), parameter_hse);
 		}
 		else if( id == XC_HYB_GGA_XC_HSE06 ) // HSE06 hybrid functional
 		{
-			add_func( XC_HYB_GGA_XC_HSE06 );	
-			double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, 
-				GlobalC::exx_info.info_global.hse_omega, 
+			add_func( XC_HYB_GGA_XC_HSE06 );
+			double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha,
+				GlobalC::exx_info.info_global.hse_omega,
 				GlobalC::exx_info.info_global.hse_omega };
 			xc_func_set_ext_params(&funcs.back(), parameter_hse);
 		}

source/module_hamilt_general/module_xc/xc_functional_libxc.h

@@ -149,8 +149,8 @@ namespace XC_Functional_Libxc
 	// the entire GGA functional, for nspin=2 case
 	extern void gcxc_spin_libxc(
 		const std::vector<int> &func_id,
-		double rhoup, double rhodw,
-		ModuleBase::Vector3<double> gdr1, ModuleBase::Vector3<double> gdr2,
+		const double rhoup, const double rhodw,
+		const ModuleBase::Vector3<double> gdr1, const ModuleBase::Vector3<double> gdr2,
 		double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud);
 
 

source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp

@@ -205,9 +205,9 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
     std::vector<double> vtau   ( nrxx * nspin            );
     std::vector<double> vlapl  ( nrxx * nspin            );
 
-    const double rho_th  = 1e-8;
-    const double grho_th = 1e-12;
-    const double tau_th  = 1e-8;
+    constexpr double rho_th  = 1e-8;
+    constexpr double grho_th = 1e-12;
+    constexpr double tau_th  = 1e-8;
     // sgn for threshold mask
     std::vector<double> sgn( nrxx * nspin);
 #ifdef _OPENMP

source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp

@@ -5,73 +5,64 @@
 #include <xc.h>
 
 void XC_Functional_Libxc::gcxc_libxc(
-        const std::vector<int> &func_id,
-    const double &rho, const double &grho, double &sxc,
-          double &v1xc, double &v2xc)
+    const std::vector<int> &func_id,
+    const double &rho, const double &grho,
+    double &sxc, double &v1xc, double &v2xc)
 {
-    const double small = 1.e-6;
-    const double smallg = 1.e-10;
-    double s,v1,v2;
     sxc = v1xc = v2xc = 0.0;
 
+    constexpr double small = 1.e-6;
+    constexpr double smallg = 1.e-10;
     if (rho <= small || grho < smallg)
     {
         return;
     }
 
     std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED);
-
-
     for(xc_func_type &func : funcs)
     {
+        double s,v1,v2;
         xc_gga_exc_vxc(&func, 1, &rho, &grho, &s, &v1, &v2);
-
         sxc += s * rho;
-        v2xc += v2 * 2.0;
         v1xc += v1;
+        v2xc += v2 * 2.0;
     }
     XC_Functional_Libxc::finish_func(funcs);
 } // end subroutine gcxc_libxc
 
+
+
 void XC_Functional_Libxc::gcxc_spin_libxc(
         const std::vector<int> &func_id,
-        double rhoup, double rhodw, 
-        ModuleBase::Vector3<double> gdr1, ModuleBase::Vector3<double> gdr2,
+        const double rhoup, const double rhodw,
+        const ModuleBase::Vector3<double> gdr1, const ModuleBase::Vector3<double> gdr2,
         double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud)
 {
-	std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED);
-    double *rho, *grho, *v1xc, *v2xc, *sgn, s;
     sxc = v1xcup = v1xcdw = 0.0;
     v2xcup = v2xcdw = v2xcud = 0.0;
-    rho = new double[2];
-    grho= new double[3];
-    v1xc= new double[2];
-    v2xc= new double[3];
-    sgn = new double[2];
-
-    rho[0] = rhoup;
-    rho[1] = rhodw;
-    grho[0] = gdr1.norm2();
-    grho[1] = gdr1 * gdr2;
-    grho[2] = gdr2.norm2();
-
-    const double rho_threshold = 1E-6;
-    const double grho_threshold = 1E-10;
+    const std::vector<double> rho = {rhoup, rhodw};
+    const std::vector<double> grho = {gdr1.norm2(), gdr1*gdr2, gdr2.norm2()};
 
+    std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED);
     for(xc_func_type &func : funcs)
     {
         if( func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)
         {
-            sgn[0] = sgn[1] = 1.0;
-            // call Libxc function: xc_gga_exc_vxc
-            xc_gga_exc_vxc( &func, 1, rho, grho, &s, v1xc, v2xc);
+            constexpr double rho_threshold = 1E-6;
+            constexpr double grho_threshold = 1E-10;
+            std::vector<double> sgn = {1.0, 1.0};
             if(func.info->kind==XC_CORRELATION)
             {
                 if ( rho[0]<rho_threshold || sqrt(std::abs(grho[0]))<grho_threshold )
                     sgn[0] = 0.0;
                 if ( rho[1]<rho_threshold || sqrt(std::abs(grho[2]))<grho_threshold )
                     sgn[1] = 0.0;
             }
+
+            double s = 0.0;
+            std::vector<double> v1xc(2), v2xc(3);
+            // call Libxc function: xc_gga_exc_vxc
+            xc_gga_exc_vxc( &func, 1, rho.data(), grho.data(), &s, v1xc.data(), v2xc.data());
             sxc += s * (rho[0] * sgn[0] + rho[1] * sgn[1]);
             v1xcup += v1xc[0] * sgn[0];
             v1xcdw += v1xc[1] * sgn[1];
@@ -80,11 +71,6 @@ void XC_Functional_Libxc::gcxc_spin_libxc(
             v2xcdw += 2.0 * v2xc[2] * sgn[1];
         }
     }
-    delete[] grho;
-    delete[] rho;
-    delete[] v1xc;
-    delete[] v2xc;
-    delete[] sgn;
     XC_Functional_Libxc::finish_func(funcs);
 }