Transcript Computer Simulation of Colloids Jürgen Horbach

Computer Simulation of Colloids
Jürgen Horbach
Institut für Materialphysik im Weltraum, Deutsches Zentrum für Luftund Raumfahrt, Linder Höhe, 51147 Köln
 lecture 1
introduction
 lecture 2
introduction to the Molecular Dynamics (MD) simulation method
 lecture 3
case study: a glassforming Yukawa mixture under shear
 lecture 4
introduction to the Monte Carlo (MC) method
 lecture 5
case study: phase behavior of colloid-polymer mixtures studied by
grandcanonical MC
 lecture 6
modelling of hydrodynamic interactions using a hybrid MD-Lattice
Boltzmann method
Computer Simulations of Colloids
Introduction
(1) colloids as model systems: different interaction potentials and
hydrodynamic interactions
(2) correlation functions: describe the structure and dynamics of
colloidal fluids
(3) outline of the forthcoming lectures
(1) colloids as model systems: different interaction potentials and
hydrodynamic interactions
(2) correlation functions: describe the structure and dynamics of
colloidal fluids
(3) outline of the forthcoming lectures
colloids as model systems
size: 10nm – 1µm, dispersed in solvent of viscosity 
 large, slow, diffusive → optical access
 can be described classically
 tunable interactions, analytically tractable
→ theoretical guidance
 soft → shear melting
phenomena studied via colloidal systems:
glass transition, crystallisation in 2d and 3d, systems under shear,
sedimentation, phase transitions and dynamics in confinement, etc.
hard sphere colloids
 confirmed theoretically
predicted crystallization
(of pure entropic origin)
 study of nucleation and
crystal growth
 test of the mode coupling
theory of the glass
transition
tunable interactions I
hard spheres
volume fraction 
soft spheres
volume fraction, hair length and density
tunable interactions II
charged spheres
number density, salt concentration, charge
super-paramagnetic spheres
B
magnetic field
tunable interactions III
entropic attraction in colloid-polymer mixtures
volume fraction of colloidal particles,
number of polymers, colloid/polymer
size ratio
 colloid-polymer mixtures may exhibit a demixing transition which is
similar to liquid-gas transition in atomistic fluids
 transition is of purely entropic origin
hydrodynamic interactions
interaction between colloids due to the momentum transport through
the solvent (with viscosity η)

D0 
vi 
Fi
k BT

1
vi 
kBT
k BT
D0 
6a
(N ) 
(N )
 

 Dij (r )  Fj with r  (r1, r2 ,, rN )
N
j 1
 Dij: hydrodynamic diffusion tensor
 decays to leading order ~ r -1 (long-range
interactions!)
simulation of a colloidal system
 choose some effective interaction between the colloidal particles
 but what do we do with the solvent?
simulate the solvent explicitly
describe the solvent on a coarse-grained level
as a hydrodynamic medium
Brownian dynamics (ignores hydrodynamic
interactions)
ignore the solvent: not important for dynamic
properties (glass transition in hard spheres),
phase behavior only depends on configurational
degrees of freedom
(1) colloids as model systems: different interaction potentials and
hydrodynamic interactions
(2) correlation functions: describe the structure and dynamics of
colloidal fluids
(3) outline of the forthcoming lectures
pair correlation function
np(r) = number of pairs (i,j) with
rij  [r , r  r ]
ideal gas:
Δr small
i
np (r ) 
N 4
N

[( r  r )3  r 3 ]   4 r 2 r
2
3
2
r
fluid with structural correlations:
np (r ) 
g (r ) 
N
 4 r 2 r g (r )
2
average number density at distance r from an arbitrary atom in the origin
the same for an ideal gas at the same density 
pair correlation function: hard spheres vs. soft spheres
soft sphere fluid:
hard sphere fluid:
 (r  rij )
  
1
1
 g (r )    (r  ri  rj )  
N i j i
N i j i 4 r 2
static structure factor
N
N

  FT
 
microscopic
 (r )    (r  ri )   q   exp( iq  rj )
density variable:
i 1
j 1
homogeneous fluid:
 (r) 
N

V
static structure factor:

  
1
1 N N
sin( qr )
2


S (q) 
 q  q   exp[ iq  (rm  rn )]  1    4r ( g (r )  1)
dr
N
N m1 n1
qr
0
low q limit: S (q  0) 
N2  N
N
2
 k BT T
 T : isothermal compressib ility
large q:
S ( q  )  1
typical structure factor for a liquid
 1st peak at ~2π/rp where rp
measures the periodicity
in g(r)
 dense liquid → low
