Topics in Molecular Modeling: II. Numerical Methods for the Time-integration

1. Topics in Molecular Modeling: II. Numerical Methods for the Time-integration Xiantao Li (xiantao.li@gmail.com) Department of Mathematics, Pennsylvania State University

2. Outline Hamiltonian system (NVE ensemble) Operator-splitting methods Energy conservation properties Stochastic models Strong and weak convergence Operator-splitting methods

3. Hamiltonian systems Hamilton’s principle Examples: mass spring model, many-body problems, pendulum model, Lotka-Volterra model etc. The mapping is symplectic A numerical method is symplectic if it defines a symplectic mapping

4. Many-Particle Systems Coordinate and Momentum The equations: Linear momentum: Angular momentum: Total energy: Conservation:

5. First Integrals (Invariants) In general, for an ODE system, is a first integral (or invariants) if Most numerical methods preserve linear invariant Only some methods preserve quadratic invariants

6. Symmetry If the function is symmetric wrt. The solutions of the corresponding ODEs will be ρ- reversible . For example, A numerical method is ρ- reversible if the solution satisfies the same properties. Explicit RK methods are not symmetric.

7. Symplectic Structure Let . If y, then the Hamiltonian ODEs can be written as A linear mapping A is symplectic if A nonlinear mapping g(x) is symplectic if is symplectic matrix. For Hamiltonian system, the Jacobian matrix is symplectic. Therefore, the mapping is symplectic

8. Symmetric and symplectic methods • Symmetric methods • Order of accuracy is always even • Very convenient for an extrapolation procedure • Symplectic methods • Very good energy conservation properties • Promising accuracy over long time integration • Provide good statistics. • Explicit RK methods are not symplectic.

9. A special class of symplectic methods Nonlinear dynamical system: . Wronskian: . Variational equation: Fundamental solution: In particular, If then Therefore,

10. Semi-group notation For any quantity of interest: Time derivative: L is a differential operator w.r.t. Time evolution: This is the foundation for operator-splitting methods.

11. Newtonian dynamics Equations of motion: Operator: Operator splitting: corresponds to corresponds to

12. Operator splitting Both equations are solvable Do we have No, unless

13. BCH (Baker-Campbell-Hausdorff) formulas Commutator: BCH formula I: First-order splitting: Implementation: Symplectic Euler:

14. Symmetric Splitting BCH formula II: Symmetric splitting: Implementation:

15. Properties of the second splitting method Second order accuracy in time Symplectic Symmetric Energy conservation: Conserve the linear and angular momentum

16. General splitting methods Introduce multiple fractional steps The coefficients are determined by maximizing the order of accuracy. Formulas are available up to eighth order.

17. Extended Hamiltonian systems • In practice, the Hamiltonian system is often connected a heat bath or pressure bath, modeled by additional equations. • Nose-Hoover thermostat • Symmetric splitting:

18. Problem with multiple scales • Hamiltonian system driven by two forces • short-range fast force • long-range slow force; More difficult to compute. • Multiple time-stepping: evaluate on the fast time scale with step size , and evaluate g on the slow time scale with step size

19. Multiple time-stepping Let be the evolution operator for the fast dynamics: Let be the evolution operator for the slow dynamics Overall algorithm: . It is symplectic. Operations: Kick + Evolve + Kick. Local accuracy:

20. More accurate splitting To eliminate the leading term in the error, we constructed the following splitting method This method also have better energy conservation property.

21. Example: dynamics of butane molecule

22. Langevin dynamics Langevin stochastic differential equation: : position, velocity : non-linear conservative potential, : friction matrix constant and diffusion matrix constant, related by the fluctuation dissipation formula : standard -dimensional Wiener process / Brownian motion Assume additive noise: is constant

23. Strong convergence Definition: A one-step method with step size is said to converge in the strong sense to a stochastic process at time provided Further, we say the convergence is of order provided there exists constants such that Reduces to usual def. of convergence when To test convergence rates numerically, we take M many samples (M>>1) of X and Y and approximate the error:

