Molecular Dynamics Simulation - Fundamentals

Molecular Simulation

1. 1
CE 530 Molecular Simulation
Lecture 18
Free-energy calculations
David A. Kofke
Department of Chemical Engineering
SUNY Buffalo
[email protected]

2. 2
Free-Energy Calculations
Uses of free energy
•Phase equilibria
•Reaction equilibria
•Solvation
•Stability
•Kinetics
Calculation methods
•Free-energy perturbation
•Thermodynamic integration
•Parameter-hopping

3. 3
Ensemble Averages
Simple ensemble averages are of the form
To evaluate:
•sample points in phase space with probability p(G)
•at each point, evaluate M(G)
•simple average of all values gives <M>
Previous example
•mean square distance from origin in region R
•sample only points in R, average r
2
Principle applies to both MD and MC
( ) ( )M d Mp= G G Gò
( )2 2
( )
( )
s
d s
r d r
G
G G
= G G
ò
ò
1 inside R
0 outside R
s=
ì
í
î
G

4. 4
Ensemble Volumes
Entropy and free energy relate to the size of the ensemble
•e.g., S = k lnW(E,V,N)
No effective way to measure the size of the ensemble
•no phase-space function that gives size
of R while sampling only R
imagine being place repeatedly at random
points on an island
what could you measure at each point to
determine the size of the island?
Volume of ensemble is numerically
unwieldy
•e.g. for 100 hard spheres
r = 0.1, W = 5 × 10
133
r = 0.5, W = 3 × 10
7
r = 0.9, W = 5 × 10
-142
Shape of important region is very complex
•cannot apply methods that exploit some simple geometric picture
W = number of states of given E,V,N
G

5. 5
Reference Systems
All free-energy methods are based on calculation of
free-energy differences
Example
•volume of R can be measured as a fraction of the total volume
sample the reference system
keep an average of the fraction of time occupying target system
•what we get is the difference
Usefulness of free-energy difference
•it may be the quantity of interest anyway
•if reference is simple, its absolute free
energy can be evaluated analytically
e.g., ideal gas, harmonic crystal
G
( )
R
s
G
G
W
= G
W
( )ln /
R R
S S k
G G
- = W W

6. 6
HS chemical potential is an entropy difference
The energy is zero or infinity
•any change in N that does not cause overlap will be
change at constant U
To get entropy difference
•simulate a system of N+1 spheres, one non-interacting
“ghost”
•occasionally see if the ghost sphere overlaps another
•record the fraction of the time it does not overlap
Here is an applet demonstrating this calculation
[ ]
,
/
( , , 1) ( , , )
U V
S k
S U V N S U V N
N
bm
¶æ ö
=- » - + -ç ¸
¶è ø
3
/ 1
non-overlap
3 3
S k N
V
N
N
V V
e e f
N N
bmD - +
L
æ ö
W
ç ¸= = =
ç ¸WL L
è ø
N+1
N
×
×
Hard Sphere Chemical Potential

