Einfluss von Vibrationsanregungen auf die Photodynamik von Bakteriorhodopsin Influence of vibrational excitation on the photodynamics of bacteriorhodopsin (Korrigierte Version) Masterarbeit

(1)

Einfluss von Vibrationsanregungen auf die Photodynamik von Bakteriorhodopsin

Influence of vibrational excitation on the photodynamics of bacteriorhodopsin

(Korrigierte Version)

Masterarbeit

Zur Erlangung des akademischen Grades Master of Science im Studiengang „Master Chemie“

an der Heinrich-Heine-Universität Düsseldorf

vorgelegt von

Dennis Ralph Dombrowski

Angefertigt am

Institut für Theoretische Chemie und Computerchemie Unter Anleitung von

PD Dr. Oliver Weingart

(2)

Eidesstattliche Erklärung

Hiermit versichere ich, dass ich die vorliegende Masterarbeit selbstständig am Institut für Theoretische Chemie und Computerchemie der Heinrich-Heine-Universität Düsseldorf verfasst und keine anderen als die angegebenen Quellen und Hilfsmittel benutzt und die aus fremden Quellen direkt oder indirekt übernommenen Gedanken als solche kenntlich gemacht habe.

Dennis Ralph Dombrowski Düsseldorf, den

Abgabedatum: 05.11.2018 Referent: PD Dr. Oliver Weingart

Korreferent: Univ.-Prof. Dr. Birgit Strodel

Korrigierte Version vom 14.02.2019

(3)

Danksagung

Als erstes möchte ich mich bei PD Dr. Oliver Weingart für das interessante Thema und die gute Betreuung bedanken. Seine Ratschläge trugen viel zum Erfolg dieser Arbeit bei. Auch für das Korrekturlesen und bei Problemen jeglicher Art möchte ich mich bei ihm bedanken.

Zudem möchte ich mich bei Univ.-Prof. Dr. Birgit Strodel dafür bedanken, dass sie das Zweitgutachten übernimmt.

Ich möchte mich zudem noch beim Arbeitskreis für die angenehme Arbeitsatmosphäre bedanken.

Zum Schluss danke ich noch meiner Familie und meinen Freunden für die Unterstützung in der Zeit während der Masterarbeit.

(4)

Inhalt

1 Einleitung ... 1

1.1 Simulation der Photodynamik in bR ... 3

1.2 Molekülschwingungen ... 6

2 Motivation und Zielsetzung ... 7

2 Theorie ... 8

2.1 Die stationäre Schrödingergleichung und die Born-Oppenheimer-Näherung ... 8

2.2 Molekülmechanik und Kraftfeldmethoden ... 9

2.3 Das TIP3P Wassermodel ... 10

2.4 Hartree-Fock ... 10

2.5 Configuration Interaction (CI) ... 13

2.6 Multiconfiguration Self-Consistent Field (MCSCF) und CASSCF ... 13

2.7 MS-CASPT2 ... 14

2.8 Semiempirische Verfahren ... 15

2.9 OMx Modelle ... 16

2.10 Multireferenz-Konfigurationswechselwirkung ... 18

2.11 Dichtefunktionaltheorie ... 18

2.12 QM/MM Methoden ... 20

2.13 Moleküldynamik ... 21

2.14 Der Velocity Verlet Algorithmus ... 22

2.15 Der Rattle-Algorithmus ... 23

2.16 Wigner-Sampling ... 23

2.17 “Tully’s fewest switches surface hopping” Algorithmus ... 24

3 Ergebnisse ... 26

3.1 Technische Details ... 26

3.2 Das minimale QM-Modell ... 27

(5)

3.2.1 Optimierung des Grundzustands ... 27

3.2.2 Charakterisierung der elektronischen Struktur des all-trans-Retinals... 30

3.2.3 Dynamik im Schwingungsgrundzustand mit OM3/MRCI ... 31

3.2.4 Schwingungsanalyse mit OM3 und betrachtete angeregte Schwingungen ... 34

3.2.5 Dynamiken mit angeregten Schwingungsmoden mit OM3/MRCI ... 35

3.2.6 Nullpunktstrajektorie ohne kinetische Energie... 38

3.2.7 Bewertung des minimalen QM-Modells ... 39

3.3 Das vollständige QM-Modell ... 40

3.3.2 Charakterisierung der elektronischen Struktur des all-trans-Retinals... 43

3.3.4 Schwingungsanalyse mit OM3 ... 46

3.3.5 Nullpunktstrajektorie der OM3 und der B3LYP optimierten Struktur ... 47

3.3.6 Bewertung des vollständigen QM-Modells ... 48

3.4 Erweiterung des β-Iononrings ... 48

3.5 Erweiterung des QM-Bereiches am Lysin ... 49

3.6 Photodynamik des Plus1_cc Modells ... 50

3.6.2 Dynamiken mit angeregten Schwingungsmoden mit OM3/MRCI ... 53

3.6.3 Bewertung des red. QM-Modells ... 57

3.7 Die Proteinstruktur 4Y9H ... 57

3.7.2 Dynamik im Schwingungsgrundzustand ... 59

3.7.3 Nullpunktstrajektorie mit OM3/MRCI aus der OM3 und B3LYP optimierten Geometrie und mit CASSCF(10,10) aus der CASSCF(12,12) optimierten Geometrie ... 61

3.7.4 Dynamik im Schwingungsgrundzustand mit OM3/MRCI aus der B3LYP optimierten Struktur ... 63

(6)

3.7.5 Zusammenfassung der Ergebnisse des 4Y9H Modells ... 64

4 Die C6-C7 Verdrehung des Retinals in Bakteriorhodopsin... 64

5 Diskussion ... 66

6 Ausblick ... 67

7 Zusammenfassung ... 68

8 Summary ... 70

9 Literaturverzeichnis ... 72

10 Anhang... 81

10.1 IR-Wasserpektrum ... 81

10.2 Populationsgraphen und Diederwinkelanalyse mit Schwingungsanregung im minimalen QM-Modell ... 82

10.2.1 HOOP Anregung mit 3 Quanten ... 82

10.2.3 In-plane Anregung mit 3 Quanten ... 84

10.2.5 Torsionsanregung mit 3 Quanten ... 86

10.3 Nullpunktstrajektorien der erweiterten Modelle ... 88

10.4 Geometrie, IR-Spektrum und elektronische Struktur in der Grundzustandsstruktur sowie Populationsgraphen und Diederwinkelanalyse mit Schwingungsanregung im Plus1_cc Modell ... 91

(7)

10.5 Nullpunktstrajektorie im Plus1_cc Modell ... 100 10.6 Diederwinkelanalyse der Simulation im Schwingungsgrundzustand aus der B3LYP Struktur in 4Y9H ... 101 10.7 CASSCF(12,12) und CASSCF(10,10) Geometrien und elektronische Struktur ... 102

(8)

1

1 Einleitung

Der Farbstoff Retinal nimmt eine wichtige Schlüsselstellung in der belebten Natur ein. Er dient als photoaktive Komponente in Rhodopsinen (Rh), einer Familie von Rezeptorproteinen, die in der Zellmembran von Bakterien und Algen sowie in den Sehzellen von Wirbeltieren zu finden sind. [1]–[3] Die Rhodopsine in den Sehzellen der Wirbeltiere ermöglichen ihnen die visuelle Wahrnehmung der Umgebung. [3] Eine ähnliche Funktion besitzen die sensorischen Rhodopsine (SRI/SRII) in Bakterien [4] oder die Channelrhodopsine (ChR) in einigen Arten von Grünalgen. [2] Das Protein Bakteriorhodopsin (bR) findet sich in der Purpurmembran von Halobacterium salinarium, einem Bakterium, das in stark salzhaltiger Umgebung leben kann.

Mit Hilfe dieses Farbstoffs, der die bewohnten Salzsehen violett erscheinen lässt, pumpen die Bakterien lichtgesteuert Protonen vom Zellinneren zum Äußeren der Zellmembran. Der so erzeugte Protonengradient wird dazu verwendet, ADP in ATP umzuwandeln und somit Energie zu gewinnen (s. Abb. 1 links). [5] Bakteriorhodopsin ist ein Transmembranprotein, welches aus insgesamt sieben α-Helizes besteht. Der Chromophor Retinal ist über eine Schiff-Base Bindung kovalent mit einer der sieben α-Helizes verbunden. [6] In der Nähe der Schiff-Base Bindung befinden sich zwei negativ geladene Gegenionen (Asp85 und Asp212), die eine wichtige Rolle sowohl beim Protonentransfer als auch bei der Absorption des Chromophors einnehmen. [7]–

[9] Nach Lichtanregung isomerisiert der Retinalchromophor von einer all-trans zu einer 13-cis Konfiguration (Abb. 1 links). Dieser Prozess ist, mit einer Quantenausbeute von ca. 65%, nach etwa 500 fs abgeschlossen. [10]

Abbildung 1: Links: Umwandlung von ADP in ATP in Halobacterium salinarum.[11] Rechts: Initialreaktion von Retinal in bR.

(9)

2 Der mit diesem Schritt eingeleitete Photozyklus (Abb. 2 a)) verläuft über mehrere, spektroskopisch detektierbare Intermediate. [12] Im K-Intermediat liegt der Chromophor in der 13-cis Form vor. Im Schritt vom L- zum M-Intermediat wird das Proton der Schiff-Base über ein komplexes Wasserstoffbrückenbindungsnetzwerk (Abb. 2 b)) an das benachbarte Gegenion Asp85 abgegeben. Diese Protonierung sorgt schließlich dafür, dass das Proton über ein erweitertes Netzwerk am Glu194 oder Glu204 lokalisiert und anschließend an das externe Medium abgegeben. Als nächstes wird die Schiffbase über mehrere Wassermoleküle vom Asp96 reprotoniert (N- und O-Intermediat, Abb. 2a) und c)). Die Übertragung eines Protons vom intrazellularen Medium auf das Asp96 sowie vom Asp85 auf die Protonenfreisetzungstelle sind noch nicht genau verstanden. Allerdings wird vermutet, dass weitere Wassermoleküle involviert sind und dass im aktiven Zustand ein Tunnel entsteht (s.

Abb. 2 rechts), der einen Austausch zwischen Asp96 und dem intrazellularen Medium erlaubt.

[13]

Abbildung 2: a) Photozyklus von Bakteriorhodopsin. Die Zahlen neben den Intermediaten stehen für die jewilige

Absorptionswellenlänge in nm. Bild angepasst nach Ref. [14] b) Wassertoffbrückenbindungsnetzwerk in der Nähe der Schiff- Base. [15] c) Protonentransferkanal innerhalb der sieben α-Helizes, gekennzeichnet durch schwarze Pfeile.[16]

570

590

550

410 560

640 a)

b)

c)

(10)

3 Das Wasserstoffbrückennetzwerk am N-Terminus des Chromophors beinhaltet mindestens drei Wassermoleküle (WA, WB und WC in Abb. 2 c)), die auch in der Kristallstrukturanalyse aufgelöst werden. [6]

1.1 Simulation der Photodynamik in bR

Die Moleküldynamik des Retinals während der Isomerisierung wurde bereits 1998 in einer Studie von Ben-Nun et al. untersucht. [17] Sie verwendeten dabei Kraftfeldverfahren zur Beschreibung der Proteinumgebung, sowie quantenmechanisch berechnete Modellpotentiale einer Schiff-Base, die aus vier Doppelbindungen bestand. Die Torsionspotentiale und nichtadiabatischen Kopplungen von insgesamt drei Singulettzuständen (S0, S1 und S2) würden dabei mit der complete active space self-consistent field Methode (CASSCF) und dem 6-31G Basissatz ermittelt. Der angeregte S1 Zustand beschreibt dabei die HOMO-LUMO-Anregung des Chromophors, ähnlich dem 1¹Bu Zustand in Polyenen. [18] In diesem spektroskopisch hellen Zustand befindet sich der Chromophor nach der Anregung. Der darüber liegende S2

Zustand beschreibt eine HOMO²-LUMO²-Doppelanregung und befindet sich in energetischer Nähe zum HOMO-LUMO Zustand. Die Doppelanregung entspricht dabei dem symmetrieverbotenen 2¹Ag Zustand der Polyene mit geringerer Oszillatorstärke. [18] In ihrem Modell gehen Ben-Nun und Mitarbeiter davon aus, dass es entlang der Torsionskoordinate zu einer Kreuzung der S1 und S2 Zustände kommt, so dass die weitere Reaktion im dunkleren HOMO²-LUMO²-Zustand stattfindet. Ein ähnliches Modell der Isomerisierung, dass die besprochenen drei Zustände enthält, wurde bereits zuvor von Anfinrud und Mitarbeitern diskutiert und mit spektroskopischen Untersuchungen belegt (Abb. 3 rechts).

Abbildung 3: Rechts: Die drei beteiligten Potentialflächen in Abhängigkeit von der Torsion um die C13=C14 Bindung. [17]

Links: Schematische Darstellung des drei Zustandsmodel von Anfinrud und Mitarbeitern. [19]

(11)

4 Der S1 HOMO²-LUMO² Zustand kreuzt schließlich, bei einem Diederwinkel von 90°, mit dem Grundzustand. In Kombination mit den errechneten QM-Potentialen und einer Kraftfeldmethode berechneten Ben-Nun et al. die klassische Moleküldynamik der Kerne mit einem Full-Multiple-Spawning Ansatz, der Übergänge zwischen den Potentialflächen erlaubt, wobei entsprechend der berechneten Übergangswahrscheinlichkeiten jeweils neue Trajektorien generiert werden. In ihrer ersten Simulation berechneten sie eine Reaktionszeit von ca. 300 fs für das Erreichen der S0/S1 Kreuzung mit einer Quantenausbeute von 0.48. In einer erweiterten Simulation [20] kamen sie auf insgesamt 71% 13-cis Ausbeute sowie auf Reaktionszeiten von 330 fs für die S2/S1 Kreuzung und 460 fs für die S1/S0 Kreuzung. Sie beobachten dabei Rotationen entlang der C13=C14 Doppelbindung, die sowohl im Uhrzeigersinn als auch gegen den Uhrzeigersinn auftraten.

Hayashi und Schulten [21] berichteten 2003 über eine erste atomistische Simulation der Photoisomerisierung mit kombinierten quantenmechanischen und molekülmechanischen (QM/MM) Verfahren. Zur Beschreibung des Retinalchromophors verwendeten sie ein noch kleineres Modell, welches nur drei Doppelbindungen beinhaltete. Die quantenmechanischen Kräfte wurden dabei während der Simulation (engl. „on-the-fly“) mit dem CASSCF Verfahren und einem double-zeta Basissatz berechnet. Der Rest des Proteins wurde mit dem Amber Kraftfeld beschrieben. Für den Chromophor wurden lediglich zwei Zustände, d.h. der Grundzustand und der erste angeregte Zustand, berücksichtigt. Für die anfänglich 11 Trajektorien der Simulation berechneten sie die zeitabhängige Schrödingergleichung, um Übergangswahrscheinlichkeiten zum Grundzustand zu erhalten. An Punkten mit hoher Übergangswahrscheinlichkeit zum Grundzustand wurden, entsprechend den berechneten Zustandspopulationen, neue Trajektorien generiert. Mit diesen Simulationen erhalten sie eine Reaktionszeit von ca. 200 fs und eine Quantenausbeute von 0.54 für das 13-cis Photoprodukt.

Im Gegensatz zu Ben-Nun et al beobachten sie lediglich eine Rotationsrichtung für die Isomerisierung, in der sich das Proton der Schiff-Base zum Asp212 hin ausrichtet.

In einer aktuelleren Arbeit von Punwong und Martinez [22] wird für die Beschreibung des Retinalchromophors die semiempirische FOMO-CASCI (floating occupation molecular orbitals CAS configuration interaction) Methode verwendet. Diese wurde so parametrisiert, dass sie die Ergebnisse des genaueren, störungskorrigierten Multi-State CAS-Verfahrens (MS-CASPT2) wiedergibt. Die mit Multiple-Spawning berechneten Simulationen ergeben eine

(12)

5 unidirektionale Torsion entlang der C13=C14 Bindung mit einer Quantenausbeute von 0.7 und einer Lebensdauer von 620 fs.

In seiner Promotionsarbeit beschreibt B. Demoulin [23] einen weiteren wichtigen Aspekt, der die Lage der Absorptionsbande in bR betrifft. MM/CASSCF//CASPT2-Berechnungen, die die Standard-Protonierung der umgebenden Aminosäuren berücksichtigen, überschätzen diesen Wert. Altoè et al. [24] z. B. berechnen die maximale Absorption zu 534 nm (gegenüber 570 nm im Experiment), wobei sie nur die Geometrie des Retinalchromophors optimieren. Der Rest der Proteinumgebung blieb bei den Kristallstrukturdaten fixiert. Bei Optimierung des komplexen Gegenions (Abb 2 b)) sinkt in Demoulins Berechnung der Wert der Absorption auf ca. 480 nm. Eine Protonierung von Asp212 führt in diesem Fall zu einer deutlichen Verbesserung der Anregungsenergie, sie liegt dann zwischen 560 und 580 nm. P. Tuchel [25]

erreichte während seiner Masterarbeit einen ähnlichen Effekt durch Verschiebung eines Protons von Arg82 zu Asp212. Einzeldynamikrechnungen mit CASSCF(12,12), gemittelt über zwei Zustände, ergeben bei Demoulin eine Reaktionszeit von ca. 200 fs für die Isomerisierung bei protoniertem Asp212. Bei dieser Reaktion sind zwei Koordinaten maßgeblich beteiligt:

Nach der Photoanregung in den HOMO-LUMO Zustand invertieren die Bindungslängen des Chromophors, d.h. Doppelbindungen werden zu Einfachbindungen und umgekehrt. Erst danach bewegt sich das Molekül entlang der Torsionskoordinate. Die Erkenntnisse werden durch Reaktionspfadberechnungen von Altoè [24] gestützt. [26]

Eine aktuelle Studie [27] verfolgt die Isomerisierung des Retinalchromophors in bR mit kristallographischen Methoden. Dabei wurden Femtosekunden-Röntgenlaser verwendet. In Kombination mit QM/MM Berechnungen wird die Zeit, die der Chromophor bis zum Erreichen einer Torsion von 90° benötigt, mit Werten >646 fs angegeben.

