View source for Rep:Mod:phys3.1

===Molecular Dynamics Simulation of liquids===

Dr. Fernando Bresme, Room 165, C2, Department of Chemistry Imperial College London

===Introduction===
Molecular Dynamics computer simulation is a theoretical approach to investigate the structure and dynamics of atoms, molecules and molecular assemblies, using microscopic information based on intermolecular potentials. In this tutorial you will investigate the structural and thermodynamic properties of liquid phases using a general purpose molecular dynamics code called DL_POLY. DL_POLY is a parallel simulation package developed at Daresbury Laboratory. Most of the documentation relevant to this package can be found in the [http://www.cse.scitech.ac.uk/ccg/software/DL_POLY/MANUALS/USRMAN2.19.pdf manual].
 
The manual contains information on the principles and implementation of the Molecular Dynamics algorithms, as well as practical information on the DL_POLY error messages. If something goes wrong during the computation the code will produce an error message, you are advice to refer for information to the DL_POLY manual.
 
The tutorial focuses on molecular dynamics simulations of water and water solutions. Water is the most abundant compound on the surface of Earth. It is the most important solvent, providing the environment where many physical processes occur, from protein folding to electron transfer processes at interfaces, mineral dissolution and corrosion processes. Water exhibits anomalous behavior in almost any property one can measure (see for instance <ref> Eisenberg and W. Kauzmann, The Structure and Properties of Water, Clarendon Press – Oxford 1969 (new edition 2005)</ref>); high dielectric constant, high heat capacity, high melting, boiling and critical temperatures. Also, the thermal conductivity exhibits an anomalous decrease with temperature, the solid phase diagram contains a very large number of polymorphs, it features a density maximum at 4 C, it gets [http://physics.aps.org/synopsis-for/10.1103/PhysRevLett.101.020602 polarized] under temperature gradients... etc
 
The computer simulation of liquid water is far from trivial. Many theoretical [http://www.lsbu.ac.uk/water/models.html models] have been developed over the years). In this tutorial you we will be using the SPC/E model<ref>H.J.C. Berendsen, J.R. Grigera and T.P. Straatsma, J. Phys. Chem., 91, 6269 (1987).</ref>, which has shown to be very successful reproducing the properties of liquid water. Using this model we will investigate the structural and thermodynamic properties of water and water solutions at ambient conditions and near critical conditions. 

<div class="references-small"></div>
<references/>

===Before you do anything else===
Before you run DL_POLY  it is a good idea to copy all you need in your local directory.

In your login shell write the following: '''cp -r /opt/chemistry/dlpoly .''' 

and press RETURN.

In the directory '''dlpoly''' you should have two subdirectories: '''execute''' and '''java'''.
We recommend you to run all the exercises in the '''execute''' directory.

To run DL_POLY write '''./DLPOLY.X''' in the '''execute''' directory.

DL_POLY comes with a java applet called GUI that we will use to create input files and analyse results. To run the GUI 
write '''./gui''' in the '''execute''' directory.

===Molecular Dynamics Simulations using DL_POLY ===

To perform a computer simulation with DL_POLY you need the executable, DLPOLY.X and three files; CONTROL, CONFIG and FIELD.  The CONTROL file defines the simulation conditions. It makes use of keywords to define the temperature, the ensemble, timestep, ...etc. An example of a typical CONTROL file is given below:
[[Image:FBresme-FIG1.jpg]]

The example above corresponds to a simulation in the microcanonical ensemble (NVE). The temperature is set to 300 K, and the timestep, used to integrate the equations of motion, is 0.0005 ps, The forces between particles are cut-off at 7.6 angstroms. This means that the forces are set to zero beyond this distance. In this exercise you will be using the Ewald method<ref>D. Frenkel and B. Smit, Understanding Molecular Simulation, From algorithms to applications, Academic Press,  1996. </ref>  (Ewald precision 1e-6 defines the level of accuracy of the Ewald sum), to compute the coulombic interactions between the partial charges  in the water SPC/E model. 

To write a comment in the CONTROL file you may use the hash symbol in front of the text (see example above).  

The CONFIG file contains the dimensions of the unit cell, defines the periodic boundary conditions, atomic labels and coordinates of the atoms we want to simulate. An example of a CONFIG file is given below:

[[Image:config.gif]]


For each atom the first row represents the coordinates in Angstrom, and the second and third rows the velocities and the forces respectively. The CONFIG file  can be created using the DL_POLY GUI. 

Finally, the FIELD file contains the force field information, i.e., it defines the intermolecular interactions. 

[[Image:FBresme-FIG3.gif]]

The example above defines the interactions for the water SPC/E model. In this example we have considered 256 water molecules. Each molecule is modelled as a rigid triangular frame. This is done by using the CONSTRAINTS command. The SPC/E model has three rigid bonds; two OH = 1 Ǻ, and one HH = 1.632880 Ǻ, which corresponds to an HOH angle of 109 degrees (see figure below). 

[[Image:FBresme-FIG4.gif]]

In addition to the molecular geometry we define in the force field the atomic mass of the atoms, oxygen and hydrogen, and the partial charges on each atom. Finally the VDW 3 command in the FIELD file defines the intermolecular interactions. These are modelled using a Lennard-Jones potential (key “lj” in the FIELD file):

[[Image:FBresme-FIG5.gif]]

where ε is the potential depth and σ the atom diameter (in Ǻ). The units for the potential depth, ε, are defined in the top of the FIELD file, “UNITS kJ”, meaning kJ/mol. In the SPC/E model the van der Waals interactions only involve oxygen pairs, the oxygen-hydrogen and hydrogen-hydrogen interactions are set to zero. Oxygen-hydrogen and hydrogen-hydrogen interact between them through coulombic interactions only. Because the molecules are considered rigid there are not intramolecular interactions.

Once the CONTROL, CONFIG and FIELD files have been define, DL_POLY can be executed using the command DLPOLY.X

==== Generating CONTROL, CONFIG, and FIELD files using the Java Graphical User Interface====

The Java Graphical User interface (GUI) provides an easy way to generate the files needed to run a molecular dynamics simulation. To start the GUI write ./gui (in the execute directory) in the Linux command line. Two windows will appear, “Monitor” and “DL_POLY GUI” (see figure below). 

[[Image:FBresme-FIG6.gif]]

To generate an initial configuration of liquid water do the following: under the “FileMaker” option select “CONFIG” and then select “Lattice”. A new window will appear. Set the A, B and C vectors to: A=(22,0,0), B=(0,22,0), C=(0,0,22) Ǻ. Click “Make”. A cubic box of length 22 Ǻ will be created. Now we will fill the box with a pre-equilibrated configuration of water. Under the “FileMaker” option, select “Tools” and then “Add Water”. A new window will appear. Set the “Min water-water distance” to 1.35 Ǻ, and click “Make”. A new configuration should appear (see below).

[[Image:FBresme-FIG7.gif]]

The GUI program will also generate a file called “CFGH20.0”. This is the CONFIG file we will use  in our first simulation of water.

Let us generate now templates for the “CONTROL” and “FIELD” files using the “FileMaker” option.  Select “FIELD” and “BLANK” to generate the FIELD file. This will generate a file called FLDBLK.0 that we will edit later to define the SPC/E water model. To create the CONTROL file select “CONTROL” under the “FileMaker” option. A new window will appear. Set the “time step” to 0.002 ps, the “Cutoff” to 9 Ǻ, the “Verlet shell width” to 1 Ǻ, the “VDW cut off” (van der waals cutoff) to 9 Ǻ.  Select the “EWALD” method to compute the Coulombic forces. Click “Make”. The GUI will create the file “CNTROL.0”.   Before we run the simulation you will need to rename<ref>To rename a file within Linux use the command “mv”: “mv oldfilename newfilename”</ref> the files CFGH20.0, FLDBLK.0 and CNTROL.0, to CONTROL, CONFIG and FIELD respectively .

You will notice in the picture above that some molecules, in the edges of the simulation box, seem broken. This is because in computer simulations we use [http://en.wikipedia.org/wiki/Periodic_boundary_conditions Periodic Boundary Conditions](PBC). PBC provide a route to simulate a system that resembles an infinite system.

====Example:  Simulation of water at ambient temperature ====

Using the files generated above we can simulate water at 298 K in the microcanonical ensemble (NVE). First, edit  your CONTROL file so that you have the same definitions given in the example below. Delete any lines not appearing in the example.

[[Image:FBresme-FIG8.gif]]

Edit your FIELD file and define the partial charges, bond distances and pair potentials of the water SPC/E model. Your FIELD file should look like:

[[Image:FBresme-FIG9.gif]]

Check that the definitions in your FIELD file are consistent with the SPC/E model parameters published by Berendsen et al <ref>H.J.C. Berendsen, J.R. Grigera and T.P. Straatsma, J. Phys. Chem., 91, 6269 (1987)</ref>.

'''IMPORTANT''':  The input file, FIELD, is read with a fixed format. This means that the columns defining the potential, charge etc, have to be separated by a pre-defined number of spaces. If you are getting error messages when trying to execute DL_POLY, check first whether you are using the right format in your FIELD file. To be on the safe side you are recommended to use the FIELD file generated by the GUI as template for all your runs. 

Now you can execute DLPOLY.X, the code will generate several files:

'''OUTPUT'''

[[Image:FBresme-FIG10.gif]]

The meaning of the labels appearing in the OUTPUT file can be found in the DL_POLY 
[http://www.cse.scitech.ac.uk/ccg/software/DL_POLY/MANUALS/USRMAN2.19.pdf manual]. The properties that  we need to consider in this tutorial are listed below:

'''eng_tot''' 		 total internal energy of the system

'''temp_tot''' 		 system temperature

'''eng_cfg''' 		 configurational energy of the system

'''volume'''         system volume

'''press''' 		          pressure

In addition to the '''OUTPUT''' file, DL_POLY generates:

o	'''STATIS''', contains a summary of statistical data.

o	'''REVCON''', contains a sample of the final configuration.

The '''STATIS''' file can be used to analyze the simulation trajectory. The DL_POLY GUI provides an interface to do this in an easy way. Open the GUI and select “Analysis / Tools / Statistics”. A new window will appear. Now you can plot the system properties as a function of the simulation time. Select “E-TOT” and click “RUN”, to plot the variation of the total energy with time (see below).

[[Image:FBresme-FIG11.gif]]


Try plotting other properties, temperature, configurational energy, van der Waals and coulombic energies. Note that each time you select one of these properties the GUI creates a file with the name STATS0.XY, STATS1.XY, …, which contains the data represented in the graphs.

The GUI can be used to draw the final configuration, REVCON.  Just select “Analysis/ Display/ REVCON”. The configuration can be rotated, and bonds can be represented using the options appearing on the right panel of the DL_POLY GUI.

====What time step should you use in the simulations?====

One of the main variables in molecular dynamics simulations is the timestep. The timestep determines the accuracy in the numerical integration of the equations of motion. A too small time step provides good accuracy in the integration of the equations of motion. The disadvantage is that very long simulations times are needed to ensure adequate sampling of the system. A too large time step enables longer simulation times but introduces inaccuracies in the integration of the equations of motion that result in incorrect simulations. As a general rule an adequate timestep ensures that the constraints in an ensemble, eg NVE, are conserved. For instance, in the microcanonical ensemble, NVE, (once the equilibration period has finished) the total energy should be conserved, or in other words, it should be constant. 

To explore the impact that the timestep has on molecular dynamics perform two simulations in the microcanonical ensemble, NVE, using as timesteps 0.002 and 0.01 picoseconds. You may consider a trajectory of  1 ps, with 0.5 ps for the equilibration period. Plot the total energy as a function of time for the two simulations. Is 0.01 ps a reasonable choice for the simulation of SPC/E water?

'''Important''': each time you run a simulation DL_POLY  overwrites any existing STATIS, OUTPUT and REVCON files;  to keep the results from previous exercises, remember to rename the files at the end of each run.

=== Exercises ===

In the following exercises you will be investigating the properties of water at different thermodynamic conditions, from ambient conditions, to high pressures and near critical conditions. The thermodynamic states we will consider are indicated by two circles (p≈1 atm, and T≈647.3 K) in the phase diagram given below <ref> (P.W. Atkins and de Paula, Oxford University Press, 7th Edition)</ref>

[[Image:FBresme-FIG12.jpg]]



<div class="references-small"></div>
<references/>

====Exercise  1: Internal energy and density of SPC/E water at ambient conditions====

In the previous example we have considered simulations of water in the microcanonical ensemble (NVE). Very often it is convenient to perform the simulation in an alternative ensemble where we can set the pressure and temperature of the system. Such ensemble is called isothermal-isobaric (NpT), and is the ensemble that corresponds to experimental conditions. In this exercise you are asked to perform a molecular dynamics simulation of SPC/E water at 298 K and 1 bar pressure in the NpT ensemble. To do this you need to introduce a minor change in the CONTROL file you have used before. Just substitute 
  
'''ensemble nve''' by  '''ensemble npt hoover 0.1 0.1''' 

and add  to the CONTROL file : '''press 0.001''' (pressure units are in kbars)

The NpT algorithm works by coupling your system to a virtual thermostat and barostat that define the temperature and pressure of the system. The two parameters after the ensemble definition, “0.1”, are time constants in picosedonds, which specify the relaxation of the thermostat and barostat. “0.1 ps” is an adequate value for the simulation of SPC/E water.
 
(a) Perform a simulation of SPC/E water in the NpT ensemble at T=298 K and p=1 bar. Analyse the equilibration of the system by plotting the time variation of the volume with time. To perform this simulation you should use typically 4000 steps, (2000 for equilibration).

(b)Calculate the configurational energy of water at normal conditions. Compare your result with the data given in references <ref>H.J.C. Berendsen, J.R. Grigera and T.P. Straatsma, J. Phys. Chem., 91, 6269 (1987)</ref> and  <ref>E. Sanz et. al., Phys. Rev. Lett., 92, 255701 (2004)</ref>.

'''Note''': DL_POLY reports the total energy for the system of N molecules. To obtain the energy per molecule you need to divide your result by the number of molecules in your simulation.

(c) What is the density predicted by the  SPC/E water model at 298 K and 1 bar pressure? Compare your result with the experimental one.

<div class="references-small"></div>
<references/>

====Exercise 2: Liquid water structure and hydrogen bonding====

The structure of gas, liquid and solid phases, can be analysed through the so called pair correlation function g(R). The pair correlation function quantifies the probability that a molecule will be found at a distance R from another molecule. Mathematically it is defined as:

[[Image:FBresme-FIG13.jpg]]

where the numerator is the local density of water molecules at a distance R from a given molecule, and the denominator is the average density of the liquid at that distance. This is the density that the liquid would have in an ideal gas, i.e., in a system without interactions, this is why we have used the "id" subscript. <ref>P. Atkins and J. de Paula, Physical Chemistry, Oxford University Press, 7th Edition.</ref>,<ref>D. Frenkel and B. Smit, Understanding Molecular Simulation, From algorithms to applications, Academic Press,  1996</ref>.

The figure below shows a typical radial distribution function for a dense liquid. The oscillations indicate the formation of shells of molecules around a given molecule (see sketch). There are regions with a higher density than the average, g(R) > 1, which are separated by low density regions  g(R) < 1. At long interparticle distances g(R) converges towards the average density of the liquid, g(R)→1.

[[Image:FBresme-FIG14.jpg]]

Liquid water has a structure that is much more complicated than that of a simple liquid. One of the reasons of this complex structure is the existence of hydrogen bonding.  

(a)	Using the SPC/E model investigate the structure of liquid water  at standard conditions, p=1 bar and T=298 K. Since water consists of  oxygen and hydrogen, there are three pair correlation functions that define the water structure, gOH, gOO and gHH. 

To compute the pair correlation functions you need to perform your simulation (NpT ensemble) adding the following commands to the CONTROL file:

rdf 5

print rdf

This tells DL_POLY to compute the pair correlation function every fifth time step and print it at the end of the run. When the run finishes DL_POLY creates a file called '''RDFDAT'''.

You can analyse the pair correlation functions using the DL_POLY GUI. (See “Analysis/Structure/RDF-Plot” option). Compare your pair correlation function with the experimental pair correlation functions for water<ref>A.K. Soper, F. Bruni and M.A. Ricci, J. Chem. Phys., 106, 247 (1997).</ref>. Comment on your results.

(b) Walrafen<ref>G.E. Walrafen, J. Chem. Phys., 40, 3249 (1964).</ref>  proposed a hydrogen bond structure for water consisting of a tetrahedral motif of water molecules (see Figure below).

[[Image:FBresme-FIG15.gif]]

Are your pair correlation functions compatible with this structure?

To answer this question you may want to investigate whether the peaks in the radial distribution functions are consistent with the oxygen-oxygen, oxygen-hydrogen distances in the tetrahedral structure suggested by Walrafen. Also try to compute the coordination numbers, nOO, nOH and nHH as a function of radial distance, R. The coordination number gives the number of atoms of a given species surrounding a specific atom, and it can be computed by calculating (numerically) the following integral:

[[Image:FBresme-FIG16.gif]]

Where α and β represent the species, O or H, and ρβ, is the density of species β. If Walrafen's model is correct the nOH coordination number at R~2.4 Ǻ should be ~2. Also the nOO coordination number at 3.4 Ǻ should be ~4, compatible with the four water molecules surrounding the central one (See figure above).

<div class="references-small"></div>
<references/>

==== Exercise 3: Water structure near the critical point ====

The temperature has an important effect on the hydrogen bonding structure of water. In the following we consider a thermodynamic state near the critical point of the water SPC/E model. The estimated critical point is Tc=651.7 K, ρc=0.326 g/cm3, and Pc=189 bars, which is good agreement with the experimental critical properties of water <ref> Y. Guissani, B. Guillot, J. Chem. Phys., 98, 8221 (1993)</ref>. Near the critical point the liquid undergoes very large fluctuations in density. This is reflected in the compressibility of the high temperature liquid which is different from that of water at normal conditions. 

(a) Perform one simulation of liquid water at P=160 bars and T=620 K. Monitor the equilibration through the variation of the volume of the simulation cell with time. What is the equilibrium density for this thermodynamic state?

(b) Compute the radial distribution function of water at these thermodynamic conditions. Compare your result with the radial distribution function at normal conditions.  Is there hydrogen bonding still present at P=160 bars and T=620 K? Justify your answer.

<div class="references-small"></div>
<references/>

==== Exercise 4: Computer simulations of aqueous solutions ====

Most processes of  interest in colloidal chemistry or biochemistry occur in aqueous solutions.  In the second year (Electrochemistry course) you studied the Debye-Huckel theory (DHT). This is the most basic approach to study ionic solutions from a theoretical perspective. This theory is powerful due to its simplicity and accuracy at low ionic concentrations. The DHT treats water as a structureless dielectric medium. This assumption works well for highly diluted solutions, but it becomes inaccurate at relatively high concentrations. Computer simulations provide an approach to investigate the structure and thermodynamics of ionic solutions including solvent effects. 

In this exercise you will  investigate the solvation structure of ions in water. The first part of the exercise involves the simulation of a NaCl solution - concentration 2.5 molal -  at ambient conditions. You will need a new configuration [http://www.ch.ic.ac.uk/bresme/lab3www/CONFIG CONFIG] and a new [http://www.ch.ic.ac.uk/bresme/lab3www/FIELD FIELD] file. The figure below shows an example of a configuration where Cl- anions and Na+ cations are represented as big spheres (yellow) and  small spheres (green) respectively.  

[[Image:FBresme-FIG17.jpg]]

(a) Perform an NpT simulation at 1 bar pressure and 298 K. Compute the pair correlation functions, and analyse the water structure around the ions. 

Try to address the following points in your investigation: 

•How many water molecules are solvating the  Na+ and Cl- ions? (to answer this question you will need to calculate the Na+ O and Cl- O coordination numbers. To do this integrate the pair correlation functions up to the first minimum). 

•How do the ions modify the structure of water? 

•How do Na+ and Cl- ions organise in the aqueous solution? Do you find any evidence for the formation of ion pairs?

(b) If you have time, perform a simulation  at P=160 bars and T=620 K. This thermodynamic state is near the  water critical point. Calculate the Na+ Cl- correlation functions and compare your result with the aqueous solution at 298 K. Explain the differences in the ionic structure of both thermodynamics states.