From f7a21799194f6626994195302ac95449158bdcd9 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 25 May 2017 13:26:38 -0400 Subject: started wolff code in new repository --- .gitignore | 1 + CMakeLists.txt | 21 ++++ lib/queue.c | 32 ++++++ lib/wolff.h | 51 +++++++++ lib/wolff_tools.c | 179 ++++++++++++++++++++++++++++++ src/pt_wolff.c | 317 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/wolff.c | 140 ++++++++++++++++++++++++ src/wolff.h | 52 +++++++++ 8 files changed, 793 insertions(+) create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 lib/queue.c create mode 100644 lib/wolff.h create mode 100644 lib/wolff_tools.c create mode 100644 src/pt_wolff.c create mode 100644 src/wolff.c create mode 100644 src/wolff.h diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a007fea --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +build/* diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..c325cb4 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,21 @@ + +cmake_minimum_required(VERSION 3.7) +project(wolff) + +include_directories(src ~/.local/include) +link_directories(~/.local/lib) + +file(GLOB SOURCES lib/*.c) +add_executable(wolff src/wolff.c ${SOURCES}) + +find_package(OpenMP) +if (OPENMP_FOUND) + set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") +endif() + +target_link_libraries(wolff gsl m jst cblas fftw3 tcmalloc profiler) +#target_link_libraries(pt_wolff gsl m jst cblas tcmalloc profiler) + +install(TARGETS wolff DESTINATION bin) + diff --git a/lib/queue.c b/lib/queue.c new file mode 100644 index 0000000..9a01741 --- /dev/null +++ b/lib/queue.c @@ -0,0 +1,32 @@ + +#include "wolff.h" + +void stack_push(ll_t **q, uint32_t x) { + ll_t *nq = malloc(sizeof(ll_t)); + nq->x = x; + nq->next = *q; + + *q = nq; +} + +uint32_t stack_pop(ll_t **q) { + ll_t *old_q = *q; + + *q = old_q->next; + uint32_t x = old_q->x; + + free(old_q); + + return x; +} + +bool stack_contains(const ll_t *q, uint32_t x) { + if (q == NULL) { + return false; + } else if (q->x == x) { + return true; + } else { + return stack_contains(q->next, x); + } +} + diff --git a/lib/wolff.h b/lib/wolff.h new file mode 100644 index 0000000..6beb03e --- /dev/null +++ b/lib/wolff.h @@ -0,0 +1,51 @@ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +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); + diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c new file mode 100644 index 0000000..72cf9f7 --- /dev/null +++ b/lib/wolff_tools.c @@ -0,0 +1,179 @@ + +#include "wolff.h" + +graph_t *graph_add_ext(const graph_t *g) { + graph_t *h = (graph_t *)calloc(1, sizeof(graph_t)); + + h->nv = g->nv + 1; + h->ne = g->ne + g->nv; + + h->ev = (uint32_t *)malloc(2 * h->ne * sizeof(uint32_t)); + h->vei = (uint32_t *)malloc((h->nv + 1) * sizeof(uint32_t)); + h->ve = (uint32_t *) malloc(2 * h->ne * sizeof(uint32_t)); + h->vx = (double *)malloc(2 * h->nv * sizeof(double)); + h->bq = (bool *)malloc(h->nv * sizeof(bool)); + + memcpy(h->ev, g->ev, 2 * g->ne * sizeof(uint32_t)); + memcpy(h->vx, g->vx, 2 * g->nv * sizeof(double)); + memcpy(h->bq, g->bq, g->nv * sizeof(bool)); + h->vx[2 * g->nv] = -1; + h->vx[2 * g->nv + 1] = -0.5; + h->bq[g->nv] = false; + + for (uint32_t i = 0; i < g->nv; i++) { + h->ev[2 * g->ne + 2 * i] = i; + h->ev[2 * g->ne + 2 * i + 1] = g->nv; + } + + for (uint32_t i = 0; i < g->nv; i++) { + h->vei[i] = g->vei[i] + i; + + for (uint32_t j = 0; j < g->vei[i + 1] - g->vei[i]; j++) { + h->ve[h->vei[i] + j] = g->ve[g->vei[i] + j]; + } + + h->ve[h->vei[i] + g->vei[i + 1] - g->vei[i]] = g->ne + i; + } + + h->vei[g->nv] = g->vei[g->nv] + g->nv; + h->vei[g->nv + 1] = h->vei[g->nv] + g->nv; + + for (uint32_t i = 0; i < g->nv; i++) { + h->ve[h->vei[g->nv] + i] = g->ne + i; + } + + return h; +} + +double get_hamiltonian(graph_t *g, double *coupling, bool *x) { + double hamiltonian = 0; + + for (uint32_t i = 0; i < g->ne; i++) { + uint32_t v1, v2; + + v1 = g->ev[2 * i]; + v2 = g->ev[2 * i + 1]; + + if (x[v1] == x[v2]) { + hamiltonian -= coupling[i]; + } else { + hamiltonian += coupling[i]; + } + } + + return hamiltonian; +} + +cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, gsl_rng *r) { + uint32_t v0; + bool x0; + cluster_t *c; + + v0 = gsl_rng_uniform_int(r, g->nv); // pick a random vertex + x0 = x[v0]; // record its orientation + + ll_t *stack = NULL; // create a new stack + stack_push(&stack, v0); // push the initial vertex to the stack + + // initiate the data structure for returning flip information + c = (cluster_t *)calloc(1, sizeof(cluster_t)); + c->nv = 0; + c->dH = 0; + + while (stack != NULL) { + uint32_t v; + uint16_t nn; + + v = stack_pop(&stack); + nn = g->vei[v + 1] - g->vei[v]; + + if (x[v] == x0) { // if the vertex hasn't already been flipped + x[v] = !x[v]; // flip the vertex + + for (uint16_t i = 0; i < nn; i++) { + uint32_t e, v1, v2, vn; + + e = g->ve[g->vei[v] + i]; // select the ith bond connected to site + v1 = g->ev[2 * e]; + v2 = g->ev[2 * e + 1]; + + vn = v == v1 ? v2 : v1; // distinguish neighboring site from site itself + + if (x[vn] == x0) { // if the neighboring site matches the flipping cluster... + if (v1 == g->nv - 1 || v2 == g->nv - 1) { + c->dH += H; + } else { + c->dH += 1; + } + + if (gsl_rng_uniform(r) < ps[e]) { // and with probability ps[e]... + stack_push(&stack, vn); // push the neighboring vertex to the stack + } + } else { + if (v1 == g->nv - 1 || v2 == g->nv - 1) { + c->dH -= H; + } else { + c->dH -= 1; + } + } + } + + if (v != g->nv - 1) { + c->nv++; + } + } + } + + if (x0) { + c->nv = -c->nv; + } + + return c; +} + +double hh(double th) { + return (th - pow(th, 3) / 1.16951) * (1 - 0.222389 * pow(th, 2) - 0.043547 * pow(th, 4) - 0.014809 * pow(th, 6) - 0.007168 * pow(th, 8)); +} + +double *get_bond_probs(double T, double H, ising_state_t *s) { + double p = 1 - exp(-2 / T); + double q = 1 - exp(-2 * H / T); + + double *ps = (double *)malloc(s->g->ne * sizeof(double)); + + for (uint32_t i = 0; i < s->g->ne; i++) { + uint32_t v1, v2; + v1 = s->g->ev[2 * i]; + v2 = s->g->ev[2 * i + 1]; + if (v1 == s->g->nv - 1 || v2 == s->g->nv - 1) { + ps[i] = q; + } else { + ps[i] = p; + } + } + + return ps; +} + +int32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps) { + if (r == NULL) { + r = gsl_rng_alloc(gsl_rng_mt19937); + gsl_rng_set(r, jst_rand_seed()); + } + + if (ps == NULL) { + ps = get_bond_probs(T, H, s); + } + + cluster_t *c = flip_cluster(s->g, ps, H, s->spins, r); + + s->M += 2 * c->nv; + s->H += c->dH; + + int32_t n_flipped = c->nv; + + free(c); + + return n_flipped; +} + 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +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); + -- cgit v1.2.3-70-g09d2