Molecular dynamics simulations have become a very important tool in computational chemistry and biology. Once primarily useful for understanding the dynamic properties of molecular systems, molecular dynamics has become invaluable for such diverse tasks as building protein models from crystallographic data[22] and determining the relative free energy of binding of different drug molecules to the same receptor[23]. As these calculations become more accurate and applicable to a wider variety of problems, researchers seek to apply them to larger and more complex systems as well as to use them to study processes occurring over increasingly longer timescales. This requires continual improvements in both computer hardware and software performance. The former challenge is being addressed by such developments as vector processing supercomputers, RISC workstations and now massively parallel supercomputers. The software problem is being tackled on several overlapping fronts: optimization of programs for new computer architectures[24], improved efficiency in calculating interatomic forces[26][25], and development of techniques for allowing larger timesteps in molecular dynamics simulations[33][27].
Molecular dynamics simulations typically involve numerical integration
of Newton's equations of motion. Timesteps for the integration must be
small enough for the fastest modes to be handled accurately.
Typically, a timestep of 1 femtosecond ( s) is necessary to
handle bond stretches involving hydrogen atoms. Although success has
been achieved lately by separating short and long-range forces and
using different timesteps for the different forces[27],
the most popular approach for increasing timesteps is to fix the
fastest degrees of freedom (bond stretches and angles) and to solve
the equations of motion for the slower (dihedral) degrees of freedom.
Such an approach is especially justified for studies of large
biological molecules, where bond lengths and angles vary little from
one structure to another and nearly all important conformational
transitions are due to dihedral angle motions.
The SHAKE algorithm[28] has become the standard approach
for doing molecular dynamics with fixed bond lengths. It can also be
used to hold angles fixed, but this is less common. SHAKE is a
modification of the Verlet algorithm for integrating the equations of
motion for the Cartesian coordinates degrees of freedom in an
-particle system. Particle velocities are first calculated for the
unconstrained system, then modified to meet each constraint. An
iterative process is required to meet all the constraints
concurrently. The SHAKE algorithm works well for timesteps up to
5 fs[31][30], thereby enabling a
five-fold speedup in computational time as long as the process of
iteratively solving the constraint equations does not consume too much
time[30].
An alternative to the SHAKE methodology of solving Cartesian
coordinate dynamics plus constraints, is to solve the equations of
motion directly for the internal degrees of freedom. Solutions to
these equations automatically fulfill the desired bond length and/or
angle constraints, so their efficiency is not limited by a secondary
constraint-solving step. Indeed, Mazur et al.[35]
were able to simulate accurately a small polypeptide, (Ala), with
timesteps as large as 13 fs. This is a significant improvement over
the SHAKE algorithm. Unfortunately, their method requires the direct
solution of matrix equations; the computational time for solving these
equations is usually proportional to
, where
is the number of degrees of freedom. This can be prohibitive for
large systems.
Recently, Jain et al.[33][32] have developed an
alternative method for solving the equations of motion for internal
coordinates. This new method, the Newton-Euler Inverse Mass Operator
(NEIMO) method, does not require direct manipulation of matrices and
is proportional to , rather than
. This method
was developed for spacecraft dynamics, but in a separate report, Jain
et al[33] described how the method could be applied
to molecular dynamics. This report presents the first implementation of
the NEIMO method for molecular systems. We have studied the dynamics of
polypeptide systems and have found that we are able to calculate the
dynamics of some systems accurately with timesteps as large as 20 fs.
Because the NEIMO method is computationally proportional to
,
it can be applied to very large systems. The method is shown to be
rigorously proportional to
for systems as large as the tomato
bushy stunt virus crystal structure[45], which has three
chains of nearly 300 residues each.