Document [original]

Atomistische Modellierung der Struktur

und Funktion von Biomolek¨

ulen

Marcus Elstner

Habilitationsschrift zum Habilitationsgesuch

Fachbereich Physik

der

Universit¨

at Paderborn

Mai 2003

Inhaltsverzeichnis

1 Einf¨

uhrung 5

1.1 Methoden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.2 Problemstellung am Beispiel Bakteriorhodopsin . . . . . . . . . . 13

2 DFTB 19

2.1 Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.2 Einbeziehung weiterer Atomtypen . . . . . . . . . . . . . . . . . . 22

2.3 Erweiterung fuer H-Bruecken . . . . . . . . . . . . . . . . . . . . 24

2.4 Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.5 Struktur von Peptiden . . . . . . . . . . . . . . . . . . . . . . . . 26

2.6 Struktur des Retinals in Bakteriorhodopsin . . . . . . . . . . . . . 29

3 Grosse Systeme 31

3.1 Kombination von Methoden:QM/MM . . . . . . . . . . . . . . . . 31

3.2 O(N) Methoden . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.3 Kombinierte QM/MM O(N) . . . . . . . . . . . . . . . . . . . . . 35

4 Stabilit¨

at von Protein Sekund¨

arstrukturen 37

5 Optische Eigenschaften 41

6 Zeitabh¨

angige Ph¨

anomene 43

7 Zusammenfassung und Ausblick 45

A Relevante Ver¨

offentlichungen 51

4INHALTSVERZEICHNIS

Kapitel 1

Einf¨

uhrung

In den letzten Jahren sind theoretische Studien realer Materialien ein wichtiger

Teil der Materialforschung geworden, sei es in der Physik, Chemie oder Biologie.

Dies liegt einerseits an den grossen Fortschritten in der Computertechnik und

andererseits an der Enwicklung effizienter numerischer Simulationsmethoden.

Untersucht werden neben kristallinen und amorphen Festk¨

orpermodifikatio-

nen auch Cluster, Fullerene, organische Molek¨

ule, Polymere und biologische Mo-

lek¨

ule. Ein Hauptschwerpunkt liegt dabei in dem atomistischen Verst¨

andnis der

Struktur und Dynamik dieser Materialien, aber auch weitergehend in der Auf-

kl¨

arung des Verh¨

altnisses von Struktur, Eigenschaft und Funktion. Ferner kann

durch die theoretische Vorhersage experimentell relevanter Daten Hilfestellung

bei der Strukturaufkl¨

arung und bei der Interpretation experimenteller Ergebnisse

gegeben werden. ¨

Uber Beitr¨

age zum Verst¨

andnis von Materialien hinaus k¨

onnen

aber auch gezielte Vorschl¨

age zur Modifikation, zum Design oder Optimierung

gemacht werden. Ein prominentes Beispiel ist hier das sogenannte ’computer ai-

ded drug design’, wo durch theoretische Simulation Struktur und Eigenschaften

neuer Medikamente gezielt untersucht werden k¨

onnen, bevor sie synthetisiert und

getestet werden. Fener k¨

onnen Informationen auf theoretischem Wege gewonnen

werden, die dem Experiment nicht, oder nicht ohne weiteres zug¨

anglich sind.

Dies betrifft etwa die Modellierung chemischer Reaktionen, die Simulation von

Prozessen der Anlagerung von Molek¨

ulen, oder auch Ladungs- und Energietrans-

portprozesse, die auf solch kurzen Zeitskalen stattfinden, dass sie vom Experiment

z.Z. nicht aufgel¨

ost werden k¨

onnen.

Die Simulation biomolekularer Prozesse, stellt eine besondere Herausforderung

f¨

ur theoretische Methoden dar. Dies liegt an der Zusammensetzung, Ausdehnung

und spezifischen Art der Wechselwirkung dieser Systeme mit ihrer Umgebung.

Ausgehend von Methoden, die in der Cluster- und Festk¨

orperphysik inzwischen

eine breite Verwendung finden, liegt der Schwerpunkt dieser Arbeit deshalb in der

Entwicklung neuer, und in der Weiterenwicklung und Kombination bestehender

6KAPITEL 1. EINF ¨

UHRUNG

Methoden. Ziel ist es, ein Methodenspektrum zu entwickeln, das den Herausfor-

derungen bei der Modellierung biochemischer Prozesse gewachsen ist.

Dazu wird im Folgenden zun¨

achst ein ¨

Uberblick ¨

uber die verschiedenen Si-

mulationsmethoden gegeben, die in der theoretischen Materialforschung einge-

setzt werden (Kapitel 1.1). Daran anschliessend wird skizziert, welche besonde-

ren Anforderungen biologische Strukturen und Prozesse an diese Methoden stellen

(Kapitel 1.2). Der Hauptteil dieser Arbeit besch¨

aftigt sich mit der Entwicklung

und dem Testen von Methoden, die eine Beschreibung komplexer biologischer

Strukturen erlauben. Ausgangspunkt der Entwicklung ist ein quantenmechnisches

Rechenverfahren, das bisher auf organische Molek¨

ule, Cluster und Festk¨

orper-

strukturen sehr erfolgreich angewendet wurde. Dieses wird in Kapitel 2 zun¨

achst

kurz eingef¨

uhrt, und erste Erweiterungen des Schemas und Anwendungen und

Tests der Methode f¨

ur biologische Modellstrukturen werden vorgestellt. Kapitel

3 besch¨

aftigt sich mit der Beschreibung grosser Systeme. Hier werden verschie-

dene Algorithmen vorgestellt, die eine quantenmechanische Behandlung grosser

biologischer Molek¨

ule erlauben. In Kapitel 4 wird eine Anwendung vorgestellt,

die Simulation der Stabilit¨

at und Dynamik von Polypeptiden und eines Proteins,

welche erst durch die vorgestellten methodischen Neuerungen m¨

oglich geworden

ist. Kapitel 5 widmet sich der Beschreibung von optischen Anregungen innerhalb

der hier verwendeten Methode, und Kapitel 6 generalisiert dieses Vorgehen auf die

Beschreibung von zeitabh¨

angigen Prozessen, die durch eine Wechselwirkung der

Molek¨

ule mit Licht ausgel¨

ost werden k¨

onnen. Kapitel 7 fasst den gegenw¨

artigen

Stand der Methode zusammen und gibt einen Ausblick auf weitere Entwicklungs-

richtungen. Die Originalarbeiten zu den Themen in den jeweiligen Kapiteln sind

in dem Anhang beigef¨

ugt.

1.1 Methoden

Zur Bestimmung der Eigenschaften und der Dynamik molekularer Systeme muss

im allgemeinen Fall das zeitabh¨

angige quantenmechanische Vielteilchenproblem

f¨

ur die elektronischen und nuklearen Freiheitsgrade gel¨

ost werden. Dieses Vor-

gehen ist jedoch sehr aufwendig. Um zu praktisch realisierbaren Algorithmen zu

kommen wird zun¨

achst der Spezialfall des zeitunabh¨

angigen Problems betrachtet.

Ferner wird angenommen, daß die Dynamik der elektronischen Freiheitsgrade auf

wesentlich kl¨

urzeren Zeitskalen stattfindet, als die Dynamik der Kernfreiheits-

grade. Im Rahmen der sogenannten Born-Oppenheimer (BO) N¨

aherung wird die

Schr¨

odingergleichung (SG) separiert in einen Teil, der die elektronischen Frei-

heitsgrade bei festgehaltenen Kernfreiheitsgraden beschreibt und einen Teil, der

die Kernwellenfunktion in einem Potential beschreibt, das aus der L¨

osung der

elektronischen SG resultiert. Die elektronische Wellenfunktion h¨

angt dann nur

1.1. METHODEN 7

noch parametrisch von den Kernkoordinaten ab (und nicht noch zus¨

atzlich etwa

von den Kernimpulsen), es entsteht das Bild von Potentialenergiefl¨

achen (PES:

’potential energy surface’), auf denen sich die Kerne bewegen. Der Fehler durch

Einf¨

uhren der BO Approximation ist i.A. klein gegen¨

uber den nachfolgend disku-

tierten weiteren N¨

aherungen. Es gibt jedoch chemische Situationen, in denen sie

’zusammenbricht’, etwa bei hochenergetischen St¨

ossen von Teilchen, Wechselwir-

kung mit starken elektromagentischen Feldern, oder in bestimmten molekularen

Konformationen, in denen sich zwei PES entlang einer Reaktionskoordinate (etwa

elektronischer Grund- und angeregter Zustand) kreuzen. In diesem Fall m¨

ussen al-

ternative Zug¨

ange gew¨

ahlt werden, wie in Kapitel 6 ausgef¨

uhrt. In einigen F¨

allen

ist auch die klassische Beschreibung der Atomkerne nicht angemessen. Die wird

akut, wenn chemische Reaktionen untersucht werden, die den Transfer von Pro-

tonen beinhalten.

In einem weiteren N¨

aherungsschritt werden die Atomkerne als klassische Punkt-

teilchen behandelt. Die Quanteneffekte der Atomkerne k¨

onnen jedoch mit Kennt-

nis der PES entlang einer Reaktionskoordinate im Nachinein ber¨

ucksichtig wer-

den. Hierf¨

ur stehen in der Literatur eine Reihe von Zug¨

angen zur Verf¨

ugung, auf

die aber in dieser Arbeit nicht weiter eingegangen wird.

Viele Eigenschaften molekularer Systeme sind demnach durch die Kenntnis

der elektronischen Gesamtenergie des Systems zug¨

anglich. Chemische Informatio-

nen k¨

onnen aus der Kenntnis der Energie in Abh¨

angigkeit von den Kernkoordi-

naten gewonnen werden. Dies betrifft etwa molekulare Strukturen, wie aber auch

Energiebarrieren f¨

ur chemische Reaktionen etc. Weitere Eigenschaften k¨

onnen

aus der Kenntnis der Energie¨

anderungen in Abh¨

angigkeit von ¨

ausseren St¨

orun-

gen gewonnen werden. Auf diese Weise k¨

onnen spektroskopische Daten, wie etwa

optische Eigenschaften, NMR/ESR Spektren, Schwingungsspektren (IR, Raman)

berechnet werden.

Zur Bestimmung der elektronischen Gesamtenergie ist aber die exakte L¨

osung

der Schr¨

odingergleichung (SG) f¨

ur Mehrelektronensysteme wegen des hohen Re-

chenaufwandes i.A. nicht praktikabel. Die L¨

osung der SG etwa durch Monte Carlo

Integration [1] ist auf die Beschreibung von Systemen mit wenigen Atomen be-

schr¨

ankt. Praktikable numerischen Rechenverfahren f¨

ur gr¨

ossere Systeme basieren

daher auf weitergehenden geeigneten N¨

aherungen der elektronischen SG.

Die Methoden, die ausgehend von der molekularen Struktur die elektronische

Gesamtenergie komplexer Systeme zu bestimmen erlauben, lassen sich in drei

Kategorien einteilen: Die wellenfunktionsbasierten ab initio (1.1.1) und Dichte-

funktional Methoden (1.1.2) , die semi-empirischen Verfahren (1.1.3) und die

klassischen Kraftfeldmethoden (1.1.4).

8KAPITEL 1. EINF ¨

UHRUNG

Wellenfunktionsbasierte ab initio Methoden

Einen quantenmechnischen Zugang zu grossen Systemen erlauben erst Wellen-

funktionsbasierte Methoden, die eine n¨

aherungsweise L¨

osung der SG durch einen

Ansatz f¨

ur die Vielteilchen-Wellenfunktion erreichen. Startpunkt ist die Hartree-

Fock (HF) Methode, die als Ansatz ein antisymmetrisiertes Produkt aus Ein-

teilchenwellenfunktionen zu Grunde legt. Davon ausgehend kann die L¨

osung der

HF Methode durch Einbeziehung der Elektronenkorrelationen verbessert werden

(siehe etwa [2]). Zum Einen durch st¨

orungstheoretische Methoden (benannt nach

M¨

oller und Plesset, die als erste die Raileigh-Schr¨

odingersche St¨

orungstheorie in

diesem Kontext anwendeten: MP2, MP4), zum Anderen durch Methoden, die aus-

gehend von der HF Wellenfunktion erweiterte Wellenfunktionsans¨

atze verwenden

(CI, CASSCF, ...).

Die HF Methode ist heute ein Standardverfahren, das insbesondere zur Be-

schreibung von Struktur und Schwingungseigenschaften von Molek¨

ulen sehr er-

folgreich eingesetzt wird. Mit ihr k¨

onnen Systeme von bis zu etwa 100 Atomen

routinem¨

assig behandelt werden. Die Energetik der chemischen und die Struk-

tur und Energetik der H-Br¨

ucken und VdW Bindungen k¨

onnen mit HF nicht,

oder nicht mit ausreichender Genauigkeit berechnet werden. Hierzu ist die Ein-

beziehung der Elektronenkorrelation unbedingt erforderlich. Die entsprechenden

Methoden sind aber wesentlich aufwendiger, mit MP2 etwa k¨

onnen nur noch

Systeme mit bis zu etwa 30 Atomen betrachtet werden.

Die Dichtefunktional Theorie (DFT)

Als Alternative zu den ab initio Verfahren hat sich in den letzten Jahren die

sogenannte Dichtefunktionaltheorie (DFT) etabliert (siehe etwa [3, 4]). Ihre At-

traktivit¨

at beruht darin, daß sie viele relvante molekulare Systeme mit der Genau-

igkeit von post-HF Methoden (etwa MP2) beschreibt, bei gleichzeitig wesentlich

geringerer Rechenzeit.

In diesem Abschnitt sollen einige grundlegende Gleichungen der Dichtefunk-

tionaltheorie eingef¨

uhrt werden ([3, 4]), die als Startpunkt f¨

ur die Entwicklung

eines approximativen Verfahrens zur Bestimmung der elektronischen Struktur

molekularer Systeme dienen.

Der (nicht relativistische) Hamiltonoperator f¨

ur ein System von N Elektronen

in dem Coulomb-Potential der Atomkerne ist (in atomaren Einheiten) durch

i

−∆i

2+VKern(ri) +

j>i

|rj−ri|

(1.1)

gegeben. VKern(ri) ist das auf die Elektronen wirkende Potential der Atomr¨

umpfe.

1.1. METHODEN 9

Die Eigenzust¨

ande Ψ(r1, ..., rN) des Systems k¨

onnen durch die L¨

osung der

Schr¨

odingergleichung erhalten werden:

HΨ(r1, ..., rN) = EΨ(r1, ..., rN) (1.2)

Die Grundidee der Dichtefunktionaltheorie ist, die komplizierte N-Elektronenwellenfunktion

Ψ(r1, ..., rN) und die dazugeh¨

orige Schr¨

odingergleichung durch die viel einfa-

chere Elektronendichte ρ(r) und ein entsprechendes Schema zu ihrer Berech-

nung zu ersetzen. Grundlegend f¨

ur die Dichtefunktionaltheorie ist ein Theorem

von Hohenberg und Kohn [5], nach dem eine eineindeutige Abbildung zwischen

einem externen Potential Vext(r) (z.B. das Coulomb-Potential der Atomkerne,

Vext(r) = VKern(r)) und der Elektronendichte ρ(r) besteht. Nach dem ersten

Hohenberg-Kohn Theorem ist ρ(r) durch das externe Potential Vext(r) eindeutig

bestimmt.

Durch das zweite Hohenberg-Kohn Theorem ist eine Beziehung der Grund-

zustandsenergie zu der N-Elektronendichte gegeben: Die exakte Grundzustands-

energie E0ist durch die Grundzustandsdichte ρG(r) eindeutig festgelegt. F¨

ur eine

Testfunktion ρ0(r), mit ρ0(r)>0 und

Zρ0(r)dr=N

gilt:

E[ρ0(r)] ≥E0.

Das Energiefunktional l¨

aßt sich in folgende Teilausdr¨

ucke zerlegen:

E[ρ(r)] = Ts[ρ(r)] + ZVKern(r)ρ(r)dr+1

2ZVH[ρ(r)]ρ(r)dr+Exc[ρ(r)] (1.3)

VKern(r) = −X

|RA−r|, VH[ρ(r)] = Zρ(r0)

|r−r0|dr0(1.4)

Ts[ρ(r)] ist die kinetische Energie eines nicht wechselwirkenden Elektronenga-

ses, Exc[ρ(r)] ist die sogenannte Austausch- und Korrelationsenergie und ZAsind

die Kernladungszahlen der Atome mit den Koordinaten RA.

Kohn und Sham [6] haben f¨

ur die Bestimmung von ρG(r) ein auf der Einf¨

uhrung

von Einelektronenorbitalen Ψi(r) basierendes Schema vorgeschlagen. Hierbei wird

ein System nicht wechselwirkender Elektronen betrachtet, die sich in einem ex-

ternen Potatial Vsbewegen, wobei sich die Dichte durch die Kohn-Sham (KS)

Orbitale Ψi(r) und Besetzungszahlen nigem¨

aß

ρ(r) =

occ

i=1

ni<Ψi(r)|Ψi(r)>

10 KAPITEL 1. EINF ¨

UHRUNG

darstellen l¨

aßt und die kinetische Energie durch

Ts[ρ(r)] =

occ

i=1

ni<Ψi(r)|−∇2|Ψi(r)>

gegeben ist.

Die Grundzustandsenergie kann durch Anwendung des Variationsprinzips er-

halten werden. Im Grundzustand des Systems gilt, daß die Variation der Energie

unter der Nebenbedingung der Orthogonalit¨

at der Wellenfunktionen

Z|Ψi(r)|2dr= 1,

verschwindet:

δΨ∗

i(r)

E[ρ(r)] −X

njjZ|Ψi(r)|2dr−1

= 0

Dies f¨

uhrt zu dem Eigenwertproblem (Kohn-Sham Gleichungen)

h−∇2+Vs[ρ(r)]iΨi(r) = iΨi(r),(1.5)

mit

Vs[ρ(r)] = VKern(r) + VH[ρ(r)] + Vxc[ρ(r)]

und

Vxc[ρ(r)] = δExc[ρ(r)]

δρ(r)

Vs[ρ(r)] ist ein effektives Einteilchenpotential, das die Wechselwirkung eines Elek-

trons mit den ¨

ubrigen Elektronen und den Atomr¨

umpfen beschreibt. Das Eigen-

wertproblem muß, beginnend mit einer Startdichte, iterativ gel¨

ost werden, da der

Hamiltonoperator die Dichte selbst enth¨

alt. Die Iteration wird bis zur Selbstkon-

sistenz durchgef¨

uhrt, d.h., bis sich die Dichten zwischen zwei Zyklen nur noch

innerhalb einer festgelegten Toleranz unterscheiden.

Die Bestimmung des Austausch- und Korrelationsfunktionals Exc ist eine

große Herausforderung im Rahmen der Dichtefunktionaltheorie. Die Hauptde-

fizite dieser Theorie liegen immer noch darin, daß eine angemessene Repr¨

asen-

tation der Austausch- und Korrelationseffekte durch entsprechende Funktionale

im großen und ganzen nicht befriedigend gel¨

ost ist. Auf einige Defizite wird auch

in diesen Ausf¨

uhrungen noch eingegangen. Die Form von Exc ist analytisch nicht

bekannt, wodurch jedes Funktional auf N¨

aherungen beruht. Die Bestimmung von

guten Approximationen f¨

ur Exc ist Gegenstand intensiver Forschung und kann

hier nicht weitergehend behandelt werden.

1.1. METHODEN 11

Die Effizienz der DFT Methoden ist in etwa mit der der HF Methode vergleich-

bar, wie oben schon erw¨

ahnt, es k¨

onnen in etwa Systeme mit 100 Atomen routi-

nem¨

assig behandelt werden. Mit Hilfe der Eigenwertgleichung kann die Gesamt-

energie auch folgendermaßen geschrieben werden:

E[ρ(r)] = X

nii−1

2ZVH[ρ(r)]ρ(r)dr−ZVxc[ρ(r)]ρ(r)dr+Exc[ρ(r)] (1.6)

A6=B

ZAZB

|RA−RB|

Diese Gleichung ist der Ausgangspunkt f¨

ur N¨

aherungen zur vereinfachten

Berechnung elektronischer Eigenschaften molekularer Systeme, von denen die

n¨

achsten Abschnitte handeln.

Die meisten heutigen Verfahren zur Bestimmung der elektronischen Struktur

molekularer Systeme basieren auf einer Einelektronenorbitalformulierung, wie sie

am Beispiel der Dichtefunktionaltheorie gerade vorgestellt wurde. In einem sol-

chen Schema wird der Grundzustand E0eines molekularen Systems durch eine

iterative L¨

osung der Eigenwertgleichungen bestimmt. Die L¨

osung beinhaltet da-

mit zwei Rechenschritte: Zun¨

achst m¨

ussen die Matrixelemtente des Hamilton

Operators und der ¨

Uberlappmatrix in einer System von Basisfunktionen berech-

net werden. Dazu werde die Wellenfunktionen Ψidurch die Basisfunktionen φµ

ΨiX

µφµ

dargestellt. Damit k¨

onnen dann die Hamilton- und ¨

Uberlappmatrixelemente in

dieser Basis berechnet werden. Im n¨

achsten Schritt werden dann die Eigenwerte

und Eigenfunktionen der Hamiltonmatrix bestimmt (Dieser Ablauf gilt auch f¨

die oben angesprochenen HF und post-HF Methoden). F¨

ur den ersten Schritt,

sind numerisch aufwendige Integrationen durchzuf¨

uhren. Die auftretenden Inte-

grale sind i. A. Mehrzentrenintegrale, d.h. sie h¨

angen von den Kernkoordinaten

mehrerer Atome ab. Daher skaliert der Aufwand zu ihrer Berechnung auch mit

der dritten (und bei den ab initio Methoden auch der vierten oder f¨

unften) Potenz

der Anzahl der im System vorhandenen Atome. Der zweite Schritt besteht in der

Diagonalisierung der Hamiltonmatrix. Der Zeitaufwand f¨

ur diesen Schritt steigt

mit der dritten Potenz der Dimension dieser Matrix. Die Dimension wiederum

h¨

angt von der Anzahl der Basisfunktionen ab. Bei Verwendung von Basiss¨

atzen,

die an den Atomkernenen lokalisiert sind, h¨

angt somit die Rechenzeit von der

Anzahl der Basisfunktionen pro Atom und der Anzahl der Atome ab. Die Anzahl

der Atome ist durch die Systemgr¨

osse gegeben. Somit l¨

asst sich die Rechenzeit

12 KAPITEL 1. EINF ¨

UHRUNG

durch die Basissatzgr¨

osse (Anzahl der Basisfunktionen pro Atom) beeinflussen.

Allerdings verschlechtert eine kleinere Basis die Qualit¨

at der Ergebnisse, so daß-

eine Balance zwischen Qualit¨

at der Rechnung und Rechenzeit gesucht werden

muß.

Es lassen sich zwei Ansatzpunkte f¨

ur schnellere Verfahren benennen. Der Eine

bezieht sich auf die Bestimmung der Matrixelemente. Hier gibt es vereinfachte,

semi-empirische und empirische Verfahren, die die Matrixelemente im Voraus

und m¨

oglichst unabh¨

angig von den zu betrachteten Systemen zu bestimmen ver-

suchen. Damit f¨

allt eine Berechnung der Matrixelemente weg, was einen signifi-

kanten Rechenzeitgewinn bedeutet, allerdings auch Abstriche bei der Pr¨

azision

der Ergebnisse nach sich zieht. Andererseits kann man versuchen, die Aufwen-

dige Diagonalisierung der Hamiltonmatrix zu umgehen. Die sogenannten O(N)

Verfahren versuchen genau dies, die Bestimmung der Grundzustandsenergie bei

einem Rechenaufwand, der linear mit der Systemgr¨

oße w¨

achst.

Semi-empirische Methoden

Semi-empirische Verfahren vernachl¨

assigen und n¨

ahern die aufwendig zu berech-

nenden Integrale der ab initio Verfahren und ersetzen sie anschließend durch zu

bestimmende Parameter (siehe etwa [7]). Ausgehend von der HF Theorie wurden

in den letzten 40 Jahren eine Reihe von semi-empirschen Verfahren entwickelt.

Bei den erfolgreichsten Vertretern dieser Klasse von Methoden, den AM1 [8] und

PM3 [9] Methoden werden die Parameter an experimentelle Daten von organi-

schen Molek¨

ulen gefittet. Eine zweite Klasse von semi-empirischen Verfahren, die

sogenannten tight-binding (TB) Methoden k¨

onnen als stationa¨

are Approximati-

on an die DFT verstanden werden [10]. Charakteristisch ist die Darstellung der

Gesamtenergie aus Gl. 1.6, welche aus einem Bandstrukturterm Ebs und einem

repulsiven Energiebeitrag, Erep hervorgeht. Ebs repr¨

asentiert den ersten Term auf

der rechten Seite von Gl. 1.6, die restlichen Terme werden in Erep zusammenge-

fasst. Die Kohn-Sham Orbitale werden dann in einer lokalen Basis dargestellt, und

die resultierenden Kohn-Sham Gleichungen werden nicht-selbstkonsistent gel¨

ost.

Da bei diesen beiden Klassen von semi-empirischen Methoden die Wechsel-

wirkungsintegrale nicht mehr berechnet werden m¨

ussen, ist der zeitlimitierende

Schritt durch die Diagonalisierung der Hamiltonmatrix gegeben. Die Rechen-

zeit dieser Methoden skaliert daher kubisch mit der Systemgr¨

osse. Mit semi-

empirischen Verfahren sind Systeme bis zu 1000 Atomen und Simulationszeiten

von 10-100 ps m¨

oglich, sie k¨

onnen damit etwa 10 fach gr¨

ossere Systeme behandeln

als DFT.

1.2. PROBLEMSTELLUNG AM BEISPIEL BAKTERIORHODOPSIN 13

Empirische Kraftfeld Methoden

Kraftfeldmethoden modellieren die Wechselwirkungen zwischen den Atomen durch

klassische Potentiale [11, 12, 13, 14]. Dies f¨

uhrt zu einer erheblichen Vereinfa-

chung in der Berechnung der interatomaren Wechselwirkungen. Der zeitlimitie-

rende Faktor besteht in der Berechnung der Coulomb-Wechselwirkung zwischen

den Atomen, die durch Punktladungen repr¨

asentiert werden. Daher w¨

achst die

Rechenzeit quadratisch mit der Systemgr¨

osse, solange nicht weitere N¨

aherun-

gen, wie etwa das Abschneiden der Coulomb-WW f¨

ur große Abst¨

ande zwischen

den Atomen oder Multipol-Ans¨

atze zur Verwendung kommen. Kraftfeldmethoden

sind allerdings aufwendig zu parametrisieren, da die Wechselwirkung der Atome

von der chemischen Umgebung abh¨

angt. So m¨

ussen beispielsweise verschiedene

Typen eines Elements unterschieden werden, je nach funktionaler Gruppe, in der

es auftritt. Eine sorgf¨

altige Kalibrierung hat jedoch zu erfolgreichen Anwendun-

gen dieser Verfahren gef¨

uhrt, so daß sie inzwischen zum Standardrepertoire bei

der Modellierung komplexer Strukturen geh¨

oren. Mit ihnen lassen sich Systeme

bis zu einigen 10000 Atomen und Simulationszeiten bis zu Nanosekunden reali-

sieren.

Empirische Kraftfelder werden sehr erfolgreich zur Modellierung der Struktur

und Dynamik von Proteinen und DNA eingesetzt. Bestimmte Probleme, wie etwa

das der Proteinfaltung k¨

onnen dennoch zur Zeit auch mit diesen nicht Methoden

behandelt werden, da die hierf¨

ur relevanten Zeitskalen die mit den Kraftfeldern

nicht erreicht werden.

Neben dem Problem ihrer Genauigkeit und ¨

Ubertragbarkeit zwischen ver-

schiedenen chemischen Umgebungen liegt die massgebliche Einschr¨

ankung dieser

Methoden darin, chemische Reaktionen nicht beschreiben zu k¨

onnen. Dies be-

trifft etwa enzymatische Katalyse, Protonentransport, Elektronentransport, die

Wechselwirkung von Materie mit elektromagnetischer Strahlung oder elektronisch

angeregte Zust¨

ande. Hier wird eine quantenmechanische Methode erforderlich

Die Herausforderungen an Simulationsmethoden bei der Modellierung solcher

Prozesse sollen stellvertretend am Beispiel des Protonentransfers in dem Mem-

branprotein Bakteriorhodopsin (bR) skizziert werden.

1.2 Problemstellung am Beispiel Bakteriorho-

dopsin

bR stellt das einfachste bekannte photosynthetische Zentrum dar. Es dient damit

als Modellsubstanz f¨

ur das Studium grundlegender photoinduzierter Prozesse in

der Biologie und wird sowohl experimentell als auch theoretisch intensiv unter-

sucht (einen zusammenfassenden ¨

Uberblick findet man in [15]). Dar¨

uber hinaus

14 KAPITEL 1. EINF ¨

UHRUNG

werden aber auch bedeutende Anwendungspotentiale dieses biologischen Systems

in der molekularen Elektronik gesehen.

Trifft ein Photon auf das Retinal Chromophor, welches das aktive Zentrum in

bR darstellt, so wird ein Protonentransferprozess in Gang gesetzt. Der dadurch

erzeugte pH Gradient ist der Ausgangspunkt f¨

