VIEWS: 3 PAGES: 18 POSTED ON: 11/25/2011 Public Domain
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 412 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 412 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 (r11 ) f (r22 )... f (rhh ) homogenous isotropic fluid: f (r11 ) / dr d 1 1 f (r11 ) f (r11 ) dr1d1 V (r11 ) N How to present the results of MC simulations? angular correlation function, g(rhωh) : f (r h h ) f (r11 ) f (r22 )... f (rhh ) g (r h h ) generally: f (r h h ) h g (r h h ) / h pair correlation function: N ( N 1)2 g (r1212 ) 2 dr ...dr 3 N d3 ...dN 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 (r12 ) 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 (r12 ) g (l1l2l; r )C (l1l2l ; m1m2 m)Yl1m1 (1 )Yl2 m2 (2 )Ylm ( ) * l1l2 l m1m2 m intermolecular frame ω=0φ : g r12 g (l1l2 m; r )Yl1m 1 Yl 2 2m l1 l2 m How to present the results of MC simulations? g r12 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