softmatter:Monte Carlo Simulation (MC)

From NSDL Materials Digital Library Soft Matter Wiki

Monte Carlo simulation is a stochastic method where particle configurations are randomly generated at each step. This differs greatly from methods such as MD ) which involve generating a series of configurations that are sequential in time. Although the MC algorithm of choice may depend heavily on the quantity of interest (usually this quantity is an ensemble average), the Metroplis MC (MMC) algorithm is commonly used. Per the Boltzman distribution law, the ensemble average of a quantity A is given by:

$\sum_i \frac{A_iN_i}{N_i} = \sum_i \frac{A_{i}e^{-\epsilon_i / kT}}{e^{-\epsilon_i / kT}} \qquad \qquad [1]$

where, $ε i$ is the energy of configuration i, T is the system temperature, and the sum is over all states accessible to the system (i.e., phase space). The MMC method is particularly efficient for calculating ensemble averages because configurations are sampled with a frequency that is proportional to their Boltzmann weight (i.e., their contribution to the ensemble average). This is much more efficient than naively generating states at random, where configurations with negligible Boltzmann weights are sampled with the same frequency as those with high Boltzmann weights. It can be proven that imposing the Metroplis constraint on an MC simulation results in the simplification:

$<A > \sim \frac{A_{sum}}{N} \qquad \qquad [2]$

where $A s u m$ is the sum of quantity A over all states sampled in an MMC simulation. While MC simulation can often be carried out using a simple algoirthm, it is a relatively complex approach to molecular simulation compared to MD, particularly for those not well-versed in statistical physics. The most basic MC simulation, which involves sampling with respect to the canonical ( NVT) ensemble, is described below.

Basic algorithm for carrying out NVT Metroplis Monte-Carlo

Initialize particle positions
Perform a trial move by choosing a particle at random and perturbing its position
Accept the trial move with a probability proportional to the ratio of Boltzmann weights of the new configuration vs. the old configuration
Calculate A
Go to step 2)

Two steps in this algorithm require further explanation. First, to carry out step 3) we must derive a formula to calculate the ratio $ρ$ of Bolmann weights of the new configuration vs. the old configuration. This is given by

$\rho = \frac{ e^{-\epsilon_{new} / kT} } { e^{-\epsilon_{old} / kT} } = e^{-(\epsilon_{new} - \epsilon_{old} ) / kT} \qquad \qquad [3]$

We can accept a trial move in proportion to this probability by drawing a random number (R) [0, 1] and accepting the move if

$R < min(1, e^{-(\epsilon_{new} - \epsilon_{old} ) / kT}) \qquad \qquad [4]$

It is important to note that we perform step 4) regardless of whether the trial move is accepted. This reflects the fact that rejecting a trial move is equivalent to accepting the current configuration.