diff --git a/source/source_pw/module_ofdft/kedf_tf.cpp b/source/source_pw/module_ofdft/kedf_tf.cpp index c5bf62173a5..7c78380733c 100644 --- a/source/source_pw/module_ofdft/kedf_tf.cpp +++ b/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; } -} \ No newline at end of file +} diff --git a/source/source_pw/module_ofdft/test/test_kedf_manager.cpp b/source/source_pw/module_ofdft/test/test_kedf_manager.cpp index c04ce916246..3ae65cde414 100644 --- a/source/source_pw/module_ofdft/test/test_kedf_manager.cpp +++ b/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 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 perturbed_rho(basis.nrxx); + std::vector 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); + } + } +} diff --git a/source/source_pw/module_ofdft/test/test_kedf_tf.cpp b/source/source_pw/module_ofdft/test/test_kedf_tf.cpp index b1b6ab6e0e6..e307ade96e3 100644 --- a/source/source_pw/module_ofdft/test/test_kedf_tf.cpp +++ b/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); + } + } +} diff --git a/tests/07_OFDFT/11_OF_TF_weight/README b/tests/07_OFDFT/11_OF_TF_weight/README index 213b0f1553c..d36bd1775b9 100644 --- a/tests/07_OFDFT/11_OF_TF_weight/README +++ b/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. diff --git a/tests/07_OFDFT/11_OF_TF_weight/result.ref b/tests/07_OFDFT/11_OF_TF_weight/result.ref index c712de7d05a..2a29d3a9608 100644 --- a/tests/07_OFDFT/11_OF_TF_weight/result.ref +++ b/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