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;
+    }
+  }
+};