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 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 ); 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 polarized under temperature gradients... etc

The computer simulation of liquid water is far from trivial. Many theoretical models have been developed over the years). In this tutorial you we will be using the SPC/E model, 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, high pressure and near critical conditions.

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:

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 (Ewald precision 1e-6), to compute the coulombic interactions between the partial charges in the water SPC/E model. In order 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 atomic coordinates. This file contains the positions of the atoms we want to simulate. An example of a CONFIG file is given below:

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.

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).

I 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):

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.

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).

To generate a 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).

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 the files CFGH20.0, FLDBLK.0 and CNTROL.0, to CONTROL, CONFIG and FIELD respectively .

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.

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:

Check that the definitions in your FIELD file are consistent with the SPC/E model parameters published by Berendsen et al .

IMPORTANT: The input file, FIELD, is read with a fixed format. This means that the variables defining the potential, charge etc, have to be separated by a pre-defined number of spaces between the different definitions. If you are getting error messages when trying to execute DLPOLY, 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

To find the meaning of the labels you are referred to the DL_POLY [ http://www.cse.scitech.ac.uk/ccg/software/DL_POLY/MANUALS/USRMAN2.19.pdf manual]. The properties we will consider in this experiment 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).

Try plotting other properties, temperature, configurational, van der Waals and coulombic energies. Note that each time you select one of these properties the GUI generates 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. This is a good way to see whether the timestep you are using is correct.

To explore the the impact that the timestep has on the computer simulation 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 solid circles in the phase diagram given below

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 maintain the pressure and temperature constant. The two parameters after the ensemble definition, “0.1 ps”, are constants that specify the relaxation of the thermostat and barostat, respectively. “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 5000 steps, (2500 for equilibration).

(b)Calculate the configurational energy of water at normal conditions. Compare your result with the data given in references and .

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.

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:

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. ,.

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.

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. Comment on your results.

(b) Walrafen proposed a hydrogen bond structure for water consisting of a tetrahedral motif of water molecules (see Figure below).

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:

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).

Exercise 3: Water at high temperatures

The temperature has an important effect on the hydrogen bonding structure of water. In the following we consider a thermodynamics 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 . 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, and that diverges at the critical point.

(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) Computer 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.

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 explicitly the 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 CONFIG and a new 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.

(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 within the first coordination shell of Na+ and Cl-? (To answer this question you will need to calculate the Na+ O and Cl- O coordination numbers. You can integrate the pair correlation function up to the first minimum).

How does the salt modify the structure of water with respect to the structure of bulk 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.