1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
|
#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 = 24;
unsigned N = 1200;
double R = 0.025;
double ε = 0.005;
unsigned steps = 1e4;
double β = 0.5;
initializeAnimation(argc, argv);
Rng r;
Model<2, SoftSphere<2, n>> m(N, 2 * R * 1.25, gContinuousPolydisperse(R), r);
for (unsigned i = 0; i < steps; i++) {
for (unsigned j = 0; j < N; j++) {
m.randomMove(β, ε, r);
m.randomSwap(β, r);
}
// std::cout << m.randomClusterSwap(β, r) << std::endl;
if (i % 3 == 0)
draw(m);
}
for (unsigned i = 0; i < steps; i++) {
for (unsigned j = 0; j < N; j++) {
m.randomMove(β, ε, r);
}
unsigned k = m.randomClusterSwap(β, r);
if (k != 0 && k != N) {
std::cout << k << std::endl;
}
if (i % 3 == 0)
draw(m);
}
return 0;
}
|
