In this paper, we develop a recursive algorithm for solving the dynamical equations of motion for molecular systems. We make use of internal variable models which have been shown to reduce the computation times of molecular dynamics simulations by an order of magnitude when compared with Cartesian models. The O(N) algorithm described in this paper for solving the equations of motion provides additional significant improvements in computational speed. We make extensive use of the spatial operator methods which have been developed recently for the analysis and simulation of the dynamics of multibody systems. The spatial operators are used to derive the equations of motion and obtain an operator expression for the system mass matrix. An alternative square factorization of the mass matrix leads to a closed form expression for its inverse. From this follows the recursive algorithm for computing the generalized accelerations. The computational cost of this algorithm grows only linearly with the number of degrees of freedom. This is in contrast to conventional constrained dynamics algorithms whose cost is a cubic function of the number of degrees of freedom. For the case of a polypeptide molecule with 400 residues, the O(N) algorithm provides computational speedup by a factor of 450 over the conventional O(N 3) algorithm. We also describe a simplified method for computing and handling the potential function gradients within the dynamics computations.