An elegant solution to molecular Monte Carlo simulations that sample from the Boltzmann distribution was published by Metropolis and co-workers of the Los Alamos team in 1953. Instead of accumulating configurations randomly, and then evaluating their probability-weighed contribution to the desired property we can accumulate configurations according to their Boltzmann probability and then take a simple arithmetic average. This is accomplished in the Metropolis Monte Carlo algorithm by first making a random move, then evaluating the Boltzmann probability of such a move and comparing the probability against a random number. If the Boltzmann probability of the move is larger than the random number, the move is accepted; otherwise the system is returned to its original configuration. At the end, a set of configurations is obtained according to Boltzmann statistics and the expectation value of a property is obtained as a simple arithmetic average of property values from individual accepted configurations:
The computer implementation of importance sampling from the Boltzmann distribution is known as the Metropolis Monte Carlo algorithm. We illustrate such implementation by calculating the average position, the root mean square displacement, and the average energy of a classical particle in harmonic potential. Harmonic potential is often used in molecular simulations to describe bond vibrations. Below, a core of the Metropolis Monte Carlo algorithm is implemented in the Python programming language. The complete program can be downloaded here.
#
# Hamiltonian: Harmonic potential
#
def hamiltonian(r):
pot = 0.5 * k * (r - r_eq)**2
return pot
#
# Monte Carlo code
#
acc = 0
mc_dst = 0.0
mc_dst2 = 0.0
mc_ener = 0.0
R = 1.9872156e-3 # Gas constant kcal/mol/degree
beta = 1 / ( R * temp)
r = r_ini
for point in range(1,points):
r_new = r + step * ( random() - 0.5 ) * 2.0 # Move particle randomly
v = hamiltonian(r)
v_new = hamiltonian(r_new)
if v_new < v: # Downhill move always accepted
r = r_new
v = v_new
acc = acc + 1
else: # Uphill move accepted with luck
A_mov = exp( -beta *(v_new - v) )
if A_mov > random():
r = r_new
v = v_new
acc = acc + 1
# Update regardless of acceptance!
mc_dst = mc_dst + r
mc_dst2 = mc_dst2 + r*r
mc_ener = mc_ener + v
One of the user-selectable parameters in the Monte Carlo simulation is the maximum amplitude of the step (the step variable above). The maximum amplitude gets multiplied by a random number, which in the example above can be between -1 and +1; the result determine how much the position of the particle or length of a bond can change in one step. If the step size is small (black) compared to the size of the molecular system the atom moves mainly within its local minimum, and the likelyhood of crossing to the other side of the potential energy barrier is small. Such Monte Carlo search is likely unable to find other conformations. If the step size is too big (red), many Monte Carlo attempts are rejected because large steps are likely to take system into a region of high potential energy (e.g. two atoms come too close). Thus, one needs to determine the optimal step size by trial and error. A good step (green) size should give the acceptance ratio between 30 and 60%.
