From 3140c22dc32ae525b491bcb9653ec755762321d2 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 17 Sep 2021 12:40:44 +0200 Subject: Some refactoring. --- spheres.hpp | 258 +++++++++++++++++++++++++++++++----------------------------- 1 file changed, 135 insertions(+), 123 deletions(-) (limited to 'spheres.hpp') diff --git a/spheres.hpp b/spheres.hpp index ef6b665..6dd5b04 100644 --- a/spheres.hpp +++ b/spheres.hpp @@ -59,10 +59,30 @@ public: } }; -template concept Particle = requires(T v, const T& w) { +template class Euclidean { +public: + Vector t; + Matrix r; + + Euclidean() { + t.setZero(); + r.setIdentity(); + } + + Euclidean(const Vector& t0, const Matrix& r0) : t(t0), r(r0) {} + + Euclidean inverse() const { return Euclidean(-r.transpose() * t, r.transpose()); } + + Vector act(const Vector& x) const { return r * x + t; } + Position act(const Position& x) const { return (Position)(r * x) + t; } + Euclidean act(const Euclidean& R) const { return {act(R.t), r * R.r}; } +}; + +template concept Particle = requires(T v, const Position& x, double r) { { v.x } -> std::same_as>; + { v.r } -> std::same_as; { v.m } -> std::same_as; - { v.U(w) } -> std::same_as; + { v.U(x, r) } -> std::same_as; }; template T> class NeighborLists { @@ -157,162 +177,97 @@ public: } }; -template class Sphere { -public: - Position x; - bool m; - double r; - - Sphere(const Position& x, double r) : x(x), m(false), r(r) {}; - - double regularizedDistanceSquared(const Sphere& s) const { - return (x - s.x).squaredNorm() / pow(r + s.r, 2); - } - - bool intersectsBoundary() const { - for (double xi : x) { - if (xi < r || 1 - xi < r) { - return true; - } - } - return false; - } -}; - -template class HardSphere : public Sphere { -public: - using Sphere::Sphere; - - double U(const HardSphere& s) const { - if (Sphere::regularizedDistanceSquared(s) < 1) { - return 1e6; - } else { - return 0; - } - } -}; - -template class SoftSphere : public Sphere { +template T> class Model { private: - const double c0 = -pow(2, 2*n-3) * pow(5, -n) * (8 + 6*n + pow(n, 2)); - const double c2 = pow(2, 2+2*n) * pow(5, -n-2) * n * (4+n); - const double c4 = -pow(2, 5+2*n) * pow(5, -n-4) * n * (2+n); -public: - using Sphere::Sphere; - - double U(const SoftSphere& s) const { - double r2 = Sphere::regularizedDistanceSquared(s); - if (r2 < pow(1.25, 2)) { - return pow(1 / r2, n / 2) + c0 + c2 * r2 + c4 * pow(r2, 2); - } else { - return 0; - } + bool metropolisAccept(double β, double ΔE, Rng& r) const { + return ΔE < 0 || exp(-β * ΔE) > r.uniform(0.0, 1.0); } -}; -template class Euclidean { -public: - Vector t; - Matrix r; - - Euclidean() { - t.setZero(); - r.setIdentity(); + bool swapAccept(const T& p1, const T& p2) const { + double rat = p1.r / p2.r; + return (rat < 1.2 && rat > 0.83); } - Euclidean(const Vector& t0, const Matrix& r0) : t(t0), r(r0) {} - - Euclidean inverse() const { return Euclidean(-r.transpose() * t, r.transpose()); } - - Vector act(const Vector& x) const { return r * x; } - Position act(const Position& x) const { return (Position)(r * x) + t; } -}; - -template T> class Model { -private: - NeighborLists nl; - public: + NeighborLists nl; std::vector particles; + Euclidean orientation; - Model(unsigned N, double l) : nl(l) { - particles.reserve(N); + Model(unsigned N, double l, std::function g, Rng& r) : nl(l), particles(N) { + for (T& p : particles) { + for (double& xi : p.x) { + xi = r.uniform(0.0, 1.0); + } + p.r = g(r); + nl.insert(p); + } }; - void insert(const T& p) { - particles.push_back(p); - nl.insert(particles.back()); - } - void move(T& p, const Position& x) { nl.move(p, x); p.x = x; } - std::set neighborsOf(const Position& x) const { - return nl.neighborsOf(x); - } - - double energyChange(const T& p, const T& pNew) const { + bool move(double β, T& p, const Vector& Δx, Rng& r) { + Position xNew = p.x + Δx; double ΔE = 0; - for (const T* neighbor : neighborsOf(p.x)) { + + for (const T* neighbor : nl.neighborsOf(p.x)) { if (neighbor != &p) { - ΔE -= neighbor->U(p); + ΔE -= neighbor->U(p.x, p.r); } } - for (const T* neighbor : neighborsOf(pNew.x)) { + + for (const T* neighbor : nl.neighborsOf(xNew)) { if (neighbor != &p) { - ΔE += neighbor->U(pNew); + ΔE += neighbor->U(xNew, p.r); } } - return ΔE; + if (metropolisAccept(β, ΔE, r)) { + move(p, xNew); + return true; + } else { + return false; + } } - double energyChangeSwap(const T& p1, const T& p2) const { + bool swap(double β, T& p1, T& p2, Rng& r) { + if (!swapAccept(p1, p2)) + return false; + double ΔE = 0; - T p1New = p1; - T p2New = p2; - p1New.r = p2.r; - p2New.r = p1.r; - for (const T* neighbor : neighborsOf(p1.x)) { + + for (const T* neighbor : nl.neighborsOf(p1.x)) { if (neighbor != &p1 && neighbor != &p2) { - ΔE -= neighbor->U(p1); - ΔE += neighbor->U(p1New); + ΔE += neighbor->U(p1.x, p2.r) - neighbor->U(p1.x, p1.r); } } - for (const T* neighbor : neighborsOf(p2.x)) { + + for (const T* neighbor : nl.neighborsOf(p2.x)) { if (neighbor != &p1 && neighbor != &p2) { - ΔE -= neighbor->U(p2); - ΔE += neighbor->U(p2New); + ΔE += neighbor->U(p2.x, p1.r) - neighbor->U(p2.x, p2.r); } } - return ΔE; - } - - bool metropolis(double β, T& p, const Vector& Δx, Rng& r) { - T pNew = p; - pNew.x += Δx; - - double ΔE = energyChange(p, pNew); - - if (exp(-β * ΔE) > r.uniform(0.0, 1.0)) { - move(p, pNew.x); + if (metropolisAccept(β, ΔE, r)) { + std::swap(p1.r, p2.r); return true; } else { return false; } } - bool swap(double β, T& p1, T& p2, Rng& r) { - double ΔE = energyChangeSwap(p1, p2); - if (exp(-β * ΔE) > r.uniform(0.0, 1.0)) { - std::swap(p1.r, p2.r); - return true; - } else { - return false; + bool randomSwap(double β, Rng& r) { + return swap(β, r.pick(particles), r.pick(particles), r); + } + + bool randomMove(double β, double δ, Rng& r) { + Vector Δx; + for (double& Δxi : Δx) { + Δxi = r.variate(0, δ); } + return move(β, r.pick(particles), Δx, r); } unsigned clusterFlip(double β, const Euclidean& R, T& p0, Rng& r) { @@ -326,19 +281,18 @@ public: if (!p->m) { p->m = true; - T pNew = *p; - pNew.x = R.act(p->x); + Position xNew = R.act(p->x); - for (T* neighbor : neighborsOf(pNew.x)) { + for (T* neighbor : nl.neighborsOf(xNew)) { if (!neighbor->m) { - double ΔE = neighbor->U(pNew) - neighbor->U(*p); + double ΔE = neighbor->U(xNew, p->r) - neighbor->U(p->x, p->r); if (1 - exp(-β * ΔE) > r.uniform(0.0, 1.0)) { q.push(neighbor); } } } - move(*p, pNew.x); + move(*p, xNew); n++; } } @@ -347,6 +301,64 @@ public: p.m = false; } + if (n > particles.size() / 2) { + orientation = R.act(orientation); + } + return n; } + + unsigned clusterSwap(double β, T& s1, T& s2, Rng& r) { + if (!swapAccept(s1, s2)) + return 0; + + Vector<2> t1 = s1.x; + Vector<2> t2 = s2.x; + Vector<2> t = (t1 + t2) / 2; + + Matrix<2> mat; + + mat(0, 0) = -1; + mat(1, 1) = -1; + mat(0, 1) = 0; + mat(1, 0) = 0; + + Euclidean<2> g(t - mat * t, mat); + + return clusterFlip(β, g, s1, r); + } + + unsigned randomClusterSwap(double β, Rng& r) { + return clusterSwap(β, r.pick(particles), r.pick(particles), r); + } +}; + +template class Sphere { +public: + Position x; + bool m; + double r; + + Sphere() : m(false) {}; + + Sphere(const Position& x, double r) : x(x), m(false), r(r) {}; + + double regularizedDistanceSquared(const Position& y, double ry) const { + return (x - y).squaredNorm() / pow(r + ry, 2); + } + + bool intersectsBoundary() const { + for (double xi : x) { + if (xi < r || 1 - xi < r) { + return true; + } + } + return false; + } }; + +std::function gContinuousPolydisperse(double Rmax) { + return [Rmax] (Rng& r) -> double { + return pow(r.uniform(pow(2.219, -2), 1.0), -0.5) * Rmax / 2.219; + }; +} -- cgit v1.2.3-54-g00ecf