5.4 Monte Carlo Simulationen

(1)

5. Computersimulationen

1

5.1 Einleitung

2

5.2 Periodische Randbedinungen

3

5.3 Molekulardynamik Simulationen

4

5.4 Monte Carlo Simulationen

G. Kahl & F. Libisch (136) Statistische Physik II – Kapitel 5 2. Juni 2015 1 / 41

(2)

5.1 Einleitung

(a) Motivation

Computerunterstz¨ utzte Physik (’computational physics’), die sich methodisch der Computersimulationen bedient, wird zunehmend als drittes Standbein der Physik neben der theoretischen und der experimentellen Physik betrachtet

Begr¨ undung, Notwendigkeit:

• es existiert nur eine sehr kleine Zahl an exakt (d.h. analytisch) l¨ osbaren Modellen (z.B., Ising Modell f¨ ur D = 2); Computersimulationen erg¨ anzen/ersetzen/testen daher theoretische (N¨ aherungs-)Verfahren

• Simulationen setzen Standards, in dem sie ’quasi-exakte’ Ergebnisse liefern

• Simulationen erm¨ oglichen Behandlung komplexerer Probleme, als dies mit Hilfe theoretischer (N¨ aherungs-)Verfahren m¨ oglich ist

• Simulationen f¨ uhren oft ¨ uber Experimente und experimentelle M¨ oglichkeiten hinaus:

◦ extreme Bedingungen (Temperatur, Druck, ...)

◦ computational materials design: Abtasten hochdimensionaler

Parameterr¨ aume, freie Variation des Modells ⇒ Suche nach neuen

Systemen (Physik, Chemie, Pharmazie, Werkstoffwissenschaften, ...)

(3)

5.1 Einleitung

Aus Ref. [4.1].

• aus dem Vergleich zwischen den Ergebnissen aus Computersimulationen und

experimentellen Ergebnissen erh¨ alt man Informationen ¨ uber

Anwendbarkeit/Verl¨ aßlichkeit eines Modells f¨ ur ein System

• aus dem Vergleich zwischen den Ergebnissen aus Computersimulationen und theoretischen Vorhersagen erh¨ alt man Informationen ¨ uber die Qualit¨ at und Anwendbarkeit des theoretischen (N¨ aherungs-)Methode

(4)

5.1 Einleitung

(b) (Kurze) Geschichte von Computersimulationen

• einfache Modelle fl¨ ussiger Strukturen:

Gelatinekugeln, Eisenkugeln – J.D. Bernal (∼ 1930)

• erste Computersimulationen von N. Metropolis et al. (1953) – Los Alamos: Monte Carlo Simulationen an harten Teilchen (Scheiben)

• erste Molekulardynamik Simulationen von Alder und Wainwright (1956) an harten Teilchen

Trajektorien harter Scheiben: Festk¨orper (links) und fl¨ussige Phase (rechts); aus Ref. [4.2].

• Monte Carlo Simulationen an einem Lennard-Jones System von Wood und Parker (1957)

• erste Molekulardynamik Simulationen an Lennard-Jones Systemen von Rahman (1964)

• Computersimulationen beginnen sich ∼ 1967 durchzusetzen

(5)

5.1 Einleitung

in der Folge:

• Simulationen an Systemen mit zunehmender Komplexit¨ at werden durchgef¨ uhrt: Molek¨ ule (ab 1968), Wasser (ab 1969), ..., ab initio Simulationen (unter Ber¨ ucksichtigung atomarer und elektronischer Freiheitsgrade), ...

weitere Anwendungen:

• Nichtgleichgewichtssimulationen

• Quantenfl¨ ussigkeiten

• biologische Systeme

• ...

andere Bereiche, wo Simulationen eine große Rolle spielen:

• Katalyse

• computational materials design

• Aerodynamik

• Methode der Finiten Elemente

• ...

(6)

5.1 Einleitung

(c) Grundlagen

• Festlegung des Modells durch (i) die Hamilton-Funktion H (ii) die Geometrie des Systems, (iii) ...

• Wahl des Ensembles legt den Phasenraum Γ = {z} = {p

^N

, r

^N

} fest

◦ N = const., V = const., E = const. ⇒ mikrokanonisches Ensemble

◦ N = const., V = const., T = const. ⇒ kanonisches Ensemble

◦ ...

• Ziel von Computersimulationen ist die Berechnung physikalischer Eigenschaften (Thermodynamik, strukturelle Eigenschaften,

elektronische Eigenschaften, Transporteigenschaften, ...) mit Hilfe einer geeigneten Mittelwertbildung im Phasenraum

sei A = A(z) eine Observable, dann gibt es (in der klassischen Statistischen Physik) zwei M¨ oglichkeiten der Mittelwertbildung;

(i) Ensemblemittelwert mit Hilfe jener Verteilungsfunktion f (z), die das Ensemble charakterisiert

Z = Z

Γ

dzf (z) Zustandssumme hAi

e

= 1

Z Z

Γ

dzA(z)f (z) Ensemblemittelwert

(7)

5.1 Einleitung

(ii) Zeitmittelwert: Konstruktion einer Trajektorie z = z(t) durch den Phasenraum, Mittelwertbildung erfolgt dann entlang dieser Trajektorie

¨

uber ein Zeitintervall T , mit (idealerweise) T → ∞

hAi

t

= lim

T→∞

1 T

Z

T

0

dt

⁰

A[z(t

⁰

)] Zeitmittelwert zentrale Frage: gilt die Gleichheit dieser Mittelwerte, also

hAi

t

= hAi

e

f¨ uhrt zur Ergodenhypothese in der Praxis

◦ k¨ onnen nur endliche Bereiche im Phasenraum Γ abgetastet werden

◦ kann nur eine endliche Beobachtungsdauer T realisiert werden praktische Umsetzung:

(i) stochastische Methode:

es wird (fast immer) nur der Konfigurationsanteil von H ben¨ otigt;

’Bewegung’ von einer Konfiguration zur anderen im

Konfigurationsraum Γ mit Hilfe einer ’probabilistischen Mechanik’ ⇒ Monte Carlo Simulation

(8)

5.1 Einleitung

(ii) deterministische Methode:

Verwendung der intrinsischen Dynamik (Newton, Lagrange, Hamilton) des Systems:

es werden die Bewegungsgleichungen aufgestellt und numerisch integriert

t = t

₀

z(t

₀

) = {p

₁

(t

₀

), · · · , p

_N

(t

₀

); r

₁

(t

₀

), · · · , r

_N

(t

₀

)}

numerische Integration w w

t > t

0

z(t) = {p

1

(t), · · · , p

N

(t ); r

1

(t), · · · , r

N

(t)}

⇒ Molekulardynamik Simulation

(9)

5.1 Einleitung

(d) Mittelwerte und Fluktuationen

Aussage aus der LVA ”Statistische Physik I”:

... im thermodynamischen Grenzwert (also, z.B., N → ∞ und V → ∞ mit

% = N/V = const.) stimmen die in verschiedenen Ensembles ermittelten Ergebnisse f¨ ur eine Observable (die mit Hilfe von Mittelwerten berechnet wurden) ¨ uberein ...

in der Praxis von Computersimulationen k¨ onnen aber nur Ensembles mit endlicher Gr¨ oße (z.B. Teilchenzahl, Volumen, ...) betrachtet werden, daher ist obige Aussage nicht mehr richtig

aber: es k¨ onnen Beziehungen zwischen den Mittelwerten in den verschiedenen Ensembles hergeleitet werden (vgl., z.B., Ref. [4.1])

Mittelwerte:

bei der Berechnung physikalischer Observablen ¨ uber Zeit- oder

Ensemblemittelwerte werden grundlegende thermodynamische Relationen ausgen¨ utzt

Beispiele:

• Temperatur: sei H = K + V , dann ist hKi = 3

2 Nk

B

T ⇒ (kinetische) Temperatur

(10)

5.1 Einleitung

• Druck: sei F

i

= −∇

i

V (r

1

, · · · , r

N

) die Kraft, die auf das Teilchen mit Index i wirkt, dann gilt

PV = Nk

B

T − 1 3 h

N

X

i=1

r

i

· F

i

i (Virial−)Druck

• statische Korrelationsfunktionen:

Γ(r, r

⁰

) = h (%(r) − %) (%(r

⁰

) − %) i geben Auskunft ¨ uber die statische Struktur

• dynamische Korrelationsfunktionen:

Geschwindigkeitsautorkorrelationsfunktion Ψ(t) Ψ(t ) = hv

i

(t ) · v

_i

(t

₀

)i ⇒ D = 1

3 Z

∞

0

dtΨ(t )

Ψ(t ) gibt Auskunft ¨ uber die dynamischen Eigenschaften des Systems

und, im weiteren, ¨ uber Transporteigenschaften (D ist hier der

Diffusionskoeffizient)

(11)

5.1 Einleitung

Fluktuationen:

• die Wahl des Systems, des Ensembles und der Randbedingungen legen (rigorose) Erhaltungsgr¨ oßen fest (Noether-Theorem)

Beispiele f¨ ur Erhaltungsgr¨ oßen:

◦ Gesamtenergie im mikrokanonischen Ensemble

◦ Gesamtimpuls in Ensembles mit fixer Teilchenzahl und periodischen Randbedingungen

◦ . . .

Erhaltungsgr¨ oßen sind wichtig, um Aufschluß ¨ uber die numerische Genauigkeit der Implementierung der Computersimulation zu erhalten

• andere Gr¨ oßen d¨ urfen/sollen fluktuieren:

sei A eine Observable und δA = A − hAi mit h(δA)

²

i = hA

²

i − hAi

²

es gilt, zum Beispiel (Beweise sind in Ref. [4.1] angegeben):

◦

h(δV)

²

i

NVE

= h(δK)

²

i

NVE

= 3 2 Nk

B

T

²

1 − 3Nk

B

2C

V

◦

h(δV)

²

i

NPT

= Vk

B

T κ

T

◦

h(δN)

²

i

µVT

= k

B

T N

²

V κ

T

◦ etc.

(12)

5.1 Einleitung

(e) Nachteile und Schwachpunkte von Computersimulationen

• begrenzte Ensemblegr¨ oße:

10

⁵

∼ 10

⁷

Teilchen im Vergleich zu ∼ 10

²³

Teilchen in realen Systemen

• verschiedene Zeitskalen innerhalb eines Modells/Problems

z.B.: Elektronen vs. Kerne, Molek¨ ule vs. Bestandteile der Molek¨ ule, ...

• begrenzte Computerleistung

(sowohl bez¨ uglich Rechenzeit als auch Speicherplatz)

• begrenzte numerische Genauigkeit

• Mittelungsprozeß

• ...

(13)

5.2 Periodische Randbedinungen

5.2 Periodische Randbedingungen

Hauptproblem von Computersimulationen: endliche Systemgr¨ oße (Teilchenzahl N) Systemgr¨ oße durch Rechenleistung begrenzt, da die Berechnung von Kr¨ aften oder Energien (im allgemeinen) wie N

²

skaliert

die (relativ) kleine Systemgr¨ oße f¨ uhrt zu erheblichen Einschr¨ ankungen:

• begrenzte Statistik

• Oberfl¨ acheneffekte

in einem W¨ urfel mit 10 × 10 × 10 = 1000 Teilchen befinden sich 488 (!) Teilchen an den sechs Oberfl¨ achen; ohne zus¨ atzliche Maßnahmen wirken auf die Teilchen an den Oberfl¨ achen (bzw. nahe den Oberfl¨ achen) andere Kr¨ afte als auf jene Teilchen, die sich im Inneren des W¨ urfels befinden

• eingeschr¨ ankte Beschreibung kritischer Ph¨ anomene

• ...

M¨ ogliche Randbedingungen:

• harte, undurchdringliche W¨ ande

• periodische Randbedingungen (Born & von Karman, 1912)

• spezielle, problemspezifische Randbedingungen (Oberfl¨ achen, Schichtstrukturen, ...)

(14)

periodische Randbedingungen

Schematische Darstellung einer (zwei dimensionalen)Simulationszelle mit periodischen Randbedingungen; aus Ref. [4.3].

• die Simulationszelle wird ’vervielf¨ altigt’; somit wird ein periodisches Gitter aufgebaut;

die Simulationszelle muß nicht notwendigerweise kubisch sein!

• die Bewegungen der Teilchen in der zentralen Zelle und in den Spiegelzellen verlaufen parallel

• der Austritt eines Teilchens aus der (zentralen) Simulationszelle entspricht dem

Eintritt seines Spiegelteilchens aus einer Kopie der Simulationszelle ⇒ die

Teilchenzahl und somit die Teilchendichte bleibt erhalten

(15)

• es werden nur die Daten (Orte, Geschwindigkeiten, ...) der Teilchen in der zentralen Simulationszelle abgespeichert

• f¨ ur die Reichweite des Potentials wird ein Abschneideradius r

c

(mit r

c

< L/2) eingef¨ uhrt

es gilt die ’minimum image convention’: ein herausgegriffenes Teilchen ist mit all jenen Teilchen in Wechselwirkung, die sich zu einem gegebenen Zeitpunkt in einer (fiktiven) Zelle befinden, die (maximal) so groß ist wie die Simulationszelle und um dieses Teilchen zentriert ist

• die Einf¨ uhrung eines Abschneideradius erfordert spezielle Methoden (z.B.

Ewald-Summen) bei langreichweitigen Potentialen (z.B. Coulomb-Wechselwirkung, Dipolwechselwirkung, ...)

Fragen, Probleme:

• ist das ’kleine’ System mit periodischen Randbedinungen ¨ aquivalent dem großen (’realistischen’) System?

• was passiert bei Ph¨ anomenen, wo langreicheweitige Korrelationen relevant sind (kritische Ph¨ anomene, ...)

(16)

5.3 Molekulardynamik Simulationen

(a) Molekulardynamik Simulationen im mikrokanonischen Ensemble System ist gegeben durch:

◦ Teilchenzahl N (fest)

◦ Volumen V (fest); z.B. W¨ urfel

◦ Annahme: es existieren keine ¨ außeren Kr¨ afte (wie, z.B., Schwerkraft)

◦ die Hamiltonfunktion ist gegeben durch

H = H(r

1

, · · · , r

N

; p

1

, · · · , p

N

) = K(p

^N

) + V(r

^N

) mit K =

N

X

i=1

p

²_i

2m Bewegungsgleichungen:

(i) Hamilton dr

_i

(t)

dt = ∂H

∂p

i

d p

_i

(t )

dt = − ∂H

∂r

i

i = 1, · · · , N (ii) Newton

d²r_i(t)

dt²

= p

_i

(t)/m

dpi(t)

dt

= −∇

i

V = F

i

)

m d

²

r

i

(t )

dt

²

= F

i

= −∇

i

V(r

^N

)

(17)

Bemerkungen:

• es existieren je nach Randbedinungen Erhaltungsgr¨ oßen

• die Bewegungsgleichungen sind formal invariant unter Zeitumkehr

• besondere Algorithmen werden bei Systemen mit unstetigen Potentialen ben¨ otigt

numerische Integration der Bewegungsgleichungen:

• Diskretisierung der Zeit: t ⇒ i∆t , i = 1, · · · , N

T

∆t ist das Zeitinkrement – wird systemspezifisch gew¨ ahlt (z.B. Argon: ∆t ∼ 10

⁻¹⁵

s)

Simulationsdauer N

T

– wird systemspezifisch gew¨ ahlt (z.B. ’einfache’ Fl¨ ussigkeiten: N

T

∼ 10

⁹

)

somit ergibt sich eine makroskopische Simulationsdauer von ∼ 10

⁻⁹

s

(18)

• durch die Diskretisierung in der Zeit werden die Differentialgleichungen der Bewegungsgleichungen zu Differenzengleichungen f¨ ur die

diskretisierten Orte, Geschwindigkeiten, Beschleunigungen, ... der Teilchen:

r

k

(t) ⇒ r

k

(t = i∆t ) dr

k

(t)

dt = v

k

(t) ⇒ v

k

(t = i∆t ) d

²

r

_k

(t )

dt

²

= a

_k

(t) ⇒ a

_k

(t = i∆t) k = 1, · · · , N numerische L¨ osungsstrategien der Differenzengleichungen:

es handelt sich um 3N gekoppelte Differenzengleichungen, die zu vorgegebenen Anfangsbedingungen

{r

^N

(t = 0); v

^N

(t = 0)} bzw. {r

^N

(t = 0); r

^N

(t = ∆t)}

numerisch gel¨ ost werden m¨ ussen

in der Literatur existieren zahlreiche numerische Verfahren, um dieses Problem zu l¨ osen (vgl. Ref. 4.1, 4.4])

Qualit¨ atskriterien: Stabilit¨ at, Genauigkeit, numerischer Aufwand, ...

(19)

Beispiele

(a) Verlet-Algorithmus (Teilchenindex weggelassen)

r(t + ∆t ) ∼ = r(t) + ∆t v(t) + (1/2) (∆t )

²

a(t) + · · · r(t − ∆t) ∼ = r(t ) − ∆t v(t ) + (1/2) (∆t)

²

a(t ) + · · · somit

r(t + ∆t) ∼ = 2r(t ) − r(t − ∆t ) + (∆t )

²

a(t) + O (∆t)

⁴

wobei

a(t) = d

²

r(t) dt

²

= − 1

m ∇V(r

^N

) weiters

v(t) ∼ = 1

2∆t [r(t + ∆t ) − r(t − ∆t)] + O (∆t)

²

(20)

(b) ’leap-frog’-Algorithmus

r(t + ∆t) ∼ = r(t ) + ∆t v(t + ∆t /2) + · · · v(t + ∆t/2) ∼ = v(t − ∆t /2) + ∆t a(t ) + · · · mit

v(t ) ∼ = 1

2 [v(t + ∆t /2) − v(t − ∆t /2)]

die Geschwindigkeiten werden zu ’halbzahligen’ Zeitschritten berechnet, w¨ ahrend Orte und Beschleunigungen zu ’ganzzahligen’

Zeitschritten ermittelt werden

Schematische Darstellung der Berechnung der Teilchenorte aus den Geschwindigkeiten und den Beschleunigungen beimVerlet(oben) und beim’leap-frog’(unten) Algorithmus;

aus Ref. [4.1].

(21)

Bemerkungen:

◦ die Wahl des Integrationsalgorithmus ist sehr heikel;

dabei m¨ ussen konzeptuell bessere Algorithmen nicht notwendigerweise bessere Ergebnisse liefern; oft liefert bereits der Verlet-Algorithmus sehr gute Ergebnisse

◦ eine formale Zeitumkehrsymmetrie ist in der Praxis nicht realisiert

◦ Lyapunov-Instabilit¨ at

(22)

Umsetzung einer Molekulardynamik Simulation:

(i) Anfangsbedingungen f¨ ur die Teilchenkoordinaten und/oder Teilchen Geschwindigkeiten werden f¨ ur t = 0 vorgegeben, also

{r

^N

(t = 0); v

^N

(t = 0)} bzw. {r

^N

(t = 0); r

^N

(t = ∆t)}

(ii) Berechnung von V (r

^N

) f¨ ur t = 0 und t = ∆t;

daraus erfolgt die (numerische) Berechnung der Kr¨ afte F

i

zum Zeitpunkt ∆t ¨ uber F

_i

= −∇

i

V

(iii) Berechnung von r

_k

(∆t + ∆t ) = r

_k

(2∆t) aus r

k

(∆t + ∆t ) = 2r

k

(∆t) − r

k

(0) + (∆t)

²

1 m F

k

(∆t ) (iv) Abspeichern und Berechnung verschiedener Observablen:

{r

^N

(2∆t)}, {v

^N

(2∆t )}, {a

^N

(2∆t )}, ... bzw. A[p

^N

(2∆t), r

^N

(2∆t)]

(v) Vorr¨ ucken der Zeit um ∆t

(vi) zur¨ uck zu Schritt (ii); etc.

(23)

Ablauf einer Molekulardynamik Simulation:

• Initialisierungsphase: Festlegung der Anfangskonfiguration bzw. der Anfangsgeschwindigkeiten (z.B. ¨ uber eine Maxwell-Boltzmann Verteilung, die durch die Temperatur T charakterisiert ist)

• Aquilibrierungsphase: ¨ System wird ’beobachtet’ (z.B. zur ¨ Uberpr¨ ufung der Erhaltungsgr¨ oßen), ...

hier besteht die M¨ oglichkeit, korrigierend einzugreifen (Modifikation von ∆t , Reskalierung der Geschwindigkeiten, ...)

• Produktionsphase: Zeit wird auf t = 0 zur¨ uckgesetzt, das System wird sich selbst ¨ uberlassen;

die Koordinaten und die Geschwindigkeiten der Teilchen, sowie weitere relevante Gr¨ oßen werden in regelm¨ aßigen Abst¨ anden abgespeichert;

sei N

_T

Zahl der Integrationsschritte, so gilt f¨ ur den (Zeit-)Mittelwert einer Observablen A

hAi = 1 N

T

∆t

NT

X

i=1

A

p

^N

(i∆t), r

^N

(i∆t)

∆t

(24)

Verlauf derkinetischen (oben), der potentiellen (Mitte)und dergesamten Energie (unten)alsFunktionen der Zeit (d.h. Molekulardynamikschritte)in einer Molekulardynamiksimulation im mikrokanonischen Ensemble; das simulierte System ist ein Lennard-Jones System zu den angegebenen

(reduzierten) Dichte- und Temperaturwerten; aus Ref. [4.5].

(25)

(b) Molekulardynamik Simulationen im kanonischen Ensemble (i) Andersen Thermostat:

die Teilchen eines (fiktiven) W¨ armebades (durch die Temperatur T charakterisiert) wechselwirken ¨ uber zuf¨ allig generierte ”St¨ oße” mit den Teilchen des Systems

◦ in regelm¨ aßigen Intervallen wird die Geschwindigkeit eines zuf¨ allig herausgegriffenen Teilchens gem¨ aß einer Maxwell-Boltzmann Verteilung ver¨ andert; dies entspricht einer fiktiven Kollision mit einem Teilchen aus dem W¨ armebad

◦ durch die Kollision wird die Energie des Systems ver¨ andert, das System wird somit auf eine andere Energiefl¨ ache ’verschoben’; zwischen zwei dieser Kollisionen wird das System im Rahmen eines mikrokanonischen Ensembles simuliert; das System ’bewegt’ sich somit auf einer Folge von Fl¨ achen konstanter Energie

◦ dadurch wird der Phasenraum so abgetastet, wie dies im Rahmen eines kanonischen Ensembles geschehen w¨ urde

• trotz des sehr einfachen Konzepts, sind einige technische Schwierigkeiten zu

¨ uberwinden (z.B. Wahl der Kollisionsrate, ..)

• das Verfahren ist zur Berechnung mancher thermodynamischer oder

struktureller Gr¨ oßen (z.B. dynamische Korrelationsfunktionen) nicht geeignet

(26)

(ii) Nos´ e-Hoover Thermostat:

• es werden eine zus¨ atzlichen Koordinate s und ihr konjugierter Impuls p

_s

eingef¨ uhrt

f¨ ur den erweiterten Satz von Variablen (virtuelle Variablen) wird eine Molekulardynamik Simulation im mikrokanonischen Ensemble durchgef¨ uhrt

• Hamiltonfunktion H

Nose

=

N

X

i=1

p

²_i

2ms

²

+ V (r

^N

) + p

²_s

2Q + g ln s β

• Q ... effektive Masse

g ... ein vorerst nicht n¨ aher definierter Parameter (sh. sp¨ ater)

• Bewegungsgleichungen f¨ ur die virtuellen Variablen d r

_i

dt = ∂H

_Nose

∂p

i

= p

_i

ms

²

dp

_i

dt = − ∂H

_Nose

∂r

i

= − ∂V

∂r

i

ds

dt = ∂H

Nose

∂p

s

= p

s

Q

dp

s

dt = − ∂H

Nose

∂s = 1

s X

i

p

²_i

ms

²

− g

β

!

(27)

• Berechnung der mikrokanonischen Zustandssumme f¨ ur das erweiterte Ensemble (basierend auf H

Nose

)

Z

N

∼ Z

ds dp

s

dr

^N

d p

^N

δ(H

Nose

− E ) =

= Z

ds dp

s

dr

^N

(dp

^0N

s

^3N

| {z }

p⁰=~p/s

)δ X

^N

i=1

p

⁰²_i

2m + V (r

^N

)

| {z }

H(~p^0N,~r^N)

+ p

²_s

2Q + g ln s β − E

mit δ[f (s)] = δ(s − s

₀

)/|f

⁰

(s

₀

)|, wobei s

₀

(einfache) Nullstelle von f (s) ist, gilt somit

Z

_N

∼ Z

dp

_s

dp

^0N

d r

^N

Z

ds β

g s

^(3N+1)

×

× δ h

s − exp n

− β g

H(p

^0N

, r

^N

) − p

²_s

2Q − E

oi

(28)

Z

_N

∼ β

g exp[E (3N + 1)/g ] Z

dp

_s

exp

−β 3N + 1 g

p

_s²

2Q

×

× Z

dr

^N

dp

^0N

exp

− β

g (3N + 1)H(p

^0N

, r

^N

)

= C

Z

dr

^N

dp

^0N

exp

−βH(p

^0N

, r

^N

)

wenn g = 3N + 1 der letzte Ausdruck entspricht aber einer kanonischen Zustandssumme f¨ ur die Hamiltonfunktion H(p

^0N

, r

^N

), also in den Variablen p

^0N

und r

^N

zusammenfassend ist die

mikrokanonische Zustandssumme

[mit Hamiltonfunktion H

Nose

(p

^N

, r

^N

, s, p

s

)]

(im wesentlichen) gleich der

kanonischen Zustandssumme

[mit Hamiltonfunktion H(p

^0N

, r

^N

)]

(29)

virtuelle Koordinaten reelle Koordinaten (mikrokanonisches Ensemble) (kanonisches Ensemble)

p p

⁰

= p/s

r r

⁰

= r

s s

⁰

= s

∆t ∆t

⁰

= ∆t

¹_s

• die Bewegungsgleichungen f¨ ur die virtuelle Koordinaten sind auf Folie 25 zusammengestellt

die Bewegungsgleichungen f¨ ur reellen Koordinaten lassen sich aus den obigen Relationen herleiten

• Wahl von Q ⇔ Kopplung an den Thermostat

(30)

5.4 Monte Carlo Simulationen

(a) Grundlagen, Monte Carlo Simulationen im kanonischen Ensemble

• in Monte Carlo Simulationen werden zur Berechnung von Observablen Ensemblemittelwerte ausgewertet; zum Beispiel im kanonischen Ensemble

hAi

_e

= R

Γ

dzA(z) exp[−β H(z)]

R

Γ

dz exp[−βH(z)] ∼ P

NK

i=0

A(z

i

) exp[−βH(z

i

)]

P

N_K

i=0

exp[−βH(z

i

)]

• die z

_i

sind ausgew¨ ahlte S¨ atze von verallgemeinerten Koordinaten

• zentrale Frage: wie werden die {z

_i

} ausgew¨ ahlt?

• in Monte Carlo Simulationen werden die {z

i

} bem¨ aß einer vorgegebenen Wahrscheinlichkeitsverteilung {p

i

} ausgew¨ ahlt;

somit gilt f¨ ur obigen N¨ aherungsausdruck des Ensemblemittelwertes hAi

e

∼

P

N_K

i=1

A(z

i

)p

i−1

exp[−βH(z

i

)]

P

NK

i=1

p

i−1

exp[−β H(z

i

)]

d.h. es werden Zust¨ ande im Phasenraum zuf¨ allig ausgew¨ ahlt, wobei die

Wahrscheinlichkeit der Auswahl konsistent ist mit der vorgegebenen

Wahrscheinlichkeitsverteilung der {p

i

}

(31)

• Beispiele f¨ ur die {p

_i

}:

(i) alle p

i

sind gleich: dann werden auch

unphysikalische/unwahrscheinliche Konfigurationen z

i

(z.B. mit

¨

uberlappenden Teilchen) mitber¨ ucksichtigt, deren Boltzmannfaktor exp[−βH(z

i

)] keinen relevanten Beitrag zu hAi liefert

(ii) ’importance sampling’: die p

i

werden gem¨ aß der Verteilungsfunktion f (z) ausgew¨ ahlt, die das Ensemble charakterisiert, also p

i

∝ f (z

i

) konkret f¨ ur das kanonische Ensemble gilt somit

p

i

∝ exp[−βH(z

i

)]

und daher

hAi

e

∼ P

N_K

i=1

A(z

i

)p

⁻¹_i

exp[−βH(z

i

)]

P

N_K

i=1

p

_i⁻¹

exp[−βH(z

i

)] ∼ 1 N

K

N_K

X

i=1

A(z

i

) (1)

• zentrale Frage: wie kann man die Zust¨ ande z

i

a priori so ausw¨ ahlen, sodaß ihre Wahrscheinlichkeiten p

i

(im kanonischen Ensemble) gem¨ aß exp[−βH(z

i

)] verteilt sind? ⇒ Realisierung in einem Markov Prozeß

(32)

• Konstruktion einer Folge von Zust¨ anden z

_i

:

◦ Markov Prozeß, Markov Kette

? ein Markov Prozeß erzeugt – ausgehend von einem Zustand z

i

– einen neuen Zustand z

j

; dieser Prozeß ist zufallsbehaftet, d.h., ausgehend von einem z

i

, wird nicht immer derselbe Zustand z

j

erzeugt

? die Wahrscheinlichkeit, daß aus z

i

der Zustand z

j

erzeugt wird, wird mit P(i → j) bezeichnet

? in einem Markov Prozeß ist (i) diese Wahrscheinlichkeit von der Zeit unabh¨ angig und (ii) sollte nur von den Zust¨ anden z

i

und z

j

abh¨ angen

? weiters gilt: (i) P(i → i ) kann von null verschieden sein;

(ii) offensichtlich ist

X

j

P(i → j) = 1

? die wiederholte Anwendung eines Markov Prozesses f¨ uhrt zu einer Markov Kette

? in Monte Carlo Simulationen werden durch speziell konstruierte Markov Ketten Folgen von Zust¨ anden erzeugt

? damit die Wahrscheinlichkeiten dieser Folge von Zust¨ anden ’im Gleichgewicht’ nach einer kanonischen Verteilungsfunktion verteilt sind, m¨ ussen zwei weitere Bedingungen eingef¨ uhrt werden:

’ergodicity’ und ’detailed balance’

(33)

◦ ’ergodicity’

ausgehend von einem beliebigen Zustand des Systems, soll jeder andere Zustand des Systems mit Hilfe eines Markov Prozesses erreicht werden k¨ onnen (egal wie lange das dauern kann)

◦ ’detailed balance’

diese Bedingung garantiert, daß im Gleichgewicht die Wahrscheinlichkeit der Zust¨ ande gem¨ aß der kanonischen Verteilungsfunktion verteilt sind

Gleichgewicht bedeutet, daß die ¨ Ubergangsrate in einen Zustand und aus diesem Zustand (Index j) gleich sind;

Ausgangspunkt ist ein Spezialfall der s.g. ’master’-Gleichung:

X

i

p

i

P(i → j) = X

i

p

j

P(j → i ) bzw. X

i

p

i

P(i → j) = p

j

ohne Beweis: diese Beziehung alleine reicht nicht aus, um zu garantieren, daß alle p

j

gem¨ aß der vorgegebenen Verteilungsfunktion verteilt sind; diese Forderung kann durch folgende zus¨ atzliche Annahme gew¨ ahrleistet werden:

p

j

P(j → i) = p

i

P(i → j)

⁰

detailed balance

⁰

(34)

wir betrachten nun ein kanonisches Ensemble und sei E

i

= H(z

i

) und E

j

= H(z

j

), dann gilt offensichtlich

P(j → i) P(i → j) = p

i

p

j

= exp[−β(E

i

− E

j

)] (2) Achtung: hier geht die Verteilungsfunktion des Ensembles ein!

erf¨ ullen die P(i → j) diese Relation, dann ist garantiert, daß die Wahrscheinlichkeiten der Zust¨ ande z

i

im Gleichgewicht nach den p

i

verteilt sind

allerdings l¨ aßt Relation (2) zwischen den P(i → j) noch einen relativ hohen Grad an Freiheit f¨ ur den Algorithmus zu

konkretes Beispiel

◦ Metropolis Algorithmus

wurde von Metropolis und Mitarbeitern ∼ 1953 eingef¨ uhrt Erzeugung der Zust¨ ande wird in zwei Schritte geteilt:

(i) ausgehend von Zustand z

i

: Auswahl eines neuen Zustandes z

j

mit Wahrscheinlichkeit g (i → j) – ’selection probability’;

(ii) Annahme des Zustandes z

j

mit Wahrscheinlichkeit A(i → j) –

’acceptance ratio’

somit

P(j → i )

P(i → j) = g (j → i)A(j → i)

g (i → j)A(i → j)

(35)

(i) konkret w¨ ahlt man aus rechentechnischen Gr¨ unden g(i → j) = const.; somit

P(j → i)

P(i → j) = A(j → i)

A(i → j) = exp[−β(E

i

− E

j

)] (3) (ii) damit die Akzeptanzrate aus Effizienzgr¨ unden m¨ oglichst groß ist, w¨ ahlt man sinnvollerweise folgende L¨ osung f¨ ur die A(i → j):

A(i → j) =

exp[−β(E

j

− E

i

)] wenn E

j

− E

i

> 0

1 sonst

aber: jede andere Akzeptanzrate, die (3) erf¨ ullt, ist de facto genauso gut

• Monte Carlo Simulation im kanonischen Ensemble (Fl¨ ussigkeit)

◦ da die zeitliche Entwicklung bei Monte Carlo Simulationen offensichtlich keine Rolle spielt, werden in den z

i

die Impulse weggelassen; somit z

i

= {r

^N

}

i

und E

i

= V({r

^N

}

i

)

◦ ausgehend von einem Zustand {r

^N

}

i

wird ein Zustand {r

^N

}

j

erzeugt;

dies geschieht, in dem man die Position eines zuf¨ allig herausgegriffenen Teilchens, r

α;i

um einen Vektor δr

α;i

verschiebt; der Betrag von δr

α;i

ist durch das ’maximal displacement’ begrenzt, seine Richtung ist durch Zufallszahlen gegeben

(36)

◦ Berechnung von E

i

= V({r

^N

}

i

) und E

j

= V({r

^N

}

j

)

◦ der neue Zustand {r

^N

}

j

wird mit der Wahrscheinlichkeit

A(i → j) =

exp[−β(E

j

− E

i

)] wenn E

j

− E

i

> 0

1 sonst

= min {1, exp[−β(E

j

− E

i

)]}

angenommen oder verworfen:

(i) wird der Zustand {r

^N

}

j

angenommen, so kann aus ihm ein weiterer Zustand, {r

^N

}

k

, erzeugt werden

(ii) wird der Zustand {r

^N

}

j

verworfen, so muß gem¨ aß obigem Algorithmus aus {r

^N

}

i

ein weiterer Zustand {r

^N

}

j⁰

erzeugt werden

◦ Fortf¨ uhrung dieses Verfahrens f¨ uhrt zu einer Folge von Zust¨ anden, deren Wahrscheinlichkeit konstruktionsgem¨ aß (!) durch die Verteilungsfunktion des kanonischen Ensembles gegeben ist

◦ Berechnung von Mittelwerten gem¨ aß Formel (1) auf Folie 31

• Monte Carlo Simulation im kanonischen Ensemble (Spinsystem) anstelle der {r

^N

}

_j

treten die Spineinstellungen {s

^N

}

_i

; ein neuer Zustand wird dadurch erzeugt, indem man an einem zuf¨ allig

herausgegriffenen Teilchen die Richtung die Spineinstellung ver¨ andert

(37)

• Ablauf einer Monte Carlo Simulation vgl. Molekulardynamiksimulationen

◦ Initialisierungsphase

◦ Aquilibrierungsphase ¨

◦ Produktionsphase

(38)

(b) Monte Carlo Simulationen in anderen thermodynamischen Ensembles

• kanonisch-harmonisches (NPT) Ensemble

neben den Teilchenverschiebungen muß auch das Volumen V des Systems in einem Monte Carlo Schritt ver¨ andert werden

V wird wie eine zus¨ atzliche Koordinate behandelt: V → V

⁰

= V + δV , wobei δV innerhalb eines sinnvoll gew¨ ahlten Intervalls variiert werden kann/soll

Auswahlregeln:

(i) Teilchenverschiebung (Zust¨ ande i und j):

A(i → j ) = min {1, exp[−β(E

j

− E

i

)]}

(ii) Volums¨ anderung (V → V

⁰

):

A(V → V

⁰

) = min

1, exp

−β[E (V

⁰

) − E(V )] + P(V

⁰

− V ) − Nβ

⁻¹

ln(V

⁰

/V )]

(39)

• großkanonisches (µVT) Ensemble

neben den Teilchenverschiebungen k¨ onnen gem¨ aß dem Wert des chemischen Potentials, µ, Teilchen in das System eingesetzt (N → N + 1) oder Teilchen aus dem System entfernt (N → N − 1) werden

Auswahlregeln:

(i) Teilchenverschiebung (Zust¨ ande i und j):

A(i → j ) = min {1, exp[−β(E

j

− E

i

)]}

(ii) Teilchenzahl¨ anderung:

A(N → N + 1) = min

1, V

Λ

³

(N + 1) exp [β[µ − E (N + 1) + E (N)]]

A(N → N − 1) = min

1, Λ

³

N

V exp [−β [µ + E (N − 1) − E (N)]]

(40)

Gibbs Ensemble Monte Carlo Simulation (vgl. Ref. 4.6])

Schematische Darstellung der verschiedenen Ver¨anderungen am System inGibbs Ensemble Monte Carlo Simulationen(oben). Unten: Teilchendichte in denTeilsystemen (= Phasen)als Funktion der Simulationsschritte in einemphasentrennenden Lennard-Jones System; aus Ref.

[4.6].

(41)