Use mapped-density FFT gradients for noncollinear gga_grad=2

source/source_hamilt/module_xc/libxc_abacus.h

@@ -10,7 +10,6 @@
 #include <xc.h>
 #include <xc_funcs.h>
 
-#include <array>
 #include <tuple>
 #include <vector>
 
@@ -30,15 +29,14 @@ namespace XC_Functional_Libxc
 
     // Complete forward data for the gga_grad=2 noncollinear Libxc graph:
     //   rho_s = N_s(x),
-    //   g_s   = sum_A (d N_s / d x_A) G_h x_A.
-    // Keeping the local map and all input gradients together lets the reverse
-    // use the exact same branch choices and radial Hessian as the forward.
+    //   g_s   = G_h N_s(x).
+    // Reverse the FFT gradient before applying the local map Jacobian;
+    // this uses the same clipping branches as the forward energy.
     struct NclSfDiscreteData
     {
         std::vector<ModuleXC::NcggaSpinMapPoint> spin_map;
         std::vector<double> rho;
         std::vector<std::vector<ModuleBase::Vector3<double>>> spin_gradient;
-        std::array<std::vector<ModuleBase::Vector3<double>>, 3> grad_m;
     };
 
 //-------------------
@@ -133,7 +131,7 @@ namespace XC_Functional_Libxc
         const Charge* const chr);
 
     // Build the exact gga_grad=2 local spin map and, when requested, its
-    // projected FFT-gradient graph.  LDA-only callers set need_gradient=false.
+    // mapped-density FFT-gradient graph.  LDA-only callers set need_gradient=false.
     extern NclSfDiscreteData make_ncl_sf_discrete_data(
         const std::size_t nrxx,
         const double tpiba,

source/source_hamilt/module_xc/libxc_pot.cpp

@@ -40,7 +40,7 @@ void XC_Functional_Libxc::gradcorr_ncgga_sf_libxc(const std::vector<int>& func_i
     }
 
     // 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
+    // FFT gradients of the mapped densities, 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);
@@ -105,10 +105,11 @@ void XC_Functional_Libxc::gradcorr_ncgga_sf_libxc(const std::vector<int>& func_i
         }
     }
 
