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

---
 src/wolff.c | 84 +++++++++++++++++++++++++++++++++++++++----------------------
 1 file changed, 54 insertions(+), 30 deletions(-)

(limited to 'src/wolff.c')

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