Molecular Simulation/Monte Carlo Methods

Monte Carlo Methods are stochastic techniques that use random numbers to sample conformation space. In this method, a random conformation is generated then it would be determined whether to reject or accept it.

A Metropolis Monte Carlo simulation of Liquid Argon at 98 K

A Monte Carlo simulation starts from a given conformation, then random numbers will generate a new trial conformation. This trial conformation will be determined whether to be accepted or rejected. If it is accepted, this conformation will become the current conformation and the next iteration will start from this conformation.

Monte Carlo Integration

Using Monte Carlo method for evaluating $\int _{0}^{1}x^{2}dx$. The random samples will be accepted with a probability of $0.{\overline {3}}$

Computing a definite integral is a simple illustrative example of using Monte Carlo method. Consider the simple function $f(x)=x^{2}$ and its definite integral from 0 to 1:

$I=\int _{0}^{1}f(x)dx$

In this case, the conformation space is two dimensional and random x and y values will be generated to make the new trial conformation. Only the trial points (a pair of x and y values) that are below the graph can be accepted. The number of accepted conformations divided by the total number of samples represents the probability for samples to be accepted. Multiplying this probability by the volume of the sampling space ($V=1\times 1$) is an approximation of the correct value of the integral.

$I\approx Q_{N}=V{\frac {1}{N}}\sum _{i=1}^{N}f(x_{i})$

A large number of samples will result a closer value to the expected value.

$\lim _{N\to \infty }Q_{N}=I$

This method is usually used for higher-dimensional integrals.

Metropolis Monte Carlo Simulation

The schematic of the Metropolis Monte Carlo Algorithm. Moves that decrease the potential energy of the system ($\Delta {\mathcal {V}}\leq 0$) are always accepted. Moves that increase the potential energy ($\Delta {\mathcal {V}}>0$) are accepted with a probability determined by the Boltzmann distribution.

In molecular simulations, Metropolis Monte Carlo Algorithm is usually used to accept or reject the trial conformation. If the energy of the trial conformation is lower than or equal to the current energy, it will always be accepted. If the energy of the trial conformation is higher than that of the current energy, then it will be accepted with a probability determined by the Boltzmann distribution.

$P_{acc}=\left\{{\begin{array}{lr}1&\Delta {\mathcal {V}}\leq 0\\\exp({\frac {-\Delta {\mathcal {V}}}{k_{B}T}})&\Delta {\mathcal {V}}>0\end{array}}\right.$

At higher temperatures, high energy increases are more likely to be accepted. Conversely, at lower temperatures there is a greater probability of rejecting the state.

It means that, if $\Delta {\mathcal {V}}>0$, a random number ($\xi$) will determine whether to accept or reject this conformation.

$\left\{{\begin{array}{lr}\xi \leq \exp({\frac {-\Delta {\mathcal {V}}}{k_{B}T}})&{\text{Accept}}\\\xi >\exp({\frac {-\Delta {\mathcal {V}}}{k_{B}T}})&{\text{Reject}}\end{array}}\right.$

Metropolis algorithm satisfies microscopic reversibility, meaning that the probability of being in a state and going to the next state is the same as being the new state and going back to the previous state.

Sampling

Configurational integral requires integration of Boltzmann distribution over all spatial variables. To employ Monte Carlo integration to estimate configurational integral (Z), a large number of configurations are randomly generated and Boltzmann weight is evaluated only at these points. With enough samples, configurational integral can be approximated.

$Z\approx {\frac {V^{N}}{N_{samples}}}\sum _{i}^{N_{samples}}\exp \left({\frac {-{\mathcal {V(\{r_{i}\})}}}{k_{B}T}}\right)$

But direct Monte Carlo sampling is extremely inefficient because a huge number of states must be sampled. Besides, there is an equal probability to sample high energy states as low energy states, but by Boltzmann distribution system is likely to be in a low energy state.

Metropolis algorithm changes the acceptance probability, so that the system samples low energy states preferentially.

Limitations of Monte Carlo Simulations

With Monte Carlo simulations, only an equilibrium distribution of states can be sampled. Moves are random and non-physical, so trajectories only corresponds to sequence of Monte Carlo moves and there is no temporal information, so dynamical properties like diffusion coefficients cannot be calculated from a Monte Carlo trajectory.

Equilibration

Metropolis Monte Carlo cycle favor transitions to lower energy and reach more probable structures. Many iterations of the metropolis Monte Carlo will move the system to an equilibrium state. Multiple properties can be plotted to determine when equilibration is complete.

Moving to equilibrium can be seen in a slow drift, equilibration is not complete until drift has stopped and properties fluctuate around average value.

Equilibrium properties

At equilibrium, the system is not undergoing a systematic change. There is no large scale gradient in physical properties of matter. If small perturbation is made to the system, the system will return to its equilibrium state. In equilibrium there is no momentum gradient, concentration gradient and temperature gradient in the system.