Neben der Umkehr der Bindungsalternanz und der Torsion entlang der C13=C14 Bindung spielen HOOP Moden an der isomerisierenden Doppelbindung eine entscheidende Rolle für die Entstehung des Photoprodukts in Schiff-Basen. Dies wird in mehreren Studien der Photodynamik von Retinalen belegt [28], [29]. Die Richtung der HOOP-Torsion beim Übergang zum Grundzustand entscheidet dabei, ob das jeweilige cis- oder aber das trans-Isomer entsteht. Die Bewegung der Kohlenstoffatome folgt dabei der Bewegung der beteiligten Wasserstoffe, d.h. für eine erfolgreiche trans-cis Isomerisierung ändert sich der Diederwinkel H-C=C-H beim Übergangspunkt in Richtung des cis-Produktes.

(13)

6 1.2 Molekülschwingungen

Jedes mehratomige Molekül ist in der Lage, auf bestimmten gequantelten Vibrationsniveaus zu schwingen. Bei Raumtemperatur befinden sie sich normalerweise im Schwingungsgrundzustand. Um diese Schwingungen mithilfe von Infrarot-Strahlung anzuregen, muss sich während der Vibration das Dipolmoment des Moleküls ändern. Dieses Verhalten kann in der Spektroskopie benutzt werden, um über spezifische Schwingungen die Struktur eines Moleküls aufzuklären und zu verstehen. Ergänzend kann Raman-Spektroskopie verwendet werden, wenn das Molekül seine Polarisierbarkeit während einer Schwingung ändert. So haben Smith et al. [30] mithilfe eines Ramanspektrums erfolgreich die Kohlenstoffstreckschwingung für Doppelbindungen bei 1500 bis 1660 cm^-1 und für eine Einfachbindung bei 1100 bis 1250 cm^-1 in Bakteriorhodopsin zuweisen können. Auch die Schwingung, die die Wasserstoffatome aus der Ebene des Kohlenstoffgerüstes auslenkt (HOOP, eng. hydrogen out-of-plane), konnte von ihnen zwischen 700 und 1000 cm^-1 beobachtet werden.

Abbildung 4: [7-(diethylamino)coumarin-4-yl]methyl-azide. [31]

Schwingungsmoden können aber auch über Anregung gezielt beeinflusst werden. So können mit sog. „chirped“ IR-Pulsen nicht nur einzelne Schwingungsmoden angeregt werden, sondern es lässt sich auch Einfluss auf die Phase der Schwingungsanregung nehmen. [32] Somit erscheint es möglich, z. B. eine spezifische Mode in einer vorbestimmten Richtung anzuregen, und damit nach Anregung mit einem Laserpuls Einfluss auf die Photodynamik einfacher Systeme zu nehmen. Experimentell wurde bereits von Shin et al. [33] nachgewiesen, dass Vibrationsanregungen die Reaktion von Dimethylsulfoxid und Methanol unter Anregung von C-H und O-H Schwingungen gezielt beeinflussen können. Des Weiteren konnten Kern-Michler et al. [31] zeigen, dass sich auch die Photochemie mit Vibrationsschwingungen beeinflussen

(14)

7 lässt. Dazu haben sie zwei DEACM-N3 ([7-(diethylamino)coumarin-4-yl]methyl-azide, s.

Abbildung 4) Isotopomere mit einer Viper (Vibrationally Promoted Electronic Resonance) Sequenz angeregt. Als Resultat konnten sie selektiv eines der Isotopomere anregen und zur Photoreaktion bringen.

2 Motivation und Zielsetzung

Obwohl die grundlegende Photochemie von Bakteriorhodopsin als bekannt vorausgesetzt wird, ergeben sich nach wie vor Unstimmigkeiten in der aktuellen Literatur. Dies betrifft vor Allem die Reaktionszeiten in Experiment und Simulation sowie die spektroskopischen Eigenschaften von Retinal in bR. Aus diesem Grund sollten moleküldynamische Berechnungen der Retinalphotoisimerisierung mit einer genügend großen Anzahl von Trajektorien durchgeführt werden, die eine statistische Analyse der Reaktion und damit das Erhalten von Reaktionszeiten und Quantenausbeuten erlaubt. Dabei sollte die semiempirische OMx/MRCI Strategie, die sich bereits bei der Beschreibung der Retinalisomerisierung in der Channelrhodopsin-Chimäre C1C2 als erfolgreich erwies, [34] verwendet werden.

Ausgangspunkt für die Untersuchungen ist ein reduziertes QM/MM Modell, welches in der Masterarbeit von Philipp Tuchel [25] erstellt wurde. Dieses betrachtet nur die Polyenkette des Retinals. Nach Erweiterung auf das vollständige System und Validierung sollte mit diesem Modell der Einfluss von Schwingungsanregungen auf die Photodynamik des Retinals in bR untersucht werden. Die Untersuchung spezifischer Schwingungsanregungen wurde vor allem durch eine neue experimentelle Methode, die zurzeit in der Arbeitsgruppe von Dwayne Miller (DESY, Hamburg) entwickelt wird, motiviert. Die Methode sieht vor, mindestens zwei Schwingungsfreiheitsgrade unabhängig voneinander anzuregen, und dann das Molekül mit einem UV-Puls in den angeregten Zustand zu verbringen. Dabei kann bei den individuellen Moden auch die Phase, d.h. die Richtung der Schwingung, vorgegeben werden.

Für die Untersuchungen von Schwingungsanregungen in bR sollte der Einfluss der Hydrogen out-of-plane (HOOP) Mode, der Torsionsmode um die C13=C14 Bindung und der Kohlenstoffstreckschwingungen entlang des Kohlenstoffgerüstes berücksichtigt werden.

(15)

8

2 Theorie

2.1 Die stationäre Schrödingergleichung und die Born-Oppenheimer-Näherung

Um das Verhalten von kleinen Teilchen wie Elektronen und Atomkernen mit quantenmechanischen Verfahren beschreiben zu können, kann die sogenannte stationäre Schrödingergleichung gelöst werden.

ℋ𝛹 = 𝐸𝛹 1

Dabei ist ℋ der Hamiltonoperator, 𝛹 die Wellenfunktion und E die entsprechende Energie des Systems.

Diese Gleichung ist für ein System mit mehreren Teilchen schwer zu lösen, weshalb sie durch die sogenannte Born-Oppenheimer-Näherung vereinfacht werden kann. Dabei geht man davon aus, dass sich die Atomkerne deutlich langsamer aufgrund ihrer Masse sind als die Elektronen. Dies hat zur Folge, dass sich die Kernbewegungen von den Elektronenbewegungen abseparieren lassen und die Elektronenbewegung somit im fixierten Feld der Kerne stattfindet. Der elektronische Hamiltonoperator ℋ_{𝑒𝑙𝑒𝑐} kann dann durch folgende Gleichung beschrieben werden:

ℋ_{𝑒𝑙𝑒𝑐} = − ∑1 2∇_𝑖²

𝑁

𝑖=1

− ∑ ∑𝑍_𝐴 𝑟_𝑖𝐴

𝑀

𝐴=1 𝑁

𝑖=1

+ ∑ ∑ 1 𝑟_𝑖𝑗

𝑁

𝑗>𝑖 𝑁

𝑖=1

2

∇² ist dabei der Laplaceoperator, der auf das i-te Elektron wirkt, ZA ist die Kernladung, riA ist der Abstand zwischen Elektron i und Atomkern A, und rij ist der Abstand zwischen den Elektronen i und j.

Damit lässt sich das elektronische Lösung der Schrödingergleichung wie folgt formulieren:

ℋ_{𝑒𝑙𝑒𝑐}𝛷_{𝑒𝑙𝑒𝑐} = 𝐸_{𝑒𝑙𝑒𝑐}𝛷_{𝑒𝑙𝑒𝑐} 3

Analog dazu kann bei der Betrachtung der Kerne die Elektronenwechselwirkung auf diese als Mittelwert angenommen werden:

ℋ_{𝑛𝑢𝑐𝑙}= − ∑ 1 2𝑀_𝐴

𝑀

𝐴=1

∇_𝐴² + 〈− ∑1 2∇_𝑖²

𝑁

𝑖=1

− ∑ ∑𝑍_𝐴 𝑟_𝑖𝐴

𝑀

𝐴=1 𝑁

𝑖=1

+ ∑ ∑ 1 𝑟_𝑖𝑗

𝑁

𝑗>𝑖 𝑁

𝑖=1

〉 + ∑ ∑𝑍_𝐴𝑍_𝐵 𝑅_𝐴𝐵

𝑀

𝐵>𝐴 𝑀

𝐴=1

4

(16)

9 Dabei ist MA die Masse des Kerns A.

Die letzten beiden Terme in Gleichung 4 beschreiben eine Potentialfläche, auf der sich die Kerne innerhalb der Born-Oppenheimer-Näherung bewegen.

In der Nähe von konischen Durchdringungen, d.h. wenn zwei Potentialflächen mit der gleichen Symmetrie sich einander annähern, bricht die Born-Oppenheimer-Näherung zusammen. Die Elektronenbewegung ist in diesen Bereichen nicht mehr unabhängig von der Kernbewegung und die Kopplung zwischen den Beiden muss berücksichtigt werden, indem Korrekturterme hinzugezogen werden müssen. Eine Möglichkeit der Korrektur ist der Born-Huang Ansatz [35].

Dabei wird die Gesamtwellenfunktion aus einer Linearkombination von Kernwellenfunktionen und elektronischen Wellenfunktionen gebildet, um eine nichtadiabatische Beschreibung mit zu berücksichtigen. [36]–[38]

2.2 Molekülmechanik und Kraftfeldmethoden

