Transcript Document

Free Energy MD and Nanoscale Polymers
Lecture notes for Computational Nanotechnogoly and Molecular
Engineering Workshop, Pan American Advanced Study
Institutes
1/14/2004
Shiang-Tai Lin, Youyong Li, Seung Soon Jang, Prabal Maiti
Tahir Çağın, Mario Blanco, and William A. Goddard III
Materials and Process Simulation Center, Caltech
Free Energy Calculation in Nanotechnology
Free Energy is a key parameter in Nanofabrication
• Equilibrium structure is determined by free energy
– formation of self-assembled monolayers (SAM)
– nanoscale patterns in liquid crystals and block copolymers
– secondary structures of DNA, RNA
• However, “Many of the ideas that are crucial to the development of this area-"molecular shape", the interplay between enthalpy and entropy, the nature of noncovalent forces that connect the particles in self-assembled molecular aggregates-are simply not yet under the control of investigators.”
George M. Whitesides
(http://www.zyvex.com/nanotech/nano4/whitesidesAbstract.html)
NanoStructures
Functions
Free Energy
Weak interactions
(vdW, Coulomb, Hb)
Outline
• How to obtain Free Energies from MD Simulations
–Test Particle method
–Perturbation method
–Nonequilibrium method
–Normal mode analysis
• A new 2PT approach for efficient Free Energy Estimation
–Basic Ideas (with Blanco and Goddard)
–Test of method with LJ fluids (with Blanco and Goddard)
• Applications of 2PT methods in the study of Dendrimers
–Zimmerman H-bond dendrimer (with Jang, Çağın and Goddard)
–PAMAM dendrimer (with Maiti and Goddard)
–Percec dendrimer (with Li and Goddard)
Free Energy Calculation from Molecular Simulations
• Test Particle Method (insertion or deletion)
•Good for low density systems
•Available in the Sorption Module of Cerius2
• Perturbation Method (Thermodynamic integration, Thermodynamic perturbation)
•Applicable to most problems
•Require long simulations to maintain “reversibility”
• Nonequilibrium Method (Jarzinski’s equality)
•Obtaining differential equilibrium properties from irreversible processes
•Require multiple samplings to ensure good statistics
• Normal model Method
•Good for gas and solids
•Fast
•Not applicable for liquids
Reference: Frenkel, D.; Smit, B. Understanding Molecular Simulation from Algorithms to Applications. Academic press: Ed., New York, 2002.
McQuarrie, A. A. Statistical Mechanics. Harper & Row: Ed., New York, 1976.
Jarzynski, C. Nonequilibrium Equality for Free Energy Differences. Phys. Rev. Lett. 1997, 78, 2690.
From Normal Modes to Free Energy
• Total number of normal modes = 3N, N=number of particles
For an isolated molecule with N atoms
3 translation, 3 rotation (or 2 for linear molecules, 0 for monoatomic molecules), 3N -6 vibration
For a crystal with N particles
3 translation (acoustic modes), 3N-3 vibrational modes
• Each vibrational mode can be treated as a harmonic oscillator
1
2
 n  (n  )h
Q
qHO
( )   exp(β n ) 
n
exp(βh/2)
1 - exp(βh/2)
• The partition function is the sum of contributions from HOs

ln Q   dS ( ) ln q HO ( )
0
• All the thermodynamic properties are defined
  ln Q 
E  V0  Tβ 1 


T

 N ,V
1
A  V0  β ln Q  V0  β
  ln Q 
S  k ln Q  β 1 


T

 N ,V
1

 dS ( )W
0
A
( )
Determine Normal Modes from Molecular Simulations
•The density of states S()
•S()d number of modes between  and +d
•


0
S ( )d  3N
•Determination of S()
• The eigenvalues of the Hessian matrix (vibrational analysis, phonon spectrum )
• Covariance matrix of atomic position fluctuations
• Fourier transform of velocity autocorrelation function
S ( υ) 

2
lim  C (t )e i 2 πυt dt
kT   
C(t )  v(0)  v(t )
The Density of States Distribution S()
Solid
Gas
Liquid
Debey crystal
S(v) ~2
S ( )
S ( )
S ( )


exponential
decay

• Singularity at zero frequency
• Strong anharmonicity at low frequency regime
The 2PT idea: Liquid  Solid+Gas
solid-like
gas-like
S ( )

• Decompose liquid S(v) to a gas and a solid contribution
• S(0) attributed to gas phase diffusion
• Gas component contains anharmonic effects
• Solid component contains quantum effects
• Two-Phase Thermodynamics Model (2PT)
The 2PT Method
Lin, S. T.; Blanco, M.; Goddard, W. A. The Two-Phase Model for
Calculating Thermodynamic Properties of Liquids from Molecular
Dynamics: Validation for the Phase Diagram of Lennard-Jones Fluids. J.
Chem. Phys. 2003, 119, 11792.
• The basic idea
•The DoS
S ( )  S gas ( )  S solid ( )

•Thermodynamic properties
P   dS ( )W
s

HO
P
( )   dS g ( )WPg ( )
0
• The gas component
0
• VAC for hard sphere gas
c HS (t ) 
3kT
exp( a t )
m
a (T ,  HS , HS ) : frictioncoefficient
Gas
• DoS for hard sphere gas
S HS ( ) 
12N a

2
2 2
a  4 
gas
s0
 s  
1  0 
 6 fN 
S0
N gas  fN
2
s0  S
HS
(0) 
• Two unknowns (a and Ngas) or (so and f)
12N gas
S ( )
f
a

Determining so and f from MD Simulation
• so (DoS of the gas component at =0)
• completely remove S(0) of the fluid
• s0  S (0), S solid (0)  0
• f (gas component fraction)
• T or 0 : f1 (all gas)
• : f0 (all solid)
•
f 
D(T ,  )
D0HS (T ,  ; HS )
• one unknown HS
D(T ,  ) 
kTS (0)
12 mN
D0HS (T ,  ; HS ) 
3 1
kT 1/ 2
(
)
8  HS 2 m
(Chapman- Enskog)
Determining HS
• HS (hard sphere radius for describing the gas molecules)
• gas component diffusivity should agree with statistical
mechanical predictions at the same T and 
• gas component diffusivity from MD simulation
D HS (T , f ) 
kTs0
12m fN
• HS diffusivity from the Enskog theory
y
D (T , f )  D (T , f ;
HS
HS
0
HS
4 fy
)
z ( fy)  1

6

HS 3
1 y  y2  y3
z ( y) 
(1  y) 3
At Last…
• A universal equation for f
29 / 2 f 15 / 2  63 f 5  3 / 2 y 7 / 2  63 / 2 f 5 / 2  2 f  2  0
normalized diffusivit y : (T ,  , m, s0 ) 
2s0 kT 1/ 2 1/ 3 6 2 / 3
(
)  ( )
9N m

• Graphical representation
1.0
f or f y
0.8
f
fy
0.6
0.4
0.2
0.0
1.E-05
1.E-03
1.E-01
1.E+01
 (normalized diffusivity)
1.E+03
Comparison of the 1PT and 2PT methods
Run a MD simulation
(trajectory information saved)
Calculate VAC
Calculate DoS (FFT of VAC)
Apply HO
approximation
To S()
1PT
thermodynamic
predictions
Apply HO statistics
To Ssolid()
Apply HS statistics
to Sgas()
2PT
thermodynamic
predictions
Calculate S(0) and 
Solve for f
Determine Sgas(), Ssolid()
An Overview of the 2PT Method
5
0
To
Phase Behaviors
-5
G*
-10
-15
T*=1.8
T*=1.4
T*=1.1
T*=0.9
2PT(Q)
2PT(C)
-20
DoS
100

90
DoS S(cm)
50
G β
40
30


10
0

0

20
40
60
80
100
frequency v(cm-1)
400

2
S ( υ) 
lim  C (t )e i 2 πυt dt
kT   
0
-400
-800
0
MD
simulations
0.5
1
1.5
time (ps)
2
2.5
C(t )  v(0)  v(t )
3
0.4
0.6
*
1
0.8
1

 dS ( )W
0
20

800
0.2


60
VAC
1200
C (t)
0

70
1600
-30

80
From
MD Simulations
-25
A
1.2
( )
Test the 2PT Method Using the LJ System
• Intermolecular potential
Lennard-Jones Potential
 

V (r )  4 ( )12  ( ) 6 
r 
 r
V(r)
r=
0
-
r
T -  diagram for Lennard Jones Fluid
• Phase diagram
●stable
– critical point
1.8
Tc*  1.316  0.006
T  0.69
●unstable
1.4
 c*  0.304  0.006
– triple point
Gas
T*
Liquid
1.0
*
tp
(T*=kT/*=3)
●metastable
Supercritical Fluid
Solid
0.6
0.0
0.4
*
0.8
1.2
VAC and DoS of LJ Fluids
Velocity Autocorrelation
Density of States

2
S ( υ) 
lim  C (t )e i 2 πυt dt
kT   
C(t )  v(0)  v(t )
1600
100
800
C (t)
90

80

70
DoS S(cm)
1200



400
0
gas
liquid
solid
60





50
40
30
gas
liquid
solid
-400
-800
0
0.5
1
1.5
time (ps)
2
2.5
20
10
0
3
0
20
40
60
frequency v(cm-1)
80
100
2PT DoS Decomposition
• Examples
liquid
gas
1200
30
1000
25
solid
35
30
25
600
gas-like
solid-like
400
0
0
5
solid-like
gas-like
15
10
5
 [cm-1]
20
10
200
0
solid-like
gas-like
15
S ( ) [cm]
20
S ( ) [cm]
S ( ) [cm]
800
10
5
0
0
50
100
 [cm-1]
150
0
50
100
 [cm-1]
150
Pressure and Energy
Total Energy
Pressure
3
18
T*=1.8
16
T*=1.4
14
T*=1.1
1
12
T*=0.9
0
MD
10
P*
-1
E*
8
-2
6
-3
4
-4
2
-5
0
-6
-2
-7
0
0.2
0.4
0.6
*
0.8
T*=1.8
T*=1.4
T*=1.1
T*=0.9
MD
2PT(Q)
2
1
1.2
0
0.2
0.4
0.6
*
Pressures and MD Energies agree with EOS values
Quantum Effect (ZPE) most significant for crystals (~2%)
0.8
1
1.2
Entropy
2PT model
1PT
20
20
T*=1.8
T*=1.4
T*=1.1
T*=0.9
2PT(Q)
2PT(C)
T*=1.8
18
gas
16
T*=1.4
18
T*=1.1
16
T*=0.9
14
S*
14
1PT(Q)
S*
12
crystal
10
8
12
10
8
6
6
liquid
4
4
0
0.2
0.4
0.6
*
0.8
1
• Overestimate entropy for low
density gases
• Underestimate entropy for liquids
• Accurate for crystals
1.2
0
0.2
0.4
0.6
*
0.8
1
1.2
• Accurate for gas, liquid, and crystal
• Accurate in metastable regime
• Quantum Effects most important for
crystals (~1.5%)
Gibbs Free Energy
1PT
2PT model
5
5
liquid
0
0
-5
G*
-5
-10
crystal
G*
-15
-10
-15
T*=1.8
-20
T*=1.4
T*=1.1
-20
-25
T*=0.9
1PT(Q)
-25
gas
-30
T*=1.8
T*=1.4
T*=1.1
T*=0.9
2PT(Q)
2PT(C)
-30
0
0.2
0.4
0.6
*
0.8
1
• Underestimate free energy for
low density gases
• overestimate entropy for liquids
• Accurate for crystals
1.2
0
0.2
0.4
0.6
*
0.8
1
1.2
• Accurate for gas, liquid, and crystal
• Accurate in metastable regime
Why does 2PT work?

P   dS ( )W
s
1200

HO
P
( )   dS g ( )WPg ( )
0
HS fy = 0.036
0
HS fy = 0.309
4
gas
1000
WS( )
800
S( ) [cm]
5
600
3
2
400
QHO
1
200
CHO
0
0
0
2
4
6
 [cm-1]
8
10
0
50
100
150
  [cm ]
-1
30
liquid
25
S( ) [cm]
20
• 1PT overestimates Wsgas for gas for
modes < 5 cm-1
• 1PT underestimates Wsgas for liquid for
modes between 5 and 100 cm-1
15
10
5
• 2PT properly corrects these errors
0
0
50
100
 [cm-1]
150
Convergence of 2PT
15.5
14.5
• For gas, the entropy
13.5
2PT(Q)
2PT(C)
12.5
converges to within 0.2% with
MBWR EOS
11.5
S*
gas (*=0.05 T*=1.8)
10.5
• For liquid, the entropy
liquid (*=0.85 T*=0.9)
9.5
2500 MD steps (20 ps)
8.5
converges to within 1.5% with
7.5
2500 MD steps (20 ps).
6.5
100
1000
10000
MD steps
100000
1000000
2PT for Melting and Solidification
• Initial amorphous structure is used in
the cooling process
• The fluid remains amorphous in
simulation even down to T*=0.8
(supercooled)
• The predicted entropy for the fluid and
supercooled fluid agree well with EOS
for LJ fluids
metatstable supercritical
solid
unstable
fluid
Simulation conditions
2.0
supercritical
fluid
1.8
1.6
1.4
T*
1.2
1.0
0.8
solid
0.6
0.00
0.40
0.80
8
1.20
starting with
amorphous liquid

• Initial fcc crystal is used in the
heating process
• The crystal appears stable in
simulation even up to T*=1.8
(superheated)
• The predicted entropies for the
crystal and superheated crystal
agree well with EOS for LJ solids
Entropy
7
6
S*
liquid (EOS)
5
solid (EOS)
4
starting with
fcc crystal
heating
cooling
classical
3
0.80
1.20
1.60
T*
2.00
Extension to Mixtures
LJ Mixtures at T*=0.85
Combination Rules:
ij= lij ( ii+ jj)/2
ij= bij ( ii+ jj)1/2
22/11=0.80
22/11=0.85
l12=0.95
b12=0.70
2PT Literature
x1(I) ~0.02
0.04
x1(II) ~0.79
0.84
Efficiency of 2PT for Mixtures
100 ps
400 ps
T*=0.95 P*=0.125
2.0E-04
1.5E-04
1.5E-04
1.0E-04
400ps
100ps
5.0E-05
Gmix/RT
1.0E-04
5.0E-05
Gmix/RT
T*=0.95 P*=0.125
2.0E-04
0.0E+00
0.0E+00
-5.0E-05
-5.0E-05
-1.0E-04
-1.0E-04
-1.5E-04
-1.5E-04
-2.0E-04
-2.0E-04
-2.5E-04
0
-2.5E-04
0
0.2
0.4
0.6
0.8
0.2
1
0.4
0.6
0.8
1
x(Ar)
x(Ar)
25 ps
12.5 ps
T*=0.95 P*=0.125
T*=0.95 P*=0.125
2.0E-04
2.0E-04
1.5E-04
1.5E-04
25ps
1.0E-04
5.0E-05
Gmix/RT
Gmix/RT
12.5ps
1.0E-04
5.0E-05
0.0E+00
-5.0E-05
-1.0E-04
0.0E+00
-5.0E-05
-1.0E-04
-1.5E-04
-1.5E-04
-2.0E-04
-2.0E-04
-2.5E-04
0
0.2
0.4
0.6
x(Ar)
0.8
1
-2.5E-04
0
0.2
0.4
0.6
x(Ar)
0.8
1
An Overview of the 2PT Method
5
0
To
Phase Behaviors
-5
G*
-10
-15
T*=1.8
T*=1.4
T*=1.1
T*=0.9
2PT(Q)
2PT(C)
-20
DoS
100

90
DoS S(cm)
50
G β
40
30


10
0

0

20
40
60
80
100
frequency v(cm-1)
400

2
S ( υ) 
lim  C (t )e i 2 πυt dt
kT   
0
-400
-800
0
MD
simulations
0.5
1
1.5
time (ps)
2
2.5
C(t )  v(0)  v(t )
3
0.4
0.6
*
1
0.8
1

 dS ( )W
0
20

800
0.2


60
VAC
1200
C (t)
0

70
1600
-30

80
From
MD Simulations
-25
A
1.2
( )
Summary of 2PT
1.
2.
3.
4.
5.
6.
7.
Thermodynamic and transport properties are determined simultaneously.
Only short simulation times (20 ps) are needed to obtain high accuracy. For a
system with N particles, we expect 2PT to be N times faster than methods
such as particle insertion and thermodynamic integration.
The efficiency of 2PT does not deteriorate with increasing density (a severe
limitation in most other techniques).
The properties are obtained under real equilibrium conditions (no
perturbation in the simulation itself).
Zero point energy and corrections for quantum effects are included.
2PT can be used to determine the properties in metastable and unstable
regimes.
2PT could also be used for nonequilibrium systems to estimate effects of
transient effects, reaction, and phase transitions, since it is only necessary to
have stabilities over time scales of ~20 ps.
Application: Zimmerman H-bond Dendrimers
Dendrimer study on Zimmerman system
Initiative: Science, 271, 1095 (1996)
R=
H
O
O
g1
g2
g3
g4
O O
O OO O
O
O
O
O
O
O H
O
O
O O
O
O
O
O
O
R
O O
O
O
O
O
O
O O
H
N
O
O
O O
O
O O
OO O OO O
O
O
O
O
O
O
O
O
O
O O
O
O
O O
O
O
O
O H
O
dendrons
core
Experimental Observations
HO O
O
OH
• Aggregates found in non-polar solvents such as CH2Cl2
and CHCl3 but not in polar solvents such as THF and
DMSO
OH O HO
O
HO O
O
OH
OH O HO
O
• Circular and ladder coexist in lower generations (gen 1)
• Circular type is dominant for generations 2, 3 and 4
• The generation dependent stability behavior is attributed to
the subtle interplay between H-bond, vdW and steric
repulsions
O HO
OH O
O H O
OH O
O HO
OH O
O H O
OH O
OH O H
O O
HO O
O
OH
OH O H
O O
HO O
O
OH
circular
O
HO
OH
O
HO
O
H
O O
OH O
O
HO
OH
O
O H O
OH
O
O
H O
linear
O
O
H
H
O O
O O H
H O O
OH
O
H
O
O
Application: Zimmerman H-bond Dendrimers
• Relative stability determined by the free energy differences is in consistent
with experimental observations
• Linear and circular forms are energetically similar
• Stability of Zimmerman H-bond mediated dendrimers are dominated by
entropic effects
• Opens up possibility to design thermodynamically stable suprastructures
• This method is potentially useful for study other supramolecules such as
Protein, DNA, etc
Energy
Entropy
Free Energy
Application: PAMAM dendrimers
NH 2
NH 2
EDA core
NH
NH 2
O
repeating (monomer) unit
Generation 3
Applications of PAMAM dendrimer
• Catalysts
• Environment applications (metal encapsulation)
• Medical applications (drug delivery, gene therapy)
• Surface active agents
• Viscosity modifier
• New electrical materials
Application: PAMAM dendrimers
Water molecules in PAMAM
Free Water Domain
Surface Water Domain
Inner Water Domain
Schematic
Atomistic
• Number of water molecules
solution
overall
inner
surface
free
SESA (Å2)
high pH
9406
376
1014
1440
10642±185
neutral pH
low pH
8948
12977
547
747
1328
1729
1806
2472
13216±198
18526±723
Application: PAMAM dendrimers
high pH
neutral
low pH
-6
E (kcal/mol)
-7
• Engetically slightly less favored for
water at the PAMAM surface
• Entropically less favored for water
inside PAMAM
• Entropically slightly less favored for
water at PAMAM surface
• Surface and Inner Waters are in a
higher free energy state compared
to the bulk
-8
-9
-10
-11
-12
inner
surface
free
water location
high pH
neutral
low pH
6
4
3
2
-13
-14
-15
1
-16
0
-17
inner
surface
water location
free
neutral
low pH
-12
A (kcal/mol)
TS (kcal/mol)
5
high pH
-11
inner
surface
water location
free
Application: Percec Dendrimer Crystals
1500
Free energy profile over volume at 277K
1450
A(kJ/mol/dendron)
1400
1350
1300
1250
a. Condense
phase
1200
Critical pressure:
0.033GPa
1150
1100
b. Isolated Micelle
phase
1050
AA15
1000
0
2000
4000
6000
8000
10000
12000
Volume (A^3)
14000
16000
18000
20000