# Calculation of Mean Square Displacement for Brownian dynamics system with Lennard Jones interactions in python3

I have a problem getting a sensible result for the Mean Square Displacement (MSD) for a simulation of $$N$$ particles under Brownian dynamics with Lennard-Jones interaction between them with or without periodic conditions applied.

According to my class notes, the result should look similar to this: Instead, with periodic boundary conditions, I get something like this: and without periodic boundary The system begins with $$N$$ particles positioned in a lattice way. One calculates the force on the particles as the derivative of the Lennard-Jones potential, for each axes $$x$$ and $$y$$. So in the code, the final force in each direction gets computed with the use of the chain rule. For the total displacement of a particle, one also adds a random number drawn from a Gaussian distribution. This process is repeated for each particle and for a series of time-repetitions.

The program saves every 1000 times, the positions of the particles in a dump lammbs like format ( for visualizing with ovito or vmd), it saves the distances between the particles for calculating in the end the radial distribution function, and also calculates the MSD by averaging over particles the sum of the particles' distance from their original position: $$ = \frac{1}{N}\sum_i ^N (|\bar(r_i (t)) -\bar(r(0) ) |)^2 ,$$ where t is the specific time of the calculation. The several MSD are appended to a list to be plotted ασ one function of time.

I thought that perhaps a problem in my code could be the sign of the forces, but I am sure it is now correct, the signs from the derivatives $$\frac{dV}{dr}$$ and $$\frac{dV}{dx}$$ and respectively for $$y$$, together with an overall minus ($$-$$) sign give the correct direction to the forces.

So, in the end, I have no idea what goes wrong with the program, and after having checked it 1000 times over the last 5 days, I am at a moment where I am just changing things in the code, only to see what will come as a result - which is really frustrating.

In my laptop (Intel i5, 8GB RAM) it takes around 30 minutes for the code to run, which I think is relatively correct for the amount of calculations.

My code is this (you can see the main program and all the function calls after MAIN PROGRAM comment line):

from pylab import *
import pylab
import random
import math
import matplotlib.pyplot as plt
import numpy as np

#------------------------ΠΡΟΣΗΜΑ (+) Frx, Fry, RDF----------------------

def dump(particles,step,n):
fileoutput = open('coord.txt', 'a')
fileoutput.write("ITEM: TIMESTEP \n")
fileoutput.write("%i \n" % step)
fileoutput.write("ITEM: NUMBER OF ATOMS \n")
fileoutput.write("%i \n" % n)
fileoutput.write("ITEM: BOX BOUNDS \n")
fileoutput.write("%e %e xlo xhi \n" % (0, L))
fileoutput.write("%e %e xlo xhi \n" % (0, L))
fileoutput.write("%e %e xlo xhi \n" % (-0.25, 0.25))
fileoutput.write("ITEM: ATOMS id type x y z \n")
i = 0
while i < n:
x = particles[i]
y = particles[i]
#fileoutput.write("%i %i %f %f %f \n" % (i, 1, x*1e10, y*1e10, z*1e10))
fileoutput.write("%i %i %f %f %f \n" % (i, 1, x, y, 0))
i += 1

fileoutput.close()

N = 49
L = 10
rcut = 4
meanz = 0
varz = 1
particles = []
x=0
y=0
epsilon =1
sigma = 1
# tau = sigma**2*ksi/(kT)

# Starting time
t_0 = 0

# Time increments
dt = 10**(-4)  # dt/tau

# Ending time
tsim = 5*10**(5)*dt # tsim/tau

numberofsaves=int(tsim/dt/1000 +1)
rdf_dist=[]

# Produce random particles and avoid overlap:

times = np.arange(t_0, tsim+dt, dt)
tindex=0
step=0
rMeansq=[]

def rdists(particles,N,L):
distances=[]
for i in range(N):
for j in range(i+1,N):
dii=distancesq(particles[i], particles[j],L)
distances.append(dii)
# print(distances)
# input()
return distances

def distancesq(particle1, particle2, L):
# Calculate distances
dx = particle1-particle2
dy = particle1-particle2

# Minimum image convention
if dx > L/2:
dx=L%dx
elif dx < -L/2:
dx=L%dx
if dy > L/2:
dy=L%dy
elif dy < -L/2:
dy=L%dy
return math.sqrt(dx**2+dy**2)

def wrap(particle):
'''Apply perodic boundary conditions.'''
if particle > L:
particle = particle%L
elif particle < 0:
particle =particle%L
if particle > L:
particle =particle%L
elif particle < 0:
particle =particle%L

return particle

def show_conf(particles, sigma, L, t):
pylab.axes()
for [x,y] in particles:
cir = pylab.Circle((x, y), radius=sigma,  fc='r')
pylab.axis('scaled')
pylab.title("title")
pylab.axis([-2*L,2*L, -2*L,2*L])
pylab.savefig('snapshot=%s' % t)
# pylab.show()
pylab.close()

