LINCS: A Linear Constraint Solver for Molecular Simulations

1. LINCS: A Linear Constraint Solver for Molecular Simulations Berk Hess, Hemk Bekker, Herman J.C.Berendsen, Johannes G.E.M.Fraaije Journal of Computational Chemistry, 1997 Ankur Dhanik

2. Outline • Introduction to Molecular Dynamics • Problem description • Some solutions • LINCS • Results

3. Introduction to Molecular Dynamics • A classical molecular simulations method • Successive configurations of the system are generated by integrating Newton’s laws of motion • Integration algorithms should strive to reduce computation and conserve energy • Various integration algorithms • Standard Verlet algorithm • Leap-frog algorithm

4. Introduction to Molecular Dynamics • Standard Verlet Algorithm • r(t+Δt) = 2r(t) - r (t-Δt) + Δt2a(t) • v(t) = [r(t+Δt) - r(t-Δt)]/2Δt • Velocities are not directly generated, and are one time step behind • Leap-frog algorithm • r(t+Δt) = r(t) + Δtv(t+Δt/2) • v(t+Δt/2) = v(t-Δt/2) + Δta(t) • v(t) = [v(t+Δt/2) + v(t-Δt/2)]/2 • Velocities are half time step behind

5. Introduction to Molecular Dynamics • Choosing the time step • One order of magnitude smaller than the shortest motion (bond vibrations) • Severe restriction as these high frequency motions have minimal effect on the overall behavior of the system • Constrained dynamics

6. Introduction to Molecular Dynamics The different types of motion present in various systems together with suggested time steps

7. Fc g Introduction to Molecular Dynamics • In constrained dynamics bonds and angles are forced to adopt specific values throughout a simulation • Constraints are categorized as holonomic and non-holonomic. • Suppose a particle on the surface of a sphere • r2 – a2 = 0 Holonomic • r2 – a2 >= 0 Non-holonomic • In a constrained system • Particles are not independent • Magnitude of constraint forces are unknown

8. Introduction to Molecular Dynamics • Integration algorithms • Standard Verlet Algorithm • Leap-frog algorithm • Choosing time step • Constrained dynamics • Holonomic and non-holonomic constraints

9. Outline • Introduction to Molecular Dynamics • Problem description • Some solutions • LINCS • Results

10. Problem Description • Design an algorithm for solving constrained molecular dynamics, the constraints being holonomic • The algorithm should strive for following features: • Numerical stability • Energy conservation • Computational efficiency

11. Some solutions • Reset coupled constraints after an unconstrained update • Non-linear problem • SETTLE • Solves analytically • Very fast, but unsuited for large molecules • SHAKE • Iterative method • Sequentially all the bonds are set to the correct length • Simple and numerically stable • No solutions may be found when displacements are large, difficult to parallelize

12. Some solutions • EEM • The second derivatives of constraint equations are set to zero • All the constraints are dealt with simultaneously • Corrections are necessary to achieve accuracy and stability

13. LINCS • Does an unconstrained update • Sets the projection of the new bonds onto the old directions of the bonds to the prescribed lengths • Similar to EEM, with some practical differences • Implements • Efficient solver for the matrix equation • A velocity correction that prevents rotational lengthening • A length correction that improves accuracy and stability

14. Newton’s equation of motion Constraint equations Lagrangian formulation Gradient matrix of constraint equations, B LINCS

15. 1st and 2nd derivatives of constraints are zero Constraint forces , where LINCS Newton’s equation of motion (I-TB) is projection matrix which sets the constrained coordinates to zero. T transforms motions in the constrained coordinates into motions in cartesian coordinates

16. Constrained leap-frog algorithm =0 is the projection matrix that sets the constrained coordinates to zero =0 Since If we set We can use LINCS

17. LINCS • Numerical errors can accumulate which leads to drifts Velocity correction Position correction • Drawback: the projection of the new bonds onto the old directions rather than the new bonds are set to prescribed lengths

18. LINCS • Correction for rotational lengthening • To correct rotation of bond i, the projection of the bond on the old direction is set to • The corrected positions are

19. Half of the CPU time goes to invert with diagonal element & A is symmetric, sparse and has zeros on the diagonal LINCS Constrained new position is given by

20. LINCS • The first power of An gives the coupling effects of neighboring bonds • The second power gives the coupling effect over a distance of two bonds • The inversion through a series expansion makes parallelization easy • In one timestep, the bonds influence each other when they are separated by fewer bonds than highest order in expansion

21. LINCS • Parallelization • Consider a linear-bond constrained molecule of 100 atoms to be simulated on a dual processor computer • Uses rotation correction and an expansion to the second power of An • Because order of expansion is two, bonds influence each other over a distance of 6 • Update of position and call of LINCS algorithm must be done for atom 1-56 and 44-100 on processors 1 and 2 respectively • 1-50 update from processor 1 and 50-100 from processor 2

22. Results

23. Problem Description • Design an algorithm for solving constrained molecular dynamics, the constraints being holonomic • The algorithm should strive for following features: • Numerical stability • Energy conservation • Computational efficiency Results • Solves constrained molecular dynamics • Numerically stable • Conserves energy • Three to four times faster than SHAKE • Can be easily parallelized