Hard+disks+and+the+liquid-solid+transition+in+two+dimensions 

toc
 * DM3 Hard disks and the liquid-solid transition in two dimensions **


 * Don't hesitate to ask questions and make remarks on this wiki page.**

=Introduction= We examine different algorithms to sample equilibrium configurations of hard disks in a box with or without periodic boundary conditions. We examine the physical properties of the system as a function of the density.

For a generic reference, see Chapter 2 of W. Krauth, //Statistical mechanics: Algorithms and Computations//.

=Hard disks: two Monte Carlo Algorithms= We start by defining a few programs using different algorithms to sample configurations of hard disks.

Direct sampling with periodic boundary conditions
Here is program direct_disks.py : code format="python" from random import random import math

def dist(x, y): d_x = abs(x[0] - y[0]) % 1.0 d_x = min(d_x, 1.0 - d_x) d_y = abs(x[1] - y[1]) % 1.0 d_y = min(d_y,1-d_y) return d_x**2 + d_y**2 N = 16 sigma = 0.1 sigma_sq = sigma**2 condition = False while condition == False: L = [(random, random)] for k in range(1, N): b = (random, random) min_dist = min(dist(b, x) for x in L)       if min_dist < 4.0 * sigma_sq: condition = False break else: L.append(b) condition = True print L code

Markov Chain. Version A: the hard box
Here is program markov_disks_box.py : 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(1000000): 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.0 - sigma if box_cond or min_dist < 4.0 * sigma ** 2: L.append(a) else: L.append(b) print L code

Markov Chain. Version B: the periodic boundary conditions
Here is program markov_disks_box_B.py : code format="python" from random import uniform as ran, choice def dist(x,y): d_x = abs(x[0] - y[0]) % 1.0 d_x = min(d_x, 1.0 - d_x) d_y = abs(x[1] - y[1]) % 1.0 d_y = min(d_y, 1.0 - d_y) return d_x ** 2 + d_y ** 2 L = [(0.25, 0.25), (0.75, 0.25), (0.25, 0.75), (0.75, 0.75)] sigma = 0.20 delta = 0.15 for iter in range(1000000): a = choice(L) L.remove(a) b = ((a[0] + ran(-delta, delta)) % 1.0, (a[1] + ran(-delta,delta)) % 1.0) min_dist = min(dist(b, x) for x in L)   if  min_dist < 4.0 * sigma ** 2: L.append(a) else: L.append(b) print L code


 * Answer the following questions:
 * 1) Discuss the three algorithms (Direct sampling, and the two Markov-chain algorithms).
 * 2) Justify precisely why one can forget about the velocities of the disks.
 * 3) Explain why (and in which limit) direct_disks.py and markov_disks_box_B.py give the same results.
 * 4) Explain why (and in which limit) these Monte Carlo algorithms give the same results as molecular dynamics.


 * Run the program Version A:
 * 1) Plot the histogram of the //x//-coordinate of the disks.
 * 2) Property: the probability density for hard spheres moving in a finite volume to visit a given point is uniform in space. Is this true?
 * 3) Do you think it is possible to test this property using the 2-dimensional histogram introduced in Tutorial 03?
 * 4) Propose a way to test numerically this property.

=The high density phase=

Crystals in two dimensions?
At moderate density //η// both algorithms above give configurations that are liquid-like. Alder and Wainwright discovered in 1962, using the molecular dynamics approach, that, surprisingly, a phase transition occurs at //η//≃0.7. At high density the configurations seem to crystalize on a hexagonal lattice. However, a rigorous theorem (Mermin and Wagner 1966) forbids positional long-range order for two-dimensional systems with short-range interactions, a class to which hard disks belong. This means that an infinitely large system cannot have endless patterns of disks nearly aligned, as shown in the figure. Nevertheless, in two dimensions, long-range order is possible for the orientation of links between neighbors. A detailed discussion of crystalline order in two dimensions would go beyond the scope of this exercise (see however, the following beautiful thesis by E. P. Bernard). Here we propose to "rediscover" the beautiful phase transition observed by Alder and Wainwright.



