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 (T), and random(0,1) is a random number between 0 and 1.

Here is the code from the Perl Monte Carlo program that implements the above test;

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

    # Pick a random atom
    $atom = int( rand($n_atoms) );

    # save the old coordinates
    @old_coords = ( $coords[$atom][0], $coords[$atom][1],
                    $coords[$atom][2] );

    # Make the move - translate by a delta in each dimension
    $delta_x = $max_translate * ( 2.0 * rand(1.0) - 1.0 );
    $delta_y = $max_translate * ( 2.0 * rand(1.0) - 1.0 );
    $delta_z = $max_translate * ( 2.0 * rand(1.0) - 1.0 );

    $coords[$atom][0] = $coords[$atom][0] + $delta_x;
    $coords[$atom][1] = $coords[$atom][1] + $delta_y;
    $coords[$atom][2] = $coords[$atom][2] + $delta_z;

    # wrap the coordinates back into the box
    $coords[$atom][0] = wrap_into_box($coords[$atom][0], $box_size[0]);
    $coords[$atom][1] = wrap_into_box($coords[$atom][1], $box_size[1]);
    $coords[$atom][2] = wrap_into_box($coords[$atom][2], $box_size[2]);

    # calculate the new energy
    $new_energy = calculate_energy();

    $accept = 0;

    # Automatically accept if the energy goes down
    if ($new_energy <= $old_energy)
    {
        $accept = 1;
    }
    else
    {
        # Now apply the Monte Carlo test - compare
        # exp( -(E_new - E_old) / kT ) >= rand(0,1)
        $x = exp( -($new_energy - $old_energy) / $kT );

        if ($x >= rand(1.0))
        {
            $accept = 1;
        }
        else
        {
            $accept = 0
        }
    }

    if ($accept == 1)
    {
        # accept the move
        $naccept = $naccept + 1;
        $total_energy = $new_energy;
    }
    else
    {
        # reject the move - restore the old coordinates
        $nreject = $nreject + 1;

        $coords[$atom][0] = $old_coords[0];
        $coords[$atom][1] = $old_coords[1];
        $coords[$atom][2] = $old_coords[2];
        $total_energy = $old_energy;
    }

Find this code in metropolis_pl (it is near the bottom of the script). The "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
    $atom = int( rand($n_atoms) );

Note that int( rand($n_atoms) ) is a function 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 coordinates of the chosen atom are saved into the array @old_coords using the line;

    # save the old coordinates
    @old_coords = ( $coords[$atom][0], $coords[$atom][1],
                    $coords[$atom][2] );

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 = $max_translate * ( 2.0 * rand(1.0) - 1.0 );
    $delta_y = $max_translate * ( 2.0 * rand(1.0) - 1.0 );
    $delta_z = $max_translate * ( 2.0 * rand(1.0) - 1.0 );

    $coords[$atom][0] = $coords[$atom][0] + $delta_x;
    $coords[$atom][1] = $coords[$atom][1] + $delta_y;
    $coords[$atom][2] = $coords[$atom][2] + $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][0] = wrap_into_box($coords[$atom][0], $box_size[0]);
    $coords[$atom][1] = wrap_into_box($coords[$atom][1], $box_size[1]);
    $coords[$atom][2] = wrap_into_box($coords[$atom][2], $box_size[2]);

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;

    $accept = 0;

    # Automatically accept if the energy goes down
    if ($new_energy <= $old_energy)
    {
        $accept = 1;
    }
    else
    {
        # Now apply the Monte Carlo test - compare
        # exp( E_new - E_old / kT ) >= rand(0,1)
        $x = exp( -($new_energy - $old_energy) / $kT );

        if ($x >= rand(1.0))
        {
            $accept = 1;
        }
        else
        {
            $accept = 0
        }
    }

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 1, 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 (rand(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 0, which means "fail").

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

    if ($accept == 1)
    {
        # accept the move
        $naccept = $naccept + 1;
        $total_energy = $new_energy;
    }

This code increases the number of accepted moves (stored in $naccept) by one. It also saves the 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 = $nreject + 1;

        $coords[$atom][0] = $old_coords[0];
        $coords[$atom][1] = $old_coords[1];
        $coords[$atom][2] = $old_coords[2];
        $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 total energy of the system is then saved into the variable $total_energy.

That is Metropolis Monte Carlo :-).