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

Don't hesitate to ask questions and make remarks onthis 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 :

fromrandomimportrandomimportmathdef 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 =Falsewhile condition ==False:
L =[(random(),random())]for k inrange(1, N):
b =(random(),random())
min_dist =min(dist(b, x)for x in L)if min_dist <4.0 * sigma_sq:
condition =Falsebreakelse:
L.append(b)
condition =Trueprint L

Markov Chain. Version A: the hard box

Here is program markov_disks_box.py :

fromrandomimport 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.15foriterinrange(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]) ** 2for x in L)
box_cond =min(b[0], b[1])< sigma ormax(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

Markov Chain. Version B: the periodic boundary conditions

Here is program markov_disks_box_B.py :

fromrandomimport 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.15foriterinrange(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

Answer the following questions:

Discuss the three algorithms (Direct sampling, and the two Markov-chain algorithms).

Justify precisely why one can forget about the velocities of the disks.

Explain why (and in which limit) direct_disks.py and markov_disks_box_B.py give the same results.

Explain why (and in which limit) these Monte Carlo algorithms give the same results as molecular dynamics.

Run the program Version A:

Plot the histogram of the x-coordinate of the disks.

Property: the probability density for hard spheres moving in a finite volume to visit a given point is uniform in space. Is this true?

Do you think it is possible to test this property using the 2-dimensional histogram introduced in Tutorial 02 (Sinai's histogram)?

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.

Initial condition for //η//≃0.7. The box has periodic boundary conditions.

Initial condition for //η//≃0.4. The box has periodic boundary conditions.

Write a version C of markov_disks_box.py for N = n 2 = 64 hard disks in a box of size

with periodic boundary conditions.

Run Version C for densities η≃0.4 and η≃0.7.

Note that η is defined as

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

Second line

and so on, with

Analyze and discuss few snapshots of these two simulations.

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 :

fromrandomimport uniform as ran, choice
importmath, 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 =[]foriterinrange(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()

Consider the following two variants:

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 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

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

## Table of Contents

DM3 Hard disks and the liquid-solid transition in two dimensionsDon'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 :## Markov Chain. Version A: the hard box

Here is program markov_disks_box.py :## Markov Chain. Version B: the periodic boundary conditions

Here is program markov_disks_box_B.py :x-coordinate of the disks.## 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.N=n2 = 64 hard disks in a box of sizeη≃0.4 andη≃0.7.Note that

ηis defined asAs 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

Second line

and so on, with

Here is program initial.py :

## 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 :

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

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.

[Print this page]