diff options
Diffstat (limited to 'lib/cluster.h')
-rw-r--r-- | lib/cluster.h | 111 |
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; -} - |