fix: handle zero-temperature rescaling

changes/user/bugfix.zero-temperature-rescaling.md

@@ -0,0 +1 @@
+- Zero-temperature rescaling remains finite and rejects positive targets that cannot be reached from zero kinetic energy.

src/resetKinetics/resetKinetics.cpp

@@ -27,6 +27,8 @@
 #include <cstddef>     // 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
@@ -127,7 +131,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(),

src/thermostat/berendsenThermostat.cpp

@@ -24,13 +24,15 @@
 
 #include <cmath>   // 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::BerendsenThermostat;
+using namespace customException;
 using namespace settings;
 using namespace simulationBox;
 using namespace physicalData;
@@ -69,11 +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))
-        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();
     const auto tempRatio = _targetTemperature / _temperature;

src/thermostat/velocityRescalingThermostat.cpp

@@ -24,15 +24,19 @@
 
 #include <cmath>   // 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
@@ -98,6 +102,17 @@ void VelocityRescalingThermostat::applyThermostat(
 
     _temperature = physicalData.getTemperature();
 
+    if (isZero(_temperature))
+    {
+        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());

tests/src/resetKinetics/testResetKinetics.cpp

@@ -26,6 +26,7 @@
 #include <memory>
 
 #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();

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
@@ -197,29 +198,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 ---------- */