Kraftfeldmethoden ermöglichen die Beschreibung großer molekularer Systeme. Die verwendeten Energiefunktionen setzen sich additiv aus einfachen physikalischen Potentialen zusammen. Dabei werden nur die Kerne, jedoch nicht die Elektronen berücksichtigt. Die Atome eines Moleküls werden als geladene Massezentren beschrieben, die über harmonische oder periodische Potentiale mit anderen Atomen verbunden sind. Nichtbindende Wechselwirkungen werden über das Coulombgesetz bzw. über Lennard-Jones Potentiale berechnet. Das Amber Kraftfeld, dass in dieser Arbeit verwendet wurde, besitzt die Gestalt:

𝑉(𝑟^𝑁) = ∑ 𝑘_𝑏(𝑙 − 𝑙₀)²

𝐵𝑖𝑛𝑑𝑢𝑛𝑔𝑒𝑛

+ ∑ 𝑘_𝑎(𝜃 − 𝜃₀)²

𝑊𝑖𝑛𝑘𝑒𝑙

+ ∑ ∑1

2𝑉_𝑛[1 + cos(𝑛𝜔 − 𝛾)]

𝑛 𝑇𝑜𝑟𝑠𝑖𝑜𝑛𝑒𝑛

+ ∑ ∑ 𝑓_𝑖𝑗{𝜖_𝑖𝑗[(𝑟_0𝑖𝑗 𝑟_𝑖𝑗)

12

− 2 (𝑟_0𝑖𝑗 𝑟_𝑖𝑗 )

6

] + 𝑞_𝑖𝑞_𝑗 4𝜋𝜖₀𝑟_𝑖𝑗}

𝑁

𝑖=𝑗+1 𝑁−1

𝑗=1

5

Dabei stehen die ersten drei Terme, unterteilt in Bindungen, Winkeländerungen und

(17)

10 Torsionen, für bindende Wechselwirkungen. Der letzte Term fasst alle nichtbindenden Wechselwirkungen zusammen. [36], [39], [40]

Das Amber Kraftfeld wurde bereits in den frühen 1980er Jahren entwickelt. In der Arbeit wurden zwei Varianten des Kraftfelds (ff99SB [41] und ff14SB [42]) eingesetzt, die sich geringfügig in der Parametrisierung unterscheiden. Die neuere ff14SB Version zeichnet sich durch verbesserte Parameter für das Proteinrückgrat aus, während bei ff99SB für die Parametrisierung nur relativ kleine organische Moleküle verwendet wurden. Die benötigten Kernladungen 𝑞_𝑖 in Formel 5 werden dabei aus einem Fit von partiellen Atomladungen an das quantenchemisch berechnete elektrostatische Potential mit der MP2//cc-pVTZ//HF/6-31G*

Methode erhalten. [43], [44]

2.3 Das TIP3P Wassermodel

Das TIP3P (Transferable Intermolecular Potential 3 Point) [45] ist ein starres Modell und behandelt Wasser als Teilchen mit drei möglichen Bindungstellen, die in Wechselwirkungen mit anderen Teilchen treten können. Dabei tragen die beiden Wasserstoffatome eine Ladung von +0.417, während das Sauerstoffatom eine Ladung von -0.834 aufweist. Nur dieses Atom verfügt über ein Lennard-Jones Potential, dessen Parameter so gewählt wurden, dass die Dichte von Wasser bei Raumtemperatur und Normaldruck sowie die Verdampfungsenthalpie bei Raumtemperatur gut wiedergegeben werden. Die Geometrie aus experimentellen Beobachtung wurde aus der Gasphase entnommen. [46] Innerhalb einer Moleküldynamik- simulation muss das TIP3P Modell starr gehalten werden, da seine Parameter keine physikalisch sinnvollen Schwingungsfrequenzen ergeben. Dies geschieht z. B. mit dem Rattle- Algorithmus. [47] (s. Kapitel 2.15)

2.4 Hartree-Fock

Das Hartree-Fock Verfahren ist eine Methode zur Berechnung der zeitunabhängigen Schrödingergleichung für Systeme mit mehreren Elektronen. Sie beruht auf dem sog.

Variationsprinzip, welches besagt, dass die Energie einer genäherten Wellenfunktion immer über der Energie der exakten Wellenfunktion liegen muss. Zur Beschreibung eines Mehrelektronensystems wird eine Slaterdeterminante verwendet. Diese hat, wie die Wellenfunktion selbst, einen antisymmetrischen Charakter bezüglich der Vertauschung zweier Elektronen. Da die meisten Systeme nur über doppelbesetzte Molekülorbitale

(18)

11 verfügen, kann in diesen Fällen die dazugehörige Slaterdeterminante vereinfacht werden. Für ein System mit N Elektronen in N/2 Orbitalen besitzt sie dann die folgende Form:

𝛹(𝑥₁, 𝑥₂, … , 𝑥_𝑛)

= (𝑁!)^−1/2

|

𝜓₁(𝑥₁)𝛼(1)𝜓₁(𝑥₁)𝛽(1) … 𝜓^𝑛

2(𝑥₁)𝛼(1)𝜓^𝑛

2(𝑥₁)𝛽(1) 𝜓₁(𝑥₂)𝛼(2)𝜓₁(𝑥₂)𝛽(2) … 𝜓^𝑛

2(𝑥₂)𝛼(2)𝜓^𝑛

2(𝑥₂)𝛽(2)

… 𝜓₁(𝑥_𝑛)𝛼(𝑛)𝜓₁(𝑥_𝑛)𝛽(𝑛) … 𝜓… ^𝑛

2(𝑥_𝑛)𝛼(𝑛)𝜓^𝑛

2(𝑥_𝑛)𝛽(𝑛)

|

6

Dabei ist (𝑁!)^−1/2 der Normierungsfaktor und 𝜓_𝑖(𝑥_𝑁) das i-te doppelbesetzte Spinorbital in dem sich das N-te Elektron 𝑥_𝑁 befindet.

Die Elektronenabstoßung wird in dieser Methode gemittelt behandelt. Dabei wird das Mehrelektronenproblem in Einelektronenprobleme unterteilt, die sich deutlich leichter lösen lassen. Die Einelektronengleichung für ein einzelnes Spinorbital in der Hartree Fock-Näherung lässt sich durch folgende Gleichung beschreiben:

𝑓(𝑖)𝜓_𝑖(𝑥_𝑖) = 𝜀𝜓_𝑖(𝑥_𝑖) 7

Dabei ist ε die Einelektronenenergie, 𝜒(𝑥_𝑖) das Spinorbital mit dem Elektron 𝑥_𝑖 und f(i) der Fockoperator, der die Wechselwirkung eines einzelnen Elektrons im gemittelten elektrostatischen Feld aller anderen Elektronen verkörpert. Der Fock-Operator setzt sich wie folgt zusammen:

𝑓(𝑖) = ℎ(𝑖) + ∑ 2𝐽_𝑎(𝑖) − 𝐾_𝑎(𝑖)

𝑁/2

𝑎

8

Dabei ist h(i) der Kernhamiltonoperator, der die kinetische und die potentielle Energie eines Elektrons im Feld der Kerne angibt. Jb(i) ist der Coulomboperator und Kb(i) ist der Austauschoperator des i-ten Elektrons. Coulomb- und Austauschoperatoren vermitteln dabei die Wechselwirkung des i-ten Elektrons mit dem gemittelten Feld aller anderen j Elektronen, gekennzeichnet durch das Integral über 𝑥_𝑗:

(19)

12 𝐽_𝑏(𝑖)𝜓_𝑎(𝑖) = [∫ 𝑑𝜓_𝑗𝜓_𝑏^∗(𝑗) 𝑟₁₂⁻¹𝜓_𝑏(𝑗)] 𝜓_𝑎(𝑖) 9

𝐾_𝑏(𝑖)𝜓_𝑎(𝑖) = [∫ 𝑑𝜓_𝑗𝜓_𝑏^∗(𝑗) 𝑟₁₂⁻¹𝜓_𝑎(𝑗)] 𝜓_𝑏(𝑖) 10

Beim Austauschoperator sind die Orbitalindizes a und b vertauscht.

Für einfache Systeme lassen sich die HF-Gleichungen numerisch lösen. Für die Beschreibung der Molekülorbitale 𝜒_𝑖 verwendet man Linearkombinationen von Basisfunktionen. Im sog.

LCAO-Verfahren (engl. Linear combination of atomic orbitals) werden die Molekülorbitale in einen Satz atomzentrierter Basisfunktionen entwickelt:

𝜓_𝑖 = ∑ 𝑐_𝜇𝑖𝜙_𝜇

𝑁

𝜇=1

11

Setzt man diese Beziehung in die HF-Gleichungen ein und integriert über die Koordinaten der Elektronen, so gelangt man zu den Roothaan-Hall-Gleichungen

∑ 𝐹_𝜇𝜈𝑐_𝜈𝑖 =

𝑁

𝜈=1

𝜖_𝑖∑ 𝑆_𝜇𝜈𝑐_𝜈𝑖

𝑁

𝜈=1

12

Diese lassen sich vereinfacht auch in Matrixform schreiben:

𝐹𝐶 = 𝑆𝐶𝜖 13