ur die ATP Synthese in der Zelle.

Sowohl das Retinal, als auch die Beschaffenheit der Proteinumgebung des Retinals

sind wesentlich f¨

ur das Funktionieren dieses Prozesses verantwortlich. Durch ge-

zielten Austausch von Aminos¨

auren des bR mit Hilfe gentechnischer Methoden

konnte dies im Detail gezeigt werden. Die optischen Eigenschaften, aber auch

der Protonentransport h¨

angen sensitiv von den elektrostatischen und sterischen

Wechselwirkungen mit der Proteinumgebung ab. Eine genaue Aufkl¨

arung die-

ser Wechselwirkungen ist somit nicht nur f¨

ur das prinzipielle Verst¨

andnis von

bR, sondern auch im Hinblick auf dessen Modifizierung f¨

ur technologische An-

wendungen wichtig. Die Proteinumgebung hat somit einen großen Einfluß auf

die r¨

aumliche Struktur und die elektronischen Eigenschaften des Retinals. Oh-

ne diese Umgebung k¨

onnte das Retinal seine Rolle in der Lichtspeicherung und

umwandlung nicht ausf¨

uhren. Daher ist die Ber¨

ucksichtigung der Proteinumge-

bung in der theoretischen Modellierung des Retinals unumg¨

anglich. Die Gr¨

osse

des Proteins inklusiver seiner Hydrath¨

ulle (etwa 10000 Atome) w¨

urde f¨

ur den

Einsatz empirischer Karftfeldmethoden sprechen. Allerdings kann das konjugier-

te π-Elektronensystem des Retinals, dessen elektronisch angeregten Zust¨

ande wie

auch der Protonentransfer nur mit Hilfe von quantenmechanischen Methoden be-

schrieben werden.

Diese zwei Anforderungen, die Notwendigkeit der quantenmechanischen Be-

schreibung des Retinals zusammen mit der Gr¨

osse des gesamten bR Molek¨

uls

machen die Anwendung neuer Techniken zwingend.

Die hier beschriebene Problemstellung ist typisch f¨

ur eine Reihe biologisch in-

teressanter Fragen, die eine besondere Herausforderung an theoretische Methoden

stellen:

•Gegen¨

uber anorganischen und organischen Systemen heben sich die biologi-

schen durch die zu behandelnde Systemgr¨

osse ab. Sie umfassen von einigen

100 bis zu einige 1000 Atomen. Die Behandlung grosser Systeme kann ¨

ubli-

cherweise durch Ausnutzung der Systemsymmetrie (etwa bei Festk¨

orperpro-

blemen) vereinfacht werden, es m¨

ussen dann nur noch kleine Aussschnitte

des Systems behandelt werden. Dieser Ausweg ist jedoch aufgrund der feh-

lenden Symmetrie bei Biomolek¨

ulen nicht gegeben.

•Biomolek¨

ule enthalten eine Vielzahl von verschiedenen Atomtypen, von

Wasserstoff bis zu Metallen und ¨

Ubergangsmetallelementen. Dies erschwert

die theoretische Beschreibung im Vergleich zu den oft homogeneren Festk¨

orper-

systemen erheblich.

1.2. PROBLEMSTELLUNG AM BEISPIEL BAKTERIORHODOPSIN 15

Abbildung 1.1: Das bR Molek¨

ul und das Retinal

•Biomolek¨

ule liegen meist nicht in der Gasphase vor. Die Wechselwirkung

mit dem sie umgebenden L¨

osungsmittel, das zudem auch teilweise in sie

eindringt, ver¨

andert Struktur, Eigenschaften und Funktion dieser Molek¨

ule

signifikant. Daher muss eine Beschreibung die Umgebung stets mit ber¨

uck-

sichtigen. Diese bedeutet sowohl methodisch, als auch im Rechenzeitbedarf

einen signifikanten Mehraufwand.

•Die Stabilit¨

at biomolekularer Strukturen beruht zum grossen Teil auf Was-

serstoffbr¨

ucken und Van der Waals Bindungen, die hohe Anforderungen an

die Genauigkeit der Simulationsmethoden stellen, wie weiter unten noch

n¨

aher diskutiert wird.

•Auch ist eine genaue Beschreibung der Struktur und Energetik der che-

mischen Bindungen erforderlich. Zum Einen, um etwaige chemische Reak-

tionen beschreiben zu k¨

onnen, zum Anderen aber, um die Stabilit¨

at der

Konformationen des Proteins zu gew¨

ahrleisten.

•Eine zus¨

atzliche Schwierigkeit bei der Beschreibung von Br liegt darin, daß

die Reaktion Photoinduziert wird. Die Beschreibung von Prozessen im elek-

tronisch angeregten Zustand ist theoretisch sehr aufwendig, viele Methoden

erreichen die hierzu erforderlichen Genauigkeitsanforderungen nicht. In den

16 KAPITEL 1. EINF ¨

UHRUNG

letzten Jahren lag daher auch ein Schwerpunkt auf der Entwicklung neuer

Methoden zur Beschreibung angeregter Zust¨

ande.

Die Systemgr¨

osse und die Anforderungen bestimmen demnach die Wahl der

Simulationsmethode. Aus dem bisher gesagten d¨

urfte ersichtlich werden, dass

keine der unter 1.1 beschriebenen Methoden per se in der Lage ist, das skizzierte

Problem zu anzugehen.

Das quantenmechanisch zu beschreibenden Problem (Chromophor und einige

Aminos¨

auren aus der Proteinmgebung) umfasst etwa 100 Atome, ist daher zu

groß um mit ab initio Verfahren behandelt zu werden. F¨

ur derartige System-

gr¨

oßen m¨

ussen quantenmechanische Methoden bereitgestellt werden, welche die

chemische Genauigkeit und ¨

Ubertragbarkeit von Post-Hartree-Fock bzw. DFT-

Rechnungen mit der Effizienz und Einfachheit von st¨

arker gen¨

aherten semiempi-

rischen Verfahren verbinden. Im Rahmen meiner Dissertation [16] habe ich eine

solche Methode, das sogenannte SCC-DFTB Verfahren (’Self-consistent charge

density-functional thight-binding’; im Folgenden wird jedoch der Lesbarkeit hal-

ber die Methode nur noch mit DFTB abgek¨

urzt) entwickelt und anhand von

Tests und ersten Anwendungen f¨

ur organische Molek¨

ule und biologische Systeme

eingef¨

uhrt [17].

In den letzten Jahren wurde die Parametrisierung der Methode um neue

Atomtypen erweitert und die Methode vor allem in Hinblick auf die beschriebe-

nen Anforderungen getestet. Desweiteren wurden verschieden M¨

oglichkeiten un-

tersucht, die Beschreibung von H-Br¨

ucken im DFTB Formalismus zu verbessern.

Van der Waals Wechselwirkungen werden aufgrund der approximativen Behand-

lung der Elektronenkorrelationen in der Dichtefunktionaltheorie nicht ber¨

uck-

sichtigt. Da das DFTB Verfahren als N¨

aherung an die DFT auch ihre Probleme

ubernimmt, wurde ein empirischer Ansatz zur Einbeziehung von VdW (Disper-

sions) Kr¨

aften in die DFTB Methode verfolgt. Damit wurden die Grundsteine

gelegt, mit Hilfe des DFTB biologische Strukturen, wie das hier skizzierte bR,

erfolgreich behandeln zu k¨

onnen (Kapitel 2).

Die Prozesse im reaktiven Zentrum sind allerdings sensitiv abh¨

angig von der

Proteinumgebung, eine isolierte Behandlung des reaktiven Zentrums macht daher

wenig Sinn. Auf Workstations k¨

onnen mit der DFTB Methode komplexe Ato-

manordnungen mit bis zu ca. 500 Atomen behandelt werden. Die Bearbeitung

gr¨

oßerer Strukturen, wie sie bei der hier anvisierten Problemstellung auftreten,

kann durch drei Vorgehensweisen erreicht werden: eine massive Parallelisierung

der Methode [19], die Verwendung von sogenannten Order-N-Verfahren oder ein

sogenanntes QM/MM Verfahren.

•Die massive Parallelisierung erlaubt zwar eine Vergr¨

osserung der Systeme,

bedeutet aber nicht eine effizientere Nutzung der Computerressourcen, da

die Gesamtrechenzeit nur auf mehrerere Prozessoren aufgeteilt wird.

1.2. PROBLEMSTELLUNG AM BEISPIEL BAKTERIORHODOPSIN 17

•Eine vielversprechende Entwicklung stellen die sogenannten O(N) Algorith-

men dar. Die L¨

osung der quantenmechanischen Eigenwertgleichungen ska-

liert mit O(N3), d.h. die Rechenzeit steigt kubisch mit wachsender System-

gr¨

osse (N: Anzahl der Basisfunktionen), w¨

ahrend unter Verwendung von

O(N) Algorithmen ein lineares Anwachsen erreicht werden kann.

•Die Idee von QM/MM Methoden besteht darin, einen Teil des Systems

quantenmechanisch zu beschreiben und den Rest des Systems in der An-

kopplung mit einem empirischen Kraftfeld zu ber¨

ucksichtigen.

In den letzten Jahren wurden alle drei Optionen implementiert und k¨

urzlich so-

gar eine Kombination der drei Strategien realisiert. Diese Entwicklungen, ihre

Implementation und Anwendung werden in Kapitel 3 ausgef¨

uhrt. In Kapitel 5

werden die Entwicklungen zur Beschreibung elektronisch angeregter Zust¨

ande

und in Kapitel 6 Methoden, wie sie etwa zur Behandlung der Wechselwirkung

von Licht mit Materie notwendig sind, dargestellt. Zun¨

achst soll jedoch kurz auf

die charakteristischen Merkmale der DFTB Methode und deren Erweiterungen

und Test eingegangen werden.

18 KAPITEL 1. EINF ¨

UHRUNG

Kapitel 2

DFTB

Eine ausf¨

urliche Diskussion der DFTB Methdode findet sich in meiner Dissertati-

on [16] und in Ref. [17, 18]. Hier sollen nur knapp die wichtigsten Charakteristika

der Methode angesprochen werden.

Die DFTB Methode basiert auf einer Entwicklung des effektiven potentials

in dem Kohn-Sham Energiefunktional nach Ladungsdichtefluktuationen δn an

einer gegebenen Referenz-Elektronendichte n0[16]. Die Gesamtenergie in zweiter

Ordnung in den Ladungsfluktuation ist:

occ

ihψi|ˆ

H0|ψii+1

2ZZ 1

|~r −~r0|+δ2Exc

δn(~r)δn(~r0)n0!δn(~r)δn(~r0)d~rd~r0

−1

2ZZ0n0(~r)n0(~r0)

|~r −~r0|d~rd~r0+Z[Exc[n0]−Vxc[n0]] n0(~r)d~r +Eii (2.1)

H0=ˆ

H[n0] ist der effektive Kohn-Sham Hamiltonian, ausgewertet an der Refe-

renzdichte, und die ψisind die sogenannten Kohn-Sham Orbitale. Exc und Vxc

bezeichnen die Austausch- und Korrelationsenergie- und Potentialbeitr¨

age und

Eii ist die Kern-Kern Abstoßung. Um den Ausdruck f¨

ur die Gesamtenergie in der

DFTB Methode zu erhalten, werden noch folgende N¨

aherungsschritte ausgef¨

uhrt:

•Die Hamiltonmatrixelemente hψi|ˆ

H0|ψiiwerden in einer minimalen Basis

optimierter atomarer Wellenfunktionen φµdargestellt,

ψi=X

µφµ.

Um die φµzu erhalten werden die atomaren Kohn-Sham Gleichungen un-

ter Verwendung eines harmonischen Zusatzpotentials gel¨

ost [20], welches

die Aufgabe hat, die Wellenfunktionen r¨

aumlich st¨

arker zu lokalisieren. Die

Hamiltonmatrixelemente in dieser Basis H0

µν werden dann wie folgt berech-

net. Als Diagonalbeitr¨

age H0

µµ werden die atomare Eigenwerte angesetzt

20 KAPITEL 2. DFTB

und die nicht-Diagonalelemente H0

µν werden durch eine Zweizentrenn¨

ahe-

rung erhalten.

µν =hφµ|ˆ

T+Veff [n0

α+n0

β]|φνiµα, νβ,

Diese werden zusammen mit den ¨

Uberlappmatrixelementen Sµν vorab be-

rechnet und in Abh¨

angigkeit vom interatomaren Abstand Rαβ tabelliert.

Veff ist das effektive Kohn-Sham Potential und die n0

αsind die Elektro-

nendicheten der neutralen Atome α;ˆ

Tist der Operator der kinetischen

Energie.

•Die Ladungsdichetfluktuationen δn werden als Superposition atomarer Bei-

tr¨

age δnαgeschrieben,

δn =X

δnα,

welche durch die Ladungsdichtefluktuationen ∆qα=qα−q0

αausgedr¨

uckt

werden. q0

αist die Anzahl der Valenelektronen des neutralen Atoms αund

die Ladungen qαwerden durch eine Mulliken Populationsanalyse erhalten.

Der Term zweiter Ordnung in δn in Gl. (2.1) wird dann durch Pαβ γαβ∆qα∆qβ

gen¨

ahert. F¨

ur α6=β, wird γαβ durch eine analytische Funktion bestimmt,

die die Coulomb Wechselwirkung sp¨

arischer Ladungsverteilungen an den

Kernorten Rαund Rβrepr¨

asentiert. Die Diagonalterme γαα werden aus

∂E2

∂2qαerhalten.

•Die ’double counting’ Beitr¨

age und die Kern-Kern Abstoßung werden in

einem Term Erep zusammengefaßt. Dieser wird als Summe von Zweizen-

trenpotentialen U[Rαβ] gen¨

ahert.

Erep =X

α6=β

U[Rαβ ].

Mit diesen N¨

aherungen und Definitionen l¨

asst sich dann die Gesamtenergie

folgendermaßen ausdr¨

ucken:

Etot =X

iµν

µci

νH0

µν +1

αβ

γαβ∆qα∆qβ+Erep (2.2)

Zur Berechnung der Hamiltonmatrixelemente wird das Austausch- und Korrela-

tionspotential von Perdew, Burcke und Ernzerhof verwendet [21].

Unter Anwendung des Variationsprinzips erh¨

alt man gen¨

aherte Kohn-Sham

Gleichungen, die nun iterativ gel¨

ost werden m¨

ussen, da die Matrixelemente selbst

von der Eletronendichte ¨

uber die Mulliken Ladungen abh¨

angen.

2.1. TESTS 21

DFTB AM1 PM3 LDA GGA

Reaktionen (33)

F(E)(kcal/mol)5.0 - - 11.9 4.6

Geometrien (65 Molek¨

ule getestet)

F(R) (˚

A)0.013 0.017 0.012 - -

F(θ)(◦) 1.9 2.0 2.4 - -

Frequenzen (33 Molek¨

ule getestet)

F(ν) 6% 3-5% 3-5%

Tabelle 2.1: Reaktionsenergien, Bindungparameter und Schwingungsfrequenzen.

Vergleich der semi-empirischen und DFT basierten Methoden mit Bezug auf ex-

perimentelle Werte. LDA bezeichnet DFT in der lokalen Dichten¨

aherung, GGA in

der gradientenkorrigierten Variante. ¯

Fbezeichnet den mitteren absoluten Fehler

in Bezug auf die experimentellen Daten.

F¨

ur die atomaren Kr¨

afte wurden analytische Ausdr¨

ucke hergeleitet, die zwei-

ten Ableitungen nach den Kernkoordinaten werden numerisch ausgef¨

uhrt.

Die Zweizentrenpotentiale U[Rαβ] erh¨

alt man durch Subtraktion der DFT Ge-

samtenergie von der elektronischen Energie des DFTB Formalismus in Abh¨

angig-

keit vom zwischenatomaren Abstand r=Rαβ.

2.1 Tests

F¨

ur Anwendungen auf organische und biologische Systeme wurde die DFTB Me-

thode zun¨

achst f¨

ur die Atomtypen H, C, N, und O parametrisiert und an einem

umfangreichen Satz organischer Molek¨

ule getestet. Zur Evaluierung wurden ins-

besondere 33 Reaktionen zwischen kleinen organischen Molek¨

ulen untersucht und

die erhaltenen Reaktionsenergien mit experimentellen Werten und den Resulta-

ten anderer Methoden verglichen (Daten von Anzelm und Wimmer [22]). Weiter

wurden die Geometrien von 65 (Molek¨

ule aus dem Satz von Ref [8]) und die

Schwingungsfrequenzen von 33 Molek¨

ulen im Vergleich mit anderen Methoden

und experimentellen Daten untersucht. Es wurde jeweils der mittlere Fehler der

verschiedenen Methoden bez¨

uglich der experimentellen Daten ermittelt, wie in

Tabelle 2.1 wiedergegeben. Bez¨

uglich der Reaktionsenergien weist DFTB eine

ahnliche Abweichung von den experimentellen Daten auf, wie DFT-GGA [22],

was als exzellentes Ergebniss zu werten ist. DFT-LDA [22] schneidet hier durch

das wohlbekannte ’¨

uberbinden’ (overbinding) wesentlich schlechter ab. F¨

ur die

22 KAPITEL 2. DFTB

semi-empirischen Methoden AM1 [8] und PM3 [9] liegen f¨

ur dieses Set von Re-

aktionen explizit keine Daten vor. Allerdings liegen die Fehler f¨

ur molekulare

Bindungsenthalpien f¨

ur diese Methoden bei 4-6 kcal/mol, der Fehler in den Re-

aktionsenergien sollte daher von der gleichen Gr¨

osse sein.

Die Abweichungen in den Bindungsl¨

angen ¯

F(R) sind bei den drei approxima-

tiven Methoden sehr klein, und sind vergleichbar mit dem Ergebniss von ab initio

MP2 und DFT Methoden, die auch Abweichungen von etwa 0.015 (˚

A) aufweisen.

Ahnliches gilt f¨

ur die Vorhersage der Bindungswinkel θ. Allerdings liegen DFT

und MP2 Daten nicht f¨

ur einen solch grossen Satz von Molek¨

ulen vor, dennoch

kann man aus den vorhandenen Daten einen Eindruck von deren Genauigkeit

gewinnen [23].

Die Schwingungsfrequenzen werden mit einem mittleren Fehler von 6% etwas

schlechter beschrieben als durch DFT-GGA und DFT-LDA, aber wesentlich bes-

ser als bei den semi-empirischen AM1 und PM3 Methoden, die mittlere Fehler

von ¨

uber 10% aufweisen (siehe etwa [24]).

In den untersuchten Molek¨

ulen treten Bindungssituationen auf, wie sie auch

in gr¨

osseren organischen und biologischen Molek¨

ulen typisch sind. Die bisher vor-

gestellten Tests lassen damit auf eine sehr gute Beschreibung der kovalenten Bin-

dungen schließen. Allerdings m¨

ussen ¨

ahnliche Tests f¨

ur jeden weiteren Atomtyp,

der in die Berechnungen einbezogen werden soll, wiederholt werden. Dies betrifft

bei biologischen Anwendungen vor allem die Elemente Schwefel und Phosphor.

Desweiteren aber auch Alkalimetalle und ¨

Ubergangsmetallelemente. Weitere Ele-

mente f¨

ur die Methode zur Verf¨

ugung zu stellen, war daher ein Schwerpunkt der

Weiterentwicklung der Methode.

Dar¨

uberhinaus geben die Tests keinen Aufschluß ¨

uber die Qualit¨

at der Metho-

de bei der Beschreibung von nicht-bindenden Wechselwirkungen (’non-bonding’

interactions), wie etwa bei Wasserstoffbr¨

ucken oder Stapelungs Wechselwirkun-

gen (stacking).

2.2 Einbeziehung weiterer Atomtypen

Zur Beschreibung von Proteinen wird die Parametrisierung von Schwefel in DFTB

notwendig, da einige Aminos¨

auren, wie etwa Cystein oder Methionin dieses Ele-

ment enthalten. Schwefelatome finden ihre fundamentale Bedeutung bei der Sta-

bilisierung der r¨

aumlichen Struktur der Proteine durch die Bildung von Disul-

fidbr¨

ucken. Bei der Parameterenwicklung f¨

ur DFTB stellte sich heraus [24], daß

es zur Beschreibung der sogenannten hypervalenten Schwefelverbindungen wich-

tig ist, die minimale sp-Basis um die unbesetzten 3d Orbitale des Schwefels zu

erweitern. Probleme treten in der minimalen Basis auch schon bei SO2auf, des-

sen Bindungswinkel ohne die d-Funktionen um 13◦untersch¨

atzt wird. Oder auch

2.2. EINBEZIEHUNG WEITERER ATOMTYPEN 23

bei der S-C Bindung in (C6H5)2SO2, die ohne d-Basisfunktionen um 0.15 ˚

ubersch¨

atzt wird. In der minimalen Basis findet im Allgemeinen ein zu star-

ker Ladungstransfer zwischen S und O statt (was zu der Bindungsverl¨

angerung

f¨

uhrt), der durch die sogenannte ’back-donation’ von Elektronen in die Schwe-

fel d-Funktionen abgeschw¨

acht wird. Daher ist die Einbeziehung von d-Orbitalen

in die Basis unbedingt notwendig. Zwei Parameter sind hierf¨

ur zu festzulegen,

den atomaren Energieeigenwert des d-Orbitals, d, und den Kompressionsradius

f¨

ur das d Orbital, rd

0. Bei der Bestimmung der Kompressionsradien ist folgendes

zu beachten: Die unbesetzten Orbitale sind sehr langreichweitig, so daß f¨

ur das

d-Orbital eine st¨

arkere Kompression als bei den besetzten Orbitalen n¨

otig sein

k¨

onnte. Andererseits k¨

onnen die d-Orbital die minimale Basis aber auch als ’dif-

fuse’ Funktionen erweitern, was f¨

ur eher gr¨

ossere Kontraktionsradien sprechen

w¨

urde (und so auch in der Parametrisierung realisiert). Festzuhalten bleibt je-

doch, daß eine vorsichtige Optimierung der Kontraktionsradien n¨

otig ist und das

sonst robuste Schema

r0= 2 ∗rcov

in diesem Fall nicht umstandslos anzuwenden ist. Der atomare Eigenwert dbe-

stimmt, wie stark die unbesetzten Orbitale in der Verbindung mit den besetz-

ten Orbitalen mischen. Da die unbesetzten Orbitale durch die DFT Rechnungen

i.A. energetisch zu niedrig liegend berechnet werden, ist gegebenenfalls auch eine

Nachjustierung des Energiewerts notwendig, wovon im Falle des Schwefels abge-

sehen wurde, z.B. aber bei Silizium wichtig zu sein scheint.

Die Parametrisieung wurde an einem großen Satz organischer Molek¨

ule, ana-

log zu der Parametrisierung f¨

ur O, N, C, H, getestet. Die Abweichungen von den

experimentellen Daten f¨

ur Reaktionsenergien, Geometrien, Schwingungsfrequen-

zen etc. sind durchaus mit den entsprechenden Werten aus DFT Rechnungen in

der Qualit¨

at vergleichbar, vor allem aber zeigt sich eine deutliche Verbesserung

gegen¨

uber den getesteten semi-empirischen quantenchemischen Verfahren.

Die Erfahrungen bei der Parametrisierung des Schwefels konnten unmittelbar

auf die Parametrisierung von Phosphor ¨

ubertragen werden, wo ¨

ahnliche Parame-

terwerte verwendet wurden und eine vergleichbar gute Beschreibung molekularer

Eigenschaften erhalten wurde [25].

Desweiteren wurden Parameter f¨

ur die Wechelwirkung mit Zink entwickelt

[26]. Zink ist zwar durch die gef¨

ullte 3d und 4s Schale nicht kompliziert zu para-

metrisieren, doch waren einige Unsicherheiten bez¨

uglich der Parametrisierungs-

prozedur zu bedenken, deren L¨

osung paradigmatisch f¨

ur die Einbeziehung von 3d

Elementen ist. Zink ist ein biologisch wichtiges Ion, es liegt meist zweifach positiv

geladen vor. Zum Einen hat es katalytische, zum Anderen wichtige strukturelle

Funktionen, indem es Proteinkonformationen durch seine Anlagerung stabilisieren

kann. Obwohl es haupts¨

achlich im doppelt positivierten Zustand vorliegt, wurde

als Startdichte n0die Dichte es neutralen Zink Atoms verwendet. Dies ist deswe-

24 KAPITEL 2. DFTB

gen empfehlenswert, da dann die Energiebeitr¨

age, die in Erep zusammengefasst

sind, keine langreichweitigen Coulombanteile aufweisen, damit Erep insgesamt al-

so eher kurzreichweitig wird. Desweiteren wurde an neutralen, fiktiven Molek¨

ulen

der Form HnZnXHm(X=O, N, C, S, P, H) angepasst und die freien Valenzen

durch eine variable Anzahl von Wassertoffatomen n und m abges¨

attigt. Ange-

passt wurde Erep demnach f¨

ur eine ’kovalente’ Bindung, was nicht repr¨

asentativ

f¨

ur die Bindungssituationen in Proteinen ist. Daher wurden umfangreiche Tests

an f¨

ur die Biologie relevanten Modellkomplexen vorgenommen, z.B. an:

[Zn(H2O)n]2+,[Zn(H2O)n(NH3)m]2+usw.

Sowohl die Struktur, als auch die Energetik (Bindungsenergien, Deprotonierungs-

energien) werden hervorragend im Vergleich zu aufwendigen DFT Rechungen wie-

dergegeben. Mit dieser Entwicklung k¨

onnen nun eine Vielzahl von katalytischen

Prozessen in Zink Proteinen beschrieben werden. In Zukunft werden weitere 3d

Elemente (Cu, Fe) und Alkali Metalle in Anlehnung an die hier beschriebene Pro-

zedur parametrisiert, was eine Beschreibung vieler der interessanten katalytischen

Prozesse in der Biologie erlaubt.

2.3 Erweiterung fuer H-Bruecken

Im Vordergrund der Entwicklung steht die Anwendbarkeit der DFTB Methode

auf Biomolek¨

ule, speziell zun¨

achst auf Proteine und DNA Fragmente. Wichtig f¨

deren Struktur und Funktion sind die langreichweitigen Wechselwirkungen, die

uber Wasserstoff Br¨

ucken (H-Br¨

ucken) und Van der Waals (VdW) Wechselwir-

kungen vermittelt werden. Bisher wurden nur die Eigenschaften der kovalenten

Bindungen, Bindungsl¨

angen, Frequenzen und die Potentialtiefen (Energetik) ge-

testet. Die DNA Doppelhelix wird durch H-Br¨

ucken zwischen den beiden Str¨

angen

der Helix stabilisiert, ebenso Proteine, die durch H-Br¨

ucken entlang des Protein

R¨

uckrats in ihrer dreidimensionalen Struktur stabilisiert werden. Daher ist aber

eine genaue Beschreibung der Struktur als auch der Energetik von Wasserstoff-

br¨

ucken (H-Br¨

ucken) gebundenen Systemen unabdingbar. Die DFTB Methode

wurde deshalb an einem grossen Satz von H-Br¨

ucken gebundenen Modellsyste-

men getestet [27, 28, 29, 30]. Die erhaltenen Geometrien aller Komplexe sind ex-

zellent, wobei die Differenzen zu ab-initio Ergebnissen lediglich innerhalb deren

Streubreite (je nach Methode und verwendetem Basissatz) liegen. Die energeti-

sche Ordnung verschiedener Konformere der H-Br¨

ucken-gebundenen Komplexe

wird ebenfalls ausnahmslos korrekt wiedergegeben. Eine leichte Untersch¨

atzung

aller Bindungsenergien auf Grund der verwendeten minimalen Basis konnte durch

eine Erweiterung der Basis f¨

ur die H-Atome beseitigt werden [27].

2.4. DISPERSION 25

Abbildung 2.1:

2.4 Dispersion

Neben der gerade angesprochenen Stabilisierung der DNA durch H-Br¨

ucken zwi-

schen den Basen auf verschiedenen Str¨

angen, wie in Abbildung 2.1 dargestellt, ist

die Wechselwirkung zwischen den Stufen der DNA durch Van der Waals (VdW)

Kr¨

afte elementar (Abbildung 2.2). Diese sogenannte Stapelungswechselwirkung

resultiert aus den rein quantenmechanischen Energiebeitr¨

agen zur Gesamtener-

gie, der Austausch und Korrelationsenergie. Die Austauschwechselwirkung f¨

uhrt

zur Repulsion, wenn sich die Elektronendichten der zwei molekularen Fragmen-

te ¨

uberlappen (Pauli Repulsion). Der attraktive Teil, die sogenannte Dispersi-

onswechelwirkung resultiert aus der Elektronenkorrelation, die auf HF Niveau

nicht ber¨

ucksichtigt ist. Zahlreiche DFT Studien zu VdW Komplexen in den

letzten Jahren zeigten eine ernstzunehmende Schw¨

ache von DFT in der Beschrei-

bung dieser Wechselwirkung, sie wird beispielsweise mit LDA Funktionalen stark

ubersch¨

atzt, mit einigen GGA Funktionalen wiederum treten starke repulsive

Effekte auf, so daß DFT zum gegenw¨

artigen Zeitpunkt nicht ohne weiteres auf

solche Komplexe anwendbar ist.

Da in der HF Methode nicht enthalten, ist die Disperison erst auf dem MP2

