# 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: $$\langle r^2 \rangle = \frac{1}{N}\sum_i ^N (|(\overline r_i (t)) -(\overline 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][0]
y = particles[i][1]
#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[0]-particle2[0]
dy = particle1[1]-particle2[1]

# 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[0] > L:
particle[0] = particle[0]%L
elif particle[0] < 0:
particle[0] =particle[0]%L
if particle[1] > L:
particle[1] =particle[1]%L
elif particle[1] < 0:
particle[1] =particle[1]%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][0]-particle2[0])
drdy = (1/dist)*(particles[fp][1]-particle2[1])
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], 0.5*(rDf[1][1:]+rDf[1][:-1])]
rDf[0]=np.multiply(rDf[0],1/(rDf[1]*( width )))

plt.figure()
plt.plot(rDf[1],rDf[0])
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][0]- pos_start[j][0] )**2 + (particles[j][1] - pos_start[j][1])**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][0]+((2*dt*sigma**2)**(1/2))*z > L or particles[j][0]+((2*dt*sigma**2)**(1/2))*z<0):   #reflecting boundary conditions
#   z=-z
particles[j][0] = particles[j][0] + ((2*dt*sigma**2)**(1/2))*z + dt*Fx/epsilon
z = random.gauss(meanz, varz)
# if(particles[j][1]+((2*dt*sigma**2)**(1/2))*z > L or particles[j][1]+((2*dt*sigma**2)**(1/2))*z<0):    #reflecting boundary conditions
#   z=-z
particles[j][1] = particles[j][1] + ((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?

• 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. Commented Jan 12, 2020 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
Commented Jan 17, 2020 at 18:51
• @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. Commented Jan 17, 2020 at 19:44
• 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. Commented Jan 17, 2020 at 19:49
• 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 Commented Jan 17, 2020 at 19:51

This will not an answer to your problem, more an excessive comment and few things you might consider, when writing such code even for self educational purpose.

Constants You asked whether your dimensionless code uses correct constants / normalization. So let's have a look. For this question it helps to write down the equations, you want to solve. Here I get from your comments, you want to solve the Ermak-McCammon equation ( I'll I show this on the standard MD equation, and leave the adaption to the Ermak-McCammon equation to you):

$$m\ddot{r}= F(r)$$

Your forces come from the particles' Lennard-Jones potentials. $$V_{LJ}(r)=4\epsilon\left(\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^6\right)$$

A short reminder of units of $$[\sigma]=m$$, $$[\epsilon]=\frac{kg m^2}{s^2}$$.

Next, you do a nondimensionalization of your equation:

• $$t=\tau t'$$
• $$r=r_0 r'$$

This gives you following equation for the variables $$r'$$ and $$t'$$: $$\frac{m r_0}{\tau^2}\ddot{r}'=-\frac{1}{r_0} \nabla' V_{LJ}(r_0 r')$$ $$\ddot{r}'=-\frac{\tau^2}{m r_0^2} \nabla' V_{LJ}(r_0 r')$$

Dropping the primes, you get the following equations for your implementation $$\ddot{r}=-4 \epsilon'\nabla \left(\left(\frac{\sigma'}{r}\right)^{12}-\left(\frac{\sigma'}{r}\right)^6\right)$$

Here, $$\epsilon'=\frac{\epsilon\tau^2}{m r_0^2}$$ and $$\sigma'=\frac{\sigma}{r_0}$$ have the dimension 1, which was the intention.

You did the nondimensionalization, which is good, when you developed your code, but you did not document in form of a couple of functions which express the scaled result in SI units again. Correct units help to determine errors if one has a bit of experience in which order of magnitude one expects results.

So, let's say you want to simulate Ar in 10 nm box, for which you can easily find some parameters, $$\sigma = 3.41 \text{ angstrom}$$, $$\epsilon=119.8 \text{ K }\cdot k_B$$. Now, you want to pick the scaling parameters $$\tau, r_0$$ such that the dimensionless constants are in the order of 1. So let's pick $$r_0=\sigma$$ and $$\tau=1\text{ ps}$$, since you would expect a $$dt=1 \text{ fs}$$ for time stepping (of course, you could also choose $$\tau$$ based on $$\epsilon$$ rather then picking some number). $$\epsilon'=0.213$$, $$\sigma '=1$$. $$L/r_0=29.326$$ will be the dimensionless box size.

From the ideal gas law, you get the right amount of particles for a given temperature and pressure in your simulation box, e.g. $$T=300 \text{ K}$$ and $$p=1013 \text{ hPa}$$ lead to $$N\approx24$$.

Let's compare to your chosen values:

N = 49
L = 10
epsilon = 1
sigma = 1


So you have a much stronger interaction and much more particle in a smaller box in 2D. This means that your system has a much much higher energy. This is per se not a problem but can lead to serious issues. e.g. strong acceleration, relativistic velocities, etc.

Error prone functions Writing code error free is very difficult. So I recommend to put your functions under unit test, which cover all edge cases. I say this because there are a couple of statements in your code which don't make sense, e.g.

        dist = distancesq(particles[fp], particle2,L)
if dist <= rcut and dist!=0:


Why is it possible in your code that two particles (after filtering particle2) have equal positions when there is a repulsive force $$\propto \frac{1}{r^{12}}$$ between them? Although it might happen this should be a unlikely event because you ask for equality of floating point numbers (you might even consider runs where this happens as invalid physics). So I suspect that there is a bug. Another example:

drdy = (1/dist)*(particles[fp][1]-particle2[1])


Although you determine the distance between two particles in distancesq correctly with their images, you ignore to do this for their difference in the derivative. So this is likely to be wrong as well.

There is probably more. Since I haven't written any code for Ermak-McCammon equation I cannot comment further on possible issues with that. I hope this helps you out.

• Really thank you for your notes and thoughts. You are correct of course, and I surely understand the second part of your answer. I have to admit, I didn't thought that much( and it seems I needed) what's going on with the constants: the problem gave me those values, so I thought I had to use them as they were and the nondimensionalization would work out. I will try to re-check the implementation with other values too, and see what I would be getting. Again, thank you a lot. Commented Jan 21, 2020 at 11:58