compressibility and therefore
small value of S(q→0)
 in general S(q) appropriate
quantity to extract information
about intermediate and large
length scales
 S(q) can be measured in
scattering experiments
van Hove correlation function
 generalization of the pair correlation function to a time-dependent
correlation function
 

1
G (r , t )    (r  [rk (t )  rl (0)])
N k l
proportional to the probability that a particle k at time t is separated by
a distance |rk(t) - rl(0)| from a particle l at time 0
 self part:
1
Gs (r , t ) 
N

k
 

 (r  [rk (t )  rk (0)])
intermediate scattering functions
coherent intermediate scattering function:
  
S (q, t )   G (r , t ) exp( iq  r )dr
1
1



 q (t )  q (0) 
N
N
S ( q, t )
F ( q, t ) 
S (q)

k
 

exp( iq  [rk (t )  rl (0)])
l
incoherent intermediate scattering function:
   1
 

Fs (q, t )   Gs (r , t ) exp( iq  r )dr   exp( iq  [rk (t )  rk (0)])
N k
can one relate these quantities to transport processes?
F(q,t) for a glassforming colloidal system
3.5
S(q)
1.0

 = 0.620
3.0

0.8
S(q,t)/S(q)
2.5
2.0
1.5
1.0

0.5
0.0
0.010
0.015
0.020
0.025
0.4
0.2
0.030

q [nm ]
 = 0.581
0.6
 = 0.513
0.0
-6
-5
-4
-3
-2
-1
0
1
2
3
10 10 10 10 10 10 10 10 10 10
Eckert, Bartsch, Faraday Discuss. 123, 51 (2003)
t/s
diffusion constants
self-diffusion
interdiffusion
tagged
atoms in I
 s
 Ds  s
t
concentration
differences
in I and II
↓
diffusion process
homogeneous distribution of
tagged particles in I+II
↓
cA
 DABcA
t
homogeneous distribution of
red and blue atoms in I+II
calculation of the self-diffusion constant
Einstein relation:
Ds  lim


[rtag (t )  rtag (0)]2
t 
2dt

Green-Kubo relation:


1
Ds   d v tag ( )  v tag (0)
30
vtag(t): velocity of tagged particle at time t
the two relations are strictly equivalent!
Computer Simulations of Colloids
Molecular Dynamics I:
How does it work?
Outline
MD – how does it work?
 numerical integration of the equations of motion
 periodic boundary conditions, neighbor lists
 MD in NVT ensemble: thermostatting of the system
 MD simulation of a 2d Lennard-Jones fluid:
structure and dynamics
Newton‘s equations of motion

 classical system of N particles at positions {ri }, i  1,..., N
U pot 


mi ri     Fi
ri
 Upot : potential function, describes interactions
between the particles
 simplest case: pairwise additive interaction between point particles
N 1 N
U pot   u (rij )
i 1 j i
 
rij  ri  rj
 solution of equations of motion yield trajectories of the particles, i.e.
positions and velocities of all the particles as a function of time
microcanonical ensemble
 consider closed system with periodic boundary conditions, no coupling
to external degrees of freedom (e.g. a heat bath, shear, electric fields)
 particle number N and volume V fixed → microcanonical ensemble
 momentum conserved, set initial conditions such that total momentum
is zero
 total energy conserved:
E  Ekin  U pot  const.
simple observables and ergodicity hypothesis
 potential energy:
 kinetic energy:
 pressure:
U
1
u 

N
N
ekin 
Ekin
N
N 1 N
 u(r )
1

N
Nk T
1
p B 
V
dV
i 1 j i
ij
mi vi2
 kBT

2
i 1
N
 
 Fij  rij
i
j i
here < ... > is the ensemble average; we assume the ergodicity hypothesis:
time average = ensemble average
a simple model potential for soft spheres
WCA potential:
u(r )  4 [( / r )12  ( / r )6 ]  
 σ=1, ε=1
 cut off at rcut=21/6σ,
corresponds to minimum
of the Lennard-Jones
potential
a simple problem
N = 144 particles in 2 dimensions:
 starting configuration: particles
sit on square lattice
 velocities from Maxwell-Boltzmann
distribution (T = 1.0)
 potential energy U = 0 due to
d > rcut
d>rcut
L=14.0
how do we integrate the equations of motion for this system?
the Euler algorithm
Taylor expansion with respect to discrete time step δt



ri (t  t )  ri (t )  v i (t )t



Fi (t )
v i (t  t )  v i (t ) 
t
mi
bad algorithm: does not recover time reversibility property of Newton‘s
equations of motion
→ unstable due to strong energy drift, very small time step required
the Verlet algorithm
consider the following Taylor expansions:




