1 files changed, 115 insertions, 0 deletions
diff --git a/ising_autocorrelation.cpp b/ising_autocorrelation.cpp
new file mode 100644
index 0000000..053be01
--- /dev/null
+++ b/ising_autocorrelation.cpp
@@ -0,0 +1,115 @@
+#include <iostream>
+#include <fstream>
+
+#include "ising.hpp"
+#include "quantity.hpp"
+
+class Autocorrelation : public measurement<signed, D, TorusGroup<signed, D>, signed> {
+public:
+  Quantity E;
+  Autocorrelation(unsigned lag) : E(lag) {}
+
+  void post_cluster(const isingModel& m) override {
+    double E_tmp = 0;
+
+    TorusGroup<signed, D> s0Inv = m.s0.inverse();
+
+    for (const isingSpin* si : m.s) {
+      E_tmp -= m.B(s0Inv.act(*si));
+
+      for (const isingSpin* sj : m.dict.neighbors(si->x)) {
+        if (si != sj) {
+          E_tmp -= m.Z(*si, *sj) / 2;
+        }
+      }
+    }
+
+    E.add(E_tmp);
+  }
+};
+
+int main(int argc, char* argv[]) {
+  unsigned L = 32;
+  unsigned N = 1000;
+  unsigned mod = 0;
+  unsigned multi = 1e4;
+  double mag = 0.5;
+  double pop = 1.0;
+  double T = 2.0 / log(1.0 + sqrt(2.0));
+  double H = 1.0;
+  double ε = 0.1;
+  unsigned lag = 1e2;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:r:p:l:")) != -1) {
+    switch (opt) {
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'L':
+      L = atoi(optarg);
+      break;
+    case 'T':
+      T = atof(optarg);
+      break;
+    case 'H':
+      H = atof(optarg);
+      break;
+    case 'e':
+      ε = atof(optarg);
+      break;
+    case 'm':
+      mod = atoi(optarg);
+      break;
+    case 'M':
+      multi = atoi(optarg);
+      break;
+    case 'r':
+      mag = atof(optarg);
+      break;
+    case 'p':
+      pop = atof(optarg);
+      break;
+    case 'l':
+      lag = (unsigned)atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  isingModel ising(L, isingZ(L), isingBMod(L, mod, H));
+  isingPopulate(ising, L, pop, mag);
+
+  auto g = isingGen(L);
+
+  measurement<signed, D, TorusGroup<signed, D>, signed> A_empty;
+
+  ising.wolff(T, {g}, A_empty, N);
+
+  Autocorrelation A(lag);
+
+  while (true) {
+    ising.wolff(T, {g}, A, N);
+    std::array<double, 2> τ = A.E.τ();
+    std::cout << A.E.num_added() << " " << A.E.avg() << " " << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n";
+    if (τ[1] / τ[0] < ε && τ[0] * multi < A.E.num_added()) {
+      break;
+    }
+  }
+
+  std::ofstream outfile;
+  outfile.open("out.dat", std::ios::app);
+
+  std::array<double, 2> act = A.E.τ();
+  std::vector<double> ρ = A.E.ρ();
+
+  outfile << L << " " << T << " " << mod << " " << H << " " << A.E.num_added() << " " <<
+  A.E.avg() << " " << A.E.serr() << " " << act[0] << " " << act[1]; for (double ρi : ρ) {
+    outfile << " " << ρi;
+  }
+  outfile << "\n";
+
+  return 0;
+}