summaryrefslogtreecommitdiff
path: root/lib/cluster.h
diff options
context:
space:
mode:
Diffstat (limited to 'lib/cluster.h')
-rw-r--r--lib/cluster.h111
1 files changed, 0 insertions, 111 deletions
diff --git a/lib/cluster.h b/lib/cluster.h
deleted file mode 100644
index f57bb68..0000000
--- a/lib/cluster.h
+++ /dev/null
@@ -1,111 +0,0 @@
-
-#pragma once
-
-#include <gsl/gsl_randist.h>
-#include <gsl/gsl_rng.h>
-#include <cmath>
-#include <vector>
-#include <stack>
-
-#include "types.h"
-#include "state.h"
-#include "graph.h"
-
-template <class R_t, class X_t>
-void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, gsl_rng *rand) {
- v_t nv = 0;
-
- std::stack<v_t> stack;
- stack.push(v0);
-
- std::vector<bool> marks(s.g.nv, false);
-
- while (!stack.empty()) {
- v_t v = stack.top();
- stack.pop();
-
- if (!marks[v]) { // don't reprocess anyone we've already visited!
- marks[v] = true;
-
- X_t si_new;
-#ifndef NOFIELD
- R_t R_new;
-
- bool v_is_ghost = (v == s.nv); // ghost site has the last index
-
- if (v_is_ghost) {
- R_new = r.act(s.R); // if we are, then we're moving the transformation
- } else
-#endif
- {
- si_new = r.act(s.spins[v]); // otherwise, we're moving the spin at our site
- }
-
- for (const v_t &vn : s.g.v_adj[v]) {
- double dE, prob;
-
-#ifndef NOFIELD
- bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost
-
- if (v_is_ghost || vn_is_ghost) { // this is a ghost-involved bond
- X_t rs_old, rs_new;
- v_t non_ghost;
-
- if (vn_is_ghost) {
- // if our neighbor is the ghost, the current site is a normal
- // spin - rotate it back!
- rs_old = s.R.act_inverse(s.spins[v]);
- rs_new = s.R.act_inverse(si_new);
- non_ghost = v;
- } else {
- /* if we're the ghost, we need to rotate our neighbor back in
- both the old and new ways */
- rs_old = s.R.act_inverse(s.spins[vn]);
- rs_new = R_new.act_inverse(s.spins[vn]);
- non_ghost = vn;
- }
-
- dE = s.H(rs_old) - s.H(rs_new);
-
-#ifdef FINITE_STATES
- prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)];
-#endif
-
- s.update_magnetization(rs_old, rs_new);
- s.update_fourierZero(non_ghost, rs_old, rs_new);
- } else // this is a perfectly normal bond!
-#endif
- {
- dE = s.J(s.spins[v], s.spins[vn]) - s.J(si_new, s.spins[vn]);
-
-#ifdef FINITE_STATES
- prob = J_probs[state_to_ind(s.spins[v])][state_to_ind(si_new)][state_to_ind(s.spins[vn])];
-#endif
- }
-
- s.update_energy(dE);
-
-#ifndef FINITE_STATES
- prob = 1.0 - exp(-dE / s.T);
-#endif
-
- if (gsl_rng_uniform(rand) < prob) {
- stack.push(vn); // push the neighboring vertex to the stack
- }
- }
-
-#ifndef NOFIELD
- if (v_is_ghost) {
- s.R = R_new;
- } else
-#endif
- {
- s.spins[v] = si_new;
- nv++;
- }
- }
- }
-
- s.last_cluster_size = nv;
-}
-