Here is my model of the molecule evolution…

```
# Argon function
function Argon(u,p,t)
np = Int(ceil(length(u)/4))
#
x = u[1:np]
y = u[np+1:2*np]
px = u[2*np+1:3*np]
py = u[3*np+1:4*np]
#
nw = 4
Xk = p[1]
Yk = p[2]
sig = p[3]
epsi = p[4]
m = p[5]
kB = p[6]
#
dfimdx=zeros(np)
dfimdy=zeros(np)
tmp = 1:np
for i in 1:np
r = zeros(np)
tm = collect(copy(tmp))
deleteat!(tm,i)
for j in 1:np
r[j] =sqrt((x[i]-x[j])^2 +(y[i]- y[j])^2)
end
for j in 1:nw
k = tm[j]
dfimdx[i] = dfimdx[i] +
24*epsi*((sig/r[k])^6 - 2*(sig/r[k])^12)*(x[i] - x[k])/r[k]^2
dfimdy[i] = dfimdy[i] +
24* epsi*((sig/r[k])^6 - 2*(sig/r[k])^12)*(y[i] - y[k])/r[k]^2
end
end
#
dfivdx = zeros(np)
dfivdy = zeros(np)
for i in 1:np
R = zeros(nw)
Xg= zeros(nw)
Yg= zeros(nw)
L2 = zeros(nw)
for k in 1:nw
L2[k] = (Xk[2,k]-Xk[1,k])^2 +(Yk[2,k]-Yk[1,k])^2
gam =((x[i]-Xk[1,k])*(Xk[2,k]-Xk[1,k]) +
(y[i]-Yk[1,k])*(Yk[2,k]-Yk[1,k]))/L2[k]
Xg[k] = Xk[1,k] + (Xk[2,k]-Xk[1,k])*gam
Yg[k] = Yk[1,k] + (Yk[2,k]-Yk[1,k])*gam
R[k] = sqrt((x[i]-Xg[k])^2 + (y[i]-Yg[k])^2)
end
for k in 1:nw
dfivdx[i] = dfivdx[i] +
24*epsi*((sig/R[k])^6 - 2*(sig/R[k])^12)*
(x[i]-Xg[k])/R[k]^2*(Yk[2,k]-Yk[1,k])^2/L2[k]
dfivdy[i] = dfivdy[i] +
24*epsi*((sig/R[k])^6 - 2*(sig/R[k])^12)*
(y[i]-Yg[k])/R[k]^2*(Xk[2,k]-Xk[1,k])^2/L2[k]
end
end
fpx = - (dfimdx + dfivdx)*kB
fpy = - (dfimdy + dfivdy)*kB
fx = px/m
fy = py/m
du = [fx...,fy...,fpx...,fpy...]
return du
end
#
```

Here is my RK4 solver:

```
# Runge Kutta 4 solver
function RK4(func,u,tspan,p)
#
#
t0 = tspan[1]
t1 = tspan[2]
h = tspan[3]
N = Int(ceil((t1-t0)/h))
nu = length(u)
#
T = zeros(N+1)
U = zeros(N+1,nu)
#
T[1] = t0
U[1,:] = u'
for i in 1:N
#
# -----------------------------------------------------------------
# Runge-Kutta-4 updating
#
# Computing vector fields and temporary quantities k1-k4
#
k1 = func(u,p,0)
k2 = func(u+0.5*h*k1,p,0)
k3 = func(u+0.5*h*k2,p,0)
k4 = func(u+h*k3,p,0)
#
# Updating states
u = u + h*(k1+2*k2+2*k3+k4)/6
# ------------------------------------------------------------------
# Storing time and states
#
T[i+1] = T[i]+h
U[i+1,:] = u'
#
# Next time
end;
#
return T,U
end
#
```

Here is my set-up for the RK4 solver:

