Docstoc

Monte-Carlo simulations of the structure of complex liquids with

Document Sample
Monte-Carlo simulations of the structure of complex liquids with Powered By Docstoc
					 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
                                                            complex liquid
simple liquid          What are complex liquids?
                 hard spheres,
                Lennard-Jones
                   particles-
                   SIMPLE
                 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

                                         S
    Angular correlations
    completely ignored!!
                                 ρ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                        quantum mechanical vibrational and
external rotational degrees of                           internal rotational degrees of freedom
freedom                                 H n E n
two sets of independent quantum states (i.e. eigenstates can be taken as a product)
                 n  ncl nqu   with energy En=Encl+Enqu

                              The partition function factorises
                                        Q = Q cl Q qu

                        En    iqu
                         qu
                                            and       Qqu  ( exp( iqu / kT )) N
                                    i                          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 )  exp( H (r N p N  N p )) / Z
                              N                            N



              Z   dr N dp N d  N dp exp(   H (r N p N  N p ))
                                       N                          N



    The classical Hamiltonian can be split into kinetic and potential energy
                        H=Kt+Kr+U(rNωN)                                Iα
                 N                              N
          Kt   p / 2m
                                                                                     I xx           0
                                         Kr           
                       2                                                                      0
                       i                                            J i2 / 2 I
                                                                                   
                                                                                  I 0      I yy    0
                                                                                                         
                i 1                           i 1   x , y , z                    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 )  P(p N ) P( p ) P(r N  N )
                                  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 ) Jac 
                                                             N       N


                                   N      ( p )
                                                N
                             Jac                    Jac1 Jac2 ...JacN
                                            (J N )
                 Statistical mechanics of complex liquids
                                      new probability density

      Jac1  Jac 
                        ( p p p )
                                                sin 
                                                                   dr N dp N d' N dJ N P (r N p N  N J N )  1
                       ( J x J y J z )

                                                            P(r N p N  N J N )  P(r N p N  N p )(sin 1...sin  N )
                                                                                                  N




P(r N p N  N J N )  P(r N p N  N p ) Jac 
                                      N       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

                                 dr          dp N d N dp P(r N p N  N J N )  1
                                          N               N



                                                  We can write
                    P(r N p N  N J N )  exp( H (r N p N  N J N )) / Z
                     Z   dr N dp N d N dJ N exp(   H (r N p N  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
                                                                        Zt   dp exp(  pi2 / 2m)
                               N
        P(p )  exp(   p / 2m) / Z t
            N                         2
                                      i
                                                                                    N

                            i 1                                                                 i 1
                                                                                             N
                                                                        Z r   dJ N exp(          
                           N
        P(J N )  exp(    J i2 / 2 I ) / Z r
                                                                                                                   J i2 / 2I )
                                                                                                                       
                          i 1   x , y , z                                                 i 1       x, y, z

        P(r N  N )  exp(U (r N  N )) / Z c                         Z c   dr N d N exp(   U (r N  N ))
                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              Λt=h/(2πmkT)1/2
                                                                                           Λr=(h/(8π2IxkT)1/2)×
         1                                                 ZC
   N !3 N  r  N 
Q                   dr N d N exp(U (r N  N ))                                        (h/(8π2IykT)1/2)(h/(8π2IzkT)1/2)
       t
             N
                                                     N !t3 N  r  N
                                                                N

                                                                                            Ω=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   
                                           2  r          r12    r12   
                                         1
                                              12                         
             electrostatic interaction                                        dispersion polarisation
                                                     exchange repulsion
                                                     (finite size of atom)
                          q2                      electrostatic interaction between two molceules=
              z                · r2
      q1                                            interaction between two charge distributions

   r1·             r                               spherical harmonic expansion of r12-1=|r+r2-r1|-1
                                                    1
                                                         Alm1m2 m (r1r2 r )Yl1m1 (1 )Yl2 m2 (2 )Ylm ( )*
                                                                     ll
                               y                   r12 l1l2l m1m2 m 1 2

                                                            potential of a charge distribution:
  x                                                       qi              4  l l 1
                                                   i   r  ri
                                                                qi  
                                                                i    l m 
                                                                                   (ri / r )Ylm (i )Ylm ( ) *
                                                                           2l  1 
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)
                 u12   Al1l2 / r
                               l 1
                                          C (l l l; m m m)Q
                                                      1 2    1    2       l1m1   Ql2 m2 Ylm ( ) *
                       l1l2          m1m2 m


                                Introduction of body-fixed coordinate frame:

        z                                     z
                                                                          z’            y’                  z’
                                              z’                 y’
             r                                                                               Ω
                                                                      y                          x’
                                                                               x’                                    y’
                       y                                x’
                                     x                                         Ylm ' ( ')   m Dmm ' ()Ylm ( )
                                                                                                  l



 x
