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

Matter Modeling Asked by user207526 on July 22, 2020

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: