[PPT] - Monte Carlo Methods Ensembles (Chapter 5) Biased Sampling (Chapter PowerPoint Presentation

SLIDE 1

Monte Carlo Methods

Ensembles (Chapter 5) Biased Sampling (Chapter 14) Practical Aspects

SLIDE 2

2

Lecture 1

SLIDE 3

3

Lecture 1&2

SLIDE 4

4

Lecture 1&3

SLIDE 5

5

Different Ensembles

Ensemble Name Constant (Imposed) Fluctuating (Measured) NVT Canonical N,V,T P NPT Isobaric-isothermal N,P,T V µVT Grand-canonical µ,V,T N

SLIDE 6

6

NVT

Liquids

Equation of State of Liquid Carbon

Imposed

SLIDE 7

7

NPT

Liquids

Equation of State of Liquid Carbon

Imposed

… if force is difficult to calculate … e.g. carbon force field

SLIDE 8

8

NPT

Non-Isotropic Systems e.g. Solids

Structure and Transformation of Carbon Nanotube Arrays

SLIDE 9

9

µVT

Adsorption

Adsorption in Carbon Nanostructues

Schwegler et al.

SLIDE 10

10

Statistical Thermodynamics

( )

NVT

Q F ln − = β

( ) ( )

[ ]

N N N N NVT NVT

U A N Q A r exp r dr ! 1 1

3

β − Λ =

∫

( )

[ ]

N N N NVT

U N Q r exp dr ! 1

3

β − Λ =

∫

Partition function Ensemble average Free energy

( ) ( ) ( ) ( )

3

1 1 N r dr' r' r exp r' exp r !

N N N N N N N NVT

U U Q N δ β β ⎡ ⎤ ⎡ ⎤ = − − ∝ − ⎣ ⎦ ⎣ ⎦ Λ

∫

Probability to find a particular configuration Lecture 2

SLIDE 11

Ensemble average

( ) ( ) [ ]

N N N N NVT NVT

U A N Q A r exp r dr ! 1 1

3

β − Λ =

∫

( ) ( )

∫

=

N N N

P A r r dr

( ) ( ) ( )

∫ ∫

=

N N N N N

P P A r dr r r dr

( ) ( ) [ ] ( ) [ ]

∫ ∫

− − =

N N N N N

U C U C A r exp dr r exp r dr β β

( ) ( ) [ ] ( ) [ ]

∫ ∫

− − =

N N N N N

U U A r exp dr r exp r dr β β

Generate configuration using MC:

( ) ( ) [ ]

N MC N MC

U C P r exp r β − =

{ }

N M N N N N

r , r , r , r , r

4 3 2 1

!

( )

N i M i

A M A r 1

1

∑

=

with

( ) ( ) ( )

∫ ∫

=

N MC N N MC N N

P P A r dr r r dr

( ) ( ) [ ] ( ) [ ]

∫ ∫

− − =

N MC N N MC N N

U C U C A r exp dr r exp r dr β β

( ) ( ) [ ] ( ) [ ]

∫ ∫

− − =

N N N N N

U U A r exp dr r exp r dr β β

Lecture 2 Weighted Distribution

SLIDE 12

Monte Carlo: Detailed balance

acc(o → n) acc(n → o) = N(n)×α(n → o) N(o)×α(o → n) = N(n) N(o)

( ) ( ) K o n K n

→

= → ( ) ( ) ( ) acc( ) K o n N o

n
n

α → = × → × →

n

( ) ( ) ( ) acc( ) K n

N n

n

n
α

→ = × → × →

Lecture 2

SLIDE 13

13

NVT-ensemble

acc( ) ( ) acc( ) ( )

n

N n n

N o

→ = →

( )

( ) exp N n U n β ⎡ ⎤ ∝ − ⎣ ⎦

( ) ( )

acc( ) exp acc( )

n

U n U o n

β

→ ⎡ ⎤ ⎡ ⎤ = − − ⎣ ⎦ ⎣ ⎦ →

Lecture 2

SLIDE 14

14

NPT ensemble

We control the

Temperature (T)
Pressure (P)
Number of particles (N)

SLIDE 15

15

Scaled coordinates

/

i i L

= s r

( )

[ ]

N N N NVT

U N Q r exp dr ! 1

3

β − Λ =

∫

Partition function Scaled coordinates This gives for the partition function

( ) ( )

3 3 3

ds exp s ; ! ds exp s ; !

N N N NVT N N N N N

L Q U L N V U L N β β ⎡ ⎤ = − ⎣ ⎦ Λ ⎡ ⎤ = − ⎣ ⎦ Λ

∫ ∫

The energy depends on the real coordinates Intermezzo

SLIDE 16

The NPT ensemble

Here they are an ideal gas Here they interact

What is the statistical thermodynamics of this ensemble?

N in volume V M-N in volume V0-V V0 is fixed V varies from 0 to V0 V0 : total volume M : total number of particles

SLIDE 17

17

The NPT ensemble: partition function

( )

3

ds exp s ; !

N N N NVT N

V Q U L N β ⎡ ⎤ = − ⎣ ⎦ Λ

∫

QMV0,NV ,T = V0 −V

( )

M −N

Λ3(M −N ) M − N

( )!

dsM −N

∫

exp −βU0 sM −N;L

( )

! " # $ V N Λ3N N! × dsN

∫

exp −βU sN;L

( )

! " # $ QMV0,NV ,T = V0 −V

( )

M −N

Λ3(M −N ) M − N

( )!

V N Λ3N N! dsN

∫

exp −βU sN;L

( )

$ % & '

Fixed V

SLIDE 18

18

QMV0,NV ,T = V0 −V

( )

M −N

Λ3(M −N ) M − N

( )!

V N Λ3N N! dsN

∫

exp −βU sN;L

( )

$ % & '

To get the Partition Function of this system, we have to integrate over all possible volumes:

QMV0,N ,T = dV V0 −V

( )

M −N

Λ3(M −N ) M − N

( )!

V N Λ3N N! dsN

∫

exp −βU sN;L

( )

$ % & '

∫

Now let us take the following limits:

constant M M V V ρ → ∞⎫ = → ⎬ → ∞ ⎭

As the particles are an ideal gas in the big reservoir we have: P ρ β =

SLIDE 19

19

We have

QMV0,N ,T = dV V0 −V

( )

M −N

Λ3(M −N ) M − N

( )!

V N Λ3N N! dsN

∫

exp −βU sN;L

( )

$ % & '

∫

( ) ( ) ( )

1 exp

M N M N M N M N

V V V V V V M N V V

− − − −

⎡ ⎤ − = − ≈ − − ⎣ ⎦

( )

[ ] [ ]

exp exp

M N M N M N

V V V V V PV ρ β

− − −

− ≈ − = − This gives:

[ ]

( )

3

d exp ds exp s ; !

N N N NPT N

P Q V PV V U L N β β β ⎡ ⎤ = − − ⎣ ⎦ Λ ∫

∫

To make the partition function dimensionless (see Frenkel/Smit for more info)

SLIDE 20

20

NPT Ensemble

Partition function:

[ ]

( )

3

d exp ds exp s ; !

N N N NPT N

P Q V PV V U L N β β β ⎡ ⎤ = − − ⎣ ⎦ Λ ∫

∫

Probability to find a particular configuration:

( )

[ ]

( )

, exp exp s ;

N N N NPT

N V V PV U L β β ⎡ ⎤ ∝ − − ⎣ ⎦ s

Sample a particular configuration:

change of volume
change of reduced coordinates

Acceptance rules ??

De Detailed balance

SLIDE 21

21

Detailed balance

acc( ) ( ) ( ) ( ) acc( ) ( ) ( ) ( )

n

N n n

N n

n

N o
n

N o α α → × → = = → × → ( ) ( ) K o n K n

→

= → ( ) ( ) ( ) acc( ) K o n N o

n
n

α → = × → × →

n

( ) ( ) ( ) acc( ) K n

N n

n

n
α

→ = × → × →

SLIDE 22

22

NPT-ensemble

acc( ) ( ) acc( ) ( )

n

N n n

N o

→ = →

( )

[ ]

( )

, exp exp s ;

N N N NPT

N V V PV U L β β ⎡ ⎤ ∝ − − ⎣ ⎦ s

Suppose we change the position of a randomly selected particle

[ ]

( )

[ ]

( )

N n N

exp

exp s ; acc( ) acc( ) exp exp s ;

N N

V PV U L

n

n

V

PV U L β β β β ⎡ ⎤ − − → ⎣ ⎦ = → ⎡ ⎤ − − ⎣ ⎦

( ) ( )

{ }

N n N

exp

s ; exp exp s ; U L U n U o U L β β β ⎡ ⎤ − ⎣ ⎦ ⎡ ⎤ = = − − ⎣ ⎦ ⎡ ⎤ − ⎣ ⎦

SLIDE 23

23

NPT-ensemble

acc( ) ( ) acc( ) ( )

n

N n n

N o

→ = →

( )

[ ]

( )

, exp exp s ;

N N N NPT

N V V PV U L β β ⎡ ⎤ ∝ − − ⎣ ⎦ s

Suppose we change the volume of the system

[ ]

( )

[ ]

( )

N N

exp exp s ; acc( ) acc( ) exp exp s ;

N n n n N

V

PV U L

n

n

V

PV U L β β β β ⎡ ⎤ − − → ⎣ ⎦ = → ⎡ ⎤ − − ⎣ ⎦

( ) ( ) ( )

{ }

exp exp

N n n

V

P V V U n U V β β ⎛ ⎞ ⎡ ⎤ ⎡ ⎤ = − − − − ⎜ ⎟ ⎣ ⎦ ⎣ ⎦ ⎝ ⎠

SLIDE 24

24

Algorithm: NPT

Randomly change the position of a particle
Randomly change the volume

SLIDE 25

25

SLIDE 26

26

SLIDE 27

27

SLIDE 28

28

NPT simulations

Equation of State of Lennard Jones System

SLIDE 29

29

Measured and Imposed Pressure

Imposed pressure P
Measured pressure <P> from virial

P = − < ∂F ∂V # $ % & ' ( >N,T = − dVV Ne−βPV

∫

dsNe−βU(sN ) ∂F ∂V # $ % & ' (

N,T

∫

dVV Ne−βPV

∫

dsNe−βU(sN )

∫

p(V) = exp −β F(V)+ PV

( )

* + ,

QNPT

QNPT = βP dV exp −β F(V)+ PV

( )

* + ,

∫

SLIDE 30

30

P = − < ∂F ∂V # $ % & ' ( >N,T = − dVV Ne−βPV

∫

dsNe−βU(sN ) ∂F ∂V # $ % & ' (

N,T

∫

dVV Ne−βPV

∫

dsNe−βU(sN )

∫

P = − βP Q(NPT) dV ∂F ∂V # $ % & ' (

N,T

exp −β F V

( )+ PV

( )

* + ,

∫

P = βP Q(NPT) dV exp −βPV

[ ]

β ∂exp −βF V

( )

* + ,

∂V

∫

SLIDE 31

31

Measured and Imposed Pressure

Partial integration
For V=0 and V=∞
Therefore,

P = βP Q(NPT) dV exp −βPV

[ ]

β ∂exp −βF V

( )

! " # $ ∂V

∫

P = βP Q(NPT) = dVPexp −β F V

( )+ PV

( )

! " # $

∫

= P

[ ]

∫ ∫

− =

b a b a b a

gdf fg fdg

exp −β F V

( )+ PV

( )

" # $ %= 0

SLIDE 32

32

Grand-canonical ensemble

What are the equilibrium conditions?

SLIDE 33

33

Grand-canonical ensemble

We impose:

– Temperature (T) – Chemical potential (µ) – Volume (V) – But NOT pressure

SLIDE 34

34

The ensemble of the total system

Here they are an ideal gas Here they interact

What is the statistical thermodynamics of this ensemble?

SLIDE 35

35

The ensemble: partition function

( )

3

ds exp s ; !

N N N NVT N

V Q U L N β ⎡ ⎤ = − ⎣ ⎦ Λ

∫

QMV0,NV ,T = V0 −V

( )

M −N

Λ3(M −N ) M − N

( )!

dsM −V

∫

exp −βU0 sM −N;L

( )

! " # $ V N Λ3N N! × dsN

∫

exp −βU sN;L

( )

! " # $

( ) ( )

( )

, , 3 3

ds exp s ; ! !

M N N N N MV NV T M N N

V V V Q U L M N N β

− −

− ⎡ ⎤ = − ⎣ ⎦ Λ − Λ

∫

SLIDE 36

36

QMV0,NV ,T = V0 −V

( )

M −N

Λ3(M −N ) M − N

( )!

V N Λ3N N! dsN

∫

exp −βU sN;L

( )

$ % & '

To get the Partition Function of this system, we have to sum over all possible number of particles

QMV0,N ,T = V0 −V

( )

M −N

Λ3(M −N ) M − N

( )!

V N Λ3N N! dsN

∫

exp −βU sN;L

( )

$ % & '

N=0 N=M

∑

Now let us take the following limits:

constant M M V V ρ → ∞⎫ = → ⎬ → ∞ ⎭

As the particles are an ideal gas in the big reservoir we have:

( )

3

ln

B

k T µ ρ = Λ

( )

3

exp ds exp s ; !

N N N N VT N N

N V Q U L N

µ

βµ β

=∞ =

⎡ ⎤ = − ⎣ ⎦ Λ

∑ ∫

SLIDE 37

37

µVT Ensemble

Partition function: Probability to find a particular configuration:

( )

3

exp , exp s ; !

N N N VT N

N V N V U L N

µ

βµ β ⎡ ⎤ ∝ − ⎣ ⎦ Λ s

Sample a particular configuration:

Change of the number of particles
Change of reduced coordinates

Acceptance rules ??

De Detailed balance

( )

3

exp ds exp s ; !

N N N N VT N N

N V Q U L N

µ

βµ β

=∞ =

⎡ ⎤ = − ⎣ ⎦ Λ

∑ ∫

SLIDE 38

38

Detailed balance

acc( ) ( ) ( ) ( ) acc( ) ( ) ( ) ( )

n

N n n

N n

n

N o
n

N o α α → × → = = → × → ( ) ( ) K o n K n

→

= → ( ) ( ) ( ) acc( ) K o n N o

n
n

α → = × → × →

n

( ) ( ) ( ) acc( ) K n

N n

n

n
α

→ = × → × →

SLIDE 39

39

µVT-ensemble

acc( ) ( ) acc( ) ( )

n

N n n

N o

→ = →

Suppose we change the position of a randomly selected particle

( )

N n 3 N

3

exp exp s ; acc( ) ! exp acc( ) exp s ; !

N N N N

N V U L

n

N N V n

U

L N βµ β βµ β ⎡ ⎤ − ⎣ ⎦ → Λ = → ⎡ ⎤ − ⎣ ⎦ Λ

( ) ( )

{ }

exp U n U β ⎡ ⎤ = − − ⎣ ⎦

( )

3

exp , exp s ; !

N N N VT N

N V N V U L N

µ

βµ β ⎡ ⎤ ∝ − ⎣ ⎦ Λ s

SLIDE 40

40

µVT-ensemble

acc( ) ( ) acc( ) ( )

n

N n n

N o

→ = →

Suppose we change the number of particles of the system

( )

1 N+1 3 3 N 3

exp 1 exp s ; 1 ! acc( ) exp acc( ) exp s ; !

N n N N

N

N V U L N

n

N V n

U

L N βµ β βµ β

+ +

+ ⎡ ⎤ − ⎣ ⎦ Λ + → = → ⎡ ⎤ − ⎣ ⎦ Λ

( )

3

exp , exp s ; !

N N N VT N

N V N V U L N

µ

βµ β ⎡ ⎤ ∝ − ⎣ ⎦ Λ s

( ) ( )

[ ]

3

exp exp 1 V U N βµ β = − Δ Λ +

SLIDE 41

41

SLIDE 42

42

SLIDE 43

43

Application: equation of state of Lennard-Jones

SLIDE 44

44

Application: adsorption in zeolites

SLIDE 45

45

Summary

Ensemble Constant (Imposed) Fluctuating (Measured) Function NVT N,V,T P βF=-lnQ(N,V,T) NPT N,P,T V βG=-lnQ(N,P,T) µVT µ,V,T N βΩ=-lnQ(µ,V,T)=-βPV

SLIDE 46

Practical Points

Boundaries
CPU saving methods
Reduced units
Long ranged forces

SLIDE 47

Boundary effects

In small systems, boundary effects are always large.
1000 atoms in a simple cubic crystal – 488 boundary

atoms.

1000000 atoms in a simple cubic crystal – still 6%

boundary atoms.

SLIDE 48

Periodic boundary conditions

SLIDE 49

ener(x,en) energies)

