Modeling multiple particles in a box

We consider \(N_{tot} = 25\) identical particles in a two-dimensional square box of size \(a = 6.25\). Below shows a configuration of the system where the particles are located at evenly distributed grid points.

ntot = 25
a = 6.25
grid_x = seq(0, 5, by = 1.25)
grid_y = seq(0, 5, by = 1.25)
mat_xy = as.matrix(expand.grid(grid_x, grid_y))
x0 = c(t(mat_xy))
xy0 = matrix(x0, ncol = 2, byrow = TRUE)
plot(xy0[,1], xy0[,2], type = 'p', xlab = "X", ylab = "Y", xlim = c(0, a), ylim = c(0,a))

There are pairwise interactions between any two particles governed by the following force function:

\[F(r) = \frac{1}{r^{13}} - \frac{1}{r^7} \tag{1}.\]

par(mar = c(4, 4, 1, 1))   #set plot margins
force <- function(r){
  r3 = r*r*r
  r6 = r3*r3
  f = (1/r6/r6 - 1/r6)/r
  return (f)
}
r_all = seq(0, 4, by = 0.01)
plot(r_all, sapply(r_all, force), type = "l", col = 1, 
     xlab = "r", ylab = "F", xlim = c(0,4), ylim = c(-0.5,1))
abline(h = 0, lty = 3, col = 2)

\(r\) is the Euclidean distance between the particles. The force in Equation (1) is attractive (i.e., \(F(r) < 0\)) for \(r > 1\) and repulsive (i.e., \(F(r) > 0\)) for \(0 < r < 1\). The masses of all particles are 1. Note, this force function correponds to the famous Lennard Jones (LJ) potential.

par(mar = c(4, 4, 1, 1))   #set plot margins
E_LJ <- function(r){
  r3 = r*r*r
  r6 = r3*r3
  e = 1/12/r6/r6 - 1/6/r6
  return (e)
}
r_all = seq(0, 4, by = 0.01)
plot(r_all, sapply(r_all, E_LJ), type = "l", col = 1, 
     xlab = "r", ylab = "LJ Potential", xlim = c(0,4), ylim = c(-1,2))
abline(h = 0, lty = 3, col = 2)

Periodic boundary condition & minimum image convention

In this simulation, we consider the periodic boundary condition in that particles travel outside of the box would go inside the box from the other side.

To update the position x (a vector of size 2) of a particle:

x = c(-5.5, 4.2); a = 2
xnew = x - floor(x/a) * a # floor: the largest integer below x/a
x; xnew
## [1] -5.5  4.2
## [1] 0.5 0.2

There is no need to update the velocities.

Using the periodic boundary condition, we simulate a small box of particles to model the behavior of a much larger system consisting of the repetitive boxes. Since each particle can now present a lattice of particles of images, there are an infinite number of interactions to consider between any two particles \(i\) and \(j\). A common way to address this is to consider the minimum image convention, where individual particle in the simulation interacts with the closest image of the remaining particles in the system.

To compute the corresponding distance between two particles \(i\) and \(j\), we

xi = c(0.7, 1.9); xj = c(1.8, 1.2); a = 2
dx = xj - xi  # a vector of size 2
r1 = sqrt(sum(dx*dx))  # distance directly from xi and xj
dx = dx - round(dx/a)*a  # round: nearest integer
r_min = sqrt(sum(dx*dx)) # distance from the minimum image convention
r1; r_min
## [1] 1.30384
## [1] 1.140175
sqrt((2.7-1.8)**2 + (1.9 - 1.2)**2)  # this is actual calculation in r_min
## [1] 1.140175

Compute the total force

Now we devise a function to compute the x and y components of the total force for each particle.

The x-component of the force from particle \(j\) to particle \(i\) is

\[f_x(j \rightarrow i) = F(r_{ij}) \frac{(x_i - x_j)}{r_{ij}} \tag{2}\] , where \(r_{ij} = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2}\). (Note that positive \(F\) means repulsive interaction)

Similarly, the y-component of the force is

\[f_y(j \rightarrow i) = F(r_{ij}) \frac{(y_i - y_j)}{r_{ij}} \tag{3}.\] The total force of particle \(i\) is the summation of the pairwise forces from any particle \(j\) (where \(j \ne i\)) to particle \(i\).

# x-y components of the force from j to i
force_j2i <- function(i,j, x_all, a){
  # x_all a vector of 2 * Ntot
  xi = x_all[2*i-1]
  yi = x_all[2*i]
  xj = x_all[2*j-1]
  yj = x_all[2*j]
  
  dx = xi - xj
  dx = dx - round(dx/a)*a 
  dy = yi - yj
  dy = dy - round(dy/a)*a 
  
  rij = sqrt(dx**2 + dy**2)
  fij = force(rij)
  fxi = fij*dx/rij
  fyi = fij*dy/rij
  return(c(fxi,fyi))
}

force_all <- function(t, x_all, a){
  #Ntot: number of particles
  # x_all: a vector of 2 * Ntot
  # a: box size
  # f: a force vector of 2*Ntot
  Ntot = length(x_all)/2
  f = numeric(2*Ntot)
  for (i in 1:(Ntot-1)){
    for (j in (i+1):Ntot){
      fij = force_j2i(i,j, x_all, a)
      f[(2*i-1):(2*i)] = f[(2*i-1):(2*i)] + fij
      f[(2*j-1):(2*j)] = f[(2*j-1):(2*j)] - fij
    }
  }
  return(f)
}

The function force_all computes the forces from all pairwise interactions between the particles. The input argument for this function is a vector of size \(2N_{tot}\) for the \(x\) and \(y\) coordinates. The output of this function is a vector of size \(2N_{tot}\) for the \(x\) and \(y\) components of the total force \(f\) for each particle.

Update positions by the periodic boundary condition

After every MD integration step, particle positions need to be adjusted according to the periodic boundary condition.

bc_periodic <- function(x_all, a){
  # X_all: a vector of 2 * Ntot for x and y positions. 
  # a: box size
  x_all_new = x_all - floor(x_all/a) * a # floor: the largest integer below x/a
  return(x_all_new)
}

Verlocity Verlet (periodic boundary condition)

We simulated the dynamics of the particles according to

\[ m_i\frac{d^2x_i}{dt^2} = \sum_{j\neq i}{f_x(j\rightarrow i)} \tag{4}\]

, where \(m_i = 1\) (thus can be omitted) and the summation from the right hand side is over all particles other than \(i\). A similar equation is needed for \(y_i\).

# Velocity Verlet method for high dimensional systems (2D, 3D, or multiple particles)
# Updated with periodic boundary condition
velocity_verlet_bc <- function(f, t0, x0, v0, t.total, dt, ...){
  # f:  2nd derivative function
  # t0: initial time
  # x0: a vector of initial position x
  # v0: a vector of initial velocity v
  # t.total: total simulation time
  # dt: time step size 
  # a: box size (input in ellipsis)
  t_all = seq(t0, t.total, by=dt)
  nt_all = length(t_all)
  nx = length(x0)
  x_all = matrix(0, nrow = nt_all, ncol = nx)
  v_all = matrix(0, nrow = nt_all, ncol = nx)
  x_all[1,] = x0
  v_all[1,] = v0 
  for (i in 1:(nt_all-1)) {
    v_half = v_all[i,] + 0.5 * dt * f(t_all[i], x_all[i,],...)
    x_all[i+1,] = x_all[i,] + dt * v_half
    v_all[i+1,] = v_half + 0.5 * dt * f(t_all[i+1], x_all[i+1,],...)
    
    x_all[i+1,] = bc_periodic(x_all[i+1,], ...)
  }
  return(cbind(t_all, x_all, v_all))   # the output is a matrix of t, x, v for all time steps
}

Initial simulation for t.total = 10

We now perform a dynamics simulation of the box of the \(N_{total}=25\) particles using the velocity Verlet method. For the initial condition, we use evenly distributed grid points for the particle positions \(x\) and \(y\) (see the first figure). We use random velocities in the range of (-0.05, 0.05) for both \(v_x\) and \(v_y\). The initial simulation takes a short simulation (t.total \(= 10\)). Time step size \(dt = 0.01\).

set.seed(11)
#force_all(0, x0, a) # if you want to check the forces for the IC
ntot = 25
a = 6.25
grid_x = seq(0, 5, by = 1.25)
grid_y = seq(0, 5, by = 1.25)
mat_xy = as.matrix(expand.grid(grid_x, grid_y))
x0 = c(t(mat_xy))
v0 = runif(n = ntot*2, min = -0.05, max = 0.05)
results_init = velocity_verlet_bc(f = force_all, t0 = 0, x0 = x0, v0 = v0, 
                                       t.total = 10, dt = 0.01, a = a)
nt_all = nrow(results_init)
init_last_step = results_init[nt_all, ]

x0 = init_last_step[2:(2*ntot+1)]
v0 = init_last_step[(2*ntot+2):(4*ntot+1)]

xy0 = matrix(x0, ncol = 2, byrow = TRUE)
plot(xy0[,1], xy0[,2], type = 'p', xlab = "X", ylab = "Y", xlim = c(0, a), ylim = c(0,a))

The plot shows the particles’ positions from the last time point of the initial simulation in a 2D scatter plot.

Final simulation for t.total = 100

After the initial simulation, the system should reach to an “equilibrium” – a more disordered configuration. Now, we continue to perform a longer simulation for t.total \(= 100\). (Before you run this simulation, check the time cost from the initial simulation to estimate the waiting time for this longer simulation.)

