# Metropolis Monte Carlo

There are many different types of Monte Carlo. The one used in molecular simulation is called “Metropolis Monte Carlo”. This is named after Metropolis, who was one of the five authors of the famous 1953 paper that first introduced the method.

Metropolis Monte Carlo is based on random numbers and has a very simple algorithm.

1. Calculate the energy of the molecular system. Save this as the `old energy`.
2. Randomly choose a molecule in the system to move, and save its coordinates (the `old coordinates`).
3. Randomly move this molecule (e.g. randomly translate or rotate it). This results in the `new coordinates`.
4. Calculate the energy of the molecular system. Save this as the `new energy`.
5. Use the difference between the `new` and `old` energies in a Monte Carlo test. If this test passes, then keep the `new coordinates`, else, restore the `old coordinates`.

For Metropolis Monte Carlo, in a NVT (canonical) ensemble, the Monte Carlo test is;

``exp( -(E_{new} - E_{old}) / kT ) >= random(0,1)``

where `E_{new}` is the `new energy`, `E_{old}` is the `old energy`, `kT` is Boltzmann’s constant (`k`) multiplied by temperature in Kelvin (`T`), and `random(0,1)` is a uniform random number between 0 and 1. If you want to know where this test comes from, take a look at these slides.

Take a look at `metropolis.py` (e.g. by opening it in `nano`). Scroll down until you can find the code that implements the above test (it is copied below here);

``````    # calculate the old energy
old_energy = calculate_energy();

# Pick a random atom (random.randint(x,y) picks a random
# integer between x and y, including x and y)
atom = random.randint(0, n_atoms-1);

# save the old coordinates
old_coords = ( coords[atom], coords[atom],
coords[atom] )

# Make the move - translate by a delta in each dimension
delta_x = random.uniform( -max_translate, max_translate )
delta_y = random.uniform( -max_translate, max_translate )
delta_z = random.uniform( -max_translate, max_translate )

coords[atom] += delta_x
coords[atom] += delta_y
coords[atom] += delta_z

# wrap the coordinates back into the box
coords[atom] = wrap_into_box(coords[atom], box_size)
coords[atom] = wrap_into_box(coords[atom], box_size)
coords[atom] = wrap_into_box(coords[atom], box_size)

# calculate the new energy
new_energy = calculate_energy();

accept = False;

# Automatically accept if the energy goes down
if (new_energy <= old_energy):
accept = True

else:
# Now apply the Monte Carlo test - compare
# exp( -(E_new - E_old) / kT ) >= rand(0,1)
x = math.exp( -(new_energy - old_energy) / kT )

if (x >= random.uniform(0.0,1.0)):
accept = True
else:
accept = False

if accept:
# accept the move
naccept += 1
total_energy = new_energy
else:
# reject the move - restore the old coordinates
nreject += 1

coords[atom] = old_coords
coords[atom] = old_coords
coords[atom] = old_coords
total_energy = old_energy``````

Find this code in `metropolis.py` (it is near the bottom of the script). You should be able to see that `old energy` is calculated using the line;

``````    # calculate the old energy
old_energy = calculate_energy()``````

Note that `calculate_energy()` is a function that calculates the total energy of the molecular system.

An atom to move is chosen at random using the line;

``````    # Pick a random atom (random.randint(x,y) picks a random
# integer between x and y, including x and y)
atom = random.randint(0, n_atoms-1)``````

Note that `random.randint(0, n_atoms-1)` is a function from the python `random` module that generates a random number between 0 and `n_atoms - 1`. This chooses the index of the random atom in the `coords` array, and saves that index in the variable `atom`.

The old coordinates of the chosen atom are saved into the array `old_coords` using the line;

``````    # save the old coordinates
old_coords = ( coords[atom], coords[atom],
coords[atom] )``````

In this program, a random move is just a random translation along the `x`, `y` and `z` axis by a random value between `-max_translate` and `max_translate` Angstroms. This is performed using these lines;

``````    # Make the move - translate by a delta in each dimension
delta_x = random.uniform( -max_translate, max_translate )
delta_y = random.uniform( -max_translate, max_translate )
delta_z = random.uniform( -max_translate, max_translate )

coords[atom] += delta_x
coords[atom] += delta_y
coords[atom] += delta_z``````

The krypton atoms are in a periodic box. This random translation may have moved them outside of the box, so the atoms have to be wrapped back into the box. This is achieved using these lines;

``````    # wrap the coordinates back into the box
coords[atom] = wrap_into_box(coords[atom], box_size)
coords[atom] = wrap_into_box(coords[atom], box_size)
coords[atom] = wrap_into_box(coords[atom], box_size)``````

The `new energy` is calculated using this line, and stored in the variable `new_energy`;

``````    # calculate the new energy
new_energy = calculate_energy()``````

The Metropolis Monte Carlo test is performed using these lines, which perform the test and store whether or not the test passed in the variable `accept`;

``````    accept = False

# Automatically accept if the energy goes down
if (new_energy <= old_energy):
accept = True

else:
# Now apply the Monte Carlo test - compare
# exp( -(E_new - E_old) / kT ) >= rand(0,1)
x = math.exp( -(new_energy - old_energy) / kT )

if (x >= random.uniform(0.0,1.0)):
accept = True
else:
accept = False``````

If the energy goes down (`new_energy` is less than or equal to `old_energy`) then the test is automatically passed (`accept` is set to `True`, which means `pass`). Otherwise, the exponential of the difference of `new_energy` and `old_energy`, divided by `kT` is calculated, and stored in the variable `x`. This is compared to a random number between `0` and `1` (`random.uniform(0.0,1.0)`). If it is greater than or equal to the random number, then the test is passed. Otherwise the test fails (and `accept` is set equal to `False`, which means `fail`).

If the Metropolis Monte Carlo test is passed, then the following code is run;

``````    if accept:
# accept the move
naccept += 1
total_energy = new_energy``````

This code increases the number of accepted moves (stored in `naccept`) by one. It also saves the new total energy of the system into the variable `total_energy`.

If the Metropolis Monte Carlo test is failed, then the following code is run;

``````    else:
# reject the move - restore the old coordinates
nreject += 1

coords[atom] = old_coords
coords[atom] = old_coords
coords[atom] = old_coords
total_energy = old_energy``````

The number of rejected moves (stored in `nreject`) is increased by one. The `old coordinates` of the atom are restored from the `old_coords` array. The old total energy of the system is then saved into the variable `total_energy`.

That is Metropolis Monte Carlo :-)

# PreviousUpNext

Quick Links | Home | Courses | Software | Contact