24. Strong Itô-Taylor expansion Analysis will involve comparison with the Itô-Taylor formula. Notation: step size , independent Brownian increments , and independent stochastic integrals Order 1 and : We derived an order 3. Higher order stochastic integrals are involved:

25. Existing methods References: for example, (1) Leimkuhler & Matthews, Molecular Dynamics (2014) and (2) Kloeden & Platen, Num. Sol. Of SDE’s (1991) Euler Maruyama (order 1): Milstein (order 1) and Itô-Taylor method (order 1): reduce to the same method as Euler Maruyama for Langevin dynamics with additive noise (i.e. constant diffusion ) Itô-Taylor method (order 2): requires computing Stochastic velocity Verlet (SVV): no and only one evaluation. No clear generalization to a third order method. Refer to Swope, Anderson, Berens & Wilson (1982). Theorem: (Telatovich and Li) SVV has strong order 2. Itô-Taylor method (order 3): requires computing and

26. Existing methods (continued) Direct splitting methods: split into Solve these simpler linear systems exactly, back-to-back Solution expansion involves Lie brackets High order convergence only in the weak sense Can do more sophisticated direct splittings Direct splitting methods lack stochastic integrals necessary for higher order convergence in strong sense Theorem:(Telatovich & Li) The splitting methods have strong order at most 1. Challenge:Higher order schemes require not only more smoothness but also more information about W!

27. Higher order operator-splitting methods Starting point: Kunita’s solution operator (Kunita, 1980). is given by: Vector fields Lie bracket Third order brackets and integrals defined similarly ODE special case (no noise ): is a simple finite sum (comprehensive exam) We can make truncations of

28. Truncations of the operator Truncation I: Truncation II: Truncation III:

29. Multiple stochastic integrals For truncations I, II, and III (Misawa, 2000): Joint symmetric covariances (T.): Computed using Itô’s isometry

30. Reduced Differential Equations Once sampled, we can treat as constants Rescale and solve from 0 to instead of 0 to Approximate, still non-linear due to f Approximate numerically by standard splitting techniques

31. Operator-splitting methods Useful integration method for ODEs. Order 1 naïve splitting: Order 2 symmetric/Strang splitting: Order 4 Neri splitting: where and

32. Baker-Campbell-Hausdorff (BCH) formula Useful for obtaining a “local” time estimate for the strong error. Theorem: If A and B are (non-commuting) matrices and , then , where is the matrix exponential, and the coefficients are given by, , , and . Proof: Follow the definitions above. For the formula in full generality, with a description of every coefficient, see Hairer, Lubich & Wanner (2005).

33. Doob’s inequality (special case) Useful for obtaining a uniform “global” time estimate on the strong error. Theorem (Doob): If is a martingale then

34. Truncation I Recall: Split where Notice: and These equations can be solved exactly. Theorem: The naïve, symmetric, and Nerisplittings give consistent methods with order 1 convergence. Proof (sketch): Use BCH formula, compare with Itô-Taylor expansions, and use Doob.

35. Truncation II Recall: Split where Notice: and These equations can be solved exactly. Theorem: (T. & Li) The naïve, symmetric, and Nerisplittings give consistent methods with convergence orders 1, 2, and 2, respectively. Proof (sketch): Use the same strategy as for truncation I.

36. Truncation III Recall: Splitwhere Notice: and Numerical evidence: truncation III with the Neri splitting gives strong order 3 Analysis is complicated

37. Numerical tests One dimensional pendulum: Lennard-Jones cluster (LJ7), function of distance:

38. Trajectory comparison One trajectory using three different operator splitting methods. Constants ; initial conditions . Same driving Wiener process with small . Exact solution: Euler Maruyama with step size . Approximate solutions: use larger step size .

39. Numerical tests