Computational Statistical Mechanics

Computational Statistical Mechanics

2 agosto 2018

2

Indice

1 Termodinamica 5

1.1 Ensemble in termodinamica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.1.1 Quantità intensive ed estensive . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.2 Gas ideale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.3 Gas non ideale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.4 Le quantità fisiche di interesse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Meccanica statistica 11

2.1 Il teorema di ricorrenza di Poincarè . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.2 Teorema ergodico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3 Ensemble microcanonico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.3.1 Gas ideale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.4 Ensemble canonico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.4.1 Gas interagente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.5 Ensemble isobaro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.6 Ensemble grancanonico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.7 Gas ideale diatomico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.8 Sviluppo del viriale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.9 Calcolo dell’energia media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.9.1 Il problema dell’energia di un gas diatomico . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.10 Funzioni di distribuzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.10.1 Modello a sfere dure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.10.2 Misura sperimentale della funzione di correlazione . . . . . . . . . . . . . . . . . . . . . . 27

3 Metodo Monte Carlo 29

3.1 Stima di una funzione dipendente dalle osservabili . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.1.1 Funzione di un solo osservabile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.1.2 Funzione di due osservabili . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.2 Metodo Jack Knife . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.2.1 Stima dell’errore . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.2.2 Eliminazione del bias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.3 Catene di Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.3.1 Il problema del campionamento indipendente . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.3.2 Proprietà delle catene di Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.3.3 Campionamento markoviano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.3.4 Differenze tra campionamento indipendente e markoviano . . . . . . . . . . . . . . . . . . 39

3.4 Analisi a blocchi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.5 Calcolo delle funzioni di autocorrelazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3.6 Algoritmo di Metropolis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.6.1 Esposizione dell’algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.6.2 Modello di Ising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.6.3 Gas interagente - ensemble canonico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3.6.4 Gas interagente - ensemble isobaro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.6.5 Gas interangente - ensemble grancanonico . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

3.7 Widom - calcolo del potenziale chimico µ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.8 Metodo del reweighting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3.9 Metodo dei multi-istogrammi di Ferremberg-Swenden . . . . . . . . . . . . . . . . . . . . . . . . . 57

3.10 Umbrella Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

3.11 Simulated Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3

4 INDICE

3.12 Parallel Tempering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4 Dinamica Molecolare 65

4.1 Algoritmo di Verlet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

4.1.1 Verlet delle velocità . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

4.1.2 Verlet delle posizioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

4.1.3 Alcune problematiche . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

4.1.4 Scelta del passo temporale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

4.1.5 Conservazione della misura microcanonica . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

4.1.6 La presenza di una quantità conservata . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.1.7 Reversibilità . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

4.2 Moto vincolato . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

4.2.1 Multiple time steps algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

4.3 Metodo di Nosè-Hoover . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

4.3.1 Trattazione teorica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

4.3.2 Implementazione dell’algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

4.4 Equazioni differenziali stocastiche . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

4.4.1 Oscillatore armonico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

5 Appendice 85

5.1 Funzione Gamma di Eulero . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

Capitolo 1

Termodinamica

Primo principio. L’energia si conserva, cioè dE = δQ + δW,

dove dE è la variazione di energia interna, δQ è il calore dato al sistema mentre δW è il lavoro compiuto

sul sistema. Si noti che la convenzione circa il segno del lavoro è opposta rispetto a quella usuale, in cui il

lavoro era quello compiuto dal sistema. Consideriamo una trasformazione quasi statica, in cui applichiamo una

forza F ad un pistone di superficie S spostandolo di un tratto dx, da cui il lavoro è δW = F dx. Essendo una

trasformazione quasi statica si ha F = P S, da cui δW = P Sdx. Adesso, se dV < 0 è la variazione di volume

−dV −P

del gas si ha Sdx = da cui δW = dV . La variazione di energia interna risulta

−

dE = δQ P dV.

Si noti che se si fosse usata la convenzione opposta circa il lavoro si sarebbe ottenuto δW = P dV , ma il risultato

−

sarebbe stato lo stesso dato che in tale circostanza dE = δQ δW . Ora siamo interessati al calore scambiato

1

δQ, il quale in presenza di un trasformazione reversibile è δQ = T dS, con S pari all’entropia. Abbiamo dunque

−

dE = T dS P dV.

Si evince però che l’energia interna può cambiare anche variando il numero di particelle. La variazione al primo

ordine dovuta ad una variazione dN può essere integrata introducendo il potenziale chimico µ, ottenendo

−P . (1.1)

dE = T dS dV + µdN

| {z } {z }

|

W W

mec. chim.

1.1 Ensemble in termodinamica

La (1.1) ci dice che l’energia interna è funzione di (S, V, N ), in cui però l’entropia non è misurabile. Possiamo

però invertire la relazione in modo da ottenere un’espressione di S

1 P µ

− →

dS = dE + dV dN S = S(E, V, N ). (1.2)

T T T

Però, nuovamente negli esperimenti non abbiamo l’energia interna. Anche se questo non è vero nel caso delle

