1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
#include "smiley.hpp"
#include "space_wolff.hpp"
#include "torus_symmetries.hpp"
const unsigned D = 2;
typedef Model<signed, D, TorusGroup<signed, D>, signed> isingModel;
typedef Spin<signed, D, signed> isingSpin;
std::function<double(const isingSpin&, const isingSpin&)> isingZ(unsigned L) {
return [L](const isingSpin& s1, const isingSpin& 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(isingSpin)> isingBFace(unsigned L, double H) {
return [L, H](const isingSpin& s) -> double {
return H * s.s * smiley[15 - s.x(1) * 16 / L][s.x(0) * 16 / L];
};
}
std::function<double(isingSpin)> isingBMod(unsigned L, unsigned mod, double H) {
return [mod, L, H](const isingSpin& s) -> double {
return H * s.s * sin(s.x(0) * mod * 2 * M_PI / L);
};
}
Gen<signed, D, TorusGroup<signed, D>, signed> isingGen(unsigned L) {
std::vector<Vector<signed, 2>> torusVectors = torus_vecs<signed, 2>(L);
std::vector<Matrix<signed, 2>> torusMatrices = torus_mats<signed, 2>();
return [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);
if (s2s.empty()) {
return new SpinFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss);
} else {
return new PairFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss, *s2s.begin());
}
};
}
void isingPopulate(isingModel& m, unsigned L, double pop, double mag) {
Rng rng;
for (unsigned i = 0; i < L; i++) {
for (unsigned j = 0; j < L; j++) {
if (rng.uniform<double>(0, 1) < pop) {
Spin<signed, D, signed>* ss = new Spin<signed, D, signed>();
ss->x = {i, j};
if (rng.uniform<double>(0, 1) < mag) {
ss->s = 1;
} else {
ss->s = -1;
}
m.s.push_back(ss);
m.dict.insert(ss);
}
}
}
}
|
