nnn using new wolff edge properties

2 files changed, 242 insertions, 265 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3ddd0d6..3363340 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -13,5 +13,5 @@ link_directories(~/.local/lib64)
 
 add_executable(new_model src/new_model.cpp)
 target_include_directories(new_model PUBLIC ~/.local/include)
-target_link_libraries(new_model wolff)
+target_link_libraries(new_model)
 
diff --git a/src/new_model.cpp b/src/new_model.cpp
index 96789d1..b96421d 100644
--- a/src/new_model.cpp
+++ b/src/new_model.cpp
@@ -1,375 +1,352 @@
 
 #include <getopt.h>
 #include <iostream>
+#include <fstream>
 
 #include "../include/randutils/randutils.hpp"
 
 #include <functional>
 #include <vector>
 #include <queue>
-
-#include <wolff/types.h>
-#include <wolff/graph.hpp>
-
-using namespace wolff;
+#include <list>
 
 namespace ising {
 #define WOLFF_SITE_DEPENDENCE
-#include <wolff.hpp>
-#include <wolff/models/ising.hpp>
+#include <wolff_models/ising.hpp> // ising includes wolff
 #undef WOLFF_SITE_DEPEDENCE
+  using namespace wolff;
 }
 
 namespace xy {
 #undef WOLFF_H
-#undef WOLFF_SYSTEM_H
-#undef WOLFF_CLUSTER_H
-#undef WOLFF_MEASUREMENTS_H
 
 #define WOLFF_NO_FIELD
 #define WOLFF_BOND_DEPENDENCE
-#include <wolff.hpp>
-#include <wolff/models/vector.hpp>
-#include <wolff/models/orthogonal.hpp>
+#include <wolff_models/vector.hpp>
+#include <wolff_models/orthogonal.hpp> // orthogonal includes wolff
 #undef WOLFF_NO_FIELD
 #undef WOLFF_BOND_DEPENDENCE
+  using namespace wolff;
 }
 
-double E;
+using ising::ising_t;
+using xy::vector_t;
+using xy::orthogonal_t;
+
+// at each ising site, we need to know which sublattice we're in and the xy
+// spins to our left and right (top and bottom)
+typedef struct _ising_site_props {
+  unsigned sublattice;
+  vector_t<2, double>* back;
+  vector_t<2, double>* front;
+} ising_site_props;
+
+// at each xy bond, we need to know whether we're nearest or next-nearest, the
+// displacement the bond makes, and the ising spin that lies on us (if any)
+typedef struct _xy_bond_props {
+  unsigned distance;
+  vector_t<2, int> dNeighborSelf;
+  ising_t* spin;
+} xy_bond_props;
+
+typedef ising::graph<ising_site_props> graph_i;
+typedef xy::graph<std::tuple<>, xy_bond_props> graph_x;
+
+double vsin(const vector_t<2, double>& v1, const vector_t<2, double>& v2) {
+  return v1[0] * v2[1] - v1[1] * v2[0];
+}
 
-using ising::wolff::ising_t;
-using xy::wolff::vector_t;
-using xy::wolff::orthogonal_t;
+double cos_xy;
+vector_t<2, double> sin_xy;
 