Niveau beschreibbar. Gr¨

ossere DNA Framente sind aber mit MP2 aufgrund des

Rechenaufwandes nicht zu behandeln. Im Prinzip w¨

are es denkbar, daß semi-

emprische Methoden wie AM1 und PM3 hier erfolgreich einsetzbar sind, da sie

uber die Parametrisiung Elektronenkorrelationen enthalten. Es stellt sich aber

heraus, daß diese Methoden anstatt einer Attraktion der Basenpaare f¨

alschlicher-

weise eine Repulsion beschreiben, was zu einer Instabilit¨

at der DNA Doppelhelix

f¨

uhrt. Damit sind diese Methoden in ihrer gegenw¨

artigen Parametrisierung f¨

diese Anwendung unbrauchbar.

Da die DFTB Methode auf der DFT basiert, hat sie auch deren Schw¨

achen.

Dies betrifft insbesondere auch die fehlerhafte Ber¨

ucksichtigung der VdW (Di-

spersions) Kr¨

afte. Innerhalb DFT liegt das an dem approximativen Charakter

der Austausch- und Korrelationsfunktionale. Mit Kenntnis der exakten Funktio-

26 KAPITEL 2. DFTB

Abbildung 2.2:

nale w¨

are auch eine korrekte Behandlung der VdW Wechselwirkung m¨

oglich.

Durch die Parametrisierungsprozedur des DFTB werden aber genau die f¨

die Dispersionswechselwirkung verantwortlichen Energiebeitr¨

age vernachl¨

assigt,

selbst wenn sie exakt bekannt w¨

aren. Somit besteht die M¨

oglichkeit, diese Wech-

selwirkung in empirischer Weise auf das DFTB Schema aufzusetzten, ohne die

Korrelationsbeitr¨

age doppelt zu ber¨

ucksichtigen, wie es bei dem analogen Vorge-

hen f¨

ur DFT der Fall w¨

are (da hier die Korrelationsenergie teilweise ber¨

ucksich-

tigt ist). Hierzu wird eine empirische Formel f¨

ur die Dispersion, die eine 1/R6

Abh¨

angigkeit vom zwischenatomaren Abstand R aufweist und nur von atoma-

ren Polarisierbarkeiten abh¨

angt, zu dem DFTB Gesamtenergieausdruck addiert

(DFTB+VdW) [31]. Allerdings wird diese empirische Formel f¨

ur kurze intera-

tomare Abst¨

ande R, bei denen die Elektronendichten der Atome zu ¨

uberlappen

beginnen, ung¨

ultig. Eine bloße Addition der Formel f¨

uhrt daher f¨

ur kleine R zu

zu großer Attraktion, weshalb ein D¨

ampfungsterm eingef¨

uhrt wurde, der den Di-

spersionsteil f¨

ur kleine R zu Null f¨

uhrt.

In dieser Weise k¨

onnen die Energetik und Struktur der Basenpaarstapelung in

hervorragender ¨

Ubereinstimmung mit MP2 Rechungen beschrieben werden [31].

Damit wird der Weg zu einer quantenmechanischen Beschreibung von Fragmenten

der DNA Doppelhelix erstmals geebnet.

2.5 Struktur von Peptiden

Die Bausteine, aus denen Proteine und Polypeptide bestehen, sind die Ami-

nos¨

auren. Die Abfolge der Amins¨

auren entlang des Proteinr¨

uckgrats, die Sequenz,

wird dessen Prim¨

arstruktur genannt. Die dreidimensionale Strukur von Polypep-

tiden und Proteinen ensteht durch die r¨

aumliche Anordnung von sogenannten

Sekund¨

arstrukturelementen, wie etwa αHelizes und β-Faltblattstruktur. Diese

Sekund¨

arstrukturen enstehen durch ein Wechselspiel zwischen den Torsionsbar-

rieren f¨

ur Drehungen um einzelne Bindungen entlang des Proteinr¨

uckgrat und

internen Wasserstoffbr¨

ucken, die zwischen den einzelnen Aminos¨

auren der jewei-

ligen Sekund¨

arstruktur bestehen. Die Stabilit¨

at der Terti¨

arstruktur resultiert aus

2.5. STRUKTUR VON PEPTIDEN 27

den Wasserstoffbr¨

ucken oder Disulfidbr¨

ucken die zwischen den Sekund¨

arstruktu-

relementen ausgebildet werden.

Die Abbildung 2.3 zeigt zwei solcher Sekund¨

arstrukturen, die αHelix und

die βFaltblattstruktur. Diese beiden Strukturen k¨

onnen durch Drehungen der

Bindungswinkel im Peptid R¨

uckgrat ineinander ¨

uberf¨

uhrt werden.

Abbildung 2.3:

Wie von anderen approximativen Methoden (etwa AM1 und PM3) bekannt,

geben Tests an kleineren Molek¨

ulen nur wenig Aufschluss ¨

uber die Zuverl¨

assigkeit

der Methode f¨

ur ausgedehnte komplexe biologische Strukturen. AM1 und PM3

reproduzieren die Energetik und Struktur kleiner organischer Molek¨

ule mit guter

Genauigkeit, k¨

onnen aber die Peptidbindung nicht ad¨

aquat beschreiben. Diese

Methoden k¨

onnen dadurch die Energetik und Struktur verschiedener Konformere

von Polypeptiden nur unzureichend wiedergeben.

Daher ist eine genaue Evaluierung approximativer Verfahren (das bezieht sich

im engeren Sinne auf semi-empirische Methoden, generell aber auch auf ab initio

und DFT Methoden) f¨

ur die Anwendung auf Polypeptide und Proteine unbedingt

erforderlich.

Um die Zuverl¨

assigkeit des DFTB in der Beschreibung von Protein-Sekund¨

arstruk-

turen zu dokumentieren, wurden poly-Alanin Ketten im Vergleich zu ab initio,

DFT und semi-empirischen (AM1, PM3) Methoden untersucht. Die atomaren

Strukturen, Dipolmomente und relativen Energien verschiedener Konformere wer-

den auf der Grundlage des DFTB Verfahren mit einer hervorragenden Genauig-

keit wiedergegeben [30, 27, 32].

28 KAPITEL 2. DFTB

In einer weiteren Arbeit [34] wurden die Geometrien und relativen Energien

zweier Helixtypen untersucht, die in Proteinen am h¨

aufigsten auftreten: der α

und der 310 Helix [34]. Die αHelix wird durch interne Wasserstoffbr¨

ucken stabi-

lisiert, die zwischen der i-ten und der i+4-ten Aminos¨

aure geformt werden, bei

der 310 Helix bestehen diese zwischen der i-ten und i+3-ten Amninos¨

aure entlang

des Proteinr¨

uckgrats. Aufgrund der unterschiedlichen Verkn¨

upfungen ist die 310

Helix d¨

unner und l¨

anger als die αHelix, wie in Bild 2.4 zu sehen ist. Entgegen

Rechnungen mit klassischen Kraftfeldern, wurden mit DFTB und DFT die 310

Helizes als energetisch stabiler beschrieben. Auch sind die αHelizes f¨

ur Ketten

mit weniger als 8 Einheiten strukturell nicht stabil, sie relaxieren in 310 Helizes.

Ebenfalls ein Widerspruch zu den Rechnungen mit empirischen Kraftfeldern. Die

DFTB und DFT Energien und Strukturen stimmen sehr gut ¨

uberein, eine her-

vorragende Best¨

atigung f¨

ur die DFTB Methode, die damit in der Lage ist, diese

feinen Energieunterschiede, wie sie bei Konformations¨

anderungen auftreten, zu

beschreiben. Die αHelizes werden erst durch die Einbeziehung von L¨

osungsmit-

teleffekten stabiler als die 310 Helizes. Diese Problemstellung wird weiter unten

noch fortgef¨

uhrt.

Abbildung 2.4:

Um die Anwendbarkeit f¨

ur Molekulardynamiksimlationen zu evaluieren, wur-

den die Schwingungsfrequenzen von kleinen Peptiden in Abh¨

angikeit von der

Konformationen mit der DFTB-Methode im Vergleich mit MP2 und DFT unter-

sucht. F¨

ur die Standardabweichungen der Frequenzen von den experimentellen

Werten mit MP2 und DFT werden 3.0%, und 4.4% ermittelt, w¨

ahrend das DFTB

Verfahren mit 6.7 % nur unwesentlich schlechter abschneidet [33]. Speziell ist die

Methode in der Lage, die sehr feinen Differenzen in den Freuqenzen zwischen den

verschiedenen Konformeren zu reprozieren.

2.6. STRUKTUR DES RETINALS IN BAKTERIORHODOPSIN 29

Desweiteren wurden IR und Raman Intensitit¨

aten mit Hilfe des DFTB auf

zwei N¨

aherungsstufen berechnet. Zum Einen wurden die Geometrien und Hes-

sische mit dem DFTB Verfahren gewonnen, und die IR und Raman Tensoren

f¨

ur die DFTB-Geometrie mit Hilfe der DFT berechnet. Dieses Hybridverfah-

ren erm¨

oglicht die Berechnung von IR und Raman Intensit¨

aten, die sehr gut

mit den vollen ab initio Werten ¨

ubereinstimmen [33]. Ein weitergehender N¨

ahe-

rungsschritt besteht darin, die IR Intensit¨

aten ganz in der DFTB N¨

aherung un-

ter Verwendung der Mulliken N¨

aherung f¨

ur die (Ableitung der) Dipolmomente

zu berechnen. Dieses Vorgehen erm¨

oglicht sehr effizient die Berechnung von IR

Spektren. Qualitativ werden die Spektren reproduziert, man findet allerdings ei-

ne ¨

Ubersch¨

atzung der Intensit¨

at der intensiven Moden relativ zu den Moden mit

geringer IR Intensit¨

at. Dennoch zeigen die Spektren, eventuell nach einer Ska-

lierung der intensiven Moden, alle qualitiv wichtigen Merkmale der ’fingerprint’

Regionen, die zur Identifikation der Konformation notwendig sind [27].

2.6 Struktur des Retinals in Bakteriorhodopsin

F¨

ur eine Vertiefung des Verst¨

andnisses der Effekte der Proteinumgebung auf die

strukturelle und elektronische Eigenschaften des Chromophors in dem Protein

Bakteriorhodopsin (siehe Einleitung) m¨

ussen alle Einflußfaktoren aus der Protei-

numgebung mit einer hinreichenden Genauigkeit in den Rechnungen ber¨

ucksich-

tigt werden. Die alleinige Einbeziehung einzelner Aminos¨

auren in der Bindungsta-

sche kann den Einfluß des Gesamtproteins einschließlich seiner L¨

osungsmolek¨

ule

auf die RSB nur unvollst¨

andig beschreiben. Die Wechselwirkung der Proteinum-

gebung mit dem Chromophor ist vielmehr durch kollektive Effekte gepr¨

agt, die

durch die Dynamik einer Vielzahl von geladenen, polaren bzw. polarisierbaren

molekularen Gruppen hervorgerufen werden. Die Einbeziehung all dieser Grup-

pen in dynamische Simulationen wird auch in der nahen Zukunft bei weitem die

Leistungsf¨

ahigkeit von ab-initio Rechnungen ¨

ubersteigen. Andererseits w¨

aren die

traditionellen semi-empirischen AM1 und PM3 Verfahren zwar aufgrund ihrer Re-

chengeschwindigkeit f¨

ur die Behandlung entsprechender Fragestellungen geeignet,

versagen aber vollst¨

andig durch zu große Defizite in der erreichbaren Genauigkeit.

Deshalb wurden ab-initio und DFT Rechnungen als Datenbasis f¨

ur umfangrei-

che Tests der DFTB Methode genutzt. Zun¨

achst wurden die Geometrien, Proto-

nenaffinit¨

aten und Isomerisierungsbarriern f¨

ur einige Modelle des Chromophors

berechnet. Desweiteren wurde der Einfluß einzelner Aminos¨

auren und Wasser-

molek¨

ule der Proteinumgebung auf die Geometrie und Ladungsverteilung des

Chromophors untersucht. Alle Ergenisse der DFT Rechnungen wurden dabei

hervorragend reproduziert. Deprotonierungsenergien, Rotationsbarrieren und Ge-

metrie¨

anderungen werden ebenso mit hoher Pr¨

azision wiedergegeben wie Trends

30 KAPITEL 2. DFTB

dieser Gr¨

oßen in Abh¨

angigkeit von der L¨

ange des betrachteten Polyenes oder von

der Einbeziehung einzelner geladener Gruppen in der Umgebung [35]. Die DFTB

Methode ist damit gegenw¨

artig das einzige N¨

aherungsschema, mit dem in einer

mit semiempirischen Verfahren vergleichbaren Effizienz die Wechselwirkung von

Retinoiden in ihrer Proteinumgebung voll quantenmechanisch mit hoher Genau-

igkeit behandelt werden kann.

Kapitel 3

Grosse Systeme

3.1 Kombination von Methoden:QM/MM

Die Idee von QM/MM Methoden besteht darin, einen Teil des Systems quan-

tenmechanisch zu beschreiben (etwa das Reaktionszentrum, wenn eine genauere

Beschreibung n¨

otig ist oder wenn chemische Reaktionen modelliert werden sol-

len), und den Rest des Systems mit einem empirischen Kraftfeld.

Die QM/MM Methodik wurde urspr¨

unglich eingef¨

uhrt, um die enzymatische

Katalyse untersuchen zu k¨

onnen [44, 45]. In den letzten zwei Jahrzehnten wur-

den diese Methodik von verschiedenen Gruppen weiterentwickelt, wobei im QM

Bereich haupts¨

achlich semi-empirische Methoden eingesetzt wurden. Erst in den

letzten Jahren werden vermehrt auch ab initioi (vornehmlich HF) und DFT Me-

thoden im QM Bereich verwendet.

Am Beispiel des bR (siehe Einf¨

uhrung) soll die QM/MM Methode kurz erl¨

autert

werden. Die Proteinumgebung hat einen großen Einfluß auf die r¨

aumliche Struk-

tur und die elektronischen Eigenschaften des Retinals. Ohne diese Umgebung

k¨

onnte das Retinal seine Rolle in der Lichtabsorption- und Umwandlung nicht

ausf¨

uhren. Daher ist die Ber¨

ucksichtigung der Proteinumgebung in der theoreti-

schen Modellierung des Retinals unumg¨

anglich. Das konjugierte π-Elektronensystem

des Retinals kann aber nur mit Hilfe von quantenmechanischen Methoden be-

schrieben werden. Diese zwei Anforderungen, die Notwendigkeit der quantenme-

chanischen Beschreibung des Retinals zusammen mit der Gr¨

osse des gesamten

bR Molek¨

uls machen die Anwendung moderner QM/MM Techniken, wie oben

vorgestellt, zwingend. Hierbei kann das Retinal zusammen mit einigen wichtigen