1 Fi (t ) 2
ri (t  t )  ri (t )  v i (t )t 
t
(1)
2 mi




1 Fi (t  t ) 2
ri (t )  ri (t  t )  v i (t  t )t 
t (2)
2
mi
addition of Eqs. (1) and (2) yields:



t 
v i (t  t )  v i (t ) 
[ Fi (t )  Fi (t  t )]
2mi
(3)
equations (1) and (2): velocity form of the Verlet algorithm
symplectic algorithm: time reversible and conserves phase space volume
implementation of velocity Verlet algorithm
do i=1,3*N ! update positions
displa=hstep*vel(i)+hstep**2*acc(i)/2
pos(i)=pos(i)+displa
fold(i)=acc(i)
enddo
do i=1,3*N ! apply periodic boundary conditions
if(pos(i).lt.0) then
pos(i)=pos(i)+lbox
else
if(pos(i).gt.lbox) pos(i)=pos(i)-lbox
endif
enddo
call force(pos,acc) ! compute forces on the particles
do i=1,3*N ! update velocities
vel(i)=vel(i)+hstep*(fold(i)+acc(i))/2
enddo
neighbor lists
force calculation is the most time-consuming part in a MD simulation
→ save CPU time by using neighbor lists
 Verlet list:
update when a particle has
made a displacement
> (rskin-rcut)/2
 cell list:
most efficient is a combination
of Verlet and cell list
simulation results I
time step t  0.0007 MD with  MD  m 2 / 
initial configuration with u=0:
after 105 time steps:
energies as a function of time
→ determine histogramm P(u) for the potential energy for different N
pressure
fluctuations of the potential energy
 dotted lines: fits with
Gaussians
 width of peaks is directly
related to specific heat cV
per particle at constant
volume V
u
2
 u
2
dk B2T 2
d

(1 
)
2N
2cV
 compare to canonical
ensemble
2
etot
 etot
2
kB2T 2

cV
N
simulations at constant temperature: a simple thermostat
idea: with a frequency ν assign new velocities to randomly selected
particles according to a Maxwell-Boltzmann (MB) distribution
with the desired temperature
simple version: assign periodically new velocities to all the particles
(typically every 150 time steps)
algorithm:
(i)
take new velocities from distribution
(ii) total momentum should be zero
(iii) scale velocities to desired temperature
P(v ,i )  exp(  v2 ,i / 2)
v , i  v , i
v , i  v , i
m v


i
,i
Nmi
3NkBT
 mi v , i
energies for the simulation at constant temperature
pair correlation function
 pair correlation function
very similar to 3d
 no finite size effects
visible
static structure factor
 S(q) also similar to 3d
 no finite size effects
 low compressibility
indicated by small
value of S(q) for q→0
mean squared displacement
D  lim
t 
r 2 (t )
4t
diffusion constant increases with increasing system size !?
long time tails in 3d

Green-Kubo relation for diffusion constant:
low density

1 
D   vi (t )  vi (0) dt
d0
high density
cage effect in
dense liquid:
 t 3 / 2
Alder, Wainwright,
PRA 1, 18 (1970)
Levesque, Verlet,
PRA 2, 2514 (1970)
→ power law decay of velocity autocorrelation function (VACF) at long times
Alder‘s argument
Stokes equation describes momentum diffusion:


jt
2
   jt
t


(  kinematic viscosity , jt  transvers e mass current)

backflow pattern around
particle:
consider momentum in volume element δV at t=0
volume at time t:
 (t ) d / 2
amount of momentum in original volume
element at time t →
VACF  V /(t ) d / 2  t  d / 2
1
 d=2: VACF  t  D diverges logarithmi cally
 refined theoretical prediction VACF  (t ln( t ) ) 1
effective diffusion constants
→ larger system sizes required to check theoretical prediction !
literature
 M. Allen, D. J. Tildesley, Computer Simulation of Liquids
(Clarendon Press, Oxford, 1987)
 D. Rapaport, The Art of Molecular Dynamics
(Cambridge University Press, Cambridge, 1995)
 D. Frenkel, B. Smit, Understanding Molecular Simulation: From
Algorithms to Applications
(Academic Press, San Diego, 1996)
 K. Binder, G. Ciccotti (eds.), Monte Carlo and Molecular Dynamics
