Transcript Monte-Carlo simulations of the structure of complex

Monte-Carlo simulations of the
structure of complex liquids with
various interaction potentials
Aljaž Godec
Advisers: prof. dr. Janko Jamnik and doc. dr. Franci Merzel
National Institute of Chemistry
Contents
1. Introduction
2. Statistical mechanics of complex liquids
3. Spherical multipole expansion of the electrostatic interaction energy
4. Monte-Carlo simulations of ensembles of anisotropic particles
5. How to present the results of MC simulations?
6. Conclusions and considerations for future work
Introduction
simple liquid
complex liquid
What are complex liquids?
hard spheres,
Lennard-Jones
particlesSIMPLE
ISOTROPIC
POTENTIALS
anisotropic particles,
Importance of complex liquids? COMPLEX POTENTIALS
ΔU = 0
Hydrophobic interactions ΔU > 0
≈ρvapour
ΔF= ΔU-TΔS
ρbulk
Introduction
hard sphere in LJ fluid
≈ρvapour
S
Angular correlations
completely ignored!!
S
ρbulk
Statistical mechanics of complex liquids
Assumption: separable Hamiltonian
(intermolecular interactions have no effect on the quantum states)
H=Hclass+Hquant
classical: centre of mass and the
external rotational degrees of
freedom
H n E n
quantum mechanical vibrational and
internal rotational degrees of freedom
two sets of independent quantum states (i.e. eigenstates can be taken as a product)
with energy En=Encl+Enqu
n  ncl nqu
The partition function factorises
Q = Q cl Q qu
Enqu    iqu
and
i
Q qu  ( exp( iqu / kT )) N
i
individual molecular energy
Consequence of the above assumption: the contributions of quantum
coordinates to physical properties are independent of density
Statistical mechanics of complex liquids
Probability density for the classical states
P(r N p N  N pN )  exp( H (r N p N  N pN )) / Z
Z   dr N dpN d N dpN exp( H (r N pN  N pN ))
The classical Hamiltonian can be split into kinetic and potential energy
Iα
H=Kt+Kr+U(rNωN)
N
N
Kt   p / 2m
i 1
Kr  
2
i

i 1   x , y , z
 I xx

I 0
 0

J i2 / 2 I
0
I yy
0
0 

0 
I zz 
In Monte-Carlo calculations we need only the configurational
probability density, but
J  f ( p )
P(r N p N  N pN )  P(p N ) P( pN ) P(r N  N )
we introduce a new distribution P'(rNpNωNJN)
P(r N p N  N J N )  P(r N p N  N pN ) Jac
Jac
N
 ( pN )

 Jac1 Jac2 ...JacN
N
 (J )
N
Statistical mechanics of complex liquids
new probability density
Jac1  Jac 
 ( p p p )
( J x J y J z )
  sin 
N
N
N
N
N N N N

d
r
d
p
d'

d
J
P
(
r
p  J ) 1

P(r N p N  N J N )  P(r N p N  N pN )(sin 1...sin  N )
P(r N p N  N J N )  P(r N p N  N pN ) Jac
N
it is convinient to introduce a new distribution P(rNpNωNJN)
P(r N p N  N J N )  P '(r N p N  N J N )(sin 1...sin  N )1
N
N
N
N
N N N N
d
r
d
p
d

d
p
P
(
r
p  J ) 1


We can write
P(r N p N  N J N )  exp( H (r N p N  N J N )) / Z
Z   dr N dpN d N dJ N exp( H (r N pN  N J N ))
Statistical mechanics of complex liquids
We can now directly factorize
P(r N p N  N J N )  P(p N ) P(J N ) P(r N  N )
Z  Zt Z r Z c
N
N
Zt   dp exp(  pi2 / 2m)
P(p )  exp(  p / 2m) / Zt
N
i 1
N
2
i
N
i 1
P(J N )  exp(   J i2 / 2I ) / Z r
Z r   dJ exp(  
P(r N  N )  exp(U (r N  N )) / Zc
Zc   dr N d N exp(U (r N  N ))
N
i 1   x , y , z
N
i 1


 x, y, z
J i2 / 2 I )
the ps and Js of different molecules are uncorrelated
P(p N )  P(p1 ) P(p2 )...P(p N )
P(J N )  P(J1 ) P(J 2 )...P(J N )
furthermore
P(p1 )  P( p1x ) P( p1 y )...P( p1z )
P(J1 )  P( J1x ) P( J1 y )...P( J1z )
thus we can directly integrate
ZC
1
N
N
N
N
Q
d
r
d

