From 0e0b2aaa4f8914f20cb52646f16fd7843257830c Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Tue, 28 Jul 2026 17:45:59 +0200 Subject: [PATCH 1/4] fix: handle zero-temperature rescaling --- src/resetKinetics/resetKinetics.cpp | 28 +++++++- src/thermostat/berendsenThermostat.cpp | 15 +++-- .../velocityRescalingThermostat.cpp | 18 +++++- tests/src/resetKinetics/testResetKinetics.cpp | 39 +++++++++++ tests/src/thermostat/testThermostat.cpp | 64 +++++++++++++++---- 5 files changed, 144 insertions(+), 20 deletions(-) diff --git a/src/resetKinetics/resetKinetics.cpp b/src/resetKinetics/resetKinetics.cpp index acc238269..6c609f85b 100644 --- a/src/resetKinetics/resetKinetics.cpp +++ b/src/resetKinetics/resetKinetics.cpp @@ -27,6 +27,8 @@ #include // for size_t #include "constants/conversionFactors.hpp" // for _FS_TO_S_, _S_TO_FS_ +#include "exceptions.hpp" // for UserInputException +#include "mathUtilities.hpp" // for isZero #include "physicalData.hpp" // for PhysicalData #include "simulationBox.hpp" // for SimulationBox #include "staticMatrix.hpp" // for operator*, operator+= @@ -38,7 +40,9 @@ using namespace linearAlgebra; using namespace physicalData; using namespace simulationBox; using namespace constants; +using namespace customException; using namespace settings; +using namespace utilities; /** * @brief Construct a new Reset Kinetics:: Reset Kinetics object @@ -129,7 +133,27 @@ void ResetKinetics::reset( void ResetKinetics::resetTemperature(SimulationBox &simBox) { const auto targetTemp = ThermostatSettings::getActualTargetTemperature(); - const auto lambda = ::sqrt(targetTemp / _temperature); + + if (isZero(targetTemp)) + { + std::ranges::for_each( + simBox.getAtoms(), + [](auto &atom) { atom->scaleVelocity(0.0); } + ); + + _temperature = simBox.calculateTemperature(); + _momentum = simBox.calculateMomentum(); + _angularMomentum = simBox.calculateAngularMomentum(_momentum); + return; + } + + if (isZero(_temperature)) + throw UserInputException( + "Cannot rescale a zero-temperature system to a positive target " + "temperature. Initialize velocities first." + ); + + const auto lambda = ::sqrt(targetTemp / _temperature); std::ranges::for_each( simBox.getAtoms(), @@ -323,4 +347,4 @@ size_t ResetKinetics::getFrequencyMomentumReset() const size_t ResetKinetics::getNStepsForcesReset() const { return _nStepsForcesReset; -} \ No newline at end of file +} diff --git a/src/thermostat/berendsenThermostat.cpp b/src/thermostat/berendsenThermostat.cpp index 5f2c51361..2377d8f0a 100644 --- a/src/thermostat/berendsenThermostat.cpp +++ b/src/thermostat/berendsenThermostat.cpp @@ -24,6 +24,7 @@ #include // for sqrt +#include "exceptions.hpp" // for UserInputException #include "mathUtilities.hpp" // for isZero #include "physicalData.hpp" // for PhysicalData #include "simulationBox.hpp" // for SimulationBox @@ -31,6 +32,7 @@ #include "timingsSettings.hpp" // for TimingsSettings using thermostat::BerendsenThermostat; +using namespace customException; using namespace settings; using namespace simulationBox; using namespace physicalData; @@ -69,13 +71,16 @@ void BerendsenThermostat::applyThermostat( _temperature = data.getTemperature(); - // If the kinetic energy is (approximately) zero, there is nothing to - // thermostat: dividing by _temperature would NaN all velocities - // (1 / 0 -> Inf, then vel * Inf = NaN when vel is 0). Skip silently. if (isZero(_temperature)) { stopTimingsSection("Berendsen"); - return; + if (isZero(_targetTemperature)) + return; + + throw UserInputException( + "Cannot apply Berendsen coupling to a zero-temperature system " + "with a positive target temperature. Initialize velocities first." + ); } const auto dt = TimingsSettings::getTimeStep(); @@ -113,4 +118,4 @@ void BerendsenThermostat::setTau(const double tau) { _tau = tau; } ThermostatType BerendsenThermostat::getThermostatType() const { return ThermostatType::BERENDSEN; -} \ No newline at end of file +} diff --git a/src/thermostat/velocityRescalingThermostat.cpp b/src/thermostat/velocityRescalingThermostat.cpp index c81ac9cb5..59500a7de 100644 --- a/src/thermostat/velocityRescalingThermostat.cpp +++ b/src/thermostat/velocityRescalingThermostat.cpp @@ -24,15 +24,19 @@ #include // for sqrt +#include "exceptions.hpp" // for UserInputException +#include "mathUtilities.hpp" // for isZero #include "physicalData.hpp" // for PhysicalData #include "simulationBox.hpp" // for SimulationBox #include "thermostatSettings.hpp" // for ThermostatType #include "timingsSettings.hpp" // for TimingsSettings using thermostat::VelocityRescalingThermostat; +using namespace customException; using namespace settings; using namespace simulationBox; using namespace physicalData; +using namespace utilities; /** * @brief Construct a new Velocity Rescaling Thermostat:: Velocity Rescaling @@ -80,6 +84,18 @@ void VelocityRescalingThermostat::applyThermostat( _temperature = physicalData.getTemperature(); + if (isZero(_temperature)) + { + stopTimingsSection("Velocity Rescaling"); + if (isZero(_targetTemperature)) + return; + + throw UserInputException( + "Cannot apply velocity rescaling to a zero-temperature system " + "with a positive target temperature. Initialize velocities first." + ); + } + const auto timeStep = TimingsSettings::getTimeStep(); const auto tempRatio = _targetTemperature / _temperature; const auto dof = double(simulationBox.getDegreesOfFreedom()); @@ -136,4 +152,4 @@ void VelocityRescalingThermostat::setTau(const double tau) { _tau = tau; } ThermostatType VelocityRescalingThermostat::getThermostatType() const { return ThermostatType::VELOCITY_RESCALING; -} \ No newline at end of file +} diff --git a/tests/src/resetKinetics/testResetKinetics.cpp b/tests/src/resetKinetics/testResetKinetics.cpp index 76b9d68dd..1458deed5 100644 --- a/tests/src/resetKinetics/testResetKinetics.cpp +++ b/tests/src/resetKinetics/testResetKinetics.cpp @@ -26,6 +26,7 @@ #include #include "atom.hpp" +#include "exceptions.hpp" #include "gtest/gtest.h" #include "molecule.hpp" #include "physicalData.hpp" @@ -124,6 +125,44 @@ TEST(TestResetKinetics, resetTemperatureRescalesVelocitiesAndStaysFinite) delete box; } +TEST(TestResetKinetics, resetTemperatureSupportsZeroKelvin) +{ + auto *box = makeBox(); + resetKinetics::ResetKinetics resetKinetics; + + for (const auto &atom : box->getAtoms()) atom->setVelocity({0.0, 0.0, 0.0}); + + settings::ThermostatSettings::setTargetTemperature(0.0); + resetKinetics.setTemperature(0.0); + resetKinetics.resetTemperature(*box); + + auto data = physicalData::PhysicalData(); + data.calculateTemperature(*box); + EXPECT_DOUBLE_EQ(data.getTemperature(), 0.0); + for (const auto &atom : box->getAtoms()) + EXPECT_EQ(atom->getVelocity(), linearAlgebra::Vec3D(0.0, 0.0, 0.0)); + + delete box; +} + +TEST(TestResetKinetics, rejectsPositiveTargetFromZeroTemperature) +{ + auto *box = makeBox(); + resetKinetics::ResetKinetics resetKinetics; + + for (const auto &atom : box->getAtoms()) atom->setVelocity({0.0, 0.0, 0.0}); + + settings::ThermostatSettings::setTargetTemperature(300.0); + resetKinetics.setTemperature(0.0); + + EXPECT_THROW( + resetKinetics.resetTemperature(*box), + customException::UserInputException + ); + + delete box; +} + TEST(TestResetKinetics, resetMomentumZerosTotalLinearMomentum) { auto *box = makeBox(); diff --git a/tests/src/thermostat/testThermostat.cpp b/tests/src/thermostat/testThermostat.cpp index ab0303b7d..e1135e77d 100644 --- a/tests/src/thermostat/testThermostat.cpp +++ b/tests/src/thermostat/testThermostat.cpp @@ -26,6 +26,7 @@ #include "berendsenThermostat.hpp" // for BerendsenThermostat #include "constants/internalConversionFactors.hpp" // for _TEMPERATURE_FACTOR_ +#include "exceptions.hpp" // for UserInputException #include "gtest/gtest.h" // for InitGoogleTest #include "langevinThermostat.hpp" // for LangevinThermostat #include "noseHooverThermostat.hpp" // for NoseHooverThermostat @@ -195,29 +196,68 @@ TEST_F(TestThermostat, velocityRescaling_applyDoesNotNaN) } } -// Regression test: starting from zero kinetic energy (T == 0) used to -// produce NaN velocities, because tempRatio = T_target / 0 = Inf and -// the velocity scaling 0 * Inf = NaN. The guard skips the scaling and -// leaves velocities at zero. -TEST_F(TestThermostat, applyBerendsen_zeroTemperatureNoNaN) +TEST_F(TestThermostat, berendsenZeroTemperatureDoesNotNaN) +{ + delete _thermostat; + _thermostat = new thermostat::BerendsenThermostat(0.0, 100.0); + settings::TimingsSettings::setTimeStep(0.1); + + for (auto &atom : _simulationBox->getAtoms()) + atom->setVelocity({0.0, 0.0, 0.0}); + + _thermostat->applyThermostat(*_simulationBox, *_data); + + EXPECT_TRUE(std::isfinite(_data->getTemperature())); + for (const auto &atom : _simulationBox->getAtoms()) + for (size_t dimension = 0; dimension < 3; ++dimension) + EXPECT_TRUE(std::isfinite(atom->getVelocity()[dimension])); +} + +TEST_F(TestThermostat, berendsenRejectsPositiveTargetFromZero) { delete _thermostat; _thermostat = new thermostat::BerendsenThermostat(300.0, 100.0); settings::TimingsSettings::setTimeStep(0.1); + for (auto &atom : _simulationBox->getAtoms()) + atom->setVelocity({0.0, 0.0, 0.0}); + + EXPECT_THROW( + _thermostat->applyThermostat(*_simulationBox, *_data), + customException::UserInputException + ); +} + +TEST_F(TestThermostat, velocityRescalingZeroTemperatureDoesNotNaN) +{ + delete _thermostat; + _thermostat = new thermostat::VelocityRescalingThermostat(0.0, 100.0); + settings::TimingsSettings::setTimeStep(0.1); + for (auto &atom : _simulationBox->getAtoms()) atom->setVelocity({0.0, 0.0, 0.0}); _thermostat->applyThermostat(*_simulationBox, *_data); - EXPECT_FALSE(std::isnan(_data->getTemperature())); - EXPECT_FALSE(std::isinf(_data->getTemperature())); + EXPECT_TRUE(std::isfinite(_data->getTemperature())); for (const auto &atom : _simulationBox->getAtoms()) - for (size_t i = 0; i < 3; ++i) - { - EXPECT_FALSE(std::isnan(atom->getVelocity()[i])); - EXPECT_FALSE(std::isinf(atom->getVelocity()[i])); - } + for (size_t dimension = 0; dimension < 3; ++dimension) + EXPECT_TRUE(std::isfinite(atom->getVelocity()[dimension])); +} + +TEST_F(TestThermostat, velocityRescalingRejectsPositiveTargetFromZero) +{ + delete _thermostat; + _thermostat = new thermostat::VelocityRescalingThermostat(300.0, 100.0); + settings::TimingsSettings::setTimeStep(0.1); + + for (auto &atom : _simulationBox->getAtoms()) + atom->setVelocity({0.0, 0.0, 0.0}); + + EXPECT_THROW( + _thermostat->applyThermostat(*_simulationBox, *_data), + customException::UserInputException + ); } /* ---------- LangevinThermostat ---------- */ From 43b253c1cda9887dc6072204952d51bf7a99f6f9 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Thu, 30 Jul 2026 12:59:47 +0200 Subject: [PATCH 2/4] chore: add required changelog fragment --- changes/zero-temperature-rescaling.user.bugfix.md | 1 + 1 file changed, 1 insertion(+) create mode 100644 changes/zero-temperature-rescaling.user.bugfix.md diff --git a/changes/zero-temperature-rescaling.user.bugfix.md b/changes/zero-temperature-rescaling.user.bugfix.md new file mode 100644 index 000000000..228975b75 --- /dev/null +++ b/changes/zero-temperature-rescaling.user.bugfix.md @@ -0,0 +1 @@ +- Zero-temperature rescaling remains finite and rejects positive targets that cannot be reached from zero kinetic energy. From ce2db3f040208b69b5010970729b54665184cfc6 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Sat, 1 Aug 2026 22:57:06 +0200 Subject: [PATCH 3/4] chore: migrate zero-temperature changelog fragment --- .../bugfix.zero-temperature-rescaling.md} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename changes/{zero-temperature-rescaling.user.bugfix.md => user/bugfix.zero-temperature-rescaling.md} (100%) diff --git a/changes/zero-temperature-rescaling.user.bugfix.md b/changes/user/bugfix.zero-temperature-rescaling.md similarity index 100% rename from changes/zero-temperature-rescaling.user.bugfix.md rename to changes/user/bugfix.zero-temperature-rescaling.md From cc9ed25330ebb28c243cba1f823bc34d9dcfd713 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Tue, 4 Aug 2026 10:03:59 +0200 Subject: [PATCH 4/4] style: format thermostat test includes --- tests/src/thermostat/testThermostat.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/tests/src/thermostat/testThermostat.cpp b/tests/src/thermostat/testThermostat.cpp index 1559ebc28..ab65908a8 100644 --- a/tests/src/thermostat/testThermostat.cpp +++ b/tests/src/thermostat/testThermostat.cpp @@ -24,17 +24,17 @@ #include // for sqrt -#include "berendsenThermostat.hpp" // for BerendsenThermostat -#include "constants/internalConversionFactors.hpp" // for _TEMPERATURE_FACTOR_ -#include "exceptions.hpp" // for UserInputException -#include "gtest/gtest.h" // for InitGoogleTest -#include "langevinThermostat.hpp" // for LangevinThermostat -#include "noseHooverThermostat.hpp" // for NoseHooverThermostat -#include "physicalData.hpp" // for PhysicalData -#include "simulationBox.hpp" // for SimulationBox -#include "thermostatSettings.hpp" // for ThermostatType -#include "timingsSettings.hpp" // for TimingsSettings -#include "velocityRescalingThermostat.hpp" // for VelocityRescalingThermostat +#include "berendsenThermostat.hpp" // for BerendsenThermostat +#include "constants/internalConversionFactors.hpp" // for _TEMPERATURE_FACTOR_ +#include "exceptions.hpp" // for UserInputException +#include "gtest/gtest.h" // for InitGoogleTest +#include "langevinThermostat.hpp" // for LangevinThermostat +#include "noseHooverThermostat.hpp" // for NoseHooverThermostat +#include "physicalData.hpp" // for PhysicalData +#include "simulationBox.hpp" // for SimulationBox +#include "thermostatSettings.hpp" // for ThermostatType +#include "timingsSettings.hpp" // for TimingsSettings +#include "velocityRescalingThermostat.hpp" // for VelocityRescalingThermostat TEST_F(TestThermostat, calculateTemperature) {