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/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/simulationBox/celllist.cpp b/src/simulationBox/celllist.cpp index e3c6d422b..a14b31ad5 100644 --- a/src/simulationBox/celllist.cpp +++ b/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 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/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) diff --git a/tests/src/thermostat/testThermostat.cpp b/tests/src/thermostat/testThermostat.cpp index ab0303b7d..87c656273 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 @@ -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) diff --git a/tests/src/utilities/CMakeLists.txt b/tests/src/utilities/CMakeLists.txt index cbfd29b18..c6ffff0c1 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}) @@ -9,6 +10,7 @@ foreach(source_file ${source_files}) 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 ) 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); +}