exp(


U
(
r

))

N ! 3t N  rN  N 
N ! 3t N  rN  N
Λt=h/(2πmkT)1/2
Λr=(h/(8π2IxkT)1/2)×
(h/(8π2IykT)1/2)(h/(8π2IzkT)1/2)
Ω=8π2 (4π)
Spherical multipole expansion of the electrostatic interaction energy
a molceule= a distribution of charges (placed in the atomic centres);
Atoms have finite sizes and also interact with polarization interactions
Pair potential energy:
qq
  12   6  
U12    1 2  412  12    12   
 r12 
1
2  r12
 r12   

electrostatic interaction
q2
z
q1
r1·
· r2
r
y
exchange repulsion
(finite size of atom)
dispersion polarisation
electrostatic interaction between two molceules=
interaction between two charge distributions
spherical harmonic expansion of r12-1=|r+r2-r1|-1
1
   Alm1l12ml 2 m (r1r2 r )Yl1m1 (1 )Yl2 m2 (2 )Ylm ( ) *
r12 l1l2 l m1m2 m
potential of a charge distribution:
x

i
qi
 4  l l 1
 qi  
 (ri / r )Ylm (i )Ylm ( ) *
r  ri
2
l

1


i
l
m
Spherical multipole expansion of the electrostatic interaction energy
mth component of the spherical multipole moment of order l :
Qlm   qi ril Ylm (i )
i
interaction between two charge distributions=
∑(interactions of components of multipole moments of charge
distributions)
l 1
u12   Al1l2 / r
l1l2
 C (l l l; m m m)Q
1 2
1
2
l1m1
Ql2 m2 Ylm ( )*
m1m2 m
Introduction of body-fixed coordinate frame:
z
z
z’
z’
y’
Ω
r
y
x’
y
x
x
z’
y’
x’
x’
l
Ylm ' ( ')   m Dmm
' ()Ylm ( )
y’
Spherical multipole expansion of the electrostatic interaction energy
z
z’
q1
r1·
y’
q2
x’’
z’’
· r2
Relation between multipole moments in the
space-fixed and body-fixed coordinate frames:
l
Qlm   Dmk
() *Qlk
y’’
r
y
x’
n
u12   Al1l2 / r l 1
l1l2
 C (l l l; m m m)Q
1 2
1
2
l1m1
Ql2 m2 Ylm ( )*
m1m2 m
xyz: space-fixed
x’y’z’ and x’’y’’z’’: body-fixed
x
u12   Al1l2  Ql1k1 Ql2 k2 / r l 1
l1l2
k1k2
 C (l l l; m m m) D
1 2
1
2
l1
m1k1
(1 ) * Dml22 k2 (2 ) * Ylm ( ) *
m1m2 m
TIP5P water model
calculated only once
What is gained?
example: molecule consisting of four charges
  12   6 
qq
U12  412  12    12     1 2
 r12   r12   i j rij
17 terms /pair
Spherical multipole expansion
5 (10) terms /pair
Monte-Carlo simulations of ensembles of anisotropic particles
What do we do in a MC calculation?
F
x2
 f ( x)dx
x1
x2
F   ( f ( x) / P( x)) P( x)dx
x1
P(x)... probability density function
Monte-Carlo: perform a number of trials τ: in each trial choose a
random number ζ from P(x) in the interval (x1,x2)
f (  )
F 
 trials
P (  )
How to choose P in a way, which allows the function evaluation to be
concentrated in the region of space that makes importatnt contributions to the
integral?
Construction of P(x) by Metropolis algorithm
P  exp( U ) / Z c generates a Markov chain of states
1. outcome of each trial depends only upon the preceding trial
2. each trial belongs to a finite set of possible outcomes
Monte-Carlo simulations of ensembles of anisotropic particles
a state of the system m is characterized by positions and orientations of all molecules
probability of moving from m to n = πmn
N possible states
πmn constitute a N×N matrix, π
each row of π sums to 1
probability that the system is in a particular state is given by the probability vector
ρ=(ρ1, ρ2, ρ3,..., ρm, ρn,..., ρN)
probability of the initial state = ρ(1)
ρ(2)  ρ(1) π
ρ(3)  ρ(1) ππ
Microscopic reversibility
(detailed balance):
Metropolis:
 mn  m   nm  n
equilibrium distribution
ρlim  lim ρ(1)πN
N 
 mn
