Transcript Παρουσίαση του PowerPoint

Simulations of Soft Matter under Equilibrium
and Non-equilibrium Conditions
VAGELIS HARMANDARIS
International Conference on Applied Mathematics
Heraklion, 16/09/2013
Outline
 Introduction: Motivation, Length-Time Scales, Simulation Methods.
 Multi-scale Particle Approaches: Atomistic and systematic coarse-grained
simulations of polymers.
 Application: Equilibrium polymeric systems.
 Application: Non-equilibrium (flowing) polymer melts.
 Conclusions – Open Questions.
COMPLEX SYSTEMS: TIME - LENGTH SCALES:
• A wide spread of characteristic times: (15 – 20 orders of magnitude!)
-- bond vibrations: ~ 10-15 sec
-- dihedral rotations: 10-12 sec
-- segmental relaxation: 10-9 - 10-12 sec
-- maximum relaxation time, τ1: ~ 1 sec (for Τ < Τm)
Modeling of Complex Systems: Molecular Dynamics
Classical mechanics: solve classical equations of motion in phase space, Γ:=Γ(r, p).
 In microcanonical (NVE) ensemble:
The evolution of system from time t=0 to time t is given by :
Liouville operator:
(t )  exp iLt  (0)
g 
 
iL  K , H    ri
 Fi


r

p
i 1 
i
i
ri
pi
ri :

t
mi
Hamiltonian (conserved quantity):
N
pi
U
pi :

 Fi
t
ri
H NVE
pi 2
 K  U r   
 U (r)
i 2mi
Modeling of Complex Systems: Molecular Dynamics
 Various methods for dynamical simulations in different ensembles.
 In canonical (NVT) ensemble:
-- Langevin (stochastic) Thermostat
-- Nose-Hoover thermostat: [Nosé 1984; Hoover, 1985]: add one more degree of freedom ζ.
H NVT
p 2
pi 2
 H NVE  Ks  Vs  
 V (r) 
 3NkBT ln 
Q
i 2mi
pi
ri 
mi
i 
p
Q
U p
pi  

pi
ri
Q
p  
i
p 2i
 3NkBT
mi
Molecular Interaction Potential (Force Field)
Molecular model: Information for the functions describing the molecular interactions
between atoms.
U
(
r
)

V

V

V

V

V

V

V
b
o
n
d
e
dn
o
n

b
o
n
d
e
d
s
t
r
b
e
n
dt
o
r
s
L
J
q
-- stretching potential
-- bending potential
-- dihedral potential
1 (
2
V
lo)
s
trk
s
trl
2
1
2
V

k
(



)
o
b
e
n
d 2
b
e
n
d
5
i()
V
c
o
s


to
r
s
ic
i
0
Van der Waals (LJ)
-- non-bonded potential
 1

2
6

 


V
4





L
J
r
 
 


 r




Coulomb
qq
i j
Vq 
εrij
-- Potential parameters are obtained from more detailed simulations or fitting to
experimental data.
MULTI-SCALE DYNAMIC MODELING OF COMPLEX SYSTEMS
Atomistic MD Simulations: Quantitative predictions of the dynamics in soft matter.
Limits of Atomistic MD Simulations (with usual computer power):
-- Length scale:
-- Time scale:
few Å – O(10 nm)
few fs - O(1 μs)
(10-15 – 10-6 sec)
~ 107 – 109 time steps
-- Molecular Length scale (concerning the global dynamics):
up to a few Me for “simple” polymers like PE, PB
much below Me for more complicated polymers (like PS)
Need:
- Simulations in larger length – time scales.
- Application in molecular weights relevant to polymer processing.
- Quantitative predictions.
Proposed method:
- Coarse-grained particle models obtained directly from the chemistry.
Systematic Coarse-Graining: Overall Procedure
1. Choice of the proper CG description.
-- Microscopic (N particles)
-- Mesoscopic (M “super particles”)
Q : q1, q2 ,..., qM 
R : r1 , r2 ,..., rN 
Q : T R
-- Usually T is a linear operator (number of particles that correspond to a ‘super-particle’
N
qi   ciri ,
j 1
i  1, 2,..., M
Systematic Coarse-Graining: Overall Procedure
2. Perform microscopic (atomistic) simulations of short chains (oligomers) (in vacuum)
for short times.
3. Develop the effective CG force field using the atomistic data-configurations.
4. CG simulations (MD or MC) with the new coarse-grained model.
Re-introduction (back-mapping) of the atomistic detail if needed.
Effective (Mesoscopic) CG Interaction Potential (Force Field)
CG Potential: In principle UCG is a function of all CG degrees of freedom in the system and
of temperature (free energy):
 CG Hamiltonian – Renormalization Group Map:
e  U
CG
( q1 ,q 2 ,...,q M ,T )

AT

....
exp


U
 r1 , r2 ,..., rN  dr1 ,...rN | q1 , q 2 ,..., q M
  
ZN
N
qi   ci ri
 Remember:
j 1
 Assumption 1:
CG
CG
U CG (Q, T )  Ubonded
(Q, T )  Unon
bonded (Q, T )
: P CG  Q, T 
Bonded CG Interaction Potential
Bonded Potential
 Degrees of freedom: bond lengths (r), bond angles (θ), dihedral angles ()
CG
CG
Ubonded
(Q, T )  Ubonded
(r, , , T )
r
Procedure:
 From the microscopic simulations we calculate the distribution functions of the
degrees of freedom in the mesoscopic representation, PCG(r,θ,).
 PCG(r,θ, ) follow a Boltzmann distribution:
 Assumption 2:
 Finally:
CG

U
(r ,  ,  ) 
PCG  r , ,    exp   bonded

kT


PCG  r, ,   PCG  r  PCG   PCG  
CG
Ubonded
( x, T )  kBT ln PCG  x, T  ,
(x  r, , )
Non-bonded CG Interaction Potential: Reversible Work
 Assumption 3: Pair-wise additivity
M 1 M
U
CG
non bonded
CG
(Q, T )  U nb
 qi , q j , T 
q
i 1 j i
Reversible work method [McCoy and Curro, Macromolecules, 31, 9362 (1998)]
 By calculating the reversible work (potential of mean force) between the centers
of mass of two isolated molecules as a function of distance:
e  U
CG
nb
( q1 ,q 2 ,T )

AT

....
exp


U
 r1 , r2 ,..., rN   dr1 ,...rN  | q1, q 2
  
ZN
U nbCG (q, T )   ln exp  U AT  r,   
U AT  r,    U AT rij 
i, j
 Average < > over all degrees of freedom Γ that are integrated out (here orientational )
keeping the two center-of-masses fixed at distance r.
CG MD DEVELOPMENT OF CG MODELS DIRECTLY FROM THE CHEMISTRY
APPLICATION: POLYSTYRENE (PS)
[Harmandaris, et al. Macromolecules, 39, 6708 (2006); Macromol. Chem. Phys. 208, 2109 (2007);
Macromolecules 40, 7026 (2007); Fritz et al. Macromolecules 42, 7579 (2009)]
1) CHOICE OF THE PROPER COARSE-GRAINED MODEL
 2:1 model: Each chemical repeat unit replaced by two CG spherical beads (PS: 16
atoms or 8 “united atoms” replaced by 2 beads).
 CG operator T: from “CHx” to “A” and “B”
description.
 Each CG bead corresponds to O(10) atoms.
σΑ = 4.25 Å
σB = 4.80 Å
 Chain tacticity is described through the effective bonded potentials.
 Relatively easy to re-introduce atomistic detail if needed.
2) ATOMISTIC SIMULATIONS OF ISOLATED PS RANDOM WALKS
CG MD Simulations: Structure in the Atomistic Level after Re-introducing the
Atomistic Detail in CG Configurations.
 Simulation data: atomistic configurations of polystyrene obtained by reinserting
atomistic detail in the CG ones.
 Wide-angle X-Ray diffraction measurements [Londono et al., J. Polym. Sc. B, 1996.]
grem: total g(r) excluding
correlations between first
and second neighbors.
CG Polymer Dynamics is Faster than the Real Dynamics
PS, 1kDa, T=463K
ATCG
Free Energy Landscape
-- CG effective interactions are softer than the
real-atomistic ones due to lost degrees of
freedom (lost forces).
Free energy
CG
Atomistic
This results into a smoother energy landscape.
 CG MD: We do not include friction forces.
