Work on the Ciamarra model.

1 files changed, 96 insertions, 23 deletions

diff --git a/glass.hpp b/glass.hpp
index fc7b186..8ad0b52 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -111,6 +111,8 @@ template <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() {
   return mats;
 }
 
+template <unsigned D> class CiamarraState;
+
 template <unsigned D> class Transformation {
 public:
   unsigned L;
@@ -146,6 +148,11 @@ public:
   Transformation<D> inverse() const {
     return Transformation<D>(L, m.transpose(), -m.transpose() * v);
   }
+
+  CiamarraState<D> apply(const CiamarraState<D>& s) const {
+    CiamarraState<D> s2(m * s);
+    return s2;
+  }
 };
 
 template <typename T> concept State = requires(T v, const T& v2) {
@@ -225,25 +232,31 @@ public:
     return iPow(L, D);
   }
 
-  void setInitialPosition() {
+  virtual void setInitialPosition() {
+    orientation = Transformation<D>(L);
     for (Vertex<D, S>& v : vertices) {
       v.initialPosition = v.position;
     }
   }
 
-  double selfIntermediateScattering() const {
+  virtual double selfIntermediateScattering(const std::vector<Matrix<D>>& ms) const {
     double F = 0;
 
-    Vector<D> k1 = {M_PI, 0};
-    Vector<D> k2 = {0, M_PI};
+    std::array<Vector<D>, D> ks;
+    for (unsigned i = 0; i < D; i++) {
+      ks[i].setZero();
+      ks[i](i) = 1;
+    }
     for (const Vertex<D, S>& v : vertices) {
       if (!v.empty()) {
-        F += cos(k1.dot(v.position - v.initialPosition));
-        F += cos(k2.dot(v.position - v.initialPosition));
+        Vector<D> Δx = v.position - orientation.apply(v.initialPosition);
+        for (const Vector<D>& k : ks) {
+          F += cos(M_PI * k.dot(Δx));
+        }
       }
     }
 
-    return F / (2 * this->N);
+    return F / (D * System::N);
   }
 };
 
@@ -371,7 +384,13 @@ public:
           std::swap(v.initialPosition, vNew.initialPosition);
         }
 
-        n += 1;
+        if (!v.empty()) {
+          n += 1;
+        }
+
+        if (!vNew.empty()) {
+          n += 1;
+        }
       }
     }
 
@@ -380,6 +399,9 @@ public:
 
   unsigned wolff(const Transformation<D>& R, double T, Rng& r) {
     unsigned n = flipCluster(R, r.pick(this->vertices), T, r);
+    if (n > SoftSystem::N / 2) {
+      SoftSystem::orientation = R.apply(SoftSystem::orientation);
+    }
     for (Vertex<D, S>& v : this->vertices) {
       v.marked = false;
     }
@@ -431,6 +453,7 @@ public:
   bool tryRandomMove(Rng& r) {
     Vertex<D, S>& v = r.pick(HardSystem::vertices);
     S oldState = v.state;
+    Vector<D> oldPos = v.initialPosition;
 
     if (!tryDeletion(v)) {
       return false;
@@ -440,6 +463,7 @@ public:
       for (HalfEdge<D, S>& e : v.adjacentEdges) {
         if (1 == e.Δx.dot(oldState)) {
           if (insert(e.neighbor, oldState.flip())) {
+            e.neighbor.initialPosition = oldPos;
             return true;
           }
           break;
@@ -451,17 +475,46 @@ public:
         newState = S(r);
       }
       if (insert(v, newState)) {
+        v.initialPosition = oldPos;
         return true;
       }
     }
+    v.initialPosition = oldPos;
+    v.state = oldState;
+    HardSystem::N++;
+    return false;
+  }
+
+  bool tryRandomNonlocalMove(Rng& r) {
+    Vertex<D, S>& v = r.pick(HardSystem::vertices);
+    S oldState = v.state;
+    Vector<D> oldPos = v.initialPosition;
+
+    if (!tryDeletion(v)) {
+      return false;
+    }
+
+    Vertex<D, S>& vNew = r.pick(HardSystem::vertices);
+    S newState(r);
+
+    if (insert(vNew, newState)) {
+      vNew.initialPosition = oldPos;
+      return true;
+    }
+
+    v.initialPosition = oldPos;
     v.state = oldState;
     HardSystem::N++;
     return false;
   }
 
   bool trySwap(Vertex<D, S>& v1, Vertex<D, S>& v2) {
+    if (v1.state == v2.state) {
+      return false;
+    }
     if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) {
       std::swap(v1.state, v2.state);
+      std::swap(v1.initialPosition, v2.initialPosition);
       return true;
     } else {
       return false;
@@ -539,34 +592,54 @@ public:
         Vector<D> xNew = R.apply(v.position);
         Vertex<D, S>& vNew = HardSystem::vertices[HardSystem::vectorToIndex(xNew)];
 
-        v.marked = true;
-        vNew.marked = true;
-
         S s = R.apply(v.state);
         S sNew = R.apply(vNew.state);
 
-        std::list<std::reference_wrapper<Vertex<D, S>>> overlaps1 = overlaps(vNew, s, true);
-        std::list<std::reference_wrapper<Vertex<D, S>>> overlaps2 = overlaps(v, sNew, true);
-        overlaps1.splice(overlaps1.begin(), overlaps2);
+        if (s != vNew.state) {
+          v.marked = true;
+          vNew.marked = true;
+
+          std::list<std::reference_wrapper<Vertex<D, S>>> overlaps1 = overlaps(vNew, s, true);
+          std::list<std::reference_wrapper<Vertex<D, S>>> overlaps2 = overlaps(v, sNew, true);
+          overlaps1.splice(overlaps1.begin(), overlaps2);
 
-        for (Vertex<D, S>& vn : overlaps1) {
-          if (!vn.marked) {
-            q.push(vn);
+          for (Vertex<D, S>& vn : overlaps1) {
+            if (!vn.marked) {
+              q.push(vn);
+            }
           }
-        }
 
-        if (!dry) {
-          v.state = sNew;
-          vNew.state = s;
-        }
+          if (!dry) {
+            v.state = sNew;
+            vNew.state = s;
+            std::swap(v.initialPosition, vNew.initialPosition);
+          }
+
+          if (!v.empty()) {
+            n += 1;
+          }
 
-        n += 1;
+          if (!vNew.empty()) {
+            n += 1;
+          }
+        }
       }
     }
 
     return n;
   }
 
+  unsigned wolff(const Transformation<D>& R, Rng& r) {
+    unsigned n = flipCluster(R, r.pick(this->vertices), r);
+    if (n > HardSystem::N / 2) {
+      HardSystem::orientation = R.apply(HardSystem::orientation);
+    }
+    for (Vertex<D, S>& v : this->vertices) {
+      v.marked = false;
+    }
+    return n;
+  }
+
   void swendsenWang(const Transformation<D>& R, Rng& r) {
     for (Vertex<D, S>& v : HardSystem::vertices) {
       if (!v.marked) {