# Demo entry 6688228

monte

Submitted by anonymous on Jan 02, 2018 at 08:59
Language: Python 3. Code size: 5.2 kB.

```import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt

def RandomVector(d):
return 2*rnd.random(d)-1

def IsotropicVector(d): #for square
while True:
u = RandomVector(d)
if(-1<u[0]<1 and -1<u[1]<1):
break
return u

def AllocateAtoms(n,d):
return np.empty([n,d])

def InitPositions(r,n,d,boxsize):
for i in range(n):
r[i,:] = IsotropicVector(d)*boxsize

def InsideSquare(p,d,boxsize):
if -boxsize<p[0]<boxsize and -boxsize<p[1]<boxsize:
return True
else:
return False

def AllInsideSquare(r,n,d,boxsize):
inside = True
for i in range(n):
inside = inside and InsideSquare(r[i,:],d,boxsize)
return inside

def GenerateMove(r,n,d,stepsize):
imov = rnd.randint(n)
oldpos = r[imov,:]
u = IsotropicVector(d)
newpos = oldpos + stepsize*u
return imov, oldpos, newpos

def AcceptMove(imov,oldpos,newpos,r,n,d):
r[imov,:] = newpos

def RejectMove(imov,oldpos,newpos,r,n,d):
r[imov,:] = oldpos

def Potential(r2): #lenard jones potential
return 4*((1.0/r2**6)-(1.0/r2**3))

def ComputeEnergy(r,n,d):
E=0.0
for i in range(n-1):
for j in range(i+1,n):
r2 = np.sum((r[i,:]-r[j,:])**2)
E += Potential(r2)
return E

def ComputeDiffEnergy(imov,oldpos,newpos,r,n,d):
DE = 0.0
for i in range(n):
if (not(i == imov)):
r2 = np.sum((r[i,:]-newpos)**2)
DE += Potential(r2)
r2 = np.sum((r[i,:]-oldpos)**2)
DE -= Potential(r2)
return DE

def MCstep(r,n,d,Energy,nacc,Temperature,boxsize,stepsize,Shape=InsideSquare):
imov, oldpos, newpos = GenerateMove(r,n,d,stepsize)
if (not Shape(newpos,d,boxsize)):
RejectMove(imov,oldpos,newpos,r,n,d)
else:
DE = ComputeDiffEnergy(imov,oldpos,newpos,r,n,d)
if DE <= 0.0:
AcceptMove(imov,oldpos,newpos,r,n,d)
Energy += DE
nacc += 1
else:
if rnd.random()<np.exp(-DE/Temperature):
AcceptMove(imov,oldpos,newpos,r,n,d)
Energy += DE
nacc += 1
else:
RejectMove(imov,oldpos,newpos,r,n,d)
return Energy, nacc

def PlotPositions(r,n,d,boxsize,Shape = InsideSquare):
f,((ax)) = plt.subplots(1)
ax.scatter(r[:,0],r[:,1])
ax.plot([ boxsize, boxsize],[ boxsize,-boxsize],'r')
ax.plot([ boxsize,-boxsize],[-boxsize,-boxsize],'r')
ax.plot([-boxsize,-boxsize],[ boxsize,-boxsize],'r')
ax.plot([ boxsize,-boxsize],[ boxsize, boxsize],'r')
plt.show()

def analyse(u,boxsize):
def circ(x):
return np.sqrt(circ.b**2.-x**2.)
circ.b = 0
nBins = u.shape[0]//1
binEdges = np.linspace(0., np.sqrt(2.)*boxsize, nBins+1)
counts = np.zeros(nBins, dtype=int)
for j in range(u.shape[0]):
for i in range(nBins):
if np.linalg.norm(u[j, :])<binEdges[i+1]:
counts[i] += 1
break
areas = np.empty(nBins)
areaError = np.zeros(nBins)
for i in range(nBins):
if binEdges[i+1]>boxsize:
alpha  = np.sqrt(max(0, binEdges[i]**2.-boxsize**2.))
beta   = np.sqrt(binEdges[i+1]**2.-boxsize**2.)
gamma  = min(boxsize,binEdges[i])
circ.b = binEdges[i+1]
circ.b = binEdges[i]
circ.b = binEdges[i+1]
areas[i] = 4.*(int1[0]-int2[0]-int3[0]+boxsize*(beta-alpha))
areaError[i] = 4.*np.sqrt(int1[1]**2+int2[1]**2+int3[1]**2)
else:
areas[i] = np.pi*(binEdges[i+1]**2.-binEdges[i]**2)
f, ((a,b),
(c,d)) = plt.subplots(2,2)
x = binEdges[:-1]+np.diff(binEdges)[0]*0.5
a.plot(x,counts)
a.set_title('Counts')
b.plot(x, areas)
b.set_title('Area')
c.plot(x, counts/areas)
c.set_title('Density')
d.plot(x, areaError)
d.set_title('Error in area')
plt.show()

d=2
boxsize = 2.5
Temperature = 10.0
stepsize = np.float(boxsize/50)
n = 15
N = 100000
Nequil = int(N/10)
r = AllocateAtoms(n,d)
InitPositions(r,n,d,boxsize)

# Equilibration:
nacc = 0
Energy = ComputeEnergy(r,n,d)
for i in range(Nequil):
Energy, nacc = MCstep(r,n,d,Energy,nacc,Temperature,boxsize,stepsize)

# Production run
nacc = 0
sum_Energy = 0.0
sum_Energy2 = 0.0
Energy = ComputeEnergy(r,n,d)

for i in range(N):
Energy, nacc = MCstep(r,n,d,Energy,nacc,Temperature,boxsize,stepsize)
sum_Energy  += Energy
sum_Energy2 += Energy**2

PlotPositions(r,n,d,boxsize)
analyse(r,boxsize)

# Calculate averages
av_Energy  = sum_Energy/float(N)
av_Energy2 = sum_Energy2/float(N)
CV = (av_Energy2-av_Energy**2)/Temperature**2

# print results
print("acceptance ratio",nacc/float(N))
print("average energy",av_Energy)
print("heat capacity",CV)
```

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.