AND I CHOSE THE TOP 3 STATE POINTS FOR THEM: describing an equation of state for
ID: 1770049 • Letter: A
Question
AND I CHOSE THE TOP 3 STATE POINTS FOR THEM:
describing an equation of state for the Lennard-Jones fluid. Choose three state points used in the paper and simulate these for argon with the Monte Carlo code provided, commenting on the accuracy and any deviations between your simulations and the published results. Clearly show all calculations for the input parameters at your chosen state points and relate how many moves you used for the equilibration and production stages. The Lennard Jones potential for argon has the parameters = 3.405 A and elk 119.8K
Explanation / Answer
This is the Monte Carlo python code I have written to produce these results :
First, run the build.py program to build your system with defined density,
#########**** _______________Supriyo Naskar, Department of Physics___________________________ ********###############
############************** Code for Montecarlo *********************#############
############************* Ensemble = NVT or NPT ********************#############
########*** System is Interacting via 12-6 Lennard-Jones Potential *******#######
#_______________Import some library__________#
from __future__ import division
import os
import numpy as np
import math as ma
import random
from pylab import *
from visual import sphere,rate, display,box,curve
#from ljenergy import TotalPairEnergy
#from metropolis import ConfProb
#__________________Inputs___________________#
N = 125 # Number of atoms for the simulation
density = 0.5 # Number of atoms/leng_volume for the simulation
epsilon = 1.0 # Lennard-Jones epsilon
sigma = 1.0 # Lennard-Jones sigma
mass = 1.0 # Mass of each particle
center = 0.0, 0.0, 0.0 # Boundary Conditions -- using hard cubic leng so need a center and leng length
leng = (N/density)**(1/3) # leng dimension
temp = 1.0 # Temperature
press = 1.0 # Pressure
kb=1.0 # Boltzmann Constant
rt = kb*temp/epsilon # if using real units this is the corresponding reduced temperature
f = open('build.txt', 'wb')
'''#__________________Visualization___________________#
def goto(linenum):
global line
line = linenum
def add(i,j,k,lengt,rad):
aa=lengt/2
s = sphere(radius=rad)
sphere(radius=rad).pos = i,j,k #i-50,j-50
cur=curve(pos=[(aa,-aa,-aa),(aa,aa,-aa),(aa,aa,aa),(aa,-aa,aa),(aa,-aa,-aa)])
cur.color=2,0,0
cur=curve(pos=[(-aa,-aa,-aa),(-aa,aa,-aa),(-aa,aa,aa),(-aa,-aa,aa),(-aa,-aa,-aa)])
cur.color=2,0,0
cur=curve(pos=[(aa,aa,aa),(-aa,aa,aa),(-aa,-aa,aa),(aa,-aa,aa),(aa,aa,aa)])
cur.color=2,0,0
cur=curve(pos=[(aa,aa,-aa),(-aa,aa,-aa),(-aa,-aa,-aa),(aa,-aa,-aa),(aa,aa,-aa)])
cur.color=2,0,0'''
#__________________Pair Energy between two atoms using LJ12-6___________________#
def LJPairEnergy(i,j):
epsilon = 1.0 # Lennard-Jones epsilon
sigma = 1.0 # Lennard-Jones sigma
eps = epsilon**2
sig = sigma**2
rcs = 2.5*sigma
dx = atom[i][0] - atom[j][0]
dy = atom[i][1] - atom[j][1]
dz = atom[i][2] - atom[j][2]
dx=dx-leng*round(dx/leng)
dy=dy-leng*round(dy/leng)
dz=dz-leng*round(dz/leng)
rs = dx**2.0 + dy**2.0 + dz**2.0
if rs<rcs:
E = 4.0*eps*( (sig/rs)**6.0 - (sig/rs)**3.0 ) - (4.0*eps*( (sig/rcs)**6.0 - (sig/rcs)**3.0 ))
else:
E = 0.0
return E
#Computes the total pair energy
def TotalPairEnergy(N):
Et = 0.0
for i in xrange(0, N-1):
for j in xrange(i+1, N):
Et += LJPairEnergy(i,j)
return Et
#computes the probability or probability ratio using Boltzmann weight
def ConfProb(dE, Beta):
prob = np.exp(-Beta*dE)
return prob
#__________________Assigning Initial Coordinates___________________#
###_____________________ Initial coordinates_______________________________###
#####-------- Assign randomly inside a cube ( overlap checking done)--###
####------------------ length**3 = V -----------------------------------###
atom=np.zeros((N,N))
for index in xrange(0, N):
atom[index][0] = center[0]+ leng * (random()-0.5)
atom[index][1] = center[1]+ leng * (random()-0.5)
atom[index][2] = center[2]+ leng * (random()-0.5)
Energy1 = TotalPairEnergy(N)
print "Assigning Initial Coordinates..."
overlap=np.zeros(N)
i_o=0
for index in range(0, N-1):
px=atom[index][0]
py=atom[index][1]
pz=atom[index][2]
for index1 in range (index+1,N):
if (abs(np.sqrt((px-atom[index1][0])**2 + (py-atom[index1][1])**2 +(pz-atom[index1][2])**2)) <= sigma):
overlap[index]=1
i_o=i_o+1
for index in range(0, N):
if (overlap[index]==1):
atom[index][0] = center[0]+ leng * (random()-0.5)
atom[index][1] = center[1]+ leng * (random()-0.5)
atom[index][2] = center[2]+ leng * (random()-0.5)
index1=0
while index1 < N :
if (index1!=index):
if (abs(np.sqrt((px-atom[index1][0])**2 + (py-atom[index1][1])**2 +(pz-atom[index1][2])**2)) <= sigma):
index1=index1
goto(127)
else :
index1=index1+1
elif(index1==index):
index1=index1+1
print index
Energy = TotalPairEnergy(N)
#for index in range(0,N):
# f.write("%f %f %f " %(atom[index][0],atom[index][1], atom[index][2]))
#f.close()
print "Press Esc to start Monte-Carlo simulation"
Then run this program to calculate your pressure and internal energy of the system,
#########**** Supriyo Naskar, Department of Physics********###############
#########**** __________________________________________ ********###############
############************** Code for Montecarlo *********************#############
############************* Ensemble = NVT or NPT ********************#############
########*** System is Interacting via 12-6 Lennard-Jones Potential *******#######
#_______________Import some library__________#
from __future__ import division
import os
import numpy as np
import math as ma
import random
from pylab import *
#from visual import sphere,rate, display,box,curve
#from ljenergy import TotalPairEnergy
#from metropolis import ConfProb
f = open('eq.txt', 'wb')
f1 = open('pro.txt', 'wb')
#__________________Inputs___________________#
N = 100 # Number of atoms for the simulation
density = 1.1 # Number of atoms/leng_volume for the simulation
epsilon = 1.0 # Lennard-Jones epsilon
sigma = 1.0 # Lennard-Jones sigma
mass = 1.0 # Mass of each particle
mc_it = 20000 # Number of trial moves
eq_it = 20000 # Number of Monte carlo iteration for equilibration
c_it = 10 # frequency for collect data after equilibration
s_it = 10 # Step interval to output running averages to screen
center = 0.0, 0.0, 0.0 # Boundary Conditions -- using hard cubic leng so need a center and leng length
leng = (N/density)**(1/3) # leng dimension
temp = 1.0 # Temperature
press = 1.0 # Pressure
kb=1.0 # Boltzmann Constant
rt = kb*temp/epsilon # if using real units this is the corresponding reduced temperature
#__________________Pair Energy and Force between two atoms using LJ12-6___________________#
def LJPairEnergy(i,j):
epsilon = 1.0 # Lennard-Jones epsilon
sigma = 1.0 # Lennard-Jones sigma
eps = epsilon**2
sig = sigma**2
rcs = 2.5*sigma
dx = atom[i][0] - atom[j][0]
dy = atom[i][1] - atom[j][1]
dz = atom[i][2] - atom[j][2]
dx=dx-leng*round(dx/leng)
dy=dy-leng*round(dy/leng)
dz=dz-leng*round(dz/leng)
rs = dx**2.0 + dy**2.0 + dz**2.0
if rs<rcs:
E = 4.0*eps*( (sig/rs)**6.0 - (sig/rs)**3.0 ) - (4.0*eps*( (sig/rcs)**6.0 - (sig/rcs)**3.0 ))
F = 48.0*eps*( (sig/rs)**6.0 - 0.5*(sig/rs)**3.0 )
else:
E = 0.0
F = 0.0
return E,F
#Computes the total pair energy and Pressure
def TotalPairEnergy(N):
Et = 0.0
pres1 = 0.0
for i in xrange(0, N-1):
for j in xrange(i+1, N):
Et += LJPairEnergy(i,j)[0]
pres1 += LJPairEnergy(i,j)[1]/3
pres = density/(kb*temp) + pres1/(leng**3)
return Et,pres
#computes the probability or probability ratio using Boltzmann weight
def ConfProb(dE, Beta):
prob = np.exp(-Beta*dE)
return prob
pos = np.loadtxt("build.txt",float)
atom_old=np.zeros((N,3))
atom=np.zeros((N,3))
#Assigning the positions
for index in xrange(0, N):
atom[index][0] = pos[index][0]
atom[index][1] = pos[index][1]
atom[index][2] = pos[index][2]
Energy = TotalPairEnergy(N)[0]
Pres = TotalPairEnergy(N)[1]
print "Energy of initial configuration : ", Energy
print "Pressure of initial configuration : ", Pres
def CalcCOM(N):
M = 0.0
sumx = sumy = sumz = 0.0
for i in xrange(0, N):
m = mass
sumx += atom[i][0] * m
sumy += atom[i][1] * m
sumz += atom[i][2] * m
M += m
Cx = sumx/M
Cy = sumy/M
Cz = sumz/M
return Cx,Cy,Cz
#print CalcCOM(N)
def RadGyr(N):
sumgyr = 0.0
M = 0.0
comx = CalcCOM(N)[0]
comy = CalcCOM(N)[1]
comz = CalcCOM(N)[2]
for i in xrange(0, N):
M += mass
mc = mass
rx = atom[i][0]
ry = atom[i][1]
rz = atom[i][2]
sgx = (rx - comx)**2.0
sgy = (ry - comy)**2.0
sgz = (rz - comz)**2.0
sumgyr += mc*(sgx+sgy+sgz)
Rsq = sumgyr / M
R = ma.sqrt(Rsq)
return R
def energy_RMSD(sc_it):
E=np.zeros(sc_it)
for j in range (0,sc_it):
E[j]=E_j
e=np.average(E)
e2=np.average(E*E)
c_v=(e2-e**2)/(kb*temp**2)
rms_eng=np.sqrt(e2)
return rms_eng, c_v
#print RadGyr(N)
#Periodic Boundary Condition
def PBC(index):
dx = atom[index][0] - 0
dy = atom[index][1] - 0
dz = atom[index][2] - 0
if atom[index][0] > 0 and dx > leng/2:
atom[index][0] = leng - atom[index][0]
elif atom[index][0] < 0 and dx < leng/2:
atom[index][0] = leng + atom[index][0]
else:
atom[index][0] = atom[index][0]
if atom[index][1] > 0 and dy > leng/2:
atom[index][1] = leng - atom[index][1]
elif atom[index][1] < 0 and dy < leng/2:
atom[index][1] = leng + atom[index][1]
else:
atom[index][1] = atom[index][1]
if atom[index][2] > 0 and dz > leng/2:
atom[index][2] = leng - atom[index][2]
elif atom[index][2] < 0 and dz < leng/2:
atom[index][2] = leng + atom[index][2]
return atom[index][0],atom[index][1],atom[index][2]
'''ball=np.zeros(N)
for index in range(0,N):
ball_index = sphere(pos=(atom[index][0],atom[index][1],atom[index][2]), radius=0.5)
add(atom[index][0],atom[index][1],atom[index][2],leng,0.5)'''
################### Start the Equilibration loop ##############################
count_eq=0
step=0.5
tot_count=0
for i in range(0, eq_it):
j=0
count_eq=0
for j1 in range (0,s_it):
#rate(0.00011)
Beta=1/kb/temp
#save energy before trial move
E_old_j = TotalPairEnergy(N)[0]
#save Pressure before trial move
P_old_j = TotalPairEnergy(N)[1]
#save COM before trial move
COM_old_j = CalcCOM(N)
#save radius of Gyration before trial move
Rad_old_j = RadGyr(N)
# Choose randomly a particle
index=randint(0,N)
#ball = sphere(pos=(atom[index][0],atom[index][1],atom[index][2]), radius=0.5)
#remove(atom[index][0],atom[index][1],atom[index][2],leng,0.5)
#save the old co-ordinate
atom_old[index][0]=atom[index][0]
atom_old[index][1]=atom[index][1]
atom_old[index][2]=atom[index][2]
#energy_old=TotalPairEnergy(N)[0]
# Give that particle a random displacement
p=step*(random()-0.5)
p1=p,p,p
atom[index][0] = p + atom_old[index][0]
atom[index][1] = p + atom_old[index][1]
atom[index][2] = p + atom_old[index][2]
# Periodic Boundary Condition Checking
atom[index][0] = PBC(index)[0]
atom[index][1] = PBC(index)[1]
atom[index][2] = PBC(index)[2]
#add(atom[index][0],atom[index][1],atom[index][2],leng,0.5)
#new energy
E_new_j = TotalPairEnergy(N)[0]
P_new_j = TotalPairEnergy(N)[1]
dE=E_new_j-E_old_j
prob=ConfProb(dE, Beta)
ranprob = random()
#New ROG
Rad_new_j = RadGyr(N)
#New COM
COM_new_j = CalcCOM(N)
if dE>0 :
if (random()<prob):
count_eq=count_eq+1
tot_count=tot_count+1
atom[index][0]=atom[index][0]
atom[index][1]=atom[index][1]
atom[index][2]=atom[index][2]
E_j=E_new_j
R_j=Rad_new_j
P_j=P_new_j
C_j=COM_new_j
else:
atom[index][0]=atom_old[index][0]
atom[index][1]=atom_old[index][1]
atom[index][2]=atom_old[index][2]
E_j=E_old_j
R_j=Rad_old_j
P_j=P_old_j
C_j=COM_old_j
else:
count_eq=count_eq+1
tot_count=tot_count+1
atom[index][0]=atom[index][0]
atom[index][1]=atom[index][1]
atom[index][2]=atom[index][2]
E_j=E_new_j
R_j=Rad_new_j
P_j=P_new_j
C_j=COM_new_j
j=j+1
if (count_eq/(j))<.45 and step>.11 :
step=step-.1
elif (count_eq/(j))>.55 and step>0 :
step=step+.1
elif (count_eq/(j))<.45 and step <.11 :
step=step/2
else:
step=step
print "===================================================================================================="
print " "
print "Equilibration is going on "
print "Statistics are given below"
print "Statistics after every =",j,"iteration and total no of iteration =",(i+1)*(j)
print "Aceptence Ratio",count_eq/(j),step
print "Energy = ",E_j
#print "rms energy =",energy_RMSD(s_it)[0]/N
print "Center of Mass =",C_j
print "Radius of gyration =", R_j
print "Pressure = ",P_j
print " "
f.write("%i %f %f %f %f " %(((i+1)*(j)),(E_j),(P_j),(energy_RMSD(s_it)[0]/N),(R_j)))
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
print "===================================================================================================="
################### Start the Production loop ##############################
count_pro=0
for i in range(0, mc_it):
j=0
count_pro=0
for j1 in range (0,c_it):
#rate(0.00011)
Beta=1/kb/temp
#save energy before trial move
E_old_j = TotalPairEnergy(N)[0]
#save Pressure before trial move
P_old_j = TotalPairEnergy(N)[1]
#save COM before trial move
COM_old_j = CalcCOM(N)
#save radius of Gyration before trial move
Rad_old_j = RadGyr(N)
# Choose randomly a particle
index=randint(0,N)
#ball = sphere(pos=(atom[index][0],atom[index][1],atom[index][2]), radius=0.5)
#remove(atom[index][0],atom[index][1],atom[index][2],leng,0.5)
#save the old co-ordinate
atom_old[index][0]=atom[index][0]
atom_old[index][1]=atom[index][1]
atom_old[index][2]=atom[index][2]
#energy_old=TotalPairEnergy(N)[0]
# Give that particle a random displacement
p=step*(random()-0.5)
p1=p,p,p
atom[index][0] = p + atom_old[index][0]
atom[index][1] = p + atom_old[index][1]
atom[index][2] = p + atom_old[index][2]
# Periodic Boundary Condition Checking
atom[index][0] = PBC(index)[0]
atom[index][1] = PBC(index)[1]
atom[index][2] = PBC(index)[2]
#add(atom[index][0],atom[index][1],atom[index][2],leng,0.5)
#new energy
E_new_j = TotalPairEnergy(N)[0]
dE=E_new_j-E_old_j
prob=ConfProb(dE, Beta)
ranprob = random()
#new pressure
P_new_j = TotalPairEnergy(N)[1]
#New ROG
Rad_new_j = RadGyr(N)
#New COM
COM_new_j = CalcCOM(N)
if dE>0 :
if (random()<prob):
count_pro=count_pro+1
tot_count=tot_count+1
atom[index][0]=atom[index][0]
atom[index][1]=atom[index][1]
atom[index][2]=atom[index][2]
E_j=E_new_j
R_j=Rad_new_j
P_j=P_new_j
C_j=COM_new_j
else:
atom[index][0]=atom_old[index][0]
atom[index][1]=atom_old[index][1]
atom[index][2]=atom_old[index][2]
E_j=E_old_j
R_j=Rad_old_j
P_j=P_old_j
C_j=COM_old_j
else:
count_pro=count_pro+1
tot_count=tot_count+1
atom[index][0]=atom[index][0]
atom[index][1]=atom[index][1]
atom[index][2]=atom[index][2]
E_j=E_new_j
R_j=Rad_new_j
P_j=P_old_j
C_j=COM_new_j
j=j+1
if (count_pro/(j))<.45 and step > 0.11:
step=step-.1
elif (count_pro/(j))>.55 and step>0:
step=step+.1
elif (count_pro/(j))<.45 and step < 0.11:
step=step/2
else:
step=step
print "===================================================================================================="
print " "
print "Production is going on "
print "Statistics are given below"
print "Statistics after every =",j,"iteration and total no of iteration =",(i+1)*(j)+(s_it*eq_it)
print "Aceptence Ratio",count_pro/(j),", Adjusting Step length=",step
print "Energy = ",E_j
print "Center of Mass =",C_j
print "Radius of gyration =", R_j
print "Pressure = ",P_j
print " "
print "===================================================================================================="
f1.write("%i %i %f %f %f %f " %(((i+1)*(j)),((i+1)*(j)+(s_it*eq_it)),(E_j),(P_j),(energy_RMSD(s_it)[0]/N),(R_j)))
print "Final Aceptence Ratio =",tot_count/(eq_it*s_it+c_it*mc_it)
#print "rms energy =",energy_RMSD((eq_it*s_it+c_it*mc_it))[0]/N,"Specific Heat =",energy_RMSD((eq_it*s_it+c_it*mc_it))[1]
Related Questions
drjack9650@gmail.com
Navigate
Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.