Configuration
CG Polymer Dynamics – Quantitative Predictions
CG dynamics is faster than the real dynamics.
Time Mapping (semi-empirical) method:
 Find the proper time in the CG description by moving the raw data in time. Choose a
reference system. Scaling parameter, τx, corresponds to the ratio between the two
friction coefficients.
Time Mapping using the mean-square displacement of the chain center of mass
Time Scaling

Check transferability of τx for different systems, conditions (ρ, T, P, …).
Polymer Melts through CG MD Simulations: Self Diffusion Coefficient
 Correct raw CG diffusion data using a time mapping approach.
[V. Harmandaris and K. Kremer, Soft Matter, 5, 3920 (2009)]
 Crossover regime: from Rouse to reptation dynamics. Include the chain end (free
volume) effect.
-- Rouse: D ~ M-1
-- Reptation: D ~ M-2
Crossover region:
-- CG MD: Me ~ 28.000-33.000 gr/mol
-- Exp.: Me ~ 30.0000-35.000 gr/mol
-- Exp. Data: NMR [Sillescu et al. Makromol. Chem., 188, 2317 (1987)]
CG Simulations – Application: Non-Equilibrium Polymer Melts
 Non-equilibrium molecular dynamics (NEMD): modeling of systems out of
equilibrium - flowing conditions.
 NEMD: Equations of motion (pSLLOD)
q i
p i :
t
simple shear flow
pi
 Fi  pi  u  mi qi  u  u
t
0

u
0

0

0
0
0
0

0

 In canonical ensemble (Nose-Hoover) [C. Baig et al., J. Chem. Phys., 122, 11403, 2005] :
p
pi
 Fi  pi  u  mi qi  u  u   pi
t
Q


t
Q
p
p
p 2i

 3NkBT
t
i mi
CG Simulations – Application: Non-Equilibrium Polymer Melts
 NEMD: equations of motion are not enough: we need the proper periodic boundary
conditions.
 Steady shear flow:
Lees-Edwards Boundary Conditions
Lt
 
y
L
ux
y
ux
Primary
x
x
simple shear flow
0

u
0

0

0
0
0
0

0

 L
CG Polymer Simulations: Non-Equilibrium Systems
 CG NEMD - Remember: CG interaction potentials are calculated as potential of mean
force (they include entropy). In principle UCG(x,T) should be obtained at each state point, at
each flow field.
U CG (Q, T )  kBT ln PCG Q, T 
Important question: How well polymer systems under non-equilibrium (flowing)
conditions can be described by CG models developed at equilibrium?
Method:
[C. Baig and V. Harmandaris, Macromolecules, 43, 3156 (2010)]
Use of existing equilibrium CG polystyrene (PS) model.
 Direct comparison between atomistic and CG NEMD simulations for various flow
fields. Strength of flow (Weissenberg number, Wi = 0.3 - 200)
Wi  
 Study short atactic PS melts (M=2kDa, 20 monomers) by both atomistic and CG NEMD
simulations.
CG Non-Equilibrium Polymers: Conformations
 Properties as a function of strength of flow (Weissenberg number)
Wi  
 Conformation tensor
3 RR
 
c  2
R
eq
R
 Atomistic cxx: asymptotic behavior at high Wi because of (a) finite chain extensibility, (b)
chain rotation during shear flow.
 CG cxx: allows for larger maximum chain extension at high Wi because of the softer
interaction potentials.
CG Non-Equilibrium Polymers: Conformation Tensor
 cyy, czz: excellent agreement between atomistic and CG configurations.
CG Non-Equilibrium Polymers: Dynamics
 Is the time mapping factor similar for different flow fields?
[C. Baig and V. Harmandaris, Macromolecules, 43, 3156 (2010)]
Translational motion
 Purely convective contributions from the applied strain rate are excluded.
R
(
t
)

R
(
0
)


c
m
c
m
2
D

l
i
m
G
t


6
t
 Very good qualitative agreement between atomistic and CG (raw) data at low and
intermediate flow fields.
CG Non-Equilibrium Polymers: Dynamics
Orientational motion



t
2
R
(
t
)
R
(
0
)
R
(
t
)
e
x
p



e
q
r

 Rotational relaxation time: small variations at low strain rates, large decrease at high
flow fields.
 Good agreement between atomistic and CG at low and intermediate flow fields.