funktionalen Gruppen der Proteinumgebung quantenmechanisch, der Rest des

Molek¨

uls aber mit den wesentlich weniger aufwendigen empirischen Kraftfeldme-

thoden behandelt werden. Auf diese Weise kann der Einfluß der Proteinumgebung

auf das Retinal, sowie der Protonentransportrozess untersucht werden.

Die Kopplung zwischen beiden Beschreibungsweisen erfolgt ¨

uber ein link-

32 KAPITEL 3. GROSSE SYSTEME

atom-Konzept und ber¨

ucksichtigt sowohl elektrostatische als auch Van-der-Waals

Wechselwirkungen.

Die Gesamtenergie f¨

ur den QM/MM Formalismus lautet:

Etot =EQM +EMM +EQM−MM ,(3.1)

wobei EQM als Energie des quantenmechnischen Teils durch die DFTB Gesamt-

energie repr¨

asentiert wird, der Energiebeitrag EMM durch die Energiefunktion des

klassischen empirischen Kraftfeldes gegeben ist. EQM−MM beschreibt die Kopp-

lung der beiden Teilgebiete.

Das Schema wurde zuerst f¨

ur den Fall implementiert [28], daß durch die

Systemgrenzen keine kovalenten Bindungen verlaufen. In diesem Fall besteht

EQM−MM nur aus Coulomb und Van der Waals Wechselwirkungen. Im DFTB

Schema l¨

asst sich dies besonders einfach ausdr¨

ucken:

EQM−MM =X

QA∆qk

Rak

+EV dW

QAsind die Ladungen auf den MM Atomen, ∆qkdie Mulliken Ladungen

aus dem DFTB und RAk ist der Richtungsvektor zwischen Atom A und Atom k.

EV dW gibt die VdW Wechselwirkung zwischen den Atomen aus den beiden Berei-

chen wieder, wobei auch f¨

ur den QM Bereich die Parameter des MM Kraftfeldes

verwendet werden.

Rechnungen f¨

ur H-Br¨

ucken gebundene Komplexe zeigen eine sehr gute ¨

Uber-

einstimmung mit ab-initio Resultaten [28] f¨

ur Geometrien, Bindungsenergien und

die energetischen Ordnung verschiedener Konformere. Es konnten auch erstmals

die oben diskutierten Polyalanin Helices in einer Umgebung von Wassermolek¨

ulen

untersucht werden, wobei die Wassermolek¨

ule in dem MM Bereich liegen und die

Helix Atome durch QM beschrieben werden. Allerding wurden nur Geometrien-

optimierungen durchgef¨

uhrt, was die Aussage der Rechnungen sehr limitiert hat.

Dennoch konnte man eine energetischer Favorisierung der αHelix feststellen. Die

Analyse ergab, daß die etwas gr¨

oßeren Dipolmomente der αHelix die Energieer-

niedrigung im Vergleich zur 310 Helix erkl¨

aren k¨

onnen.

Durchschneidet die QM-MM Grenze eine kovalente Bindung, wird das Pro-

blem etwas komplexer. Die einfachste L¨

osung besteht in dem sogenannten ’link

atom’ Konzept. Hierbei wird die QM Region durch ein Wasserstoffatom ab-

ges¨

attigt, das entlang der Achse der durchgeschnittenen Bindung im entsprechen-

den Abstand von dem QM Atom plaziert wird. Die Wechselwirkung zwischen den

Atomen der durchgeschnittenen Bindung wird durch das empirische Kraftfeld ver-

mittelt. Das ’link atom’ Wasserstoff ist dabei ein Fremdk¨

orper, dessen Einfluß auf

das System in vielen F¨

allen vernachl¨

assigt werden kann. Dies wurde auch f¨

ur die

DFTB Methode intensiv getestet [29]. Es wurden eine Reihe von Protonentrans-

ferprozessen in dem Enzym Triosephosphate Isomerase (TIM) durchgef¨

uhrt [29].

3.2. O(N) METHODEN 33

Im Vergleich mit DFT-Rechnungen ergibt sich eine mittlere Abweichung von 2-

4 kcal/mol, w¨

ahrend die semi-empirischen Methoden AM1 und PM3 mittlere

Abweichungen von 11 kcal/mol aufweisen. Zus¨

atzlich hervorzuheben ist, daß die

Geometrien der ¨

Ubergangszust¨

ande sehr gut mit den DFT-Geometrien ¨

uberein-

stimmen. In der gleichen Arbeit wurde auch die Stabilit¨

at der Polyalanin Helices

in L¨

osung untersucht, diesmal allerdings mit Hilfe einer QM/MM MD Simulati-

on. Hier zeigte sich, daß die 310 Helix mit acht Aminos¨

aureeinheiten in L¨

osung

dynamisch instabil ist, sie geht innerhalb von 10 ps Simulationszeit in eine α

Helix ¨

uber.

Leztere Simulationen dokumentieren die Begrenztheit der QM/MM Metho-

den. Simulationen mit einem QM Bereich dieser Gr¨

oße werden sehr rechenzei-

tintensiv, wobei die meiste Rechenzeit im QM Teil bei der Diagonalisierung ver-

braucht wird. F¨

ur gr¨

ossere Systeme empfiehlt sich daher die Verwendung von

O(N) Algorithmen in der QM Region.

3.2 O(N) Methoden

Die Entwicklung von linear skalierenden Methoden stellt einen Durchbruch in

der Modellierung großer biologischer Systeme dar. W. Yang entwickelte als erster

einen solchen Algorithmus, das sogenannte ’divide-and-conquer’ (DAC) Schema

[38, 39]. In den folgenden Jahren wurden weitere O(N) Verfahren vorgeschlagen,

die auf ¨

ahnlichen Prinzipien beruhen. Einen zusammenfassenden ¨

Uberblick ¨

uber

die Methoden geben Yang und Perez-Jorda [40]. Geodecker [41] und Galli [42].

Kombiniert man O(N) mit approximativen Verfahren wie AM1, PM3 oder DFTB,

wird die Behandlung von Systemen mit mehreren 1000 Atomen m¨

oglich [43].

In den meisten O(N) Methoden wird die r¨

aumliche Lokalisierung der Elektro-

nendichte ausgenutzt, um zu einer linearen Abh¨

angigkeit der Rechenzeit von der

Systemgr¨

osse zu kommen. In diesem Fall l¨

aßt sich die elektronische Struktur des

Systems aus lokalen Gr¨

oßen, etwa aus den lokalisierten molekularen Orbitalen

oder einer lokalisierten Dichtematrix aufbauen. Eine Dichtematrixformulierung

der Methode von Yang und Lee vereinfacht den Algorithmus weitergehend [39].

Wenn die Dichtematrix

ρµν =X

nici

µci

ν(3.2)

r¨

aumlich lokalisiert ist (nicht in Metallen), l¨

asst diese sich aus Beitr¨

agen von lo-

kal ermittleten Dichtematritzen aufbauen. Dies geschieht folgendermassen: Man

teilt das zu behandelnde System in M Untersysteme, die sich ¨

uberlappen. Dann

wird das Eigenwertproblem durch Diagonalisierung f¨

ur die Untersysteme gel¨

ost,

man erh¨

alt die ci

µf¨

ur jedes Untersystem. Anstatt also das Eigenwertproblem f¨

das Gesamtsystem zu l¨

osen, l¨

ost man es f¨

ur die M Untersysteme. Aus den lokal

34 KAPITEL 3. GROSSE SYSTEME

berechneten Koeffizienten ci

µwird dann die Dichtematrix Gl. (3.2) konstruiert.

Um die Besetzungszahlen der Eigenzusta¨

ande der Subsysteme zu finden, wird

das chemische Potential in die Rechnung eingef¨

uhrt. Die Eigenzust¨

ande werden

bis zu diesem Energiewert besetzt, und das chemische Potential wird solange ite-

riert, bis die berechnete Elektronenzahl mit der Gesamtelektronezahl des Systems

identisch ist. Da mit wachsender Systemgr¨

osse die Anzahl der Untersysteme line-

ar anw¨

achst, w¨

achst somit der Gesamtrechenaufwand linear. F¨

ur kleine Systeme

ist allerdings die DAC Methode aufwendiger als die direkte Diagonalisierung, da

bei DAC mehrere Subsysteme diagonalisiert werden m¨

ussen. Erst ab einem ’cross

over’ Punkt wird die DAC Methode effizienter. Dieser liegt bei einer Kombination

mit DFTB bei etwa 200 Atomen (bei etwa 40% Wasserstoffatomen, der Rest der

Atome wird mit minimaler sp Basis beschrieben). Ein Problem entsteht nur an

den Grenzen der Untersysteme, wenn die ¨

Uberlappungen zu klein gew¨

ahlt sind,

so daß Beitr¨

age zur Dichtematrix vernachl¨

assigt werden, die noch signifikant f¨

deren Darstellung sind. Die Qualit¨

at de Methode wird somit ¨

uber die Gr¨

osse

der ¨

Uberlappungsregion bestimmt. Ist diese sehr groß so wird die Methode sehr

genau, aber die Rechenzeit steigt. Ist sie zu klein, werden die Energien und ato-

maren Kr¨

afte ungenau F¨

ur die DFTB Methode wurde der ¨

Uberlappbereich durch

ausprobieren bestimmt. Dabei wurden Alanin Polypeptide mit 20 Aminos¨

aure-

einheiten (212 Atome) verwendet und die Energien und Kr¨

afte f¨

ur verschieden

Konformationen mit verschieden großen ¨

Uberlappbereichen ermittelt. Ein Unter-

system wurde als eine Aminos¨

aure definiert. Zur Bestimmung der ¨

Uberlappbe-

reiche wurden zwei Kriterien verwendet: Zum Einen wurden N Nachbarn jedes

Atoms in Subsystem M ber¨

ucksichtigt, zum Anderen alle Atome im Abstand Rb

eines jeden Atoms aus dem Subsystem M. Die Vereinigung dieser beiden Mengen

ergab die ¨

Uberlappmenge. In Tabelle 3.1 sind die Ergebnisse der Testrechnungen

dargestellt. Zwischen den Konformationen , wie durch Diagonalisierung erhalten,

bestehen Energieunterschiede von etwa 10 kcal/mol. Eine Kombination Rb= 3.5

Aund N=3 ergibt einen Fehler, der in der Gr¨

ossenordnung der Energeidifferen-

zen zwischen den Konformeren liegt. Erst die Kombination Rb= 4.0 ˚

Aund N=

4 f¨

uhrt zu akzeptabel kleinen Fehlern. Ein Fehler von 5 kcal/mol enspricht ei-

nem Fehler von 1meV pro Atom im System, was 0.03 % der koh¨

asiven Energie

entspricht.

Die verschiedenen Konformationen ben¨

otigen unterschiedliche Rechenzeiten

im DAC Formalismus, da je nach Kompaktheit des Molek¨

uls die ¨

Uberlappregio-

nen verschieden groß ausfallen. Bei der gew¨

ahlten Gr¨

oße der ¨

Uberlappregionen

erh¨

alt man eine Halbierung der Rechenzeit im Vergleich zur Diagonalisierung.

Als erste Anwendung wurde das Protein Crambin betrachtet (Abb. 3.1). Die-

ses Molek¨

ul enth¨

alt 46 Aminos¨

auren (=46 Subsysteme in der Rechnung) und

insgesamt 639 Atome.

Die volle Diagonalisierung ben¨

otigt etwa 15 Minuten (auf einer Dec alpha

3.3. KOMBINIERTE QM/MM O(N) 35

Konformation Ceq

7Cext

5310 αRechenzeit[s]

Diagonalisierung 0.0 34.5 -11.9 -32.0 85

Rb(˚

A) N

3.5 3 16.2 13.4 7.1 5.3 16-37

4.0 4 3.1 2.3 1.6 3.2 22-57

5.0 4 0.4 2.1 0.7 1.4 35-110

Tabelle 3.1: Die relativen Energien (kcal/mol) der verschiedenen Konformere,

wie sie durch Diagonalisierung erhalten werden (2. Zeile). In den folgenen Zeilen

werden die Fehler in Bezug auf die exakten Energien in Abh¨

angigkeit von Rbund

N wiedergegeben. Die letzte Spalte gibt die Rechenzeiten (s) wieder.

ev6), w¨

ahrend die DAC Rechnung mit den obigen Parametern nur eine Minute

be¨

otigt. Die Geometrie des Molek¨

uls wurde einer Geometrieoptimierung unterzo-

gen, wobei die Startstruktur, die der Kristallstruktur entspricht, sehr gut erhalten

wird.

3.3 Kombinierte QM/MM O(N)

Da die biologischen Molek¨

ule meist in (w¨

assriger) L¨

osung oder in der Zellumge-

bung sind, machen Rechnungen an isolierten Molek¨

ulen wenig Sinn: Umgebungs-

einfl¨

usse k¨

onnen die elektronische und geometrische Struktur stark beeinflussen.

Um auch Simulationen zu Systemen mit grossen QM Regionen durchf¨

uhren zu

k¨

onnen, wurde das DAC Verfahren in einen QM/MM Algorithmus implementiert.

Als erstes Testbeispiel wurde das Protein Crambin, umgeben von etwa 3000

Wassermolek¨

ulen betrachtet. Das Crambin Molek¨

ul wurde hierbei mit der linear

skalierenden DFTB-Methode behandelt, die Wassermolek¨

ule mit einem empiri-

schen Kraftfelder. F¨

ur 2 Picosekunden Molekulardynamiksimulation dieses Sy-

stems von 7000 Atomen, ben¨

otigt eine DEC/alpha Workstation ca. 3 Tage [46].

Dieses Programmsystem wurde auch paralellisiert, so daß nun f¨

ur Systeme der

oben genannten Gr¨

ossenordung MD Simulationen ¨

uber 100ps in etwa 2 Wochen

auf 16 Prozessoren durchgef¨

uhrt werden k¨

onnen.

Die Simulation ergab, daß Crambin strukturell nicht stabil ist. Crambin enth¨

alt

zwei αHelices, und beide gehen w¨

ahrend der MD in 310 Helices ¨

uber. Nachfolgend

wurden Simulationen durchgef¨

uhrt, in denen die Dispersionswechselwirkung mit

in die Rechnung einbezogen wurde. Damit war die Struktur von Crambin sta-

bil. Im Detail wird dieses Ergbnis weiter unten diskutiert. F¨

ur die geometrischen

Parameter wurde eine hervorragende ¨

Ubereinstimmung mit den experimentellen

Daten gefunden, eine bessere ¨

Ubereinstimmung sogar als die empirischen Kraft-

feldmethoden erzielten, die zum Vergleich herangezogen wurden. Damit wurde

36 KAPITEL 3. GROSSE SYSTEME

Abbildung 3.1: Das Protein Crambin

