SHAKE is a widely used algorithm to impose general holonomic constraints during molecular simulations. By imposing constraints on stiff degrees of freedom that require integration with small time steps (without the constraints) we are able to calculate trajectories with time steps larger by approximately a factor of two. The larger time step makes it possible to run longer simulations. Another approach to extend the scope of Molecular Dynamics is parallelization. Parallelization speeds up the calculation of the forces between the atoms and makes it possible to compute longer trajectories with better statistics for thermodynamic and kinetic averages. A combination of SHAKE and parallelism is therefore highly desired. Unfortunately, the most widely used SHAKE algorithm (of bond relaxation) is inappropriate for parallelization and alternatives are needed. The alternatives must minimize communication, lead to good load balancing, and offer significantly better performance than the bond relaxation approach. The algorithm should also scale with the number of processors. We describe the theory behind different implementations of constrained dynamics on parallel systems, and their implementation on common architectures.