The Lennard-Jones Potential
Contents
3.3. The Lennard-Jones Potential#
In this exercise you will study the configuration of a collection of gaseous particles using the Metropolis Monte Carlo algorithm. The system includes \(N\) particles within a cubic box of volume \(V\) at a given temperature \(T\), in any configuration permitted by the Lennard-Jones potential:
This potential traditionally has an infinite range, however, the potential decays rapidly with separation distance and can be effectively ignored at large \(|r|\), resulting in a faster calculation. In practical applications it is customary to establish a cutoff \(r_c\) and disregard pairwise interactions separated beyond this radius. This truncation leads to a discontinuity in the pairwise potential energy function; large numbers of these events are likely to spoil energy conservation thus an improvement is to shift the potential such that the energy continuously approaches zero at \(r_c\):
This approach results in a potential that produces discontinuities in the first and higher order derivatives. To compensate, switching functions are often employed to smoothly and continuously taper the pair potential to zero between two cutoff limits.
Truncating pair interactions systematically removes a non-trivial contribution to the net potential energy and pressure. For interactions that are cut but not shifted, one can approximately add the interactions beyond \(r_c\) to the total energy and pressure, assuming the radial distribution function \(g(r > r_c) \approx 1\):
where
To sample configurational space using the Lennard-Jones potential, a randomly selected particle is first randomly translated to generate a new system configuration. Whether the new configuration is accepted depends on the acceptance probability discussed in the preceding section. This procedure repeats iteratively such that classical phase space is directly sampled and ensemble averages of physical properties become arithmetic averages over their sampled values.
3.3.1. Configurational Sampling Code#
3.3.1.1. Imports and Helper Functions#
Below are some necessary imports and useful helper functions to make the code run, just run the following cells once
# Necessary imports
import numpy as np
import MDAnalysis as md
import pytraj as pt
import sys
import os
sys.path.append("..")
import helpers
from helpers import *
helpers.set_style()
def remove_file(filename):
"Delete trajectory file with previous runs"
try:
os.remove(filename)
except:
pass
def printXYZ(coords, filename, step, mode='vmd'):
""" Write Coordinates to trajectory file """
nPart=coords.shape[0]
with open(filename+'.pdb', 'a') as xyz:
xyz.write('MODEL\n')
for i in range(nPart):
xyz.write(f"ATOM {i+1:3} Ar X 1 {coords[i][0]: 3.3f} {coords[i][1]: 3.3f} {coords[i][2]: 3.3f} 0.00 0.00 Ar\n")
xyz.write('ENDMDL\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.3.1.2. Main Code NVT Ensemble#
Here you can find the main code to run MC in the NVT ensemble
#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
nSteps = 1000; #Total simulation time (in integration steps)
printFreq = 100 #Printing frequency
# Remove trajectory file
remove_file('trajectory.pdb')
# 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)
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
3.3.1.3. Visualization of the Results#
Executing the code below you can visualize the MC trajectory for the particles
traj = md.Universe('trajectory.pdb')
show_trajectory(traj, only_spheres=True)
Exercise 8
Perform a simulation at \(T = 2.0\) and \(\rho \in [0.05, 0.5, 2, 10.0]\). What do you observe?
Exercise 9
The program produces a trajectory. Look at it and explain the behavior of the system over time for \(\rho\) 0.85
Exercise 10
Instead of performing a trial move in which only one particle is displaced, one can do a trial move in which more particles are
displaced simultaneously. You can find below a modified version of the code, where for each step nParticlesMove are displaced.
From the given code, explain how a trial move is now performed and what changes with respect to the previous case.
Hint: It can be helpful to perform different tests, e.g. using nParticlesMove = nPart or, for example, nParticlesMove = nPart/2. Describe the changes you see in the energies and the trajectories in this new version of the code.
#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
nSteps = 1000; #Total simulation time (in integration steps)
printFreq = 100 #Printing frequency
nParticlesMove = int(nPart/1.0) # Number of particles moved simultaneously
# Remove trajectory file
remove_file('trajectory.pdb')
# Set initial configuration
coords, L = placeParticlesOnGrid(nPart,density)#
for step in range(nSteps):
TrialCoords = coords.copy()
if step%printFreq==0:
print(f"MC Step {step:2} Energy {LJ_Energy(coords,L):10.5f}")
printXYZ(coords, 'trajectory', step)
for i in np.random.randint(0,nPart,nParticlesMove): # the particles to move are extracted randomly between all the particles
# np.random.randint(low, high, N) returns N random integers from low to high
rTrial = TrialCoords[i]+maxDr*(np.random.rand(1,3).squeeze()-0.5)
rTrial = putParticlesInBox(rTrial, L)
TrialCoords[i] = rTrial
deltaE = LJ_Energy(TrialCoords,L) - LJ_Energy(coords,L)
if (np.random.rand() < np.exp(-beta*deltaE)):
coords = TrialCoords
traj = md.Universe('trajectory.pdb')
show_trajectory(traj, only_spheres=True)
3.3.1.4. Main Code NPT Ensemble#
These are the modifcations to run in the NPT ensemble
Exercise 11 - Bonus
The code to sample the NPT ensemble is provided below. What needs to be changed in the code to sample from the isothermic-isobaric ensemble (NpT) instead of the microcanonical ensemble (NVT)?
#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.001 # Maximal displacement
press = 1.0
maxDv = 0.1 # Change in volume
nSteps = 1000; # Total simulation time (in integration steps)
printFreq = 100 # Printing frequency
sampleFreq = 100
# Remove trajectory file
remove_file('trajectory.pdb')
# Set initial configuration
coords, L = placeParticlesOnGrid(nPart,density)
# Calculate initial energy
energy = LJ_Energy(coords,L)
avgL = 0
sampleCounter = 0 # Set counter to 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):
# Displacement move
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:
# Volume move
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}")