Handle simulation edge cases safely

changes/simulation-edge-cases.user.bugfix.md

@@ -0,0 +1 @@
+- PQ now rejects invalid temperature ramps and handles zero-temperature, cell-list, and single-step simulation edge cases safely.

docs/sphinx/src/userGuide/inputFile.rst

@@ -815,6 +815,8 @@ Temperature Ramp Frequency
 
 With the ``temp_ramp_frequency`` keyword the user can specify the frequency of the temperature ramping from the ``start_temp`` to the ``temp`` value. If no starting temperature is given the keyword will be ignored. If a starting temperature is given and this keyword is omitted the temperature ramping will be performed, so that each step the temperature is increased by the same value.
 
+If the ramp length is not divisible by this frequency, the temperature increments are scaled by the number of scheduled updates so that the final update reaches the requested target temperature exactly.
+
 .. centered:: *default value* = 1 step
 
 .. _thermostatKey:

external/progressbar/include/progressbar.hpp

@@ -145,7 +145,7 @@ inline void progressbar::update() {
     int perc = 0;
 
     // compute percentage, if did not change, do nothing and return
-    perc = progress*100./(n_cycles-1);
+    perc = n_cycles == 1 ? 100 : progress*100./(n_cycles-1);
     if (perc < last_perc) return;
 
     // update percentage each unit

include/input/parameterFileReader/parameterFileReader.hpp

@@ -52,6 +52,7 @@ namespace input::parameterFile
 
        public:
         ParameterFileReader(const std::string &filename, pq::Engine &engine);
+        ~ParameterFileReader();
 
         void read();
         void deleteSection(const pq::ParamFileSection *section);
@@ -72,4 +73,4 @@ namespace input::parameterFile
 
 }   // namespace input::parameterFile
 
-#endif   // _PARAMETER_FILE_READER_HPP_
+#endif   // _PARAMETER_FILE_READER_HPP_

src/input/parameterFileReader/parameterFileReader.cpp

@@ -68,6 +68,8 @@ ParameterFileReader::ParameterFileReader(
     _parameterFileSections.push_back(make_unique<NonCoulombicsSection>());
 }
 
+ParameterFileReader::~ParameterFileReader() = default;
+
 /**
  * @brief determines which section of the parameter file the header line belongs
  * to
@@ -219,4 +221,4 @@ std::vector<std::unique_ptr<ParameterFileSection>> &ParameterFileReader::
 const std::string &ParameterFileReader::getFilename() const
 {
     return _fileName;
-}
+}

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
@@ -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;
-}
+}

src/setup/thermostatSetup.cpp

@@ -225,36 +225,47 @@ void ThermostatSetup::setupTemperatureRamp()
         return;
 
     /*************************************************************
-     * resetting the target temperature to the start temperature *
+     * If steps is 0, set the steps to the total number of steps *
      *************************************************************/
 
-    const auto startTemp = ThermostatSettings::getStartTemperature();
+    auto       steps = ThermostatSettings::getTemperatureRampSteps();
+    const auto useFullSimulation = steps == 0;
 
-    _engine.getThermostat().setTargetTemperature(startTemp);
-    ThermostatSettings::setActualTargetTemperature(startTemp);
+    if (useFullSimulation)
+        steps = TimingsSettings::getNumberOfSteps();
+
+    if (steps == 0)
+        throw InputFileException(
+            "Temperature ramp requires at least one simulation step"
+        );
+
+    const auto frequency = ThermostatSettings::getTemperatureRampFrequency();
 
-    auto steps = ThermostatSettings::getTemperatureRampSteps();
+    if (frequency == 0)
+        throw InputFileException(
+            "Temperature ramp frequency must be greater than zero"
+        );
+
+    if (useFullSimulation)
+        ThermostatSettings::setTemperatureRampSteps(steps);
 
     /*************************************************************
-     * If steps is 0, set the steps to the total number of steps *
+     * resetting the target temperature to the start temperature *
      *************************************************************/
 
-    if (steps == 0)
-    {
-        steps = TimingsSettings::getNumberOfSteps();
-        ThermostatSettings::setTemperatureRampSteps(steps);
-    }
+    const auto startTemp = ThermostatSettings::getStartTemperature();
 
+    _engine.getThermostat().setTargetTemperature(startTemp);
+    ThermostatSettings::setActualTargetTemperature(startTemp);
     _engine.getThermostat().setTemperatureRampingSteps(steps);
-
-    const auto frequency = ThermostatSettings::getTemperatureRampFrequency();
+    _engine.getThermostat().setTemperatureRampingFrequency(frequency);
 
     const auto targetTemp   = ThermostatSettings::getTargetTemperature();
     const auto tempDelta    = targetTemp - startTemp;
-    const auto tempIncrease = tempDelta / double(steps) * frequency;
+    const auto updates      = steps / frequency + (steps % frequency != 0);
+    const auto tempIncrease = tempDelta / double(updates);
 
     _engine.getThermostat().setTemperatureIncrease(tempIncrease);
-    _engine.getThermostat().setTemperatureRampingFrequency(frequency);
 }
 
 void ThermostatSetup::writeSetupInfo() const

src/simulationBox/celllist.cpp

@@ -313,7 +313,27 @@ Vec3Dul CellList::getCellIndexOfAtom(
  * @brief resize cells
  *
  */
-void CellList::resizeCells() { _cells.resize(prod(_nCells)); }
+void CellList::resizeCells()
+{
+    auto numberOfCells = size_t{1};
+
+    for (size_t dimension = 0; dimension < 3; ++dimension)
+    {
+        if (0 == _nCells[dimension])
+            throw CellListException(
+                "Number of cells must be positive"
+            );   // GCOVR_EXCL_BR_LINE
+
+        if (_nCells[dimension] > _cells.max_size() / numberOfCells)
+            throw CellListException(
+                "Number of cells exceeds the supported size"
+            );   // GCOVR_EXCL_BR_LINE
+
+        numberOfCells *= _nCells[dimension];
+    }
+
+    _cells.resize(numberOfCells);
+}
 
 /**
  * @brief add cell to cell list

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,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;
-}
+}

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
@@ -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;
-}
+}

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();