Kraftfeldmethoden (Teil 1)

Marcus Weber

Vorlesung: Algorithmische Bioinformatik

WS 05/06

1. Kraftfelder

1.1 Kraftfelder (mathematisch)

1.2 Wie real sind Kraftfelder?

TEIL 1

2. Minimierung

2.1 lokale Minimierunng

2.2 globale Minimierung

3. Moleküldynamik

3.1 Hamilton-Dynamik

3.2 Numerische Lösung

TEIL 2

4. Statistische Thermodynamik

4.1 Ensemble-Begriff

4.2 Monte-Carlo Methoden

5. Metastabile Konformationen

1. Kraftfelder

1.1 Die analytische Form von Kraftfeldern

1.2 Wie real sind Kraftfelder?

1.1 Kraftfelder (mathematisch)

x

z

y

Atompositionen (Orte)

pro Atom (x,y,z) Position

3 x Anzahl Atome

viele Positionen

Wechselwirkungspotential (Kraftfeld)

Abstoßungskräfte

Bindungsschwingung

elektrostatische

Kräfte

Rotation um

Bindungen

z.B.

Verbreitete Programmpakete und Kraftfelder

AMBER; Assisted Model Building and Energy Refinement)

CHARMM; Chemistry at HARvard Macromolecular Mechanics)

GROMOS (GROningen MOlecular Simulation package)

MMFF94, the Merck Molecular Force Field

Tripos, the force field of the Sybyl molecular modeling program

http://msdlocal.ebi.ac.uk/docs/mmrefs.html

Unterschiede: Modellierung von Wasserstoffbrücken- bindungen, des Lösungsmittels, physikalische Terme, Anwendungsgebiete (Proteine, Nukleinsäuren, etc.)

Die Entwicklung brauchbarer Kraftfelder ist „eine Wissenschaft für sich“!

… um nur einige Programmpakete/Kraftfelder zu nennen.

P.J. Steinbach „Introduction to Macromolecular Simulation“

source: http://cmm.info.nih.gov/intro_simulation/node15.html

… und so sieht’s in Formeln gegossen aus:

Zur Kraftfeldberechnung werden die kartesischen Atom-

Koordinanten q in internen Freiheitsgrade/Variablen

(Bindungslängen l, Winkel , Abstände rij.) umgerechnet.

(n bezeichnet Anzahl der Minima)

1.2 Wie real sind Kraftfelder?

Born-Oppenheimer Näherung: Ort q und Impulse p der Atomkerne

sind zur Berechnung der Gesamtenergie gleichzeit scharf zu bestimmen.

Das widerspricht der Heisenberg'schen Unschärferelation.

H(q,p) = V(q) + K(p)

Gesamtenergie potentielle Energie kinetische Energie

„Die Elektronen richten sich nach den Atomkernen“

Beste Approximation für Elektronenanordnung z und Energie

V(q) in dieser Näherung ist die Lösung der Schrödinger-Gleichung abhängig von den gegebenen Kernpositionen q.

L(q) z = V(q) z

Operator Zustandsfunktion ergibt Elektronenverteilung

diskrete Energiewerte

Die Lösung der Schrödinger-Gleichung (als Eigenwertgleichung)

ergibt ein Spektrum von möglichen Energiewerten V(q) und

dazu gehörenden Elektronendichten, sprich Bindungsmodi.

Im Kraftfeld jedoch gibt es zu jeder Kernanordnung q genau einen

Energiewert V(q), eine elektronische Anregung ist nicht möglich.

Grundzustands-Näherung

Im Kraftfeld resultiert der Bindungsmodus und die Ladungsdichte

nicht aus den Kernpositionen, sondern wird a priori „geraten“.

Deshalb gilt das Kraftfeld auch nur, solange der Bindungsmodus und

die Partialladungen für gegebenes q richtig geraten wurden.

Lokalitäts-Näherung

In klassischen Kraftfeldern gibt es im allgemeinen keine

Bindungsbrüche, Ladungsverschiebungen, Komplex-Bildung...

Verbreitete quantenchemische Programmpakete

Gaussian 03; sehr populäres ab-initio Programmpaket, Basisfunktionen

sind Gauss-Glockenförmig

Jaguar; ein auf Geschwindigkeit optimiertes ab-initio Programmpaket

PSI 3; ein auf Genauigkeit optimiertes ab-initio Paket

MOPAC; ein semi-empirisches QC-Programmpaket, das gut zur Analyse chemischer Reaktionen geeignet ist

