Hard_Spheres_MD_MC

=The hard-sphere model (short history)=

The study of hard-sphere systems goes back a long time. Many people see a precursor in the Roman poet and philosopher [|Lucretius]. Another mile-stone was Daniel Bernoulli's (1700 – 1782) discussion of the pressure dependence of a hard-sphere gas (1738). Boltzmann worked on hard spheres, so did Maxwell, and many researchers since then. The phase transition in hard disks, discovered by Alder and Wainwright in 1962 is an important discovery made by numerical simulations. The true understanding of this transition dates from 2011.

=Molecular Dynamics (animation)=



The Event-driven algorithm (scheme)


Here is a cartoon of the time-evolution of four hard disks in a square box. At time t=0, each disk starts at a given position and with a given velocity. The entire time-evolution is simply the solution of Newton's equations: Each disk moves freely until it collides either with another disk or with a wall (see ). In this simulation, time is continuous, and the time for the next "event" can be computed exactly: it is given by the minimum of //N//(//N//-1)/2 pair collision times for isolated pairs of particles, without walls, and N wall-collision times (for isolated disks). The figure on the left is SMAC fig. 2.1. Molecular dynamics was invented by Alder and Wainwright, in 1957.

Wall collisions, pair collisions
Wall collisions are trivial. Pair collisions (with periodic boundary conditions) will be treated and programmed in this week's practical session.

Chaos
The event-driven molecular dynamics algorithm has no time-step error, and the only source of error comes from the finite precision of the arithmetic. These errors are magnified from iteration to iteration. The cartoons on the left illustrates the influence of tiny rounding errors on the dynamics. This is SMAC fig. 2.5. This extreme influence on the initial conditions is called chaos. In this simulation, chaos has two consequences:
 * Molecular dynamics cannot really give rigorous results for the dynamics: How can we compute the exact positions of particles at //t// = 100000 if for //t// = 33, we already need double precision?
 * The hard-sphere system is chaotic (unlike the Fermi-Pasta-Ulam chain, subject of last week's homework) for all radii. Chaos implies that statistical mechanics applies to a finite system of hard spheres and hard disks. This was proven in fundamental theorems by Sinai and Simanyi . The mathematical problems in these works are formidable, however, and the situation does not seem to be fully under control.

Complexity of hard-sphere MD: Local Cells, Priority queues
The naive calculation of the next collision in the system does not take into account that upcoming events can be computed ahead of time: The event at time 5.84 (in the above figure of the event-driven algorithm) was already computed at time 0. In addition, especially at appreciable density, collisions are local.

In modern implementations of the event-driven algorithm, one thus uses a socalled "Priority queue", an abstract data structure often implemented as a "heap". This allows to compute the minimum element very easily (in //O//(1) operations, and to add new ones, in log //N// operations for a heap of //N// elements. The initialization of the heap takes log //N// operations. [|Here's a video lecture on priority queues] by Berkeley Computer Scientist J. Shewchuk on priority queues. It nicely explains the logarithmic complexity of the insertion step.

We thus have the following complexity of the algorithm: > For the special case of hard-disk systems, there is even an //O//(1) per event algorithm although it is currently not faster than the priority queue algorithm.
 * naive event-driven molecular dynamics algorithm, as in SMAC Algorithm 2.1 event-disks on the order of //N// 2 operations per event
 * event-driven molecular dynamics algorithm using local cells, as discussed in SMAC sect. 2.4.1 on the order of //N// operations per event
 * event-driven molecular dynamics algorithm using priority queues, of the order of log //N// per event (!)

Here is an illustration of priority queues (heaps) in Python, in a context similar to an event-driven MD simulation. Some events in the heap are "False". This corresponds to collisions which no longer exist. It would be more expensive to take them out of the heap than to simply mark them as "invalid".

code format="python" import random import heapq def heap_simulation: a = random.random, 1] for i in xrange(100000)]   heapq.heapify(a)    for iter in xrange(100000):        for i in range(2):            j  =random.randint(0, len(a) - 1)            a[j][1] = 0        for i in range(2): heapq.heappush(a, [random.random, 1])        while True:            d = heapq.heappop(a)            if d[1] == 1: break    return cProfile.run('heap_simulation') [[code =Monte Carlo Algorithms=

Equiprobability


The equal-probability principle translates the Boltzmann distribution for a system whose energy is either zero or infinity: math \pi(x_1 \dots x_N) = \begin{cases} 1\quad \text{if configuration legal} \\ 0 \quad \text{otherwise} \end{cases} math

The concept of "equiprobability" often poses problems. Notice that:
 * In the above picture, equiprobability is a statement in 18 dimensions (//d// //N//, with //d// the dimension and //N// the number of spheres): The 18-dimensional probability distribution, in Cartesian coordinates, is uniform, but not its lower-dimensional projections.
 * The equal-probability principle (= the Boltzmann distribution for hard spheres) is valid independent of the boundary conditions.
 * The equal-probability principle neither contradicts the presence of indirect attraction between particles nor the existence of phase transitions.

Direct-sampling algorithm (Python code)
code format="python" from random import uniform

def dist(x,y): return (x[0] - y[0])**2 + (x[1] - y[1]) ** 2 N = 4 sigma = 0.2; sigma_sq = sigma**2 condition = False while condition == False: L = [(uniform(sigma, 1 - sigma), uniform(sigma, 1 - sigma))] for k in range(1,N): b = (uniform(sigma, 1 - sigma), uniform(sigma, 1 - sigma)) min_dist = min(dist(b, x) for x in L)       if min_dist < 4 * sigma_sq: condition = False break else: L.append(b) condition = True print L code

Local Metropolis algorithm (Python code)
code format="python" from random import uniform as ran, choice

L = [(0.25, 0.25), (0.75, 0.25), (0.25, 0.75), (0.75, 0.75)] sigma = 0.20; sigma_sq = sigma**2 delta = 0.15 for iter in range(100000): a = choice(L) L.remove(a) b = (a[0] + ran(-delta, delta), a[1] + ran(-delta, delta)) min_dist = min((b[0] - x[0]) ** 2 + (b[1] - x[1]) ** 2 for x in L)   box_cond = min(b[0],b[1]) < sigma or max(b[0], b[1]) > 1 - sigma if box_cond or min_dist < 4 * sigma ** 2: L.append(a) else: L.append(b) print L code

Variant of the Markov-chain Monte Carlo algorithm
In the Markov-chain Monte Carlo algorithm, for two disks, we can restrict ourselves to always move in one direction. This is explained in the following picture.

=Equivalence between Monte Carlo algorithm and molecular dynamics= The Boltzmann distribution describes physical reality, as we all know. But for any decades it was unclear whether this could actually be proven, or was more like an axiom of statistical mechanics. The mentioned mathematical papers by Sinai (1970) and by Simanyi (2003) have proven the exact validity of equipartition for finite hard-disk systems. =References=