math 1 \times \sqrt{3} / 2 ~ math > with periodic boundary conditions.
 * Write a version C of markov_disks_box.py for //N// = //n// 2 = 64 hard disks in a box of size
 * Run Version C for densities //η//≃0.4 and //η//≃0.7.

Note that //η// is defined as math \eta= \frac{\text{total area of the disks}}{\text{ total area of the box}} math As initial condition for your simulation (which might be rather long: more than 10 minutes) prepare the disks on a hexagonal lattice (the program initial.py given below generates the correct initial condition and produces snapshots) :

First line math (0,0), (\Delta_x,0), (2\Delta_x,0) \ldots (L-\Delta_x,0) math Second line math (\tfrac{1}{2} \Delta_x,\Delta_y), (\tfrac{3}{2} \Delta_x,\Delta_y) \ldots (L - \tfrac{1}{2} \Delta_x, \Delta_y ) ~ math and so on, with

math \Delta_x = \frac{L}{n} \quad\text{and}\quad \Delta_y= \Delta_x \sqrt{3} / 2 ~ math


 * Analyze and discuss few snapshots of these two simulations.

Here is program initial.py :

code format="python" import math,random,pylab N=8; Ntot=N*N; eta=.4; sigma=math.sqrt(math.sqrt(3.)*eta/math.pi/128.) Dx=1./8. Dy=math.sqrt(3.)*Dx/2. positions=(k1+(k2%2)/2.)*Dx,k2*Dy] for k1 in range(N) for k2 in range(N)] pylab.axes for [x,y] in positions:  cir=pylab.Circle((x,y), radius=sigma,  fc='r')   pylab.gca.add_patch(cir) pylab.axis('scaled') pylab.show [[code

=Variant of Monte Carlo Algorithms: true or false?=

In the following Version B of the Markov chain we implemented the computation of the Pair Correlation Function. Here is program markov_disks_box_B.py : code format="python" from random import uniform as ran, choice import math, pylab def dist(x,y): d_x = abs(x[0]-y[0])%1 d_x = min(d_x,1-d_x) d_y = abs(x[1]-y[1])%1 d_y = min(d_y,1-d_y) return d_x**2 + d_y**2 L=[(0.25,0.25),(0.75,0.25),(0.25,0.75),(0.75,0.75)] sigma = 0.20 delta = 0.15 data = [] for iter in range(1000000): a = choice(L) L.remove(a) b = (a[0] + ran(-delta,delta),a[1] + ran(-delta,delta))%1 #bug corrected 01-OCT-2012, WK   min_dist = min(dist(b,x) for x in L)    if  min_dist < 4*sigma**2: L.append(a) else: L.append(b) a = choice(L) b = choice(L) if a != b: data.append( math.sqrt(dist(a,b))) pylab.hist(data,bins=40,normed=True,facecolor='green') pylab.show code Consider the following two variants: > b = (a[0] + ran(-delta,delta),a[1] + ran(-delta,delta))%1 > with > b = (a[0] + ran(0,delta),a[1] + ran(0,delta))%1 > Is detailed balance respected? Test with the pair correlation function if the result is correct
 * replace a = choice(L) with a=L[iter%N] . Is detailed balance respected? Test with the pair-correlation function if the result is correct
 * replace


 * Super Bonus**: the variant 2 does not satisfy the global balance condition. Can you show explicitly why? [Try the case of three particles first.]

=References=

W. Krauth Statistical mechanics: Algorithms and Computations (Oxford University Press, 2006). See also the webpage of the book. B. J. Alder and T. E. Wainwright, Phase Transition in Elastic Disks Phys. Rev. 127, 359 (1962). E. P. Bernard, PhD Thesis, 2011 (UPMC/ENS) Algorithms and applications of the Monte Carlo method: Two-dimensional melting and perfect sampling.

url}?f=print| [Print this page]