MOLECULAR MECHANICS AND DYNAMICS

04/13/16

G.  Molecular Dynamics

This section comes directly from the NIH tutorial:

"In the broadest sense, MD is concerned with molecular motion . Motion is inherent to all chemical processes. Simple vibrations, like bond stretching, and angle bending, give rise to IR spectra. Chemical reactions, hormone-receptor binding, and other complex processes are associated with many kinds of intra- and intermolecular motions.

The driving force for chemical processes is described by thermodynamics. The mechanism by which chemical processes occurs is described by kinetics. Thermodynamics dictates the energetic relationships between different chemical states, whereas the sequence or rate of events that occur as molecules transform between their various possible states is described by kinetics.

Conformational transitions and local vibrations are the usual subjects of molecular dynamics studies. MD alters the intermolecular degrees of freedom in a step-wise fashion, analogous to energy minimization. The individual steps in energy minimization are merely directed at establishing a down-hill direction to a minimum. The steps in MD, on the other hand, meaningfully represent the changes in atomic position, ri, over time (i.e. velocity).

Newton's equation (Fi = mi ai) is used in the MD formalism to simulate atomic motion. The rate and direction of motion (velocity) are governed by the forces that the atoms of the system exert on each other as described by Newton's equation. In practice, the atoms are assigned initial velocities that conform to the total kinetic energy of the system, which in turn, is dictated by the desired simulation temperature. This is carried out by slowly heating the system (initially at absolute zero) and then allowing the energy to equilibrate among the constituent atoms. The basic ingredients of MD are the calculation of the force on each atom, and form that information, the position of each atom through a s specified period of time (typically on the order of picoseconds = 10-12 seconds).

The force on an atom can be calculated from the change in energy between its current position and its position a small distance away. This can be recognized as the derivative of the energy with respect to the change in the atom's position: -dE/dri = Fi

Energies can be calculated using either MM or quantum mechanics methods. MM energies are limited to applications that do not involve drastic changes in electronic structure such as bond making/breaking. Quantum mechanical energies can be used to study dynamic processes involving chemical changes. The latter technique is extremely novel and of limited availability.

Knowledge of the atomic forces and masses can then be used to solve for the positions of each atom along a series of extremely small time steps (on the order of femtoseconds). The resulting series of snapshots of structural changes over time is called a trajectory. The use of this method to compute trajectories can be more easily seen when Newton's equation is expressed in the following form

-dE/dr_i = m_ia= md²r_i/dt²

In practice, trajectories are not directly obtained from Newton's equation due to lack of an analytical solution. First the atomic accelerations are computed from the forces and masses. The velocities are next calculated from the accelerations based on the following relationship:
ai = dvi/dt. Lastly, the positions are calculated from the velocities: vi = dr_i/dt. A trajectory between two states can be subdivided into a series of sub-states separated by a small time step, delta t (e.g. 1 fs).

The initial atomic positions at time t are used to predict the atomic positions at time t = delta t. The positions at t = delta t are used to predict the positions at t = 2delta t, and so on.

The leapfrog method is a common numerical approach to calculating trajectories based on Newton's equation. The method derives its name from the fact that the velocity and position information successively alternate at ½ time step intervals.

MD has no defined point of termination other than the amount of time that can be practically covered. Unfortunately, the current ps order of magnitude limit is often not long enough to follow many kinds of state to state transformations, such as large conformational transitions in proteins."

A quick-time movie of a molecular dynamics simulations of 86 water molecules, run for 1.6 ps, created using Quanta/Charmm!

QuickTime:  PhiPsi movement

MD simlulations can be used to obtained theoretical values for ΔG and Keq values for conformational changes, binding of small  ligands and changes in protonation states for side chains.   This process is based on the idea that the conformations sampled in in silico MD simulations reflect those found in vitro (i.e. they are part of the thermodynamically expected and available conformations available to the molecules during its normal conformational shifts).  This is called the Ergodic Hypothesis.  Given the short time spans for MD simulations (limited by computer power) this hypothesis can't apply to the dynamic results unless the sample conformations are close in energy without a large activation energy barrier between them.     If it is then the following equation could apply:

ΔG0 = -RTlnK_eq

ΔG0 = -RTlnP₂/P₁ = -RTlnf₂/f₁ where

Pn is the probability  of being in a given state and fn is the fraction in a given state.

You will note that none of the potential energy function use quantum mechanical parameters.  This is due, in part, to the complexity of the systems studied.   This is beginning to change as more effort is devoted to understanding quantum mechanical aspects of complex bonded systems.  New advances in even simple systems can illustrate this point.  Take for example ethane.  The conformational analysis of this simple molecule is discussed in all organic chemistry books.  Plots of energy vs dihedral angle (viewing the molecule down the C-C bond and measuring the angle between C-H bonds on adjacent C atoms) oscillates every 120^o.  The E is at a maximum when the dihedral angle is 0, 120, 240 and 360^o, each representing the eclipsed conformation.  It reaches minimums at the staggered (gauche) conformations at 60, 180, and 270^o.  Why is the eclipsed form higher in energy than the staggered form.  All organic books would state that there is greater steric repulsion (of the electron clouds) in the eclipsed forms, which raises their energy compared to the staggered forms?  However,  Pophristic shows that to be incorrect.  For the correct answer, you must turn to quantum mechanics and the phenomena of hyperconjugation.  The staggered conformation is energetically favored not since it is less sterically restricted, but since its is a lower energy form due to resonance like stabilization of the σ CH molecular orbitals.  There is greater correct phase overlap of σ CH and s* CH molecular orbitals on the adjacent Cs when they are in the staggered conformation than in the eclipsed form.

Navigation

Return to Molecular Mechanics and Dynamics Contents

Return to Biochemistry Online Table of Contents

Archived version:  Molecular Mechanics and Dynamics