of Condensed Matter Systems
(Societa Italiana di Fisica, Bologna, 1996)
Computer Simulations of Colloids
Molecular Dynamics II:
Application to Systems Under Shear
with Jochen Zausch (Universität Mainz)
shear viscosity of glassforming liquids
1013 Poise →
10-3 Poise →
dramatic slowing down of dynamics in a relatively small temperature range
F(q,t) for a glassforming colloidal system
3.5
S(q)
1.0

 = 0.620
3.0

0.8
S(q,t)/S(q)
2.5
2.0
1.5
1.0

0.5
0.0
0.010
0.015
0.020
0.025
0.4
0.2
0.030

q [nm ]
 = 0.581
0.6
 = 0.513
0.0
-6
-5
-4
-3
-2
-1
0
1
2
3
10 10 10 10 10 10 10 10 10 10
Eckert, Bartsch, Faraday Discuss. 123, 51 (2003)
t/s
glassforming fluids under shear
hard sphere colloid glass: confocal microscopy
apply constant shear rate 
drastic acceleration of dynamics
Weissenberg number:
W   (T )
τ(T): relaxation time in equilibrium
Besseling et al., PRL (2007)
 central question: transient dynamics towards steady state for W >>1
 problem amenable to experiment, mode coupling theory and simulation
outline
 simulation details: interaction model, thermostat, boundary conditions
 dynamics from equilibrium to steady state?
 dynamics from steady state back to equilibrium?
technical requirements
(1) model potential
(2) external shear field produces heat: thermostat necessary
(3) can we model the colloid dynamics realistically?
(4) how can we shear the system without using walls?
model potential: a binary Yukawa mixture
requirements: → binary mixture (no crystallization or phase separation)
→ colloidal system (more convenient for experiments)
→ possible coupling to solvent: density not too high
U  (r ) 
   
exp[  (r    )]
r
 1.2 AA ,  AB  1.1 AA
 BB
 BB  2 AA ,  AB  1.4 AA
  6 /  AA , cut - off U  (rc )  10 7
equations of motion and thermostat
 use of thermostat necessary + colloid dynamics → solve Langevin equation
 
r
 equations of motion: i  pi / mi
 D R

pi  [ Fij  Fij  Fij ]
j i
 dissipative particle dynamics (DPD):
D
R
rˆij

2
Fij  w (rij )( rˆij  v ij )rˆij , Fij  2k BT w(rij )ij
t
w(rij )  1 
rij
1.25 AA
for rij  1.25 AA
ij   ji : uniform random numbers with zero mean and unit varia nce
→ Galilean invariant, local momentum conservation, no bias on shear profile
 low ξ (ξ=12): Newtonian dynamics, high ξ (ξ=1200): Brownian dynamics
Lees-Edwards boundary conditions
 no walls, instead modified periodic
boundary conditions
 shear flow in x direction,
velocity gradient in y direction,
“free“ z direction
 shear rate:
  2v x / L
linear shear profile
vs ( y)   ( y  L / 2)
further simulation details
2
 generalized velocity Verlet algorithm, time step t  0.0083 mA AA
/  AA
 system size: equimolar AB mixture of 1600 particles (NA=NB=800)
 L=13.3 AA , density ρ=0.675 AA , volume fraction Φ≈48%
3
 at each temperature at least 30 independent runs
 temperature range 1.0 ≥ T ≥ 0.14 (40 million time steps for production
runs at T=0.14)
self diffusion
 B particles slower than
A particles
 weak temperature
dependence of D for
sheared systems if
W   2 / D  1
 along T = 0.14:
W ~ 100 - 1000
shear stress: equilibrium to steady state
definition:  xy 
1
L3
 [m ( v
i
i
i,x
 vs ( yi )) vi , y   rij, x Fij, y ]
j i
averaging over 250
independent runs
 maximum around t  0.1 marks
transition to plastic flow regime
 steady state reached on
time scale t  1
 time scale on which linear
profile evolves?
shear profile
 velocity profile becomes almost linear in the elastic regime
 maximum in stress not related to evolution of linear profile
structural changes: pair correlation function
pair correlation function not
very sensitive to structural
changes
→ consider projections of g(r)
onto spherical harmonics
structure and stresses
2

 U (r ) xy 
U (r )
 xy   xykin   dr
g (r )   xykin  A dr r 3
Im[ g 22 (r )]
2
r r
r
0
how is the stress overshoot reflected in dynamic correlations?
mean squared displacement
 occurrence of superdiffusive regime in the transient plastic flow regime
 superdiffusion less pronounced in hard sphere experiment (→ Joe Brader)
EQ to SS: superdiffusive transient regime
 y: gradient direction
z: free direction
 measure effective exponent:
s
d log( r 2 (t , t w ) )
d log( t  t w )
EQ to SS: Incoherent intermediate scattering function
 transients can be described by
