4/26/2005 CHEM 408 - Sp05L16 - 1Molecular Dynamics SimulationsMolecular Dynamics Simulations• compute thermally averaged properties by sampling a representative phase-space trajectory via numerical solution of the classical equations of motion• for a system of N interacting particles the coupled e.o.m. are:(Leach 7.1-7.5)iiiFmdtrdrr122=xiiiFmdtxd122=Ni →=1)(RFirirrrrΦ∇−=where the force on particle i is related to the 1st derivative of the intermolecular potential Φ by: all 3N particlecoordinates)(RxFixirΦ∂∂−=4/26/2005 CHEM 408 - Sp05L16 - 2• some simple examples: km1d harmonic oscillatorx-x0-1.0 -0.5 0.0 0.5 1.0v(x)/k, F(x)/(2k)-0.20.00.20.4r / σ1.0 1.5 2.0 2.5 3.0Φ(r)/ε, F(r)/(2ε/σ)-1.5-1.0-0.50.02021)()( xxkx −=Φ)()(0xxkxF−−=⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎟⎟⎠⎞⎜⎜⎝⎛−⎟⎟⎠⎞⎜⎜⎝⎛=Φ6124)(ijijijrrrσσεijrrij⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎟⎟⎠⎞⎜⎜⎝⎛−⎟⎟⎠⎞⎜⎜⎝⎛=7132ˆ24ijijijirrrFσσσεr2 interacting LJ particlesFΦΦF4/26/2005 CHEM 408 - Sp05L16 - 3• trajectory is propagated using finite difference methods, dt→δt• use Taylor expansion ...))((!31))((!21))(()()(32+−′′′+−′′+−′+= axafaxafaxafafxfon r(t) to find r(t+δt):...)(!31)(!21)()()(32++++=+ ttbttattvtrttrδδδδrrrrrdtrdtrtvrrr=′= )()(where22)()(dtrdtrtarrr=′′=33)()(dtrdtrtbrrr=′′′=44)()()(dtrdtrtcivrrr==velocityaccelerationetc.jerk4/26/2005 CHEM 408 - Sp05L16 - 4• consider expanding both r(t+δt) and r(t-δt):...)()()()()(361221++++=+ ttbttattvtrttrδδδδrrrrr...)()()()()(361221+−+−=− ttbttattvtrttrδδδδrrrrr(a)(b)(a+b))()()(2)()(42tOttatrttrttrδδδδ++=−++rrrr• these result is the original Verlet algorithm: 2)()()(2)( ttattrtrttrδδδrrrr+−−=+• the Verlet algorithm is correct to 4thorder in positions; velocities are implicit but can be approximated by: tttrttrtvδδδ2)()()(−−+≅rr• a(t) is determined by F{r(t)}; so calculation sequence is:• calc a(t)=F{r(t)}/m• calc r(t+δt) from a(t) and r(t-δt) • calc v(t) if desired from r(t+δt) and r(t-δt)4/26/2005 CHEM 408 - Sp05L16 - 5• variants to the original Verlet algorithm provide for more accurate velocities:Verlet Leap-Frog algorithm: ttattvttvtttvtrttrδδδδδδ)()()()()()(212121+−=+++=+rrrrvelocity Verlet algorithm: tttattvttvtttvtrttrttatvttvδδδδδδδδδ)()()()()()()()()(2121212121+++=+++=++=+rrrrrrr and v known at the same set of times• closely related to these Verlet methods is the Beeman algorithm• more complicated Gear predictor-corrector methods also usedboth r and v known with good accuracy4/26/2005 CHEM 408 - Sp05L16 - 6originalleap-frogvelocitySchematic Summary of Verlet VariantsFigure from Allen & Tildesley, Computer Simulation of Liquids (Oxford, 1987)4/26/2005 CHEM 408 - Sp05L16 - 7• choosing an integrator is partly a matter of taste; key point is to minimize # force evaluations (the costly part) and maintain accuracy • more complicated algorithms involve more computation (often moreF evaluations) but allow for larger δt for given accuracy• for a given integrator, choice of time step is importantAr - Ar collision “just right”small δtlargeδtFigures from Leach4/26/2005 CHEM 408 - Sp05L16 - 8Digression on Energies)()(),(3333 NNNNpKqpqH +Φ=(classical) Hamiltonianpotential energykinetic energydependence of system energy on phase space coordinate><+>Φ>=<=<KHEwhere <x> denotes a phase space average of x either over a dynamic trajectory (MD) or a static ensemble (MC)• the energy of a classical system can be described in terms of its phase space coordinates by:• the thermodynamic (internal) energy E is an average quantity:• in simple (N,V,E) MD simulations H is a constant of the motion whereas Φ and Kfluctuate in timeTime0 1020304050Energy-2-1012V(t)H(t)K(t)4/26/2005 CHEM 408 - Sp05L16 - 9“Coordinate” q3NEnergyΦ(t)K(t)q3N(t)Φ(q3N)EEnergyCoordinate q0Φ(q)K(t)Φ(t)q(t)E1d Harmonic OscillatorAr Liquid2 Schematic Examples4/26/2005 CHEM 408 - Sp05L16 - 10Figure from Leach2. Energy Conservation in MDLiquid Ar (N=256; T~100 K) .01K(t)E(t)Φ(t)~8Effect of Step Size δtEnergy FluctuationTime Step δtFigure from Allen & Tildesley, Computer Simulation of Liquids (Oxford, 1987)velocity VerletG4G6G54/26/2005 CHEM 408 - Sp05L16 - 11• calculating exact trajectories is impossible Liquid Ar N~108(?)positions ∆renergies ∆K/K2012|)()(|1|)(| trtrNtriNiirr∑=−=∆• compare reference system “0” to ones in which particles are displaced by 10-3, 10-6, 10-9σ• settle for accurate trajectories for times > molecular correlation times (ps)Figure from Allen & Tildesley, Computer Simulation of Liquids (Oxford, 1987)• energy conservation such that 22KE ∆<<∆4/26/2005 CHEM 408 - Sp05L16 - 12• Leach’s recommendations on step size:• efficient simulation is made difficult by the presence of motions having different characteristic time constants; one example encountered in many problems is the fact that H-atom vibrations are much faster than most other vibrations• two solutions are:multiple time step algorithms: • calculate forces at different δt depending on how quickly they change with time• r-RESPAconstraint dynamics: • freeze the motion of fastest degrees of freedom• SHAKE, RATTLE4/26/2005 CHEM 408 - Sp05L16 - 133. Initializing & Running Simulations• initial coordinates from:- prior simulation- lattice configuration- crystal structure (proteins)• velocities sampled from Maxwell-Boltzmann distriution at proper T:⎭⎬⎫⎩⎨⎧−⎟⎟⎠⎞⎜⎜⎝⎛=TkvmTkmvPBixiBiix22/121exp2)(π• equilibrate (evolve) system, perhaps changing K until properties have reached apparent equilibriumAr(l) near TP started from fcc lattice (N=108)Figure from Allen & Tildesley, Computer Simulation of Liquids (Oxford, 1987)rms displacementtranslational orderpressureKVPH∑=⋅=NiirkNk1)cos(1)(rrρ4/26/2005 CHEM 408 - Sp05L16 - 14• the thermodynamic pressure P is related to the “virial”∑===NiiiBmpNkNKT1231324. Calculating Some Simple Properties (N,V,E)• internal energy E is a constant • thermodynamic temperature T is proportional to the average kinetic energy <K> by:∑∑>=−=NijNiijijBFrVVTNkP131• the dielectric constant εis related to the fluctuations in the net sample dipole moment Σµby:)2()12(319421εεµπ++−=⎟⎠⎞⎜⎝⎛∑=ssNiiBDDTVkr(expression depends on electrostatic BCs; this equation applies to reaction field BCs with