Molekulardynamik Franziska Dilger 27. Januar 2011

Inhaltsverzeichnis Einleitung

1

1

2

2

3

4

Einführung klassische Molekulardynamik

1.1 Historisches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Hard-Sphere-Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 N-Körper-Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Potentiale und Kraftauswertung

2.1 Kurzreichweitige Potentiale . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Langreichweitige Potentiale . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Algorithmen zur ezienten Kraftberechnung . . . . . . . . . . . . . . 2.3.1 Linked-Cell Verfahren für kurzreichweitige Potentiale . . . . . . 2.3.2 Barnes-Hut Verfahren und schnelles Multipolverfahren für langreichweitige Potentiale . . . . . . . . . . . . . . . . . . . . . . . Zeitintegrationsverfahren

3.1 Anforderungen . . . . . . . . . . . . . . 3.2 Verwendete Verfahren . . . . . . . . . . 3.3 Fehler der Zeitintegration . . . . . . . . 3.3.1 lokaler Fehler . . . . . . . . . . . 3.3.2 globaler Fehler . . . . . . . . . . 3.3.3 Experiment Lyapunov Stabilität Ausblick

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . . . . . . .

2 2 2

4

4 6 6 7 7

10

10 11 11 11 12 12

13

4.1 Parallelisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4.2 Thermostate und Ensembles . . . . . . . . . . . . . . . . . . . . . . . . . 13

Franziska Dilger

Molekulardynamik

TUM

Einleitung Computersimulationen spielen heutzutage in beinahe allen Zweigen der Wissenschaft eine wichtige Rolle. Mithilfe von Computersimulationen können Experimente ergänzt und schrittweise ersetzt werden. So ist es möglich, Experimente unter Extrembedingungen zu simulieren, die im Labor nur schwer herzustellen sind. Betrachtet man Computersimulationen im Molekulardesign so führt dies zur Molekulardynamik (MD). Im folgenden beschränke ich mich auf die klassische Molekulardynamik. In der klassischen Molekulardynamik werden quantenmechanische Aspekte vernachlässigt, die Teilchendynamik wird somit nicht durch die Schrödinger-Gleichung beschrieben. Das Lösen von Schrödinger-Gleichung ist bei hoher Dimension nur für sehr kleine Systeme numerisch möglich, analytisches Lösen ist nicht möglich und ihre Lösungen liefern lediglich Wahrscheinlichkeitsverteilungen für die Postition und Geschwindigkeit der Teilchen. Die klassische MD basiert auf den Newtonschen Bewegungsgleichungen. Es werden Wechselwirkungen zwischen Partikeln, sogenannten interagierenden Teilchen wie zum Beispiel Atomen, Molekülen, Sterne oder Teile von Galaxien betrachtet und über einen kurzen Zeitraum hinweg dargestellt. Damit lässt sich die zeitliche Veränderung von Position, Geschwindigkeit und Orientierung beschreiben. Die Partikel tragen die Eigenschaften von physikalischen Teilchen wie Masse, Geschwindigkeit, Ort und Ladung. Man betrachtet die Energieerhaltung der kinetischen Energie der Teilchen und der potentiellen Energie des Systems aufgrund der Wechselwirkungen. Hauptsächlich verwendet wird die MD Simulation zur Untersuchung biologischer Moleküle, deren Struktur und Dynamik. In den Materialwissenschaften werden mithilfe der MD Simulation neue Stoe entwickelt, aber auch in der Biologie, Chemie, Physik und Astrophysik spielt die MD eine wichtige Rolle. Gerade die rasche Entwicklung paralleler Rechensysteme macht die Molekulardynamik immer populärer.

1

Franziska Dilger

Molekulardynamik

TUM

1 Einführung klassische Molekulardynamik

1.1 Historisches

Da die historische Entwicklung der Computersimulation von Partikelmodellen eng mit der Entwicklung des Computers verbunden ist, bewegte man sich in der Vergangenheit immer an der Grenze des Machbaren. Der erste Artikel über MD Simulationen erschien im Jahre 1957 von Alder und Wainwright. Sie entwickelten ein Modell mit einigen hundert Partikeln, harten Kugeln, die untereinander elastische Stöÿe ausführen. Durch Phasendiagramme auf Basis ihrer Ergebnisse endeckten sie die Existenz eines Phasenwechsels von üssig nach fest. 1964 untersuchte Rahman die Eigenschaften von üssigem Argon. Er verwendete dabei das Lennard-Jones Potential, erstmals ein realistisches Potential. Verlet stellte 1967 die Nachbarschaftslisten vor, ein Verfahren zur ezienten Verwaltung der Daten während einer MD Simulation und entwickelte ein Verfahren zur Zeitintegration, das heute noch Standard in der MD Simulation ist. 1974 führten Stillinger und Rahman erste MD Simulationen von Wasser durch und drei Jahre später, 1977 erschien der SHAKE-Algorithmus. Josh Barnes und Piet Hut entwickelten 1986 den Barnes-Hut Algorithmus für langreichweitige Potentiale. In den 90er Jahren wurde die MD Simulation dann auf parallelen Prozessoren durchgeführt, dadurch wurde es möglich gröÿere Teilchensysteme und längere Simulationzeiten zu betrachten.

1.2 Hard-Sphere-Model

Im Modell harter Kugeln werden zur Vereinfachung die Teilchen im zweidimensionalen Raum als ÿchwebende runde Kreiseünd im dreidimensionalen Raum als ÿchwebende runde Kugeln"mit festem Radius betrachtet. Dabei führt eine eher dichte Anordnung in einem Kristall zu einem Festkörper. Will man Gase simulieren, so wählt man die durchschnittlichen Teilchenabstände sehr groÿ. Die Atome von Flüssigkeiten sind ohne feste Anordnung, aber im Vergleich zu Gasen eher dicht angeordnet. Dieses Modell ist nicht ohne Schwächen. So sind Moleküle meist keine Kugeln. Betrachtet man zum Beispiel Kohlenwasserstoe, beispielsweise Decan so sieht man, dass Moleküle auch sehr in die Länge gezogen sein können. Ein weiterer Nachteil ist, dass in diesem Modell quantenmechanische Eekte vernachlässigt werden. Im Ausgleich kann man eektive Potentiale betrachten.

1.3 N-Körper-Problem

Gegeben seien N Partikel in einem Simulationsgebeit Ω ⊆ Rd . Jedes Teilchen besitzt dabei eine bestimmte Masse mi ∈ R, einen Ort xi ∈ Rd , Geschwindigkeit vi ∈ Rd und damit einen Impuls pi = mi vi ∈ Rd . Jeder dieser Partikel hat ein bestimmtes Potential, besitzt also ein Kraftfeld. Die Kräfte der einzelnen Partikel überlagern sich zu einem gemeinsamen Kraftfeld Fij , das wiederum auf die einzelnen Körper wirkt.

2

Franziska Dilger

Molekulardynamik

TUM

Die Partikel erfahren damit eine Beschleunigung. Das Kraftfeld ist somit abhängig von der Position der Partikel und die Position und Bewegung der Partikel ist jeweils abhängig von dem Gesamtkraftfeld. Betrachte man nun die Dynamik eines Teilchens mithilfe der Newtonschen Bewegungsgleichung i = 1, ..., N Fi = −∇xi V (x1 , ..., xN ) P wobei die Kraft auf ein Teilchen Fi = j6=i Fij aus dem negativen Gradienten des Potentials V (x1 , ..., xN ) berechnet wird. Schreibt man diese Gleichung nun um in mi x¨i = Fi

Hamilton Formalismus so erhält man die Hamiltonschen Dierentialgleichungen x˙ i = ∇pi H

p˙i = −∇xi H

mit der Hamilton Funktion H(x1 , ..., xN , p1 , ...pN ) =

N X p2i + V (x1 , ..., xN ) 2mi

.

k=1

Und damit ein System von 2dN Dierentialgleichungen erster Ordnung: x˙ i = m−1 i pi

p˙i = −∇xi V (x1 , ..., xN )

3

Franziska Dilger

Molekulardynamik

TUM

2 Potentiale und Kraftauswertung Ein Potential U (rij ) wird beschrieben durch ein skalares Feld mit Fij = −∇U (rij ). In der Physik ist das Potential die Fähigkeit eines Kraftfeldes Fij Arbeit zu verrichten. Die potentielle Energie lässt sich berechnen durch V (x1 , ...xN ) =

N N X X

U (rij )

rij := kxj − xi k

i=1 j=1,j>i

Die potentielle Gesamtenergie Vtot ergibt sich durch Summation der im System beteiligten potentiellen Energien: Vtot = Vbond + Vangle + ... + VvdW

Im Allgemeinen muss immer unterschieden werden zwischen einem kurzreichweitigen, schnell abfallenden Potential und einem langreichweitigen, langsam abfallenden Potential. Denition: Ein Potential heiÿt kurzreichweitig, wenn es für d > 2 mit zunehmendem Abstand rij zwischen zwei Teilchen schneller als r−d abfällt, andernfalls heiÿt das Potential langreichweitig.

2.1 Kurzreichweitige Potentiale

Man betrachte hier zunächst Wechselwirkungen zwischen gebundenen Atomen. Das sind Potentiale für die die Kraftberechnung relativ einfach ist und in O(N ) liegt, da pro Atom nur eine konstante Anzahl von Nachbarn betrachtet werden muss. Modelliert durch das Hook'sche Gesetz, der Vorstellung, dass Atompaare durch elastische Federn verbunden sind, ist die Bindungslängendeformation Vbond gegeben durch Vbond =

X (i,j),i>j

1 bij (∆rij )2 2

∗ mit bij als experimentell bestimmte Kopplungskonstante und ∆rij = rij − rij . Ebenso betrachte man die Bindungswinkeldeformation Vangle

Vangle =

X (i,j,k),i>k

1 wijk (∆Θijk )2 2

mit Bindungswinkel Θijk , ∆Θijk = Θijk − Θ∗ijk mit Θ∗ijk als natürlicher Bindungswinkel und wijk als Kraftkonstante. Als klassisches Beispiel für ein kurzreichweitiges Potential zwischen nicht gebundenen Atomen nun das Van-der-Waals Potential mit potentieller Energie VvdW . Dieses Potential modelliert elektrostatische Wechselwirkungen zwischen zwei ungeladenen Atomen oder Molekülen. Die Elektronen in der Hülle sind leicht beweglich, dadurch 4

Franziska Dilger

Molekulardynamik

TUM

ergeben sich temporäre Ladungsverschiebungen und man erhält temporäre Dipole. Die Kräfte haben eine geringe Reichweite, da sich die Dipolladungen im Fernfeld schnell auslöschen. Das Van-der-Waals Potential ist gegeben durch 

1 UvdW (rij ) = −a rij

6

Dieses Potential tritt als anziehender Anteil im Lennard-Jones-(12,6) Potential auf. Das Lennard-Jones-(n,m) Potential beschreibt ebenfalls Wechselwirkungen zwischen ungeladenen, nicht chemisch aneinander gebundenen Atomen. Es hat allerdings zusätzlich einen repulsiven Potentialanteil, der bei Annäherung zweier Teilchen überwiegt. Bei Annäherung zweier Teilchen überlappen sich die Elektronenwolken und die positiven Kerne wechselwirken miteinander. Für groÿe Entfernungen zwischen zwei Teilchen überwiegen die anziehenden Kräfte. Das Lennard-Jones-(n,m) Potential ist gegeben durch  ULJ (rij ) = αε

σ rij

n

 −

σ rij

m ! m rcut für das Potential gleich Null gesetzt. Wenn man nun die Seitenlänge der Zellen ≥ rcut wählt, so kann man sich bei der Kraftauswertung auf die Partikel benachbarter Zellen beschränken.

In der rechten Abbildung sehen wir das Simulationsgebiet mit den unterteilten Zellen im Fall von periodischen Randbedingungen. MD-Simulationen nden meist unter periodischen Randbedingungen statt. Jedes Teilchen, das das Simulationsgebiet auf einer Seite verlässt, taucht auf der gegenüber liegenden wieder auf und wie man in der Abbildung sieht, nden alle Wechselwirkungen auch über diese Grenzen hinweg statt. Damit gilt für die Kraftberechnung Fi ≈

X

X

Fij

Zelle kc, j∈{P artikel der Zelle kc} kc∈N (ic) j6=i

Der Aufwand für die Kraftberechnung verringert sich also mit dem Linked-Cell Verfahren auf O(N ). 2.3.2 Barnes-Hut Verfahren und schnelles Multipolverfahren für langreichweitige Potentiale

Das Barnes-Hut Verfahren für langreichweitige Potentiale gehört zu der Klasse der Baumverfahren und ist das älteste, einfachste und meist verbreiteste Verfahren für 7

Franziska Dilger

Molekulardynamik

TUM

langsam abfallende Potentiale. Dabei verringert man die Anzahl der zu berechnenden Kräfte zwischen zwei Teilchen durch geeignetes Zusammenfassen von Teilchengruppen zu Pseudoteilchen. Hauptsächlich wird dieses Verfahren in der Astrophysik verwendet, wobei dort das Gravitationspotential betrachtet wird. Die Idee des Verfahrens ist, die Kraft einer Partikelgruppe auf einen Einzelpartikel durch die Wirkung einer einzelnen Masse, dem Pseudoteilchen zugehörig, im Massenschwerpunkt der Partikelgruppe anzunähern. Der Algorithmus untergliedert sich in drei Schritte. Im ersten Schritt werden die Partikel im Simulationsgebiet in eine hierarchische Struktur, im zweidimensionalen in einen Quadtree, im dreidimensionalen in einen Octtree einsortiert. Man zerlegt das Simulationsgebiet solange rekursiv in gleich groÿe Teilgebiete, sogenannte Zellen, bis höchstens ein Partikel pro Zelle vorhanden ist. Die Partikel stehen dann in den Blättern des Baumes. Da sich die Partikel im System mit der Zeit bewegen, muss der Baum für jeden Zeitschritt neu berechnet werden. Die Komplexität ist also abhängig von der Anordnung und der Dichte der Partikel.

Im zweiten Schritt des Algorithmus wird die Masseverteilung berechnet. Bis zu diesem Zeitpunkt enthält der Baum nur Informationen über die räumliche Verteilung der Teilchen. Man berechne jetzt die Gesamtmasse aller unter einem Knoten enthaltenen Teilchen, sowie den Massenschwerpunkt für alle Knoten des Baumes. Die Knoten des Baumes repräsentieren somit die, mit einer bestimmten Masse und einem Ort, gegeben durch den Massenschwerpunkt der Partikelgruppe, bestückten Pseudopartikel des Systems. Im dritten Schritt werden nun die Kräfte berechnet. Dieser Teil des Algorithmus ist für die Ezienzierung zuständig. Die Kraft mehrerer Teilchen j auf ein Teilchen i wird approximiert durch die Kraft des Pseudopartikels auf das Teilchen i, wobei das Pseudoteilchen die Zelle repräsentiert in dem die Teilchen j liegen. Damit diese Näherung zulässig ist, muss der Abstand r der Gruppe zum Einzelteilchen groÿ sein und der Gruppendurchmesser bzw. der Zellendurchmesser d im Verhältnis zum Abstand klein sein. Es muss also gelten, dass Θ≤

d r

Θ dr so muss der Quadrant erneut unterteilt werden und das Θ-Kriterium muss erneut für alle Knoten überprüft werden, anderfalls werte man die Kraft wie oben beschrieben aus. Der Aufwand der Kraftberechnung für langreichweitige Potentiale verringert sich mit dem Barnes-Hut Verfahren auf O(N log(N )). In dem Barnes-Hut Verfahren werden Patikel-Partikel-Wechselwirkungen durch PartikelCluster-WW ersetzt. Im Schnellen Multipolverfahren will man diese Strategie nochmals weiter führen und ersetzt nun die Partikel-Cluster-WW durch Cluster-ClusterWW. Als Cluster werden die Zellen mit ihren enthaltenen Partikeln bezeichnet, die durch die Pseudopartikel repräsentiert werden. Für alle Partikel eines Clusters können aus den Cluster-Cluster-Wechselwirkungen die Partikel-Cluster-Wechselwirkungen abgeleitet und berechnet werden.

Durch dieses Vorgehen reduziert sich der Aufwand der Kraftberechnung nochmals auf O(N ) und man ist optimal, da für jede Kraftberechnung jeder Partikel mindestens einmal in Betracht gezogen werden muss.

9

Franziska Dilger

Molekulardynamik

TUM

3 Zeitintegrationsverfahren Für die eziente Integration gewöhnlicher Dierentialgleichungen denkt man sofort an Verfahren höherer Ordnung, es spielen aber zusätzliche Eigenschaften eine wichtige Rolle. Viele davon wurden bereits in vorherigen Vorträgen behandelt.

3.1 Anforderungen

Zum einen muss das gewählte Integrationsverfahren exakt sein, da die numerisch berechnete Trajektorie der Partikel nach einem Zeitschritt so wenig wie möglich von der exakten Trajektorie der Partikel abweichen soll. Weitere einzuhaltende Eigenschaften sind die Stabilität des Verfahrens und die Energieerhaltung. Bei Hamiltonschen Systemen die nicht explizit zeitabhängig sind wird die Energie entlang der Trajektorie der Partikel erhalten. Die numerisch berechnete Trajektorie kann von der exakten Trajektorie abweichen und einen leichten Drift in der Energie aufweisen. Symplektische Verfahren besitzen ausgezeichnete Eigenschaften der Energieerhaltung. Auÿerdem sollte das Verfahren wenn möglich Reversibilität aufweisen, das bedeutet in Abwesenheit von numerischen Rundungsfehlern soll man die berechnete Trajektorie durch Vorzeichenwechsel in der Zeit exakt zurücklaufen können. Ganz wichtig ist auch wie bei der Kraftauswertung, dass das Verfahren rechnerisch ezient ist. Es sollte unter anderem auch möglich sein, groÿe Schrittweiten zu wählen ohne dabei die Stabilität zu verlieren. Da die Kräfteberechnung in jedem Zeitschritt auf einen einzelnen Partikel durchgeführt werden muss, ist daurauf zu achten, dass nur eine Kraftauswertung pro Zeitschritt im Zeitintegrationsverfahren fällig ist. Implizite Verfahren sind demnach ausgeschlossen, da viele Kraftauswertungen bei der iterativen Lösung des nichtlinearen Gleichungssystem anfallen. Wie bereits erwähnt, kann man um die Zeitintegration zu beschleunigen, Iterationsverfahren höherer Ordnung verwenden. Durch die bessere Konvergenzordnung wird der Fehler der Simulation verkleinert, wir erhalten einen besseren Stabilitätsbereich und können damit gröÿere Zeitschritte verwenden. In der Molekulardynamik hat man im Gesamtpotential durch die unterschiedlichen Zeitskalen der Bindungs-,Winkel-, Lennard-Jones und Coulombanteile verschiedene Zeitskalen. Am besten verwendet man Multiple Zeitschrittverfahren, trennt also die Potentiale auf und verwendet jeweils unterschiedliche Zeitschrittweiten. Durch Einfrieren der Bindungs- und Winkelkräfte lassen sich direkt gröÿere Zeitschrittweiten verwenden, ohne die Stabilität des Verfahrens zu gefährden. Man xiert somit hochfrequente Moden durch das Erfüllen zusätzlicher Nebenbedingungen. Bekannt ist dieses Verfahren auch unter dem Namen Shake- und Rattle Verfahren.

10

Franziska Dilger

Molekulardynamik

TUM

3.2 Verwendete Verfahren

Das meist verwendete Verfahren zur Integration der Newtonschen Bewegungsgleichung ist das Störmer-Verlet Verfahren. xn+1 = xn + dtvn vn+1 = vn +

dt F (xn+1 ) m

Es werden aber auch Mehrschrittverfahren, Splittingverfahren, Kompositionsverfahren, Partitionierte Runge-Kutta Verfahren, Multiple Zeitschrittverfahren und Verfahren mit Zwangsbedingungen wie der Rattle Algorithmus verwendet. Bei den Splittingverfahren zählt die Idee des Operator-Splittings. Man zerlegt die Hamiltonfunktion in zwei Teile, wendet auf diese getrennt voneinander Integrationsverfahren an und führt die Teilverfahren dann hintereinander aus. Damit erhält man dann wieder ein Verfahren erster Ordnung. Bei den Kompositionsverfahren gewinnnt man aus symplektischen Verfahren zeitreversible. Aus einem Integrationsverfahren Ψ entsteht durch Umkehrung der Zeit t ein adjungiertes Verfahren Ψ∗ . Durch Hintereinanderausführung der Verfahren Ψ∗ ◦ Ψ Ψ ◦ Ψ∗ kann man Zeitintegrationsverfahren höherer Ordnung basteln. Symplektische Runge-Kutta Verfahren sind meist aufwendig und im allgemeinen nicht explizit. Durch Aufspaltung der Hamiltonfunktion bekommt man symplektische partitionierte Runge-Kutta Verfahren, die explizit sein können. Multiple Zeitschrittverfahren werden verwendet, da die gebundenen und ungebundenen Atome auf unterschiedlichen Zeitskalen wechselwirken. Die Bindungslängen und Bindungswinkel oszillieren sehr schnell, es sind sehr kleine Zeitschritte nötig. Das Zeitintegrationsverfahren muss also aufgesplittet werden für die verschiedenen Potentiale. Beim Rattle Algorithmus werden die Freiheitsgrade der Winkel- und Bindungspotentiale eingefroren, da der Energietransport zwischen hoch- und niederfrequenten Freiheitsgraden sehr langsam ist. Man betrachtet hier also die Längen von chemischen Bindungen oder die Gröÿen der Bindungswinkel in Molekülen als xiert.

3.3 Fehler der Zeitintegration 3.3.1 lokaler Fehler

Für Einschrittverfahren der Ordnung p zur Zeitintegration der Newtonschen Bewegungsgleichung gilt bei geeigneter Glattheit der Lösung (x(t), p(t)) k Ψ(xn , pn , h) − Φ(xn , pn , h) k= O(hp+1 )

wobei Ψ die numerische
und Φ die exakte Verfahrensfunktion ist mit n) Φ(xn , pn , h) = x(h+t . p(h+tn )

Damit gilt für den lokalen Fehler in einem Zeitschritt:

k xn+1 − x(h + tn ) k= O(hp+1 )

11

Franziska Dilger

Molekulardynamik

TUM

k pn+1 − p(h + tn ) k= O(hp+1 )

Die Newtonschen Bewegungsgleichunge sind in einem Zeitschritt h somit sehr genau berechenbar. Die Genauigkeit lässt sich zusätzlich durch noch kleinere Schrittweiten h oder durch Integrationsverfahren höherer Ordnung verbessern. 3.3.2 globaler Fehler

Betrachtet man ein festes Zeitintervall von t0 bis tend = t0 + nh mit M als LipschitzKonstante der Verfahrensfunktion Ψ und (x(t), p(t)) als exakte Lösung durch (x0 , p0 , t0 ), so gilt für den globalen Fehler der Zeitintegration von t0 bis tend : k xn − x(tend ) k≤ C · hp ·

eM (tend −t0 ) − 1 M

k pn − p(tend ) k≤ C · hp ·

eM (tend −t0 ) − 1 M

Der Fehler wächst exponentiell in der Zeit tend an. In der Molekulardynamik hat man immer das Interesse an numerischer Simulation über lange Zeitintervalle mit tend groÿ im Verhältnis zur Lipschitz-Konstanten M . Eine kleine Störung der Anfangsdaten wird also exponentiell in der Zeit tend verstärkt. Man erhält ein chaotisches Verhalten. 3.3.3 Experiment Lyapunov Stabilität

Zur Untersuchung des exponentiellen Fehlerwachstums betrachte man zwei einfache Simulationen mit jeweils tausend Partikeln und mit unterschiedlichen Anfangsbedingungen zweier Partikel. Man erwartet dann, dass sich diese kleine Störung der Anfangsbedingung eines Partikels in der Simulation exponentiell in der Zeit verstärkt.

Wie zu erwarten war, verstärkt sich der Änfangsfehlerëines Partikels im Ort um 0.001 in der Dierenz der Trajektorien der zwei Simulationen exponentiell. → Partikelsimulation in C, Fehlerberechnung in Matlab s. les

12

Franziska Dilger

Molekulardynamik

TUM

4 Ausblick Die MD-Simulation ist eine etablierte und mittlerweile sehr robuste Simulationstechnik. Es gibt ständig neue Anwendungsgebiete da sich die Algorithmen immer verbessern, die Prozessorgeschwindigkeiten immer besser werden und sich aufgrund der Parallelisierung der Programme immer mehr Teilchen über längere Zeiträume hinweg simulieren lassen.

4.1 Parallelisierung

Um die Gesamtrechenzeit zu verkürzen werden die zu leistenden Berechnungen auf mehrere Prozessoren verteilt und somit bis zu einem bestimmten Grad gleichzeitig ausgeführt. Die leistungsfähigen Prozessoren werden zusammengeschalten, rechnen unabhängig voneinander aber es ndet Kommunikation zwischen ihnen statt. Bei der Aufteilung der Rechenarbeit auf mehrere Prozessoren ist darauf zu achten, dass man die Kommunikation zwischen den Prozessoren so gering wie möglich hält. Durch aufwendige Kommunikation zwischen den Prozessoren verliert man viel von der gewonnenen Leistung durch die Parallelisierung. Man muss somit bei der Kraftberechnung eine sinnvolle Abbildung der Zellen, in denen die Partikel liegen, auf die Prozessoren nden, sodass so geringe Grenzächen wie möglich zwischen den verschiedenen Prozessoren existieren. Diese Abbildungen nennt man raumfüllende Kurven.

4.2 Thermostate und Ensembles

Ein statistisches Ensemble ist die Menge der möglichen Zustände eines Systems. Das meist verwendete Ensemble ist das (kanonische) NVT- Ensemble, was nicht viel mehr aussagt, als dass wir ein System betrachten in dem die Zahl der Moleküle N , das Volumen des Simulationsgebiet V und die Temperatur T konstant gehalten werden soll. Die kinetische Energie ist durch die Geschwindigkeit der Moleküle gegeben. Ekin =

1X mi vi2 2 i

Die Temperatur T ist deniert über T =

2 Ekin 3N kB

13

Franziska Dilger

Molekulardynamik

TUM

wobei kB die Boltzmann-Konstante ist. Bei Fluktuationen während der Simulation reguliert ein Thermostat die Temperatur (kinetische Energie). Dafür gibt es verschiedene Verfahren wie beispielsweise das NoséHoover Thermostat oder das Berendsen Thermostat.

14

Franziska Dilger

Molekulardynamik

Literatur [1] Michael Griebel, Stephan Knapek, Gerhard Zumbusch, Attila Caglar Numerische Simulation in der Molekulardynamik Springer, 2004 [2] Peter Deuhard, Folkmar Bornemann Numerische Mathematik 2 [3] Vorlesung Algorithms of Scientic Computing 2 Winter 10 www5.in.tum.de

15

TUM