diff options
Diffstat (limited to 'glass.hpp')
| -rw-r--r-- | glass.hpp | 209 |
1 files changed, 205 insertions, 4 deletions
@@ -180,6 +180,7 @@ public: const unsigned L; unsigned N; std::vector<Vertex<D, S>> vertices; + Transformation<D> orientation; unsigned vectorToIndex(const Vector<D>& x) const { unsigned i = 0; @@ -197,7 +198,7 @@ public: return x; } - System(unsigned L, unsigned N) : L(L), N(N), vertices(iPow(L, D)) { + System(unsigned L) : L(L), N(0), vertices(iPow(L, D)), orientation(L) { for (unsigned i = 0; i < iPow(L, D); i++) { vertices[i].position = indexToVector(i); vertices[i].initialPosition = vertices[i].position; @@ -216,8 +217,14 @@ public: } } + unsigned size() const { return vertices.size(); } + double density() const { return N / pow(L, D); } + unsigned maxOccupation() const { + return iPow(L, D); + } + void setInitialPosition() { for (Vertex<D, S>& v : vertices) { v.initialPosition = v.position; @@ -240,7 +247,7 @@ public: } }; -template <unsigned D, State S> class FiniteEnergySystem : public System<D, S> { +template <unsigned D, State S> class SoftSystem : public System<D, S> { public: using System<D, S>::System; @@ -326,7 +333,7 @@ public: bool dry = false) { std::queue<std::array<std::reference_wrapper<Vertex<D, S>>, 2>> q; Vector<D> x0New = R.apply(v0.position); - Vertex<D, S>& v0New = this->vertices[this->vectorToIndex(x0New)]; + Vertex<D, S>& v0New = SoftSystem::vertices[SoftSystem::vectorToIndex(x0New)]; q.push({v0, v0New}); unsigned n = 0; @@ -346,7 +353,7 @@ public: Vertex<D, S>& vn = e.neighbor; Vector<D> xnNew = R.apply(vn.position); - Vertex<D, S>& vnNew = this->vertices[this->vectorToIndex(xnNew)]; + Vertex<D, S>& vnNew = SoftSystem::vertices[SoftSystem::vectorToIndex(xnNew)]; if (!vn.marked && !vnNew.marked) { double E0 = pairEnergy(v.state, vn.state) + pairEnergy(vNew.state, vnNew.state); @@ -393,3 +400,197 @@ public: } }; +template <unsigned D, State S> class HardSystem : public System<D, S> { +public: + + HardSystem(unsigned L) : System<D, S>(L) {} + + virtual std::list<std::reference_wrapper<Vertex<D, S>>> overlaps(Vertex<D, S>&, const S&, + bool = false) { return {}; } + + bool insert(Vertex<D, S>& v, const S& s) { + if (overlaps(v, s).empty()) { + v.state = s; + HardSystem::N++; + return true; + } else { + return false; + } + } + + bool tryDeletion(Vertex<D, S>& v) { + if (v.empty()) { + return false; + } else { + v.state.remove(); + HardSystem::N--; + return true; + } + } + + bool tryRandomMove(Rng& r) { + Vertex<D, S>& v = r.pick(HardSystem::vertices); + S oldState = v.state; + + if (!tryDeletion(v)) { + return false; + } + + if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) { + for (HalfEdge<D, S>& e : v.adjacentEdges) { + if (1 == e.Δx.dot(oldState)) { + if (insert(e.neighbor, oldState.flip())) { + return true; + } + break; + } + } + } else { + S newState(r); + while (newState == oldState) { + newState = S(r); + } + if (insert(v, newState)) { + return true; + } + } + v.state = oldState; + HardSystem::N++; + return false; + } + + bool trySwap(Vertex<D, S>& v1, Vertex<D, S>& v2) { + if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) { + std::swap(v1.state, v2.state); + return true; + } else { + return false; + } + } + + bool tryRandomSwap(Rng& r) { + Vertex<D, S>& v1 = r.pick(HardSystem::vertices); + Vertex<D, S>& v2 = r.pick(HardSystem::vertices); + + return trySwap(v1, v2); + } + + + bool compatible() { + for (Vertex<D, S>& v : HardSystem::vertices) { + if (overlaps(v, v.state, true).size() > 0) { + return false; + } + } + + return true; + } + + void sweepGrandCanonical(double z, Rng& r) { + for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) { + if (0.5 < r.uniform(0.0, 1.0)) { + double pIns = HardSystem::maxOccupation() * z / (HardSystem::N + 1); + + if (pIns > r.uniform(0.0, 1.0)) { + while (true) { + Vertex<D, S>& v = r.pick(HardSystem::vertices); + if (v.empty()) { + insert(v, S(r)); + break; + } + } + } + } else { + + double pDel = HardSystem::N / (z * HardSystem::maxOccupation()); + + if (pDel > r.uniform(0.0, 1.0)) { + tryDeletion(r.pick(HardSystem::vertices)); + } + } + + tryRandomMove(r); + } + } + + void sweepLocal(Rng& r) { + for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) { + tryRandomMove(r); + } + } + + void sweepSwap(Rng& r) { + for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) { + tryRandomSwap(r); + } + } + + unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, Rng& r, bool dry = false) { + std::queue<std::reference_wrapper<Vertex<D, S>>> q; + q.push(v0); + + unsigned n = 0; + + while (!q.empty()) { + Vertex<D, S>& v = q.front(); + q.pop(); + + if (!v.marked) { + Vector<D> xNew = R.apply(v.position); + Vertex<D, S>& vNew = HardSystem::vertices[HardSystem::vectorToIndex(xNew)]; + + v.marked = true; + vNew.marked = true; + + S s = R.apply(v.state); + S sNew = R.apply(vNew.state); + + std::list<std::reference_wrapper<Vertex<D, S>>> overlaps1 = overlaps(vNew, s, true); + std::list<std::reference_wrapper<Vertex<D, S>>> overlaps2 = overlaps(v, sNew, true); + overlaps1.splice(overlaps1.begin(), overlaps2); + + for (Vertex<D, S>& vn : overlaps1) { + if (!vn.marked) { + q.push(vn); + } + } + + if (!dry) { + v.state = sNew; + vNew.state = s; + } + + n += 1; + } + } + + return n; + } + + void swendsenWang(const Transformation<D>& R, Rng& r) { + for (Vertex<D, S>& v : HardSystem::vertices) { + if (!v.marked) { + bool dry = 0.5 < r.uniform(0.0, 1.0); + unsigned n = flipCluster(R, v, r, dry); + if (n > pow(HardSystem::L, D) / 4 && !dry) { + orientation = R.apply(orientation); + } + } + } + + for (Vertex<D, S>& v : HardSystem::vertices) { + v.marked = false; + } + } + + int overlap(const System<D, S>& s) const { + int o = 0; + + for (unsigned i = 0; i < HardSystem::vertices.size(); i++) { + S s2 = orientation.apply(s.vertices[HardSystem::vectorToIndex(orientation.inverse().apply(HardSystem::indexToVector(i)))].state); + o += HardSystem::vertices[i].state.dot(s2); + } + + return o; + } +}; |