A Second-Order Strong Method for the Langevin Equations with Holonomic Constraints

A numerical predictor-corrector method is presented for the integration of the Langevin equations with holonomic constraints, specifically Kramers' ball-rod model for long-chain polymers. In the predictor, constraint forces are expressed as Lagrange multipliers and evaluated as implicit functions of particle coordinates. The resulting expressions, with and without the Duhamel form, are developed as a strong second-order Itô-Taylor method by the expansion technique of Wagner and Platen. The corrector is a projection that enforces constraints to machine precision. The numerical evaluation of the stochastic integrals, as well as the “coarsening” relations used to recursively construct integral sets when the time step is doubled, are described in detail. We present numerical simulations that demonstrate the strong and weak orders of convergence. With parameters approximating $\lambda$-phage DNA, the Duhamel formulation is stable with time steps exceeding the relaxation time. As a validation study we evolve initially linear polymers and observe that they evolve to configurations with mean square end-to-end distance predicted by analytical theory. We also simulate a velocity correlation function in agreement with classical theory.