gezeigt, daß die DFTB + Dispersion Methode sehr gut geeignet ist, Proteinstruk-

turen zu beschrieben. Weiter wurde ein nichtlokaler Ladungstransfer innerhalb

des Proteins gefunden, der von empirischen Potentialen nicht beschreiben werden

kann. Inwieweit dieser Ladungstransfer f¨

ur die Struktur wichtig ist, oder ob es

auch Ladungstransfer zwischen Protein und Wasser gibt, bleibt nachfolgenden

Untersuchungen vorbehalten.

Kapitel 4

Stabilit¨

at von Protein

Sekund¨

arstrukturen

Eine Reihe von Arbeiten haben sich mit der Stabilit¨

at von Sekund¨

arstrukturele-

menten kurzer Polypeptide besch¨

aftigt [34, 32, 27, 28, 29], insbesondere mit der

relativen Stabilit¨

at der 310 vs. der αRHelix.

Unter Verwendung von DFT und DFTB Methoden wurden die Geometrien

von 310 und αRHelices optimiert. Dabei stellte sich heraus, daß zwar die 310

Strukturen stabil f¨

ur alle untersuchten Kettenl¨

angen sind, die αRHelices aber

erst ab einer bestimmten L¨

ange von etwa acht Aminos¨

aureeinheiten. F¨

ur N=8

(N: Anzahl der Aminos¨

aureeinheiten) ist nur eine αRtypische i →i+4 H-Br¨

ucke

zu finden, so daß der Großteil der Helix, n¨

amlich die beiden Enden, in der 310

Konformation vorliegt.

In Proteinen, wie etwa Crambin, findet man aber αRFragmente von der L¨

ange

N=6-8. Daher wurde die relative Stabilit¨

at der zwei Helices f¨

ur Dihedralwinkel

in dem Proteinr¨

uckgrat untersucht, die typischerweise in Proteinen aufzufinden

sind. W¨

ahrend der Geometrieoptimierung wurden die Dihedralwinkel bei den

idealen Werten festgeghalten. Dabei stellte sich heraus, daß die αRerst f¨

ur N>>

10 stabiler als die 310 Helix wird. Die Rechnungen beziehen sich auf Differenzen

von Gesamtenergien, nicht von freien Energien. Allerdings scheinen die Entropie-

beitr¨

age zus¨

atzlich die 310 Helix zu favorisieren. Dies steht im Gegensatz zu den

experimentellen Ergebnissen.

Die Stabilisierung der αRgegen¨

uber der 310 Helix findet durch Umgebungsein-

fl¨

usse statt. So fanden wir durch die oben beschriebenen QM/MM Rechnungen,

daß f¨

ur N= 8 die αRHelix in w¨

assriger L¨

osung stabiler wird. Der Wert der dielek-

trischen Konstante von Wasser liegt bei =78, die Vakuumrechnungen unterstel-

len einen Wert von = 1, w¨

ahrend in der Proteinumgebung ein Wert von etwa

= 4-8 angenommen werden kann. Unterstellt man eine monotone Abh¨

angigkeit

der relativen Stabilit¨

aten von dem Wert des Dielektrikums, so bewirkt die Pro-

38 KAPITEL 4. STABILIT ¨

AT VON PROTEIN SEKUND ¨

ARSTRUKTUREN

teinumgebung in jedem Fall eine Stabilisierung der αRKonformation relativ zur

310 Helix. Die Simulation an dem Crambin Molek¨

ul hat jedoch ergeben, daß der

Einfluß der Proteinumgebung nicht ausreicht, um die zwei αRHelices in Crambin

zu stabilisieren. Diese gehen im Verlauf der Dynamik innerhalb von 30 ps in 310

Konformationen ¨

uber, im Gegensatz zu den QM/MM Rechnungen in Wasser [29].

Auf der anderen Seite ergab eine MD Simulation unter Einbeziehung der Di-

spersionskr¨

afte [46] in dem DFTB Formalismus strukturell stabile αHelizes in

Crambin. Die Dispersion (Korrelation) scheint daher die αHelix gegen¨

uber der 310

Helix zu favoriserieren. Die Dispersionskr¨

afte haben eine reine 1/R6Abh¨

angig-

keit, kompaktere Molek¨

ule, wie die αRStrukturen weisen daher eine gr¨

ossere

Dispersionswechselwirkung auf als die eher l¨

angliche und d¨

unnere 310 Konforma-

tion.

Um dieses Ergebnis etwas genauer zu verstehen, wurden abermals Rechnun-

gen an Polyalanin Helices durchgef¨

uhrt. Wiederum wurden die Geomtrien beider

Helices f¨

ur N=3,5,8 bei festgehaltenen Dihedralwinkeln des Peptidr¨

uckrades op-

timiert. Es wurden HF, MP2, DFTB und DFTB+Dis Energien f¨

ur die jeweilige

Konformation bestimmt. Mit MP2 sind wegen der langen Rechenzeit Geome-

trieoptimierungen schwer, bzw. z. Z. nicht durchf¨

uhrbar. Die Ergebnisse sind in

Tabelle 4 wiedergegeben.

N MP2 HF B3LYP DFTB DFTB+vdW

2 3.5 3.3 3.6 2.9 2.9

3 4.8 5.4 4.9 3.8 3.3

5 5.1 7.6 6.3 5.2 3.3

8 2.6 8.9 6.1 5.9 1.4

Tabelle 4.1: Die relativen Stabilit¨

aten von αRund 310 (in Kcal/mole) f¨

ur verschie-

dene Kettenl¨

angen N. Die positiven Werte bezeichnen den Energieunterschied

zwischen der stabileren 310 Helix zu der αRHelix.

MP2 und DFTB+vdW beschreiben eine wesentlich schnellere Stabilisierung

der αRHelix, ein crossover ist f¨

ur N=9,10 zu erwarten. Auffallend sind die gr¨

oße-

ren Energieunterschiede in der MP2 Methode. In vorhergehenden Arbeiten wurde

aber eine generelle Tendenz des DFTB gefunden, kleinere Energiedifferenzen zwi-

schen den Konformationen vorherzusagen als die ab initio und DFT Methoden.

Dies muß nicht ein Defizit der DFTB Methode sein, es gibt Hinweise, daß die

Energiedifferenzen aufgrund des internen ’basis set superposition error’ (BSSE)

Effektes ¨

ubersch¨

atzt werden [32]. Die 310 Helix hat eine interne Wasserstoffbr¨

ucke

mehr als die αRHelix, und BSSE f¨

uhrt zu einer ¨

Ubersch¨

atzung der Bindungsener-

gie von typischerweise 1-2 kcal/mole f¨

ur Basiss¨

atze mittlerer Qualit¨

at. wie den

hier verwendeten. Daher kann man annehmen, daß die Energie der 310 Helix etwa

um diesen Betrag untersch¨

atzt wird, die Energiedifferenz also um diesen Betrag

zu groß ist. Die DFTB+VdW Werte st¨

unden dann aber in sehr guter ¨

Uberein-

stimmung mit solcherart ’BSSE-korrigierten’ ab initio Werten. Leider ist zur Zeit

keine Methode verf¨

ugbar, die den internen BSSE korrigieren kann. Man kann nur

auf Rechnungen zu kleinen Modellpeptiden zur¨

uckgreifen, die diese Diskussion zu

unterst¨

utzen scheinen [32].

Den Werten in Tabelle 4 liegen die Geometrien aus DFT Rechnungen mit

festgehaltenen Dihedralwinkeln des Peptidr¨

uckgrades zu Grunde. Mit dieser Vor-

gehensweise ist der absolute Kreuzungspunkt der relativen Stabilit¨

at nicht genau

bestimmbar, dazu sind volle Geometrieoptimierungen n¨

otig. Um dies zu bewerk-

stelligen, wurden die Geometrien f¨

ur N=4-12 mit der DFTB+VdW Methode

optimiert. F¨

ur N=4 findet noch ein ¨

Ubergang der αRin die 310 Konformation

statt. Aber schon f¨

ur N= 5 bilden sich stabile i→i+4 H-Br¨

ucken aus, also par-

tiell stabile αRHelices, wobei aber erst bei N=7-8 der Kreuzungspunkt in den

Stabilit¨

aten liegt. Im Vergleich zu der Situation ohne Dispersionswechselwirkung

bedeutet dies, daß der Kreuzungspunkt bei sehr viel kleineren Peptiden zu finden

ist, als vormals angenommen.

Dieser Wert l¨

asst sich sehr gut mit dem Wert N=6 aus der Analyse der Struk-

turen in den Proteindatenbanken vergleichen. Bis N=6 ist die 310 Konformation

stabil, f¨

ur l¨

angere Ketten findet man αR. Die Umgebungseinfl¨

usse werden die αR

Konformation zudem zus¨

atzlich beg¨

unstigen.

Zur Berechnung von Proteinstrukturen ist also eine korrekte Einbeziehung

der Dispersionskr¨

afte unumg¨

anglich. HF, aber auch DFT Rechnungen sind damit

nicht in der Lage, Proteinstrukturen i. A. vorherzusagen, es ist zu erwarten, daß

auch bei diesen Methoden αRim Protein nicht stabil sind.

Dagegen scheint die Einbeziehung der Dispersionskr¨

afte in den DFTB sehr

vielversprechend zu sein, die gegenw¨

artigen Ergebnisse lassen eine zuverl¨

assi-

ge Beschreibung von Proteinstrukturen erwarten. Da HF und DFT Methoden

z.Z. nicht uneingeschr¨

ankt anwendbar scheinen, andere semi-empirische Metho-

den wie AM1 und PM3 Proteinstrukturen sehr schlechte beschreiben, scheint

DFTB+VdW die z.Z einzige Alternative zu empirischen Kraftfeldmethoden zu

sein.

40 KAPITEL 4. STABILIT ¨

AT VON PROTEIN SEKUND ¨

ARSTRUKTUREN

Kapitel 5

Optische Eigenschaften

Die theoretische Beschreibung optischer Eigenschaften ist sehr aufwendig. Wellen-

funktions basierte Methoden sowie Propagatortechniken erzielen zwar eine sehr

gute ¨

Ubereinstimmung mit dem Experiment, sind gleichzeitig aber extrem re-

chenintensiv und k¨

onnen nur f¨

ur sehr kleine Systeme angewendet werden.

Die DFT, urspr¨

unglich nur streng g¨

ultig f¨

ur den Grundzustand, erm¨

oglicht in

ihrer Erweiterung auf zeitabh¨

angige Dichten (’time dependent’ DFT - TD-DFT)

einen robusten und einfachen Zugang zur Berechnung von Anregungsenergien im

Rahmen der ’linear response’ (LR) Theorie.

Dazu ist in einem ersten Schritt eine DFT Grundzustandsrechnung zur Be-

stimmung der Kohn-Sham Orbitale ψi(r) und Kohn-Sham Energien in¨

otig. In

einem zweiten Schritt wird die sogenannte Kopplungsmatrix Kijσ,klτ berechnet:

Kijσ,klτ =ZZ ψi(r)ψj(r) 1

|r−r0|+δ2Exc

δnσ(r)δnτ(r0)n0!ψk(r0)ψl(r0)drdr0

(5.1)

Die Differenzen von Kohn-Sham Energien ωij =i−jrepr¨

asentieren nicht die

Anregungsenergien des Systems. Um diese zu erhalten, muss die im Angeregten

Zustand modifizierte Elektron-Elektron Wechselwirkung, wie sie durch Kijσ,klτ

berechnet wird, ber¨

ucksichtigt werden. Dies geschieht durch L¨

osung des folgenden

Eigenwertproblems:

ijσ hω2

ijδikδjlδστ + 2√ωijKijσ,klτ √ωijiFI

ijσ =ω2

IFI

klτ (5.2)

Die ωIsind die Anregungsenergien des Systems, aus den Eigenvektoren Fijσ

k¨

onnen Oszillatorst¨

arken berechnet werden.

Obwohl die TD-DFT wesentlich weniger rechenzeitintensiv ist als die herk¨

omm-

lichen Wellenfunktions basierten Methoden, ist dennoch die Untersuchung von

42 KAPITEL 5. OPTISCHE EIGENSCHAFTEN

Nanomaterialien mit mehreren 100 Atomen nicht m¨

oglich. Hierzu ist ein Forma-

lismus n¨

otig, der in etwa die Rechenzeiteffizienz von semi-empirischen Methoden

hat, also zwei bis drei Gr¨

ossenordnungen schneller ist, als DFT. Deshalb wurde

der TD-DFT Formalismus auch im Rahmen des DFTB Schmas implementiert.

Die DFTB Methode stellt die Kohn-Sham Orbitale und Energien zur Verf¨

ugung,

ein Problem besteht in der effizienten Bestimmung von Kijσ,klτ .

Im Rahmen der DFTB Methode wurde jedoch schon eine N¨

aherung f¨

ur ein

Integral gefunden, das die gleiche Form wie Kijσ,klτ aufweist: Der zweite Term auf

der rechten Seite in Gl. 2.1, die zweite Ableitung der Coulomb- und Austausch-

Korrelationsenergie nach der Elektronendichte wurde durch die Funktion γαβ ap-

proximmiert.

Diese N¨

aherung kann ebenfalls zur Bestimmung von Kijσ,klτ verwendet wer-

den. Daraus resultiert ein Schema, das sehr effizient die Anregungsenergien von

Molek¨

ulen und Clustern im Rahmen der DFTB Methoden bestimmen kann [47,

48]. Damit steht mit der TD-DFTB Methode eine sehr schnelle und genaue Me-

thode zur Berechnung von Anregungsenergie biologischer Molek¨

ule zur Verf¨

ugung.

Durch Implementation von numerischen Energiegradienten sind in Zukunft Geo-

metrieoptimierungen und MD Simulationen von Molek¨

ulen im elektronisch ange-

regten Zustand m¨

oglich.

Kapitel 6

Zeitabh¨

angige Ph¨

anomene

Alle bisher diskutierten Methoden beruhen auf der Born Oppenheimer N¨

aherung.

In einigen Anwendungsf¨

allen kann jedoch die Entkopplung der Dyanmik der Elek-

tronen von der der Atomkerne nicht mehr aufrecht erhalten werden. Dies betrifft

zum Einen Stossprozesse, zum Anderen die Wechselwirkung von Materie mit elek-

tromagnetischer Strahlung aber auch strahlungslose ¨

Uberg¨

ange beispielsweise in

photochemischen Reaktionen, wie am Beispiel bR in der Einleitung beschrieben.

Gemeinsam ist allen drei Prozessen, daß im Verlauf der Dynamik der Atomkerne

elektronische ¨

Uberg¨

ange stattfinden. Um diese in die Beschreibung einbeziehen zu

k¨

onnen, ist die Ber¨

ucksichtigung der Zeitabh¨

angigkeit des quantenmechanischen

Zustandes notwendig. Hierzu kann man von einer modifizierten Lagrangefunktion