Dabei ist F die Matrix der Fock-Operatoren, C die Matrix der Orbitalkoeffizienten, S ist die Überlappungsmatrix und 𝜖 die Diagonalmatrix der Orbitalenergien. Im HF-Verfahren sucht man dabei die optimalen Koeffizienten C, die die Energie des Systems minimieren. Da die entsprechenden Orbitale von den anderen Elektronen abhängig sind, muss Gleichung 13 mithilfe des Variationsprinzips iterativ gelöst werden. Die Verwendung einer vollständigen Basis führt zum sogenannten Hartree-Fock Limit. In diesem Limit ist die beste Wellenfunktion erreicht, die mit nur einer Slater-Determinanten beschrieben werden kann. [36], [48]

(20)

13 2.5 Configuration Interaction (CI)

Im HF-Verfahren wird die Elektronenkorrelation durch die gemittelte Beschreibung vernachlässigt. Eine Verbesserung lässt sich über die Beschreibung der Wellenfunktionen mit mehreren Determinanten erreichen. Die Konfigurationswechselwirkungsmethode (CI, eng.

Configuration interaction) berücksichtigt neben der Grundzustandsdeterminante auch angeregte Konfigurationen, die aus dieser durch Tausch von besetzten Orbitalen mit virtuellen Orbitalen hervorgehen:

𝛹_𝐶𝐼 = 𝑎₀𝛷_𝐻𝐹+ ∑ 𝑎_𝑆𝛷_𝑆

𝑆

+ ∑ 𝑎_𝐷𝛷_𝐷

𝐷

+ ∑ 𝑎_𝑇𝛷_𝑇

𝑇

+ ⋯ = ∑ 𝑎_𝑖𝛷_𝑖

𝑖=0

14

Dabei beschreibt 𝛷_𝐻𝐹 die Slater-Determinante aus der Hartree-Fock Rechnung und 𝛷_𝑖 die i- fach angeregten Slater-Determinanten, in denen eine bestimmte Anzahl virtueller und besetzter Orbitale vertauscht sind. S steht in diesem Fall für einfach angeregte, D für doppelt angeregte, und T für dreifach angeregte Determinanten in Bezug auf die Hartree-Fock Konfiguration. Bei Berücksichtigung aller möglichen (d.h. bis zu N-fach) Anregungen spricht man von einer „Full CI“ gesprochen. Das Ergebnis einer „Full CI” Rechnung beinhaltet die exakte Elektronenkorrelation innerhalb der verwendeten Basis.

Da eine „Full CI“ Rechnung nur für die kleinsten Systeme möglich ist, beschränkt man sich oft nur auf bestimmte Arten von Anregungen. Eine CIS Rechnung würde nur Einfachanregungen mitbetrachten, während eine CISDT auch Doppel- und Dreifachanregungen berücksichtigen würde. [36]

2.6 Multiconfiguration Self-Consistent Field (MCSCF) und CASSCF

Die Multikonfigurations-SCF Methode optimiert in einem iterativen Verfahren zusätzlich zu den Koeffizienten 𝑎_𝑖 der Determinanten auch die Orbitalkoeffizienten 𝑐_𝑖 und liefert damit, insbesondere für angeregte Zustände, bessere Ergebnisse. Die Optimierung der Koeffizienten erfordert dabei deutlich mehr Computerressourcen. Für die Beschreibung der Wellenfunktion in diesem Ansatz werden häufig Konfigurationszustandsfunktionen (engl. Configuration State Functions, CSFs) verwendet. Diese Funktionen bestehen aus mehreren, linearkombinierten Slater-Determinanten und einer Funktion, die die korrekte Spinsymmetrie gewährleistet. Um auch größere Moleküle untersuchen zu können, verringert man in vielen Ansätzen die Anzahl der betrachteten Elektronen und Anregungen. In der Complete Active Space Self-Consistent

(21)

14 Field (CASSCF) Methode unterteilt man den Orbitalraum in aktive und inaktive Orbitale (Abbildung 5). [36] Im aktiven „Fenster“ werden dabei alle Möglichkeiten der Anregung, d.h.

alle Verteilungen von Elektronen auf Orbitale, berücksichtigt. Die inaktiven Orbitale ändern ihre Besetzung während der Rechnung nicht.

Abbildung 5: Schematische Darstellung eines aktiven Fensters in CASSCF. Bearbeitet nach Ref. [36]

Die Auswahl der aktiven Orbitale beeinflusst dabei maßgeblich das Ergebnis. In der Regel wählt man diejenigen Orbitale aus, die für die betrachtete Reaktion von chemischer Bedeutung sind. Will man z.B. einen Bindungsbruch beschreiben, so muss das aktive Fenster bindende und antibindende Orbitale beinhalten, die diese Bindung charakterisieren. Für die Beschreibung konjugierter π-Systeme wie dem Retinalchromophor würde man entsprechend alle Orbitale mit π-Symmetrie einbeziehen.

Soll die Energie mehrere Zustände gleichwertig beschrieben werden, verwendet man das state-averaging Verfahren. Dabei wird die Energie der Wellenfunktion bezüglich des Mittelwerts der Energie der beteiligten Zustände optimiert:

𝐸_{𝑎𝑣𝑒𝑟} = ∑ 𝜔_𝑖𝐸_𝑖

𝑖

, ∑ 𝜔_𝑖 = 1

𝑖

15

𝐸_𝑖 ist dabei die Energie des Zustands i und 𝜔_𝑖 der Wichtungsfaktor für den entsprechenden Zustand. [36]

2.7 MS-CASPT2

Das CASSCF-Verfahren bietet eine bessere Beschreibung angeregter Zustände, jedoch treten nur die Orbitale im aktiven Fenster miteinander in Wechselwirkung. Vernachlässigt wird bei

(22)

15 CASSCF vor allem die dynamische Korrelation. Dies führt häufig zu schlechten absoluten Werten für Anregungsenergien. Auch die Reihenfolge der Zustände kann davon betroffen sein. Die Störungstheorie bietet einen Ansatz, um die Elektronenkorrelation besser zu beschreiben. Diese Verfahren unterteilen den Hamiltonoperator in einen Term nullter Ordnung sowie in einen Störterm, der die Korrelationseffekte beschreibt. Im CASPT2 Verfahren wird eine CASSCF Wellenfunktion verwendet um mithilfe der Störung zweiter Ordnung die dynamische Elektronenkorrelation zu berücksichtigen. Diese Methode versagt jedoch bei der Beschreibung von Potentialbereichen, in denen eine starke Mischung verschiedener Zustände auftritt. Um dieses Problem zu beheben, wird ein mehrdimensionaler Referenzraum verwendet, der durch zwei oder mehr state-average CASSCF Zustände beschrieben wird. [49]

2.8 Semiempirische Verfahren

Semiempirische Methoden verwenden Integralnäherungen, um die Anzahl der zu berechnenden Terme zu reduzieren. Der dabei auftretende Fehler wird durch eine Parametrisierung mit empirischen Daten verringert. Die zentrale Näherung, die bei semiempirischen Verfahren verwendet wird, ist die ZDO Näherung (engl. zero differential overlap). Dabei werden alle Produkte von Funktionen, die an unterschiedlichen Zentren lokalisiert sind und sich auf dieselbe Elektronenkoordinate beziehen, gleich null gesetzt. Das sorgt dafür, dass die Überlappungsmatrix zur Einheitsmatrix wird, alle Einelektronen- Dreizentrenintegrale null werden und alle Zweielektronen-Dreizentren- und -Vierzentrenintegrale vernachlässigt werden. Damit ergibt sich die Energie aus der Diagonalisierung der Gleichung:

𝐹𝐶 = 𝐶𝜖 16

Als Integralnäherung kann unter anderem CNDO (engl. complete neglect of differential overlap) [50], INDO (engl. intermediate neglect of differential overlap) [51] und NDDO (engl.

neglect of diatomic differential overlap) [51] verwendet werden. Diese Näherungen unterscheiden sich in der Anzahl der Überlappungsterme, die vernachlässigt werden. NDDO unterscheidet sich von den anderen beiden Näherungen unter anderem darin, dass durch eine bessere Behandlung der Zweizentren-Wechselwirkungen ein gewisser Grad an Anisotropie in

(23)

16 diesen Wechselwirkungen erhalten bleibt. Einer der bekannteren Methoden, die die NDDO- Näherung verwenden, ist die MNDO-Methode. [52] Aus der MNDO-Methode sind die Orthogonalisierungsmodelle OM1, OM2 und das in dieser Arbeit hauptsächlich verwendete OM3 [53], [54] hervorgegangen. [36], [52]

2.9 OMx Modelle

Die OMx Modelle [53] (x = 1, 2, 3) sind semiempirische Verfahren, die auf HF-SCF Molekülorbitalen basieren und die Valenzelektronen mit einer minimalen Basis behandeln.

Die inneren Elektronen werden dabei durch eine reduzierte Kernladung und ein effektives Kernpotential dargestellt. Die dynamischen Elektronenkorrelationseffekte werden gemittelt durch eine bestimmte Auswahl an Zweielektronenintegralen und passender Parametrisierung dargestellt. Die Drei- und Vierzentren-Zweielektronenintegrale werden vernachlässigt. In den OMx Methoden wird ein zentrales Problem der semiempirischen Verfahren korrigiert. Durch die Vernachlässigung der Überlappmatrix S werden, in einem homonuklearen Molekül, die HOMO bzw. LUMO-Orbitale im Vergleich zu den Atomorbitalen gleichermaßen stabilisiert bzw. destabilisiert (s. Abbildung 6). Im HF-Verfahren jedoch ist die Aufspaltung unsymmetrisch. Versucht man, eine semiempirische Methode auf die Eigenschaften eines angeregten Zustands zu parametrisieren, wird so zugleich die Beschreibung des Grundzustands ungenauer.

