Hi, I have been writing a program that simulates the behaviour of a system of hard disks using the monte carlo method. I ran the program upto around 2 million steps and it's really slow(takes around 3 1/2 to 4 hours).

I first wrote it in vpython as I'm a complete beginner and I needed to see where I was going wrong!

However I'm trying to move away from the visual side now as that should help speed things up. I tested psyco over 5,000 steps and found that it was actually 4 seconds slower than without it!

I was just wondering if anyone could give me any pointers as to things that could obviously be optimised.

Here's the full code, Ignore the PairCorrelationFunction function as I didn't write that part.

```
from visual import * #
from math import * #IMPORT LIBRARIES
from random import * #
from numpy import array
import time
start = time.time()
def PairCorrelationFunction_2D(x,y,S,rMax,dr):
"""Compute the two-dimensional pair correlation function, also known
as the radial distribution function, for a set of circular particles
contained in a square region of a plane. This simple function finds
reference particles such that a circle of radius rMax drawn around the
particle will fit entirely within the square, eliminating the need to
compensate for edge effects. If no such particles exist, an error is
returned. Try a smaller rMax...or write some code to handle edge effects! ;)
Arguments:
x an array of x positions of centers of particles
y an array of y positions of centers of particles
S length of each side of the square region of the plane
rMax outer diameter of largest annulus
dr increment for increasing radius of annulus
Returns a tuple: (g, radii, interior_x, interior_y)
g(r) a numpy array containing the correlation function g(r)
radii a numpy array containing the radii of the
annuli used to compute g(r)
interior_x x coordinates of reference particles
interior_y y coordinates of reference particles
"""
from numpy import zeros, sqrt, where, pi, average, arange, histogram
# Number of particles in ring/area of ring/number of reference particles/number density
# area of ring = pi*(r_outer**2 - r_inner**2)
# Find particles which are close enough to the box center that a circle of radius
# rMax will not cross any edge of the box
bools1 = x>1.1*rMax
bools2 = x<(S-1.1*rMax)
bools3 = y>rMax*1.1
bools4 = y<(S-rMax*1.1)
interior_indices, = where(bools1*bools2*bools3*bools4)
num_interior_particles = len(interior_indices)
if num_interior_particles < 1:
raise RuntimeError ("No particles found for which a circle of radius rMax\
will lie entirely within a square of side length S. Decrease rMax\
or increase the size of the square.")
edges = arange(0., rMax+1.1*dr, dr)
num_increments = len(edges)-1
g = zeros([num_interior_particles, num_increments])
radii = zeros(num_increments)
numberDensity = len(x)/S**2
# Compute pairwise correlation for each interior particle
for p in range(num_interior_particles):
index = interior_indices[p]
d = sqrt((x[index]-x)**2 + (y[index]-y)**2)
d[index] = 2*rMax
(result,bins) = histogram(d, bins=edges, normed=False, new=True)
g[p,:] = result/numberDensity
# Average g(r) for all interior particles and compute radii
g_average = zeros(num_increments)
for i in range(num_increments):
radii[i] = (edges[i] + edges[i+1])/2.
rOuter = edges[i+1]
rInner = edges[i]
g_average[i] = average(g[:,i])/(pi*(rOuter**2 - rInner**2))
return (g_average, radii, x[interior_indices], y[interior_indices])
#my code starts here
nlat=20 #HOW MANY DISKS ON THE SIDES OF THE BOX
r = 0.5 #DISK RADIUS
d = 2.0*r #DIAMETER
s = 1.22 #GAPS, ALTER TO CHANGE DENSITY
S = s - 1
L = (nlat*d)+(nlat-1)*S #LENGTH OF SIDES OF BOX
rMax = 5.0
dr=0.08 # Radius for increasing the size of the bins
S = L
area = L**2
pi = 3.141592654
diskarea= pi*(r**2) * nlat**2
density = diskarea/area
cpd = 0.9069
cpdf = (density / cpd)*100
print 'density = %f' % (density)
print 'which is %f percent of the closest packed density' % (cpdf)
scene.autoscale=0 #KEEP CAMERA STILL
scene = display(title='Monte Carlo Simulation of a Hard Disk Fluid',width=600, height=600,center=(L/2,L/2,0), background=(1,1,1))
#################################DRAW THE BOX#########################################
thk=0.3
s2 = L - thk
s3 = L + thk
wallR = box (pos=vector(L, L/2, 0), length=thk, height=s2, width=d, color = color.yellow)
wallL = box (pos=vector(0, L/2, 0), length=thk, height=s2, width=d, color = color.yellow)
wallB = box (pos=vector(L/2, 0, 0), length=s3, height=thk, width=d, color = color.yellow)
wallT = box (pos=vector(L/2, L, 0), length=s3, height=thk, width=d, color = color.yellow)
###################################################################################
#######CREATE SQUARE ARRAY OF PARTICLES (NLAT BY NLAT)############
alldisks = []
for x in range(0,nlat):
for y in range(0,nlat):
disk=sphere(pos=(x*s,y*s), radius = r, color=color.red)
alldisks.append(disk)
disk.id = (x,y)
####################################################
nstep=0.0 #SET FIRST STEP = 0
reject=0.0 #COUNT THE NUMBER OF TIMES A STEP IS REJECTED
maxpos=L #DEFINE BOUNDARY CONDITIONS
while (nstep<5000):####HOW MANY STEPS?####
nstep=nstep+1
p = choice(alldisks) #CHOOSE A PARTICLE,P, AT RANDOM FROM THE ARRAY
pnewx=p.x + 0.9*(random() -0.5) #GENERATE A STEP IN X CHANGE COEFFICIENT OF RANDOM TO GET 30% ACCEPTANCE
pnewy=p.y + 0.9*(random() -0.5) #GENERATE A STEP IN Y
######PERIODIC BOUNDARY CONDITIONS############
if pnewx > L: #IF X GOES BEYOND THE RIGHT WALL
pnewx=pnewx - L #MOVE X POSITION TO OTHER SIDE
if pnewx < 0: #IF X FALLS BELOW THE LEFT WALL
pnewx=pnewx + L #MOVE X POSITION TO OTHER SIDE
if pnewy > L: #IF Y IS ABOVE BOX
pnewy=pnewy - L #MOVE TO THE BOTTOM
if pnewy < 0: #IF Y BELOW
pnewy=pnewy + L #TO TOP
#################################################
rejected = False
for disk in alldisks:
distance=mag((vector(pnewx,pnewy))-disk.pos)#WORK OUT THE DISTANCE BETWEEN THE NEW POSITION FOR THE DISK AND ANY OTHER DISK
if (distance <= 2.0 * r) and p.id != disk.id: #if p overlaps any other disk
rejected = True
reject=reject+1 # COUNT HOW MANY TIMES THE STEP IS REJECTED
break
if not rejected:
p.x=pnewx #
p.y=pnewy #ACCEPT THE STEPS AS THE NEW POSITION FOR THE PARTICLE, P
if (nstep % 1800) == 0: #DUMP AFTER SO MANY STEPS
x_list = []
y_list = []
for disk in alldisks:
x_list.append(disk.x)
y_list.append(disk.y)
output = open('output''%s' % (nstep), 'w') #open file for outputting
(g, radii, newx, newy) = PairCorrelationFunction_2D(array(x_list), array(y_list), S, rMax, dr) #CALCULATE PAIR CORRELATION FUNCTION FOR EACH DUMP
for index in range(len(g)):
output.write( "%s %s\n" % (radii[index], g[index]))
output.close() #close file
continue
rejection = (reject/nstep) * 100 #percentage of steps rejected
acceptance = 100 - rejection
print 'acceptance rate was %f per cent' % (acceptance)
end = time.time()
elapsed= end - start
print "The simulation took", elapsed, "seconds to run"
```