fix(ofdft): apply TF weight consistently to single-spin energy

source/source_pw/module_ofdft/kedf_tf.cpp

@@ -14,7 +14,7 @@ void KEDF_TF::set_para(int nx, double dV, double tf_weight)
 
 /**
  * @brief Get the energy of TF KEDF,
- * \f[ E_{TF} = c_{TF} * \int{\rho^{5/3} dr} \f]
+ * \f[ E_{TF} = tf_weight * c_{TF} * \int{\rho^{5/3} dr} \f]
  *
  * @param prho charge density
  * @return the energy of TF KEDF
@@ -28,7 +28,7 @@ double KEDF_TF::get_energy(const double* const* prho)
         {
             energy += std::pow(prho[0][ir], 5. / 3.);
         }
-        energy *= this->dV_ * this->c_tf_;
+        energy *= this->dV_ * this->c_tf_ * this->tf_weight_;
     }
     else if (PARAM.inp.nspin == 2)
     {
@@ -139,4 +139,4 @@ void KEDF_TF::get_stress(double cell_vol)
     {
         this->stress(i, i) = temp;
     }
-}
+}

source/source_pw/module_ofdft/test/test_kedf_manager.cpp

@@ -201,3 +201,63 @@ TEST_F(TestKEDFManager, InitWTAndReinitializeWithChangedParameters)
     inp.of_vw_weight = 0.8;
     check_wt(inp);
 }
+
+TEST_F(TestKEDFManager, WeightedEnergyMatchesConservedDensityDerivative)
+{
+    const char* kinetics[] = {"tf", "wt"};
+    const double weights[] = {1.0, 1.1};
+    const double step = 1e-4;
+    std::vector<double> direction(basis.nrxx);
+    double particle_change = 0.0;
+    for (int ir = 0; ir < basis.nrxx; ++ir)
+    {
+        const int iz = ir % basis.nplane + basis.startz_current;
+        const double phase = ModuleBase::TWO_PI * iz / basis.nz;
+        direction[ir] = std::cos(phase);
+        particle_change += direction[ir] * dV;
+    }
+    ASSERT_NEAR(particle_change, 0.0, 1e-10);
+
+    for (const auto kinetic : kinetics)
+    {
+        SCOPED_TRACE(kinetic);
+        for (const double weight : weights)
+        {
+            SCOPED_TRACE(weight);
+            Input_para inp;
+            inp.of_kinetic = kinetic;
+            inp.of_tf_weight = weight;
+            manager.init(inp, &basis, dV, nelec);
+            const double* density[] = {rho.data()};
+            const double* wave[] = {phi.data()};
+            ModuleBase::matrix gradient(1, basis.nrxx);
+            manager.get_potential(density, wave, &basis, gradient);
+            double product = 0.0;
+            for (int ir = 0; ir < basis.nrxx; ++ir)
+            {
+                // The manager returns dT/dphi; convert to dT/drho.
+                product += gradient(0, ir) * direction[ir] / (2.0 * phi[ir]) * dV;
+            }
+
+            double energies[2];
+            std::vector<double> perturbed_rho(basis.nrxx);
+            std::vector<double> perturbed_phi(basis.nrxx);
+            for (int sign_index = 0; sign_index < 2; ++sign_index)
+            {
+                const double sign = 1.0 - 2.0 * sign_index;
+                for (int ir = 0; ir < basis.nrxx; ++ir)
+                {
+                    perturbed_rho[ir] = rho[ir] + sign * step * direction[ir];
+                    perturbed_phi[ir] = std::sqrt(perturbed_rho[ir]);
+                }
+                const double* perturbed_density[] = {perturbed_rho.data()};
+                const double* perturbed_wave[] = {perturbed_phi.data()};
+                ModuleBase::matrix perturbed_gradient(1, basis.nrxx);
+                manager.get_potential(perturbed_density, perturbed_wave, &basis, perturbed_gradient);
+                energies[sign_index] = manager.get_energy();
+            }
+            const double derivative = (energies[0] - energies[1]) / (2.0 * step);
+            EXPECT_NEAR(derivative, product, 1e-5);
+        }
+    }
+}

source/source_pw/module_ofdft/test/test_kedf_tf.cpp

@@ -77,3 +77,73 @@ TEST_F(TestKEDFTF, StressIsIsotropicAndScalesWithVolume)
         EXPECT_NEAR(kedf.stress(i, i), expected * 0.5, 1e-12);
     }
 }
+
+TEST_F(TestKEDFTF, WeightedEnergyMatchesPointwiseIntegralAndStress)
+{
+    const double density[] = {1.0, 8.0};
+    const double* rho[] = {density};
+    const double coefficient = 5.742468000376382;
+    const double dV = 0.25;
+    const double volume = 2.0 * dV;
+    const double weights[] = {0.0, 0.4, 1.0, 1.1};
+    for (const double weight : weights)
+    {
+        SCOPED_TRACE(weight);
+        kedf.set_para(2, dV, weight);
+        const double expected = weight * coefficient * 33.0 * dV;
+        const double energy = kedf.get_energy(rho);
+        EXPECT_NEAR(energy, expected, 1e-12);
+
+        double integral = 0.0;
+        for (int ir = 0; ir < 2; ++ir)
+        {
+            integral += kedf.get_energy_density(rho, 0, ir) * dV;
+        }
+        EXPECT_NEAR(energy, integral, 1e-12);
+
+        ModuleBase::matrix potential(1, 2);
+        kedf.tf_potential(rho, potential);
+        EXPECT_NEAR(kedf.tf_energy, expected, 1e-12);
+        kedf.get_stress(volume);
+        const double expected_stress = 2.0 * expected / (3.0 * volume);
+        for (int i = 0; i < 3; ++i)
+        {
+            EXPECT_NEAR(kedf.stress(i, i), expected_stress, 1e-12);
+        }
+    }
+}
+
+TEST_F(TestKEDFTF, WeightedPotentialMatchesConservedDensityDerivative)
+{
+    const double density[] = {1.0, 8.0};
+    const double direction[] = {1.0, -1.0};
+    const double* rho[] = {density};
+    const double dV = 0.25;
+    const double weights[] = {0.0, 0.4, 1.0, 1.1};
+    const double steps[] = {1e-3, 1e-4};
+    for (const double weight : weights)
+    {
+        SCOPED_TRACE(weight);
+        kedf.set_para(2, dV, weight);
+        ModuleBase::matrix potential(1, 2);
+        kedf.tf_potential(rho, potential);
+        const double product = dV * (potential(0, 0) - potential(0, 1));
+        for (const double step : steps)
+        {
+            SCOPED_TRACE(step);
+            double plus_density[2];
+            double minus_density[2];
+            for (int ir = 0; ir < 2; ++ir)
+            {
+                plus_density[ir] = density[ir] + step * direction[ir];
+                minus_density[ir] = density[ir] - step * direction[ir];
+            }
+            const double* plus_rho[] = {plus_density};
+            const double* minus_rho[] = {minus_density};
+            const double plus_energy = kedf.get_energy(plus_rho);
+            const double minus_energy = kedf.get_energy(minus_rho);
+            const double derivative = (plus_energy - minus_energy) / (2.0 * step);
+            EXPECT_NEAR(derivative, product, 1e-6);
+        }
+    }
+}

tests/07_OFDFT/11_OF_TF_weight/README

@@ -1 +1,9 @@
-Test the energy, force, and stress of `WT + 1.1 * TF + vW`, with of_tf_weight=1.1, symmetry=on
+Test WT with of_tf_weight=1.1 and symmetry enabled, including energy, force, and stress.
+
+The nonlocal WT kernel uses the same TF/vW weights as the local contributions.
+The reference includes the correction for issue #8243: the single-spin TF energy
+must include of_tf_weight, consistently with the potential. The earlier reference
+used an unweighted TF energy and consequently an inconsistent cached TF stress.
+With the corrected energy, this Al case converges to -57.4900280759 eV and
+40.538472 kbar pressure. The stress collector sums absolute tensor components,
+so totalstressref is 121.615416 kbar. Forces remain zero by symmetry.

tests/07_OFDFT/11_OF_TF_weight/result.ref

@@ -1,7 +1,7 @@
-etotref -59.7090017316024841
-etotperatomref -59.7090017316
+etotref -57.4900280758884890
+etotperatomref -57.4900280759
 totalforceref 0.000000
-totalstressref 332.896200
+totalstressref 121.615416
 pointgroupref O_h
 spacegroupref O_h
 nksibzref 1