P-SHAKE: A quadratically convergent SHAKE in O(n2)O(n2)

An algorithm for solving arbitrary linear constraints in molecular dynamics simulations of rigid and semi-rigid molecules is presented. The algorithm – P-SHAKE – is a modified version of the SHAKE [J.-P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n  -alkanes, J. Comput. Phys. 23 (1977) 327–341.] algorithm with a preconditioner applied which effectively de-couples the constraint equations. It achieves quadratic convergence, as does M-SHAKE [V. Kräutler, W.F. van Gunsteren, P.H. Hünenberger, A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J. Comput. Chem. 22 (5) (2001) 501–508.], yet at a cost of only O(n2)O(n2) operations per iteration, as opposed to O(n3)O(n3) per iteration for M-SHAKE. The algorithm is applied to simulations of rigid water, DMSO, chlorophorm and non-rigid ethane and cyclohexane and is shown to be faster than M-SHAKE by up to a factor of three for relatively small error tolerances.