From aa81a529be992b323257cadb5f7b9d4eb4dfb211 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 30 Jul 2019 13:38:14 -0400 Subject: added measurement tools --- ising.cpp | 35 ++++++++--- space_wolff.hpp | 177 ++++++++++++++++++++++++++++++++++++-------------------- 2 files changed, 142 insertions(+), 70 deletions(-) diff --git a/ising.cpp b/ising.cpp index 6939749..645c759 100644 --- a/ising.cpp +++ b/ising.cpp @@ -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& m, unsigned i0, spin s1) -> std::set { std::set nn; if (i0 < m.s.size()) { - std::set os1 = m.dict.on_site(s1.x); - std::set nn0 = m.dict.nearest_neighbors(m.s[i0].x); - std::set nn1 = m.dict.nearest_neighbors(s1.x); + std::set nn0 = m.dict.neighbors(m.s[i0].x, 1); + std::set 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 τ = ising.Eq.τ(); + std::cout << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n"; + if (τ[1] / τ[0] < ε) { + break; + } + } std::vector 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 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 #include #include +#include #include #include #include @@ -66,7 +67,7 @@ class euclidean { } template - spin act(spin s) { + spin act(const spin& s) const { spin 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 tnew = r * x.t + t; matrix rnew = r * x.r; @@ -88,7 +89,7 @@ class euclidean { return pnew; } - euclidean inverse() { + euclidean inverse() const { vector tnew = - r.transpose() * t; matrix rnew = r.transpose(); @@ -98,6 +99,7 @@ class euclidean { } }; +// to-do: clean up nnn... code template class dictionary { private: @@ -108,7 +110,7 @@ class dictionary { dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {}; template - unsigned dictionary_index(vector x) { + unsigned dictionary_index(vector x) const { unsigned pos_ind = 0; for (unsigned i = 0; i < D; i++) { @@ -129,42 +131,25 @@ class dictionary { }; template - std::set on_site(vector x) { - return d[dictionary_index(x)]; + std::set neighbors(vector x, unsigned depth) const { + return nearest_neighbors_of(dictionary_index(x), depth, {}); }; - template - std::set nearest_neighbors(vector x) { - unsigned ind = dictionary_index(x); - std::set ns; + std::set nearest_neighbors_of(unsigned ind, unsigned depth, std::list ignores) const { + std::set 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 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 - std::set next_nearest_neighbors(vector x) { - unsigned ind = dictionary_index(x); - std::set ns; + std::set 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 - std::set next_next_nearest_neighbors(vector x) { - unsigned ind = dictionary_index(x); - std::set ns; +class quantity { + private: + double total; + double total2; + std::vector 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 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 ρ() const { + double C0 = C.front() / (n - wait); + double avg2 = pow(total / (n - wait), 2); + + std::vector ρtmp; + + for (double Ct : C) { + ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2)); + } + + return ρtmp; + } + + std::array τ() const { + double τtmp = 0.5; + unsigned M = 1; + double c = 8.0; + + std::vector ρ_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 @@ -216,15 +249,29 @@ class model { std::function(model&, unsigned, spin)> neighbors; std::function, spin)> Z; std::function)> B; - double E; + long double E; + quantity Eq; + quantity Cq; model(unsigned L, std::function, spin)> Z, std::function)> B, std::function(model&, unsigned, spin)> 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 r, std::mt19937& rng) { + unsigned cluster_size = 0; std::uniform_real_distribution dist(0.0, 1.0); std::queue 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 g(L, t, m); this->step(T, ind_dist(rng), g, rng); + + this->update_energy(); + Eq.add(E); } } }; -- cgit v1.2.3-70-g09d2