From 0c827b1eee42d88ea42d5830518480199170b27e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefanie=20Kr=C3=B6ll?= Date: Mon, 19 May 2025 10:27:10 +0200 Subject: [PATCH 1/4] calculateMu isotropy divide all terms by 3 --- src/manostat/stochasticRescalingManostat.cpp | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/manostat/stochasticRescalingManostat.cpp b/src/manostat/stochasticRescalingManostat.cpp index 8f0fd602b..fc8f63038 100644 --- a/src/manostat/stochasticRescalingManostat.cpp +++ b/src/manostat/stochasticRescalingManostat.cpp @@ -56,7 +56,9 @@ StochasticRescalingManostat::StochasticRescalingManostat( : Manostat(other), _tau(other._tau), _compressibility(other._compressibility), - _dt(other._dt) {}; + _dt(other._dt) +{ +} /** * @brief Construct a new Stochastic Rescaling Manostat:: Stochastic Rescaling @@ -78,7 +80,9 @@ SemiIsotropicStochasticRescalingManostat:: ) : StochasticRescalingManostat(targetPressure, tau, compressibility), _2DAnisotropicAxis(anisotropicAxis), - _2DIsotropicAxes(isotropicAxes) {}; + _2DIsotropicAxes(isotropicAxes) +{ +} /** * @brief Construct a new Stochastic Rescaling Manostat:: Stochastic Rescaling @@ -155,7 +159,9 @@ tensor3D StochasticRescalingManostat::calculateMu(const double volume) const auto deltaP = _targetPressure - _pressure; - return diagonalMatrix(::exp(-compress * (deltaP) + stochasticFactor / 3.0)); + return diagonalMatrix( + ::exp(-compress * (deltaP) / 3.0 + stochasticFactor / 3.0) + ); } /** From 854c5af4fac13a09b0edb9347df5e07899c0d79c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefanie=20Kr=C3=B6ll?= Date: Mon, 19 May 2025 10:28:42 +0200 Subject: [PATCH 2/4] add pbc for imaging the positions back to new box --- src/simulationBox/molecule.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/simulationBox/molecule.cpp b/src/simulationBox/molecule.cpp index 6e5859b16..f69edccc9 100644 --- a/src/simulationBox/molecule.cpp +++ b/src/simulationBox/molecule.cpp @@ -40,14 +40,14 @@ using namespace settings; * * @param name */ -Molecule::Molecule(const std::string_view name) : _name(name){}; +Molecule::Molecule(const std::string_view name) : _name(name) {} /** * @brief Construct a new Molecule:: Molecule object * * @param moltype */ -Molecule::Molecule(const size_t moltype) : _moltype(moltype){}; +Molecule::Molecule(const size_t moltype) : _moltype(moltype) {} /** * @brief finds number of different atom types in molecule @@ -124,6 +124,9 @@ void Molecule::scale(const tensor3D &shiftTensor, const Box &box) if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) position = box.toSimSpace(position); + // add pbc to the position + box.applyPBC(position); + atom->setPosition(position); }; From 0606bcf976a44e155be08c202bf9a71a30135fc9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefanie=20Kr=C3=B6ll?= Date: Thu, 5 Jun 2025 09:28:26 +0200 Subject: [PATCH 3/4] rescaling atomwise --- src/manostat/stochasticRescalingManostat.cpp | 29 ++++++-- src/simulationBox/molecule.cpp | 7 +- tests/src/simulationBox/testMolecule.cpp | 73 +++++++++++++++++--- 3 files changed, 87 insertions(+), 22 deletions(-) diff --git a/src/manostat/stochasticRescalingManostat.cpp b/src/manostat/stochasticRescalingManostat.cpp index fc8f63038..59ac4a995 100644 --- a/src/manostat/stochasticRescalingManostat.cpp +++ b/src/manostat/stochasticRescalingManostat.cpp @@ -119,13 +119,6 @@ void StochasticRescalingManostat::applyManostat( const auto mu = calculateMu(simBox.getVolume()); - simBox.scaleBox(mu); - - physicalData.setVolume(simBox.getVolume()); - physicalData.setDensity(simBox.getDensity()); - - simBox.checkCoulRadiusCutOff(ExceptionType::MANOSTATEXCEPTION); - auto scalePositions = [&mu, &simBox](auto &molecule) { molecule.scale(mu, simBox.getBox()); }; @@ -133,8 +126,30 @@ void StochasticRescalingManostat::applyManostat( { atom->scaleVelocityOrthogonalSpace(inverse(mu), simBox.getBox()); }; std::ranges::for_each(simBox.getMolecules(), scalePositions); + std::ranges::for_each(simBox.getAtoms(), scaleVelocities); + simBox.scaleBox(mu); + + physicalData.setVolume(simBox.getVolume()); + physicalData.setDensity(simBox.getDensity()); + + simBox.checkCoulRadiusCutOff(ExceptionType::MANOSTATEXCEPTION); + + // auto pbcPositions = [&simBox](auto &molecule) + // { + // molecule.calculateCenterOfMass(simBox.getBox()); + // for (std::size_t i = 0; i < molecule.getNumberOfAtoms(); ++i) + // { + // auto position = molecule.getAtomPosition(i); + // position += molecule.getCenterOfMass(); + // simBox.applyPBC(position); + // molecule.setAtomPosition(i, position); + // } + // }; + + // std::ranges::for_each(simBox.getMolecules(), pbcPositions); + stopTimingsSection("Stochastic Rescaling"); } diff --git a/src/simulationBox/molecule.cpp b/src/simulationBox/molecule.cpp index f69edccc9..d90a55442 100644 --- a/src/simulationBox/molecule.cpp +++ b/src/simulationBox/molecule.cpp @@ -110,7 +110,7 @@ void Molecule::scale(const tensor3D &shiftTensor, const Box &box) if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) centerOfMass = box.toOrthoSpace(_centerOfMass); - const auto shift = shiftTensor * centerOfMass - centerOfMass; + const auto shift = shiftTensor; // * centerOfMass - centerOfMass; auto scaleAtomPosition = [&box, shift](auto atom) { @@ -119,14 +119,11 @@ void Molecule::scale(const tensor3D &shiftTensor, const Box &box) if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) position = box.toOrthoSpace(position); - position += shift; + position = shift * position; if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) position = box.toSimSpace(position); - // add pbc to the position - box.applyPBC(position); - atom->setPosition(position); }; diff --git a/tests/src/simulationBox/testMolecule.cpp b/tests/src/simulationBox/testMolecule.cpp index 3dacabd72..61d164449 100644 --- a/tests/src/simulationBox/testMolecule.cpp +++ b/tests/src/simulationBox/testMolecule.cpp @@ -22,6 +22,8 @@ #include "testMolecule.hpp" +#include // for std::ranges::for_each + #include "gtest/gtest.h" // for Message, TestPartResult #include "orthorhombicBox.hpp" // for OrthorhombicBox #include "staticMatrix.hpp" // for diagonalMatrix @@ -39,26 +41,77 @@ TEST_F(TestMolecule, calculateCenterOfMass) TEST_F(TestMolecule, scaleAtoms) { - const linearAlgebra::tensor3D scale = - diagonalMatrix(linearAlgebra::Vec3D{1.0, 2.0, 3.0}); - const linearAlgebra::Vec3D atomPosition1 = _molecule->getAtomPosition(0); - const linearAlgebra::Vec3D atomPosition2 = _molecule->getAtomPosition(1); - const linearAlgebra::Vec3D atomPosition3 = _molecule->getAtomPosition(2); + // const linearAlgebra::tensor3D scale = + // diagonalMatrix(linearAlgebra::Vec3D{1.0, 2.0, 3.0}); + // const linearAlgebra::Vec3D atomPosition1 = _molecule->getAtomPosition(0); + // const linearAlgebra::Vec3D atomPosition2 = _molecule->getAtomPosition(1); + // const linearAlgebra::Vec3D atomPosition3 = _molecule->getAtomPosition(2); + + // simulationBox::OrthorhombicBox box; + // box.setBoxDimensions({10.0, 10.0, 10.0}); + + // _molecule->calculateCenterOfMass(box); + + // const auto centerOfMassBeforeScaling = _molecule->getCenterOfMass(); + // const linearAlgebra::Vec3D shift = + // centerOfMassBeforeScaling * (diagonal(scale) - 1.0); + + // _molecule->scale(scale, box); + + // EXPECT_EQ(_molecule->getAtomPosition(0), atomPosition1 + shift); + // EXPECT_EQ(_molecule->getAtomPosition(1), atomPosition2 + shift); + // EXPECT_EQ(_molecule->getAtomPosition(2), atomPosition3 + shift); + + // Test atoms scaling on the boundaries simulationBox::OrthorhombicBox box; - box.setBoxDimensions({10.0, 10.0, 10.0}); + simulationBox::OrthorhombicBox box_control; + box.setBoxDimensions({1.0, 2.0, 3.0}); + box_control.setBoxDimensions({1.0, 2.0, 3.0}); + _molecule->setNumberOfAtoms(3); + _molecule->setAtomPosition(0, linearAlgebra::Vec3D(0.5, 0.0, 0.0)); + _molecule->setAtomPosition(1, linearAlgebra::Vec3D(0.0, 1.5, 0.0)); + _molecule->setAtomPosition(2, linearAlgebra::Vec3D(0.0, 0.0, 2.5)); + + EXPECT_EQ( + _molecule->getAtomPosition(0), + linearAlgebra::Vec3D(0.5, 0.0, 0.0) + ); + EXPECT_EQ( + _molecule->getAtomPosition(1), + linearAlgebra::Vec3D(0.0, 1.5, 0.0) + ); + EXPECT_EQ( + _molecule->getAtomPosition(2), + linearAlgebra::Vec3D(0.0, 0.0, 2.5) + ); + + const linearAlgebra::tensor3D scale = + diagonalMatrix(linearAlgebra::Vec3D{-2.0, -2.0, -2.0}); + + const linearAlgebra::Vec3D atomPosition_1 = _molecule->getAtomPosition(0); + const linearAlgebra::Vec3D atomPosition_2 = _molecule->getAtomPosition(1); + const linearAlgebra::Vec3D atomPosition_3 = _molecule->getAtomPosition(2); _molecule->calculateCenterOfMass(box); const auto centerOfMassBeforeScaling = _molecule->getCenterOfMass(); + const linearAlgebra::Vec3D shift = - centerOfMassBeforeScaling * (diagonal(scale) - 1.0); + diagonal(scale); //* centerOfMassBeforeScaling - centerOfMassBeforeScaling; + + linearAlgebra::Vec3D pos1 = shift * atomPosition_1; + linearAlgebra::Vec3D pos2 = shift * atomPosition_2; + linearAlgebra::Vec3D pos3 = shift * atomPosition_3; _molecule->scale(scale, box); - EXPECT_EQ(_molecule->getAtomPosition(0), atomPosition1 + shift); - EXPECT_EQ(_molecule->getAtomPosition(1), atomPosition2 + shift); - EXPECT_EQ(_molecule->getAtomPosition(2), atomPosition3 + shift); + // scale box + box_control.scaleBox(scale); + + EXPECT_EQ(_molecule->getAtomPosition(0), pos1); + EXPECT_EQ(_molecule->getAtomPosition(1), pos2); + EXPECT_EQ(_molecule->getAtomPosition(2), pos3); } TEST_F(TestMolecule, setAtomForceToZero) From dbc7d0d07062af3bb59a1fc0aa3faf627a0ddef5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefanie=20Kr=C3=B6ll?= Date: Fri, 6 Jun 2025 15:59:34 +0200 Subject: [PATCH 4/4] try of center of mass rescaling --- include/simulationBox/molecule.hpp | 1 + src/manostat/stochasticRescalingManostat.cpp | 19 ++---- src/simulationBox/molecule.cpp | 65 +++++++++++++++++++- tests/src/simulationBox/testMolecule.cpp | 48 +++++++++++++-- 4 files changed, 114 insertions(+), 19 deletions(-) diff --git a/include/simulationBox/molecule.hpp b/include/simulationBox/molecule.hpp index 178a38dda..ecf6eb2c1 100644 --- a/include/simulationBox/molecule.hpp +++ b/include/simulationBox/molecule.hpp @@ -67,6 +67,7 @@ namespace simulationBox void calculateCenterOfMass(const Box &); void scale(const pq::tensor3D &, const Box &); + void decenter(const Box &box); [[nodiscard]] size_t getNumberOfAtomTypes(); [[nodiscard]] std::vector getExternalGlobalVDWTypes() const; diff --git a/src/manostat/stochasticRescalingManostat.cpp b/src/manostat/stochasticRescalingManostat.cpp index 59ac4a995..5966fee06 100644 --- a/src/manostat/stochasticRescalingManostat.cpp +++ b/src/manostat/stochasticRescalingManostat.cpp @@ -136,19 +136,12 @@ void StochasticRescalingManostat::applyManostat( simBox.checkCoulRadiusCutOff(ExceptionType::MANOSTATEXCEPTION); - // auto pbcPositions = [&simBox](auto &molecule) - // { - // molecule.calculateCenterOfMass(simBox.getBox()); - // for (std::size_t i = 0; i < molecule.getNumberOfAtoms(); ++i) - // { - // auto position = molecule.getAtomPosition(i); - // position += molecule.getCenterOfMass(); - // simBox.applyPBC(position); - // molecule.setAtomPosition(i, position); - // } - // }; - - // std::ranges::for_each(simBox.getMolecules(), pbcPositions); + simBox.calculateCenterOfMass(); + + auto decenterPositions = [&simBox](auto &molecule) + { molecule.decenter(simBox.getBox()); }; + + std::ranges::for_each(simBox.getMolecules(), decenterPositions); stopTimingsSection("Stochastic Rescaling"); } diff --git a/src/simulationBox/molecule.cpp b/src/simulationBox/molecule.cpp index d90a55442..f685f1f57 100644 --- a/src/simulationBox/molecule.cpp +++ b/src/simulationBox/molecule.cpp @@ -103,6 +103,32 @@ void Molecule::calculateCenterOfMass(const Box &box) * * @param shiftFactors */ +// void Molecule::scale(const tensor3D &shiftTensor, const Box &box) +// { +// auto centerOfMass = _centerOfMass; + +// if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) +// centerOfMass = box.toOrthoSpace(_centerOfMass); + +// const auto shift = shiftTensor; // * centerOfMass - centerOfMass; + +// auto scaleAtomPosition = [&box, shift](auto atom) +// { +// auto position = atom->getPosition(); + +// if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) +// position = box.toOrthoSpace(position); + +// position = shift * position; + +// if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) +// position = box.toSimSpace(position); + +// atom->setPosition(position); +// }; + +// std::ranges::for_each(_atoms, scaleAtomPosition); +// } void Molecule::scale(const tensor3D &shiftTensor, const Box &box) { auto centerOfMass = _centerOfMass; @@ -110,7 +136,7 @@ void Molecule::scale(const tensor3D &shiftTensor, const Box &box) if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) centerOfMass = box.toOrthoSpace(_centerOfMass); - const auto shift = shiftTensor; // * centerOfMass - centerOfMass; + const auto shift = shiftTensor * centerOfMass - centerOfMass; auto scaleAtomPosition = [&box, shift](auto atom) { @@ -119,7 +145,7 @@ void Molecule::scale(const tensor3D &shiftTensor, const Box &box) if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) position = box.toOrthoSpace(position); - position = shift * position; + position += shift; if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) position = box.toSimSpace(position); @@ -129,6 +155,41 @@ void Molecule::scale(const tensor3D &shiftTensor, const Box &box) std::ranges::for_each(_atoms, scaleAtomPosition); } +/** + * @brief decenter the positions of the molecule after bringing molecule into + * center of mass + * + * @details decenter has to be done in orthogonal space since pressure scaling + * is done in orthogonal space + * + * @param box + */ +void Molecule::decenter(const Box &box) +{ + auto centerOfMass = _centerOfMass; + + if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) + centerOfMass = box.toOrthoSpace(_centerOfMass); + + auto decenterAtomPosition = [&box, centerOfMass](auto atom) + { + auto position = atom->getPosition(); + + if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) + position = box.toOrthoSpace(position); + + position += centerOfMass; + + if (ManostatSettings::getIsotropy() != Isotropy::FULL_ANISOTROPIC) + position = box.toSimSpace(position); + + box.applyPBC(position); + + atom->setPosition(position); + }; + + std::ranges::for_each(_atoms, decenterAtomPosition); +} /** * @brief returns the external global vdw types of the atoms in the molecule diff --git a/tests/src/simulationBox/testMolecule.cpp b/tests/src/simulationBox/testMolecule.cpp index 61d164449..47fd4a66a 100644 --- a/tests/src/simulationBox/testMolecule.cpp +++ b/tests/src/simulationBox/testMolecule.cpp @@ -98,11 +98,11 @@ TEST_F(TestMolecule, scaleAtoms) const auto centerOfMassBeforeScaling = _molecule->getCenterOfMass(); const linearAlgebra::Vec3D shift = - diagonal(scale); //* centerOfMassBeforeScaling - centerOfMassBeforeScaling; + diagonal(scale) * centerOfMassBeforeScaling - centerOfMassBeforeScaling; - linearAlgebra::Vec3D pos1 = shift * atomPosition_1; - linearAlgebra::Vec3D pos2 = shift * atomPosition_2; - linearAlgebra::Vec3D pos3 = shift * atomPosition_3; + linearAlgebra::Vec3D pos1 = atomPosition_1 + shift; + linearAlgebra::Vec3D pos2 = atomPosition_2 + shift; + linearAlgebra::Vec3D pos3 = atomPosition_3 + shift; _molecule->scale(scale, box); @@ -114,6 +114,46 @@ TEST_F(TestMolecule, scaleAtoms) EXPECT_EQ(_molecule->getAtomPosition(2), pos3); } +TEST_F(TestMolecule, decenterPositions) +{ + simulationBox::OrthorhombicBox box; + box.setBoxDimensions({1.0, 2.0, 3.0}); + _molecule->setNumberOfAtoms(3); + _molecule->setAtomPosition(0, linearAlgebra::Vec3D(0.5, 0.0, 0.0)); + _molecule->setAtomPosition(1, linearAlgebra::Vec3D(0.0, 1.5, 0.0)); + _molecule->setAtomPosition(2, linearAlgebra::Vec3D(0.0, 0.0, 2.5)); + const auto centerOfMassBeforeScaling = _molecule->getCenterOfMass(); + const linearAlgebra::tensor3D scale = + diagonalMatrix(linearAlgebra::Vec3D{-2.0, -2.0, -2.0}); + const linearAlgebra::Vec3D atomPosition_1 = _molecule->getAtomPosition(0); + const linearAlgebra::Vec3D atomPosition_2 = _molecule->getAtomPosition(1); + const linearAlgebra::Vec3D atomPosition_3 = _molecule->getAtomPosition(2); + + _molecule->scale(scale, box); + _molecule->calculateCenterOfMass(box); + const auto centerOfMassAfterScaling = _molecule->getCenterOfMass(); + + _molecule->decenter(box); + const linearAlgebra::Vec3D shift = + diagonal(scale) * centerOfMassBeforeScaling - centerOfMassBeforeScaling; + + linearAlgebra::Vec3D pos1 = + atomPosition_1 + shift + centerOfMassAfterScaling; + linearAlgebra::Vec3D pos2 = + atomPosition_2 + shift + centerOfMassAfterScaling; + linearAlgebra::Vec3D pos3 = + atomPosition_3 + shift + centerOfMassAfterScaling; + + box.applyPBC(pos1); + box.applyPBC(pos2); + box.applyPBC(pos3); + + EXPECT_EQ(_molecule->getAtomPosition(0), pos1); + EXPECT_EQ(_molecule->getAtomPosition(1), pos2); + EXPECT_EQ(_molecule->getAtomPosition(2), pos3); + EXPECT_EQ(_molecule->getCenterOfMass(), centerOfMassAfterScaling); +} + TEST_F(TestMolecule, setAtomForceToZero) { _molecule->setAtomForcesToZero();