diff options
Diffstat (limited to 'space_wolff.hpp')
-rw-r--r-- | space_wolff.hpp | 333 |
1 files changed, 333 insertions, 0 deletions
diff --git a/space_wolff.hpp b/space_wolff.hpp new file mode 100644 index 0000000..aaef673 --- /dev/null +++ b/space_wolff.hpp @@ -0,0 +1,333 @@ + +#include <vector> +#include <list> +#include <set> +#include <iostream> +#include <functional> +#include <random> +#include <queue> +#include <eigen3/Eigen/Dense> +#include "randutils/randutils.hpp" + +const std::array<std::array<unsigned, 16>, 16> smiley = + {{ + {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}}, + {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}}, + {{0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0}}, + {{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}}, + {{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1}}, + {{1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1}}, + {{1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1}}, + {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}, + {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}, + {{1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1}}, + {{1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}}, + {{1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1}}, + {{0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0}}, + {{0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0}}, + {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}}, + {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}} + }}; + +template <class T, int D> +using vector = Eigen::Matrix<T, D, 1>; + +template <class T, int D> +using matrix = Eigen::Matrix<T, D, D>; + +template <class T, int D, class state> +class spin { + public: + vector<T, D> x; + state s; +}; + +template <class T, int D> +class euclidean { + private: + unsigned L; + + public: + vector<T, D> t; + matrix<T, D> r; + euclidean(unsigned L) : L(L) { + for (unsigned i = 0; i < D; i++) { + t(i) = 0; + r(i, i) = 1; + for (unsigned j = 1; j < D; j++) { + r(i, (i + j) % D) = 0; + } + } + } + + euclidean(unsigned L, vector<T, D> t0, matrix<T, D> r0) : L(L) { + t = t0; + r = r0; + } + + template <class state> + spin<T, D, state> act(spin<T, D, state> s) { + spin<T, D, state> s_new; + + s_new.x = t + r * s.x; + s_new.s = s.s; + + for (unsigned i = 0; i < D; i++) { + s_new.x(i) = fmod(L + s_new.x(i), L); + } + + return s_new; + } + + euclidean act(const euclidean& x) { + vector<T, D> tnew = r * x.t + t; + matrix<T, D> rnew = r * x.r; + + euclidean pnew(this->L, tnew, rnew); + + return pnew; + } + + euclidean inverse() { + vector<T, D> tnew = - r.transpose() * t; + matrix<T, D> rnew = r.transpose(); + + euclidean pnew(this->L, tnew, rnew); + + return pnew; + } +}; + +template <int D> +class dictionary { + private: + unsigned L; + std::vector<std::set<unsigned>> d; + + public: + dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {}; + + template <class T> + unsigned dictionary_index(vector<T, D> x) { + unsigned pos_ind = 0; + + for (unsigned i = 0; i < D; i++) { + pos_ind += pow(L, i) * (unsigned)x(i); + }; + + return pos_ind; + } + + template <class T> + void record(vector<T, D> x, unsigned ind) { + d[dictionary_index<T>(x)].insert(ind); + }; + + template <class T> + void remove(vector<T, D> x, unsigned ind) { + d[dictionary_index<T>(x)].erase(ind); + }; + + template <class T> + std::set<unsigned> on_site(vector<T, D> x) { + return d[dictionary_index<T>(x)]; + }; + + template <class T> + std::set<unsigned> nearest_neighbors(vector<T, D> x) { + unsigned ind = dictionary_index<T>(x); + std::set<unsigned> ns; + + for (unsigned i = 0; i < D; i++) { + for (signed j : {-1, 1}) { + unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1)); + for (unsigned nii : d[ni]) { + ns.insert(nii); + } + } + } + + return ns; + }; + + template <class T> + std::set<unsigned> next_nearest_neighbors(vector<T, D> x) { + unsigned ind = dictionary_index<T>(x); + std::set<unsigned> ns; + + for (unsigned i = 0; i < D; i++) { + for (signed j : {-1, 1}) { + unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1)); + for (unsigned k = 0; k < D; k++) { + if (k != i) { + for (signed l : {-1, 1}) { + unsigned nni = pow(L, k + 1) * (ni / ((unsigned)pow(L, k + 1))) + fmod(pow(L, k + 1) + ni + l * pow(L, k), pow(L, k + 1)); + for (unsigned nnii : d[nni]) { + ns.insert(nnii); + } + } + } + } + } + } + + return ns; + }; + + template <class T> + std::set<unsigned> next_next_nearest_neighbors(vector<T, D> x) { + unsigned ind = dictionary_index<T>(x); + std::set<unsigned> ns; + + for (unsigned i = 0; i < D; i++) { + for (signed j : {-1, 1}) { + unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1)); + for (unsigned k = 0; k < D; k++) { + if (k != i) { + for (signed l : {-1, 1}) { + unsigned nni = pow(L, k + 1) * (ni / ((unsigned)pow(L, k + 1))) + fmod(pow(L, k + 1) + ni + l * pow(L, k), pow(L, k + 1)); + for (unsigned m = 0; m < D; m++) { + if (m != i && m != k) { + for (signed n : {-1, 1}) { + unsigned nnni = pow(L, m + 1) * (nni / ((unsigned)pow(L, m + 1))) + fmod(pow(L, m + 1) + nni + n * pow(L, m), pow(L, m + 1)); + for (unsigned nnnii : d[nnni]) { + ns.insert(nnnii); + } + } + } + } + } + } + } + } + } + + return ns; + }; +}; + +template <class U, int D, class state> +class model { + public: + unsigned L; + euclidean<U, D> s0; + std::vector<spin<U, D, state>> s; + dictionary<D> dict; + std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> neighbors; + std::function<double(spin<U, D, state>, spin<U, D, state>)> Z; + std::function<double(spin<U, D, state>)> B; + double E; + + model(unsigned L, std::function<double(spin<U, D, state>, spin<U, D, state>)> Z, + std::function<double(spin<U, D, state>)> B, + std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> ns) : + L(L), s0(L), dict(L), neighbors(ns), Z(Z), B(B) { + } + + void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) { + std::uniform_real_distribution<double> dist(0.0, 1.0); + + std::queue<unsigned> queue; + queue.push(ind); + + std::vector<bool> visited(s.size() + 1, false); + + while (!queue.empty()) { + unsigned i = queue.front(); + queue.pop(); + + if (!visited[i]) { + visited[i] = true; + + bool we_are_ghost = i == s.size(); + + spin<U, D, state> si_new; + euclidean<U, D> s0_new(L); + + if (we_are_ghost) { + s0_new = r.act(s0); + } else { + si_new = r.act(s[i]); + } + + for (unsigned j : neighbors(*this, i, si_new)) { + if (j != i) { + double dE; + bool neighbor_is_ghost = j == s.size(); + + if (we_are_ghost || neighbor_is_ghost) { + spin<U, D, state> s0s_old, s0s_new; + unsigned non_ghost; + + if (neighbor_is_ghost) { + non_ghost = i; + s0s_old = s0.inverse().act(s[i]); + s0s_new = s0.inverse().act(si_new); + } else { + non_ghost = j; + s0s_old = s0.inverse().act(s[j]); + s0s_new = s0_new.inverse().act(s[j]); + } + + dE = B(s0s_old) - B(s0s_new); + } else { + dE = Z(s[i], s[j]) - Z(si_new, s[j]); + } + + double p = 1.0 - exp(-dE / T); + + if (dist(rng) < p) { + queue.push(j); + } + } + } + + if (we_are_ghost) { + s0 = s0_new; + } else { + dict.remove(s[i].x, i); + s[i] = si_new; + dict.record(s[i].x, i); + } + } + } + } + + void wolff(double T, unsigned N, std::mt19937& rng) { + std::uniform_real_distribution<double> t_dist(0, L); + std::uniform_int_distribution<unsigned> r_dist(0, D - 1); + std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1); + + for (unsigned i = 0; i < N; i++) { + vector<U, D> t; + matrix<U, D> m; + for (unsigned j = 0; j < D; j++) { + t(j) = (U)t_dist(rng); + } + + unsigned flip_D1 = r_dist(rng); + unsigned flip_D2 = r_dist(rng); + + for (unsigned j = 0; j < D; j++) { + for (unsigned k = 0; k < D; k++) { + if ((j == flip_D1 && k == flip_D2) || (j == flip_D2 && k == flip_D1)) { + if (flip_D1 <= flip_D2) { + m(j, k) = -1; + } else { + m(j, k) = 1; + } + } else if ((j == k && j != flip_D1) && j != flip_D2) { + m(j, k) = 1; + } else { + m(j, k) = 0; + } + } + } + + euclidean<U, D> g(L, t, m); + + this->step(T, ind_dist(rng), g, rng); + } + } +}; + |