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

from scipy.integrate import quad 

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]
            int1   = quad(circ, alpha, boxsize)
            circ.b = binEdges[i]
            int2   = quad(circ, alpha, gamma)
            circ.b = binEdges[i+1]
            int3   = quad(circ, alpha, beta)
            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.

Delete this entry (admin only).