diff options
Diffstat (limited to 'glass.hpp')
| -rw-r--r-- | glass.hpp | 230 |
1 files changed, 229 insertions, 1 deletions
@@ -1,4 +1,5 @@ #include <list> +#include <queue> #include <eigen3/Eigen/Dense> #include <eigen3/Eigen/src/Core/CwiseNullaryOp.h> @@ -11,6 +12,14 @@ template <int D> using Matrix = Eigen::Matrix<int, D, D>; using Rng = randutils::random_generator<pcg32>; +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 iPow(int x, unsigned p) { if (p == 0) return 1; @@ -139,9 +148,10 @@ public: } }; -template <typename T> concept State = requires(T v) { +template <typename T> concept State = requires(T v, const T& v2) { { v.empty() } -> std::same_as<bool>; {v.remove()}; + { v.operator==(v2) } -> std::same_as<bool>; }; template <unsigned D, State S> class HalfEdge; @@ -165,3 +175,221 @@ public: HalfEdge(Vertex<D, S>& n, const Vector<D>& d) : neighbor(n), Δx(d) {} }; +template <unsigned D, State S> class System { +public: + const unsigned L; + unsigned N; + std::vector<Vertex<D, S>> vertices; + + 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, unsigned N) : L(L), N(N), vertices(iPow(L, D)) { + for (unsigned i = 0; i < iPow(L, D); i++) { + vertices[i].position = indexToVector(i); + vertices[i].initialPosition = vertices[i].position; + 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, S>(vertices[j], Δx)); + vertices[j].adjacentEdges.push_back(HalfEdge<D, S>(vertices[i], -Δx)); + } + } + } + + double density() const { return N / pow(L, D); } + + void setInitialPosition() { + for (Vertex<D, S>& v : vertices) { + v.initialPosition = v.position; + } + } + + double selfIntermediateScattering() const { + double F = 0; + + Vector<D> k1 = {M_PI, 0}; + Vector<D> k2 = {0, M_PI}; + for (const Vertex<D, S>& v : vertices) { + if (!v.empty()) { + F += cos(k1.dot(v.position - v.initialPosition)); + F += cos(k2.dot(v.position - v.initialPosition)); + } + } + + return F / (2 * this->N); + } +}; + +template <unsigned D, State S> class FiniteEnergySystem : public System<D, S> { +public: + using System<D, S>::System; + + virtual double pairEnergy(const S&, const S&) const { return 0; } + + double siteEnergy(const Vertex<D, S>& v) const { + double E = 0; + + for (const HalfEdge<D, S>& e : v.adjacentEdges) { + E += pairEnergy(v.state, e.neighbor.state); + } + + return E; + } + + double energy() const { + double E = 0; + + for (const Vertex<D, S>& v : this->vertices) { + E += siteEnergy(v); + } + + return E / 2; + } + + bool trySwap(Vertex<D, S>& v1, Vertex<D, S>& v2, double T, Rng& r) { + double E0 = siteEnergy(v1) + siteEnergy(v2); + std::swap(v1.state, v2.state); + double E1 = siteEnergy(v1) + siteEnergy(v2); + double ΔE = E1 - E0; + + if (exp(-ΔE / T) > r.uniform(0.0, 1.0)) { + /* Accept the swap. */ + std::swap(v1.initialPosition, v2.initialPosition); + return true; + } else { + /* Revert the swap. */ + std::swap(v1.state, v2.state); + return false; + } + } + + bool tryRandomMove(double T, Rng& r) { + Vertex<D, S>& v1 = r.pick(this->vertices); + + if (v1.empty()) { + return false; + } + + Vertex<D, S>& v2 = (r.pick(v1.adjacentEdges)).neighbor; + + if (!v2.empty()) { + return false; + } + + return trySwap(v1, v2, T, r); + } + + bool tryRandomSwap(double T, Rng& r) { + Vertex<D, S>& v1 = r.pick(this->vertices); + Vertex<D, S>& v2 = r.pick(this->vertices); + + if (v1.state != v2.state) { + return trySwap(v1, v2, T, r); + } else { + return false; + } + } + + void sweepLocal(Rng& r) { + for (unsigned i = 0; i < iPow(this->L, D); i++) { + tryRandomMove(r); + } + } + + void sweepSwap(Rng& r) { + for (unsigned i = 0; i < iPow(this->L, D); i++) { + tryRandomSwap(r); + } + } + + unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, double T, Rng& r, + 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)]; + q.push({v0, v0New}); + + unsigned n = 0; + + while (!q.empty()) { + auto [vR, vNewR] = q.front(); + q.pop(); + + Vertex<D, S>& v = vR; + Vertex<D, S>& vNew = vNewR; + + if (!v.marked && !vNew.marked) { + v.marked = true; + vNew.marked = true; + + for (HalfEdge<D, S>& e : v.adjacentEdges) { + Vertex<D, S>& vn = e.neighbor; + + Vector<D> xnNew = R.apply(vn.position); + Vertex<D, S>& vnNew = this->vertices[this->vectorToIndex(xnNew)]; + + if (!vn.marked && !vnNew.marked) { + double E0 = pairEnergy(v.state, vn.state) + pairEnergy(vNew.state, vnNew.state); + double E1 = pairEnergy(vNew.state, vn.state) + pairEnergy(v.state, vnNew.state); + double ΔE = E1 - E0; + + if (exp(-ΔE / T) < r.uniform(0.0, 1.0)) { + q.push({vn, vnNew}); + } + } + } + + if (!dry) { + std::swap(v.state, vNew.state); + std::swap(v.initialPosition, vNew.initialPosition); + } + + n += 1; + } + } + + return n; + } + + unsigned wolff(const Transformation<D>& R, double T, Rng& r) { + unsigned n = flipCluster(R, r.pick(this->vertices), T, r); + for (Vertex<D, S>& v : this->vertices) { + v.marked = false; + } + return n; + } + + void swendsenWang(const Transformation<D>& R, double T, Rng& r) { + for (Vertex<D, S>& v : this->vertices) { + if (!v.marked) { + bool dry = 0.5 < r.uniform(0.0, 1.0); + flipCluster(R, v, T, r, dry); + } + } + + for (Vertex<D, S>& v : this->vertices) { + v.marked = false; + } + } +}; + |