Work on the Ciamarra model.

3 files changed, 167 insertions, 53 deletions

diff --git a/ciamarra.cpp b/ciamarra.cpp
index 7d7c8af..d1087a4 100644
--- a/ciamarra.cpp
+++ b/ciamarra.cpp
@@ -36,6 +36,7 @@ public:
 
 template <unsigned D>
 class CiamarraSystem : public HardSystem<D, CiamarraState<D>> {
+  public:
   using HardSystem<D, CiamarraState<D>>::HardSystem;
 
   std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>>
@@ -83,6 +84,38 @@ class CiamarraSystem : public HardSystem<D, CiamarraState<D>> {
       }
     }
   }
+
+  /* For the Ciamarra system, position within sites must be specially acconted
+   * for in the scattering function. We therefore expand the lattice and add
+   * the state when setting positions, so as later to be able to take particle
+   * sublattice positions into account.
+   * */
+  void setInitialPosition() override {
+    CiamarraSystem::orientation = Transformation<D>(CiamarraSystem::L);
+    for (Vertex<D, CiamarraState<D>>& v : CiamarraSystem::vertices) {
+      v.initialPosition = 4 * v.position + v.state;
+    }
+  }
+
+  double selfIntermediateScattering(const std::vector<Matrix<D>>& ms) const override {
+    double F = 0;
+
+    std::array<Vector<D>, D> ks;
+    for (unsigned i = 0; i < D; i++) {
+      ks[i].setZero();
+      ks[i](i) = 1;
+    }
+    for (const Vertex<D, CiamarraState<D>>& v : CiamarraSystem::vertices) {
+      if (!v.empty()) {
+        Vector<D> Δx = ((4 * CiamarraSystem::orientation.inverse().apply(v.position)) + CiamarraSystem::orientation.inverse().apply(v.state)) - v.initialPosition;
+        for (const Vector<D>& k : ks) {
+          F += cos((M_PI / 4) * k.dot(Δx));
+        }
+      }
+    }
+
+    return F / (D * CiamarraSystem::N);
+  }
 };
 
 
@@ -123,37 +156,35 @@ int main() {
 
   double z = exp(1 / 0.08);
 
-  while (s.density() < 0.818) {
-    s.sweepGrandCanonical(z, r);
-  }
-
-  if (!s.compatible()) {
-    std::cerr << "Net compatible!" << std::endl;
-    return 1;
-  } 
-
-  std::cerr << "Found state with appropriate density." << std::endl;
-
-  CiamarraSystem<D> s0 = s;
-
   std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
 
-  unsigned n = 1;
-  for (unsigned i = 0; i < 1e5; i++) {
-    if (n < 20 * log(i + 1)) {
-      n++;
-      std::cout << i << " " << s.overlap(s0) << std::endl;
+  for (double ρ : {0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.76, 0.78, 0.8, 0.81, 0.815, 0.818, 0.821, 0.825}) {
+    while (s.density() < ρ) {
+      s.sweepGrandCanonical(z, r);
+    }
+    std::cout << s.density() << " ";
+    s.setInitialPosition();
+    auto start = std::chrono::high_resolution_clock::now();
+    for (unsigned i = 0; i < 1e4; i++) {
+      /*
+      for (unsigned j = 0; j < s.vertices.size(); j++) {
+        s.tryRandomNonlocalMove(r);
+      }
+      */
+      Matrix<D> m = r.pick(ms);
+      Vertex<D, CiamarraState<D>>& v = r.pick(s.vertices);
+      unsigned n = s.wolff(Transformation<D>(L, m, v.position - m * v.position), r);
+      if (i % 1 == 0) {
+        auto stop = std::chrono::high_resolution_clock::now();
+        auto duration = duration_cast<std::chrono::microseconds>(stop - start);
+        std::cout << duration.count() << " " << s.selfIntermediateScattering(ms) << " " << n << " ";
+      }
     }
-    Matrix<D> m = r.pick(ms);
-    Vertex<D, CiamarraState<D>>& v = r.pick(s.vertices);
-    unsigned nC = s.flipCluster(Transformation<D>(L, m, v.position - m * v.position), v, r);
-    std::cerr << nC << std::endl;
-    for (Vertex<D, CiamarraState<D>>& v : s.vertices) {
-      v.marked = false;
+    if (!s.compatible()) {
+      std::cerr << "Bad configuration" << std::endl;
+      return 1;
     }
-//    s.sweepLocal(r);
-//    s.sweepSwap(r);
-//    s.swendsenWang(Transformation<D>(L, ms, r), r);
+    std::cout << std::endl;
   }
 
   return 0;
diff --git a/distinguishable.cpp b/distinguishable.cpp
index 22f228a..b739927 100644
--- a/distinguishable.cpp
+++ b/distinguishable.cpp
@@ -83,13 +83,23 @@ int main() {
     */
     /*
      */
-    s.setInitialPosition();
     std::cout << T << " ";
-    for (unsigned i = 0; i < 1e4; i++) {
+    s.setInitialPosition();
+    auto start = std::chrono::high_resolution_clock::now();
+    for (unsigned i = 0; i < 1e5; i++) {
+      for (unsigned j = 0; j < s.vertices.size(); j++) {
+        s.tryRandomSwap(T, r);
+      }
+      /*
       Matrix<D> m = r.pick(ms);
       Vertex<D, DistinguishableState>& v = r.pick(s.vertices);
       s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
-      std::cout << s.selfIntermediateScattering() << " ";
+      */
+      if (i % 10 == 0) {
+        auto stop = std::chrono::high_resolution_clock::now();
+        auto duration = duration_cast<std::chrono::microseconds>(stop - start);
+        std::cout << duration.count() << " " << s.selfIntermediateScattering(ms) << " ";
+      }
     }
     /*
     for (unsigned i = 0; i < 1e5; i++) {
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) {