Running metropolis.py
Now you understand the algorithm, please run the metropolis.py
Monte Carlo program again using the command;
python metropolis.py
You will see printed out the move number, then the total energy (in kcal mol-1), then the number of accepted moves, then the number of rejected moves, printed every 10 Monte Carlo moves. It will look something like this…
10 389339.37706 3 7
20 389245.260273 10 10
30 17523.5326617 18 12
40 5672.37548777 22 18
50 5595.74183933 26 24
60 2581.84585808 32 28
70 2539.87664826 36 34
80 721.890052807 44 36
90 701.680122453 53 37
100 701.710146729 58 42
etc. etc.
Copy the output into a spreadsheet (e.g. by either copying and pasting from the terminal, or by running the simulation again and writing the output into a text file, using `python metropolis.py > metropolis.dat’. You can then load the output text file into a spreadsheet).
In the spreadsheet plot the energy versus move number. Does the energy appear to equilibrate? If it does, calculate the average energy and standard deviation of the average over the equilibrated region. Next, calculate the acceptance ratio of the simulation. This is the ratio of the number of accepted move over the total number of attempted moves. For a good Monte Carlo simulation, this should be over 40%. If you are stuck, take a look here for an example spreadsheet.
As well as printing out the energy and move information, the program will also write a PDB of the coordinates every 100 moves. These output PDB files are called output??????.pdb
, where ??????
is the move number. You can a movie of this simulation using VMD. The process to do this depends on whether you are using Linux/OS X, or if you are using windows…
##If you are using Linux or OS X View these files in VMD. You can make a movie of these files by typing;
vmd output*.pdb
The movie is of a periodic box of krypton. To see the periodic cells, open the Graphical Representations
window in VMD and select the Periodic
tab. Under Select periodic images to draw
click the +X
, -X
, +Y
, -Y
, +Z
and -Z
boxes. This will display those periodic images. Click on the Play
icon in the VMD Main
window to play the movie. What do you see? Are the atoms in the solid, liquid or gaseous state? Does the structure appear to be equilibrated?
##If you are using windows If you are using Windows, then the process to make the movies is a little harder. First, you have to create a VMD animation script. You can do this by then you have to create a VMD animation script. You can do this by downloading the script animatepdb.py. Then type;
python animatepdb.py output*.pdb > movie.vmd
Open VMD, then use the Tcl shell to change into the directory that contains movie.vmd
. In the Tcl shell type;
play movie.vmd
The movie is of a periodic box of krypton. To see the periodic cells, open the Graphical Representations
window in VMD and select the Periodic
tab. Under Select periodic images to draw
click the +X
, -X
, +Y
, -Y
, +Z
and -Z
boxes. This will display those periodic images. Click on the Play
icon in the VMD Main
window to play the movie. What do you see? Are the atoms in the solid, liquid or gaseous state? Does the structure appear to be equilibrated?
Running metropolis.cpp
Now that you have run metropolis.py
, repeat the above using metropolis.cpp
. First, remove the previous output using;
rm output*.pdb
Then, compile metropolis.cpp
using
g++ -O3 metropolis.cpp -o metropolis
Finally, run metropolis.cpp
using
./metropolis
Again, you will see printed out the move number, then the total energy (in kcal mol-1), then the number of accepted moves, then the number of rejected moves, this time printed every 1000 Monte Carlo moves. It will look something like this…
1000 -16.509815 743 257
2000 -15.274366 1516 484
3000 -17.819591 2302 698
4000 -16.643222 3101 899
5000 -17.111819 3848 1152
6000 -17.186787 4593 1407
7000 -16.200351 5371 1629
8000 -15.643928 6132 1868
9000 -16.994626 6904 2096
10000 -16.582194 7635 2365
etc. etc.
In the spreadsheet plot the energy versus move number. Does the energy appear to equilibrate? If it does, calculate the average energy and standard deviation of the average over the equilibrated region. Next, calculate the acceptance ratio of the simulation. This is the ratio of the number of accepted move over the total number of attempted moves. How do these values compare to the short simulation?
As before, visualise the movie of the trajectory using VMD. How does the movie from metropolis.cpp
compare to the one you saw using metropolis.py
?