Implementierung und Vergleich weitererZeit-Integrationsverfahren

IDP „Numerische Aspekte beiMolekulardynamik-Simulationen“

Jakob Kummerow

15. September 2010

Technische Universität MünchenFakultät für InformatikLehrstuhl V: Wissenschaftliches Rechnen

Inhaltsverzeichnis1 Einleitung 2

1.1 Überblick über die Molekulardynamiksimulation . . . . . . . . . . . . . . . . . . . 21.2 Numerische Aspekte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 Theoretische Überlegungen und Herleitung der Integrationsstrategien 32.1 Integratorenvergleich . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 Euler-Methode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.3 Störmer-Verlet-Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.4 Leapfrog-Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.5 Predictor-Corrector-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.6 Runge-Kutta-Methoden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3 Implementierung 8

4 Auswertung 94.1 Betrachtung der einzelnen Integratoren . . . . . . . . . . . . . . . . . . . . . . . . . 94.2 Vergleich aller Integratoren für ausgewählte Zeitschrittweiten . . . . . . . . . . . . 134.3 Extra-lange Zeitschritte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

5 Fazit 16

1

1 Einleitung

1.1 Überblick über die MolekulardynamiksimulationFür viele Simulationsanwendungen genügt es, von kontinuierlichen Modellen der simulierten Ma-terialien auszugehen, und diese relativ grob zu diskretisieren. Bestimmte Fragestellungen, wiezum Beispiel das genaue Verhalten an Grenzflächen verschiedener Stoffe sowie Vermischungs-,Kristallisations- und Kondensationsprozesse, lassen sich so aber nicht untersuchen: In diesen Fäl-len ist es zielführender, die untersuchten Materialien Molekül für Molekül zu simulieren, was fürgrößere Simulationsgebiete zwar prinzipiell ebenso funktionieren würde, aber wegen der enormenAnzahl an Molekülen viel zu viel Rechenzeit benötigen würde. Interessiert man sich auf der anderenSeite für extrem kleine Maßstäbe, z.B. subatomare Prozesse, Quanteneffekte und chemische Reak-tionen, ist die Betrachtung einzelner Moleküle zu grob, sodass wiederum andere Herangehensweisenbenötigt werden. Der Anwendungsbereich der Molekulardynamik ist dadurch auf Interessensgebieteeiner bestimmten Größenordnung festgelegt.

Um physikalische Vorgänge im Bereich einzelner Moleküle korrekt zu simulieren, müsste maneigentlich Quanteneffekte berücksichtigen. Dies würde jedoch zu so komplizierten Modellen führen,dass nur sehr kleine Modelle numerisch lösbar wären. Eine weit verbreitete Vereinfachung bestehtdaher darin, sich auf Newton’sche Mechanik zu beschränken, d.h. die Moleküle als Kugeln im Raumzu modellieren, die jeweils eine gewisse Masse und Geschwindigkeit haben und in Abhängigkeit ihrerEntfernung miteinander interagieren. Die simulierte Zeit wird diskretisiert, d.h. die Simulationarbeitet mit Zeitschritten einer bestimmten Länge (die nicht notwendigerweise global konstantist). Im Kern besteht ein MD-Algorithmus also aus folgenden Schritten, die in einer Schleife fürjeden Zeitschritt wiederholt werden:

1. Zeitschritt initialisieren

2. Aus den aktuellen Positionen der Moleküle deren Interaktionskräfte berechnen

3. Unter Verwendung der wirkenden Kräfte die Position der Moleküle aktualisieren

Schritt 1 beinhaltet in Abhängigkeit des verwendeten Verfahrens diverse interne Operationen, hataber nichts mit der Physiksimulation als solcher zu tun.

Für Schritt 2 definiert man in der Regel eine Potenzialfunktion für das in einer bestimmtenMolekülkonfiguration enthaltene Potenzial. Häufig wird das so genannte Lennard-Jones-PotenzialULJ(r) = 4ε

((σr

)12 −(σr

)6) verwendet, das sowohl Abstoßung extrem nah beieinander befindli-cher Moleküle (erster Summand) als auch Anziehung aufgrund von van-der-Waals-Kräften (subtra-hierter Summand) modelliert, wobei σ die Größe der Moleküle beschreibt, r den aktuellen Abstandder beiden betrachteten Moleküle, und ε eine stoffspezifische Konstante ist. Die resultierende Kraftzwischen den beiden Molekülen ergibt sich als F (r) = − ∂

∂rU(r). Um den Rechenaufwand niedrigzu halten, werden typischerweise folgende Maßnahmen ergriffen:

• Es werden nur Interaktionen zwischen jeweils zwei Molekülen betrachtet. In der Realitätdarüber hinaus vorhandene Wechselwirkungen zwischen drei und mehr Molekülen werden

vernachlässigt. Auf diese Weise müssen bei n Molekülen im System nur(n2

)∈ O

(n2)

Molekülpaare betrachtet werden, statt∑nk=2

(nk

)≈ 2n Gruppen zu je k Molekülen.

• Das Lennard-Jones-Potenzial ist ein kurzreichweitiges Potenzial, da es mit zunehmender Ent-fernung „schnell“ gegen 0 geht (d.h. im d-dimensionalen Raum gilt ULJ ∈ o

(r−d)). Folglich

kann ohne großen Präzisionsverlust das Potenzial zwischen weit voneinander entfernt liegen-den Molekülen (r > rcutoff für einen definierten Abschneideradius rcutoff ) vernachlässigtwerden. Anhand dieser Erkenntnis wird die Anzahl zu betrachtender Moleküle weiter redu-ziert, indem die Moleküle in geeigneten, Nachbarschaft berücksichtigenden Datenstrukturenverwaltet werden, so dass im Erwartungswert nur eine konstante Anzahl Nachbarmolekülefür jedes Molekül angesehen werden müssen.

2

• Unter Ausnutzung des dritten Newton’schen Gesetzes („actio = reactio“) muss für jedesMolekülpaar die herrschende Kraft nur einmal berechnet werden, und kann (mit invertiertemVorzeichen) für beide beteiligten Moleküle verwendet werden.

Schritt 3 wird vom so genannten Integrator erledigt. Die Grundideen sind Newtons GleichungenF = m ·a, ∂

∂tx = v, ∂∂tv = a, oder zusammengefasst als eine einzige Differenzialgleichung: F = m · x

(wobei x die Position, v die Geschwindigkeit, a die Beschleunigung undm die Masse eines Teilchenssind). Da diese Gleichung mit vertretbarem Aufwand nur für extrem einfache Systeme analytischlösbar ist, besteht die Aufgabe des Integrators darin, eine möglichst gut annähernde Diskretisierungzu finden.

Um beliebige Moleküle korrekt simulieren zu können, muss man neben der Translation des Mo-leküls auch dessen Rotation berücksichtigen. Für diese Arbeit wird dieser Aspekt jedoch ignoriertund alle Moleküle werden als kugelförmig angenommen, da die Formeln und Lösungsmethoden fürDrehmoment/-geschwindigkeit nicht grundlegend anders sind als für die einfache Bewegung.

1.2 Numerische AspekteEine Diskretisierung hat im Wesentlichen mit zwei Schwierigkeiten, d.h. Ursachen für Fehlerhaf-tigkeit des berechneten Ergebnisses, zu kämpfen:

Zum einen entsteht bei der Diskretisierung notwendigerweise ein Diskretisierungsfehler, da dertatsächliche Verlauf der jeweiligen Funktion nur angenähert wird. Typischerweise ist der Diskretisie-rungsfehler mit der Zeitschrittweite korreliert; daher ist es Aufgabe des Integrators, eine möglichstgute Annäherung zu liefern, die möglichst große Zeitschritte zulässt, da der Gesamtrechenaufwandin etwa proportional zur Anzahl simulierter Zeitschritte ist.

Kleine Zeitschrittweiten führen zu einem anderen Problem: Beim Rechnen mit endlicher Präzi-sion treten beim Addieren stark unterschiedlich großer Zahlen Genauigkeitsverluste auf – addiertman beispielsweise eine sehr geringe zurückgelegte Wegstrecke auf die Koordinaten eines Moleküls,geht die Änderung im schlimmsten Fall beim notwendigen Runden auf die von den Koordinatendarstellbare Genauigkeit wieder komplett verloren. Einen ähnlichen Genauigkeitsverlust erleidetman bei der Subtraktion ähnlich großer Zahlen (sofern solche Operationen verwendet werden).Infolgedessen kann man (selbst unter Vernachlässigung des höheren Rechenaufwandes) nicht einenschlechten Integrator durch Wahl einer geringen Zeitschrittweite zu einem beliebig guten Integratormachen. Statt dessen ist dies ein zweiter Grund, weshalb es sich lohnt, einen Integrator zu verwen-den, der möglichst große Zeitschritte und damit möglichst signifikante Änderungen der internenModellvariablen zulässt.

Inhalt dieser Arbeit ist es, verschiedene Integratoren in den existierenden Molekulardynamik-Simulator MarDyn zu integrieren und miteinander zu vergleichen.

2 Theoretische Überlegungen und Herleitung der Integrati-onsstrategien

2.1 IntegratorenvergleichWill man verschiedene Integratoren vergleichen, stellt sich die offensichtliche Frage, wie man diesenVergleich anstellen kann.

Gewisse Aussagen kann man mittels rein theoretischer Betrachtung treffen. Theoretische Ana-lyse allein ist aber nicht genug, da beispielsweise zwei Näherungsverfahren mit Fehlertermen derGrößenordnung ∆t2 deswegen noch lange nicht in der Praxis gleich gut sein müssen1. Um das tat-sächliche Verhalten der Integratoren beurteilen zu können, muss man sie empirisch untersuchen,d.h. die von ihnen berechneten Ergebnisse vergleichen.

Aber wie entscheidet man, welches das bessere von zwei verschiedenen Ergebnissen ist? Im Fallvon einfachen Problemen (z.B. numerische Integration einer analytisch integrierbaren Funktion)

1Konvention: Der Übersichtlichkeit halber habe der Operator „∆“ stärkere Assoziativität als Potenzierungen wie„2“, d.h. ∆t2 = (∆t)2

3

kann man recht einfach das numerische Ergebnis mit dem theoretischen vergleichen; Molekulardy-namik mit mehreren tausend Molekülen ist für diesen Ansatz jedoch nicht geeignet, da die Trajek-torien der Moleküle analytisch nicht berechenbar sind. Die Integratoren einfach anhand einfacherProblemstellungen zu testen und die Gültigkeit der Ergebnisse auch für komplexe Simulationenanzunehmen, ist ebenfalls kein gangbarer Weg, da gerade in der Komplexität und dem Chaos derzu simulierenden Systeme eine besondere Herausforderungen für die Integratoren liegt, und diesedort möglicherweise deutlich anders abschneiden als bei einfach gelagerten Problemen.

Ein Ausweg ist, die Genauigkeit der Integratoren nur indirekt zu messen, nämlich über dieEnergieerhaltung. Bei exakter Rechnung bleibt die Gesamtenergie des Systems, d.h. die Summeaus der Bewegungsenergie der Moleküle und der durch das aufsummierte Potenzial repräsentiertenin der Konstellation enthaltenen Energie, konstant. Ändert sich dagegen die Gesamtenergie, hatder Integrator offenbar nicht exakt gerechnet; der Betrag der Abweichung gibt einen Hinweis aufdie Genauigkeit der erfolgten Rechnung.

Typischerweise simuliert man so genannte NVT-Ensembles, d.h. die Anzahl Moleküle, das Vo-lumen und die Temperatur werden konstant gehalten. Die Temperatur ist, vereinfacht gesagt, diedurchschnittliche Geschwindigkeit der Moleküle (Brown’sche Bewegung; für Details und den genau-en Zusammenhang sei auf die einschlägige Literatur verwiesen2), steht also in direktem Zusammen-hang mit der Gesamtenergie. Um beurteilen zu können, ob die Energie erhalten wird, wurde daherfür die Tests das Thermostat, das normalerweise mittels Skalierung aller Molekülgeschwindigkeitendie Temperatur konstant hält, im Programmcode deaktiviert.

Molekulardynamik-Simulationen sind chaotische Systeme, d.h. die Trajektorien der Molekülein zwei Systemen mit nur minimal unterschiedlichen Ausgangsbedingungen können schon nach re-lativ kurzer Zeit beliebig verschieden werden, bzw. anders betrachtet: schon kleinste Rechen- oderRundungsfehler können zu dramatisch unterschiedlichen Ergebnissen führen. Man könnte daherargumentieren, dass es gar nicht so wichtig ist, wie genau der Integrator nun rechnet, da es ohne-hin ein aussichtsloses Unterfangen ist, die Bewegung jedes einzelnen Moleküls exakt zu simulieren,sondern lediglich beispielhafte Verläufe berechnet werden können (die in vielen Fällen auch völ-lig ausreichende Informationen liefern). Ungeachtet dessen ist eine möglichst genaue Berechnungtrotzdem wünschenswert, vor allem bei ähnlichem Rechenaufwand. Außerdem können Ungenau-igkeiten nicht nur zu „gutartigen“ Verfälschungen des Ergebnisses führen, sondern beispielsweiseauch zu explosionsartigen Zuständen, wenn durch zu grobe Annäherung der Molekülpositionen dasPotenzial der Molekülkonstellation von einem Zeitschritt zum nächsten ins Unermessliche steigt.

2.2 Euler-MethodeDer wohl simpelste Algorithmus zur numerischen Integration der Molekülbewegung ist die so ge-nannte explizite Euler-Methode, die die wirkenden Kräfte sowie die daraus resultierende Geschwin-digkeit jeweils am Anfang des Zeitschritts betrachtet und stückweise konstant interpoliert:

x(t+ ∆t) = x(t) + v(t) ·∆tv(t+ ∆t) = v(t) + a(t) ·∆t

Diese Formeln lassen sich direkt aus der Taylor-Entwicklung von x(t) bzw. v(t) herleiten, indemdiese nach dem zweiten Term abgebrochen wird. Alternativ kann man sie auch als Umformung derFinite-Differenzen-Methode auffassen:

v(t) = x(t) ≈ x(t+ ∆t)− x(t)

∆t

Die Euler-Methode ist relativ ungenau: Da alle Taylor-Terme ab einschließlich dem ∆t2-Termabgeschnitten werden und in den zu erwartenden Fehler einfließen, handelt es sich um ein Verfahrenerster Ordnung. Anders gesagt, aufgrund der stückweise konstanten Interpolation sind sehr kleineZeitschritte erforderlich, um akzeptable Genauigkeit zu erzielen; aber selbst dann ist zu erwarten,dass die Euler-Methode den anderen hier vorgestellten Algorithmen unterlegen sein wird. DieVorteile der Euler-Methode seien trotzdem nicht unerwähnt: Der Algorithmus ist sehr einfach undbenötigt bei der Ausführung sehr wenig Speicher, nämlich nur Koordinaten und Geschwindigkeit(und temporär Kräfte bzw. Beschleunigung) jedes Moleküls.

2Eine Erklärung ist beispielsweise in [Griebel04], S. 91ff. zu finden.

4

Implizite Euler-Methode

In der bisher dargestellten Form wird die Euler-Methode „explizit“ genannt, weil sie die Kräfteam Anfang des jeweiligen Zeitschritts als Grundlage der Berechnungen verwendet. Zwar nichtdie theoretische Genauigkeit, wohl aber die akzeptablen Zeitschrittweiten lassen sich signifikantverbessern, indem man die Kräfte am Ende des Zeitschritts zugrunde legt:

v(t+ ∆t) = v(t) + a(t+ ∆t) ·∆tx(t+ ∆t) = x(t) + v(t+ ∆t) ·∆t

Das Problem an diesen Formeln ist allerdings, dass a(t+ ∆t) nur berechnet werden kann, wennx(t+ ∆t) bereits bekannt ist; in anderen Worten: diese Formeln stellen ein Gleichungssystem dar,in dem nicht (wie bei der expliziten Euler-Methode) jede Zeile nacheinander in wenigen Rechen-operationen gelöst werden kann. Unter den verschiedenen Lösungsverfahren für Gleichungssystemebietet sich hier ein iterativer Löser an. Will man den Rechenaufwand verglichen mit der explizitenEuler-Methode ungefähr konstant halten, kann man genau eine Iteration durchführen:

(1) : x(t+ ∆t) = x(t) + v(t) ·∆t(2) : a(t+ ∆t) aus x(t+ ∆t) berechnen(3) : v(t+ ∆t) = v(t) + a(t+ ∆t) ·∆t(4) : x(t+ ∆t) = x(t) + v(t+ ∆t) ·∆t

In diesem Fall ist die erreichte Genauigkeit aber nicht besser als bei der expliziten Methode, derFehler hat lediglich ein umgekehrtes Vorzeichen. Führt man dagegen mehrere Iterationen durch(d.h. verwendet das in der letzten Zeile berechnete x(t + ∆t) als neues x(t + ∆t) und wiederholtdamit den Algorithmus ab der zweiten Zeile), werden mehrere Kräfte-Berechnungen pro Zeitschrittbenötigt (nämlich eine Kraftberechnung pro Iteration). Da das Berechnen der Kräfte in einerMolekulardynamik-Simulation jedoch sehr teuer ist, nämlich ca. 90% der gesamten Rechenzeit inAnspruch nimmt, ist das keine attraktive Option: Ein Algorithmus, der k Kräfteberechnungen proZeitschritt durchführt, braucht fast k-mal so viel Gesamtrechenzeit wie einer, der pro Zeitschrittnur einmal die Kräfte berechnet.

2.3 Störmer-Verlet-AlgorithmusAus der Taylor-Entwicklung für x(t) lässt sich ein besseres Integrationsverfahren ableiten, nämlichVerlets Algorithmus. Man beginnt mit einer Taylor-Entwicklung in jede Richtung und addiertanschließend beide Gleichungen, wobei alle Terme ungerader Ordnung wegfallen:

x(t+ ∆t) = x(t) + x(t)∆t+1

2x(t)∆t2 +

1

3!x(3)(t)∆t3 +O(∆t4)

x(t−∆t) = x(t)− x(t)∆t+1

2x(t)∆t2 − 1

3!x(3)(t)∆t3 +O(∆t4)

x(t+ ∆t) + x(t−∆t) = 2 · x(t) + x(t) ·∆t2 +O(∆t4)

Daraus ergibt sich nun Verlets Algorithmus:

x(t+ ∆t) = 2 · x(t)− x(t−∆t) + a(t) ·∆t2

wobei wie immer a = x und der zu erwartende Fehler in jedem Zeitschritt vierter Ordnung ist;Verlets Algorithmus ist also lokal gesehen ein Verfahren dritter Ordnung. Betrachtet man dage-gen den globalen Fehler (bei Ausführung mehrerer Zeitschritte nacheinander), treten doch wiederFehlerterme dritter Ordnung auf, die globale Ordnung des Verfahrens ist daher 2.

Verglichen mit Eulers Methode steigt die erwartete Genauigkeit aufgrund der höheren Ord-nung. Ebenfalls steigt der Speicherplatzbedarf, da für jedes Molekül x(t−∆t) gespeichert werdenmuss. Wohlgemerkt muss das Verfahren für den allerersten Zeitschritt einer Simulation angepasstwerden, da x(t −∆t) in diesem Fall natürlich noch nicht existiert; beispielsweise kann man einenrückwärtsgewandten Euler-Schritt durchführen, um die Molekülpositionen zu einem hypothetischen

5

Zeitpunkt t = −∆t zu definieren. Die numerischen Eigenschaften des Störmer-Verlet-Algorithmussind nicht ganz optimal, da zum einen subtrahiert wird und zum anderen der ∆t2-Term bei kleinem∆t sehr klein wird.

Außerdem fällt auf, dass die Geschwindigkeit der Moleküle nicht mehr direkt vorkommt. Oftwird sie allerdings benötigt, z.B. um die Gesamtenergie des Systems zu überwachen oder die Tem-peratur zu kontrollieren. Sie muss aus den Positionen der Moleküle berechnet werden:

v(t) =x(t+ ∆t)− x(t−∆t)

2∆toder

v

(t+

∆t

2

)=

x(t+ ∆t)− x(t)

∆t

Geschwindigkeits-Störmer-Verlet

Man kann Verlets Algorithmus so umformen, dass die Geschwindigkeit der Moleküle explizit ver-wendet wird. Die genauen Schritte der Herleitung seien dem Leser zur Übung überlassen. Nach derUmformung lauten die Formeln wie folgt:

x(t+ ∆t) = x(t) + v(t) ·∆t+a(t)

2·∆t2

v(t+ ∆t) = v(t) +a(t) + a(t+ ∆t)

2·∆t

Bei exakter Arithmetik ist diese Variante identisch zum normalen Verlet-Algorithmus; entspre-chend ist auch der zu erwartende Fehler derselbe. Die Werte für x(t − ∆t) müssen nicht mehrgespeichert werden. Die numerischen Eigenschaften sind etwas besser, da immerhin keine Sub-traktion mehr vorkommt, allerdings werden nach wie vor potenziell sehr kleine ∆t2-Terme zu denKoordinaten addiert.

2.4 Leapfrog-AlgorithmusDer Leapfrog-Algorithmus ist eine weitere arithmetisch äquivalente Umformung des Verlet-Algo-rithmus, ist aber von den Formeln her so einfach wie die Euler-Methode. Der Trick ist, dass Orte undBeschleunigungen der Moleküle auf der einen und Geschwindigkeiten auf der anderen Seite nichtimmer zum gleichen Zeitpunkt berechnet werden, sondern die Geschwindigkeiten um einen halbenZeitschritt versetzt angegeben werden. Die Formeln für Position und Geschwindigkeit lauten:

x(t+ ∆t) = x(t) + v

(t+

∆t

2

)·∆t

v

(t+

∆t

2

)= v

(t− ∆t

2

)+ a(t) ·∆t

Der Beweis, dass diese Formeln sich in diejenigen des (Geschwindigkeits-)Störmer-Verlet-Algo-rithmus umformen lassen, sei wiederum dem Leser überlassen. Benötigt man die Geschwindigkeitenzum Zeitpunkt t, ist es praktikabel, die Aktualisierung des Geschwindigkeitswerts in zwei gleichgroße Schritte aufzuteilen, und v(t) zwischendurch abzulesen:

v(t) = v

(t− ∆t

2

)+

1

2a(t) ·∆t

v

(t+

∆t

2

)= v(t) +

1

2a(t) ·∆t

Der Leapfrog-Algorithmus vereint die Vorteile aller bisher vorgestellten Methoden: Er ist einVerfahren dritter Ordnung wie Verlets Algorithmus, benötigt nur so wenig Speicher wie die Euler-Methode, erfordert keine spezielle Behandlung des ersten Zeitschritts und verwendet weder Sub-traktionen noch ∆t2-Terme.

6

2.5 Predictor-Corrector-VerfahrenEinen ganz anderen Ansatz als Verlets Algorithmus verfolgen die Predictor-Corrector-Verfahren.Der Name selbst bezeichnet eigentlich nur eine recht generelle Idee bzw. eine ganze Familie vondarauf aufbauenden Algorithmen; im Kontext von Molekulardynamik meint man in der Regel dennach seinem Erfinder benannten Gear-Algorithmus. Das grundsätzliche P-C-Schema besteht ausdrei Schritten:

1. Vorhersage (prediction): Aus den bereits bekannten Informationen zu einem gegebenen Zeit-punkt t werden die Molekülpositionen (und ggfs. Geschwindigkeiten) zum Zeitpunkt t+ ∆tgeschätzt.

2. Berechnung (evaluation): Die an den geschätzten Positionen auf die Moleküle wirkendenKräfte sowie die daraus resultierende Beschleunigung werden berechnet.

3. Korrektur (correction): Anhand der berechneten Beschleunigungen werden die geschätztenWerte korrigiert.

Man kann die implizite Euler-Methode als ein sehr einfaches P-C-Verfahren interpretieren. GearsAlgorithmus dagegen verwendet höhere Ableitungen der Molekülpositionen, um bessere Genauig-keit zu erzielen, d.h. er basiert auf der Taylor-Entwicklung von x(t), schneidet diese aber späterab als Verlets Algorithmus. Dementsprechend lautet der Prediction-Schritt (beispielhaft für einenGear-Algorithmus vierter Ordnung):

xP (t+ ∆t) = x(t) + x(t)∆t+ x(t)∆t2

2!+ x(3)(t)

∆t3

3!+ x(4)(t)

∆t4

4!

xP (t+ ∆t) = x(t) + x(t)∆t+ x(3)(t)∆t2

2!+ x(4)(t)

∆t3

3!

xP (t+ ∆t) = x(t) + x(3)(t)∆t+ x(4)(t)∆t2

2!

x(3)P (t+ ∆t) = x(3)(t) + x(4)(t)∆t

x(4)P (t+ ∆t) = x(4)(t)

Im Correction-Schritt wird zunächst der Fehler in der vorhergesagten Beschleunigung xP ver-glichen mit der aus den Kräften am vorhergesagten Ort bestimmten Beschleunigung x berechnetund in einen Korrekturterm c umgewandelt:

∆x = x(t+ ∆t)− xP (t+ ∆t)

c =∆x ·∆t2

2!

Mit Hilfe dieses Korrekturterms werden nun die Vorhersagen korrigiert:

x(t+ ∆t) = xP (t+ ∆t) + α0c

x(t+ ∆t)∆t = xP (t+ ∆t)∆t+ α1c

x(t+ ∆t)∆t2

2!= xP (t+ ∆t)

∆t2

2!+ α2c

x(3)(t+ ∆t)∆t3

3!= x(3)P (t+ ∆t)

∆t3

3!+ α3c

x(4)(t+ ∆t)∆t4

4!= x(4)P (t+ ∆t)

∆t4

4!+ α4c

Dabei sind α0, . . . , α4 ∈ (0, 1] „magische“ Konstanten, die entsprechend der Ordnung des Verfah-rens passend bestimmt werden müssen, um numerische Stabilität zu erhalten. Natürlich gilt dabeiimmer α2 = 1, um die Beschleunigung (deren Soll-Wert ja bekannt ist) vollständig zu korrigieren.Für bereits ermittelte, sinnvolle Werte für gängige Ordnungen sei auf die entsprechende Literaturverwiesen3.

3z.B. [Haile92] S. 162

7

Der essenzielle Unterschied zu den bisher vorgestellten Integratoren besteht im Correction-Schritt, der eine Art Feedback-Mechanismus darstellt. Dadurch kann etwaiges instabiles Verhaltendes Prediction-Schritts gedämpft werden. Wie beim impliziten Euler-Verfahren wäre es theoretischmöglich (und bezüglich der erreichbaren Genauigkeit vorteilhaft), mehrere Evaluation-Correction-Runden durchzuführen, d.h. das Ergebnis eines Correction-Schritts als neue Vorhersage zu verwen-den; aber wiederum würde dies in der Praxis zu viel zusätzlichen Rechenaufwand bedeuten, undempirische Ergebnisse zeigen, dass Gears P-C-Verfahren auch ohne wiederholte Correction-Schrittegut funktioniert.

P-C-Verfahren sind, was theoretische Herleitung, erforderlichen Programmcode und auch Re-chenaufwand betrifft, deutlich aufwändiger als beispielsweise der Leapfrog-Algorithmus. Letzteresfällt aber nicht ins Gewicht, da der eigentliche Integrator ohnehin nur einen kleinen Bruchteilder Gesamtrechenzeit einer Molekulardynamik-Simulation in Anspruch nimmt. Wenn ein Integra-tor höhere Genauigkeit bietet (oder, noch besser, größere Zeitschritte zulässt und somit wenigerKräfteberechnungen erfordert) als ein anderer, so ist es unerheblich, ob er dafür ein paar mehrOperationen benötigt. Die interessante Frage für die empirische Untersuchung ist damit, ob dasP-C-Verfahren dank seiner höhergradigen internen Schätzwerte tatsächlich bessere Ergebnisse lie-fert.

2.6 Runge-Kutta-MethodenRunge-Kutta-Methoden werden gern für numerische Integration eingesetzt, da sie in Bezug auf Ge-nauigkeit und Robustheit sehr gut sind. Je nach Ordnung erfordern sie allerdings mehrere Berech-nungen der auf die Moleküle wirkenden Kräfte pro Zeitschritt, daher sind sie für Molekulardynamikzu aufwändig. Eine Ausnahme ist das Runge-Kutta-Schema erster Ordung: dieses ist identisch mitder expliziten Euler-Methode (und stellt somit eine dritte Möglichkeit dar, selbige herzuleiten).

3 ImplementierungIn MarDyns Code war bislang nur der Leapfrog-Integrator implementiert, und dessen Arbeitsweisewar nicht sauber gekapselt. Um möglichst effizient alternative Integratoren einbinden zu können,war einiges Refactoring nötig.

Die Klasse Molecule, die eigentlich als Datenstruktur für die simulierten Moleküle dient, ent-hielt einige Leapfrog-spezifische Funktionen. Diese wurden in die Leapfrog-Klasse verschoben. ImGegenzug mussten einige neue Zugriffsmethoden eingeführt werden, damit alle Integratoren aufdie für sie jeweils interessanten Eigenschaften der Moleküle zugreifen können.

Die existierende Leapfrog-Klasse wiederum enthielt zum Großteil Integrator-unabhängigen Co-de, der daher in der Superklasse Integrator zusammengefasst wurde und somit von allen Integra-toren genutzt werden kann, ohne dupliziert werden zu müssen. Die Klasse Integrator existiertezuvor schon, fungierte aber lediglich als Interface und enthielt keine eigene Funktionalität.

Schließlich wurde das Framework (Simulation, Domain) so erweitert, dass erstens der zu ver-wendende Integrator an der Kommandozeile als Parameter spezifiziert werden kann und dann zurLaufzeit entsprechend initialisiert wird, und dass zweitens die Gesamtenergie des Systems in jedemZeitschritt berechnet und ausgegeben wird.

Als besondere Schwierigkeit stellte sich heraus, dass an mehreren Stellen in MarDyns Code An-nahmen über das Verhalten und die Arbeitsweise des Integrators getroffen werden, so dass fehler-hafte Ergebnisse eines neuen Integrators nicht unbedingt auf Fehler im neuen Code begründet seinmüssen, sondern sich auch aus falschem Zusammenspiel mit anderen Komponenten von MarDyn er-geben können. Beispielsweise wurden Variablenwerte, die der Leapfrog-Algorithmus nicht benötigt,unerwarteterweise einmal pro Zeitschritt auf 0 zurückgesetzt; oder neu eingeführte Eigenschaftender Molecule-Klasse wurden bei bestimmten Operationen (Herumkopieren des Moleküls wie vonder LinkedCells-Datenstruktur gefordert) nicht berücksichtigt. Eine Reihe solcher Hürden wurdeim Verlauf dieser Arbeit gefunden und beseitigt, aber es ist nicht auszuschließen, dass eventuellekünftig implementierte weitere Integratoren auf ähnliche Hindernisse stoßen werden.

Im Unterschied zu ihrer obigen Beschreibung mussten alle Integratoren aufgeteilt werden in eineMethode upd_preF(), die zu Beginn des jeweiligen Zeitschritts abgearbeitet wird, bevor die neuen

8

Kräfte berechnet wurden, sowie eine Methode upd_postF(), die entsprechend nach der Kräftebe-rechnung ausgeführt wird. Die bisherige Implementierung des Leapfrog-Algorithmus kapselte diesebeiden Schritte in einen Zustandsautomaten, um ihre abwechselnde Bearbeitung in der richtigenReihenfolge sicherzustellen. Dies war in frühen Phasen von MarDyns Entwicklung sicher hilfreichzur Fehlerfindung/-vermeidung, wurde nun aber zwecks Vereinfachung des Codes entfernt.

Speziell zum Predictor-Corrector-Verfahren ist noch anzumerken, dass dieses einmal mit dritterund einmal mit fünfter Ordnung implementiert wurde. Dabei werden intern nicht die Werte von x(t)

usw. verwaltet, sondern direkt die Taylor-Terme x(t)∆t2

2! usw. Die Formeln des Prediction-Schrittslauten dadurch etwas anders, sind aber arithmetisch identisch.

4 Auswertung

4.1 Betrachtung der einzelnen IntegratorenZunächst wird untersucht, wie sich die verschiedenen Integratoren bei unterschiedlichen Zeitschritt-weiten verhalten. Dafür wird mit jedem der Integratoren jeweils dieselbe Simulation eine bestimmte,konstante (simulierte) Zeit T = 3.2 lang berechnet (also in Abhängigkeit von ∆t unterschiedlichviele Zeitschritte lang). Die dabei verwendete Simulation ist nicht frisch initialisiert, sondern hatbereits eine gewisse Zeit damit verbracht, sich in eine natürliche Konstellation zu bewegen (diesenProzess nennt man Equilibrierung). Während dem Test wird in jedem Zeitschritt die Gesamtenergiedes Systems berechnet und ausgegeben. Die folgenden Diagramme zeigen, wie die Gesamtenergiesich im Lauf der Simulation verändert (folglich ist der Wert zu Beginn immer 0). Was die Vergleich-barkeit der Diagramme betrifft, ist zu beachten, dass die Skalierung der y-Achse jeweils passendzu den Daten gewählt ist.

Die voreingestellte Zeitschrittweite für die verwendete Simulation war ∆t = 0.01. Um dasVerhalten der Integratoren zu untersuchen, wurde diese zweimal halbiert und dreimal verdoppelt;getestet wird also mit

∆t ∈ {0.0025, 0.005, 0.01, 0.02, 0.04, 0.08}.

Die Zeit ist eine „dimensionslose Größe“, d.h. sie ist passend zu den übrigen physikalischenGrößen in der Simulation normalisiert4. Übersetzt in reale physikalische Größen sind die Zeit-schrittweiten in der Molekulardynamik-Simulation typischerweise im Femtosekundenbereich. Dergenaue Umrechnungsfaktor ist abhängig vom simulierten Stoff, im Fall von Argon gilt beispielsweise∆t = „0.01“ ≈ 22 · 10−15s.

Die Energie ist ebenso dimensionslos, wiederum im Fall von Argon gilt ∆E = „1“ ≈ 1.65 ·10−21J . Das ist aber nicht weiter von Interesse, da wir primär die Genauigkeit der Integratorenmiteinander vergleichen wollen, ohne uns Gedanken um absolute Werte zu machen.

4Eine detaillierte Erklärung findet sich beispielsweise in [Griebel04], S. 102f.

9

Euler-Methode

0

50

100

150

200

250

300

350

400

450

500

0 0.5 1 1.5 2 2.5 3 3.5

Energ

iediffe

renz

Zeit t

Euler

Δt=0.0025Δt=0.005

Δt=0.01Δt=0.02

Δt=0.04Δt=0.08

Auffällig bei der Euler-Methode ist, dass die Gesamtenergie in allen Experimenten beinaheexakt linear zunimmt. Dabei wird die Zunahme erwartungsgemäß geringer, je kleiner der Zeitschrittgewählt wird. Für noch größere Zeitschritte als hier dargestellt ist die Energiezunahme übrigensnicht mehr annähernd linear.

Störmer-Verlet-Algorithmus

−1

−0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

0 0.5 1 1.5 2 2.5 3 3.5

Energ

iediffe

renz

Zeit t

Störmer−Verlet

Δt=0.0025Δt=0.005

Δt=0.01Δt=0.02

Δt=0.04Δt=0.08

Der Störmer-Verlet-Algorithmus funktioniert erwartungsgemäß deutlich besser als die Euler-Methode. Für kleine Zeitschritte ist die Energieerhaltung recht zufriedenstellend; für ∆t ≥ 0.04wird das Verhalten jedoch deutlich schlechter.

10

Geschwindigkeits-Störmer-Verlet-Algorithmus

−0.18

−0.16

−0.14

−0.12

−0.1

−0.08

−0.06

−0.04

−0.02

0

0.02

0.04

0 0.5 1 1.5 2 2.5 3 3.5

Energ

iediffe

renz

Zeit t

Geschwindigkeits−Störmer−Verlet

Δt=0.0025Δt=0.005

Δt=0.01Δt=0.02

Δt=0.04Δt=0.08

Im direkten Vergleich zum normalen Störmer-Verlet- zeigt der Geschwindigkeits-Störmer-Verlet-Algorithmus, was für numerische Vorteile eine arithmetische Umformung des Algorithmus habenkann. Nicht nur für die groben, sondern auch für die feineren Zeitschritte ist die Rechengenauigkeitoffenbar deutlich besser. Auffällig ist zudem, dass die Zeitschrittweite wesentlich weniger Einflussauf das Ergebnis hat. Die Kurven für ∆t ∈ {0.0025, 0.005, 0.01} sind fast identisch, und auch∆t = 0.02 liegt noch sehr nah an diesen. Für noch gröbere Werte werden zwar die Fehler etwasgrößer, aber wenn man bedenkt, dass der Gesamtrechenaufwand für die Simulation sich mit jederVerdoppelung von ∆t halbiert, sind vielleicht sogar diese noch akzeptabel.

Dass der Fehler nicht mehr so regelmäßig aussieht wie noch bei der Euler-Methode, hat nichtviel zu bedeuten; bestünde der Energieverlauf bei der Euler-Methode aus einer perfekten linearenFunktion plus einer unregelmäßigen Schwankung wie dem hier dargestellten Verlauf, wäre das imDiagramm von einem rein linearen Verlauf nicht zu unterscheiden.

Leapfrog-Algorithmus

−0.14

−0.12

−0.1

−0.08

−0.06

−0.04

−0.02

0

0.02

0.04

0 0.5 1 1.5 2 2.5 3 3.5

Energ

iediffe

renz

Zeit t

Leapfrog

Δt=0.0025Δt=0.005

Δt=0.01Δt=0.02

Δt=0.04Δt=0.08

Die Ergebnisse des Leapfrog-Algorithmus sind denen der Geschwindigkeits-Störmer-Verlet-Methode sehr ähnlich, dementsprechend gelten auch dieselben Kommentare. Für gröbere Zeit-schritte liefert der Leapfrog etwas bessere Ergebnisse; wiederum sind alle Kurven sehr dicht bei-

11

einander. Die Kurve für ∆t = 0.02 liegt so nah an den (beinahe ununterscheidbaren) Kurven fürkleinere Zeitschritte, dass die Empfehlung naheliegt, zumindest für die vorliegende Simulation dieZeitschrittweite zu erhöhen, um den Rechenaufwand zu reduzieren.

Predictor-Corrector (3. Ordnung)

−0.14

−0.12

−0.1

−0.08

−0.06

−0.04

−0.02

0

0.02

0.04

0 0.5 1 1.5 2 2.5 3 3.5

Energ

iediffe

renz

Zeit t

PredictorCorrector−3

Δt=0.0025Δt=0.005

Δt=0.01Δt=0.02

Δt=0.04Δt=0.08

Das Diagramm für das Predictor-Corrector-Verfahren dritter Ordnung ist von dem des Leapfrog-Algorithmus fast nicht zu unterscheiden. Für ∆t ≤ 0.02 sind die Verläufe der Gesamtenergie prak-tisch identisch, für die beiden gröberen Zeitschritte sind die Ergebnisse sogar etwas besser als beimLeapfrog-Algorithmus.

Predictor-Corrector (5. Ordnung)

−1.2

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0 0.5 1 1.5 2 2.5 3 3.5

Energ

iediffe

renz

Zeit t

PredictorCorrector−5

Δt=0.0025Δt=0.005

Δt=0.01Δt=0.02

Δt=0.04Δt=0.08

Das Predictor-Corrector-Verfahren fünfter Ordnung liefert ein – zumindest auf den ersten Blick– deutlich anderes und signifikant schlechteres Bild. Bei genauerem Hinsehen zeigt sich jedoch,dass die Kurven wiederum sehr ähnlich wie bei den drei letztgenannten Algorithmen sind. DerUnterschied ist, dass bei gröberen Zeitschrittweiten der PC-5 Startschwierigkeiten an den Tag legt:In den ersten zwei Zeitschritten treten sehr starke Fehler in der Energieerhaltung auf, danach sinddie Verläufe aber so konstant, wie man sich das wünschen würde. Dies ist recht einfach zu erklären:Wie oben beschrieben, verwaltet das PC-5-Verfahren intern Schätzwerte für die 0. bis 5. Terme

12

der Taylorentwicklung der Molekülkoordinaten. Zu Beginn des Testlaufs sind diese Daten nochnicht vorhanden, da anfangs nur Positionen und Geschwindigkeiten bekannt sind. Während denersten Zeitschritten müssen daher die mit 0 initialisierten Schätzungen der höheren Ableitungenkorrigiert werden, erst danach erreicht das Verfahren seine optimale Leistungsfähigkeit. Würdedie Simulation von Anfang an mit dem PC-5-Verfahren durchgeführt, sähe das Ergebnis dahermöglichweise deutlich besser aus.

4.2 Vergleich aller Integratoren für ausgewählte ZeitschrittweitenUm den Vergleich der Integratoren zu erleichtern, werden nun dieselben Daten anders visualisiert:Jeweils eine feste Zeitschrittweite, dafür aber alle Integratoren in einem Diagramm. Wohlgemerktwurden die Ergebnisse der Euler-Methode nicht in die Diagramme aufgenommen, da diese die y-Achsenskalierung so gestaucht hätten, dass zwischen den übrigen Integratoren keine Unterschiedemehr zu sehen gewesen wären.

−0.2

−0.15

−0.1

−0.05

0

0.05

0 0.5 1 1.5 2 2.5 3 3.5

Energ

iediffe

renz

Zeit t

∆t=0.0025

Störmer−VerletGeschwindigkeits−Störmer−Verlet

LeapfrogPredictorCorrector−3PredictorCorrector−5

Für den kleinsten untersuchten Zeitschritt, ∆t = 0.0025, liefern die Geschwindigkeits-Störmer-Verlet-, Leapfrog-, Predictor-Corrector-3- und Predictor-Corrector-5-Algorithmen praktisch gleicheErgebnisse, während der Störmer-Verlet-Algorithmus ein gutes Stück schlechter abschneidet.

−0.35

−0.3

−0.25

−0.2

−0.15

−0.1

−0.05

0

0.05

0 0.5 1 1.5 2 2.5 3 3.5

Energ

iediffe

renz

Zeit t

∆t=0.01

Störmer−VerletGeschwindigkeits−Störmer−Verlet

LeapfrogPredictorCorrector−3PredictorCorrector−5

Mit ∆t = 0.01 fällt das Ergebnis ein kleines bisschen differenzierter aus: Leapfrog und PC-3teilen sich den ersten Platz, während Geschwindigkeits-Störmer-Verlet und PC-5 knapp dahinter

13

liegen; letzterer allerdings nur wegen seiner Startschwierigkeiten. Verlets Algorithmus ist deutli-cher abgeschlagen als zuvor und bildet wiederum das Schlusslicht (von der im Diagramm nichtenthaltenen Euler-Methode abgesehen).

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0 0.5 1 1.5 2 2.5 3 3.5

Energ

iediffe

renz

Zeit t

∆t=0.04

Störmer−VerletGeschwindigkeits−Störmer−Verlet

LeapfrogPredictorCorrector−3PredictorCorrector−5

Für ∆t = 0.04 ergibt sich dieselbe Rangfolge. Das PC-5-Verfahren scheint weiter zurückzufallen,was aber nur an einem größeren Fehler in den ersten beiden Zeitschritten liegt. Nimmt man denWert im dritten Zeitschritt als Basis und betrachtet den weiteren Verlauf der Energieerhaltung, istdie (auf diese Weise geschönte) Kurve für PC-5 praktisch identisch mit der für PC-3.

4.3 Extra-lange ZeitschritteDie bisher untersuchten Zeitschrittweiten sind für reale Simulationsanwendungen vernünftig, wiedie in Kapitel 4.1 besprochenen Ergebnisse zeigen. Da die besten der untersuchten Integratoren dortaber sehr ähnliche Ergebnisse liefern, wurden sie zusätzlich unter etwas widrigeren Bedingungen,sprich gröberen Zeitschritten, getestet, um ihre Robustheit besser beurteilen zu können. Wiederumwurde mit verschieden feinen Zeitschritten (∆t ∈ {0.1, 0.2, 0.4}) eine konstante Gesamtzeit (T =64) simuliert. Auf die Euler-Methode wurde wie bei den bisherigen Vergleichen schon verzichtet,da sie um Größenordnungen schlechter arbeitet als die übrigen Verfahren.

−5

0

5

10

15

20

0 10 20 30 40 50 60 70

Energ

iediffe

renz

Zeit t

∆t=0.1

Störmer−VerletGeschwindigkeits−Störmer−Verlet

LeapfrogPredictorCorrector−3PredictorCorrector−5

Das Ergebnis für ∆t = 0.1 ist nicht unerwartet: Der Störmer-Verlet-Algorithmus führt zu relativstarken Energieschwankungen, das Predictor-Corrector-5-Verfahren hat Startschwierigkeiten, die

14

verbleibenden drei Integratoren liegen sehr dicht beieinander. Um ihre Differenzen besser sehenzu können, entfernen wir für den nächstgröberen Zeitschritt die Kurve für den Störmer-Verlet-Algorithmus aus dem Diagramm:

−7

−6

−5

−4

−3

−2

−1

0

1

0 10 20 30 40 50 60 70

Energ

iediffe

renz

Zeit t

∆t=0.2

Geschwindigkeits−Störmer−VerletLeapfrog

PredictorCorrector−3PredictorCorrector−5

Für ∆t = 0.2 kann der Leapfrog-Algorithmus leicht in Führung gehen, wobei das PC-5-Verfahren wegen der sehr großen Fehler in den ersten zwei Zeitschritten wieder einmal sehr schwerzu beurteilen ist.

Um die Predictor-Corrector-Verfahren angemessen bewerten zu können, werden für die gröbsteZeitschrittweite ihre Fehler in den ersten Zeitschritten ignoriert, und statt dessen für PC-3 derEnergiewert nach dem 1. Zeitschritt, für PC-5 derjenige nach dem 3. Zeitschritt als Basiswertverwendet. Dass diese Manipulation legitim ist, lässt sich wie folgt begründen: Die Initialisierungdes Integrationsalgorithmus aus den wenigen zu Beginn des Testlaufs vorhandenen Daten (nämlichnur Koordinaten und Geschwindigkeiten der Moleküle) stellt eine besondere Schwierigkeit dar, dienicht repräsentativ für die Leistungsfähigkeit eines Verfahrens „auf offener Strecke“ ist. Ignoriertman also die ersten paar Zeitschritte für diejenigen Algorithmen, die offenbar einige Zeit brauchen,bis sie ihre maximale Ergebnisqualität erreichen, so erhält man Beobachtungen, die aussagekräftigersind für das typische Verhalten des Algorithmus während einer längeren Simulation.

−4

−3.5

−3

−2.5

−2

−1.5

−1

−0.5

0

0.5

1

0 10 20 30 40 50 60 70

Energ

iediffe

renz

Zeit t

Δt=0.4

Geschwindigkeits−Störmer−VerletLeapfrog

PredictorCorrector−3 ab Zeitschritt 1Predictor−Corrector−5 ab Zeitschritt 3

Damit zeigt sich für die besten vier hier besprochenen Integratoren bei der gröbsten geteste-ten Zeitschrittlänge folgendes Bild: Das Predictor-Corrector-5-Verfahren erreicht die beste Ener-gieerhaltung, rechnet also anscheinend am genauesten. Die einfachere Variante PC-3 zeigt ähn-

15

liche durchschnittliche Änderungen pro Zeitschritt, aber mehr langfristigen Drift. Der Leapfrog-Algorithmus scheint zwar zumindest für diesen Simulationslauf im langfristigen Mittel mit demPC-5-Verfahren mitzuhalten, führt aber zu wesentlich stärkeren kurzfristigen Schwankungen derEnergiemenge. Der Geschwindigkeits-Störmer-Verlet-Algorithmus ist von diesen vier Verfahren dasschlechteste, da sowohl die hochfrequenten Schwankungen als auch der langfristige Drift vergleichs-weise hoch sind.

5 FazitDer Leapfrog-Algorithmus als bisher einziger in MarDyn implementierter Integrator ist eine guteWahl und liefert insbesondere für üblicherweise verwendete Zeitschrittweiten sehr gute Ergebnis-se. Eine signifikante Verbesserung durch Verwendung eines anderen Integrators ist hier nicht zuerwarten.

Gleichwohl ist das Predictor-Corrector-Verfahren von Grad 3 dem Leapfrog äquivalent bis leichtüberlegen, allerdings auch etwas aufwändiger zu implementieren. Anhand der vorliegenden Ergeb-nisse fällt es schwer, zwischen diesen beiden Alternativen eine konkrete Empfehlung auszusprechen.

Das PC-Verfahren mit Grad 5 braucht zwei bis drei Zeitschritte, während denen es relativ großeFehler macht, um seine internen Werte korrekt zu initialisieren, liefert danach aber hervorragendeErgebnisse. Insbesondere für extrem lange Zeitschritte (die je nach Anwendungsfall möglicherweisegrundsätzlich nicht in Frage kommen, da sie zu grob wären) ist das PC-5-Verfahren allen ande-ren hier besprochenen Algorithmen deutlich überlegen. Ist man bereit, auf etwas Genauigkeit zuverzichten, ist es im Interesse eines möglichst geringen Gesamtrechenaufwandes daher eine nahe-liegende Idee, das PC-5-Verfahren und eine sehr lange Zeitschrittweite zu wählen.

Die Euler-Methode ist deutlich ungenauer als die ausgefeilteren Algorithmen. Sie hat lediglicheinen Vorteil: da sie extrem einfach zu implementieren ist, ist sie sehr unanfällig für Bugs (sei es inihrem eigenen Code oder im Zusammenspiel mit den übrigen Teilen des Simulationsprogramms)und eignet sich daher gut, um die grundsätzliche Funktion des Programms zu verifizieren, sowie alsuntere Grenze für den Erwartungswert, wie gut ein anderer Algorithmus mindestens funktionierenmuss, wenn er korrekt implementiert ist.

Um noch aussagekräftigere Ergebnisse zu erhalten, könnte man noch weitere Simulationendurchführen; insbesondere könnte man erstens verschiedene Ausgangssituationen verwenden imHinblick auf Molekülanzahl, -Dichte, -Geschwindigkeit und -Masse, sowie zweitens die Simulationdeutlich länger laufen lassen, um die langfristige Stabilität der Algorithmen besser beurteilen zukönnen. Aufgrund von beschränkter zur Verfügung stehender Zeit und Prozessorleistung ist diesim Rahmen dieser Arbeit nicht erfolgt.

Bevor einer der alternativen Integratoren voll einsatzbereit ist, müsste er zudem für die Dreh-bewegungen nicht-kugelförmiger Moleküle erweitert werden. Es ist zu erwarten, dass die anhandder Translation gewonnenen Beobachtungen sich auf die Rotation übertragen.

Literatur[Griebel04] Griebel, M., Knapek, S., Zumbusch, G., Caglar, A.: Numerische Simulation in der

Moleküldynamik. Springer 2004

[Haile92] Haile, J.M.: Molecular Dynamics Simulation: Elementary Methods. Wiley 1992

[Frenkel96] Frenkel, D., Smit, B.: Understanding Molecular Simulation. From Algorithms to App-lications. Academic Press / Elsevier, 1996

16