summaryrefslogtreecommitdiff
path: root/spheres.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'spheres.hpp')
-rw-r--r--spheres.hpp258
1 files changed, 135 insertions, 123 deletions
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;
+ };
+}