From 4c3a0d01a979ef186b9dbec87323dc282fb3cdc0 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sat, 20 Mar 2021 19:05:44 +0100 Subject: Initial commit. --- glass.cpp | 341 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 341 insertions(+) create mode 100644 glass.cpp diff --git a/glass.cpp b/glass.cpp new file mode 100644 index 0000000..e0c5e25 --- /dev/null +++ b/glass.cpp @@ -0,0 +1,341 @@ +#include +#include +#include +#include +#include + +#include "pcg-cpp/include/pcg_random.hpp" +#include "randutils/randutils.hpp" + +using Rng = randutils::random_generator; + +template using Vector = Eigen::Matrix; +template using Matrix = Eigen::Matrix; + +int iPow(int x, unsigned p) { + if (p == 0) return 1; + if (p == 1) return x; + + int tmp = iPow(x, p / 2); + if (p % 2 == 0) return tmp * tmp; + else return x * tmp * tmp; +} + +unsigned mod(signed a, unsigned b) { + return ((a %= b) < 0) ? a + b : a; +} + +template class State { + public: + Vector σ; + + State(unsigned a, signed b) : σ(Vector::Zero()) { + σ[a] = b; + } + + State apply(const Matrix& m) const { + return {m * σ}; + } +}; + +template class Vertex; + +template class Particle { + public: + Vertex& position; + State state; + bool marked; + + Particle(Vertex& p, const State& s) : position(p), state(s) {} +}; + +template State randomState(Rng& r) { + return State(r.uniform((unsigned)0, D - 1), r.pick({-1, 1})); +} + +template class HalfEdge { + public: + Vertex& neighbor; + Vector Δx; + + HalfEdge(Vertex& n, const Vector& d) : neighbor(n), Δx(d) {} +}; + +template class Vertex { + public: + Vector position; + std::vector> adjacentEdges; + Particle* occupant; + + bool empty() const { + return occupant == NULL; + } +}; + +template class System { + public: + const unsigned L; + std::vector> vertices; + std::set*> particles; + + unsigned vectorToIndex(const Vector& x) const { + unsigned i = 0; + for (unsigned d = 0; d < D; d++) { + i += x[d] * iPow(L, d); + } + return i; + } + + Vector indexToVector(unsigned i) const { + Vector x; + for (unsigned d = 0; d < D; d++) { + x[d] = (i / iPow(L, d)) % L; + } + return x; + } + + 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 = NULL; + } + + for (unsigned d = 0; d < D; d++) { + Vector Δx = Vector::Zero(); + Δx[d] = 1; + for (signed i = 0; i < iPow(L, D); i++) { + unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1)); + vertices[i].adjacentEdges.push_back(HalfEdge(vertices[j], Δx)); + vertices[j].adjacentEdges.push_back(HalfEdge(vertices[i], -Δx)); + } + } + } + + std::list*> overlaps(const Particle& p, bool excludeSelf = false) const { + std::list*> o; + + if (!p.position.empty() && !excludeSelf) { + o.push_back(p.position.occupant); + } + + for (const HalfEdge& e : p.position.adjacentEdges) { + if (!e.neighbor.empty()) { + if (p.state.σ.dot(e.Δx) == 1 || e.neighbor.occupant->state.σ.dot(e.Δx) == -1) { + o.push_back(e.neighbor.occupant); + } + } + } + + return o; + } + + bool insert(const Particle& p) { + if (overlaps(p).empty()) { + Particle* pN = new Particle(p); + particles.insert(pN); + p.position.occupant = pN; + return true; + } else { + return false; + } + } + + void remove(Particle* p) { + p->position.occupant = NULL; + particles.erase(p); + delete p; + } + + /* + bool randomMove(Rng& r) { + Particle& p = *(r.pick(particles)); + Particle pN = p; + + if (0.2 < r.uniform(0.0, 1.0)) { + pN.position = r.pick(p.position.adjacentEdges).neighbor; + } + pN.state = randomState(r); + + if (overlaps(pN, true).empty()) { + remove(&p); + insert(pN); + return true; + } else { + return false; + } + } + */ + + bool swap(Particle& p1, Particle& 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 randomSwap(Rng& r) { + return swap(*r.pick(particles), *r.pick(particles)); + } + + void setGroundState() { + for (Particle* p : particles) { + remove(p); + } + + for (Vertex& v : vertices) { + unsigned a = 0; + for (unsigned d = 0; d < D; d++) { + a += (d + 1) * v.position(d); + } + a %= 2 * D + 1; + + if (a > 0) { + if (a <= D) { + State s(a - 1, -1); + insert(Particle(v, s)); + } else if (D < a) { + State s(2 * D - a, 1); + insert(Particle(v, s)); + } + } + } + } + + bool compatible() { + for (Particle* p : particles) { + if (!overlaps(*p, true).empty()) { + return false; + } + } + + return true; + } + + double density() const { + return N() / pow(L, D); + } + + unsigned maxOccupation() const { +// return (2 * D * iPow(L, D)) / (2 * D + 1); + return iPow(L, D); + } + + 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); + + if (pIns > r.uniform(0.0, 1.0)) { + while (true) { + Vertex& v = r.pick(vertices); + if (v.empty()) { + insert(Particle(v, randomState(r))); + break; + } + } + } + } else { + + double pDel = N() / (z * maxOccupation()); + + if (pDel > r.uniform(0.0, 1.0)) { + remove(r.pick(particles)); + } + } + +// randomMove(r); + randomSwap(r); + } + } + + void sweepLocal(double z, Rng& r) { + for (unsigned i = 0; i < iPow(L, D); i++) { +// randomMove(r); + } + } + + void sweepSwap(double z, Rng& r) { + for (unsigned i = 0; i < iPow(L, D); i++) { + randomSwap(r); + } + } + + void flipCluster(const Matrix& R, Vertex& v0, Rng& r, bool dry = false) { + std::list>> q = {v0}; + + while (!q.empty()) { + Vertex& v = q.front(); + q.pop_front(); + + Vector xNew = R * v.position; + State sNew = v.state.apply(R); + + for (Vertex& vn : overlaps(vertices[vectorToIndex(xNew)], sNew)) { + if (!vn.marked) { + vn.marked = true; + q.push_back(vn); + } + } + } + } +}; + +template Vector randomVector(unsigned L, Rng& r) { + Vector x; + for (unsigned i = 0; i < D; i++) { + x[i] = r.uniform((unsigned)0, L - 1); + } + return x; +} + +int main() { + const unsigned D = 3; + unsigned L = 28; + unsigned Nmin = 2e2; + unsigned Nmax = 2e5; + double Tmin = 0.04; + double Tmax = 0.2; + double δT = 0.02; + + System s(L); + + 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); + + for (unsigned n = 0; n < N; n++) { + s.sweepGrandCanonical(z, r); + } + + std::cout << L << " " << T << " " << N << " " << δT << " " << 1 / s.density() << std::endl; + } + + for (double T = Tmax; T > Tmin; T /= 1 + δT) { + double z = exp(1 / T); + + for (unsigned n = 0; n < N; n++) { + s.sweepGrandCanonical(z, r); + } + + std::cout << L << " " << T << " " << N << " " << δT << " " << 1 / s.density() << std::endl; + } + } + + return 0; +} + + -- cgit v1.2.3-70-g09d2