1 files changed, 58 insertions, 27 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;