simulazioni, dove si ha una Hamiltoniana. In presenza della conservazione dell’energia possiamo studiare la

dinamica del nostro sistema mediante la legge di Newton. Questa circostanza coincide con il regime di validità

dell’ensemble microcanonico, che è utilizzato in particolare per lo studio della dinamica molecolare. A

partire da (1.2) si ottengono le relazioni

∂S P ∂S µ ∂S

1

−

= , = , = . (1.3)

T ∂E T ∂V T ∂N

T,V E,V E,V

In generale negli esperimenti non abbiamo una energia fissata, mentre lo sono (T, V, N ). Questa circostanza

corrisponde all’ensemble canonico. Vogliamo ottenere delle relazioni analoghe a (1.3) in tale contesto. A tal

fine −

dE = T dS P dV + µdN

− −

= d(T S) SdT P dV + µdN

− −SdT −

d(E T S) = P dV + µdN,

1 Dato che la variazione di energia è una funzione di stato possiamo sempre considerare un processo reversibile.

5

6 CAPITOLO 1. TERMODINAMICA

−

dove definiamo l’energia libera di Helmotz A = E T S. Dalla conoscenza di A si possono ricavare le

seguenti quantità termodinamiche ∂A ∂A ∂A

− −

S = , P = , µ = . (1.4)

∂T ∂V ∂N

V,N T,N T,V

L’ensemble canonico è particolarmente utile se abbiamo un sistema in contatto termico con una sorgente.

Ciononostante, in alcune circostanze anche le variabili (T, V, N ) non sono buone. Ad esempio, consideriamo un

2

liquido in un contenitore privo di coperchio. Il volume del sistema non è ben definito, mentre lo è la pressione.

Abbiamo dunque l’ensemble isobaro in cui le variabili buone sono (T, N, P ). Vogliamo ancora una volta

ottenere delle relazioni con le quantità termodinamiche a partire da un’energia. Seguiamo la stessa procedura

di prima da cui −SdT −

dA = P dV + µdN

−SdT −

= d(P V ) + V dP + µdN

−SdT

d(A + P V ) = + V dP + µdN,

dove definiamo l’energia libera di Gibbs G = A + P V . A partire da questa si ottengono

∂G ∂G

∂G

− , V = , µ = . (1.5)

S = ∂T ∂P ∂N

P,N T,N T,P

In alcuni casi, di interesse prettamente teorico, si hanno sistemi in cui il numero di particelle non è fissato

(ensemble grancanonico). L’unica eccezione in cui si ha sperimentalmente tale circostanza è in presenza di

sistemi in cui abbiamo la coesistenza tra due fasi diverse. Le variabili utili sono in questo caso (T, V, µ), dove si

è sostituito µ a N essendo quest’ultimo non fissato. A partire dall’energia libera di Helmotz si ha

−SdT −

dA = P dV + µdN

−SdT − −

= P dV + d(µN ) N dµ

− −SdT −

d(A µN ) = + V dP N dµ,

−

dove definiamo il potenziale grancanonico Ω = A µN . A partire da questo si ottengono

∂Ω ∂Ω ∂Ω

− − −

S = , P = , N = . (1.6)

∂T ∂V ∂µ

V,µ T,µ T,V

1.1.1 Quantità intensive ed estensive

Supponiamo di lavorare nel canonico, con variabili (N, T, V ), e vogliamo cambiarle rimanendo nel limite termo-

dinamico. Come farlo? Mandiamo sia N che V all’infinito mantenendo la densità ρ = N/V costante. Questo

risulta vero nel caso si escludano interazioni a lungo raggio, quali quella coulombiana, dato che in loro presenza

la trattazione si complica. Vogliamo adesso capire come le varie quantità cambino a seguito dell’aumento di N

e V a parità di ρ. Abbiamo che le seguenti quantità sono estensive, essendo proporzionali a N : E, S, A, G e

Ω. Mentre sono intensive: P , ρ, µ e T . Come possiamo sfruttare questa loro caratteristica? Ad esempio, dato

che G(P, T, N ) è estensiva allora il rapporto G/N sarà intensivo. Possiamo assumere ragionevolmente che tale

rapporto sarà una funzione f (P, T, N ) da cui, data l’intensività

G ∀λ,

= f (P, T, N ) = f (P, T, λN ) ,

N

che è verificata solo se f (P, T, N ) non dipende da N . Abbiamo dunque G = N f (P, T ), cioè sappiamo che la

dipendenza di G da N è data da un semplice fattore. Dalla (1.5) discende dunque

∂G

→

µ = = f (P, T ) G = µN.

∂N P,T

−

Dato che G = A + P V e A = E T S si ottengono −

A = µN P V (1.7)

−

E = µN P V + T S.

Si faccia attenzione al fatto che l’espressione di E ottenuta non può essere data da una banale integrazione della

(1.1). Per ottenerla si è infatti sfruttato l’estensività e l’intensività.

