z e−x2dx geschrieben
Der erste Term wird noch in den Raum der reziproken Gittervektoren k transfor-miert, so dass man schließlich
EC = 1
erh¨alt (V: Volumen der Elementarzelle). Der zweite Term ergibt sich durch die Ber¨ucksichtigung i 6=j f¨ur n= 0 im reziproken Raum. Der Parameter β (Einheit:
m−1) bestimmt das Verh¨altnis zwischen reziprokem und direktem Anteil bei der Berechnung vonEC, und wird so gew¨ahlt, dass beide Teile schnell konvergieren. Mit abnehmendemβ nimmt dabei der im direkten Raum bestimmte Anteil zu, f¨urβ = 0 geht Gleichung (4.10) in Gleichung (4.7) ¨uber (erfc(0) = 1).
Bei der Particle-Mesh-Ewald-(PME)-Methode [130, 131] werden die Ladungs-punkte zur Berechnung des reziproken Anteils n¨aherungsweise durch LadungsLadungs-punkte auf einem ¨aquidistanten Gitter ersetzt; damit wird es m¨oglich, diesen Anteil durch schnelle Fourier-Transformation (FFT) auszurechnen.
Bei der praktischen Anwendung von PME in GROMACS wird der Wert von β nicht vom Anwender angegeben, sondern stattdessen eine Cut-Off-Distanz rcoulomb, bis zu der die Berechnung im direkten Raum erfolgen soll (siehe oben).
Der Parameter β wird dann programmintern so bestimmt, dass der Abschirmfaktor erfc(β|rij+n|) bei dieser Entfernung kleiner als der Eingabeparameterewald rtol ist, also erfc(β rcoulomb) < ewald rtol, so dass der direkte Anteil f¨ur gr¨oßere Abst¨ande vernachl¨assigt werden kann. Somit bestimmt ewald rtol die gew¨unschte Genauigkeit der Wechselwirkung im direkten Raum.
4.2 Molekulardynamik
Wir kommen nun zur n¨aherungsweisen L¨osung der Newtonschen Bewegungsglei-chungen Fi =miai; die Kr¨afte Fi werden dabei vom vorher besprochenen Kraftfeld bestimmt.
Orte
Zeit
Geschwindigkeiten
t
r(t) r(t+ ∆t)
t+ 3∆t/2 t+ ∆t
t−∆t/2
v(t−∆t/2) v(t+ ∆t/2) v(t+ 3∆t/2)
t+ ∆t/2
Abbildung 4.2:Leap-Frog-Algorithmus zur Integration der Bewegungsgleichungen, nach [128].
4.2.1 Der Leap-Frog-Algorithmus
GROMACS benutzt den Leap-Frog-Algorithmus, zu deutsch manchmal auch Bock-sprung-Algorithmus genannt, zur Integration der Bewegungsgleichungen [128, 132].
Durch Taylor-Entwicklung nach ∆t2 mit Abbruch nach dem Glied zweiter Ordnung lassen sich folgende Gleichungen herleiten (Indexi weggelassen):
v(t+ ∆t
2 ) = v(t− ∆t
2 ) + F(t)
m ∆t+O(∆t3) (4.11) r(t+ ∆t) = r(t) +v(t+ ∆t
2 )∆t+O(∆t3) (4.12) Wie man sieht, werden Orte und Geschwindigkeiten zu verschiedenen Zeitpunk-ten berechnet, was auch namensgebend war, vgl. Abbildung 4.2. Der Leap-Frog-Algorithmus zeichnet sich anderen Methoden gegen¨uber durch seinen minimalen Rechen- und Speicheraufwand aus.
4.2.2 Nebenbedingungen
Da die Gr¨oße des Zeitschritts ∆tumso kleiner gew¨ahlt werden muss, je h¨oher die ma-ximale Frequenz des Systems ist, ist es sinnvoll, hochfrequente Freiheitsgrade, das heißt Streckschwingungen, einzufrieren, was einen etwa zwei bis viermal l¨angeren Zeitschritt zul¨asst. Solche Nebenbedingungen k¨onnen auf verschiedene Weise im-plementiert werden, in dieser Arbeit wurde der SETTLE-Algorithmus [133] f¨ur die Wassermolek¨ule und LINCS (LINear Constraint Solver) [134] f¨ur den Rest verwen-det. Die Auswirkungen von Abstandsfixierungen wurden in der Literatur f¨ur Gleich-gewichtssituationen untersucht, wobei sich herausgestellt hat, dass die berechneten physikalischen Eigenschaften dadurch nicht verzerrt werden [135,136]. Bindungswin-kel d¨urfen hingegen nicht ohne weiteres festgehalten werden [136].
Der LINCS-Algorithmus funktioniert, kurz gesagt, wie folgt (Abbildung 4.3):
Zun¨achst wird die Bindungs¨anderung ohne Nebenbedingung berechnet; dann werden
4.2 Molekulardynamik
d d
d
d
Abbildung 4.3: LINCS-Verfahren zur Ber¨ucksichtigung von Nebenbedingungen bez¨uglich der Bindungsl¨ange, nach [134]. d bezeichnet die urspr¨ungliche L¨ange der Bindung, die Pfeile deuten die Verschiebungsvektoren der beiden Atome an. Zur Erl¨auterung der Schritte, siehe Text.
die Anteile der Verschiebungsvektoren entlang der alten Bindung herausprojiziert, und schließlich wird nochmals entlang dieser Richtung gestaucht, um die durch die Rotation der Bindung verbleibende L¨angen¨anderung zu kompensieren.
F¨ur starre Wassermolek¨ule, wie Single-Point-Charge-(SPC)-Wasser [137], exi-stieren analytische Gleichungen, um die feste Molek¨ulgeometrie in einem Schritt zu ber¨ucksichtigen (SETTLE-Algorithmus). Das SPC-Wasser-Modell ist gekennzeich-net durch einen O–H-Abstand von 1,0 ˚A und einer Ladung von−0,82ebzw. +0,41e auf dem O-Atom bzw. den H-Atomen. Der H–H-Abstand betr¨agt 1,6333 ˚A (die H–
Atome sind somit in Richtung der Ecken eines Tetraeders mit O im Schwerpunkt ausgerichtet). Van-der-Waals-Wechselwirkungen treten nur mit dem Sauerstoffatom auf, wof¨ur ein Lennard-Jones-Potential verwendet wird.
4.2.3 Temperatur- und Druckkopplung
Aufgrund numerischer Fehler, Unzul¨anglichkeiten des Kraftfelds, etc., kann es im Laufe einer MD-Simulation zu einem unphyskalischen Aufheizen/Abk¨uhlen des Sy-stems kommen. Um dem entgegenzuwirken ist es notwendig, die Temperatur des Systems zu kontrollieren.
Die ”instantane“ Temperatur ˜T eines N-Teilchensystems kann man definieren
¨
uber 1
2NFkT˜=Ekin, (4.13)
mit NF als Anzahl der Freiheitsgrade und Ekin = 1
2 XN
i=1
miv2i. (4.14)
Nach dem Gleichverteilungssatz gilt n¨amlich im Gleichgewicht, dass hEkini= 1
2NFkT, (4.15)
wobei hEkini der Ensemble-Mittelwert der kinetischen Energie im thermodynami-schen Gleichgewicht ist. Wie man sieht, geht Gleichung (4.13) in (4.15) ¨uber, wenn man den Mittelwert von hT˜i uber ein Ensemble (wie es beispielsweise eine MD-¨ Simulation liefert) berechnet, und mit der klassischen Temperatur T gleichsetzt.
Insofern kann Gleichung (4.13) als eine m¨ogliche Verallgemeinerung des Tempe-raturbegriffs angesehen werden; im folgenden wird nicht mehr zwischen ˜T und T unterschieden.
F¨ur die Regelung der Temperatur eines MD-Systems gibt es verschiedene M¨oglichkeiten. In dieser Arbeit wurde der sogenannte Berendsen-Thermostat [138]
verwendet. Dabei wird vorausgesetzt, dass die zeitliche Temperatur¨anderung pro-portional zur Differenz zwischen Soll-Temperatur T0 und Ist-Temperatur T(t) ist (schwache Kopplung):
dT(t)
dt = T0−T(t) τT
. (4.16)
Das heißt, die Temperaturabweichung verschwindet exponentiell mit einer Zeitkon-stanteτT. Die Nachregelung der Temperatur erfolgt durch Skalierung der Geschwin-digkeiten mit einem gemeinsamen Faktorλ, f¨ur den sich durch Einsetzen von Glei-chung (4.13) in (4.16)
λ =
Analog wie die Temperatur l¨asst sich auch der Druck des Systems durch schwache Kopplung nachregeln. Dies geschieht entsprechend zur Temperaturkopplung durch eine Skalierung der Simulationsbox und somit aller Ortskoordinaten [128].
Zur Temperaturkopplung ist noch zu sagen, dass man in Systemen wie
” Pro-tein in Wasser“ oder ¨ahnlichem, das ProPro-tein und das L¨osungsmittel getrennt an das W¨armebad der Temperatur T0 koppelt, das heißt, zwei verschiedene Parameter λP und λLM verwendet. Der Grund daf¨ur ist, dass wegen nicht perfekter Kopplung zwischen Protein und L¨osungsmittel, Cut-Off-Effekten, etc., der Energieaustausch zwischen Protein und L¨osungsmittel nicht gleichm¨aßig ist. Koppelt man ein solches System insgesamt an ein Temperaturbad, wird im allgemeinen ein Abk¨uhlen des Proteins und ein Aufheizen des L¨osungsmittels beobachtet. Der Temperaturunter-schied liegt dabei typischerweise um die 100 K (f¨ur eine Simulation bei ca. 300 K).