Bugfix/rescale

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<size_t> getExternalGlobalVDWTypes() const;

src/manostat/stochasticRescalingManostat.cpp

@@ -56,7 +56,9 @@
     : 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 @@
     )
     : StochasticRescalingManostat(targetPressure, tau, compressibility),
       _2DAnisotropicAxis(anisotropicAxis),
-      _2DIsotropicAxes(isotropicAxes) {};
+      _2DIsotropicAxes(isotropicAxes)
+{
+}
 
 /**
  * @brief Construct a new Stochastic Rescaling Manostat:: Stochastic Rescaling
@@ -115,21 +119,29 @@
 
     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());
     physicalData.setDensity(simBox.getDensity());
 
     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 @@
 
     const auto deltaP = _targetPressure - _pressure;
 
-    return diagonalMatrix(::exp(-compress * (deltaP) + stochasticFactor / 3.0));
+    return diagonalMatrix(
+        ::exp(-compress * (deltaP) / 3.0 + stochasticFactor / 3.0)
+    );
 }
 
 /**

src/simulationBox/molecule.cpp

@@ -40,14 +40,14 @@
  *
  * @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 @@
  *
  * @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 @@
 
     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

tests/src/simulationBox/testMolecule.cpp

@@ -22,6 +22,8 @@
 
 #include "testMolecule.hpp"
 
+#include <algorithm>   // 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)