summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitmodules3
m---------pcg-cpp0
-rw-r--r--space_wolff.hpp241
-rw-r--r--spheres_infinite.cpp30
4 files changed, 189 insertions, 85 deletions
diff --git a/.gitmodules b/.gitmodules
index 0a97926..15de333 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,3 +1,6 @@
[submodule "randutils"]
path = randutils
url = https://gist.github.com/imneme/540829265469e673d045/
+[submodule "pcg-cpp"]
+ path = pcg-cpp
+ url = https://github.com/imneme/pcg-cpp
diff --git a/pcg-cpp b/pcg-cpp
new file mode 160000
+Subproject 5b5cac8d61339e810c5dbb4692d868a1d7ca1b2
diff --git a/space_wolff.hpp b/space_wolff.hpp
index 4c935db..8ea7dfb 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -7,6 +7,7 @@
#include <queue>
#include <vector>
+#include "pcg-cpp/include/pcg_random.hpp"
#include "randutils/randutils.hpp"
#include "euclidean.hpp"
@@ -14,6 +15,8 @@
#include "quantity.hpp"
#include "spin.hpp"
+using Rng = randutils::random_generator<pcg32>;
+
template <class U, int D, class R, class S> class Model;
template <class U, int D, class R, class S> class measurement {
@@ -31,42 +34,171 @@ public:
virtual void post_cluster(const Model<U, D, R, S>&){};
};
-template <class U, int D, class R, class S> using Gen = std::function<R(const Model<U, D, R, S>&, randutils::mt19937_rng&)>;
-
-template <class U, int D, class R, class S> class Model;
template <class U, int D, class R, class S> class Transformation {
public:
- virtual void ready() {};
- virtual double p(const Transformation&) const { return 0; }
- virtual void apply() {};
- virtual std::set<Transformation*> to_consider() const {};
+ 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>& s;
- const R& r;
- Spin<U, D, S> s_new;
+ 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:
- SpinFlip(Model<U, D, R, S>& M, Spin<U, D, S>& s, const R& r) : M(M), s(s), r(r) {}
+ 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;
- void ready() override {
- s_new = r.act(s);
- };
+ for (Spin<U, D, S>& s : M.s) {
+ neighbors.insert(&s);
+ }
- double p(const SpinFlip& sf, double T) const override {
- return 1.0 - exp(-(M.Z(s, sf.s) - M.Z(s_new, sf.s)) / T);
+ return neighbors;
}
- double p(const SpinFlip& sf, double T) const override {
- return 1.0 - exp(-(M.Z(s, sf.s) - M.Z(s_new, sf.s)) / T);
+ 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;
@@ -75,70 +207,36 @@ public:
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;
- randutils::mt19937_rng rng;
+ Rng rng;
Model(U L, std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z,
std::function<double(const Spin<U, D, S>&)> B)
: L(L), s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {}
- void step(double T, unsigned ind, R& r, measurement<U, D, R, S>& A) {
- std::queue<Spin<U, D, S>*> queue;
- queue.push(&(s[ind]));
+ void step(double T, Transformation<U, D, R, S>* t0, measurement<U, D, R, S>& A) {
+ std::queue<Transformation<U, D, R, S>*> queue;
+ queue.push(t0);
std::set<Spin<U, D, S>*> visited;
while (!queue.empty()) {
- Spin<U, D, S>* si = queue.front();
+ Transformation<U, D, R, S>* t = queue.front();
queue.pop();
- if (!visited.contains(si)) {
- visited.insert(si);
-
- if (si == NULL) {
- R s0_new = r.act(s0);
- for (Spin<U, D, S>& ss : s) {
- Spin<U, D, S> s0s_old = s0.inverse().act(ss);
- Spin<U, D, S> s0s_new = s0_new.inverse().act(ss);
- double p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T);
- A.ghost_bond_visited(*this, s0s_old, s0s_new, p);
- if (rng.uniform(0.0, 1.0) < p) {
- queue.push(&ss);
- }
- }
- A.ghost_site_transformed(*this, s0_new);
- s0 = s0_new;
- } else {
- Spin<U, D, S> si_new = r.act(*si);
- std::set<Spin<U, D, S>*> all_neighbors;
- std::set<Spin<U, D, S>*> current_neighbors = dict.neighbors(si->x);
- std::set<Spin<U, D, S>*> new_neighbors = dict.neighbors(si_new.x);
-
- all_neighbors.insert(current_neighbors.begin(), current_neighbors.end());
- all_neighbors.insert(new_neighbors.begin(), new_neighbors.end());
- all_neighbors.insert(NULL);
- for (Spin<U, D, S>* sj : all_neighbors) {
- if (sj != si) {
- double p;
- if (sj == NULL) {
- Spin<U, D, S> s0s_old = s0.inverse().act(*si);
- Spin<U, D, S> s0s_new = s0.inverse().act(si_new);
- p = 1.0 - exp(-(B(s0s_new) - B(s0s_old)) / T);
- A.ghost_bond_visited(*this, s0s_old, s0s_new, p);
- } else {
- p = 1.0 - exp(-(Z(*si, *sj) - Z(si_new, *sj)) / T);
- A.plain_bond_visited(*this, si, sj, si_new, p);
- }
- if (rng.uniform(0.0, 1.0) < p) {
- queue.push(sj);
- }
- }
+ std::set<Spin<U, D, S>*> c = t->current();
+ if (!std::any_of(c.begin(), c.end(), [&visited](Spin<U, D, S>* s) { return visited.contains(s); })) {
+ visited.insert(c.begin(), c.end());
+
+ for (Spin<U, D, S>* s : t->toConsider()) {
+ double ΔE = t->ΔE(s);
+ if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform<double>(0, 1)) {
+ queue.push(t->createNew(s));
}
- A.plain_site_transformed(*this, si, si_new);
- dict.remove(si);
- *si = si_new;
- dict.insert(si);
}
+
+ t->apply();
}
+ delete t;
}
}
@@ -148,12 +246,9 @@ public:
measurement<U, D, R, S>& A, unsigned N) {
for (unsigned i = 0; i < N; i++) {
Gen<U, D, R, S>& g = rng.pick(gs);
- R r = g(*this, rng);
- unsigned ind = rng.uniform((unsigned)0, (unsigned)(s.size() - 1));
-
- A.pre_cluster(*this, ind, r);
+ Transformation<U, D, R, S> *t = g(*this, rng);
- this->step(T, ind, r, A);
+ this->step(T, t, A);
A.post_cluster(*this);
}
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index a9cd7d0..3505d6d 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -68,7 +68,7 @@ public:
};
Gen<double, D, Euclidean<double, D>, double> eGen(double ε) {
- return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+ return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
Vector<double, D> t;
Matrix<double, D> m;
@@ -85,12 +85,12 @@ Gen<double, D, Euclidean<double, D>, double> eGen(double ε) {
}
Euclidean<double, D> g(t - m * t, m);
- return g;
+ return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, M.s[f_ind]);
};
}
Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
- return [ε](const model& M, randutils::mt19937_rng& rng) -> Euclidean<double, 2> {
+ 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());
@@ -110,7 +110,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 g;
+ return new PairFlip<double, D, Euclidean<double, D>, double>(M, g, M.s[f_ind1], M.s[f_ind2]);
};
}
@@ -123,9 +123,12 @@ int main(int argc, char* argv[]) {
double H = 1.0;
unsigned n = 25;
+ double k = 1e2;
+ double a = 0.1;
+
int opt;
- while ((opt = getopt(argc, argv, "n:N:L:T:H:")) != -1) {
+ while ((opt = getopt(argc, argv, "n:N:L:T:H:a:k:")) != -1) {
switch (opt) {
case 'n':
n = (unsigned)atof(optarg);
@@ -142,14 +145,17 @@ int main(int argc, char* argv[]) {
case 'H':
H = atof(optarg);
break;
+ case 'a':
+ a = atof(optarg);
+ break;
+ case 'k':
+ k = atof(optarg);
+ break;
default:
exit(1);
}
}
- double k = 1000;
- double a = 0.05;
-
std::function<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
[L, a, k](const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
Vector<double, D> d = s1.x - s2.x;
@@ -170,19 +176,19 @@ int main(int argc, char* argv[]) {
return H * s.x.norm();
};
- auto g1 = eGen(0.25);
- auto g2 = mGen(0.1);
+ auto g1 = eGen(0.5);
+ auto g2 = mGen(0.2);
animation A(L, 750, argc, argv);
model sphere(1, Z, B);
- randutils::mt19937_rng rng;
+ Rng rng;
sphere.s.reserve(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.45, 0.45)});
+ sphere.s.push_back({pos, rng.uniform<double>(0.25, 0.45)});
sphere.dict.insert(&sphere.s.back());
}