compressed exponential decay:
Fs (q, t )  exp[ (t /  q )  ]
with   1
tw=0: β ≈ 1.8 → different for
Brownian dynamics?
EQ to SS: Newtonian vs. overdamped dynamics
change friction coefficient
in DPD forces:
  12 : Newtonian dynamics
  1200 : " Brownian" dynamics
→ same compressed exponential decay also for overdamped dynamics
shear stress: from steady state back to equilibrium
 decay can be well described by stretched exponential with β≈0.7
 stressed have completely decayed at t
1
SS to EQ: mean squared displacement
 stresses relax on time scale of
the order of 300 (“crossover“ in
tw=0 curve)
 then slow aging dynamics toward
equilibrium
conclusions
phenomenology of transient dynamics in glassforming liquids under shear:
 binary Yukawa mixture:
model for system of charged colloids, using a DPD thermostat and
Lees-Edwards boundary conditions
 EQ→SS: superdiffusion on time scales between the occurrence of the
maximum in the stress and the steady state
 SS→EQ: stresses relax on time scale 1/ , followed by slow aging
dynamics
Computer Simulations of Colloids
Monte Carlo Simulation I:
How does it work?
Outline
MC – how does it work?
 calculation of π via simple sampling and Markov chain sampling
 Metropolis algorithm for the canonical ensemble
 MC at constant pressure
 MC in the grandcanonical ensemble
calculation of π
 first method: integrate over unit circle
 second method: use uniform random numbers
→ direct sampling
→ Markov chain sampling
Werner Krauth, Statistical Mechanics: Algorithms and Computations
(Oxford University Press, Oxford, 2006)
calculation of π : direct sampling
Acircle 

 estimate π from
Asquare 4
 algorithm:
Nhits=0
do i=1,N
x=ran(-1,1)
y=ran(-1,1)
if(x2+y2<1) Nhits=Nhits+1
enddo
 estimate of π from ratio Nhits/N
calculation of π : Markov chain sampling
 algorithm:
Nhits=0; x=0.8; y=0.9
do i=1,N
Δx=ran(-δ,δ)
Δy=ran(-δ,δ)
if(|x+Δx|<1.and.|y+Δy|<1) then
x=x+Δx; y=y+Δy
endif
if(x2+y2<1) Nhits=Nhits+1
enddo
 optimal choice for δ range:
not too small and not too large
 here good choice δmax=0.3
calculation of π with simple and Markov chain sampling
simple sampling (N = 4000):
Markov chain sampling
(N = 4000, δmax = 0.3)
run
Nhits
estimate of π
run
Nhits
estimate of π
1
2
3
4
5
3156
3150
3127
3171
3148
3.156
3.150
3.127
3.171
3.148
1
2
3
4
5
3123
3118
3040
3066
3263
3.123
3.118
3.040
3.066
3.263
simple sampling vs. Markov chain sampling
N hits 1 N
  Ai  A 
N
N i 1
observable: A  1 inside circle
0 elsewhere
1
1
1
1
1
1
1
1
 dx dy  ( x, y) A( x, y)
 dx dy  ( x, y)
1 inside square
probability
 ( x, y )  
density:
elsewhere
0
 simple sampling: A sampled directly through π
 Markov sampling: acceptance probability
p(old→new) given by
 (new)
p(old  new)  min [1,
]
 (old)
more on Markov chain sampling
 Markov chain: probability of generating configuration i+1 depends
only on the preceding configuration i
 important: loose memory of initial condition
 consequence of algorithm: points pile up at the boundaries
shaded region:
piles of sampling points
δmax
detailed balance
 consider discrete system: configurations
a, b, c should be generated with equal
probability
{ (a),  (b),...} : stationary probabilit ies of the system to be at a, b,...
{ p(a  b), p(a  c),...} : transitio n probabilit ies
 detailed balance condition holds
 (a) p(a  b)   (b) p(b  a)
 ( a ) p ( a  c )   ( c ) p (c  a )
detailed balance in the calculation of π
 acceptance probability for a trial move
 (new)
p(old  new)  min[ 1,
]
 (old)
 one can easiliy show that
 (old) p(old  new)   (new) p(new  old)
detailed balance: analog to the time-reversibility property of Newton‘s
equations of motion
a priori probabilities
 moves Δx and Δy in square of linear size δ defines a priori
probability pap(old→new)
 different choices possible
P(old  new) 
pap (old  new) p(old  new)
 generalized acceptance criterion:
