分子模拟第一章
《分子模拟方法》课件
加速研发进程
分子模拟可以大大缩短药 物研发、材料合成等领域 的实验周期，降低研发成 本。
揭示微观机制
通过模拟，可以揭示分子 间的相互作用机制和反应 过程，有助于深入理解物 质的性质和行为。
分子模拟的发展历程
经典力学模拟
基于牛顿力学，适用于 较大分子体系，但精度
较低。
量子力学模拟
适用于小分子体系，精 度高，但计算量大，需
详细描述
利用分子模拟方法，模拟小分子药物与生物大分子（如蛋白质、核酸等）的相 互作用过程，探究药物的作用机制和药效，为新药研发提供理论支持。
高分子材料的模拟研究
总结词
研究高分子材料的结构和性能，优化 材料的设计和制备。
详细描述
通过模拟高分子材料的结构和性能， 探究高分子材料的物理和化学性质， 优化材料的设计和制备过程，为新材 料的研发提供理论指导。
分子动力学方法需要较高的计算资源和 精度，但可以获得较为准确的结果，因 此在计算化学、生物学、材料科学等领
域得到广泛应用。
介观模拟的原理
介观模拟是一种介于微观和宏观之间的模拟方 法，通过模拟一定数量的粒子的相互作用和演 化来研究介观尺度的结构和性质。
介观模拟方法通常采用格子波尔兹曼方法、粒 子流体动力学等方法，适用于模拟流体、表面 、界面等介观尺度的问题。
分子模拟基于量子力学、经典力 学、蒙特卡洛等理论，通过建立 数学模型来描述分子间的相互作
用和运动。
分子模拟可以用于药物研发、材 料科学、环境科学等领域，为实 验研究和工业应用提供重要支持
。
分子模拟的重要性
01
02
03
预测分子性质
通过模拟，可以预测分子 的性质，如稳定性、溶解 度、光谱等，为实验设计 和优化提供指导。
分子模拟课件
ME346–Introduction to Molecular Simulations–Stanford University–Winter2007 Handout1.An Overview of Molecular SimulationWei Caic All rights reservedSeptember26,2005Contents1Examples of molecular simulations1 2What consists of a molecular dynamics?2 3Newton’s equation of motion2 4A simple numerical integrator:Verlet algorithm5 1Examples of molecular simulationsPictures and movies of modern large scale molecular simulations can be downloaded from the course web site /∼me346(click“Course Materials”).Exam-ples include plastic deformation in metals containing a crack,phase transformation induced by shock waves,nano-indentation of metalfilm,laser ablation on metal surface,etc.In each of these examples,molecular simulations not only produce beautiful pictures,but they also offer new understanding on the behavior of materials at the nano-scale.Many times such understanding provides new insights on how the materials would behave at the micro-scale and nano-scale.With unique capabilities to offer,molecular simulation has become a widely used tool for scientific investigation that complements our traditional way of doing science: experiment and theory.The purpose of this course is to provide thefirst introduction to molecular level simulations—how do they work,what questions can they answer,and how much can we trust their results—so that the simulation tools are no longer black boxes to us.Students will develop a set of simulation codes in Matlab by themselves.(A simple Matlab tutorial written by Sigmon[5]may help you refresh your memory.)We will also give tutorial and exercise problems based on MD++,a molecular simulation package used in my research group.MD++will be used for illustration in class as well as tools to solve homework problems.2What consists of a molecular dynamics?In a molecular simulation,we view materials as a collection of discrete atoms.The atoms interact by exerting forces on each other and they follow the Newton’s equation of motion. Based on the interaction model,a simulation compute the atoms’trajectories numerically. Thus a molecular simulation necessarily contains the following ingredients.•The model that describes the interaction between atoms.This is usually called the interatomic potential:V({r i}),where{r i}represent the position of all atoms.•Numerical integrator that follows the atoms equation of motion.This is the heart of the ually,we also need auxiliary algorithms to set up the initial and bound-ary conditions and to monitor and control the systems state(such as temperature) during the simulation.•Extract useful data from the raw atomic trajectory pute materials properties of interest.Visualization.We will start our discussion with a simple but complete example that illustrates a lot of issues in simulation,although it is not on the molecular scale.We will simulate the orbit of the Earth around the Sun.After we have understood how this simulation works(in a couple of lectures),we will then discuss how molecular simulations differ from this type of “macro”-scale simulations.3Newton’s equation of motionThroughout this course,we limit our discussions to classical dynamics(no quantum mechan-ics).The dynamics of classical objects follow the three laws of Newton.Let’s review the Newton’s laws.Figure1:Sir Isaac Newton(1643-1727England).First law:Every object in a state of uniform motion tends to remain in that state of motion unless an external force is applied to it.This is also called the Law of inertia.Second law:An object’s mass m,its acceleration a,and the applied force F are related byF=m a.The bold font here emphasizes that the acceleration and force are vectors.Here they point to the same direction.According to Newton,a force causes only a change in velocity(v)but is not needed to maintain the velocity as Aristotle claimed(erroneously).Third law:For every action there is an equal and opposite reaction.The second law gives us the equation of motion for classical particles.Consider a particle (be it an atom or the Earth!depending on which scale we look at the world),its position is described by a vector r=(x,y,z).The velocity is how fast r changes with time and is also a vector:v=d r/dt=(v x,v y,v z).In component form,v x=dx/dt,v y=dy/dt,v z=dz/dt. The acceleration vector is then time derivative of velocity,i.e.,a=d v/dt.Therefore the particle’s equation of motion can be written asd2r dt2=Fm(1)To complete the equation of motion,we need to know how does the force F in turn depends on position r.As an example,consider the Earth orbiting around the Sun.Assume the Sun isfixed at origin(0,0,0).The gravitational force on the Earth is,F=−GmM|r|2·ˆe r(2)whereˆe r=r/|r|is the unit vector along r direction.|r|=√x2+y2+z2is the magnitudeof r.m is the mass of the Earth(5.9736×1024kg),M is mass of the Sun(1.9891×1030kg), and G=6.67259×10−11N·m2/kg2is gravitational constant.If we introduce a potential field(the gravitationalfield of the Sun),V(r)=−GmM|r|(3)Then the gravitational force can be written as minus the spatial derivative of the potential,F=−dV(r)d r(4)In component form,this equation should be interpreted as,F x=−dV(x,y,z)dxF y=−dV(x,y,z)dyF z=−dV(x,y,z)dzNow we have the complete equation of motion for the Earth:d2r dt2=−1mdV(r)d r(5)in component form,this reads,d2x dt2=−GMx|r|3d2y dt2=−GMy|r|3d2z dt2=−GMz|r|3In fact,most of the system we will study has a interaction potential V(r),from which the equation of motion can be specified through Eq.(5).V(r)is the potential energy of the system.On the other hand,the kinetic energy is given by the velocity,E kin=12m|v|2(6)The total energy is the sum of potential and kinetic energy contributions,E tot=E kin+V(7) When we express the total energy as a function of particle position r and momentum p=m v, it is called the Hamiltonian of the system,H(r,p)=|p|22m+V(r)(8)The Hamiltonian(i.e.the total energy)is a conserved quantity as the particle moves.To see this,let us compute its time derivative,dE tot dt =m v·d vdt+dV(r)d r·d rdt=md rdt·−1mdV(r)d r+dV(r)d r·d rdt=0(9)Therefore the total energy is conserved if the particle follows the Newton’s equation of motion in a conservative forcefield(when force can be written in terms of spatial derivative of a potentialfield),while the kinetic energy and potential energy can interchange with each other.The Newton’s equation of motion can also be written in the Hamiltonian form,d r dt =∂H(r,p)∂p(10)d p dt =−∂H(r,p)∂r(11)Hence,dH dt =∂H(r,p)∂r·d rdt+∂H(r,p)∂p·d pdt=0(12)Figure2:Sir William Rowan Hamilton(1805-1865Ireland).4A simple numerical integrator:Verlet algorithmA dynamical simulation computes the particle position as a function of time(i.e.its trajec-tory)given its initial position and velocity.This is called the initial value problem(IVP). Because the Newton’s equation of motion is second order in r,the initial condition needs to specify both particle position and velocity.To solve the equation of motion on a computer,thefirst thing we need is to discretize the time.In other words,we will solve for r(t)on a series of time instances t ually the time axis is discretized uniformly:t i=i·∆t,where∆t is called the time step.Again,the task of the simulation algorithm is tofind r(t i)for i=1,2,3,···given r(0)and v(0).The Verlet algorithm begins by approximating d2r(t)/dt2,d2r(t) dt2≈r(t+∆t)−2r(t)+r(t−∆t)∆t2(13)Thusr(t+∆t)−2r(t)+r(t−∆t)∆t2=−1mdV(r(t))dr(14)r(t+∆t)=2r(t)−r(t−∆t)−∆t·1m·dV(r(t))dr(15)Therefore,we can solve r(t+∆t)as along as we know r(t)and r(t−∆t).In other words,if we know r(0)and r(−∆t),we can solve for all r(t i)for i=1,3,4,···.Question:Notice that to start the Verlet algorithm we need r(0)and r(−∆t).However, the initial condition of a simulation is usually given in r(0)and v(0).How do we start the simulation when this initial condition is specified?References1.Alan and Tildesley,Computer Simulation of Liquids,(Oxford University Press,1987)pp.71-80.2.Frenkel and Smit,Understanding Molecular Simulation:From Algorithms to Applica-tions,2nd ed.(Academic Press,2002).3.Bulatov and Cai,Computer Simulation of Dislocations,(Oxford University Press,2006)pp.49-51.ndau and Lifshitz,Mechanics,(Butterworth-Heinemann,1976).5.Kermit Sigmon,Matlab Primer,Third Edition,/Matlab/matlab-primer.pdfME346–Introduction to Molecular Simulations–Stanford University–Winter2007 Handout2.Numerical Integrators IWei Caic All rights reservedJanuary12,2007Contents1Which integrator to use?1 2Several versions of Verlet algorithm2 3Order of accuracy3 4Other integrators4 1Which integrator to use?In the previous lecture we looked at a very simple simulation.The Newton’s equation of motion is integrated numerically by the original Verlet algorithm.There are many algorithms that can be used to integrate an ordinary differential equation(ODE).Each has a its own advantages and disadvantages.In this lecture,we will look into more details at the numerical integrators.Our goal is to be able to make sensible choices when we need to select a numerical integrator for a specific simulation.What questions should we ask when choosing an integrator?Obviously,we would like something that is both fast and accurate.But usually there are more issues we should care about.Here are a list of the most important ones(Allen and Tildesley1987):•Order of accuracy•Range of stability•Reversibility•Speed(computation time per step)•Memory requirement(as low as possible)•Satisfying conservation law•Simple form and easy to implementItem1and2(especially2)are important to allow the use of long time steps,which makes the simulation more efficient.2Several versions of Verlet algorithmThe Verlet algorithm used in the previous lecture can be written as,r(t+∆t)=2r(t)−r(t−∆t)+a(t)·∆t2(1)wherea(t)=−1mdV(r(t))d r(t)(2)While this is enough to compute the entire trajectory of the particle,it does not specify the velocity of the particle.Many times,we would like to know the particle velocity as a function of time as well.For example,the velocity is needed to compute kinetic energy and hence the total energy of the system.One way to compute velocity is to use the followingapproximation,v(t)≡d r(t)dt≈r(t+∆t)−r(t−∆t)2∆t(3)Question:if the initial condition is given in terms of r(0)and v(0),how do we use the Verlet algorithm?Answer:The Verlet algorithm can be used to compute r(n∆t)for n=1,2,···provided r(−∆t)and r(0)are known.The task is to solve for r(−∆t)given r(0)and v(0).This is achieved from the following two equations,r(∆t)−r(−∆t)2∆t=v(0)(4)r(∆t)−2r(0)+r(−∆t)∆t2=−1mdV(r(t))d r(t)t=0(5)Exercise:express r(−∆t)and r(∆t)in terms of r(0)and v(0).Therefore,r(−∆t)and r(∆t)are solved together.During the simulation,we only know v(t)after we have computed(t+∆t).These are the(slight)inconveniences of the(original) Verlet algorithm.As we will see later,in this algorithm the accuracy of velocity calculation is not very high(not as high as particle positions).Several modifications of the Verlet algorithms were introduced to give a better treatment of velocities.The leap-frog algorithm works as follows,v(t+∆t/2)=v(t−∆t/2)+a(t)·∆t(6)r(t+∆t)=r(t)+v(t+∆t/2)·∆t(7) Notice that the velocity and position are not stored at the same time slices.Their values are updated in alternation,hence the name“leap-frog”.For a nice illustration see(Allen and Tildesley1987).This algorithm also has the draw back that v and r are not known simultaneously.To get velocity at time t,we can use the following approximation,v(t)≈v(t−∆t/2)+v(t+∆t/2)2(8)The velocity Verlet algorithm solves this problem in a better way.The algorithm reads,r(t+∆t)=r(t)+v(t)∆t+12a(t)∆t2(9)v(t+∆t/2)=v(t)+12a(t)∆t(10)a(t+∆t)=−1mdV(r(t+∆t))d r(t+∆t)(11)v(t+∆t)=v(t+∆t/2)+12a(t+∆t)∆t(12)Exercise:show that both leap-frog and velocity Verlet algorithms gives identical results as the original Verlet algorithm(except forfloating point truncation error).3Order of accuracyHow accurate are the Verlet-family algorithms?Let’s take a look at the original Verletfirst. Consider the trajectory of the particle as a continuous function of time,i.e.r(t).Let us Taylor expand this function around a given time t.r(t+∆t)=r(t)+d r(t)dt∆t+12d2r(t)dt2∆t2+13!d3r(t)dt3∆t3+O(∆t4),(13)where O(∆t4)means that error of this expression is on the order of∆t4.In other words, (when∆t is sufficiently small)if we reduce∆t by10times,the error term would reduce by 104times.Replace∆t by−∆t in the above equation,we have,r(t−∆t)=r(t)−d r(t)dt∆t+12d2r(t)dt2∆t2−13!d3r(t)dt3∆t3+O(∆t4),(14)Summing the two equations together,we have,r(t+∆t)+r(t−∆t)=2r(t)+d2r(t)dt2∆t2+O(∆t4)≡2r(t)+a(t)∆t2+O(∆t4)(15)r(t+∆t)=2r(t)−r(t−∆t)+a(t)∆t2+O(∆t4)(16) Therefore,the local accuracy for the Verlet-family algorithms is to the order of O(∆t4). This is quite accurate for such a simple algorithm!The trick here is that,by adding t+∆t and t−∆t terms,the third order term cancels out.Notice that the accuracy analysis here deals with the local error.That is,assuming r(t)is exact,how much error we make when computing r(t+∆t).In practice,the total error of computed r at time t is the accumulated (and possibly amplified)effect of all time steps before time t.Because for a smaller timestep ∆t more steps is required to reach a pre-specified time t,the global order of accuracy is always local order of accuracy minus one.For example,the global order of accuracy for Verlet algorithms is3.Since all Verlet-family algorithms generate the same trajectory,they are all O(∆t4)in terms of particle positions.But they differ in velocity calculations.How accurate is the velocity calculation in the original Verlet algorithm?Let us subtract Eq.(13)and(14),r(t+∆t)−r(t−∆t)=2v(t)∆t+13d3r(t)dt3∆t3+O(∆t4),(17)v(t)=12∆tr(t+∆t)−r(t−∆t)−16d3r(t)dt2∆t3+O(∆t3)(18)v(t)=r(t+∆t)−r(t−∆t)2∆t+O(∆t2)(19)Thus the velocity is only accurate to the order of O(∆t2)!This is much worse than O(∆t4) of the positions.A reason for the lose of accuracy is that the velocities are computed by taking the difference between two large numbers,i.e.r(t+∆t)and r(t−∆t),and then divide by a small number.1Intuitively,higher order algorithms allows the use of bigger time steps to achieve the same level of accuracy.However,the more important determining factor for the choice of time steps is usually not the order of accuracy,but instead is the algorithm’s numerical stability (see next section).4Other integratorsBefore closing this lecture,let us briefly look at two other well known integrators.Gear’s predictor-corrector method is a high-order integrator.It keeps track of several time deriva-tives of the particle position.Definer0(t)≡r(t)(20) 1Strictly speaking,Eq.(17)is valid when we have exact values for r(t+∆t),r(t)and r(t−∆t).But in doing the order-of-accuracy analysis,we can only assume that we have exact values for r(t)and r(t−∆t) and we know that we are making an error in r(t+∆t).Fortunately,the error we make in r(t+∆t)is also O(∆t4).Thus Eq.(17)still holds.r1(t)≡∆t d r(t)dt(21)r2(t)≡∆t22!d2r(t)dt2(22)r3(t)≡∆t33!d3r(t)dt3(23)···Sincer(t+∆t)=r(t)+∆t d r(t)dt+∆t22!d2r(t)dt2+∆t33!d3r(t)dt3+O(∆t4)=r(t)+r1(t)+r2(t)+r3(t)+O(∆t4)r1(t+∆t)=r1(t)+2r2(t)+3r3(t)+O(∆t4)r2(t+∆t)=r2(t)+3r3(t)+O(∆t4)r3(t+∆t)=r3(t)+O(∆t4)we can predict the values of r and r i at time t+∆t simply based on their values at time t,r P(t+∆t)=r(t)+r1(t)+r2(t)+r3(t)r P1(t+∆t)=r1(t)+2r2(t)+3r3(t)r P2(t+∆t)=r2(t)+3r3(t)r P3(t+∆t)=r3(t)Notice that we have not even calculated the force at all!If we compute the force at time t+∆t based on the predicted position r P(t+∆t),we will obtain the second derivative of r at time t+∆t—let it be r C(t+∆t).Mostly likely,r C(t+∆t)will be different from the predicted value r P(t+∆t).Thefinal results are obtained by adding a correction term to each predictor,r n(t+∆t)=r Pn (t+∆t)+C n[r C2(t+∆t)−r P2(t+∆t)](24)The Gear’s6th order scheme(uses r0through r5)takes the following C n values,C0=3 16C1=251 360C2=1C3=11 18C4=1 6C5=1 60The global accuracy is of the order O(∆t)6,hence the method is called the Gear’s6th order algorithm.Exercise:show that the local truncation error of the above scheme is of the order of O (∆t 7).In a way,the Runge-Kutta method is similar to predictor-corrector methods in that the final answer is produced by improving upon previous predictions.However,the Runge-Kutta method requires multiple evaluation of forces.The Runge-Kutta method is most conveniently represented when the equation of motion is written as a first order differential equation,˙η=g (η)≡ω∂H (η)∂η(25)whereη≡(r ,p )(26)andω= 0I−I 0 (27)and I is the 3×3identity matrix and H is the Hamiltonian.The fourth-order Runge-Kutta method goes as (from ),k 1=∆t ·g (η(t ))k 2=∆t ·g η(t )+12k 1 k 3=∆t ·g η(t )+12k 2 k 4=∆t ·g (η(t )+k 3)η(t +∆t )=η(t )+16k 1+13k 2+13k 3+16k 4(28)Exercise:show that the local truncation error of the above scheme is of the order of O (∆t 5)and the global accuracy is of the order of O (∆t 4).Suggested reading1.(Alan and Tildesley 1987)pp.71-84.2.J.C.Butcher,The numerical analysis of ordinary differential equations:Runge-Kutta and general linear methods ,(Wiley,1987),pp.105-151.ME346–Introduction to Molecular Simulations–Stanford University–Winter2007 Handout3.Symplectic ConditionWei Caic All rights reservedJanuary17,2007Contents1Conservation laws of Hamiltonian systems1 2Symplectic methods4 3Reversibility and chaos5 4Energy conservation7 1Conservation laws of Hamiltonian systemsThe numerical integrators for simulating the Newton’s equation of motion are solvers of ordinary differential equations(ODE).However,the Newton’s equation of motion that can be derived from a Hamiltonian is a special type of ODE(not every ODE can be derived from a Hamiltonian).For example,a Hamiltonian system conserves it total energy,while an arbitrary ODE may not have any conserved quantity.Therefore,special attention is needed to select an ODE solver for Molecular Dynamics simulations.We need to pick those ODE solvers that obeys the conservation laws of a Hamiltonian system.For example,if an ODE solver predicts that the Earth will drop into the Sun after100cycles around the Sun,it is certainly the ODE solver to blame and we should not immediately start worrying about the fate of our planet based on this prediction.A conservation law is a(time invariant)symmetry.A Hamiltonian system possess more symmetries than just energy conservation.Ideally the numerical integration scheme we pick should satisfy all the symmetries of the true dynamics of the system.We will look at these symmetries in this and following sections.When discussing a Hamiltonian system,it is customary to use letter q to represent the particle position(instead of r).Recall that the equation of motion for a system with Hamiltonian H(q,p)is,˙q=∂H(q,p)∂p(1)˙p=−∂H(q,p)∂q,(2)where˙q≡d q/dt,˙p≡d p/dt.If the system has N particles,then q is a3N-dimensional vector specifying its position and p is a3N-dimensional vector specifying its momentum. (In the above Earth orbiting example,N=1.)Notice that in terms of q and p,the equation of motion is a set of(coupled)first-order differential equations(whereas F=m a is second order).The equation can be reduced to an even simpler form if we define a6N-dimensional vectorη≡(q,p)(3) Then the equation of motion becomes,˙η=ω∂H(η)∂η(4)whereω=0I−I0(5)and I is the3×3identity matrix.The6N-dimensional space whichη≡(q,p)belongs to is called the phase space.The evolution of system in time can be regarded as the motion of a pointηin the6N-dimensional phase space following thefirst-order differential equation(4).(a)(b)Figure1:(a)Motion of a point(q,p)in phase space is confined to a shell with constant energy.(b)Area of a phase space element(shaded)is conserved.Because the Newton’s equation of motion conserves total energy,the motion of a point in phase space must be confined to a subspace(or hyperplane)with constant energy,i.e.a constant energy shell,as illustrated in Fig.1(a).Think about all the points in the constant energy shell.As time evolves,the entire volume that these points occupy remains the same —the points simply rearrange within the constant energy shell.Let us now consider the evolution of a small element in phase space over time,as illus-trated in Fig.1(b).Think of this as the ensemble of trajectories of many Molecular Dynamics simulations with very similar initial conditions.We will show that the area enclosed by this element remains constant,even though the element inevitably experiences translation anddistortion.1Let the element at time t be a hypercube around pointη,whose area is,|dη|=|d q|·|d p|=dq1···dq3N dp1···dp3N(6) For a pointηat time t,letξbe its new position at time t+δt.In the limit ofδt→0,wehave,ξ=η+˙ηδt=η+ω∂H(η)∂ηδt(7)Because the above equation is true only in the limit ofδt→0,we will ignore any higher order terms(e.g.δt2)in the following discussions.Let M be the Jacobian matrix of the transformation fromηtoξ,i.e.,M≡∂ξ∂η=1+ω∂2H(η)∂η∂ηδt(8)For clarity,we write M explicitly in terms of q and p,M=1+0I−I0·∂2H∂q∂q∂2H∂q∂p∂2H∂p∂q∂2H∂p∂p·δt(9)=1+∂2Hδt∂2Hδt−∂2H∂q∂qδt1−∂2H∂q∂pδt(10)The area of the new element is related to the the determinant of the Jacobian matrix,|dξ|=|dη|·|det M|=|dη|·(1+O(δt2))(11) Because thefirst order term ofδt in det M vanishes,we can show that the area of the element remains constant for an arbitrary period of time(i.e.forever).2This is an important property of a Hamiltonian system.Because the area of any element in phase space always remains constant,the evolution of phase space points is analogous to that in an incompressiblefluid.Exercise:Write down the explicit expression for the Jacobian matrix M for problem of the Earth orbiting around the Sun.The Hamiltonian dynamics has even more symmetries.Notice that the transpose of M is,M T=1−∂2H∂η∂ηωδt(12)1It then follows that the area enclosed by an arbitrarily shaped element in the phase space also remains constant.2To show that area remains constant after afinite period∆t,we only need to divide the time interval into N subintervals,each withδt=∆t/N.The area change per subinterval is of the order of1/N2.The accumulated area change over time period∆t is then of the order of1/N,which vanishes as N→∞.we haveM ωM T = 1+ω∂2H ∂η∂ηδt ω 1−∂2H ∂η∂ηωδt =ω+ω∂2H ∂η∂ηωδt −ω∂2H ∂η∂ηωδt +O (δt 2)(13)=ω+O (δt 2)(14)(15)Therefore the Jacobian matrix M thus satisfies the so called symplectic condition (up to O (δt 2)),M ωM T =ω(16)Again,we can show that the symplectic condition is satisfied for an arbitrary period of time.To be specific,let ξbe the point in phase space at time t ;at time 0this point was at η.ξcan then be regarded as a function of t and η(initial condition),i.e.ξ=ξ(η,t ).Define M as the Jacobian matrix between ξand η,we can show that M satisfies the symplectic condition Eq.(16).3Obviously,the Hamiltonian dynamics is also time reversible.This means that if ηevolves to ξafter time t .Another phase space point η that has the same q as ηbut with reversed p (momentum)will exactly come to point ξafter time t .To summarize,the dynamics of a system that has a Hamiltonian has the following symmetry properties:•Conserves total energy•Conserves phase space area (incompressible flow in phase space)•Satisfies symplectic condition Eq.(16)•Is time reversibleIdeally,the numerical integrator we choose to simulate Hamiltonian system should satisfy all of these symmetries.In the following,let us take a look at how specific integrators are doing in this regard.In fact,if an algorithm satisfies the symplectic condition,it automatically conserves phase space area 4,is reversible,and should have good long-term energy conserva-tion properties.This is why symplectic methods are good choices for simulating Hamiltonian systems.2Symplectic methodsWhat do we mean if we ask whether a numerical integrator satisfies the symplectic condition or not?For a numerical integrator,the system moves forward in time by discrete jumps.Let3Again,we devide the time period t into N subintervals,each with δt =∆t/N .Let M i be the Jacobian matrix from time (i −1)δt to iδt .We have M = N i =1M i .We can show that M T ωM =ωin the limit of N →∞.4Since M T ωM =ω,det ω=det(M T ωM )=(det M )2det ω.Thus det M =±1,which means area conservation.η(i)=(q(i),p(i))be the place of the system in the phase space at time step i,i.e.t i=i∆t. At the next time step,the algorithm brings the system to pointη(i+1)=(q(i+1),p(i+1)). Let M be the Jacobian matrix that connectsη(i)withη(i+1).The method is symplectic if MωM T=ω.The good news is that the Verlet-family algorithms are symplectic.To show this is the case,let us consider the leap-frog algorithm.Let q(i)correspond to q(i·∆t)and p(i) correspond to p((i−12)·∆t).Then the leap-frog algorithm is,p(i+1)=p(i)−∂V∂q(i)∆t(17)q(i+1)=q(i)+p(i+1)m∆t=q(i)+1mp(i)−∂V∂q(i)∆t∆t(18)Thus,M=∂q(i+1)∂q(i)∂q(i+1)∂p(i)∂p(i+1)(i)∂q(i+1)(i)=1−1m∂2V∂q(i)∂q(i)∆t2∆tm−∂2V∂q(i)∂q(i)∆t1(19)Thus M is symplectic.Exercise:show that for M defined in Eq.(19),MωM T=ω(symplectic)and det M=1 (area preserving).Exercise:derive the Jacobian matrix M for velocity Verlet algorithm and show that it is symplectic and area preserving.3Reversibility and chaosBecause the Verlet-family algorithms(original Verlet,leap-frog and velocity Verlet)are sym-plectic,they are also time reversible.Nonetheless,it is instructive to show the time reversibil-ity explicitly.If we run the leap-frog algorithm in reverse time,then the iteration can be written as,˜q(i+1)=˜q(i)+˜p(i)m∆t(20)˜p(i+1)=˜p(i)−∂V∂q(i+1)∆t(21)(22)where i is the iteration number,˜q(i)=q(−i∆t),˜p(i)=p((−i−12)∆t).It is not difficult toshow that if˜q(i)=q(i+1)and˜p(i)=−p(i+1),then˜q(i+1)=q(i)and˜p(i+1)=−p(i).Exercise:show that the leap-frog algorithm is reversible.Exercise:show that the velocity Verlet algorithm is reversible.Exercise:show that the Forward Euler method:q(i+1)=q(i)+∆t·p(i)/m,p(i+1)=p(i)−∆t·∂V/∂q(i),is not reversible and not symplectic.Since every integration step using the Verlet algorithm is time reversible,then in principle an entire trajectory of Molecular Dynamics simulation using the Verlet algorithm should be time reversible.To be specific,let the initial condition of the simulation be q(0)and p(0). After N steps of integration,the system evolves to point q(N)and p(N).Time reversibility means that,if we run the reverse Verlet iteration starting from q(N)and−p(N),after N steps we should get exactly to point q(0)and−p(0).In practice,however,we never have perfect reversibility.This is due to the combined effect of numerical round offerror and the chaotic nature of trajectories of many Hamiltonian systems.The chaotic nature of many Hamiltonian system can be illustrated by considering the two trajectories with very close initial conditions.Let the distance between the two points in phase space at time0be d(0).At sufficiently large t,the two phase space points will diverge exponentially,i.e.d(t)∼d(0)eλt,whereλis called the Lyapunov exponent and this behavior is called Lyapunov instability.Given that we can only represent a real number on a computer withfinite precision,we are in effect introducing a small(random) perturbation at every step of the integration.Therefore,sooner or later,the numerical trajectory will deviate from the analytical trajectory(if the computer had infinite precision) significantly5.While the analytical trajectory is reversible,the one generated by a computer is not.This is why there exist a maximum number of iterations N,beyond which the original state cannot be recovered by running the simulation backwards.This by itself does not present too big a problem,since we seldom have the need to recover the initial condition of the simulation in this way.However,this behavior does mean that we will not be able to follow the“true”trajectory of the original system forever on a computer.Sooner or later, the“true”trajectory and the simulated one will diverge significantly.This is certainly a problem if we would like to rely entirely on simulations to tell us where is our satellite or space explorer.However,this usually does not present a problem in Molecular Simulations, in which we are usually not interested in the exact trajectories of every atom.We will return to this point in future lectures when we discuss Molecular Dynamics simulations involving many atoms.5Notice that the analytical trajectory here is not the same as the“true”trajectory of the original system. The analytical trajectory is obtained by advancing the system in discrete jumps along the time axis.。
分子模拟教案—
教案课程名称分子模拟授课题目(章、节) 第一章授课教师李慎敏授课对象化学081-2 授课时间9月2日教学方法讲授选用教具多媒体本教案以讲授一个单元(2—4学时)或一次实验（实习）为单位填写。
填写要用钢笔，字迹要清晰、工整。
按表项目逐一填写。
教案课程名称分子模拟授课题目(章、节) 第二章授课教师李慎敏授课对象化学081-2 授课时间9月9日教学方法讲授选用教具多媒体本教案以讲授一个单元(2—4学时)或一次实验（实习）为单位填写。
填写要用钢笔，字迹要清晰、工整。
按表项目逐一填写。
教案课程名称分子模拟授课题目(章、节) 第三章授课教师李慎敏授课对象化学081-2 授课时间9月16 日教学方法讲授选用教具多媒体本教案以讲授一个单元(2—4学时)或一次实验（实习）为单位填写。
填写要用钢笔，字迹要清晰、工整。
按表项目逐一填写。
教案课程名称分子模拟授课题目(章、节) 第三章授课教师李慎敏授课对象化学081-2 授课时间9月23 日教学方法讲授选用教具多媒体本教案以讲授一个单元(2—4学时)或一次实验（实习）为单位填写。
填写要用钢笔，字迹要清晰、工整。
按表项目逐一填写。
教案课程名称分子模拟授课题目(章、节) 第四章授课教师李慎敏授课对象化学081-2 授课时间10月9日教学方法讲授选用教具多媒体本教案以讲授一个单元(2—4学时)或一次实验（实习）为单位填写。
填写要用钢笔，字迹要清晰、工整。
按表项目逐一填写。
教案课程名称分子模拟授课题目(章、节) 第四章授课教师李慎敏授课对象化学081-2 授课时间10月14日教学方法讲授选用教具多媒体本教案以讲授一个单元(2—4学时)或一次实验（实习）为单位填写。
填写要用钢笔，字迹要清晰、工整。
按表项目逐一填写。
教案课程名称分子模拟授课题目(章、节) 第五章授课教师李慎敏授课对象化学081-2 授课时间10月21 日教学方法讲授选用教具多媒体本教案以讲授一个单元(2—4学时)或一次实验（实习）为单位填写。
《分子模拟设计》课件
目录
• 分子模拟设计概述 • 分子模拟设计的基本方法 • 分子模拟设计的应用领域 • 分子模拟设计的挑战与展望 • 分子模拟设计案例分析
01
CATALOGUE
分子模拟设计概述
定义与特点
定义
分子模拟设计是指利用计算机模 拟技术，对分子结构和性质进行 预测和设计的过程。
蒙特卡洛模拟
总结词
基于概率统计的模拟方法
详细描述
蒙特卡洛模拟是一种基于概率统计的模拟方法，通过随机抽样和统计计算来获 得系统的性质。该方法适用于模拟复杂系统，能够考虑系统的随机性和不确定 性。
分子力学模拟
总结词
基于势能函数的模拟方法
详细描述
分子力学模拟是一种基于势能函数的模拟方法，通过势能函数来描述分子间的相互作用和分子结构。该方法适用 于快速计算分子的结构和性质，常用于药物设计和材料科学等领域。
材料的界面行为等多个方面。
高分子材料的模拟设计有助于缩短新材料研发周期、 降低研发成本，提高新材料开发的成功率。
高分子材料的模拟设计是利用分子模拟技术对 高分子材料的结构和性质进行预测和优化的一 种方法。
通过模拟高分子材料的结构和性质，可以预测材 料的力学性能、热性能、电性能等，从而优化材 料的设计和制备工艺。
在生物大分子模拟中，研究人员可以使用分子模拟设计来研究蛋白质、 核酸和糖等生物大分子的结构和动力学性质。
这有助于理解这些大分子在细胞中的功能和相互作用的机制，以及与疾 病相关的生物大分子的异常行为。
04
CATALOGUE
分子模拟设计的挑战与展望
计算资源的限制
计算资源不足
高性能计算机和计算集群的资源有限，难以满足 大规模分子模拟的需求。
分子动力学模拟I
Gromacs中文教程淮海一粟分子动力学（MD）模拟分为三步：首先，要准备好模拟系统；然后，对准备好的系统进行模拟；最后，对模拟结果进行分析。
虽然第二步是最耗费计算资源的，有时候需要计算几个月，但是最耗费体力的步骤在于模拟系统准备和结果分析。
本教程涉及模拟系统准备、模拟和结果分析。
一、数据格式处理准备好模拟系统是MD最重要的步骤之一。
MD模拟原子尺度的动力学过程，可用于理解实验现象、验证理论假说，或者为一个待验证的新假说提供基础。
然而，对于上述各种情形，都需要根据实际情况对模拟过程进行设计；这意味着模拟的时候必须十分小心。
丢失的残基、原子和非标准基团本教程模拟的是蛋白质。
首先需要找到蛋白质序列并选择其起始结构，见前述；然后就要检查这个结构是否包含所有的残基和原子，这些残基和原子有时候也是模拟所必需的。
本教程假定不存在缺失，故略去。
另一个需要注意的问题是结构文件中可能包含非标准残基，被修饰过的残基或者配体，这些基团还没有力场参数。
如果有这些基团，要么被除去，要么就需要补充力场参数，这牵涉到MD的高级技巧。
本教程假定所有的蛋白质不含这类残基。
结构质量对结构文件进行检查以了解结构文件的质量是一个很好的练习。
例如，晶体结构解析过程中，对于谷氨酰胺和天冬酰胺有可能产生不正确的构象；对于组氨酸的质子化状态和侧链构象的解析也可能有问题。
为了得到正确的结构，可以利用一些程序和服务器(如WHATIF)。
本教程假定所用的结构没有问题，我们只进行数据格式处理。
二、结构转换和拓扑化一个分子可以由各个原子的坐标、键接情况与非键相互作用来确定。
由于.pdb结构文件只含有原子坐标，我们首先必须建立拓扑文件，该文件描述了原子类型、电荷、成键情况等信息。
拓扑文件对应着一种力场，选择何种力场对于拓扑文件的建立是一个值得仔细考虑的问题。
这里我们用的是GROMOS96 53a6连接原子力场，该力场对于氨基酸侧链的自由能预测较好，并且与NMR试验结果较吻合。
《分子模拟教程》课件
人工智能和机器学习技术将在分子模拟中发挥越 来越重要的作用，例如用于优化模拟参数、预测 性质等。
多尺度模拟
目前分子模拟主要集中在原子或分子级别，未来 将进一步发展多尺度模拟方法，将微观尺度和宏 观尺度相结合，以更全面地理解物质性质和行为 。
跨学科融合
分子模拟将与生物学、医学、材料科学等更多学 科领域进行交叉融合，为解决实际问题提供更多 可能性。
环境科学
在环境科学领域，分子模拟可用于研究污 染物在环境中的迁移转化机制，为环境保 护提供理论依据。
THANKS.
分子动力学模拟的常见算法
Verlet算法
一种基于离散时间步长的算法，用于计算分子位置和速度。
leapfrog算法
一种常用的分子动力学模拟算法，具有数值稳定性和计算效率高的特 点。
Parrinello-Rahman算法
一种基于分子力场的算法，可以用于模拟大尺度分子体系的运动。
Langevin动力学算法
材料科学
通过模拟材料中分子的运动和相互作 用，可以研究材料的力学、热学和电 学等性质，为材料设计和优化提供依 据。
03
Monte Carlo模拟
Monte Carlo模拟的基本概念
随机抽样
Monte Carlo模拟基于随 机抽样的方法，通过大量 随机样本的统计结果来逼 近真实结果。
概率模型
Monte Carlo模拟建立概 率模型，模拟系统的状态 变化和行为。
通过模拟药物分子与靶点分子的相互作用，预测 药物活性并优化药物设计。
材料科学
研究材料中分子的结构和性质，预测材料的物理 和化学性质。
生物大分子模拟
模拟生物大分子的结构和动力学行为，如蛋白质 、核酸等，有助于理解其功能和性质。
分子模拟PPT—第一章 概论
对1998 年诺贝尔化学奖 划时代的评价
瑞典皇家科学院的评价：
“ ··量子化学已发展成为广大化学家都能使用 · 的工具，将化学带入一个新时代 — 实验 与理论能携手协力揭示分子体系的性质。 化学不再是一门纯实验科学了”
I2 在 光 滑 球 模 型
不 Ar 同溶 半剂 径中 光的 滑振 球动 内弛 豫 时 间
Radius / nm 1.2 1.5 1.8 2.2 Shift / cm-1 1.20 1.79 1.76 1.75
受限于不同半径的光滑球内I2在Ar溶剂中的振动光谱位移
单 壁 碳 纳 米 管 模 型
Radius / nm Shift / cm-1
0.68 (10,10)
1.02 (15,15) 1.36 (20,20) 1.70 (25,25) 2.04 (30,30)
3.42
3.53 3.54 3.55 3.55
受限于不同半径的碳纳米管内I2在 Ar溶剂中的振动光谱位移密度为0.5 g/cm3
纳米反应器
自然界生命体系中的化学变化 绝佳的反应环境
R
+
R
R
R R
product shape selectivity
• “三十年前，如果说并非大多数化学家，那末至少 是有许多化学家嘲笑量子化学研究，认为这些工 作对化学用处不大，甚至几乎完全无用。现在的 情况却是完全两样了…。当90年代行将结束之际， 我们看到化学理论和计算研究的巨大进展，导致 整个化学正在经历一场革命性的变化。Kohn和 Pople是其中的两位最优秀代表”
1986：李远哲：“ 在十五年前，如果理论 结果与实验有矛盾，那么经常证明是理论结果错 了。但是最近十年则相反，常常是实验错了。… 量子力学有些结果是实验工作者事先未想到的，
结构生物学大分子模拟课件
现在学习的是第12页，共69页
到2007年3月，已知的蛋白质序列超过60万 条，而测定了三维结构的蛋白质仅为4万多个
蛋白质三维结构的测定已经成为生命科学 发展的“瓶颈”
现在学习的是第13页，共69页
蛋白质结构的基本概念
蛋白质的二级结构单元 – a螺旋 – b折叠股（b-strand）和b折叠片（b-sheet）
Turns Random coil (its neither!)
现在学习的是第18页，共69页
酰胺平面
肽键的键长介于单键和双键之间，具有部分双键的性质，不能自由旋转。
现在学习的是第19页，共69页
a螺旋
α－螺旋结构
最常见、含量最丰富的二级结构
现在学习的是第20页，共69页
现在学习的是第21页，共69页
现在学习的是第35页，共69页
蛋白质分子表面的环区域
现在学习的是第36页，共69页
蛋白质分子表面的环区域（LOOP）
大多数蛋白质分子是由a-螺旋和/或b回折构成的。这些a-螺旋和/
或b回折通常由不同长度和不规则形状的环区域相连接。 螺旋与回折的组合形成了分子的稳定的疏水内核，环区
域位于分子的表面。主链上的这环区域的C＝O和NH基团一
特点：不依赖于已知的结构模式，是一种普适的解 决方法，目前尚处于探索阶段
现在学习的是第42页，共69页
同源和比较模建方法
根据蛋白质结构的相似性，以已知蛋白质的结构为 模板构建未知蛋白质的三维结构
– 在进化过程中蛋白质三维结构的保守性远大 于序列的保守性，当两个蛋白质的序列同源 性/相似性高于35%时，一般情况下它们的三 维结构基本相同
针对局域能量极小，不是整个系统 能够计算含有大量原子的体系 简单有效，目前应用得最广泛
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
➢ 分子模拟——更广的一个定义，包括上述全部
精品课件
9
分子模拟历史及意义
分子模拟的思想溯源于上世纪初量子力学的兴起
1925年 1926年 1927年
1930年 1933年
1946年
……
Heisenberg发表了第一篇量子力学的文章， Schrodinger发表了著名的波动方程 采用量子力学的方法计算了氢分子的轨道，
Molecular modeling is
a technique for deriving, representing , manipulating the structures and reactions of molecules, and those properties that are dependent on these three dimensional structures
简单的说：计算机图形学和计算理论的结合
特点：更抽象的模型有关
精品课件
4
常见的分子模型
CPK
线
框 模 型
模 型
球
棍
棍
状
模
模
型
型
精品课件
5
常见的生物分子模型
碳飘片状精骨带品管课架模件状模型模型12型
6
应用举例——分子动力学
精品课件
7
分子模拟计算理论
根据基本原理的不同，分子模拟主要有量子力学 模型和分子力学模型两类
精品课件
13
ห้องสมุดไป่ตู้
分子模拟的常用软件
1，常用的平台软件
1，SYBYL 2，InsightII 3，HyperChem 4，PCmodel
2， 常用的分子动力学软件
1，AMBER（Assisted Model Building with Energy
Refinement )
2，GROMACS，免费
/pdb/Welcome.do
精品课件
15
分子模拟相关杂志
➢ J. Compt. Chem. ➢ J. Molecular Modeling ➢ J. Comp. Aid. Mol. Des. ➢ J. Chem. Inf. & Comp. Sci. ➢ J. Molecular Graphics and Modeling ➢ Network Science ➢ Computers and Chemsitry ➢ 计算化学邮件列表 (/chemistry/)
量子化学由此诞生； Andrews提出分子力学的基本思想； Bernal, J. D.和Flowler, R. H提出了水的
原子模型； Frank Westheiner完成了第一次分子力学计
算；
精品课件
10
分子模拟历史及意义
1974年 Allinger等人发布了分子力场MM
1981年 Peter Kollman发表了AMBER力场；
精品课件
11
分子模拟的主要方法
➢ 从头计算 ➢ 半经验计算 ➢ 固体计算 ➢ 分子力学 ➢ 分子动力学 ➢ 统计力学
精品课件
12
分子模拟的主要步骤
➢ 确定描述体系所用的模型 1. 量子力学 2. 分子力学
➢ 计算 1.能量优化， 2.分子动力学 3.蒙特卡罗模拟， 4.构象搜索等
➢ 分析 1，计算体系性质 2，检验结果
Model的特点：只显示分子的三维结构
常见的model如下面所示,
“stick”模型（Dreiding发明）
“CPK” （ “space filling” ）模型 （ Corey,Pauling和Koltun），
“Ribbon飘带模型 精品课件
3
The definition of molecular modelling
精品课件
14
分子模拟的常用软件
5，常用的分子对接软件
1，flex，SYBYL的一个模块
2，affinity，InsightII的一个模块
3，dock，免费 /
4，autodock，免费 /mb/olson/doc/autodock/
考虑电子运动 状态
考虑核运动状 态（电子运动 作近似假设）
量子力学模型
用波函数表示
分子力学模型
用质点运动轨 迹表示
解薛定谔方程
解力学方程 （经典牛顿力 学方程）
精品课件
8
几个容易混淆的术语
➢ 理论化学——利用数学模型来描述化学，常常 认为与量子化学同义
➢ 计算化学——包含量子力学和分子力学、极小 化、构象分析和基于计算机的其他方法，或者 说数学模型的计算机化
6，常用的分子显示软件
1，Viewlite，免费
2，VMD，免费 /Research/vmd/
3，Pymol，免费 /
7，常用的数据库
1，MDL，
2，ACD，可用晶体结构数据库
3，CSDC，剑桥晶体结构数据库
4，PDB库，美国布鲁克海文国家重点实验室蛋白质数据库，免 费
分子模拟原理及应用
molecular modelling
principles and appilications
计明娟 2006.9
精品课件
1
第一章 序言
化学就其物质的根本组成来说，是由电子及原子核 构成的。
1964年，加州大学的Walter Kohn教授和Hohenberg 博士从薛定谔方程出发，严格证明了一个重要性 质： 分子结构决定了分子性质
Ps：上述二人共同获得1998年诺贝尔化学奖
因此分子模拟就是要用理论方法（模拟分子结构） 去实现以往用实验才能证明的东西（性质）。
精品课件
2
The definition of model
The definition of model in Oxford English Dictionary：a simplified or idealized description of a system or process, often in mathematical terms, devised to facilitate calculations and predictions.
1983年 Martin Karplus发布了分子力学与分子动力学模拟程序CHARMM； …..
分子模拟曾经限制在部分科学家手中，现在随着理论和计算机 技术的发展，任何有兴趣的人都可以应用。而且应用范围越来越广， 尤其在生命科学，高分子，材料科学等多个领域取得了重要成果。
21世纪也正好是这些领域大展身手的世界，分子模拟的作用也 将更加重要。