Phase Changes

Now you have an understanding of how the Metropolis Monte Carlo program works, it is time to make it do something useful. Try to change the temperature ($temperature, box size (@box_size) and maximum translation ($max_translate) variables so that you can make solid, then gas, then liquid krypton. You will need to change the variables, and then run the metropolis_pl script. For Monte Carlo to be efficient, you need about 40%-60% of the moves to be accepted. If too many moves are rejected, then reduce the value of $max_translate. If too many moves are accepted, then increase $max_translate.

Once you have working parameters, copy them into the C++ program metropolis.cpp. You can set the variables in the lines near the top of this program;

// Set the number of atoms in the box
const int n_atoms = 25;

// Set the number of Monte Carlo moves to perform
const int num_moves = 500000;

// Set the size of the box (in Angstroms)
const double box_size[3] = { 10.0, 10.0, 10.0 };

// The maximum amount that the atom can be translated by
const double max_translate = 0.5;  // angstroms

// Simulation temperature
const double temperature = 298.15;   // kelvin

The variables have the same name as in Perl.

Compile the C++ program and run the program by typing;

$ g++ -O3 metropolis.cpp -o metropolis
$ ./metropolis

Just like with the Perl program, you will get a print out of the moves, energy, and number of accepted and rejected steps, e.g.

1000: 19.876323  264  736
2000: -1.941197  420  1580
3000: -8.801116  565  2435
4000: -10.333685  688  3312
5000: -10.266989  838  4162
6000: -13.508316  976  5024
7000: -3.880514  1141  5859
8000: -11.871016  1300  6700
9000: -9.966295  1455  7545
10000: -7.182279  1618  8382

Copy these into a spreadsheet program (like Excel) and graph the total energy versus move number. Also get the spreadsheet to calculate the acceptance ratio of the Monte Carlo moves. This is the number of accepted moves divided by the total number of moves (e.g. $naccept / ($naccept + $nreject)). Is the acceptance ratio constant during the simulation? Does the simulation need to be equilibrated, or is the energy stable after 500,000 moves?