diff options
| -rw-r--r-- | animation.hpp | 8 | ||||
| -rw-r--r-- | hard_spheres.cpp | 99 | ||||
| -rw-r--r-- | soft_spheres.cpp | 67 | ||||
| -rw-r--r-- | spheres.hpp | 258 |
4 files changed, 234 insertions, 198 deletions
diff --git a/animation.hpp b/animation.hpp index 26b3d55..b1ff815 100644 --- a/animation.hpp +++ b/animation.hpp @@ -28,10 +28,12 @@ void initializeAnimation(int argc, char** argv, unsigned window_size = 1000) { template <Particle<2> T> void draw(const Model<2, T>& m) { glClear(GL_COLOR_BUFFER_BIT); for (const Sphere<2>& p : m.particles) { - draw(p); - if (p.intersectsBoundary()) { + Sphere<2> pTmp = p; + pTmp.x = m.orientation.inverse().act(p.x); + draw(pTmp); + if (pTmp.intersectsBoundary()) { for (Vector<2> Δy : {(Vector<2>){1,0}, {-1,0}, {0,1}, {0,-1}, {-1,1}, {1,1}, {1,-1},{-1,-1}}) { - draw({(Vector<2>)p.x + Δy, p.r}); + draw({(Vector<2>)pTmp.x + Δy, p.r}); } } } diff --git a/hard_spheres.cpp b/hard_spheres.cpp index 6a14953..b8be1cf 100644 --- a/hard_spheres.cpp +++ b/hard_spheres.cpp @@ -1,53 +1,82 @@ #include "animation.hpp" +template <int D> class HardSphere : public Sphere<D> { +public: + using Sphere<D>::Sphere; + + double U(const Position<D>& y, double ry) const { + if (Sphere<D>::regularizedDistanceSquared(y, ry) < 1) { + return 1e6; + } else { + return 0; + } + } +}; + int main(int argc, char* argv[]) { unsigned N = 1200; - double R = 0.021; - double ε = 0.01; + double R = 0.024; + double ε = 0.005; unsigned steps = 1e6; initializeAnimation(argc, argv); Rng r; - Model<2, HardSphere<2>> m(N, 2 * R); + Model<2, HardSphere<2>> m(N, 2 * R, gContinuousPolydisperse(R), r); - for (unsigned i = 0; i < N; i++) { - Position<2> x = {r.uniform(0.0, 1.0), r.uniform(0.0, 1.0)}; - double rad = pow(r.uniform(pow(2.219, -2), 1.0), -0.5) * R / 2.219; - m.insert(HardSphere<2>(x, rad)); - } + for (unsigned i = 0; i < steps; i++) { + for (unsigned j = 0; j < N; j++) { + HardSphere<2>& p = r.pick(m.particles); + Position<2> xNew = p.x; + Vector<2> Δx; + for (double& xi : Δx) { + xi = r.variate<double, std::normal_distribution>(0, ε); + } + xNew += Δx; + bool reject = false; + for (const HardSphere<2>* neighbor : m.nl.neighborsOf(xNew)) { + if (neighbor != &p) { + if (neighbor->regularizedDistanceSquared(xNew, p.r) < 1 && neighbor->regularizedDistanceSquared(p.x, p.r) >= 1) { + reject = true; + break; + } + } + } + if (!reject) { + m.move(p, xNew); + } + + HardSphere<2>& p1 = r.pick(m.particles); + HardSphere<2>& p2 = r.pick(m.particles); - for (unsigned i = 0; i < steps * N; i++) { - HardSphere<2>& s = r.pick(m.particles); - Vector<2> Δx; - for (double& Δxi : Δx) { - Δxi = r.variate<double, std::normal_distribution>(0, ε); + double rat = p1.r / p2.r; + if (rat < 1.2 && rat > 0.83) { + reject = false; + for (const HardSphere<2>* neighbor : m.nl.neighborsOf(p2.x)) { + if (neighbor != &p1 && neighbor != &p2) { + if (neighbor->regularizedDistanceSquared(p2.x, p1.r) < 1 && neighbor->regularizedDistanceSquared(p1.x, p1.r) >= 1) { + reject = true; + break; + } + } + } + for (const HardSphere<2>* neighbor : m.nl.neighborsOf(p1.x)) { + if (neighbor != &p1 && neighbor != &p2) { + if (neighbor->regularizedDistanceSquared(p1.x, p2.r) < 1 && neighbor->regularizedDistanceSquared(p2.x, p2.r) >= 1) { + reject = true; + break; + } + } + } + if (!reject) { + std::swap(p1.r, p2.r); + } + } } - m.metropolis(1, s, Δx, r); - - HardSphere<2>& s1 = r.pick(m.particles); - HardSphere<2>& s2 = r.pick(m.particles); - - 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); - - std::cout << m.clusterFlip(1, g, s1, r) << std::endl; - - if (i % (N) == 0) { + if (i % 10 == 0) draw(m); - } } return 0; diff --git a/soft_spheres.cpp b/soft_spheres.cpp index cdd6f31..7944aa0 100644 --- a/soft_spheres.cpp +++ b/soft_spheres.cpp @@ -1,55 +1,48 @@ #include "animation.hpp" +template <int D, int n> class SoftSphere : public Sphere<D> { +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<D>::Sphere; + + double U(const Position<D>& y, double ry) const { + double r2 = Sphere<D>::regularizedDistanceSquared(y, ry); + if (r2 < pow(1.25, 2)) { + return pow(1 / r2, n / 2) + c0 + c2 * r2 + c4 * pow(r2, 2); + } else { + return 0; + } + } +}; + int main(int argc, char* argv[]) { - const unsigned n = 12; + const unsigned n = 24; unsigned N = 1200; - double R = 0.023; - double ε = 0.01; + double R = 0.025; + double ε = 0.005; unsigned steps = 1e6; - double β = 1; + double β = 0.5; initializeAnimation(argc, argv); Rng r; - Model<2, SoftSphere<2, n>> m(N, 2 * R * 1.25); + Model<2, SoftSphere<2, n>> m(N, 2 * R * 1.25, gContinuousPolydisperse(R), r); - for (unsigned i = 0; i < N; i++) { - Position<2> x = {r.uniform(0.0, 1.0), r.uniform(0.0, 1.0)}; - double rad = pow(r.uniform(pow(2.219, -2), 1.0), -0.5) * R / 2.219; - m.insert(SoftSphere<2, n>(x, rad)); - } - - for (unsigned i = 0; i < steps * N; i++) { - SoftSphere<2, n>& s = r.pick(m.particles); - Vector<2> Δx; - for (double& Δxi : Δx) { - Δxi = r.variate<double, std::normal_distribution>(0, ε); + for (unsigned i = 0; i < steps; i++) { + for (unsigned j = 0; j < N; j++) { + m.randomMove(β, ε, r); + m.randomSwap(β, r); } - for (unsigned j = 0; j < N; j++) - m.metropolis(β, s, Δx, r); - - SoftSphere<2, n>& s1 = r.pick(m.particles); - SoftSphere<2, n>& s2 = r.pick(m.particles); - - 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); - - std::cout << m.clusterFlip(1, g, s1, r) << std::endl; +// std::cout << m.randomClusterSwap(β, r) << std::endl; - draw(m); + if (i % 3 == 0) + draw(m); } return 0; diff --git a/spheres.hpp b/spheres.hpp index ef6b665..6dd5b04 100644 --- a/spheres.hpp +++ b/spheres.hpp @@ -59,10 +59,30 @@ public: } }; -template <typename T, int D> concept Particle = requires(T v, const T& w) { +template <int D> class Euclidean { +public: + Vector<D> t; + Matrix<D> r; + + Euclidean() { + t.setZero(); + r.setIdentity(); + } + + Euclidean(const Vector<D>& t0, const Matrix<D>& r0) : t(t0), r(r0) {} + + Euclidean inverse() const { return Euclidean(-r.transpose() * t, r.transpose()); } + + Vector<D> act(const Vector<D>& x) const { return r * x + t; } + Position<D> act(const Position<D>& x) const { return (Position<D>)(r * x) + t; } + Euclidean<D> act(const Euclidean<D>& R) const { return {act(R.t), r * R.r}; } +}; + +template <typename T, int D> concept Particle = requires(T v, const Position<D>& x, double r) { { v.x } -> std::same_as<Position<D>>; + { v.r } -> std::same_as<double>; { v.m } -> std::same_as<bool>; - { v.U(w) } -> std::same_as<double>; + { v.U(x, r) } -> std::same_as<double>; }; template <int D, Particle<D> T> class NeighborLists { @@ -157,162 +177,97 @@ public: } }; -template <int D> class Sphere { -public: - Position<D> x; - bool m; - double r; - - Sphere(const Position<D>& x, double r) : x(x), m(false), r(r) {}; - - double regularizedDistanceSquared(const Sphere<D>& 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 <int D> class HardSphere : public Sphere<D> { -public: - using Sphere<D>::Sphere; - - double U(const HardSphere<D>& s) const { - if (Sphere<D>::regularizedDistanceSquared(s) < 1) { - return 1e6; - } else { - return 0; - } - } -}; - -template <int D, int n> class SoftSphere : public Sphere<D> { +template <int D, Particle<D> 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<D>::Sphere; - - double U(const SoftSphere<D, n>& s) const { - double r2 = Sphere<D>::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 <int D> class Euclidean { -public: - Vector<D> t; - Matrix<D> 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<D>& t0, const Matrix<D>& r0) : t(t0), r(r0) {} - - Euclidean inverse() const { return Euclidean(-r.transpose() * t, r.transpose()); } - - Vector<D> act(const Vector<D>& x) const { return r * x; } - Position<D> act(const Position<D>& x) const { return (Position<D>)(r * x) + t; } -}; - -template <int D, Particle<D> T> class Model { -private: - NeighborLists<D, T> nl; - public: + NeighborLists<D, T> nl; std::vector<T> particles; + Euclidean<D> orientation; - Model(unsigned N, double l) : nl(l) { - particles.reserve(N); + Model(unsigned N, double l, std::function<double(Rng& r)> 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<D>& x) { nl.move(p, x); p.x = x; } - std::set<T*> neighborsOf(const Position<D>& x) const { - return nl.neighborsOf(x); - } - - double energyChange(const T& p, const T& pNew) const { + bool move(double β, T& p, const Vector<D>& Δx, Rng& r) { + Position<D> 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<D>& Δ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<D> Δx; + for (double& Δxi : Δx) { + Δxi = r.variate<double, std::normal_distribution>(0, δ); } + return move(β, r.pick(particles), Δx, r); } unsigned clusterFlip(double β, const Euclidean<D>& 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<D> 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 <int D> class Sphere { +public: + Position<D> x; + bool m; + double r; + + Sphere() : m(false) {}; + + Sphere(const Position<D>& x, double r) : x(x), m(false), r(r) {}; + + double regularizedDistanceSquared(const Position<D>& 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<double(Rng&)> gContinuousPolydisperse(double Rmax) { + return [Rmax] (Rng& r) -> double { + return pow(r.uniform(pow(2.219, -2), 1.0), -0.5) * Rmax / 2.219; + }; +} |