《纺织复合材料》课程参考文献（Computational Materials Science，From Basic Principles to Material Properties）09 Molecular Dynamics Simulations in Biology, Chemistry and Physics

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

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 time￾reversal invariant trajectories. Using a different approach Nos´e reproduced the canonical phase-space distribution, so that the kinetic energy can fluctuate with a distribution pro￾portional 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 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 “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 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

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

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 dynam￾ics 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 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)

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 co￾efficients. 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 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 Γ 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