After step (8b).a force evaluation is carried out,to give fi(t+ot)for step (8c).This scheme advances the coordinates and momenta over a timestep 6t.A piece of pseudo-code illustrates how this works: do step =1,nstep pp 0.5*dt*f r =r dt*p/m f force(r) p =p +0.5*dt*f enddo As we shall see shortly there is an interesting theoretical derivation of this version of the algorithm.Important features of the Verlet algorithm are:(a)it is exactly time reversible; (b)it is symplectic (to be discussed shortly);(c)it is low order in time,hence permitting long timesteps:(d)it requires just one(expensive)force evaluation per step:(e)it is easy to program. 3.2 Constraints It is quite common practice in classical computer simulations not to attempt to represent intramolecular bonds by terms in the potential energy function,because these bonds have very high vibration frequencies (and arguably should be treated in a quantum mechanical way rather than in the classical approximation).Instead,the bonds are treated as being constrained to have fixed length.In classical mechanics,constraints are introduced through the Lagrangian5 or Hamiltonian46 formalisms.Given an algebraic relation between two atomic coordinates,for example a fixed bond length b between atoms 1 and 2,one may write a constraint equation,plus an equation for the time derivative of the constraint X(r1,T2)=(r1-T2)·(r1-T2)-b2=0 (9a) X(T1,T2)=2(v1-v2)·(r1-T2)=0. (9b) In the Lagrangian formulation,the constraint forces acting on the atoms will enter thus: miri=fi+Agi where A is the undetermined multiplier and 0y=-2(r1-r2) g1=-r1 _0X=2(r1-r2) 92=-0r2 It is easy to derive an exact expression for the multiplier A from the above equations; if several constraints are imposed,a system of equations(one per constraint)is obtained. However,this exact solution is not what we want:in practice,since the equations of motion are only solved approximately.in discrete time steps,the constraints will be increasingly violated as the simulation proceeds.The breakthrough in this area came with the pro- posal to determine the constraint forces in such a way that the constraints are satisfied exactly at the end of each time step47-49.For the original Verlet algorithm,this scheme is called SHAKE.The appropriate version of this scheme for the velocity Verlet algorithm is called RATTLES0.Formally,we wish to solve the following scheme,in which we combine 8After step (8b), a force evaluation is carried out, to give f i (t + δt) for step (8c). This scheme advances the coordinates and momenta over a timestep δt. A piece of pseudo-code illustrates how this works: do step = 1, nstep p = p + 0.5*dt*f r = r + dt*p/m f = force(r) p = p + 0.5*dt*f enddo As we shall see shortly there is an interesting theoretical derivation of this version of the algorithm. Important features of the Verlet algorithm are: (a) it is exactly time reversible; (b) it is symplectic (to be discussed shortly); (c) it is low order in time, hence permitting long timesteps; (d) it requires just one (expensive) force evaluation per step; (e) it is easy to program. 3.2 Constraints It is quite common practice in classical computer simulations not to attempt to represent intramolecular bonds by terms in the potential energy function, because these bonds have very high vibration frequencies (and arguably should be treated in a quantum mechanical way rather than in the classical approximation). Instead, the bonds are treated as being constrained to have fixed length. In classical mechanics, constraints are introduced through the Lagrangian45 or Hamiltonian46 formalisms. Given an algebraic relation between two atomic coordinates, for example a fixed bond length b between atoms 1 and 2, one may write a constraint equation, plus an equation for the time derivative of the constraint χ(r1, r2) = (r1 − r2) · (r1 − r2) − b 2 = 0 (9a) χ˙(r1, r2) = 2(v1 − v2) · (r1 − r2) = 0 . (9b) In the Lagrangian formulation, the constraint forces acting on the atoms will enter thus: mir¨i = fi + Λgi where Λ is the undetermined multiplier and g1 = − ∂χ ∂r1 = −2(r1 − r2) g2 = − ∂χ ∂r2 = 2(r1 − r2) It is easy to derive an exact expression for the multiplier Λ from the above equations; if several constraints are imposed, a system of equations (one per constraint) is obtained. However, this exact solution is not what we want: in practice, since the equations of motion are only solved approximately, in discrete time steps, the constraints will be increasingly violated as the simulation proceeds. The breakthrough in this area came with the pro￾posal to determine the constraint forces in such a way that the constraints are satisfied exactly at the end of each time step47–49 . For the original Verlet algorithm, this scheme is called SHAKE. The appropriate version of this scheme for the velocity Verlet algorithm is called RATTLE50 . Formally, we wish to solve the following scheme, in which we combine 8