# I wrote a Python code to do MD for a Lennard-Jones fluid but the VACF is wrong: What might be the problem? [closed]

I am trying to write a molecular dynamics simulation for a Lennard-Jones fluid in a box with periodic boundary conditions. The box has no net momentum. I am writing this in python.

I have written a library of functions to set up my box of particles. Then I am implementing it in a separate script. Here is dynamics.py:

import itertools
import numpy as np
import random
import time
import math

random.seed(time.time())

#create a box of particles
#make object Box which will hold all the particles
class Box:
def __init__(self, numberOfParticles, boxLength, dimension, sigma, epsilon, temperature, dt):
self.numberOfParticles = numberOfParticles
self.boxLength = boxLength
self.dimension = dimension
self.sigma = sigma
self.epsilon = epsilon
self.temperature = temperature
self.dt = dt #time step
##### non given quantities
self.nrho = numberOfParticles/(boxLength**(dimension)) #number density
self.particlePositions = np.zeros((numberOfParticles, dimension)) #do a cubic lattice
self.particleVelocities = self.boxLength*(np.random.rand(numberOfParticles, dimension)-0.5) #assign randomly
self.particleForces = np.zeros((numberOfParticles, dimension))
#
#
#now to evaluate energy of configuration
#evaluating kinetic energy of the system

def latticePositions(self):

pointsInLattice = math.ceil(self.numberOfParticles**(1/3))

spots = np.linspace(0, self.boxLength, num=pointsInLattice, endpoint=False)
count = 0

for p in itertools.product(spots, repeat=3):
p = np.asarray(list(p))
self.particlePositions[count, :] = p
count += 1
if count>self.numberOfParticles-1:
break
#
return self
#

def evaluateKineticEnergy(self):
#square every element, add up the elements of each row
kineticEnergy = 0.5*np.sum(np.square(self.particleVelocities))
return kineticEnergy
#
#I will be selecting a particle, and summing up all the potential energy arising
#due to interactions with every other particle

def evaluatePotentialEnergy(self):
energy = 0
for i in range(self.numberOfParticles):
for j in range(i+1, self.numberOfParticles):
displacement = self.particlePositions[i,:]-self.particlePositions[j,:]
for k in range(self.dimension):
if abs(displacement[k])>self.boxLength/2:
displacement[k] -= self.boxLength*np.sign(displacement[k])
r = np.linalg.norm(displacement,2) #finding euclidean distance between two particles
energy += (4*self.epsilon*((self.sigma/r)**12-(self.sigma/r)**6)) #evaluating potential energy, multiply by 2?
return energy
#sum of potential and kinetic energy is equal to the total energy

def evaluateTotalEnergy(self):
totalEnergy = self.evaluatePotentialEnergy()+self.evaluateKineticEnergy()
#end of energy calculations
#
#
#find the force each particle is experiencing due to the other particles
#force = - gradient of potential
def evaluateForce(self):
self.particleForces = np.zeros((self.numberOfParticles, self.dimension))
def LJForce(displacement):
r = np.linalg.norm(displacement, 2)
force = 48/(r**2)*(1/(r**12)-0.5*1/(r**6))*displacement
return force
for i in range(self.numberOfParticles):
for j in range(i+1, self.numberOfParticles):
rij = self.particlePositions[i,:]-self.particlePositions[j,:]
for k in range(self.dimension):
if abs(rij[k])>self.boxLength/2:
rij[k] -= self.boxLength*np.sign(rij[k])
rji = -rij
self.particleForces[i,:] += LJForce(rij)
self.particleForces[j,:] += -self.particleForces[i,:]
return self
#end of force evaluations
#make sure total momentum of box is zero
def stationaryCenterOfMass(self):
#ensure center of mass is stationary
v_cm = np.mean(self.particleVelocities, axis=0)
self.particleVelocities = self.particleVelocities - v_cm
return self
#
#
#
def VelocityVerletTimeStepping(currentBox):
previousParticleForces = currentBox.particleForces
currentBox.particlePositions = (currentBox.particlePositions + currentBox.particleVelocities*currentBox.dt + 0.5*currentBox.particleForces*(currentBox.dt)**2)%(currentBox.boxLength)
currentBox = currentBox.evaluateForce()
currentBox.particleVelocities = currentBox.particleVelocities + 0.5*(previousParticleForces + currentBox.particleForces)*currentBox.dt
return currentBox
#


