fix(xc): close gga_grad 2 builtin and LibXC stress derivatives

Based-on: #7758
Co-authored-by: dyzheng <zhengdy@aisi.ac.cn>

Author: aboys-cb

source/source_hamilt/module_xc/libxc_pot.cpp

@@ -16,6 +16,170 @@
 #include <vector>
 #include <complex>
 
+void XC_Functional_Libxc::gradcorr_ncgga_sf_libxc(const std::vector<int>& func_id,
+                                                  const std::size_t nrxx,
+                                                  const double tpiba,
+                                                  const Charge* const chr,
+                                                  const std::map<int, double>* scaling_factor,
+                                                  const double hybrid_alpha,
+                                                  const double hse_omega,
+                                                  std::vector<double>& stress_gga)
+{
+    constexpr int nspin = 2;
+    stress_gga.assign(9, 0.0);
+
+    std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED, hybrid_alpha, hse_omega);
+    bool has_gga = false;
+    for (const xc_func_type& func: funcs)
+    {
+        has_gga = has_gga || func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA;
+    }
+    if (!has_gga)
+    {
+        XC_Functional_Libxc::finish_func(funcs);
+        return;
+    }
+
+    // This is the same forward graph used by v_xc_libxc: the local spin map,
+    // its projected FFT gradients, and the sigma invariants are constructed
+    // once and shared by all Libxc components.
+    const XC_Functional_Libxc::NclSfDiscreteData sf_data
+        = XC_Functional_Libxc::make_ncl_sf_discrete_data(nrxx, tpiba, chr, true);
+    const std::vector<double>& rho = sf_data.rho;
+    const std::vector<double> sigma = XC_Functional_Libxc::convert_sigma(sf_data.spin_gradient);
+    std::vector<double> aggregate_dsigma(3 * nrxx, 0.0);
+
+    for (xc_func_type& func: funcs)
+    {
+        if (func.info->family != XC_FAMILY_GGA && func.info->family != XC_FAMILY_HYB_GGA)
+        {
+            continue;
+        }
+
+        constexpr double rho_threshold = 1.0e-6;
+        constexpr double grho_threshold = 1.0e-10;
+        xc_func_set_dens_threshold(&func, rho_threshold);
+        const std::vector<double> sgn
+            = XC_Functional_Libxc::cal_sgn(rho_threshold, grho_threshold, func, nspin, nrxx, rho, sigma);
+        std::vector<double> exc(nrxx);
+        std::vector<double> vrho(nspin * nrxx);
+        std::vector<double> vsigma(3 * nrxx);
+        constexpr int nr_batch_size = 1024;
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, nr_batch_size)
+#endif
+        for (int ir_start = 0; ir_start < static_cast<int>(nrxx); ir_start += nr_batch_size)
+        {
+            const int ir_end = std::min(ir_start + nr_batch_size, static_cast<int>(nrxx));
+            const int nrxx_thread = ir_end - ir_start;
+            xc_gga_exc_vxc(&func,
+                           nrxx_thread,
+                           rho.data() + ir_start * nspin,
+                           sigma.data() + ir_start * 3,
+                           exc.data() + ir_start,
+                           vrho.data() + ir_start * nspin,
+                           vsigma.data() + ir_start * 3);
+        }
+
+        double factor = 1.0;
+        if (scaling_factor != nullptr)
+        {
+            const std::map<int, double>::const_iterator entry = scaling_factor->find(func.info->number);
+            if (entry != scaling_factor->end())
+            {
+                factor = entry->second;
+            }
+        }
+        const XC_Functional_Libxc::LibxcWeightedDerivatives weighted
+            = XC_Functional_Libxc::make_libxc_weighted_derivatives(func,
+                                                                   nspin,
+                                                                   nrxx,
+                                                                   sgn,
+                                                                   rho,
+                                                                   sigma,
+                                                                   exc,
+                                                                   vrho,
+                                                                   vsigma);
+        for (std::size_t index = 0; index < aggregate_dsigma.size(); ++index)
+        {
+            aggregate_dsigma[index] += factor * weighted.dsigma[index];
+        }
+    }
+
+// For g_s=sum_A J_sA G_h x_A, a reciprocal deformation changes G_h but
+// not the pointwise map J.  Therefore the metric derivative is exactly
+// sum_s h_s,l g_s,m, with h_s=dE/dg_s built from the sanitizer-reversed,
+// component-scaled aggregate above.
+#ifdef _OPENMP
+#pragma omp parallel
+    {
+        std::vector<double> local_stress(9, 0.0);
+#pragma omp for schedule(static, 512)
+        for (std::size_t ir = 0; ir < nrxx; ++ir)
+        {
+            const std::size_t sigma_index = 3 * ir;
+            const ModuleBase::Vector3<double>& grad_up = sf_data.spin_gradient[0][ir];
+            const ModuleBase::Vector3<double>& grad_down = sf_data.spin_gradient[1][ir];
+            const ModuleBase::Vector3<double> h_up
+                = ModuleBase::e2
+                  * (2.0 * aggregate_dsigma[sigma_index] * grad_up + aggregate_dsigma[sigma_index + 1] * grad_down);
+            const ModuleBase::Vector3<double> h_down
+                = ModuleBase::e2
+                  * (2.0 * aggregate_dsigma[sigma_index + 2] * grad_down + aggregate_dsigma[sigma_index + 1] * grad_up);
+            const double grad_up_component[3] = {grad_up.x, grad_up.y, grad_up.z};
+            const double grad_down_component[3] = {grad_down.x, grad_down.y, grad_down.z};
+            const double h_up_component[3] = {h_up.x, h_up.y, h_up.z};
+            const double h_down_component[3] = {h_down.x, h_down.y, h_down.z};
+            for (int l = 0; l < 3; ++l)
+            {
+                for (int m = 0; m <= l; ++m)
+                {
+                    local_stress[l * 3 + m]
+                        += h_up_component[l] * grad_up_component[m] + h_down_component[l] * grad_down_component[m];
+                }
+            }
+        }
+#pragma omp critical(libxc_ncgga_stress_reduce)
+        {
+            for (int l = 0; l < 3; ++l)
+            {
+                for (int m = 0; m <= l; ++m)
+                {
+                    stress_gga[l * 3 + m] += local_stress[l * 3 + m];
+                }
+            }
+        }
+    }
+#else
+    for (std::size_t ir = 0; ir < nrxx; ++ir)
+    {
+        const std::size_t sigma_index = 3 * ir;
+        const ModuleBase::Vector3<double>& grad_up = sf_data.spin_gradient[0][ir];
+        const ModuleBase::Vector3<double>& grad_down = sf_data.spin_gradient[1][ir];
+        const ModuleBase::Vector3<double> h_up
+            = ModuleBase::e2
+              * (2.0 * aggregate_dsigma[sigma_index] * grad_up + aggregate_dsigma[sigma_index + 1] * grad_down);
+        const ModuleBase::Vector3<double> h_down
+            = ModuleBase::e2
+              * (2.0 * aggregate_dsigma[sigma_index + 2] * grad_down + aggregate_dsigma[sigma_index + 1] * grad_up);
+        const double grad_up_component[3] = {grad_up.x, grad_up.y, grad_up.z};
+        const double grad_down_component[3] = {grad_down.x, grad_down.y, grad_down.z};
+        const double h_up_component[3] = {h_up.x, h_up.y, h_up.z};
+        const double h_down_component[3] = {h_down.x, h_down.y, h_down.z};
+        for (int l = 0; l < 3; ++l)
+        {
+            for (int m = 0; m <= l; ++m)
+            {
+                stress_gga[l * 3 + m]
+                    += h_up_component[l] * grad_up_component[m] + h_down_component[l] * grad_down_component[m];
+            }
+        }
+    }
+#endif
+
+    XC_Functional_Libxc::finish_func(funcs);
+}
+
 std::tuple<double, double, ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc( // Peize Lin update for nspin==4 at
                                                                                 // 2023.01.14
     const std::vector<int>& func_id,

source/source_hamilt/module_xc/test/test_xc_functional_ncgga_sf.cpp

@@ -1491,13 +1491,92 @@ TEST_F(RealPwNcgga, LibxcGgaGrad2DifferentiatesMixedLdaAndGgaComponents)
     expect_vtxc_matches_returned_potential("libxc_gga2_mixed_lda_gga", evaluate);
 }
 
