diff options
| -rw-r--r-- | distinguishable.cpp | 246 | ||||
| -rw-r--r-- | glass.hpp | 230 |
2 files changed, 240 insertions, 236 deletions
diff --git a/distinguishable.cpp b/distinguishable.cpp index 4b66989..611bf50 100644 --- a/distinguishable.cpp +++ b/distinguishable.cpp @@ -1,5 +1,4 @@ #include <iostream> -#include <queue> #include "glass.hpp" @@ -12,262 +11,39 @@ public: DistinguishableState(unsigned t) { type = t; } bool empty() const { return type == 0; } + bool operator==(const DistinguishableState& s) const { return type == s.type; } void remove() { type = 0; } }; -template <unsigned D> class System { +template <unsigned D> +class DistinguishableSystem : public FiniteEnergySystem<D, DistinguishableState> { public: - const unsigned L; - unsigned N; - std::vector<Vertex<D, DistinguishableState>> vertices; std::vector<double> interaction; - 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, Rng& r) : L(L), N(N), vertices(iPow(L, D)), interaction(N * N) { - 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, DistinguishableState>(vertices[j], Δx)); - vertices[j].adjacentEdges.push_back(HalfEdge<D, DistinguishableState>(vertices[i], -Δx)); - } - } - + DistinguishableSystem(unsigned L, unsigned N, Rng& r) + : FiniteEnergySystem<D, DistinguishableState>(L, N), interaction(N * N) { for (double& V : interaction) { V = r.uniform(-0.5, 0.5); } - for (unsigned i = 0; i < N; i++) { - vertices[i].state.type = i + 1; + for (unsigned i = 0; i < this->N; i++) { + this->vertices[i].state.type = i + 1; } } - double pairEnergy(const DistinguishableState& s1, const DistinguishableState& s2) const { + double pairEnergy(const DistinguishableState& s1, const DistinguishableState& s2) const override { if (s1.empty() || s2.empty()) { return 0; } else { if (s1.type < s2.type) { - return interaction[(s1.type - 1) * N + (s2.type - 1)]; + return interaction[(s1.type - 1) * this->N + (s2.type - 1)]; } else { - return interaction[(s2.type - 1) * N + (s1.type - 1)]; + return interaction[(s2.type - 1) * this->N + (s1.type - 1)]; } } } - - double siteEnergy(const Vertex<D, DistinguishableState>& v) const { - double E = 0; - - for (const HalfEdge<D, DistinguishableState>& e : v.adjacentEdges) { - E += pairEnergy(v.state, e.neighbor.state); - } - - return E; - } - - double energy() const { - double E = 0; - - for (const Vertex<D, DistinguishableState>& v : vertices) { - E += siteEnergy(v); - } - - return E / 2; - } - - bool trySwap(Vertex<D, DistinguishableState>& v1, Vertex<D, DistinguishableState>& 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, DistinguishableState>& v1 = r.pick(vertices); - - if (v1.empty()) { - return false; - } - - Vertex<D, DistinguishableState>& 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, DistinguishableState>& v1 = r.pick(vertices); - Vertex<D, DistinguishableState>& v2 = r.pick(vertices); - - if (v1.state.type != v2.state.type) { - return trySwap(v1, v2, T, r); - } else { - return false; - } - } - - double density() const { return N / pow(L, D); } - - void sweepLocal(Rng& r) { - for (unsigned i = 0; i < iPow(L, D); i++) { - tryRandomMove(r); - } - } - - void sweepSwap(Rng& r) { - for (unsigned i = 0; i < iPow(L, D); i++) { - tryRandomSwap(r); - } - } - - unsigned flipCluster(const Transformation<D>& R, Vertex<D, DistinguishableState>& v0, double T, - Rng& r, bool dry = false) { - std::queue<std::array<std::reference_wrapper<Vertex<D, DistinguishableState>>, 2>> q; - Vector<D> x0New = R.apply(v0.position); - Vertex<D, DistinguishableState>& v0New = vertices[vectorToIndex(x0New)]; - q.push({v0, v0New}); - - unsigned n = 0; - - while (!q.empty()) { - auto [vR, vNewR] = q.front(); - q.pop(); - - Vertex<D, DistinguishableState>& v = vR; - Vertex<D, DistinguishableState>& vNew = vNewR; - - if (!v.marked && !vNew.marked) { - v.marked = true; - vNew.marked = true; - - for (HalfEdge<D, DistinguishableState>& e : v.adjacentEdges) { - Vertex<D, DistinguishableState>& vn = e.neighbor; - - Vector<D> xnNew = R.apply(vn.position); - Vertex<D, DistinguishableState>& vnNew = vertices[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(vertices), T, r); - for (Vertex<D, DistinguishableState>& v : vertices) { - v.marked = false; - } - return n; - } - - void swendsenWang(const Transformation<D>& R, double T, Rng& r) { - for (Vertex<D, DistinguishableState>& v : vertices) { - if (!v.marked) { - bool dry = 0.5 < r.uniform(0.0, 1.0); - flipCluster(R, v, T, r, dry); - } - } - - for (Vertex<D, DistinguishableState>& v : vertices) { - v.marked = false; - } - } - - int overlap(const System<D>& s) const { - int o = 0; - - for (unsigned i = 0; i < vertices.size(); i++) { - // DistinguishableState<D> s2 = - // orientation.apply(s.vertices[vectorToIndex(orientation.inverse().apply(indexToVector(i)))].state); - // o += vertices[i].state.dot(s2); - } - - return o; - } - - void setInitialPosition() { - for (Vertex<D, DistinguishableState>& 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, DistinguishableState>& v : vertices) { - if (!v.empty()) { - F += cos(k1.dot(v.position - v.initialPosition)); - F += cos(k2.dot(v.position - v.initialPosition)); - } - } - - return F / (2 * N); - } }; -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 = 2; unsigned L = 40; @@ -278,7 +54,7 @@ int main() { Rng r; - System<D> s(L, N, r); + DistinguishableSystem<D> s(L, N, r); std::vector<Matrix<D>> ms = generateTorusMatrices<D>(); @@ -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; + } + } +}; + |