Working Biroli-Mezard clusters.

1 files changed, 84 insertions, 0 deletions

diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
new file mode 100644
index 0000000..ac68776
--- /dev/null
+++ b/biroli-mezard.cpp
@@ -0,0 +1,84 @@
+#include <iostream>
+#include <fstream>
+
+#include "glass.hpp"
+
+void print(const BiroliSystem<2>& s) {
+  for (const Vertex<2, BiroliState>& v : s.vertices) {
+      std::cerr << v.state.type;
+
+    if (v.position(0) == s.L - 1) {
+      std::cerr << std::endl;
+    }
+  }
+}
+
+int main() {
+  const unsigned D = 3;
+  unsigned L = 30;
+  unsigned Nmin = 2e2;
+  unsigned Nmax = 2e5;
+  double Tmin = 0.04;
+  double Tmax = 0.2;
+  double δT = 0.02;
+
+  BiroliSystem<D> s(L);
+
+  Rng r;
+
+  double z = exp(1 / 0.2);
+
+  if (!s.compatible()) {
+    std::cerr << "Storted incompatible!" << std::endl;
+    return 1;
+  } 
+
+  while (s.density() < 0.57) {
+    s.sweepGrandCanonical(z, r);
+  }
+
+  if (!s.compatible()) {
+    std::cerr << "Not compatible!" << std::endl;
+    return 1;
+  } 
+
+  std::cerr << "Found state with appropriate density." << std::endl;
+
+  BiroliSystem<D> s0 = s;
+
+  std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
+
+  std::vector<unsigned> clusterDist(s.size() + 1);
+
+  unsigned n = 1;
+  unsigned i = 0;
+  double nC = 0;
+  while (nC / s.size() < 1e5) {
+    if (n < 20 * log(i + 1)) {
+      n++;
+      std::cout << nC / s.size() << " " << (double)s.overlap(s0) / s.size() << std::endl;
+    }
+    unsigned nn = s.flipCluster(Transformation<D>(L, ms, r), r.pick(s.vertices));
+    nC += nn;
+    clusterDist[nn]++;
+//    s.sweepLocal(r);
+//    nC += s.size();
+//    s.sweepSwap(r);
+//    s.swendsenWang(Transformation<D>(L, ms, r), r);
+    i++;
+  }
+
+  if (!s.compatible()) {
+    std::cerr << "Not compatible!" << std::endl;
+    return 1;
+  } 
+
+  std::ofstream file("dist.dat");
+  for (unsigned i : clusterDist) {
+    file << i << " ";
+  }
+  file.close();
+
+  return 0;
+}
+