fix: guard simulation edge cases

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

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/simulationBox/celllist.cpp

@@ -313,7 +313,25 @@ 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");
+
+        if (_nCells[dimension] > _cells.max_size() / numberOfCells)
+            throw CellListException(
+                "Number of cells exceeds the supported size"
+            );
+
+        numberOfCells *= _nCells[dimension];
+    }
+
+    _cells.resize(numberOfCells);
+}
 
 /**
  * @brief add cell to cell list

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

tests/src/simulationBox/testCelllist.cpp

@@ -22,6 +22,7 @@
 
 #include "testCelllist.hpp"
 
+#include <limits>   // for numeric_limits
 #include <memory>   // for make_shared, __shared_ptr_access
 #include <string>   // for allocator, basic_string
 #include <vector>   // for vector
@@ -265,6 +266,28 @@ TEST_F(TestCellList, activateDeactivateToggles_isActive)
     EXPECT_TRUE(_cellList->isActive());
 }
 
+TEST_F(TestCellList, resizeCellsRejectsOverflow)
+{
+    _cellList->setNumberOfCells(std::numeric_limits<int>::max());
+
+    EXPECT_THROW_MSG(
+        _cellList->resizeCells(),
+        customException::CellListException,
+        "Number of cells exceeds the supported size"
+    );
+}
+
+TEST_F(TestCellList, resizeCellsRejectsZeroDimensions)
+{
+    _cellList->setNumberOfCells(0);
+
+    EXPECT_THROW_MSG(
+        _cellList->resizeCells(),
+        customException::CellListException,
+        "Number of cells must be positive"
+    );
+}
+
 /* ---------- clone() copies the configured cell counts ---------- */
 
 TEST_F(TestCellList, clone_preservesNumberOfCellsAndNeighbourCells)

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
@@ -220,6 +221,38 @@ TEST_F(TestThermostat, applyBerendsen_zeroTemperatureNoNaN)
         }
 }
 
+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_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, 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 ---------- */
 
 TEST_F(TestThermostat, langevin_constructorComputesSigma)

tests/src/utilities/CMakeLists.txt

@@ -2,13 +2,15 @@ set(source_files
     testStringUtilities.cpp
     testMathUtilities.cpp
     testCollectionUtilities.cpp
+    testProgressbar.cpp
 )
 
 foreach(source_file ${source_files})
     get_filename_component(test_name ${source_file} NAME_WE)
     add_executable(${test_name} ${source_file})
     target_include_directories(${test_name}
         PRIVATE
+        ${PROJECT_SOURCE_DIR}/external/progressbar/include
         ${PROJECT_SOURCE_DIR}/tests/include/utilities
         ${PROJECT_SOURCE_DIR}/tests/include/macros
     )

tests/src/utilities/testProgressbar.cpp

@@ -0,0 +1,38 @@
+/*****************************************************************************
+<GPL_HEADER>
+
+    PQ
+    Copyright (C) 2023-now  Jakob Gamper
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+<GPL_HEADER>
+******************************************************************************/
+
+#include <gtest/gtest.h>
+
+#include <sstream>
+#include <string>
+
+#include "progressbar.hpp"
+
+TEST(TestProgressbar, singleIterationCompletes)
+{
+    auto output = std::ostringstream();
+    auto bar    = progressbar(1, true, output);
+
+    bar.update();
+
+    EXPECT_NE(output.str().find("100%"), std::string::npos);
+}