1. No category

1 Molekulardynamische Simulation eines Argon-Fluids
Aufgabenstellung
• Wir simulieren ein Argon-Fluid auf Basis eines Lennard-Jones Modells f¨ur Argon mit den
Modellparametern Ar = 0.34 nm und "Ar /kB = 120 K mit Hilfe der Methode der Molekulardynamischen Simulation.
• F¨ur zwei Temperaturen und bei kleiner Dichte (Gasphasen-Dichte) werden die Trajektorien
eines Systems bestehend aus 32 Teilchen unter Anwendung periodischer Randbedingungen
berechnet.
– Es werden ausgew¨ahlte Trajektorien stoßender Teilchen visualisiert und graphisch dargestellt.
– Die mittlere freie Wegl¨ange zwischen den St¨oßen wird anhand der Zahl der St¨oße und
des zur¨uckgelegten Weges berechnet und mit den theoretischen Werten verglichen.
– Anhand der mittleren quadratischen Verschiebung der Teilchen soll der Selbstdiffusionskoeffizient berechnet werden.
• Systeme von 256 Argon-Teilchen werden f¨ur vorgegebene Temperaturen und Dichten, variierend u¨ ber einen großen Dichtebereich, simuliert. Es sollen die Geschwindigkeitsverteilungen und radiale Verteilungsfunktionen f¨ur jeden dieser (T, ⇢) Zustandspunkte berechnet
und graphisch dargestellt werden.
• Anhand des berechneten Drucks und des berechneten chemischen Potentials entlang einer
Isotherme eines aus 256 Argon-Teilchen bestehenden Systems ist die Lage des Phasen¨ubergangs fl¨ussig/gas zu ermitteln.
Grundlagen
Allgemeines
Durch die hohe Rechengeschwindigkeit moderner Computer ist es m¨oglich geworden, das Verhalten physikalischer und chemischer Systeme durch “Computerexperimente” nachzuvollziehen.
Mittlerweile sind Computersimulationen zu einem neuen und wichtigen Werkzeug zur Untersuchung der Eigenschaften molekularer Systeme, insbesondere von Fl¨ussigkeiten, geworden. In den
Computerexperimenten wird das Verhalten einer begrenzten Zahl von Teilchen unter vorgegebenen Bedingungen simuliert. Dabei unterscheidet man zwei Methoden: Die statische M ONTE C ARLO-Methode (MC), bei der eine große Zahl von zuf¨allig erzeugten (daher der Name) Konfigurationen ausgewertet werden, und die Molekulardynamik-Simulation (molecular dynamics, MD),
bei der das dynamische Verhalten (die zeitliche Entwicklung) eines Ensembles von Teilchen entsprechend der Newton’schen Bewegungsgesetze berechnet wird. Beide Methoden behandeln das
Verhalten statistisch-thermodynamischer Ensembles und es k¨onnen daher die thermodynamischen
1
1 Computersimulation
und physikalischen Eigenschaften der simulierten Systeme als Ensemble-Mittelwerte bestimmt
werden.
F¨ur die Molekulardynamische Simulation ben¨otigt man ein Computerprogramm, das, ausgehend von gegebenen Anfangsbedingungen (Orte und Geschwindigkeiten der Teilchen), die Bewegung (Flugbahnen = Trajektorien) einer Anzahl von Teilchen unter Ber¨ucksichtigung der Kr¨afte
zwischen den Teilchen berechnet.
W¨ahrend der Simulation werden bestimmte Eigenschaften des Systems (Temperatur, Druck,
Geschwindigkeiten und Energien) ermittelt. Von Interesse ist auch die radiale Verteilungsfunktion, die Einblicke in die Struktur des Systems erlaubt. Der Diffusionskoeffizient der Molek¨ule
l¨asst sich aus der mittleren quadratischen Verschiebung der Teilchen bestimmen. Aus der Berechnung des chemischen Potentials l¨asst sich schließlich die Lage von Phasen¨uberg¨angen (hier:
fl¨ussig/gasf¨ormig) ermitteln.
Eine wesentliche Voraussetzung f¨ur die Simulationen ist die Kenntnis der Potentialfunktion, d.h.
einer Funktion, welche die Wechselwirkungsenergie eines Systems von Teilchen als Funktion ihrer
Ortskoordinaten beschreibt. F¨ur die Beschreibung der Molekularen Wechselwirkung der ArgonTeilchen verwenden wir ein sogenanntes “effektives” L ENNARD -J ONES-12-6-Potentialmodell,
dessen Parameter an den Verlauf der Fl¨ussig/Gas-Koexistenzlinie angepasst worden sind. Das Adjektiv “effektiv” deutet hierbei an, dass der Effekt von Mehrk¨orperwechselwirkungen durch die
gew¨ahlten Parameter des Zweik¨orper-Lennard-Jones Modells Ber¨ucksichtigung findet.
¨
Wechselwirkungspotential und Krafte
Zur Beschreibung der Potentialfunktion verwenden wir einen Paarpotentialansatz, d.h. die Gesamtwechselwirkungsenergie des Systems ergibt sich aus der Summe aller Paarwechselwirkungen. Jede Paarwechselwirkung wird hierbei durch eine Lennard-Jones-Funktion (LJ) definiert,
welche in guter N¨aherung die Wechselwirkung zwischen zwei nichtpolaren Teilchen als Funktion
ihres Abstands r beschreibt:
U (r) = 4 "
"✓ ◆
12
r
✓ ◆6 #
(1.1)
r
Der Potentialverlauf wird durch zwei charakteristische Gr¨oßen definiert: " gibt den Betrag der
Tiefe des Potentialminimums an der Stelle r" an und beschreibt den “Stoßdurchmesser”, d.h.
den Abstand, bei dem das Potential gerade Null ist. Es gilt:
✓
d U (r)
dr
◆
=0
r" = 21/6
U (r" ) =
"
U( ) = 0
(1.2)
r"
In Gleichung 1.1 repr¨asentiert der negative r 6 -Term weitreichende anziehende Wechselwirkungen (sogenannte Dispersionswechselwirkungen). Kurzfristige Fluktuationen der Elektronendichte
bewirken das Auftreten von zeitweiligen Dipolmomenten in einem Molek¨ul, welches dann seinerseits in einem zweiten Molek¨ul ein induziertes Dipolmoment erzeugt, woraus im Mittel eine elektrostatische Anziehung erfolgt. Bei gr¨oßerer Ann¨aherung der Teilchen kommt es zur Abstoßung,
die in diesem Fall durch den r 12 -Term ausgedr¨uckt wird. In Abbildung 1.1 ist das “effektive”
LJ-12-6-Potential f¨ur Argon dargestellt.
2
80
-12
~r
σ = 0.34 nm
-1
U(r) kB / K
40
1/6
σ 2
0
ε/kB = 120 K
σ
-40
-80
-6
~r
-120
0.3
0.4
0.5
0.6
r / nm
0.7
0.8
Abbildung 1.1: “Effektives” Lennard-Jones Wechselwirkungspotential U (r) f¨ur Argon. Die Parameter des verwendeten Lennard-Jones-Modells lauten: Ar = 0.34 nm und
"Ar /kB = 120 K [Rahman, A].
F¨ur die MD-Simulation werden neben den Wechselwirkungsenergien insbesondere auch die
Kr¨afte ben¨otigt. Die Kraft F~ij , die auf das Teilchen i durch das Teilchen j ausge¨ubt wird, ergibt sich aus der Ableitung des Potentials nach dem Abstand der Teilchen, multipliziert mit dem
normierten Abstandsvektor (um die Richtung der Kraftwirkung zu ber¨ucksichtigen):
F~ij =
ri U (rij ) =
d U (rij ) ~rij
·
d rij
rij
mit
rij = |~rij | = |~ri
~rj |
(1.3)
Die Gesamtkraft, die auf ein Teilchen i wirkt, ergibt sich dann als Summe u¨ ber alle Kr¨afte, die auf
das Teilchen i durch die umgebenden Teilchen j ausge¨ubt werden:
F~i =
X
F~ij
(1.4)
j
Berechnung der Teilchenbewegungen
Wir betrachten ein System aus N Teilchen. Zu einem Zeitpunkt t kennen wir f¨ur dieses System
s¨amtliche Orte ~ri (t) der Teilchen, ihre Geschwindigkeiten ~vi (t) und die auf die Teilchen wirkenden Kr¨afte F~i (t). Zur Berechnung der Orte und der Geschwindigkeiten zu einem kurz darauf
folgende Zeitpunkt t + t wird nun der “Velocity-Verlet”-Algorithmus verwendet: entsprechend
~ri (t +
t) = ~ri (t) + ~vi (t) ·
~vi (t +
t) = ~vi (t) +
t+
1
~ai (t) ·
2
1
[~ai (t) + ~ai (t +
2
t2
t)] ·
(1.5)
t
Die Beschleunigung des Teilchens i ergibt sich aus der auf das Teilchen wirkenden Kraft gem¨aß
~ai (t) = F~i (t)/mi . Das diskrete Zeitintervall t muss sehr klein sein, da sich die intermolekularen
¨
Kr¨afte sehr schnell a¨ ndern k¨onnen. Ublicherweise
w¨ahlt man Intervalle im Bereich von 10 15 bis
3
1 Computersimulation
Abbildung 1.2: Schematische Illustration zu den periodischen Randbedingungen. Die Basiszelle
(mitte, stark umrandet) ist umgeben von ihren virtuellen Abbildern. Aus Gr¨unden
¨
der Ubersichtlichkeit
sind hier nur die direkt benachbarten Zellen dargestellt.
10 14 s. F¨ur die Simulationen der Argon-Systeme verwenden wir hier t = 5 fs. Die Berechnung der gesamten Teilchenflugbahnen u¨ ber ein l¨angeres Zeitintervall erfolgt nun durch iterative
Wiederholung der Rechenschritte in Gleichung 1.5.
Das “Computerexperiment” l¨auft insgesamt nach folgendem Schema ab: Zun¨achst werden den
beteiligten Teilchen Startpositionen (z.B. auf einem Kristallgitter) und Anfangsgeschwindigkeiten zugewiesen. Letztere sind in der Regel Zufallsgr¨oßen, die allerdings so gew¨ahlt sind, dass
der Impuls des Schwerpunktes des Systems Null ist und die Geschwindigkeitsverteilung einer
Maxwell-Boltzmann-Verteilung bei einer vorzugebenden Temperatur entspricht. Die Flugbahnen
(Trajektorien) der Teilchen werden dann entsprechend Gleichung 1.5 schrittweise berechnet. Die
Geschwindigkeiten der Teilchen werden dabei allerdings zun¨achst so skaliert, das eine vorgegebene Temperatur eingestellt wird. F¨ur die hier untersuchten Systeme ist es ausreichend, insgesamt
etwa 104 · · · 106 Zeitschritte zu betrachten. Typischerweise teilt man die Simulation auf in einen
“Equilibrierungslauf”, bei welchem das System zun¨achst von einer beliebigen Startkonfiguration
in das thermodynamische Gleichgewicht gebracht wird. Daran schließt sich der eigentliche “Produktionslauf” an, f¨ur welchen die letzte Konfiguration (Orte und Geschwindigkeiten) des Equilibrierungslaufs als Startkonfiguration verwendet wird. S¨amtliche Simulationen und Auswertungen
werden mit Hilfe MD-Simulationspakets MOSCITO [MOSCITO] durchgef¨uhrt.
Die Basiszelle
Eine wesentliche Beschr¨ankung jeder Computersimulation ist die Anzahl der Teilchen, die man
einbeziehen kann. Bekanntlich enth¨alt 1 Mol 6.02214 ⇥ 1023 Partikel. Um f¨ur ein solches System auch nur einen einzigen Zeitschritt zu berechnen, w¨urde man selbst mit den heute schnellsten
Computern noch Jahrhunderte brauchen. Anstelle der Simulation einer makroskopischen Anzahl
von Teilchen reicht es jedoch in den meisten F¨allen aus, einen Ausschnitt aus der Fl¨ussigkeit
zu betrachten und diesen, analog zu einer Elementarzelle in einem Kristall, unendlich oft in alle
Raumrichtungen zu duplizieren. In Abh¨angigkeit vom betrachteten physikalischen Problem werden molek¨uldynamische Rechnungen (MD) mit etwa 102 bis 107 Teilchen pro Volumeneinheit
4
2000
-6
2
2
<|r(0)-r(t)| > / (nm)
2
D=1.6 ·10 m s
-1
1500
linear
1000
parabolisch
500
0
0
50
100
150
t / ps
200
250
Abbildung 1.3: Mittlere quadratische Verschiebung der Argon-Teilchen bei einer Temperatur von
T = 150 K und einer Dichte von 0.008 g cm 3 .
(Zelle oder Box) durchgef¨uhrt. In unserem Fall sind es 32 bzw. 256 Teilchen. Das Verfahren, die
Zelle oder Box mit den Teilchen nach allen Seiten unendlich oft fortzusetzen, beruht auf periodischen Randbedingungen. Das ist in Abbildung 1.2 f¨ur den zweidimensionalen Fall demonstriert.
Verl¨asst ein Teilchen die Box auf der einen Seite, so tritt auf der gegen¨uberliegenden Seite das
virtuelle Abbild des Teilchen aus der benachbarten Box mit gleicher Geschwindigkeit ein. F¨ur
die Berechnung der Wechselwirkung eines Teilchens mit seiner Umgebung wird jeweils immer
der Abstand zu dem am n¨achsten befindlichen Abbild der Teilchen betrachtet (siehe gestrichelte
Linien in Abbildung 1.2).
Diffusionskoeffizient
Die Bewegung von Teilchen in der Gasphase und in Fl¨ussigkeiten ist “diffusiv”. Die Ursache
hierf¨ur ist, dass die Trajektorien der Teilchen durch St¨oße mit umgebenden Teilchen abgelenkt
werden und die Bewegungsrichtung nach der Kollision quasi zuf¨allig verteilt ist. Die Trajektorie
eines Argon-Teilchens entspricht daher schon nach einem sehr kurzen Zeitintervall dem Verhalten eines sogenannten dreidimensionalen “random walkers”. F¨ur diffundierende Teilchen gilt die
“Einstein-Beziehung”, nach welcher sich die mittlere quadratische Verschiebung der Teilchen proportional zur Zeit verh¨alt:
D
|~ri (0)
E
~ri (t)|2 = 6D · t
(1.6)
D bezeichnet hierbei den Selbstdiffusionskoeffizienten der Molek¨ule. Im Zeitbereich unterhalb
des mittleren Kollisionsintervalls weist die mittlere quadratische Verschiebung einen parabelf¨ormigen Verlauf auf, da sich die Teilchen zwischen den St¨oßen mit konstanter Geschwindigkeit bewegen (sog. ballistischer Bereich). Im angegebenen Beispiel in Abbildung 1.3 geht der ballistische
Bereich etwa ab 120 ps in den linearen, diffusiven Bereich u¨ ber. Aus der Steigung der Geraden
(f¨ur t > 120 ps) l¨asst sich entsprechend Gleichung 1.6 der Selbstdiffusionskoeffizient des Argons
bei der gegebenen Temperatur und Dichte ermitteln.
5
1 Computersimulation
4
4
ρ=0.05 g cm
-3
ρ=1.0 g cm
3
T=90 K
g(r)
g(r)
3
-3
2
1
T=90 K
2
1
exp(-ULJ(r)/kBT)
0
0.4
0.6
0.8
r / nm
1
1.2
0
0.4
0.6
0.8
1
1.2
r / nm
Abbildung 1.4: Simulierte radiale Verteilungsfunktionen f¨ur Argon bei T = 90 K. Im Falle der
Gasphasendichte (links) wird die radiale Verteilungsfunktion im wesentlichen
durch den Verlauf der Zweiteilchen-(LJ)-Potentialfunktion bestimmt. Im Falle
der Fl¨ussigkeitsdichte (rechts) sind mehrere Maxima erkennbar, die auf eine “packungsinduzierte” Schalenstruktur in der Fl¨ussigkeit hinweisen.
Radiale Verteilungsfunktion
Auch Gase und Fl¨ussigkeiten besitzen eine Ordnung. Im Gegensatz zu Festk¨orpern existiert in diesen Phasen allerdings keine Fernordnung, sondern nur eine lokale Nahordnung. Gase und Fl¨ussigkeiten sind zudem isotrop, d.h., die Nahordnung ist unabh¨angig davon, welche Raumrichtung
betrachtet wird. Die Nahordnung in Gasen und Fl¨ussigkeiten besteht daher in der Korrelation von
Teilchenabst¨anden. Die radiale Verteilungsfunktion oder Paarverteilungsfunktion g(r) beschreibt
diese Teilchenkorrelationen. Sie ist ein Maß f¨ur die Dichte, die wir um ein zentrales Teilchen
herum im Abstand r finden, bezogen auf die mittlere Dichte des Fluids. Dazu bestimmt man die
Zahl der Teilchen in einer Kugelschale der Dicke r mit dem Volumen V = 4⇡r2 r im Abstand zwischen r und r + r relativ zu einem Referenzteilchen und normiert durch die mittlere
Teilchenzahldichte ⇢ = N/V :
g(r) =
⇢(r)
1
1
= ·
· h n(r)i
⇢
⇢ 4⇡r2 r
(1.7)
Der Wert der Paarverteilungsfunktion strebt f¨ur kleine Abst¨ande gegen Null, da sich am Ort r = 0
das Referenzteilchen befindet: lim r!0 g(r) = 0. F¨ur große Abst¨ande konvergiert die Dichte ⇢(r)
gegen die mittlere Teilchenzahldichte und die Paarverteilungsfunktion strebt entsprechend gegen
Eins: lim r!1 g(r) = 1.
F¨ur Fl¨ussigkeiten weist g(r) einen relativ steilen “Peak” etwa bei r" auf, der sich, schw¨acher
werdend, periodisch wiederholt. Je gr¨oßer r wird, desto breiter werden die Peaks. Letztlich n¨ahert
sich g(r) dem Wert 1. Die Ursache f¨ur die alternierenden Maxima und Minima ist die schalenf¨ormige Anordnung von Teilchen um ein zentrales Teilchen herum, wenn die Dichte hinreichend groß ist. Aber auch in der Gasphase finden wir eine Paarverteilungsfunktion mit einem
ausgepr¨agten ersten Peak (siehe Abbildung 1.4). Dieser r¨uhrt daher, dass in der Gasphase die
Verteilung der Teilchenabst¨ande durch das Boltzmann-gewichtete Paarpotential bestimmt wird.
6
-11.6
-11.8
-1
µeq = -12.19±0.05 kJ mol
-12
-12.2
peq= 19.8±1.2 bar
Gasphas
e
-12.4
-12.6
-12.8
-20
0
µ / kJ mol
-1
µ / kJ mol
ρg = 0.100±0.009 g cm
sig
ues
-12
-1
Fl
keit
-1
µeq=-12.19 kJ mol
-12.4
-12.8
20 40
p / bar
60
-3
ρl = 1.108±0.013 g cm
0
80
0.2 0.4 0.6 0.8
-3
ρ / g cm
1
-3
1.2
Abbildung 1.5: Links: Berechnetes chemisches Potential µ als Funktion des Drucks entlang einer Isotherme. Aus dem Schnittpunkt von Gasphasen-Ast und Fl¨ussigkeits-Ast
k¨onnen der Druck und das chemische Potential am Gleichgewichtspunkt abgelesen werden. Rechts: Chemisches Potential als Funktion der Dichte. Aus den
Schnittpunkten mit dem chemischen Potential des Gleichgewichtszustandes l¨asst
sich die Dichte der Gasphase und der Fl¨ussigkeit ermitteln.
Abst¨ande in der N¨ahe von r" werden demnach bevorzugt und es gilt g(r) ⇡ exp[ U (r)/kB T ] f¨ur
⇢ ! 0.
Chemisches Potential und Flussig/Gas-Phasen
¨
ubergang
¨
Die Bestimmung von Phasen¨uberg¨angen mit der Computersimulation ist mit Hilfe des chemischen
Potentials m¨oglich. Am Siedepunkt beispielsweise sind der Druck und das chemische Potential in
der fl¨ussigen Phase und der Gasphase gleich. Dies wollen wir ausnutzen, um den Phasen¨ubergang
zwischen der fl¨ussigen Phase und der Gasphase entlang einer Isotherme zu lokalisieren.
Das chemische Potential µ besteht formal aus einem Anteil µid f¨ur das ideale Gas, der analytisch
berechnet werden kann, und einem wechselwirkungsbasierten Anteil µres , der numerisch ermittelt
werden muss:
µ = µid + µres
(1.8)
F¨ur den Anteil des idealen Gases gilt hierbei:
µid =
p
kB T ln
V
⇤3 · N
.
(1.9)
Dabei bezeichnet ⇤ = h2 /(2⇡mkB T ) die thermische de Broglie Wellenl¨ange. Man kann formal
beweisen, dass sich der wechselwirkungsbasierte Anteil des chemischen Potentials µres als ein
Ensemble-Mittelwert entsprechend
µres =
kB T ln hexp[
/(kB T )]i
(1.10)
berechnen l¨asst. ist hierbei die Energie eines in die N -Teilchen-Simulation zus¨atzlich eingef¨ugten (N + 1)-ten Teilchens gem¨aß = UN +1 UN . Die Positionen des zus¨atzlich eingef¨ugten Teilchens werden nach dem Zufallsprinzip ermittelt und sind u¨ ber das Volumen der Box
7
1 Computersimulation
v¨ollig gleichverteilt. Das (N + 1)-te Teilchen beeinflusst dabei die Trajektorie des N -TeilchenSystems nicht, weshalb man auch von einem “Testteilchen” spricht.
Durch eine Auftragung des chemischen Potentials µ gegen den Druck p l¨asst sich der Phasen¨ubergang durch den Schnittpunkt von µg (pg ) des Gasphasen-Astes und µl (pl ) des Fl¨ussigkeitsAstes ermitteln (siehe Abbildung 1.5). Am Schnittpunkt gilt jeweils pl = pg und µl = µg .
Versuchsapparatur
Die Rechnungen erfolgen auf einem Linux/Unix-PC, der Ihnen vom Assistenten zugewiesen wird.
Es findet eine Einweisung in die ben¨otigten Unix-Kommandos durch den Assistenten statt.
Versuchsdurchfuhrung
¨
Zu Beginn des Versuches bekommen Sie die genaue Aufgabenstellung f¨ur die Rechnungen vom
Assistenten mitgeteilt.
Hinweis: Das Simulationsprogramm erzeugt zahlreiche unterschiedliche Ausgabedateien. Um
Probleme durch gleich lautende Dateinamen zu vermeiden, legen Sie f¨ur jede Simulation am besten ein separates Unterverzeichnis an und starten die Simulation dort.
Simulation mit 32 Teilchen
8
Kommandos:
Aufbau der Simulation (Startkonfiguration):
MD Simulation bei Temperatur:
Visualisierung der Systemtrajektorie mittels VMD
Visualisierung einzelner Trajektorien
Zur¨uckgelegte Wegstrecke pro Atom (nm)
buildlj32.pl -rho ⇢
mdsimlj32.pl -T T
argonview.pl
xmgrace sample 01.xy sample 01.xz
tralen.pl
Dateien:
Trajektorien (xy-Projektion, in nm)
Trajektorien (xz-Projektion, in nm)
Systemtrajektorie (bin¨ar)
Geschwindigkeiten (bin¨ar)
Simulationsprotokoll
Mittlere quadratische Verschiebung (ps, nm2 )
Mittlere quadratische Verschiebung (ps, nm2 )
sample 00.xy ...
sample 00.xz ...
sim.xtc
sim.vel
sim.log
msd.dat
msd.csv
sample 09.xy
sample 09.xz
Simulation mit 256 Teilchen
Kommandos:
Aufbau der Simulation (Startkonfiguration):
MD Simulation bei Temperatur:
Visualisierung der Systemtrajektorie mittels VMD
buildlj256.pl -rho ⇢
mdsimlj256.pl -T T
argonview.pl
Dateien:
Paarverteilungsfunktion
Geschwindigkeitsverteilung
Systemtrajektorie (bin¨ar)
Geschwindigkeiten (bin¨ar)
Simulationsprotokoll
gofr Ar Ar.dat, gofr Ar Ar.csv
v histogram.dat, v histogram.csv
sim.xtc
sim.vel
sim.log
Isotherme
Kommandos:
Berechnung der Isotherme (Dauer: ca. 1 h)
Aufbereitung der Daten
run isotherm.pl -T T
format data
Dateien:
Thermodynamische Daten (Klartext)
Thermodynamische Daten (Excel)
Thermodynamische Daten Protokoll
Isotherm.txt
Isotherm.csv
Isotherm.pdf
Versuchsauswertung
1. T RAJEKTORIEN
Analysieren Sie mit Hilfe des Visualisierungsprogrammes VMD (argonview.pl) die
Trajektorien einzelner Argon-Teilchen. Erzeugen Sie zwei Abbildungen, die jeweils charakteristisch einen direkten Stoß als auch einen streifenden Stoß abbilden. Als Letzteren
verstehen wir hierbei eine Richtungs¨anderung, die wesentlich durch die attraktive Wechselwirkung zweier Teilchen bestimmt ist.
¨
2. M ITTLERE FREIE W EGL ANGE
Die Trajektorien der Teilchen im verd¨unnten Gas werden genutzt, um die mittlere freie
Wegl¨ange l zu berechnen. Bestimmen Sie dazu zun¨achst mit dem Programm tralen.pl
die L¨ange des zur¨uckgelegten Weges der Argon-Partikel. Anschließend plotten Sie die Projektionen der Trajektorien einzelner Teilchen in die xz- und xy-Ebene mit Hilfe des Programmes xmgrace:
z.B. xmgrace sample 01.xy sample 01.xz
Drucken Sie die Trajektorien aus und ermitteln Sie anhand der markierten Richtungs¨anderungen die Zahl der erfolgten St¨oße. Berechnen Sie so die mittlere freie Wegl¨ange f¨ur f¨unf
ausgew¨ahlte Argon-Teilchen und bestimmen Sie den Mittelwert. Vergleichen Sie den so
erhaltenen Wert mit dem Resultat nach der kinetischen Gastheorie: l = 1/(21/2 ⇡⇢ 2 ).
Erl¨autern Sie, warum die erhaltene mittlere freie Wegl¨ange typischerweise kleiner ist als
von der kinetischen Gastheorie vorhergesagt, wenn als Stoßdurchmesser der Lennard-JonesParameter verwendet wird.
9
1 Computersimulation
3. D IFFUSIONSKOEFFIZIENTEN
Die u¨ ber alle Teilchen gemittelte mittlere quadratische Verschiebung ist in den Dateien
“msd.dat” bzw. “msd.csv”
enthalten. Die erste Spalte beschreibt die Zeit in ps, die zweite Spalte bezeichnet die mittlere
quadratische Verschiebung in (nm)2 . Aus der Steigung des linearen Bereichs kann der Diffusionskoeffizient ermittelt werden. Aus dem Diffusionskoeffizienten kann dann mittels des
Zusammenhangs D = 1/3 · hvi · l gem¨aß der kinetischen Gastheorie ebenfalls die mittlere
freie Wegl¨ange l berechnet werden.
4. G ESCHWINDIGKEITSVERTEILUNG
Die berechneten Geschwindigkeitsverteilungen der Teilchen sind in den Dateien
“v histogram.dat” bzw. “v histogram.csv”
abgelegt. Die erste Spalte beschreibt die Geschwindigkeit (in m s 1 ), die zweite Spalte die
entsprechende Wahrscheinlichkeitsdichte (in (m s 1 ) 1 ). Vergleichen Sie die so erhalte Geschwindigkeitsverteilung mit der Maxwell-Boltzmann-Geschwindigkeitsverteilung gem¨aß
✓
M
P (v) = 4⇡
2⇡RT
◆3/2
· v 2 · exp[ M v 2 /(2RT )]
(1.11)
5. R ADIALE V ERTEILUNGSFUNKTION
Die Daten f¨ur die radiale Verteilungsfunktion g(r) sind in den Dateien
“gofr Ar Ar.dat” bzw. “gofr Ar Ar.csv”
abgespeichert. Die erste Spalte beschreibt den Abstand r (in nm) und die zweite Spalte
enth¨alt g(r) (dimensionslos). Zeichnen Sie die radiale Verteilungsfunktion zun¨achst bei
gleichen Temperaturen T f¨ur unterschiedliche Dichten ⇢ in ein Diagramm. Entsprechend
werden dann Diagramme bei konstantem ⇢ f¨ur verschiedene T erstellt. Wie liegt das erste
Maximum verglichen zu r" = 21/6 ?
¨
¨
6. F L USSIG
/G AS -P HASEN UBERGANG
In den Ausgabedateien
“Isotherm.txt” bzw. “Isotherm.csv”
werden die errechneten Daten f¨ur T [K], ⇢ [g cm 3 ], Vm [m3 mol 1 ], p [bar], µid. [kJ mol 1 ],
µres. [kJ mol 1 ] und µ [kJ mol 1 ] abgespeichert. Ein vollst¨andiges Protokoll der IsothermenBerechnung finden Sie weiterhin in der Datei
“Isotherm.pdf”
Die csv-Dateien k¨onnen zur Auswertung problemlos in ein Tabellenkalkulationsprogramm
importiert werden. Zun¨achst ist nun f¨ur die Isotherme µ als Funktion von p in einem Dia¨
gramm darzustellen. Hier erhalten Sie f¨ur die unterkritischen Isothermen zwei Aste,
zwischen denen das f¨ur die Betrachtung unbrauchbare metastabile Zwei-Phasen-Gebiet liegt.
¨
Um den eigentlichen Phasen¨ubergang zu erfassen, ermitteln Sie den Schnittpunkt der Aste
und bestimmen den Gleichgewichtsdruck peq (Dampfdruck) und das chemische Potential
µeq im Gleichgewicht.
10
Erstellen Sie nun ein Diagramm von µ als Funktion von ⇢ und tragen Sie µeq ein. An den
Schnittpunkten von µeq mit der Funktion µ(⇢) lassen sich die Gasphasendichte und die
Fl¨ussigkeitsdichte am fl¨ussig/gas-Koexistenzpunkt ablesen.
Literatur
Rahman, A. A. Rahman, “Correlations in the Motions of Liquid Argon”, Physical Review 136,
A405-A411 (1964).
MOSCITO MD-Simulationsprogramm: http://www.moscitomd.org
11