Three new algorithms are presented for incorporating nonholonomic constraints into molecular dynamics (MD) simulations, along with any additional holonomic constraints. The advantages of these algorithms over the commonly used Gaussian approach are discussed. Of the three algorithms presented, the optimal one can efficiently ensure satisfaction of large numbers of nonholonomic and holonomic constraints at every MD time step, without introducing additional numerical errors in the coordinate or velocity trajectories. Numerical results from MD simulations of Lennard-Jones particles, rigid water molecules, and partially rigid methane molecules are given, illustrating the advantages of this algorithm. In addition, this algorithm is suggested as a more advantageous alternative to velocity scaling, for maintaining fixed temperature during equilibration of constant energy MD simulations.