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 8f0fd602b..5966fee06 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 @@ -115,6 +119,16 @@ void StochasticRescalingManostat::applyManostat( const auto mu = calculateMu(simBox.getVolume()); + auto scalePositions = [&mu, &simBox](auto &molecule) + { molecule.scale(mu, simBox.getBox()); }; + + auto scaleVelocities = [&mu, &simBox](auto &atom) + { 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()); @@ -122,14 +136,12 @@ void StochasticRescalingManostat::applyManostat( simBox.checkCoulRadiusCutOff(ExceptionType::MANOSTATEXCEPTION); - auto scalePositions = [&mu, &simBox](auto &molecule) - { molecule.scale(mu, simBox.getBox()); }; + simBox.calculateCenterOfMass(); - auto scaleVelocities = [&mu, &simBox](auto &atom) - { atom->scaleVelocityOrthogonalSpace(inverse(mu), simBox.getBox()); }; + auto decenterPositions = [&simBox](auto &molecule) + { molecule.decenter(simBox.getBox()); }; - std::ranges::for_each(simBox.getMolecules(), scalePositions); - std::ranges::for_each(simBox.getAtoms(), scaleVelocities); + std::ranges::for_each(simBox.getMolecules(), decenterPositions); stopTimingsSection("Stochastic Rescaling"); } @@ -155,7 +167,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) + ); } /** diff --git a/src/simulationBox/molecule.cpp b/src/simulationBox/molecule.cpp index 6e5859b16..f685f1f57 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 @@ -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; @@ -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 3dacabd72..47fd4a66a 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,117 @@ 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 = atomPosition_1 + shift; + linearAlgebra::Vec3D pos2 = atomPosition_2 + shift; + linearAlgebra::Vec3D pos3 = atomPosition_3 + shift; _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, 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)