7

Monte Carlo-Simulationen

Aufgabe 7.1: Das Lennard-Jones Potential, "  12

u(r) = 4ε

σ r

−

 6 #

σ r

,

hat sich als ein recht brauchbares einfaches Modell f¨ ur die effektive Paarwechselwirkung zwischen Edelgasatomen bzw. sph¨arisch symmetrischen Molek¨ ulen erwiesen. Als Paradebeispiel wird gern Argon angef¨ uhrt, wobei die konventionelle (aber nicht unbedingt optimale) Wahl der Parameter ε/kB = 119.8 K und σ = 3.405 ˚ A ist. F¨ uhre eine N V T -MC Simulation von fl¨ ussigem Argon am Tripelpunkt (T = 83.80 K, V = 28.24 cm3 /mol) mit einem System von N = 256 Teilchen durch und berechne die Autokorrelationsfunktionen f¨ ur die Momentanwerte der Kontrollgr¨oßen “potentielle Energie” und “Druck”. Was l¨aßt sich daraus u ¨ber die Korrelation der Stichprobe und den statistischen Fehler der Mittelwerte von U und p aussagen? Zum Vergleich: Die experimentell bestimmten Werte f¨ ur den Druck und die potentielle Energie der fl¨ ussigen Phase von realem Argon am Tripelpunkt sind p = 0.689 bar bzw. U = −5.97 kJ/mol. Der Konfigurationsanteil der spezifischen W¨arme (d.h. abz¨ uglich des Beitrags der kinetischen Energie) ist Cv 0 = 6.9 J/(K mol). Aufgabe 7.2: In der Praxis wird die Absch¨atzung der statistischen Fehler meist mit Hilfe sogenannter “Blockmittelwerte” durchgef¨ uhrt. Dazu werden die Momentanwerte einer Observablen A zu Bl¨ocken der L¨ange n zusammengefaßt, A1 , A2 , . . ., An , An+1 , An+2 , . . ., A2n , . . .,

|

{z

} |

1. Block

{z

}

2. Block

und die Mittelwerte (A¯(n) )i = [¨ uber den i-ten Block der L¨ange n] berechnet. Wenn die Anzahl der Einzelbeobachtungen insgesamt M = n·m ist, dann k¨onnen m solche Blockmittelwerte gebildet werden. Da die Einzelbeobachtungen jedes Blocks i.a. eine korrelierte Stichprobe darstellen, wird f¨ ur die Varianz der Blockmittelwerte gelten  

σA2¯(n) = σA2¯(1)

s , n

wobei σA2¯(1) ≡ σA2 die Varianz von A ist und s die “statistische Ineffizienz” (d.h. im wesentlichen die Korrelationsl¨ange). Man kann daher so vorgehen, daß man Blockungen verschiedener L¨ange n vornimmt und s aus dem Plateauwert n σA2¯(n) n→∞ σ 2 ¯(1) A

s = lim

2

bestimmt. F¨ ur die Varianz (den Standardfehler) des Gesamtmittelwertes A¯ ergibt sich dann σA2¯ =

σA2 . M/s

Wende dieses Verfahren auf die Daten der vorigen Aufgabe an. Aufgabe 7.3: Berechne f¨ ur das Lennard–Jones und das exp–6 Potential, u(r) = A e−br − C r−6 , die “Hintergrundkorrekturen” zur potentiellen Energie, Z

∆Uc /N = 2πρ

∞

dr r2 u(r),

rc

und zum Virial, Z

∆Wc /N = 2πρ

∞

dr r2 [r u0 (r)].

rc

Diese Beitr¨age sollen n¨aherungsweise die auf Grund des Cutoffs vernachl¨assigten Wechselwirkungen ber¨ ucksichtigen. Wie groß sind die Korrekturen zu potentieller Energie, Virial und Druck der fl¨ ussigen Phase des Lennard–Jones Systems am Tripelpunkt, wenn in der Simulation N = 256 Teilchen und ein Cutoff von rc = L/2 verwendet werden? (L ist die Kantenl¨ange des w¨ urfelf¨ormigen Simulationsvolumens.) Aufgabe 7.4: Bei der Besch¨aftigung mit einfachen Modellsystemen verwendet man gern reduzierte oder “*”-Einheiten, d.h. daß alle dimensionsbehafteten Gr¨oßen durch Multiplikation mit einer geeigneten Kombination der Parameter des Modells dimensionslos gemacht werden. Im Fall des Lennard–Jones Potentials sind das die Teilchenmasse m, der Teilchendurchmesser σ und die Potentialtiefe ε. Massen, L¨angen und Energien werden also in Einheiten von m, σ und ε gemessen. Was sind die entsprechenden Einheiten f¨ ur Zeit, Ladung, Viskosit¨at und Wirkung? Welchem Absolutwert des Drucks entspricht die Aussage p∗ = 1, wenn zur Interpretation die Parameter von Argon (s. Aufgabe 7.1) verwendet werden? Im Gegensatz dazu ist im Rahmen “realistischer” Simulationen vielfach ein Einheitensystem in Gebrauch, das auf der atomaren Masseneinheit, einem Nanometer als L¨angen- und einer Picosekunde als Zeiteinheit aufbaut. Daraus ergibt sich z.B. die f¨ ur physikalisch– chemische Fragestellungen u ¨bliche Energieeinheit von 1 kJ/mol. Wie groß sind in diesem Einheitensystem ein Druck von 1 atm, die elektrischen Elementarladung und h ¯?

3

Aufgabe 7.5: In statistisch-mechanischen Simulationen interessiert man sich meist f¨ ur den sogenannten “thermodynamischen Limes”, d.h. die Eigenschaften eines unendlich großen Systems. F¨ uhre eine Reihe von MC-Simulationen des Systems (fl¨ ussiges) LJ-Argon am Tripelpunkt mit N = 108, 256, 500, . . . Teilchen durch und untersuche die N -Abh¨angigkeit von potentieller Energie U , Druck p und Paarverteilungsfunktion g(r). Aufgabe 7.6: Die mittlere Anzahl von Nachbarteilchen in einer Kugel vom Radius R um ein Referenzteilchen, Z

n(R) = 4πρ

R

dr r2 g(r),

0

wird als “laufende Koordinationszahl” bezeichnet. Integriert man insbesondere bis R = Rmin , d.h. bis zum Ort des ersten Minimums von g(r), so wird das Resultat gern als “Koordinationszahl” schlechthin oder als Anzahl der n¨achsten Nachbarn interpretiert (obwohl diese Definition keineswegs eindeutig ist). Berechne f¨ ur einen oder mehrere Zust¨ande auf den Koexistenzlinien des LJ-Systems aus der folgenden (aus verschiedenen Quellen zusammengestellten) Tabelle n(R) und bestimme daraus die Koordinationszahl. T∗ ρ∗g 2.74 1.35 1.15 0.073 0.75 0.0035 0.694 0.0020

ρ∗` ρ∗` ρ∗s 1.167 1.230 0.982 1.056 0.605 0.947 1.028 0.822 0.856 0.967 0.84 0.84 0.96

(Die letzte Zeile enth¨alt die Dichten f¨ ur die am Tripelpunkt koexistierenden drei Phasen. Da es sich hier um die Daten f¨ ur das LJ-Modell handelt, m¨ ussen sie nicht genau mit den reduzierten Werten f¨ ur Argon von Aufgabe 7.1 u ¨bereinstimmen.) Aufgabe 7.7: Zwischen der Paarverteilungsfunktion g(r) und dem experimentell mit R¨ontgenoder Neutronenstreuung meßbaren Strukturfaktor S(k) besteht im Prinzip der eindeutige Zusammenhang 4πρ Z ∞ S(k) = 1 + dr r sin kr [g(r) − 1] , k 0 g(r) = 1 +

1 Z∞ dk k sin kr [S(k) − 1] . 2π 2 ρr 0

Berechne f¨ ur einen der in den vorhergehenden Aufgaben simulierten Zust¨ande des LJSystems den Strukturfaktor aus g(r) und versuche, aus diesem fiktiven “experimentellen” S(k) durch R¨ ucktransformation wieder die Paarverteilungsfunktion zu gewinnen. [In der Praxis werden bei der numerischen Transformation in beiden Richtungen die Integrale an Stellen abgebrochen, wo g(r) bzw. S(k) gleich 1 sind. Bei g(r) ergibt sich der gr¨oßtm¨ogliche Cutoff f¨ ur das Integral aus der Dimension der bei der Simulation verwendeten Box; bei der Inversion experimenteller Daten ist ein Cutoff von der Gr¨oßenordnung kσ ≈ 20 typisch.] 4

R¨ontgenstrukturfaktor Hd (k) = S(k) − 1 von CH4 bei T = 92 ± 1 K, ρ = −3 0.01702 ˚ A , aus A. Habenschuss, E. Johnson & A. H. Narten, J. Chem. Phys. 74, −1 5234 (1981); k in ˚ A . k 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0

Hd (k) -0.969 -0.969 -0.969 -0.969 -0.966 -0.960 -0.950 -0.935 -0.911 -0.880 -0.855 -0.834 -0.811 -0.776 -0.712 -0.579 -0.286 0.362 1.392 1.635 0.867 0.232 -0.114 -0.297 -0.387 -0.421 -0.419 -0.366 -0.284 -0.180 -0.047 0.114 0.207 0.238 0.226 0.180 0.143 0.090 0.028 -0.043 -0.091

k

Hd (k)

4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7.0 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 8.0

-0.112 -0.124 -0.110 -0.090 -0.065 -0.021 0.021 0.048 0.064 0.067 0.067 0.047 0.027 0.015 -0.001 -0.019 -0.026 -0.039 -0.034 -0.035 -0.019 -0.021 -0.003 0.006 0.018 0.011 0.022 0.011 0.016 0.002 0.002 -0.002 -0.002 -0.005 -0.009 -0.009 -0.007 -0.001 -0.0085 -0.003

k

Hd (k)

k

Hd (k)

8.1 0.001 12.1 0.002 8.2 0.002 12.2 -0.003 8.3 0.002 12.3 -0.003 8.4 0.009 12.4 -0.007 8.5 0.007 12.5 -0.006 8.6 0.008 12.6 -0.007 8.7 0.002 12.7 -0.001 8.8 -0.002 12.8 0.010 8.9 0.001 12.9 0.002 9.0 -0.001 13.0 0.006 9.1 -0.005 13.1 0.003 9.2 -0.009 13.2 0.004 9.3 0.002 13.3 0.003 9.4 -0.004 13.4 0.004 9.5 -0.004 13.5 -0.001 9.6 0.004 13.6 -0.005 9.7 -0.004 13.7 -0.001 9.8 0.001 13.8 -0.001 9.9 0.004 13.9 -0.001 10.0 0.003 14.0 -0.005 10.1 0.008 14.1 -0.003 10.2 -0.003 14.2 0.001 10.3 -0.001 14.3 -0.000 10.4 -0.004 14.4 0.000 10.5 -0.002 14.5 0.001 10.6 -0.000 14.6 -0.002 10.7 0.001 14.7 0.001 10.8 -0.004 14.8 0.005 10.9 -0.004 14.9 0.003 11.0 0.002 15.0 0.002 11.1 -0.004 15.1 0.002 11.2 0.004 15.2 0.001 11.3 0.004 15.3 -0.001 11.4 0.001 15.4 -0.003 11.5 0.001 15.5 -0.000 11.6 0.005 15.6 -0.004 11.7 0.003 15.7 -0.001 11.8 -0.000 15.8 -0.001 11.9 0.000 15.9 0.002 12.0 0.000 16.0 0.005

Aufgabe 7.8: Methan ist ein nahezu kugelf¨ormiges, bei nicht zu tiefen Temperaturen auch in den kondensierten Phasen relativ frei rotierendes Molek¨ ul. Im OPLS-”Baukasten” (Optimized Potentials for Liquid Simulations), einem f¨ ur die Simulation von Fl¨ ussigkeiten optimierten Satz von einfachen Modellpotentialen, finden sich daher f¨ ur Methan sowohl ein 5-Zentren- (Atom-Atom) als auch ein 1-Zentrenmodell (United Atom Model). Im ersten Fall ist jedes Atom Tr¨ager eines Lennard-Jones-Zentrums und einer Partialladung mit den Parametern εCC /kB = 33.2 K, σCC = 3.50 ˚ A, qC = −0.24 |e| εHH /kB = 15.1 K, σHH = 2.50 ˚ A, qH = +0.06 |e| (der intramolekulare CH-Abstand ist 1.09 ˚ A), im zweiten Fall wird das ganze Molek¨ ul durch ein einziges LJ-Zentrum mit ε/kB = 147.9 K,

σ = 3.73 ˚ A

ersetzt [G. Kaminski, E. M. Duffy, T. Matsui & W. L. Jorgensen, J. Phys. Chem. 98, 13077 (1994)]. F¨ uhre mit beiden Modellen eine Monte Carlo-Simulation am Siedepunkt von CH4 (T = 111.6 K, V = 37.97 cm3 /mol) durch und vergleiche die Atom–AtomVerteilungsfunktion gCC (r), die in diesem Fall mit der Paarverteilungsfunktion g(r) identisch ist, mit der inversen Fouriertransformierten des in der obigen Tabelle gegebenen Strukturfaktors aus einem R¨ontgenstreuexperiment. Zum Vergleich: Die potentielle Energie von Methan am Siedepunkt ist U = −7.37 kJ/mol, der Druck p = 1 atm. Aufgabe 7.9: Die folgende Abbildung zeigt die Atom–Atom-Verteilungsfunktionen, die mit Hilfe von Reverse Monte Carlo aus Neutronenstreuexperimenten an u ¨berkritischem deuteriertem Methan (CD4 ) bei hohem Druck gewonnen wurden [G. Strauß, A. Bassen, H. Zweier, H. Bertagnolli, K. T¨odheide, A. K. Soper & J. Turner, Phys. Rev. E 53, 3505 (1996)]. F¨ uhre eine Monte Carlo-Simulation mit dem 5-Zentrenmodell f¨ ur Methan bei einem der experimentell untersuchten Zust¨ande durch und vergleiche die Verteilungsfunktionen sowohl mit dem Experiment als auch mit der im vorigen Beispiel bestimmten Struktur von Methan am Siedepunkt.

6

7

Aufgabe 7.10: Die spezifische W¨arme (W¨armekapazit¨at) bei konstantem Druck, Cp = T

∂S ∂T

!

, p

und die isotherme Kompressibilit¨at, 1 χT = − V

∂V ∂p

!

, T

h¨angen beide von T und p als unabh¨angigen Variablen ab und lassen sich daher—sowohl theoretisch als auch in Simulationen—am einfachsten im isotherm-isobaren Ensemble berechnen. Ben¨ utze den Zusammenhang zwischen dem Gibbs-Potential, dG = −S dT + V dp + µ dN, und der zugeh¨origen Zustandssumme, e−βG = QT pN , sowie die Definition von QT pN , um Fluktuationsformeln herzuleiten, mit deren Hilfe Cp und χT in isotherm-isobaren Monte Carlo-Simulationen berechnet werden k¨onnen. Aufgabe 7.11: F¨ ur viele Substanzen ist der Tripelpunktsdruck so klein, daß man in guter N¨aherung pt = 0 setzen kann. Die reduzierte Tripelpunktstemperatur f¨ ur das Lennard– ∗ Jones System liegt bei T ≈ 0.694 (vergleiche Aufgabe 7.6). Bestimme die Dichten der fl¨ ussigen und festen Phase bei dieser Temperatur mit Hilfe von N pT -Monte Carlo-Simulationen. Als Startkonfiguration f¨ ur die Simulation der Fl¨ ussigkeit kann im Prinzip eine beliebige, also z.B. eine etwas gr¨oßere oder geringere als die erwartete Dichte gew¨ahlt werden; die korrekte Dichte stellt sich durch die Vorgabe des Drucks im Lauf der Rechnung dann von selbst ein. (Ist die Anfangsdichte allerdings zu klein, kann die Simulation auch in der Gasphase enden!) Bei der festen Phase wird man hingegen mit einem Kristallgitter in der N¨ahe der richtigen Dichte beginnen, um zu verhindern, daß das System versehentlich “schmilzt”. Aufgabe 7.12: Das Phasengleichgewicht von Reinsubstanzen ist dadurch gekennzeichnet, daß Temperatur T , Druck p und chemisches Potential µ der koexistierenden Phasen u ¨bereinstimmen. Da T , V und µ im großkanonischen Ensemble die unabh¨angigen Variablen sind, kann man Monte Carlo-Simulationen in diesem Ensemble dazu ben¨ utzen, um z.B. das Phasengleichgewicht gasf¨ormig–fl¨ ussig zu bestimmen. Dazu gibt man sich f¨ ur jede der beiden Phasen ein fixes Volumen vor und ermittelt jeweils in einer Reihe von Simulationen den Druck als Funktion des chemischen Potentials (die Temperatur ist ebenfalls fix 8

und f¨ ur beide Phasen dieselbe). Jenes chemische Potential, bei dem die beiden Druckkurven einander schneiden, entspricht dann dem Phasengleichgewicht, und aus den mittleren Teilchenzahlen in den beiden Volumina ergeben sich die Dichten der koexistierenden Phasen. Bestimme auf diese Weise das Phasengleichgewicht gasf¨ormig–fl¨ ussig des Lennard–JonesModells bei T ∗ = 1.15 oder einer anderen sinnvollen Temperatur. [Die Koexistenzdichten bei T ∗ = 1.5 sind ρg ≈ 0.07 und ρ` ≈ 0.61. Es ist vorteilhaft, diese Dichten im vorhinein wenigstens ungef¨ahr zu kennen, da man die Simulationsvolumina so groß w¨ahlen sollte, daß in beiden Phasen im Mittel etwa gleich viel (und mindestens 250–500) Teilchen sind.] Statt des chemischen Potentials µ verwendet man als Simulationsgr¨oßen gern andere Kombinationen der Parameter, in denen die (f¨ ur beide Phasen identische und daher f¨ ur das Gleichgewicht irrelevante) thermische de Broglie-Wellenl¨ange Λ subsumiert ist, z.B. in der Form z = eβµ /Λ3 . F¨ ur das Lennard–Jones-Modell bei T ∗ = 1.15 ist dann das gesuchte z ∗ ≈ 0.04. Aufgabe 7.13: Bestimme das Phasengleichgewicht gasf¨ormig-fl¨ ussig (d.h. den Gleichgewichtsdruck, das chemische Potential und die Dichten der Phasen) f¨ ur das Lennard-Jones∗ System bei der reduzierten Temperatur T = 1.15 mit Hilfe einer Gibbs-Ensemble Monte Carlo-Simulation mit insgesamt N = 250 und/oder 500 Teilchen. Da Gibbs-Ensemble-Simulationen wesentlich aufwendiger und empfindlicher sind als gew¨ohnliche MC-Simulationen, ist es vorteilhaft, schon eine ungef¨ahre Vorstellung von den Gleichgewichtsdichten zu haben (s. Tabelle zu Aufgabe 7.6). Eine Faustregel f¨ ur GEMCSimulationen besagt n¨amlich, daß die Teilchenzahl in der Gasphase wenigstens halb so groß sein sollte wie in der fl¨ ussigen Phase. Man wird daher das Gesamtvolumen und die Aufteilung von Volumen und Teilchenzahlen auf die beiden Boxen schon f¨ ur die Startkonfiguration so w¨ahlen, daß sie in der N¨ahe der erwarteten Endergebnisse liegen. Weiters sollten Z¨ uge vom Typ “Teilchentransfer” h¨aufig genug versucht werden, z.B. im Verh¨altnis Teilchenverschiebung:Teilchentransfer:Volumsaustausch wie (gr¨oßenordnungsm¨aßig) 100 : 10 : 1.

9