p(old  new)  min[ 1,
 (new)
pap (new  old)
pap (old  new)
 (old)
]
canonical ensemble
 problem: compute expectation value
(N )
(N )
(N )
dr exp[ H (r ) /( k BT )] A(r )
(N )

A(r ) 
(N )
(N )
 dr exp[  H (r ) /( kBT )]
(N )
 

(N )
r  (r1 , r2 ,......, rN ); H (r ) : Hamiltonia n
 direct evaluation of the integrals: does not work!
(N )
 simple sampling: sample random configurations rl
(N )
(N )
exp[

H
(
r
)
/(
k
T
)]
A
(
r
)

l
B
l
M
(N )
A(r ) 
l 1
(N )
exp[

H
(
r
) /( k BT )]

l
M
l 1
samples mostly states in the tails of the Boltzmann distribution
idea of importance sampling
(N )
(N )
r
 sample configurations l
according to probability P(rl )
(N )
(N )
(N )
 exp[  H (rl ) /( kBT )] A(rl ) / P(rl )
M

A(r ( N ) ) 
l 1
(N )
(N )
 exp[  H (rl ) /( kBT )] / P(rl )
M
l 1
(N )
(N )
 choose P (rl )  exp[  H (rl ) /( k BT )], then
(N )
1
A(r ) 
M
(N )
 A(rl )
M
l 1
 this can be realized by Markov chain sampling: random walk through
regions in phase space where the Boltzmann factor is large
importance sampling
(N )
(N )
 choose transition probability W (rl  rl ) such that
(N )
1
Peq  exp[  H (rl ) /( k BT )]
Z
for
M 
 to achieve this use detailed balance condition