def Forces(fp,particles,L):   #fp: i-particle
Frx=0
Fry=0
for particle2 in particles:
if particle2 != particles[fp]:
dist = distancesq(particles[fp], particle2,L)
if dist <= rcut and dist!=0:
dVdr = 4*epsilon*(-12*dist**(-13) + 6*dist**(-7))
drdx = (1/dist)*(particles[fp]-particle2)
drdy = (1/dist)*(particles[fp]-particle2)
Frx = Frx -dVdr*drdx
Fry = Frx -dVdr*drdy
return Frx, Fry

def RDF(N,particles, L,dist):
minb=0
maxb=10
nbin=200
width=(maxb-minb)/(nbin)

rings=np.linspace(minb, maxb, nbin)
d=np.asarray(dist).flatten()
rDf = np.histogram(d, rings ,density=True)
prefactor = (1/( np.pi*(N/L**2)))
rDf = [prefactor*rDf, 0.5*(rDf[1:]+rDf[:-1])]
rDf=np.multiply(rDf,1/(rDf*( width )))

plt.figure()
plt.plot(rDf,rDf)
plt.xlabel("r")
plt.ylabel("g(r)")
plt.savefig("g(r)")

def msd_dists(particles,pos_start, N):
msd_distances=[]
for j in range (0,N):
msd_r = ((particles[j]- pos_start[j] )**2 + (particles[j] - pos_start[j])**2 )**(1/2)
msd_distances.append(msd_r)
return msd_distances

#--------------------------------------------MAIN PROGRAMM----------------------------------

# Create the lattice
N_sqrt=math.sqrt(N)
N_sqrt=int(N_sqrt)
delxy=L/(2*N_sqrt)
for k in range(0,N):
particles = [[delxy + i * 2*delxy, delxy + j * 2*delxy] for i in range(N_sqrt) for j in range(N_sqrt)]
pos_start = [[delxy + i * 2*delxy, delxy + j * 2*delxy] for i in range(N_sqrt) for j in range(N_sqrt)]

# Simulation

for t in times:
# distances = []
for j in range(0,N):
Fx,Fy = Forces(j, particles,L)
z = random.gauss(meanz, varz)
# if(particles[j]+((2*dt*sigma**2)**(1/2))*z > L or particles[j]+((2*dt*sigma**2)**(1/2))*z<0):   #reflecting boundary conditions
#   z=-z
particles[j] = particles[j] + ((2*dt*sigma**2)**(1/2))*z + dt*Fx/epsilon
z = random.gauss(meanz, varz)
# if(particles[j]+((2*dt*sigma**2)**(1/2))*z > L or particles[j]+((2*dt*sigma**2)**(1/2))*z<0):    #reflecting boundary conditions
#   z=-z
particles[j] = particles[j] + ((2*dt*sigma**2)**(1/2))*z + dt*Fy/epsilon
particles[j]=wrap(particles[j])

if (step%1000 == 0):
gen_dists=rdists(particles,N,L)
rdf_dist.append(gen_dists)
msd_distances=msd_dists(particles,pos_start,N)
rmeansq=sum(np.square(msd_distances))/N
rMeansq.append(rmeansq)
dump(particles,t,N)
print(step)
print()
# save distances-positions of all the particles (every 1000 steps) in file
step+=1

# Distribution Function
# rdf_dist = rdists(particles,N,L)
RDF(N,particles,L,rdf_dist)

# ---------Plot MSD---------
plt.figure()
plt.plot(np.linspace(0,times[-1],num=numberofsaves),rMeansq)
plt.title('MSD')
plt.savefig("MSD ")


I am not sure what is the problem here; I don't think I have missed some normalization nor I thing that there is a mistake to the parameters of the problem- and even then, I don't think it would be so severe as to produce such a different result.

Any ideas and references are truly welcomed.

Edit:

If I plot the MSD for small values I get something like this: Can it be that there is a problem with the dimensional analysis of my implementation of that some constants are not correct?

New contributor
Constantine Black is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
• If someone is downvote, could he give a reason for doing so? This is really not helpful. If the post can be written better, then comment, don't downvote. Otherwise, I see no reason to have this as bad post, since it shows effort, code, and tries to explain the context of the problem as good as possible. – Constantine Black Jan 12 at 18:22
• I guess people downvote, because this is a "please debug my code"-type of question. Anyway, it seems to be a simple particle simulation but I don't see any information about particle velocities. You seem to update the particle positions with the forces. That is certainly wrong. – Bort yesterday
• @Bort So, if you are working constantly with some program over a week, and you don' t seem able to find what's going wrong with the implementation, and you visit a question and answer site to get some insight from other that might have faced a similar problem, you are not trying to get help in understanding, but you wan't others to debug your code... Ok. The second part of your comment though is helpful. If you are correct, then the updating part is wrong, but this is what was given as certain from the tutors of the tutorial: the equation was given. Interesting. Thanks for your time. – Constantine Black yesterday
• And if I check my notes on the Ermak-McCammon implementation of brownian dynamics, there is a non-velocity term dependent on the forces and the diffusion coefficient for updating the position of the particles, that is no velocity is needed. – Constantine Black yesterday
• One can check for example the last page of this pdf, where you af course begin with langevin equation, but you reach a non-velocity equation solving for the position variable. dasher.wustl.edu/bio5476/reading/stochastic.pdf – Constantine Black yesterday