Reflection move sampler

src/core/move_samplers/ReflectionSampler.cpp

@@ -0,0 +1,48 @@
+//
+// Created by ciesla on 1/1/25.
+//
+
+#include "ReflectionSampler.h"
+
+ReflectionSampler::ReflectionSampler(std::size_t flipEvery, const Vector<3> &plane) : flipEvery{flipEvery}, planeAxis{plane.normalized()}{
+    constexpr double EPSILON = 1e-12;
+    Expects(plane.norm2() > EPSILON * EPSILON);
+    Expects(flipEvery > 0);
+}
+
+std::size_t ReflectionSampler::getNumOfRequestedMoves(std::size_t numParticles) const {
+    Expects(numParticles > this->flipEvery);
+    return numParticles / this->flipEvery;
+}
+
+Matrix<3, 3, double> ReflectionSampler::getRotationMatrix(const Vector<3, double> &axis, double cosangle)
+{
+    double sina = -2*cosangle*std::sqrt(1-cosangle*cosangle);
+    double cosa = 1-2*cosangle*cosangle;
+    Matrix<3, 3> K = {{      0 , -axis[2],  axis[1],
+                             axis[2],        0, -axis[0],
+                             -axis[1],  axis[0],        0}};
+
+    // Rodrigues' rotation formula
+    return Matrix<3, 3>::identity() + K*sina + (1 - cosa)*K*K;
+}
+
+
+MoveSampler::MoveData ReflectionSampler::sampleMove(const Packing &packing, const std::vector<std::size_t> &particleIdxs, std::mt19937 &mt){
+    Expects(this->geometry != nullptr);
+
+    MoveData moveData;
+    std::uniform_int_distribution<std::size_t> particleDistribution(0, particleIdxs.size() - 1);
+    moveData.particleIdx = particleIdxs[particleDistribution(mt)];
+    Shape shape = packing[moveData.particleIdx];
+
+    Vector<3> shapeAxis = geometry->getPrimaryAxis(shape);
+    Vector<3> rotationAxis = (shapeAxis^this->planeAxis).normalized();
+
+    double cosangle = shapeAxis*this->planeAxis;
+    moveData.rotation = this->getRotationMatrix(rotationAxis, cosangle);
+    moveData.moveType = MoveType::ROTATION;
+    return moveData;
+}
+
+

src/core/move_samplers/ReflectionSampler.h

@@ -0,0 +1,56 @@
+//
+// Created by ciesla on 1/1/25.
+//
+
+#ifndef RAMPACK_REFLECTIONSAMPLER_H
+#define RAMPACK_REFLECTIONSAMPLER_H
+
+#include <variant>
+
+#include "core/MoveSampler.h"
+
+/**
+ * @brief ReflectionSampler performing the reflection relative to a given plane described by its normal vector.
+ * @details Particles are sampled at random. Reflection is performed in respect to a given plane which can be relative to
+ * simulation box coordinate system (@global=True) or to particle coordinate system.
+ * Internally it consists of a single move named @a reflection. The group name is @a reflection.
+ */
+class ReflectionSampler : public MoveSampler{
+
+private:
+    std::size_t flipEvery{};
+    const ShapeGeometry *geometry = nullptr;
+    const Vector<3> planeAxis{};
+
+    Matrix<3, 3, double> getRotationMatrix(const Vector<3, double> &axis, double cosangle);
+
+
+public:
+    /**
+     * @brief Constructs the class specifying how often to perform a reflection (i.e. how many moves should be requested,
+     * calculated by dividing the number of molecules by @a flipEvery).
+     */
+
+    explicit ReflectionSampler(std::size_t flipEvery, const Vector<3> &plane);
+
+    [[nodiscard]] std::string getName() const override { return "reflection"; }
+
+    [[nodiscard]] std::size_t getNumOfRequestedMoves(std::size_t numParticles) const override;
+
+    MoveData sampleMove(const Packing &packing, const std::vector<std::size_t> &particleIdxs,
+                        std::mt19937 &mt) override;
+
+    bool increaseStepSize() override { return false; }
+
+    bool decreaseStepSize() override { return false; }
+
+    [[nodiscard]] std::vector<std::pair<std::string, double>> getStepSizes() const override { return {{"rotation", 0}};}
+
+    void setStepSize(const std::string &, double) override { }
+
+    void setupForShapeTraits(const ShapeTraits &shapeTraits) override {
+        this->geometry = &shapeTraits.getGeometry();
+    }
+};
+
+#endif //RAMPACK_REFLECTIONSAMPLER_H

