From e7e466c6ac611d638a778938f06eb3f267c00eab Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Tue, 28 Jul 2026 17:34:18 +0200 Subject: [PATCH 1/7] fix: define parameter reader destructor out of line --- include/input/parameterFileReader/parameterFileReader.hpp | 3 ++- src/input/parameterFileReader/parameterFileReader.cpp | 4 +++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/include/input/parameterFileReader/parameterFileReader.hpp b/include/input/parameterFileReader/parameterFileReader.hpp index 257626e09..bab66a76e 100644 --- a/include/input/parameterFileReader/parameterFileReader.hpp +++ b/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_ \ No newline at end of file +#endif // _PARAMETER_FILE_READER_HPP_ diff --git a/src/input/parameterFileReader/parameterFileReader.cpp b/src/input/parameterFileReader/parameterFileReader.cpp index 116132157..85e056730 100644 --- a/src/input/parameterFileReader/parameterFileReader.cpp +++ b/src/input/parameterFileReader/parameterFileReader.cpp @@ -68,6 +68,8 @@ ParameterFileReader::ParameterFileReader( _parameterFileSections.push_back(make_unique()); } +ParameterFileReader::~ParameterFileReader() = default; + /** * @brief determines which section of the parameter file the header line belongs * to @@ -219,4 +221,4 @@ std::vector> &ParameterFileReader:: const std::string &ParameterFileReader::getFilename() const { return _fileName; -} \ No newline at end of file +} From ecf1bbd613f1ddea90d4dcea3a561f64065dad4c 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 2/7] fix: validate cell-list dimensions --- src/simulationBox/celllist.cpp | 22 +++++++++++++++++++++- tests/src/simulationBox/testCelllist.cpp | 23 +++++++++++++++++++++++ 2 files changed, 44 insertions(+), 1 deletion(-) diff --git a/src/simulationBox/celllist.cpp b/src/simulationBox/celllist.cpp index e3c6d422b..380d0c699 100644 --- a/src/simulationBox/celllist.cpp +++ b/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 diff --git a/tests/src/simulationBox/testCelllist.cpp b/tests/src/simulationBox/testCelllist.cpp index 56e929b69..9ee21fe87 100644 --- a/tests/src/simulationBox/testCelllist.cpp +++ b/tests/src/simulationBox/testCelllist.cpp @@ -22,6 +22,7 @@ #include "testCelllist.hpp" +#include // for numeric_limits #include // for make_shared, __shared_ptr_access #include // for allocator, basic_string #include // for vector @@ -265,6 +266,28 @@ TEST_F(TestCellList, activateDeactivateToggles_isActive) EXPECT_TRUE(_cellList->isActive()); } +TEST_F(TestCellList, resizeCellsRejectsOverflow) +{ + _cellList->setNumberOfCells(std::numeric_limits::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) From 4a9cb42a93fd64fee7f3ea223fece3d3a6ff6870 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 3/7] fix: complete single-step progress bars --- external/progressbar/include/progressbar.hpp | 2 +- tests/src/utilities/CMakeLists.txt | 6 ++++ tests/src/utilities/testProgressbar.cpp | 38 ++++++++++++++++++++ 3 files changed, 45 insertions(+), 1 deletion(-) create mode 100644 tests/src/utilities/testProgressbar.cpp diff --git a/external/progressbar/include/progressbar.hpp b/external/progressbar/include/progressbar.hpp index 50ab6d089..eb154c66c 100644 --- a/external/progressbar/include/progressbar.hpp +++ b/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 diff --git a/tests/src/utilities/CMakeLists.txt b/tests/src/utilities/CMakeLists.txt index cbfd29b18..95af86ef6 100644 --- a/tests/src/utilities/CMakeLists.txt +++ b/tests/src/utilities/CMakeLists.txt @@ -2,6 +2,7 @@ set(source_files testStringUtilities.cpp testMathUtilities.cpp testCollectionUtilities.cpp + testProgressbar.cpp ) foreach(source_file ${source_files}) @@ -27,6 +28,11 @@ foreach(source_file ${source_files}) set_property(TEST ${test_name} PROPERTY LABELS utilties) endforeach() +target_include_directories(testProgressbar + PRIVATE + ${PROJECT_SOURCE_DIR}/external/progressbar/include +) + if(${BUILD_WITH_GCOVR}) include(CodeCoverage) setup_target_for_coverage_gcovr_html( diff --git a/tests/src/utilities/testProgressbar.cpp b/tests/src/utilities/testProgressbar.cpp new file mode 100644 index 000000000..531cd548d --- /dev/null +++ b/tests/src/utilities/testProgressbar.cpp @@ -0,0 +1,38 @@ +/***************************************************************************** + + + 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 . + + +******************************************************************************/ + +#include + +#include +#include + +#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); +} From 34cf0eea7818b14661fb6e94647ae80d0bc17155 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 4/7] 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 5c60f5d7b814693a501da9a1d04a8493cd9aa30d Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Tue, 28 Jul 2026 17:47:15 +0200 Subject: [PATCH 5/7] fix: complete partial temperature ramps --- src/setup/thermostatSetup.cpp | 3 ++- tests/src/setup/testThermostatSetup.cpp | 27 +++++++++++++++++++++++++ 2 files changed, 29 insertions(+), 1 deletion(-) diff --git a/src/setup/thermostatSetup.cpp b/src/setup/thermostatSetup.cpp index c9f4277c7..04585e341 100644 --- a/src/setup/thermostatSetup.cpp +++ b/src/setup/thermostatSetup.cpp @@ -251,7 +251,8 @@ void ThermostatSetup::setupTemperatureRamp() const auto targetTemp = ThermostatSettings::getTargetTemperature(); const auto tempDelta = targetTemp - startTemp; - const auto tempIncrease = tempDelta / double(steps) * frequency; + const auto updates = (steps + frequency - 1) / frequency; + const auto tempIncrease = tempDelta / double(updates); _engine.getThermostat().setTemperatureIncrease(tempIncrease); _engine.getThermostat().setTemperatureRampingFrequency(frequency); diff --git a/tests/src/setup/testThermostatSetup.cpp b/tests/src/setup/testThermostatSetup.cpp index 16e377f0b..413990b73 100644 --- a/tests/src/setup/testThermostatSetup.cpp +++ b/tests/src/setup/testThermostatSetup.cpp @@ -99,7 +99,34 @@ TEST_F(TestSetup, setupThermostat_temp_ramping) thermostatSetup.getEngine().getThermostat().getRampingFrequency(), 2 ); +} + +TEST_F(TestSetup, temperatureRampReachesEndWithPartialFinalInterval) +{ + ThermostatSetup thermostatSetup(*_mdEngine); + + settings::TimingsSettings::setNumberOfSteps(10); + settings::ThermostatSettings::setThermostatType("berendsen"); + settings::ThermostatSettings::setTargetTemperature(300); + settings::ThermostatSettings::setStartTemperature(200); + settings::ThermostatSettings::setTemperatureRampSteps(10); + settings::ThermostatSettings::setTemperatureRampFrequency(3); + settings::ThermostatSettings::setEndTemperatureSet(false); + + thermostatSetup.setup(); + + EXPECT_DOUBLE_EQ( + thermostatSetup.getEngine().getThermostat().getTemperatureIncrease(), + 25.0 + ); + + for (size_t step = 0; step < 10; ++step) + thermostatSetup.getEngine().getThermostat().applyTemperatureRamping(); + EXPECT_DOUBLE_EQ( + thermostatSetup.getEngine().getThermostat().getTargetTemperature(), + 300.0 + ); } TEST_F(TestSetup, setupThermostat_only_end_temp_defined) From 6f437f6e9091968070983db9db61e50ec9c551cf Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Tue, 28 Jul 2026 21:48:27 +0200 Subject: [PATCH 6/7] fix: validate temperature ramp schedules --- docs/sphinx/src/userGuide/inputFile.rst | 2 + src/setup/thermostatSetup.cpp | 40 +++++++++++-------- tests/src/setup/testThermostatSetup.cpp | 51 +++++++++++++++++++++---- 3 files changed, 70 insertions(+), 23 deletions(-) diff --git a/docs/sphinx/src/userGuide/inputFile.rst b/docs/sphinx/src/userGuide/inputFile.rst index 30e04b017..4a960b795 100644 --- a/docs/sphinx/src/userGuide/inputFile.rst +++ b/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: diff --git a/src/setup/thermostatSetup.cpp b/src/setup/thermostatSetup.cpp index 04585e341..605eeee94 100644 --- a/src/setup/thermostatSetup.cpp +++ b/src/setup/thermostatSetup.cpp @@ -225,37 +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 updates = (steps + frequency - 1) / 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 diff --git a/tests/src/setup/testThermostatSetup.cpp b/tests/src/setup/testThermostatSetup.cpp index 413990b73..d92fca40d 100644 --- a/tests/src/setup/testThermostatSetup.cpp +++ b/tests/src/setup/testThermostatSetup.cpp @@ -28,7 +28,9 @@ #include "berendsenThermostat.hpp" // for BerendsenThermostat #include "constants/conversionFactors.hpp" // for _FS_TO_S_, _KG_TO_GRAM_ #include "constants/natureConstants.hpp" // for _UNIVERSAL_GAS_CONSTANT_ +#include "exceptions.hpp" // for InputFileException #include "gtest/gtest.h" // for Message, TestPartResult +#include "inputFileReader.hpp" // for InputFileReader #include "langevinThermostat.hpp" // for LangevinThermostat #include "noseHooverThermostat.hpp" // for NoseHooverThermostat #include "testSetup.hpp" // for TestSetup @@ -103,15 +105,15 @@ TEST_F(TestSetup, setupThermostat_temp_ramping) TEST_F(TestSetup, temperatureRampReachesEndWithPartialFinalInterval) { - ThermostatSetup thermostatSetup(*_mdEngine); + ThermostatSetup thermostatSetup(*_mdEngine); + input::InputFileReader reader("input.in", *_mdEngine); - settings::TimingsSettings::setNumberOfSteps(10); - settings::ThermostatSettings::setThermostatType("berendsen"); - settings::ThermostatSettings::setTargetTemperature(300); - settings::ThermostatSettings::setStartTemperature(200); - settings::ThermostatSettings::setTemperatureRampSteps(10); - settings::ThermostatSettings::setTemperatureRampFrequency(3); - settings::ThermostatSettings::setEndTemperatureSet(false); + reader.process({"nstep", "=", "10"}); + reader.process({"thermostat", "=", "berendsen"}); + reader.process({"temp", "=", "300"}); + reader.process({"start_temp", "=", "200"}); + reader.process({"temp_ramp_steps", "=", "10"}); + reader.process({"temp_ramp_frequency", "=", "3"}); thermostatSetup.setup(); @@ -127,6 +129,39 @@ TEST_F(TestSetup, temperatureRampReachesEndWithPartialFinalInterval) thermostatSetup.getEngine().getThermostat().getTargetTemperature(), 300.0 ); + + settings::ThermostatSettings::setTemperatureRampSteps(0); + settings::ThermostatSettings::setTemperatureRampFrequency(1); +} + +TEST_F(TestSetup, rejectsEmptyTemperatureRamp) +{ + ThermostatSetup thermostatSetup(*_mdEngine); + + settings::TimingsSettings::setNumberOfSteps(0); + settings::ThermostatSettings::setThermostatType("berendsen"); + settings::ThermostatSettings::setTargetTemperature(300); + settings::ThermostatSettings::setStartTemperature(200); + settings::ThermostatSettings::setTemperatureRampSteps(0); + + EXPECT_THROW(thermostatSetup.setup(), customException::InputFileException); +} + +TEST_F(TestSetup, rejectsZeroTemperatureRampFrequency) +{ + ThermostatSetup thermostatSetup(*_mdEngine); + + settings::TimingsSettings::setNumberOfSteps(10); + settings::ThermostatSettings::setThermostatType("berendsen"); + settings::ThermostatSettings::setTargetTemperature(300); + settings::ThermostatSettings::setStartTemperature(200); + settings::ThermostatSettings::setTemperatureRampSteps(10); + settings::ThermostatSettings::setTemperatureRampFrequency(0); + + EXPECT_THROW(thermostatSetup.setup(), customException::InputFileException); + + settings::ThermostatSettings::setTemperatureRampSteps(0); + settings::ThermostatSettings::setTemperatureRampFrequency(1); } TEST_F(TestSetup, setupThermostat_only_end_temp_defined) From 58c98bd53f2a40ce16ee88264b3c6b43637ef82d Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Thu, 30 Jul 2026 08:42:55 +0200 Subject: [PATCH 7/7] chore: add changelog fragment for #385 --- changes/simulation-edge-cases.user.bugfix.md | 1 + 1 file changed, 1 insertion(+) create mode 100644 changes/simulation-edge-cases.user.bugfix.md diff --git a/changes/simulation-edge-cases.user.bugfix.md b/changes/simulation-edge-cases.user.bugfix.md new file mode 100644 index 000000000..1e5f194e9 --- /dev/null +++ b/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.