summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--biroli-mezard.cpp84
1 files changed, 59 insertions, 25 deletions
diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
index c1ed606..6f390e3 100644
--- a/biroli-mezard.cpp
+++ b/biroli-mezard.cpp
@@ -145,6 +145,7 @@ public:
template <int D> class Vertex {
public:
+ Vector<D> initialPosition;
Vector<D> position;
std::vector<std::reference_wrapper<Vertex<D>>> neighbors;
unsigned occupiedNeighbors;
@@ -170,6 +171,9 @@ public:
return o == occupiedNeighbors;
}
+ void setInitial() {
+ initialPosition = position;
+ }
};
template <int D> class System {
@@ -293,6 +297,12 @@ public:
return true;
}
+ void setInitial() {
+ for (Vertex<D>& v : vertices) {
+ v.setInitial();
+ }
+ }
+
bool insert(Vertex<D>& v, unsigned t, bool force = false) {
if (force || (v.empty() && canReplaceWith(v, t))) {
v.maximumNeighbors = t;
@@ -340,6 +350,7 @@ public:
remove(v1);
insert(v2, t, true);
}
+ std::swap(v1.initialPosition, v2.initialPosition);
}
}
@@ -352,20 +363,41 @@ public:
}
}
- int overlap(const System<D>& s) const {
+ std::complex<double> selfIntermediateScattering(unsigned k, unsigned t = 2) const {
+ std::complex<double> F = 0;
+
+ using namespace std::complex_literals;
+
+ for (const Vertex<D>& v : vertices) {
+ if (v.maximumNeighbors == t) {
+ for (unsigned d = 0; d < D; d++) {
+ F += exp(1i * 2.0 * M_PI * (double)k * (double)(v.position(d) - v.initialPosition(d)) / (double)L);
+ }
+ }
+ }
+
+ F /= D * N[t];
+
+ return F;
+ }
+
+ double overlap(const System<D>& s, unsigned t = 2) const {
int o = 0;
- for (unsigned i = 0; i < size(); i++) {
- unsigned t2 =
- s.vertices[vectorToIndex(orientation.inverse().apply(indexToVector(i)))].maximumNeighbors;
- if (t2 == vertices[i].maximumNeighbors) {
+ for (const Vertex<D>& v1 : vertices) {
+ const Vertex<D>& v2 = s.vertices[vectorToIndex(orientation.inverse().apply(v1.position))];
+ if (v1.maximumNeighbors == v2.maximumNeighbors) {
o++;
- } else {
- o--;
}
}
- return o;
+ unsigned nn = 0;
+
+ for (unsigned n : N) {
+ nn += iPow(n, 2);
+ }
+
+ return ((double)o / (double)size() - (double)nn / pow((double)size(), 2)) / (1 - (double)nn / pow((double)size(), 2));
}
bool randomExchange(Rng& r) {
@@ -403,13 +435,13 @@ public:
}
void sweepLocal(Rng& r) {
- for (unsigned i = 0; i < size(); i++) {
+ for (unsigned i = 0; i < size() / 2; i++) {
randomLocalExchange(r);
}
}
void sweepSwap(Rng& r) {
- for (unsigned i = 0; i < size(); i++) {
+ for (unsigned i = 0; i < size() / 2; i++) {
randomExchange(r);
}
}
@@ -420,7 +452,7 @@ public:
Vertex<D>& v0New = vertices[vectorToIndex(R.apply(v0.position))];
if (&v0New != &v0) {
- std::queue<std::reference_wrapper<Vertex<D>>> q;
+ std::list<std::reference_wrapper<Vertex<D>>> q;
v0.marked = true;
v0New.marked = true;
@@ -429,18 +461,18 @@ public:
for (Vertex<D>& vn : overlaps(v0)) {
if (!vn.marked) {
- q.push(vn);
+ q.push_front(vn);
}
}
for (Vertex<D>& vn : overlaps(v0New)) {
if (!vn.marked) {
- q.push(vn);
+ q.push_front(vn);
}
}
while (!q.empty()) {
Vertex<D>& v = q.front();
- q.pop();
+ q.pop_front();
if (!v.marked && !overlaps(v).empty()) {
v.marked = true;
@@ -453,13 +485,13 @@ public:
for (Vertex<D>& vn : overlaps(v)) {
if (!vn.marked) {
- q.push(vn);
+ q.push_front(vn);
}
}
for (Vertex<D>& vn : overlaps(vNew)) {
if (!vn.marked) {
- q.push(vn);
+ q.push_front(vn);
}
}
@@ -469,11 +501,13 @@ public:
}
}
- // For some reason, the process above misses frustrated states. This is an inelegant stopgap.
+ // For some reason, the process above misses frustrated states. This is
+ // an inelegant stopgap.
if (q.empty()) {
for (Vertex<D>& vv : vertices) {
if (!vv.marked && !overlaps(vv).empty()) {
- q.push(vv);
+ std::cerr << "Found bad vertex." << std::endl;
+ q.push_front(vv);
}
}
}
@@ -528,7 +562,7 @@ int main() {
return 1;
}
- while (s.density() < 0.58) {
+ while (s.density() < 0.56) {
s.sweepGrandCanonical(z, r);
}
@@ -545,22 +579,22 @@ int main() {
std::vector<unsigned> clusterDist(s.size() + 1);
+ s.setInitial();
+
unsigned n = 1;
unsigned i = 0;
double nC = 0;
while (nC / s.size() < 1e5) {
if (n < 20 * log(i + 1)) {
n++;
- std::cout << nC / s.size() << " " << (double)s.overlap(s0) / s.size() << std::endl;
+ std::cout << nC / s.size() << " " << s.overlap(s0) << std::endl;
}
unsigned nn = s.flipCluster(Transformation<D>(L, ms, r), r.pick(s.vertices));
nC += nn;
clusterDist[nn]++;
- /*
- s.sweepLocal(r);
- nC += s.size();
- s.sweepSwap(r);
- */
+ //s.sweepLocal(r);
+ //nC += s.size();
+ //s.sweepSwap(r);
i++;
}