src/frontend/matchers/MoveSamplerMatcher.cpp

@@ -9,6 +9,7 @@
 #include "core/move_samplers/RotationSampler.h"
 #include "core/move_samplers/AxialRotationSampler.h"
 #include "core/move_samplers/FlipSampler.h"
+#include "core/move_samplers//ReflectionSampler.h"
 
 using namespace pyon::matcher;
 
@@ -19,6 +20,7 @@ namespace {
     MatcherDataclass create_rotation();
     MatcherDataclass create_axis_rotation();
     MatcherDataclass create_flip();
+    MatcherDataclass create_reflection();
 
 
     MatcherDataclass create_rototranslation() {
@@ -125,9 +127,20 @@ namespace {
                 return std::make_shared<FlipSampler>(every);
             });
     }
-}
 
+    MatcherDataclass create_reflection() {
+        return MatcherDataclass("reflection")
+                .arguments({{"every", MatcherInt{}.positive().mapTo<std::size_t>(), "10"},
+                            {"planeAxis", MatcherArray(MatcherFloat{}.mapTo<double>(),3), "[0,0,1]"}})
+                .mapTo([](const DataclassData &reflection) -> std::shared_ptr<MoveSampler> {
+                    auto every = reflection["every"].as<std::size_t>();
+                    auto planeAxisData = reflection["planeAxis"].as<pyon::matcher::ArrayData>();
+                    auto planeAxis = Vector<3, double>({planeAxisData[0].as<double>(), planeAxisData[1].as<double>(), planeAxisData[2].as<double>()});
+                    return std::make_shared<ReflectionSampler>(every, planeAxis);
+                });
+    }
+}
 
 MatcherAlternative MoveSamplerMatcher::create() {
-    return create_rototranslation() | create_translation() | create_rotation() | create_axis_rotation() | create_flip();
+    return create_rototranslation() | create_translation() | create_rotation() | create_axis_rotation() | create_flip() | create_reflection();
 }

test/unit_tests/core/move_samplers/ReflectionSamplerTest.cpp