+TEST_F(RealPwNcgga, LibxcGgaGrad2StressProductionDispatchClosesSmoothMetricDerivative)
+{
+    const BranchMargins margins = report_branch_margins("libxc_stress_smooth");
+    EXPECT_GT(margins.min_abs_total_density, 1.0);
+    EXPECT_GT(margins.min_signed_saturation_gap, 0.8);
+    EXPECT_GT(margins.min_eta_distance, 0.3);
 
+    // The analytic tensor must come through the existing public production
+    // dispatch.  Keeping the new helper out of the test-only patch lets the
+    // matched parent build and fail on the numerical identity, not at link
+    // time because the child-only helper does not exist yet.
+    expect_gradient_stress_metric_derivative(
+        "libxc_smooth_dispatch",
+        [this]() { return evaluate_libxc_gga(); },
+        [this]() { return evaluate_libxc_pbe_gradient_stress_dispatch(); },
+        3.0e-4,
+        1.0e-11);
+}
 
+TEST_F(RealPwNcgga, LibxcGgaGrad2StressProductionDispatchClosesLocalMapBranches)
+{
+    const Evaluator energy = [this]() { return evaluate_libxc_gga(); };
+    const StressEvaluator stress = [this]() { return evaluate_libxc_pbe_gradient_stress_dispatch(); };
 
+    set_negative_gga_state();
+    const BranchMargins negative = report_branch_margins("libxc_stress_negative_abs");
+    EXPECT_GT(negative.min_abs_total_density, 1.3);
+    EXPECT_GT(negative.min_signed_saturation_gap, 0.8);
+    expect_gradient_stress_metric_derivative("libxc_negative_abs_dispatch", energy, stress, 3.0e-4, 1.0e-11);
 
+    set_saturated_gga_state();
+    const BranchMargins saturated = report_branch_margins("libxc_stress_saturated");
+    EXPECT_LT(saturated.max_signed_saturation_gap, -0.1);
+    expect_gradient_stress_metric_derivative("libxc_saturated_dispatch", energy, stress, 3.0e-4, 1.0e-11);
 
+    set_inside_eta_state();
+    const BranchMargins radial = report_branch_margins("libxc_stress_inside_eta");
+    EXPECT_LT(radial.max_magnitude, 6.0e-4);
+    EXPECT_GT(radial.min_eta_distance, 4.0e-4);
+    expect_gradient_stress_metric_derivative("libxc_inside_eta_dispatch", energy, stress, 3.0e-4, 1.0e-11);
+}
 
