From 3b470b46f600ee56babb71a337f2acfe7b9ae4b3 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 24 Jun 2021 14:08:57 +0200 Subject: More refactoring. --- distinguishable.cpp | 246 +++------------------------------------------------- 1 file changed, 11 insertions(+), 235 deletions(-) (limited to 'distinguishable.cpp') diff --git a/distinguishable.cpp b/distinguishable.cpp index 4b66989..611bf50 100644 --- a/distinguishable.cpp +++ b/distinguishable.cpp @@ -1,5 +1,4 @@ #include -#include #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 class System { +template +class DistinguishableSystem : public FiniteEnergySystem { public: - const unsigned L; - unsigned N; - std::vector> vertices; std::vector interaction; - 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; - } - - 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 Δ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)); - } - } - + DistinguishableSystem(unsigned L, unsigned N, Rng& r) + : FiniteEnergySystem(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& v) const { - double E = 0; - - for (const HalfEdge& e : v.adjacentEdges) { - E += pairEnergy(v.state, e.neighbor.state); - } - - return E; - } - - double energy() const { - double E = 0; - - for (const Vertex& v : vertices) { - E += siteEnergy(v); - } - - return E / 2; - } - - bool trySwap(Vertex& v1, Vertex& 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& v1 = r.pick(vertices); - - if (v1.empty()) { - return false; - } - - Vertex& v2 = (r.pick(v1.adjacentEdges)).neighbor; - - if (!v2.empty()) { - return false; - } - - return trySwap(v1, v2, T, r); - } - - bool tryRandomSwap(double T, Rng& r) { - Vertex& v1 = r.pick(vertices); - Vertex& 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& R, Vertex& v0, double T, - Rng& r, bool dry = false) { - std::queue>, 2>> q; - Vector x0New = R.apply(v0.position); - Vertex& v0New = vertices[vectorToIndex(x0New)]; - q.push({v0, v0New}); - - unsigned n = 0; - - while (!q.empty()) { - auto [vR, vNewR] = q.front(); - q.pop(); - - Vertex& v = vR; - Vertex& vNew = vNewR; - - if (!v.marked && !vNew.marked) { - v.marked = true; - vNew.marked = true; - - for (HalfEdge& e : v.adjacentEdges) { - Vertex& vn = e.neighbor; - - Vector xnNew = R.apply(vn.position); - Vertex& 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& R, double T, Rng& r) { - unsigned n = flipCluster(R, r.pick(vertices), T, r); - for (Vertex& v : vertices) { - v.marked = false; - } - return n; - } - - void swendsenWang(const Transformation& R, double T, Rng& r) { - for (Vertex& v : vertices) { - if (!v.marked) { - bool dry = 0.5 < r.uniform(0.0, 1.0); - flipCluster(R, v, T, r, dry); - } - } - - for (Vertex& v : vertices) { - v.marked = false; - } - } - - int overlap(const System& s) const { - int o = 0; - - for (unsigned i = 0; i < vertices.size(); i++) { - // DistinguishableState s2 = - // orientation.apply(s.vertices[vectorToIndex(orientation.inverse().apply(indexToVector(i)))].state); - // o += vertices[i].state.dot(s2); - } - - return o; - } - - void setInitialPosition() { - for (Vertex& v : vertices) { - v.initialPosition = v.position; - } - } - - double selfIntermediateScattering() const { - double F = 0; - - Vector k1 = {M_PI, 0}; - Vector k2 = {0, M_PI}; - for (const Vertex& 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 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 = 2; unsigned L = 40; @@ -278,7 +54,7 @@ int main() { Rng r; - System s(L, N, r); + DistinguishableSystem s(L, N, r); std::vector> ms = generateTorusMatrices(); -- cgit v1.2.3-54-g00ecf