AMSOL, Quantenchemie in Lösung (semi-empirisch)

http://www.asc.edu/supercomputing/software.shtml

Unterschiede: Es gibt Programme, die ganz auf der Basis der

Schrödinger-Gleichung arbeiten (ab-initio) und solche, die schon Expertenwissen über Atomtypen benutzen (semi-empirisch)

Quantenchemische Berechnungen sind oft nur für kleine

Molekülsysteme möglich, dann liefern sie jedoch die

genauesten Ergebnisse zu Ladungsdichte und Bindungsmodi!

2. Minimierung

2.1 lokale Minimierung

2.2 globale Minimierung

2.1 Lokale Minimierung

1-dimensional

mehrdimensional

Hinreichende Bedingung für ein lokales Minimum bei

x*=[x1*,...,xn*]:

Gradient und Hesse-Matrix

1. Der Gradient gradV(x*) ist der Null-Vektor.

Mehrdimensionales Äquivalent zu V'(x*)=0

2. Die Hesse-Matrix hessV(x*) ist positiv-definit

(für alle nichtnull Vektoren v ist (vT hessV(x*) v) >0).

Mehrdimensionales Äquivalent zu V''(x*)>0

Berechnung des Gradienten-Vektor: Die Komponente i des

Gradientenvektors erhält man, indem man die Funktion V nach xi ableitet,

wobei alle anderen Komponenten als konstant angesehen werden.

Berechnung der Hesse-Matrix: Das Matrix-Element (i,j) der Hesse-Matrix

erhält man, indem man die Komponente i des Gradienten nach xj ableitet.

Taylorentwicklung 2.Ordnung

Die meisten lokalen Optimierungsverfahren basieren auf einer

Approximation der Potentialfunktion mittels einer Taylor-Entwicklung

(lokale Näherung des Potentials als quadratische Funktion):

V(x) = V(x*) + (x-x*)T gradV(x*) + 0.5 (x-x*)T hessV(x*) (x-x*) + r(x)

Je nachdem, wieviele Terme dieser Entwicklung in der

Minimierungsroutine verwendet werden, spricht man von Verfahren

nullter, erster oder zweiter Ordnung.

Grundalgorithmus lokaler Verfahren

Gegeben ist ein Startpunkt x0 und der Schrittzähler k=0.

1. Prüfe, ob ein Abbruchkriterium erfüllt ist, z.B. ist die Norm

des Vektors gradV(xk) kleiner als ein vorgegebener Schwellwert?

Wenn Abbruchbedingung erfüllt, dann STOP.

2. Berechne als „Suchrichtung“ einen n-dimensionalen Vektor

dk , z.B. aus einer Taylor-Approximation des Potentials.

3. Bestimme von xk aus in Suchrichtung einen Punkt mit minimalem

Funktionswert, d.h. bestimme einen optimalen Wert s>0, so dass V

in xk+1= xk + s dk minimal ist. (1-dimensionales „line-search“)

4. k=k+1, gehe zu 1.

Suchrichtung (1.Ordnung)

V(xk+1) = V(xk) + (xk+1-xk)T gradV(xk) + r(xk+1)

In dieser Formel stellt (xk+1-xk) die Suchrichtung dar. Man kann ablesen,

dass kleinster Wert V(xk+1) zu erwarten ist, wenn als Suchrichtung der

Vektor dk=-gradV(xk) gewählt wird (steepest descend).

Erinnerung: Skalarprodukt

Suchrichtung (2.Ordnung)

V(xk+1) = V(xk) + (dk)T gradV(xk) + 0.5 (dk)T hessV(xk) dk + r(xk+1)

Bestimme Minimum:

Obige quadratische Funktion in dk hat folgende Ableitung:

damit ist die Lösung dk des linearen Gleichungssystems

hessV(xk) dk = - gradV(xk) eine Nullstelle dieser Ableitung und

auch ein lokales Minimum der obigen Funktion, solange hessV(xk)

positiv definit ist (trust region, Konvergenzradius, Einzugsgebiet des

Minimums). Das so gewonnene Verfahren heißt Newton-Verfahren.

gradV(xk) + hessV(xk) dk

Quasi-Newton Verfahren

Da die Berechnung der Hessematrix von Kraftfeldern im Allgemeinen

als Programmcode nicht vorliegt, verwenden die sogenannten

Quasi-Newton Verfahren eine Approximation der Hesse-Matrix

(DFP-Algorithmus) oder eine Approximation der inversen