+TEST_F(RealPwNcgga, LibxcGgaGrad2StressAggregatesPublicFunctionalScaling)
+{
+    const std::vector<int> gga = {XC_GGA_X_ITYH, XC_GGA_C_LYPR, XC_GGA_X_B88, XC_GGA_C_LYP};
+    const std::map<int, double> scaling
+        = {{XC_GGA_X_ITYH, -1.0}, {XC_GGA_C_LYPR, -1.0}, {XC_GGA_X_B88, 1.0}, {XC_GGA_C_LYP, 1.0}};
+    const std::vector<double> exchange_short = evaluate_gradient_stress_dispatch("GGA_X_ITYH");
+    const std::vector<double> correlation_short = evaluate_gradient_stress_dispatch("GGA_C_LYPR");
+    const std::vector<double> exchange_full = evaluate_gradient_stress_dispatch("GGA_X_B88");
+    const std::vector<double> correlation_full = evaluate_gradient_stress_dispatch("GGA_C_LYP");
+    const std::vector<double> scaled = evaluate_gradient_stress_dispatch("BLYP_LR");
+    ASSERT_EQ(exchange_short.size(), 9U);
+    ASSERT_EQ(correlation_short.size(), 9U);
+    ASSERT_EQ(exchange_full.size(), 9U);
+    ASSERT_EQ(correlation_full.size(), 9U);
+    ASSERT_EQ(scaled.size(), 9U);
+    for (int row = 0; row < 3; ++row)
+    {
+        for (int column = 0; column <= row; ++column)
+        {
+            const int index = row * 3 + column;
+            const double expected
+                = -exchange_short[index] - correlation_short[index] + exchange_full[index] + correlation_full[index];
+            EXPECT_NEAR(scaled[index], expected, 3.0e-12 * std::max(1.0, std::abs(expected)));
+        }
+    }
+    expect_gradient_stress_metric_derivative(
+        "libxc_blyp_lr_dispatch",
+        [this, &gga, &scaling]() { return evaluate_libxc(gga, &scaling); },
+        [this]() { return evaluate_gradient_stress_dispatch("BLYP_LR"); },
+        3.0e-4,
+        1.0e-11);
+}
 
