Constant and inconstant temperature molecular dynamics

1
Constant and inconstant temperature molecular
dynamics

• Ben Leimkuhler
• Centre for Mathematical Modelling
• University of Leicester
• work supported by EPSRC

CIMMS/IPAM Workshop 2002
2
Outline

• Statistical Mechanics from Pseudo-
• Microcanonical Trajectories
• Nosé Thermostats and Variants

• Generalized Thermostatting Baths

• Generalized Distributional Dynamics
• New approaches based on modified ensembles

• Applications work in Progress

3
Classical Molecular Models
Cyclosporin -- an immunosuppresant used in organ
transplants -- changes shape during delivery -
simulation may help reveal unlikely candidates
for new drugs...
Gigatom MD simulation atomic accuracy helps to
clarify the ductile-brittle transition
4
Key Simulation Challenges

• Dynamics

2. Sampling
3. Minimization
5
Improving Sampling, Accelerating Dynamics
Bond Stretch potential Stiff spring with rest
length ? replace with rigid rod
constraint. removes the highest frequency
components from the dynamics.
Multiple Scale Dynamics slow and fast variables
or slow and fast forces
Constant temperature and non-constant temperature
dynamics - heating and cooling, annealing,
barrier melts, Tsallis statistics
Voter dynamics transition state theory, reaction
paths, kinetic monte-carlo, Onsager-Machlup
action
6
Characteristic Features of MD Simulations
Very long time intervals
Importance of rare events
Simulation stepsize dominated by stability
demandsnot by accuracy demands
Work dominated by force computation
Computable quantites of interest are averages
with respect to some statistical mechanical
ensemble (especially the canonical ensemble)
7
Why Does MD Work?
Standard error analysis produces a bound of the
form error lt K exp(? T) h p where ? ? largest
Lyapunov exponent gt0. Since the time interval T
is very large, this bound is also large.
Traditional numerical analysis has little to say
about the validity of molecular simulation.
Instead What is needed is a backward error
analysis based on the fact that MD integrators
are symplectic they mimic classical mechanics
8
Properties of SymplecticDiscretization
Symplectic discretization of system with
Hamiltonian H is equivalent to solving a
perturbed Hamiltonian system with Hamiltonian Hh.
Backward Error Analysis Symplectic treatment
of classical mechanics ? Discrete Classical
Mechanics
9
Dynamic Sampling
Microcanonical sampling dynamics average
F(q,p) along a constant energy trajectory
Relies on
Ergodicity assumption time-average
spatial-average
Backward error analysis average along discrete
trajectory spatial average w.r.t. to
perturbed ensemble r dHh E
10
Canonical Sampling
Canonical sampling dynamics average F(q,p)
with respect to rr(q,p)exp(-bH)

• Periodic Boundary Conditions
• (allows for constant Volume simulation),
• and
• 2. Nosé Dynamics

Simulation of a small system exchanging energy
with a larger system at a given fixed temperature
T Temperature time-average of the k.e. per
particle
11
Nosés Dynamical Thermostat
Introduce a new variable s with momentum ps, and
define
HNose pTM -1p/2s2 ps2 /2Q V(q) g kBT log
s
Nosé showed that, for g N1
const ? òò exp(-b H) dq dp òòòò dHNose - E
dq dp ds dps  
canonical ensemble from pseudo-microcanonical
dynamics
12
Separating TransformationL. 2001
By using time and coordinate transformations we
can separate kinetic and potential parts --- Nosé
with a flat metric.
HNose pTM -1p/2s2 ps2 /2Q V(q) g kBT log
s
1. s eq , ps e-q pq Symplectic
coordinate change 2. dt/dt e2q
Poincaré time transformation
G pTM-1p/2 pq2/2Q e2q(V(q) gkTq E)
T(p, pq) V(q,q T )
13
Harmonic Oscillator
Nosé Potentials
14
Double Well
low T
low T
high T
high T
15
Lennard-JonesOscillator Chain
A
A
Velocity Autocorrelation Function
Recovered by binning the data after interpolation
of the solution to invert the effect of the
incorrect time transformation dt/dt s2
16
TP
HTP pTM -1p/(2s2 m2/3) ps2 /(2Q m2/3 ) pm2
/2Qm s2 V(m1/3 q) g kBT log s Pm - kBT
logm/3

• (s, ps) (eq , e-q pq ) coordinate change
• (m, pm) (y3/2, (2/3) y-1/2 py) coordinate
  change
• dt/dt exp(2q) y time transformation

HTP ? Separated form
17
Symplectic Time Rescaling Method
Nosé Poincaré Bond, Laird L. 1998 Based on
Poincaré time transformation dt/dt
s H(q,p)?s (H - E) and direct semi-explicit
symplectic discretization. Nosé 2000
splitting methods based on this... NP is as
efficient as methods based on Nosé-Hoover NP is
more stable than Nosé-Hoover Bond 98, Laird and
Sturgeon 99 and....NP preserves momentum and
angular momentum...and provides a non-stochastic
alternative to dissipative particle dynamics (DPD)
18
Generalized Baths Laird and L. 2002
Nearly any system can be used as a heat bath
for any other system, given the right coupling
term H(q,p) Primary System G(q,pq) Arbitr
ary Bath System s,ps Coupling
variables Generalized real-time Hamiltonian for
canonical sampling with an arbitrary bath
HGN s (H(q,p/x) G(q,pq) a(q)ps2/2 gkT ln
s - E)
19
Generalized Baths
Theorem on Generalized Baths

• Assuming mild conditions of regularity and
  convergence of the integrals, with gN, define
