Case Study 2: MC of Hard Disks

Download the code hdisk.c by clicking here. This code simulates disks confined to a circle. The Hamiltonian for this system may be expressed as

$$\mathscr{H} = \sum_{i=1}^N H_1({\bf r}_i) + \sum_{i=1}^{N}\sum_{j=i+1}^{N} H_2({\bf r}_i,{\bf r}_j)$$ (85)

where

$$H_1({\bf r}_i) = \left\{\begin{array}{ll} 0 & r_i < R\\ \infty & r_i > R \end{array}\right.$$ (86)

and

$$H_2({\bf r}_i,{\bf r}_j) = \left\{\begin{array}{ll} 0 & \sq... ...eft({\bf r}_i-{\bf r}_j\right)^2} < \sigma \end{array}\right.$$ (87)

$H_1$ acts to keep the particles confined, and $H_2$ prevents them from overlapping. One nice thing about using hard-disk Hamiltonians is that there is never a reason to evaluate a Boltzmann factor. Any trial move that results in an overlap or a particle crossing the boundary gives an ``infinite'' $\Delta\mathscr{U}$, so $e^{-\beta\Delta\mathscr{U}}$ is identically 0 and the trial is unconditionally rejected.

hdisk.c accepts as user input any two of the following three parameters: $R$, the radius of the circle; $rho$, the areal number density of particles (# per square $\sigma$); and $N$, the number of particles. The user may also specify $\delta r$, the scalar displacement, and nCycle, the number of MC cycles, where one cycle is $N$ attempted particle displacements. The only output of the code currently is the acceptance ratio:

$$\left[\begin{array}{c} \mbox{acceptance}\\ \mbox{ratio} \en... ...{number of successful trials}}{\mbox{number of trials}}\right]$$

One important aspect of any MC simulation code is how the particle positions are initialized. Here, it is best to assign initial positions to the particles such that the initial energy is 0 (i.e., there are no overlaps nor particles out of bounds.) Try to figure out how the function init() in the program hdisk.c accomplishes this.

As a suggested further exercise, use hdisk.c to determine a reasonable displacement to achieve a 30% acceptance ratio at a density of 0.5. Compare your results across differently sized systems and runs with different numbers of cycles. For fewer than 10$^6$ cycles, you will have large acceptance ratios because the initial condition is not yet fully destroyed.

Below is a plot of acceptance ratio vs. $\Delta R$ for densities $rho$ of 0.2, 0.4, 0.6, from a simulations of 200 particles. 10,000 cycles were performed for each run. Are your results consistent with this data?

Acceptance ratio vs. $\Delta R$ for various densities in a 200-particle hard-disk system from 100,000-cycle MC simulations.

As a third suggested exercise, consider generating a trial move in the following way:

1. Randomly select a particle $i$.
2. Randomly choose a direction. In 2D, this is an angle $\theta$ chosen uniformly from the interval $[0,2\pi]$. In 3D, this is two angles $\theta$ and $\phi$, where $\cos\theta$ is chosen uniformly from the interval $[-1,1]$, and $\phi$ is chosen uniformly from the interval $[0,2\pi]$.
3. Compute the displacement vector components. In 2D: $dx = \Delta R\cos\theta$ and $dy = \Delta R\sin\theta$. In 3D: $dx = \Delta R\sin\theta\cos\phi$, $dy = \Delta R\sin\theta\sin\phi$, $dz = \Delta R\cos\theta$.
4. Execute the move by adding the displacement vector components to the position of particle $i$.

The main difference of this scheme with the first one is that, in 2D, only one random number must be chosen per displacement attempt. Explore the acceptance ratio vs. $\Delta R$ relationship with this new trial move. Does this scheme generate a 30% acceptance ratio for a higher or lower maximum displacement that does the original scheme?