Abbildung 6: Aufspaltung des HOMO und LUMO Orbitals mit (links) und ohne Vernachlässigung der Überlappmatrix S (rechts). [55]

Die OMx-Methoden führen eine Orthogonalisierungskorrektur ein, die wiederum eine antisymetrische Aufspaltung erlaubt. Die Orthogonalisierungskorrekturen betreffen dabei vor allem die Einelektronenintegrale. [53]

Die drei Methoden OM1, OM2 und OM3 unterscheiden sich in der Anzahl der Terme, auf die die Orthogonalisierungskorrektur angewandt wird. In OM1 wird die Korrektur nur für den Ein-

(24)

17 zentren Anteil des Kernhamiltonoperators verwendet, während sie im OM2 auch in den Zwei- zentren Termen bzw. Resonanzintegralen verwendet wird.

Der Einelektronhamiltonoperator kann in diesen Methoden wie folgt ausgedrückt:

ℋ_𝜇𝜈 = 𝑈_𝜇𝜇𝛿_𝜇𝜈+ ∑[𝑉_{𝜇𝜈,𝐵}^𝑠 + 𝑉_{𝜇𝜈,𝐵}^𝑂𝑅𝑇 + 𝑉_{𝜇𝜈,𝐵}^𝑃𝐼 + 𝑉_{𝜇𝜈,𝐵}^𝐸𝐶𝑃]

𝐵

17 ℋ_𝜇𝜆 = 𝛽_𝜇𝜆+ ∑ 𝑉_{𝜇𝜆,𝐶}^𝑂𝑅𝑇

𝐶

18

Dabei sind μ und υ Orbitale am selben Atom A, λ ein Orbital am Atom B und C entspricht dem dritten Atom. 𝑈_𝜇𝜇 sind anpassbare Einelektronenergien, 𝑉_{𝜇𝜈,𝐵}^𝑠 ist die Kern-Elektron Anziehung, 𝑉_{𝜇𝜈,𝐵}^𝑂𝑅𝑇 ist die Orthogonalisierungskorrektur, 𝑉_{𝜇𝜈,𝐵}^𝑃𝐼 ist das Durchdringungsintegral, dass die Unterschiede zwischen der Kern-Elektron Anziehung und der Elektron-Elektron Abstoßung wiedergeben, 𝑉_{𝜇𝜈,𝐵}^𝐸𝐶𝑃 ist die Kraft, die von den inneren Kernelektronen ausgeht und 𝛽_𝜇𝜆 ist das Resonanzintegral. Der Ausdruck für die Zwei-Zentren (𝑉_{𝜇𝜈,𝐵}^𝑂𝑅𝑇) und Drei-Zentren Orthogonalisierungskorrektur ergibt sich dann wie folgt:

𝑉_{𝜇𝜈,𝐵}^𝑂𝑅𝑇 = −1

2𝐹₁^𝐴∑(𝑆_𝜇𝜌𝛽_𝜌𝜈+ 𝛽_𝜇𝜌𝑆_𝜌𝜈)

𝜌∈𝐵

+1

8𝐹₂^𝐴∑ 𝑆_𝜇𝜌𝑆_𝜌𝜈(ℋ_{𝜇𝜇,𝐵}^𝑙𝑜𝑐 + ℋ_{𝜈𝜈,𝐵}^𝑙𝑜𝑐 − 2ℋ_{𝜌𝜌,𝐴}^𝑙𝑜𝑐 )

𝜌∈𝐵

19

𝑉_{𝜇𝜈,𝐶}^𝑂𝑅𝑇 = −1

2𝐺₁^𝐴𝐵∑(𝑆_𝜇𝜌𝛽_𝜌𝜆+ 𝛽_𝜇𝜌𝑆_𝜌𝜆)

𝜌∈𝐶

+1

8𝐺₂^𝐴𝐵∑ 𝑆_𝜇𝜌𝑆_𝜌𝜆(ℋ_{𝜇𝜇,𝐶}^𝑙𝑜𝑐 + ℋ_{𝜆𝜆,𝐶}^𝑙𝑜𝑐 − ℋ_{𝜌𝜌,𝐴}^𝑙𝑜𝑐 − ℋ_{𝜌𝜌,𝐵}^𝑙𝑜𝑐 )

𝜌∈𝐶

20

Dabei geben 𝐹₁^𝐴, 𝐹₂^𝐴, 𝐺₁^𝐴𝐵, 𝐺₂^𝐴𝐵 die Größe der Orthogonalisierungskorrektur an und sind anpassbar. ℋ_{𝜇𝜇,𝑋}^𝑙𝑜𝑐 ist ein Matrixelement eines lokalen Hamiltonoperators, der auf das Atompaar A und X wirkt.

OM3 unterscheidet sich von OM2 nur durch Vernachlässigung von den Korrekturtermen höherer Ordnung, d.h. in Gleichung 19 und 20 wird der letzte Term gleich null gesetzt.

Außerdem besitzen OM2 und OM3 im Vergleich zu OM1 Drei-Zentren Korrekturen für die

(25)

18 Resonanzintegrale, die dafür sorgen, dass Elektronenpaare in Konformeren richtig ausgerichtet sind. [37], [53]

2.10 Multireferenz-Konfigurationswechselwirkung

Häufig genügt die Betrachtung einer einzigen Referenzdeterminante nicht, um die statische Elektronenkorrelation mit zu berücksichtigen. Diese sind aber wichtig um Übergänge in angeregte Zustände beschreiben zu können. Für solche Fälle kann die Multireferenzkonfigurationswechselwirkungsmethode (MRCI, engl. Multireference Configuration Interaction) benutzt werden. Diese Methode verwendet mehrere Determinanten als Referenz. In Kombination mit semiempirischen Verfahren wie der OM3- Methode, kann dieser, ansonsten rechenaufwändige Prozess, in relativ kurzer Zeit durchgeführt werden. [36]

2.11 Dichtefunktionaltheorie

Die Grundzüge der Dichtefunktionaltheorie basieren auf der Feststellung von Hohenberg und Kohn [56], dass der Grundzustand allein durch die Elektronendichte beschrieben werden kann. Im Vergleich zur Beschreibung mit einer Wellenfunktion ist die Elektronendichte nur von drei Raumkoordinaten abhängig. Das Funktional, das eine solche Beschreibung vermittelt, ist jedoch unbekannt. Aufgabe der Dichtefunktionaltheorie ist es, ein solches Funktional zu finden. Dabei wird die Gesamtenergie zunächst mit drei einzelnen Termen beschrieben:

𝐸_𝐷𝐹𝑇 = 𝑇[𝜌] + 𝐸_𝑛𝑒[𝜌] + 𝐸_𝑒𝑒[𝜌] 21

𝑇[𝜌] ist die kinetische Energie, 𝐸_𝑛𝑒[𝜌] die Kernanziehung, die auf die Elektronen wirkt und 𝐸_𝑒𝑒[𝜌] ist die Elektronenabstoßung. Da die kinetische Energie mit den zuvor verwendeten Ansätzen, die ein gleichverteiltes Elektronengas beschrieben, nur unzureichend berechnet werden konnte, führten Kohn und Sham unter der Annahme nicht wechselwirkender Elektronen wieder Orbitale in die Dichtefunktionaltheorie ein. Durch diese lässt sich die kinetische Energie in einen exakt lösbaren Teil und einen Korrekturterm unterteilen:

𝐻_𝜆 = 𝑇 + 𝑉_𝑒𝑥𝑡(𝜆) + 𝜆𝑉_𝑒𝑒 22

Dabei ist 𝑉_𝑒𝑥𝑡(𝜆) der Operator für ein externes Potential. Er entspricht dem Potential der Kernanziehung 𝑉_𝑛𝑒, wenn 𝜆 = 1 ist. 𝑉_𝑒𝑒 entspricht dem Potential der

(26)

19 Elektronenwechselwirkung. Bei 𝜆 = 0 würde in diesem System keine Wechselwirkung zwischen den Elektronen stattfinden. Die Energiegleichung bei Verwendung von Kohn-Sham Orbitalen lässt sich dann beschreiben als:

𝐸_𝐷𝐹𝑇[𝜌] = 𝑇_𝑆[𝜌] + 𝐸_𝑛𝑒[𝜌] + 𝐽[𝜌] + 𝐸_𝑥𝑐[𝜌] 23

Dabei ist 𝐽[𝜌] der Coulombteil analog zum Hartree-Fock Coulomboperator und 𝐸_𝑥𝑐[𝜌] ist ein Austauschkorrelationsterm. [36] Das Auffinden eines geeigneten Austauschkorrelations- funktionals ist ein zentrales Anliegen der DFT-Methodik. Es gibt mehrere Ansätze, mit denen diese Problematik behandelt werden kann. In dem lokalisierte-Dichte-Näherung (LDA, eng.

Local Density Approximation) wird die Elektronendichte als gleichverteiltes Gas angenommen.

Der generalisierte-Gradienten-Ansatz (GGA, eng. Generalized Gradient Approximation) erweitert LDA, sodass die Elektronenkorrelation nicht nur von der Elektronendichte eines gleich verteilten Elektronengases abhängt, sondern auch noch von ihrer Änderung abhängt.

