diff options
| -rw-r--r-- | glass.cpp | 341 |
1 files changed, 341 insertions, 0 deletions
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 <iostream> +#include <vector> +#include <list> +#include <set> +#include <eigen3/Eigen/Dense> + +#include "pcg-cpp/include/pcg_random.hpp" +#include "randutils/randutils.hpp" + +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>; + +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<unsigned D> class State { + public: + Vector<D> σ; + + State(unsigned a, signed b) : σ(Vector<D>::Zero()) { + σ[a] = b; + } + + State<D> apply(const Matrix<D>& m) const { + return {m * σ}; + } +}; + +template <unsigned D> class Vertex; + +template<unsigned D> class Particle { + public: + Vertex<D>& position; + State<D> state; + bool marked; + + Particle(Vertex<D>& p, const State<D>& s) : position(p), state(s) {} +}; + +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 { + public: + Vertex<D>& neighbor; + Vector<D> Δx; + + 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; + Particle<D>* occupant; + + bool empty() const { + return occupant == NULL; + } +}; + +template<unsigned D> class System { + public: + const unsigned L; + std::vector<Vertex<D>> vertices; + std::set<Particle<D>*> particles; + + unsigned vectorToIndex(const Vector<D>& x) const { + unsigned i = 0; + for (unsigned d = 0; d < D; d++) { + i += x[d] * iPow(L, d); + } + return i; + } + + Vector<D> indexToVector(unsigned i) const { + Vector<D> 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<D> Δx = Vector<D>::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<D>(vertices[j], Δx)); + vertices[j].adjacentEdges.push_back(HalfEdge<D>(vertices[i], -Δx)); + } + } + } + + std::list<Particle<D>*> overlaps(const Particle<D>& p, bool excludeSelf = false) const { + std::list<Particle<D>*> o; + + if (!p.position.empty() && !excludeSelf) { + o.push_back(p.position.occupant); + } + + for (const HalfEdge<D>& 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<D>& p) { + if (overlaps(p).empty()) { + Particle<D>* pN = new Particle<D>(p); + particles.insert(pN); + p.position.occupant = pN; + return true; + } else { + return false; + } + } + + void remove(Particle<D>* p) { + p->position.occupant = NULL; + particles.erase(p); + delete p; + } + + /* + 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 (overlaps(pN, true).empty()) { + remove(&p); + insert(pN); + return true; + } else { + return false; + } + } + */ + + 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 randomSwap(Rng& r) { + return swap(*r.pick(particles), *r.pick(particles)); + } + + void setGroundState() { + for (Particle<D>* p : particles) { + remove(p); + } + + for (Vertex<D>& 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<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 (Particle<D>* 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<D>& v = r.pick(vertices); + if (v.empty()) { + insert(Particle<D>(v, randomState<D>(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<D>& R, Vertex<D>& v0, Rng& r, bool dry = false) { + std::list<std::reference_wrapper<Vertex<D>>> q = {v0}; + + while (!q.empty()) { + Vertex<D>& v = q.front(); + q.pop_front(); + + Vector<D> xNew = R * v.position; + State<D> sNew = v.state.apply(R); + + for (Vertex<D>& vn : overlaps(vertices[vectorToIndex(xNew)], sNew)) { + if (!vn.marked) { + vn.marked = true; + q.push_back(vn); + } + } + } + } +}; + +template<unsigned D> Vector<D> randomVector(unsigned L, Rng& r) { + Vector<D> 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<D> 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; +} + + |