diff options
| -rw-r--r-- | lib/wolff.h | 1 | ||||
| -rw-r--r-- | lib/wolff_tools.c | 13 | ||||
| -rw-r--r-- | src/wolff.c | 19 |
3 files changed, 11 insertions, 22 deletions
diff --git a/lib/wolff.h b/lib/wolff.h index 6beb03e..ca7635a 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -29,6 +29,7 @@ typedef struct ll_tag { typedef struct { int32_t nv; double dH; + int32_t dM; } cluster_t; double get_hamiltonian(graph_t *g, double *coupling, bool *x); diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c index adef40f..d361344 100644 --- a/lib/wolff_tools.c +++ b/lib/wolff_tools.c @@ -123,10 +123,7 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, g } c->dH = n_bonds + H * n_h_bonds; - - if (x0) { - c->nv = -c->nv; - } + c->dM = n_h_bonds; return c; } @@ -137,7 +134,7 @@ double hh(double th) { 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 q = 1 - exp(-2 * fabs(H) / T); double *ps = (double *)malloc(s->g->ne * sizeof(double)); @@ -167,13 +164,11 @@ int32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps) cluster_t *c = flip_cluster(s->g, ps, H, s->spins, r); - s->M += 2 * c->nv; + s->M += -2 * c->dM; s->H += 2 * c->dH; - int32_t n_flipped = c->nv; - free(c); - return n_flipped; + return c->nv; } diff --git a/src/wolff.c b/src/wolff.c index 22e6857..5f06665 100644 --- a/src/wolff.c +++ b/src/wolff.c @@ -56,7 +56,7 @@ int main(int argc, char *argv[]) { s->g = graph_add_ext(h); s->spins = (bool *)calloc(h->nv + 1, sizeof(bool)); - s->M = -h->nv; + s->M = h->nv; s->H = -(1.0 * h->ne) - H * h->nv; double *bond_probs = get_bond_probs(T, H, s); @@ -82,24 +82,17 @@ int main(int argc, char *argv[]) { n_clust++; } - int32_t MM; - if (s->spins[h->nv]) { - MM = s->M; - } else { - MM = -s->M; - } - E1 = E1 * (n_runs / (n_runs + 1.)) + s->H * 1. / (n_runs + 1.); - M1 = M1 * (n_runs / (n_runs + 1.)) + MM * 1. / (n_runs + 1.); + M1 = M1 * (n_runs / (n_runs + 1.)) + 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(MM, 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(MM - M1, 2) * 1. / (n_runs + 1.); - Mmu4 = Mmu4 * (n_runs / (n_runs + 1.)) + pow(MM - M1, 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){ + 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; |