summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--biroli-mezard.cpp84
-rw-r--r--glass.hpp241
2 files changed, 289 insertions, 36 deletions
diff --git a/biroli-mezard.cpp b/biroli-mezard.cpp
new file mode 100644
index 0000000..ac68776
--- /dev/null
+++ b/biroli-mezard.cpp
@@ -0,0 +1,84 @@
+#include <iostream>
+#include <fstream>
+
+#include "glass.hpp"
+
+void print(const BiroliSystem<2>& s) {
+ for (const Vertex<2, BiroliState>& v : s.vertices) {
+ std::cerr << v.state.type;
+
+ if (v.position(0) == s.L - 1) {
+ std::cerr << std::endl;
+ }
+ }
+}
+
+int main() {
+ const unsigned D = 3;
+ unsigned L = 30;
+ unsigned Nmin = 2e2;
+ unsigned Nmax = 2e5;
+ double Tmin = 0.04;
+ double Tmax = 0.2;
+ double δT = 0.02;
+
+ BiroliSystem<D> s(L);
+
+ Rng r;
+
+ double z = exp(1 / 0.2);
+
+ if (!s.compatible()) {
+ std::cerr << "Storted incompatible!" << std::endl;
+ return 1;
+ }
+
+ while (s.density() < 0.57) {
+ s.sweepGrandCanonical(z, r);
+ }
+
+ if (!s.compatible()) {
+ std::cerr << "Not compatible!" << std::endl;
+ return 1;
+ }
+
+ std::cerr << "Found state with appropriate density." << std::endl;
+
+ BiroliSystem<D> s0 = s;
+
+ std::vector<Matrix<D>> ms = generateTorusMatrices<D>();
+
+ std::vector<unsigned> clusterDist(s.size() + 1);
+
+ 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;
+ }
+ 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.swendsenWang(Transformation<D>(L, ms, r), r);
+ i++;
+ }
+
+ if (!s.compatible()) {
+ std::cerr << "Not compatible!" << std::endl;
+ return 1;
+ }
+
+ std::ofstream file("dist.dat");
+ for (unsigned i : clusterDist) {
+ file << i << " ";
+ }
+ file.close();
+
+ return 0;
+}
+
diff --git a/glass.hpp b/glass.hpp
index e67bb51..b35964b 100644
--- a/glass.hpp
+++ b/glass.hpp
@@ -1,3 +1,4 @@
+#include <limits>
#include <list>
#include <queue>
#include <vector>
@@ -130,7 +131,7 @@ public:
unsigned occupiedNeighbors;
BiroliState() {
- type = 0;
+ type = std::numeric_limits<unsigned>::max();
occupiedNeighbors = 0;
}
@@ -139,9 +140,13 @@ public:
occupiedNeighbors = o;
}
- bool isEmpty() const { return type == 0; }
+ BiroliState(Rng& r) {
+ type = r.pick({1,2,2,2,2,2,3,3,3,3});
+ }
+
+ bool isEmpty() const { return type == std::numeric_limits<unsigned>::max(); }
- void remove() { type = 0; };
+ void remove() { type = std::numeric_limits<unsigned>::max(); };
};
template <int D> class Transformation {
@@ -246,10 +251,9 @@ public:
}
}
- virtual std::list<std::reference_wrapper<Vertex<D, State>>> overlaps(Vertex<D, State>&,
- const State&, bool = false) { return {}; };
+ virtual std::list<std::reference_wrapper<Vertex<D, State>>> overlaps(Vertex<D, State>&) { return {}; };
- virtual bool insert(Vertex<D, State>& v, const State& s) { return false; };
+ virtual bool insert(Vertex<D, State>& v, const State& s, bool force = false) { return false; };
virtual bool remove(Vertex<D, State>& v) { return false;};
@@ -268,7 +272,7 @@ public:
bool compatible() {
for (Vertex<D, State>& v : vertices) {
- if (overlaps(v, v.state, true).size() > 0) {
+ if (overlaps(v).size() > 0) {
return false;
}
}
@@ -319,43 +323,78 @@ public:
}
}
- unsigned flipCluster(const Transformation<D>& R, Vertex<D, State>& v0, Rng& r, bool dry = false) {
- std::queue<std::reference_wrapper<Vertex<D, State>>> q;
- q.push(v0);
-
+ unsigned flipCluster(const Transformation<D>& R, Vertex<D, State>& v0, bool dry = false) {
unsigned n = 0;
- while (!q.empty()) {
- Vertex<D, State>& v = q.front();
- q.pop();
+ Vertex<D, State>& v0New = vertices[vectorToIndex(R.apply(v0.position))];
- if (!v.marked) {
- Vertex<D, State>& vNew = vertices[vectorToIndex(R.apply(v.position))];
+ if (&v0New != &v0) {
+ std::queue<std::reference_wrapper<Vertex<D, State>>> q;
+
+ v0.marked = true;
+ v0New.marked = true;
+
+ swap(v0, v0New);
+
+ for (Vertex<D, State>& vn : overlaps(v0)) {
+ if (!vn.marked) {
+ q.push(vn);
+ }
+ }
+ for (Vertex<D, State>& vn : overlaps(v0New)) {
+ if (!vn.marked) {
+ q.push(vn);
+ }
+ }
- v.marked = true;
- vNew.marked = true;
+ while (!q.empty()) {
+ Vertex<D, State>& v = q.front();
+ q.pop();
- State s = R.apply(v.state);
- State sNew = R.apply(vNew.state);
+ if (!v.marked && !overlaps(v).empty()) {
+ v.marked = true;
+ Vertex<D, State>& vNew = vertices[vectorToIndex(R.apply(v.position))];
- std::list<std::reference_wrapper<Vertex<D, State>>> overlaps1 = overlaps(vNew, s, true);
- std::list<std::reference_wrapper<Vertex<D, State>>> overlaps2 = overlaps(v, sNew, true);
- overlaps1.splice(overlaps1.begin(), overlaps2);
+ if (&vNew != &v) {
+ vNew.marked = true;
+
+ swap(v, vNew);
+
+ for (Vertex<D, State>& vn : overlaps(v)) {
+ if (!vn.marked) {
+ q.push(vn);
+ }
+ }
+ for (Vertex<D, State>& vn : overlaps(vNew)) {
+ if (!vn.marked) {
+ q.push(vn);
+ }
+ }
- for (Vertex<D, State>& vn : overlaps1) {
- if (!vn.marked) {
- q.push(vn);
+ n += 2;
+ } else {
+ n += 1;
}
}
-
- if (!dry) {
- swap(v, vNew);
- v.state = sNew;
- vNew.state = s;
+ if (q.empty()) {
+ for (Vertex<D, State>& vv : vertices) {
+ if (!vv.marked && !overlaps(vv).empty()) {
+// std::cerr << "Found bad state at end" << std::endl;
+ q.push(vv);
+ }
+ }
}
-
- n += 1;
}
+
+
+ }
+
+ if (n > size() / 4) {
+ orientation = R.apply(orientation);
+ }
+
+ for (Vertex<D, State>& v : vertices) {
+ v.marked = false;
}
return n;
@@ -365,7 +404,7 @@ public:
for (Vertex<D, State>& v : vertices) {
if (!v.marked) {
bool dry = 0.5 < r.uniform(0.0, 1.0);
- unsigned n = flipCluster(R, v, r, dry);
+ unsigned n = flipCluster(R, v, dry);
if (n > size() / 4 && !dry) {
orientation = R.apply(orientation);
}
@@ -408,8 +447,8 @@ template <int D> class CiamarraSystem : public System<D, CiamarraState<D>> {
return o;
}
- bool insert(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s) override {
- if (overlaps(v, s).empty()) {
+ bool insert(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s, bool force = false) override {
+ if (force || overlaps(v, s).empty()) {
v.state = s;
this->N++;
return true;
@@ -506,3 +545,133 @@ template <int D> class CiamarraSystem : public System<D, CiamarraState<D>> {
return o;
}
};
+
+template <int D> class BiroliSystem : public System<D, BiroliState> {
+public:
+ using System<D, BiroliState>::System;
+
+ bool canReplaceWith(const Vertex<D, BiroliState>& v, const BiroliState& s) {
+ if (s.isEmpty()) {
+ return true;
+ }
+ if (s.type < v.state.occupiedNeighbors) {
+ return false;
+ } else {
+ int Δn = 0;
+ if (v.isEmpty()) {
+ Δn += 1;
+ }
+ for (const HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
+ if (e.neighbor.state.type < e.neighbor.state.occupiedNeighbors + Δn) {
+ return false;
+ }
+ }
+ }
+ return true;
+ }
+
+ std::list<std::reference_wrapper<Vertex<D, BiroliState>>>
+ overlaps(Vertex<D, BiroliState>& v) override {
+ std::list<std::reference_wrapper<Vertex<D, BiroliState>>> o;
+
+ if (v.isEmpty()) { // an empty site cannot be frustrating anyone
+ return o;
+ }
+
+ bool selfFrustrated = v.state.occupiedNeighbors > v.state.type;
+ bool anyNeighborFrustrated = false;
+
+ for (const HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
+ if (!e.neighbor.isEmpty()) {
+ bool thisNeighborFrustrated = e.neighbor.state.type < e.neighbor.state.occupiedNeighbors;
+
+ if (thisNeighborFrustrated) {
+ anyNeighborFrustrated = true;
+ }
+
+ if (selfFrustrated || thisNeighborFrustrated) {
+ o.push_back(e.neighbor);
+ }
+ }
+ }
+
+ if (selfFrustrated || anyNeighborFrustrated) {
+ o.push_back(v);
+ }
+
+ return o;
+ }
+
+ bool insert(Vertex<D, BiroliState>& v, const BiroliState& s, bool force = false) override {
+ if (force || (v.isEmpty() && canReplaceWith(v, s))) {
+ v.state.type = s.type;
+ for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
+ e.neighbor.state.occupiedNeighbors++;
+ }
+ this->N++;
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ bool remove(Vertex<D, BiroliState>& v) override {
+ if (v.isEmpty()) {
+ return false;
+ } else {
+ v.state.remove();
+ for (HalfEdge<D, BiroliState>& e : v.adjacentEdges) {
+ e.neighbor.state.occupiedNeighbors--;
+ }
+ this->N--;
+ return true;
+ }
+ }
+
+ bool randomMove(Rng& r) override {
+ Vertex<D, BiroliState>& v = r.pick(this->vertices);
+
+ return exchange(v, r.pick(v.adjacentEdges).neighbor);
+ }
+
+ void swap(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) override {
+ if (v1.state.type != v2.state.type) {
+ if (!v1.isEmpty() && !v2.isEmpty()) {
+ std::swap(v1.state.type, v2.state.type);
+ } else if (v1.isEmpty()) {
+ unsigned t = v2.state.type;
+ remove(v2);
+ insert(v1, BiroliState(t, 0), true);
+ } else {
+ unsigned t = v1.state.type;
+ remove(v1);
+ insert(v2, BiroliState(t, 0), true);
+ }
+ }
+ }
+
+ bool exchange(Vertex<D, BiroliState>& v1, Vertex<D, BiroliState>& v2) override {
+ if (canReplaceWith(v1, v2.state) && canReplaceWith(v2, v1.state)) {
+ swap(v1, v2);
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ int overlap(const System<D, BiroliState>& s) const override {
+ int o = 0;
+
+ for (unsigned i = 0; i < this->vertices.size(); i++) {
+ BiroliState s2 = this->orientation.apply(
+ s.vertices[this->vectorToIndex(this->orientation.inverse().apply(this->indexToVector(i)))].state);
+ if (s2.type > 0 && s2.type == this->vertices[i].state.type) {
+ o++;
+ } else {
+ o--;
+ }
+ }
+
+ return o;
+ }
+};