Broken: attempt to implement dynamic matching of pair-flipped spheres.

2 files changed, 46 insertions, 12 deletions

diff --git a/space_wolff.hpp b/space_wolff.hpp
index 9d529a3..141396a 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -62,7 +62,7 @@ public:
         for (Spin<U, D, S>* s : t->toConsider()) {
           double ΔE = t->ΔE(s);
           if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform<double>(0, 1)) {
-            queue.push(t->createNew(s));
+            queue.push(t->createNew(T, visited, s, rng));
           }
         }
 
diff --git a/transformation.hpp b/transformation.hpp
index 0731ec9..fc1833e 100644
--- a/transformation.hpp
+++ b/transformation.hpp
@@ -12,7 +12,7 @@ public:
   virtual std::set<Spin<U, D, S>*> current() const { return {NULL}; }
   virtual std::set<Spin<U, D, S>*> toConsider() const { return {}; }
   virtual double ΔE(const Spin<U, D, S>*) const { return 0; }
-  virtual Transformation* createNew(Spin<U, D, S>*) const { return new Transformation(); }
+  virtual Transformation* createNew(double T, const std::set<Spin<U, D, S>*>&, Spin<U, D, S>*, Rng& rng) const { return new Transformation(); }
   virtual void apply(){};
 };
 
@@ -56,7 +56,7 @@ public:
     }
   }
 
-  Transformation<U, D, R, S>* createNew(Spin<U, D, S>* s) const override;
+  Transformation<U, D, R, S>* createNew(double, const std::set<Spin<U, D, S>*>&, Spin<U, D, S>* s, Rng&) const override;
 
   void apply() override {
     M.dict.remove(sOld);
@@ -112,8 +112,8 @@ public:
     }
   }
 
-  Transformation<U, D, R, S>* createNew(Spin<double, D, S>* s) const;
-  Transformation<U, D, R, S>* createNew(Spin<signed, D, S>* s) const;
+  Transformation<U, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const;
+  Transformation<U, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<signed, D, S>* s, Rng&) const;
 
   void apply() override {
     M.dict.remove(s1Old);
@@ -150,11 +150,11 @@ public:
     return M.B(s0s_new) - M.B(s0s_old);
   }
 
-  Transformation<double, D, R, S>* createNew(Spin<double, D, S>* s) const {
+  Transformation<double, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const {
     return new SpinFlip<U, D, R, S>(M, r, s);
   }
 
-  Transformation<signed, D, R, S>* createNew(Spin<signed, D, S>* s) const {
+  Transformation<signed, D, R, S>* createNew(double, const std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
     Vector<signed, 2> v = r.act(*s).x;
     std::set<Spin<U, D, S>*> on_site = M.dict.at(v);
     if (on_site.empty()) {
@@ -168,25 +168,59 @@ public:
 };
 
 template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(Spin<U, D, S>* s) const {
+Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<U, D, S>* s, Rng& rng) const {
   if (s == NULL) {
     return new FieldFlip<U, D, R, S>(M, r);
   } else {
-    return new SpinFlip(M, r, s);
+    SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, r, s);
+    Spin<U, D, S>* sMax = NULL;
+    double ΔEMax = 0.0;
+    for (Spin<U, D, S>* ss : t->toConsider()) {
+      if (ss != NULL && ss != s && !v.contains(ss)) {
+        double ΔE = t->ΔE(ss);
+        if (ΔE > ΔEMax) {
+          sMax = ss;
+          ΔEMax = ΔE;
+        }
+      }
+    }
+    if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1)) {
+      return t;
+    } else {
+      delete t;
+      return new PairFlip<U, D, R, S>(M, r, s, sMax);
+    }
   }
 }
 
 template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<double, D, S>* s) const {
+Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng& rng) const {
   if (s == NULL) {
     return new FieldFlip<U, D, R, S>(M, r);
   } else {
-    return new SpinFlip<U, D, R, S>(M, r, s);
+    SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, r, s);
+    Spin<U, D, S>* sMax = NULL;
+    double ΔEMax = 0.0;
+    for (Spin<U, D, S>* ss : t->toConsider()) {
+      if (ss != NULL && ss != s && !v.contains(ss)) {
+        double ΔE = t->ΔE(ss);
+        if (ΔE > ΔEMax) {
+          sMax = ss;
+          ΔEMax = ΔE;
+        }
+      }
+    }
+    if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1)) {
+      return t;
+    } else {
+      delete t;
+      return new PairFlip<U, D, R, S>(M, r, s, sMax);
+    }
   }
 }
 
 template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<signed, D, S>* s) const {
+Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double, const std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
   if (s == NULL) {
     return new FieldFlip<U, D, R, S>(M, r);
   } else {