Assignment 6: Off-lattice simulation of the hard-disk fluid in 2-d

Write a Monte Carlo code for a system of 225 disks confined within a box of size 40x40. Each disk has radius 1, which means they overlap if their center-to-center separation is less than 2.

Option 1: use periodic boundary conditions, and calculate both g(r) and mean square displacement vs time. For g(r), use concentric circles with spacing of 0.2 between them.

Option 2: use hard boundaries; calculate density as a function of x near the hard wall, using x-increments of size 0.2. Calculate also the mean square displacement as a function of time.

Note: the fact that rmax is fluctuating means that your mean square displacement may not be linear. However its continued growth is a sign that the system is in a fluid state and is not a glass.

Here is a meta-code to get you started. Details are included for calculating the pair correlation function.

Set up data structures

n=225 number of particles

nsteps=100000 number of monte carlo steps

define arrays x(n),y(n) these are particle positions

define arrays dx(n),dy(n) these are particle displacements from initial

positions

Data structures for the cell list

define array num(20,20) counts how many particles are in

cell number icell,jcell

define array clist(20,20,10) lists the particles contained in

cell number icell,jcell

define some real variables: dxcell=2.,dycell=2. these are cell sizes

define some integer variables: nxcells=20,nycells=20 number of cells in

x,y, directions

define some real variables:

hx=40,hy=40 these are my cell dimensions

rmax=1.0 this is my initial max step size for mc steps

theta this is a random angle for mc steps

rho0 this is the average number density of particles, = n/(hx*hy)

pi=4.*atan(1.0), twopi=pi*2.0

define variables for calculating g(r)

dr=0.2 spacing of concentric circles

if we go out to a distance of r=10 that gives us 50 circles

define real array g(50), this is g(r) for the current configuration

define real array gsave(50), a running average of g(r) over time

define integer array kount(50) counts how many particle pairs are

separated by a given distance

define integer kountg counts how many configs of g(r) are included in our running average

initialize the random number generator with a seed of your choice

iseed=5628901

call srand(iseed)

open a data file called r2ms.dat to store mean square displacement vs time

zero out dx(n),dy(n) arrays

zero out g(50),gsave(50), and set kountg=0

set up initial configuration by putting particles at the sites of a 15x15 lattice with spacing 2.66

k=0

do 1 ii=1,15

do 2 jj=1,15

k=k+1

x(k)=ii*2.66

y(k)=jj*2.66

2 continue

1 continue

Loop on istep=1 to nsteps

refresh the cell list

do 50 i=1,nxcells

do 51 j=1,nycells

num(i,j)=0

do 52 k=1,10

clist(i,j,k)=0

52 continue

51 continue

50 continue

do 60 i=1,n

icell=1+x(i)/dxcell

jcell=1+y(i)/dycell

num(icell,jcell)=num(icell,jcell)+1

clist(icell,jcell,num(icell,jcell))=i

60 continue

set naccept=0

Loop itry=1,n

pick a particle ip at random between 1 and n, save its x and y positions

figure out what cell it's in: icellold, jcellold

pick a random step size r between 0 and rmax

pick a random angle theta between 0 and twopi

move particle ip

if it went outside the periodic boundaries, add or subtract the boxsize to get it back in bounds

figure out what cell it's in now, call that icell,jcell

if not the same cell, add the particle to its new cell

num(icell,jcell)=num(icell,jcell)+1

clist(icell,jcell,num(icell,jcell))=ip

CHECK for overlaps:

step 1: loop on ic0=icell-1 to icell+1

ic=ic0

if ic<1 then ic=ic+nxcells

if ic>nxcells then ic=ic-nxcells

step 2: loop on jc0=jcell-1 to jcell+1

jc=jc0

if jc<1 then jc=jc+nycells

if jc>nycells then jc=jc-nycells

are there any particles in cell ic,jc?

if so, loop over them, k=1,num(ic,jc)

jp=clist(ic,jc,k)

if ip and jp aren't the same particle, calculate the distance squared between them using periodic bc's

if the distance squared is less than 4, there is an overlap, so put particle ip back where it was before, and go to next itry.

next k

next jc0

next ic0

made it all the way through with no overlaps, it's a success, hurray! naccept=naccept+1

add new step to the displacement arrays

dx(ip)=dx(ip)+r*cos(theta)

dy(ip)=dy(ip)+r*sin(theta)

next itry

check acceptance rate and adjust rmax as needed

lower rate threshold is 15%, upper rate threshold is 25%

if istep=100,200,300.... then let's print some data and make some

measurements

open a data file disk1.dat, dump x(i),y(i) for the current config

find the mean square displacement from your arrays dx(n),dy(n),

and write it into the data file r2ms.dat

calculate the pair correlation function in the present

configuration and add it to the running average. Here is my fortran code for this part. First, zero out the kount(50) and g(50) arrays

do 499 ir=1,50

kount(ir)=0

g(ir)=0.d0

499 continue

Now, loop on all pairs, calculate distance between, convert it to an integer ir, and increment kount(ir) by one.

do 500 i=1,n-1

do 501 j=i+1,n

ddx=dabs(x(i)-x(j))

ddx=min(ddx,hx-ddx)

ddy=dabs(y(i)-y(j))

ddy=min(ddy,hy-ddy)

r=dsqrt(ddx*ddx+ddy*ddy)

ir=1+int(r/dr)

kount(ir)=kount(ir)+1

501 continue

500 continue

Now normalize g(ir) by the area between the outer and inner circles. There are several other factors which we can talk about in class.

Accumulate a running sum of g(50) in the array gsave(50). Increment

kountg by one for every configthat goes into the running average.

g.dat is the g(r) for the current config only

gsave.dat is the average g(r) over the whole run

open(unit=16,file='g.dat',status='unknown')

open(unit=17,file='gsave.dat',status='unknown')

kountg=kountg+1

do 510 ir=1,50

r1=(ir-1)*dr

r2=ir*dr

area=pi*(r2*r2-r1*r1)

g(ir)=float(kount(ir))/area

g(ir)=g(ir)/rho0

g(ir)=g(ir)*2.d0

g(ir)=g(ir)/float(n)

write(16,*)r2,g(ir)

gsave(ir)=gsave(ir)+g(ir)

write(17,*)r2,gsave(ir)/float(kountg)

510 continue

close(unit=16)

close(unit=17)

next istep

close data files

all done!