From 496dcbd9960677db246a84bcd3e4b4230ee28e0a Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 16 Sep 2021 16:47:05 +0200 Subject: Initial commit. --- .gitignore | 8 ++ .gitmodules | 6 ++ animation.hpp | 39 ++++++++ compile_flags.txt | 3 + hard_spheres.cpp | 41 ++++++++ pcg-cpp | 1 + randutils | 1 + soft_spheres.cpp | 44 +++++++++ spheres.hpp | 289 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 9 files changed, 432 insertions(+) create mode 100644 .gitignore create mode 100644 .gitmodules create mode 100644 animation.hpp create mode 100644 compile_flags.txt create mode 100644 hard_spheres.cpp create mode 160000 pcg-cpp create mode 160000 randutils create mode 100644 soft_spheres.cpp create mode 100644 spheres.hpp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..fb16dee --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +* +!*.txt +!*.hpp +!*.cpp +!.gitignore +!.gitmodules +*/ + diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..d9c0c06 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,6 @@ +[submodule "randutils"] + path = randutils + url = https://gist.github.com/makokal/3bf4f1f66b9686384f75 +[submodule "pcg-cpp"] + path = pcg-cpp + url = https://github.com/imneme/pcg-cpp diff --git a/animation.hpp b/animation.hpp new file mode 100644 index 0000000..26b3d55 --- /dev/null +++ b/animation.hpp @@ -0,0 +1,39 @@ +#pragma once + +#include "spheres.hpp" +#include + +void draw(const Sphere<2>& s, std::array color = {0, 0, 0}, + unsigned nPoints = 20) { + glColor3d(color[0], color[1], color[2]); + glBegin(GL_POLYGON); + for (unsigned i = 0; i < nPoints; i++) { + glVertex2d(s.x(0) + s.r * cos(2 * i * M_PI / nPoints), s.x(1) + s.r * sin(2 * i * M_PI / nPoints)); + } + glEnd(); +} + +void initializeAnimation(int argc, char** argv, unsigned window_size = 1000) { + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(window_size, window_size); + glutCreateWindow("spheres"); + glClearColor(0.0, 0.0, 0.0, 0.0); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(0, 1.0, 0, 1.0); + glClearColor(1.0f, 1.0f, 1.0f, 1.0f); +} + +template 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()) { + 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}); + } + } + } + glFlush(); +} diff --git a/compile_flags.txt b/compile_flags.txt new file mode 100644 index 0000000..02d9f07 --- /dev/null +++ b/compile_flags.txt @@ -0,0 +1,3 @@ +-std=c++20 +-O3 +-I/usr/include/eigen3 diff --git a/hard_spheres.cpp b/hard_spheres.cpp new file mode 100644 index 0000000..b53f7cf --- /dev/null +++ b/hard_spheres.cpp @@ -0,0 +1,41 @@ +#include "animation.hpp" + +int main(int argc, char* argv[]) { + unsigned N = 1200; + double R = 0.021; + double ε = 0.01; + unsigned steps = 1e6; + + initializeAnimation(argc, argv); + + Rng r; + + Model<2, HardSphere<2>> m(N, 2 * 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 * N; i++) { + HardSphere<2>& s = r.pick(m.particles); + Vector<2> Δx; + for (double& Δxi : Δx) { + Δxi = r.variate(0, ε); + } + + m.metropolis(1, s, Δx, r); + + HardSphere<2>& s1 = r.pick(m.particles); + HardSphere<2>& s2 = r.pick(m.particles); + + m.swap(1, s1, s2, r); + + if (i % (N * 3) == 0) { + draw(m); + } + } + + return 0; +} diff --git a/pcg-cpp b/pcg-cpp new file mode 160000 index 0000000..ffd522e --- /dev/null +++ b/pcg-cpp @@ -0,0 +1 @@ +Subproject commit ffd522e7188bef30a00c74dc7eb9de5faff90092 diff --git a/randutils b/randutils new file mode 160000 index 0000000..8486a61 --- /dev/null +++ b/randutils @@ -0,0 +1 @@ +Subproject commit 8486a610a954a8248c12485fb4cfc390a5f5f854 diff --git a/soft_spheres.cpp b/soft_spheres.cpp new file mode 100644 index 0000000..a71a472 --- /dev/null +++ b/soft_spheres.cpp @@ -0,0 +1,44 @@ +#include "animation.hpp" + +int main(int argc, char* argv[]) { + const unsigned n = 12; + + unsigned N = 1200; + double R = 0.024; + double ε = 0.01; + unsigned steps = 1e6; + double β = 1; + + initializeAnimation(argc, argv); + + Rng r; + + Model<2, SoftSphere<2, n>> m(N, 2 * R * 1.25); + + 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(0, ε); + } + + m.metropolis(β, s, Δx, r); + + SoftSphere<2, n>& s1 = r.pick(m.particles); + SoftSphere<2, n>& s2 = r.pick(m.particles); + + m.swap(β, s1, s2, r); + + if (i % (N * 3) == 0) { + draw(m); + } + } + + return 0; +} diff --git a/spheres.hpp b/spheres.hpp new file mode 100644 index 0000000..5ecad39 --- /dev/null +++ b/spheres.hpp @@ -0,0 +1,289 @@ +#pragma once + +#include +#include + +#include + +#include "pcg-cpp/include/pcg_random.hpp" +#include "randutils/randutils.hpp" + +using Rng = randutils::random_generator; + +template using Vector = Eigen::Matrix; + +template class Position : public Vector { +private: + void recenter() { + for (double& vi : *this) { + if (vi >= 0) { + vi = fmod(vi, 1); + } else { + vi = fmod(vi + ceil(-vi), 1); + } + } + } + +public: + using Vector::Vector; + + Position& operator+=(const Vector& u) { + Vector::operator+=(u); + recenter(); + return *this; + } + + friend Position operator+(Position v, const Vector& u) { + v += u; + return v; + } + + Position& operator-=(const Vector& u) { + Vector::operator-=(u); + recenter(); + return *this; + } + + friend Vector operator-(Position v, const Position& u) { + v -= (Vector)u; + + for (double& vi : v) { + if (vi > 0.5) + vi = 1 - vi; + } + + return (Vector)v; + } +}; + +template concept Particle = requires(T v, const T& w) { + { v.x } -> std::same_as>; + { v.U(w) } -> std::same_as; +}; + +template T> class NeighborLists { +private: + unsigned n; + std::vector> lists; + + unsigned mod(signed a, unsigned b) const { + if (a >= 0) { + return a % b; + } else { + return (a + b * ((b - a) / b)) % b; + } + } + + int iPow(int x, unsigned p) const { + if (p == 0) + return 1; + if (p == 1) + return x; + + int tmp = iPow(x, p / 2); + if (p % 2 == 0) { + return tmp * tmp; + } else { + return x * tmp * tmp; + } + } + + unsigned index(Position x) const { + unsigned ind = 0; + for (unsigned i = 0; i < D; i++) { + ind += pow(n, i) * std::floor(n * x(i)); + } + return ind; + } + + std::pair::iterator, bool> insert(T& p, unsigned ind) { + return lists[ind].insert(&p); + } + + size_t erase(T& p, unsigned ind) { + return lists[ind].erase(&p); + } + + void neighborsOfHelper(signed i0, unsigned d, std::set& ns) const { + if (d == 0) { + const std::set& a = lists[i0]; + ns.insert(a.begin(), a.end()); + } else { + for (signed j : {-1, 0, 1}) { + signed i1 = iPow(n, d) * (i0 / iPow(n, d)) + mod(i0 + j * iPow(n, d-1), iPow(n, d)); + neighborsOfHelper(i1, d - 1, ns); + } + } + } + +public: + NeighborLists(unsigned m) : n(m), lists(pow(n, D)) {} + NeighborLists(double l) : NeighborLists((unsigned)floor(1 / l)) {} + + std::pair::iterator, bool> insert(T& p, const Position& x) { + return insert(p, index(x)); + } + + std::pair::iterator, bool> insert(T& p) { + return insert(p, p.x); + } + + size_t erase(T& p, const Position& x) { + return erase(p, index(x)); + } + + size_t erase(T& p) { + return erase(p, p.x); + } + + void move(T& p, const Position& xNew) { + unsigned iOld = index(p.x); + unsigned iNew = index(xNew); + if (iNew != iOld) { + erase(p, iOld); + insert(p, iNew); + } + } + + std::set neighborsOf(const Position& x) const { + std::set neighbors; + neighborsOfHelper(index(x), D, neighbors); + return neighbors; + } +}; + +template class Sphere { +public: + Position x; + double r; + + Sphere(const Position& x, double r) : x(x), 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 { +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; + } + } +}; + +template T> class Model { +private: + NeighborLists nl; + +public: + std::vector particles; + + Model(unsigned N, double l) : nl(l) { + particles.reserve(N); + }; + + 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 { + double ΔE = 0; + for (const T* neighbor : neighborsOf(p.x)) { + if (neighbor != &p) { + ΔE -= neighbor->U(p); + } + } + for (const T* neighbor : neighborsOf(pNew.x)) { + if (neighbor != &p) { + ΔE += neighbor->U(pNew); + } + } + + return ΔE; + } + + double energyChangeSwap(const T& p1, const T& p2) const { + double ΔE = 0; + for (const T* neighbor : neighborsOf(p1.x)) { + if (neighbor != &p1 && neighbor != &p2) { + ΔE -= neighbor->U(p1); + ΔE += neighbor->U({p1.x, p2.r}); + } + } + for (const T* neighbor : neighborsOf(p2.x)) { + if (neighbor != &p1 && neighbor != &p2) { + ΔE -= neighbor->U(p2); + ΔE += neighbor->U({p2.x, p1.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); + 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; + } + } +}; -- cgit v1.2.3-70-g09d2