diff options
| -rw-r--r-- | ising.cpp | 35 | ||||
| -rw-r--r-- | space_wolff.hpp | 177 |
2 files changed, 142 insertions, 70 deletions
@@ -8,10 +8,11 @@ int main(int argc, char* argv[]) { unsigned N = 1000; double T = 2.0 / log(1.0 + sqrt(2.0)); double H = 1.0; + double ε = 0.1; int opt; - while ((opt = getopt(argc, argv, "N:L:T:H:")) != -1) { + while ((opt = getopt(argc, argv, "N:L:T:H:e:")) != -1) { switch (opt) { case 'N': N = (unsigned)atof(optarg); @@ -25,6 +26,9 @@ int main(int argc, char* argv[]) { case 'H': H = atof(optarg); break; + case 'e': + ε = atof(optarg); + break; default: exit(1); } @@ -67,12 +71,10 @@ int main(int argc, char* argv[]) { [] (model<signed, D, signed>& m, unsigned i0, spin<signed, D, signed> s1) -> std::set<unsigned> { std::set<unsigned> nn; if (i0 < m.s.size()) { - std::set<unsigned> os1 = m.dict.on_site(s1.x); - std::set<unsigned> nn0 = m.dict.nearest_neighbors(m.s[i0].x); - std::set<unsigned> nn1 = m.dict.nearest_neighbors(s1.x); + std::set<unsigned> nn0 = m.dict.neighbors(m.s[i0].x, 1); + std::set<unsigned> nn1 = m.dict.neighbors(s1.x, 1); nn.insert(nn0.begin(), nn0.end()); nn.insert(nn1.begin(), nn1.end()); - nn.insert(os1.begin(), os1.end()); nn.insert(m.s.size()); } else { for (unsigned i = 0; i < m.s.size(); i++) { @@ -119,7 +121,16 @@ int main(int argc, char* argv[]) { } */ - ising.wolff(T, N, rng); + ising.update_energy(); + + while (true) { + ising.wolff(T, N, rng); + std::array<double, 2> τ = ising.Eq.τ(); + std::cout << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n"; + if (τ[1] / τ[0] < ε) { + break; + } + } std::vector<signed> output(pow(L, D)); @@ -128,13 +139,21 @@ int main(int argc, char* argv[]) { output[L * rs.x(1) + rs.x(0)] = s.s; } + std::ofstream outfile; + outfile.open("snap.dat"); + for (unsigned i = 0; i < L; i++) { for (unsigned j = 0; j < L; j++) { unsigned out = output[L * i + j] == 1 ? 1 : 0; - std::cout << out; + outfile << out << " "; } - std::cout << "\n"; + outfile << "\n"; } + outfile.close(); + + std::array<double, 2> act = ising.Eq.τ(); + + std::cout << act[0] << " " << act[1] << "\n"; return 0; } diff --git a/space_wolff.hpp b/space_wolff.hpp index aaef673..b60e274 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -3,6 +3,7 @@ #include <list> #include <set> #include <iostream> +#include <fstream> #include <functional> #include <random> #include <queue> @@ -66,7 +67,7 @@ class euclidean { } template <class state> - spin<T, D, state> act(spin<T, D, state> s) { + spin<T, D, state> act(const spin<T, D, state>& s) const { spin<T, D, state> s_new; s_new.x = t + r * s.x; @@ -79,7 +80,7 @@ class euclidean { return s_new; } - euclidean act(const euclidean& x) { + euclidean act(const euclidean& x) const { vector<T, D> tnew = r * x.t + t; matrix<T, D> rnew = r * x.r; @@ -88,7 +89,7 @@ class euclidean { return pnew; } - euclidean inverse() { + euclidean inverse() const { vector<T, D> tnew = - r.transpose() * t; matrix<T, D> rnew = r.transpose(); @@ -98,6 +99,7 @@ class euclidean { } }; +// to-do: clean up nnn... code template <int D> class dictionary { private: @@ -108,7 +110,7 @@ class dictionary { dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {}; template <class T> - unsigned dictionary_index(vector<T, D> x) { + unsigned dictionary_index(vector<T, D> x) const { unsigned pos_ind = 0; for (unsigned i = 0; i < D; i++) { @@ -129,42 +131,25 @@ class dictionary { }; template <class T> - std::set<unsigned> on_site(vector<T, D> x) { - return d[dictionary_index<T>(x)]; + std::set<unsigned> neighbors(vector<T, D> x, unsigned depth) const { + return nearest_neighbors_of(dictionary_index<T>(x), depth, {}); }; - template <class T> - std::set<unsigned> nearest_neighbors(vector<T, D> x) { - unsigned ind = dictionary_index<T>(x); - std::set<unsigned> ns; + std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth, std::list<unsigned> ignores) const { + std::set<unsigned> ns = d[ind]; - 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; - }; + if (depth > 0) { + for (unsigned i = 0; i < D; i++) { + if (std::none_of(ignores.begin(), ignores.end(), [i](unsigned k){return i==k;})) { + std::list<unsigned> ignores_new = ignores; + ignores_new.push_back(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)); - template <class T> - std::set<unsigned> next_nearest_neighbors(vector<T, D> x) { - unsigned ind = dictionary_index<T>(x); - std::set<unsigned> ns; + std::set<unsigned> nns = nearest_neighbors_of(ni, depth - 1, ignores_new); - 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); - } + for (unsigned guy : nns) { + ns.insert(guy); } } } @@ -173,37 +158,85 @@ class dictionary { 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; +class quantity { + private: + double total; + double total2; + std::vector<double> C; + unsigned wait; - 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); - } - } - } - } - } - } + public: + unsigned n; + std::list<double> hist; + + quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) { + n = 0; + total = 0; + total2 = 0; + } + + void add(double x) { + if (n > wait) { + hist.push_front(x); + if (hist.size() > C.size()) { + hist.pop_back(); + unsigned t = 0; + for (double a : hist) { + C[t] += a * x; + t++; } + total += x; + total2 += pow(x, 2); } } + n++; + } - return ns; - }; + double avg() const { + return total / (n - wait); + } + + double avg2() const { + return total2 / (n - wait); + } + + std::vector<double> ρ() const { + double C0 = C.front() / (n - wait); + double avg2 = pow(total / (n - wait), 2); + + std::vector<double> ρtmp; + + for (double Ct : C) { + ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2)); + } + + return ρtmp; + } + + std::array<double, 2> τ() const { + double τtmp = 0.5; + unsigned M = 1; + double c = 8.0; + + std::vector<double> ρ_tmp = this->ρ(); + + while (c * τtmp > M && M < C.size()) { + τtmp += ρ_tmp[M]; + M++; + } + + return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)}; + } + + double σ() const { + return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2)); + } + + double serr() const { + return sqrt(this->σ()); + } }; template <class U, int D, class state> @@ -216,15 +249,29 @@ class model { 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; + long double E; + quantity Eq; + quantity Cq; 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) { + L(L), s0(L), dict(L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) { } + void update_energy() { + E = 0; + for (unsigned i = 0; i < s.size(); i++) { + for (unsigned j = 0; j < i; j++) { + E -= Z(s[i], s[j]); + } + + E -= B(s0.inverse().act(s[i])); + } + } + void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) { + unsigned cluster_size = 0; std::uniform_real_distribution<double> dist(0.0, 1.0); std::queue<unsigned> queue; @@ -288,9 +335,12 @@ class model { dict.remove(s[i].x, i); s[i] = si_new; dict.record(s[i].x, i); + cluster_size++; } } } + + Cq.add(cluster_size); } void wolff(double T, unsigned N, std::mt19937& rng) { @@ -327,6 +377,9 @@ class model { euclidean<U, D> g(L, t, m); this->step(T, ind_dist(rng), g, rng); + + this->update_energy(); + Eq.add(E); } } }; |