Particles in a box revisited

We again model a box of particles with pairwise Lennard-Jones potential.

For the pairwise interactions:

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

The corresponding potential function is

\[U(r) = \int{-F(r)dr} = \frac{1}{12r^{12}}-\frac{1}{6r^{6}}. \tag{2}\] We can then compute the total potential for the particles in the box. We consider the periodic boundary condition, so the minimum image convention is also applied.

# pairwise force
potential <- function(r2){
  # r2: r square (not r to save a square root operation)
  r4 = r2*r2
  r6 = r4*r2
  r12 = r6*r6
  u = 1/12/r12 - 1/6/r6
  return (u)
}

r_all = seq(0, 4, by = 0.01)
plot(r_all, sapply(r_all, potential), type = "l", col = 1, 
     xlab = "r", ylab = "U", xlim = c(0,4), ylim = c(-0.2,0.3))
abline(h = 0, lty = 3, col = 2)

# potential between particles i and j
potential_ij <- function(xi, xj, a){
  # xi: position of particle i (a vector of size 2)
  # xj: position of particle j (a vector of size 2)
  # a: box size
  
  dx = xi - xj
  dx = dx - round(dx/a)*a  # minimum image convention 
  
  rij2 = sum(dx*dx)
  return(potential(rij2))
}

potential_tot <- function(x_all, a){
  # Ntot: number of particles
  # x_all: a vector of 2 * Ntot
  # a: box size
  Ntot = length(x_all)/2
  p_tot = 0
  for (i in 1:(Ntot-1)){
    xi = x_all[(2*i-1):(2*i)]
    for (j in (i+1):Ntot){
      xj = x_all[(2*j-1):(2*j)]
      p_ij = potential_ij(xi,xj, a)
      p_tot = p_tot + p_ij
    }
  }
  return(p_tot)
}

Interestingly, we don’t need to consider the kinetic energy in the Boltzmann distribution, as the kinetic energy is only velocity dependent and can be integrated out in the partition function (we will not discuss this Statistical Mechanics topic in details). Thus, in the Metropolis-Hastings algorithm, we just need to consider the particle positions and the potential energy and just omit the velocities and the kinetic energy.

Local move & energy update

To efficiently sample another configuration of the system, we propose a local move, which involves a randomly selected particle and a small displacement. Such a local move would allow a small change in the total energy; therefore, the proposed move is more likely to be accepted.

Below shows the calculation of the energy differences after a proposed local move. Only potential terms related to the moved particle \(i\) are computed.

dp <- function(x_all, a, x_new, i){
  # x_all: a vector f 2 *Ntot
  # a: box size
  # x_new: proposed position of a particle i
  # i: index of the moved particle
  xi = x_all[(2*i-1):(2*i)]
  
  Ntot = length(x_all)/2
  p_old = 0
  p_new = 0
  for (j in 1:Ntot){
    if(i != j){
      xj = x_all[(2*j-1):(2*j)]
      p_old = p_old + potential_ij(xi, xj, a)
      p_new = p_new + potential_ij(x_new, xj, a)
    }
  }
  return(p_new - p_old)
}

Q: How much improvement is this method in terms of computational efficiency? Is there a way to further improve the efficiency?

Metroplis-Hastings algorithm

Specifying the initial configuration. Note that, in MCMC, only positions are needed. We also make the box size slightly larger (why?).

set.seed(10)
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))

potential_tot(x0, a) # initial potential
## [1] -2.24004

Use Metropolis-Hastings algorithm to sample particle configurations.

metropolis_particle <- function(ntot, cal_e, cal_de, x0, dxmax, nstep, temp, a){
  # ntot: number of particles
  # cal_e: name of the function to compute the total energy
  # cal_de: name of the function to compute the energy differences
  # x0: initial configuration (a vector of size 2*ntot)
  # dxmax: maximum displacement in a proposed move (a vector of size 2)
  # nstep: number of MCMC steps
  # temp: temperature (scaled)
  # a: box size
  
  x = matrix(0, nrow = (nstep + 1), ncol = 2*ntot)  # positions
  e = numeric(nstep + 1) # total energy
  
  x[1, ] = x0
  e[1] = cal_e(x0, a) 
  num_accept = 0
  for(i in 1:nstep){
    k = sample(x = ntot, size = 1) # randomly select a particle
    dx = runif(n = 2, min = -dxmax, max = dxmax)
    x_new = x[i, (2*k-1):(2*k)] + dx
    x_new = x_new - floor(x_new/a) * a  # periodic boundary condition
    
    de = cal_de(x[i,], a, x_new, k)
    if((de < 0) || (runif(1) <  exp(-de/temp))){
      x[i+1,] = x[i,]
      x[i+1, (2*k-1):(2*k)] = x_new
      e[i+1] = e[i] + de
      num_accept = num_accept + 1
    }else{
      x[i+1,] = x[i,]
      e[i+1] = e[i]
    }
  }
  print(paste0("Acceptance rate ", num_accept/nstep))
  return(list(x = x, e = e))
}

MCMC sampling - the initial phase. This is needed as the initial configuration has a higher energy state for \(T = 0.05\). The parameter dxmax was selected to make the acceptance rate around 40% - 50%.

set.seed(10)
dxmax = c(0.25, 0.25)
temp = 0.05
nstep = 5*10^3

results_init = metropolis_particle(ntot = ntot, cal_e = potential_tot, cal_de = dp, x0 = x0, 
                    dxmax = dxmax, nstep = nstep, temp = temp, a = a)
## [1] "Acceptance rate 0.4512"
plot(0:nstep, results_init$e, type = "l", col = 2, lty = 1, xlab = "Step", ylab = "Total Energy", xlim = c(0,nstep))

x1 = results_init$x[nstep + 1,]

A particle configuration after the initial phase of the MCMC sampling.

xy_1 = matrix(x1, ncol = 2, byrow = TRUE)
plot(xy_1[,1], xy_1[,2], type = 'p', xlab = "X", ylab = "Y")

A second phase of MCMC sampling, prepared for the RDF calculation later.

nstep = 1*10^4

results_final = metropolis_particle(ntot = ntot, cal_e = potential_tot, cal_de = dp, x0 = x1, 
                    dxmax = dxmax, nstep = nstep, temp = temp, a = a)
## [1] "Acceptance rate 0.4079"
plot(0:nstep, results_final$e, type = "l", col = 2, lty = 1, xlab = "Step", ylab = "Total Energy", xlim = c(0,nstep))

Q: Can we propose a move that involves displacements of multiple particles?

Radial distribution function

The RDF is computed by an ensemble average from the second phase of the MCMC sampling.

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, ]
    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$x, a = a, dr = 0.05)

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

Q: Compare the results with those from an MD simulation. Which method is more efficient?