7. 7
Free-Energy Perturbation
Widom method is an example of a free-energy perturbation
(FEP) technique
FEP gives free-energy difference between two systems
•labeled 0, 1
Working equation
1 0 0
0
1 0
1
1 0
0
0
1
( )
( )
0
( ) 1
0
( )
0
)
U
A A
U U U
U
U U
U U
U
d e e
d e
Q
e
Q
d e
d e
e
e
d
b b
b
b
b
b
b
b
p
-
- -
- - -
-
- -
-
-
-
G
= =
G
G
G
= G (G
=
=
ò
ò
ò
ò
ò
G
0
1
Free-energy difference
is a ratio of partition
functions

8. 8
Free-Energy Perturbation
Widom method is an example of a free-energy perturbation
(FEP) technique
FEP gives free-energy difference between two systems
•labeled 0, 1
Working equation
1 0
1 0
1
1 0
1
0
0
0
0
( )
0
( )
( ) 1
0
( )
0
)
U U
U U
U
A A
U
U
U
U U
d e
Q
e
Q d e
d e e
d e
d e
e
b
b
b
b b
b
b
b
p
- -
-
- -
-
- -
-
-
-
-
G
= =
G
G
=
G
G (= G
=
ò
ò
ò
ò
ò
G
0
1
Add and subtract
reference-system
energy

9. 9
Free-Energy Perturbation
Widom method is an example of a free-energy perturbation
(FEP) technique
FEP gives free-energy difference between two systems
•labeled 0, 1
Working equation
1
1 0
0
1 0 0
0
1
1
0
0
( ) 1
0
( )
( )
0
( )
0
)
U
A A
U
U U U
U
U
U U
U
d e
Q
e
Q d e
d e e
e
e
d
d e
b
b
b
b b
b
b
b
p
-
- -
-
- - -
-
- -
- -
G
= =
G
G
=
G
= G (G
=
ò
ò
ò
ò
ò
G
0
1
Identify
reference-system
probability
distribution

10. 10
Free-Energy Perturbation
Widom method is an example of a free-energy perturbation
(FEP) technique
FEP gives free-energy difference between two systems
•labeled 0, 1
Working equation
Sample the region important to 0 system, measure properties
of 1 system
1
1 0
0
1 0 0
0
1 0
1 0
( ) 1
0
( )
( )
0
( )
0
)
U
A A
U
U U U
U
U U
U U
d e
Q
e
Q d e
d e e
d e
d e
e
b
b
b
b b
b
b
b
p
-
- -
-
- - -
-
- -
- -
G
= =
G
G
=
G
= G (G
=
ò
ò
ò
ò
ò
G
0
1
Write as
reference-system
ensemble average

11. 11
Chemical potential
For chemical potential, U1 - U0 is the energy of turning on the
ghost particle
•call this ut, the “test-particle” energy
•test-particle position may be selected at random in simulation volume
•for hard spheres, e
-bu
t is 0 for overlap, 1 otherwise
then (as before) average is the fraction of configurations with no overlap
This is known as Widom’s insertion method
1 0
3
( )
0
t
A A
uV
N
e e
e
b bm
b
- - -
-
L
=
=

12. 13
Deletion Method
The FEP formula may be used also with the roles of the
reference and target system reversed
•sample the 1 system, evaluate properties of 0 system
Consider application to hard spheres
•e
bu
t is infinity for overlap, 0 otherwise
•but overlaps are never sampled
•true average is product of 0 × ∞
technically, formula is correct
•in practice simulation average is always zero
method is flawed in application
many times the flaw with deletion is not as obvious as this
1 0 1 0( ) ( )
0
A A U U
e e
b b- - - -
=
1 0 1 0( ) ( )
1
A A U U
e e
b b+ - + -
=
G
0
1
Original: 0  1 Modified: 1  0
1
tuN
V
e e
bbm ++
=

13. 14
Other Types of Perturbation
Many types of free-energy differences can be computed
Thermodynamic state
•temperature, density, mixture composition
Hamiltonian
•for a single molecule or for entire system
•e.g., evaluate free energy difference for hard spheres with and
without electrostatic dipole moment
Configuration
•distance/orientation between two solutes
•e.g, protein and ligand
Order parameter identifying phases
•order parameter is a quantity that can be used to identify the
thermodynamic phase a system is in
•e.g, crystal structure, orientational order, magnetization

14. 15
General Numerical Problems
Sampling problems limit range of FEP calculations
Target system configurations must be encountered when
sampling reference system
Two types of problem arise
•first situation is more common
deletion FEP provides an avoidable example of the latter
G
0
1
target-system space very small
G
0
1
target system outside of reference

15. 16
Staging Methods
Multistage FEP can be used to remedy the sampling problem
•define a potential Uw intermediate between 0 and 1 systems
•evaluate total free-energy difference as
Each stage may be sampled in
either direction
•yielding four staging schemes
•choose to avoid deletion calculation
1 0 1 0
( ) ( )
w w
A A A A A A- = - + -
G
0
1
W
0 1 Umbrella sampling
0 1 Bennett's method
0 1 Staged deletion
0 1 Staged insertion
W
W
W
W
Ø ®
® Ø
Ø Ø
® ®
Use staged insertion Use umbrella sampling
G
0
1
W
G
0
1
Use Bennett’s method
W

16. 17
Example of Staging Method
Hard-sphere chemical potential
Use small-diameter sphere as intermediate
( ) ( )
w N tA A u
N
e e
b b s- D - -
=
In first stage, measure fraction
of time random insertion of
small sphere finds no overlap
( )
1
( ) ( )( ) tN w
u uA A
w
e e
b s sb
+
- -- D -
=
In second stage, small sphere moves around with
others. Measure fraction of time no overlap is
found when it is grown to full-size sphere
( )
1
( ) ( )( ) ( ) tN N t
u uA A u
N w
e e e
b s sb b s
+
- -- D - -
=
G

17. 18
Multiple Stages
G
0
1
W
1
W
2
W
3
G
0
2
W
1
Multistage insertionMultistage umbrella samplingMultistage Bennett’s method
W
0
1
2
G
0
2
W
01
1
W
12
0W
01
1
1W
12
2
0W
1
W
2
W
3
1

18. 20
Overlap Sampling Optimization
When using overlap sampling, we need to
define our intermediate (W).
W is important only where both 0 and 1 are
important.
This for was first proposed by Bennett,
often called Bennett's acceptance ratio
(BAR).
Requires A to be known before
measuring data.
•Guess A
•Iterate to refine estimate
G
0
1
Bennett’s method
W
e
W
=
e
−βU
0
e
−βU
1
e
−βU
0
e
−β(A
1
−A
0
)
+e
−βU
0

19. 21
Umbrella Sampling Optimization
When using umbrella sampling, we need to
define a different intermediate (W).
W is important only where either 0 or 1 are
important
Also requires A, iterating requires
sampling new configurations.
Other intermediate systems might work
even better
•Include regions important to neither 0 nor 1
•While sampling W, simulation can hop between
0 and 1 more easily.
•Not easy to determine W
G
1
Umbrella sampling
e
−βU
W
=e
−βU
0
e
−β(A
1−A
0)
+e
−βU
1
G
0
1
W

20. 22
Thermodynamic Integration 1.
Thermodynamics gives formulas for variation of free energy
with state
These can be integrated to obtain a free-energy difference
•derivatives can be measured as normal ensemble averages
•this is usually how free energies are “measured” experimentally
( )d A Ud PdV dNb b b bm= - +
,, T NV N
A A
U P
V
b b
b
b
¶ ¶æ ö æ ö
= = -ç ¸ç ¸
¶ ¶ è øè ø
2
1
2 1
( ) ( ) ( )
V
V
A V A V P V dVb b= - ò

21. 23
Thermodynamic Integration 2.
TI can be extended to follow uncommon (or unphysical)
integration paths
•much like FEP, can be applied for any type of free-energy change
Formalism
•Let l be a parameter describing the path
•the potential energy is a function of l
•ensemble formula for the derivative
•then
3
3
( ; )1
!
( ; )1
!
ln 1
1
( ; )
N
N
N
N
N U r
N
N U r N
N
A Q
dr e
Q
dr e U r
Q
U
b l
b l
b
l l l
b l
l
b
l
-
L
-
L
¶ ¶ ¶ é ù
= - = -
ê úë û¶ ¶ ¶
¶é ù
= +
ê ú
¶ë û
¶
=
¶
ò
ò
2
1
2 1
( ) ( )
U
A A d
l
l
b
b l b l l
l
¶
= +
¶
ò

22. 24
Thermodynamic Integration Example
The soft-sphere pair potential is
given by
Softness and range varies with n
•large n limit leads to hard spheres
•small n leads to Coulombic
behavior
Thermodynamic integration can
be used to measure free energy as
a function of softness s = 1/n
•Integrand is
( )
n
u r
r
s
e
æ ö
=ç ¸
è ø
Exhibits simplifying behavior because
es
n
is the only potential parameter
2
( )ln( / )
A
u r r
s s
b be
s
¶
= -
¶
3.0
2.5
2.0
1.5
1.0
0.5
0.0
U/e
2.52.01.51.0
Separation, r/s
n = 20
n = 12
n = 6
n = 3

23. 25
Parameter Hopping. Theory
View free-energy parameter l as another dimension in phase
space
Partition function
Monte Carlo trials include changes in λ
Probability that system has λ = λ
0 or λ = λ
1
( , , )
N N
E E l=p r
G
0
l
1
l
0 1
( , , )
( , , ) ( , , )
0 1
N N
N N N N
N N E
E EN N N N
Q d d e
d d e d d e
Q Q
b l
l
b l b l
-
- -
=
= +
= +
åò ò
ò ò ò ò
p r
p r p r
p r
p r p r
00
0 1 0 1
11
0 1 0 1
( , , )1
0
( , , )1
1
( )
( )
N N
N N
QEN N
Q Q Q Q
QEN N
Q Q Q Q
d d e
d d e
b l
b l
p l
p l
-
+ +
-
+ +
= =
= =
ò ò
ò ò
p r
p r
p r
p r

24. 26
Parameter Hopping. Implementation
Monte Carlo simulation in which l-change trials are attempted
Accept trials as usual, with probability min[1,e
-bDU]
Record fractions f
0, f
1 of configurations spent in l = l
0 and l = l
1
Free energy is given by ratio
In practice, system may spend almost no time in one of the values
•Can apply weighting function w(l) to encourage it to sample both
•Accept trials with probability min[1,(w
n/w
o) e
-bDU]
•Free energy is
•Good choice for w has f
1 = f
0
Multivalue extension is particularly effective
•l takes on a continuum of values
1 0( ) 1 0 11 1
0 0 0 1 0
( )
( )
A A Q Q QQ f
e
Q Q Q Q f
b- - +
= = =
+
1 0( ) 0 1
1 0
A A w f
e
w f
b- -
=

25. 27
Summary
Free energy calculations are needed to model the most
interesting physical behaviors
•All useful methods are based on computing free-energy difference
Four general approaches
•Free-energy perturbation
•Thermodynamic integration
•Parameter hopping
FEP is asymmetric
•Deletion method is awful
Four approaches to basic multistaging
•Umbrella sampling, Bennett’s method, staged insertion/deletion