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
An overview of SDE based thermostats Gentle thermostats & NEMD Adaptive thermostats, polymer models and noisy gradients Ensemble preconditioning
Ben Leimkuhler
Stochastic Dynamics out of Equilibrium Research School, CIRM, Marseille, 2017
From L., Generalized Bulgak-Kusnezov Thermostats, PRE, 2010
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
Gibbs distribution Overdamped Langevin Nosé-Hoover Langevin Dynamics Nosé-Hoover Langevin Generalized Bulgac-Kusnezov Preconditioned Methods Ensemble Quasi-Newton
dx = [J(x) + S(x)]r log π(x) + r · [J(x) + S(x)] + p 2S(x)dW
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?
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
in order to achieve thermal equilibriation.
example with momentum constraints for shear flow applications and their use in machine learning.
use of an ensemble of “particles” (walkers).
The thermostats can be combined in most cases without altering their effectiveness (often improving it).
f ρ = 0
g ρ = 0
f+gρ = 0
Works for SDEs too…
Many schemes make good use of the concept of extension
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
Conservative dynamical system
Preserves Ex: Langevin dynamics
Ornstein-Uhlenbeck SDE
ρβ = e−βpT M −1pe−βU
Ergodic for
The ergodicity of the OU process implies ergodicity
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
Seminar, Cambridge University 1997: Nosé Dynamics
Too fast: balls move to outside,
reducing pressure, reducing speed Too slow: balls fall to inside, closing valve, leading to an increase in pressure, increasing speed James Watt’s Engine Nose-Hoover dynamics - a “Gibbs Governor” Preserves
e−β[p2/2+U(q)] × e−βξ2/2
Undergraduate research project of Josiah Willard Gibbs
µ = 2 µ = 1 µ = 1/2 µ = 4 A l l W r
g !
Stochastic version: Nosé-Hoover-Langevin dynamics ‘Histograms’ matches theoretical behavior
scalar OU process
˙ q = M −1p ˙ p = rU ˙ p = −ξp ˙ ξ = pT M −1p − nkBT e−βpT M −1p/2e−βU(q) e−βpT M −1p/2e−βξ2/2 e−βξ2/2
preserves preserves preserves ergodically
dξ = −γξdt + p 2β−1γdWt
dq = M −1pdt dp = rU(q)dt ξpdt dξ = [pT M −1p nkBT]dt γξdt + p 2β−1γdWt
NHL is clearly compatible with an extended Gibbs distribution meaning that
NHL[ρβe−βξ2/2] = 0
We can also prove it is ergodic by using the theory developed for Langevin dynamics and explained in the previous lectures. [L., Noorizadeh, Theil 2009]
Example: clamped harmonic spring chain [Mukamel’s talk]
q, p ∈ Rd
[L., Noorizadeh, Theil 2009]
NH NHL butane molecule
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]
Similar (but less smooth) for “Stochastic Velocity Rescaling”
0.1 0.004 Langevin NHL VAF Error VAF Error
Stochastic Velocity Rescaling [Bussi, Donadio, Parinello] dqi = −∂H ∂pi dt kinetic energy uses multiplicative noise
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
β = −0.006 β = 0.01
real ‘bath’ of 100 weak vortices a 1 variable bath
Dubinkina, Frank, L. MMS 2010
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.
Superficially similar to thermostatted MD: Newton’s equations + perturbation
graining of MD MD - Lennard-Jones DPD
(George Karniadakis) tDPD,cDPD, aDPF, sDPD, mDPD, eDPD
[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.
Fluctuation-dissipation: Momentum-conserving Langevin 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.
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.
OBAB (Shardlow) Others: DPD-Trotter = A(B+O)A (Coveney et al) Also Lowe-Anderson, NHLA, …
Nose-Hoover-Like (+ “gentle noise”)
PNHL L. & Shang, JCP, 2015 kinetic energy control A method or DPD at low friction or for NEMD with momentum conservation.
PNHL-S = ABCDODCBA PNHL-N = ABCDODCAB not symmetric! due to lack of symmetry, this method requires TWO force evaluations per timestep, instead of one.
All the schemes are convergent, giving the same RDF at small stepsize.
As the stepsize is increased, significant differences appear among the different methods
500 DPD particles
PNHL-N exhibits 2nd order convergence, but it is a non-symmetric method! Observed 2nd Order Convergence for PNHL-N
Due to cancellations, PNHL-N leads to: Invariant distribution: For large , the aux variable plays no role.
Reduce the perturbation equation using the projection operator Allowing to show: Proposition: PNHL-N is 2nd order accurate for all observables of the form
Critical Stepsize: Stepsize to achieve a 10% relative error in config temp
Posterior probability density (from Bayes’ Theorem):
Understand choice of parameters q given observations X Use Maximum Likelihood Estimate/“Subsampling”: ˜ N << N log p(X|q) ≈ N ˜ N
˜ N
X
i=1
log p(xi|q) X = {x1, x2, . . . xN} p(X|q)p(q) = p(q|X) model prior
D a t a
Alan Turing Institute to research big data
Nate Silver’s Nowcast (fivethirtyeight.com) Bayesian Election Forecasting The ultra-nerd October 2016
Variational Auto-Encoders and other “Deep” Generative Models in machine learning define smooth maps
(random) inputs sampled from prior parameters
[Graham & Storkey, Asymptotically exact inference in deep generative models and differentiable simulators, Arxiv, July 2016] problem: explore input space consistent with
sample the parameters.
=> Stochastic, constrained model & geodesic integrator
1719
Student Roll, University of Edinburgh
To sample a parameter distribution, apply Langevin dynamics with the potential computed from log p(X|q) ≈ N ˜ N
˜ N
X
i=1
log p(xi|q) This suggests a force polluted by a random field Special case (but usually assumed)
σ unknown.
Bayesian Learning via Stochastic Gradient Langevin Dynamics, M. Welling, & Y.-W. Teh, International Conference on Machine Learning, 2010
Like Euler-Maruyama for the modified SDE
This is problemmatic - If we don’t know the perturbation we can’t control the distribution. Solution of Welling and Teh: h = hn → 0 Not a super idea! As the sampling proceeds we need to shrink the stepsize, slowing exploration. How to decide the timestep schedule?
Noise Perturbation Negative Feedback Control
Use negative feedback loop control to stabilize the system against force perturbation (even unknown) Jones & L., J. Chemical Physics, 2011
driving thermostat steady state L † (ρβ(q, p) × ˆ ρ(ξ)) = 0, ˆ ρ(ξ) = exp
Nose-Hoover can automatically thermostat a simple stochastic driving perturbation. Jones & L., J. Chem. Phys., 2011
2014 NIPS Dear Bob…. A victim of google!
noise due to force error additional driving noise boxed term must be taken as a unit ergodicity by Hörmander, etc.
σ2
F = σ2∆t
large thermal mass limit BADODAB Under strong driving
superconvergence Role of thermal mass
f: logistic function
data e.g. voting intention covariates e.g. age, income, …
posterior parameter distribution
Gaussian prior
GOOGLE O u r M e t h
revenge of the nerd
Logistic Regression with 100 parameters Shang, Zhu, Leimkuhler, Storkey, NIPS, 2015 CCAdL (covariance controlled adaptive Langevin)
Dissipative Particle Dynamics simulation under shear flow Understanding particle margination in blood flow —- a step toward optimised drug delivery systems, K Mueller, D. Fedosov, and G. Gompper, 2015 Predicting human blood viscosity in silico, D. Fedosov et al, PNAS, 2011
Implemented using Lees-Edwards Boundary Conditions
“Adaptive variant of DPD” PAdL L. & Shang, JCP 2016 Adaptively parameterizes DPD, i.e. selects the right value of so that the fluctuation-dissipation relation is satisfied Idea: apply this for non-equilibrium DPD simulation, viewing non-equilibrium forces as continual perturbation of conservative ones.
PAdL-S = ABODOBA Better than alternative orderings….
PAdL mimics DPD VAFs
low friction high friction
Kremer-Grest Melts: joint with Xiaocheng Shang and Martin Kröger
Lees-Edwards Boundary Conditions Shear rates between 0.01 and 1 (above this bonds break) Possible Application: e.g. Shear Banding simulation M = 20 chains of length N = 30, density = 0.84. box size = 8.939 (unentangled).
PAdL: much faster convergence (to similar accuracy) than Langevin in equilibrium setting (PAdL 5-6 times more efficient!) DPD - Much much slower than these in our tests.
Langevin PAdL
Orientation relaxation: PAdL is independent of friction
Low Shear Rate High Shear Rate
figure of merit = Integrated Autocorrelation Time (IAT)
We would like to have as small as possible
Motivating Example: Eigenvalues: For MCMC schemes like Euler-Maruyama or Leimkuhler-Matthews, stability requires But for Poor Scaling Slow Convergence
More generally, we wish to sample problems with complicated energy functions, where each basin or local approximation may be very poorly scaled. Related Concepts BFGS Method Stochastic Newton schemes MC Hammer (Goodman and Weare) Compare to: Riemannian Manifold HMC (Girolami et al)
Use local information to estimate inverse Hessian matrix; precondition (rescale) dynamics to enhance convergence (reduce IAT)
Wishlist:
approximation
Idea: Use a collection of “walkers” to generate local covariance information and use this to estimate the inverse Hessian adaptively
Use an ensemble of L walkers:
each walker samples the same target distribution
We construct an dynamics in the extended space and compute ensemble averages by marginalisation
collects covariance info
walkers blend with identity (robustness) Basing Bi on other walkers only eliminates the problems of multiplicative noise
Discretization of SDEs: similar to BAOAB
“Adventures in Stamp Collecting”
But the eigensystems are different in different basins, so the localized covariance is needed…
Integrated Autocorrelation Times of Different Schemes
Ensemble Quasi-Newton python package http://bitbucket.org/c_matthews/ensembleqn