-class i_measurement : public ising::wolff::measurement<ising_t, ising_t> {
+class i_measurement : public ising::wolff::measurement<ising_t, ising_t, graph_i> {
   private:
-    N_t nc;
-    N_t ns;
+    double a;
 
-    int M;
+    unsigned n;
     int A;
-    v_t C;
-
-    double totalaA;
-    double totalA2;
-    double totalaM;
-    double totalM2;
-    double totalC;
-    double totalC2;
 
   public:
-    i_measurement(const ising::wolff::system<ising_t, ising_t>& S) {
-      nc = 0;
-      ns = 0;
-      M = 0;
+    double totalA2;
+    i_measurement(const ising::wolff::system<ising_t, ising_t, graph_i>& S, double a) : a(a) {
       A = S.nv;
 
-      totalaA = 0;
       totalA2 = 0;
-      totalaM = 0;
-      totalM2 = 0;
-      totalC = 0;
-      totalC2 = 0;
     }
 
-    void pre_cluster(N_t, N_t, const ising::wolff::system<ising_t, ising_t>&, v_t, const ising_t&) override {
-      C = 0;
-    }
-
-    void plain_bond_visited(const ising::wolff::system<ising_t, ising_t>&, v_t, const ising_t&, v_t, double dE) override {
-      E += dE;
-    }
-
-    void ghost_bond_visited(const ising::wolff::system<ising_t, ising_t>& S, v_t i, const ising_t& s_old, const ising_t& s_new, double dE) override {
-      E += dE;
-      M += s_new - s_old;
-
-      if (i >= S.nv / 2) {
-        A -= s_new - s_old;
-      } else {
+    void ghost_bond_visited(const ising::wolff::system<ising_t, ising_t, graph_i>& S, const graph_i::vertex& v, const ising_t& s_old, const ising_t& s_new, double dE) override {
+      if (v.prop.sublattice == 0) {
         A += s_new - s_old;
+      } else {
+        A -= s_new - s_old;
       }
-    }
-
-    void plain_site_transformed(const ising::wolff::system<ising_t, ising_t>& S, v_t i, const ising_t& si_new) override {
-      C++;
-    }
-
-    void post_cluster(N_t, N_t, const ising::wolff::system<ising_t, ising_t>&) override {
-      totalaA += abs(A);
-      totalA2 += pow(A, 2);
-      totalaM += abs(M);
-      totalM2 += pow(M, 2);
-      totalC += C;
-      totalC2 += pow(C, 2);
-
-      nc++;
-
-      ns++;
-    }
-
-    double avgaA() {
-      return totalaA / ns;
-    }
 
-    double avgA2() {
-      return totalA2 / ns;
-    }
-
-    double avgaM() {
-      return totalaM / ns;
-    }
+      sin_xy[v.prop.sublattice] += a * (s_new - s_old) * vsin(*(v.prop.back), *(v.prop.front));
 
-    double avgM2() {
-      return totalM2 / ns;
+      cos_xy -= dE;
     }
 
-    double avgC() {
-      return totalC / nc;
-    }
-
-    double avgC2() {
-      return totalC2 / nc;
+    void post_cluster(unsigned, unsigned, const ising::wolff::system<ising_t, ising_t, graph_i>&) override {
+      totalA2 += pow(A, 2);
     }
 };
 
-class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, vector_t<2, double>> {
+class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, vector_t<2, double>, graph_x> {
   private:
-    N_t ns;
-    N_t nc;
-    orthogonal_t<2, double> R;
-
-    vector_t<2, double> M;
-    v_t C;
-    vector_t<2, double> CM;
-
-    vector_t<2, double> totalaM;
-    vector_t<2, double> totalM2;
-    double totalM;
-    double totalCM;
-    double totalC;
-    double totalC2;
+    double Jxy;
+    double Jxyn;
+    double a;
+    const ising_t* si0;
 
   public:
-    x_measurement(const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>>& S) {
-      ns = 0;
-      nc = 0;
-      M = S.nv * S.s[0];
-
-      totalaM.fill(0);
-      totalM2.fill(0);
-      totalC = 0;
-      totalC2 = 0;
-      totalM = 0;
-      totalCM = 0;
-    }
-
-    void pre_cluster(N_t, N_t, const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>>&, v_t, const orthogonal_t<2, double>& r) override {
-      R = r;
-      C = 0;
-      CM.fill(0);
-    }
-
-    void plain_bond_visited(const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>>&, v_t, const vector_t<2, double>&, v_t, double dE) override {
-      E += dE;
-    }
-
-    void plain_site_transformed(const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>>& S, v_t i, const vector_t<2, double>& si_new) override {
-      C++;
-      CM += S.s[i];
-      M += si_new - S.s[i];
-    }
-
-    void post_cluster(N_t, N_t, const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>>&) override {
-      totalaM[0] += fabs(M[0]);
-      totalaM[1] += fabs(M[1]);
-      totalM2[0] += pow(M[0], 2); 
-      totalM2[1] += pow(M[1], 2); 
-      totalM += (pow(M[0], 2) + pow(M[1], 2));
-      totalCM += pow(CM[0] * R[0][0] + CM[1] * R[0][1], 2) / C;
-      totalC += C;
-      totalC2 += pow(C, 2);
-
-      nc++;
-
-      ns++;
-    }
-
-    vector_t<2, double> avgaM() {
-      return totalaM / ns;
-    }
-
-    vector_t<2, double> avgM2() {
-      return totalM2 / ns;
+    x_measurement(const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>, graph_x>& S, double Jxy, double Jxyn, double a, const ising_t *si0) : Jxy(Jxy), Jxyn(Jxyn), a(a), si0(si0) {
     }
 
-    double avgM() {
-      return totalM / ns;
-    }
+    void plain_bond_visited(const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>, graph_x>& S, const graph_x::halfedge& e, const vector_t<2, double>& s1_new, double dE) override {
+      cos_xy -= dE;
 
-    double avgCM() {
-      return totalCM / ns;
-    }
+      vector_t<2, double> dsin = (vsin(s1_new, S.s[e.neighbor.ind]) - vsin(S.s[e.self.ind], S.s[e.neighbor.ind])) * e.prop.dNeighborSelf;
 
-    double avgC() {
-      return totalC / nc;
-    }
-
-    double avgC2() {
-      return totalC2 / nc;
+      if (e.prop.distance == 0) {
+        sin_xy += (Jxy + (*si0) * a * *(e.prop.spin)) * dsin;
+      } else if (e.prop.distance == 1) {
+        sin_xy += Jxyn * dsin;
+      }
     }
 };
 
-using namespace wolff;
-
 int main (int argc, char *argv[]) {
 
   // set defaults
-  N_t N = (N_t)1e4;
-  const D_t D = 2;
-  L_t L = 128;
-  double T = 2.26918531421;
-  double J = 1.0;
+  unsigned N = (unsigned)1e4;
+  const unsigned D = 2;
+  unsigned L = 128;
+  double T = 1.0;
+  double Ji = 1.0;
+  double Jx = 1.0;
+  double Jxn = 0.0;
   double a = 1.5;
 
   int opt;
 
   // take command line arguments
-  while ((opt = getopt(argc, argv, "N:L:T:J:a:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:L:T:I:X:a:n:")) != -1) {
     switch (opt) {
-    case 'N': // number of steps
-      N = (N_t)atof(optarg);
-      break;
-    case 'L': // linear size
-      L = atoi(optarg);
-      break;
-    case 'T': // temperature
-      T = atof(optarg);
-      break;
-    case 'J': // temperature
-      J = atof(optarg);
-      break;
-    case 'a': // temperature
-      a = atof(optarg);
-      break;
-    default:
-      exit(EXIT_FAILURE);
+      case 'N': // number of steps
+        N = (unsigned)atof(optarg);
+        break;
+      case 'L': // linear size
+        L = atoi(optarg);
+        break;
+      case 'T': // temperature
+        T = atof(optarg);
+        break;
+      case 'I': // temperature
+        Ji = atof(optarg);
+        break;
+      case 'X': // temperature
+        Jx = atof(optarg);
+        break;
+      case 'a': // temperature
+        a = atof(optarg);
+        break;
+      case 'n': // temperature
+        Jxn = atof(optarg);
+        break;
+      default:
+        exit(EXIT_FAILURE);
     }
   }
 
+  // initialize random numbers
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
 
+  // antiferromagnetic Ising coupling
   std::function <double(const ising_t&, const ising_t&)> Zi = [=] (const ising_t& s1, const ising_t s2) -> double {
     if (s1.x == s2.x) {
-      return -J;
+      return -Ji;
     } else {
-      return J;
+      return Ji;
     }
   };
 
-  std::vector<vector_t<2, double>> *xy_spins;
-
-  std::function <double(v_t, const ising_t&)> Bi = [L, D, a, &xy_spins] (v_t v, const ising_t& s) -> double {
-    v_t xyv1 = v % (v_t)pow(L, D);
-    v_t xyv2;
-    if (v / (v_t)pow(L, D) == 0) {
-      xyv2 = L * (xyv1 / L) + ((xyv1 % L) + 1) % L;
-    } else {
-      xyv2 = L * (((xyv1 / L) + 1) % L) + (xyv1 % L);
-    }
+  // ising field depends on dot product of surrounding xy spins
+  std::function <double(const graph_i::vertex&, const ising_t&)> Bi =
+    [a] (const graph_i::vertex& v, const ising_t& s) -> double {
+      return a * s * (*(v.prop.back) * *(v.prop.front));
+    };
 
-    double value = a * ((xy_spins->at(xyv1)) * (xy_spins->at(xyv2)));
+  // initialize Ising diagonal-lattice graph
+  graph_i Gi;
 
-    if (s.x) {
-      return -value;
-    } else {
-      return value;
-    }
-  };
+  Gi.L = L;
+  Gi.D = 2;
 
-  graph Gi;
-  
   Gi.nv = 2 * pow(L, 2);
   Gi.ne = 2 * Gi.nv;
 
-  Gi.adjacency.resize(Gi.nv);
-  Gi.coordinate.resize(Gi.nv);
+  Gi.vertices.resize(Gi.nv);
 
-  for (std::vector<v_t> adj_i : Gi.adjacency) {
-    adj_i.reserve(4);
-  }
+  for (unsigned i = 0; i < 2; i++) {
+    unsigned sb = i * pow(L, 2);
 
-  for (D_t i = 0; i < 2; i++) {
-    v_t sb = i * pow(L, 2);
+    for (unsigned j = 0; j < pow(L, 2); j++) {
+      unsigned vc = sb + j;
+      Gi.vertices[vc].ind = vc;
+      Gi.vertices[vc].prop.sublattice = i;
 
-    for (v_t j = 0; j < pow(L, 2); j++) {
-      v_t vc = sb + j;
+      graph_i::halfedge e1(Gi.vertices[vc], Gi.vertices[((i + 1) % 2) * pow(L, 2) + j]);
+      graph_i::halfedge e2(Gi.vertices[vc], Gi.vertices[((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * (j / (unsigned)pow(L, 2 - 1)) + (j + 1 - 2 * (i % 2)) % L]);
+      graph_i::halfedge e3(Gi.vertices[vc], Gi.vertices[((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * ((L + (j/ (unsigned)pow(L, 2 - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L]);
+      graph_i::halfedge e4(Gi.vertices[vc], Gi.vertices[((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * ((L + (j/ (unsigned)pow(L, 2 - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L]);
 
-      Gi.adjacency[vc].push_back(((i + 1) % 2) * pow(L, 2) + j);
-      Gi.adjacency[vc].push_back(((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * (j / (v_t)pow(L, 2 - 1)) + (j + 1 - 2 * (i % 2)) % L);
-      Gi.adjacency[vc].push_back(((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * ((L + (j/ (v_t)pow(L, 2 - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L);
-      Gi.adjacency[vc].push_back(((i + 1) % 2) * pow(L, 2) + pow(L, 2 - 1) * ((L + (j/ (v_t)pow(L, 2 - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L);
+      Gi.vertices[vc].edges.push_back(e1);
+      Gi.vertices[vc].edges.push_back(e2);
+      Gi.vertices[vc].edges.push_back(e3);
+      Gi.vertices[vc].edges.push_back(e4);
     }
   }
-  
-  ising::wolff::system<ising_t, ising_t> si(Gi, T, Zi, Bi);
 
-  std::function <double(v_t, const vector_t<2, double>&, v_t, const vector_t<2, double>&)> Zxy = [L, D, a, &si] (v_t v1, const vector_t<2, double>& s1, v_t v2, const vector_t<2, double>& s2) -> double {
-    v_t iv;
-    if (v1 / L == v2 / L) {
-      v_t vs = v1 < v2 ? v1 : v2;
-      v_t vl = v1 < v2 ? v2 : v1;
+  ising::wolff::system<ising_t, ising_t, graph_i> si(Gi, T, Zi, Bi);
 
-      if (vl == L - 1 && vs == 0) {
-        iv = vl;
+  // ferromagnetic XY coupling, with nearest-neighbor energy that deponds on
+  // state of Ising spin on the bond
+  std::function <double(const graph_x::halfedge&, const vector_t<2, double>&, const vector_t<2, double>&)>
+    Zxy = [Jx, Jxn, a, &si] (const graph_x::halfedge& e, const vector_t<2, double>& s1, const vector_t<2, double>& s2) -> double {
+      if (e.prop.distance == 0) {
+        return (Jx + si.s0 * a * *(e.prop.spin)) * (s1 * s2);
+      } else if (e.prop.distance == 1) {
+        return Jxn * (s1 * s2);
       } else {
-        iv = vs;
+        return 0;
       }
-    } else {
-      v_t vs = v1 < v2 ? v1 : v2;
-      v_t vl = v1 < v2 ? v2 : v1;
+  };
+
+  // initialize square-lattice xy graph
+  graph_x Gxy(2, L);
+
+  // assign edge properties to the XY nearest-neighbor bonds
+  for (graph_x::vertex& v : Gxy.vertices) {
+    for (graph_x::halfedge& e : v.edges) {
+      e.prop.distance = 0;
+
+      int v1 = e.self.ind;
+      int v2 = e.neighbor.ind;
+
+      e.prop.dNeighborSelf[0] = (v2 % L) - (v1 % L);
+      e.prop.dNeighborSelf[1] = v2 / L - v1 / L;
 
-      if (vl / L == L - 1 && vs / L == 0) {
-        iv = vl;
+      for (unsigned i = 0; i < 2; i++) {
+        if (e.prop.dNeighborSelf[i] == L - 1) {
+          e.prop.dNeighborSelf[i] = -1;
+        } 
+        if (e.prop.dNeighborSelf[i] == -(L - 1)) {
+          e.prop.dNeighborSelf[i] = 1;
+        }
+      }
+
+      unsigned vs = v1 < v2 ? v1 : v2;
+      unsigned vl = v1 < v2 ? v2 : v1;
+
+      if (v1 / L == v2 / L) {
+        if (vl % L == L - 1 && vs % L == 0) {
+          e.prop.spin = &(si.s[vl]);
+        } else {
+          e.prop.spin = &(si.s[vs]);
+        }
       } else {
-        iv = vs;
+        if (vl / L == L - 1 && vs / L == 0) {
+          e.prop.spin = &(si.s[pow(L, 2) + vl]);
+        } else {
+          e.prop.spin = &(si.s[pow(L, 2) + vs]);
+        }
       }
     }
+  }
 
-    if (si.s[iv].x) {
-      return (1 + a) * (s1 * s2);
-    } else {
-      return (1 - a) * (s1 * s2);
-    }
-  };
+  // add the next-nearest-neighbor bonds
+  for (graph_x::vertex &v : Gxy.vertices) {
+    int i = v.ind / L;
+    int j = v.ind % L;
+
+    graph_x::halfedge e1(v, Gxy.vertices[((i - 1) % L) * L + (j + 1) % L]);
+    graph_x::halfedge e2(v, Gxy.vertices[((i - 1) % L) * L + (j - 1) % L]);
+    graph_x::halfedge e3(v, Gxy.vertices[((i + 1) % L) * L + (j + 1) % L]);
+    graph_x::halfedge e4(v, Gxy.vertices[((i + 1) % L) * L + (j - 1) % L]);
+
+    e1.prop.distance = 1;
+    e2.prop.distance = 1;
+    e3.prop.distance = 1;
+    e4.prop.distance = 1;
+
+    e1.prop.dNeighborSelf[0] =  1;
+    e1.prop.dNeighborSelf[1] = -1;
+
+    e2.prop.dNeighborSelf[0] = -1;
+    e2.prop.dNeighborSelf[1] = -1;
 
-  graph Gxy(2, L);
+    e3.prop.dNeighborSelf[0] =  1;
+    e3.prop.dNeighborSelf[1] =  1;
 
-  xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>> sxy(Gxy, T, Zxy);
+    e4.prop.dNeighborSelf[0] = -1;
+    e4.prop.dNeighborSelf[1] =  1;
 
-  xy_spins = &(sxy.s);
+    v.edges.push_back(e1);
+    v.edges.push_back(e2);
+    v.edges.push_back(e3);
+    v.edges.push_back(e4);
+  }
+
+  xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>, graph_x> sxy(Gxy, T, Zxy);
+
+  // now that the xy spins have been initialized, go back to the ising graph
+  // and add a reference to the ones on each bond
+  for (unsigned i = 0; i < 2; i++) {
+    unsigned sb = i * pow(L, 2);
 
-  for (v_t i = 0; i < pow(L, D); i++) {
+    for (unsigned j = 0; j < pow(L, 2); j++) {
+      unsigned vc = sb + j;
+
+      unsigned xyv1 = j;
+      unsigned xyv2;
+      if (i == 0) {
+        xyv2 = L * (xyv1 / L) + ((xyv1 % L) + 1) % L;
+      } else {
+        xyv2 = L * (((xyv1 / L) + 1) % L) + (xyv1 % L);
+      }
+
+      si.G.vertices[vc].prop.back = &(sxy.s[xyv1]);
+      si.G.vertices[vc].prop.front = &(sxy.s[xyv2]);
+    }
+  }
+
+  // start in the ground state of the Ising model
+  for (unsigned i = 0; i < pow(L, D); i++) {
     si.s[pow(L, D) + i].x = true;
   }
 
-  E = - J * 4 * pow(L, D) - 2 * pow(L, D);
+  cos_xy = Jx * 2 * pow(L, 2) + Jxn * 2 * pow(L, 2);
+  sin_xy.fill(0);
+
+  i_measurement mi(si, a);
+  x_measurement mxy(sxy, Jx, Jxn, a, &(si.s0));
 
-  i_measurement mi(si);
-  x_measurement mxy(sxy);
+  double total_cos = 0;
+  double total_sin2 = 0;
 
-  double total_E;
-  double total_E2;
+  for (unsigned i = 0; i < N; i++) {
+    si.run_wolff(1, ising::wolff::gen_ising<graph_i>, mi, rng);
+    sxy.run_wolff(1, xy::wolff::generate_rotation_uniform<2, graph_x>, mxy, rng);
+    total_cos += cos_xy;
+    total_sin2 += sin_xy * sin_xy;
 
-  for (N_t i = 0; i < N; i++) {
-    si.run_wolff(1, ising::wolff::gen_ising, mi, rng);
-    sxy.run_wolff(1, xy::wolff::generate_rotation_uniform<2>, mxy, rng);
-    total_E += E;
-    total_E2 += pow(E, 2);
+//    std::cout << cos_xy / sxy.ne << " " << sin_xy / pow(sxy.ne, 2) << " " << E << "\n";
   }
 
-  std::cout << mi.avgC() / si.nv << " " << mi.avgA2() / pow(si.nv, 2) << "\n";
-  std::cout << mxy.avgM() / pow(sxy.nv, 1) << " " << 2 * mxy.avgCM() / pow(sxy.nv, 0) << "\n";
+  std::ofstream outfile;
+  outfile.open("out.dat", std::ios::app);
+  outfile << L << " " << Ji << " " << Jx << " " << a << " " << mi.totalA2 / N / pow(si.nv, 2) << " " << (total_cos - total_sin2 / T) / N / sxy.ne << " " <<(total_cos + total_sin2 / T) / N / sxy.ne << "\n";
+  outfile.close();
 }