-// 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.
+    // For g_s = G_h N_s(x), a reciprocal-metric variation at fixed real-grid
+    // density values changes G_h while holding N_s(x) fixed. The gradient-metric
+    // contribution is sum_s h_s,l g_s,m, with h_s = df/dg_s built from the
+    // sanitizer-reversed, component-scaled aggregate above. This is only the
+    // gradient-metric contribution, not the full cell-strain derivative.
 #ifdef _OPENMP
 #pragma omp parallel
     {
@@ -200,7 +201,7 @@ std::tuple<double, double, ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(
     const int nspin = (nspin_in == 1 || (nspin_in == 4 && !domag && !domag_z)) ? 1 : 2;
 
     // For nspin=4 with noncollinear magnetism, gga_grad=2 selects the
-    // regularized projected local-collinear graph; gga_grad=0/1 keeps the
+    // regularized mapped-density local-collinear graph; gga_grad=0/1 keeps the
     // original collinear algorithm.
     const bool has_mag = domag || domag_z;
     const bool use_lca = (nspin_in == 4) && has_mag && gga_grad == 2;
@@ -246,7 +247,7 @@ std::tuple<double, double, ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(
     else if (use_lca)
     {
         // gga_grad=2 uses one complete local map for both the Libxc density
-        // input and the projected FFT-gradient graph.  LDA-only functionals
+        // input and the mapped-density FFT-gradient graph. LDA-only functionals
         // need the same local map but do not pay for gradients.
         sf_data = XC_Functional_Libxc::make_ncl_sf_discrete_data(nrxx, tpiba, chr, is_gga);
         rho = sf_data.rho;
@@ -372,7 +373,7 @@ std::tuple<double, double, ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(
 
         // Keep the established energy accumulation and reduction order.  In
         // gga_grad=2, reverse every sanitizer now, apply the component scaling,
-        // and aggregate before traversing the shared projected graph once.
+        // and aggregate before reversing the shared mapped-density graph once.
         etxc += XC_Functional_Libxc::convert_etxc(nspin, nrxx, energy_mask, rho, exc) * factor;
         if (use_lca)
         {

source/source_hamilt/module_xc/libxc_tools.cpp

@@ -87,50 +87,9 @@ XC_Functional_Libxc::make_ncl_sf_discrete_data(
 		return data;
 	}
 
-	std::vector<ModuleBase::Vector3<double>> grad_total(nrxx);
-	std::vector<double> real_field(nrxx);
-	std::vector<std::complex<double>> reciprocal(chr->rhopw->npw);
-	#ifdef _OPENMP
-	#pragma omp parallel for schedule(static, 1024)
-	#endif
-	for (std::size_t ir = 0; ir < nrxx; ++ir)
-	{
-		real_field[ir] = chr->rho[0][ir] + chr->rho_core[ir];
-	}
-	chr->rhopw->real2recip(real_field.data(), reciprocal.data());
-	XC_Functional::grad_rho(
-		reciprocal.data(), grad_total.data(), chr->rhopw, tpiba);
-
-	for (int mu = 0; mu < 3; ++mu)
-	{
-		data.grad_m[mu].resize(nrxx);
-		chr->rhopw->real2recip(chr->rho[mu + 1], reciprocal.data());
-		XC_Functional::grad_rho(
-			reciprocal.data(), data.grad_m[mu].data(), chr->rhopw, tpiba);
-	}
-
-	data.spin_gradient.resize(nspin);
-	for (int spin = 0; spin < nspin; ++spin)
-	{
-		data.spin_gradient[spin].resize(nrxx);
-	}
-	#ifdef _OPENMP
-	#pragma omp parallel for schedule(static, 512)
-	#endif
-	for (std::size_t ir = 0; ir < nrxx; ++ir)
-	{
-		for (int spin = 0; spin < nspin; ++spin)
-		{
-			ModuleBase::Vector3<double> gradient
-				= data.spin_map[ir].jacobian(spin, 0) * grad_total[ir];
-			for (int mu = 0; mu < 3; ++mu)
-			{
-				gradient += data.spin_map[ir].jacobian(spin, mu + 1)
-						* data.grad_m[mu][ir];
-			}
-			data.spin_gradient[spin][ir] = gradient;
-		}
-	}
+	// Differentiate the mapped spin densities with the actual FFT operator.
+	// The reverse must apply its adjoint before the local map Jacobian.
+	data.spin_gradient = cal_gdr(nspin, nrxx, data.rho, tpiba, chr);
 	return data;
 }
 
@@ -176,12 +135,9 @@ ModuleBase::matrix XC_Functional_Libxc::reverse_ncl_sf_discrete(
 	{
 		assert(data.spin_gradient[spin].size() == nrxx);
 	}
-	for (int mu = 0; mu < 3; ++mu)
-	{
-		assert(data.grad_m[mu].size() == nrxx);
-	}
 
-	std::vector<ModuleBase::Vector3<double>> h_up(nrxx), h_down(nrxx);
+	std::vector<ModuleBase::Vector3<double>> h_up(nrxx);
+	std::vector<ModuleBase::Vector3<double>> h_down(nrxx);
 	#ifdef _OPENMP
 	#pragma omp parallel for schedule(static, 512)
 	#endif
@@ -198,43 +154,20 @@ ModuleBase::matrix XC_Functional_Libxc::reverse_ncl_sf_discrete(
 				 + dsigma[sigma_index + 1] * data.spin_gradient[0][ir]);
 	}
 
-	std::vector<ModuleBase::Vector3<double>> flux(nrxx);
-	std::vector<double> divergence(nrxx);
-	for (int channel = 0; channel < nchannel; ++channel)
+	std::vector<double> divergence_up(nrxx);
+	std::vector<double> divergence_down(nrxx);
+	XC_Functional::grad_dot(h_up.data(), divergence_up.data(), chr->rhopw, tpiba);
+	XC_Functional::grad_dot(h_down.data(), divergence_down.data(), chr->rhopw, tpiba);
+	#ifdef _OPENMP
+	#pragma omp parallel for schedule(static, 512)
+	#endif
+	for (std::size_t ir = 0; ir < nrxx; ++ir)
 	{
-		#ifdef _OPENMP
-		#pragma omp parallel for schedule(static, 512)
-		#endif
-		for (std::size_t ir = 0; ir < nrxx; ++ir)
-		{
-			flux[ir]
-				= data.spin_map[ir].jacobian(0, channel) * h_up[ir]
-				  + data.spin_map[ir].jacobian(1, channel) * h_down[ir];
-		}
-		XC_Functional::grad_dot(
-			flux.data(), divergence.data(), chr->rhopw, tpiba);
-
-		#ifdef _OPENMP
-		#pragma omp parallel for schedule(static, 512)
-		#endif
-		for (std::size_t ir = 0; ir < nrxx; ++ir)
+		for (int channel = 0; channel < nchannel; ++channel)
 		{
-			potential(channel, ir) -= divergence[ir];
-			if (channel == 0 || data.spin_map[ir].saturated)
-			{
-				continue;
-			}
-
-			const ModuleBase::Vector3<double> spin_flux
-				= 0.5 * (h_up[ir] - h_down[ir]);
-			double local_response = 0.0;
-			for (int nu = 0; nu < 3; ++nu)
-			{
-				local_response
-					+= data.spin_map[ir].radial.jacobian(nu, channel - 1)
-					   * (spin_flux * data.grad_m[nu][ir]);
-			}
-			potential(channel, ir) += local_response;
+			potential(channel, ir)
+				-= data.spin_map[ir].jacobian(0, channel) * divergence_up[ir]
+				   + data.spin_map[ir].jacobian(1, channel) * divergence_down[ir];
 		}
 	}
 	return potential;

source/source_hamilt/module_xc/test/test_xc_functional_ncgga_sf.cpp

@@ -274,6 +274,75 @@ TEST_P(RealPwNcgga, ZeroMagnetizationIsFinite)
             EXPECT_NEAR(std::get<2>(result)(c, ir), 0.0, 1e-13);
 }
 
+// A local rotation preserves both mapped spin densities pointwise. This is
+// stronger than global spin covariance and fails for J(x) G_h(x) at finite cutoff.
+TEST_P(RealPwNcgga, LocalSpinRotationPreservesEnergyAndCovariesPotential)
+{
+    const VxcResult reference = evaluate();
+    std::vector<double> cosine(pw.nrxx);
+    std::vector<double> sine(pw.nrxx);
+    for (int ir = 0; ir < pw.nrxx; ++ir)
+    {
+        const double angle = 0.7 * std::sin(0.31 * ir);
+        cosine[ir] = std::cos(angle);
+        sine[ir] = std::sin(angle);
+        const double mx = density[1][ir];
+        const double my = density[2][ir];
+        density[1][ir] = cosine[ir] * mx - sine[ir] * my;
+        density[2][ir] = sine[ir] * mx + cosine[ir] * my;
+    }
+    const VxcResult rotated = evaluate();
+    EXPECT_NEAR(std::get<0>(reference), std::get<0>(rotated), 1e-10);
+    const ModuleBase::matrix& before = std::get<2>(reference);
+    const ModuleBase::matrix& after = std::get<2>(rotated);
+    for (int ir = 0; ir < pw.nrxx; ++ir)
+    {
+        const double expected_x = cosine[ir] * before(1, ir) - sine[ir] * before(2, ir);
+        const double expected_y = sine[ir] * before(1, ir) + cosine[ir] * before(2, ir);
+        EXPECT_NEAR(after(0, ir), before(0, ir), 1e-11);
+        EXPECT_NEAR(after(1, ir), expected_x, 1e-11);
+        EXPECT_NEAR(after(2, ir), expected_y, 1e-11);
+        EXPECT_NEAR(after(3, ir), before(3, ir), 1e-11);
+    }
+}
+
+TEST_P(RealPwNcgga, ClippingBoundaryEnergyIsContinuous)
+{
+    // Each point reaches |m|=|n+core| together, while the spin direction varies.
+    const std::array<std::vector<double>, 4> original = density;
+    double previous_jump = 0.0;
+    for (int level = 0; level < 4; ++level)
+    {
+        const double epsilon = 1e-5 / (1 << level);
+        double energies[2];
+        for (int side = 0; side < 2; ++side)
+        {
+            const double factor = side == 0 ? 1.0 - epsilon : 1.0 + epsilon;
+            for (int ir = 0; ir < pw.nrxx; ++ir)
+            {
+                const double mx = original[1][ir];
+                const double my = original[2][ir];
+                const double mz = original[3][ir];
+                const double magnitude = std::sqrt(mx * mx + my * my + mz * mz);
+                const double target = factor * std::abs(density[0][ir] + core[ir]);
+                for (int c = 1; c < noncollinear_spin; ++c)
+                {
+                    density[c][ir] = original[c][ir] * target / magnitude;
+                }
+            }
+            const VxcResult result = evaluate();
+            energies[side] = std::get<0>(result);
+        }
+        const double jump = std::abs(energies[0] - energies[1]);
+        if (level > 0)
+        {
+            EXPECT_LE(jump, 0.6 * previous_jump + 1e-10);
+        }
+        previous_jump = jump;
+    }
+}
+
+
 INSTANTIATE_TEST_SUITE_P(ProductionBranches, RealPwNcgga,
                         testing::Combine(testing::Bool(), testing::Values(Smooth, Negative, Saturated, Radial)));
 } // namespace