Hesse-Matrix (BFGS-Algorithmus).

Diese Approximationen werden in jedem Schritt des Grundalgorithmus

verbessert und dabei immer positiv-definit gewählt.

Das BFGS-Verfahren wird sehr häufig in der Optimierung von

Molekülpotentialfunktionen verwendet. Bei der Approximation der

inversen Hesse-Matrix entfällt das Lösen eines linearen Gleichungssystems

zur Bestimmung der Suchrichtung.

2.2 globale Minimierung

Multi-Start-Ansätze: Lokale Optimierung von verschiedenen Startpunkten aus.

Branch-and-Bound: Wird jetzt behandelt...

Statistische Verfahren: Zum Beispiel „simulated annealing“ (globales Minimum=Zustand eines Moleküles, der bei 0 Kelvin angenommen wird), siehe Abschnitt 4.

Andere Verfahren: Homotopie-Verfahren, CGU-Verfahren, Tabu search, genetische Algorithmen etc. (werden nicht besprochen)

Grundalgorithmus (branch & bound)

Globale Variablen: Bisher bester Funktionswert V* und die zugehörigen

Koordinaten x* lassen sich mit Ergebnis einer lokalen Minimierung

initialisieren. Branch & Bound ist eine rekursive Methode. Eingabe ist

ein (mehrdimensionales) Koordinatenintervall [xu, xo] (=Box).

Prozedur bbStep(xu,xo)

BEGIN

Berechne V(x) für einen beliebigen Punkt x aus [xu,xo].

Ist V(x) besser als V*, dann ersetze V* durch V(x) und x* durch x.

Berechne eine untere Schranke u für V auf dem Bereich [xu,xo].

Ist u größer als V*, dann lohnt eine weitere Untersuchung des Intervals [xu,xo]

nicht, da das Optimum hier nicht zu finden ist: Verlasse Prozedur!

Ist das Intervall [xu,xo] kleiner als eine vorgegebene Mindestgroesse, dann

lohnt eine weitere Unterteilung nicht: Verlasse Prozedur!

Zerlege [xu,xo] (entlang der längsten Ausdehnung) in zwei Unterintervalle [xu1,xo1]

und [xu2,xo2]. Rufe bbStep(xu1,xo1) und danach bbStep(xu2,xo2) auf.

END.

Problemstelle im Grundalgorithmus ist die Bestimmung der unteren

Schranke für V auf einem Intervall [xu,xo].

Ist die untere Schranke zu schlecht bestimmt, dann untersucht der

Algorithmus im schlimmsten Fall das gesamte Ausgangsintervall mit

einer Auflösung gemäß der vorgegebenen Intervallmindestgröße ab.

Ist die untere Schranke falsch bestimmt, dann wird wohlmöglich ein

Intervall verworfen, obwohl es das globale Optimum enthält.

WICHTIG: Zur Bestimmung der unteren Schranke auf einem Intervall

benötigt man auf jeden Fall eine „globale Information“ über das Potential

V, z.B. so etwas wie „Was ist der größtmögliche Gradient von V?“ oder

„Was ist die maximale Krümmung von V?“. Oft können die Antworten

auf diese Fragen nur geraten werden.

a- Branch and Bound

Androulakis, Maranas and Floudas (1995)

Problem der globalen Optimierung: hessV(x) ist nicht überall positiv

definit, lokale Methoden wie das Newton-Verfahren konvergieren nicht

für jeden Startwert. Wäre hessV(x) überall positiv definit, dann gäbe es auch

nur ein eindeutiges lokales Minimum (das damit auch global minimal wäre).

Idee: Zwinge hessV auf positiv definit. Finde zunächst eine (negative)

Zahl a, so dass

Auf dem Intervall [xu,xo] konstruiere eine Hilfsfunktion F mittels:

Eigenschaften dieser Funktion:

F(x) ist auf dem Intervall [xu,xo] stets kleiner als V(x). Das Minimum
von F auf dem Intervall [xu,xo] ist also auch eine untere Schranke
für V, die im B&B-Algorithmus verwendet werden kann.

Die Hesse-Matrix von F lautet hessF(x)=hessV(x)-a E, wobei E die
Einheitsmatrix ist. Diese ist per Konstruktion von a stets positiv
definit. Damit kann man das Newton-Verfahren auf F(x) anwenden,
um durch restringierte Optimierung auf [xu,xo] eine untere Schranke
von V zu berechnen.

Marcus Weber

2006-01-17T08:36:28