2 In presenza di un gas il volume è ben definito essendo questo confinato in una scatola chiusa, dato che altrimenti sfugge.

1.2. GAS IDEALE 7

1.2 Gas ideale

L’equazione dei gas ideali è →

P V = N kT P = ρkT, (1.8)

la quale ne enfatizza l’intensività. Come possiamo trattarlo lavorando nel canonico? Abbiamo dalla (1.7) e

(1.8) −

A = N (µ kT ),

in cui µ = µ(N, T, V ). Dato che è una quantità intensiva allora ∀λ,

µ(N, T, V ) = µ(λN, T, λV ) ,

implicando che la dipendenza di µ da N e V può essere derivante solo dalla densità ρ = N/V , da cui A =

−

N (µ(ρ, T ) kT ). Sfruttiamo la (1.4) ottenendo ∂µ

∂A ∂ρ ∂µ

− −

= µ kT + N

µ = = µ kT + ρ , (1.9)

∂N ∂ρ ∂N ∂ρ

T,V

da cui si evince che affinchè l’uguaglianza sia verificata bisogna avere

∂µ ρ

−kT →

+ ρ = 0 µ(ρ, T ) = kT ln ρ + h(T ) = kT ln ,

∂ρ ρ (T )

0

−kT

con h(T ) = ln ρ (T ). L’espressione cosı̀ ottenuta è più corretta essendo l’argomento del logaritmo adimen-

0

sionale. Sostituendola in A ρ

− 1 . (1.10)

A = N kT ln ρ (T )

0

Ci mancano delle informazioni su ρ (T ), che possiamo però ottenere dall’energia interna E = A + T S. Al fine

0

di calcolarla ricaviamo l’entropia A 1 ∂ρ (T )

∂A

0

− −

− = N kT ,

S = ∂T T ρ (T ) ∂T

0

da cui 1 ∂ρ (T )

0

2

E = N kT . (1.11)

ρ (T ) ∂T

0

L’energia interna è ragionevolmente indipendente dalla densità, dato che le particelle non interagiscono tra loro

essendo un gas perfetto. Ora, se conosciamo E(T ) possiamo ottenere a meno di una costante ρ (T ). Supponiamo

0

di ottenere dal principio di equipartizione E = αN kT da cui

1 ∂ρ (T )

0

2 α

→

αN kT = N kT ρ (T ) = CT ,

0

ρ (T ) ∂T

0

che sostituita in (1.10) ρ

−

A = N kT ln 1 .

α

CT

Il gas ideale risulta interessante essendo ragionevolmente il sistema a cui un arbitrario, sempre rimanendo nelle

→

interazioni a corto raggio, sistema interagente deve tendere per ρ 0. Sia data un’equazione di stato

P (N, T, V ) = P (ρ, T ),

→

se facciamo tendere ρ 0 abbiamo che il sistema tende a divenire un gas ideale, implicando

2 3 →

P (ρ, T ) = kT ρ + α (T )ρ + α (T )ρ + ... , ρ 0,

1 2

che costituisce lo sviluppo del viriale.

1.3 Gas non ideale

Alla luce di quanto detto nella sezione precedente vogliamo ottenere un’espressione relativa ai gas non ideali,

anche detti interagenti. L’equazione di stato risulta particolarmente importante dato che, come si mostrerà, se

essa è nota è possibile descrivere completamente il sistema.

8 CAPITOLO 1. TERMODINAMICA

Nel caso ideale si è trovata una funzione di stato P = ρkT , la quale deve essere il limite a cui l’equazione di un

3

→

gas non ideale deve tendere per ρ 0. Relativamente alla funzione di stato P = P (N, V, T ), da argomentazioni

generali legate alle proprietà delle quantità intensive si ha ⇒

P (N, V, T ) = P (λN, λV, T ) P = P (ρ, T ),

cioè funzione di N e V solo mediante la densità di particelle. Al fine di mostrare che data l’equazione di stato

si può trovare qualsiasi quantità , consideriamo l’energia libera di Helmotz nel canonico. A partire da (1.4) si

può ottenere, dalle derivate miste, la seguente identità

2 2

∂ A ∂µ ∂ A ∂P ∂µ

∂P

− ⇒ −

= , = = . (1.12)

∂N ∂N ∂V ∂V ∂V ∂N ∂N ∂V

V,T T N,T T V,T N,T

Queste relazioni sono del tutto generali, ma possiamo andare oltre sfruttando la proprietà di intensività della

pressione e del potenziale chimico da cui

∂ ∂ρ 1 ∂P

∂P

=

P (ρ, T ) =

∂N ∂ρ ∂N V ∂ρ

T T T

T ∂ρ N ∂P

∂ ∂µ

−

= ,

µ(ρ, T ) = 2

∂V ∂ρ ∂V V ∂ρ

T T T

T

che unite a (1.12) ∂P ∂µ

= ρ . (1.13)

∂ρ ∂ρ

T T

Abbiamo dunque una relazione che lega i tassi di variazione della pressione e del potenziale chimico rispetto

alla densità. Se è n