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


#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:
        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()
        return totalEnergy
    #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()



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.ylim(np.amin(energy)*0.999, np.amax(energy)*1.001)

plt.plot(timepoints, ACF)
plt.title("Normalized VACF plot")

Here are the results:

enter image description here

enter image description here

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

enter image description here

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: enter image description here

enter image description here

  • 1
    $\begingroup$ “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. $\endgroup$ Commented Jul 21, 2020 at 4:06
  • 1
    $\begingroup$ 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. $\endgroup$
    – TAR86
    Commented Jul 21, 2020 at 4:12
  • 2
    $\begingroup$ 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. $\endgroup$
    – Godzilla
    Commented Jul 21, 2020 at 8:29
  • 2
    $\begingroup$ 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. $\endgroup$
    – TAR86
    Commented Jul 21, 2020 at 17:45
  • 2
    $\begingroup$ 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 $\endgroup$
    – B. Kelly
    Commented Jul 21, 2020 at 18:33