1 files changed, 27 insertions, 48 deletions

diff --git a/ciamarra.cpp b/ciamarra.cpp
index d3acc38..1d4538f 100644
--- a/ciamarra.cpp
+++ b/ciamarra.cpp
@@ -84,38 +84,6 @@ 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);
-  }
 };
 
 void print(const CiamarraSystem<2>& s) {
@@ -161,27 +129,38 @@ int main() {
     while (s.density() < ρ) {
       s.sweepGrandCanonical(z, r);
     }
-    std::cout << s.density() << " ";
-    s.setInitialPosition();
+
+    std::cout << s.N << std::endl;
+
+    CiamarraSystem<D> s0 = s;
+    double last_overlap = 1;
     auto start = std::chrono::high_resolution_clock::now();
-    for (unsigned i = 0; i < 1e4; i++) {
-      /*
+
+    while (last_overlap > 0.1) {
       for (unsigned j = 0; j < s.vertices.size(); j++) {
-        s.tryRandomNonlocalMove(r);
+        s.tryRandomMove(r);
       }
-      */
+
+      auto stop = std::chrono::high_resolution_clock::now();
+      auto duration = duration_cast<std::chrono::microseconds>(stop - start);
+      last_overlap = ((double)s.overlap(s0) / s.N - s.density() / 6) / (1 - s.density() / 6);
+      std::cout << duration.count() << " " << last_overlap << " ";
+    }
+    std::cout << std::endl;
+
+    CiamarraSystem<D> s1 = s;
+    last_overlap = 1;
+    auto start1 = std::chrono::high_resolution_clock::now();
+
+    while (last_overlap > 0.1) {
       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 << " ";
-      }
-    }
-    if (!s.compatible()) {
-      std::cerr << "Bad configuration" << std::endl;
-      return 1;
+      s.wolff(Transformation<D>(L, m, v.position - m * v.position), r.pick(s.vertices));
+
+      auto stop = std::chrono::high_resolution_clock::now();
+      auto duration = duration_cast<std::chrono::microseconds>(stop - start1);
+      last_overlap = ((double)s.overlap(s1) / s.N - s.density() / 6) / (1 - s.density() / 6);
+      std::cout << duration.count() << " " << last_overlap << " ";
     }
     std::cout << std::endl;
   }