Molecular Dynamics (MD) Simulation

As already noted, MD simulation generally begins where experimental structure determination leaves off, if not during the structure refinement itself. It is generally not used to predict structure from sequence or to model the protein folding pathway. MD simulation can fold extended sequences to `global' potential energy minima for very small systems (peptides of length ten, or so, in vacuum), but it is most commonly used to simulate the dynamics of known structures.

An initial velocity is assigned to each atom, and Newton's laws are applied at the atomic level to propagate the system's motion through time (see `Classical and Quantum Mechanics - in a Nutshell' above). Thus, dynamical properties such as time correlation functions and transport coefficients (e.g., diffusion constants, bulk viscosities) can be calculated from a sufficiently long MD trajectory.

Once again, Newton's second law is: $\vec{F}_i = m_i \vec{a}_i$, where $\vec{F}_i$ is the sum of all forces acting on atom i that results in its acceleration $\vec{a}_i$. The acceleration is the second derivative of the position with respect to time: $\vec{a}_i = d \vec{v}_i/dt = d^2 \vec{r}_i/dt^2$. In words, it is the rate of change of the velocity $\vec{v}_i$, which in turn, is the rate of change of the position $\vec{r}_i$.

The `Leap Frog' algorithm is one method commonly used to numerically integrate Newton's second law. We obtain all atomic positions $\vec{r}_i$ at all times $t_n$ and all atomic velocities $\vec{v}_i$ at intermediate times $t_{n+1/2}$. This method gets its name from the way in which positions and velocities are calculated in an alternating sequence, `leaping' past each other in time:

$$\vec{v}_i(t_{n+1/2}) = \vec{v}_i(t_{n-1/2}) + \frac{\vec{F}_i(t_n)}{m_i} \Delta t,$$

$$\vec{r}_i(t_{n+1}) = \vec{r}_i(t_n) + \vec{v}_i(t_{n+1/2}) \Delta t.$$

Initial velocities are assigned so as to reflect equilibrium at the desired temperature T (a Maxwellian distribution), without introducing a net translation or rotation of the system.

The energy of an isolated system (as opposed to, for example, one in contact with a thermal bath) is conserved in nature, but it may not be in simulations. Energy conservation can be violated in simulations because of an insufficiently short integration time step $\Delta t$, an inadequate cutoff method applied to long-range (electrostatic and Lennard-Jones) forces, or even bugs in the program. Of course, energy conservation alone is not sufficient to ensure a realistic simulation. The realism of the dynamics trajectory depends on the empirical potential energy function $V(\vec{R})$, the treatment of long-range forces, the value of $\Delta t$, etc.