Transcript Document

General considerations
Monte Carlo methods (I)
Averages
A 
   AX f XdV
X
x1
xN



f
X
dV

1
X
 
x1
xN
X = (x1, x2, …, xN) – vector of variables
P – probability density function
Random variables and their distributions
A random variable: a number assigned to an event
Distribuand:
F  X   PY  X 
dP X 
Probability density: f  X  
dX
Any function of a random variable is a random
variable.
Moments of a distribution
Continuous variables:

n
E({x}) 
 x P x
i 1
EH(x) 
 xi 
i

n
 H x Px
i 1
E {x}  xˆ 
i
 xi  EH x  
 xf x dx

 H x  f x dx

If H(x)=(x-xc)n then E{H(X)} is called the nth moment
of x with respect to c; if c= x̂ then E is the nth central
moment, mn({x}).
Some useful central moments
Variance


2
x  m 2 x   x  x̂  f x dx
2

Skewness
m 3 x
1
3
x  x̂  f x dx
 x  3 / 2
 3

m 2 x  x 

Kurtosis
m 4 x
1
4
x  x̂  f x dx  3
x  2
3  4

m 2 x
 x 

For a set of points:
1 n
x̂   x i
n i 1
2
n
n


1
1


2
2




2 
x

x̂

n
x


i
i    xi 
n  1 i 1
n  1  i 1
 i 1  
n
n

 x
i 1
i
(n  1)3
n

 x̂ 
3
 x
i 1
 x̂ 
4
i
(n  1)
4
3
Some examples of central moments
Boltzmann distribution
 Ei 
1

Pi  exp  
Q
 k BT 
 Ei 

Q   exp  
i
 k BT 
 E q, p   N N
1
Pq, p d qd p  exp 
 d qd p
Q
 k BT 
 E q, p  N N
Q    exp 
d qd p
 k BT 
q p
 E 
1
dE
PE dE   E  exp  
Q
 k BT 
 E 
dE
Q    E  exp  
 k BT 
E
N
N
kB (Boltzmann constant) = 1.3810-23 J/K
NAkB = R (universal gas constant) = 8.3143 J/(molK)
q – positions; p – momenta; (E) – density of states with energy E
Thermodynamic quantities
 i 
1
2   ln Q 
  k BT 
U  E T 
 i exp  


Q i
 T  N ,V
 k BT 
 i 
1
  i 
  ln Q 
 k BT 
p  p T 

 exp  


Q i  V 
 V  N ,T
 k BT 
  ln Q 
S  k BT  Pi ln Pi k BT 
 k B ln Q

 T  N ,V
i
F  k BT ln Q
2
 U 
2
2
CV  k BT E
 E T  

T
 T  N ,V


Ergodicity
t
1
A  lim
t  t

A d
0
1
A  lim
N  N
ergodic
theorem
A 
A
N
A
i 1
sampled
acc. to ρ
i
Trajectories in the phase space. For ergodic sampling, sub-trajectories should be
part of a trajectory that passes through all points.
Monte Carlo methods: use of random-number generators to
follows the evolution of a system according to a given
probability-distribution function.
Beginnings: Antiquity (?); estimation of the results of dice
game.
First documented use: G. Comte de Buffon (1777): estimation
of the number p by throwing a needle on a sheet of paper with
parallel lines and counting the number of hits.
First large-scale (wartime) application: J. von Neumann, S.
Ulam, N. Metropolis, R.P. Feynman i in. (1940’s; the
Manhattan project) calculations of neutron creation and
scattering. For security, the calculations were disguised as
,,Monte Carlo” calculations.
Kinds of Monte Carlo methods
 The von Neumann (rejection) sampling
 The Metropolis (in general: Markov chain)
sampling. Also known as „importance
sampling”
Simple Monte Carlo averaging
f(x)
Compute f(xi)
x
1
A 
N
N
Sample a point on x from a uniform distribution
 Ax  f x 
i 1
i
i
Rejection or hit-and-miss sampling
f(x)
reject
Sample a point on f
accept
x
Sample a point on x from a uniform distrubution
A 
1
N accepted
 Ax 
iaccepted
i
Algorithm
Generate a random
point X in the
configurational space
Generate a random
point y in [0,1]
Reject X
P(X)>y
Accept X
A:=A+A(X)
Application of the rejection sampling to
compute the number p
1
S

S
  p 1 2  p