SLIDE 50

Energy evaluation costs!

The most time-consuming part of any simulation is

the evaluation of all the interactions between the molecules.

In general: N(N-1)/2 = O(N2)
But often, intermolecular forces have a short range:
Therefore, we do not have to consider interactions

with far- away atoms.

SLIDE 51

Saving CPU

Cell list

Verlet List

Han sur Lesse

SLIDE 52

Application: Lennard Jones potential

u LJ r

( ) = 4ε

σ r # $ % & ' (

12

− σ r # $ % & ' (

6

* + , ,

.

/ /

u r

( ) = u LJ r ( )

r ≤ r

c

r > r

c

" # $ u r

( ) = u LJ r ( ) − u LJ r

c

( )

r ≤ r

c

r > r

c

# $ %

The truncated and shifted Lennard-Jones potential
The truncated Lennard-Jones potential
The Lennard-Jones potential

SLIDE 53

ener(x,en) energies)

SLIDE 54

Phase diagrams of Lennard Jones fluids

SLIDE 55

Long ranged interactions

Long-ranged forces require special techniques.

– Coulomb interaction (1/r in 3D) – Dipolar interaction (1/r3 in 3D)

...and, in a different context:

– Interactions through elastic stresses (1/r in 3D) – Hydrodynamic interactions (1/r in 3D) –