Damit können Ungleichmäßigkeiten besser berücksichtigt werden. Insgesamt liefert die Berechnung der Austauschenergie in der Dichtefunktionaltheorie schlechte Ergebnisse, weshalb oft sogenannte Hybridfunktionale verwendet werden. Diese mischen einen bestimmten Anteil der exakten Hartree-Fock Austauschenergie mit der Energie, die durch LDA oder anderer Ansätze berechnet wird. Zu dieser Gruppe an Funktionalen gehört auch das hier verwendete B3LYP-Funktional.

Um angeregte Zustände zu beschreiben reicht die normale DFT-Methode nicht aus. Deshalb wird zeitabhängige DFT (TDDFT) verwendet. Dabei geht man von einem störungstheoretischen Ansatz aus, um kleine Anregungen in höhere Zustände mithilfe der linearen Antworttheorie (eng. linear Response) zu beschreiben. Die berechneten Anregungsenergien werden auf die Energie des Grundzustands addiert, sodass eine Potentialfläche des betrachteten angeregten Zustands entsteht. Da allerdings einfach angeregte Slaterdeterminanten verwendet werden, erzielt die Berechnung von Zuständen, die durch andere Konfigurationen beschrieben werden, schlechte Ergebnisse.

Ein Nachteil der DFT-Methodik ist, dass die Funktionale nicht systematisch verbessert werden können und damit viele unterschiedliche Funktionale für viele unterschiedliche Aufgabenstellungen entwickelt wurden. [36]

(27)

20 2.12 QM/MM Methoden

Wenn das betrachtete System zu anspruchsvoll für eine reine ab initio Methode ist, können kombinierte quantenmechanische und molekülmechanische (QM/MM) Methoden zum Einsatz kommen. Diese unterteilen das System in mindestens zwei Teile: eine innere Region (I, s. Abb. 7 links), die mit einem möglichst genauen Verfahren betrachtet werden soll, und einen äußeren Bereich (O). Dieser spielt für das betreffende System eine untergeordnete Rolle und kann deshalb mit weniger aufwendigen Verfahren beschrieben werden. Zusätzlich dazu gibt es in den verwendeten Programmpaket COBRAMM [57] eine weitere Unterteilung in drei Bereiche, die die Beweglichkeit einzelner Systemteile angibt. In dem obersten Teil (H in Abb.

7 rechts) befinden sich alle Atome, die mit QM-Methoden beschrieben werden und sich während der Berechnung bewegen können. Die mittlere Ebene (M) besitzt die mobilen MM- Atome und auf der untersten Ebene (L) werden alle Atome fixiert. Eine Stoffklasse, bei der diese Methoden häufig eingesetzt werden, sind Proteine. Diese beinhalten oft ein aktives Zentrum, an dem die eigentliche Reaktion stattfindet, während die Proteinumgebung das aktive Zentrum nur indirekt beeinflusst, ohne selbst zu reagieren.

Abbildung 7: Schematische Darstellung der unterschiedlichen Ebenen einer QM/MM-Methode(links) und der beweglichen Region in COBRAMM (rechts). [57]

Das aktive Zentrum wird dann mit QM-Verfahren und das Proteinrückgrat mit MM-Methoden beschrieben. Das Proteinrückgrat kann normalerweise nicht ignoriert werden, da es die Form und Reaktivität des aktiven Zentrums beeinflusst.

Eine wichtige Entscheidung betrifft die Behandlung der Bindungsstellen zwischen den beiden Teilen (L in Abb. 7) und ihre Interaktion untereinander. Damit im QM-Bereich keine ungesättigten Bindungen verbleiben, werden diese oft mit sogenannten Verbindungsatomen abgesättigt. Dazu kommen meistens Wasserstoffatome zum Einsatz. Die Verbindungsatome

(28)

21 in dem, in dieser Arbeit verwendeten Programmpaket, besitzen keine MM Beschreibung, weshalb sie keine van der Waals Wechselwirkungen zu anderen MM Atomen aufweisen. Um die elektrostatische Wechselwirkung der Verbindungsatome mit anderen MM Atomen zu beschränken, muss die Ladung am eigentlich verbundenen MM Atom (M¹in Abb. 8) auf null gesetzt werden und diese unter den anderen MM Atomen gewichtet umverteilt werden.

In dieser Arbeit wurde die elektrostatische Einbettung (engl: Electrostatic embedding) verwendet, bei der zusätzlich zu den mechanischen Kräften auch das Ladungsfeld, welches durch die MM Atome erzeugt wird, mit den QM Atomen interagieren kann. Dabei gehen die Punktladungen der MM Atome als Ein-Elektron Terme in die Berechnung der QM Region mit ein. [36], [57]

Abbildung 8: QM und MM Regionen. In grün ist die Grenze eingezeichnet. Das Verbindungsatom ist mit L beschriftet. [57]

Cobramm verwendet eine subtraktive Formel für die Berechnung der Energie:

𝐸_{𝑄𝑀/𝑀𝑀}(𝑆) = 𝐸_𝑀𝑀(𝑆) + 𝐸_𝑄𝑀(𝐼 + 𝑙𝑖𝑛𝑘) − 𝐸_𝑀𝑀(𝐼 + 𝑙𝑖𝑛𝑘) 24

Dabei steht S für das gesamte System, I für den High-Layer, link für die Verbindungsatome.

[57]

2.13 Moleküldynamik

Eine Moleküldynamik Methode betrachtet ein System mit gegebenen Anfangskoordinaten und -geschwindigkeiten über einen bestimmten Zeitraum mit Newtons zweitem Gesetz. Für Moleküldynamiken werden die Kräfte, also die ersten Ableitungen der Kernkoordinaten nach

(29)

22 der Energie, benötigt. Werden parametrisierte Kraftfeldmethoden für die Dynamik verwendet, können die Kräfte fast genauso schnell berechnet werden wie die Energie.

Die Bewegung der Atomkerne kann in guter Näherung durch Newtons zweites Gesetz beschrieben werden:

−𝑑𝑉

𝑑𝒓 = 𝑚𝑑²𝒓 𝑑𝑡²

25

Dabei ist V die Geschwindigkeit der Kerne, r bezeichnet die Kernkoordinaten, m ist die Masse und t die Zeit. [36]

2.14 Der Velocity Verlet Algorithmus

Ein oft verwendeter Algorithmus in Moleküldynamiken ist die velocity Verlet Routine. [58] Er beruht auf der Bestimmung der Halbschrittgeschwindigkeiten:

𝑣_𝑡+1

2Δ𝑡 = 𝑣_𝑡+𝑓_𝑡𝛥𝑡 2𝑚

26

Zuerst wird die Geschwindigkeit 𝑣_𝑡+¹

2Δ𝑡 einen halben Zeitschritt nach dem Zeitpunkt t betrachtet. Dabei ist 𝑓_𝑡 die Kraft am aktuellen Zeitpunkt t und m die Masse.

Anschließend wird mit dieser Geschwindigkeit die Position der Kerne zum Zeitpunkt 𝑡 + Δ𝑡 berechnet:

𝑟_𝑡+Δ𝑡 = 𝑟_𝑡+ 𝑣_𝑡+1

2Δ𝑡𝛥𝑡 27

Diese neue Geometrie kann dann verwendet werden, um die Kräfte zum Zeitpunkt t zu bestimmen. Die korrigierte Geschwindigkeit ergibt sich dann zu:

𝑣_𝑡+Δ𝑡 = 𝑣_𝑡+1

2Δ𝑡+𝑓_𝑡+Δ𝑡𝛥𝑡 2𝑚

28

Man erhält so Informationen über die zeitliche Änderung der Kernbewegung, d.h. eine sog.

Trajektorie. Dieser Algorithmus ist selbststartend, d.h. er benötigt jeweils nur die Geschwindigkeiten, Kräfte und die Geometrie zum Zeitpunkt t. Er ist zeitreversibel und der

(30)

23 Fehler durch die Halbschrittgeschwindigkeiten ist vom gewählten Zeitschritt selbst abhängig.

Dies führt, bei Wahl eines entsprechenden 𝛥𝑡, zu einer besseren Energieerhaltung auch bei langen Simulationszeiten. [36], [59]

2.15 Der Rattle-Algorithmus

Der Zeitschritt in MD-Simulationen muss in der Regel so gewählt werden, dass die schnellste Schwingung des Systems ausreichend beschrieben wird. Um einen größeren Zeitschritt zu verwenden, können schnelle Freiheitsgrade, wie z.B C-H Schwingungen, eingefroren werden.

Der Rattle-Algorithmus [47] bietet die Möglichkeit, Bindungslängen zu fixieren, ohne dabei an Genauigkeit zu verlieren. Diese Eigenschaft ist besonders beim TIP3P Wassermodell, welches in der Simulation starr gehalten werden muss, von Nutzen. Dabei wird der Abstand zwischen den beiden Wasserstoffatomen sowie die O-H Abstände festgesetzt und die Geschwindigkeiten werden so umverteilt, dass sich das Wassermolekül als starre Einheit in der Simulation bewegt. Dies wird erreicht, indem die entsprechenden Restriktionen im Velocity Verlet Algorithmus (s. Kapitel 2.14) mitberücksichtigt werden. Die Berechnung der korrigierten Geschwindigkeiten geschieht dabei iterativ. [47]

2.16 Wigner-Sampling

