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 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


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
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_BT}) & \Delta \mathcal{V} > 0 \end{array} \right. $$

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_BT}) & \text{Accept} \\ \xi > \exp(\frac{-\Delta \mathcal{V}}{k_BT}) & \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^{N_{samples}}_{i} \exp\left(\frac{-\mathcal{V(\{r_i\})}}{k_BT} \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.