summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--.gitignore8
-rw-r--r--.gitmodules6
-rw-r--r--animation.hpp39
-rw-r--r--compile_flags.txt3
-rw-r--r--hard_spheres.cpp41
m---------pcg-cpp0
m---------randutils0
-rw-r--r--soft_spheres.cpp44
-rw-r--r--spheres.hpp289
9 files changed, 430 insertions, 0 deletions
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 <GL/glut.h>
+
+void draw(const Sphere<2>& s, std::array<double, 3> 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 <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()) {
+ 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<double, std::normal_distribution>(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
+Subproject ffd522e7188bef30a00c74dc7eb9de5faff9009
diff --git a/randutils b/randutils
new file mode 160000
+Subproject 8486a610a954a8248c12485fb4cfc390a5f5f85
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<double, std::normal_distribution>(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 <set>
+#include <vector>
+
+#include <eigen3/Eigen/Dense>
+
+#include "pcg-cpp/include/pcg_random.hpp"
+#include "randutils/randutils.hpp"
+
+using Rng = randutils::random_generator<pcg32>;
+
+template <int D> using Vector = Eigen::Matrix<double, D, 1>;
+
+template <int D> class Position : public Vector<D> {
+private:
+ void recenter() {
+ for (double& vi : *this) {
+ if (vi >= 0) {
+ vi = fmod(vi, 1);
+ } else {
+ vi = fmod(vi + ceil(-vi), 1);
+ }
+ }
+ }
+
+public:
+ using Vector<D>::Vector;
+
+ Position<D>& operator+=(const Vector<D>& u) {
+ Vector<D>::operator+=(u);
+ recenter();
+ return *this;
+ }
+
+ friend Position<D> operator+(Position<D> v, const Vector<D>& u) {
+ v += u;
+ return v;
+ }
+
+ Position<D>& operator-=(const Vector<D>& u) {
+ Vector<D>::operator-=(u);
+ recenter();
+ return *this;
+ }
+
+ friend Vector<D> operator-(Position<D> v, const Position<D>& u) {
+ v -= (Vector<D>)u;
+
+ for (double& vi : v) {
+ if (vi > 0.5)
+ vi = 1 - vi;
+ }
+
+ return (Vector<D>)v;
+ }
+};
+
+template <typename T, int D> concept Particle = requires(T v, const T& w) {
+ { v.x } -> std::same_as<Position<D>>;
+ { v.U(w) } -> std::same_as<double>;
+};
+
+template <int D, Particle<D> T> class NeighborLists {
+private:
+ unsigned n;
+ std::vector<std::set<T*>> 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<D> 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<typename std::set<T*>::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<T*>& ns) const {
+ if (d == 0) {
+ const std::set<T*>& 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<typename std::set<T*>::iterator, bool> insert(T& p, const Position<D>& x) {
+ return insert(p, index(x));
+ }
+
+ std::pair<typename std::set<T*>::iterator, bool> insert(T& p) {
+ return insert(p, p.x);
+ }
+
+ size_t erase(T& p, const Position<D>& x) {
+ return erase(p, index(x));
+ }
+
+ size_t erase(T& p) {
+ return erase(p, p.x);
+ }
+
+ void move(T& p, const Position<D>& xNew) {
+ unsigned iOld = index(p.x);
+ unsigned iNew = index(xNew);
+ if (iNew != iOld) {
+ erase(p, iOld);
+ insert(p, iNew);
+ }
+ }
+
+ std::set<T*> neighborsOf(const Position<D>& x) const {
+ std::set<T*> neighbors;
+ neighborsOfHelper(index(x), D, neighbors);
+ return neighbors;
+ }
+};
+
+template <int D> class Sphere {
+public:
+ Position<D> x;
+ double r;
+
+ Sphere(const Position<D>& x, double r) : x(x), 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> {
+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;
+ }
+ }
+};
+
+template <int D, Particle<D> T> class Model {
+private:
+ NeighborLists<D, T> nl;
+
+public:
+ std::vector<T> 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<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 {
+ 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<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);
+ 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;
+ }
+ }
+};