From ffa85e71c09a9603479f87f782734ce926239e29 Mon Sep 17 00:00:00 2001 From: Armin Penz Date: Wed, 19 Aug 2026 15:48:39 +0200 Subject: [PATCH] feat: add fixed_axis keyword to enable 2d pressure coupling --- .../user/enhancement.2d-pressure-couling.md | 1 + .../inputFileParser/manostatInputParser.hpp | 1 + include/manostat/berendsenManostat.hpp | 17 ++- .../manostat/stochasticRescalingManostat.hpp | 17 +-- include/settings/manostatSettings.hpp | 13 +++ include/setup/manostatSetup.hpp | 2 + .../inputFileParser/manostatInputParser.cpp | 42 +++++++ src/manostat/berendsenManostat.cpp | 92 +++++++++++++-- src/manostat/stochasticRescalingManostat.cpp | 105 ++++++++++++++++-- src/settings/manostatSettings.cpp | 31 ++++++ src/setup/manostatSetup.cpp | 67 +++++++---- tests/data/inputFileReader/keywordList.txt | 1 + tests/src/manostat/testManostat.cpp | 73 +++++++++--- 13 files changed, 395 insertions(+), 67 deletions(-) create mode 100644 changes/user/enhancement.2d-pressure-couling.md diff --git a/changes/user/enhancement.2d-pressure-couling.md b/changes/user/enhancement.2d-pressure-couling.md new file mode 100644 index 000000000..96dd95a2f --- /dev/null +++ b/changes/user/enhancement.2d-pressure-couling.md @@ -0,0 +1 @@ +- Enable isotropic, anisotropic and full-anisotropic 2d pressure coupling by adding the "fixed_axis" keyword diff --git a/include/input/inputFileParser/manostatInputParser.hpp b/include/input/inputFileParser/manostatInputParser.hpp index b792551a5..07c65cd76 100644 --- a/include/input/inputFileParser/manostatInputParser.hpp +++ b/include/input/inputFileParser/manostatInputParser.hpp @@ -56,6 +56,7 @@ namespace input ); void parseIsotropy(const std::vector &, const size_t); + void parseFixedAxis(const std::vector &, const size_t); }; } // namespace input diff --git a/include/manostat/berendsenManostat.hpp b/include/manostat/berendsenManostat.hpp index e3e87d336..7f7dc0a75 100644 --- a/include/manostat/berendsenManostat.hpp +++ b/include/manostat/berendsenManostat.hpp @@ -40,12 +40,18 @@ namespace manostat class BerendsenManostat : public Manostat { protected: - double _tau; - double _compressibility; - double _dt; + double _tau; + double _compressibility; + double _dt; + settings::FixedAxis _fixedAxis; public: - explicit BerendsenManostat(const double, const double, const double); + explicit BerendsenManostat( + const double, + const double, + const double, + const settings::FixedAxis + ); void applyManostat( simulationBox::SimulationBox &simBox, @@ -83,7 +89,8 @@ namespace manostat const double, const double, const size_t, - const std::vector & + const std::vector &, + const settings::FixedAxis ); [[nodiscard]] linearAlgebra::tensor3D calculateMu() const override; diff --git a/include/manostat/stochasticRescalingManostat.hpp b/include/manostat/stochasticRescalingManostat.hpp index 8dbdf345a..4dfba5873 100644 --- a/include/manostat/stochasticRescalingManostat.hpp +++ b/include/manostat/stochasticRescalingManostat.hpp @@ -40,16 +40,18 @@ namespace manostat protected: randomNumberGenerator::RandomNumberGenerator _randomNumberGenerator{}; - double _tau; - double _compressibility; - double _dt; + double _tau; + double _compressibility; + double _dt; + settings::FixedAxis _fixedAxis; public: StochasticRescalingManostat() = default; explicit StochasticRescalingManostat( - const double targetPressure, - const double tau, - const double compressibility + const double targetPressure, + const double tau, + const double compressibility, + const settings::FixedAxis fixedAxis ); ~StochasticRescalingManostat() override = default; @@ -97,7 +99,8 @@ namespace manostat const double tau, const double compressibility, const size_t anisotropicAxis, - const std::vector &isotropicAxes + const std::vector &isotropicAxes, + const settings::FixedAxis fixedAxis ); [[nodiscard]] diff --git a/include/settings/manostatSettings.hpp b/include/settings/manostatSettings.hpp index 503b51739..43b6a5a97 100644 --- a/include/settings/manostatSettings.hpp +++ b/include/settings/manostatSettings.hpp @@ -60,6 +60,14 @@ namespace settings FULL_ANISOTROPIC }; + enum class FixedAxis + { + NONE, + X, + Y, + Z + }; + [[nodiscard]] std::string string(const ManostatType &manostatType); [[nodiscard]] std::string string(const Isotropy &isotropy); @@ -74,6 +82,7 @@ namespace settings private: static inline ManostatType _manostatType = ManostatType::NONE; static inline Isotropy _isotropy = Isotropy::ISOTROPIC; + static inline FixedAxis _fixedAxis = FixedAxis::NONE; static inline double _targetPressure; @@ -99,6 +108,9 @@ namespace settings static void setIsotropy(const std::string_view &isotropy); static void setIsotropy(const Isotropy &isotropy); + static void setFixedAxis(const std::string_view &fixedAxis); + static void setFixedAxis(const FixedAxis &fixedAxis); + static void setTargetPressure(const double targetPressure); static void setTauManostat(const double tauManostat); static void setCompressibility(const double compressibility); @@ -113,6 +125,7 @@ namespace settings [[nodiscard]] static ManostatType getManostatType(); [[nodiscard]] static Isotropy getIsotropy(); + [[nodiscard]] static FixedAxis getFixedAxis(); [[nodiscard]] static double getTargetPressure(); [[nodiscard]] static double getTauManostat(); [[nodiscard]] static double getCompressibility(); diff --git a/include/setup/manostatSetup.hpp b/include/setup/manostatSetup.hpp index c6416c297..afac79a56 100644 --- a/include/setup/manostatSetup.hpp +++ b/include/setup/manostatSetup.hpp @@ -53,6 +53,8 @@ namespace setup void setupBerendsenManostat(); void setupStochasticRescalingManostat(); + void validateIsotropyFixedAxisCombination() const; + void writeSetupInfo() const; void writeManostatSelection() const; void writeBerendsenSetup() const; diff --git a/src/input/inputFileParser/manostatInputParser.cpp b/src/input/inputFileParser/manostatInputParser.cpp index 3f0c154aa..1574d5610 100644 --- a/src/input/inputFileParser/manostatInputParser.cpp +++ b/src/input/inputFileParser/manostatInputParser.cpp @@ -86,6 +86,12 @@ ManostatInputParser::ManostatInputParser(Engine &engine) bindMember(&ManostatInputParser::parseIsotropy, this), false ); + + addKeyword( + std::string("fixed_axis"), + bindMember(&ManostatInputParser::parseFixedAxis, this), + false + ); } /** @@ -288,3 +294,39 @@ void ManostatInputParser::parseIsotropy( ); } } + +void ManostatInputParser::parseFixedAxis( + const std::vector &lineElements, + const size_t lineNumber +) +{ + checkCommand(lineElements, lineNumber); + + const auto fixed_axis = toLowerAndReplaceDashesCopy(lineElements[2]); + + using enum FixedAxis; + + if (fixed_axis == "none") + ManostatSettings::setFixedAxis(NONE); + + else if (fixed_axis == "x") + ManostatSettings::setFixedAxis(X); + + else if (fixed_axis == "y") + ManostatSettings::setFixedAxis(Y); + + else if (fixed_axis == "z") + ManostatSettings::setFixedAxis(Z); + + else + { + throw InputFileException( + std::format( + "Invalid fixed_axis \"{}\" at line {} in input file.\n" + "Possible options are: none, x, y, z", + lineElements[2], + lineNumber + ) + ); + } +} diff --git a/src/manostat/berendsenManostat.cpp b/src/manostat/berendsenManostat.cpp index c7ef0d6a0..3374a17b0 100644 --- a/src/manostat/berendsenManostat.cpp +++ b/src/manostat/berendsenManostat.cpp @@ -45,16 +45,19 @@ using namespace physicalData; * @param targetPressure * @param tau * @param compressibility + * @param fixedAxis */ BerendsenManostat::BerendsenManostat( - const double targetPressure, - const double tau, - const double compressibility + const double targetPressure, + const double tau, + const double compressibility, + const FixedAxis fixedAxis ) : Manostat(targetPressure), _tau(tau), _compressibility(compressibility), - _dt(TimingsSettings::getTimeStep()) + _dt(TimingsSettings::getTimeStep()), + _fixedAxis(fixedAxis) { } @@ -64,15 +67,19 @@ BerendsenManostat::BerendsenManostat( * @param targetPressure * @param tau * @param compressibility + * @param anisotropicAxis + * @param isotropicAxes + * @param fixedAxis */ SemiIsotropicBerendsenManostat::SemiIsotropicBerendsenManostat( const double targetPressure, const double tau, const double compressibility, const size_t anisotropicAxis, - const std::vector &isotropicAxes + const std::vector &isotropicAxes, + const FixedAxis fixedAxis ) - : BerendsenManostat(targetPressure, tau, compressibility), + : BerendsenManostat(targetPressure, tau, compressibility, fixedAxis), _2DAnisotropicAxis(anisotropicAxis), _2DIsotropicAxes(isotropicAxes) { @@ -118,14 +125,47 @@ void BerendsenManostat::applyManostat( /** * @brief calculate mu as scaling factor for Berendsen manostat (isotropic) * + * @details If a fixed axis is specified, that axis is not scaled (mu = 1.0) + * and the remaining axes are scaled isotropically + * * @return tensor3D */ tensor3D BerendsenManostat::calculateMu() const { - const auto p = trace(_pressureTensor) / 3.0; + using enum FixedAxis; + const auto preFactor = _compressibility * _dt / _tau; - return diagonalMatrix(::cbrt(1.0 - preFactor * (_targetPressure - p))); + // 2D pressure coupling + if (_fixedAxis != NONE) + { + const auto fixedAxisIndex = static_cast(_fixedAxis) - 1; + const auto p_xyz = diagonal(_pressureTensor); + + // Calculate average pressure of non-fixed axes + double p_avg = 0.0; + for (size_t i = 0; i < 3; ++i) + if (i != fixedAxisIndex) + p_avg += p_xyz[i]; + + p_avg /= 2.0; + + // Scale factor for non-fixed axes + const auto mu_2D = ::sqrt(1.0 - preFactor * (_targetPressure - p_avg)); + + Vec3D mu = {1.0, 1.0, 1.0}; + for (size_t i = 0; i < 3; ++i) + if (i != fixedAxisIndex) + mu[i] = mu_2D; + + return diagonalMatrix(mu); + } + + // 3D pressure coupling + const auto p = trace(_pressureTensor) / 3.0; + const auto mu_scalar = ::cbrt(1.0 - preFactor * (_targetPressure - p)); + + return diagonalMatrix(mu_scalar); } /** @@ -162,30 +202,64 @@ tensor3D SemiIsotropicBerendsenManostat::calculateMu() const /** * @brief calculate mu as scaling factor for Berendsen manostat (anisotropic) * + * @details If a fixed axis is specified, that axis is not scaled (mu = 1.0) + * and the other axes are scaled independently + * * @return tensor3D */ tensor3D AnisotropicBerendsenManostat::calculateMu() const { + using enum FixedAxis; + const auto pxyz = diagonal(_pressureTensor); const auto preFactor = _compressibility * _dt / _tau; - return diagonalMatrix(1.0 - preFactor * (_targetPressure - pxyz)); + auto mu = 1.0 - preFactor * (_targetPressure - pxyz); + + // 2D pressure coupling + if (_fixedAxis != NONE) + { + const auto fixedAxisIndex = static_cast(_fixedAxis) - 1; + mu[fixedAxisIndex] = 1.0; + } + + return diagonalMatrix(mu); } /** * @brief calculate mu as scaling factor for Berendsen manostat (full * anisotropic including angles) * + * @details If a fixed axis is specified, the corresponding row and column + * are zeroed (no coupling with other axes) and the diagonal is set to 1.0 + * * @return tensor3D */ tensor3D FullAnisotropicBerendsenManostat::calculateMu() const { + using enum FixedAxis; + const auto pTarget = diagonalMatrix(_targetPressure); const auto preFactor = _compressibility * _dt / _tau; const auto kronecker = kroneckerDeltaMatrix(); auto mu = kronecker - preFactor * (pTarget - _pressureTensor); + // 2D full anisotropic: fix one axis and remove its coupling + if (_fixedAxis != NONE) + { + const auto fixedAxisIndex = static_cast(_fixedAxis) - 1; + + // Zero out the row and column of the fixed axis + for (size_t i = 0; i < 3; ++i) + { + mu[fixedAxisIndex][i] = 0.0; + mu[i][fixedAxisIndex] = 0.0; + } + // Set diagonal to 1.0 (no scaling) + mu[fixedAxisIndex][fixedAxisIndex] = 1.0; + } + rotateMu(mu); return mu; diff --git a/src/manostat/stochasticRescalingManostat.cpp b/src/manostat/stochasticRescalingManostat.cpp index 59390a244..f0e3bc9f4 100644 --- a/src/manostat/stochasticRescalingManostat.cpp +++ b/src/manostat/stochasticRescalingManostat.cpp @@ -55,7 +55,8 @@ StochasticRescalingManostat::StochasticRescalingManostat( : Manostat(other), _tau(other._tau), _compressibility(other._compressibility), - _dt(other._dt) + _dt(other._dt), + _fixedAxis(other._fixedAxis) { } @@ -75,6 +76,7 @@ StochasticRescalingManostat &StochasticRescalingManostat::operator=( _tau = other._tau; _compressibility = other._compressibility; _dt = other._dt; + _fixedAxis = other._fixedAxis; } return *this; } @@ -87,6 +89,7 @@ StochasticRescalingManostat &StochasticRescalingManostat::operator=( * @param compressibility * @param anisotropicAxis * @param isotropicAxes + * @param fixedAxis * @return SemiIsotropicStochasticRescalingManostat:: */ SemiIsotropicStochasticRescalingManostat:: @@ -95,9 +98,15 @@ SemiIsotropicStochasticRescalingManostat:: const double tau, const double compressibility, const size_t anisotropicAxis, - const std::vector &isotropicAxes + const std::vector &isotropicAxes, + const FixedAxis fixedAxis ) - : StochasticRescalingManostat(targetPressure, tau, compressibility), + : StochasticRescalingManostat( + targetPressure, + tau, + compressibility, + fixedAxis + ), _2DAnisotropicAxis(anisotropicAxis), _2DIsotropicAxes(isotropicAxes) { @@ -110,16 +119,19 @@ SemiIsotropicStochasticRescalingManostat:: * @param targetPressure * @param tau * @param compressibility + * @param fixedAxis */ StochasticRescalingManostat::StochasticRescalingManostat( - const double targetPressure, - const double tau, - const double compressibility + const double targetPressure, + const double tau, + const double compressibility, + const FixedAxis fixedAxis ) : Manostat(targetPressure), _tau(tau), _compressibility(compressibility), - _dt(TimingsSettings::getTimeStep()) + _dt(TimingsSettings::getTimeStep()), + _fixedAxis(fixedAxis) { } @@ -168,19 +180,57 @@ void StochasticRescalingManostat::applyManostat( * @brief calculate mu as scaling factor for Stochastic Rescaling manostat * (isotropic) * + * @details If a fixed axis is specified, that axis is not scaled (mu = 1.0) + * and the remaining axes are scaled isotropically with stochastic coupling + * * @param volume * @return Vec3D */ tensor3D StochasticRescalingManostat::calculateMu(const double volume) { + using enum FixedAxis; + const auto compress = _compressibility * _dt / _tau; const auto kb = BOLTZMANN_CONSTANT_IN_KCAL_PER_MOL; const auto kT = kb * ThermostatSettings::getActualTargetTemperature(); const auto random = _randomNumberGenerator.getNormalDistribution(0.0, 1.0); + // 2D pressure coupling + if (_fixedAxis != NONE) + { + const auto fixedAxisIndex = static_cast(_fixedAxis) - 1; + const auto p_xyz = diagonal(_pressureTensor); + + // Calculate average pressure of non-fixed axes + double p_avg = 0.0; + for (size_t i = 0; i < 3; ++i) + if (i != fixedAxisIndex) + p_avg += p_xyz[i]; + p_avg /= 2.0; + + auto stochasticFactor = 2.0 * kT * compress / volume; + stochasticFactor *= PRESSURE_FACTOR; + stochasticFactor = ::sqrt(stochasticFactor) * random; + + const auto deltaP = _targetPressure - p_avg; + + // 2D isotropic scaling + constexpr auto dimension = 2.0; + const auto mu_2D = + ::exp((-compress * deltaP + stochasticFactor) / dimension); + + Vec3D mu = {1.0, 1.0, 1.0}; + for (size_t i = 0; i < 3; ++i) + if (i != fixedAxisIndex) + mu[i] = mu_2D; + + return diagonalMatrix(mu); + } + + // 3D pressure coupling auto stochasticFactor = 2.0 * kT * compress / volume; - stochasticFactor *= PRESSURE_FACTOR; + stochasticFactor *= PRESSURE_FACTOR; stochasticFactor = ::sqrt(stochasticFactor) * random; const auto deltaP = _targetPressure - _pressure; @@ -243,6 +293,9 @@ tensor3D SemiIsotropicStochasticRescalingManostat::calculateMu( * @brief calculate mu as scaling factor for Stochastic Rescaling manostat * (anisotropic) * + * @details If a fixed axis is specified, that axis is not scaled (mu = 1.0) + * and the other axes are scaled independently with stochastic coupling + * * @param volume * @return Vec3D */ @@ -250,6 +303,8 @@ tensor3D AnisotropicStochasticRescalingManostat::calculateMu( const double volume ) { + using enum FixedAxis; + const auto compress = _compressibility * _dt / _tau; const auto kb = BOLTZMANN_CONSTANT_IN_KCAL_PER_MOL; @@ -262,15 +317,26 @@ tensor3D AnisotropicStochasticRescalingManostat::calculateMu( const auto deltaP = _targetPressure - diagonal(_pressureTensor); - return diagonalMatrix( - exp(-compress * (deltaP) / _pressureTensor.size + stochasticFactor) - ); + auto mu = + exp(-compress * (deltaP) / _pressureTensor.size + stochasticFactor); + + // 2D anisotropic: fix one axis + if (_fixedAxis != NONE) + { + const auto fixedAxisIndex = static_cast(_fixedAxis) - 1; + mu[fixedAxisIndex] = 1.0; + } + + return diagonalMatrix(mu); } /** * @brief calculate mu as scaling factor for Stochastic Rescaling manostat (full * anisotropic including angles) * + * @details If a fixed axis is specified, the corresponding row and column + * are zeroed (no coupling with other axes) and the diagonal is set to 1.0 + * * @param volume * @return tensor3D */ @@ -278,6 +344,8 @@ tensor3D FullAnisotropicStochasticRescalingManostat::calculateMu( const double volume ) { + using enum FixedAxis; + const auto compress = _compressibility * _dt / _tau; const auto kb = BOLTZMANN_CONSTANT_IN_KCAL_PER_MOL; @@ -292,6 +360,21 @@ tensor3D FullAnisotropicStochasticRescalingManostat::calculateMu( auto mu = expPade(-compress * deltaP / _pressureTensor.size + stochasticFactor); + // 2D full anisotropic: fix one axis and remove its coupling + if (_fixedAxis != NONE) + { + const auto fixedAxisIndex = static_cast(_fixedAxis) - 1; + + // Zero out the row and column of the fixed axis + for (size_t i = 0; i < 3; ++i) + { + mu[fixedAxisIndex][i] = 0.0; + mu[i][fixedAxisIndex] = 0.0; + } + // Set diagonal to 1.0 (no scaling) + mu[fixedAxisIndex][fixedAxisIndex] = 1.0; + } + rotateMu(mu); return mu; diff --git a/src/settings/manostatSettings.cpp b/src/settings/manostatSettings.cpp index 906bc3ebb..f7577e2e2 100644 --- a/src/settings/manostatSettings.cpp +++ b/src/settings/manostatSettings.cpp @@ -140,6 +140,30 @@ void ManostatSettings::setIsotropy(const Isotropy &isotropy) _isotropy = isotropy; } +void ManostatSettings::setFixedAxis(const std::string_view &fixedAxis) +{ + using enum FixedAxis; + const auto FixedAxisToLower = + utilities::toLowerAndReplaceDashesCopy(fixedAxis); + + if (FixedAxisToLower == "none") + _fixedAxis = NONE; + + else if (FixedAxisToLower == "x") + _fixedAxis = X; + + else if (FixedAxisToLower == "y") + _fixedAxis = Y; + + else if (FixedAxisToLower == "z") + _fixedAxis = Z; +} + +void ManostatSettings::setFixedAxis(const FixedAxis &fixedAxis) +{ + _fixedAxis = fixedAxis; +} + /** * @brief sets the targetPressure to double in settings * @@ -222,6 +246,13 @@ ManostatType ManostatSettings::getManostatType() { return _manostatType; } */ Isotropy ManostatSettings::getIsotropy() { return _isotropy; } +/** + * @brief get the FixedAxis + * + * @return FixedAxis + */ +FixedAxis ManostatSettings::getFixedAxis() { return _fixedAxis; } + /** * @brief get the target pressure * diff --git a/src/setup/manostatSetup.cpp b/src/setup/manostatSetup.cpp index cb5ca1744..e1dd6969d 100644 --- a/src/setup/manostatSetup.cpp +++ b/src/setup/manostatSetup.cpp @@ -27,6 +27,7 @@ #include "berendsenManostat.hpp" // for BerendsenManostat #include "constants/conversionFactors.hpp" // for _PS_TO_FS_ +#include "exceptions.hpp" #include "manostat.hpp" // for BerendsenManostat, Manostat, manostat #include "manostatSettings.hpp" // for ManostatSettings #include "mdEngine.hpp" // for Engine @@ -38,6 +39,7 @@ using namespace engine; using namespace settings; using namespace manostat; using namespace constants; +using namespace customException; /** * @brief wrapper for setupManostat @@ -84,6 +86,7 @@ void ManostatSetup::setup() else _engine.makeManostat(Manostat()); + validateIsotropyFixedAxisCombination(); writeSetupInfo(); } @@ -95,12 +98,13 @@ void ManostatSetup::setup() */ void ManostatSetup::setupBerendsenManostat() { - const auto isotropy = ManostatSettings::getIsotropy(); - const auto pTarget = ManostatSettings::getTargetPressure(); - const auto tau = ManostatSettings::getTauManostat() * PS_TO_FS; - const auto compress = ManostatSettings::getCompressibility(); - const auto aniso = ManostatSettings::get2DAnisotropicAxis(); - const auto iso = ManostatSettings::get2DIsotropicAxes(); + const auto isotropy = ManostatSettings::getIsotropy(); + const auto pTarget = ManostatSettings::getTargetPressure(); + const auto tau = ManostatSettings::getTauManostat() * PS_TO_FS; + const auto compress = ManostatSettings::getCompressibility(); + const auto aniso = ManostatSettings::get2DAnisotropicAxis(); + const auto iso = ManostatSettings::get2DIsotropicAxes(); + const auto fixedAxis = ManostatSettings::getFixedAxis(); switch (isotropy) { @@ -108,20 +112,20 @@ void ManostatSetup::setupBerendsenManostat() // clang-format off case SEMI_ISOTROPIC: - _engine.makeManostat(SemiIsotropicBerendsenManostat(pTarget, tau, compress, aniso, iso)); + _engine.makeManostat(SemiIsotropicBerendsenManostat(pTarget, tau, compress, aniso, iso, fixedAxis)); break; case ANISOTROPIC: - _engine.makeManostat(AnisotropicBerendsenManostat(pTarget, tau, compress)); + _engine.makeManostat(AnisotropicBerendsenManostat(pTarget, tau, compress, fixedAxis)); break; case FULL_ANISOTROPIC: - _engine.makeManostat(FullAnisotropicBerendsenManostat(pTarget, tau, compress)); + _engine.makeManostat(FullAnisotropicBerendsenManostat(pTarget, tau, compress, fixedAxis)); break; case NONE: // fall through case ISOTROPIC: - _engine.makeManostat(BerendsenManostat(pTarget, tau, compress)); + _engine.makeManostat(BerendsenManostat(pTarget, tau, compress, fixedAxis)); // clang-format on } @@ -135,12 +139,13 @@ void ManostatSetup::setupBerendsenManostat() */ void ManostatSetup::setupStochasticRescalingManostat() { - const auto isotropy = ManostatSettings::getIsotropy(); - const auto pTarget = ManostatSettings::getTargetPressure(); - const auto tau = ManostatSettings::getTauManostat() * PS_TO_FS; - const auto compress = ManostatSettings::getCompressibility(); - const auto aniso = ManostatSettings::get2DAnisotropicAxis(); - const auto iso = ManostatSettings::get2DIsotropicAxes(); + const auto isotropy = ManostatSettings::getIsotropy(); + const auto pTarget = ManostatSettings::getTargetPressure(); + const auto tau = ManostatSettings::getTauManostat() * PS_TO_FS; + const auto compress = ManostatSettings::getCompressibility(); + const auto aniso = ManostatSettings::get2DAnisotropicAxis(); + const auto iso = ManostatSettings::get2DIsotropicAxes(); + const auto fixedAxis = ManostatSettings::getFixedAxis(); switch (isotropy) { @@ -149,25 +154,47 @@ void ManostatSetup::setupStochasticRescalingManostat() // clang-format off case SEMI_ISOTROPIC: - _engine.makeManostat(SemiIsotropicStochasticRescalingManostat(pTarget, tau, compress, aniso, iso)); + _engine.makeManostat(SemiIsotropicStochasticRescalingManostat(pTarget, tau, compress, aniso, iso, fixedAxis)); break; case ANISOTROPIC: - _engine.makeManostat(AnisotropicStochasticRescalingManostat(pTarget, tau, compress)); + _engine.makeManostat(AnisotropicStochasticRescalingManostat(pTarget, tau, compress, fixedAxis)); break; case FULL_ANISOTROPIC: - _engine.makeManostat(FullAnisotropicStochasticRescalingManostat(pTarget, tau, compress)); + _engine.makeManostat(FullAnisotropicStochasticRescalingManostat(pTarget, tau, compress, fixedAxis)); break; case NONE: // fall through case ISOTROPIC: - _engine.makeManostat(StochasticRescalingManostat(pTarget, tau, compress)); + _engine.makeManostat(StochasticRescalingManostat(pTarget, tau, compress, fixedAxis)); // clang-format on } } +/** + * @brief validate isotropy and fixed_axis combination + * + * @throws SetupException if semi-isotropic mode conflicts with fixed_axis + */ +void ManostatSetup::validateIsotropyFixedAxisCombination() const +{ + using enum Isotropy; + + const auto isotropy = ManostatSettings::getIsotropy(); + const auto fixedAxis = ManostatSettings::getFixedAxis(); + + if (isotropy == SEMI_ISOTROPIC && fixedAxis != FixedAxis::NONE) + { + throw UserInputException( + "Invalid combination: semi-isotropic pressure coupling cannot " + "be used while any axis is fixed. For isotropic 2D pressure " + "coupling, use isotropy = isotropic with fixed_axis instead." + ); + } +} + /** * @brief write setup info * diff --git a/tests/data/inputFileReader/keywordList.txt b/tests/data/inputFileReader/keywordList.txt index 8fed9ea4f..eede282cd 100644 --- a/tests/data/inputFileReader/keywordList.txt +++ b/tests/data/inputFileReader/keywordList.txt @@ -64,6 +64,7 @@ pressure false p_relaxation false compressibility false isotropy false +fixed_axis false thermostat false temp false diff --git a/tests/src/manostat/testManostat.cpp b/tests/src/manostat/testManostat.cpp index d61a9ccea..75215c80e 100644 --- a/tests/src/manostat/testManostat.cpp +++ b/tests/src/manostat/testManostat.cpp @@ -190,7 +190,12 @@ TEST_F(TestManostat, testApplyBerendsenManostat) _box->addMolecule(molecule); settings::TimingsSettings::setTimeStep(0.5); - _manostat = new manostat::BerendsenManostat(1.0, 0.1, 4.5); + _manostat = new manostat::BerendsenManostat( + 1.0, + 0.1, + 4.5, + settings::FixedAxis::NONE + ); const auto scaleFactors = linearAlgebra::Vec3D( ::pow( @@ -223,7 +228,12 @@ TEST_F(TestManostat, testApplyBerendsenManostatPreservesCutMoleculeGeometry) { setupCutMolecule(*_box, *_data); - auto manostat = manostat::BerendsenManostat(1.0, 1.0, 0.058808); + auto manostat = manostat::BerendsenManostat( + 1.0, + 1.0, + 0.058808, + settings::FixedAxis::NONE + ); manostat.applyManostat(*_box, *_data); expectCutMoleculeScaled(*_box); @@ -241,8 +251,12 @@ TEST_F( setupCutMolecule(*_box, *_data); settings::ThermostatSettings::setActualTargetTemperature(0.0); - auto manostat = - manostat::StochasticRescalingManostat(-3.0 * ::log(0.98), 1.0, 1.0); + auto manostat = manostat::StochasticRescalingManostat( + -3.0 * ::log(0.98), + 1.0, + 1.0, + settings::FixedAxis::NONE + ); manostat.applyManostat(*_box, *_data); expectCutMoleculeScaled(*_box); @@ -276,8 +290,12 @@ TEST_F( _box->addMolecule(molecule); _box->setTotalMass(4.0); - auto manostat = - manostat::StochasticRescalingManostat(-3.0 * ::log(0.98), 1.0, 1.0); + auto manostat = manostat::StochasticRescalingManostat( + -3.0 * ::log(0.98), + 1.0, + 1.0, + settings::FixedAxis::NONE + ); manostat.applyManostat(*_box, *_data); const auto cutDistance = getMinimumImageDistance(*_box, 0); @@ -305,7 +323,8 @@ TEST_F( _manostat = new manostat::BerendsenManostat( 3.0 * constants::PRESSURE_FACTOR, 0.1, - 4.5 + 4.5, + settings::FixedAxis::NONE ); EXPECT_THROW_MSG( @@ -332,7 +351,12 @@ TEST_F(TestManostat, stochasticRescalingMuUsesLengthScaling) settings::ThermostatSettings::setActualTargetTemperature(0.0); settings::TimingsSettings::setTimeStep(0.5); - auto manostat = TestableStochasticRescalingManostat(7.0, 0.25, 0.12); + auto manostat = TestableStochasticRescalingManostat( + 7.0, + 0.25, + 0.12, + settings::FixedAxis::NONE + ); manostat.setPressure(1.0); const auto mu = manostat.calculateMu(10.0); @@ -379,7 +403,12 @@ TEST_F(TestManostat, stochasticRescalingPreservesInternalMolecularVelocities) molecule.calculateCenterOfMass(_box->getBox()); _box->addMolecule(molecule); - _manostat = new manostat::StochasticRescalingManostat(7.0, 0.25, 0.12); + _manostat = new manostat::StochasticRescalingManostat( + 7.0, + 0.25, + 0.12, + settings::FixedAxis::NONE + ); const auto mu = ::exp(-(0.12 * 0.5 / 0.25) * (7.0 - 0.0) / 3.0); const auto expectedCenterOfMassVelocity = @@ -435,20 +464,23 @@ TEST_F(TestManostat, testRotateMu) TEST_F(TestManostat, berendsenTauAndCompressibilityGetters) { - auto bm = manostat::BerendsenManostat(1.0, 0.1, 4.5); + auto bm = + manostat::BerendsenManostat(1.0, 0.1, 4.5, settings::FixedAxis::NONE); EXPECT_DOUBLE_EQ(bm.getTau(), 0.1); EXPECT_DOUBLE_EQ(bm.getCompressibility(), 4.5); } TEST_F(TestManostat, berendsenManostatType) { - auto bm = manostat::BerendsenManostat(1.0, 0.1, 4.5); + auto bm = + manostat::BerendsenManostat(1.0, 0.1, 4.5, settings::FixedAxis::NONE); EXPECT_EQ(bm.getManostatType(), settings::ManostatType::BERENDSEN); } TEST_F(TestManostat, berendsenIsotropy) { - auto bm = manostat::BerendsenManostat(1.0, 0.1, 4.5); + auto bm = + manostat::BerendsenManostat(1.0, 0.1, 4.5, settings::FixedAxis::NONE); EXPECT_EQ(bm.getIsotropy(), settings::Isotropy::ISOTROPIC); } @@ -459,7 +491,8 @@ TEST_F(TestManostat, semiIsotropicBerendsenIsotropy) 0.1, 4.5, 2U, - std::vector{0U, 1U} + std::vector{0U, 1U}, + settings::FixedAxis::NONE ); EXPECT_EQ(bm.getIsotropy(), settings::Isotropy::SEMI_ISOTROPIC); EXPECT_EQ(bm.getManostatType(), settings::ManostatType::BERENDSEN); @@ -467,14 +500,24 @@ TEST_F(TestManostat, semiIsotropicBerendsenIsotropy) TEST_F(TestManostat, anisotropicBerendsenIsotropy) { - auto bm = manostat::AnisotropicBerendsenManostat(1.0, 0.1, 4.5); + auto bm = manostat::AnisotropicBerendsenManostat( + 1.0, + 0.1, + 4.5, + settings::FixedAxis::NONE + ); EXPECT_EQ(bm.getIsotropy(), settings::Isotropy::ANISOTROPIC); EXPECT_EQ(bm.getManostatType(), settings::ManostatType::BERENDSEN); } TEST_F(TestManostat, fullAnisotropicBerendsenIsotropy) { - auto bm = manostat::FullAnisotropicBerendsenManostat(1.0, 0.1, 4.5); + auto bm = manostat::FullAnisotropicBerendsenManostat( + 1.0, + 0.1, + 4.5, + settings::FixedAxis::NONE + ); EXPECT_EQ(bm.getIsotropy(), settings::Isotropy::FULL_ANISOTROPIC); EXPECT_EQ(bm.getManostatType(), settings::ManostatType::BERENDSEN); }