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 = 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
1 Asked on January 5, 2022 by quantumx
1 Asked on January 5, 2022
1 Asked on January 3, 2022
3 Asked on January 1, 2022 by anyon
1 Asked on January 1, 2022 by bnd
1 Asked on December 25, 2021
1 Asked on December 22, 2021
2 Asked on December 22, 2021
1 Asked on December 20, 2021
1 Asked on December 20, 2021 by c-alexander
1 Asked on December 20, 2021
1 Asked on December 18, 2021 by kfc
1 Asked on December 16, 2021
1 Asked on December 16, 2021 by luphys
1 Asked on December 16, 2021
1 Asked on December 13, 2021 by roman-korol
1 Asked on December 13, 2021
9 Asked on December 11, 2021
4 Asked on December 9, 2021 by georgios-karaiskakis
Get help from others!