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()
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()
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
```

Matter Modeling Asked by user207526 on July 22, 2020

0 Answers1 Asked on January 5, 2022 by quantumx

computational chemistry coupled cluster gaussian quantum chemistry

1 Asked on January 5, 2022

ab initio calculations density functional theory functional many body perturbation theory quantum chemistry

1 Asked on January 3, 2022

density functional theory mathematical modeling quantum chemistry spectroscopy td dft

3 Asked on January 1, 2022 by anyon

1 Asked on January 1, 2022 by bnd

1 Asked on December 22, 2021

lammps molecular dynamics molecular interactions molecular modeling software

2 Asked on December 22, 2021

1 Asked on December 20, 2021

electronic structure fcidump integrals one topic per answer software

1 Asked on December 20, 2021 by c-alexander

density functional theory from phys se quantum chemistry spectroscopy td dft

1 Asked on December 20, 2021

1 Asked on December 18, 2021 by kfc

1 Asked on December 16, 2021

density functional theory from phys se mathematical modeling statistical thermodynamics variational principle

1 Asked on December 16, 2021 by luphys

alloys density functional theory from phys se relativistic software

1 Asked on December 16, 2021

1 Asked on December 13, 2021 by roman-korol

atomic physics full configuration interaction high precision molecular modeling periodic system

1 Asked on December 13, 2021

band structure density functional theory from phys se quantum espresso

9 Asked on December 11, 2021

high performance computing matter modeling specific technology one topic per answer open source software

4 Asked on December 9, 2021 by georgios-karaiskakis

computational chemistry density functional theory functional history

2 Asked on December 9, 2021 by koroma

band structure lattice dynamics perovskites phonons symmetry

Get help from others!

Recent Answers

- Philipp on How do i draw a ray in unity
- DMGregory on MouseLook Script “Pops” back to the last value when the script is enabled after being disabled or destroyed
- Justin Markwell on Unity app crashes when using unmodified custom Android manifest (didn’t find class “UnityPlayerActivity”)
- eric_kernfeld on How to test consistency of responses?
- kjetil b halvorsen on How to test consistency of responses?

Recent Questions

- MouseLook Script “Pops” back to the last value when the script is enabled after being disabled or destroyed
- Unity app crashes when using unmodified custom Android manifest (didn’t find class “UnityPlayerActivity”)
- How do i draw a ray in unity
- How to test consistency of responses?
- How can I understand these variograms?

© 2022 AnswerBun.com. All rights reserved.