From 9db6ee734df8477a2529f56e4a6f4b1784bf941b Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 22 Jun 2017 16:05:05 -0400 Subject: many changes, simplification of some functions, removal of unneeded ones --- CMakeLists.txt | 4 +- lib/queue.c | 21 ++-- lib/queue.h | 18 ++++ lib/wolff.h | 44 +++----- lib/wolff_tools.c | 61 ++++------- src/pt_wolff.c | 317 ------------------------------------------------------ src/wolff.c | 140 ++++++++++++++---------- src/wolff.h | 52 --------- 8 files changed, 149 insertions(+), 508 deletions(-) create mode 100644 lib/queue.h delete mode 100644 src/pt_wolff.c delete mode 100644 src/wolff.h diff --git a/CMakeLists.txt b/CMakeLists.txt index c325cb4..66d9e2c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,8 +1,8 @@ -cmake_minimum_required(VERSION 3.7) +cmake_minimum_required(VERSION 3.0) project(wolff) -include_directories(src ~/.local/include) +include_directories(lib ~/.local/include) link_directories(~/.local/lib) file(GLOB SOURCES lib/*.c) diff --git a/lib/queue.c b/lib/queue.c index 9a01741..aebbb2a 100644 --- a/lib/queue.c +++ b/lib/queue.c @@ -1,23 +1,23 @@ -#include "wolff.h" +#include "queue.h" void stack_push(ll_t **q, uint32_t x) { - ll_t *nq = malloc(sizeof(ll_t)); - nq->x = x; - nq->next = *q; + ll_t *nq = malloc(sizeof(ll_t)); + nq->x = x; + nq->next = *q; - *q = nq; + *q = nq; } uint32_t stack_pop(ll_t **q) { - ll_t *old_q = *q; + ll_t *old_q = *q; - *q = old_q->next; - uint32_t x = old_q->x; + *q = old_q->next; + uint32_t x = old_q->x; - free(old_q); + free(old_q); - return x; + return x; } bool stack_contains(const ll_t *q, uint32_t x) { @@ -29,4 +29,3 @@ bool stack_contains(const ll_t *q, uint32_t x) { return stack_contains(q->next, x); } } - diff --git a/lib/queue.h b/lib/queue.h new file mode 100644 index 0000000..c7acec1 --- /dev/null +++ b/lib/queue.h @@ -0,0 +1,18 @@ + +#pragma once + +#include +#include +#include +#include + +typedef struct ll_tag { + uint32_t x; + struct ll_tag *next; +} ll_t; + +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); diff --git a/lib/wolff.h b/lib/wolff.h index cec9ee3..d3bc412 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -1,19 +1,23 @@ -#include -#include -#include +#pragma once + +#include +#include #include -#include -#include +#include #include #include +#include +#include #include -#include -#include +#include +#include #include #include +#include "queue.h" + typedef struct { graph_t *g; bool *spins; @@ -21,32 +25,18 @@ typedef struct { double H; } ising_state_t; -typedef struct ll_tag { - uint32_t x; - struct ll_tag *next; -} ll_t; - typedef struct { uint32_t nv; - double dH; - int32_t dM; + int32_t dJb; + int32_t dHb; } 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); +cluster_t *flip_cluster(const graph_t *g, const double *ps, 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); -uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps); - +uint32_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 index 016e19b..f45e1b5 100644 --- a/lib/wolff_tools.c +++ b/lib/wolff_tools.c @@ -1,4 +1,5 @@ +#include "queue.h" #include "wolff.h" graph_t *graph_add_ext(const graph_t *g) { @@ -9,7 +10,7 @@ graph_t *graph_add_ext(const graph_t *g) { 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->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)); @@ -45,42 +46,21 @@ graph_t *graph_add_ext(const graph_t *g) { 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) { +cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, + gsl_rng *r) { uint32_t v0; int32_t n_h_bonds, n_bonds; bool x0; cluster_t *c; - + v0 = gsl_rng_uniform_int(r, g->nv); // pick a random vertex - x0 = x[v0]; // record its orientation + x0 = x[v0]; // record its orientation - ll_t *stack = NULL; // create a new stack + 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; - - n_h_bonds = 0; n_bonds = 0; while (stack != NULL) { uint32_t v; @@ -90,7 +70,7 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, g 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 + x[v] = !x[v]; // flip the vertex for (uint32_t i = 0; i < nn; i++) { uint32_t e, v1, v2, vn; @@ -102,9 +82,11 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, g vn = v == v1 ? v2 : v1; // distinguish neighboring site from site itself - bond_counter = (v1 == g->nv - 1 || v2 == g->nv - 1) ? &n_h_bonds : &n_bonds; + bond_counter = + (v1 == g->nv - 1 || v2 == g->nv - 1) ? &(c->dHb) : &(c->dJb); - if (x[vn] == x0) { // if the neighboring site matches the flipping cluster... + if (x[vn] == + x0) { // if the neighboring site matches the flipping cluster... (*bond_counter)++; if (gsl_rng_uniform(r) < ps[e]) { // and with probability ps[e]... @@ -115,22 +97,15 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, g } } - if (v != g->nv - 1) { + if (v != g->nv - 1) { // count the number of non-external sites that flip c->nv++; } } } - c->dH = n_bonds + H * n_h_bonds; - c->dM = n_h_bonds; - 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 * fabs(H) / T); @@ -151,7 +126,8 @@ double *get_bond_probs(double T, double H, ising_state_t *s) { return ps; } -uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps) { +uint32_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()); @@ -161,10 +137,10 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps ps = get_bond_probs(T, H, s); } - cluster_t *c = flip_cluster(s->g, ps, H, s->spins, r); + cluster_t *c = flip_cluster(s->g, ps, s->spins, r); - s->M += -2 * c->dM; - s->H += 2 * c->dH; + s->M += -2 * c->dHb; + s->H += 2 * (c->dJb + H * c->dHb); uint32_t n_flips = c->nv; @@ -172,4 +148,3 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps return n_flips; } - diff --git a/src/pt_wolff.c b/src/pt_wolff.c deleted file mode 100644 index 128b601..0000000 --- a/src/pt_wolff.c +++ /dev/null @@ -1,317 +0,0 @@ - -#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 index 0c369b7..7675889 100644 --- a/src/wolff.c +++ b/src/wolff.c @@ -1,5 +1,5 @@ -#include "wolff.h" +#include int main(int argc, char *argv[]) { int opt; @@ -22,55 +22,56 @@ int main(int argc, char *argv[]) { while ((opt = getopt(argc, argv, "N:L:T:H:m:e:oq:D")) != -1) { switch (opt) { - case 'N': - N = (uint64_t)atof(optarg); + 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; + case 'D': + use_dual = true; + break; + case 'q': + lattice_i = atoi(optarg); + switch (lattice_i) { + case 0: + lat = SQUARE_LATTICE; break; - case 'L': - L = atoi(optarg); + case 1: + lat = DIAGONAL_LATTICE; break; - case 'T': - T = atof(optarg); + case 2: + lat = TRIANGULAR_LATTICE; break; - case 'H': - H = atof(optarg); + case 3: + lat = VORONOI_HYPERUNIFORM_LATTICE; break; - case 'm': - min_runs = atoi(optarg); + case 4: + lat = VORONOI_LATTICE; break; - case 'e': - eps= atof(optarg); - break; - case 'o': - output_state = true; - break; - case 'D': - use_dual = true; - break; - case 'q': - lattice_i = atoi(optarg); - switch (lattice_i) { - case 0: - lat = SQUARE_LATTICE; - break; - case 1: - lat = DIAGONAL_LATTICE; - break; - case 2: - lat = TRIANGULAR_LATTICE; - break; - case 3: - lat = VORONOI_HYPERUNIFORM_LATTICE; - break; - case 4: - lat = VORONOI_LATTICE; - break; - default: - printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 (DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n"); - exit(EXIT_FAILURE); - } - break; default: + printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 " + "(DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n"); exit(EXIT_FAILURE); + } + break; + default: + exit(EXIT_FAILURE); } } @@ -94,12 +95,26 @@ int main(int argc, char *argv[]) { 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; + 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. || n_runs < min_runs) { - printf("\033[F\033[JWOLFF: sweep %lu, 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); + printf("\033[F\033[JWOLFF: sweep %" PRIu64 + ", 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; @@ -114,19 +129,27 @@ int main(int argc, char *argv[]) { 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(s->M - M1, 2) * 1. / (n_runs + 1.); - Mmu4 = Mmu4 * (n_runs / (n_runs + 1.)) + pow(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.); + Mmu2 = Mmu2 * (n_runs / (n_runs + 1.)) + + pow(s->M - M1, 2) * 1. / (n_runs + 1.); + Mmu4 = Mmu4 * (n_runs / (n_runs + 1.)) + + pow(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; + 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; + 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); @@ -134,14 +157,20 @@ int main(int argc, char *argv[]) { diff = fabs(dX / X); } - clust_per_sweep = clust_per_sweep * (n_runs / (n_runs + 1.)) + (n_clust * 1. / N) * 1. / (n_runs + 1.); + 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 %lu, 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); + printf("\033[F\033[JWOLFF: sweep %" PRIu64 + ", 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); + 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); @@ -164,4 +193,3 @@ int main(int argc, char *argv[]) { return 0; } - diff --git a/src/wolff.h b/src/wolff.h deleted file mode 100644 index cec9ee3..0000000 --- a/src/wolff.h +++ /dev/null @@ -1,52 +0,0 @@ - -#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 { - uint32_t nv; - double dH; - int32_t dM; -} 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); - -uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps); - -- cgit v1.2.3-54-g00ecf