E. Barth, K. Kuczera, B. Leimkuhler and R.D. Skeel

In molecular dynamics simulations, the fastest components of the potential field impose severe restrictions on the stability and hence the speed of computational methods. This difficulty can be avoided by replacing these fast components with algebraic constraints, with little effect on the underlying dynamics. In this paper, the resulting systems of mixed differential and algebraic equations are studied. Commonly used discretization schemes for constrained Hamiltonian systems are discussed. The form of the nonlinear equations is examined in detail and used to give convergence results for the traditional nonlinear solution technique SHAKE iteration and for a modification based on Successive OverRelaxation (SOR). A simple adaptive algorithm for finding the optimal relaxation parameter is presented. Alternative direct methods using sparse matrix techniques are discussed. Finally, numerical results are given for the new and existing techniques implemented in the molecular modeling software package CHARMM .