SLIDE 56

Reduced units Example: Particles with mass m and pair potential: Unit of length: σ Unit of energy: ε Unit of time:

SLIDE 57

Beyond standard MC

Non Boltzmann Sampling – Lecture 2
Parallel tempering

More to Come Tomorrow (Daan Frenkel)

SLIDE 58

Ergodicity problems can occur, especially in glassy systems: biomolecules, molecular glasses, gels, etc. The solution: go to high temperature

Parallel tempering/Replica Exchange

E

phase space

E

phase space

low T high T

High barriers in energy landscape: difficult to sample Barriers effectively low: easy to sample

SLIDE 59

Parallel tempering/Replica Exchange

Simulate two systems simultaneously

e−β1U1(rN) e−β1U1(rN)e−β2U2(rN) e−β2U2(rN)

system 1 temperature T1 system 2 temperature T2 total Boltzmann weight:

SLIDE 60

Swap move

Allow two systems to swap

system 2 temperature T1 system 1 temperature T2

e−β1U2(rN) e−β1U2(rN)e−β2U1(rN) e−β2U1(rN)

total Boltzmann weight:

SLIDE 61

Acceptance rule

The ratio of the new Boltzmann factor over the old one is The swap acceptance ratio is

N(n) N(o) = e(β2−β1)[U2(rN)−U2(rN)]

1

acc(1 ↔ 2) = min

1, e(β2−β1)[U2(rN)−U1(rN)]⇥

SLIDE 62

Consider M replica’s in the NVT ensemble at a different temperature. A swap between two systems of different temperatures (Ti,Tj) is accepted if their potential energies are near.

More replicas

ther parameters can be used: Hamiltonian exchange

SLIDE 63

63

Questions Lunch …