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