3.5. Appendix: Configurational sampling python program¶
Note that the mdmc environment needs to be enabled (or the package py3Dmol needs to be installed). In a local environment or on Google colab you can execute !pip install py3Dmol. On noto.epfl.ch this won’t work and you need to follow the instructions using the venv package.
# Necessary imports
import numpy as np
import py3Dmol
3.5.1. Helper Functions¶
Below are some helper functions to make the code run
def printXYZ(coords, filename, step, mode='vmd'):
""" Write Coordinates to trajectory file """
nPart=coords.shape[0]
if mode=='vmd':
with open(filename+'.xyz', 'a') as xyz:
xyz.write(f"{nPart}\n\n")
for i in range(nPart):
xyz.write(f"Ar {coords[i][0]:10.5f} {coords[i][1]:10.5f} {coords[i][2]:10.5f}\n")
if mode=='python':
with open(filename+'.'+str(step)+'.xyz', 'a') as xyz:
xyz.write(f"{nPart}\n\n")
for i in range(nPart):
xyz.write(f"Ar {coords[i][0]:10.5f} {coords[i][1]:10.5f} {coords[i][2]:10.5f}\n")
def placeParticlesOnGrid(nPart=100, density=0.85):
"""Initialize LJ Particles on grid"""
coords=np.zeros((nPart,3))
L=(nPart/density)**(1.0/3)
nCube=2
while (nCube**3 < nPart):
nCube+=1
index = np.array([0,0,0])
for part in range(nPart):
coords[part]=index+np.array([0.5,0.5,0.5])*(L/nCube)
index[0]+=1
if index[0]==nCube:
index[0]=0
index[1]+=1
if index[1]==nCube:
index[1]=0
index[2]+=1
return coords, L
def LJ_Energy(coords, L):
"""Calculate energy according to Lennard Jones Potential"""
energy=0
nPart=coords.shape[0]
for partA in range(nPart-1):
for partB in range(partA+1,nPart):
dr = coords[partA]-coords[partB]
dr = distPBC3D(dr, L)
dr2=np.sum(np.dot(dr,dr))
#Lennard-Jones potential:
# U(r) = 4*epsilon* [(sigma/r)^12 - (sigma/r)^6]
# Here, we set sigma = 1, epsilon = 1 (reduced distance and energy units). Therefore:
# U(r) = 4 * [(1/r)^12 - (1/r)^6]
# For efficiency, we will multiply by 4 only after summing
# up all the energies.
invDr6= 1.0/(dr2**3) # 1/r^6
energy = energy +(invDr6*(invDr6-1))
return energy*4
def distPBC3D(dr, L):
"""Calculate distance according to minimum image convention"""
hL=L/2.0
# Distance vector needs to be in [-hLx, hLX], [-hLy, hLy] [-hLz,hLz]
for dim in range(3):
if dr[dim]>hL:
dr[dim]-=L
elif dr[dim]<-hL:
dr[dim]+=L
return dr
def putParticlesInBox(vec, L):
"""Put coordinates back into primary periodic image"""
# Coord needs to be in [0, L], [0,L] [0,L]
for dim in range(3):
if vec[dim]>L:
vec[dim]-=L
elif vec[dim]<-L:
vec[dim]+=L
return vec
def deltaLJ(coords, trialPos, part, L):
"""Calculate Energy change of Move"""
deltaE=0
npart=coords.shape[0]
for otherPart in range(nPart):
if otherPart == part:
continue
# Particle Particle dist
drNew = coords[otherPart]-trialPos
drOld = coords[otherPart]-coords[part]
drNew = distPBC3D(drNew, L)
drOld = distPBC3D(drOld, L)
dr2_New = np.sum(np.dot(drNew, drNew))
dr2_Old = np.sum(np.dot(drOld, drOld))
# Lennard-Jones potential:
# U(r) = 4*epsilon* [(sigma/r)^12 - (sigma/r)^6]
# with sigma = 1, epsilon = 1 (reduced units).
# =>
# U(r) = 4 * [(1/r)^12 - (1/r)^6]
# We multiply by 4 only in the end
invDr6_New = 1.0/(dr2_New**3)
invDr6_Old = 1.0/(dr2_Old**3)
eNew = (invDr6_New*(invDr6_New-1))
eOld = (invDr6_Old*(invDr6_Old-1))
deltaE = deltaE + eNew-eOld
return deltaE*4
3.5.2. Main code to run MC in NVT ensemble¶
#Set configuration parameters
nPart = 100; # Number of particles
density = 0.4; # Density of particles
#Set simulation parameters
Temp = 2.0; # Simulation temperature
beta = 1.0/Temp; # Inverse temperature
maxDr = 0.1; # Maximal displacement
nSteps = 10000; #Total simulation time (in integration steps)
printFreq = 1000 #Printing frequency
# Set initial configuration
coords, L = placeParticlesOnGrid(nPart,density)#
# Calculate initial energy
energy = LJ_Energy(coords,L)
for step in range(nSteps):
if step%printFreq==0:
print(f"MC Step {step:2} Energy {energy:10.5f}")
printXYZ(coords, 'trajectory', step, mode='python') # mode VMD
for i in range(nPart):
rTrial = coords[i]+maxDr*(np.random.rand(1,3).squeeze()-0.5)
rTrial = putParticlesInBox(rTrial, L)
deltaE = deltaLJ(coords, rTrial, i, L)
if np.random.rand() < np.exp(-beta*deltaE):
coords[i]=rTrial
energy +=deltaE
MC Step 0 Energy -235.93251
MC Step 1000 Energy -393.70441
MC Step 2000 Energy -1056.16208
MC Step 3000 Energy -1152.87243
MC Step 4000 Energy -1309.31165
MC Step 5000 Energy -1413.22709
MC Step 6000 Energy -1550.79566
MC Step 7000 Energy -1561.81937
MC Step 8000 Energy -1868.14282
3.5.3. Visualization of results¶
Note that you can also use VMD for this.
If you use python the printing function needs to use the python mode, for vmd the vmd mode should be used.
printXYZ(coords, 'trajectory', step, mode='python') # to print single xyz files for py3Dmol
printXYZ(coords, 'trajectory', step, mode='vmd') # to trajectory xyz files for vmd
view = py3Dmol.view(data="",linked=False,viewergrid=(2,2))
with open('trajectory.0.xyz', 'r') as file:
xyz = file.read()
view.addModel(xyz, viewer=(0,0))
view.setStyle({'sphere':{'scale':0.2}}, viewer=(0,0))
view.zoomTo(viewer=(0,0))
with open('trajectory.1000.xyz', 'r') as file:
xyz = file.read()
view.addModel(xyz, viewer=(0,1))
view.setStyle({'sphere':{'scale':0.2}}, viewer=(0,1))
view.zoomTo(viewer=(0,1))
with open('trajectory.7000.xyz', 'r') as file:
xyz = file.read()
view.addModel(xyz, viewer=(1,0))
view.setStyle({'sphere':{'scale':0.2}}, viewer=(1,0))
view.zoomTo(viewer=(1,0))
with open('trajectory.9000.xyz', 'r') as file:
xyz = file.read()
view.addModel(xyz, viewer=(1,1))
view.setStyle({'sphere':{'scale':0.2}}, viewer=(1,1))
view.zoomTo(viewer=(1,1))
view.render()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
<py3Dmol.view at 0x7f3b72867ca0>
3.5.4. MC NPT Program¶
Modifcations to run in the NPT ensemble
#Initialize
#===================
#Set configuration parameters
nPart = 100 # Number of particles
density = 0.85 # Density of particles
#Set simulation parameters
Temp = 2.0 # Simulation temperature
beta = 1.0/Temp # Inverse temperature
maxDr = 0.1 # Maximal displacement
press = 1.0
maxDv = 0.01
nSteps = 1000; #Total simulation time (in integration steps)
printFreq = 100 #Printing frequency
sampleCounter = 0
sampleFreq = 100
# Set initial configuration
coords, L = placeParticlesOnGrid(nPart,density)
# Calculate initial energy
energy = LJ_Energy(coords,L)
avgL = 0
for step in range(nSteps):
# Choose whether to do a volume move or a displacement move
if (np.random.rand()*(nPart+1)+1 < nPart):
for i in range(nPart):
rTrial = coords[i]+maxDr*(np.random.rand(1,3).squeeze()-0.5)
rTrial = putParticlesInBox(rTrial, L)
deltaE = deltaLJ(coords, rTrial, i, L)
if np.random.rand() < np.exp(-beta*deltaE):
coords[i]=rTrial
energy +=deltaE
else:
oldV= L**3
lnvTrial = np.log(oldV)+(np.random.rand()-0.5)*maxDv
vTrial = np.exp(lnvTrial)
newL = vTrial**(1.0/3)
# Rescale coords
coordsTrial=coords*(newL/L)
eTrial = LJ_Energy(coordsTrial, newL)
weight = (eTrial - energy) + press*(vTrial-oldV) - (nPart+1)*Temp*np.log(vTrial/oldV)
if (np.random.rand()<np.exp(-beta*weight)):
coords = coordsTrial
energy = eTrial
L = newL
if step%printFreq==0:
print(f"MC Step {step:2} Energy {energy:10.5f}")
printXYZ(coords, 'npt_trajectory', step, mode='python')
if step%sampleFreq==0:
sampleCounter+=1
avgL+=L
print(f"Average Box Size {avgL/sampleCounter:3.3f}")
MC Step 0 Energy -107.17941
Average Box Size 4.900
MC Step 100 Energy -138.38271
Average Box Size 4.900
MC Step 200 Energy -165.94861
Average Box Size 4.900
MC Step 300 Energy -298.44862
Average Box Size 4.900
MC Step 400 Energy -305.89813
Average Box Size 4.900
MC Step 500 Energy -313.93907
Average Box Size 4.900
MC Step 600 Energy -314.25791
Average Box Size 4.900
MC Step 700 Energy -328.15041
Average Box Size 4.900
MC Step 800 Energy -421.19435
Average Box Size 4.900
MC Step 900 Energy -455.50053
Average Box Size 4.900