if  n   m (U n  U m )

 mn  
 mn (  n /  m ) if  n   m (U n  U m )
 mm  1    mn
m n
n
 exp(  (U n  U m ))
m
Monte-Carlo simulations of ensembles of anisotropic particles
How to accept trial moves?
exp(-βΔU)
Metropolis:
- allways accept if Unew ≤ Uold
- if Unew>Uold choose a random
number ζ from the interval [0,1]
How to generate trial moves?
translation
xi  x1    0,5  rmax
yi  y1    0,5  rmax
zi  z1    0,5  rmax
rotation
 '    (  0,5)max
cos  '  cos   (  0,5) cos  max
 '    (  0,5) max
sampling efficiency down by a
factor 1/N
1
reject
allways accept
×ζ 1
accept × ζ
2
0
ΔUnm
Unew-Uold
How many particles should be moved?
sampling efficiency:  ri2 / CPU time
U ~kT
U 
reasonable acceptance
U
1  2U
a

r

 ri a rib  ...
i
a
a
b
2 ri ri
ri
 0  f (U ab ) ri 2 O ( 4 )
 ri2  kT / f (Uab )
1. N particles, one at a time: CPU time ~ nN
 ri2  NkT / f (Uab )
 ri2 /(CPU )  kT / nf Uab 
2. N particles in one move: CPU time ~ nN
 ri2  kT / f (Uab )
 ri2 /(CPU )  kT / Nnf Uab 
Monte-Carlo simulations of ensembles of anisotropic particles
How to represent results (especially angular correlations)?
we introduce a generic distribution function:
f (r N  N )  N ! P(r N  N )  N !exp(U (r N  N )) / Zc
N
N
N N
d
r
d

f
(
r
 )  N!

we further introduce a reduced generic distribution function:
N!
1
N h
N h
N
N
f (r h h ) 
d
r
d

exp(


U
(
r

))
/
Z

d r N  h d N  h f (r N  N )
c


( N  h)!
( N  h)!
ideal gas:
f (r h h )  f (r11 ) f (r22 )... f (rhh )
homogenous isotropic fluid:
 dr d f (r  )  f (r  ) dr d
1
1
1 1
1 1
1
1
 V (r11 )  N
f (r11 )   / 
How to present the results of MC simulations?
angular correlation function, g(rhωh) : f (r h h )  f (r11 ) f (r22 )... f (rhh ) g (r h h )
generally: f (r h h )   h g (r h h ) / h
pair correlation function:
g (r1212 ) 

N ( N  1)2

2
2
 dr ...dr d ...d
3
2
N
3
N
exp(U (r N  N )) / Z c
  (r  r ) (r  r ) (    ) (    )
i j
i
1
j
2
i
1
j
2
δ(ω)=δ(φ) δ(cosθ) δ(χ)
spherical harmonic expansion of the pair correlation function in a space fixed frame:
g (r12 )  
  g (l l l; n n ; r )C (l l l; m m m)D
1 2
1 2
1 2
1
2
l1
m1n1
(1 )* Dml22 n2 (2 )* Ylm ( )*
l1l2 l m1m2 m n1n2
linear molecules:
g (r12 )  

g (l1l2l ; r )C (l1l2l ; m1m2 m)Yl1m1 (1 )Yl2 m2 (2 )Ylm ( ) *
l1l2 l m1m2 m
intermolecular frame ω=0φ :
g  r12    g (l1l2 m; r )Yl1m 1  Yl
l1
l2
m
2m
2 
How to present the results of MC simulations?
g  r12    g (l1l2 m; r )Yl1m 1  Yl
l1
l2
2m
m
2 
removing the m dependence:
 2l  1 
g (l1l2 m; r )   

l  4 
1/ 2
C (l1l2 l ; mm0) g (l1l2l ; r )
θ2
EXAMPLE: dipoles in LJ spheres
r
reconstruction
φ
Conclusions and considerations for the future
- we have briefly reviewed the statistical mechanics of complex liquids
- in order to reduce the number of interaction terms that have to be
evaluated in each simulation step a spherical multipole expansion of the
electrostatic interaction energy was made
- the basics of the Monte-Carlo method for simulation of ensembles of
anisotropic particles were provided along with useful methods for
representing the results of such simulations.
- finally results of MC simulations of dipoles embedded in LennardJones spheres were briefly presented.
- employ such simulations to study biophysical processes, such as the
hydrophobic effect
- possibility of including polarization effects  basis for developing a
polarizable water model for biomolecular simulations