(N )
(N )
(N )
(N )
(N )
(N )
Peq (rl )W (rl  rl 1 )  Peq (rl 1 )W (rl 1  rl )
(N )
(N )
W (rl  rl 1 )
H
or
)
(N )
 ( N )  exp( 
W (rl 1  rl )
k BT
 possible choice
(N )
 ( N ) exp[ H /( k BT )] H  0
W (rl  rl 1 )  
1
otherwise

 yields no information about the partition sum
Metropolis algorithm
N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, JCP 21, 1087 (1953)
old configuration
 displacement move
rtrial  rold  
trial configuration
ΔE = E(trial) - E(old)
yes
new config. =
trial config.
ΔE ≤ 0 ?
 acceptance probability
W (old  new) 
min(1, exp[( E ( trial)  E (old))/ (k BT )]
no
r < exp[-ΔE/(kBT)] ?
no
new config. = old config.
 r uniform random
number 0 < r < 1
implementation of the Metropolis algorithm
do icycl=1,ncycl
old=int[ran(0,1)*N]+1
call energy(x(old),en_old)
xn=x(old)+(ran(0,1)-0.5)*disp_x
call energy(x(new),en_new)
! select a particle at random
! energy of old conf.
! random displacement
! energy of new conf.
if(ran(0,1).lt.exp[ -beta*(en_new-en_old) ] ) x(old)=x(new)
if( mod(icycl,nsamp).eq.0) call sample
enddo
! sample averages
dynamic interpretation of the Metropolis precedure
(N )
(N)
P(r , t ) : probability that at time t a configuration r occurs during
the Monte Carlo simulation
rate equation (or master equation):






dP(r , t )

  W ( r  r ' ) P (r , t )   W (r '  r ) P (r ' , t )


dt
r'
r'


t   : P(r , t )  Peq (r )
equilibrium:






W (r  r ' ) Peq (r )  W (r '  r ) Peq (r ' )

r'

r'



W (r  r ' ) Peq (r ' )
detailed balance condition one possible solution:

 

W (r '  r ) Peq (r )
remarks
 Monte Carlo measurements:
first equilibration of the system (until Peq is reached), before
measurement of physical properties
 random numbers:
real random numbers are generated by a deterministic algorithm and
thus they are never really random (be aware of correlation effects)
 a priori probability:
can be used to get efficient Monte Carlo moves!
Monte Carlo at constant pressure
system of volume V in contact with an ideal gas reservoir via a piston
acceptance probability for volume
moves:
V0 - V
V
W (old  new ) 


1
min[ 1, exp{ 
( E ( s ( N ) ,V ' )  E ( s ( N ) , V ))
k BT
 P(V 'V ) 
N
V'
ln( )}]
k BT
V
trial move:
V '  V  V with V uniform random number from interval [Vmax ,Vmax ]
grand-canonical Monte Carlo
system of volume V in contact with an ideal gas reservoir with
fugacity z  exp[  ]
k BT
V0 - V
V
insertion:
W ( N  N  1) 
Vz
1
min[ 1,
exp{ 
( E ( N  1)  E ( N ))}]
N 1
k BT
removal:
W ( N  N  1) 
N
1
min[ 1, exp{ 
( E ( N  1)  E ( N ))}]
Vz
k BT
textbooks on the Monte Carlo simulation method
 D. P. Landau, K. Binder, A Guide to Monte Carlo Simulations in
Statistical Physics
(Cambridge University Press, Cambridge, 2000)
 W. Krauth, Statistical Mechanics: Algorithms and Computations
(Oxford University Press, Oxford, 2006)
 M. Allen, D. J. Tildesley, Computer Simulation of Liquids
(Clarendon Press, Oxford, 1987)
 D. Frenkel, B. Smit, Understanding Molecular Simulation: From
Algorithms to Applications
(Academic Press, San Diego, 1996)
 K. Binder, G. Ciccotti (eds.), Monte Carlo and Molecular Dynamics
of Condensed Matter Systems
(Societa Italiana di Fisica, Bologna, 1996)
Computer Simulations of Colloids
Monte Carlo II:
Phase Behaviour of Colloid-Polymer
Mixtures
collaborations: Richard Vink, Andres De Virgiliis, Kurt Binder
depletion interactions in colloid-polymer mixtures
polymers cannot move into
depletion zones of colloids
with their center of mass
 gain of free volume if there is
overlap of the depletion zones
of two colloids
 effective attraction between
colloids (of entropic origin)
→ demixing transition possible into colloid-rich phase (liquid) and
colloid-poor phase (gas)
Asakura-Oosawa (AO) model
 r  2 Rc
ucc (r )  
 0 otherwise
 r  Rc  Rp
ucp (r )  
 0 otherwise
upp (r )  0
in the following AO model for Rp/Rc = 0.8:
→ exhibits demixing transition into a colloid-poor phase (gas) and a
colloid-rich phase (liquid)
→ temperature not relevant; what is the corresponding variable here?
Outline
phase diagram:
 pr  zP
4 3
4 3 N C
RP , C 
RC
3
3
V
AO model for RP/RC = 0.8:
 binodal
 critical point (FSS)
 interfacial tension
 interfacial width
 capillary waves
grandcanonical MC for the AO model
 volume V, polymer fugacity zP and colloid fugacity zC fixed → removal
and insertion of particles
 problem 1: high free energy barrier between coexisting phases in the
two-phase region (far away from the critical point)
→ use successive umbrella sampling
 problem 2: low acceptance rate for trial insertion of colloidal particle, high
probability that it overlaps with a polymer particles
→ use cluster move
successive umbrella sampling
problem: high free energy barrier between coexisting phases
cluster move
problem: low acceptance rate for trial insertion of a colloid → cluster move
sphere of radius δ and volume Vδ around randomly selected point
try to replace nr ≤ nP polymers (A)
by a colloid (B):
A  min [1,
zCV
(nP )!
 E
]
n e
N C  1 (nP  nr )!( zPV )
r
reverse move (C+D):
n
N C (nP )!( zPV ) r  E
A  min [1,
e
]
(nP  nr )!
zCV
R.L.C. Vink, J. H., JCP 121, 3253 (2004)
phase diagram
critical point from finite size scaling
→ 3D Ising universality class
→ estimate interfacial tension by γ = ΔF/(2A) (Binder 1982)
finite size scaling
cumulant ratio
M  m / | m|
2
2
with m  C  C
→ critical value
ηrP,cr≈0.766
interfacial tension
γ*≡(2RC)2γ as a function of the difference in the colloid packing
fraction of the liquid (L) and the vapor (V) phase at coexistence
→ DFT (M. Schmidt et al.)
yields accurate result
interfacial profile
mean field result:
 z  z0 

( z )  A  B tanh 
 w0 
, Lx,y=31.3
, Lx,y=23.1
→ strong finite size effects
(not due to critical fluctuations)
capillary wave theory
 free energy cost of spatial interfacial fluctuations (long
wavelength limit)
2
2

1
 h   h  
Fs    dxdy     
2
 x   y  
 result for mean square amplitude of the interfacial thickness
w  2  ln( L / a)
2
cw
1
 combination with mean-field result by convolution approximation
w2  w02  (4 ) 1 ln L  (4 ) 1 ln a
→ determination of intrinsic width w0 by simulation not possible
is CWT valid?
logarithmic divergence?
lateral dimension Lx = Ly large enough?
finite size scaling II
 map on universal 3D Ising
distribution
 account for field mixing
→ AO model belongs to 3D Ising
universality class
conclusions
discrepancies to mean-field results

AO model belongs to 3D Ising universality class

interfacial broadening by capillary waves
much more on the AO model in
R.L.C. Vink, J.H., JCP 121, 3253 (2004);
R.L.C. Vink, J.H., JPCM 16, S3807 (2004);
R.L.C. Vink, J.H., K. Binder, JCP 122, 134905 (2005);
R.L.C. Vink, J.H., K. Binder, PRE 71, 011401 (2005);
R.L.C. Vink, M. Schmidt, PRE 71, 051406 (2005);
R.L.C. Vink, K. Binder, J.H., Phys. Rev. E 73, 056118 (2006);
R.L.C. Vink, K. Binder, H. Löwen, PRL 97, 230603 (2006);
R.L.C. Vink, A. De Virgiliis, J.H. K. Binder, PRE 74, 031601 (2006);
A. De Virgiliis, R.L.C. Vink, J.H., K. Binder, Europhys. Lett. 77, 60002 (2007);
K. Binder, J.H., R.L.C. Vink, A. De Virgiliis, Softmatter (2008).
Computer Simulations of Colloids
Modelling of hydrodynamic interactions
collaboration: Apratim Chatterji (FZ Jülich)
simulation of colloids
colloids:
 size 10 nm – 1 μm
 mesoscopic time scales
solvent:
 atomistic time (~ 1 ps) and
length scales (~ Å)
simulation difficult due to different time and length scales
→ coarse graining of solvent‘s degrees of freedom
Langevin equation
 many collisions with solvent particles
on typical time scale of colloid
 systematic friction force on colloid
 


d ri
M 2  Fc,i   0 vi (t )  f r ,i
dt
2
Brownian particles of mass M :
fr,i uncorrelated random force with zero mean:

f (r , t )  0 with   x, y, z
 
( ) 
( ) 
f r ,i (r , t ) f r ,i (r ' , t ' )  A   (r  r ' ) (t  t ' )
( )
r,i
fluctuation-dissipation theorem:
A  2k BT 0
decay of velocity correlations
consider single Brownian particle of mass M :
d v x (t )
d
M
  0 v x (t )  f r , x  M
v x (t ) v x (0)   0 v x (t ) v x (0)
dt
dt
→ exponential decay of velocity autocorrelation function (VACF)
0
v x (t ) v x (0)  exp(  t )
M
correct behavior: power law decay of VACF due to local momentum
conservation
v x (t )vt (0)  t
3 / 2
generalized Langevin equation
 idea: couple velocity of colloid to the
local fluid velocity field via frictional force
 point particle:


 
Ffric (r , t )  0 [vi (t )  u (r , t )]
problems:
 determination of hydrodynamic velocity field
of Navier-Stokes equation
 
u (r , t ) from solution
 rotational degrees of freedom → colloid-fluid coupling?
lattice Boltzmann method I



ni (r , t ): number
of particles at a lattice node r , at a time t, with

a velocity
ci
 discrete analogue of kinetic equation:
 


ni (r  ci , t  1)  ni (r , t )  i (r , t )
Δi: collision operator
 velocity space from projection of 4D
FCHC lattice onto 3D (isotropy)
 moments of ni :
18

 18 

   ni , j  u   ni ci ,    ni ci ci
18
i 1
i 1
i 1
lattice Boltzmann method II

 linearized collision operator:  i   Lij [n j (r , t )  n eq
j ]
18
j 1
 equilibrium distr.:
1   1   
n  a [   2 u  ci  4 uu : (ci ci  cs21)]
cs
cs
eq
i
ci
 recover linearized Navier-Stokes equations:

 j

 2
   j ,
 p   j
t
t

 thermal fluctuations: add noise terms to stress tensor
lattice Boltzmann algorithm III
 collision step: compute


ni (r , t )  i (r , t )
 propagation step: compute
 
ni (r  ci , t  1)
colloid fluid momentum exchange
sphere represented by (66)
uniformly distributed points
on its surface
each point exchanges
momentum with surrounding
fluid nodes
velocity correlations at short times
give particle a kick and determine decay of its velocity
blue curves:
g (t )  exp( 
0
M
red curves:
fluid velocity at
the surface of the
sphere
t)
Cv(t) and Cω(t)
M0
A
( ) 3 / 2
12 
independent of R and ξ0
I
B  (4 ) 5 / 2

I: moment of inertia
thermal fluctuations ok?
linear response theory:
velocity relaxation of kicked particle in a fluid at rest =
velocity autocorrelation function of particle in a thermal bath
charged colloids
colloids in electric field