summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--space_wolff.hpp8
-rw-r--r--spheres_infinite.cpp27
-rw-r--r--transformation.hpp46
3 files changed, 49 insertions, 32 deletions
diff --git a/space_wolff.hpp b/space_wolff.hpp
index 141396a..63bf118 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -50,19 +50,23 @@ public:
queue.push(t0);
std::set<Spin<U, D, S>*> visited;
+ std::set<Spin<U, D, S>*> fused;
while (!queue.empty()) {
Transformation<U, D, R, S>* t = queue.front();
queue.pop();
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); })) {
+ if (!std::any_of(c.begin(), c.end(), [&visited](Spin<U, D, S>* s) { return visited.contains(s); }) && !(c.size() == 1 && std::any_of(c.begin(), c.end(), [&fused](Spin<U, D, S>* s) { return fused.contains(s); }))) {
visited.insert(c.begin(), c.end());
+ fused.insert(c.begin(), c.end());
for (Spin<U, D, S>* s : t->toConsider()) {
+ if (!fused.contains(s)) {
double ΔE = t->ΔE(s);
if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform<double>(0, 1)) {
- queue.push(t->createNew(T, visited, s, rng));
+ queue.push(t->createNew(T, fused, s, rng));
+ }
}
}
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index 5287a84..b39f164 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -117,6 +117,28 @@ Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
};
}
+Gen<double, D, Euclidean<double, D>, double> rGen(double ε) {
+ return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
+ Vector<double, D> t;
+ Matrix<double, D> m;
+
+ double θ = rng.uniform((double)0.0, 2 * M_PI);
+ m(0, 0) = -cos(2 * θ);
+ m(1, 1) = cos(2 * θ);
+ m(0, 1) = -2 * cos(θ) * sin(θ);
+ m(1, 0) = -2 * cos(θ) * sin(θ);
+
+ Spin<double, D, double>* s = rng.pick(M.s);
+ t = M.s0.t;
+ 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, s);
+ };
+}
+
int main(int argc, char* argv[]) {
const unsigned D = 2;
@@ -179,8 +201,9 @@ int main(int argc, char* argv[]) {
return H * s.x.norm();
};
- auto g1 = eGen(1);
+ auto g1 = eGen(2);
auto g2 = mGen(0.5);
+ auto g3 = rGen(5);
animation A(L, 750, argc, argv);
model sphere(1.0, Z, B);
@@ -197,7 +220,7 @@ int main(int argc, char* argv[]) {
sphere.dict.insert(ss);
}
- sphere.wolff(T, {g1, g2}, A, N);
+ sphere.wolff(T, {g1, g2, g3}, A, N);
std::ofstream outfile;
outfile.open("test.dat");
diff --git a/transformation.hpp b/transformation.hpp
index fc1833e..ad3b49f 100644
--- a/transformation.hpp
+++ b/transformation.hpp
@@ -12,7 +12,7 @@ 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(double T, const std::set<Spin<U, D, S>*>&, Spin<U, D, S>*, Rng& rng) const { return new Transformation(); }
+ virtual Transformation* createNew(double T, std::set<Spin<U, D, S>*>&, Spin<U, D, S>*, Rng& rng) const { return new Transformation(); }
virtual void apply(){};
};
@@ -43,6 +43,8 @@ public:
neighbors.insert(new_neighbors.begin(), new_neighbors.end());
neighbors.insert(NULL);
+ neighbors.erase(sOld);
+
return neighbors;
}
@@ -56,7 +58,7 @@ public:
}
}
- Transformation<U, D, R, S>* createNew(double, const std::set<Spin<U, D, S>*>&, Spin<U, D, S>* s, Rng&) const override;
+ Transformation<U, D, R, S>* createNew(double, std::set<Spin<U, D, S>*>&, Spin<U, D, S>* s, Rng&) const override;
void apply() override {
M.dict.remove(sOld);
@@ -97,6 +99,9 @@ public:
neighbors.insert(new_neighbors_2.begin(), new_neighbors_2.end());
neighbors.insert(NULL);
+ neighbors.erase(s1Old);
+ neighbors.erase(s2Old);
+
return neighbors;
}
@@ -112,8 +117,8 @@ public:
}
}
- Transformation<U, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const;
- Transformation<U, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<signed, D, S>* s, Rng&) const;
+ Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const;
+ Transformation<U, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<signed, D, S>* s, Rng&) const;
void apply() override {
M.dict.remove(s1Old);
@@ -150,11 +155,11 @@ public:
return M.B(s0s_new) - M.B(s0s_old);
}
- Transformation<double, D, R, S>* createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const {
+ Transformation<double, D, R, S>* createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng&) const {
return new SpinFlip<U, D, R, S>(M, r, s);
}
- Transformation<signed, D, R, S>* createNew(double, const std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
+ Transformation<signed, D, R, S>* createNew(double, std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) 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()) {
@@ -168,33 +173,16 @@ public:
};
template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(double T, const std::set<Spin<U, D, S>*>& v, Spin<U, D, S>* s, Rng& rng) const {
+Transformation<U, D, R, S>* SpinFlip<U, D, R, S>::createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<U, D, S>* s, Rng& rng) const {
if (s == NULL) {
return new FieldFlip<U, D, R, S>(M, r);
} else {
- SpinFlip<U, D, R, S>* t = new SpinFlip<U, D, R, S>(M, r, s);
- Spin<U, D, S>* sMax = NULL;
- double ΔEMax = 0.0;
- for (Spin<U, D, S>* ss : t->toConsider()) {
- if (ss != NULL && ss != s && !v.contains(ss)) {
- double ΔE = t->ΔE(ss);
- if (ΔE > ΔEMax) {
- sMax = ss;
- ΔEMax = ΔE;
- }
- }
- }
- if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1)) {
- return t;
- } else {
- delete t;
- return new PairFlip<U, D, R, S>(M, r, s, sMax);
- }
+ 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(double T, const std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng& rng) const {
+Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, std::set<Spin<U, D, S>*>& v, Spin<double, D, S>* s, Rng& rng) const {
if (s == NULL) {
return new FieldFlip<U, D, R, S>(M, r);
} else {
@@ -210,17 +198,19 @@ Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double T, const std:
}
}
}
- if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1)) {
+ if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform<double>(0, 1) || v.contains(sMax)) {
return t;
} else {
delete t;
+ v.insert(s);
+ v.insert(sMax);
return new PairFlip<U, D, R, S>(M, r, s, sMax);
}
}
}
template <class U, int D, class R, class S>
-Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double, const std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
+Transformation<U, D, R, S>* PairFlip<U, D, R, S>::createNew(double, std::set<Spin<U, D, S>*>&, Spin<signed, D, S>* s, Rng&) const {
if (s == NULL) {
return new FieldFlip<U, D, R, S>(M, r);
} else {