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

Question

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
radius=0.4
# 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.gca().add_patch(cir)
    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)]


snapshot_before=show_conf(particles, radius, L, t='start')


# 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()
        show_conf(particles, radius, L, step)
                        # save distances-positions of all the particles (every 1000 steps) in file
    step+=1

snapshot_after=show_conf(particles, radius, L,t='end')

# 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?

Answer 1

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.