diff options
| -rw-r--r-- | glass.cpp | 184 |
1 files changed, 94 insertions, 90 deletions
@@ -1,6 +1,7 @@ #include <iostream> #include <vector> #include <list> +#include <set> #include <eigen3/Eigen/Dense> #include "pcg-cpp/include/pcg_random.hpp" @@ -11,6 +12,7 @@ using Rng = randutils::random_generator<pcg32>; template<int D> using Vector = Eigen::Matrix<int, D, 1>; template<int D> using Matrix = Eigen::Matrix<int, D, D>; +// Integer power function int iPow(int x, unsigned p) { if (p == 0) return 1; if (p == 1) return x; @@ -32,37 +34,22 @@ template<unsigned D> class State { σ[a] = b; } - State() : σ(Vector<D>::Zero()) {} - - bool isEmpty() const { - return σ.squaredNorm() == 0; - } - - void remove() { - σ = Vector<D>::Zero(); - } + State(Rng& r) : State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})) {} State<D> apply(const Matrix<D>& m) const { return {m * σ}; } }; -template<unsigned D> State<D> randomState(Rng& r) { - return State<D>(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})); -} - -template<unsigned D> class HalfEdge; +template <unsigned D> class Vertex; -template<unsigned D> class Vertex { +template<unsigned D> class Particle { public: - Vector<D> position; + Vertex<D>& position; State<D> state; - std::vector<HalfEdge<D>> adjacentEdges; bool marked; - bool isEmpty() const { - return state.isEmpty(); - } + Particle(Vertex<D>& p, const State<D>& s) : position(p), state(s) {} }; template<unsigned D> class HalfEdge { @@ -73,11 +60,18 @@ template<unsigned D> class HalfEdge { HalfEdge(Vertex<D>& n, const Vector<D>& d) : neighbor(n), Δx(d) {} }; +template<unsigned D> class Vertex { + public: + Vector<D> position; + std::vector<HalfEdge<D>> adjacentEdges; + typename std::list<Particle<D>>::iterator occupant; +}; + template<unsigned D> class System { public: const unsigned L; - unsigned N; std::vector<Vertex<D>> vertices; + std::list<Particle<D>> particles; unsigned vectorToIndex(const Vector<D>& x) const { unsigned i = 0; @@ -95,10 +89,19 @@ template<unsigned D> class System { return x; } - System(unsigned L) : L(L), N(0), vertices(iPow(L, D)) { + bool isEmpty(const Vertex<D>& v) const { + return v.occupant == particles.end(); + } + + unsigned N() const { + return particles.size(); + } + + System(unsigned L) : L(L), vertices(iPow(L, D)) { for (unsigned i = 0; i < iPow(L, D); i++) { vertices[i].position = indexToVector(i); vertices[i].adjacentEdges.reserve(2 * D); + vertices[i].occupant = particles.end(); } for (unsigned d = 0; d < D; d++) { @@ -112,21 +115,17 @@ template<unsigned D> class System { } } - std::list<std::reference_wrapper<Vertex<D>>> overlaps(Vertex<D>& v, const State<D>& s, bool excludeSelf = false) { - std::list<std::reference_wrapper<Vertex<D>>> o; + std::list<typename std::list<Particle<D>>::iterator> overlaps(const Particle<D>& p, bool excludeSelf = false) const { + std::list<typename std::list<Particle<D>>::iterator> o; - if (s.isEmpty()) { - return o; + if (!excludeSelf && !isEmpty(p.position)) { + o.push_back(p.position.occupant); } - if (!v.isEmpty() && !excludeSelf) { - o.push_back(v); - } - - for (const HalfEdge<D>& e : v.adjacentEdges) { - if (!e.neighbor.isEmpty()) { - if (s.σ.dot(e.Δx) == 1 || e.neighbor.state.σ.dot(e.Δx) == -1) { - o.push_back(e.neighbor); + for (const HalfEdge<D>& e : p.position.adjacentEdges) { + if (!isEmpty(e.neighbor)) { + if (p.state.σ.dot(e.Δx) == 1 || e.neighbor.occupant->state.σ.dot(e.Δx) == -1) { + o.push_back(e.neighbor.occupant); } } } @@ -134,65 +133,59 @@ template<unsigned D> class System { return o; } - bool tryInsertion(Vertex<D>& v, const State<D>& s) { - if (overlaps(v, s).empty()) { - v.state = s; - N++; + bool insert(const Particle<D>& p) { + if (overlaps(p).empty()) { + particles.push_back(p); + p.position.occupant = std::prev(particles.end()); return true; } else { return false; } } - bool tryDeletion(Vertex<D>& v) { - if (v.isEmpty()) { - return false; - } else { - v.state.remove(); - N--; - return true; - } + void remove(typename std::list<Particle<D>>::iterator ip) { + ip->position.occupant = particles.end(); + particles.erase(ip); } - bool tryRandomMove(Rng& r) { - Vertex<D>& v = r.pick(vertices); - State<D> oldState = v.state; - if (!tryDeletion(v)) { - return false; + /* + bool randomMove(Rng& r) { + Particle<D>& p = *(r.pick(particles)); + Particle<D> pN = p; + + if (0.2 < r.uniform(0.0, 1.0)) { + pN.position = r.pick(p.position.adjacentEdges).neighbor; } + pN.state = randomState<D>(r); - if (0.2 > r.uniform(0.0, 1.0)) { - if (tryInsertion(v, randomState<D>(r))) { - return true; - } + if (overlaps(pN, true).empty()) { + remove(&p); + insert(pN); + return true; } else { - if (tryInsertion(r.pick(v.adjacentEdges).neighbor, randomState<D>(r))) { - return true; - } + return false; } - v.state = oldState; - N++; - return false; } + */ - bool trySwap(Vertex<D>& v1, Vertex<D>& v2) { - if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) { - std::swap(v1.state, v2.state); + bool swap(Particle<D>& p1, Particle<D>& p2) { + std::swap(p1.state, p2.state); + if (overlaps(p1, true).empty() && overlaps(p2, true).empty()) { return true; } else { + std::swap(p1.state, p2.state); return false; } } - bool tryRandomSwap(Rng& r) { - Vertex<D>& v1 = r.pick(vertices); - Vertex<D>& v2 = r.pick(vertices); - - return trySwap(v1, v2); + bool randomSwap(Rng& r) { + return swap(r.pick(particles), r.pick(particles)); } void setGroundState() { - N = 0; + for (typename std::list<Particle<D>>::iterator ip = particles.begin(); ip != particles.end(); ip++) { + remove(ip); + } for (Vertex<D>& v : vertices) { unsigned a = 0; @@ -201,21 +194,21 @@ template<unsigned D> class System { } a %= 2 * D + 1; - v.state.σ = Vector<D>::Zero(); - - if (0 < a && a <= D) { - v.state.σ(a - 1) = -1; - N++; - } else if (D < a) { - v.state.σ(2 * D - a) = 1; - N++; + if (a > 0) { + if (a <= D) { + State<D> s(a - 1, -1); + insert(Particle<D>(v, s)); + } else if (D < a) { + State<D> s(2 * D - a, 1); + insert(Particle<D>(v, s)); + } } } } bool compatible() { - for (Vertex<D>& v : vertices) { - if (overlaps(v, v.state, true).size() > 0) { + for (Particle<D> p : particles) { + if (!overlaps(p, true).empty()) { return false; } } @@ -224,7 +217,7 @@ template<unsigned D> class System { } double density() const { - return N / pow(L, D); + return N() / pow(L, D); } unsigned maxOccupation() const { @@ -235,40 +228,40 @@ template<unsigned D> class System { void sweepGrandCanonical(double z, Rng& r) { for (unsigned i = 0; i < iPow(L, D); i++) { if (0.5 < r.uniform(0.0, 1.0)) { - double pIns = maxOccupation() * z / (N + 1); + double pIns = maxOccupation() * z / (N() + 1); if (pIns > r.uniform(0.0, 1.0)) { while (true) { Vertex<D>& v = r.pick(vertices); - if (v.isEmpty()) { - tryInsertion(v, randomState<D>(r)); + if (isEmpty(v)) { + insert(Particle<D>(v, State<D>(r))); break; } } } } else { - double pDel = N / (z * maxOccupation()); + double pDel = N() / (z * maxOccupation()); if (pDel > r.uniform(0.0, 1.0)) { - tryDeletion(r.pick(vertices)); + remove(r.choose(particles)); } } - tryRandomMove(r); - tryRandomSwap(r); +// randomMove(r); + randomSwap(r); } } void sweepLocal(double z, Rng& r) { for (unsigned i = 0; i < iPow(L, D); i++) { - tryRandomMove(r); +// randomMove(r); } } void sweepSwap(double z, Rng& r) { for (unsigned i = 0; i < iPow(L, D); i++) { - tryRandomSwap(r); + randomSwap(r); } } @@ -313,8 +306,12 @@ int main() { Rng r; + /* for (unsigned N = Nmin; N <= Nmax; N *= 10) { s.setGroundState(); + if (!s.compatible()) { + return 1; + } for (double T = Tmin; T < Tmax; T *= 1 + δT) { double z = exp(1 / T); @@ -335,6 +332,13 @@ int main() { std::cout << L << " " << T << " " << N << " " << δT << " " << 1 / s.density() << std::endl; } } + */ + + s.setGroundState(); + + for (unsigned n = 0; n < 10; n++) { + s.sweepGrandCanonical(exp(1/0.15), r); + } return 0; } |