B5 COMPUTER SIMULATION OF A LIQUID USING THE MONTE CARLO METHOD
There is no simple theory for the liquid state as there is for the gaseous and solid states. Until the advent of fast computing in the mid 60's there was therefore a whole area of important phenomena for which it was extremely difficult to gain any real insight. Over the past 25 years or so computers have revolutionised this situation completely - although there are still many important problems (e.g moderately dilute aqueous electrolytes) which are only just coming into reach. There are now theoretical approaches, such as perturbation theory and integral equation theory, which can be used to calculate liquid properties with the aid of a computer and which are successful when applied to relatively simple problems, but computer simulation remains a major avenue of approach, and the only viable one in many situations.
There are two basic methods for computer simulation viz: Monte Carlo (MC) and molecular dynamics (MD). Under these headings there are numerous variations on the methodology, with fresh inventions appearing almost by the week. References [1-3] at the end give accounts of the basic procedures. A brief description of MC is given in Appendix 1.
The only input to a simulation is the model of the potential energy
of interaction between the particles in the system. The simple model used
here was originally developed by Lennard-Jones in the 20's. Although it
is an approximate potential it has proved to be remarkably successful,
and has the added advantage of being fairly cheap to incorporate into a
simulation. In this model the interaction potential u between two
spherical particles, separated by a distance r, can be written
where u is positive for repulsion and negative for attraction.
The function passes through a minimum of depth e (the well depth) and is
zero at a separation s, the hard sphere diameter, and at infinity. The
parameters e and s are sufficient to characterise any particular molecule
in this model. It is customary in MC and MD programs to express quantities
in terms of reduced variables. Thus s is used as a convenient unit of length,
and e as a unit of energy.
This experiment (a 'numerical experiment'!) gives you the opportunity to run a Monte Carlo simulation program to calculate the properties of a liquid. The simulation uses a program (SIMAPOS) developed for the simulation of adsorbate mixtures in different adsorbent environments. However, in this exercise you will use an input file that specifies a bulk isotropic fluid. There is no solid adsorbent phase present. The simulations proceed in two stages; first a starting configuration is created and stored, then the main simulation is run from this configuration.
The main simulation is run in the canonical ensemble. The canonical ensemble is one in which the number of particles (N), the volume of the system (V) and the temperature (T) are held constant - it is often simply called the (N,V,T) ensemble. The computer program is designed to simulate a bulk liquid whose molecules, or the atoms comprising the molecules, are modelled as 'Lennard-Jones' (or '12-6') spheres, and to compare the results with experimental data.
The initial simulation is run in the grand canonical ensemble. In this ensemble the chemical potential m is specified, and its conjugate variable N, the number of particles in the simulation, is varied by creating and destroying particles in the simulation box. As for the canonical ensemble, Vand T are held constant.
The source code is written in FORTRAN77 and will run satisfactorily on a suitable PC. You should be able to obtain reasonable results with runs of a few minutes each in length. When you carry out production runs, keep a record of some of the timings. It is interesting to note that when this project was first introduced in 1990, typical run lengths for 64 particles on the College mainframe computer were of the order of 2 hours. Simulation times scale roughly as the square of the number of particles!
To get started you should read the next section and Appendix
2 (which gives a few tips on using the computer) then follow the instructions
in section 3 (RUNNING THE PROGRAM). A brief summary of the way in which
a simulation is set up and run is given in Appendix 1. Once you have gained
confidence in the basic procedures, start on the standard exercises in
section 4.
2. THE PROGRAMS
The compiled executable file (SIMAPOS.EXE) is stored in a directory \mcsim\
The input and output files are summarised below, the parameters can be altered by editing the file (see Appendix 2):
When the program starts it looks for a file named input.dat. The file is divided into four sections:
1. The first two lines specify the temperature and number of molecular species (single species or a mixture of 2).
2. The section, headed ADSORBATES, gives a specification of the molecules in the simulation. Specification for multi-atom and single atom molecular models can be inserted here. For example, methane is treated as a single united atom with interactions given by the Lennard-Jones potential.
3. The section, headed ADSORBENTS, specifies a model for a containing solid. The different models that can be used are identified by an integer KTYPE. In the case of a simple isotropic fluid, to be studied here, KTYPE=0.
4. The final section is headed RUN
PARAMETERS. These include several variables that can be set to vary
the run length etc.
When you open this file under the editor (Appendix 2) it will look something
like this (The line numbers do not appear in the actual file and have been
inserted here for convenience):
Line
1. 222.3 !Temperature in Kelvin
2. 1 !number of adsorbate molecule types in simulation
3. ADSORBATE(S) !Edit in adsorbate specs here:
4. METHANE
5. 10. !ZACT Absolute activity
1. 1,0 ! Number of LJ sites, number of Coulomb sites
7. C1 ! Symbol for atom type
8. 148.2,0.3812,0.,0. ! (epsilon/k)/K, sigma/nm, atom position in the molecule (x,y)
9. ADSORBENT
10. 0 !KTYPE Type of adsorbent.
11. Isotropic fluid !Name; up to 50 characters
12. 2.0 2.0 2.0 !box size nm
13. RUN PARAMETERS
14. 1.5 !CUTOFF in nm
15. 200 !NPLANES Number of distribution function boxes (<=200)
16. 0.5 !HASPCUT Hard sphere cutoff
17. .05 ! DMOVE Maximum distance moved in a trial(nm)
18. .05 ! DPHI Maximum rotation (rad) in a trial (if non spherical adsorbate)
19. 500000,100000 ! NRUN # configs, NREJ # rejects
20. 500,10000,200 ! NDUMP dump frequency, NWRITE screen dump frequency, NDIST distribution functions
21. 2 ! ISIMTYPE, GCMC=1, CEMC=2
22. 1 ! ISTART =0 No restart file. Starts empty; GCMC, or 1 mol; CEMC.
23. 0 ! IDISTEN=1, Calculate energy distribution with single mol.
24. 0 ! IBIAS =0 no bias.
25. 0 !IACCUM=1, results accumulated to accum.dat accum2.dat,
etc
The text after the exclamation marks is annotation and is not read by the program. The names in upper case, following the exclamations, refer to variables used by the program. Many of the options are not required in the present exercise.
The input file is set up to simulate an isotropic fluid of simple Lennard-Jones molecules. The input variables that will be altered at different stages are discussed below. In some cases detailed editing is not required, eg files containing molecule parameters are already present and all that is needed is to replace one section of 'input.dat' by a new section.
Line 1: Temperature in K.
Line 4-8: The specification of the spherical Lennard-Jones molecule:
Line 4: Molecule name (methane in the above example).
Line 5: The absolute activity. This is only used in grand canonical runs (see Introduction and Appendix 1) and is proportional to exp(-m/kT) where m is the chemical potential. This is discussed in more detail below (see 3.1 Creating a startup configuration).
Line 6: The number of Lennard-Jones and electrostatic sites on the molecule.
Some molecules (eg N2, CO2) can be modelled as multi-atom centres, interacting via Lennard-Jones potentials, plus point charges positioned in the molecule. The latter are chosen so as to represent any permanent dipoles or quadrupoles in the molecule (see section 4.2 below). In the simple model above, methane is represented as a single Lennard-Jones centre with no charges.
Line 7: This a symbol consisting of up to 2 characters, that will be recognised by the simple molecular modelling program.
Line 8: There are 4 entries: (i) e/k (in K) and (ii) s (in nm) used by the Lennard-Jones model for the atom in the molecule. (iii), (iv), the location of the atom in molecule fixed coordinates. Since there is only one atom in the above example this is at the centre ie (0,0).
When a multi-centre molecular model is used, more lines are needed to specify all its properties. Look at the file N2KE.dat to see an example of a model for 2-centre nitrogen.
Lines 9-12. These specify the adsorbent.
Since, in this example there is no solid adsorbent phase, the specification is very simple:
Line 9: ADSORBENT is a keyword to help input file editing.
Line 10: KTYPE identifies the adsorbent type. 0 stands for isotropic bulk fluid.
Line 11: A name of up to 50 characters which is sent to the output file.
Line 12: Box size (nm). A cubic box is conventional for an isotropic
fluid.
Lines 13-25.Various parameters that control the run.
Line 13. Header.
Line 14. The long-range cut off (nm). Outside a sphere of this radius, no detailed molecule-molecule interactions are calculated. An approximate long range correction is calculated by the program for isotropic fluids. It does not make sense to have a long range cut off larger than half the box side length (why?), and values larger than this are flagged as an error.
Line 15,16. Not used.
Line 17. Maximum distance that a molecule can move in a simulation step.
Line 18. Maximum rotation (for non spherical molecular models).
Line 19. NRUN is the number of configurations you wish to run. In a production run, the simulation is allowed to settle for an initial number of configurations (NREJ) before averages are calculated.
Line 20. Three integers: NDUMP, NWRITE, NDIST. NDUMP determines the frequency of writing data to the file control.dat. NWRITE determines the frequency of writing to the screen, and NDIST is the number of calls to the singlet distribution routine. The last is not used in this project.
Line 21. ISIMTYPE is a switch that selects either grand canonical ensemble (ISIMTYPE= 1) or canonical ensemble (ISIMTYPE= 2). See the Introduction for the meaning of these terms.
Line 22. ISTART. If this switch is =0, the simulation begins from an empty box. This option only makes sense when GCMC is being run, then molecules are created as the run proceeds. If ISTART=1, the program reads a configuration from the file config.dat. This has the same format as dump.dat, so that a run can be continued from the final configuration of some previous run.
Line 23,24,25 Not used.
3. RUNNING THE PROGRAM
3.1 Creating a startup configuration
The following table gives some data calculated by molecular dynamics
for spherical molecules with a slightly modified Lennard-Jones potential
[4]
TABLE 1 Some state points for cut and shifted Lennard-Jones molecules
|
|
T* | U* |
| 0.835 | 0.902 | -5.85 |
| 0.710 | 0.977 | -4.96 |
| 0.650 | 3.670 | -3.32 |
| 0.550 | 1.350 | -3.70 |
Reduced variables: r* = Ns3/V
T* = kT/e
U* = U/Ne
The reduced variables express density, temperature and internal energy in terms of the Lennard Jones parameters; e, s.
Your first objective is to choose one of these points and study the effects on energy of altering some of the simulation variables. It is immaterial which Lennard-Jones molecule we use, so the methane model, as in the above example of 'input.dat', will do perfectly well. [If this is not already in the file 'input.dat'. It can be edited in by inserting the file C11C.dat using the editor - see Appendix 2.]
To run simulations in the canonical (N,V,T) ensemble, the first
step is to create a suitable starting configuration by filling the simulation
box with the required number of molecules. This is done by running a GCMC
simulation from an empty box until the required density has been reached.
On examining the model input file, you will see that e/k=148.2K and s=0.3812nm
for this methane. Choose a box size of 2.0nm x 2.0nm x 2.0nm to begin with
(Line 12 in the input file) and use the box volume and the equations for
the reduced variables (Table 1) to calculate the temperature in K and the
number of molecules corresponding to your targeted density.
Now set the run parameters to run about 5000 configurations (NRUN on
Line 19) and set the number of rejects to zero (NREJ on line 19). In order
to print at every configuration, set NWRITE=1 on Line 20. Set ISIMTYPE
to 1 (for grand canonical simulation), and ISTART to 0 (start from an empty
box) on Lines 21 and 22.
To run a simulation, open an MS-DOS window, use the DOS command
cd ..
to take you to the top of the C: drive, then
cd mcsim
to change to the directory containing the executable SIMAPOS file. At the prompt type
simapos
You should see a list of numbers flowing up the screen. [If this doesn't happen it may be time to call for help!]. The first column is the configuration number, the second column is the energy. Watch the third column. This is the number of molecules in the simulation box. It should increase gradually as the simulation proceeds. The final number is the number of molecules in the simulation box after 5000 configurations. To have more information about the simulation open the file monitor.dat under the editor. This will confirm the final number you saw on the screen.
To obtain a dump.dat file with the required number of molecules you can do two things:
1. Change the absolute activity (ZACT on Line 5 in the input file). Imagine this as the pressure of an external gas phase. The higher the value of ZACT, the higher the density of the dense phase in the simulation box. Consequently the box tends to fill more rapidly when ZACT is large.
2. Reset the number of configurations (Line 19) so that the simulation terminates with the exact number of molecules required. A little experimentation should produce the desired result.
3.2 Running a canonical ensemble simulation
To make a production run in the canonical ensemble you need to do the following:
1. Reset NRUN to the desired run length, and NREJ to a large enough value to ensure that the simulation has settled (see below). 500,000 and 100,000 are reasonable choices.
2. Change IWRITE to a larger value such as 1000 or 10000 to reduce the number of screen dumps printed.
3. Set the simulation type to ISIMTYPE=2 for canonical ensemble and ISTART =1, so that the program picks up the file config.dat as a startup configuration.
4. To make this file, copy dump.dat to config.dat (see Appendix 2).
Now rerun the simulation program. You should see a line printed every
NWRITE configurations. The screen dumps are really just to reassure you
that the program is running (in fact they are rather a waste of computer
time). The results you need are found in monitor.dat
on completion of the run. Remember that this file is overwritten every
time the program is run. NB. Energies are given in units of kT in
the output file and may need to be converted in order to compare with other
data.
4. STANDARD EXERCISES
4.1 Investigation of the run parameters
The value U/NkT appears in the output file >monitor.dat=. After
conversion to U/Ne it can be compared with the data in Table 1.
Note that this is the configurational part of the internal energy - ie
the part coming from potential energy of molecule-molecule interaction.
The kinetic energy part is easily calculated, see section
4.2 below.
Choose TWO of the (NVT) points in the above table and examine
the effect of the following:
(a) Changing the number of particles. Choose about 3 values up to about 300. An interesting question here is how few particles you can get away with, bearing in mind that the molecular mass contains 6 x 1023.
(b) Length of run. (Line 19) i.e. number of configurations. Remember that a run is split into a settling part and a production part. Make some experiments to decide optimum procedures. It is preferable practice to express this in terms of the number of passes through the system - i.e. each configuration is an attempt to move a single particle, so on average it takes N configurations to make a trial on each of N particles. A pass is the number of configurations divided by the (mean) number of particles in the system.
(c) Long range cutoff. (Line 14) Choose some different values for the largest number of particles which you used in (a), and observe the relative contribution from exact calculation and long range correction.
(d) Maximum step length. (Line 17) Not every trial move in a MC simulation is accepted (see Appendix 2 and references). In fact it is necessary to have some rejections in order to ensure that the system has entropy. Step length influences % acceptance. Check the monitor.dat file.
The file control.dat can be used to help decide how many configurations are needed for convergence - i.e to obtain a stable value of the running mean energy (see Appendix 1 and refs [1-3]. Each line in this file is written every NDUMP configurations (Line 20). If you make NDUMP small this will be a huge file. The columns are: (i) Configuration number. (ii) Mean energy over the last NDUMP configurations. (iii) The mean energy over the whole simulation so far (the running mean energy). You can import this file into a graphics package and plot the partial sum and running means versus configuration number. Such graphs are known as control plots. You should produce some sample control plots and discuss these in your report. Be selective!
N.B. There is no distinct criterion for convergence - this is largely a matter of judgement, see refs [1,2].
The programs include an elementary graphics package. To obtain a picture of the molecules in the simulation box type
picture
at the mcsim> prompt. See Appendix 3.
4.2. An investigation of liquid nitrogen
Liquid nitrogen boils at 77.4K. Its heat of vaporisation is 5.59 kJ
mol-1 [see the Tables at the back of Atkins]. The density of
liquid nitrogen is 28.88 mol dm-3. From the heat of vaporisation,
you can obtain a good estimate for the configuration part of the internal
energy Uconf of the liquid at its boiling point. Uconf
is the quantity calculated in a simulation. It is that part of the internal
energy that comes from the potential energy of interaction between the
molecules. The remainder comes from the kinetic energy, ie
U = Uconf + Ukin = <PE + KE>
Assume that pV for the gas phase =RT, and is approximately
0 for the liquid phase. Uconf for the gas is very small;
it can be found from a knowledge of the second virial coefficient, and
is -0.063 kJ mol-1. Finally assume that Ukin
is the same in the gas phase as in the liquid phase. With this information
you can calculate Uconf for liquid nitrogen at 77.4K.
Note that Uconf is negative unless the molecules are
packed very closely together. To understand this make a sketch (or a plot)
of the Lennard-Jones potential.
In the files provided, there are two models for the nitrogen molecule:
(i) A spherical model, in the file N21C.dat.
This models N2 as a single molecule. The parameters have been
used in several calculations in the literature eg refs [5-7].
(ii) A two-centre model that includes a quadrupole [7]. Each atom in
the molecule is represented by a Lennard-Jones potential and the quadrupole
is modelled by placing 4 partial point charges on the molecule. This model
is to be found in the file N2KE.dat.
Use the editor to replace the methane model with each of these models
in turn, and compare the internal energies obtained from simulation with
those from experiment. Note that you cannot use the single centre config.dat
file for the two centre calculations, since this contains no information
about molecular orientation, and will lead to erroneously large values
of U. You may find snapshots from these two runs instructive - see
Appendix 3
Optional extra It is clear that some adjustment to the
and parameters might result in excellent agreement between experimental
and theoretical internal energies, in either of these models. So why has
this not been done? There are several reasons, and you will have to do
a little reading to dig them out. As an (optional) extra exercise, make
a copy of the file N21C.dat and
see if you can arrive at a parameter set that gives exact agreement with
your experimental value of Uconf for this state point.
APPENDIX 1 - The Monte Carlo method
A simulation is performed by defining a space (the simulation box).
This is usually cubic. The box is filled with molecules by specifying the
location and orientation of each particle. In canonical ensemble Monte
Carlo (CEMC), new configurations are generated by moving one (or more)
molecule(s) to a randomly chosen position and orientation. In the grand
ensemble method (GCMC), the simulation box is filled by running creation/destruction
trials as well as move trials. Any new configuration of lower total energy
is accepted, and new configurations with a higher energy are accepted at
random according to a Boltzmann criterion. All new configurations add to
the averages, whether accepted or rejected. In this way a Boltzmann distribution
is maintained. Since a relatively tiny number of molecules is used in the
simulation, it is necessary to avoid creating spurious edge effects at
the surfaces of the simulation box. This is done by: (a) having >periodic
boundaries= so that a molecule passing through one wall re-enters the simulation
at the opposite wall. (b) Using a >minimum image convention=. The simulation
box is imagined to be surrounded by images of itself. For the purpose of
calculating intermolecular interactions, each molecule in turn is considered
to be at the centre of a sphere of surrounding molecules. The sphere has
the radius of the long range cutoff (see section 2.1)
and the surrounding molecules are always chosen to be the nearest image
from the surrounding boxes. Have a look at refs [1] or [2] to get a clearer
idea of these concepts and other details of simulation technique.
A simulation is run for a large number of trials (configurations), between105
to 107 . Obviously it is not easy to choose "realistic" starting
configuration, and the system takes some time to settle. For example in
this project, the initial energy may be too high or too low. Gradually
the mean value over the whole simulation will become steady, and averages
over a few hundred or thousand configurations will fluctuate about this
"running mean". The simulation can then claim to be "converged". See section
4.1. Note that a simulation may be converged with respect to one property
but not necessarily converged with respect to others. Fluctuation properties
(such as heat capacity) are notoriously difficult to converge.
APPENDIX 2 - file manipulation
Copying files
You will want to keep copies of some of your input and output files, so that you can refer back to them for writing up etc. You should create a new folder under mcsim with a personalised name myfolder ['start' (left mouse button), 'explore', then go to 'file >new']. Copies can be made using explore and dragging and dropping the file into your folder. Then click on the file name and rename. For example a set of input files might be named inp001.dat, inp002.dat etc (assuming you won't need more than 1000!). Copying can also be done under DOS using
copy filename1 filename2
eg:
copy dump.dat config.dat
This is useful in batch files where you can keep a copy of the standard
output file, or copy an input file to the standard name input.dat,
that is recognised by the program.
Editing
To edit the input file, go to 'start', 'programs' and click on 'MS-DOS prompt'. A window opens in the WINDOWS directory, type
edit
and click on 'file', 'open'. Use the right hand window to search for the required directory [ mcsim], then click on the required file [eg input.dat] in the left hand window. Alternatively you can transfer to the directory you wish to work in by using
cd ..
to get to the top directory and
cd mcsim
to transfer to the simulation directory. Then type
edit
The files present can be examined using 'explore' [go to 'start' and click left mouse button] or by typing
dir
at the MS-DOS prompt (or dir/p to prevent file names running out of the window).
In part B you need to replace a section of the input file by a new section.
One way to do this is to open both files under the editor and copy and
paste between them.
Killing a job
If you want to stop a job before it is completed click on the simulation window and use <CTRL-C>.
Batch files
Running in batch mode means that you can stack up jobs and run them overnight, whilst you are down at the pub etc. An example of a batch file is given in mcrun.bat. More elaborate batch files can be constructed that pick up data files one after the other and copy output to a corresponding series of report files. You may wish to experiment with using batch files.
A batch file can be run from 'start', 'run', or by typing the name of the batch file (e.g. mcrun) at the mcsim prompt in an MS-DOS window.
To check that a batch file is performing correctly, you can watch a few configurations, then
type <CTRL-C> to abort the simulation and skip to the next batch
command.
Leaving MS-DOS
To close the MS-DOS window, type
exit
at the prompt.
APPENDIX 3 - Molecular pictures
The molecular modeller can be run by typing
picture
at the MS-DOS prompt. Then go to 'file', 'open', you should see the names of two .mol files in the right hand window. Click on mcsnap2.mol. (the other file is the same snapshot without the box). A side view of the simulation box containing the molecules should appear. Use the up and down arrows, or the mouse to rotate the box. Under 'style' there are options to display the picture as 'space fill' or 'shade fill'. For the latter delete the contents of the box next to 'screen fill', and type in 1000, then hit return. The box is not drawn with the last two options, but you will be able to see its corners as tiny black points. If you prefer a picture without these, use the file mcsnap3.mol.
[1] Allen, M.P. and Tildesley, D.J. Computer Simulation of Liquids,OUP (1984).
[3] Frankel, D and Smit, B., Understanding computer simulation, Academic Press, 1996.
[3] Nicholson, D. and Parsonage, N.G Computer Simulation and the Statistical Mechanics of Adsorption,1982, Academic Press.
[4] Nicolas, J. B., Gubbins, K., Street, W and Tildesley, D. J. Molecular Physics,37,1429,(1979).
[5] Hirschfelder, J.,Curtiss, C.F. and Bird, R.B., Molecular Theory of Gases and Liquids, Wiley (1961).
[6] Lastoskie, C., Gubbins, K. E. and Quirke, N. APore size heterogeneity and the carbon slit pore: A density functional theory model.@, Langmuir 9, 2693 (1993).
[7] Kaneko, K., Cracknell, R. F. and Nicholson, D. ANitrogen adsorption
in slit pores at ambient temperatures: Comparison of simulation and experiment.@,
Langmuir 10,4606 (1994).
David Nicholson, October 1998. d.nicholson@ic.ac.uk