This paper is a review of the method of constraints. The method was devised to carry out Molecular Dynamics simulations of complex molecular systems with some internal degrees of freedom frozen, in terms of atomic Cartesian coordinates. The method has been subsequently generalized to treat all kinds of holonomic constraints and has been adapted to the more recent dynamical simulations in ensembles different from the microcanonical one. We start by deriving the statistical-mechanical formalism for these systems. We then proceed to derive the equations of motion. We conclude with a detailed discussion of the relevant MD algorithms. Some technical details on the FORTRAN code are presented in an appendix.