Sampling Free Energy Surfaces by MD Marcella Iannuzzi Department - - PowerPoint PPT Presentation
3rd CP2K tutorial: Enabling the Power of Imagination in MD Simulations June 17-21 2013, Zrich Sampling Free Energy Surfaces by MD Marcella Iannuzzi Department of Chemistry, University of Zurich http://www.cp2k.org OUTLINE Free energy
3rd CP2K tutorial: Enabling the Power of Imagination in MD Simulations June 17-21 2013, Zürich
2
3
1.5 2 2.5 3 3.5 4 4.5 5
Predictive power frustrated by sampling fast degrees of freedom with time-steps from < 0.1 fs (CPMD) up to 1 fs (MM)
Choose a suitable model of the system Determine the thermodynamic conditions ⇒ Ensemble averages Equilibrium sampling of physical quantities
4
Phase Transitions, Conformational Rearrangements, Chemical Reactions, Nucleation, Diffusion, Growth, etc.
Activation Energies Minimum energy pathways
Complex and high dimensional configurational space Multitude of unknown intermediates and products Unforeseen processes, many irrelevant transition states Intrinsically multidimensional
Entropic bottlenecks
Exploration of configurational space
5
A system of N particles in a volume V is completely determined through its Hamiltonian
Choice of ensemble: portion of phase space, microstates Difference in averages and fluctuations
NVE-P total energy and linear momentum are constant of motion
H({RI}, {PI}) =
P2
I
2MI + U(RI})
Counters, blue on one side and green on the other, 6×6 checkerboard Each pattern is a microstate Subset belonging to macrostate “15 green”
6
The Laplace transform of the density of state
Q(N, V, T) = Z exp(−βE)Ω(N, V, E)dE Probability of the macrostate at a given T
Q(N, V, T) = 1 N!h3N Z exp ⇥ −βH(rN, pN) ⇤ drNdpN = 1 Λ(β)3NN!Z(N, V, T)
analytic kinetic part is integrated out
Z(N, V, T) =
configurational partition function
7
Helmholtz free energy or thermodynamic potential
A = − 1 β ln Q(N, V, T)
Thermodynamics Statistical Mechanics
∆A = − 1 β ln Z1 Z0 ⇥ Q0 ∝
e−βH(r,p)drdp
Macroscopic state 0 corresponds to a portion of the phase space : Γ0
Q0 ∝
e−βH0(r,p)drdp
Macroscopic state 0 corresponds to H0
Q0 ∝
e−β0H(r,p)drdp
Macroscopic state 0 corresponds to a value of a macroscopic parameter, e.g T
entropic and enthalpic contributions
8
Reference (0) and target system (1)
H1(r, p) = H0(r, p) + ∆H(r, p)
Probability of finding (0) in configuration (r,p)
P0(r, p) = e−βH0(r,p)
∆A = − 1 β ln R e−βH1drNdpN R e−βH0drNdpN = − 1 β ln R e−βH0e−β∆HdrNdpN R e−βH0drNdpN
Free energy difference
∆A = − 1 β ln ⌦ exp ⇥ −β∆H(rN, pN) ⇤↵
Integrating out the analytic kinetic part
⇥F(r, p)⇤1 = ⇥Fe−β∆U⇤0 ⇥e−β∆U⇤0 ∆A0,1 = 1 β ln⇥e−β∆U⇤0
9
Accuracy ⇒ target and reference systems are similar⇒ overlapping regions insufficient statistics or incomplete overlap ⇒ enhanced sampling
Γ 1 Γ 1 Γ 1
∆A = − 1 β ln Z exp [−β∆U] P0(∆U)d∆U
−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 DU 0.0 0.4 0.8 1.2 1.6 2.0 2.4 P(DU )
P0(DU )x
exp(−bDU ) P0(DU )
exp(−bDU ) }
Shifted function Low-𝚬U tail is poorly sampled low statistical accuracy but important contribution to 𝚬A
10
Protein ligand binding, host-guest chemistry, solvation properties, ...
Transformation (0)→(1) as series of non-physical, intermediate states along a pathway characterized by the “coupling parameter” λ Free energy as continuous function of λ through H(λ)
H(λ) = Henv + λH0 + (1 − λ)H1
Point mutation of alanine into serine: coexistence without seeing each other
L 1 L 2 L1 + R L 2 + R D A1 D A4 D A2 D A3
11
Variables chosen to describe changes in the system Reaction coordinate : the order parameter corresponds to the pathway along which the transformation occur in nature
Different possible definitions of OP
Set up of system with desired values of OP
annihilation non-bonded torsion angle
12
Select parameters, continuous functions of coordinates
Ωξ(N, V, E, ξ) = Z δ[U(rN) − E] ⇣ Πiδ[ˆ ξi(rN) − ξi] ⌘ drN
ξ = {ξi} ˆ ξi(rN) Qξ(N, V, T, ξ) = Z e−βU(rN) ⇣ Πiδ[ˆ ξi(rN) − ξi] ⌘ drN Aξ = − 1 β ln Qξ
Density of States Canonical Partition Function Free Energy
must distinguish among metastable states
ˆ ξi(rN)
13
1 2 3 4 5 6 7
60 120 180 240
A (kcal/mol) Θ (degrees)
Free energy butane isomerization C-C-C-C
Probability distribution of the order parameter
∆A(ξ) = A(ξ1) − A(ξ0) = −β−1 ln P(ξ1) P(ξ0)
Histogram of M bins
δξ = (ξ1 − ξ0)/M
P(ξ0 + (i − 0.5)δξ) = fi P
j fj
Not uniform sampling Restrain the system within a window by harmonic potential
τ = Lτw ∝ (ξ1 − ξ0)2 LDξ
14
Non-Boltzmann sampling to enhance the probability of important regions
P(w)(ξ, T) = R w(ˆ ξ(rN)) exp[−βU(rN)]δ[ξ − ˆ ξ(rN)]drN R w(ˆ ξ(rN)) exp[−βU(rN)]drN
positive bias w(ξ) Free energy differences ∆A(w)(ξ) = −β−1 ln w(ξ0)P(w)(ξ, T) w(ξ)P(w)(ξ0, T) = −β−1 ln P(w)(ξ, T) P(w0)(ξ, T) + ln w(ξ0) − ln w(ξ)
w(ξ) = exp[−βV (ξ)]
U(b)(rN) = U(rN) + V (ξ(rN))
Umbrella Potential connecting different regions of the phase space
15
Modify the underlying potential to obtain a uniform sampling : Vb(s) = -ΔA(s) First guess & iterative improvement
P(b)(s) = 1 Z(b) Z dx exp (−β(U(x) + Vb(s(x)))) δ(s − s(x)) = Z Z(b) exp (−βVb(s)) 1 Z Z dx exp (−βU(x)) δ(s − s(x)) = Z Z(b) exp (−β(Vb(s) + A(s)))
Probability in the presence of the bias & un-biasing
A(s) = − 1 β ln P(b)(s) − Vb(s) − 1 β ln ✓ Z Z(b) ◆ U (p)(x) = U(x) + Vb(s(x))
Free Energy Free Energy
good Vb bad Vb
as flat as possible
∆A(b)(s)
∆A(s)
−Vb(s)
16
Capture the essential physics include all relevant DoF and properly describes the dynamics
q distinguishes between A and B but might fail in describing essential aspects of the transition Discriminate configurations of reactants and products and intermediates
mechanisms of transition Reversibility
DoF (no relevant bottlenecks)
17
(b) (a)
x x* x′ x A(x ) B A
Not including important DoF leads to wrong interpretation
18
Distance Angle Dihedral Difference of distances
Derivable function of the degrees of freedom to which a given value can be assigned
⇧ ⌃
Coordination number
D[klm]
L1L2 =
1 NL1
di · ˆ v[klm] − 1 NL2
dj · ˆ v[klm] |RI − RJ|
θ(RI, RJ, Rk)
Θ(RI, RJ, Rk, RL) |RI − RJ| −| RJ − RK| CL1L2 = 1 NL1
NL1 j=1
⇤ ⌥ ⇧
NL2 i=1
1 −
r0
⇥n 1 −
r0
⇥m ⌅
RMSD
19
Knowing end-points or a full approximate path (NEB)
B A
s(R) =
[Rk(R)]nkck s(R) = 1 P − 1 P
k=1 k e−λ||S(R)−S(k)||2
P
k=1 e−λ||S(R)−S(k)||2
z(R) = 1 λ ln P ⇤
k=1
e−λ||S(R)−S(k)||2 ⇥
Rk(R) = ⇥Ndist
j
(dj(R) − d(k)
j (R))2
Ndist Rk(R) = ⇥
i(Ri − R(k) i
)2 N
20
Ramachandran plot competitive paths linear interpolation for the initial path
21
In SUBSYS add one section per CV
&COLVAR &COORDINATION KINDS_FROM N KINDS_TO O R_0 [angstrom] 1.8 NN 8 ND 14 &END COORDINATION &END COLVAR &COLVAR &DISTANCE_FUNCTION ATOMS 4 6 6 1 COEFFICIENT -1.00000 # distance 1 = ( 4 - 6 ) # distance 2 = ( 6 - 1 ) &END DISTANCE_FUNCTION &END COLVAR &COLVAR &DISTANCE AXIS X ATOMS 1 4 &END DISTANCE &END COLVAR &COLVAR &RMSD &FRAME COORD_FILE_NAME planar.xyz &END &FRAME COORD_FILE_NAME cage.xyz &END SUBSET_TYPE LIST ATOMS 1 3 5 6 8 9 ALIGN_FRAMES T &END &END
22
&CONSTRAINT &COLLECTIVE COLVAR 1 INTERMOLECULAR TARGET 5. TARGET_GROWTH 1.1 TARGET_LIMIT 10. &END COLLECTIVE &END CONSTRAINT
&COLLECTIVE TARGET [deg] 0.0 MOLECULE 1 COLVAR 1 &RESTRAINT K [kcalmol] 4.90 &END &END COLLECTIVE
In MOTION add one section per constraint
23
To freeze fast degrees of freedom and increase the time step: e.g., intra-molecular bonds by distance constraints To explore only a sub-region of the conformational space As reaction coordinates : constrained dynamics and thermodynamic integration To prevent specific events or reactions
Implicit functions of the degrees of freedom of the system Lagrange formulation for simple constraints, functions of RI
σ({RI}, h, Ψ) = 0 ˙ σ({RI}, h, Ψ) = 0 L({RI}, {PI}) = L({RI}, {PI}) −
λασ({RI})
24
First update of velocities (first half step) and positions
Calculation of the new forces FI(t+∂t) Update of velocity (second half step) Correction by the constraint forces Modified velocity Verlet scheme by additional constraint forces
H.C. Andersen, J. Comp. Phys., 52, 24 (1983)
V
I = VI(t) +
δt 2MI FI(t)
R
I = RI(t) + δtV I
RI(t + δt) = R
I + δt2
2MI g(p)
I (t)
VI(t + δt) = V
I +
δt 2MI [FI(t + δt) + g(v)
I (t + δt)]
Constraint Forces
g(v)
I (t) = −
λ(v)
α
∂σα({RI}) ∂RI
g(p)
I (t) = −
λ(p)
α
∂σα({RI}) ∂RI
Set of non-linear equations solved iteratively
∂σα ∂RI VI =
∂σα ∂RI · V′
I +
δt2 2MI ∂σα ∂RI ∂σβ ∂RI
β = 0
eα({λγ}) = −
J−1
αβ σβ({λγ})
Jαβ = ∂σα({λγ}) ∂λβ
25
A(ξ1) − A(ξ0) = ξ1
ξ0
dA dξ dξ
along a one dimensional generalized coordinate ξ(x) Path-independent
dA dξ = ∂H
∂ξ e−βHdq1..dqN−1dpξ..dpq N−1
N−1
(x, p) ⇒ (ξ, q1..qN−1, pξ..pq
N−1)
generalized coordinate to simplify derivative
= ∂H ∂ξ ⇥
ξ force at ξ, averaged over fluctuations of other DoF instantaneous force on ξ
Potential of Mean Force exerted on ξ
∂H ∂ξ ⇥
ξ
= ∂U ∂ξ − 1 β ∂ ln |J| ∂ξ ⇥ [J]ij = ∂xi ∂qj
mechanical + entropic
26
Hλ = H + λ(ξ − ξ(r))
F(S) S tMD1 tMD2 SMD1 SMD2 SMD3 MD1 MD1
fSMDn = −∂F(S) ∂S
ξ A(ξ) ξ1 ξ2 ξ3
λ⇥(ξ ξ(r))
Fixman Potential Series of constrained MD simulations
˙ ξ = 0 : pξ(q, pq)
un-constrained <..> = constrained corrected <..>F
Hλ
F = Hλ + 1
2β ln Zξ Zξ = ⇤
i
1 mi ∂ξ ∂xi ⇥2
independent MD replica
mean force acting on the system to hold ξ constant
dA dξ = λF ⇥F
ξ ˙ ξ
mean force acting on ξ related to external force to hold ξ constant
27
Many estimates at ξ to reduce statistical errors Many windows to get accurate integrals May not be easy to prepare by hand the system at given ξ Different possible pathways: the starting configuration selects one path, but crossing is rare, ξ(r)=ξ partially sampled
Multidimensionality (more coordinates) too expensive MD performed at fixed ξ, collecting statistics of the force acting
28
Order Parameter : coordination of a specific donor phosphorane O (r0=1.3 A) axial equa. gradually reducing n corrected MF ∆A(n) de-protonation and formation of Zundel ion H3O+ breaks loose chain of proton transfers protonation of axial OH formation of H3PO4 and subsequent de-protonation equivalent PO in solvated P(OH)5 Pronounced equatorphilicity of O- in PO5H4 -
29
Running M replica at different T high T to explore large part of the phase space low T to sample precisely local minima
T1 < T2 < .. < TM
exchange ensures access to all local minima at low T
Equilibration + swap attempt + rescaling of velocities Swapping acceptance Energy histogram
σ ∝ (N)1/2
likelihood that two replicas have overlap in phase space M and ∆T strongly affect sampling efficiency and computational costs
αIJ = min
30
Penta-peptide system in gas phase(Try-Gly-Gly-Phe-Met), all-atom AMBER Eight temperatures replica exchange simulation, over 1 ns
Regular NVT MD Replica exchange 200 K 700 K
31
Statistical, reaction-coordinate free description of pathways connecting long-lived stable states z(τ)≣{z0,z∆t,...zτ} The probability distribution in the pathway ensemble depends on ρ(z0) and the dynamics
definition of multidimensional ξ(r) to identify A and B large enough to contain equilibrium fluctuations no overlap with opposite basin of attraction
P[z(τ)] = ρ(z0)
τ/∆t−1
p(zi∆t → z(i+1)∆t)
r ρ(z0)
Markovian processes ξ1 ξ2
32
Focused on the sub-ensemble of pathways containing barrier-crossing events
detailed balance criterion
identified and transition mechanisms are characterised Transition-path-ensemble Shooting Shifting pacc[zn(τ) → zo(τ)] = min
PAB[zo(τ)]pg[zn → zo] ⇥
Efficiency:
as possible pacc[zn ➔ zo] as large as possible
pacc[zn(τ) → zo(τ)] = hA(zn
0)hB(zn τ ) min
0)
ρ(zo
0)
⇥
PAB[z(τ)] = hA(z0)P[z(τ)]hB(zτ)/ZAB(τ)
33
Canonical equilibrium distribution under potential V(r) Choose a set of relevant Collective Variables S(r):{Sα(r)}, such that the process is well defined in the reduced space Σ(S)
into the configuration space Σ(S): meta-trajectory S(r(t)) Enhance the exploration by adding a penalty potential that discourages the system to visit already explored macro-states
P(S) = e−βA(S)
A(S) = − 1 β ln ⇤ dr e−βV (r)δ(S − S(r)) ⇥ VMTD(S(r), t) =
Wt⇥e− [S(r)S(rG(t⇥))]2
2∆S2
A Laio et al. Proc. Natl. Acad. Sci. U.S.A., 99, 12562 (2002) M Iannuzzi et al, PRL, 90, 238302 (2003)
34
Non Markovian Coarse-grained MD
Escape S F(S) V(t) t1 t2 t3 t4
Fill-up by successive contributions Keeping track of the explored space Direct estimate of FES topology Eliminate metastability and reconstruct A(S) within Σ(S)
AG(S, t) = −VMTD(S(r), t)
Flattening of free energy surface
W/τG → 0
δA(S) = A(S) − AG(S, t)
P(S) ∝ e−β[A(S)−AG(S,t)]
δA(S)⇥
C Micheletti et al, PRL, 92, 170601 (2004)
δtMD τG τs
35
Histogram of a trajectory generated by V+VMTD provides a direct estimate of free energy
36
M Iannuzzi et al, PRL, 90, 238302 (2003)
Introduction of auxiliary variables s : {sα} one for each Sα(r) with large enough M Enforcing adiabatic separation, other time scales and memory effects thermalization coupling to system DoF sampling enhancement
L = K − V (r) +
1 2Mα ˙ sα −
1 2kα(sα − Sα(r))2 − VG(s, t)
VG(s, t) =
Wt⇥e− (ssG(t⇥))2
2∆s2
Ak(s) = − 1 β ln ⇤ dr e−β[V (r)+ 1
2
P
α kα(sα−Sα(r))2]
⇥
For large t and slow deposition rate, VG approximates the free energy and the meta-trajectory s(t) follows the MEP
lim
k→∞ Ak(s) = A(s)
τs
37
Smooth and continuous meta-trajectories that follow the minimum energy pathway are obtained by a modified Velocity Verlet algorithm
Modified force on ions due to the coupling with the dynamics of the meta-variables Force on the meta-variable Vharm and VMTD generate opposite contributions The Metadynamics Lagrangian generates fictitious dynamics describing the most probable pathways
fα(t) = kα(Sα({RI(t)}) − sα(t)) − ∂ ∂sα VMTD(s, t) fI(t) = f (0)
I
(t) − kα(Sα({RI(t)}) − sα(t))∂Sα ∂RI
38
Different CVs, anisotropic A(s), different diffusion coefficients
τsα = Lα/Dα σα = sα/Lα
VG(s, t) =
Wt⇥e
− P
α (sαsα(t⇥))2 2(∆sα)2
Δs⊥ (input) Δs// t1 t2 t3 t4 tn Accumulate slices along the trajectory
MTD trajectory re-shaped Gaussian tube
VG(s, t) =
Wt⇥e− (ssG(t⇥))2
2∆s2
× e
− [(ssG(t⇥))·dG(t⇥)]2
2||dG(t⇥)||4
dG
39
C Micheletti et al, PRL, 92, 170601 (2004) A Laio et al, JPC-B 109, 6714 (2005)
Dynamics generating the equilibrium distribution associated with A(s)-AG(s,t)
(s, t) =
Averaging over many independent trajectories
¯ (t) =
100 300 500 700 900 0.2 0.6 1 number of Gaussians error
different A(s) profiles
¯ ⇥ = C(d)
D⇤G ttot = τG ⇥
Ω:A(s)<Amax ds(Amax − A(s))
(2π)d/2W
α ∆sα
too large ∆s would smear out A(s) details : ∆s/L<0.1
small Gaussians more frequently is better Empirical error estimate
¯ ⇥ ⇥ ⇥⇤s⇥ ¯ A ttot VΩ
40
The selected CV must discriminate among the relevant states (reactants, products, TS) The number of hills required to fill the well is proportional to 1/(∆)NCV The sampling of large variations of the CV over almost flat energy regions is expensive: diffusive behavior MTD is not the true dynamics. Reaction rates are derived a posteriori from the estimated FES The analysis of the trajectory is needed to isolate the TS With proper choices of CV and parameters, the MTD trajectory describes the most probable pathway taking into account also possible kinetic effects (lager and shallower channels are preferred) The accuracy in the evaluation of the FES depends on hills’ shape and size, and on the deposition rate. The ideal coverage VG ({S})= -A({S}) (flat surface)
42
Crystalline Silicon can assume different lattice structures. Phase transitions induced by external pressure are known from experiment Metadynamics is able to explore all the accessible metastable states, without requiring any over-pressurization
h = [a1, a2, a3]
Dia βtin SHex
Evolution of the cell parameters during the metadynamics run
43
Structural analysis by radial correlation function
44
Ar+ at 50 eV striking the rim
Time [fs]
Temperature [K] Ar N B Rh Partial Charge (Mulliken) B (662) N (829)
45
Metadynamics to accelerate the simulation of the healing process Coordination B def N def
Time [ps] 2.5 2 1.5 1 6 5 4 3 2 1
0.4 0.8 Partial Charge 6 5 4 3 2 1
46
Spacer Iter-pair Intra-pair Defect in the Chain [001]
Generalized relative displacement Good conductivity in doped samples
bond patterns
G.R. Gowers, J. Phys. Chem. B, 106, 9322 (2002)
Coordination numbers and/or rotation angles
D[001]
H,N ({RI}) =
1 NH
di · ˆ v[001] − 1 NN
dj · ˆ v[001]
a b
DFTPT in CPMD
H-Doped Optimized Structure
Fluctuating hydrogen-bond network: structure diffusion Mobility induced by excess protons Monitoring defect diffusion by activating 5 Collective Variables Facile intermolecular proton hopping Proton transferred over 4 molecules Complex pathway with several intermediates Rearrangement of the h-bond network Stabilizing role of O...H interactions
SCN SD SCN
1 1.2 1.4 1.8 1.6
FES(kcal/mol)
10
E F B C SD A
M Iannuzzi, Parrinello, PRL, 93, 025901(2004);
48
49
Three CN as CVs : N vs O, O vs C, N vs C
The first event is OH dissociation and hydroxilation of the surface (~35 kcal/mol), with formation of bi- radical system. System stabilized by the formation of H-bond and delocalization of the unpaired electron transferred to the surface. sp3 hybridization of C’ hydroxilation
H transfer reactive NO2 system is slightly stabilized; closed shell, only 1.9 kcal/mol more stable than
Epoxide structure migrates
MEP energy profile TS Epoxide migration TS back transfer The direct epoxidation is a much less probable process O...H 1.93 A
50 Early interaction between the incipient radical and the surface. N-OH 2.1 A OH-graphite 2.5 A interaction with π el. As soon as the OH radical is adsorbed on the surface, its unpaired electron is delocalized among the atoms of the graphite layer. MD snapshot during the transition MD snapshot at dissociation completed
The OH-graphite interaction starts before the dissociation is complete The dissociation barrier is 12 kcal/mol lower than in gas phase
if the partial pressure is such as to guarantee a high probability of having HNO3 in the vicinity of the surface.
51
Simplest type of defects, predominant in carbon layers σ type sp2 dangling bonds at the unsaturated C’ 4 CN as CVs: N->O; O->C; N->C;H->C
H transfer Oxidation Full saturation
Most stable
unpaired H el. saturates σ type dangling, π spin density delocalized C=O, unpaired π el. engaged no π spin density Ether group
Important effects of translational entropy on the free energy barriers
52
Hydroxilation Direct oxidation Replacement The estimated activation energies are much smaller than what
pristine graphite
MEP: ABC
C-O σ type dangling, π spin density delocalized weak hydrogen bond (1.88A); closed shell sp3 character, not symmetric, electropositive O, no H-bond the direct oxidation leads to more stable configurations but requires higher energy barriers E (kcal/mol)