2 files changed, 59 insertions, 8 deletions

diff --git a/soft_spheres.cpp b/soft_spheres.cpp
index 10e90c7..e438362 100644
--- a/soft_spheres.cpp
+++ b/soft_spheres.cpp
@@ -21,10 +21,10 @@ public:
 int main(int argc, char* argv[]) {
   const unsigned n = 24;
 
-  unsigned N = 1200;
+  unsigned N = 1000;
   double R = 0.025;
   double ε = 0.005;
-  unsigned steps = 1e2;
+  unsigned steps = 1e3;
   double β = 0.5;
 
   initializeAnimation(argc, argv);
@@ -45,15 +45,14 @@ int main(int argc, char* argv[]) {
       draw(m);
   }
 
-  for (unsigned i = 0; i < steps; i++) {
+  for (unsigned i = 0; i < 10 * steps; i++) {
     for (unsigned j = 0; j < N; j++) {
       m.randomMove(β, ε, r);
     }
 
-    unsigned k = m.randomClusterSwap(β, r);
-    std::cout << k << std::endl;
+    unsigned k = m.randomCluster(β, r);
 
-    if (k > 0 && k < N) {
+    if (k > 2 && k < N - 2) {
       draw(m);
       getchar();
     }
diff --git a/spheres.hpp b/spheres.hpp
index 161ef9c..dd2a189 100644
--- a/spheres.hpp
+++ b/spheres.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <assert.h>
+#include <random>
 #include <set>
 #include <vector>
 #include <queue>
@@ -224,7 +225,6 @@ public:
   }
 };
 
-
 template <int D, Particle<D> T> class Model {
 private:
   bool metropolisAccept(double β, double ΔE, Rng& r, double a = 1.0) const {
@@ -356,7 +356,7 @@ public:
         }
 
         if (pMax != NULL) {
-          if (1 - 1e2 * exp(-β * ΔEmax) > r.uniform(0.0, 1.0)) {
+          if (1 - exp(-β * ΔEmax) > r.uniform(0.0, 1.0)) {
             m.insert(pMax);
           }
         }
@@ -404,6 +404,58 @@ public:
     return clusterFlip(β, g, m, r);
   }
 
+  unsigned randomCluster(double β, Rng& r) {
+    Vector<2> t = {r.uniform(0.0, 1.0), r.uniform(0.0, 1.0)};
+
+    Matrix<2> mat;
+
+    mat(0, 0) = -1;
+    mat(1, 1) = -1;
+    mat(0, 1) = 0;
+    mat(1, 0) = 0;
+
+    Euclidean<2> g(t - mat * t, mat);
+
+    Move<D, T> m(g, {&r.pick(particles)});
+
+    return clusterFlip(β, g, m, r);
+  }
+
+  unsigned clusterNudge(double β, T& s, Rng& r) {
+    Vector<2> t = s.x + (Vector<2>){r.variate<double, std::normal_distribution>(0, s.r / 100),
+        r.variate<double, std::normal_distribution>(0, s.r / 100)}
+      ;
+
+    Matrix<2> mat;
+
+    mat(0, 0) = -1;
+    mat(1, 1) = -1;
+    mat(0, 1) = 0;
+    mat(1, 0) = 0;
+
+    Euclidean<2> g(t - mat * t, mat);
+
+    Move<D, T> m(g, {&s});
+
+    return clusterFlip(β, g, m, r);
+  }
+
+  unsigned clusterNeighborSwap(double β, T& s, Rng& r) {
+    T* s2 = &s;
+    while (s2 == &s) {
+      s2 = r.pick(nl.neighborsOf(s.x));
+    }
+    return clusterSwap(β, s, *s2, r);
+  }
+
+  unsigned randomClusterNeighborSwap(double β, Rng& r) {
+    return clusterNeighborSwap(β, r.pick(particles), r);
+  }
+
+  unsigned randomClusterNudge(double β, Rng& r) {
+    return clusterNudge(β, r.pick(particles), r);
+  }
+
   unsigned randomClusterSwap(double β, Rng& r) {
     return clusterSwap(β, r.pick(particles), r.pick(particles), r);
   }