summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--.gitignore7
-rw-r--r--ising.cpp250
-rw-r--r--octree.hpp19
-rw-r--r--random.hpp6
-rw-r--r--space_wolff.hpp172
-rw-r--r--spheres_infinite.cpp46
-rw-r--r--spin.hpp2
-rw-r--r--torus_symmetries.hpp2
-rw-r--r--transformation.hpp201
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
+*/
+
diff --git a/ising.cpp b/ising.cpp
index 1987a8f..38dce3b 100644
--- a/ising.cpp
+++ b/ising.cpp
@@ -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;
}
diff --git a/octree.hpp b/octree.hpp
index 1be6e0a..bb13242 100644
--- a/octree.hpp
+++ b/octree.hpp
@@ -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();
diff --git a/spin.hpp b/spin.hpp
index d4d39cc..bb8369b 100644
--- a/spin.hpp
+++ b/spin.hpp
@@ -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());
+ }
+ }
+}