From 12a3623728dcbe2ef1e8082310c86cec4e4578d8 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sun, 4 Feb 2018 14:18:44 -0500
Subject: added support for arbitrary spin interactions

---
 lib/wolff.h       |  3 +-
 lib/wolff_tools.c | 19 ++++++-------
 src/wolff.c       | 84 +++++++++++++++++++++++++++++++++++--------------------
 3 files changed, 65 insertions(+), 41 deletions(-)

diff --git a/lib/wolff.h b/lib/wolff.h
index dae907a..44e5926 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -24,8 +24,9 @@ typedef struct {
   graph_t *g;
   q_t *spins;
   double T;
+  double *J;
   double *H;
-  double T_prob;
+  double *J_probs;
   double *H_probs;
   double E;
   v_t *M;
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index ddce40f..20bc529 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -1,10 +1,6 @@
 
 #include "wolff.h"
 
-q_t delta(q_t s0, q_t s1) {
-  return s0 == s1 ? 1 : 0;
-}
-
 v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) {
   q_t s0 = s->spins[v0];
   v_t nv = 0;
@@ -36,12 +32,13 @@ v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) {
 
         bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1);
 
+        q_t M_ind_0;
+        q_t M_ind_1;
+
         if (is_ext) {
-            q_t M_ind_0;
-            q_t M_ind_1;
           if (vn == s->g->nv - 1) {
-            M_ind_0 = (s_old + s->q - s->spins[s->g->nv - 1]) % s->q;
-            M_ind_1 = (s_new + s->q - s->spins[s->g->nv - 1]) % s->q;
+            M_ind_0 = (s_old + s->q - sn) % s->q;
+            M_ind_1 = (s_new + s->q - sn) % s->q;
           } else {
             M_ind_0 = (sn + s->q - s_old) % s->q;
             M_ind_1 = (sn + s->q - s_new) % s->q;
@@ -51,8 +48,10 @@ v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) {
           s->M[M_ind_1]++;
           s->E += - s->H[M_ind_1] + s->H[M_ind_0];
         } else {
-          prob = s->T_prob * delta(s0, sn);
-          s->E += - ((double)delta(s->spins[v], sn)) + ((double)delta(s0, sn));
+          M_ind_0 = (s_old + s->q - sn) % s->q;
+          M_ind_1 = (s_new + s->q - sn) % s->q;
+          prob = s->J_probs[M_ind_1 * s->q + M_ind_0];
+          s->E += - s->J[M_ind_1] + s->J[M_ind_0];
         }
 
         if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
diff --git a/src/wolff.c b/src/wolff.c
index 67c05a5..fee796f 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -10,14 +10,18 @@ int main(int argc, char *argv[]) {
   q_t q = 2;
   D_t D = 2;
   double T = 2.26918531421;
+  double *J = (double *)calloc(MAX_Q, sizeof(double));
+  J[0] = 1.0;
   double *H = (double *)calloc(MAX_Q, sizeof(double));
   double eps = 0;
   bool pretend_ising = false;
+  bool planar_potts = false;
 
   int opt;
+  q_t J_ind = 0;
   q_t H_ind = 0;
 
-  while ((opt = getopt(argc, argv, "N:n:D:L:q:T:H:m:e:I")) != -1) {
+  while ((opt = getopt(argc, argv, "N:n:D:L:q:T:J:H:m:e:Ip")) != -1) {
     switch (opt) {
     case 'N':
       N = (count_t)atof(optarg);
@@ -37,6 +41,10 @@ int main(int argc, char *argv[]) {
     case 'T':
       T = atof(optarg);
       break;
+    case 'J':
+      J[J_ind] = atof(optarg);
+      J_ind++;
+      break;
     case 'H':
       H[H_ind] = atof(optarg);
       H_ind++;
@@ -50,6 +58,9 @@ int main(int argc, char *argv[]) {
     case 'I':
       pretend_ising = true;
       break;
+    case 'p':
+      planar_potts = true;
+      break;
     default:
       exit(EXIT_FAILURE);
     }
@@ -60,9 +71,14 @@ int main(int argc, char *argv[]) {
 
   if (pretend_ising) {
     q = 2;
-    T /= 2;
-    H[0] /= 2;
     H[1] = -H[0];
+    J[1] = -J[0];
+  }
+
+  if (planar_potts) {
+    for (q_t i = 0; i < q; i++) {
+      J[i] = cos(2 * M_PI * i / ((double)q));
+    }
   }
 
   ising_state_t *s = (ising_state_t *)calloc(1, sizeof(ising_state_t));
@@ -76,8 +92,14 @@ int main(int argc, char *argv[]) {
 
   s->T = T;
   s->H = H;
+  s->J = J;
 
-  s->T_prob = 1.0 - exp(- 1.0 / T);
+  s->J_probs = (double *)calloc(pow(q, 2), sizeof(double));
+  for (q_t i = 0; i < q; i++) {
+    for (q_t j = 0; j < q; j++) {
+      s->J_probs[q * i + j] = 1.0 - exp((s->J[i] - s->J[j]) / T);
+    }
+  }
   s->H_probs = (double *)calloc(pow(q, 2), sizeof(double));
   for (q_t i = 0; i < q; i++) {
     for (q_t j = 0; j < q; j++) {
@@ -87,7 +109,7 @@ int main(int argc, char *argv[]) {
 
   s->M = (v_t *)calloc(q, sizeof(v_t));
   s->M[0] = h->nv;
-  s->E = - ((double)h->nv) * s->H[0];
+  s->E = - ((double)h->ne) * s->J[0] - ((double)h->nv) * s->H[0];
 
   double diff = 1e31;
   count_t n_runs = 0;
@@ -137,35 +159,35 @@ int main(int argc, char *argv[]) {
 
   FILE *outfile = fopen("out.m", "a");
 
-  if (pretend_ising) {
-    fprintf(outfile,
-            "%" PRIL " %.15f %" PRIcount "  %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n",
-            L, 2 * T, n_runs,
-            (2 * E->x - h->ne) / h->nv, 2 * E->dx / h->nv, (M[0]->x - M[1]->x) / h->nv, M[0]->dx / h->nv, 4 * E->c / h->nv, 4 * E->dc / h->nv, 4 * M[0]->c / h->nv, 4 * M[0]->dc / h->nv); 
-  } else {
-    fprintf(outfile, "{\"D\"->%" PRID ",\"L\"->%" PRIL ",\"q\"->%" PRIq ",\"T\"->%.15f,\"H\"->{", D, L, q, T);
-    for (q_t i = 0; i < q; i++) {
-      fprintf(outfile, "%.15f", H[i]);
-      if (i != q-1) {
-        fprintf(outfile, ",");
-      }
+  fprintf(outfile, "{\"D\"->%" PRID ",\"L\"->%" PRIL ",\"q\"->%" PRIq ",\"T\"->%.15f,\"J\"->{", D, L, q, T);
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "%.15f", J[i]);
+    if (i != q-1) {
+      fprintf(outfile, ",");
     }
-    fprintf(outfile, "},\"E\"->{%.15f,%.15f},\"C\"->{%.15f,%.15f},\"M\"->{", E->x / h->nv, E->dx / h->nv, E->c / h->nv, E->dc / h->nv);
-    for (q_t i = 0; i < q; i++) {
-      fprintf(outfile, "{%.15f,%.15f}", M[i]->x / h->nv, M[i]->dx / h->nv);
-      if (i != q-1) {
-        fprintf(outfile, ",");
-      }
+  }
+  fprintf(outfile, "},\"H\"->{");
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "%.15f", H[i]);
+    if (i != q-1) {
+      fprintf(outfile, ",");
     }
-    fprintf(outfile, "},\"X\"->{");
-    for (q_t i = 0; i < q; i++) {
-      fprintf(outfile, "{%.15f,%.15f}", M[i]->c / h->nv, M[i]->dc / h->nv);
-      if (i != q-1) {
-        fprintf(outfile, ",");
-      }
+  }
+  fprintf(outfile, "},\"E\"->{%.15f,%.15f},\"C\"->{%.15f,%.15f},\"M\"->{", E->x / h->nv, E->dx / h->nv, E->c / h->nv, E->dc / h->nv);
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "{%.15f,%.15f}", M[i]->x / h->nv, M[i]->dx / h->nv);
+    if (i != q-1) {
+      fprintf(outfile, ",");
+    }
+  }
+  fprintf(outfile, "},\"X\"->{");
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "{%.15f,%.15f}", M[i]->c / h->nv, M[i]->dc / h->nv);
+    if (i != q-1) {
+      fprintf(outfile, ",");
     }
-    fprintf(outfile, "},\"n\"->{%.15f,%.15f}}\n", clust->c / h->nv, clust->dc / h->nv);
   }
+  fprintf(outfile, "},\"n\"->{%.15f,%.15f}}\n", clust->c / h->nv, clust->dc / h->nv);
 
   fclose(outfile);
 
@@ -176,11 +198,13 @@ int main(int argc, char *argv[]) {
   }
   free(M);
   free(s->H_probs);
+  free(s->J_probs);
   free(s->M);
   free(s->spins);
   graph_free(s->g);
   free(s);
   free(H);
+  free(J);
   graph_free(h);
   gsl_rng_free(r);
 
-- 
cgit v1.2.3-70-g09d2