1 / 45

V5 Moleküldynamik

V5 Moleküldynamik. “Es gibt “nur” zwei Probleme bei MD-Simulationen: das Kraftfeld und ausreichendes Sampling .” Pflicht V2 Elektrostatik V3 Kraftfelder V4 Energieminimierung + Sampling V6 Ensembles und Thermodynamik Heute (V5): MD- Algorithmen

butch
Download Presentation

V5 Moleküldynamik

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. V5 Moleküldynamik • “Es gibt “nur” zwei Probleme bei MD-Simulationen: • das Kraftfeld und ausreichendes Sampling.” • Pflicht • V2 Elektrostatik • V3 Kraftfelder • V4 Energieminimierung + Sampling • V6 Ensembles und Thermodynamik • Heute (V5): • MD- Algorithmen • - Behandlung der langreichweitigen elektrostatischen WW • - Simulation bei konstanter Temperatur • Kür • V12 Freie Energie-Rechnungen,Berechnung von Bindungskonstanten • V13 Anwendung von MD-Simulationen zur Proteinfaltung Computational Chemistry

  2. Der Gasdruck (I) Impuls p ist das Produkt aus Masse m mal Geschwindigkeit v: p = m ∙ v Kraft F ist das Produkt aus Masse m mal Beschleunigung A: F = m ∙ a, also ist die Kraft = Impuls pro Zeit. Wenn ein Gasmolekül mit dem Impuls m∙vx auf die Wand trifft, wird es mit diesem Impuls in negativer Richtung elastisch reflektiert. Damit wird auf die Wand ein Gesamtimpuls von 2m∙vx ausgeübt. Für einen Würfel der Kantenlänge l ergibt sich die Kraft Fxzu Impuls pro Zeit: 2l / v ist die Zeit, die das Gasmolekül für die Strecke 2l benötigt.. Computational Chemistry

  3. Der Gasdruck (II) Druck P ist definiert als Kraft F pro Flächeneinheit A : P = F / A Die Kraft pro Flächeneinheit bzw. der Druck Px ist dann mit l2 als der Fläche und V = l3 als dem Volumen des Kastens. Gleiches gilt entsprechend für die y- und z-Flächen des Kastens. Wenn der Kasten nun N Moleküle enthält, so ergibt sich mit dem Mittelwert (Erwartungswert) der quadratischen Geschwindigkeiten, wobei Computational Chemistry

  4. Der Gasdruck (III) Da die Bewegungen der Moleküle in jede Richtung zufällig und unabhängig voneinander sind, wird auch der Mittelwert der quadratischen Geschwindigkeiten in jede Richtung identisch sein: Damit ist auch der Druck an allen Wänden gleich. Durch Umformulieren von erhalten wir für den Gesamtdruck, bzw. Computational Chemistry

  5. Der Gasdruck (IV) Die mittlere Geschwindigkeit der Moleküle kann also nur noch von der Temperatur und nicht vom Druck oder Volumen abhängen. Die kinetische Energie der Moleküle wird daher durch den Druck wiedergegeben: pro Molekül, und pro Mol, mit der Avogadrozahl N = 6.022∙1023 mol-1 Der Vergleich mit liefert Oder umformuliert Mit der allgemeinen Gaskonstante R = 8.314 JK-1 mol-1 Computational Chemistry

  6. Freiheitsgrade der Bewegung (I) Für ein Edelgas (einatomig) kommen nur translatorische Bewegungen als Ursache für die kinetische Energie in Frage. Teilen wir die gesamte kinetische Energie auf jede der drei Raumrichtungen (x,y,z) auf, so erhalten wir pro Richtung, sprich Freiheitsgrad pro Mol, bzw. pro einzelnem Molekül mit der Boltzmannkonstante k = R / N = 1.3807∙10-23 J K-1 Computational Chemistry

  7. Freiheitsgrade der Bewegung (II) Für mehratomige Moleküle gibt es noch weitere Freiheitsgrade: 3 Freiheitsgrade der Rotation (im allgemeinen) Symmetrische und anti-symmetrische Schwingungen/Vibrationen 3 Freiheitsgrade der Vibration (im allgemeinen) Computational Chemistry

  8. Freiheitsgrade der Bewegung (III) Ein N-atomiges Moleküle hat damit (in der Regel) 3N - 6 Freiheitsgrade der Bewegung. (Die 3 Translationen und die 3 Rotationen des Schwerpunkts sind invariante Symmetrieoperationen). Bei Raumtemperatur (T = 298 K) sind meist nur die Freiheitsgrade der Translation und Rotation (durch die zur Verfügung stehende thermische Energie) angeregt. Durch entsprechende Energiezufuhr (z.B. IR-Spektroskopie) können Freiheitsgrade selektiv angeregt werden. Computational Chemistry

  9. Anwendungsbereich von Newton-Dynamik Die Newtonsche Physik verteilt die Energie gleichmäßig über alle vibrationelle Moden (Gleichverteilungs-Theorem), wogegen die Eigenzustände des harmonischen Oszillator jedoch gequantelte Energieabstände h haben. MD Simulationen sind daher gut geeignet für Simulationen bei hohen Temperaturen (z.B. Raumtemperatur) und zum Sampling von Moden mit Moden mit relativ großen Frequenzen wie Bindungsschwingungen (siehe Tabelle) sind bei Raumtemperatur jedoch fast immer im Grundzustand. Daher ist klassische MD nicht gut geeignet um die Vibrationen solcher Bindungen und Bindungswinkeln zu simulieren. Hierzu sind quantenmechanische Methoden notwendig. [Schlick] Computational Chemistry

  10. Maxwell-Boltzmannverteilung der Geschwindigkeiten Wenn es eine mittlere Geschwindigkeit der Moleküle gibt, wie sieht dann die Verteilung der Geschwindigkeit aus ? Bei höherer Temperatur verschiebt sich die mittlere Geschwindigkeit daher zu höheren Geschwindigkeiten. Computational Chemistry

  11. Boltzmannverteilung der Energien Analog zu den Geschwindigkeiten sind auch die möglichen Energieniveaus Ei für alle N Teilchen je nach Temperatur unterschiedlich besetzt. Diese Verteilung wird als Boltzmannverteilung bezeichnet, und gilt für Systeme, in denen ein bestimmtes Energieniveau beliebig oft besetzt werden kann. Für welche Teilchen gilt dies nicht ? Computational Chemistry

  12. Moleküldynamik (MD) MD-Simulation ist eine der Hauptmethoden für die Untersuchung biologischer Moleküle (Struktur, Dynamik und Thermodynamik) und ihrer Komplexe. MD berechnet das zeitliche Verhalten eines molekularen Systems. MD-Simulationen enthalten detaillierte Informationen über Fluktuationen und Konformationsänderungen von Proteinen und Nukleinsäuren. MD-Simulationen werden auch zur Strukturaufklärung in Röntgenkristallographie und NMR-Experimenten eingesetzt. Computational Chemistry

  13. Dynamische Prozesse in Biomolekülen • Lokale Bewegungen (0.01 to 5 Å, 10-15 to 10-1 s) • Atomare Fluktuationen • Bewegungen der Seitenketten • Bewegung der Loops • Bewegung starrer Körper (1 to 10 Å, 10-9 to 1s) • Bewegung von Helices • Domänenbewegung (hinge bending) • Bewegung von Untereinheiten • Weiträumige Bewegungen (> 5Å, 10-7 to 104 s) • Helix coil Übergänge • Dissoziation/Assoziation • Faltung und Entfaltung Computational Chemistry

  14. Historischer Hintergrund: immer am Abgrund des gerade noch „Machbaren“ 1957/1959: Einführung der MD-Methode für die Wechselwirkung von harten Kugeln als Modell von Flüssigkeiten (Alder & Wainwright) 1964 erste Simulation mit realistischem Potential für flüssiges Argon (Rahman, 1964) 1974 erste MD-Simulation von Wasser (Stillinger and Rahman, 1974) 1977 erste MD-Simulation des Proteins BPTI (McCammon, et al, 1977). 1977 SHAKE-Algorithmus 1983 CHARMM, AMBER-Kraftfeld 1984 Freie-Energie-Rechnungen 1985 Car-Parrinello MD 1987 Gromos87-Kraftfeld seit 1990 QM/MM Methode, z.B. für enzymatische Reaktionen 1992 Ewald-Methode ´90er Jahre: MD-Simulation auf parallelen Prozessoren  größere Systeme Alder, B. J. and Wainwright, T. E. J. Chem. Phys.27, 1208 (1957) Alder, B. J. and Wainwright, T. E. J. Chem. Phys.31, 459 (1959) Rahman, A. Phys. Rev. A136, 405 (1964) Stillinger, F. H. and Rahman, A. J. Chem. Phys.60, 1545 (1974) McCammon, J. A., Gelin, B. R., and Karplus, M. Nature (Lond.)267, 585 (1977) Computational Chemistry

  15. Ein einfaches Moleküldynamik-Programm program md einfaches MD program call init Initialisierung t = 0 do while (t.lt.tmax) MD Schleife call force(f,en) berechne die Kräfte call integrate(f,en) integriere Bewegungsgleichungen t = t + delta call sample berechne Mittelwerte enddo stop end nach [Frenkel, Smit 1996] Computational Chemistry

  16. Initialisierung eines MD-Programms subroutine init Initialisierung sumv=0 sumv2=0 do i=1,npart x(i)= lattice_pos(i) platziere Teilchen auf Gitter (oder lese Koordinate ein) v(i)=(ranf()-0.5) ordne Zufallsgeschwindigkeiten zu sumv=sumv+v(i) Geschwindigkeit des CMS sumv2=sumv2+v(i)**2 kinetische Energie enddo sum=sumv/npart Geschwindigkeit des CMS sumv2=sumv2/npart mittlere quadrierte Geschwindigkeit fs=sqrt(3*temp/sumv2) Skalierungsfaktor für Geschwindigkeiten – warum? do i=1,npart setze gewünschte kinetische Energie und setze v(i)=(v(i)-sumv)*fs - Geschwindigkeit des Massenschwerpunkts = 0 xm(i)=x(i)-v(i)*dt Position bei vorherigem Zeitschritt enddo - (für Integrationsalgorithmus) return end nach [Frenkel, Smit 1996] Computational Chemistry

  17. periodische Randbedingungen (PBC) Verwendung von PBC erlaubt es, Simulationen mit relativ kleiner Teilchenzahl durchzuführen, wobei die Teilchen Kräfte wie in Lösung erfahren. Die Koordinaten der „Kopien“ gehen aus denen der Teilchen in der zentralen Box durch einfache Translationen hervor. Kräfte auf die zentralen Teilchen werden durch Wechselwirkungen innerhalb der Box sowie mit den Nachbarboxen berechnet. Der cut-off wird so gewählt, daß ein zentrales Teilchen nicht gleichzeitig mit einem anderen zentralen Teilchen und dessen Kopie wechsel-wirkt. Für eine kubische Box: cut-off < box/2 Bsp. 2-dimensionale Box. Zentrale Box ist von 8 Nachbarn umgeben. nach [Frenkel, Smit 1996] Computational Chemistry

  18. Wiederholung: Molekülmechanik-Kraftfeld Computational Chemistry

  19. Berechnung der Kräfte subroutine force(f,en) Initialisierung en=0 do i=1,npart f(i)= 0 setze Kräfte auf 0 enddo do i=1,npart-1 do j=i+1,npart Schleife über alle Paare xr=x(i)-x(j) berechne Distanz zwischen Teilchen i und j xr=xr-box*nint(xr/box) periodische Randbedingungen r2=xr**2 if (r2.lt.rc2) then überprüfe, ob Distanz < cut-off r2i=1/r2 r6i=r2i**3 ff=48*r2i*r6i*(r6i-0.5) in diesem Beispiel:Lennard-Jones Potential f(i)=f(i)+ff*xr Update der Kräfte f(j)=f(j)-ff*xr en=en+4*r6i*(r6i-1)-ecut Update der Energie endif enddo enddo return end nach [Frenkel, Smit 1996] Computational Chemistry

  20. Nachbarlisten bringen Speedup Es gibt N*(N-1) Paar-Wechselwirkungen. Nicht alle müssen bei jedem Iterationsschritt berechnet werden! Spare Zeit mit Nachbarliste, die Buch führt, welche Teilchen in benachbarten Zellen liegen (links) bzw. in äußerer Kugelschale (rechts). Zell-Liste Verlet-Liste [Frenkel,Smit 1996] Bei Schritt i+1 braucht nur überprüft zu werden, ob Teilchen aus diesen Zellen nun innerhalb des cut-off liegen. Skaliert mit N. Alle n Schritte ist Neuberechnung der Nachbarliste nötig. Skaliert mit N2. Computational Chemistry

  21. Satz von Taylor Eine Funktion f sei in einem Intervall [x0-, x0+] (n+1)-mal differenzierbar. Dann gilt in diesem Intervall die Taylorentwicklung: Computational Chemistry

  22. Integration der Bewegungsgleichungen Taylor-Entwicklung der Teilchenkoordinate r(t) Dies ist der Verlet-Algorithmus. Die Geschwindigkeiten kommen nicht vor, können aber aus der Trajektorie berechnet werden: Computational Chemistry

  23. Leap-frog Algorithmus Leap-frog (“Bocksprung-“) Algorithmus Für den update der Geschwindigkeiten gilt: Computational Chemistry

  24. Integriere die Bewegungsgleichungen subroutine integrate(f,en) integriere Bewegungsgleichungen sumv=0 sumv2=0 do i=1,npart MD Schleife xx=2*x(i)-xm(i)+delt**2*f(i) Verlet-Algorithmus v(i)=(xx-xm(i))/(2*delt) Geschwindigkeit sumv=sumv+v(i) Geschwindigkeit des CMS sumv2=sumv2*v(i)**2 totale kinetische Energie xx(i)=xx update Positionen enddo temp=sumv2/(3*npart) aktuelle Temperatur etot=(en+sumv2)/(2*npart) Gesamtenergie pro Teilchen return end nach [Frenkel, Smit 1996] Computational Chemistry

  25. Algorithmen für MD Leap-frog und Verlet-Algorithmen sind am gebräuchlichsten in typischen MD-Simulationspaketen. Sie benötigen lediglich die ersten Ableitung der Energie nach den Teilchenkoordinaten, also die Kräfte auf die Atome. In Simulationen im NVE-Ensemble (also bei konstanter Teilchenzahl N, konstantem Volumen V und konstanter Energie E) bleibt die Energie bei Verwendung von Simulationszeitschritten t  0.5 fs erhalten. Bei Einfrieren der Bindungslängen (“SHAKE-Algorithmus”) kann man sogar bis t  2.0 fs verwenden. Computational Chemistry

  26. Algorithmen für MD Es existieren kompliziertere Algorithmen, die höhere Ableitungen der Energie verwenden. Diese Algorithmen erlauben es, längere Zeitschritte bis ca. t  5.0 fs zu verwenden. Allerdings ist der Aufwand, die höheren Ableitungen zu berechnen, gewöhnlich grösser als die Einsparung durch die kleinere Anzahl an Schritten. Computational Chemistry

  27. Moleküldynamik MD-Simulationen haben mehrere Aufgaben * Generierung von Konformationen, Suche im Konformationsraum * Dynamische Simulation von Systemen im Gleichgewicht * Nicht-Gleichgewichts-Simulationen Computational Chemistry

  28. MD-Simulation für Protein-Assoziation Computational Chemistry

  29. Langreichweitige elektrostatische Wechselwirkung • Empirische Kraftfelder beschreiben die nichtgebundene Wechselwirkung • meistens durch 2 Terme: • (1) ein Lennard-Jones-Potential C12/r12 - C6/r6 • Der Lennard-Jones Term fällt sehr schnell mit zunehmender Entfernung ab und • kann guten Gewissens für r 1 nm vernachlässigt werden. • die elektrostatische Coulombsche Wechselwirkung aller Teilchen • miteinander mittels ihrer Punktladungen qi • Die elektrostatische Wechselwirkung fällt nur sehr langsam mit r-1ab. • Die Verwendung einer Abschneidefunktion (cut-off) führt daher oft zu Artefakten. Computational Chemistry

  30. Cut-off Artefakte • Aufheizung des Systems durch Austritt dipolarer Gruppen mit Vorzugsrichtung und Wiedereintritt mit beliebiger Orientierung • strukturelle + dynamische Abweichung von experimentellen Daten! • Zu große “Ordnung” des Systems, da jede cut-off-Sphäre quasi von Vakuum umgeben ist. • Probleme bei Equilibrierung von kritischen Systemen, die z.B. viele polare Gruppen enthalten, deren Ladungen in der Lösung durch Gegenionen abgesättigt werden, wie DNA-Lösungen oder Biomembranen. • In MD-Simulationen von DNA mit cut-off wird z.B. während ns Simulation keine stabile Konformation gefunden; große Abweichung von Kristallstrukturen. Computational Chemistry

  31. Ewald-Summe Bei dieser Methode stellt man sich das System unendlich weit in 3 Dimensionen fortgesetzt, also periodisch, vor. Die Summationsmethode der gesamten Wechselwirkung in Kristallen wurde in den ‘20er Jahren von Ewald entwickelt. Der entscheidende Schritt ist, den kritischen Term r-1 in zwei rasch konvergierende Terme zu zerlegen. Der linke Termin konvergiert schnell im reellen Raum, der rechte schnell im reziproken Raum. Computational Chemistry

  32. NVE Bedingung Bislang wurden die Zustände des klassischen Systems unter den Bedingungen N (Teilchenzahl) = konstant V (Volumen) = konstant E (Gesamtenergie) = konstant simuliert. Wenn man annimmt, daß zeitliche Mittelwerte = Ensemblemittelwerte sind, dann entsprechen Mittelwerte über eine konventionelle MD-Simulationen den Ensemble-Mittelwerten eines mikrokanonischen Ensembles. Jedoch ist es oft wünschenswert, Simulationen in anderen Ensembles durchzuführen, z.B. NVT oder NPT, in denen die Temperatur (T) und/oder der Druck (p) konstant gehalten werden. Computational Chemistry

  33. Bedingung konstanter Temperatur: Kopplung an Wärmebad Wenn man ein System mit einem großen Wärmereservoir koppelt, bekommt das System eine Temperatur. D.h. die Anzahl der Zustände bei Energie E gehorcht Boltzmann-Verteilung, und die Verteilung der Geschwindigkeiten gehorcht für ein klassisches System der Maxwell-Boltzmann-Verteilung: Computational Chemistry

  34. Grundsätze In einem klassischen Vielteilchensystem ist die Energie über all diejenigen Freiheitsgrade Nfgleichverteilt, die quadratisch in den Hamiltonian (Energiefunktion) des Systems eingehen, also z.B. die Geschwindigkeiten. Für die mittlere kinetische Energie pro Freiheitsgrad gilt: Die Klammer ... bedeutet eine zeitliche Mittelung über die Konfigurationen des Systems z.B. während einer MD-Simulation. In einer MD-Simulation kann man diese Beziehung bequem zur Berechnung der Temperatur des Systems heranziehen. Computational Chemistry

  35. Andersen-Thermostat Zu bestimmten Intervallen wird die Geschwindigkeit eines zufällig ausgewählten Teilchens neu aus einer Maxwell-Boltzmann-Verteilung der Geschwindigkeiten zugeteilt, entsprechend einer stochastischen Kollision dieses Teilchens mit einem Solvensmolekül eines imaginären Bads. Problem: die stochastischen Kollisionen stören die dynamischen Eigenschaften in unrealistischer Weise. Die Geschwindigkeits-Autokorrelationsfunktion fällt zu schnell ab. Dieser Thermostat ist nicht geeignet um dynamische Eigenschaften wie Diffusionskoeffizienten zu simulieren. Alle statischen Eigenschaften werden aber korrekt berechnet. Andere Verfahren (Nosé-Hoover) erlauben es, sowohl thermodynamische wie dynamische Eigenschaften methodisch sauber zu berechnen. Computational Chemistry

  36. Protokoll für Moleküldynamik-Simulationen generiere Koordinaten minimiere Struktur (EM) ordne Anfangsgeschwindigkeiten zu Aufwärmphase (MD) Equilibrierung (MD) Produktionslauf (MD) Analyse Skaliere Geschwindigkeiten Temperatur OK? Nein Ja Computational Chemistry

  37. Analyse von Trajektorien Mittlere Energie RMS-Unterschied zwischen 2 Strukturen RMS-Fluktuationen (diese hängen mit den kristallographischen B-Faktoren zusammen) Computational Chemistry

  38. Zusammenfassung • MD-Simulation ist eine etablierte und mittlerweile sehr robuste Simulationstechnik. • Ständig neue Anwendungsgebiete • durch Beschleunigung der Algorithmen, • durch Zunahme der Prozessorgeschwindigkeit, • durch Parallelisierung der Programme. • Ausgangspunkt ist gewöhnlich Röntgenkristallstruktur eines Biomoleküls • oder eine modellierte Struktur. • Wichtig sind • sorgfältige Konstruktion der Startkonfiguration • sinnvolle Wahl von Kraftfeld und Simulationsbedingungen • sorgfältige Equilibrierung des Systems • analysiere nur Eigenschaften des Systems für die Sampling OK ist. Computational Chemistry

  39. Zusätzliche Folien Computational Chemistry

  40. Ewald-Summe Die Fehlerfunktionerf(x) ist definiert Für die entsprechende als das Integral über eine Gaußverteilung Fehlerfunktion erfc(x) gilt: mit dem Mittelwert 0 und der Varianz ½. Mit erfc(x) + erf(x) =1 zerlegt man nun 1/r in zwei Terme: Für r fällt erfc(r) schnell ab. Der erste Term konvergiert daher im reellen Raum schnell. erf(r)/r ist eine stetige Funktion. Ihre Fouriertransformierte fällt daher schnell ab. Computational Chemistry

  41. Ewald (2) • Der Coulomb-Term wird damit in drei Terme zerlegt: • Eine Summe im reellen Raum, die die Coulomb-Wechselwirkungen zwischen benachbarten Atomen innerhalb eines cut-off Radius beschreibt • Eine Summe im reziproken Raum, die die Wechselwirkung mit dem restlichen System beschreibt: • Hierbei ist V das Volumen der Simulationsbox, m = (l, j, k) ein Vektor im reziproken Raum, und Computational Chemistry

  42. Ewald (2) (c) Die Selbst-Energie der artifiziell hinzugefügten Gegenladungen: Mittels des Wertes von  kann man den Anteil von direkter und reziproker Summe verschieben. In der Praxis benutzt man einen cut-off von ca. 1nm für den reellen Raum. Computational Chemistry

  43. Ewald (3) Eine physikalische Interpretation dieses Vorgehens ist oben dargestellt: Die Punktladungen werden in der direkten Summe ergänzt durch Gaußfunktionen mit umgekehrtem Vorzeichen, so dass sich die integrierte Ladungsdichte aufhebt. Die Wechselwirkung dieser hinzugefügten Gaußfunktionen wird dann im reziproken Raum berechnet. Computational Chemistry

  44. Langreichweite Elektrostatik • Es gibt weitere Techniken • Particle-Mesh-Ewald, eine vor allem für größere Systeme sehr effiziente Variante der ursprünglichen Ewald-Summe • Fast-Multipole-Methode • - Reaktionsfeld-Methoden, nur für isotrope Systeme wie Flüssigkeiten... • die alle in bestimmten Bereichen Vorteile bieten. • Es stehen damit eine Reihe von zuverlässigen Methoden zur Berücksichtigung langreichweitiger elektrostatischer Effekte in Simulationen zur Verfügung. • Wichtig: • Durch Verwendung einer künstlichen Periodizität werdendie strukturellen, dynamischen und vermutlich thermodynamischen Eigenschaften der Systeme verändert! • Dies muß jedoch nicht von Nachteil sein, sofern man die relative Größe dieser Effekte kennt und evtl. bei der (Neu-) Parametrisierung von Kraftfeldern berücksichtigt. Computational Chemistry

  45. Elektrostatisches Potential an Grenzflächen Wasser: Flüssigkeit - Dampf Wasser - Lipid Computational Chemistry

More Related