+TEST_F(RealPwNcgga, LibxcGgaGrad2StressProductionDispatchClosesFullXcDiagonalWithoutCore)
+{
+    set_zero_core_density();
+    const BranchMargins margins = report_branch_margins("libxc_full_xc_no_core");
+    EXPECT_GT(margins.min_abs_total_density, 1.0);
+    EXPECT_GT(margins.min_signed_saturation_gap, 0.8);
+    expect_full_xc_diagonal_stress(
+        "libxc_full_xc_dispatch_no_core",
+        [this]() { return evaluate_libxc_gga(); },
+        [this]() { return evaluate_libxc_pbe_gradient_stress_dispatch(); });
+}
 #endif
 
 TEST_F(RealPwNcgga, BuiltinGgaGrad2VtxcEqualsFinalValencePotentialInnerProduct)
@@ -1653,13 +1732,64 @@ TEST_F(RealPwNcgga, BuiltinGgaGrad2IsCovariantUnderGlobalSpinRotation)
     }
 }
 
+TEST_F(RealPwNcgga, BuiltinGgaGrad2StressClosesSmoothSixComponentMetricDerivative)
+{
+    const BranchMargins margins = report_branch_margins("stress_smooth");
+    EXPECT_GT(margins.min_abs_total_density, 1.0);
+    EXPECT_GT(margins.min_signed_saturation_gap, 0.8);
+    EXPECT_GT(margins.min_eta_distance, 0.3);
+    expect_builtin_gradient_stress_metric_derivative("smooth");
 
+    set_zero_core_density();
+    const BranchMargins zero_core = report_branch_margins("stress_smooth_zero_core");
+    EXPECT_GT(zero_core.min_abs_total_density, 1.0);
+    EXPECT_GT(zero_core.min_signed_saturation_gap, 0.8);
+    expect_builtin_full_xc_diagonal_stress("smooth");
+}
 
+TEST_F(RealPwNcgga, BuiltinGgaGrad2StressClosesNegativeAbsSixComponentMetricDerivative)
+{
+    set_negative_gga_state();
+    const BranchMargins margins = report_branch_margins("stress_negative_abs");
+    EXPECT_GT(margins.min_abs_total_density, 1.3);
+    EXPECT_GT(margins.min_signed_saturation_gap, 0.8);
+    expect_builtin_gradient_stress_metric_derivative("negative_abs");
 
+    set_zero_core_density();
+    const BranchMargins zero_core = report_branch_margins("stress_negative_abs_zero_core");
+    EXPECT_GT(zero_core.min_abs_total_density, 1.3);
+    EXPECT_GT(zero_core.min_signed_saturation_gap, 0.8);
+    expect_builtin_full_xc_diagonal_stress("negative_abs");
+}
 
+TEST_F(RealPwNcgga, BuiltinGgaGrad2StressClosesSaturatedSixComponentMetricDerivative)
+{
+    set_saturated_gga_state();
+    const BranchMargins margins = report_branch_margins("stress_saturated");
+    EXPECT_LT(margins.max_signed_saturation_gap, -0.1);
+    expect_builtin_gradient_stress_metric_derivative("saturated");
 
+    set_zero_core_density();
+    const BranchMargins zero_core = report_branch_margins("stress_saturated_zero_core");
+    EXPECT_LT(zero_core.max_signed_saturation_gap, -0.1);
+    expect_builtin_full_xc_diagonal_stress("saturated");
+}
 
+TEST_F(RealPwNcgga, BuiltinGgaGrad2StressClosesRadialEtaSixComponentMetricDerivative)
+{
+    set_inside_eta_state();
+    const BranchMargins margins = report_branch_margins("stress_inside_eta");
+    EXPECT_GT(margins.min_abs_total_density, 0.025);
+    EXPECT_LT(margins.max_magnitude, 6.0e-4);
+    EXPECT_GT(margins.min_eta_distance, 4.0e-4);
+    expect_builtin_gradient_stress_metric_derivative("inside_eta");
 
+    set_zero_core_density();
+    const BranchMargins zero_core = report_branch_margins("stress_inside_eta_zero_core");
+    EXPECT_GT(zero_core.min_abs_total_density, 0.025);
+    EXPECT_LT(zero_core.max_magnitude, 6.0e-4);
+    expect_builtin_full_xc_diagonal_stress("inside_eta");
+}
 
 } // namespace