Molecular Dynamics Lecture 4 Ben Leimkuhler An overview of SDE - PowerPoint PPT Presentation

Molecular Dynamics Lecture 4 Ben Leimkuhler An overview of SDE based thermostats Gentle thermostats & NEMD Adaptive thermostats, polymer models and noisy gradients Ensemble preconditioning Stochastic Dynamics out of Equilibrium Research School, CIRM, Marseille, 2017

From L., Generalized Bulgak-Kusnezov Thermostats, PRE, 2010

Finding the “Right” Dynamics for the Job There are many different stochastic models that can be used in MD, but they can have very different efficiencies for a particular task. Overdamped Langevin Dynamics great for sampling well scaled multivariate Gaussian distribution, lousy for a highly corrugated landscape Nosé-Hoover gentle - good for autocorrelation functions in systems with strong internal mixing properties… not ergodic - for nucleic acid simulations in implicit solvent

Thermostats Overdamped Langevin Gibbs distribution Langevin Dynamics Nosé-Hoover Langevin Nosé-Hoover Generalized Bulgac-Kusnezov Preconditioned Methods Ensemble Quasi-Newton

A generic sampling dynamics p d x = [ J ( x ) + S ( x )] r log π ( x ) + r · [ J ( x ) + S ( x )] + 2 S ( x )dW J(x) antisymmetric S(x) symmetric Includes, e.g., SDEs like Brownian and Langevin dynamics non-reversible perturbation methods various ensemble sampling schemes Questions: Which approach converges most rapidly ? (small IAT) What is the sampling bias under discretization ? How to effectively combine with extension ?

Generalized Sampling Up to know we have assumed the situation of a known distribution with invariant density What if we don’t know U or cannot exactly resolve the force? Multiscale models, e.g. ab initio MD Methods and QM/MM methods (heating due to force mismatch) Nonequilibrium MD (e.g Shear Flows ) Applications in Bayesian Inference & Big Data Analytics

Problems for today 1. How to gently perturb Hamiltonian dynamics in order to achieve thermal equilibriation. 2. How to handle driven systems efficiently, for example with momentum constraints for shear flow applications and their use in machine learning . 3. How to accelerate convergence to equilibrium by use of an ensemble of “particles” (walkers).

Additivity The thermostats can be combined in most cases without altering their effectiveness (often improving it). x = f ( x ) + g ( x ) ˙ L † f ρ = 0 ⇒ L † f + g ρ = 0 L † g ρ = 0 Works for SDEs too…

Extension Many schemes make good use of the concept of extension Z π ( x )˜ π ( y )d y ∝ π ( x ) This looks banal but the key point is that although x and y decouple in the invariant distribution, they may be tightly coupled in the associated SDEs. Example: Langevin dynamics Z e − β p 2 / 2 e − β U ( x ) d p ∝ e − β U ( x )

Remote control of thermal equilibration Ex: Langevin dynamics Conservative Ornstein-Uhlenbeck dynamical SDE system Preserves Ergodic for ρ = e − β p T M − 1 p ρ β = e − β p T M − 1 p e − β U ¯ • The two systems are both compatible with ρ β • Sufficient mixing The ergodicity of the OU process implies ergodicity of the full system

Gentle thermostats & NEMD

Where I learned about Nosé Dynamics Seminar, Cambridge University 1997: Nosé Dynamics Sir Henry Peter Francis Swinnerton-Dyer, 16th Baronet KBE FRS Number Theorist, Student of Littlewood, Polya and Sylvester Prizeholder Vice-Chancellor of Cambridge University 1973-83

James Watt’s Engine Too fast: balls move to outside, opening valve, releasing steam, reducing pressure, reducing speed Too slow: balls fall to inside, closing valve, leading to an increase in pressure, increasing speed Nose-Hoover dynamics - a “ Gibbs Governor ” q = p ˙ Preserves p = �r U ( q ) � ξ p ˙ e − β [ p 2 / 2+ U ( q )] × e − βξ 2 / 2 ξ = p 2 � kT ˙

Problems with the Gibbs Governor a. It doesn’t actually work. b. It’s not the Gibbs Governor. This is: Undergraduate research project of Josiah Willard Gibbs

Nosé-Hoover Dynamics for Harmonic Oscillator µ = 1 µ = 2 ! g n o r W l l A µ = 4 µ = 1 / 2

Nosé-Hoover Chains are not ergodic..

Stochastic version: Nosé-Hoover-Langevin dynamics q = p ˙ p = �r U ( q ) � ξ p ˙ ξ = p 2 � kT � γξ + ση ( t ) ˙ scalar OU process ‘Histograms’ matches theoretical behavior

Nosé-Hoover-Langevin d q = M − 1 p d t d p = �r U ( q )d t � ξ p d t d ξ = [ p T M − 1 p � nk B T ]d t � γξ d t + p 2 β − 1 γ dW t q = M − 1 p e − β p T M − 1 p/ 2 e − β U ( q ) ˙ preserves p = �r U ˙ p = − ξ p ˙ e − β p T M − 1 p/ 2 e − βξ 2 / 2 preserves ξ = p T M − 1 p − nk B T ˙ e − βξ 2 / 2 preserves p 2 β − 1 γ dW t d ξ = − γξ d t + ergodically