N tot
1
4
N
lim
n
2
Applications
• For one-dimensional integrals classical
quadratures (Newton-Cotes, Gauss) are better
than that.
• For multi-dimensional integrals sampling the
integrations space is somewhat better; however
most points have zero contributions.
• We cannot compute ensemble averages of
molecular systems that way. The positions of
atoms would have to generated at random, this
usually leading to HUGE energies.
Illustration of the difference between the direct- and importancesampling methods to measure the depth of the river of Nile
Von Neumann: all points are visited
Metropolis: The walker stays in the river
Configuration Xo, energy Eo
Perturb Xo: X1 = Xo + DX
Compute the new energy (E1)
E1<Eo ?
N
Draw Y from U(0,1)
Compute W=exp[-(E1-Eo)/kT]
A:=A+A(Xo)
W>Y?
Y
N
Xo=X1, Eo=E1
Y
E1
E0
Accept
with probability
exp[-(E2-E1)/kBT]
Accept
E1
Space representation in MC
simulations
• Lattice (discrete). The particles are on
lattice nodes
• Continuous. The particles move in the 3D
space.
On- and off-lattice representations
Application of Metropolis Monte Carlo
• Determination of mechanical and
thermodynamic properties(density, average
energy, heat capacity, conductivity, virial
coefficients).
• Simulations of phase transitions.
• Simulations of polymer properties.
• Simulations of biopolymers.
• Simulations of ligand-receptor binding.
• Simulations of chemical reactions.
Computing averages with Metropolis Monte
Carlo
1
A
N
N
A
i 1
i
It should be noted that all MC steps are considered, including
those which resulted in the rejection of a new configuration.
Therefore, if a configuration has a very low energy, it will be
counted multiple times.
Detailed balance (Einstein’s theorem)
N (o)p (o  n)  N (n)p (n  o)
old
new
In real life the detailed balance
condition is rarely satisfied…
A famous Russian proverb states:
„A ruble to get in, ten rubles to get out”
These gates at a
Seoul subway
station do not satisfy
the detailed-balance
condition: you can
go through but you
cannot go back…
Nature teaches us the
hard way that detailed
balance is not something
to meet in the macroworld..
I thought there was
a way out….
I was sssoooo
busy working
for my Queen
and
Community….
It is only too easy to violate the
detailed-balance conditions
Little chance to
get from C to B
A
B
C
Monte Carlo with minimization: energy is minimized
after each move. Transition probability is proportional
to basin size.
Translational perturbations (straightforward)
D
x[o]
o : int N  ranf ()   1
xnew : xo  D  ranf ()  0.5
ynew : yo  D  ranf ()  0.5
z new : zo  D  ranf ()  0.5
Orientational perturbations
Euler angles. We rotate the
system first about the z
axis by f so that the new x
axis is axis w, then about
the w axis by q so that the
new z axis is z’ and finally
about the z’ axis by y so
that the new x axis is x’
and the new y axis is y’.
Uniform sampling the Euler
angles would result in a
serious bias.
Rotational perturbations: rigid linear molecules
Genarate a random
unit vector v
Orientational perturbation
u’
u
Compute t:=u+v
Compute u’:= t/||t||
Rigid non-linear molecules
Sample a unit vector on a 4D sphere (quaternion)
(E. Veseley, J. Comp. Phys., 47:291-296, 1982), then
compute the Euler angles from the following formulas:
Q  q0 , q1 , q2 , q3 ; q02  q12  q22  q32  1
q0 
q1 
q2 
q3 
q
  y 
cos cos

2
 2 
q
  y 
sin cos

2
 2 
q
  y 
sin sin 

2
 2 
q
  y 
cos sin 

2
 2 
Rotation matrix
 q02  q12  q22  q32

R   2q1q2  q0 q3 

 2q1q3  q0 q2 

r '  Rr
2q1q2  q0 q3 
q02  q12  q22  q32
2q2 q3  q0 q1 
2q1q3  q0 q2 


2q2 q3  q0 q1  

2
2
2
2
q0  q1  q2  q3 
Choosing the initial configuration and step size
• Generally: avoid clashes.
• For rigid liquids molecules start from a
configuration on a lattice.
• For flexible molecules: build the chain to
avoid overlap with the atoms already
present.
MC step size
• Too small: high acceptance rate but poor
ergodicity (can’t get out of a local minimum).
• Too large: low acceptance rate.
• Avoid accepting „good” advices that the
acceptance rate should be 10/20/50%, etc. Do
pre-production simulations to select optimal
step size
• This can help:
– Configurational-bias Monte Carlo,
– Parallel tempering.
Reference algorithms for MC/MD
simulations (Fortran 77)
M.P. Allen, D.J. Tildesley, „Computer Simulations of
Liquids” , Oxford Science Publications, Clardenon Press,
Oxford, 1987
http://www.ccp5.ac.uk/software/allen_tildersley.shtml
F11: Monte Carlo simulations of Lennard-Jones fluid.