• Then
• òòòòòò ds(HGN - E) dq dp ds dps dq dpq
• const ? òò exp(-bH(q,p)) dq dp

HGN H(q,p/s) G(q,p) a(q)ps2/2 gkT ln s

• Numerics Various approaches based on splitting,
  e.g.
• sH(q,p/s) sG(q,p) s(a(q)ps2/2 gkT ln s
  - E)
• And semi-explicit treatment of resulting systems

20
Proof of Theorem

• òòòòòò dHGN - E dq dp ds dps dqdpq
• òòòòòò dH(q,p/s)a(q)ps2 /2QG(q,pq )g
  kBT log s-E dq dpq
• òòòòòò dH(q,p)a(q)ps2 /2QG(q,pq )g kBT
  log s-EsNdqdpq
• òòòòò (ò dF(sq,p, q,pq, ps)sN ds ) dqdp
  dsdps dqdpq
• òòòòò F-1(0)N / F(F-1(0)) dq dp dps
  dqdpq
• òòòòò exp-(N1) (H(q,p)a(q)ps2 /2QG(q,pq
  ))/ g kBT dqdp dps dqdpq
• const ? òò exp-(N1) H(q,p)/ g kBT dq
  dp
• so gN ? canonical ensemble partition function

F(s q,p, ps)
21
Harmonic Oscillator
The hardest system to sample using dynamics is
the harmonic oscillator H(q,p) q2/2 p2/2
Configurational (q) distribution using a handful
of oscillators as a bath (with a bit of help from
the large timestep perturbed Hamiltonian)
Nosé
22
Generalized Baths
An important potential application for
generalized baths improved ergodicity Standard
approach Nose-Hoover chains, a cascade of
Nose-Hoover thermostats. The new approach is
potentially more general and more powerful,
allowing a broad class of ergodicity-enhancing
devices.
23
Generalized Density Dynamics
Barth, Laird, L. 2002
There are many situations where it would be
useful to be able to sample with respect to some
arbitrary given density F(q,p) with
? ? F(q,p)dqdp 1 and F(q,p)0
24
Motivation Potential Smoothing

• Many approaches to global minimization are based
  on some sort of potential modification replace
  the potential function by a smoothed alternative.
  Such potential modifications are also at the
  heart of Voters accelerated dynamics.
• If properly chosen, a potential smoothing V ?
  f(V) can sometimes simplify the approach to
  minimum energy configuration.
• We can view any such smoothing as being
  associated to a corresponding energetic
  transformation H ? f(H)
• Canonical ensemble sampling of f(H) represents
  sampling in a modified distribution for H

25
Tsallis Entropy
Tsallis-Straub potential
Tsallis Statistics ??(q,p)(1?H(q,p)/kT)1/?

?? ? exp(-bH) as ? ? 0
26
Potential Modification
An Alternative
27
Effective Energy
For any given density expressed in terms of
energy, we can write r(H) exp(-bf(H)) for
some f. Here f(H) can be viewed as a modified
energy function sampling f(H) from the canonical
ensemble is equivalent to sampling H in the
density r(H). Question Can we find an efficient
numerical scheme for dynamics in the r(H)
ensemble?
cf. Plastino and Antaneodo 97
28
Generalized Ensemble Dynamics
If we start with an arbitrary density F(p,q), we
can define Heff -b-1 ln F(p,q)
then we can apply Nose dynamics to sample,
using HeffNose -b-1 ln F(p/s,q) p2/2Q
gkT ln s or using a generalized bath. But the
straightforward approach to symplectic
integration leads to an implicit method...
29
Generalized Ensemble Dynamics
Two tractable cases 1. F(p,q) A(q)B(p) In
this case the method is like standard
canonical dynamics with a modified
separated potential and kinetic energy 2.
F(p,q) F(H(p,q)) In this case we need a
trick...
30
Generalized Ensemble Dynamics
If f is a smooth monotone function, the constant
E dynamics of the Hamiltonian f(H1(q,p))
H2(q,p) are from the zero energy dynamics of
H1(q,p) - f-1(E - H2(q,p)) together with a
time-transformation, dt/dt f An elegant tool
for developing effective numerical methods
31
Generalized Ensemble Dynamics
Example Nosé Dynamics
Equations of Motion
32
Semi-Explicit Numerical Method
33
Canonical Samplingand Simulation
We can use a modified ensemble as a tool for
canonical sampling. We first compute the
modified ensemble dynamics and then reweight
averages by a factor exp(b(f(H)-H)) If a
time-transformation is used as part of the
scheme, an additional factor of f must also be
included in the reweighting.
Example sampling a double-well potential using
a barrier melt
Coordinate q
34
Designer Ensembles
full Hamiltonian potential only
? near perfect recovery of canonical ensemble
Differences --- Hamiltonian approach vs.
potential smoothing significantly more time is
spent in high energy regions.
35
Designer Ensembles
36
Designer Ensembles
37
Current Future Work

• Minimization via Dynamical Annealing with E.
  Barth, B. Laird
• Generalized Baths with B. Laird and C. Sweet
• Local Adaptive Thermostatting and Barrier Melts
• Application to protein simulations, drug
  transport/absorbtion w./A. Laaksonen (Stockholm),
  E. Barth (Kalamazoo) Pharmacia Corp.
• Applications to nucleation and cavitation
  processes w./A. Cocks and S. Gill

38
Alanine Dipeptide
39
Alanine Dipeptide Torsion Sampling
Modify only the flexible parts of the potential
40
Dihedral Sampling
41
Dihedral Sampling