@@ -0,0 +1,85 @@
+//
+// Created by Michal Ciesla on 14.01.25.
+//
+
+#include <catch2/catch.hpp>
+
+#include "matchers/VectorApproxMatcher.h"
+
+#include "core/move_samplers/ReflectionSampler.h"
+#include "core/move_samplers/RotationSampler.h"
+#include "core/shapes/PolyspherocylinderBananaTraits.h"
+#include "core/lattice/Lattice.h"
+#include "core/PeriodicBoundaryConditions.h"
+
+namespace {
+    void test_rotationAroundAxis_move(const ShapeTraits &traits) {
+        Lattice lattice(UnitCell(TriclinicBox(2), {Shape({0.5, 0.5, 0.5})}), {2, 2, 2});
+        auto pbc = std::make_unique<PeriodicBoundaryConditions>();
+
+        Vector<3, double> zeroVector({0, 0, 0});
+        Packing packing(lattice.getLatticeBox(), lattice.generateMolecules(), std::move(pbc), traits.getInteraction());
+        RotationSampler rotation(M_PI_2);
+        ReflectionSampler reflectionX(1, Vector<3, double>({1, 0, 0}));
+        reflectionX.setupForShapeTraits(traits);
+        ReflectionSampler reflectionY(1, Vector<3, double>({0, 1, 0}));
+        reflectionY.setupForShapeTraits(traits);
+        ReflectionSampler reflectionZ(1, Vector<3, double>({0, 0, 1}));
+        reflectionZ.setupForShapeTraits(traits);
+
+        const auto &interaction = traits.getInteraction();
+
+        std::vector<std::size_t> particleIdxs(packing.size());
+        std::iota(particleIdxs.begin(), particleIdxs.end(), 0);
+        std::mt19937 mt(1234);
+
+        for (auto i=0; i<100; i++){
+            auto move = rotation.sampleMove(packing, particleIdxs, mt);
+            packing.tryRotation(move.particleIdx, move.rotation, interaction);
+            packing.acceptRotation();
+        }
+
+        for (auto i = 0; i < 100; i++) {
+            auto moveX = reflectionX.sampleMove(packing, particleIdxs, mt);
+
+            CHECK(moveX.particleIdx < packing.size());
+            CHECK(moveX.moveType == MoveSampler::MoveType::ROTATION);
+            auto axis =  traits.getGeometry().getPrimaryAxis(packing[moveX.particleIdx]);
+            packing.tryRotation(moveX.particleIdx, moveX.rotation, interaction);
+            packing.acceptRotation();
+            auto axisX = traits.getGeometry().getPrimaryAxis(packing[moveX.particleIdx]);
+            REQUIRE_THAT(axis[0]+axisX[0], Catch::Matchers::WithinAbs(0, 1e-12));
+            REQUIRE_THAT(axis[1]-axisX[1], Catch::Matchers::WithinAbs(0, 1e-12));
+            REQUIRE_THAT(axis[2]-axisX[2], Catch::Matchers::WithinAbs(0, 1e-12));
+
+            auto moveY = reflectionY.sampleMove(packing, particleIdxs, mt);
+            axis =  traits.getGeometry().getPrimaryAxis(packing[moveY.particleIdx]);
+            packing.tryRotation(moveY.particleIdx, moveY.rotation, interaction);
+            packing.acceptRotation();
+            auto axisY = traits.getGeometry().getPrimaryAxis(packing[moveY.particleIdx]);
+            REQUIRE_THAT(axis[0]-axisY[0], Catch::Matchers::WithinAbs(0, 1e-12));
+            REQUIRE_THAT(axis[1]+axisY[1], Catch::Matchers::WithinAbs(0, 1e-12));
+            REQUIRE_THAT(axis[2]-axisY[2], Catch::Matchers::WithinAbs(0, 1e-12));
+
+            auto moveZ = reflectionZ.sampleMove(packing, particleIdxs, mt);
+            axis = traits.getGeometry().getPrimaryAxis(packing[moveZ.particleIdx]);
+            packing.tryRotation(moveZ.particleIdx, moveZ.rotation, interaction);
+            packing.acceptRotation();
+            auto axisZ = traits.getGeometry().getPrimaryAxis(packing[moveZ.particleIdx]);
+            REQUIRE_THAT(axis[0]-axisZ[0], Catch::Matchers::WithinAbs(0, 1e-12));
+            REQUIRE_THAT(axis[1]-axisZ[1], Catch::Matchers::WithinAbs(0, 1e-12));
+            REQUIRE_THAT(axis[2]+axisZ[2], Catch::Matchers::WithinAbs(0, 1e-12));
+
+            CHECK_THAT(moveX.translation, IsApproxEqual(zeroVector, 1e-12));
+            CHECK_THAT(moveY.translation, IsApproxEqual(zeroVector, 1e-12));
+            CHECK_THAT(moveZ.translation, IsApproxEqual(zeroVector, 1e-12));
+        }
+    }
+}
+
+TEST_CASE("ReflectionSampler") {
+    PolyspherocylinderBananaTraits traits(5, 2 * M_PI / 3, 2, 0.1, 1);
+    SECTION("move basic test") {
+        test_rotationAroundAxis_move(traits);
+    }
+}