1 files changed, 94 insertions, 18 deletions

diff --git a/spheres.hpp b/spheres.hpp
index 6dd5b04..161ef9c 100644
--- a/spheres.hpp
+++ b/spheres.hpp
@@ -177,9 +177,57 @@ public:
   }
 };
 
+template <int D, Particle<D> T> class Move {
+public:
+  const Euclidean<D> g;
+  std::vector<T*> ps;
+  std::vector<Position<D>> xNews;
+
+  Move(const Euclidean<D>& g, std::vector<T*> ps) : g(g), ps(ps) {
+    xNews.reserve(ps.size());
+    for (T* p : ps) {
+      xNews.push_back(g.act((Position<D>)p->x));
+    }
+  }
+
+  void insert(T* p) {
+    ps.push_back(p);
+    xNews.push_back(g.act(p->x));
+  }
+
+  bool anyMarked() const {
+    for (const T* p : ps) {
+      if (p->m) {
+        return true;
+      }
+    }
+
+    return false;
+  }
+
+  unsigned size() const {
+    return ps.size();
+  }
+
+  double ΔE(const T& n) const {
+    double e = 0;
+
+    for (unsigned i = 0; i < ps.size(); i++) {
+      if (&n == ps[i]) {
+        return 0;
+      } else {
+        e += n.U(xNews[i], ps[i]->r) - n.U(ps[i]->x, ps[i]->r);
+      }
+    }
+
+    return e;
+  }
+};
+
+
 template <int D, Particle<D> T> class Model {
 private:
-  bool metropolisAccept(double β, double ΔE, Rng& r) const {
+  bool metropolisAccept(double β, double ΔE, Rng& r, double a = 1.0) const {
     return ΔE < 0 || exp(-β * ΔE) > r.uniform(0.0, 1.0);
   }
 
@@ -208,6 +256,13 @@ public:
     p.x = x;
   }
 
+  void execute(Move<D, T>& m) {
+    for (unsigned i = 0; i < m.size(); i++) {
+      move(*(m.ps[i]), m.xNews[i]);
+      m.ps[i]->m = true;
+    }
+  }
+
   bool move(double β, T& p, const Vector<D>& Δx, Rng& r) {
     Position<D> xNew = p.x + Δx;
     double ΔE = 0;
@@ -270,37 +325,56 @@ public:
     return move(β, r.pick(particles), Δx, r);
   }
 
-  unsigned clusterFlip(double β, const Euclidean<D>& R, T& p0, Rng& r) {
+  unsigned clusterFlip(double β, const Euclidean<D>& R, Move<D, T>& m0, Rng& r) {
     unsigned n = 0;
-    std::queue<T*> q;
-    q.push(&p0);
+    std::queue<Move<D, T>> q;
+    q.push(m0);
 
     while (!q.empty()) {
-      T* p = q.front();
+      Move<D, T> m = q.front();
       q.pop();
 
-      if (!p->m) {
-        p->m = true;
-        Position<D> xNew = R.act(p->x);
+      if (!m.anyMarked()) {
+        std::set<T*> neighbors = {};
+        for (unsigned i = 0; i < m.ps.size(); i++) {
+          std::set<T*> pNeighbors = nl.neighborsOf(m.ps[i]->x);
+          std::set<T*> pNeighborsNew = nl.neighborsOf(m.xNews[i]);
+          neighbors.insert(pNeighbors.begin(), pNeighbors.end());
+          neighbors.insert(pNeighborsNew.begin(), pNeighborsNew.end());
+        }
+        double ΔEmax = 0;
+        T* pMax = NULL;
 
-        for (T* neighbor : nl.neighborsOf(xNew)) {
+        for (T* neighbor : neighbors) {
           if (!neighbor->m) {
-            double ΔE = neighbor->U(xNew, p->r) - neighbor->U(p->x, p->r);
+            double ΔE = m.ΔE(*neighbor);
+            if (ΔE > ΔEmax) {
+              ΔEmax = ΔE;
+              pMax = neighbor;
+            }
+          }
+        }
+
+        if (pMax != NULL) {
+          if (1 - 1e2 * exp(-β * ΔEmax) > r.uniform(0.0, 1.0)) {
+            m.insert(pMax);
+          }
+        }
+        
+        for (T* neighbor : neighbors) {
+          if (!neighbor->m) {
+            double ΔE = m.ΔE(*neighbor);
             if (1 - exp(-β * ΔE) > r.uniform(0.0, 1.0)) {
-              q.push(neighbor);
+              q.push(Move<D, T>(R, {neighbor}));
             }
           }
         }
 
-        move(*p, xNew);
-        n++;
+        execute(m);
+        n += m.size();;
       }
     }
 
-    for (T& p : particles) {
-      p.m = false;
-    }
-
     if (n > particles.size() / 2) {
       orientation = R.act(orientation);
     }
@@ -325,7 +399,9 @@ public:
 
     Euclidean<2> g(t - mat * t, mat);
 
-    return clusterFlip(β, g, s1, r);
+    Move<D, T> m(g, {&s1, &s2});
+
+    return clusterFlip(β, g, m, r);
   }
 
   unsigned randomClusterSwap(double β, Rng& r) {