diff options
| -rw-r--r-- | biroli-mezard.cpp | 84 | ||||
| -rw-r--r-- | glass.hpp | 241 |
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; +} + @@ -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; + } +}; |