summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/pt_wolff.c317
-rw-r--r--src/wolff.c140
-rw-r--r--src/wolff.h52
3 files changed, 509 insertions, 0 deletions
diff --git a/src/pt_wolff.c b/src/pt_wolff.c
new file mode 100644
index 0000000..128b601
--- /dev/null
+++ b/src/pt_wolff.c
@@ -0,0 +1,317 @@
+
+#include "wolff.h"
+
+#define TC 2. / log(1. + sqrt(2.))
+
+int main(int argc, char *argv[]) {
+ int opt;
+ bool use_scho, use_rand_ini;
+ lattice_t lat;
+ uint16_t L;
+ uint32_t n_temps;
+ uint64_t N;
+ double Tin, Hin, eps, dT, dH;
+
+ L = 128;
+ N = 1000;
+ lat = SQUARE_LATTICE;
+ Tin = 2.3;
+ Hin = 0;
+ eps = 1e30;
+ use_scho = false;
+ use_rand_ini = false;
+ dT = 0;
+ dH = 1;
+ n_temps = 1;
+
+ while ((opt = getopt(argc, argv, "N:L:q:T:H:e:St:h:n:r")) != -1) {
+ switch (opt) {
+ case 'S':
+ use_scho = true;
+ break;
+ case 'r':
+ use_rand_ini = true;
+ break;
+ case 'N':
+ N = (uint64_t)atof(optarg);
+ break;
+ case 'L':
+ L = atoi(optarg);
+ break;
+ case 'T':
+ Tin = atof(optarg);
+ break;
+ case 'H':
+ Hin = atof(optarg);
+ break;
+ case 't':
+ dT = atof(optarg);
+ break;
+ case 'h':
+ dH = atof(optarg);
+ break;
+ case 'n':
+ n_temps = atoi(optarg);
+ break;
+ case 'e':
+ eps= atof(optarg);
+ break;
+ default:
+ exit(EXIT_FAILURE);
+ }
+ }
+
+ gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+ gsl_rng_set(r, jst_rand_seed());
+
+ graph_t *g = graph_create(lat, TORUS_BOUND, L, false);
+
+ double *Ts = (double *)malloc(n_temps * sizeof(double));
+ double *Hs = (double *)malloc(n_temps * sizeof(double));
+ double *Tins = (double *)malloc(n_temps * sizeof(double));
+ double *Hins = (double *)malloc(n_temps * sizeof(double));
+
+ for (uint32_t i = 0; i < n_temps; i++) {
+ double T, H;
+
+ // if Schofield coordinates are used, T and H are treated as R and theta.
+ if (use_scho) {
+ double R = Tin; double th = Hin;
+ double t = R * (1 - pow(th, 2));
+ T = TC / (1 - t);
+ double h0 = 0.940647;
+ double h = h0 * pow(R, 15. / 8.) * hh(th);
+ H = h * T;
+ } else {
+ T = Tin;
+ H = Hin;
+ }
+
+ Tins[i] = Tin;
+ Hins[i] = Hin;
+
+ Ts[i] = T;
+ Hs[i] = H;
+
+ Tin += dT;
+ if (use_scho) {
+ Hin = 1.08144 - (1.08144 - Hin) / dH;
+ } else {
+ Hin = Hin / dH;
+ }
+ }
+
+ bool **xs = (bool **)malloc(n_temps * sizeof(bool *));
+ for (uint32_t i = 0; i < n_temps; i++) {
+ xs[i] = (bool *)calloc(g->nv, sizeof(bool));
+ }
+
+ if (use_rand_ini) {
+ for (uint32_t j = 0; j < n_temps; j++) {
+ for (uint32_t i = 0; i < g->nv; i++) {
+ xs[j][i] = gsl_rng_uniform_int(r, 2);
+ }
+ }
+ }
+
+ int32_t *Ms = (int32_t *)calloc(n_temps, sizeof(int32_t));
+ for (uint32_t j = 0; j < n_temps; j++) {
+ if (use_rand_ini) {
+ for (uint32_t i = 0; i < g->nv; i++) {
+ if (xs[j][i]) {
+ Ms[j]++;
+ } else {
+ Ms[j]--;
+ }
+ }
+ } else {
+ Ms[j] = -g->nv;
+ }
+ }
+
+ double *Es = (double *)malloc(n_temps * sizeof(double));
+ for (uint32_t i = 0; i < n_temps; i++) {
+ Es[i] = get_hamiltonian(g, Hs[i], xs[i]);
+ }
+
+ double *ps = (double *)malloc(n_temps * sizeof(double));
+ for (uint32_t i = 0; i < n_temps; i++) {
+ ps[i] = 1 - exp(-2 / Ts[i]);
+ }
+
+ double *E1s = (double *)calloc(n_temps, sizeof(double));
+ double *E2s = (double *)calloc(n_temps, sizeof(double));
+ double *E4s = (double *)calloc(n_temps, sizeof(double));
+
+ double *M1s = (double *)calloc(n_temps, sizeof(double));
+ double *M2s = (double *)calloc(n_temps, sizeof(double));
+ double *M4s = (double *)calloc(n_temps, sizeof(double));
+
+ double *tE1s = (double *)calloc(n_temps, sizeof(double));
+ double *tE2s = (double *)calloc(n_temps, sizeof(double));
+ double *tE4s = (double *)calloc(n_temps, sizeof(double));
+
+ double *tM1s = (double *)calloc(n_temps, sizeof(double));
+ double *tM2s = (double *)calloc(n_temps, sizeof(double));
+ double *tM4s = (double *)calloc(n_temps, sizeof(double));
+
+ uint64_t n_runs = 0;
+
+ double diff = 1e30;
+
+ printf("\n");
+ while (diff > eps) {
+ printf("\033[F\033[JWOLFF: sweep %llu, %g\n", n_runs, diff);
+
+ for (uint32_t j = 0; j < n_temps; j++) {
+ tE1s[j] = 0;
+ tE2s[j] = 0;
+ tE4s[j] = 0;
+ tM1s[j] = 0;
+ tM2s[j] = 0;
+ tM4s[j] = 0;
+ }
+
+ for (uint64_t i = 0; i < N; i++) {
+#pragma omp parallel for
+ for (uint32_t j = 0; j < n_temps; j++) {
+ cluster_t *c = get_cluster(g, xs[j], ps[j], r);
+ double dHex = 0;
+ double s;
+
+ if (xs[j][c->vi->x]) {
+ s = 1;
+ } else {
+ s = -1;
+ }
+
+ dHex = s * Hs[j] * c->nv;
+
+ if (gsl_rng_uniform(r) < exp(-2 * dHex / Ts[j])) {
+ while (c->vi != NULL) {
+ uint32_t v = queue_del(&(c->vi));
+ xs[j][v] = !xs[j][v];
+ }
+ Es[j] += 2 * c->nb + dHex;
+ Ms[j] -= 2 * s * c->nv;
+ } else {
+ while (c->vi != NULL) {
+ queue_del(&(c->vi));
+ }
+ }
+
+ tE1s[j] += Es[j];
+ tM1s[j] += abs(Ms[j]);
+ tE2s[j] += pow(Es[j], 2);
+ tM2s[j] += pow(Ms[j], 2);
+ tE4s[j] += pow(Es[j], 4);
+ tM4s[j] += pow(Ms[j], 4);
+
+ free(c);
+ }
+
+ for (uint32_t j = 0; j < n_temps - 1; j++) {
+ if (gsl_rng_uniform(r) < 0.5) {
+ if (gsl_rng_uniform(r) < exp((Es[j + 1] - Es[j]) * (1 / Ts[j + 1] - 1 / Ts[j]))) {
+ bool *tmp_x = xs[j];
+ double tmp_E = Es[j];
+ int32_t tmp_M = Ms[j];
+
+ xs[j] = xs[j + 1];
+ xs[j + 1] = tmp_x;
+
+ Es[j] = Es[j + 1];
+ Es[j + 1] = tmp_E;
+
+ Ms[j] = Ms[j + 1];
+ Ms[j + 1] = tmp_M;
+ }
+ }
+ }
+ }
+
+ for (uint32_t j = 0; j < n_temps; j++) {
+ tE1s[j] /= N;
+ tM1s[j] /= N;
+ tE2s[j] /= N;
+ tM2s[j] /= N;
+ tE4s[j] /= N;
+ tM4s[j] /= N;
+
+ if (n_runs > 0) {
+ E1s[j] = E1s[j] * ((n_runs - 1.) / n_runs) + tE1s[j] * 1. / n_runs;
+ M1s[j] = M1s[j] * ((n_runs - 1.) / n_runs) + tM1s[j] * 1. / n_runs;
+ E2s[j] = E2s[j] * ((n_runs - 1.) / n_runs) + tE2s[j] * 1. / n_runs;
+ M2s[j] = M2s[j] * ((n_runs - 1.) / n_runs) + tM2s[j] * 1. / n_runs;
+ E4s[j] = E4s[j] * ((n_runs - 1.) / n_runs) + tE4s[j] * 1. / n_runs;
+ M4s[j] = M4s[j] * ((n_runs - 1.) / n_runs) + tM4s[j] * 1. / n_runs;
+ }
+ }
+
+ if (n_runs > 0) {
+ diff = 0;
+ for (uint32_t j = 0; j < n_temps; j++) {
+ double dd = sqrt(M2s[j] - pow(M1s[j], 2)) / sqrt(N * n_runs);
+ double t_diff = dd / M1s[j];
+ if (t_diff > diff) diff = t_diff;
+ }
+ } else {
+ diff = 1e30;
+ }
+
+ n_runs++;
+ }
+
+ FILE *outfile = fopen("out.dat", "a");
+
+ for (uint32_t j = 0; j < n_temps; j++) {
+ double C = (E2s[j] - pow(E1s[j], 2)) / pow(Ts[j], 2);
+ double X = (M2s[j] - pow(M1s[j], 2)) / Ts[j];
+
+ double dE1 = sqrt(E2s[j] - pow(E1s[j], 2)) / sqrt(N * n_runs);
+ double dM1 = sqrt(M2s[j] - pow(M1s[j], 2)) / sqrt(N * n_runs);
+ double dE2 = sqrt(E4s[j] - pow(E2s[j], 2)) / sqrt(N * n_runs);
+ double dM2 = sqrt(M4s[j] - pow(M2s[j], 2)) / sqrt(N * n_runs);
+
+ double dC = sqrt(pow(dE2, 2) + pow(2 * E1s[j] * dE1, 2)) / pow(Ts[j], 2);
+ double dX = sqrt(pow(dM2, 2) + pow(2 * M1s[j] * dM1, 2)) / Ts[j];
+
+ fprintf(outfile, "%u %u %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", L, use_scho, Tins[j], Hins[j], E1s[j] / g->nv, dE1 / g->nv, M1s[j] / g->nv, dM1 / g->nv, C / g->nv, dC / g->nv, X / g->nv, dX / g->nv);
+ }
+
+ fclose(outfile);
+
+ for (uint32_t i = 0; i < n_temps; i++) {
+ free(xs[i]);
+ }
+
+ free(xs);
+
+ gsl_rng_free(r);
+ graph_free(g);
+
+ free(Ts);
+ free(Hs);
+ free(Es);
+ free(Ms);
+
+ free(E1s);
+ free(E2s);
+ free(E4s);
+
+ free(M1s);
+ free(M2s);
+ free(M4s);
+
+ free(tE1s);
+ free(tE2s);
+ free(tE4s);
+
+ free(tM1s);
+ free(tM2s);
+ free(tM4s);
+
+
+ return 0;
+}
+
diff --git a/src/wolff.c b/src/wolff.c
new file mode 100644
index 0000000..80f8278
--- /dev/null
+++ b/src/wolff.c
@@ -0,0 +1,140 @@
+
+#include "wolff.h"
+
+int main(int argc, char *argv[]) {
+ int opt;
+ bool output_state;
+ lattice_t lat;
+ uint16_t L;
+ uint32_t min_runs;
+ uint64_t N;
+ double T, H, eps;
+
+ L = 128;
+ N = 1000;
+ lat = SQUARE_LATTICE;
+ T = 2.3;
+ H = 0;
+ eps = 1e30;
+ output_state = false;
+ min_runs = 10;
+
+ while ((opt = getopt(argc, argv, "N:L:T:H:m:e:o")) != -1) {
+ switch (opt) {
+ case 'N':
+ N = (uint64_t)atof(optarg);
+ break;
+ case 'L':
+ L = atoi(optarg);
+ break;
+ case 'T':
+ T = atof(optarg);
+ break;
+ case 'H':
+ H = atof(optarg);
+ break;
+ case 'm':
+ min_runs = atoi(optarg);
+ break;
+ case 'e':
+ eps= atof(optarg);
+ break;
+ case 'o':
+ output_state = true;
+ break;
+ default:
+ exit(EXIT_FAILURE);
+ }
+ }
+
+ gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+ gsl_rng_set(r, jst_rand_seed());
+
+ graph_t *h = graph_create(lat, TORUS_BOUND, L, false);
+
+ ising_state_t *s = (ising_state_t *)calloc(1, sizeof(ising_state_t));
+
+ s->g = graph_add_ext(h);
+ s->spins = (bool *)calloc(h->nv + 1, sizeof(bool));
+ s->M = -h->nv;
+ s->H = -(1.0 * h->ne) - H * h->nv;
+
+ double *bond_probs = get_bond_probs(T, H, s);
+
+ double diff = 1e31;
+ uint64_t n_runs = 0;
+
+ double E1, E2, dE1, M1, M2, dM1, C, dC, X, dX, Mmu2, Mmu4, Emu2, Emu4;
+ double clust_per_sweep = 0;
+
+ E1 = 0; E2 = 0; M1 = 0; M2 = 0; C = 0; dC = 0; X = 0; dX = 0;
+ dE1 = 0; dM1 = 0; Mmu2 = 0; Mmu4 = 0; Emu2 = 0; Emu4 = 0;
+
+ printf("\n");
+ while (diff > eps || diff == 0. || diff != diff || n_runs < min_runs) {
+ printf("\033[F\033[JWOLFF: sweep %llu, dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n", n_runs, fabs(dE1 / E1), dM1 / M1, dC / C, dX / X, clust_per_sweep);
+
+ uint32_t n_flips = 0;
+ uint32_t n_clust = 0;
+
+ while (n_flips / h->nv < N) {
+ n_flips += abs(wolff_step(T, H, s, r, bond_probs));
+ n_clust++;
+ }
+
+ E1 = E1 * (n_runs / (n_runs + 1.)) + s->H * 1. / (n_runs + 1.);
+ M1 = M1 * (n_runs / (n_runs + 1.)) + abs(s->M) * 1. / (n_runs + 1.);
+ E2 = E2 * (n_runs / (n_runs + 1.)) + pow(s->H, 2) * 1. / (n_runs + 1.);
+ M2 = M2 * (n_runs / (n_runs + 1.)) + pow(s->M, 2) * 1. / (n_runs + 1.);
+
+ Mmu2 = Mmu2 * (n_runs / (n_runs + 1.)) + pow(abs(s->M) - M1, 2) * 1. / (n_runs + 1.);
+ Mmu4 = Mmu4 * (n_runs / (n_runs + 1.)) + pow(abs(s->M) - M1, 4) * 1. / (n_runs + 1.);
+ Emu2 = Emu2 * (n_runs / (n_runs + 1.)) + pow(s->H - E1, 2) * 1. / (n_runs + 1.);
+ Emu4 = Emu4 * (n_runs / (n_runs + 1.)) + pow(s->H - E1, 4) * 1. / (n_runs + 1.);
+
+ if (n_runs > 1){
+ double Msigma2 = n_runs / (n_runs - 1) * (M2 - pow(M1, 2));
+ X = Msigma2 / T;
+ dX = sqrt((Mmu4 - (n_runs - 3.) / (n_runs - 1.) * pow(Mmu2, 2)) / n_runs) / T;
+
+ double Esigma2 = n_runs / (n_runs - 1) * (E2 - pow(E1, 2));
+ C = Esigma2 / T;
+ dC = sqrt((Emu4 - (n_runs - 3.) / (n_runs - 1.) * pow(Emu2, 2)) / n_runs) / T;
+
+ dE1 = sqrt(Esigma2 / n_runs);
+ dM1 = sqrt(Msigma2 / n_runs);
+
+ diff = fabs(dX / X);
+ }
+
+ clust_per_sweep = clust_per_sweep * (n_runs / (n_runs + 1.)) + (n_clust * 1. / N) * 1. / (n_runs + 1.);
+
+ n_runs++;
+ }
+ printf("\033[F\033[JWOLFF: sweep %llu, dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n", n_runs, fabs(dE1 / E1), dM1 / M1, dC / C, dX / X, clust_per_sweep);
+
+ FILE *outfile = fopen("out.dat", "a");
+ fprintf(outfile, "%u %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", L, T, H, E1 / h->nv, dE1 / h->nv, M1 / h->nv, dM1 / h->nv, C / h->nv, dC / h->nv, X / h->nv, dX / h->nv);
+ fclose(outfile);
+
+ free(bond_probs);
+
+ if (output_state) {
+ FILE *state_file = fopen("state.dat", "a");
+ for (uint32_t i = 0; i < h->nv; i++) {
+ fprintf(state_file, "%d ", s->spins[i]);
+ }
+ fprintf(state_file, "\n");
+ fclose(state_file);
+ }
+
+ gsl_rng_free(r);
+ graph_free(s->g);
+ free(s->spins);
+ free(s);
+ free(bond_probs);
+ graph_free(h);
+
+ return 0;
+}
+
diff --git a/src/wolff.h b/src/wolff.h
new file mode 100644
index 0000000..9619a4d
--- /dev/null
+++ b/src/wolff.h
@@ -0,0 +1,52 @@
+
+#include <string.h>
+#include <math.h>
+#include <getopt.h>
+#include <float.h>
+#include <sys/types.h>
+#include <inttypes.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_rng.h>
+#include <stdbool.h>
+#include <assert.h>
+#include <omp.h>
+#include <fftw3.h>
+
+#include <jst/graph.h>
+#include <jst/rand.h>
+
+typedef struct {
+ graph_t *g;
+ bool *spins;
+ int32_t M;
+ double H;
+} ising_state_t;
+
+typedef struct ll_tag {
+ uint32_t x;
+ struct ll_tag *next;
+} ll_t;
+
+typedef struct {
+ int32_t nv;
+ double dH;
+} cluster_t;
+
+double get_hamiltonian(graph_t *g, double *coupling, bool *x);
+
+void stack_push(ll_t **q, uint32_t x);
+
+uint32_t stack_pop(ll_t **q);
+
+bool stack_contains(const ll_t *q, uint32_t x);
+
+cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, gsl_rng *r);
+
+graph_t *graph_add_ext(const graph_t *g);
+
+double hh(double th);
+
+double *get_bond_probs(double T, double H, ising_state_t *s);
+
+int32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps);
+