diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2021-06-09 16:32:27 +0200 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2021-06-09 16:32:27 +0200 |
commit | 7548abd44c37d4495dd498193ea15b0df757b904 (patch) | |
tree | ddca95300a8b323192f2f6a8f4459bb08733d3f5 /glass.hpp | |
parent | 2bf71fe50bf2fb0a256b8f1511aafdd5031c47d9 (diff) | |
download | lattice_glass-7548abd44c37d4495dd498193ea15b0df757b904.tar.gz lattice_glass-7548abd44c37d4495dd498193ea15b0df757b904.tar.bz2 lattice_glass-7548abd44c37d4495dd498193ea15b0df757b904.zip |
Lots of work.
Diffstat (limited to 'glass.hpp')
-rw-r--r-- | glass.hpp | 677 |
1 files changed, 0 insertions, 677 deletions
diff --git a/glass.hpp b/glass.hpp deleted file mode 100644 index b35964b..0000000 --- a/glass.hpp +++ /dev/null @@ -1,677 +0,0 @@ -#include <limits> -#include <list> -#include <queue> -#include <vector> - -#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 < 0) ? (a + (1 - a / (signed)b) * b) : a) % b; } - -template <int D, typename Derived> Vector<D> mod(const Eigen::MatrixBase<Derived>& v, unsigned b) { - Vector<D> u; - for (unsigned i = 0; i < D; i++) { - u(i) = mod(v(i), b); - } - return u; -} - -template <int D> void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) { - if (level > 0) { - unsigned new_level = level - 1; - for (std::array<double, D>& sequence : sequences) { - std::array<double, D> new_sequence = sequence; - new_sequence[new_level] = -1; - sequences.push_front(new_sequence); - } - one_sequences<D>(sequences, new_level); - } -} - -template <int D> std::vector<Matrix<D>> generateTorusMatrices() { - std::vector<Matrix<D>> mats; - - std::array<double, D> ini_sequence; - ini_sequence.fill(1); - std::list<std::array<double, D>> sequences; - sequences.push_back(ini_sequence); - - one_sequences<D>(sequences, D); - - sequences.pop_back(); // don't want the identity matrix! - - for (std::array<double, D> sequence : sequences) { - Matrix<D> m; - for (unsigned i = 0; i < D; i++) { - for (unsigned j = 0; j < D; j++) { - if (i == j) { - m(i, j) = sequence[i]; - } else { - m(i, j) = 0; - } - } - } - - mats.push_back(m); - } - - for (unsigned i = 0; i < D; i++) { - for (unsigned j = 0; j < D; j++) { - if (i != j) { - Matrix<D> m; - for (unsigned k = 0; k < D; k++) { - for (unsigned l = 0; l < D; l++) { - if ((k == i && l == j) || (k == j && l == i)) { - if (i < j) { - m(k, l) = 1; - } else { - m(k, l) = -1; - } - } else if (k == l && (k != i && k != j)) { - m(k, l) = 1; - } else { - m(k, l) = 0; - } - } - } - mats.push_back(m); - } - } - } - - return mats; -} - -template <int D> class CiamarraState : public Vector<D> { -public: - CiamarraState() : Vector<D>(Vector<D>::Zero()) {} - - CiamarraState(const Vector<D>& v) : Vector<D>(v) {} - - CiamarraState(unsigned a, signed b) : Vector<D>(Vector<D>::Zero()) { this->operator()(a) = b; } - - CiamarraState(Rng& r) : CiamarraState(r.uniform(0, D - 1), r.pick({-1, 1})) {} - - bool isEmpty() const { return this->squaredNorm() == 0; } - - void remove() { this->setZero(); } - - CiamarraState<D> flip() const { - CiamarraState<D> s; - for (unsigned i = 0; i < D; i++) { - s(i) = -this->operator()(i); - } - return s; - } -}; - -class BiroliState { -public: - unsigned type; - unsigned occupiedNeighbors; - - BiroliState() { - type = std::numeric_limits<unsigned>::max(); - occupiedNeighbors = 0; - } - - BiroliState(unsigned t, unsigned o) { - type = t; - occupiedNeighbors = o; - } - - 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 = std::numeric_limits<unsigned>::max(); }; -}; - -template <int D> class Transformation { -public: - unsigned L; - Matrix<D> m; - Vector<D> v; - - Transformation(unsigned L) : L(L) { - m.setIdentity(); - v.setZero(); - } - - Transformation(unsigned L, const Matrix<D>& m, const Vector<D>& v) : L(L), m(m), v(v) {} - - Transformation(unsigned L, const std::vector<Matrix<D>>& ms, Rng& r) : m(r.pick(ms)), L(L) { - for (unsigned i = 0; i < D; i++) { - v[i] = r.uniform((unsigned)0, L - 1); - } - - v = v - m * v; - } - - Transformation<D> inverse() const { - return Transformation<D>(L, m.transpose(), -m.transpose() * v); - } - - Vector<D> apply(const Vector<D>& x) const { return mod<D>(v + m * x, L); } - - Transformation<D> apply(const Transformation<D>& t) const { - Transformation<D> tNew(L); - - tNew.m = m * t.m; - tNew.v = apply(t.v); - - return tNew; - } - - CiamarraState<D> apply(const CiamarraState<D>& s) const { return CiamarraState<D>(m * s); } - - BiroliState apply(const BiroliState& s) const { return s; } -}; - -template <int D, typename State> class HalfEdge; - -template <int D, typename State> class Vertex { -public: - Vector<D> position; - std::vector<HalfEdge<D, State>> adjacentEdges; - State state; - bool marked; - - bool isEmpty() const { return state.isEmpty(); } -}; - -template <int D, class State> class HalfEdge { -public: - Vertex<D, State>& neighbor; - Vector<D> Δx; - - HalfEdge(Vertex<D, State>& n, const Vector<D>& d) : neighbor(n), Δx(d) {} -}; - -template <int D, class State> class System { -public: - const unsigned L; - unsigned N; - std::vector<Vertex<D, State>> vertices; - Transformation<D> orientation; - - 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; - } - - 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].adjacentEdges.reserve(2 * D); - vertices[i].marked = false; - } - - 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, State>(vertices[j], Δx)); - vertices[j].adjacentEdges.push_back(HalfEdge<D, State>(vertices[i], -Δx)); - } - } - } - - virtual std::list<std::reference_wrapper<Vertex<D, State>>> overlaps(Vertex<D, State>&) { return {}; }; - - virtual bool insert(Vertex<D, State>& v, const State& s, bool force = false) { return false; }; - - virtual bool remove(Vertex<D, State>& v) { return false;}; - - virtual bool randomMove(Rng& r) { return false;}; - - virtual void swap(Vertex<D, State>& v1, Vertex<D, State>& v2) {}; - - virtual bool exchange(Vertex<D, State>& v1, Vertex<D, State>& v2) { return false;}; - - bool randomExchange(Rng& r) { - Vertex<D, State>& v1 = r.pick(vertices); - Vertex<D, State>& v2 = r.pick(vertices); - - return exchange(v1, v2); - } - - bool compatible() { - for (Vertex<D, State>& v : vertices) { - if (overlaps(v).size() > 0) { - return false; - } - } - - return true; - } - - double density() const { return N / pow(L, D); } - - unsigned size() const { return vertices.size(); } - - 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 = size() * z / (N + 1); - - if (pIns > r.uniform(0.0, 1.0)) { - while (true) { - Vertex<D, State>& v = r.pick(vertices); - if (v.isEmpty()) { - insert(v, State(r)); - break; - } - } - } - } else { - - double pDel = N / (z * size()); - - if (pDel > r.uniform(0.0, 1.0)) { - remove(r.pick(vertices)); - } - } - - randomMove(r); - } - } - - void sweepLocal(Rng& r) { - for (unsigned i = 0; i < iPow(L, D); i++) { - randomMove(r); - } - } - - void sweepSwap(Rng& r) { - for (unsigned i = 0; i < iPow(L, D); i++) { - randomExchange(r); - } - } - - unsigned flipCluster(const Transformation<D>& R, Vertex<D, State>& v0, bool dry = false) { - unsigned n = 0; - - Vertex<D, State>& v0New = vertices[vectorToIndex(R.apply(v0.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); - } - } - - while (!q.empty()) { - Vertex<D, State>& v = q.front(); - q.pop(); - - if (!v.marked && !overlaps(v).empty()) { - v.marked = true; - Vertex<D, State>& vNew = vertices[vectorToIndex(R.apply(v.position))]; - - 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); - } - } - - n += 2; - } else { - n += 1; - } - } - 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); - } - } - } - } - - - } - - if (n > size() / 4) { - orientation = R.apply(orientation); - } - - for (Vertex<D, State>& v : vertices) { - v.marked = false; - } - - return n; - } - - void swendsenWang(const Transformation<D>& R, Rng& r) { - for (Vertex<D, State>& v : vertices) { - if (!v.marked) { - bool dry = 0.5 < r.uniform(0.0, 1.0); - unsigned n = flipCluster(R, v, dry); - if (n > size() / 4 && !dry) { - orientation = R.apply(orientation); - } - } - } - - for (Vertex<D, State>& v : vertices) { - v.marked = false; - } - } - - virtual int overlap(const System<D, State>& s) const { return 0; }; -}; - -template <int D> class CiamarraSystem : public System<D, CiamarraState<D>> { - public: - - using System<D, CiamarraState<D>>::System; - std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>> - overlaps(Vertex<D, CiamarraState<D>>& v, const CiamarraState<D>& s, - bool excludeSelf = false) override { - std::list<std::reference_wrapper<Vertex<D, CiamarraState<D>>>> o; - - if (s.isEmpty()) { - return o; - } - - if (!v.isEmpty() && !excludeSelf) { - o.push_back(v); - } - - for (const HalfEdge<D, CiamarraState<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); - } - } - } - - return o; - } - - 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; - } else { - return false; - } - } - - bool remove(Vertex<D, CiamarraState<D>>& v) override { - if (v.isEmpty()) { - return false; - } else { - v.state.remove(); - this->N--; - return true; - } - } - - bool randomMove(Rng& r) override { - Vertex<D, CiamarraState<D>>& v = r.pick(this->vertices); - CiamarraState<D> oldState = v.state; - - if (!remove(v)) { - return false; - } - - if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) { - for (HalfEdge<D, CiamarraState<D>>& e : v.adjacentEdges) { - if (1 == e.Δx.dot(oldState)) { - if (insert(e.neighbor, oldState.flip())) { - return true; - } - break; - } - } - } else { - CiamarraState<D> newState(r); - while (newState.dot(oldState) == 1) { - newState = CiamarraState<D>(r); - } - if (insert(v, newState)) { - return true; - } - } - v.state = oldState; - this->N++; - return false; - } - - void swap(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override { - std::swap(v1.state, v2.state); - } - - bool exchange(Vertex<D, CiamarraState<D>>& v1, Vertex<D, CiamarraState<D>>& v2) override { - if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) { - swap(v1, v2); - return true; - } else { - return false; - } - } - - void setGroundState() { - this->N = 0; - - for (Vertex<D, CiamarraState<D>>& v : this->vertices) { - unsigned a = 0; - for (unsigned d = 0; d < D; d++) { - a += (d + 1) * v.position(d); - } - a %= 2 * D + 1; - - v.state.setZero() = Vector<D>::Zero(); - - if (0 < a && a <= D) { - v.state(a - 1) = -1; - this->N++; - } else if (D < a) { - v.state(2 * D - a) = 1; - this->N++; - } - } - } - - int overlap(const System<D, CiamarraState<D>>& s) const override { - int o = 0; - - for (unsigned i = 0; i < this->vertices.size(); i++) { - CiamarraState<D> s2 = this->orientation.apply( - s.vertices[this->vectorToIndex(this->orientation.inverse().apply(this->indexToVector(i)))].state); - o += this->vertices[i].state.dot(s2); - } - - 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; - } -}; |