Now I am calling these functions to perform time-evolution. I am evaluating energy and ACF and plotting them against time to see if I have done this right.

import dynamics
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import math

#parameters of simulation
numberOfParticles = 100
dimension = 3
sigma = 1
epsilon = 1
temperature = 1
#max_iterations = 100
boxLength = 10*sigma
dt = 0.001 #time step
kb = 1 #boltzmann
##
#
#

nmax = 100 #number of time steps to take

#set up the box

currentBox = dynamics.Box(numberOfParticles, boxLength, dimension, sigma, epsilon, temperature, dt)

#placing particles in a lattice
currentBox = currentBox.latticePositions()

#ensuring box has no net momentum
currentBox = currentBox.stationaryCenterOfMass()

#calculating forces on particles in the box
currentBox = currentBox.evaluateForce()

#making a list of particle positions and velocities for acf and what not
particlePositionsList = [currentBox.particlePositions]
particleVelocityList = [currentBox.particleVelocities]

#making a list of energies at various time steps to plot later
energy = np.zeros(nmax+1,)
energy[0] = currentBox.evaluateTotalEnergy()

timepoints = np.arange(nmax+1)*currentBox.dt

#start time stepping routine
#time points = 0, 1, 2, ..., nmax
for i in range(nmax):

#do the time step, knowing that currentBox already knows the particle forces at the moment
currentBox = dynamics.VelocityVerletTimeStepping(currentBox) #evaluates forces on particles, updates particle positions and velocities
energy[i+1] = currentBox.evaluateTotalEnergy()

particlePositionsList.append(currentBox.particlePositions)
particleVelocityList.append(currentBox.particleVelocities)

#
#print(energy)

ACF = np.zeros(nmax+1,)

for i in range(nmax+1):
for j in range(nmax+1-i):
ACF[i] = ACF[i] + np.sum(particleVelocityList[j]*particleVelocityList[j+i])
#ACF[j] = ACF[j] + np.sum(particleVelocityList[i]*particleVelocityList[j+i]) #this one works
#
ACF[i] = ACF[i]/(nmax+1-i)
#

plt.plot(timepoints, energy)
plt.title("Energy over time")
plt.xlabel("time")
plt.ylabel("energy")
plt.ylim(np.amin(energy)*0.999, np.amax(energy)*1.001)
plt.show()

plt.plot(timepoints, ACF)
plt.title("Normalized VACF plot")
plt.xlabel("time")
plt.ylabel("VACF")
plt.show()



Here are the results:

Energy looks good, but ACF looks really wrong. This is how it is supposed to look like:

I am unsure where I am going wrong here. I am a lone software engineer thrust into the world of physics and molecular modelling, so any advice you have would be appreciated!!

Edit 1: After initializing velocity according to a Gaussian distribution,

self.particleVelocities = self.boxLength*(np.random.normal(0, 1, (numberOfParticles, dimension))) #assign velocity as per normal distr


These are the figures I get:

• “Energy looks good” no that does not look good. MD with with LJ force field should show an energy drift at the beginning and if you use NVT ensemble (what’s your thermostat?) then it should reach a near constant values after the defined relaxation in your thermostat. So a constant energy from beginning is sign of a something wrong fundamentally in your code. Commented Jul 21, 2020 at 4:06
• It's been a long time for me, but I remember my instructor insisting that we initialize the velocities according to a Gauss distribution. Does not look like this is the case here to me. Commented Jul 21, 2020 at 4:12
• I would strongly suggest spending some time on Frenkel and Smit's "Understanding Molecular Simulation" textbook. They walk you through the physics and they also give pseudocode for many applications, so this will help you debug (they also have an example on diffusion). Debugging is already hard enough to do when you know the physics, so you can't expect much progress if you don't understand the algorithms behind it. Commented Jul 21, 2020 at 8:29
• I also noticed that you appear to take 100 steps? That's never enough. I looked at an old input file of mine and saw 30k steps. It may be instructive to plot $E_\text{kin}$ and $E_\text{pot}$ in one plot, not just their sum. I expect the latter to drop over time. Commented Jul 21, 2020 at 17:45
• I recommend looking at Allen & Tildesley - it is one of the two bibles, Frenkel & Smit being the other. However, A & T came out with a version only a couple years ago and have all their codes running in both python and fortran: github.com/Allen-Tildesley/examples/tree/master/python_examples Commented Jul 21, 2020 at 18:33