Spherical multipole expansion of the electrostatic interaction energy

                z                                      z’’
                                    q2                           Relation between multipole moments in the
                                   x’’      · r2               space-fixed and body-fixed coordinate frames:
     q1
z’            y’                                   y’’                                 Qlm   Dmk () *Qlk
                                                                                                l

     r1·                r                                                                         n


                                                                 u12   Al1l2 / r l 1     C (l l l; m m m)Q                    Ql2 m2 Ylm ( ) *
                                                   y                                                   1 2     1   2       l1m1
     x’                                                                 l1l2              m1m2 m

                   xyz: space-fixed
           x’y’z’ and x’’y’’z’’: body-fixed
 x
                   u12   Al1l2  Ql1k1 Ql2 k2 / r l 1       C (l l l; m m m) D
                                                                      1 2      1   2
                                                                                        l1
                                                                                        m1k1   (1 ) * Dm2 k2 (2 ) * Ylm ( ) *
                                                                                                        l2

                            l1l2     k1k2                    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                      Spherical multipole expansion
                      r12 
                               r12   i j rij
                                       
                 17 terms /pair                                                    5 (10) terms /pair
Monte-Carlo simulations of ensembles of anisotropic particles
                        What do we do in a MC calculation?
               x2                                x2

          F    f ( x)dx
               x1
                                           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 ) / Zc 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)                equilibrium distribution
            ρ(2)  ρ(1)π                                          ρlim  lim ρ(1) π N
                                                                         N 
             ρ(3)  ρ(1)ππ
                                                              mn       if n  m (U n  U m )
                                                 mn  
      Microscopic reversibility                         mn ( n / m ) if n  m (U n  U m )
         (detailed balance):                     mm  1    mn
                               Metropolis:
           mn m   nm n                                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:                              1           reject
                                                                               ×ζ 1
 - allways accept if Unew ≤ Uold        allways accept

 - if Unew>Uold choose a random                                      accept × ζ
                                                                                2
 number ζ from the interval [0,1]                        0                    ΔUnm        Unew-Uold
         How to generate trial moves?       How many particles should be moved?
      translation                              sampling efficiency:  ri 2 / CPU time
 xi  x1    0,5   rmax
                                                U ~kT                  reasonable acceptance
 yi  y1    0,5   rmax
                                                        U           1  2U
 zi  z1    0,5   rmax                      U         ri 
                                                                  a
                                                                                  ri a rib  ...
                                                        ri a
                                                                     2 ri ri
                                                                          a    b


        rotation                                              0  f (U ab ) ri 2 O ( 4 )
 '    (  0,5)max
                                          ri 2  kT / f (U ab )
cos '  cos   (  0,5) cos max
                                        1. N particles, one at a time: CPU time ~ nN
 '    (  0,5)max
                                         ri 2  NkT / f (U ab )           ri 2 /(CPU )  kT / nf U ab 

    sampling efficiency down by a       2. N particles in one move: CPU time ~ nN
             factor 1/N                   ri 2  kT / f (U ab )         ri 2 /(CPU )  kT / Nnf U ab 
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
                                         dr       d N f (r N  N )  N !
                                               N




                 we further introduce a reduced generic distribution function:
                 N!                                                        1
f (r h h )             dr N  h d N  h exp(U (r N  N )) / Zc             dr d f (r  )
                                                                                     N h N h N N

              ( N  h)!                                                ( N  h)!
                                                  ideal gas:
                                f (r h h )  f (r11 ) f (r22 )... f (rhh )
                                    homogenous isotropic fluid:
                                                                                     f (r11 )   / 
                      dr d 
                         1   1   f (r11 )  f (r11 )  dr1d1  V (r11 )  N
                        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:
                 N ( N  1)2
g (r1212 ) 
                            2        dr ...dr
                                           3          N   d3 ...dN exp(   U (r N  N )) / Z c

               2
           
                 2
                        (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 ) * Dm2 n2 (2 ) * Ylm ( ) *
                                                                                                                             l2

                                 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
                                                                                                                          2 
                                                                                                                      2m
                                                                           l1   l2   m
        How to present the results of MC simulations?

                  g  r12    g (l1l2 m; r )Yl1m 1  Yl
                                                                              2 
                                                                          2m
                                 l1   l2   m

                           removing the m dependence:
                                      2l  1 
                                                 1/ 2

                  g (l1l2 m; r )                    C (l1l2l ; mm0) g (l1l2l ; r )
                                   l  4 



                                                                                             θ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 Lennard-
Jones 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

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:3
posted:11/25/2011
language:English
pages:18