Volume Moves

The Monte Carlo simulation you have run sampled configurations from the canonical ensemble. This meant that the volume of the periodic box was constant, and it was not possible to specify the pressure acting on the krypton atoms. We live in a constant pressure world, not a constant volume world. To be able to specify pressure, we need to sample from the isothermal-isobaric (NPT) ensemble. This can be achieved by making two changes to the metropolis_pl script;

1. The acceptance test must be changed from the NVT test to the NPT test.
2. Monte Carlo moves that change the volume of the box need to be added.

The first part is easy. Currently, the script uses the canonical MC acceptance test;

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

This must be changed to the isothermal-isobaric test;

        # Now apply the Monte Carlo test - compare
        # exp( -(E_new - E_old + P(V_new - V_old)) / kT
        #             +  N ln (V_new - V_old) ) >= rand(0,1)
        $x = exp( -($new_energy - $old_energy + $pressure * ($V_new - $V_old)) / $kT
                      + $n_atoms * (log($V_new) - log($V_old) ) );

(note that natural logarithm in Perl in "log")

This uses the old and new volumes of the box, that are calculated and stored in the V_old and V_new variables via;

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

    # calculate the old volume of the box
    $V_old = $box_size[0] * $box_size[1] * $box_size[2];

and

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

    # calculate the new volume of the box
    $V_new = $box_size[0] * $box_size[1] * $box_size[2];

(note here that you can see the old and new volumes calculated after the old and new energies, so that you can work out where in the script these lines should be.).

As well as the volumes, the test also requires the pressure, which can go at the top of the script underneath the temperature;

# Simulation temperature
$temperature = 298.15;   # kelvin

# Simulation pressure
$pressure = 1;   # atmospheres

# Convert pressure to internal units (kcal mol-1 A-3)
$pressure = $pressure * 1.458397506863647E-005;

(note that pressure is normally given in units of atmospheres, while the internal units of this program are kcal mol-1 for energy, and angstrom for length. This is why the pressure has to be converted to internal units above).

Once you have made the changes and run the simulation, you should have found that adding in the volume test on its own has not changed the simulation at all. This is because the simulation has only one type of move - the translation of a single atom of krypton. This move does not change the volume of the periodic box, so V_new will always equal V_old. This means that V_new - V_old is always zero, so the NPT acceptance test becomes equal to the NVT acceptance test.

To have an NPT simulation, we need to add in a MC move that changes the volume of the box. This can be achieved by adding in a move that changes the volume of the box by a random amount. This move should only occur once for every $n_atoms moves. This can be achieved using this code;

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

    # save the old box dimensions
    @old_box_size = @box_size;

    # Decide if we are performing an atom move, or a volume move
    if (rand(1.0) <= 1.0 / $n_atoms)
    {
        $is_volume_move = 1;

        # 1 in $n_atoms chance of being here. Perform a volume move
        # by changing the volume for a random amount
        $delta_vol = rand($max_volume_change);

        # Volume is the cube of the box length, so add the cube root
        # of this change onto the box size
        $delta_vol = $delta_vol**(1.0/3.0);

        # Are we making the box bigger or smaller?
        if (rand(1.0) >= 0.5)
        {
            $box_size[0] = $box_size[0] + $delta_vol;
            $box_size[1] = $box_size[1] + $delta_vol;
            $box_size[2] = $box_size[2] + $delta_vol;
        }
        else
        {
            $box_size[0] = $box_size[0] - $delta_vol;
            $box_size[1] = $box_size[1] - $delta_vol;
            $box_size[2] = $box_size[2] - $delta_vol;
        }

        # wrap the atoms back into the box
        wrap_all_atoms();
    }
    else
    {
        $is_volume_move = 0;

        # Make the atom 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]);
    }

We start saving the old coordinates as before, but now we must save the old box size too, using the code;

    # save the old box dimensions
    @old_box_size = @box_size;

We then have to choose between performing an atom move, or a volume move. Each atom should move once between each volume move, so we choose a volume move with a probability of 1 / n_atoms, using the code;

    # Decide if we are performing an atom move, or a volume move
    if (rand(1.0) <= 1.0 / $n_atoms)
    {
        $is_volume_move = 1;
        ...
    }
    else
    {
        $is_volume_move = 0;
        ...
    }

The atom move code is the same as before. The volume move code is straightforward. This adds or subtracts a random volume chosen to have a maximum value of $max_vol_change to each dimension of the box, using the code;

        # 1 in $n_atoms chance of being here. Perform a volume move
        # by changing the volume for a random amount
        $delta_vol = rand($max_volume_change);

        # Volume is the cube of the box length, so add the cube root
        # of this change onto the box size
        $delta_vol = $delta_vol**(1.0/3.0);

        # Are we making the box bigger or smaller?
        if (rand(1.0) >= 0.5)
        {
            $box_size[0] = $box_size[0] + $delta_vol;
            $box_size[1] = $box_size[1] + $delta_vol;
            $box_size[2] = $box_size[2] + $delta_vol;
        }
        else
        {
            $box_size[0] = $box_size[0] - $delta_vol;
            $box_size[1] = $box_size[1] - $delta_vol;
            $box_size[2] = $box_size[2] - $delta_vol;
        }

Because the box has changed size, it may be necessary to wrap some atoms back into the box. This is achieved by calling the wrap_all_atoms() function, which has been added using;

sub wrap_all_atoms
{
     for ($i=0; $i<$n_atoms; $i = $i + 1)
     {
        # wrap the coordinates back into the box
        $coords[$i][0] = wrap_into_box($coords[$i][0], $box_size[0]);
        $coords[$i][1] = wrap_into_box($coords[$i][1], $box_size[1]);
        $coords[$i][2] = wrap_into_box($coords[$i][2], $box_size[2]);
     }      
}                        

The flag $is_volume_move is used to say whether a volume move has been performed. It is equal to one if the move was a volume move, or zero if the move was an atom move. This is used in the code run when a move is accepted or rejected to increment the number of accepted and rejected volume moves, and also to restore the old box dimensions if the move failed, e.g.

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

        if ($is_volume_move){ $nvolaccept = $nvolaccept + 1; }

        $total_energy = $new_energy;
    }
    else
    {
        # reject the move - restore the old coordinates and box size
        $nreject = $nreject + 1;

        if ($is_volume_move){ $nvolreject = $nvolreject + 1; }

        $coords[$atom][0] = $old_coords[0];
        $coords[$atom][1] = $old_coords[1];
        $coords[$atom][2] = $old_coords[2];

        @box_size = @old_box_size;

        if ($is_volume_move){ wrap_all_atoms(); }

        $total_energy = $old_energy;
    }

The parameter for the move is the maximum volume change. This is normally 10% of the number of atoms, e.g. the script should also be modified to add in the maximum volume change;

# The maximum amount that the atom can be translated by
$max_translate = 0.5;  # angstroms
        
# The maximum amount to change the volume - the
# best value is 10% of the number of atoms
$max_volume_change = 0.1 * $n_atoms;   # Angstroms**3

Finally, for completeness, it is a good idea to print out information about the box size and number of accepted and rejected volume moves during a simulation, e.g.

    # print the energy every 10 moves
    if ($move % 10 == 0)
    {
        print "$move: $total_energy  ($naccept : $nreject)  ".
              "($nvolaccept : $nvolreject)\n";

        $vol = $box_size[0] * $box_size[1] * $box_size[2];

        print "    Box size = @box_size. Volume = $vol A^3\n";
    }

2. ADVANCED. Make the same changes to the C++ metropolis.cpp program. Run the program and see what happens during a longer simulation. Note that you can download the program if you get stuck.