Fix #6286: avoid redundant libxc init_func calls in per-point evaluation paths

source/source_hamilt/module_xc/libxc_abacus.h

@@ -73,7 +73,6 @@ namespace XC_Functional_Libxc
 
     extern void finish_func(std::vector<xc_func_type> &funcs);
 
-
 //-------------------
 //  libxc_pot.cpp
 //-------------------
@@ -282,6 +281,17 @@ namespace XC_Functional_Libxc
         const double hybrid_alpha,
         const double hse_omega);
 
+    // Overload accepting an already-initialized functional vector. The caller
+    // is responsible for init_func/finish_func; this overload does not touch
+    // the lifetime of funcs. Useful for per-thread reuse inside OpenMP loops.
+    extern void gcxc_libxc(
+        const std::vector<xc_func_type>& funcs,
+        const double &rho,
+        const double &grho,
+        double &sxc,
+        double &v1xc,
+        double &v2xc);
+
     // the entire GGA functional, for nspin=2 case
     extern void gcxc_spin_libxc(
         const std::vector<int> &func_id,
@@ -298,6 +308,20 @@ namespace XC_Functional_Libxc
         const double hybrid_alpha,
         const double hse_omega);
 
+    // Overload accepting an already-initialized functional vector.
+    extern void gcxc_spin_libxc(
+        const std::vector<xc_func_type>& funcs,
+        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);
+
 
 //-------------------
 //  libxc_mgga_wrap.cpp
@@ -318,6 +342,21 @@ namespace XC_Functional_Libxc
         const double &hybrid_alpha,
         const double &hse_omega);
 
+    // Overload accepting an already-initialized functional vector.
+    // hybrid_alpha scales the semilocal SCAN exchange for SCAN0.
+    extern void tau_xc(
+        const std::vector<xc_func_type>& funcs,
+        const double &rho,
+        const double &grho,
+        const double &lapl_rho,
+        const double &atau,
+        double &sxc,
+        double &v1xc,
+        double &v2xc,
+        double &v3xc,
+        double &vlaplxc,
+        const double hybrid_alpha);
+
     extern void tau_xc_spin(
         const std::vector<int> &func_id,
         double rhoup,
@@ -341,6 +380,30 @@ namespace XC_Functional_Libxc
         const double &hybrid_alpha,
         const double &hse_omega);
 
+    // Overload accepting an already-initialized functional vector.
+    // hybrid_alpha scales the semilocal SCAN exchange for SCAN0.
+    extern void tau_xc_spin(
+        const std::vector<xc_func_type>& funcs,
+        double rhoup,
+        double rhodw,
+        ModuleBase::Vector3<double> gdr1,
+        ModuleBase::Vector3<double> gdr2,
+        double laplup,
+        double lapldw,
+        double tauup,
+        double taudw,
+        double &sxc,
+        double &v1xcup,
+        double &v1xcdw,
+        double &v2xcup,
+        double &v2xcdw,
+        double &v2xcud,
+        double &v3xcup,
+        double &v3xcdw,
+        double &vlaplxcup,
+        double &vlaplxcdw,
+        const double hybrid_alpha);
+
 } // namespace XC_Functional_Libxc
 
 #endif // __LIBXC

source/source_hamilt/module_xc/libxc_gga_wrap.cpp

@@ -8,34 +8,35 @@
 #include <xc.h>
 #include <array>
 
+namespace
+{
+// gcxc_libxc returns zero below these thresholds.
+bool gcxc_is_negligible(const double rho, const double grho)
+{
+    constexpr double small = 1.e-6;
+    constexpr double smallg = 1.e-10;
+    return rho <= small || grho < smallg;
+}
+} // namespace
+
 void XC_Functional_Libxc::gcxc_libxc(
-    const std::vector<int>& func_id,
+    const std::vector<xc_func_type>& funcs,
     const double& rho,
     const double& grho,
     double& sxc,
     double& v1xc,
-    double& v2xc,
-    const double hybrid_alpha,
-    const double hse_omega)
+    double& v2xc)
 {
     sxc = 0.0;
     v1xc = 0.0;
     v2xc = 0.0;
 
-    constexpr double small = 1.e-6;
-    constexpr double smallg = 1.e-10;
-    if (rho <= small || grho < smallg)
+    if (gcxc_is_negligible(rho, grho))
     {
         return;
     }
 
-    std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(
-        /* func_id = */ func_id,
-        /* xc_polarized = */ XC_UNPOLARIZED,
-        /* hybrid_alpha = */ hybrid_alpha,
-        /* hse_omega = */ hse_omega);
-
-    for (xc_func_type& func : funcs)
+    for (const xc_func_type& func : funcs)
     {
         double s = 0.0;
         double v1 = 0.0;
@@ -45,14 +46,42 @@ void XC_Functional_Libxc::gcxc_libxc(
         v1xc += v1;
         v2xc += v2 * 2.0;
     }
+}
+
+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 double hybrid_alpha,
+    const double hse_omega)
+{
+    // Skip the libxc setup entirely for points the evaluation would zero out.
+    if (gcxc_is_negligible(rho, grho))
+    {
+        sxc = 0.0;
+        v1xc = 0.0;
+        v2xc = 0.0;
+        return;
+    }
+
+    std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(
+        /* func_id = */ func_id,
+        /* xc_polarized = */ XC_UNPOLARIZED,
+        /* hybrid_alpha = */ hybrid_alpha,
+        /* hse_omega = */ hse_omega);
+
+    gcxc_libxc(funcs, rho, grho, sxc, v1xc, v2xc);
 
     XC_Functional_Libxc::finish_func(funcs);
 } // end subroutine gcxc_libxc
 
 
 
 void XC_Functional_Libxc::gcxc_spin_libxc(
-    const std::vector<int>& func_id,
+    const std::vector<xc_func_type>& funcs,
     const double rhoup,
     const double rhodw,
     const ModuleBase::Vector3<double> gdr1,
@@ -62,9 +91,7 @@ void XC_Functional_Libxc::gcxc_spin_libxc(
     double& v1xcdw,
     double& v2xcup,
     double& v2xcdw,
-    double& v2xcud,
-    const double hybrid_alpha,
-    const double hse_omega)
+    double& v2xcud)
 {
     sxc = 0.0;
     v1xcup = 0.0;
@@ -75,13 +102,7 @@ void XC_Functional_Libxc::gcxc_spin_libxc(
     const std::array<double, 2> rho = {rhoup, rhodw};
     const std::array<double, 3> grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()};
 
-    std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(
-        /* func_id = */ func_id,
-        /* xc_polarized = */ XC_POLARIZED,
-        /* hybrid_alpha = */ hybrid_alpha,
-        /* hse_omega = */ hse_omega);
-
-    for (xc_func_type& func : funcs)
+    for (const xc_func_type& func : funcs)
     {
         if (func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)
         {
@@ -113,6 +134,32 @@ void XC_Functional_Libxc::gcxc_spin_libxc(
             v2xcdw += 2.0 * v2xc[2] * sgn[1];
         }
     }
+}
+
+void XC_Functional_Libxc::gcxc_spin_libxc(
+    const std::vector<int>& func_id,
+    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,
+    const double hybrid_alpha,
+    const double hse_omega)
+{
+    std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(
+        /* func_id = */ func_id,
+        /* xc_polarized = */ XC_POLARIZED,
+        /* hybrid_alpha = */ hybrid_alpha,
+        /* hse_omega = */ hse_omega);
+
+    gcxc_spin_libxc(funcs, rhoup, rhodw, gdr1, gdr2,
+                    sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud);
+
     XC_Functional_Libxc::finish_func(funcs);
 }