182 P.Entel et al. Neutron diffraction scattering dsta Car-Parrinello MD 2 0 0 4567 8 r(A) Fig.9.1.The oxygen-oxygen (top),oxygen-hydrogen (middle)and hydrogen- hydrogen (bottom)radial distribution functions obtained from a Car-Parrinello molecular dynamics simulation of a system of 64 H2O molecules(solid lines)[9.27] compared to results of neutron diffraction scattering 9.26(dashed lines).The total simulation time was 11 ps and the average ionic temperature was 307 K.The figure was adapted from [9.27]. Of all the RDFs,goo(r)is perhaps the most informative.Calculation of the coordination number defined by Tmin Nc=4πp drr2goo(r), (9.9) (where rmin is the location of the first minimum in goo(r))gives the number of nearest neighbor water molecules.The recent X-ray data give Ne =4.7 9.25 while the CP simulations yield an average number of H2O molecules in the first coordination shell of about 4.0.A coordination number N<5 indicates that liquid water preserves much of its tetrahedral,ice-like structuring,while differences from the ice-like structuring can be found from the hydrogen-bonding patterns,see,for example,9.33.The ratio of the dis- tances of the first two peak positions in goo(r)is about 1.64 close to the ideal
182 P. Entel et al. 0 1 2 3 gOO(r) Neutron diffraction scattering data Car-Parrinello MD 0 1 2 3 gOH(r) 0 12345678 r (Å) 0 1 2 3 gHH(r) Fig. 9.1. The oxygen-oxygen (top), oxygen-hydrogen (middle) and hydrogenhydrogen (bottom) radial distribution functions obtained from a Car-Parrinello molecular dynamics simulation of a system of 64 H2O molecules (solid lines) [9.27] compared to results of neutron diffraction scattering [9.26] (dashed lines). The total simulation time was 11 ps and the average ionic temperature was 307 K. The figure was adapted from [9.27]. Of all the RDFs, gOO(r) is perhaps the most informative. Calculation of the coordination number defined by Nc = 4πρ r min 0 dr r2 gOO(r), (9.9) (where rmin is the location of the first minimum in gOO(r)) gives the number of nearest neighbor water molecules. The recent X-ray data give Nc = 4.7 [9.25] while the CP simulations yield an average number of H2O molecules in the first coordination shell of about 4.0. A coordination number Nc < 5 indicates that liquid water preserves much of its tetrahedral, ice-like structuring, while differences from the ice-like structuring can be found from the hydrogen-bonding patterns, see, for example, [9.33]. The ratio of the distances of the first two peak positions in gOO(r) is about 1.64 close to the ideal
9 Molecular Dynamics Simulations in Biology,Chemistry and Physics 183 value of 2(2/3)1/2=1.633 in an ice Ih-type lattice.The increasing width with increasing distance is an indication of increasing disorder in the liquid.The first two peaks in gon(r)and gHH(r)correspond to the average intramolecu- lar distances of O-H and H-H,respectively.For further discussions we refer to[9.27]. The RDFs of water and ice from 220 to 673 K and at pressures up to 400 MPa have recently been discussed on the basis of neutron scattering data 9.26.It is interesting to note that in the ice formation there is still substantial disorder in the hydrogen bonding pattern as can be checked from the width of the RDFs.MD simulations of the phase transition,i.e.,freezing of water to ice,are more difficult to achieve than melting of ice.There have only been a few successful MD runs of free (i.e.,not confined)water which show ice nucleation and subsequent percolation of the nucleus throughout the simulation box containing 512 water molecules [9.34].Due to the complex global potential-energy surface,a large number of possible network config- urations are possible.This causes large structural fluctuations showing up in the simulations hindering the system to find an easy pathway from the liquid to the frozen state (in spite of the fact that water molecules forming tiny ice-like clusters with four-coordinated hydrogen bonds have by 2 kJ/mol lower potential energy than that of other water molecules 9.34).Results of MD simulations of ice nucleation are shown in Fig.9.2. The constant-temperature MD simulations have been done for 512 mole- cules in the simulation box with a time step of 1 fs.The TIP4P model for water has been employed,which is a flat 4-center model with a potential energy consisting of Coulomb and Lennard-Jones terms,whereby the oxygen charge is shifted to a fourth site located closer but equidistant from the two hydrogen atoms 9.35. The essential four stages in the freezing process are discussed in detail in [9.34].The simulations show that the number of six-member rings N6R fuctuates during the simulations but only increases after an initial nucleus has formed and the ice-growth process has set in,see Fig.9.2b and following. The authors investigated also the so-called Q6 parameter [9.36,which may serve to characterize in how far the hydrogen bonds are coherently (icosa- hedrally)ordered.From the simulations it seems to follow that neither N6R nor Q6 are suitable order parameters to describe the entire freezing process. Obviously coherent icosahedral orientational correlations are an imperfect so- lution to characterize tetrahedral packing of water molecules as in ice [9.36. The reverse process,melting of ice,is easier to simulate.In the following we briefly present results of melting of water clusters simulated by our group. Numerous studies have been devoted to understand the dynamics of small clusters of water since the beginning of simulation studies in the 1970s and to elucidate the nature of the pseudo-first order melting transition. Our aim is to simulate the melting of water clusters (H2O)n of selected sizes (shown in Fig.9.3)and how their properties evolve with size from ab initio type of calculations using the DFTB method [9.37].Some of the bigger
9 Molecular Dynamics Simulations in Biology, Chemistry and Physics 183 value of 2(2/3)1/2 = 1.633 in an ice Ih-type lattice. The increasing width with increasing distance is an indication of increasing disorder in the liquid. The first two peaks in gOH(r) and gHH(r) correspond to the average intramolecular distances of O–H and H–H, respectively. For further discussions we refer to [9.27]. The RDFs of water and ice from 220 to 673 K and at pressures up to 400 MPa have recently been discussed on the basis of neutron scattering data [9.26]. It is interesting to note that in the ice formation there is still substantial disorder in the hydrogen bonding pattern as can be checked from the width of the RDFs. MD simulations of the phase transition, i.e., freezing of water to ice, are more difficult to achieve than melting of ice. There have only been a few successful MD runs of free (i.e., not confined) water which show ice nucleation and subsequent percolation of the nucleus throughout the simulation box containing 512 water molecules [9.34]. Due to the complex global potential-energy surface, a large number of possible network configurations are possible. This causes large structural fluctuations showing up in the simulations hindering the system to find an easy pathway from the liquid to the frozen state (in spite of the fact that water molecules forming tiny ice-like clusters with four-coordinated hydrogen bonds have by 2 kJ/mol lower potential energy than that of other water molecules [9.34]). Results of MD simulations of ice nucleation are shown in Fig. 9.2. The constant-temperature MD simulations have been done for 512 molecules in the simulation box with a time step of 1 fs. The TIP4P model for water has been employed, which is a flat 4-center model with a potential energy consisting of Coulomb and Lennard-Jones terms, whereby the oxygen charge is shifted to a fourth site located closer but equidistant from the two hydrogen atoms [9.35]. The essential four stages in the freezing process are discussed in detail in [9.34]. The simulations show that the number of six-member rings N6R fluctuates during the simulations but only increases after an initial nucleus has formed and the ice-growth process has set in, see Fig. 9.2b and following. The authors investigated also the so-called Q6 parameter [9.36], which may serve to characterize in how far the hydrogen bonds are coherently (icosahedrally) ordered. From the simulations it seems to follow that neither N6R nor Q6 are suitable order parameters to describe the entire freezing process. Obviously coherent icosahedral orientational correlations are an imperfect solution to characterize tetrahedral packing of water molecules as in ice [9.36]. The reverse process, melting of ice, is easier to simulate. In the following we briefly present results of melting of water clusters simulated by our group. Numerous studies have been devoted to understand the dynamics of small clusters of water since the beginning of simulation studies in the 1970s and to elucidate the nature of the pseudo-first order melting transition. Our aim is to simulate the melting of water clusters (H2O)n of selected sizes (shown in Fig. 9.3) and how their properties evolve with size from ab initio type of calculations using the DFTB method [9.37]. Some of the bigger
184 P.Entel et al. 208ns ns Fig.9.2.Change of the hydrogen bond network structure of water with the sim- ulation time.During the freezing (for times t 320 ns)the gradual formation of a perfect honeycomb structure becomes visible(lower panel)accompanied by a considerable decrease of the potential energy or loss of kinetic energy of the water molecules.The MD simulation of water freezing is performed for 512 molecules (with density 0.96 g/cm)in the simulation box with periodic boundary conditions and involves thermalization at a high temperature followed by quenching(at time t=0)to a low temperature of 230 K(supercooled water).After 256 ns a polyhedral structure consisting of long-lasting hydrogen bonds forms spontaneously acting as an initial nucleus,see the circled region in (b).The lines indicate hydrogen bonds among water molecules.The brightest blue lines are those with hydrogen-bond lifetimes of more than 2 ns.Reprinted by permission from Nature copyright 2002 Macmillan Publishers Ltd.(http://www.nature.com/nature)[9.34]. clusters were modeled from the existing smaller units of n 3,4,5 and 6.The preference is given to those with lowest energy.A fair search for the minima of these structures was carrried out before the structures were used as the starting geometry in the MD simulation runs.The results of our simulation are compared with some of the results available for the existing clusters from simulation with a model classical pairwise additive potential (CPAP)9.38. We have,however,to point out that results to be found in the literature for the apparent global minimum structures of water clusters differ for some of
184 P. Entel et al. Fig. 9.2. Change of the hydrogen bond network structure of water with the simulation time. During the freezing (for times t > 320 ns) the gradual formation of a perfect honeycomb structure becomes visible (lower panel) accompanied by a considerable decrease of the potential energy or loss of kinetic energy of the water molecules. The MD simulation of water freezing is performed for 512 molecules (with density 0.96 g/cm3) in the simulation box with periodic boundary conditions and involves thermalization at a high temperature followed by quenching (at time t = 0) to a low temperature of 230 K (supercooled water). After 256 ns a polyhedral structure consisting of long-lasting hydrogen bonds forms spontaneously acting as an initial nucleus, see the circled region in (b). The lines indicate hydrogen bonds among water molecules. The brightest blue lines are those with hydrogen-bond lifetimes of more than 2 ns. Reprinted by permission from Nature copyright 2002 Macmillan Publishers Ltd. (http://www.nature.com/nature) [9.34]. clusters were modeled from the existing smaller units of n = 3, 4, 5 and 6. The preference is given to those with lowest energy. A fair search for the minima of these structures was carrried out before the structures were used as the starting geometry in the MD simulation runs. The results of our simulation are compared with some of the results available for the existing clusters from simulation with a model classical pairwise additive potential (CPAP) [9.38]. We have, however, to point out that results to be found in the literature for the apparent global minimum structures of water clusters differ for some of
9 Molecular Dynamics Simulations in Biology,Chemistry and Physics 185 n=20 n=27 n=30 n=36 Fig.9.3.Initial configuration for some selected (H2O)n clusters.The number of molecules in each of the clusters is indicated by n. the clusters (compare,for example,the global minmum structures of water clusters calculated using the TIP3P and TIP4P potentials [9.39]).Our final relaxed structures for some of the smaller clusters such as n =3,4,5,6,8 and 20 are in agreement with some of the relaxed geometries of Wales and Lee[9.38,9.40,9.41]. The melting temperature can then be obtained from the inflexion point of the calorific curve (energy versus temperature).However,it is sometimes difficult to see the abrupt change in the energy.Hence,instead Lindemann's criteria of melting (along with the former condition)is used,which is obtained from the root-mean-square bond-length fluctuations ooo of oxygen in each of the clusters as given by 2 (r》-(r》2 800= N(N-1) (9.10) (r》 where the brackets denote time averages,and rij is the distance between the oxygen atoms i and j.The summation is over all N molecules.According to this criterion,melting is caused by a vibrational instability when ooo reaches a critical value.We paid attention of how to define the melting temperature Tm for a particular cluster.A cluster of water may have different isomeriza-
9 Molecular Dynamics Simulations in Biology, Chemistry and Physics 185 n = 8 n = 10 n = 15 n = 27 n = 30 n = 18 n = 3 n = 4 n = 5 n = 12 n = 6 n = 20 n = 24 n = 36 Fig. 9.3. Initial configuration for some selected (H2O)n clusters. The number of molecules in each of the clusters is indicated by n. the clusters (compare, for example, the global minmum structures of water clusters calculated using the TIP3P and TIP4P potentials [9.39]). Our final relaxed structures for some of the smaller clusters such as n = 3, 4, 5, 6, 8 and 20 are in agreement with some of the relaxed geometries of Wales and Lee [9.38, 9.40, 9.41]. The melting temperature can then be obtained from the inflexion point of the calorific curve (energy versus temperature). However, it is sometimes difficult to see the abrupt change in the energy. Hence, instead Lindemann’s criteria of melting (along with the former condition) is used, which is obtained from the root-mean-square bond-length fluctuations δOO of oxygen in each of the clusters as given by δOO = 2 N(N − 1) i<j 7 r2 ij −rij 2 rij , (9.10) where the brackets denote time averages, and rij is the distance between the oxygen atoms i and j. The summation is over all N molecules. According to this criterion, melting is caused by a vibrational instability when δOO reaches a critical value. We paid attention of how to define the melting temperature Tm for a particular cluster. A cluster of water may have different isomeriza-
186 P.Entel et al. 0 (H,O)cluster (HO)cluste -0. 0.2 0.2 0. 0.3 0.8 4 035 0.25 0. 0 02 0.15 0.1 0. 0.05 0.05 00 50.100150200250300350 50 100150200250300350 T(K) T(K Fig.9.4.Calculated energy per(H2O)n cluster together with the root mean square fluctuations in the O-O bond lengths versus temperature (n=3,4).The vertical dotted lines mark the approximate melting tmperatures. tion in which there is interconversion of a water cluster from one form of isomer to another when the hydrogen-bond breaking and reforming is sub- stantial,leading to new structures before melting.This behavior shows up as fluctuations in oo(T).Figure 9.4 shows the behavior of the energy and of oo versus temperature for the case of two small clusters. The size dependence of the melting temperature obtained with the DFTB method compared to calculations using model potentials (9.38,9.42,a pair- 350 SPC 300 DFTB CPAP 250 200 150 100 5 15 0253035 40 N Fig.9.5.The melting temperature against the number of water molecules.Both results of our DFTB calculations and simulations using the CPAP 9.38 classical potential are plotted.Also shown is the Ti value for the pentamer cluster ob- tained from a point charge model [9.42].Remarkable are the rather high melting temperatures for the small clusters
186 P. Entel et al. -0.5 -0.4 -0.3 -0.2 -0.1 0 Energy (eV) -1 -0.8 -0.6 -0.4 -0.2 0 0 50 100 150 200 250 300 350 T (K) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 δOO (Å) 0 50 100 150 200 250 300 350 T (K) 0 0.05 0.1 0.15 0.2 0.25 (H2 O)3 cluster (H2 O)4 cluster Fig. 9.4. Calculated energy per (H2O)n cluster together with the root mean square fluctuations in the O–O bond lengths versus temperature (n = 3, 4). The vertical dotted lines mark the approximate melting tmperatures. tion in which there is interconversion of a water cluster from one form of isomer to another when the hydrogen-bond breaking and reforming is substantial, leading to new structures before melting. This behavior shows up as fluctuations in δOO(T). Figure 9.4 shows the behavior of the energy and of δOO versus temperature for the case of two small clusters. The size dependence of the melting temperature obtained with the DFTB method compared to calculations using model potentials ( [9.38,9.42], a pair- 0 5 10 15 20 25 30 35 40 N 100 150 200 250 300 350 Tm (K) SPC DFTB CPAP Fig. 9.5. The melting temperature against the number of water molecules. Both results of our DFTB calculations and simulations using the CPAP [9.38] classical potential are plotted. Also shown is the Tm value for the pentamer cluster obtained from a point charge model [9.42]. Remarkable are the rather high melting temperatures for the small clusters