Im Wigner-Sampling verwendet man eine Pseudoverteilung, um eine passende quantenmechanische Beschreibung von Koordinaten und Impulsen eines Systems zu erhalten.

Dabei wird von folgender Gleichung ausgegangen:

𝜌_𝑊(𝑸_𝑖, 𝑷_𝑖) = (𝜋ħ)⁻¹exp (−𝜇_𝑖𝜔_𝑖𝑸_𝑖²

ħ ) ∗ exp ( 𝑷_𝑖²

𝜇_𝑖𝜔_𝑖ħ) 29

Dabei ist 𝜌_𝑊(𝑸_𝑖, 𝑷_𝑖) die Pseudoverteilung, 𝜇_𝑖 die reduzierte Masse, 𝜔_𝑖 die harmonische Frequenz der Normalmode i und Pi der Impuls, der mit einer bestimmten Normalkoordinate Qi verbunden wird.

Mithilfe dieser Gleichung kann für jede Normalmode i ein Paar von zufälligen Auslenkungen und Impulsen erhalten werden. Dieser Prozess wird für jede Normalmode wiederholt.

Schließlich werden die Normalkoordinaten in kartesische Koordinaten transformiert. Nach dem Aufaddieren aller einzelnen Komponenten für die Impulse und Koordinaten kann das Ergebnis dieses Verfahrens als Startbedingung für Moleküldynamiken verwendet werden.

Durch Umrechnen erhält man dabei eine definierte Geometrie und einen entsprechenden

(31)

24 Geschwindigkeitsvektor. Durch Wiederholung dieser Prozedur lässt sich ein Ensemble von Molekülen generieren, die im Nullpunktspotenzial schwingen. Man spricht dann vom sog.

Sampling. Bei einer zusätzlichen Schwingungsanregung wird in den entsprechenden Moden die Energie um ein oder mehrere Quanten erhöht. Damit vergrößert sich auch die kinetische und die potentielle Energie, wodurch beim Sampling die Auslenkungen und die entsprechenden Impulse für eine Vibrationsanregung erhalten werden. [60], [61]

2.17 “Tully’s fewest switches surface hopping” Algorithmus

Dieser Algorithmus wird verwendet, um elektronische Übergänge im Verlauf von Trajektorien zu berechnen. Er baut auf der „surface-hopping“ Methode auf. Das Besondere an dieser Methode ist, dass sie den elektronischen Übergang zwischen zwei Zuständen beschreiben und für beliebig viele elektronische Zustände angewandt werden kann.

Namensgebend für diesen Algorithmus ist, dass es entlang einer Trajektorie so wenig Übergänge wie möglich geben soll. Das erlaubt es der Trajektorie, dass sie auf einer Potentialfläche wandert und nicht durch häufige Übergänge einen Mittelwert der beiden Energiepotentiale beschreiben soll. Ein weiterer zentraler Ansatz dieser Methode ist, dass Tully eine Verbindung zwischen dem Verhalten der Gesamtzahl aller Trajektorien und dem Verhalten der einzelnen Trajektorien herstellt. Dadurch, dass er die Anzahl der Übergänge so gering wie möglich hält, wird die richtige statistische Verteilung der Populationen über alle Trajektorien beibehalten. Das Quadrat der Koeffizienten 𝑐_𝑖(𝑡) des Zustands i zum Zeitpunkt t repräsentiert dabei immer die Anzahl der Trajektorien 𝑁_𝑖(𝑡) im i-ten Zustand:

𝑁_𝑖(𝑡) = 𝑐_𝑖^∗(𝑡)𝑐_𝑖(𝑡)𝑁 30

Laut Tully ergeben sich damit für die Wahrscheinlichkeiten eines Sprungs von Zustand i zu j die Wahrscheinlichkeit:

𝑃_𝑖𝑗 = 2∫_𝑡^𝑡+Δ𝑡𝑅𝑒(𝑐_𝑖^∗(𝑡)𝑐_𝑗(𝑡)𝑅̇𝑑_𝑖𝑗)𝑑𝑡 𝑐_𝑖^∗(𝑡)𝑐_𝑖(𝑡)

31

Dabei ist 𝑑_𝑖𝑗 der nichtadiabatische Kopplungsvektor zwischen den Zuständen i und j.

Diese Wahrscheinlichkeit wird dann mit einer Zufallszahl verglichen, um zu entscheiden, ob ein Zustandswechsel stattfindet. Abbildung 9 stellt den Unterschied zwischen „surface-

(32)

25 hopping“ ohne und mit „Tully’s fewest switches“ da. Zu sehen ist, dass das „surface-hopping“

ohne den Tullys Algorithmus (siehe Abbildung 9, Bild c)) sehr oft zwischen den beiden Potentialkurven wechselt. Daher wird dieses Verhalten eher eine gemittelte Potentialkurve wiedergeben (Abbildung 9, Bild b)) gestrichelte Linie). [62]

Ein Problem in diesen Algorithmus betrifft die Kohärenz. Eine Anzahl an Trajektorien kann kein quantenmechanisches System komplett beschreiben, da jede individuelle Trajektorie bei einem Wechsel zwischen zwei Potentialflächen in einer von Beiden bleiben muss, ohne dass ein entsprechender Teil auf der jeweiligen anderen Fläche bleibt. In einem richtigen quantenmechanischen System würde sich das Wellenpaket, entsprechend den Wahrscheinlichkeiten, auf die beiden Flächen aufteilen. Um diesem Problem zu begegnen, gibt es ein empirisches Verfahren, welches dämpfend auf die berechneten Amplituden wirkt und die Wahrscheinlichkeiten in Gleichung 31 entsprechend korrigiert. In dem verwendeten Programm Cobramm wurde das Dekohärenzkorrektur-Schema verwendet, dass von Truhlar et al.[63] vorgeschlagen und später durch Persico et al. [64] umgesetzt wurde. [37], [57], [62]

Abbildung 9: Schematische Darstellung der Unterschiede zwischen einem Algorithmus der häufiger die Potentiale wechselt und "Tully' fewest switches" Algorithmus. a) Darstellung der zeitlichen Entwicklung zweier Populationen. b) Die dazugehörigen Potentiale (durchgezogene Linien) und dem effektiven Potential (gestrichelte Linie). c) Darstellung eines Algorithmus, der oft das Potential wechselt und somit eher das effektive Potential beschreibt. d) Zwei mögliche Trajektorien die durch „Tully’s fewest switches„ Algorithmus entstehen können. [62]

(33)

26

3 Ergebnisse

3.1 Technische Details

Insgesamt wurde die Moleküldynamik von vier Modellsystemen untersucht. Dabei wurden die Proteinstrukturen 1C3W [6] und 4Y9H [65] für bR verwendet. Für die 1C3W Struktur wurde die zunächst Parametrisierung von P. Tuchel [25], basierend auf dem Amber ff99SB [41]

Kraftfeld, übernommen. Die Struktur 4Y9H wurde, unter Verwendung des Amber ff14SB [42]

Kraftfelds, für einen pH Wert von 7.0 mit H⁺⁺[66] protoniert, und die Atompositionen der Wasserstoffe wurden bei festgehaltenem Proteinrückgrat optimiert. Für die QM/MM Rechnungen wurde das Programm Cobramm [57] verwendet, welches die in dieser Arbeit verwendeten QM-Programme Gaussian 09 [67], MNDO 7.0 [68] und MOLCAS 8.1 [69] sowie das Programmpaket Amber17 [70] kombiniert. Für die Berechnung der Grundzustandsstrukturen wurden die OM3 Methode [53-54], DFT/B3LYP [71–74] mit der Basis 6-31G* [75] und CASSCF(12,12) [36] in gleicher Basis verwendet. Für die Analyse der elektronischen Struktur wurde TD-CAM-B3LYP [76] für die B3LYP Struktur verwendet. MS- CASPT2 [49] kam für die CASSCF Struktur zum Einsatz. Dabei wurde die Wellenfunktion über fünf Zustände gemittelt. Die Berechnung erfolgte ohne IPEA Shift, aber mit einem imaginären Shift von 0.2 in der 6-31G* Basis.

Die Moleküldynamiken wurden, wenn nicht anders angegben, mit OM3/MRCI berechnet.

Dabei wurde die OM3/UHF Wellenfunktion als Referenzfunktion verwendet. Als aktives Fenster für das MRCI-Verfahren wurden die 10 höchsten besetzten und die 10 niedrigsten unbesetzten Orbitale im Grundzustand gewählt (d.h. 20 Orbitale und 20 Elektronen). Für die Bestimmung der Übergangswahrscheinlichkeiten kam der "Tully's fewest switch surface hopping“ [62] Algorithmus mit Dekohärenzkorrektur zum Einsatz. Die Berechnung der Dynamik wurde beendet, wenn die Moleküle mindestens 50 fs im Grundzustand verbracht haben und das resultierende Photoprodukt eindeutig erkennbar war. Der Rattle [47]

Algorithmus wurde für die Fixierung des TIP3P Wassermodells verwendet. Startgeometrien für die Dynamiken wurden mit Hilfe des Wignersamplings [60-61] erstellt. Für die Auswertung der Trajektorien und Geometrien wurde das Programm Visual Molecular Dynamics (kurz:

VMD) [77] oder Molden [78] verwendet. Um die Reaktionszeiten zu bestimmen, wurde der Zerfall der S1 Population mit folgender Formel exponentiell gefittet: