Studentische Forschung an der Fakultät für Mathematik


Approximative Moleküldynamik mit transformierten Phasenraumdichten

Bachelorarbeit (SS 2011) und Masterarbeit (WS 2012/13) von Johannes Keller
Betreuerin: Prof. Dr. Caroline Lasser

Überblick: In diesem Artikel stelle ich den Inhalt meiner Bachelor- und Masterarbeiten vor, die ich im Sommersemester 2011 beziehungsweise im Wintersemester 2012/13 geschrieben habe. Nachdem ich für mein Bachelorstudium das Nebenfach Wirtschaft gewählt hatte, kam ich erst durch mein Hauptseminar im 5. Semester näher mit mathematischer Physik und Numerik in Berührung. Im Anschluss an dieses Seminar habe ich meine Bachelorarbeit geschrieben.

Ich habe Schrödingers kanonische Formulierung der molekularen Quantendynamik mit einer alternativen Formulierung verglichen, die in der theoretischen Chemie sehr verbreitet ist, da sie einfachere statistische Interpretationen zulässt. Während ich in meiner Bachelorarbeit gute Ergebnisse für den zeitunabhängigen Vergleich der beiden Formulierungen herleiten konnte, waren die analytischen und numerischen Ergebnisse für die Zeitentwicklung noch nicht zufriedenstellend.

Nach meinem Bachelorabschluss bewarb ich mich erfolgreich für die Aufnahme in das Promotionsprogramm von TopMath. Mit Frau Prof. Lasser als meiner Betreuerin durfte ich weiter an dem Thema forschen, und schließlich konnten wir gemeinsam einen Approximationssatz für die Zeitentwicklung in der alternativen Formulierung beweisen. Aus diesem Satz haben wir einen neuen Algorithmus entwickelt, mit dem man quantenphysikalische Größen in der Moleküldynamik schnell und genau berechnen kann. Der Approximationssatz und der Algorithmus bilden den Kern meiner Masterarbeit.

Wir haben einen wissenschaftlichen Aufsatz, basierend auf den Forschungsergebnissen aus meiner Masterarbeit und weiteren numerischen Experimenten, zur Publikation in der Fachzeitschrift SIAM Journal of Applied Mathematics eingereicht. Die Publikation wurde im Sommersemester 2013 zur Veröffentlichung angenommen.

Seit Juli 2012 habe ich außerdem an der TU München eine Forschungsstelle in dem Projekt ,,Potential Energy Surfaces'' im Rahmen des Sonderforschungsbereichs/Transregio 109 "Discretization in Geometry and Dynamics". Im Rahmen dieses Projekts promoviere ich über Diskretisierungen von molekularer Quantendynamik und untersuche insbesondere die internuklearen Kräfte in Molekülen. Meine Arbeit im Rahmen der Bachelor- und Masterprojekte hat mich gut auf dieses Forschungsgebiet vorbereitet, in dem viele verschiedene mathematische Disziplinen - von Hamiltonschen Systemen und Geometrischer Integration bis hin zu mikrolokaler Analysis und nichtkommutativer Geometrie - zusammentreffen.

Einführung: Motiviert von einer numerischen Problemstellung aus dem Bereich der Moleküldynamik und theoretischen Chemie habe ich in meiner Bachelorarbeit Verfahren zur Berechnung von Integralen der Form $$ \int_{\mathbb{R}^{2d}} W(q,p) \cdot a(q,p)~dq~dp$$ untersucht. Wenn \( W \) eine Wahrscheinlichkeitsdichte ist, dann können wir das obige Integral mit Monte-Carlo Integration auswerten. Dafür erzeugen wir Punkte, die gemäß der Dichte \( W \) verteilt sind, und bilden dann das arithmetische Mittel über die Auswertungen von \( a \) auf diesen Punkten. Formal lässt sich diese Quadraturmethode also folgendermaßen schreiben:$$ \int_{\mathbb{R}^{2d}} W(q,p) \cdot a(q,p)~dq~dp~ \approx ~\frac{1}{n} \sum_{j=1}^n a(z_j) ~~~~~~{\rm mit}~~~~~z_1,\ldots,z_n \stackrel{i.i.d.} {\sim}W~.$$Ich habe mich mit sogenannten Wignerfunktionen \( W \) beschäftigt, die auch negativ werden dürfen, also keine echten Wahrscheinlichkeitsdichten mehr sind. Allerdings haben diese Funktionen die Eigenschaft, dass sie echte Wahrscheinlichkeitsdichten (das heißt, nichtnegativ) werden, wenn man sie mit einer zentrierten Gaussdichte \( G_{\varepsilon/2} \) mit Varianz \( \frac{\varepsilon}{2} \) (\( \varepsilon \) ist klein) faltet. Diese Dichte hat die Form $$ G_{\varepsilon/2}(q,p) = (\pi\varepsilon)^{-d}\cdot e^{-\frac{1}{\varepsilon}(|q|^2 + |p|^2)}~~~~.$$ Die einzigen Wignerfunktionen \( W \), die schon von sich aus nichtnegativ (also Wahrscheinlichkeitsdichten) sind, sind Gaussdichten. Wenn \( W \) also zum Beispiel eine Gaussdichte mit Varianz \( \frac{\varepsilon}{2} \) ist, so ist \( H:=G_{\varepsilon/2} \star W \) (die sogenannte Husimi-Transformierte) eine Gaussdichte mit Varianz \( \varepsilon \).

Ein weiteres, aus der Moleküldynamik motiviertes Beispiel für eine Wignerfunktion \( W \), die lokal negativ ist, ist eine normierte Summe ("Überlagerung") von zwei Gaussdichten mit Varianz \( \frac{\varepsilon}{2} \) und einer mit einer Kosinusfunktion multipliziertem Gaussdichte.

wigner5-zugeschnitten-klein.jpg
Abb. 1 (zum Vergrößern klicken): Überlagerung von zwei Gaussdichten und einer Gaussdichte mit oszillierendem Vorfaktor.
Den Graphen einer solchen Funktion für den Fall \( d=1 \) und \( \varepsilon=0.1 \) kann man in Abbildung 1 sehen, das mit \( G_{\varepsilon/2} \) gefaltete Gegenstück in Abbildung 2. In Abbildung 2 kann man außerdem sehr gut sehen, dass die Oszillationen von \( W \) sehr stark geglättet wurden, was neben der garantierten Nichtnegativität ein weiterer positiver Effekt des Übergangs von \( W \) zu \( H \) ist.

Es ist nun die Frage, wie groß der Fehler für den Wert des Integrals in Abhängigkeit von \( \varepsilon \) ist, wenn wir die Wignerfunktion \( W \) im Integral durch die Husimifunktion \( H=G_{\varepsilon/2} \star W \) (also eine echte Wahrscheinlichkeitsdichte) ersetzen, sodass wir wieder Monte-Carlo Integration machen können.

Mit Hilfe einer geeigneten Taylorentwicklung kann man für diese Strategie des Ersetzens von \( W \) durch \( H \) einen Fehler der Ordnung \( \mathcal{O}(\varepsilon) \) berechnen. Ich habe in meiner Bachelorarbeit die explizite Fehlerkorrektur erster Ordnung $$ \int_{\mathbb{R}^{2d}} W(q,p) \cdot a(q,p)~dq~dp = \int_{\mathbb{R}^{2d}} H(q,p) \cdot \left( a(q,p)-\frac{\varepsilon}{4}\Delta a(q,p) \right)~dq~dp + \mathcal{O}(\varepsilon^2) ~~~~~(*) $$ hergeleitet, also ist die Fehlerabschätzung für alle nichtlinearen Funktionen \( a \) scharf. Korrekturen höherer Ordnung von \( a \) bestehen aus Termen der Form \( ~C(k)\cdot\varepsilon^k\Delta^k a \).


husimi-zugeschnitten-klein.jpg
Abb. 2 (zum Vergrößern klicken): Die mit \( G_{\varepsilon/2} \) geglättete Dichte der Überlagerung aus Abbildung 1.
Hintergrund: Integrale des obigen Typs treten auf, wenn man Erwartungswerte von quantenphysikalischen Observablen, also zum Beispiel den Impuls oder die kinetische Energie eines Moleküls, berechnet. Der Zustand eines Moleküls wird in der Quantendynamik durch eine Funktion \( \psi:\mathbb{R}^{3 N}\to \mathbb{C} \) beschrieben, wobei \( N \) die Anzahl der Atomkerne im betrachteten Molekül ist (im Folgenden setzen wir \( 3\cdot N =d \)). Observablen sind lineare Funktionen auf dem Raum aller Zustände.

Der sogenannte Erwartungswert \( \langle {\hat{A}} \rangle_{\psi} \) einer quantenphysikalischen Observablen \( {\hat{A}} \) im Zustand \( \psi \) lässt sich mit Hilfe der Identität $$ \langle {\hat{A}} \rangle_{\psi}=\int_{\mathbb{R}^{d}}\bar{\psi}(q)\cdot {\hat{A}}\cdot \psi(q)~dq = \int_{\mathbb{R}^{2d}} W_{\psi}(q,p) \cdot A(q,p)~dq~dp ~~~~~~~~~~~~~~(**)$$ berechnen, wobei die Funktion \( W_{\psi} \) eine Wignerfunktion ist. Man kann \( W_{\psi} \) aus dem Zustand \( \psi \) des Systems berechnen, die Funktion \( A:{\mathbb{R}}^{2n}\to{\mathbb{R}} \) erhält man aus dem Operator \( \hat{A} \) (die sogenannte Weyl-Korrespondenz).

Ersetzen wir nun in der Gleichung \( (**) \) \( W_{\psi} \) durch die Wahrscheinlichkeitsdichte \( H_{\psi} \) und berechnen wir das Integral mit Monte-Carlo Integration, so berechnen wir damit also den quantenphysikalischen Erwartungswert bis auf Fehler der Ordnung \( \mathcal{O}(\varepsilon) \). In meiner Bachelorarbeit habe ich mehrere numerische Experimente gemacht, in denen man sieht, dass diese numerische Methode zur Berechung von Erwartungswerten wirklich Fehler der Ordnung \( \mathcal{O}(\varepsilon) \) erzeugt.

Vertiefung: Die zeitliche Entwicklung eines Molekülzustands \( t\mapsto\psi(t) \) wird in der Quantendynamik durch die semiklassische Schrödingergleichung, eine lineare partielle Differentialgleichung, beschrieben. Diese Gleichung ist numerisch allerdings nur in sehr niedrigen Dimensionen (\( N\leq3 \)), und selbst dann nur mit sehr großem Aufwand lösbar.

In der Chemie ist man jedoch meistens weniger an der zeitlichen Entwicklung des Molekülzustandes \( \psi(t) \) selber, als vielmehr an der zeitlichen Entwicklung von Erwartungswerten einer Observablen \( \hat{A} \) interessiert. Umgeschrieben mit Identität \( (**) \) will man also die Abbildung $$ t\mapsto \int_{\mathbb{R}^{2d}} W_{\psi(t)}(q,p) \cdot A(q,p)~dq~dp=\langle {\hat{A}} \rangle_{\psi(t)} $$ numerisch auswerten.

Ausgehend von meiner Bachelorarbeit habe ich im Rahmen meiner Masterarbeit einen neuen Approximationsansatz für diese Abbildung entwickelt. Ich habe numerisch einfach lösbare gewöhnliche Differentialgleichungen hergeleitet, deren Lösungsfluss \( F^t:\mathbb{R}^{2d}\to \mathbb{R}^{2d} \) die statische Fehlerordnung \( \mathcal{O}(\varepsilon^2) \) aus \( (*) \) in der Zeitentwicklung erhält. Es gilt $$\langle \hat A\rangle_{\psi(t)}= \int_{\mathbb{R}^{2d}} H_{\psi(0)}(q,p) \cdot \left( a-\frac{\varepsilon}{4}\Delta a \right)(F^t(q,p))~dq~dp + \mathcal{O}(\varepsilon^2).$$ Aus diesem Ergebnis habe ich einen Algorithmus zur näherungsweisen Berechnung der Dynamik von molekularen Erwartungswerten entwickelt, für den man nicht die Schrödingergleichung, sondern nur wesentlich einfacher zu behandelnde gewöhnliche Differentialgleichungen lösen muss. Numerische Experimente, die ich im Rahmen meiner Masterarbeit durchgeführt habe, bestätigen diese hohe Genauigkeit, die für viele Anwendungen schon ausreichend ist.

Link zur vollständigen Bachelor-Arbeit (500 KB) und Master-Arbeit (700 KB) von Johannes Keller und zum arXiv-Preprint Pfeil.

Zurück zur Übersicht über studentische Forschung




Steckbrief des Autors

keller-bild-120px.jpg

  • Johannes Keller
  • Student an der TU seit WS 2008/09
  • B.Sc. in Mathematik mit Nebenfach Wirtschaft, SS 2011
  • nach dem Bachelor: Promotionsphase im Studienprogramm TopMath
  • M.Sc. in Mathematik, WS 2012/13
  • Forscher im Projekt Potential Energy Surfaces des SFB Discretization in Geometry and Dynamics (seit Sommer 2012)
  • Mathematische Interessen: numerische Moleküldynamik, Funktionalanalysis, mathematische Chemie
  • Webseite (M3, Wissenschaftl. Rechnen)
 

TUM Mathematik Rutschen TUM Logo TUM Schriftzug Mathematik Logo Mathematik Schriftzug Rutsche

picture math department

Impressum  |  Datenschutzerklärung  |  AnregungenCopyright Technische Universität München