ausgehen, die ¨

uber das Variationsprinzip gekoppelte Differentialgleichungen f¨

die Elektronen und Kerne ergibt [48]. Die resultierende Schr¨

odingergleichung f¨

die Kohn-Sham-Orbitale innerhalb des DFTB-Schemas wird dann mit Hilfe eines

modifizierten Cayley-Algorithmus’ gel¨

ost.

Bisher wurde das beschriebene Verfahren an sys-anti-Isomerisierungsreaktionen

verschiedener Imine getestet und ergab sowohl bez¨

uglich des Anregungsmecha-

nismus’ als auch der Beschreibung der angeregten Potentialhyperfl¨

ache gute Re-

sultate. Da die benutzte Ankopplung des elektromagnetischen Feldes hinsichtlich

der eingestrahlten Intensit¨

at keine Annahmen macht, erscheint eine weitere An-

wendung der Methode auf Photanregungen von Schwingungsmoden und Photo-

fragmentationsprozesse aussichtsreich. So wurde mit dieser Feldankopplung das

Fulleren C60 untersucht [49]. In Zukunft sollen mit dieser Methode die photoin-

duzierte Dynamik des Retinals in dem bR Protein untersucht werden. Mit dieser

Methode k¨

onnen dann sowohl der Anregungsprozess, als auch der ¨

Ubergang vom

angeregten Zustand in den Grundzustand beschrieben werden.

44 KAPITEL 6. ZEITABH ¨

ANGIGE PH ¨

ANOMENE

Kapitel 7

Zusammenfassung und Ausblick

Ausgehend von der Enwicklung des DFTB wurden hier wesentliche Erweiterun-

gen der Methode vorgestellt, die eine Anwendung auf biologische Probleme erst

erm¨

oglicht hat. Zum Einen wurden

•Parameter f¨

ur neue Atomtypen entwickelt, die im biologischen Kontext

wichtig sind. Dies betrifft etwa Schwefel, Phosphor, Zink und weitere Me-

talle und ¨

Ubergansmetalle.

•Auch wurde die Qualit¨

at der Beschreibung der schwachen Wechselwirkun-

gen, der Wasserstoffbr¨

ucken und der VdW Wechselwirkung wesentlich ver-

bessert, was erst eine stabile Beschreibung vieler wichtiger biologischer

Strukturen erm¨

oglicht

Desweiteren wurden zur Behandlung großer Strukturen, die die M¨

oglichkeiten

der DFTB Methode ¨

ubersteigen, neue Algorithmen implementiert. DFTB wurde

•in einen QM/MM Algorithmus implementiert und damit mit einem empir-

schen Kraftfeld kombiniert

•in einen linear skalierenden (O(N)) Algorithmus implementiert

•in einen kombinierten O(N)-QM/MM implementiert, der zudem paralleli-

siert wurde, was die Anwendung auf sehr große Strukturen und lange Zeits-

kalen erlaubt.

Damit sind alle z.Z. g¨

angigen Methoden zur Verbesserung der Rechenzeiteffizienz

f¨

ur den DFTB verf¨

ugbar. Auch wurde die Methode zur Berechnung spektrosko-

pischer Daten erweitert:

•IR Spektren k¨

onnen sehr effizient ermittelt werden. Die Qualit¨

at der Ergeb-

nisse ist als semi-quantitativ zu bezeichnen, eine Skalierung der Ergebnisse

w¨

urde die ¨

Ubereinstimmung mit dem Experiment wesentlich verbessern

46 KAPITEL 7. ZUSAMMENFASSUNG UND AUSBLICK

•Die Methode wurde f¨

ur die Berechnung von optischen Anregungsspektren

erweitert.

Es wurden erfolgreiche Anwendungen auf Protontransfer- und Katalysepro-

zesse, auf die Stabilit¨

at von Proteinsekund¨

arstrukturen und die Dynamik eines

Proteins in w¨

assriger L¨

osung vorgestellt. Die Vorhersage von Geometrien und

Schwingungsfrequenzen mit Hilfe des DFTB ist excellent, in Bezug auf die Ener-

getik k¨

onnte man sich eine Verbesserung insofern vorstellen, da das z. Z. vorhan-

dene ’overbinding’ eine genaue Ermittlung von ¨

Uberganszustandsenergien nicht

erlaubt. Desweitern w¨

are die Entwicklung weiterer spektroskopischer Methoden,

wie etwa die Vorhersage des Spektrums des Raman- oder resonanten Ramanef-

fekts, sehr hilfreich.

Insgesamt ist mit der hier vorgestellten Enwicklung ein Methodenspektrum

enstanden, das eine echte Erweiterung der Methoden der Quantenchemie im

Blick auf Nanostrukturen darstellt. Die DFTB Methode ist den ab initio und

DFT Verfahren an Rechenzeitaufwand ¨

uberlegen, den quantenchemischen semi-

empirischen Methoden an Genauigkeit. Sie f¨

ullt damit eine L¨

ucke im Methoden-

spektrum, die zur Behandlung nanoskaliger Materialien dringend ben¨

otigt wird.

Literaturverzeichnis

[1] Foulkes W M C, Mitas L, Needs R J and Rajagopal G 2001 Rev. Mod. Phys.

73 33.

[2] F. Jensen, Introduction to Computational Chemistry, John Wiley and Sons,

1999

[3] R. G. Parr, W. Yang, Density Functional Theory of Atoms and Molecules,

Oxford University Press, Oxford 1989.

[4] R. M. Dreizler, E. K. U. Gross, Density Functional Theory, Springer-Verlag,

Berlin, 1990.

[5] P. Hohenberg, W. Kohn, Phys. Rev. B 136 (1964) 864.

[6] W. Kohn, L. J. Sham, Phys. Rev. A 140 (1965) 1133.

[7] M. C. Zerner, Semiempirical Molecular Orbital Methods, in: Lipkowitz,Boyd

(Hrg.), Reviews in Computational Chemistry, Vol.2, VCH Publishers, New

York 1991.

[8] M. J. S. Dewar et al., J. Am. Chem. soc. 107 (1985) 3902

[9] J. J. P. Stewart, J. Comp. Chem. 10 (1989).

[10] W. M. Foulkes, R. Haydock, Phys. Rev. B 39 (1989) 12520.

[11] W. D. Cornell et al., J. Am. Chem. Soc., 117 (1995) 5179.

[12] A. D. MacKerell et al., J. Am. Chem. Soc., 117 (1995) 11946.

[13] J. R. Maple et al., J. Am. Chem. Soc., 115 (1994) 162.

[14] W. L. Jorgensen, J. T. Rives, J. Am. Chem. Soc., 110 (1988) 1657.

48 LITERATURVERZEICHNIS

[15] Janos K. Lanyi, Bacteriorhodopsin. International Journal of Cytology vol.

187, (1999) 161-202.

Janos K. Lanyi, Molecular mechanism of ion transport in bacteriorhodop-

sin: insights from crystallographic, spectroscopic, and mutational studies. J.

Phys. Chem. B vol. 194(48), (2000) 11441-11448.

[16] M .Elstner, Ph.D. Thesis, University Paderborn: Paderborn, Germany 1998.

[17] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim,

S.Suhai, G. Seifert, Phys. Rev. B 58 (1998) 7260.

[18] M. Elstner, D. Porezag, G. Jungnickel, T. Frauenheim, S. Suhai, and G.

Seifert, in Tight-Binding Approach to Computational Materials Science (P.

E. A. Turchi, A. Gonis, and L. Colombo, eds.), vol. 491, p. 131, MRS, 1998.

[19] M. Haugk, J. Elsner, Th. Heine, Th. Frauenheim, G. Seifert, J.Comp.Mat.

Science 13, 239 (1999).

[20] D. Porezag, Th. Frauenheim, Th. K¨

ohler, G. Seifert, R. Kaschner, Phys. Rev.

B51, 12947 (1995).

[21] J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865

[22] J. Andzelm, E.Wimmer, J. Chem. Phys. 96 (1992) 1280 221.

[23] B. G. Johnson et al., J. Chem. Phys. 98 (1993) 5612

[24] T. Niehaus, M. Elstner, T. Frauenheim, and S. Suhai, J. Mol. Struc., (THEO-

CHEM) 541 (2001) 185.

[25] M. Elstner, Q. Cui, T. Frauenheim, S. Suhai, Manuskript.

[26] M. Elstner, C. Qiang, T. Frauenheim, T. Kaxiras, M. Karplus, J. Chem.

Phys., zur Publikation eingereicht.

[27] M. Elstner, T. Frauenheim, E. Kaxiras, G. Seifert, and S. Suhai, phys. stat.

sol.(b), 217/1 (2000) 375.

[28] W. Han, M. Elstner, K. Jalkanen, T. Frauenheim, S. Suhai, Int. J. Quant.

Chem., 78 ( 2000) 459

[29] Q. Cui, M. Elstner, E. Kaxiras, T. Frauenheim, M. Karplus, J. Phys. Chem.

B 105 (2001) 569

[30] M. Elstner, D. Porezag, G. Seifert, T. Frauenheim, and S. Suhai, In: Multis-

cale Modelling of Materials Symp., pp. 541–546, MRS, 1999.

LITERATURVERZEICHNIS 49

[31] M. Elstner, P. Hobza, T. Frauenheim, S. Suhai, and E. Kaxiras, J. Chem.

Phys. 114 (2001) 5149

[32] M. Elstner, K. Jalkanen, M. Knapp-Mohammady, T. Frauenheim, S. Suhai,

Chem. Phys. 263 (2001) 203

[33] H. G. Bohr, K. J. Jalkanen, M. Elstner, K. Frimand, S. Suhai, Chem. Phys.

246 (1999) 13.

[34] M. Elstner, K. Jalkanen, M. Knapp-Mohammady, T. Frauenheim, S. Suhai,

Chem. Phys. 256 (2000) 15

[35] H. Zhou, E. Taikhorshid, T. Frauenheim, S. Suhai, M. Elstner, J. Mol. Struc.

(THEOMCHEM), 277 (2002) 91.

[36] C. Alem´an, R. Roca, F.J. Luque, and M. Orozco, Proteins 28 (1997) 83.

[37] W. Han, K. Jalkanen, M. Elstner S. Suhai, J. Phys. Chem. B 102 (1998)

2587

[38] W. Yang, Phys. Rev. Lett. 66, 1438 (1991).

[39] W. Yang, T. S. Lee, J. Chem. Phys. 163 (1995) 5674.

[40] W. Yang,J. M. P´erez-Jord´a, Linear Scaling Methods for Electronic Structure

Calculations. In Encyclopedia of Computational Chemistry; P. Schleyer, Ed.;

John Wiley & Sons: New York, 1998.

[41] j S. Goedecker Rev. Mod. Phys. 71 (1999) 1085.

[42] G. Galli, Phys. Status Solidi B-Basic Res.217 (2000) 231.

[43] T. S. Lee, D. York,W. Yang, J. Chem. Phys. 105 (1996) 2744

[44] A. Warshel, M. Levitt, J. Mol. Bio. 103 (1976) 227.

[45] M. J. Field, P. A. Bash, M. Karplus, J. Comp. Chem. 11 (1990) 700.

[46] H. Y. Liu, M. Elstner, E. Kaxiras, T. Frauenheim, J. Hermans, W. Yang,

PROTEINS 44 (2001) 484.

[47] T. Niehaus, S. Suhai, F. Della Sala, P. Lugli, M. Elstner, G. Seifert, T.

Frauenheim, Phys. Rev. B 63 (2001) 085108.

[48] T. Niehaus, PhD Thesis, Universtit¨

at Paderborn 2001.

[49] B. Torralva, T. Niehaus, M. Elstner, S. Suhai, T. Frauenheim, R. E. Allen,

Phys. Rev. B 64 (2001) 153105

50 LITERATURVERZEICHNIS

Anhang A

Relevante Ver¨

offentlichungen

•T. Niehaus, M. Elstner, T. Frauenheim, S´andor Suhai, Application of an

approximate density functional method to sulfur containing compounds, J.

Mol. Struc. (THEOCHEM), 541 (2001) 185.

•M. Elstner, Th. Frauenheim, E. Kaxiras, G. Seifert, S. Suhai, A self-consistent

carge density-functional based tight-binding scheme for large biomolecules,

phys. stat. sol. (b) 217 (2000) 357.

•M. Elstner, P. Hobza, T. Frauenheim, S. Suhai, and E. Kaxiras, Hydro-

gen bonding and stacking interactions of nucleic acid base pairs: a density-

functional-theory based treatment, J. Chem. Phys. 114 (2001) 5149.

•M. Elstner, K. J. Jalkanen, M. Knapp-Mohammady, T. Frauenheim and S.

Suhai, Energetics and structure of glycine and alanine based model pepti-

des: Approximate SCC-DFTB, AM1 and PM3 methods in comparison with

DFT, HF and MP2 calculations, Chem. Phys. 263 (2001) 203.

•M. Elstner, K. Jalkanen, M. Knapp-Mohammady, T. Frauenheim, S. Suhai,

DFT studies on Helix Formation in N-Acetyle-(L-Alanyle)n-N’-Methylamide

for n = 1-20 Chem. Phys. 256 (2000) 15.

•H. Zhou, E. Tajkhorshid, T. Frauenheim, S. Suhai, M. Elstner, Applicability

of SCC-DFTB method for studying the conjugated Schiff base molecules,

J. Mol. Struc. (THEOMCHEM), 277 (2002) 91.

•W. Han, M. Elstner, K. J Jalkanen, T. Frauenheim and S. Suhai, Hybrid

SCC-DFTB/Molecular Mechanical Studies of H-Bonded Systems and of N-

Acetyl-(L-Ala)n-N’-Methylamide Helices in Water Solution Int. J. Quant.

Chem., 78 (2000) 459.

52 ANHANG A. RELEVANTE VER ¨

OFFENTLICHUNGEN

•Q. Cui, M. Elstner, T. Frauenheim, E. Kaxiras, M. Karplus, Combined

self-consistent charge density functional tight-binding (SCC-DFTB) and

CHARMM,

J. Phys. Chem. B 105 (2001) 569.

•H. Liu, M. Elstner, E. Kaxiras, T. Frauenheim, J. Hermans, W. Yang,

Quantum Mechanics Simulation of Protein Dynamics on Long Time Scale,

PROTEINS 44 (2001) 484.

•T. Niehaus, S. Suhai, F. Della-Sala, P. Luigli, M. Elstner, G. Seifert T

Frauenheim, A simple tight-binding approach to Time-Dependent Density-

Functional response-Theory, Phys. Rev. B 63 (2001) 5108

•B. Torralva, T. Niehaus, M. Elstner, S. Suhai, T. Frauenheim, and R.E.

Allen, Response of C-60 and C-n to ultrashort laser pulses, Phys. Rev. B

64 (2001) 3105.