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
107
108
109
|
#include <iostream>
#include "glass.hpp"
class DistinguishableState {
public:
unsigned type;
DistinguishableState() { type = 0; }
DistinguishableState(unsigned t) { type = t; }
bool empty() const { return type == 0; }
bool operator==(const DistinguishableState& s) const { return type == s.type; }
void remove() { type = 0; }
};
template <unsigned D>
class DistinguishableSystem : public FiniteEnergySystem<D, DistinguishableState> {
public:
std::vector<double> interaction;
DistinguishableSystem(unsigned L, unsigned N, Rng& r)
: FiniteEnergySystem<D, DistinguishableState>(L, N), interaction(N * N) {
for (double& V : interaction) {
V = r.uniform(-0.5, 0.5);
}
for (unsigned i = 0; i < this->N; i++) {
DistinguishableSystem::vertices[i].state.type = i + 1;
}
}
double pairEnergy(const DistinguishableState& s1, const DistinguishableState& s2) const override {
if (s1.empty() || s2.empty()) {
return 0;
} else {
if (s1.type < s2.type) {
return interaction[(s1.type - 1) * this->N + (s2.type - 1)];
} else {
return interaction[(s2.type - 1) * this->N + (s1.type - 1)];
}
}
}
};
int main() {
const unsigned D = 2;
unsigned L = 40;
unsigned N = 1560;
double Tmin = 0;
double Tmax = 0.4;
double δT = 0.02;
Rng r;
DistinguishableSystem<D> s(L, N, r);
std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
for (unsigned j = 0; j < 10 * s.vertices.size(); j++) {
// unsigned nC = s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
s.tryRandomSwap(1000, r);
}
for (unsigned i = 0; i < 1e5; i++) {
for (unsigned j = 0; j < s.vertices.size(); j++) {
s.tryRandomMove(Tmax, r);
}
}
unsigned n = 1;
for (double T = Tmax; T > Tmin; T -= δT) {
/*
for (unsigned i = 0; i < 1e5; i++) {
Matrix<D> m = r.pick(ms);
Vertex<D, DistinguishableState>& v = r.pick(s.vertices);
s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
}
*/
/*
*/
s.setInitialPosition();
std::cout << T << " ";
for (unsigned i = 0; i < 1e4; i++) {
Matrix<D> m = r.pick(ms);
Vertex<D, DistinguishableState>& v = r.pick(s.vertices);
s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
std::cout << s.selfIntermediateScattering() << " ";
}
/*
for (unsigned i = 0; i < 1e5; i++) {
for (unsigned j = 0; j < s.vertices.size(); j++) {
s.tryRandomMove(T, r);
}
std::cout << s.selfIntermediateScattering() << " ";
}
*/
std::cout << std::endl;
std::cerr << T << " " << s.energy() / N << std::endl;
// s.sweepLocal(r);
// s.sweepSwap(r);
// s.swendsenWang(Transformation<D>(L, ms, r), r);
}
return 0;
}
|
