《Simulations Moléculaires》 Cours04VII

Monte Carlo Method (3)

Monte Carlo Method (3)

Application of mr algorithm in statistical mechanics Canonical ensemble Distribution function p(r=exp(-e(rkT)/Q Q=drN exp(e(rn)/k 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

Application of MR2T2 algorithm in statistical mechanics Canonical ensemble: Distribution function: r(r N)=exp(- E(r N)/kT)/Q Q = dr N exp(-E(r N)/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 Mr t algorithm for a molecular simulation? Configuration evolution in two steps. ) attempt move(pii 2)acceptation check (t /T;) Practical recipes for configuration evolution the position of a particle, either It is evolved by displacement of Start with an initial configuration ) chosen randomly or 2)in sequence(from 1 to N)

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(D+(2.0*RANFI-1)*DRMAX RYINEW=RY+(2.0*RANF2-I)*DRMAX RZINEW=RZ(+(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(-A, 4) (△= DRMAX) In this case,pi=1(2△)

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△E=E;-E △E=2)分 No need to calculate the energy of the whole system 2)IfAE 0, accept the move 3)IfAE>O, calculate exp (-AE/kT)(exp (AE/kT=I /T and accept the move with a probability equal to exp(-AE/kT) How to do this in practice? Generate a uniform random number on(0, 1),. If s<exp(-AE/kt) accept the move To understand this easily, make an analogy with the dice throwing

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 = − j i old ij j i new rij r E u( ) u( )

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

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 P(r, v)=exp( -[U(r)+ PVkT)/ 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 ensembl An important remark A straightforward extension of Metropolis algorithm seems to take T, exp(-[Ui+ PVil/kT) WRONG/ The correct Markov chain is generated by using T, exp( -U(S)+ PVi+ Nnvil/kt where s=r/L (L: box length

Isothermal-isobaric ensemble: Distribution function: r(r N, V) = exp( - [U(r N) + 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 (s N)+ 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 4dmF)风PHF dry ds AS expl-apvauc

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 N N N NPT r r r Q A dV d A( )exp (PV U( ) 1 0         = − +  0 1 ( )exp ( ( ) 1 V s s s Q N N N N NPT dV d A PV U    

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 old S:=s;+A(251) =sy+A(2321) S:=:+12 The displacement is accepted with a probability equal to min(1,exp(-B△U)

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: (2 1) 1 = +D x − s s old x new x (2 1) 2 = +D x − s s old y new y (2 1) 3 = +D x − s s old z new z The displacement is accepted with a probability equal to min(1, exp(-DU))

Volume change L=L+21) Important note. in the case of only a volume change, one must calculate also a0 9 Volume change causes necessarily potential energy change. So, eve The volume change is accepted with a probability equal to minf l, exp (-B[AU+P(Vi-VODV/VA)

Volume change: (2 1) 4 L =L + x − new old 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}