Monte Carlo Method (3)
Application of MR2T2 algorithm in statistical
mechanics
Canonical ensemble:
Distribution function:
r(rN)=exp(- E(rN)/kT)/Q
Q =?drN exp(-E(rN)/kT)
Different states of the Markov chain:
In the application of Metropolis algorithm to a simulation of molecular
system,the states of the Markov chain correspond to different
configurations of the molecular system,
How to implement the MR2T2 algorithm for a
molecular simulation?
Configuration evolution in two steps:
1) attempt move (p*ij);
2) acceptation check (pj/pi).
Practical recipes for configuration evolution:
Start with an initial configuration,It is evolved by displacement of
the position of a particle,either
1) chosen randomly
or
2) in sequence (from 1 to N).
Attempt move in practice:
RXINEW=RX(I) + (2.0*RANF1-1)*DRMAX
RYINEW=RY(I) + (2.0*RANF2-1)*DRMAX
RZINEW=RZ(I) + (2.0*RANF3-1)*DRMAX
Important note:
The displacement must be centered at the old position,i.e.,
xnew-xold,ynew-yold,znew-zold are uniformly distributed in (-D,D)
(D=DRMAX).
In this case,p*ij=1/(2D).
Acceptation check in practice:
1) calculate DE=Ej-Ei:
No need to calculate the energy of the whole system!
2) If DE?0,accept the move.
3) If DE>0,calculate exp(-DE/kT) (exp(-DE/kT)=pj/pi) and accept
the move with a probability equal to exp(-DE/kT).
How to do this in practice?
Generate a uniform random number on (0,1),x,If x<exp(-DE/kT),
accept the move.
To understand this easily,make an analogy with the dice throwing.
D ij o l dijij n e wij rr uuE )()(
Particular case of hard spheres
1) Attempt move:
As in the general case.
2) Acceptation check (simplified)
If the displaced particle has no overlapping with any other particle,
accept the move.
Isothermal-isobaric ensemble:
Distribution function:
r(rN,V) = exp( - [U(rN) + PV]/kT)/ QTPN
QTPN,configuration integral
Different states of the Markov chain:
Now,the states of the Markov chain are characterized not only by
system’s configurations but also by the system size since the volume
fluctuates in this ensemble.
An important remark:
A straightforward extension of Metropolis algorithm seems to take
pi? exp( - [Ui + PVi]/kT).
WRONG!
The correct Markov chain is generated by using
pi? exp( - [Ui(sN)+ PVi + NlnVi]/kT)
where s=r/L (L,box length).
Why?
In this ensemble,the volume fluctuates and becomes a random
variable,To identify correctly the distribution function of all the
random variables,one must express all of them as explicitly as
possible.
Practical procedure:
V
NNN
N P T
rrrQ UPVAddVA )((e x p)(1
0
10
)((e x p)(1 sssVQ NNNN
N P T
UPVAddV
Implementation
The Monte Carlo moves in an isothermal-isobaric ensemble include
particle displacement and volume change.
Practical recipe for MC moves:
Separate particle displacements and volume change;
A volume change every a few displacement cycles.
Particle displacement:
)12( 1?D xss ol dxne wx
)12( 2?D xss o ldyn e wy
)12( 3?D xss o ldzn e wz
The displacement is accepted with a probability equal to
min(1,exp(-?DU)).
Volume change:
)12( 4 xLL o ldn e w
Important note:
Volume change causes necessarily potential energy change,So,even
in the case of only a volume change,one MUST calculate also DU.
The volume change is accepted with a probability equal to
min{1,exp(-?[DU+P(Vj-Vi)])Vj/Vi}
Grand canonical ensemble
Distribution function:
r(rN,N) = exp( - [U(rN) - Nm]/kT -3Nln?)/?(T,V,m)
(T,V,m),partition function
= (2p 2/m)1/2,thermal wave length
Is it correct to generate the Markov chain by using
pi? exp( - [Ui + mNi]/kT - 3Niln?)?
Answer:
No!
The correct method is
pi? exp( - [Ui + mNi]/kT - 3Niln? + NlnV).
The reason is
13
))((e x p)(!1 NVAdNA sssV NNNN N m
Implementation
The Monte Carlo moves in a grand canonical ensemble include
particle displacement,particle insertion and particle removal.
Particle displacement:
The same procedure as in a canonical ensemble.
Particle insertion:
Choose a random position and insert a particle at this position.
The insertion is accepted with a probability equal to
min{1,exp[-?(DU-m)]V/(N?3)}
Particle removal:
Choose randomly a particle and remove it from the system.
The removal is accepted with a probability equal to
min{1,exp[-?(DU+m)]N?3/V}
Pitfall:
If care is not taken,a non symmetric transition probability can be
generated with respect to insertion and removal.
Remedy:
Carry out the simulation so that the ratio of successful insertions is
equal to that of successful removals.
Practical recipe:
pd=pi=pr=1/3 => fast convergence,
Limit of the Grand Canonical Monte Carlo (GCMC) Method:
Break down at high densities,rs3>0.7.
Application of MR2T2 algorithm in statistical
mechanics
Canonical ensemble:
Distribution function:
r(rN)=exp(- E(rN)/kT)/Q
Q =?drN exp(-E(rN)/kT)
Different states of the Markov chain:
In the application of Metropolis algorithm to a simulation of molecular
system,the states of the Markov chain correspond to different
configurations of the molecular system,
How to implement the MR2T2 algorithm for a
molecular simulation?
Configuration evolution in two steps:
1) attempt move (p*ij);
2) acceptation check (pj/pi).
Practical recipes for configuration evolution:
Start with an initial configuration,It is evolved by displacement of
the position of a particle,either
1) chosen randomly
or
2) in sequence (from 1 to N).
Attempt move in practice:
RXINEW=RX(I) + (2.0*RANF1-1)*DRMAX
RYINEW=RY(I) + (2.0*RANF2-1)*DRMAX
RZINEW=RZ(I) + (2.0*RANF3-1)*DRMAX
Important note:
The displacement must be centered at the old position,i.e.,
xnew-xold,ynew-yold,znew-zold are uniformly distributed in (-D,D)
(D=DRMAX).
In this case,p*ij=1/(2D).
Acceptation check in practice:
1) calculate DE=Ej-Ei:
No need to calculate the energy of the whole system!
2) If DE?0,accept the move.
3) If DE>0,calculate exp(-DE/kT) (exp(-DE/kT)=pj/pi) and accept
the move with a probability equal to exp(-DE/kT).
How to do this in practice?
Generate a uniform random number on (0,1),x,If x<exp(-DE/kT),
accept the move.
To understand this easily,make an analogy with the dice throwing.
D ij o l dijij n e wij rr uuE )()(
Particular case of hard spheres
1) Attempt move:
As in the general case.
2) Acceptation check (simplified)
If the displaced particle has no overlapping with any other particle,
accept the move.
Isothermal-isobaric ensemble:
Distribution function:
r(rN,V) = exp( - [U(rN) + PV]/kT)/ QTPN
QTPN,configuration integral
Different states of the Markov chain:
Now,the states of the Markov chain are characterized not only by
system’s configurations but also by the system size since the volume
fluctuates in this ensemble.
An important remark:
A straightforward extension of Metropolis algorithm seems to take
pi? exp( - [Ui + PVi]/kT).
WRONG!
The correct Markov chain is generated by using
pi? exp( - [Ui(sN)+ PVi + NlnVi]/kT)
where s=r/L (L,box length).
Why?
In this ensemble,the volume fluctuates and becomes a random
variable,To identify correctly the distribution function of all the
random variables,one must express all of them as explicitly as
possible.
Practical procedure:
V
NNN
N P T
rrrQ UPVAddVA )((e x p)(1
0
10
)((e x p)(1 sssVQ NNNN
N P T
UPVAddV
Implementation
The Monte Carlo moves in an isothermal-isobaric ensemble include
particle displacement and volume change.
Practical recipe for MC moves:
Separate particle displacements and volume change;
A volume change every a few displacement cycles.
Particle displacement:
)12( 1?D xss ol dxne wx
)12( 2?D xss o ldyn e wy
)12( 3?D xss o ldzn e wz
The displacement is accepted with a probability equal to
min(1,exp(-?DU)).
Volume change:
)12( 4 xLL o ldn e w
Important note:
Volume change causes necessarily potential energy change,So,even
in the case of only a volume change,one MUST calculate also DU.
The volume change is accepted with a probability equal to
min{1,exp(-?[DU+P(Vj-Vi)])Vj/Vi}
Grand canonical ensemble
Distribution function:
r(rN,N) = exp( - [U(rN) - Nm]/kT -3Nln?)/?(T,V,m)
(T,V,m),partition function
= (2p 2/m)1/2,thermal wave length
Is it correct to generate the Markov chain by using
pi? exp( - [Ui + mNi]/kT - 3Niln?)?
Answer:
No!
The correct method is
pi? exp( - [Ui + mNi]/kT - 3Niln? + NlnV).
The reason is
13
))((e x p)(!1 NVAdNA sssV NNNN N m
Implementation
The Monte Carlo moves in a grand canonical ensemble include
particle displacement,particle insertion and particle removal.
Particle displacement:
The same procedure as in a canonical ensemble.
Particle insertion:
Choose a random position and insert a particle at this position.
The insertion is accepted with a probability equal to
min{1,exp[-?(DU-m)]V/(N?3)}
Particle removal:
Choose randomly a particle and remove it from the system.
The removal is accepted with a probability equal to
min{1,exp[-?(DU+m)]N?3/V}
Pitfall:
If care is not taken,a non symmetric transition probability can be
generated with respect to insertion and removal.
Remedy:
Carry out the simulation so that the ratio of successful insertions is
equal to that of successful removals.
Practical recipe:
pd=pi=pr=1/3 => fast convergence,
Limit of the Grand Canonical Monte Carlo (GCMC) Method:
Break down at high densities,rs3>0.7.
