Some refactoring.

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);
-
-  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));
-  }
+  Model<2, SoftSphere<2, n>> m(N, 2 * R * 1.25, gContinuousPolydisperse(R), r);
 
-  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;
+  };
+}