9 Molecular Dynamics Simulations in Biology, Chemistry and Physics P.Entell,W.A.Adeagbol,M.Sugiharal,G.Rollmann',A.T.Zayak, M.Kreth',and K.Kadau2 1 Institute of Physics,University of Duisburg-Essen,47048 Duisburg,Germany 2 Los Alamos National Laboratory,T-11,MS B262,Los Alamos,NM 87545,USA Abstract.We review recent progress in understanding fundamental processes in biology,chemistry and physics on the basis of ab initio and molecular dynamics simulations.The first step of the visual process involving the excitation of bovine rhodopsin after absorption of light is taken as an example from biochemistry to demonstrate what is nowadays possible to simulate numerically.The act of freezing of water has recently been simulated,for the first time successfully,by scientists from chemistry.Martensitic transformation in bulk and nanophase materials,a typical and hitherto not completely solved problem from solid state physics,is used to illustrate the achievements of multimillion atoms simulations. 9.1 Molecular Dynamics as a Multidisciplinary Numerical Tool Molecular dynamics (MD)has proved to be an optimum numerical recipe applicable to problems with many degrees of freedom from quite different fields of science.The knowledge of the energy or potential landscape of inter- acting particles,like electrons and atoms,enables one to calculate the forces acting on the particles and to study the evolution of the system with time. As long as classical mechanics is appropriate to describe the dynamics of the individual constituents (i.e.atoms or molecules),the Newtonian equations of motion can be related to the statistical mechanics of the(classical)particles by using the equipartition theorem,i.e.by combining the equations mifi Fi, (9.1) 0=NkzT-(∑m》 (9.2) for i=1...N particles.Although the equipartition theorem holds only for classical particles,to the authors'knowledge the combination of classical and statistical physics has also been used to simulate small molecules at low tem- peratures without critically discussing so far the limitations of 9.1 and 9.2 when applying to very small quantum mechanical systems.The forces in 9.1 are then simply related to the gradients of the potential energy surface(PES) of either the classical or the quantum mechanical N-particle system.For ped- agogical reasons let us recall the classical constant-temperature case.Omit- ting the statistical average in 9.2 means that the temperature is a measure P.Entel,W.A.Adeagbo,M.Sugihara,G.Rollmann,A.T.Zayak,M.Kreth,and K.Kadau, Molecular Dynamics Simulations,Lect.Notes Phys.642,177-206(2004) http://www.springerlink.com/ C Springer-Verlag Berlin Heidelberg 2004
9 Molecular Dynamics Simulations in Biology, Chemistry and Physics P. Entel1, W.A. Adeagbo1, M. Sugihara1, G. Rollmann1, A.T. Zayak1, M. Kreth1, and K. Kadau2 1 Institute of Physics, University of Duisburg–Essen, 47048 Duisburg, Germany 2 Los Alamos National Laboratory, T–11, MS B262, Los Alamos, NM 87545, USA Abstract. We review recent progress in understanding fundamental processes in biology, chemistry and physics on the basis of ab initio and molecular dynamics simulations. The first step of the visual process involving the excitation of bovine rhodopsin after absorption of light is taken as an example from biochemistry to demonstrate what is nowadays possible to simulate numerically. The act of freezing of water has recently been simulated, for the first time successfully, by scientists from chemistry. Martensitic transformation in bulk and nanophase materials, a typical and hitherto not completely solved problem from solid state physics, is used to illustrate the achievements of multimillion atoms simulations. 9.1 Molecular Dynamics as a Multidisciplinary Numerical Tool Molecular dynamics (MD) has proved to be an optimum numerical recipe applicable to problems with many degrees of freedom from quite different fields of science. The knowledge of the energy or potential landscape of interacting particles, like electrons and atoms, enables one to calculate the forces acting on the particles and to study the evolution of the system with time. As long as classical mechanics is appropriate to describe the dynamics of the individual constituents (i.e. atoms or molecules), the Newtonian equations of motion can be related to the statistical mechanics of the (classical) particles by using the equipartition theorem, i.e. by combining the equations mi¨ri = Fi, (9.1) 0 = 3 2N kBT − 2 i 1 2mir˙2 i 3 (9.2) for i = 1 ...N particles. Although the equipartition theorem holds only for classical particles, to the authors’ knowledge the combination of classical and statistical physics has also been used to simulate small molecules at low temperatures without critically discussing so far the limitations of 9.1 and 9.2 when applying to very small quantum mechanical systems. The forces in 9.1 are then simply related to the gradients of the potential energy surface (PES) of either the classical or the quantum mechanical N-particle system. For pedagogical reasons let us recall the classical constant-temperature case. Omitting the statistical average in 9.2 means that the temperature is a measure P. Entel, W.A. Adeagbo, M. Sugihara, G. Rollmann, A.T. Zayak, M. Kreth, and K. Kadau, Molecular Dynamics Simulations, Lect. Notes Phys. 642, 177–206 (2004) http://www.springerlink.com/ c Springer-Verlag Berlin Heidelberg 2004
178 P.Entel et al. of the instantaneous kinetic energy in the system (isokinetic simulation)and the equations of motions to be solved are 9.1 mifi=F miti, (9.3) which are derived from 9.1 and 9.2 by applying the Gaussian principle of least mean square action,hence the name "Gaussian thermostat".It is to note that the friction-like term in 9.3 is deterministic and gives rise to time- reversal invariant trajectories. Using a different approach Nose reproduced the canonical phase-space distribution,so that the kinetic energy can fluctuate with a distribution pro- portional to (9.4) 2m The canonical distribution is generated with deterministic and time-reversal invariant trajectories of the Hamiltonian W-∑器+i+1e7ns+ +Vq (9.5) 20 with the time-scale variable s,its conjugate momentum ps and the"heat bath mass"Q;N is the number of degrees of freedom.Nose proved that the mi- crocanonical distribution in the augmented set of variables is equivalent to a canonical distribution of the variables q,p,where p=p/s are the scaled mo- menta 9.2].Although a single harmonic oscillator is not sufficiently chaotic to reproduce the canonical distribution,it is interesting to see the result if the method is applied to it.In this case the result depends on Q:For larger values of Q the trajectories fill in a region between two limiting curves [9.3. For an N-particle system the results (expectation values,i.e.time averages,of observables)are independent of Q and hence of H,the corresponding method of controlling the temperature during the simulations is called "Nose-Hoover thermostat".The extension of the method with respect to isobaric(constant pressure)simulations by making use of the virial theorem [9.2]allows to deal with many experimental set-ups,where the experiments can be simulated by using the so-called (NPT)ensemble.Finally for the description of struc- tural changes (for example,of condensed matter materials)fluctuations of the size and geometry of the simulation box can be taken into account by us- ing the Parinello-Rahman scheme [9.4.For further details like,for example, algorithms to solve the equations of motion or a discussion of the statistical nature of the trajectories,we refer to the literature [9.5]. We close the introductory remarks by mentioning some new trends in connection with molecular dynamics.One advantage of molecular dynamics
178 P. Entel et al. of the instantaneous kinetic energy in the system (isokinetic simulation) and the equations of motions to be solved are [9.1] mi¨ri = Fi − j Fjr˙j j mjr˙2 j mir˙ i, (9.3) which are derived from 9.1 and 9.2 by applying the Gaussian principle of least mean square action, hence the name “Gaussian thermostat”. It is to note that the friction-like term in 9.3 is deterministic and gives rise to timereversal invariant trajectories. Using a different approach Nos´e reproduced the canonical phase-space distribution, so that the kinetic energy can fluctuate with a distribution proportional to exp − 1 β i p2 i 2mi . (9.4) The canonical distribution is generated with deterministic and time-reversal invariant trajectories of the Hamiltonian H = i p2 i 2mis2 + (N4 + 1)kBT ln s + p2 s 2Q + V (q) (9.5) with the time-scale variable s, its conjugate momentum ps and the “heat bath mass” Q; N4 is the number of degrees of freedom. Nos´e proved that the microcanonical distribution in the augmented set of variables is equivalent to a canonical distribution of the variables q, p , where p = p/s are the scaled momenta [9.2]. Although a single harmonic oscillator is not sufficiently chaotic to reproduce the canonical distribution, it is interesting to see the result if the method is applied to it. In this case the result depends on Q: For larger values of Q the trajectories fill in a region between two limiting curves [9.3]. For an N-particle system the results (expectation values, i.e. time averages, of observables) are independent of Q and hence of H, the corresponding method of controlling the temperature during the simulations is called “Nos´e-Hoover thermostat”. The extension of the method with respect to isobaric (constant pressure) simulations by making use of the virial theorem [9.2] allows to deal with many experimental set-ups, where the experiments can be simulated by using the so-called (NPT) ensemble. Finally for the description of structural changes (for example, of condensed matter materials) fluctuations of the size and geometry of the simulation box can be taken into account by using the Parinello-Rahman scheme [9.4]. For further details like, for example, algorithms to solve the equations of motion or a discussion of the statistical nature of the trajectories, we refer to the literature [9.5]. We close the introductory remarks by mentioning some new trends in connection with molecular dynamics. One advantage of molecular dynamics
9 Molecular Dynamics Simulations in Biology,Chemistry and Physics 179 is that the simulations involve the physical time in contrast to other tech- niques like Monte Carlo.This,however,limits the time that can be spanned to the pico-or nanosecond range if an integration step of the order of one femtosecond is used.This means that many interesting events like conforma- tional changes of large biomolecules,protein folding,solid state reactions and transformations cannot be followed directly.One possibility to circumvent this difficulty is to rely on hybrid methods like a combination of molecular dynamics and kinetic Monte Carlo,for example,in the simulation of growth processes in solid state physics 9.6.Another way is to freeze in the rapid oscillations,for example,of light atoms like hydrogen etc.and to retain only the dynamics involved with the heavier constituents.Finally there are new developments under way in many places to find optimum tools for the simu- lations of rare events which likewise can be done by a subtle modification of the potential function [9.7].For a separation of different time scales in view of biological systems see also 9.8,9.9. Apart from the time problem,the extension to larger length scales is another important challenge which has recently been addressed by coupling the region of atomistic simulations to some effective elastic medium 9.10. Of particular interest is the problem to extend the size of biological systems playing a role in signal transduction.The present state-of-the-art with respect to such simulations is to describe important reactive parts of the molecules quantum mechanically while the rest is treated by molecular mechanics with empirical potentials.This so-called QM/MM technique,which is based on partitioning the Hamiltonian into H=HQM HMM+HQM-MM, (9.6) allows,for example,to deal with the chromophore of rhodopsin on a quantum mechanical basis,whereas the whole sequence of amino acids of the protein as well as part of the lipid membrane into which the latter is embedded,can be described by model potentials (here HoM is the Hamiltonian to be treated by ab initio methods,HMM is the purely classical Hamiltonian to be treated by force field methods and HoM-MM is the Hamiltonian describing the inter- action between the quantum mechanical and classical region).Conceptually this procedure is straightforward,difficulties may arise when trying a rigor- ous treatment of the coupling of the two regions [9.11].This mixed QM/MM approach can be extended to describe the dynamics of the first excited sin- glet state.This was done recently on the basis of using pseudopotentials in the quantum mechanical region and the simulation package AMBER in the MM region [9.12].We have done alike simulations,in particular for rhodopsin, which will be discussed below.Instead of using pseudopotentials,our method relies on a density-functional based tight-binding approach 9.13,which al- lows to treat a much larger region quantum mechanically with less compu- tational efforts.In the latter approach the MM region is described by using CHARMM [9.14].Compared to ab initio quantum chemical approaches,the
9 Molecular Dynamics Simulations in Biology, Chemistry and Physics 179 is that the simulations involve the physical time in contrast to other techniques like Monte Carlo. This, however, limits the time that can be spanned to the pico- or nanosecond range if an integration step of the order of one femtosecond is used. This means that many interesting events like conformational changes of large biomolecules, protein folding, solid state reactions and transformations cannot be followed directly. One possibility to circumvent this difficulty is to rely on hybrid methods like a combination of molecular dynamics and kinetic Monte Carlo, for example, in the simulation of growth processes in solid state physics [9.6]. Another way is to freeze in the rapid oscillations, for example, of light atoms like hydrogen etc. and to retain only the dynamics involved with the heavier constituents. Finally there are new developments under way in many places to find optimum tools for the simulations of rare events which likewise can be done by a subtle modification of the potential function [9.7]. For a separation of different time scales in view of biological systems see also [9.8, 9.9]. Apart from the time problem, the extension to larger length scales is another important challenge which has recently been addressed by coupling the region of atomistic simulations to some effective elastic medium [9.10]. Of particular interest is the problem to extend the size of biological systems playing a role in signal transduction. The present state-of-the-art with respect to such simulations is to describe important reactive parts of the molecules quantum mechanically while the rest is treated by molecular mechanics with empirical potentials. This so-called QM/MM technique, which is based on partitioning the Hamiltonian into H = HQM + HMM + HQM−MM, (9.6) allows, for example, to deal with the chromophore of rhodopsin on a quantum mechanical basis, whereas the whole sequence of amino acids of the protein as well as part of the lipid membrane into which the latter is embedded, can be described by model potentials (here HQM is the Hamiltonian to be treated by ab initio methods, HMM is the purely classical Hamiltonian to be treated by force field methods and HQM−MM is the Hamiltonian describing the interaction between the quantum mechanical and classical region). Conceptually this procedure is straightforward, difficulties may arise when trying a rigorous treatment of the coupling of the two regions [9.11]. This mixed QM/MM approach can be extended to describe the dynamics of the first excited singlet state. This was done recently on the basis of using pseudopotentials in the quantum mechanical region and the simulation package AMBER in the MM region [9.12]. We have done alike simulations, in particular for rhodopsin, which will be discussed below. Instead of using pseudopotentials, our method relies on a density-functional based tight-binding approach [9.13], which allows to treat a much larger region quantum mechanically with less computational efforts. In the latter approach the MM region is described by using CHARMM [9.14]. Compared to ab initio quantum chemical approaches, the
180 P.Entel et al. latter two methods have the advantage that a relatively large QM region can be treated with sufficient accuracy. The last remark concerns ground state and excited states properties of a quantum mechanical system in view of the normally used Born-Oppenheimer approximation to separate the dynamics of electrons and atoms.The dynam- ics of the electrons is then determined by the zero-temperature solution of the corresponding Schrodinger equation for the momentary distribution of atoms,while the atoms can be handled at finite temperatures in the molec- ular dynamics simulations.While this may be a reasonable approximation when simulating ground-state properties of close-packed metallic systems,its accuracy might be questioned when dealing with excited-state properties of biological systems.In additon the crossover from excited electronic states back to the ground sate requires a complicated mixing of basis sets.Both problems are so far not solved unambigiously.For a discussion of the exten- sion of conventional molecular dynamics for simulations not confined to the Born-Oppenheimer surface we refer to [9.15]. 9.2 Simulation of Biochemical Systems In the following we discuss two molecular systems,cyclodextrin interacting with bulk water and the protein rhodopsin with internal water,for which ex- emplary the molecular simulations help to deepen our understanding of the experimental data being available.We start the discussion with water which is the most abundant substance on earth and the only naturally occuring in- organic liquid.Due to the distinctive properties of H2O,its interaction with larger molecules,proteins or even the DNA is in many cases decisive for their properties:In the absence of water (and its hydrogen bonding network)pro- teins lack activity and the DNA loses its double-helical structure.For recent reviews concerning fundamental aspects of water (from bulk water to bio- logical,internal and surface waters,the phase diagram of ice,etc.)we refer to 9.16-9.19.The different dynamics of biological water compared to bulk water is,for example,determined by the protein-solvent interaction:Typi- cal consequences are multiple-time-scale processes,time-dependent diffusion coefficients and excess low-frequency vibrational modes (glassy behavior of water,see 9.19).In the context of methods concentrating on ab initio MD of liquids and applied to water in particular see,for example,[9.20.Much recent work is connected to the simulation of proton transfer and solvation of ions in water or high pressure ice [9.21-9.24].This is also of interest with respect to the pH value of different waters affecting organisms living in wa- ter (pH =-log10 H+],where H+]is the molar concentration (mol/L);pH ranges from 0 to 14,with 7 being neutral;pHs less than 7 are acidic while pHs greater than 7 are basic;normal rainwater has a pH value of 5.7)
180 P. Entel et al. latter two methods have the advantage that a relatively large QM region can be treated with sufficient accuracy. The last remark concerns ground state and excited states properties of a quantum mechanical system in view of the normally used Born-Oppenheimer approximation to separate the dynamics of electrons and atoms. The dynamics of the electrons is then determined by the zero-temperature solution of the corresponding Schr¨odinger equation for the momentary distribution of atoms, while the atoms can be handled at finite temperatures in the molecular dynamics simulations. While this may be a reasonable approximation when simulating ground-state properties of close-packed metallic systems, its accuracy might be questioned when dealing with excited-state properties of biological systems. In additon the crossover from excited electronic states back to the ground sate requires a complicated mixing of basis sets. Both problems are so far not solved unambigiously. For a discussion of the extension of conventional molecular dynamics for simulations not confined to the Born-Oppenheimer surface we refer to [9.15]. 9.2 Simulation of Biochemical Systems In the following we discuss two molecular systems, cyclodextrin interacting with bulk water and the protein rhodopsin with internal water, for which exemplary the molecular simulations help to deepen our understanding of the experimental data being available. We start the discussion with water which is the most abundant substance on earth and the only naturally occuring inorganic liquid. Due to the distinctive properties of H2O, its interaction with larger molecules, proteins or even the DNA is in many cases decisive for their properties: In the absence of water (and its hydrogen bonding network) proteins lack activity and the DNA loses its double-helical structure. For recent reviews concerning fundamental aspects of water (from bulk water to biological, internal and surface waters, the phase diagram of ice, etc.) we refer to [9.16–9.19]. The different dynamics of biological water compared to bulk water is, for example, determined by the protein-solvent interaction: Typical consequences are multiple-time-scale processes, time-dependent diffusion coefficients and excess low-frequency vibrational modes (glassy behavior of water, see [9.19]). In the context of methods concentrating on ab initio MD of liquids and applied to water in particular see, for example, [9.20]. Much recent work is connected to the simulation of proton transfer and solvation of ions in water or high pressure ice [9.21–9.24]. This is also of interest with respect to the pH value of different waters affecting organisms living in water (pH = −log10[H+], where [H+] is the molar concentration (mol/L); pH ranges from 0 to 14, with 7 being neutral; pHs less than 7 are acidic while pHs greater than 7 are basic; normal rainwater has a pH value of 5.7)
9 Molecular Dynamics Simulations in Biology,Chemistry and Physics 181 9.2.1 Molecular Dynamics Simulation of Liquid Water Useful tools to sample the structure and dynamics of water are the radial distribution function(RDF),the time correlation functions and diffusion co- efficients.Here we briefly discuss results for the RDFs defined by gaB(r)= ()( (9.7) which is related by Fourier transformation to the (intermolecular)structure factor, SaB(k)=1+p dr eikr[gaa(r)-1刂， (9.8) where a,B=O,H;i,j refer to the molecules (note that in the case of goH and gHn we have,in addition,intramolecular contributions from i=j); the X-ray data give the O-O pair correlations 9.25 while the neutron data provide all pair correlations including the intramolecular correlations 9.26); p=N/V,N is the number of water molecules.Sos(k)is directly related to X-ray or neutron scattering intensities. Figure 9.1 shows calculated RDFs of bulk water obtained from Car- Parrinello (CP)MD simulations 9.27 in comparison to the most recent results of X-ray 9.25 and neutron scattering experiments 9.26.In the simulations the standard CP method 9.28 has been employed within a plane-wave-basis density functional theory,using the BLYP [9.29]gradient- corrected exchange-correlation energy functional and the norm-conserving Troullier-Martins pseudopotential 9.30(with a plane-wave cutoff of 80 Ry for the "good"results and I point only in the Brillouin sampling;for further details and discussion of previous simulations see [9.27).There is remarkable agreement of the calculated RDFs with most recent both X-ray [9.25]and neutron [9.26]scattering experiments. We have checked the RDFs by simulations using two other methods, VASP [9.31]and DFTB [9.13].For the bulk density of water and room tem- perature we obtain the same good agreement with experiment as in the CP simulations when using VASP while the DFTB calculations yield a smaller and higher first peak in goo(r).However,the simulations also show that for slightly different densities than bulk density or insufficient system sizes and simulation times position and height of the first peak in goo(r)are not so well reproduced.One should also note that with respect to simulations using empirical water models in many cases a relatively high and sharp first peak, contrary to experiments,is obtained [9.25,9.26](for a survey of molecular models of water see [9.32).Thus ab initio MD simulations are very useful to obtain correct information about the structure(and its temperature depen- dence)of water,they are also useful for the interpretation of the experimental data
9 Molecular Dynamics Simulations in Biology, Chemistry and Physics 181 9.2.1 Molecular Dynamics Simulation of Liquid Water Useful tools to sample the structure and dynamics of water are the radial distribution function (RDF), the time correlation functions and diffusion coefficients. Here we briefly discuss results for the RDFs defined by gαβ(r) = 1 ρ2 5 i =j δ(riα)δ(rjβ − r) 6 , (9.7) which is related by Fourier transformation to the (intermolecular) structure factor, Sαβ(k)=1+ ρ d3r eikr [gαβ(r) − 1] , (9.8) where α, β = O, H; i, j refer to the molecules (note that in the case of gOH and gHH we have, in addition, intramolecular contributions from i = j); the X-ray data give the O–O pair correlations [9.25] while the neutron data provide all pair correlations including the intramolecular correlations [9.26]); ρ = N/V , N is the number of water molecules. Sαβ(k) is directly related to X-ray or neutron scattering intensities. Figure 9.1 shows calculated RDFs of bulk water obtained from CarParrinello (CP) MD simulations [9.27] in comparison to the most recent results of X-ray [9.25] and neutron scattering experiments [9.26]. In the simulations the standard CP method [9.28] has been employed within a plane-wave-basis density functional theory, using the BLYP [9.29] gradientcorrected exchange-correlation energy functional and the norm-conserving Troullier-Martins pseudopotential [9.30] (with a plane-wave cutoff of 80 Ry for the “good” results and Γ point only in the Brillouin sampling; for further details and discussion of previous simulations see [9.27]). There is remarkable agreement of the calculated RDFs with most recent both X-ray [9.25] and neutron [9.26] scattering experiments. We have checked the RDFs by simulations using two other methods, VASP [9.31] and DFTB [9.13]. For the bulk density of water and room temperature we obtain the same good agreement with experiment as in the CP simulations when using VASP while the DFTB calculations yield a smaller and higher first peak in gOO(r). However, the simulations also show that for slightly different densities than bulk density or insufficient system sizes and simulation times position and height of the first peak in gOO(r) are not so well reproduced. One should also note that with respect to simulations using empirical water models in many cases a relatively high and sharp first peak, contrary to experiments, is obtained [9.25, 9.26] (for a survey of molecular models of water see [9.32]). Thus ab initio MD simulations are very useful to obtain correct information about the structure (and its temperature dependence) of water, they are also useful for the interpretation of the experimental data
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