```
#
using Plots, LaTeXStrings, Statistics, DifferentialEquations; pyplot()
LS1 = :solid
LS2 = :dash
LS3 = :dot
LW1 = 2
LW2 = 1
#
#
# Setting up Parameters
m = 6.6359e-23 # mass per molecule [g]
kB = 1.3806e-23 # Boltzmann's constant - converts Kelvin to Joule
#
X01 = 0. # wall position, Angstrom
X11 = 0 # --"--
X12 = 50 # --"--
X04 = 50 # --"--
X14 = X01
X02 = X11
X03 = X12
X13 = X04
Xk = [X01 X02 X03 X04;X11 X12 X13 X14]*1e-10
#
Y01 = 0. # Angstrom
Y11 = 50 # --"--
Y12 = 50 # --"--
Y04 = 0 # --"--
Y14 = Y01
Y02 = Y11
Y03 = Y12
Y13 = Y04
Yk = [Y01 Y02 Y03 Y04;Y11 Y12 Y13 Y14]*1e-10
#
epsi = 119.8 # Lennard-Jones parameter, Kelvin
sig = 3.405 # Lennard Jones parameter, Angstrom
sig = sig*1e-10 # sig in m
m = m*1e-3 # mass in kg
#
p = [Xk,Yk,sig,epsi,m,kB]
#
# Initial values for positions
#
x = [5.,15,25,35,45] # Ångström
y = [5.,5,5,5,5] # Ångstrom
x = x*1e-10
y = y*1e-10
#
# Initial values for momenta
px = m*[3.,3,-3,-3,3] # kg*Angstrom/pico second
py = m*[3.,-3,3,-3,3] # --"--
px = px*1e2 # conversion to kg*m/s
py = py*1e2 # conversion to kg*m/s
#
u0 = [x...,y...,px...,py...]
# Time interval
#
t0 = 0. # picosecond
t1 = 500. # picosecond
h = 0.01 # picosecond
t1 = t1*1e-12 # in second
h = h*1e-12 # in second
tspan1 = (t0,t1,h)
#
@time T,U = RK4(Argon,u0,tspan1,p)
N = length(T)
# Unpacking results
X = U[:,1:5]
Y = U[:,6:10]
PX = U[:,11:15]
PY = U[:,16:20]
# Convert position units from m to Angstrom
X=X*1e10
Y=Y*1e10
```

Plotting…

```
# Plotting results
plot(X[:,1],Y[:,1],lw=LW2,ls=LS1,lc=:black,label="#1")
plot!(X[:,2],Y[:,2],lw=LW2,ls=LS3,lc=:red,la=0.7,label="#2")
plot!(X[:,3],Y[:,3],lw=LW2,ls=LS3,lc=:gold,la=0.7,label="#3")
plot!(X[:,4],Y[:,4],lw=LW2,ls=LS3,lc=:cornflowerblue,la=0.7,label="#4")
plot!(X[:,5],Y[:,5],lw=LW2,ls=LS3,lc=:green,la=0.7,label="#5")
plot!(xlim=(0,50),ylim=(0,50))
plot!(xlabel=L"$x$ [Ångstrom]", ylabel=L"$y$ [Ångstrom]")
plot!(title=L"5 Argon molecules during $500$ ps")
#
Ns = Int(ceil(N/20))
plot(X[1:Ns,1],Y[1:Ns,1],lw=LW1,ls=LS1,lc=:black,label="#1")
plot!(X[1:Ns,2],Y[1:Ns,2],lw=LW1,ls=LS3,lc=:red,label="#2")
plot!(X[1:Ns,3],Y[1:Ns,3],lw=LW1,ls=LS3,lc=:gold,label="#3")
plot!(X[1:Ns,4],Y[1:Ns,4],lw=LW1,ls=LS3,lc=:cornflowerblue,label="#4")
plot!(X[1:Ns,5],Y[1:Ns,5],lw=LW1,ls=LS3,lc=:green,label="#5")
plot!(xlim=(0,50),ylim=(0,50))
plot!(xlabel=L"$x$ [Ångstrom]", ylabel=L"$y$ [Ångstrom]")
plot!(title=L"5 Argon molecules during $25$ ps")
scatter!((X[1,1],Y[1,1]),mc=:black,ms=6,label="")
scatter!((X[1,2],Y[1,2]),mc=:red,ms=6,label="")
scatter!((X[1,3],Y[1,3]),mc=:gold,ms=6,label="")
scatter!((X[1,4],Y[1,4]),mc=:cornflowerblue,ms=6,label="")
scatter!((X[1,5],Y[1,5]),mc=:green,ms=6,label="")
```