summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/wolff.h3
-rw-r--r--lib/wolff_tools.c19
-rw-r--r--src/wolff.c84
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);