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 and determining the relative free energy of binding of different drug molecules to the same receptor. 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, improved efficiency in calculating interatomic forces, and development of techniques for allowing larger timesteps in molecular dynamics simulations.
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, 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 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, 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.
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. 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. 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 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, which has three chains of nearly 300 residues each.