CG Non-Equilibrium Polymers: Dynamics
 Time mapping parameter as a function of the strength of flow.
 Strong flow fields: smaller time mapping parameter effective CG bead friction
decreases less than the atomistic one.
Reason: CG chains become more deformed than the atomistic ones.
Conclusions
 Hierarchical systematic CG models, developed from isolated atomistic chains, correctly
predict polymer structure and dimensions.
 Time mapping using dynamical information from atomistic description allow for
quantitative dynamical predictions from the CG simulations, for many cases.
 Overall speed up of the CG MD simulations, compared with the atomistic MD, is ~ 3-5
orders of magnitude.
 System at non-equilibrium conditions can be accurately studied by CG NEMD
simulations at low and medium flow fields.
 Deviations between atomistic and CG NEMD data at high flow fields due to softer CG
interaction potentials.
Challenges – Current Work
 Estimation of CG interaction potential (free energies): Check – improve all assumptions
Ongoing work with M. Katsoulakis, D. Tsagarogiannis, A. Tsourtis
 Quantitative prediction of dynamics based on statistical mechanics
e.g. Mori-Zwanzig formalism (Talk by Rafael Delgado-Buscalioni)
 Parameterizing CG models under non-equilibrium conditions
e.g. Information-theoretic tools (Talk by Petr Plechac)
 Application of the whole procedure in more complex systems
e.g. Multi-component biomolecular systems, hybrid polymer based nanocomposites
Ongoing work with A. Rissanou
ACKNOWLEDGMENTS
Prof. C. Baig
[School of Nano-Bioscience and Chemical Engineering, UNIST
University, Korea]
Funding:
ACMAC UOC [Regional Potential Grant FP7]
DFG [SPP 1369 “Interphases and Interfaces ”, Germany]
Graphene Research Center, FORTH [Greece]
EXTRA SLIDES
APPLICATION: PRIMITIVE PATHS OF LONG POLYSTYRENE MELTS
 Describe the systems in the levels of primitive paths
[V. Harmandaris and K. Kremer, Macromolecules, 42, 791, (2009)]
 Entanglement Analysis using the Primitive Path Analysis (PPA) method
[Evereaers et al., Science 2004, 303, 823].
CG PS configuration (50kDa)
PP PS configuration (50kDa)
 Calculate directly PP contour length Lpp,, tube diameter:
appNmon
N
e
R2(N)
app
R2(N)
Lpp
-- PP CG PS: Ne ~ 180 ± 20 monomers
CALCULATION OF Me in PS: Comparison Between Different Methods
 Several methods to calculate Me: broad spread of different estimates
[V. Harmandaris and K. Kremer, Macromolecules, 42, 791 (2009)]
Method
T(K)
Ne (mers)
Reference
4RT
0
Rheology GN 
5 Me
423
140 ± 15
Liu et al., Polymer, 47, 4461
(2006)
Self-diffusion coefficient
458
280-320
Antonieti et al., Makrom.
Chem., 188, 2317 (1984)
Self-diffusion coefficient
463
240-300
This work
Segmental dynamics
463
110 ± 30
This work
Entanglement analysis
463
180 ± 20
This work
Entanglement analysis
413
124
Spyriouni et al.,
Macromolecules, 40, 3876 (2007)
MESOSCOPIC BOND ANGLE POTENTIAL OF PS
Distribution function PCG(θ,T)
CG Bending potential UCG(θ,T)
CG Simulations – Applications: Equilibrium Polymer Melts
 Systems Studied: Atactic PS melts with molecular weight from 1kDa (10 monomers)
up to 50kDa (1kDa = 1000 gr/mol).
 NVT Ensemble.
 Langevin thermostat (T=463K).
 Periodic boundary conditions.
STATIC PROPERTIES : Radius of Gyration
N
RG
2


R
R

R

i
C
M


i

1
2
G
SMOOTHENING OF THE ENERGY LANDSCAPE
Qualitative prediction: due to lost degrees of freedom (lost forces) in the local
level 
Local friction coefficient in CG mesoscopic description is
smaller than in the microscopic-atomistic one
ATCG
kBT
N
kT
a2
B
Reptation: D

N
3R2
Rouse:
D
CG diffusion coefficient is larger than the atomistic one
DAT DCG
Time Mapping Parameter: Translational vs Orientational Dynamics