Methodology

In Cartesian coordinate molecular dynamics calculations, the degrees of freedom are uncoupled, and Newton's equations of motion can be solved independently for each degree of freedom, :

Here, the mass of particle , , is used for that particle's three degrees of freedom, is the acceleration, and is the force. In internal coordinates, similar equations of motion can be written for systems with a tree-topology:

Here, is the -length vector of internal coordinates, is the -length vector of nonlinear forces and is the -length vector of generalize forces, or torques in the case of dihedral-angle degrees of freedom. In internal coordinates, the degrees of freedom are not uncoupled. The mass matrix, , has off-diagonal elements and has a nonlinear dependence on , so the equations cannot be solved independently for each degree of freedom, as is done for the Cartesian degrees of freedom in Equation (). The accelerations can be determined by solving the matrix equation in Equation (), as was done for the Lagrange equations of motion by Mazur et al.[35]. The computational cost of these matrix equations is proportional to , so they are prohibitive for large molecules. Recently, however, Jain et al.[33], have developed a recursive algorithm for solving the equations of motion which does not require explicit calculation of the mass matrix, , or direct solution of the matrix equation . The method uses spatial operator algebra in a recursive approach which is computationally proportional to . This linear dependence on opens up a new class of molecular systems to study by internal-coordinate molecular dynamics.

The details of the recursive spatial operator equations for solving the equations of motion can be found in Reference [33]. The methodology has been developed for general multibody systems configured as serial chains, topological trees, or closed-loop systems[32]. The implementation for molecular systems reported here has not yet been extended to closed-loop topologies, so only the serial chain and tree topology will be discussed here. The NEIMO method uses the concepts of ``clusters'' and ``hinges'' to describe a mechanical system. A cluster is a body which moves as a unit; in a molecule, this could be a single atom, a multiple-atom group such as a methylene group, a phenyl ring, or even an entire domain of a protein. A ``hinge'' describes the relative orientation between two connected clusters; in a molecular system, the hinges correspond to the bonds connecting adjacent clusters. There are six possible degrees of freedom for each hinge, including such modes as bond stretching, ball and socket rotation, and torsional rotational. Our interest is primarily in the torsional degrees of freedom, so each hinge in our molecular implementation is limited to a single torsional degree of freedom. The one exception per molecule is the ``base'' cluster, which is connected to the reference coordinate system by a hinge with the full six degrees of freedom, thereby enabling each molecule to be oriented correctly in space.

The relationship between adjacent clusters is described in terms of ``parents'' and ``children.'' The base cluster can have one or more child clusters; it is the parent cluster of each of these children. Each of these clusters, in turn, may have one or more children, branching outward from the base. In a topological tree, outward branching can continue with each cluster having zero, one, or more children, but each child having only one parent cluster. Clusters at the far extent of each branch are termed ``tips'' and have no children. In a serial chain, such as a linear polymer, there is a single tip cluster; between the base and tip clusters, each cluster has a unique parent and unique child. In a protein system, most of the C atoms are branch points with two children, and the outermost cluster of every sidechain (other than alanine and glycine, which have no sidechain clusters) is a tip. These concepts are visualized in Figure , where the pentapeptide Met-enkephalin is shown with the hinges numbered. This numbering system is not equivalent to the numbering used in Reference [33], but is used here to facilitate analysis of the dihedral angles. Hinge 0, which connects the base cluster to the reference frame, is not shown. The clusters can be assigned the same number as the hinge connecting it to its parent cluster and are listed that way in Table . In addition, the torsional degree of freedom for each hinge, other than hinge 0, can be defined to correspond to a specific dihedral angle. These dihedrals are also listed in Table , using the standardize nomenclature for protein dihedrals.

The NEIMO method uses spatial operator algebra to formulate the equations of motion in terms of relationships between parent and child clusters, and to solve the equations during several recursions from the base cluster to the tips, and the tips to the base. These Newton-Euler recursive equations of motion are described in detail in Reference [33], but can be summarized as:

1. A base to tips recursion, during which each cluster's spatial velocity is determined from the torsional motion of its associated hinge as well as the motion of its parent cluster.

2. A tips to base recursion, during which the effective forces acting on each cluster are derived from the explicit Cartesian forces, hinge torques, and Coriolis and related forces acting on the cluster, as well as the forces imparted to it from its children. Several other quantities related to the inverse of the mass matrix are also calculated during this recursion.

3. A final base to tips recursion, during which the acceleration of each cluster and its related hinge acceleration, are determined from the cluster's acceleration within its own frame of reference plus acceleration due to the acceleration of its parent cluster.

Currently, we use a ``leapfrog'' Verlet algorithm[29] to integrate the equations of motion, using the accelerations determined by the recursive spatial operator equations. The Verlet algorithm is a very successful and popular method in Cartesian-space dynamics, but may be less suitable for NEIMO dynamics because it calculates accelerations and velocities at alternating half-timesteps. In NEIMO dynamics, accelerations are not independent of velocities, so the half-timestep separation of accelerations and velocities must be modified. This can be done by iteratively solving for the velocities at integer timesteps. Fortunately, this iteration is very fast in practice. A major advantage of the Verlet algorithm is that it requires only a single calculation of the forces at each timestep. In simulations of large systems, the force calculation consumes the vast majority of computational time, so methods which require only a single force calculation are preferable to methods which require two or more force calculations per timestep, such as the Gear predictor-corrector algorithm[36]. Currently, other integration schemes are being investigated for use in NEIMO dynamics, but all results presented here use the leapfrog Verlet algorithm. Details of this integration method are given in Appendix .

The NEIMO calculations presented here were performed using a version of the program written to work with the BIOGRAF/POLYGRAF program from Molecular Simulations, Inc[37]. All calculations were performed on Iris PowerSeries and Iris Indigo workstations from Silicon Graphics, Inc.