Ergodicity of NHL [L., Noorizadeh, Theil 2009] NHL is clearly compatible with an extended Gibbs distribution meaning that NHL [ ρ β e − βξ 2 / 2 ] = 0 L † We can also prove it is ergodic by using the theory developed for Langevin dynamics and explained in the previous lectures.

Harmonic system w/o resonance [L., Noorizadeh, Theil 2009] H = p T M − 1 p + q T Bq q, p ∈ R d 2 2 A = M − 1 B, A ϕ k = ω k ϕ k ω k ̸ = ω l , k ̸ = l R 2 d +1 Theorem: NHL is ergodic on Example : clamped harmonic spring chain [Mukamel’s talk]

Harmonic system w/o resonance

butane molecule NH NHL

“Gentle” property of NH/NHL We can show that NHL is a “gentle” thermostat: dynamical properties are mildly perturbed for a given rate of convergence of kinetic energy. [ L., Noorizadeh and Penrose, J. Stat. Phys., 2011] VAF Error VAF Error 0.1 0.004 Langevin NHL Similar (but less smooth) for “Stochastic Velocity Rescaling” of G. Bussi, D. Donadio and M. Parinello

Stochastic Velocity Rescaling [Bussi, Donadio, Parinello] d q i = − ∂ H d t uses multiplicative noise ∂ p i kinetic energy

Autocorrelation functions (LJ System) Even the deterministic method seems to work for sufficiently complicated systems.

Temperature gradients and NEMD simulation using Nosé-Hoover [NONEQUILIBRIUM MOLECULAR DYNAMICS METHODS FOR LATTICE HEAT CONDUCTION CALCULATIONS Junichiro Shiomi, Ann. Rev. Heat Transfer , 2014] Theoretical Methods for Calculating the Lattice Thermal Conductivity of Minerals, Reviews in Mineralogy & Geochemistry Vol. 71 pp. 253-269, 2010

A thermostat for the Vortex Discretization Dubinkina, Frank, L. MMS 2010 β = 0 . 01 β = − 0 . 006 Radial density at fixed β Buhler ’02 real ‘bath’ of 100 weak vortices GBK Thermostat a 1 variable bath

Gentle Momentum Conserving Thermostat for NEMD In many examples it is crucial to preserve the momentum . This is of particular relevance in fluid dynamics, or in any situations where aggregate inertial effects are likely to play an important role. The problem is that momentum is distorted by Langevin and Nose Hoover dynamics. How to control momentum while gently restoring canonical equilibrium? The simple idea is to use pairwise forces satisfying Newton’s 3rd law that are projected onto the radial axis between particles.

Dissipative Particle Dynamics Superficially similar to thermostatted MD: Newton’s equations + perturbation • DPD typically relies on “softened” potential energy functions often ad-hoc but sometimes derived by systematic coarse- graining of MD MD - Lennard-Jones DPD • DPD typically viewed as a coarse-graining technique. • Always involves a thermostat. • Variants: DPD-e, QDPD, (George Karniadakis) tDPD,cDPD, aDPF, sDPD, mDPD, eDPD

Dissipative Particle Dynamics Momentum-conserving Langevin dynamics Fluctuation-dissipation: [1] P. Hoogerbrugge and J. Koelman. Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhysics Letters , 19(3):155, 1992. [2] P. Espanol and P. Warren. Statistical mechanics of dissipative particle dynamics. Europhysics Letters , 30(4):191, 1995.

Ergodicity of Dissipative Particle Dynamics Ergodicity: System has a unique smooth distribution which is a universal attractor. For an ergodic system, memory of initial conditions is eventually lost, and any path can be used to calculate averages. Open Question: Is DPD ergodic? Only proof is for a 1D system with high particle density (Shardlow and Yan) More general answer: Probably not! Simulation study of Pastorino et al 2007: Appears to contradict ergodicity in case of soft DPD potentials and reduced interactions.

DPD alternatives Lowe-Andersen: Allows control of the Schmidt number = ratio of kinematic viscosity to diffusion constant Momenta are updated according to conservative forces. Subsequently, each pair is (with fixed probability) updated with an added random kick. Peters: All particles perform a random step after conservative updating, with collision coefficient chosen to mimic the Lowe-Andersen collision rate. Nosé-Hoover-Lowe-Andersen: Ad-hoc & does not reproduce the canonical ensemble.

Integrators for DPD OBAB (Shardlow) Others: DPD-Trotter = A(B+O)A (Coveney et al) Also Lowe-Anderson, NHLA, …

Pairwise NHL A gentle momentum-conserving thermostat A method or DPD at low friction or for NEMD with momentum conservation. PNHL L. & Shang, JCP, 2015 Nose-Hoover-Like (+ “gentle noise”) kinetic energy control