results_final = velocity_verlet_bc(f = force_all, t0 = 0, x0 = x0, v0 = v0, 
                                       t.total = 100, dt = 0.01, a = a)
nt_all = nrow(results_final)
final_last_step = results_final[nt_all, ]
x_f = final_last_step[2:(2*ntot+1)]
xy_f = matrix(x_f, ncol = 2, byrow = TRUE)
plot(xy_f[,1], xy_f[,2], type = 'p', xlab = "X", ylab = "Y", xlim = c(0,a), ylim = c(0,a))

Below is the script to generate an animation of this simulation

{r,fig.width = 3.5, fig.height = 4, fig.align='center'}
library(ggplot2)
library(gganimate)
library(gifski)
library(png)
npoints = nrow(results_final)
ind_x = seq(2, 2*ntot, 2)
ind_y = seq(3, 2*ntot+1, 2)

results_md_dframe = data.frame(t = results_final[,1], x = c(results_final[,ind_x]),
                                  y = c(results_final[,ind_y]),
                                  particle = as.character(c(sapply(1:ntot, function(x)
                                    return(rep(x,npoints))))))

p = ggplot(results_md_dframe, aes(x = x, y = y, colour = particle)) + 
  geom_point(aes(colour = particle), show.legend = FALSE, size = 7, alpha = 0.7) + labs(x = "x", y = "y") 
p
anim = p + transition_time(t, range = c(0, 100)) + labs(title = "t: {frame_time}")
anim = animate(anim, fps = 20, nframes = 200, end_pause = 0, rewind = F, renderer = gifski_renderer())
anim_save("md.gif", animation = anim)

Radial distribution function

After the simulation has completed, we use the whole time trajectories of all particles to compute the radial distribution function (RDF), \(g(r)\). \(g(r)\) measures the probability of finding a particle at distance \(r\) from a reference particle, relative to the probability of finding the same particle at \(r\) in a uniform distribution.

\[g(r) = \frac{<\rho(r)>}{\rho_0} = \frac{<\Delta N(r)>}{2\pi r \Delta r\rho_0}. \tag{5}\]

\(\rho_0\) is the particle density:

\[\rho_0 = \frac{N_{tot}}{a^2}.\]

\(\rho(r)\) is the local density of the particles at distance \(r\) away from a reference particle. Thus, \(<\Delta N(r)>\) is the average number of particles at a shell at distance \(r\) and of thickness \(\Delta r\) (red particles) from a reference particle (orange particle).

\(<>\) is the operation of averaging over all particles and all time points during this current simulation. Write a function to compute \(g(r)\) for a series of \(r\) values until \(a/2\) (why?) with an interval \(\Delta r = 0.05\). Here, you will need to count the number of particles in a shell between \(r\) to \(r + \Delta r\) and then perform the averaging. Plot \(g(r)\). The minimum image convention needs to be considered when computing the distances.

par(mar = c(4, 4, 1, 1))   #set plot margins
cal_rdf <- function(ntot, results, a, dr){
  # ntot: number of particles
  # results: simulation output matrix
  # a: box size
  # dr: bin size
  nt_all = nrow(results)
  r_all = seq(dr, a/2, by = dr) 
  nr = length(r_all)
  counts = numeric(nr)
  for(nt in c(1:nt_all)){
    x_all = results[nt, 2:(2*ntot+1)]
    for(i in 1:ntot){
      xi = x_all[2*i-1]
      yi = x_all[2*i]
      for(j in 1:ntot){
        if(i != j){
          xj = x_all[2*j-1]
          yj = x_all[2*j]
          
          dx = xi - xj
          dx = dx - round(dx/a)*a 
          dy = yi - yj
          dy = dy - round(dy/a)*a 
          
          rij = sqrt(dx**2 + dy**2)
          ind = ceiling(rij/dr)
          if(ind <= nr){
            counts[ind] = counts[ind] + 1
          }
        }
      }
    }
  }

  r_mean = r_all-0.5*dr
  rho = counts/nt_all/ntot/2/pi/dr/r_mean
  rho_0 = ntot/a**2
  rdf = rho/rho_0
  return(list(r = r_mean, rdf = rdf))
}

results_rdf = cal_rdf(ntot = ntot, results = results_final, a = a, dr = 0.05)

plot(results_rdf$r, results_rdf$rdf, type = "l", col = 2, xlab = "r", ylab = "rdf")

The shape of the RDF is typical for a liquid system. RDF ~ 0 below \(r = 0.8\) due to strong repulsive interactions from the LJ force. There is a sharp peak near \(r = 1\), representing the first coordinate shell. The second peak at \(r ~ 2\) is supposed to correspond to the second coordinate shell. The third peak is very close to 1. The second and third peaks are much lower than the first peak due to short-range order and long-range disorder.

Please note: