Module 3 - Molecular Dynamics Simulation of Liquids

The Molecular Dynamics of water and an aqueous solution was studied using a Molecular Dynamics simulation code, DL_POLY. The model for water that we are using is the extended single point charge or SPC/E model. The SPC/E model consists of three point masses with a distance of 0.1nm (OH distance) and an angle of 109.5° (tetrahedral). The point charges are -0.8476e and 0.4238e and it has a dipole moment of 2.36D.^[1] Water is a very interesting material to study seeing as it is the most abundant compound available on the planet and it has a multitude of anomalous physical properties. Primarily we will look at water at ambient conditions, as it is most commonly found.

Example - Water at Ambient Temperature

The CONTROL file was altered so that 500 steps were carried out. All 500 of these steps were equilibration steps. The timestep was 0.002 ps so the whole simulation would take 1 ps. The calculation would be carried out at 298K in a microcanonical ensemble (NVE). The FIELD file was altered to the SPC/E parameters. DL_POLY was run but an error came up: "SHAKE algorithm failed to converge." The error was due to the incorrect formatting of the FIELD file. This problem was overcome and when DL_POLY was run again the OUTPUT, STATIS and REVCON files appeared. Using the GUI (Java Graphical User Interface) the STATIS file was used to plot system properties as a function of time. The REVCON file was opened in GUI also to give the final configuration of the molecules.

Total Energy fluctuating with Time	Temperature Fluctuation with Time	Configuration Energy Fluctuation with Time

Van der Waals fluctuation with Time	Coulombic Energy Fluctuating with Time	Final Configuration of Water Molecules

The Total Energy obviously decreases with time because it's components (Coulombic, Configuration and Van der Waals energy contributions) decrease with time as an energetically favourable configuration is reached. Van der Waals contributions are due to oxygen pair repulsions only , therefore energy is positive. Conformation and coulombic interactions are negative. Coulombic interactions are due to O-H and H-H interactions. O-H interactions are strongly negative due to hydrogen bonding, which probably leads to the overall negative total energy of the system. Hydrogen bonding will be investigated further in Exercise 2.

Temperature fluctuates between 298K and 294K over time. This is due to heat energy dissipating on collision or energy used up in changing conformation etc.

The OUTPUT file gave the following information which was compared to the parameters given in reference 1. The information was taken from the last rolling average (Data point 450) in the OUTPUT file:

Physical Property	From our calculation
eng_tot	-8.9103 kJ mol-1
temp_tot	296.14K
eng_cfg	-10.828E4 kJ mol-1
volume	1.0648E4 Å³
pressure	1.5741 (bar?)

The Final Conformation shows us that water molecules are not equally spaced at ambient temperature. They seem to be 'bunched' together in the centre of the cube. This is due to energetically favourable hydrogen bonding. There seem to be some 'broken' molecules on the outer edge of the 'hydrogen bonded sphere' which is due to the Periodic Boundary Conditions set by the DL_POLY code. These conditions approximate the system as infinte.

Time Step Study

The time step was altered to 0.01ps and 0.0020ps to see the effect of the time step on the results. Initially the number of steps was kept at 1000 from the previous exercise but this caused an error with the SHAKE algorithm on using the 0.01ps time step because the time step was too small which did not allow convergence. The number of steps was altered so the entire molecular dynamics calculation takes 1ps. The graphs below show the Total Energy Fluctuation with Time at different time steps.

The time step of 0.01ps is far too long. Although the minimum energy is obtained successfully, there is no fine detail on the curve which could be useful as more information can probably be obtained from the detail. The time step of 0.002 shows fine detail on the curve. Using more steps does not give us more information because all the information we need is within the 1ps timeframe.

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

Previously we have used the microcanonical ensemble, but we would now like to be able to set the temperature and pressure of our system. In this calculation at ambient conditions we will use the isobaric-isothermal ensemble (nPT). A SPC/E calculation was performed at 298K and 1bar. 4000 steps were used, half of which were equilibration steps. A trajectory of 1ps was used so the time step used was 0.00025ps.

The volume fluctuation of the system with time shows can tell us about the equilibration of the system. From the graph shown we can see that volume decreases overall with time. During equilibration volume decreases steadily as the water molecules find their equilibrium positions and hydrogen bonds form. After equilibration the volume fluctuates at regular 0.1ps intervals. This could be synonymous with the temperature fluctuation that occurred in the example calculation.

The E_config of the system was calculated to be -11.422kJ/mol. This is in agreement with literature (-11.5kJ/mol^[2]) There are 256 molecules of water in the system so this works out as -0.0446 kJ/molecule. Using the final volume given on the output (8.83896x10^-18cm³) and the weight of 256 water molecules (7.652x10^-18g) the density was worked out as 0.8632 g/cm³. The literature value for an SPC/E calculation is 1g/cm³ (reference 2). Although the magnitude is correct, there is an error of 14% between my calculation and the one in literaure. This could be due to the trajectory of 1ps not being long enough. In the first reference, the equilibration runs for 15ps before the calculation starts. Seeing as the volume is decreasing with time, the longer the trajectory the closer the calculated density of water would be.

Exercise 2: Liquid water structure and hydrogen bonding

In this exercise the structure of water will be established using the pair correlation function g(R). This is defined as the probability of finding a molecule within a distance, R from another molecule. The function should tend to one as distance increases as density becomes close to that of bulk water. The CONTROL file was altered so that the pair correlation functions could be calculated every 5 steps:

rdf 5
print rdf

The following g(OO), g(OH), g(HH) pair correlation functions were obtained from DL_POLY using the SPC/E model at ambient temperature.

g(OO)	g(OH	g(HH)

Max. (2.775Å, 3.49913)	Max. (1.775Å, 1.824128) (3.175Å, 1.635579)	Max. (2.375Å, 1.483233) (3.925Å, 1.199994)

All the pair correlation functions calculated seem to tend to 1 (the bulk density of water) as distance from the central molecule increases. This is not in agreement with previous density calculations (ρ_calc=0.86g/cm³) but is in agreement with literature (see previous exercise). The first peak on the g(OH) pair correlation function is in agreement with literature (lit. 1.8Å ^[3]), and the other pair correlation functions seem to fit the results in literature also. The first peak on the g(OH) plot relates to the first coordination shell of water. The coordination values calculated (see below) show us there are 3 molecules in the first coordination shell. This does not quite relate to a tetrahedral structure but the tetrahedral structure might be an average, not a perfect structure as in ice. At the time the g(R) is measured the structure may be away from the tetrahedral structure. This applies to the other n(R) values also. Seeing as the regions where the function lies above 1 are the areas of above average density (compared to water as an ideal gas), we can assume these areas show us the 'shells' of hydrogen bonded water molecules. The peak with the shortest distance is obviously gOH because hydrogen bond forms between oxygen and hydrogen. Areas around a central molecule where hydrogen bonding occur are likely to be closest to the central molecule and have the highest density of molecules. 1.775Å seems to be the length of a hydrogen bond acquired from this calculation. The first peak in the gOO plot should be an OH sigma bond [in water] added to the length of a hydrogen bond, presuming linearity. The length of an OH bond in the CONTROL file is 1Å the angle O-H-O (the shortest distance between oxygens) is 180°. The density decreases as we get further from the centre and the density tends to 1. The highest g(R) occurred at g(OO). This could be due to the corresponding distance being equivalent to the centre of the 1st hydration shell and because oxygen is hydrogen bonded to several other water molecules at once. The calculated pair correlation functions seem to be compatible with the tetrahedral structure for hydrogen bonding suggested by Walrafen^[4]. To verify the compatability the coordination of OH and OO was determined using the equation given.

	g(OH)	g(OO)	g(HH)
n_αβ(R)	3.04; ~3 coordination (lit. 2 coordination) outer coordination shell: n(R)=7.95; 8 coordination	2.2; ~2 coordination (lit. 4 coordination)	7.64; ~8 coordination (from Walrafren diagram, coordination ~6)

There are errors associated with this calculation due to the numerical nature of the integration method used. However, the numbers are at the correct magnitude for coordination numbers which is satisfying!

Exercise 3: Water structure near the critical point

We have studied some physical properties of water at ambient temperature. Now we will look at water when it is near it's critical point. The estimated conditions for water at the critical point is T_c=651.7 K, ρ_c=0.326 g/cm3, and P_c=189 bars but we will run calculations at near the critical point: P=160bars, T=620K.

The volume fluctuation of the system is completely contrary to the volume v t plot obtained from the previous exercise. Volume is increasing with time and therefore density is decreasing with time. The trajectory is constant throughout both calculations but we can observe there are fewer small fluctuations of volume after equilibration. This could be due to the higher pressure making it physically harder for volume to increase and decrease as fast as it could at ambient conditions. Using the last data point on the vol/t plot, the density is equal to 0.66 g/cm³. This is not in agreement with literature (ρ_lit=0.326g/cm³^[5] but this 50% error could be due to the trajectory of 1ps not being sufficient to let the system equilibrate sufficiently. According to Guissani and coworkers, a higher density can be attributed to a higher number of hydrogen bonds i.e. 3 formed at the critical point instead of 2.

g(OO)	g(OH)	g(HH)

Max. (2.925Å, 2.026753)	Max. (1.875Å, 0.817051) (3.325Å, 1.326993)	Max. (2.575Å, 1.08429) (3.925Å, 1.132529)
n(OO)=4.18; ~4 coordination	n(OH)=0.69; ~1 coordination n(OH)=4.25; ~4 coordination	n(HH)=3.47; ~3 coordination n(OH)=4.21; ~4 coordination

From the coordination calculations we can see that the lower density compared to ambient temperatures can be attributed to fewer hydrogen bonds and longer hydrogen bonds. From the n(OH) function ,the first coordination shell only contains one molecule of water. The outer shell contains 4 water molecules, then from there the g(R) function tends to bulk water density. It is clear from our calculations that as water approaches the critical point, water starts to move away from the tetrahedral structure that was observed at ambient conditions. Hydrogen bonding may still be present in the first coordination sphere but it is not as pronounced as at ambient conditions as there is only one water molecule to hydrogen bond to. The g(OH) plot shows the hydrogen bonds at these

conditions are lengthening, and I presume they will break completely at the critical point. The REVCON file shown opposite shows the water molecules are more diffuse than the REVCON file at ambient temperature due to transition to the gas phase.

Exercise 4: Computer simulations of aqueous solutions

In this exercise the solvation of ions in water will be investigated. The pair correlation functions gained will be compared to the pure water calculations already run to determine the structure of an aqueous solution. The CONFIG and FIELD files were altered to include a 2.5 molal solution of NaCl. As can be seen from the REVCON file opposite, the ions look randomly dispersed among the water molecules but there is a slight pairing effect of ions where the big yellow spheres (Na ions) are in close proximity to a small green sphere (Cl ions).

g(ClO)	g(NaO)	g(OH)	g(OO)

Max. (2.925Å, 2.026753)	Max. (2.42Å, 6.676297) (4.575Å, 1.424192)	Max. (1.775Å, 1.392551) (3.225Å, 1.554763)	Max. (2.775Å, 2.796648)
n(ClO)=1.79; ~2 coordination	n(NaO)=2.44; ~2 coordination n(NaO)=33.9; ~34 coordination	n(OH)=2.38; ~2 coordination

The g(R) plots for water are similar in shape to those calculated in exercise 2, apart from the slightlydecreased density of hydrogen bonding. The length of hydrogen bonds are constant in both calculations, also. The number of water molecules in the outer shell is 2, which is less than the calculation in exercise 2 but is in agreement with literature. Seeing as these calculations give us a value that is agreement with literature, we could conclude that ions in solution aid the tetrahedral structure to be maintained but the water molecules are not as dense as in pure water. Now we know that the structure of water is not much altered on addition of a salt we can look at the coordination of the ions in solution. The 'length' of the NaO and ClO bonds in literature are 2.758Å and 3.783Å ^[6] which do correspond with our results in terms of magnitude and the fact that the Cl-O interaction is longer than the Na-O interaction. This is obviously to do with the charge on the ions and that the H atom on the water molecule will coordinate to the Cl, not the O atom. There is not q difference of an O-H bond (1Å) in the calculation between the two peaks like there is in literature which argues against the tetrahedral structure of water in an aqueous solution. From integrating under the peaks it has been calculated that on average there are more water molecules solvating the Na⁺ than the Cl^-. This is in agreement with literature ^[7].

Structure of Ions in Solution

RDF NaCl)	RDF NaNa	RDF ClCl

This RDF plot is not the same as the other two plots which suggests pairs are formed of NaCl in solution. The first peak is at approx 3Å which could be the distance between the Van der Waals radius of one ion plus the radius of a hydration shell of the other ion pair. The other peaks in the plot are similar to the NaNa and ClCl plots in that there are peaks which relate to the random dispersion of ions.	After 4Å the random peaks plotted show that the ions are dispersed randomly in solution.	Similarly to Na-Na, the random, equally sized peaks show that the ions are dispersed randomly in soltuion. If the RDF was plotted for a larger distance it is probable that the peaks would continue in this way.

Aqueous Solution at Critical Point

The properties of the same solution at critical point was investigated briefly. Similarly to the previous exercise, it was noted that the hydration spheres were more diffuse and less dense. The ion pairs are also more diffuse and if their number was calculated, I presume there would be fewer ion pairs.

References

↑ H.J.C. Berendsen, J.R. Grigera and T.P. Straatsma, J. Phys. Chem., 91, 6269 (1987) DOI:10.1021/j100308a038
↑ E. Sanz et. al., Phys. Rev. Lett., 92, 255701 (2004) DOI:10.1103/PhysRevLett.92.255701
↑ A.K. Soper, F. Bruni and M.A. Ricci, J. Chem. Phys., 106, 247 (1997).DOI:10.1063/1.473030
↑ G.E. Walrafen, J. Chem. Phys., 40, 3249 (1964).
↑ Y. Guissani, B. Guillot, J. Chem. Phys., 98, 8221 (1993)DOI:10.1063/1.464527
↑ E. Sanz and C. Vega, J. Chem. Phys. 126, 014507 (2007)DOI:10.1063/1.2397683
↑ P.Wilde, Y2 Electrochemistry Lecture Course, 2007

[1] H.J.C. Berendsen, J.R. Grigera and T.P. Straatsma, J. Phys. Chem., 91, 6269 (1987) DOI:10.1021/j100308a038

[2] E. Sanz et. al., Phys. Rev. Lett., 92, 255701 (2004) DOI:10.1103/PhysRevLett.92.255701

[3] A.K. Soper, F. Bruni and M.A. Ricci, J. Chem. Phys., 106, 247 (1997).DOI:10.1063/1.473030

[4] G.E. Walrafen, J. Chem. Phys., 40, 3249 (1964).

[5] Y. Guissani, B. Guillot, J. Chem. Phys., 98, 8221 (1993)DOI:10.1063/1.464527

[6] E. Sanz and C. Vega, J. Chem. Phys. 126, 014507 (2007)DOI:10.1063/1.2397683

[7] P.Wilde, Y2 Electrochemistry Lecture Course, 2007

[1]

[2]

[3]

[4]

[5]

[6]

[7]