diff options
| -rw-r--r-- | .gitignore | 7 | ||||
| -rw-r--r-- | ising.cpp | 250 | ||||
| -rw-r--r-- | octree.hpp | 19 | ||||
| -rw-r--r-- | random.hpp | 6 | ||||
| -rw-r--r-- | space_wolff.hpp | 172 | ||||
| -rw-r--r-- | spheres_infinite.cpp | 46 | ||||
| -rw-r--r-- | spin.hpp | 2 | ||||
| -rw-r--r-- | torus_symmetries.hpp | 2 | ||||
| -rw-r--r-- | transformation.hpp | 201 |
9 files changed, 412 insertions, 293 deletions
diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..dac4c22 --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +* +!*.hpp +!*.cpp +!.gitignore +!.gitmodules +*/ + @@ -1,14 +1,53 @@ +#include "smiley.hpp" #include "space_wolff.hpp" +#include "torus_symmetries.hpp" +#include <GL/glut.h> + +const unsigned D = 2; +typedef Model<signed, D, TorusGroup<signed, D>, signed> model; + +class animation : public measurement<signed, D, TorusGroup<signed, D>, signed> { +private: + uint64_t t1; + uint64_t t2; + unsigned n; + unsigned tmp; + +public: + animation(double L, unsigned w, int argc, char* argv[]) { + t1 = 0; + t2 = 0; + n = 0; + + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(w, w); + glutCreateWindow("wolff"); + glClearColor(0.0, 0.0, 0.0, 0.0); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(0, L, 0, L); + } -std::function<double(spin<signed, 2, signed>)> B_sin(unsigned L, unsigned n, double H) { - return [n, H, L] (spin<signed, 2, signed> s) -> double { - return H * s.s * cos(2 * M_PI * n * s.x(0) / ((double)L)); - }; -} + void post_cluster(const model& m) override { + glClearColor(1.0f, 1.0f, 1.0f, 1.0f); + glClear(GL_COLOR_BUFFER_BIT); + glClearColor(0.0, 1.0, 0.0, 0.0); + for (const Spin<signed, 2, signed>* s : m.s) { + if (s->s == 1) { + glColor3f(0.0, 0.0, 0.0); + } else if (s->s == -1) { + glColor3f(1.0, 1.0, 1.0); + } + Vector<signed, 2> xx = m.s0.inverse().act(*s).x; + glRecti(xx(0), xx(1), xx(0) + 1, xx(1) + 1); + } + glFlush(); + } +}; int main(int argc, char* argv[]) { - const unsigned D = 2; unsigned L = 32; unsigned N = 1000; @@ -23,120 +62,133 @@ int main(int argc, char* argv[]) { while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:n:")) != -1) { switch (opt) { - case 'N': - N = (unsigned)atof(optarg); - break; - case 'L': - L = atoi(optarg); - break; - case 'T': - T = atof(optarg); - break; - case 'H': - H = atof(optarg); - break; - case 'e': - ε = atof(optarg); - break; - case 'm': - mod = atoi(optarg); - break; - case 'M': - multi = atoi(optarg); - break; - case 'n': - mag = atof(optarg); - break; - default: - exit(1); + case 'N': + N = (unsigned)atof(optarg); + break; + case 'L': + L = atoi(optarg); + break; + case 'T': + T = atof(optarg); + break; + case 'H': + H = atof(optarg); + break; + case 'e': + ε = atof(optarg); + break; + case 'm': + mod = atoi(optarg); + break; + case 'M': + multi = atoi(optarg); + break; + case 'n': + mag = atof(optarg); + break; + default: + exit(1); } } double pZ = 1.0 - exp(-1.0 / T); - std::function<double(spin<signed, D, signed>, spin<signed, D, signed>)> Z = - [L, pZ] (spin<signed, D, signed> s1, spin<signed, D, signed> s2) -> double { - bool old_one_one = false; - bool old_many_ones = false; - bool old_any_two = false; - - vector<signed, D> old_diff = diff<signed, D>(L, s1.x, s2.x); - - for (unsigned i = 0; i < D; i++) { - if (old_diff(i) == 1 && !old_one_one) { - old_one_one = true; - } else if (old_diff(i) == 1 && old_one_one) { - old_many_ones = true; - } else if (old_diff(i) > 1) { - old_any_two = true; - } + std::function<double(const Spin<signed, D, signed>&, const Spin<signed, D, signed>&)> Z = + [L, pZ](const Spin<signed, D, signed>& s1, const Spin<signed, D, signed>& s2) -> double { + bool old_one_one = false; + bool old_many_ones = false; + bool old_any_two = false; + + Vector<signed, D> old_diff = diff<signed, D>(L, s1.x, s2.x); + + for (unsigned i = 0; i < D; i++) { + if (old_diff(i) == 1 && !old_one_one) { + old_one_one = true; + } else if (old_diff(i) == 1 && old_one_one) { + old_many_ones = true; + } else if (old_diff(i) > 1) { + old_any_two = true; } + } - bool were_on_someone = !old_one_one && !old_any_two; - bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two; - - if (were_on_someone) { - return -std::numeric_limits<double>::infinity();; - } else if (were_nearest_neighbors) { - return s1.s * s2.s; - } else { - return 0.0; - } - }; - - std::function<double(spin<signed, D, signed>)> B_face = - [L, H] (spin<signed, D, signed> s) -> double { - return H * s.s * smiley[s.x(0) * 16 / L][s.x(1) * 16 / L]; - }; - - std::function<double(spin<signed, D, signed>)> B; + bool were_on_someone = !old_one_one && !old_any_two; + bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two; - if (mod > 0) { - B = B_sin(L, mod, H); - } else { - B = B_face; - } + if (were_on_someone) { + return -std::numeric_limits<double>::infinity(); + ; + } else if (were_nearest_neighbors) { + return s1.s * s2.s; + } else { + return 0.0; + } + }; - std::function<std::set<unsigned>(model<signed, D, signed>&, unsigned, spin<signed, D, signed>)> neighbors = - [] (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> 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(m.s.size()); - } else { - for (unsigned i = 0; i < m.s.size(); i++) { - nn.insert(i); - } - } - return nn; - }; + std::function<double(const Spin<signed, D, signed>&)> B_face = + [L, H](const Spin<signed, D, signed>& s) -> double { + return H * s.s * smiley[15 - s.x(1) * 16 / L][s.x(0) * 16 / L]; + }; - model<signed, D, signed> ising(L, L, Z, B, neighbors); + std::function<double(Spin<signed, D, signed>)> B = B_face; - randutils::auto_seed_128 seeds; - std::mt19937 rng{seeds}; + Model<signed, D, TorusGroup<signed, D>, signed> ising(L, Z, B); - std::uniform_int_distribution<unsigned> coin(0, 1); + Rng rng; unsigned n = 0; unsigned up = 0; unsigned down = 0; for (unsigned i = 0; i < L; i++) { for (unsigned j = 0; j < L; j++) { - if ((coin(rng) && up < pow(L, 2) * mag) || down >= pow(L, 2) * mag) { - ising.s.push_back({{i, j}, 1}); + Spin<signed, D, signed>* ss = new Spin<signed, D, signed>(); + ss->x = {i, j}; + if (rng.uniform<double>(0, 1) < mag) { + ss->s = 1; up++; } else { - ising.s.push_back({{i, j}, -1}); + ss->s = -1; down++; } - ising.dict.record({i, j}, n); + ising.s.push_back(ss); + ising.dict.insert(ss); n++; } } + + std::vector<Vector<signed, 2>> torusVectors = torus_vecs<signed, 2>(L); + std::vector<Matrix<signed, 2>> torusMatrices = torus_mats<signed, 2>(); + + Gen<signed, D, TorusGroup<signed, D>, signed> g = + [L, &torusVectors, + &torusMatrices](Model<signed, D, TorusGroup<signed, D>, signed>& M, + Rng& r) -> Transformation<signed, D, TorusGroup<signed, D>, signed>* { + Matrix<signed, 2> m; + Vector<signed, 2> t; + + m = r.pick(torusMatrices); + t(0) = r.uniform<signed>(0, L); + t(1) = r.uniform<signed>(0, L); + t = t - m * t; + + TorusGroup<signed, 2> g = TorusGroup<signed, 2>({(signed)L, t, m}); + + Spin<signed, 2, signed>* ss = r.pick(M.s); + Spin<signed, 2, signed> s2 = g.act(*ss); + + while (ss->x(0) == s2.x(0) && ss->x(1) == s2.x(1)) { + ss = r.pick(M.s); + s2 = g.act(*ss); + } + + std::set<Spin<signed, 2, signed>*> s2s = M.dict.at(s2.x); + + return new PairFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss, *s2s.begin()); + }; + + animation A(L, 750, argc, argv); + + ising.wolff(T, {g}, A, N); + /* for (unsigned i = 0; i < L; i++) { for (unsigned j = 0; j < L; j++) { @@ -151,8 +203,7 @@ int main(int argc, char* argv[]) { } */ - ising.update_energy(); - + /* while (true) { ising.wolff(T, N, rng); std::array<double, 2> τ = ising.Eq.τ(); @@ -168,8 +219,8 @@ int main(int argc, char* argv[]) { std::array<double, 2> act = ising.Eq.τ(); std::vector<double> ρ = ising.Eq.ρ(); - outfile << L << " " << T << " " << mod << " " << H << " " << ising.Eq.num_added() << " " << ising.Cq.avg() << " " << ising.Cq.serr() << " " << act[0] << " " << act[1]; - for (double ρi : ρ) { + outfile << L << " " << T << " " << mod << " " << H << " " << ising.Eq.num_added() << " " << + ising.Cq.avg() << " " << ising.Cq.serr() << " " << act[0] << " " << act[1]; for (double ρi : ρ) { outfile << " " << ρi; } outfile << "\n"; @@ -195,6 +246,7 @@ int main(int argc, char* argv[]) { } snapfile.close(); + */ return 0; } @@ -5,6 +5,8 @@ #include <set> #include <unordered_map> +#include <iostream> + #include "spin.hpp" #include "vector.hpp" @@ -31,16 +33,16 @@ private: std::unordered_map<std::array<signed, D>, std::set<Spin<T, D, S>*>> data; public: - Octree(T L, unsigned N) { - L = L; - N = N; + Octree(T L_tmp, unsigned N_tmp) { + L = L_tmp; + N = N_tmp; } std::array<signed, D> ind(Vector<T, D> x) const { std::array<signed, D> ind; for (unsigned i = 0; i < D; i++) { - ind[i] = (signed)std::floor(N * x(i) / L); + ind[i] = std::floor(N * x(i) / L); } return ind; @@ -55,6 +57,15 @@ public: } }; + std::set<Spin<T, D, S>*> at(const Vector<T, D>& x) const { + auto it = data.find(ind(x)); + if (it == data.end()) { + return {}; + } else { + return it->second; + } + } + std::set<Spin<T, D, S>*> neighbors(const Vector<T, D>& x) const { std::array<signed, D> i0 = ind(x); std::set<Spin<T, D, S>*> ns; diff --git a/random.hpp b/random.hpp new file mode 100644 index 0000000..9f1d33d --- /dev/null +++ b/random.hpp @@ -0,0 +1,6 @@ +#pragma once + +#include "pcg-cpp/include/pcg_random.hpp" +#include "randutils/randutils.hpp" + +using Rng = randutils::random_generator<pcg32>; diff --git a/space_wolff.hpp b/space_wolff.hpp index 8ea7dfb..c4a7a21 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -7,15 +7,12 @@ #include <queue> #include <vector> -#include "pcg-cpp/include/pcg_random.hpp" -#include "randutils/randutils.hpp" - #include "euclidean.hpp" #include "octree.hpp" #include "quantity.hpp" #include "spin.hpp" - -using Rng = randutils::random_generator<pcg32>; +#include "random.hpp" +#include "transformation.hpp" template <class U, int D, class R, class S> class Model; @@ -35,175 +32,12 @@ public: }; -template <class U, int D, class R, class S> class Transformation { - public: - virtual std::set<Spin<U, D, S>*> current() const {return {NULL};} - virtual std::set<Spin<U, D, S>*> toConsider() const {return {};} - virtual double ΔE(const Spin<U, D, S>*) const {return 0;} - virtual Transformation* createNew(Spin<U, D, S>*) const {return new Transformation();} - virtual void apply(); -}; - -template <class U, int D, class R, class S> using Gen = std::function<Transformation<U, D, R, S>*(Model<U, D, R, S>&, Rng&)>; - -template <class U, int D, class R, class S> class SpinFlip : public Transformation<U, D, R, S> { - private: - Model<U, D, R, S>& M; - Spin<U, D, S>& sOld; - const R r; - Spin<U, D, S> sNew; - public: - SpinFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>& s) : M(M), r(r), sOld(s) { - sNew = r.act(s); - } - - std::set<Spin<U, D, S>*> current() const override { - return {&sOld}; - } - - std::set<Spin<U, D, S>*> toConsider() const override { - std::set<Spin<U, D, S>*> neighbors; - std::set<Spin<U, D, S>*> current_neighbors = M.dict.neighbors(sOld.x); - std::set<Spin<U, D, S>*> new_neighbors = M.dict.neighbors(sNew.x); - - neighbors.insert(current_neighbors.begin(), current_neighbors.end()); - neighbors.insert(new_neighbors.begin(), new_neighbors.end()); - neighbors.insert(NULL); - - return neighbors; - } - - double ΔE(const Spin<U, D, S>* s) const override { - if (s == NULL) { - Spin<U, D, S> s0s_old = M.s0.inverse().act(sOld); - Spin<U, D, S> s0s_new = M.s0.inverse().act(sNew); - return M.B(s0s_new) - M.B(s0s_old); - } else { - return M.Z(sOld, *s) - M.Z(sNew, *s); - } - } - - Transformation<U, D, R, S>* createNew(Spin<U, D, S>* s) const override; - - void apply() override { - M.dict.remove(&sOld); - sOld = sNew; - M.dict.insert(&sOld); - } -}; - -template <class U, int D, class R, class S> class PairFlip : public Transformation<U, D, R, S> { - private: - Model<U, D, R, S>& M; - Spin<U, D, S>& s1Old; - Spin<U, D, S>& s2Old; - const R r; - Spin<U, D, S> s1New; - Spin<U, D, S> s2New; - public: - PairFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>& s1, Spin<U, D, S>& s2) : M(M), r(r), s1Old(s1), s2Old(s2) { - s1New = r.act(s1); - s2New = r.act(s2); - } - - std::set<Spin<U, D, S>*> current() const override { - return {&s1Old, &s2Old}; - } - - std::set<Spin<U, D, S>*> toConsider() const override { - std::set<Spin<U, D, S>*> neighbors; - std::set<Spin<U, D, S>*> current_neighbors_1 = M.dict.neighbors(s1Old.x); - std::set<Spin<U, D, S>*> current_neighbors_2 = M.dict.neighbors(s2Old.x); - std::set<Spin<U, D, S>*> new_neighbors_1 = M.dict.neighbors(s1New.x); - std::set<Spin<U, D, S>*> new_neighbors_2 = M.dict.neighbors(s2New.x); - - neighbors.insert(current_neighbors_1.begin(), current_neighbors_1.end()); - neighbors.insert(current_neighbors_2.begin(), current_neighbors_2.end()); - neighbors.insert(new_neighbors_1.begin(), new_neighbors_1.end()); - neighbors.insert(new_neighbors_2.begin(), new_neighbors_2.end()); - neighbors.insert(NULL); - - return neighbors; - } - - double ΔE(const Spin<U, D, S>* s) const override { - if (s == NULL) { - Spin<U, D, S> s0s1_old = M.s0.inverse().act(s1Old); - Spin<U, D, S> s0s1_new = M.s0.inverse().act(s1New); - Spin<U, D, S> s0s2_old = M.s0.inverse().act(s2Old); - Spin<U, D, S> s0s2_new = M.s0.inverse().act(s2New); - return M.B(s0s1_new) + M.B(s0s2_new) - M.B(s0s1_old) - M.B(s0s2_old); - } else { - return M.Z(s1Old, *s) + M.Z(s2Old, *s) - M.Z(s1New, *s) - M.Z(s2New, *s); - } - } - - Transformation<U, D, R, S>* createNew(Spin<U, D, S>* s) const override; - - void apply() override { - M.dict.remove(&s1Old); - M.dict.remove(&s2Old); - s1Old = s1New; - s2Old = s2New; - M.dict.insert(&s1Old); - M.dict.insert(&s2Old); - } -}; - -template <class U, int D, class R, class S> class FieldFlip : public Transformation<U, D, R, S> { - private: - Model<U, D, R, S>& M; - const R r; - R s0New; - public: - FieldFlip(Model<U, D, R, S>& M, const R& r) : M(M), r(r), s0New(r.act(M.s0)) {} - - std::set<Spin<U, D, S>*> toConsider() const override { - std::set<Spin<U, D, S>*> neighbors; - - for (Spin<U, D, S>& s : M.s) { - neighbors.insert(&s); - } - - return neighbors; - } - - double ΔE(const Spin<U, D, S>* s) const override { - Spin<U, D, S> s0s_old = M.s0.inverse().act(*s); - Spin<U, D, S> s0s_new = s0New.inverse().act(*s); - return M.B(s0s_new) - M.B(s0s_old); - } - - Transformation<U, D, R, S>* createNew(Spin<U, D, S>* s) const override { - return new SpinFlip<U, D, R, S>(M, r, *s); - } - - void apply() override { - M.s0 = s0New; - } -}; - -template<class U, int D, class R, class S> Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(Spin<U, D, S>* s) const { - if (s == NULL) { - return new FieldFlip<U, D, R, S>(M, r); - } else { - return new SpinFlip(M, r, *s); - } -} - -template<class U, int D, class R, class S> Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<U, D, S>* s) const { - if (s == NULL) { - return new FieldFlip<U, D, R, S>(M, r); - } else { - return new SpinFlip<U, D, R, S>(M, r, *s); - } -} template <class U, int D, class R, class S> class Model { public: U L; R s0; - std::vector<Spin<U, D, S>> s; + std::vector<Spin<U, D, S>*> s; Octree<U, D, S> dict; std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z; std::function<double(const Spin<U, D, S>&)> B; diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp index 3505d6d..026486e 100644 --- a/spheres_infinite.cpp +++ b/spheres_infinite.cpp @@ -41,13 +41,13 @@ public: void post_cluster(const model& m) override { glClearColor(1.0f, 1.0f, 1.0f, 1.0f); glClear(GL_COLOR_BUFFER_BIT); - for (const Spin<double, 2, double>& s : m.s) { + for (const Spin<double, 2, double>* s : m.s) { glBegin(GL_POLYGON); unsigned n_points = 50; glColor3f(0.0f, 0.0f, 0.0f); for (unsigned i = 0; i < n_points; i++) { - glVertex2d(m.s0.inverse().act(s).x(0) + s.s * cos(2 * i * M_PI / n_points), - m.s0.inverse().act(s).x(1) + s.s * sin(2 * i * M_PI / n_points)); + glVertex2d(m.s0.inverse().act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), + m.s0.inverse().act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); } glEnd(); } @@ -78,14 +78,14 @@ Gen<double, D, Euclidean<double, D>, double> eGen(double ε) { m(0, 1) = -2 * cos(θ) * sin(θ); m(1, 0) = -2 * cos(θ) * sin(θ); - unsigned f_ind = rng.uniform((unsigned)0, (unsigned)M.s.size()); - t = M.s[f_ind].x; + Spin<double, D, double>* s = rng.pick(M.s); + t = s->x; for (unsigned j = 0; j < D; j++) { t(j) += rng.variate<double, std::normal_distribution>(0.0, ε); } Euclidean<double, D> g(t - m * t, m); - return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, M.s[f_ind]); + return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, s); }; } @@ -93,12 +93,15 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) { return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* { Matrix<double, D> m; - unsigned f_ind1 = rng.uniform((unsigned)0, (unsigned)M.s.size()); - unsigned f_ind2 = rng.uniform((unsigned)0, (unsigned)M.s.size() - 1); - if (f_ind2 >= f_ind1) f_ind2++; + Spin<double, D, double>* s1 = rng.pick(M.s); + Spin<double, D, double>* s2 = rng.pick(M.s); - Vector<double, D> t1 = M.s[f_ind1].x; - Vector<double, D> t2 = M.s[f_ind2].x; + while (s1 == s2) { + s2 = rng.pick(M.s); + } + + Vector<double, D> t1 = s1->x; + Vector<double, D> t2 = s2->x; Vector<double, D> t12 = t1 - t2; Vector<double, D> t = (t1 + t2) / 2; @@ -110,7 +113,7 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) { m(1, 0) = -2 * cos(θ) * sin(θ); Euclidean<double, D> g(t - m * t, m); - return new PairFlip<double, D, Euclidean<double, D>, double>(M, g, M.s[f_ind1], M.s[f_ind2]); + return new PairFlip<double, D, Euclidean<double, D>, double>(M, g, s1, s2); }; } @@ -179,17 +182,19 @@ int main(int argc, char* argv[]) { auto g1 = eGen(0.5); auto g2 = mGen(0.2); animation A(L, 750, argc, argv); - model sphere(1, Z, B); + model sphere(1.0, Z, B); Rng rng; - sphere.s.reserve(n); + sphere.s.resize(n); unsigned nx = floor(sqrt(n)); - for (unsigned i = 0; i < n; i++) { - Vector<double, D> pos = {(i / nx) * L / nx, (i % nx) * L / nx}; - sphere.s.push_back({pos, rng.uniform<double>(0.25, 0.45)}); - sphere.dict.insert(&sphere.s.back()); + for (unsigned i = 0; i < sphere.s.size(); i++) { + Spin<double, 2, double>* ss = new Spin<double, 2, double>(); + ss->x = {(i / nx) * L / nx, (i % nx) * L / nx}; + ss->s = rng.uniform<double>(0.25, 0.45); + sphere.s[i] = ss; + sphere.dict.insert(ss); } sphere.wolff(T, {g1, g2}, A, N); @@ -212,9 +217,10 @@ int main(int argc, char* argv[]) { std::ofstream snapfile; snapfile.open("sphere_snap.dat"); - for (Spin<double, D, double> s : sphere.s) { - Spin<double, D, double> rs = sphere.s0.inverse().act(s); + for (Spin<double, D, double>* s : sphere.s) { + Spin<double, D, double> rs = sphere.s0.inverse().act(*s); snapfile << rs.x.transpose() << "\n"; + delete s; } snapfile.close(); @@ -7,5 +7,7 @@ template <class T, int D, class S> class Spin { public: Vector<T, D> x; S s; + + Spin() { x = Vector<T, D>::Zero(); } }; diff --git a/torus_symmetries.hpp b/torus_symmetries.hpp index 17772de..7bcb70b 100644 --- a/torus_symmetries.hpp +++ b/torus_symmetries.hpp @@ -56,7 +56,7 @@ template <class U, int D> std::vector<Matrix<U, D>> torus_mats() { one_sequences<D>(sequences, D); - sequences.pop_front(); // don't want the identity matrix! + sequences.pop_back(); // don't want the identity matrix! for (std::array<double, D> sequence : sequences) { Matrix<U, D> m; diff --git a/transformation.hpp b/transformation.hpp new file mode 100644 index 0000000..0731ec9 --- /dev/null +++ b/transformation.hpp @@ -0,0 +1,201 @@ +#pragma once + +#include <set> + +#include "random.hpp" +#include "spin.hpp" + +template <class U, int D, class R, class S> class Model; + +template <class U, int D, class R, class S> class Transformation { +public: + virtual std::set<Spin<U, D, S>*> current() const { return {NULL}; } + virtual std::set<Spin<U, D, S>*> toConsider() const { return {}; } + virtual double ΔE(const Spin<U, D, S>*) const { return 0; } + virtual Transformation* createNew(Spin<U, D, S>*) const { return new Transformation(); } + virtual void apply(){}; +}; + +template <class U, int D, class R, class S> +using Gen = std::function<Transformation<U, D, R, S>*(Model<U, D, R, S>&, Rng&)>; + +template <class U, int D, class R, class S> class SpinFlip : public Transformation<U, D, R, S> { +private: + Model<U, D, R, S>& M; + Spin<U, D, S>* sOld; + const R r; + Spin<U, D, S> sNew; + +public: + SpinFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s) : M(M), r(r) { + sOld = s; + sNew = r.act(*s); + } + + std::set<Spin<U, D, S>*> current() const override { return {sOld}; } + + std::set<Spin<U, D, S>*> toConsider() const override { + std::set<Spin<U, D, S>*> neighbors; + std::set<Spin<U, D, S>*> current_neighbors = M.dict.neighbors(sOld->x); + std::set<Spin<U, D, S>*> new_neighbors = M.dict.neighbors(sNew.x); + + neighbors.insert(current_neighbors.begin(), current_neighbors.end()); + neighbors.insert(new_neighbors.begin(), new_neighbors.end()); + neighbors.insert(NULL); + + return neighbors; + } + + double ΔE(const Spin<U, D, S>* s) const override { + if (s == NULL) { + Spin<U, D, S> s0s_old = M.s0.inverse().act(*sOld); + Spin<U, D, S> s0s_new = M.s0.inverse().act(sNew); + return M.B(s0s_new) - M.B(s0s_old); + } else { + return M.Z(*sOld, *s) - M.Z(sNew, *s); + } + } + + Transformation<U, D, R, S>* createNew(Spin<U, D, S>* s) const override; + + void apply() override { + M.dict.remove(sOld); + *sOld = sNew; + M.dict.insert(sOld); + } +}; + +template <class U, int D, class R, class S> class PairFlip : public Transformation<U, D, R, S> { +private: + Model<U, D, R, S>& M; + Spin<U, D, S>* s1Old; + Spin<U, D, S>* s2Old; + const R r; + Spin<U, D, S> s1New; + Spin<U, D, S> s2New; + +public: + PairFlip(Model<U, D, R, S>& M, const R& r, Spin<U, D, S>* s1, Spin<U, D, S>* s2) : M(M), r(r) { + s1Old = s1; + s2Old = s2; + s1New = r.act(*s1); + s2New = r.act(*s2); + } + + std::set<Spin<U, D, S>*> current() const override { return {s1Old, s2Old}; } + + std::set<Spin<U, D, S>*> toConsider() const override { + std::set<Spin<U, D, S>*> neighbors; + std::set<Spin<U, D, S>*> current_neighbors_1 = M.dict.neighbors(s1Old->x); + std::set<Spin<U, D, S>*> current_neighbors_2 = M.dict.neighbors(s2Old->x); + std::set<Spin<U, D, S>*> new_neighbors_1 = M.dict.neighbors(s1New.x); + std::set<Spin<U, D, S>*> new_neighbors_2 = M.dict.neighbors(s2New.x); + + neighbors.insert(current_neighbors_1.begin(), current_neighbors_1.end()); + neighbors.insert(current_neighbors_2.begin(), current_neighbors_2.end()); + neighbors.insert(new_neighbors_1.begin(), new_neighbors_1.end()); + neighbors.insert(new_neighbors_2.begin(), new_neighbors_2.end()); + neighbors.insert(NULL); + + return neighbors; + } + + double ΔE(const Spin<U, D, S>* s) const override { + if (s == NULL) { + Spin<U, D, S> s0s1_old = M.s0.inverse().act(*s1Old); + Spin<U, D, S> s0s1_new = M.s0.inverse().act(s1New); + Spin<U, D, S> s0s2_old = M.s0.inverse().act(*s2Old); + Spin<U, D, S> s0s2_new = M.s0.inverse().act(s2New); + return M.B(s0s1_new) + M.B(s0s2_new) - M.B(s0s1_old) - M.B(s0s2_old); + } else { + return M.Z(*s1Old, *s) + M.Z(*s2Old, *s) - M.Z(s1New, *s) - M.Z(s2New, *s); + } + } + + Transformation<U, D, R, S>* createNew(Spin<double, D, S>* s) const; + Transformation<U, D, R, S>* createNew(Spin<signed, D, S>* s) const; + + void apply() override { + M.dict.remove(s1Old); + M.dict.remove(s2Old); + *s1Old = s1New; + *s2Old = s2New; + M.dict.insert(s1Old); + M.dict.insert(s2Old); + } +}; + +template <class U, int D, class R, class S> class FieldFlip : public Transformation<U, D, R, S> { +private: + Model<U, D, R, S>& M; + const R r; + R s0New; + +public: + FieldFlip(Model<U, D, R, S>& M, const R& r) : M(M), r(r), s0New(r.act(M.s0)) {} + + std::set<Spin<U, D, S>*> toConsider() const override { + std::set<Spin<U, D, S>*> neighbors; + + for (Spin<U, D, S>* s : M.s) { + neighbors.insert(s); + } + + return neighbors; + } + + double ΔE(const Spin<U, D, S>* s) const override { + Spin<U, D, S> s0s_old = M.s0.inverse().act(*s); + Spin<U, D, S> s0s_new = s0New.inverse().act(*s); + return M.B(s0s_new) - M.B(s0s_old); + } + + Transformation<double, D, R, S>* createNew(Spin<double, D, S>* s) const { + return new SpinFlip<U, D, R, S>(M, r, s); + } + + Transformation<signed, D, R, S>* createNew(Spin<signed, D, S>* s) const { + Vector<signed, 2> v = r.act(*s).x; + std::set<Spin<U, D, S>*> on_site = M.dict.at(v); + if (on_site.empty()) { + return new SpinFlip<U, D, R, S>(M, r, s); + } else { + return new PairFlip<U, D, R, S>(M, r, s, *on_site.begin()); + } + } + + void apply() override { M.s0 = s0New; } +}; + +template <class U, int D, class R, class S> +Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(Spin<U, D, S>* s) const { + if (s == NULL) { + return new FieldFlip<U, D, R, S>(M, r); + } else { + return new SpinFlip(M, r, s); + } +} + +template <class U, int D, class R, class S> +Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<double, D, S>* s) const { + if (s == NULL) { + return new FieldFlip<U, D, R, S>(M, r); + } else { + return new SpinFlip<U, D, R, S>(M, r, s); + } +} + +template <class U, int D, class R, class S> +Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(Spin<signed, D, S>* s) const { + if (s == NULL) { + return new FieldFlip<U, D, R, S>(M, r); + } else { + Vector<signed, 2> v = r.act(*s).x; + std::set<Spin<U, D, S>*> on_site = M.dict.at(v); + if (on_site.empty()) { + return new SpinFlip<U, D, R, S>(M, r, s); + } else { + return new PairFlip<U, D, R, S>(M, r, s, *on_site.begin()); + } + } +} |