Simulated Annealing angewandt auf das Problem der Kristallstrukturvorhersage

Diplomarbeit – Wilfried Bollweg

25. Juli 1994

INSTITUT FÜR NUMERISCHE UND INSTRUMENTELLE MATHEMATIK DER WESTFÄLISCHEN WILHELMS UNIVERSITT MNSTER

THEMENVERGABE: PROF. DR. H. MAURER

Inhaltsverzeichnis

Vorwort ii

1 Das Optimierungsverfahren 1

1.1 Simulated Annealing: Versuch einer Standortbestimmung … 1 1.1.1 Vorstellung des allgemeinen Optimierungsproblems . . 1 1.1.2 Klassifizierung globaler Optimierungsprobleme . . . . . 3 1.1.3 Klassifizierung der Optimierungsstrategien . . . . . . . 4 1.1.4 Beschreibung des Simulated Annealing . . . . . . . . . 5

1.2 Die Arbeitsweise des Simulated Annealing . . . . . . . . . . . 6 1.2.1 Darstellung des Grundalgorithmus . . . . . . . . . . . 6 1.2.2 Konvergenzbeeinfluung durch Parametermodifikation . 9 1.2.3 Diskussion . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.3 Die Konvergenz des Simulated Annealing . . . . . . . . . . . . 11 1.3.1 Simulated Annealing und die stochastischen Prozesse . 11 1.3.2 Konvergenztheorem . . . . . . . . . . . . . . . . . . . . 14 1.3.3 Diskussion . . . . . . . . . . . . . . . . . . . . . . . . . 15

1.4 Anwendung des Simulated Annealing auf globale Optimierungsprobleme mit nichtdiskretem Konfigurationenraum …. 16 1.4.1 Die Algorithmen in SIMULA . . . . . . . . . . . . . . 16 1.4.2 Erluterung der Abschluverfeinerung . . . . . . . . . . . 21 1.4.3 Beispiel fr die Anwendung einer Abschluverfeinerung . 24 1.4.4 Die Parameter der Simulated-Annealing-Algorithmen . 26 1.4.5 Lnge der Markovketten . . . . . . . . . . . . . . . . . . 35

1.5 Das Ma fr die Qualitt von Optimierungsverfahren: die Testfunktionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 1.5.1 Die Testfunktionen als Qualittsurteil fr das Simulated

Annealing . . . . . . . . . . . . . . . . . . . . . . . . . 35 1.5.2 Beispiele fr Zielfunktionen . . . . . . . . . . . . . . . . 36 1.5.3 Diskussion . . . . . . . . . . . . . . . . . . . . . . . . . 41

2 Kristallstrukturvorhersage 42

2.1 Kristalle als Objekte der Mathematik . . . . . . . . . . . . . . 42 2.1.1 Eigenschaften kristalliner Festkrper . . . . . . . . . . . 42 2.1.2 Bindungsenergien bei der Kristallstruktur-Vorhersage . 44 2.1.3 Mathematisches Modell der Elementarzelle . . . . . . . 48 2.2 Anpassung des Simulated Annealing - Verfahrens . . . . . . . 49

2.2.1 Zustandsbergnge auf dreidimensionalen Kugelsphren . 49 2.2.2 Unbewegliche Teilchen . . . . . . . . . . . . . . . . . . 50 2.2.3 Das ”Winkel–Abstands–Programm“ . . . . . . . . . . . 51 2.2.4 ”Koordinatenkollisionen“ . . . . . . . . . . . . . . . . . 51 2.2.5 Berechnung der Koordinationspolyeder . . . . . . . . . 52

2.3 Formulierung der verschiedenen Zielfunktionen . . . . . . . . . 52 2.3.1 Die Zielfunktion nach Pannetier et al. . . . . . . . . . . 53 2.3.2 Die modifizierte Zielfunktion und der erweiterte Ansatz 56 2.3.3 Die konvergenzbeschleunigende Methode von Ewald,

Bertaud und Williams . . . . . . . . . . . . . . . . . . 58

2.4 Zusammenfassung der Zielfunktionen . . . . . . . . . . . . . . 62

3 Berechnete Kristalle 65 3.1 Csiumchlorid . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.2 Kochsalz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.3 Zinksulfit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.4 Calciumfluorid . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.5 Rutil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.6 Perovskit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 3.7 Quarz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 3.8 Diskussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.9 Ausblick . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

Anhang 87

A.1 Die Optimierungsroutine SIMULA . . . . . . . . . . . . . . . 87 A.1.1 Verwendung der Optimierungsroutine SIMULA …. 87 A.1.2 Vom Benutzer zu erstellende Routinen . . . . . . . . . 89

A.2 Das Kristallstruktur-Vorhersage -Programm ANNEAL …. 91 A.2.1 Verfgbare Versionen . . . . . . . . . . . . . . . . . . . . 91 A.2.2 Maximale Anzahl von Atomen . . . . . . . . . . . . . . 91 A.2.3 Syntax des Programmaufrufes . . . . . . . . . . . . . . 91 A.2.4 Aufbau eines Datenfiles fr ANNEAL . . . . . . . . . . 92 A.2.5 Erluterung der Parameter . . . . . . . . . . . . . . . . 93

A.3 Beschreibung einer Ausgabedatei von ANNEAL . . . . . . . . 100

Literaturverzeichnis 105

Vorwort

Die dramatische Entwicklung der Computer bedeutete für viele Bereiche der Naturwissenschaften und der Mathematik in den letzten Jahren einen entscheidenden Schritt nach vorn. War es vor einigen Jahren nur unter Aufwendung von viel Geld und unter Einsatz von Großrechenanlagen möglich, Optimierungsprobleme mit vielen Variablen zu bearbeiten, reicht heute oft schon ein einfacher Personalcomputer für deren Lösung aus. Dies hatte zur Folge, daß immer komplexere Probleme in Angriff genommen wurden und auch die Algorithmen trugen der Entwicklung Rechnung, daß Vorgänge in der Natur simuliert werden konnten. Es entstanden neue Strategien zur Lösung numerischer Probleme. Ein wichtiger Vertreter dieser Simulationsalgorithmen ist das Simulated Annealing. Die Methode läßt sich aus dem Abkühlvorgang (→ to anneal: engl. für abkühlen) von Festkörpern herleiten. Bringt man einen Festkörper durch Energiezufuhr in den Zustand der Schmelze und läßt ihn langsam genug abkühlen, so ordnen sich seine Atome nach dem Prinzip minimaler innerer Energie an. Diese Beobachtung machten sich als erste S. Kirkpatrik, C.D. Gelatt und M.P. Vecci (1982) zunutze [14], um neue Ansätze zur Lösung von nichtlinearen Optimierungsproblemen zu entwickeln. Der größte Bereich von Anwendungen, Weiterentwicklungen und Modifikationen des Simulated Annealing ist immer noch der Bereich der kombinatorischen Optimierung. Hier wird in diskreten, endlichen Räumen versucht, nichtlineare Optimierungsprobleme zu lösen. Ein immer wieder auftauchendes Problem aus diesem Bereich ist das sogenannte Travelling-Salesman-- Problem, bei dem es gilt, den Reiseweg eines Handlungsreisenden durch eine fest vorgegebene Anzahl von Städten zu optimieren. Ein Beispiel für den Erfolg des Simulated Annealing in diesem Bereich findet sich bei Aarts und van Laarhoven [2]. Ein weiterer großer Anwendungsbereich für das Simulated Annealing fand sich in den letzten Jahren in der Optimierung stetiger Funktionen. Litten vor etwa 15 Jahren alle Algorithmen der nichtlinearen Optimierung unter dem Defizit, daß ausgehend von einem Startwert immer wieder nur das nächste Minimum in Abstiegsrichtung gefunden wurde, so erhielt man mit dem Simulated Annealing eine zweifelsohne sehr rechenintensive, aber auch sehr stabile und erfolgreiche Methode, die in der Lage war, globale Minima zu finden. Darüber hinaus war man mit Mitteln der Wahrscheinlichkeitstheorie sogar in der Lage, die überraschende Eigenschaft zu zeigen, daß globale Minima mit der Wahrscheinlichkeit Eins angenommen werden, vorausgesetzt, daß alle Parameter des Verfahrens optimal gewählt sind. Anwendungsbeispiele für das Simulated Annealing auf Optimierungsprobleme stetiger Funktionen findet man in so unterschiedlichen Bereichen wie der Bildverarbeitung, der Plazierung und Verdrahtung elektronischer Bauelemente auf Leiterplatinen, der Entwicklung von Abfangmechanismen von Raketen und im Umweltschutz. Nun hat aber ein Verfahren mit den Referenzen des Simulated Annealing, versucht man es in die Tat umzusetzen, letztlich auch immer seine Schwachstellen. Im ersten Kapitel werde ich deshalb nach einer Erläuterung des Konvergenztheorems auf die Eigenschaften von verschiedenen Implementierungen des Simulated Annealing eingehen. Weitere wesentliche Inhalte des ersten Kapitels bestehen in der Darstellung des physikalischen Grundgedankens des Simulated Annealing sowie der Vorstellung der verschiedenen Implementationsmöglichkeiten, die sich aus den Diplomarbeiten von Frank Krimphove [17] und Jens Radzek [26] ergeben. Das zweite Kapitel geht genauer auf das Problem der Kristallstrukturvorhersage ein. Hierbei handelt es sich um eine Anwendung des Simulated Annealing, die direkt auf den physikalischen Grundgedanken des Verfahrens basiert. Es zeigt sich, daß auch vermeintliche Nachteile des Verfahrens wie zum Beispiel die Eigenschaft, daß bei Implementierungen immer auch lokale Minima möglich sind, durchaus physikalisch eine Daseinsberechtigung besitzen, gerade wenn es darum geht, einen Prozeß des Abkühlens von Festkörpern zu betrachten. Desweiteren ergibt sich aus der Arbeit mit dem Kristallgittermodell eine Verbesserung des Algorithmus von Radzek, die sich bewährt hat. Auch diese soll vorgestellt werden. Sie wurde nachträglich in das Simulated-Annealing-Programm mit aufgenommen. Der Abschluß der Arbeit besteht aus der Beschreibung von Ein- und Ausgabedateien für die Kristallstrukturvorhersage sowie einer Erläuterung, wie die aktuelle Simulated-Annealing-- Routine zu benutzen ist. Auf die Listings des Kristallstrukturvorhersageprogramms ANNEAL und der Optimierungsroutine SIMULA wurde aufgrund des Umfangs verzichtet. Sowohl die Optimierungsroutine als auch das Programm zur Kristallstrukturvorhersage befinden sich auf der beiliegenden Diskette.

Kapitel 1

Simulated Annealing als Optimierungsverfahren

1.1 Simulated Annealing: Versuch einer Standortbestimmung

Viele praktische Probleme der Natur– und Ingenieurwissenschaften lassen sich als globale Optimierungsprobleme formulieren. Die Menge der zur Verfügung stehenden Verfahren ist fast unüberschaubar. Viele dieser Verfahren scheinen auf fast jedes Problem anwendbar zu sein. Selten wird jedoch versucht, die Grundsätze der Verhältnismäßigkeit zu beachten. Es reicht nicht, das Problem zu formulieren und nach einem funktionierenden Verfahren zu suchen, sondern es gilt, auch über alternative Verfahren nachzudenken. So ist es beispielsweise unnötig, eine Zielfunktion mit vielen lokalen Minima und einer unübersichtlichen Topologie durch das Anwenden einer Newtonmethode auf globale Minima zu untersuchen. Besteht keine genauere Kenntnis über die Lage der globalen Minima, würde ein solches Vorgehen in wenig erfolgreiches Raten eines Startwertes ausarten. Um dies zu verhindern, möchte ich nach der Formulierung des allgemeinen, globalen Optimierungsproblems eine Klassifizierung der Zielfunktionen sowie der zur Verfügung stehenden Optimierungsstrategien vorstellen. Diese hält sich weitestgehend an die von Törn und Ziliskas [29] erarbeitete Klassifizierung und hat sich bei der Beurteilung von Problemen als äußerst hilfreich erwiesen.

1.1.1 Vorstellung des allgemeinen Optimierungsproblems

Betrachtet man ein wirtschaftswissenschaftliches Modell, so stellt sich dieses im allgemeinen als abgeschlossenes System dar, welches durch eine Kombination seiner Parameter geprägt ist. Diese Parameter selbst bewegen sich wieder innerhalb festgesetzter Grenzen. Um ein Mittel zu besitzen, die Wirtschaftlichkeit des Systems beurteilen zu können, führt man eine Qualitätsfunktion, die sogenannte Zielfunktion, ein. Diese hängt direkt oder indirekt von den Parametern des Systems ab. Will man das System optimieren, gilt es, die Parameter so zu wählen, daß die Zielfunktion des System minimal wird. Geht man vom wirtschaftswissenschaftlichen Modell aus, so hat man die Kosten des Systems minimiert. Der Begriff der ”Kosten“ eines Systems hat sich beim Umgang mit Optimierungsproblemen etabliert. Der Begriff wird im Folgenden häufig verwandt. Mathematisch stellt sich das Problem folgendermaßen dar: Es seien

f: IR^{n} ⟶ IR stetig differenzierbar g: IR^{n} ⟶ IR^{m} stetig differenzierbar.

Dann heißt das Problem

P: min f(x) unter den Nebenbedingungen g_{i}(x) = 0 für i ∈ {1, …, m_{e}} g_{i}(x) ≤ 0 für i ∈ {m_{e+1}, …, m}

ein nichtlineares Optimierungsproblem unter den Nebenbedingungen g_{i}. Die g_{i} werden auch Restriktionen genannt. Die Menge der Werte, die durch die Restriktionen beschränkt wird, nennt man zulässiger Bereich oder auch Konfigurationenraum. Jeder Wert im zulässigen Bereich heißt Konfiguration oder Zustand. Für die im Folgenden betrachteten Probleme ist diese Darstellung jedoch zu allgemein. Es werden nur Probleme der Form

P(x): min f(x) unter den Nebenbedingungen a_{i} ≤ x_{i} ≤ b_{i} für i ∈ {1, …, m}

behandelt. Wendet man lineare Transformationen der Form

T(x) := p · x + q mit p, q ∈ IR

auf die x_{i} an, sind diese Probleme äquivalent zu Problemen der Form

P’(x) = P(T(x)): min f(T (x)) unter den Nebenbedingungen

a ≤ x_{i} ≤ b für i ∈ {1, …, m}.

Beispiel:

betrachte P : min f(x_{1}, x_{2}) = 2 · x_{1}^{2} + 3 · x_{2}
0 ≤ x_{1} ≤ 2
2 ≤ x_{2} ≤ 3.

Wendet man die linearen Transformationen der Form:

T_{1}(x_{1}) := \frac{1}{2} · x_{1}

T_{2}(x_{2}) := x_{2} − 2

auf die x_{i} des Problems P an, so ergibt sich das äquivalente Problem

P′ : min f(T_{1}(x_{1}), T_{2}(x_{2})) = 2\frac{x_{1}}{2}^{2} + 3 · (x_{2} − 2)

0 ≤ x_{1}, x_{2} ≤ 1.

1.1.2 Klassifizierung globaler Optimierungsprobleme

Aufwendige Zielfunktionen

Viele Probleme unterscheiden sich erheblich im Aufwand, der zur Berechnung einer Funktionsauswertung notwendig ist. Natürlich beeinflußt dieser Charakterisierungspunkt direkt die Wahl der Optimierungsstrategie. Kommt es darauf an, zwischen zwei Verfahren zu wählen, sollten Vor- und Nachteile der Verfahren gegeneinander abgewogen werden. Das Klassifizierungsmerkmal des Aufwands zur Berechnung einer Funktionsauswertung wird vor allem bei der Bearbeitung der kristallographischen Zielfunktionen von entscheidender Bedeutung sein, da die verwendeten Funktionen sehr rechenintensiv sind. Im Gegensatz dazu stehen Zielfunktionen, die in geschlossener Form vorliegen und mit einigen wenigen grundlegenden Operationen zu berechnen sind. Beispielhaft hierfür sind die später noch genauer zu behandelnden Testfunktionen.

Anzahl der lokalen Minima

Die Anzahl der lokalen Minima ist ein weiterer Klassifizierungspunkt. Probleme, die sehr wenige lokale Minima im zulässigen Bereich besitzen, sind offensichtlich leichter zu bearbeiten als Probleme mit sehr vielen Minima, da im ersten Fall die Anzahl der möglichen Ergebnisse kleiner ist. Viele unter Umständen gleichwertige Minima bedeuten jedoch nicht, daß ein Problem nicht lösbar ist. Beispielhaft hierfür mag das physikalische Modell der Kristalle sein. Jedes Atom im Kristall wird mit drei Koordinaten identifiziert. Diese Koordinaten bilden ein Tripel im Zustandsvektor des Optimierungsproblems. Befinden sich mehrere gleichartige Atome im Kristall, so macht es Sinn, diese auszutauschen. Hierdurch ergeben sich häufig viele gleichwertige Minima. Ferner existieren beim Austausch von Atomen mit ähnlicher formaler Ladung noch metastabile Strukturen, die sogenannten ”Nebenstrukturen“. Sie unterscheiden sich vom Funktionswert her nur unwesentlich vom globalen Minimum. Jede Lösung hat aufgrund des zugrundeliegenden physikalischen Modells jedoch ihre Bedeutung und kann als solche akzeptiert werden. Dieses Argument erlaubt eine befriedigende Arbeit mit den kristallographischen Zielfunktionen und rechtfertigt eine Anwendung des Simulated Annealing.

Streuung der Minima

Minima, die sich in der Nähe von anderen Minima befinden, haben einen relativ kleinen ”Einzugsbereich“, in dem die Zielfunktion monoton fallend ist. Ist die Zielfunktion stetig mit vielen Minima auf engem Raum und liegen nur wenige Minima weit verstreut, so ist die Wahrscheinlichkeit größer, daß die weit verstreuten Minima gefunden werden. Dieses Charakterisierungsmerkmal wird im Folgenden nicht von entscheidender Wichtigkeit sein, da es für das physikalische Modell des Kristalls nicht zutrifft.

Problemdimension

Eine entscheidende Charakteristik ist die Problemdimension. Es läßt sich einfach einsehen, daß ein höherdimensionales Problem erheblich rechenzeit– und speicherplatzaufwendiger ist und damit erheblich höhere Ansprüche an eine Optimierungsroutine stellt. Aufgrund der oft sehr großen Zahl von Atomen pro Kristall wird diese Charakteristik bei der Bearbeitung der Probleme einen hohen Stellenwert einnehmen. Sind Probleme mit etwa 50 Variablen in noch akzeptabler Zeit auch auf einem Personalcomputer zu bearbeiten, so sind Probleme mit etwa 100 Variablen selbst auf einem Großrechner fast nicht mehr durchführbar. Daß auch das Simulated Annealing trotz relativ großer Rechenzeiten unter gewissen Einschränkungen diese Probleme lösen kann, soll im Folgenden gezeigt werden.

1.1.3 Klassifizierung der Optimierungsstrategien

Theoretische Lösungsansätze

Basis für alle theoretischen Behandlungen von differenzierbaren Optimierungsproblemen bildet der Satz von Kuhn und Tucker. Diese Methode findet sich in annähernd jedem Buch über nichtlineare Optimierung und ist zum Beispiel genauer von Maurer [21] beschrieben. Die Suche nach Minima erfolgt durch das Lösen von Gleichungssystemen und ist aufgrund des erheblichen Aufwands nur für kleine Dimensionen praktikabel.

Praktische Lösungsansätze

1. Deterministische Verfahren: Im wesentlichen umfaßt diese Gruppe die Intervall-verfahren, Gradienten-Verfahren oder Quasi–Newton-Ver-fahren. Als Vertreter der deterministischen Verfahren sind zu nennen:

(a) OPRQP (siehe Bartolomäus & Biggs, [?])

(b) E04JAF (siehe Beschreibung NAG Bibliothek [?])

(c) DBCPOL (siehe Beschreibung IMSL Bibliothek [?])

2. Stochastische Verfahren: Zu dieser Gattung gehört auch das Simulated Annealing Verfahren. Die Gattung der stochastischen Verfahren läßt sich unterteilen in:

(a) einfache Zufallsalgorithmen

(b) zusammengesetzte Zufallsalgorithmen

(c) stochastische Approximationen der Zielfunktion.

Die deterministischen Verfahren sind hierbei die ältere der beiden Gattungen. Die Verfahren sind numerisch stabil und erprobt, haben jedoch einen Fehler: erreichen die Verfahren einmal ein lokales Minimum, so sind sie nicht mehr in der Lage, es zu verlassen. Hierdurch wurde der Begriff des ”Hängenbleibens“ geprägt. Ob das globale Minimum tatsächlich erreicht wird, hängt also stark von der Wahl des Startwertes ab. Aus dieser Schwäche entwickelte sich in den letzten Jahrzehnten die Gattung der stochastischen Algorithmen dadurch, daß man versuchte, das Phänomen des ”Hängenbleibens“ durch zufällige Sprünge zu überbrücken. 1982 entwickelte S. Kirkpatrick dann ein iteratives Optimierungsverfahren zur Lösung von Problemen der obigen Art [14]. Aus diesem Verfahren entstand das Simulated Annealing. In der Klassifikation der Lösungsansätze läßt sich das Simulated Annealing am besten in die Klasse der zusammengesetzten Zufallssalgorithmen einordnen, da einfache Zufallsstrategien am Anfang des Algorithmus mit deterministischen Verfahren am Ende des Annealingprozesses kombiniert werden.

1.1.4 Beschreibung des Simulated Annealing

Die Diskretisierung des Ausgangsproblems

Die grundsätzliche Idee des Simulated Annealing besteht darin, Probleme in diskreten Mengen zu lösen. Es lassen sich also Probleme der folgenden Form bearbeiten:

Sei A eine endliche Menge und sei C: A ⟶ IR eine Kostenfunktion

P’ : min_{i∈A} C(i).

Es wird klar, daß das neue Optimierungsproblem P’ nichts anderes als einen diskretisierten Spezialfall des Optimierungsproblems P von Seite 2 darstellt. Die Diskretisierung von P ermöglicht eine Bearbeitung des Optimierungsproblems auf konkreten Rechenmaschinen. Die Schwierigkeit, nicht diskrete Optimierungsprobleme mit diskreten Algorithmen zu behandeln, wird im Folgenden noch genauer behandelt.

Physikalischer Grundgedanke

Die Idee des Simulated Annealing findet ihren Ursprung in der Festkörperphysik. Ein fester Stoff wird solange erhitzt, bis er sich im flüssigen Zustand befindet. Interpretiert man die Bausteine des Stoffes als Massepunkte, die in einem dreidimensionalen Quader angeordnet sind, so ergibt sich nach dem Aufheizen des Stoffes ein chaotisches Bild. Die dem System durch das Erhitzen zugeführte Energie bewirkt, daß die Bausteine sich ”frei“ im Quader bewegen und alle Aufenthaltsorte sowie alle Bewegungsrichtungen annähernd gleich wahrscheinlich sind. Während des Abkühlens ordnen sich die Bausteine dann nach dem Prinzip der geringsten freien Energie an. Zustandsübergänge zu einem energetisch schlechteren Zustand sind, je weiter der Abkühlungsvorgang fortgeschritten ist, um so unwahrscheinlicher. Die Wahrscheinlichkeit, daß sich ein System mit der aktuellen Energie E_{int} zu einem Zeitpunkt T in einem Zustand der Energie E zu befinden, wird durch die sogenannte Boltzmann-Verteilung beschrieben:

P (\{E_{int} = E\}) = \frac{1}{Z(T)} · exp (−\frac{E_{int}}{k_{B}·T} ) (1.1)

Z(T): Normierung k_{B}: Boltzmann-Konstante.

Verdeutlichen läßt sich dieses physikalische Gesetz anhand folgender Graphik:

Abbildung 1.1: Boltzmann-Verteilung

1.2 Die Arbeitsweise des Simulated Annealing

1.2.1 Darstellung des Grundalgorithmus

Nachdem jetzt der physikalische Grundgedanke des Annealing hinreichend erläutert ist, werde ich detaillierter auf die Algorithmen des Simulated Annealing als Analogon zum physikalischen Prozeß des Abkühlens eingehen. Zur Veranschaulichung betrachte man folgendes Programm, welches den Algorithmus in seiner grundlegenden Form skizziert:

PROCEDURE  SIMULATED_ANNEALING
BEGIN
   INITIALISIERE(Starttemperatur,Startzustand,M)
   t := Starttemperatur
   i := Startzustand
   Stopkriterium := FALSE
   REPEAT
      FOR  z := 1  TO  M  DO BEGIN
         j :=  ERMITTELE_BENACHBARTEN_ZUSTAND(i)
         delta_ij :=  KOSTEN(j) - KOSTEN(i)
         (* Ermittelung der Uebergangswahrscheinlichkeiten *)
         IF (delta_ij <= 0) THEN BEGIN
            i := j
         ELSE BEGIN
            IF (EXP(-delta_ij/t) > RANDOM[0,1[) THEN BEGIN
               i := j
            ENDIF
         ENDIF
      END
      Stopkriterium :=  TEST_STOP(i)
   UNTIL  Stopkriterium
END

Anhand eines Startzustands i_{0} wird eine Folge von Zustandsketten der Länge M, sogenannte Markovketten erzeugt, die eine Konfiguration mit minimalen Kosten ansteuert. Tritt für eine Anzahl von Zustandsketten keine weitere Verbesserung mehr ein, so ist der Zustand minimaler Kosten erreicht. Betrachtet man die innere Zählschleife, welche eine Zustandskette erzeugt, so gehört zu jedem Index k ∈ IN dieser Kette eine Konfiguration i ∈ A. Außerdem ist für alle k ∈ IN und jedes Paar (i, j) ∈ A × A eine sogenannte ”bergangswahrscheinlichkeit“ gegeben. Hier findet man die bertragung der Boltzmann-Verteilung aus dem physikalischen Modell in die Optimierung. Diese gibt in Abhängigkeit vom Kontrollparameter T an, mit welcher Wahrscheinlichkeit man zum Zeitpunkt k vom Zustand i zum Zustand j gelangt. Mathematisch läßt sich die Wahrscheinlichkeit, daß ein bergang mit der Kostendifferenz △c_{ij} = C(j) − C(i) für C(j) ≥ C(i) akzeptiert wird, folgendermaßen formulieren:

P \{△c_{ij} wird akzeptiert\} = exp (−\frac{△c_{ij}}{T} )

Um ein unkontrolliertes ”Springen“ des Algorithmus im Zustandsraum zu verhindern, wird im Algorithmus zusätzlich verlangt, daß von einem Zustand i ∈ A nur zu einem Zustand aus der Menge seiner Nachbarn gewechselt werden darf. Das verlangt die Definition des Begriffs ”Nachbar eines Zustands“. Grundsätzlich ist der Begriff des ”Nachbarn eines Zustands“ in der Literatur nicht eindeutig festgelegt, sondern problemabhängig zu definieren. Für affine n-dimensionale Gitter jedoch, und mit diesen lassen sich schon viele Optimierungsprobleme bearbeiten, liegt es nahe, die 2 · n nächsten Gitterpunkte entlang der Gittergeraden als Nachbarn zu definieren. Dies soll an der Abbildung 1.2 verdeutlicht werden. Eine solche Definition wird beispielsweise von Radzek [26] in seiner Implementierung des Simulated Annealing verwendet. Es gibt jedoch auch andere Definitionen für den Begriff ”Nachbar“. Krimphove [17] beispielsweise definiert den Begriff des ”Nachbarn“ dadurch, daß er für einen möglichen Schritt im Konfigurationenraum alle Zustände zuläßt, die

Abbildung 1.2: Nachbarn auf Gitterpunkten im IR^{2}

entlang der Achsen innerhalb einer festgelegten Maximalschrittweite liegen (siehe Abbildung 1.3). Die dritte in meinem Annealingprogramm verwen-

Abbildung 1.3: Nachbarn im Gitter mit variabler Schrittweite im IR^{2}

dete Nachbarschaftsrelation ermittelt den nächsten Zustand aus der Menge aller Punkte auf einer n–dimensionalen Sphäre, die den aktuellen Zustand umgibt (siehe Abbildung 1.4). Befindet sich das System im Zustand i, so

Abbildung 1.4: Nachbarn auf einer Kugelsphäre im IR^{2}

wird durch ”Würfeln“ eine der möglichen Nachbarkonfigurationen ermittelt.

Sind die Kosten der Nachbarkonfiguration kleiner als die bisherigen Kosten, geht das System sofort in den neu ermittelten Zustand über. Sind die Kosten der Nachbarkonfiguration allerdings höher als die bisherigen Kosten, so wird anhand eines Vergleichs von einem Term ähnlich dem Boltzmann-Term mit einer Zufallszahl ∈ [0, 1] ermittelt, ob in den neuen Zustand gewechselt wird. Man nennt dieses Verfahren ”Metropolis Kriterium“. Es soll die Boltzmann-Verteilung modellieren und wird solange wiederholt, bis sich ein thermisches Gleichgewicht eingestellt hat. In der Regel geschieht dies durch eine vorgegebene Anzahl von Zustandsübergängen. Eine Sequenz dieser Zustandsübergänge nennt man ”Markovkette“. Dieser Begriff wird im nächsten Kapitel noch von entscheidender Bedeutung sein. Danach wird die Temperatur verringert. Der Algorithmus endet bei einer kleinen Temperatur T > 0. Die Verteilungskurve ist dann für △c_{ij} → 0 so steil, daß der Algorithmus nur noch Wechsel zu Zuständen mit geringeren Kosten zuläßt. Das verdeutlicht Abbildung 1.1 unter der Voraussetzung, daß E durch △c_{ij} ersetzt wird. Die Kostenentwickung, die sich während des Simulated Annealing ergibt, ist in Abbildung 1.5 skizziert. Man sieht, daß bei optimalem Verlauf des

Abbildung 1.5: Kostenentwicklung des Simulated Annealing bei fallender Temperatur

Annealings der letzte für T → 0 erreichte Zustand theoretisch mit der Wahrscheinlichkeit 1 der optimale Zustand ist. Auf die theoretischen Grundlagen dieser, bisher vermuteten Aussage wird im nächsten Kapitel noch genauer eingegangen.

1.2.2 Konvergenzbeeinflußung durch Parametermodifikation

In diesem Abschnitt gehe ich genauer auf die Möglichkeiten ein, die zur Verfügung stehen, um die Ausführungszeit einer Implementierung des Simulated Annealing zu verringern. Diese ergeben sich weitestgehend aus dem Grundalgorithmus, werden aber bei der Verwendung des Simulated Annealing häufig vernachlässigt. Es wird sich zeigen, daß sie bei richtiger Anwendung und je nach Problemstellung viel Zeit sparen können. Da der Algorithmus stochastisch arbeitet und auf dasselbe Problem durchaus häufiger angewendet werden muß, gibt es mehrere Ansatzpunkte, an denen man viel Rechenzeit sparen kann, wenn man sich nicht scheut, die Parameter von Durchlauf zu Durchlauf zu ändern.

Starttemperatur

Im Simulated Annealing betrachtet man die Temperatur als Kontrollparameter. Durch Absenken der Temperatur von einem hohen Wert (der ”Schmelztemperatur“), hofft man, daß das System in einem Zustand minimaler Energie ”einfriert“. Mathematisch hat man viel Aufwand getrieben, um diese ”Starttemperatur“ zu ermitteln. Sie ist eng mit dem Mittelwert der Kostenfunktion verbunden und muß hoch genug angesetzt werden, damit anfangs alle Konfigurationen des Systems mit derselben Wahrscheinlichkeit nahe Eins akzeptiert werden. Die optimale Wahl des Startwertes T_{0} ist immer noch nicht geklärt und die Implementationen der theoretischen Ansätze zur Bestimmung von T_{0} sind, da sie meist mit statistischen Argumenten arbeiten, sehr zeitaufwendig. Vielfach wird jedoch übersehen, daß dieser Aufwand nicht bei jedem Algorithmusdurchlauf erneut betrieben werden muß. Arbeitet man mit einer festen Zielfunktion und hat man die Starttemperatur schon einige Male ermittelt, so ist es durchaus legitim, den Algorithmus mit einem Mittelwert aus den vorherigen Starttemperaturen zu starten.

Schrittweite im Gitter

Verwendet man ein konstantes Gitter, so sollte man besonderen Wert auf die Wahl der Schrittweite legen. Bei einer zu groß gewählten Schrittweite besteht die Gefahr, daß einige Minima nicht erreicht werden können. Wählt man die Schrittweite zu klein, so müssen die Markovketten zu lang gewählt werden, um ein brauchbares Ergebnis zu erzielen. Hier besteht die Gefahr, daß die Eigenschaft, ausgehend von jeder Startkonfiguration das globale Minimum erreichen zu können, verloren geht.

Stopkriterium

Während des Arbeitens mit dem Algorithmus kann man verfolgen, wie sinnvoll die Wahl des Abbruchkriteriums und seiner Parameter ist. Wählt man beispielsweise ein Abbruchkriterium, welches die Ergebnisse der letzten 20 Markovketten miteinander vergleicht, so bedeutet dies, daß unter Umständen nutzlos 20 ganze Markovketten durchgerechnet werden, bevor der Algorithmus abbricht. Beobachtet man die Stopketten während der ersten Algorithmusdurchläufe, so kann man mit wenig Mühe geeignetere Werte angeben.

1.2.3 Diskussion

Natürlich ist es möglich, ein Optimierungsproblem immer mit den gleichen Parametern für das Simulated Annealing zu lösen, nur ist der Aufwand dafür im Allgemeinen zu groß. Generell gilt: es ist besser, Parameter des Verfahrens geeignet zu wählen und so überflüssige Berechnungen zu vermeiden. Man sollte sich auch keinesfalls scheuen, offensichtlich zu groß oder klein gewählte Werte in späteren Durchläufen anzupassen. Gerade bei einem recht rechenzeitintensiven Verfahren wie dem Simulated Annealing macht sich das Vorgehen bezahlt und schlägt sich auch bei aufwendigeren Problemen in durchaus brauchbaren Ergebnissen nieder.

1.3 Die Konvergenz des Simulated Annealing

Mathematisch läßt sich das Simulated Annealing als stochastischer Prozeß darstellen. Soll ein Konvergenzbeweis für das Simulated Annealing geführt werden, ist man daher auf Argumente der Wahrscheinlichkeitstheorie angewiesen. Aus den Eigenschaften der in diesem Zusammenhang wesentlichen Begriffe der Markovkette und der stationären Verteilung einer Markovkette resultieren grundlegende Eigenschaften einer jeden Implementierung. Ich möchte das folgende Kapitel daher nicht als alleinstehende mathematische Basis für das Simulated Annealing Verfahren sehen, sondern als Anlaß benutzen, um die Voraussetzungen für eine Implementierung zu klären. Doch zuerst sollen die Grundlagen des Algorithmus dargestellt werden.

1.3.1 Simulated Annealing und die stochastischen Prozesse

Ein stochastischer Prozeß läßt sich wie folgt definieren:

Definition 1.1 (Stochastischer Prozeß) Sei (i_{j})_{j∈J} eine Menge von Zuständen. i_{k} sei der mit dem k. Zufallsexperiment X(k) erzeugte Zustand. Ein stochastischer Prozeß ist eine bedingte Wahrscheinlichkeit der Form:

(1.2) P {X(k + 1)} := P {X(k + 1) = i_{k+1}|X(0) = i_{0},
X(1) = i_{1}, …, X(k) = i_{k}}

Für die allgemeinen stochastischen Prozesse hängt P {k} von einer endlichen Anzahl vorhergehender Zufallsexperimente ab. Um den Simulated Annealing-Algorithmus zu beschreiben, benötigt man jedoch keine so spezielle Definition. Jeder Schritt des Simulated Annealing hängt lediglich vom aktuellen Zustand des Systems ab. So kommen wir von der speziellen Definition der stochastischen Prozesse zu der allgemeineren Beschreibung der Markovprozesse.

Definition 1.2 (Markovprozesse) Ein Markovprozeß ist ein stochastischer Prozeß, der sich wie folgt definieren läßt:

P {X(k + 1)} := P {X(k + 1) = i_{k+1}|X(k) = i_{k}} (1.3)

Auch hierbei ist i_{k} der aus dem k. Zufallsexperiment gewonnene Zustand.

Die bedingte Wahrscheinlichkeit P {X(k + 1)} wird also lediglich vom Zustand i_{k} und dem k+1. Zufallsexperiment X(k + 1) bestimmt. Wie paßt nun das Simulated Annealing in dieses Modell? Der Zustandsübergang von einem Systemzustand zum nächsten unterliegt beim Simulated Annealing einem Zufallsexperiment. Es ist möglich, den Systemübergang von einem Systemzustand X(k) zu einem Systemzustand X(k + 1) mittels einer genau definierten Wahrscheinlichkeit, der Übergangswahrscheinlichkeit P, zu beschreiben. Diese beruht auf der Bolzmann-Ver-teilung und ist, wie bei Markovprozessen nur vom direkt vorhergehenden Zustand X(k) abhängig. Der Zusammenhang zwischen Simulated Annealing und Markovprozessen ergibt sich damit offensichtlich. Man erhält einen Konvergenzbeweis für das Simulated Annealing aus berlegungen zu aneinander gereihten Markovprozessen, sogenannten Markovketten.

Definition 1.3 (Markovketten) Eine Markovkette ist ein stochastischer Prozeß P {X(k), k ∈ Z^{+}} mit abzählbarem Zustandsraum (I_{j})_{j∈J} und den Eigenschaften:

P {X(0) = i_{0}, X(1) = i_{2}, …, X(k) = i_{k}} > 0 (1.4) P {X(k + 1) = i_{k+1}|X(k) = i_{k}} ∀i_{1}, …, i_{k+1} ∈ I (1.5)

Man sieht: Jede Folge von Zustandsübergängen des Simulated Annealing, setzt man eine feste Temperatur voraus, entspricht einer Kette von Markovprozessen. Es gilt nun herauszufinden, unter welchen Bedingungen aufeinander folgende Markovketten in abzählbaren Zustandsräumen konvergieren und wie sich diese Konvergenz auf das Simulated Annealing übertragen läßt.

Der physikalische Gleichgewichtszustand

Das entscheidende Element eines Konvergenzbeweises besteht darin nachzuweisen, daß tatsächlich für jede Temperatur T des Simulated Annealing eine sogenannte stationäre Verteilung erreicht wird. Das bedeutet, daß es für jede Markovkette einen Vektor q gibt, der folgende Eigenschaft erfüllt:

q_{i} = lim_{k→∞} P (X(k) = i) , i ∈ A. (1.6)

A bezeichnet hierbei eine abzählbare Menge von Zuständen. Dieses entspricht in Analogie zum physikalischen Prozeß des Abkühlens dem Erreichen eines Zustands minimaler Energie. Man spricht in diesem Zusammenhang auch vom ”Gleichgewichtszustand“ des Systems bei einer Temperatur T. Der Vektor q repräsentiert alle Zustände des Systems und die Wahrscheinlichkeit, daß sie angenommen werden. Für T → 0 nimmt die Wahrscheinlichkeit ab, daß energetisch schlechtere Zustände angenommen werden. Kann man die Existenz eines Gleichgewichtszustands für alle T erreichen, so erhält man für T → 0 eine Folge von Gleichgewichtszuständen, die gegen den Systemzustand minimaler Energie konvergiert. An dieser Stelle zeigt sich ein erstes Problem bei der Implementierung des Simulated Annealing: Wie lang müssen die Markovketten der Zufallsexperimente bei den Zustandsübergängen gewählt werden, damit eine existierende Gleichgewichtsverteilung genügend gut approximiert wird? Nur bei einer genügend guten Approximation der stationären Verteilung ist eine Konvergenz gegen die globalen Minima gesichert. Es wird sich im Folgenden herausstellen, daß dieses fast nie der Fall ist. Je schlechter die stationäre Verteilung approximiert wird (das ist besonders bei sehr rechenintensiven Problemen mit vielen Variablen fast immer der Fall), desto größer ist die Wahrscheinlichkeit, daß der Algorithmus die Umgebung des nächsten Minimums nicht mehr verläßt. Die Tatsache, daß dennoch auch für diese Probleme ein globales Minimum erreichbar bleibt, macht den Algorithmus auch bei solchen Problemen noch interessant.

Der mathematische Gleichgewichtszustand

Unter welchen Bedingungen läßt sich nun für das Simulated Annealing zeigen, daß zu jeder Temperatur T ein Gleichgewichtszustand erreicht wird? Um die Markovketten des Simulated Annealing mathematisch beschreiben zu können, führt man die Begriffe der Generierungsmatrix, der Akzeptanzmatrix und der bergangsmatrix ein. Mit diesen Mitteln läßt sich jeder Zustandsübergang im Simulated Annealing als Markovprozeß darstellen und man erhält über eine Grenzwertbetrachtung die Gleichgewichtsverteilung. Es sei A eine abzählbare Menge von Zuständen und i, j ∈ A seien Zustände in A.

Definition 1.4 (Generierungsmatrix) G_{ij}(T_{k}) gibt die Wahrscheinlichkeit an, mit der man als möglichen bergang denjenigen von i nach j wählt. Die Matrix G_{ij}(T_{k}) gibt also an, daß j zum Zeitpunkt T_{k} eine zu i benachbarte Konfiguration ist. G heißt Generierungsmatrix.

In den meisten Implementierungen des Simulated Annealing sind die Erzeugungswahrscheinlichkeiten der Nachbarkonfigurationen gleichverteilt. Für mögliche Zustandsübergänge von i nach j gilt dann:

g_{ij} = \frac{1}{|Anzahl der Nachbarn von i|} > 0 ∀ i ∈ A.

Es gibt auch Beispiele für nicht gleichverteilte Erzeugungswahrscheinlichkeiten. Diese findet man z.B. bei Aarts und van Laarhoven [2] oder bei Radzek [26]. Ich gehe jedoch nicht darauf ein, da sie für die weiteren berlegungen nicht von Bedeutung sind.

Definition 1.5 (Akzeptanzmatrix) A_{ij}(T_{k}) definiert die Wahrscheinlichkeit, mit der dieser mögliche bergang akzeptiert wird. A_{ij}(T_{k}) beinhaltet den ausgewerteten Boltzmann-Ausdruck für alle möglichen bergänge i nach j, das heißt:

A_{ij}(T_{k}) = exp (−\frac{f(j) − f(i)}{T_{k}}) .

A heißt Akzeptanzmatrix. Die einzelnen Elemente der Matrix heißen Akzeptanzwahrscheinlichkeiten.

Mit diesen Definitionen ist es nun möglich, die Übergangswahrscheinlichkeiten des Algorithmus folgendermaßen zu beschreiben:

Definition 1.6 (bergangsmatrix) ∀i, j ∈ A sei :

P_{ij}(T_{k}) := \{ \begin{matrix}G_{ij}(T_{k}) · A_{ij}(T_{k}) ∀i ≠ j \\ 1 − ∑_{l∈A,l≠i} P_{ij}(T_{k}) ∀i = j \end{matrix} (1.7)

Dann heißt P_{ij}(T_{k}) Übergangsmatrix. Die einzelnen Elemente der Matrix werden Übergangswahrscheinlichkeiten genannt.

Mit diesen Definitionen kann man ein Kriterium für die Konvergenz des Simulated Annealing Verfahren formulieren. Ein genauer Beweis des Theorems läßt sich bei Aarts und van Laarhoven [2] nachlesen.

1.3.2 Konvergenztheorem

Theorem 1.1 (Existenz einer stationären Verteilung) Sei A = (i_{j})_{j∈J} eine abzählbare Menge von Zuständen und C : A → IR eine Kostenfunktion. Angenommen, die Matrizen G_{ij}(T_{k}) und A_{ij}(T_{k}) genügen den folgenden Bedingungen:

1. ∀T_{k} > 0 ∀i, j ∈ A ∃p ≥ 1 ∃λ_{1}, …, λ_{p} ∈ A mit i = λ_{1} und j = λ_{p} : G_{λ_{l}λ_{l+1}}(T_{k}) > 0 ∀l = 0, …, p − 1 Das heißt: alle Zustände müssen erreichbar bleiben.

2. ∀T_{k} > 0 ∀i, j ∈ A : G_{ij}(T_{k}) = G_{ji}(T_{k})

G_{ij}(T_{k}) ist symmetrisch:

Ist j ein Nachbar von i, so ist auch i ein Nachbar von j.

3. ∀T_{k} > 0 ∀i, j, l ∈ A :

C(i) ≤ C(j) ≤ C(l) ⇒ A_{il}(T_{k}) = A_{ij}(T_{k}) · A_{jl}(T_{k})

Verschiedene Wege zwischen zwei Zuständen haben die gleiche Akzeptanzwahrscheinlichkeit.

4. ∀T_{k} > 0 ∀i, j ∈ A : C(i) ≥ C(j) ⇒ A_{ij}(T_{k}) = 1 Ein bergang, bei dem die Kosten abnehmen oder gleich bleiben, wird mit der Wahrscheinlichkeit 1 akzeptiert.

5. ∀T_{k} > 0 ∀i, j ∈ A : C(i) < C(j) ⇒ 0 < A_{ij}(T_{k}) < 1 Ein bergang, bei dem die Kosten zunehmen, wird mit einer Wahrscheinlichkeit größer als 0, aber kleiner als 1 akzeptiert.

6. ∀i, j ∈ A : C(i) < C(j) ⇒ lim_{T_{k}→0} A_{ij}(T_{k}) = 0 Mit abnehmender ”Temperatur“ nimmt auch die Wahrscheinlichkeit, daß ein bergang zu höheren Kosten akzeptiert wird, immer weiter ab und erreicht für den Grenzübergang T_{k} → 0 ebenfalls 0.

Dann gilt

i. lim_{l→∞} P_{T_{k}}\{X(l) = i\} =: q_{i}(T_{k}) = _{\frac{A_{i_{0}i}(T_{k})}{∑_{j∈A} A_{i_{0}j}(T_{k})}}

Man erhält einen Ausdruck für die stationäre Verteilung und

ii. lim_{T_{k}→0} q_{i}(T_{k}) = \{ \begin{matrix}|\frac{1}{A_{opt}} | für i ∈ A_{opt} \\ 0 sonst \end{matrix}

Für T → 0 ist die Wahrscheinlichkeit, einem optimalen Zustand zu erreichen, gleich 1 und jedes globale Optimum wird gleichwahrscheinlich angenommen.

1.3.3 Diskussion

So erstaunlich die Aussage des Theorems auch ist, so schwierig gestaltet sich numerisch dessen Realisierung. Optimale Konvergenz wird nur für unendlich lange Markovketten und unendlich langsames Absenken der Temperatur erreicht. Diese Eigenschaften werden jedoch von keinem Simulated-Annealing-Programm erreicht. Was eine Implementierung des Algorithmus dennoch rechtfertigt, ist eine Eigenschaft, die erst beim zweiten Hinsehen aus dem Konvergenztheorem hervorgeht: Erfüllt eine Implementierung die Voraussetzungen des Theorems, so wird bei abbrechender Markovkette zwar nicht mehr ein globales Minimum mit der Wahrscheinlichkeit 1 gefunden, aber man erreicht bei hinreichend langer Markovkette ein globales Minimum mit einer Wahrscheinlichkeit, die deutlich größer als 0 ist. Man besitzt also ein Verfahren, welches unabhängig vom Startwert in der Lage ist, ein globales Minimum zu finden. Je mehr Mühe man in die Approximation der stationären Verteilung investiert, um so mehr nimmt die Wahrscheinlichkeit zu, eines der globalen Minima zu erreichen. Stationäre Verteilungen zu zwei nur geringfügig unterschiedlichen Temperaturen unterscheiden sich nur wenig. Deshalb wirkt sich ein vorsichtiges Herabsetzen der Temperatur, hat man einmal eine stationäre Verteilung approximiert, positiv auf die Approximation der nächsten stationären Verteilung aus. Die Wahl dieser beiden Parameter sowie die Wahl einer ausreichend hohen Starttemperatur beeinflussen direkt die Qualität der erzeugten Ergebnisse. Die Einschränkung, daß der zugrundeliegende Konfigurationenraum abzählbar diskret sein muß, wird sich als nicht wesentlich herausstellen, da eine interne Zahlendarstellung im Rechner den Konfigurationenraum immer auf ein zwar feines, aber diskretes und sogar endliches Gitter reduziert. Es ist also nicht notwendig, die benachbarten Zustände nur auf einem ndimensionalen Gitter mit fester Schrittweite zu suchen, wie es Radzek [26] in seinem Simulated-Annealing-Programm realisiert. Andere Methoden, wie zum Beispiel eine Methode mit variabler Schrittweite, die bei Krimphove [17] vorgestellt wird, oder eine Methode mit konstanter Schrittweite, bei der der nächstmögliche Zustand auf der Sphäre einer n-dimensionalen Kugel um den aktuellen Zustand liegt, sind bei geeigneter Parameteranpassung ebenso erlaubte Implementierungen des Simulated Annealing. Es wird sich herausstellen, daß sie je nach Problem sogar die besseren Ergebnisse liefern.

1.4 Anwendung des Simulated Annealing auf globale Optimierungsprobleme mit nichtdiskretem Konfigurationenraum

In diesem Abschnitt gehe ich genauer auf die verschiedenen Implementierungen des Simulated Annealing ein. Zugrunde liegen die beiden Algorithmen von Radzek [26] und Krimphove [17]. Der Algorithmus von Krimphove hat sich für kleine Optimierungsprobleme mit wenigen Variablen bewährt, ist jedoch aufgrund der immensen Zahl von Funktionsauswertungen bei größeren Problemen nicht brauchbar. Um größer dimensionierte Probleme bearbeiten zu können, eignet sich der Algorithmus von Radzek. Da er jedoch einige Schwächen zeigte, was die Art der Zustandsübergänge sowie die Qualität der Approximation des gefundenen Optimums betrifft, wurde dieser Algorithmus von mir zu einer dritten Form der Simulated-Annealing-Implementation weiter entwickelt, die sich für die Probleme der Kristallstruktur-Vorhersage sowie für einige Testfunktionen ausgezahlt hat.

1.4.1 Die Algorithmen in SIMULA

Die Optimierungsroutine SIMULA umfaßt zwei unterschiedliche Arten von Annealing-Routinen.

Der Algorithmus von Radzek

Die erste Annealing-Routine geht auf eine Diplomarbeit von Radzek [26] zurück und stellt eine Simulated-Annealing-Routine mit konstanter Schrittweite dar, die in Pascal implementiert ist. Die Stärken des Algorithmus liegen in seiner relativ geringen Laufzeit, was ihn auch auf komplexere Probleme anwendbar macht. Nachteilig wirkt sich jedoch die Tatsache aus, daß ein fest vorgegebenes, sehr grobes Gitter nicht verlassen wird. Wird der Gitterabstand zu klein gewählt, verliert der Algorithmus seine günstige Laufzeit, da die Markovketten sehr lang gewählt werden müssen, um zu garantieren, daß jeder Zustand im Zustandsraum innerhalb der Markovkette erreichbar bleibt. Um den Nachteil des sehr groben Gitters zu kompensieren, wird durch eine

Abbildung 1.6: Algorithmus von Radzek

Abschlußverfeinerung eine bessere Approximation des gefundenen Minimums vorgenommen. Der Algorithmusablauf des Simulated Annealing ohne Verfeinerung wird in dem Struktogramm 1.6 wiedergeben.

Der Algorithmus von Krimphove

Die zweite Annealing-Routine stützt sich auf eine objektorientierte Implementierung von Krimphove [17] in der Programmiersprache Smalltalk. In dieser Routine wird versucht, das Problem des fest vorgegebenen Gitters mittels einer variablen Schrittweitensteuerung, die auf Corona et al. [8] zurückgeht, zu beheben. Leider gestaltet sich eine angemessene Schrittweitensteuerung meist als sehr aufwendig, auch wenn man in diesem Fall eine Abschlußverfeinerung vermeiden kann. Das macht die Routine nur auf Probleme mit verhältnismäßig wenig Variablen anwendbar. In diesen Fällen hat sich das Verfahren allerdings bewährt. Das Struktogramm 1.7 gibt den Verlauf des Simulated Annealing mit variabler Schrittweite wieder.

Abbildung 1.7: Algorithmus von Krimphove

Die Schrittweitensteuerung im Algorithmus von Krimphove

Schritte in alle Koordinatenrichtungen sind häufig nicht gleich günstig. Für Laufrichtungen mit hoher Akzeptanzwahrscheinlichkeit bietet es sich an, programmgesteuert größere Schritte zu unternehmen. Anders ist das Verfahren bei weniger erfolgversprechenden Richtungen. Bei Schritten in diese Richtung ist Vorsicht geboten. Die Schrittweite wird herabgesetzt. Eine Schrittweitensteuerungsroutine versucht, diese Idee umzusetzen. Die Methode ist einfach: Um lokale Zusatzinformationen über den Kurvenverlauf der Zielfunktion zu gewinnen, testet die Routine über eine vorgegebene Zahl von Durchläufen für alle Koordinatenrichtungen, ob ein Schritt von zufälliger Länge in diese Richtung akzeptiert wird. Je nach Prozentsatz der akzeptierten Schritte wird die Schrittlänge in diese Richtung entweder verlängert, verkürzt oder gleich belassen. Krimphove schlägt vor, dabei folgendermaßen vorzugehen:

Satz 1.1 (Schrittweitensteuerung) Für jede Koordinatenrichtung wird in der folgenden Weise ein Multiplikator v_{aktuell} ermittelt:

v_{i_{aktuell}} := \{ \begin{matrix}v_{i_{vorher}} · (1 + c_{h} · ^{\frac{\frac{n_{akzep}}{n_{durchl}}}{0.4}} ) , \\ wenn mehr als 60 % der getesteten \\ Übergänge akzeptiert wurden. \\ v_{i_{vorher}} · (1 + c_{h} · ^{\frac{\frac{−n_{akzep}}{n_{durchl}}}{0.4}} )^{−1} , \\ wenn weniger als 40 % der getesteten \\ Übergänge akzeptiert wurden. \\ v_{i_{vorher}} sonst. \end{matrix}

Hierbei sind:

v_{i_{aktuell}} : aktueller Schrittweitenmultiplikator
v_{i_{vorher}} : letzter Schrittweitenmultiplikator
n_{akzep} : Anzahl der akzeptierten Schritte in diese Richtung
n_{durchl} : Gesamtzahl der versuchten Schritte in diese Richtung
c_{h} : problemabhängig gewählte Konstante

Ein versuchter Zustandsübergang in diese Richtung erfolgt dann auf die folgende Weise:

x_{neu} := x_{alt} + RANDOM() · v_{i_{aktuell}} · e_{i}

mit:

x_{neu} : neu ermittelter Zustand
x_{alt} : bisheriger Zustand
RANDOM() : erzeugt Zufallszahl im Intervall [0,1]
e_{i} : i-ter Einheitsvektor

Die Größe c_{h} ist hierbei nicht festgelegt, sondern problemabhängig zu wählen. Mit ihr läßt sich steuern, wie schnell der Schrittweitenmultiplikator größer oder kleiner wird. Die Wahl sollte davon abhängig gemacht werden, wie groß das Intervall zwischen maximaler und minimaler Schrittweite ist.

Beispiel:

• Für c_{h} = 1 wird die Schrittweite

– verdoppelt (falls mehr als 60 % der getesteten bergänge akzeptiert wurden).

– halbiert (falls weniger als 40 % der getesteten bergänge akzeptiert wurden).

• Für c_{h} = 2 wird die Schrittweite

– verdreifacht (falls mehr als 60 % der getesteten bergänge akzeptiert wurden).

– gedrittelt (falls weniger als 40 % der getesteten bergänge akzeptiert wurden).

Der modifizierte Algorithmus

Eine Auswertung der Zielfunktion nimmt im Falle der Kristallstruktur-Vor-hersage viel Zeit in Anspruch. Deshalb zeigte sich schnell, daß der Algorithmus mit variabler Schrittweite nach Krimphove auf diese Probleme nicht anwendbar war. Leider lieferte auch die Implementation von Radzek nur ausreichende Ergebnisse. Der Grund dafür war schnell gefunden: Da der Algorithmus den aktuellen Zustand nur auf einem festen Gitter entlang der Achsen veränderte, konnte das Optimum aufgrund der Beschaffenheit der Zielfunktion nicht erreicht werden. Beispiel im IR^{3} (siehe dazu Abbildung 1.8): Der optimale Zustand liegt

Abbildung 1.8: Beispiel im IR^{3}:

auf einem Gitterpunkt, der aktuelle Zustand befindet sich auf einem abseits liegenden Gitterpunkt. Die direkte Verbindungsgerade der beiden Punkte verläuft diagonal zum Gitter. Auf den Kanten ausgehend vom aktuellen Zustand steigt der Wert der Zielfunktion an. Man erhält lokale Maxima im Gitter, die keine echten Maxima darstellen. Um das Optimum zu erreichen, wäre ein Weg durch die Mitte der Raute ideal. Um dieses zu erreichen, wurde der Algorithmus von Radzek so abgeändert, daß der nächstmögliche Zustand auf einer n-dimensionalen Sphäre um den aktuellen Zustand gewählt wird. Damit erhält man nun einen Algorithmus, der wie der Algorithmus von Krimphove nicht mehr auf ein festes Gitter fixiert ist und trotzdem mit in etwa der gleichen Anzahl von Funktionsauswertungen auskommt wie der Algorithmus von Radzek (siehe dazu auch die Tabelle Seite 36, welche die Laufzeiten der verschiedenen Algorithmen für das Beispiel CsCl wiedergibt). Das Problem der lokalen Maxima bei konstantem Gitter wird im Kapitel über die Kristallstruktur-Vorhersage noch genauer erläutert.

1.4.2 Erläuterung der Abschlußverfeinerung

Benutzt man einen Algorithmus mit variabler Schrittweite, so ist man in den meisten Fällen nicht auf eine Verfeinerung der Ergebnisse des Simulated-Annealing-Algorithmus angewiesen, da ein Minimum bis auf Minimalschrittweite approximiert wird. Nachteilig wirkt sich jedoch der Aufwand aus, der für eine Schrittweitensteuerung betrieben wird. Die Lösung von vielen Optimierungsproblemen erfordert deshalb das Arbeiten mit fester Schrittweite. Ist diese geeignet gewählt, kann man sicher sein, sich in einer akzeptablen Umgebung um das Minimum zu befinden. Eine befriedigende Approximation der Lösung wird dieses Ergebnis jedoch nicht darstellen. Dieser Aspekt wurde bisher leider fast völlig vernachlässigt. Der einzige Ansatz zu einer Abschlußverfeinerung findet sich bei Radzek. Hierbei handelt es sich jedoch nur um ein Verfahren, das mit einer Temperatur nahe 0 den Simulated Annealing Algorithmus nochmals durchläuft. Ist das Abbruchkriterium erfüllt (und das dauert aufgrund der vielen Markovketten, die berechnet werden, meist sehr lange), wird die Schrittweite herabgesetzt. Was sich ergibt, ist ein sehr zeitaufwendiges, ”straight forward“ programmiertes stochastisches Abstiegsverfahren. Hier besteht aufgrund der aufwendigen Zielfunktion der Kristallstruktur-Vorhersage der Bedarf, für Abhilfe zu sorgen. Da das Problem der Abschlußverfeinerung ein lokales Optimierungsproblem darstellt, bieten sich mehrere bekannte und erprobte Verfahrensweisen an.

Die einfachen Abstiegsverfahren

Einfache Abstiegsverfahren sind leicht zu implementieren und lassen sich auch auf nicht differenzierbare Probleme anwenden. Nachteil: Die Rechenzeit ist zu groß. Dennoch haben sie sich bei ersten Rechnungen mit der Kristallstruktur-Vorhersage bewährt. Ausgehend vom aktuellen Zustand werden sukzessive alle Koordinaten mit einer festen Schrittweite in Abstiegsrichtung verändert. Ist mit dieser Schrittweite kein Abstieg mehr möglich, so wird die Schrittweite herabgesetzt.

Die Quasi-Newton-Verfahren

Die Gruppe der Quasi-Newton-Verfahren ist eine der effizientesten Methoden zur Lösung von lokalen Optimierungsproblemen. Leider versagen sie bei den Problemen der Kristallstruktur-Vorhersage, da sie die Anforderung an die Zielfunktion stellen, zweimal stetig differenzierbar zu sein. Dies ist im Allgemeinen nicht der Fall. Dennoch haben sich die Verfahren bei der Bearbeitung der Testfunktionen bewährt. Die Quasi-Newton-Verfahren realisieren die Suche nach dem Minimum unter Verwendung von Krümmungsinformationen. Ausgehend vom aktuellen Zustand x^{k} wird der nächste Punkt x^{k+1} durch einen Line-Search nach der folgenden Vorschrift ermittelt:

Satz 1.2 (Die Quasi-Newton Vorschrift)

x^{k+1} := x^{k} − λ^{k} · M_{k}∇F (x^{k}) (1.8)

λ wird durch einen Line-Search so bestimmt, daß gilt:

F (x^{k+1}) ≤ F (x^{k}). (1.9)

F stellt die dem Optimierungsproblem zugrundeliegende Zielfunktion dar und M_{k} approximiert die inverse Hessematrix im k-ten Schritt.

Da beim Start des Verfahrens noch keine Informationen über die Kurvenkrümmung von F vorliegen, begnügt man sich in der Regel damit, daß man für M_{k} die Einheitsmatrix wählt. In diesem Schritt entspricht das Verfahren dem Verfahren des steilsten Abstiegs. Im weiteren Verlauf erfolgt in jedem Schritt nach einer festen Formel eine Veränderung von M_{k}. Diese Veränderung, auch ”Update“ genannt, arbeitet die im Iterationsschritt neu gewonnenen Krümmungsinformationen in die Matrix M_{k} mit ein. Es gibt verschiedene Methoden, dieses ”Update“ durchzuführen. Jedes Update definiert ein eigenes Verfahren. Gemeinsam gilt jedoch für alle Verfahren:

M_{k} muß in jedem Schritt positiv definit sein.

Als erfolgsversprechendste Verfahren haben sich die Verfahren mit Rang-2-Korrekturmatrix erwiesen. Die beiden wesentlichsten Vertreter dieser Gattung von Quasi-Newton-Verfahren sind das Verfahren von Davidon, Fletcher und Powel (DFP-Verfahren) und das Verfahren von Broyden, Fletcher, Goldfarb und Shanno (BFGS-Verfahren). Beide Verfahren werden ausführlich bei Schittkowsky [28] hergeleitet und behandelt. Ich habe mich für das bei Krabs [16] vorgestellte DFP-Verfahren entschieden.

Satz 1.3 (Die DFP Update-Formel) Die Update-Formel für das DFP-Verfahren lautet wie folgt:

M_{k+1}^{DF P} := M_{k}^{DF P} + \frac{v^{k}·v^{k^{T}}}{v^{k^{T}} ·w^{k}} − \frac{(M_{k}·w^{k})·(M_{k}·w^{k})^{T}}{w^{k^{T}} ·M_{k}·w^{k}} (1.10)

mit
v^{k} := x^{k+1} − x^{k}
w^{k} := ∇F (x^{k+1}) − ∇F (x^{k})

Sollte der Bedarf bestehen, ein BFGS-Verfahren zu implementieren, so läßt sich die Implementation dadurch realisieren, daß die Updatefunktion der inversen Hessematrix neu implementiert wird. Das Quasi-Newton-Verfahren an sich muß nicht geändert werden.

Die Verfahren des steilsten Abstiegs

Die Verfahren des steilsten Abstiegs sind die am häufigsten verwendeten Verfahren in der Optimierung. Sie sind einfach zu implementieren, stabil, und stellen wenige Anforderungen an die Zielfunktion. Zur Realisierung benötigt man lediglich die Existenz des Gradienten der Zielfunktion. Da sich diese Verfahren ohne großen Aufwand auch für die im Allgemeinen nicht differenzierbaren Funktionen der Kristallstruktur-Vorhersage modifizieren ließen, erwiesen sie sich als geeignete Kandidaten zur Bearbeitung dieser Funktionen.

Ein Verfahren des steilsten Abstiegs bestimmt im aktuellen Punkt den negativen Gradienten der Zielfunktion. Dieser kennzeichnet die Richtung des steilsten Abstiegs der Funktion. In dieser Richtung sucht man mit einem Line-Search nach dem Minimum und startet das Verfahren dann neu. Mathematisch stellt sich das Verfahren für einfach stetig differenzierbare Funktionen F dann folgendermaßen dar:

Satz 1.4 (Die Vorschrift des steilsten Abstiegs)

x^{k+1} := x^{k} − λ∇F (x^{k}) (1.11)

Hierbei wird λ ≥ 0 in jedem Schritt so gewählt, daß G(λ) := F (x^{k}−λ∇F (x^{k})) minimiert wird.

Leider besteht aber nicht immer die Möglichkeit, Differenzierbarkeit für die Zielfunktion F (x) vorauszusetzen, was aber nicht sofort eine Benutzung der Methode des steilsten Abstiegs ausschließt. Für stetige und stückweise differenzierbare Funktionen lassen sich Abstiegsrichtungen dadurch errechnen, daß man eine ”Schätzung“ für den Gradienten angibt: Man variiert den aktuellen Zustand x^{k} für jede Komponente um ±ϵ (ϵ ist eine kleine reelle Zahl) und berechnet den Zielfunktionswert. Würde man nun für eine stetig differenzierbare Funktion den Gradienten durch den Mittelwert KAPITEL 1. DAS OPTIMIERUNGSVERFAHREN 24

aus links– und rechtsseitigen Differenzquotienten ermitteln, so benutzt man hier nur die Abstiegsrichtung zur Ermittelung des Gradienten:

\frac{∂F (x^{k})}{∂x_{i}} := \{ \begin{matrix}\frac{F (x^{k}+ϵ·e_{i})−F (x^{k})}{ϵ} , falls F (x^{k} + ϵ · e_{i}) < F (x^{k}) \\ \frac{F (x^{k})−F (x^{k}−ϵ·e_{i})}{ϵ} sonst. \end{matrix} (1.12)

Man erhält für den negativen ”geschätzten“ Gradienten auf jeden Fall eine Abstiegstiegsrichtung. Für differenzierbare Teilstücke wird sich die Qualität des Verfahrens kaum von der Qualität des echten Verfahrens des steilsten Abstiegs unterscheiden. Für lediglich stetige und im Allgemeinen nicht differenzierbare Zielfunktionen liefert es im Gegensatz zum echten Verfahren jedoch die gewünschten Ergebnisse.

1.4.3 Beispiel für die Anwendung einer Abschlußverfeinerung

Nach der Ermittelung der Starttemperatur befindet sich der Kristall (in diesem Fall Halit,NaCl) im ungeordneten Zustand (siehe Abbildung 1.9).

Abbildung 1.9: Die Halitstruktur zum Zeitpunkt hoher Temperatur

**********************************************************************
*   HOMOGENER SIMULATED ANNEALING ALGORITHMUS !!!                    *
**********************************************************************
Eingabedatei                     : .\nacl.ann
Anzahl der Variablen             :   24
Komponenten des Vektors zwischen  0.000 und  1.000
Die Startschrittweite war        : 0.1000
Bewegungen im Gitter             : auf 3 dimensionalen Kugelsphaeren
                                   mit konstanten Radien
**********************************************************************
Laenge einer Markovkette         :    20
Bestimmung des Startwerts der Temperatur T nach  fester Vorgabe
Startwert war                    : 800.00000000
Absenkung der Temperatur durch Multiplikation mit einer Konstanten
Konstante                        : 0.950
Bestimmung Stopkriteriums nach Kirkpatrick
**********************************************************************
Startvektor                          :
 x(  1)=   0.00000 x(  2)=   0.00000 x(  3)=   0.00000
 x(  4)=   0.20000 x(  5)=   0.80000 x(  6)=   0.30000
 x(  7)=   0.70000 x(  8)=   0.50000 x(  9)=   0.70000
 x( 10)=   0.30000 x( 11)=   0.30000 x( 12)=   0.50000
 x( 13)=   0.70000 x( 14)=   0.60000 x( 15)=   0.20000
 x( 16)=   0.40000 x( 17)=   0.90000 x( 18)=   0.80000
 x( 19)=   0.60000 x( 20)=   0.80000 x( 21)=   0.10000
 x( 22)=   0.90000 x( 23)=   0.80000 x( 24)=   0.60000
**********************************************************************

KAPITEL 1. DAS OPTIMIERUNGSVERFAHREN 25

Nach dem Durchlauf der Optimierung ist eine gewisse Nähe zum Minimum, das heißt zur optimalen Struktur, erreicht (siehe Abbildung 1.10). Um diesen Abstand noch weiter zu verringern und die Struktur noch besser zu approximieren, wendet man eine Abschlußverfeinerung an.

Abbildung 1.10: Die Halitstruktur nach dem Annealing

**********************************************************************
Stopwert von t war               : 0.0468354817638843
**********************************************************************
Funktionswert Simulated Annealing: 733.95947067
Funktionsauswertungen S. A.      : 3801
Anzahl der Zustandsaenderungen   : 680
**********************************************************************
Ergebnisvektor Sim. Annealing    :
 x(  1)=   0.00000 x(  2)=   0.00000 x(  3)=   0.00000
 x(  4)=   0.51630 x(  5)=   0.05409 x(  6)=   0.47989
 x(  7)=   0.44998 x(  8)=   0.53201 x(  9)=   0.00106
 x( 10)=   0.98801 x( 11)=   0.56105 x( 12)=   0.54208
 x( 13)=   0.02104 x( 14)=   0.96515 x( 15)=   0.53013
 x( 16)=   0.48019 x( 17)=   0.52106 x( 18)=   0.47994
 x( 19)=   0.49028 x( 20)=   0.05501 x( 21)=   0.03017
 x( 22)=   0.02103 x( 23)=   0.47899 x( 24)=   0.02017
**********************************************************************

Nach der Abschlußverfeinerung hat man sich dem Idealkristall soweit genähert, daß im Bild keine Unregelmäßigkeit mehr festzustellen ist. Das Bild des Halitkristalls in Abbildung 1.11 entspricht dem Idealkristall.

Abbildung 1.11: Die Halitstruktur nach der Abschlußverfeinerung

**********************************************************************
Anzahl der Funktionsauswertungen : 5944
Anzahl der Zustandsaenderungen   : 680
**********************************************************************

     Zusammensetzung der Zielfunktion :
     Modifizierte Zielfunktion :
     Term[1] = elektroterm=  Wichtung  :  1.00000
     Term[2] = valenzterm+  Wichtung  :  1.00000
     **********************************************************************
     Der optimale Funktionswert ist   : 627.69048035
     Valenzterm                       :   0.00105889
     **********************************************************************
     Der Ergebnisvektor lautet        :
      x(  1)=   0.00000 x(  2)=   0.00000 x(  3)=   0.00000
      x(  4)=   0.50000 x(  5)=   0.00009 x(  6)=   0.49989
      x(  7)=   0.49998 x(  8)=   0.50001 x(  9)=  -0.00006
      x( 10)=   0.00000 x( 11)=   0.50005 x( 12)=   0.50008
      x( 13)=   0.00004 x( 14)=   0.00015 x( 15)=   0.50003
      x( 16)=   0.50009 x( 17)=   0.50006 x( 18)=   0.49994
      x( 19)=   0.49998 x( 20)=   0.00001 x( 21)=  -0.00007
      x( 22)=  -0.00003 x( 23)=   0.49999 x( 24)=   0.00017
     **********************************************************************

1.4.4 Die Parameter der Simulated-Annealing-Algo-rithmen

Nach der Erläuterung der verschiedenen zur Verfügung stehenden Algorithmen soll nun genauer auf die Steuerungsparameter der Algorithmen eingegangen werden. Der Simulated Annealing Algorithmus ”lebt“ davon, daß diese Parameter richtig gesetzt werden. Beispielsweise bedeutet eine zu hoch angesetzte Starttemperatur einen berhang an nutzlos vergeudeter Rechenzeit. Wählt man jedoch den Parameter ”Starttemperatur“ zu niedrig, so wird der Algorithmus sich nicht von einem einfachen Abstiegsverfahren unterscheiden. Er findet immer das ”nächste“ Minimum. Ebenso wichtig ist die richtige Wahl der Länge der Markovkette und des Stopkriteriums. Falsch gewählte Werte verlängern unnötig die Rechenzeit oder beschneiden den Algorithmus so, daß dieser nicht mehr in der Lage ist, mit einer akzeptablen Wahrscheinlichkeit brauchbare Ergebnisse zu liefern. Im Folgenden werde ich nun einige Methoden vorstellen, welche die Parameterwahl vereinfachen sollen. Sie wurden von Radzek zusammengetragen und von mir ergänzt. Alle diese Ansätze verwenden relativ einsichtige statistische Argumente und haben sich bewährt. Um die manchmal sehr umfangreichen Formeln etwas transparenter darzustellen, werden sie ausführlich kommentiert. Die eigentliche Wahl der Parameter sollte sich jedoch nicht nur auf die Anwendung der Formeln verlassen, sondern vor allem auch eigene Beobachtungen des Algorithmusverhaltens mit berücksichtigen.

Startkriterien

Das Startkriterium von Kirkpatrick. Das Startkriterium von Kirkpatrick et al. [14] ist das einfachste der zur Verfügung stehenden Kriterien. Es stützt sich auf die berlegung, daß die bergangswahrscheinlichkeit nach relativ kurzer Zeit den Wert 1 erreicht, wenn man ausgehend von einem kleinen Wert > 0 die Temperatur sukzessive verdoppelt. Zu jeder Temperatur berechnet man eine Anzahl von Zustandsübergängen. Wird ein bestimmter Prozentsatz von berechneten Zustandsübergängen, die sogenannte Akzeptanzrate, überschritten, so hat man die Starttemperatur bestimmt. Leider überschreitet man mit dieser berlegung häufig gerade bei hohen Temperaturen die ausreichende Starttemperatur deutlich. Das schlägt sich dann später,während des Abkühlprozesses, in der Anzahl der Funktionsauswertungen nieder. Trotzdem ist das Startkriterium von Kirkpatrik immer dann empfehlenswert, wenn man das erste Mal mit einer neuen Zielfunktion arbeitet. Die Wahrscheinlichkeit, den richtigen Startwert zu unterbieten, ist sehr gering.

Das Startkriterium von Johnson et al. Johnson et al. [12] zeigten, daß man das langwierige Verdoppeln der Temperatur im Kriterium von Kirkpatrick et al. [14] einsparen kann, indem man für eine Menge von Nachbarn nur die durchschnittliche Erhöhung der Kosten ermittelt. Berechnet man deren mittlere Akzeptanzwahrscheinlichkeit, so läßt sich dieser Wert dazu verwenden, einen Startwert für die Temperatur zu berechnen. Ein Defizit des Kriteriums besteht darin, daß zu einer repräsentativen Vorhersage der Starttemperatur bei einem sehr großen Zustandsraum auch sehr viele Zustandsauswertungen nötig sind. Wie bei Kirkpatrik wird eine bestimmte Akzeptanzrate χ_{0} nahe 1 vorgegeben. Hiermit läßt sich eine Starttemperatur berechnen, wenn die Wahrscheinlichkeit, daß ein Zustand mit mittelhohen Kosten akzeptiert wird, den Schwellenwert χ_{0} erreicht. So werden genügend benachbarte Zustandsübergänge akzeptiert. Die Starttemperatur ist ermittelt.

Satz 1.5 (Johnson et. al.) Sei im Folgenden für die Zustände i und j und die Kostenfunktion F: C(j) ≥ C(i). Die Anzahl der ermittelten Zustandsübergänge sei mit m bezeichnet. Dann läßt sich das Kriterium von Johnson folgendermaßen darstellen:

Sei \overline{△C} := \frac{1}{m} · ∑_{(i,j)}(C(j) − C(i))

Ein Prozentsatz χ_{0} von Zustandsübergängen wird akzeptiert, wenn

χ_{0} = e^{\frac{−\overline{△C}}{T_{0}}}

⇒ T_{0} = _{\frac{\overline{△C}}{ln(\frac{1}{χ_{0}} )}}. (1.13)

Das Startkriterium von White et al. Sucht man ein Kriterium für eine geeignete Starttemperatur, so wird man fast immer auf den Ausdruck der Standardabweichung stoßen. Auch White et al. [30] verwenden dieses Argument in seiner Bestimmung des Startkriteriums. Eine geeignete Starttemperatur ist für alle Probleme T = ∞. Das ist jedoch nicht realisierbar. Bestimmt man daher eine niedrigere Starttemperatur, so sollte die mittlere Abweichung der Akzeptanz eines Zustandsübergangs jedoch nicht bedeutend von der Akzeptanz desselben bergangs bei T = ∞ abweichen. Dazu wird im Folgenden die stationäre Verteilung q(T) zu einer Temperatur T vorausgesetzt.

Die stationäre Verteilung beeinhaltet für alle möglichen Zustände die Wahrscheinlichkeiten ihres Auftretens im Gleichgewichtszustand. Soll jetzt eine Prognose über die in der stationären Verteilung erwarteten mittleren Kosten in einer diskreten Menge A abgeben werden, so ist dies beispielsweise mit dem folgenden Ausdruck möglich.

(1.14) Erw(C(T)) := \frac{1}{|A|} · ∑_{i∈A} C(i) · q_{i}(T)

Für T = ∞ gilt sofort:

Erw(C(∞)) := lim_{T →∞} Erw(C(T)).

Im Allgemeinen gilt, da mit abnehmender Temperatur immer weniger Zustandsübergänge zu Zuständen mit schlechteren Kosten akzeptiert werden:

Erw(C(∞)) ≥ Erw(C(T)) für T < ∞. (1.15)

Das gilt also auch insbesondere für T_{0}. T_{0} soll jetzt so groß gewählt werden, daß die Differenz von Erw(C(∞)) und Erw(C(T_{0})) kleiner ist als die Standardabweichung σ(∞). Mathematisch bedeutet das:

σ(∞) ≥ Erw(C(∞)) − Erw(C(T)). (1.16)

White et al. [30] zeigten nun, daß für große T gilt:

Erw(C(T)) ≈ Erw(C(∞)) − \frac{σ^{2}(∞)}{T}. (1.17)

Aus den Formeln (1.16) und (1.17) ergibt sich für T_{0}:

T_{0} ≥ σ(∞) (1.18)

und man erhält mit der Standardformel für die Varianz das Startkriterium von White et al. [30].

Satz 1.6 (Das Startkriterium von White et. al.)

T_{0} ≥ σ(∞) = \sqrt{Erw(C^{2}(∞)) − (Erw(C(∞)))^{2}} (1.19)

Erw(C^{2}(∞)) und Erw^{2}(C(∞)) lassen sich nun über die Kostenentwicklung einer Markovkette bei einer genügend großen Temperatur T_{∞} approximativ ermitteln. Das Interessante an diesem Kriterium ist die Tatsache, daß zur Ermittelung von T_{0} nur die Berechnung einer Markovkette notwendig ist. Das gestaltet sich aber immer dann problematisch, wenn die Markovketten relativ kurz sind. Soll dieses Verfahren verwandt werden, so empfielt es sich immer, während der Ermittelung der Starttemperatur die Länge der Markovkette zu vergrößern. Nur so erhält man repräsentative Werte für die Kostenvarianz.

Die fest vorgegebene Starttemperatur. Ist T_{0} in mehreren Versuchen mit einer festen Zielfunktion bestimmt worden, so daß man Erfahrung mit der Wahl von T_{0} besitzt, sollte man für weitere Tests die Starttemperatur fest vorgeben. Es spart erheblich an Rechenaufwand.

Abkühltabellen

Die für den Abstieg der Temperatur berechneten Werte werden als ”Abkühltabellen“ bezeichnet. Ebenso wie bei der Bestimmung des Startwertes basieren die meisten Kriterien für die Abkühltabelle mehr auf Empirie als auf theoretischen berlegungen. Im anschließenden Abschnitt wird genauer auf vier unterschiedliche Abkühltabellen eingegangen. In der Simulated-Annealing-Prozedur betrachtet man die Temperatur als Kontrollparameter. In Analogie zum physikalischen Prozeß des Abkühlens von Festkörpern hofft man durch langsames Absenken der Temperatur ab einem ”hohen“ Wert, der sogenannten ”Schmelze“, daß das System in einem Zustand minimaler Energie ”einfriert“. Die Suche nach einer angemessenen Abkühltabelle, die möglichst effizient das Simulated Annealing bei der Suche nach einem globalen Minimum unterstützt, wird in vielen Veröffentlichungen behandelt. Problematisch gestaltet sich dabei die Eigenschaft des Simulated Annealing, daß bei zu schnellem Senken der Temperatur das System in einem ungünstigen, lokalen Minimum einfrieren kann. Der Grund für dieses Verhalten liegt nahe. Je schneller die Absenkung der Temperatur durchgeführt wird, um so länger muß die Markovkette zu jeder Temperatur gewählt werden, damit das thermische Gleichgewicht genügend gut approximiert wird. Dieses ist jedoch, so zeigt der Konvergenzbeweis für das Simulated Annealing, Grundvoraussetzung für die Bestimmung eines globalen Minimums. Ein mathematischer Ansatz zur Bestimmung der Abkühltabelle nutzt die Voraussetzung, daß stationäre Verteilungen ”nahe“ beieinander liegen sollen, etwa folgendermaßen:

Vorüberlegung 1.1 Eine stationäre Verteilung q(T_{k+1}) liegt ”nahe“ bei der vorhergehenden Verteilung q(T_{k}), wenn für eine kleine vorgegebene Konstante δ > 0 gilt:

1 + δ ≥ \frac{q_{i}(T_{k})}{q_{i}(T_{k+1})} ≥ \frac{1}{1 + δ} (1.20)

Kennt man nun die Kostendifferenz aller Zustände zu den optimalen Kosten, so kann man ein Kriterium für die Absenkung der Temperatur angeben:

Satz 1.7 (Hilfssatz) Sei △c_{i} := C(i) − C_{opt}. Gelten die Bedingungen:

1. T_{k+1} ≥ \frac{T_{k}}{1 + \frac{ln(1+δ)·T_{k}}{△c_{i}}} ∀i ∈ A, k ∈ IN_{0}

2. T_{k+1} ≤ T_{k} ∀k ∈ IN_{0}

So folgt daraus, daß

1 + δ ≥ \frac{q_{i}(T_{k})}{q_{i}(T_{k+1})} ≥ \frac{1}{1 + δ}

Beweis (Hilfssatz): Seien k ∈ IN_{0} und i ∈ A beliebig.

Sei Ann(T_{k}) := ∑_{i∈A} e^{− \frac{△c_{i}}{T_{k}}} (1.21)

Durch Umformung von 1. folgt:

(1 + δ) · e^{− \frac{△c_{i}}{T_{k+1}}} ≥ e^{− \frac{△c_{i}}{T_{k}}}

⇒ Ann(T_{k+1}) ≥ \frac{1}{1 + δ} · Ann(T_{k})

Offensichtlich gilt für T_{k} ≥ T_{k+1} : Ann(T_{k+1}) ≥ Ann(T_{k})

⇒ \frac{q_{i}(T_{k})}{q_{i}(T_{k+1})} ≤ 1 + δ.

Da nach Aussage i. des Konvergenztheorems (1.1) gilt:

q_{i}(T_{k}) = e^{− \frac{△c_{i}}{T_{k}}} .

Ferner ergibt sich aus der Definition von Ann(T_{k}) (siehe (1.21)):

\frac{q_{i}(T_{k})}{q_{i}(T_{k+1})} ≥ \frac{e^{− \frac{△c_{i}}{T_{k}}}}{e^{− \frac{△c_{i}}{T_{k+1}}}} · \frac{1}{1 + δ} ≥ \frac{1}{1 + δ} (1.22)

und man erhält die Behauptung.

Wählt man also eine Absenkungsbedingung, die sich aufgrund der Voraussetzungen des Hilfssatzes 1.7 ergibt, so kann man sicher gehen, daß sich die stationären Verteilungen zu aufeinanderfolgenden Temperaturen nur geringfügig unterscheiden. Das entscheidende Problem der berlegung offenbart sich jedoch bei der genaueren Betrachtung der Formeln: Möchte man eine Absenkungsbedingung angeben, so ist man darauf angewiesen, einen Wert für den Ausdruck △c_{i} anzugeben. Hierfür gibt es verschiedene Ansätze, die zu unterschiedlichen Absenkungsbedingungen führen.

Die konstante Absenkung. Die einfachste Methode der Absenkung ist die Verwendung einer fest vorgegebenen Konstanten, die, an die vorherige Temperatur ”anmultipliziert“, die nächste Temperatur ergibt.

(1.23) T_{neu} := c · T_{alt} mit c geeignet gewählt

Obwohl dieses Kriterium sehr grob ist, lassen sich schon ansehnliche Ergebnisse erzielen. Die Größe der Konstanten wird in der Literatur oft mit einem Wert > 0.85 angegeben. Ich selbst habe gute Erfahrungen mit c ≈ 0.95 gemacht.

Die Absenkung nach Lundy und Mees. Das Absenkungskriterium von Lundy und Mees [19] verwendet direkt das Ergebnis der Vorüberlegung 1.1. Eine Approximation für △c_{i} erhält man, indem man:

△c_{i} := C_{max} − C_{min} (1.24)

setzt. C_{max} und C_{min} sind hier die maximalen bzw. minimalen bisher bestimmten Kosten. Diese Abschätzung ist sicherlich sehr grob. Viele theoretische berlegungen versuchen deshalb, für das Problem mit anderen Mitteln eine bessere Abschätzung zu erhalten. Mit der eingangs erwähnten Vorüberlegung 1.1 ergibt sich sofort:

Satz 1.8 (Absenkung nach Lundy und Mees)

T_{k+1} := T_{k} · \frac{1}{1 + T_{k} · \frac{ln(1+δ)}{U}} (1.25)

mit

U := C_{max} − C_{min}

Je besser der Wert von U, also die Approximation von △c_{i}, den tatsächlichen Wert von △c_{i} erreicht, desto weniger weichen die stationären Verteilungen von zwei aufeinanderfolgenden Temperaturen voneinander ab. Da mit zunehmender Laufzeit des Programms immer mehr Zustände untersucht werden, nimmt erst nach und nach die statistische Wahrscheinlichkeit zu, daß dies der Fall ist. Es ist also eine zu lange Vorlaufzeit nötig, um brauchbare Absenkungsgeschwindigkeiten zu berechnen. Das hingegen spricht gegen eine Anwendung des Kriteriums, da gerade ein zu schnelles Absenken im hohen Temperaturbereich das Simulated Annealing empfindlich stört.

Die Absenkung nach Aarts und van Laarhoven. Aarts und van Laarhoven [1] umgehen das Problem, eine Approximation für △c_{i} anzugeben, auf die folgende Art und Weise: Für große Werte von T ist zu erwarten, daß alle Zustände mit in etwa der gleichen Wahrscheinlichkeit akzeptiert werden. Daher ist △c_{i} ebenfalls gleichverteilt. Sei der Mittelwert der △c_{i} zur Temperatur T mit μ(T) bezeichnet und sei σ^{2}(T) die zugehörige Varianz. Für große Werte von T leitet man mit Mitteln der Statistik eine Aussage über die wahrscheinliche Größe von △c_{i} her:

P (△c_{i} ≤ μ(T) + 3 · σ^{2}(T)) ≈ 0.998 (1.26)

Kennt man die nötigen Größen μ(T) und σ^{2}(T), so kann man sich relativ sicher sein, eine geeignete obere Schranke für △c_{i} angeben zu können. Eine Abkühlvorschrift ergäbe sich dann mit der Vorüberlegung folgendermaßen:

T_{k+1} := \frac{T_{k}}{1 + \frac{ln(1+δ)·T_{k}}{μ(T_{k})+3·σ^{2}(T_{k})}} (1.27)

σ^{2}(T_{k}) läßt sich einfach dadurch approximieren, daß man die Kostenvarianz während einer Markovkette bei einer hinreichend hohen Temperatur betrachtet. Die Länge einer Markovkette wird im Folgenden mit LMK angegeben:

(1.28) σ^{2}(T_{k}) ≈ \frac{1}{LMK + 1} LMK ; ∑ ; i:=0 C^{2}(X_{i}) − (\frac{1}{LMK + 1} LMK ; ∑ ; i:=0 C(X_{i}))^{2}

X_{i} ist hierbei der i. Zustand der Kette

Da jedoch C_{min} nicht bekannt ist, gestaltet sich die Bestimmung von μ(T) als problematisch. Um dieses Problem zu umgehen greifen Aarts und van Laarhoven auf den freien Parameter δ zurück. Wählt man δ klein genug, so kann man μ(T_{k}) vernachlässigen. Aarts und van Laarhoven zeigten weiter, daß auch für kleine T die Wahrscheinlichkeit für eine günstige obere Schranke mit ≈ 0.95 noch sehr hoch liegt. Es ergibt sich das Kriterium:

Satz 1.9 (Abkühlkriterium von Aarts und van Laarhoven)

T_{k+1} := \frac{1}{1 + \frac{ln(1+δ)·T_{k}}{3·σ(T_{k})}} (1.29)

Und somit erhält man einen deutlich besseren Ansatz als bei Lundy und Mees.

Die Absenkung nach Huang. Huang [11] benutzt einen etwas anderen Ansatz als Lundy und Mees sowie Aarts und van Laarhoven. Dieser basiert direkt auf dem eigentlichen Problem der Temperaturabsenkung und fordert, daß die Differenz der Erwartungswerte der Kosten zweier aufeinanderfolgender Markovketten kleiner als eine vorgegebene Standardabweichung ist. Ist dies der Fall, so differieren die stationären Verteilungen nur wenig und man hat eine geeignete Temperaturabsenkung gefunden. Mathematisch ausgedrückt erhält man:

(1.30) Erw(C(T_{k})) − Erw(C(T_{k+1})) = λ · σ(T_{k}) mit λ ≤ 1

Ohne das genauer auf die Details der Umformungen eingegangen werden soll, sei bemerkt, daß man zeigen kann:

\frac{∂Erw(C(T))}{∂ ln(T)} = \frac{σ^{2}(T)}{T} (1.31)

Die genauen Umformungen findet man bei Huang [11]. Löst man die partielle Ableitung approximativ durch einen Differenzenquotienten auf, so erhält man für T_{k} nahe bei T_{k+1}:

(1.32) \frac{σ^{2}(T_{k})}{T_{k}} = \frac{Erw(C(T_{k+1})) − Erw(C(T_{k}))}{ln(T_{k+1}) − ln(T_{k})}

Verwendet man zusätzlich den Ausdruck (1.30), so ergibt sich die Absenkungsbedingung:

Satz 1.10 (Absenkungsbedingung nach Huang) Für λ ≤ 1 gilt:

T_{k+1} = T_{k} · e^{(\frac{−λ·T_{k}}{σ(T_{k})} )}. (1.33)

Die Bestimmung von σ^{2}(T_{k}) wird dabei in der gleichen Art und Weise durchgeführt, wie im Ausdruck (1.28).

Stopkriterien

Die Stopkriterien werden häufig auch als ”Gefrierpunkt– Bedingungen“ bezeichnet. Sie bestimmen die kleinste Temperatur, bei welcher das Simulated Annealing versucht, ein thermisches Gleichgewicht zu erreichen. Zu diesem Zeitpunkt ist die Wahrscheinlichkeit, in einen Zustand mit höheren Kosten zu wechseln, so gering, daß sich das Simulated Annealing in ein einfaches Abstiegsverfahren verwandelt hat. Man hat sich dem letzten Minimum genähert und eine Verbesserung der Kosten läßt sich nun, aufgrund der vielen nicht mehr akzeptierten Zustandsübergänge, kaum noch erzielen. Das System, in dem optimiert werden soll, ist ”eingefroren“. Sollen noch weitere Verbesserungen erzielt werden, so kommt eine lokale Optimierungsstrategie sinnvoller zur Anwendung als die rechenintensive ”Trial and Error“ Methode des Simulated Annealing. Die effektivsten Verfahren zur Bestimmung der Stoptemperatur basieren auf der Beobachtung der Kostenfunktion während des Simulated Annealing. Im Folgenden wird auf einige von ihnen eingegangen werden.

Das Stopkriterium von Kirkpatrick et al. Das Stopkriterium von Kirkpatrick et al. [14] basiert auf der berlegung, daß man nur dann entscheiden kann, ob ein Stopwert für die Temperatur erreicht ist, wenn man lange genug beobachtet hat, daß die erreichten Kosten nicht mehr variieren. Für jede Temperatur werden die Kosten des am Schluß der Markovkette erreichten Zustands in eine Warteschlange, die sogenannten Stopkette, eingetragen. Unterscheiden sich eine Anzahl der letzten Stopkettenelemente um einen Wert kleiner als eine vorgegebene Konstante von den aktuellen Kosten, so ist die Stoptemperatur erreicht.

Satz 1.11 (Stopkriterium von Kirkpatrick et al.) Seien k die aktuellen Kosten und δ eine kleine reelle Zahl ≥ 0. STOPKETTE enthalte die abschließenden Kosten von n aufeinanderfolgenden Markovketten. Eine Stoptemperatur ist dann erreicht, wenn gilt:

(1.34) |k − STOPKETTE[i]| ≤ δ ∀ i = 1, …, n

Das Stopkriterium von Huang. Der Ansatz von Huang unterscheidet sich nicht wesentlich vom letzten Ansatz. Huang vergleicht die einzelnen Elemente der Stopkette miteinander. Unterscheiden sich diese nur unwesentlich, so ist eine Stoptemperatur erreicht.

Satz 1.12 (Stopkriterium von Huang) Sei δ eine kleine reelle Zahl ≥ 0. STOPKETTE enthalte die abschließenden Kosten von n aufeinanderfolgenden Markovketten. Eine Stoptemperatur ist dann erreicht, wenn gilt:

(1.35) |STOPKETTE[i] − STOPKETTE[j]| ≤ δ ∀ i, j = 1, …, n mit i ≠ j

Das Stopkriterium von Lundy und Mees. Die bisher vorgestellten Verfahren reichen in der Regel durchaus, um Ergebnisse zu produzieren, die sich mit denen der zwei folgenden, etwas theoretischeren Verfahren messen können. Ihr entscheidender Nachteil liegt auf der Hand. Es müssen immer n überflüssige Markovketten berechnet werden, bis man entscheiden kann, ob eine Stoptemperatur erreicht ist. An diesem Punkt setzen die Verfahren von Lundy und Mees sowie Aarts und van Laarhoven an. In beiden Verfahren basiert die Theorie auf folgender berlegung:Ist das Annealing optimal durchgeführt und hat die Temperatur T einen kleinen Wert erreicht, so ist die Wahrscheinlichkeit, daß die letzte Konfiguration der zugehörigen Markovkette sich um mehr als ϵ vom globalen Minimum unterscheidet, kleiner als ein kleines, vorgegebenes δ. Mathematisch ausgedrückt muß gelten:

∑_{\begin{matrix}C(i)>C_{opt}+ϵ \\ i∈A \end{matrix}} q_{i}(T_{S}) = P (X(k) = i ∧ C(i) > C_{opt} + ϵ|T_{S} = T) < δ (1.36)

Benutzt man nun Aussage i. des Konvergenztheorems 1.1 , so erhält man:

∑_{\begin{matrix}C(i)>C_{opt}+ϵ \\ i∈A \end{matrix}} \frac{e^{\frac{−C(i)+C_{opt}}{T_{S}}}}{∑_{j∈A} e^{\frac{−C(j)+C_{opt}}{T_{S}}}} < δ

Schätzt man nun C(i) − C(j) durch ϵ ab und verwendet folgende Formel, die für genügend kleine T_{S} gilt:

∑_{j∈A} e^{\frac{−C(j)+C_{opt}}{T_{S}}} ≥ 1

So erhält man:

∑_{\begin{matrix}C(i)>C_{opt}+ϵ \\ i∈A \end{matrix}} \frac{e^{\frac{−C(i)+C_{opt}}{T_{S}}}}{∑_{j∈A} e^{\frac{−C(j)+C_{opt}}{T_{S}}}} < ∑_{\begin{matrix}C(i)>C_{opt}+ϵ \\ i∈A \end{matrix}} e^{\frac{−ϵ}{T_{S}}} ≤ (|A| − 1) · e^{\frac{−ϵ}{T_{S}}} (1.37)

Und somit ergibt sich das Kriterium:

Satz 1.13 (Stopkriterium von Lundy und Mees) Eine geeignete Stoptemperatur T_{S} ist dann erreicht, wenn:

(|A| − 1) · e^{\frac{−ϵ}{T_{S}}} < δ ⇔ T_{S} ≤ \frac{ϵ}{ln(|A| − 1) − ln δ} (1.38)

Das Stopkriterium von Aarts und van Laarhoven. Der Ansatz von Aarts und van Laarhoven ist ähnlich. Sie bestimmen zu jedem Wert T_{k} den Mittelwert der Kosten der Markovkette. Man faßt den Mittelwert als Funktion

¯ : IR ⟶ IR auf.

Der Stopwert T_{S} ist dann erreicht, wenn C̄(T_{k}) − C_{opt} hinreichend klein ist. Es ergibt sich durch eine Taylorentwicklung mit Entwicklungspunkt T_{k} :

C̄(0) ≈ C̄(T_{k}) − \frac{∂C̄(T)}{∂T}T |_{T =T_{k}}

Für T = 0 sind allerdings nach Aussage ii. des Konvergenztheorems 1.1 nur noch globale Minima möglich. Daher gilt:

C̄(0) = ∑_{i∈A_{opt}} \frac{C_{opt}}{|A_{opt}|} = C_{opt}

Insgesamt erhält man daraus für kleine T_{k}:

C̄(T_{k}) − C_{opt} ≈ \frac{1}{T_{k}} (\overline{C^{2}(T_{k})} − (C̄(T_{k}))^{2}) (1.39)

C^{2} ist hier Problematisch ist hierbei wieder, daß keine Information über C_{opt} verfügbar ist. Aarts und van Laarhoven schätzen nun für geeignetes ϵ C̄(T_{S}) − C_{opt} llϵ · C̄(T_{0}) und erhalten damit das Kriterium:

Satz 1.14 (Stopkriterium von Aarts und van Laarhoven)

\frac{\overline{C^{2}(T_{S})} − (C̄(T_{S})^{2}}{T_{S} · C̄(T_{0})} ≤ ϵ (1.40)

1.4.5 Länge der Markovketten

So groß die Auswahl an Kriterien für Startwerte und Stopbedingungen ist, eine fundierte Angabe über die Länge einer Markovkette findet man kaum. Alle Angaben, die gemacht werden, stützen sich ausschließlich auf Beobachtungen bei eigenen Experimenten. Die Schätzungen gehen von Werten, die an die Anzahl der Variablen gekoppelt sind, bis zu Werten, die durch die bisher ermittelten Maximal– bzw. Minimalkosten (Romeo, Sangiovanni-Vincentelli, [27]) bestimmt werden. Ich habe die besten Erfahrungen mit fest vorgegebenen Längen für Markovketten gemacht. Die Probleme der Kristall-Struktur-Vorhersage reagierten ingesamt nicht besonders empfindlich auf nderungen der Länge der Markovketten, wenn diese nicht allzu kurz gewählt werden. Man sieht leicht, daß für eine größere Zahl von Variablen auch längere Markovketten benötigt werden, um ein thermisches Gleichgewicht zu erreichen. Senkt man den Kontrollparameter ”Temperatur“ vorsichtig genug, so kann man die Länge der Markovketten durch einen Wert proportional zur Anzahl der Variablen angeben. Dieses Kriterium findet vor allem bei Anwendungen des Simulated Annealing in der Physik und Kristallographie immer wieder Anklang. Wie man später sieht, haben sich bei der hier programmierten Kristallstruktur-Vorhersage-Routine sogar Werte im Bereich der einfachen Anzahl von Variablen bewährt.

1.5 Das Maß für die Qualität von Optimierungsverfahren: die Testfunktionen

Gerade bei Algorithmen wie dem Simulated-Annealing ist es wichtig, eine Methode zu finden, um die Qualität der Ausgabe des Algorithmus, also die Genauigkeit der erlangten Ergebnisse, zu beurteilen. Nur so ist man in der Lage, Eigenschaften des Verfahrens richtig zu erkennen und auf diese Weise zu entscheiden, wie sinnvoll eine Anwendung des Simulated Annealing-Algorithmus auf ein bestimmtes Problem ist.

1.5.1 Die Testfunktionen als Qualitätsurteil für das Simulated Annealing

Die Tatsache, daß das Simulated-Annealing-Verfahren ein stochastisches Verfahren ist, bewirkt, daß trotz gleicher Parameterwahl und gleicher Anfangswerte durchaus nicht immer dieselben Resultate zu erzielen sind. Durch die eingebaute Zufallsrichtung führt der Algorithmus unter Umständen immer wieder zu unterschiedlichen Ergebnissen. Nur im optimalen Fall ist gewährleistet, daß ein globales Minimum mit der Wahrscheinlichkeit 1 erreicht wird. Bei Zielfunktionen, die nur wenige lokale Minima besitzen, läßt sich durch häufiges Anwenden des Algorithmus die Qualität noch heuristisch bestimmen. Viele Ausdrücke besitzen jedoch eine sehr viel größere Anzahl von lokalen Minima. Sind diese ferner sogar fast gleichwertig, wie das globale Minimum, so scheint dieses Unterfangen aussichtslos. Um trotzdem ein Verfahren zu besitzen, das die Qualität der Ausgabe eines Algorithmus wiederspiegelt, haben Törn und Zilinskas [29] in ihrem Buch eine Auswahl von Testfunktionen zusammengetragen, deren Aussehen und deren Minima bekannt sind. Sie dienen auch in anderen Bereichen der Optimierung zur Beurteilung von Algorithmen und sollen verschiedene Zielfunktionen modellieren. Liefert die Implementation des Simulated Annealing für eine Testfunktion mit einem bestimmten Parametersatz häufig ”gute“ Ergebnisse, so ist die Wahrscheinlichkeit groß, daß für ähnliche Zielfunktionen genauso gute Ergebnisse zu erzielen sind. Von Interesse ist also, wie sensibel das Verfahren auf geschätzte Parameter und unterschiedliche Zielfunktionen reagiert.

1.5.2 Beispiele für Zielfunktionen

Die Rosenbrock-Funktion

RB_{n} (n=2,4,7,10,13)

RB_{n} = ∑_{k=2}^{n} (100 · (x_{k} − x_{k−1}^{2})^{2} + (1 − x_{k−1})^{2}) (1.41)

min{RB_{n}(x)| − 100 ≤ x_{i} ≤ 100, i = 1, …, n} = (1, …, 1) ∈ IR^{n} (1.42)

RB_{n}(1, …, 1) = 0

Die Rosenbrock-Funktion (siehe Abbildung 1.12) hat bloß ein Minimum,

Abbildung 1.12: Rosenbrock-Funktion (n=2)

ist aber in einer Umgebung um das Minimum sehr flach. Ein brauchbarer Optimierungsalgorithmus sollte in der Lage sein, dieses zu finden. Das Simulated Annealing lieferte auch für alle Parametersätze das gewünschte Ergebnis. Einzig die Laufzeit differierte je nach Startkriterium und Algorithmusart gehörig.

Konstantes n-dimensionale variable
Gitter  Sphäre  Schritte
Kirkpatrick  7756  7941  14971
White  2764  2804  7570
Johnson  3970  4115  11967

Mittlere Laufzeiten in Abhängigkeit von den Startkriterien (30 Durchläufe). Für die Schrittweitenbestimmung wurden 10 mal 10 Schrittweitentests pro

Richtung durchgeführt.

Die Rastringin-Funktion

R (n=2)

R(x_{1}, x_{2}) = 10(x_{1}^{2} + x_{2}^{2} − cos(18 · x_{1}) − cos(18 · x_{2})) (1.43) min{R(x_{1}, x_{2})| − 2 ≤ x_{i} ≤ 2, i = 1, 2} = (0, 0) R(0, 0) = −20 (1.44)

R(x_{1}, x_{2}) (Abbildung 1.13) hat im zulässigen Bereich noch etwa weitere

Abbildung 1.13: Rastringin-Funktion

50 lokale Minima. An der Rastringin-Funktion lassen sich sehr schön die unterschiedlichen Eigenschaften des Simulated Annealing demonstrieren. Als Startkriterium wurde eine feste Starttemperatur vorgegeben, um in etwa gleiche Voraussetzungen für die Verfahren zu schaffen. Als Beispiel wurden als Absenkungsbedingung das Verfahren der konstanten Absenkung und die Absenkungsbedingung von Aarts und van Laarhoven untersucht. Es ergibt sich die folgende Tabelle.

Konstante Absenkung Absenkung nach Aarts und
(k=0,85)  v. Laarhoven (λ=0,01)
max. Funktionswert  4,21  0,58
min. Funktionswert  -20  -20
ø Funktionswert  -9,34  -15,76

Absenkungsbetrachtung bei der Rastringin-Funktion (jeweils 30 Durchläufe, konstante Schritte auf n-dim. Sphären).

Ein ähnliches Ergebnis ergab sich für die verschiedenen Startkriterien und für den Fall, daß Markovketten zu kurz gewählt wurden. Man sieht: Die Qualität des Ergebnisses des Annealing-Prozesses hängt direkt mit der Wahl der Parameterwahl für den Abkühlprozeß zusammen. Es ist jedoch trotz schlechter Parameterwahl immer noch eine Restwahrscheinlichkeit vorhanden, mit der ein globales Minimum gefunden wird.

Kirkpatrick White Johnson
χ_{0}=0,95  chi_{0}=0,95
ø Startwert  128,44  55,15  194,1
maximaler Funktionswert  10,27  21,17  -7,89
minimaler Funktionswert  -18,789 -13,94 -18,789
mittlerer Funktionswert -14,55 0,71 -15,330

Qualität des Ergebnisses in Abhängigkeit vom Startwert (Rastringin Funktion, 30 Durchläufe, konstante Absenkung mit k=0,95).

Die Grienwank-Funktion

GN_{n} (n=2,10)

GN_{n} = ∑_{i=1}^{n} \frac{x_{i}^{2}}{d} − ∏_{i=1}^{n} cos(\frac{x_{i}}{\sqrt{i}}) + 1 (1.45)

min{GN_{n}| − 100 ≤ x_{i} ≤ 100, i = 1, 2} = (0, …, 0) (1.46)

F ür n = 2 , d = 200 , −100 ≤ x_{i} ≤ 100 für i = 1, 2 hat (1.47)

GN etwa 500 weitere Minima im zulässigen Bereich.

F ür n = 10 , d = 4000 , −600 ≤ x_{i} ≤ 600 und i = 1, …, 10 hat (1.48)

GN_{n} einige tausend weitere Minima im zulässigen Bereich.

Die Grienwank-Funktion (Abbildung 1.14) ist ein wahrer Test für jede Optimierungsstrategie. An dieser Funktion sieht man sehr schön, daß ein Verfahren wie das Simulated Annealing für die Optimierung zwar ein Fortschritt, aber kein Allheilmittel für alle Probleme ist. Es ist mir nicht gelungen, das globale Minimum mit dem Verfahren zu bestimmen. Die nächsten Tabelle zeigt jedoch, daß der Algorithmus mit variabler Schrittweite für solche Probleme die erheblich höhere Rechenzeit mit deutlich besseren Ergebnissen belohnt.

Konstantes n-dimensionale variable
Gitter  Sphäre  Schritte
Startwert  800  800  800
maximaler Funktionswert  64,63  64,63  9,22
minimaler Funktionswert  53,71  47,94  -2,83
mittlerer Funktionswert  60,67  58,18  3,23

Erlangte Ergebnisse bei der Grienwank Funktion (30 Durchläufe, Schrittweitenbestimmung mit 10· 10 Tests pro Richtung).

KAPITEL 1. DAS OPTIMIERUNGSVERFAHREN 40

Abbildung 1.14: Grienwang Funktion (eingeschränkter Bereich) Abbildung 1.15: Grienwang Funktion (gesamter Bereich)

1.5.3 Diskussion

Als Ergebnis dieses Abschnitts kann festgehalten werden: Testfunktionen bieten zwar nicht die letzte Sicherheit bei der Beurteilung der Qualität von Algorithmen, doch dienen sie als nützliches Hilfsmittel, um zu entscheiden, wie sinnvoll die Anwendung eines speziellen Verfahrens auf eine Funktion mit ähnlicher Topologie wie bei einer der Zielfunktionen ist. Des weiteren lassen sich über die Testfunktionen Parametersätze ausloten. Dies spart viel Zeit im Umgang mit zeitaufwendig zu berechnenden Zielfunktionen. Bei vernünftigem Umgang bietet dies eine gewisse Sicherheit für die Anwendung von stochastischen Algorithmen auf Optimierungsprobleme. Eine umfangreichere Aufzählung von Testfunktionen ist bei Törn und Zilinskas [29] zu finden.

Simulated Annealing als Mittel zur Kristallstruktur-Vorhersage

2.1 Kristalle als Objekte der Mathematik

Nach der Erläuterung des Simulated Annealing-Verfahrens soll es im folgenden für eine Vorhersage von Kristallstrukturen eingesetzt werden. Um ein besseres Verständnis für die untersuchten Objekte zu erreichen, werde ich im ersten Abschnitt einen berblick über das Kristallmodell geben. Es gibt mehrere Ansätze, Kristallstrukturen in diesem Modell zu beschreiben. Die im zweiten Abschnitt vorgestellten berlegungen basieren auf einem Artikel von Pannetier et al. [25]. Der dort vorgestellte Ansatz stützt sich auf empirische Regeln. Er wurde von mir verbessert und erweitert. In dem hier entwickelten Programm ANNEAL sind diese Modifikationen implementiert. ANNEAL stützt sich auf die Optimierungsroutine SIMULA, deren Verwendung im Anhang genauer beschrieben wird. Im letzten Abschnitt wird auf die erzeugten Strukturen genauer eingegangen und gezeigt, daß es mit dem erweiterten Ansatz gelingt, eine größere Klasse von Kristallstrukturen vorherzusagen. Doch zunächst zum Kristallmodell.

2.1.1 Eigenschaften kristalliner Festkörper

Kristalline Festkörper sind wesentliche Bestandteile der Natur. Viele chemische Verbindungen treten in kristalliner Form auf. Kennt man ihren inneren Aufbau, so kann man die physikalischen Eigenschaften der Kristalle gezielt nutzen. Ein Effekt, den man in vielen Bereichen der Technik verwendet, ist der sogenannte ”Piezoelektrische Effekt“. Er entsteht immer dann, wenn man in bestimmte Richtungen Druck auf den Kristall ausübt. Es liegt nahe, daß der innere Aufbau von Kristallen sich grundlegend vom inneren Aufbau anderer Festkörper unterscheidet.

Die kristallinen Festkörper

Kristalle zeichnen sich dadurch aus, daß sie aus einem dreidimensionalen Gerüst ihrer Bausteine (Gruppen von Atomen, Molekülen oder Ionen) bestehen. Man bezeichnet diese Anordnungen als Raumgitter oder als Translationsgitter. Dabei bilden die Gitterbausteine äquidistante Punktreihen oder Liniengitter. Diese liegen in festen Abständen nebeneinander und bilden Netzebenen oder Flächengitter, die ihrerseits wieder in gleichen Abständen angeordnet sind und die Raumgitter definieren. Die Abstände der Gitterbausteine können mit dem Verfahren der röntgenographischen Untersuchung, der sogenannten Röntgendifraktometrie, bestimmt werden. Durchstrahlt man einen Kristall mit Röntgenstrahlen, so treten aufgrund von Spaltbeugungen zwischen den Bausteinen des Kristalls Interferenzerscheinungen auf, die eine Berechnung der Gitterabstände ermöglichen. Die bei der röntgenographischen Untersuchung für Kristalle ermittelten Gitterabstände liegen in der Größenordnung von etwa 0, 1Å, das entspricht 10^{−10}m.

Symmetrien in kristallinen Festkörpern

Die regelmäßige Anordnung der Gitterbausteine bewirkt, daß Kristalle auffallende Symmetrieeigenschaften aufweisen. Sie sind von ebenen Flächen begrenzt, die je nach Art des Kristalles unter ganz bestimmten Winkeln aufeinanderstoßen. Dabei treten Symmetrieachsen und Symmetrieebenen auf. Mathematisch gesehen erfüllen diese Symmetrieelemente für jeden Kristall die Eigenschaften einer algebraischen Gruppe. Es gibt ein kleinstes Arrangement von Bausteinen im Kristall, die sogenannte asymmetrische Einheit oder Basis. Wendet man nacheinander alle Symmetrieoperationen auf die asymmetrische Einheit an, so erhält man eine gedachte Zelle von endlichen Ausmaßen, die Elementarzelle. Reiht man in alle drei Koordinatenrichtungen die Elementarzellen aneinander, so erhält man das Gerüst des makroskopischen Kristalls. Mittels ihrer Elementarzellen lassen sich alle Kristalle beschreiben und in sieben Gittertypen unterteilen:

kubisch, tetragonal, rhomboedrisch, hexagonal, rhombisch, monoklin und triklin.

Die zu diesen Grundtypen gehörigen Symmetrieelemente können in verschiedenen Kombinationen auftreten, wobei sich 32 verschiedene Kristallklassen ergeben. Versuchte man bisher, Kristallstrukturen mathematisch zu erzeugen, so setzte man die Kenntnis der Kristall-Symmetrie voraus. Die einzigen freien Parameter waren die Koordinaten der Atome in der asymmetrischen Einheit. Die Elementarzelle wurde mathematisch dadurch erzeugt, daß Symmetrieoperationen auf die Koordinaten der Atome der asymmetrischen Einheit angewendet werden. Ein Beispiel hierfür findet man im Programm DVLS [18]. Die Vorgabe von Symmetrieelementen bedeutet jedoch, daß sehr detaillierte Informationen über die Kristallstruktur vorliegen müssen. Neue Versuche, Kristallstrukturen zu ”berechnen“, wollen hingegen bei der Vorgabe von Informationen auf die Symmetrien verzichten und dadurch mehr zu einer ganzheitlichen Betrachtung des Kristalls über die Elementarzelle gelangen. Man spricht von einer Kristallstruktur-Vorhersage. Eine umfassende Lösung existiert für dieses Problem bisher nicht, doch gibt es eine Vielzahl von Ansätzen.

Bindungen in kristallinen Festkörpern

Kristalle unterscheiden sich durch die Art ihrer Bausteine und der zwischen den Bausteinen wirkenden Kräfte. Man unterscheidet zwischen Ionenbindungen, Atombindungen, Van–der–Waals–Bindungen und Metallbindungen. Daneben existieren zwischen den genannten Bindungstypen verschiedene bergangsformen. Alle Bindungsarten haben ihre eigenen Charakteristika. Diese sind bei der Vorgabe von Informationen für die Kristallstruktur-Vorhersage von ausschlaggebender Wichtigkeit, denn je nach Art der wirkenden Kräfte können sich die Hüllen unterschiedlicher, benachbarter Atome mehr oder weniger überlappen. Man spricht von der sogenannten Kompressibilität.

2.1.2 Bindungsenergien bei der Kristallstruktur-Vor-hersage

Das Programm ANNEAL verwendet für die Zielfunktion Terme, die sich aus elektrostatischen Überlegungen ableiten lassen. Zur Erläuterung reicht weitestgehend das Ionenmodell. Eine weitere Möglichkeit, Bindungsverhältnisse in Kristallen zu beschreiben, erhält man mit einer empirischen Regel über den Valenzausgleich in Bindungspolyedern. Diese wurde 1929 von Pauling [?] eingeführt.

Die Idee der Paulingschen Valenz-Regel.

Betrachtet man ein Ionengitter, so stellt es sich als Gitter von positiv und negativ geladenen Ionen dar. Der Zusammenhalt der Bausteine kommt im wesentlichen dadurch zustande, daß sich ungleichnamig geladene Ionen elektrostatisch anziehen. Als Standardbeispiel für ein Ionengitter sei das Gitter des Kochsalzes (NaCl) in Abbildung 2.1 schematisch wiedergegeben. Die positiven Natriumionen und negativen Chloridionen sind in den Ecken einer kubischen Elementarzelle so angeordnet, daß jedes Natriumion von sechs Chlorionen umgeben ist. Verbindet man die das Natriumion umgebenden Chlorionen, so ergibt sich ein Polyeder: Man spricht von einem Koordinationspolyeder. Die Anzahl der direkten Nachbarn des Zentralatoms im Koordinationspolyeder wird als Koordinationszahl bezeichnet. Betrachtet man die Abbildung 2.1, so ergibt sich für jedes Natrium- und jedes Chloridion jeweils eine Koordinationszahl von 6.

Abbildung 2.1: Elementarzelle eines Kochsalzkristalls

Die Paulingsche Valenzregel. Pauling beobachtete, daß die Anzahl der benachbarten Atome im Translationsgitter direkt mit der Valenz der beteiligten Atome zusammenhängt (siehe Tabelle 2.1.2).

Valenz 1. Ion Valenz 2. Ion Koordination 1. Ion Koordination 2. Ion
Na^{+}  Cl^{−}  6  6
Si^{4+}  O^{2−}  4  2
Al^{3+}  O^{2−}  6  4

Beispiel: Valenzen und Koordinationszahlen von einigen Kristallen.

Aus dieser Beobachtung entwickelte er 1929^{1} das sogenannte Valenzprinzip als Relation zwischen der formalen Valenz eines Ions und der Koordinationszahl.

Satz 2.1 (Paulingsche Valenzregel) Die Summe der elektrostatischen Valenzen von Bindungen für beliebige Ionen ist numerisch äquivalent zu der formalen Valenz der Ionen.

(2.1) q_{i} = ∑_{j=1}^{CN_{i}} S_{ij},
q_{i} : formale Valenz des Ions i,
CN_{i} = Koordinationszahl des Ions i,
S_{ij} : Bindungsstärke zwischen den Ionen i und j.

S_{ij} steht in einem direkten Zusammenhang zum interatomaren Abstand r_{ij} der Ionen i und j. S_{ij} wird durch einen empirischen Ausdruck der folgenden Form aus r_{ij} hergeleitet (Brown und Altermatt [5]):

1Linus Pauling. Jounal of the American Chemistry Society 51, 1010–1026.

(2.2) S_{ij} = e^{\frac{r{0 \atop ij} −r_{ij}}{B}} ,
r_{ij}^{0} : bindungsabhängige Konstante,
r_{ij} : aktueller Abstand zwischen den Ionen i und j,
B : Bolzmann-Konstante (0.37).

Die Koordinationszahlen CN_{i} ergeben sich aus der Anzahl der nächsten Nachbarn j des Atoms i, die mit dem Atom i eine Bindung eingehen. Sie lassen sich wie die Werte für r_{ij}^{0} aus den Listen von Brown und Altermatt [5] oder bei O’Keeffe und Navrotzky [24] ablesen. Durch Umformung und Summierung über alle Atome der Elementarzelle ergibt sich ein Term der Form

∑_{i∈ EZ} (q_{i} − ∑_{j=1}^{CN_{i}} S_{ij}) . (2.3)

Dieser hat den Wert 0, falls die Paulingsche Regel in allen Bindungspolyedern erfüllt ist. EZ wird im Folgenden als Kürzel für den Begriff ”Elementarzelle“ verwendet.

Elektrostatische Wechselwirkungen.

Elektrostatische Wechselwirkungen sind langreichweitig. Das hat zur Folge, daß der ganze Kristall als großes Molekül aufgefaßt werden muß. Es entsteht ein Gitter. Alle Bausteine werden durch elektrostatische Kräfte, die sogenannten Coulombschen Kräfte, zusammengehalten. Hierbei handelt es sich um eine einfache Erkenntnis der Elektrostatik. Sie besagt, daß sich der Betrag der Kraft zwischen zwei Ladungen im Raum umgekehrt proportional zum Abstand der Ladungen verändert. Da im Translationsgitter eines Kristalls positive Ionen in der Regel von negativen Ionen umgeben sind, und umgekehrt negative Ionen positive Ionen als nächste Nachbarn besitzen, ergibt sich eine unendliche Summe, in der die Anziehungskräfte überwiegen.

Satz 2.2 (Coulombkräfte im Kristall)

(2.4) E_{coul} = \frac{1}{2} · ∑_{i,j} \frac{q_{i} · q_{j}}{r_{ij}} ,
q_{i} : formale Ladung des Ions i,
r_{ij} : aktueller Abstand zwische den Ionen i und j.

Abstoßende Wechselwirkungen.

Coulombkräfte reichen allein nicht aus, um den Aufbau einer Kristallstruktur zu erklären. Weitere Kräfte, die verhindern, daß ungleichnamige Atome denselben Gitterplatz einnehmen, müssen in die berlegungen mit aufgenommen werden. Nähern sich zwei Atome bis auf einen bestimmten Abstand, den sogenannten Ionenradius, so beginnen sich ihre Ladungsverteilungen teilweise zu überlappen. Nimmt der Abstand weiter ab, so wirkt dieser Effekt abstoßend, da zwei Elektronen nicht in allen ihren Quantenzahlen übereinstimmen. Der Bezug zwischen berlappungsenergie und Coulombkraft wird in Abbildung 2.2 wiedergegeben. Der berlappungseffekt erhöht die Gesamtenergie des Systems und ergibt einen abstoßenden Beitrag zur Wechselwirkung. Im Optimum kompensieren sich abstoßende und anziehende Wechselwirkungen. Es gibt

Abbildung 2.2: Abstoßende Wechselwirkung

verschiedene Ansätze, um die abstoßenden Wechselwirkungen zu beschreiben. Einen davon findet man bei Gilbert [10]. Eine Auflistung aller Ansätze findet sich bei Burnham [6].

Satz 2.3 (Gilbertscher Abstoßungsterm)

(2.5) E_{rep} = \frac{1}{2} ∑_{i∈ EZ} ∑_{j}(B_{i} + B_{j}) · e^{\frac{A_{i}+A_{j} −r_{ij}}{B_{i}+B_{j}}} ,
A_{i} : Radiusparameter Ion i,
B_{i} : Kompressibilitätsparameter Ion i.

Das abstoßende Potential ist von geringer Reichweite. Es reicht daher, für die Ionen j in der Summe nur die nächsten Nachbarn des Ions i zu betrachten. Die Gilbertsche Form des Abstoßungsterms wurde schon bei Busing [7] erfolgreich zur Beschreibung von Kristallstrukturen verwendet. Leider sind numerische Werte für die Parameter A_{i}, A_{j}, B_{i} und B_{j} in der Literatur nur schwer zu finden. Sie lassen sich durch röntgenographische Untersuchungen ermitteln oder durch Schätzungen aus anderen Strukturen ableiten. Geschätzte Werte liefern allerdings nicht optimale Ergebnisse, doch reichen sie in der Regel aus, um die Struktur in grober Näherung zu erzeugen. Die Anziehungs– und Abstoßungskräfte und das Paulingsche Valenzprinzip werden für die weiteren Betrachtungen von entscheidender Bedeutung sein.

2.1.3 Mathematisches Modell der Elementarzelle

Um ein mathematisches Modell für die Betrachtung einer Elementarzelle zu erhalten, interpretiert man die Bausteine des Kristalls als Punkte in einem dreidimensionalen, triklinen Quader, welcher die Elementarzelle symbolisiert. Aufgrund der engen Verwandtschaft zwischen dem Modellquader und der Elementarzelle werde ich im Folgenden nur noch die Bezeichnung ”Elementarzelle“ verwenden. Der Modell-Kristall ergibt sich durch Translation der Elementarzelle in alle drei Koordinatenrichtungen. Diese Translation macht sich bei der Berechnung der energetischen Zusammenhänge im Kristall folgendermaßen bemerkbar: Berechnet man für ein Atom der Elementarzelle das zugehörige Koordinationspolyeder, so beinhaltet dieses unter Umständen auch Atome aus benachbarten Zellen (siehe Abbildung 2.3). Man erhält diese Atome durch

Abbildung 2.3: Beispiel für die Lage eines Bindungspolyeders

Translation von Atomen der Elementarzelle. berschreitet ein Atom eine Begrenzungsfläche der Elementarzelle, so bedeutet das, daß durch Translation auf der gegenüberliegenden Seite ein Atom in die Zelle hineingeschoben wird. Auf diese Weise bleibt die Anzahl der Atome in der Elementarzelle konstant und man erhält ein Modell, das den vollständigen Kristall beschreibt (siehe Abbildung 2.4).

Abbildung 2.4: Translationsmodell

2.2 Anpassung des Simulated Annealing - Verfahrens

Aufgrund des Kristallmodells ist die Optimierungsroutine SIMULA ohne weiteres nicht auf das Problem der Kristallstruktur-Vorhersage anwendbar. Folgende Anpassungsroutinen erweisen sich als notwendig.

2.2.1 Zustandsübergänge auf dreidimensionalen Kugelsphären

Die Auswertung der Zielfunktion nimmt im Falle der Kristallstruktur-Vor-hersage viel Zeit in Anspruch. Wie schon in 1.4.1 erwähnt, erwies sich der Algorithmus mit variabler Schrittweite nach Krimphove auf diese Probleme als nicht anwendbar (siehe Tabelle der Laufzeiten, Seite 67). Eine Verwendung des Algorithmus von Radzek ist aufgrund der erheblich geringeren Laufzeit sinnvoller. Leider liefert dieser in seiner ursprünglichen Fassung unzureichende Ergebnisse, da der aktuelle Zustand nur auf einem festen Gitter entlang der Achsen verändert wird. Das Optimum kann aufgrund der Beschaffenheit der Zielfunktion nicht gut genug approximiert werden (siehe Tabelle der beiden Simulated Annealing Implementationen auf Seite 50). Es ist günstiger, die Bewegungen der Atome in der Elementarzelle mehr an die natürlichen Bewegungen von Atomen im Raum anzupassen. Aus dieser berlegung habe ich das Prinzip der dreidimensionalen Kugelsphären entwickelt. Hierbei erfolgen die konstanten Schritte im Gitter nicht mehr nur entlang der festen Achsen im Raum, sondern jedes Atom (also jeweils drei Koordinaten) kann sich mit fester Schrittweite frei in der Elementarzelle bewegen. Die Anzahl der Funktionsauswertungen bleibt in etwa gleich, doch es ergibt sich eine Verbesserung in der Approximation des Minimums. Die Ergebnisse werden von der folgenden Tabelle dokumentiert. Es wurden fünf Durchläufe des Annealing-Algorithmus für NaCL mit der Zielfunktion von Pannetier, λ = 0.95 und der Summation nach der ”Minimum Image Convention for Truncation of the Potential“ (siehe dazu Seite 54) gerechnet.

Funktionswerte bei Schrit- Funktionswerte bei Schrit-
ten auf Kugelsphären  ten entlang der Achsen
1. Versuch  2.925  3.497
2. Versuch  2.885  2.947
3. Versuch  2.911  3.341
4. Versuch  2.901  3.231
5. Versuch  2.931  3.388
Mittelwert  2.911  3.281
Optimalwert  2.870  2.870

Die Tabelle zeigt, daß sich vor allem die Ausgabe des Simulated Annealing-Algorithmus durch die Einbeziehung der Technik der dreidimensionalen Kugelsphären verbessert. Das ist besonders dann wichtig, wenn die Zielfunktion in der Umgebung um das globale Minimum sehr flach ist und lokale Nebenminima aufweist. Diese Topologie entsteht aufgrund der Verwendung der ”Minimum Image Convention“ (siehe dazu Seite 54). Darüber hinaus hat sich das Prinzip der dreidimensionalen Kugelsphären durch eine bessere Approximation des Minimums der Zielfunktion bewährt. Eine abschließende Verfeinerung der Ausgabe ist nicht so aufwendig. Das Prinzip der n-dimensionalen Kugelsphären stellt eine Anwendung des Prinzips der dreidimensionalen Kugelsphären auf allgemeine Optimierungsprobleme dar (siehe Seite 49). Dieses läßt sich als ”Bewegung“ des Gesamtsystems deuten. Für die Vorhersage von Kristallstrukturen ergibt das Prinzip der n-dimensionalen Kugelsphären keine weiteren Verbesserungen in der Approximation des Minimums. Ich habe daher darauf verzichtet. Beim Umgang mit allgemeinen, nichtlinearen Optimierungsproblemen macht hingegen das Prinzip der dreidimensionalen Kugelsphären keinen Sinn, da Gruppen von jeweils drei Variablen für Zustandsänderungen betrachtet werden. Hier macht es Sinn, das Prinzip der n-dimensionalen Kugelsphären anzuwenden. Es hat sich gegenüber den einfachen Zustandsübergängen von Radzek bewährt.

2.2.2 Unbewegliche Teilchen

Bei der Betrachtung der Elementarzelle kann aufgrund der Tatsache, daß sie in alle drei Koordinatenrichtungen unendlich oft translatiert wird, ohne Einschränkung eine Normierung vorgenommen werden. Im Allgemeinen liegt kein Atom genau in einer Ecke der Zelle. Verschiebt man die Zelle (und mit ihr alle translatierten Zellen) so, daß sich ein Atom in der ”linken“ unteren Ecke befindet, so wird die Lage der Atome zueinander davon offensichlich nicht beeinflußt. Der Wert der Zielfunktion ist unverändert. Das Programm ANNEAL betrachtet das erste Atom des Datensatzes ohne Einschränkung als fest in (0, 0, 0). Dadurch ergibt sich eine Reduktion der Problemdimension während des Simulated Annealing. Erst während der Abschlußverfeinerung wird diese Konvention aufgehoben, da es für die direkte Umgebung um ein Minimum sinnvoller ist, wenn sich die Atome in der Elementarzelle frei gegeneinander ausrichten können. Nach der Abschlußverfeinerung wird die Elementarzelle wieder normiert.

2.2.3 Das ”Winkel–Abstands–Programm“

Neben der Ausgabe des Zustands optimaler Kosten besteht der Bedarf, sich ein Bild von der Qualität der berechneten Struktur zu machen. Das geschieht in der Regel graphisch mit Programmen zur graphischen Darstellung von Kristallstrukturen zum Beispiel mit dem Programm ATOMS^{2} dar. Oft aber reichen dreidimensionale Darstellungen allein nicht aus, um beurteilen zu können, wie genau Bindungswinkel und Bindungslängen in der errechneten Struktur eingehalten werden. Zu diesem Zweck habe ich eine Routine mit dem Optimierungsprogramm gekoppelt, die für jedes Atom in der Elementarzelle den vollständigen Bindungspolyeder mit allen Bindungslängen und Bindungswinkeln numerisch generiert. Die Ausgabe orientiert sich an der Ausgabe der Routine ORFFE3 [?] orientiert. Diese wurde schon mit Erfolg im DVLS-Programmpaket [18] verwendet. Durch eine Analyse der erzeugten Werte erhält man einen vollständigen berblick über die numerischen Abweichungen vom Sollwert der Bindungen.

2.2.4 ”Koordinatenkollisionen“

Gegenüber allgemeinen nichtlinearen Optimierungsproblemen gibt es bei der Kristallstruktur-Vorhersage das Problem der ”Koordinatenkollisionen“. Da die elektrostatischen Abstoßungskräfte im Kristall unendlich groß werden, wenn sich zwei Atome auf dem gleichen Gitterplatz im Translationsgitter befinden, entstehen an diesen Stellen Singularitäten, die aus dem Konfigurationenraum ausgeschlossen werden müssen. Dies ist vor allem für diejenigen Prozeduren im Programm wichtig, die Zustandsübergänge realisieren und in denen zufällige Vektoren erzeugt werden. Neu erzeugte Zustände des Systems werden aus diesem Grund, falls sie eine zu große Annäherung zweier Atome zur Folge haben, sofort verworfen. Die Grenze für die Annäherung ist im Programm auf die 0.001-fache Länge des zugehörigen Basisvektors der Elementarzelle festgesetzt. Die Lösung des Optimierungsproblems wird von dieser Einschränkung des Konfigurationenraums nicht beeinflußt.

2ATOMS – A Computer Program for Displaying Atomic Structures by Eric Dowty, 521 Hidden Valley Road, Kingsport TN37663, USA.

2.2.5 Berechnung der Koordinationspolyeder

Die Berechnung der Energieterme kann sich bei der Kristallstruktur-Vor-hersage nicht auf die Elementarzelle beschränken, da Koordinationspolyeder für einzelne Atome der Elementarzelle im Allgemeinen über deren Grenzen hinweg in andere Zellen hineinreichen (siehe Abbildung 2.5). Im Falle des Valenzterms genügt es, die umliegenden 26 Elementarzellen nach benachbarten Atomen, sogenannten Liganden, zu durchsuchen. Im Fall der elektrostatischen Coulomb–Wechselwirkungen gestaltet sich das

Abbildung 2.5: Koordinationspolyeder in verschieden Elementarzellen

Ermitteln der beteiligten Atome aufwendiger. Versuche, auch hier nur eine begrenzte Zahl von Nachbarn zu untersuchen, wie zum Beispiel bei der ”Minimum Image Convention for Truncation of the Potential“ (siehe dazu Seite 54), approximieren den Coulomb-Term nicht hinreichend gut. Eine Verbesserung dieser von Pannetier et al. [25] verwendeten Methode ergibt sich durch die sogenannten Ewald-Summen. Die Vorstellung des durch Bertaud und Williams verbesserten Verfahrens zur Berechnung der Coulomb-Summe wird ein zentraler Bestandteil des nächsten Abschnitts sein.

2.3 Formulierung der verschiedenen Zielfunktionen

Das Problem der Bestimmung einer Zielfunktion für die Kristallstruktur-Vorhersage ist immer noch nicht vollständig gelöst. Erste Ansätze lieferte Khachaturyan [13] schon 1979. Er benutzte als Zielfunktion den sogenannten R–Wert, der sich aus röntgenographischen Untersuchungen am Kristall herleitet. Mittels dieser Zielfunktion gelang es, einfache Kristalle in orthorhombischen Elementarzellen zu berechnen. Schon bald darauf entstanden unterschiedliche Strategien, um sich der Beschreibung von Kristallstrukturen zu nähern. Geometrischen Ansätze, wie zum Beispiel in DVLS [18] oder in WMIN [7], verwendeten Symmetrieoperationen zur Reduktion der Problemdimension und kamen so mit lokalen Optimierungsverfahren aus. Später versuchte man, den Aufbau von Kristallen ohne die Vorgabe von Symmetrieinformationen zu ermitteln. Ein wesentlicher Durchbruch bei der möglichst voraussetzungsarmen Vorhersage von Kristallstrukturen gelang Pannetier et al. 1990 [25]. Mit einer relativ einfachen Zielfunktion gelang es ihnen, eine ganze Anzahl von Kristallstrukturen vorherzusagen. Unter Anwendung eines Simulated Annealing-Algorithmus waren sie in der Lage, Elementarzellen von einfachen Kristallen zu erzeugen.

2.3.1 Die Zielfunktion nach Pannetier et al.

Die von Pannetier et al. [25] verwendete Zielfunktion setzt sich im wesentlichen aus zwei Termen zusammen. Der erste Term modelliert die Bindungen zwischen unterschiedlich geladenen Ionen nach dem Paulingschen Valenzprinzip (siehe 2.1) und wird im Folgenden ”Valenzterm“ (E_{val}) genannt. Der zweite Term modelliert die elektrostatischen Abstoßungskräfte zwischen gleichnamig geladenen Ionen und wird daher ”gleichnamiger Coulomb-Term“ (E_{elec}) genannt. Addiert man die beiden Terme und ergänzt sie durch eine Wichtung λ, so ergibt sich folgende Zielfunktion:

C(x) = E_{val} + λ · E_{elec}, (2.6)

x ∈ IR^{3·n} mit n = Anzahl der Atome in der Elementarzelle.

Die Berechnung des Paulingschen Valenzterms

Der Paulingsche Valenzterm bestimmt die Summe der Abweichungen der Koordinationspolyeder von den optimalen Polyedern und gibt so die gesamte Unausgeglichenheit der Bindungsvalenzen wieder. Die Formulierung der Valenzsumme für die Zielfunktion unterscheidet sich von (??) durch die hinzugefügten Betragstriche:

(2.7) E_{val} = ∑_{i∈EZ} | q_{i} − ∑_{j=1}^{CN_{i}} S_{ij} |,
q_{i} = formale Valenz des Ions i,
CN_{i} = Koordinationszahl des Ions i,
S_{ij} = Bindungstärke zwischen den Ionen i und j (siehe Satz 2.1).

Im Gegensatz zu den meisten zur Strukturbeschreibung verwendeten Energiefunktionen bezieht sich der obige Ausdruck nicht auf ein einzelnes Paar von Wechselwirkungen, sondern beinhaltet die Gesamtheit des Bindungspolyeders. Diese Regel des lokalen Gleichgewichts kontrolliert die Anordnung von Kationen durch Anionen und umgekehrt. Der Pauling’sche Valenzterm kontrolliert nur Kation–Anion Abstände, schränkt aber nicht Kation–Kationoder Anion–Anion-Konstellationen ein.

Die ”Minimum Image Convention“ zur Berechnung des Coulomb-Terms

Um Kation–Kation und Anion–Anion Abstände in der Zielfunktion zu berücksichtigen, schlagen Pannetier et al. [25] eine Methode vor, die eine Wechselwirkung entsprechend der elektrostatischen Coulomb–Wechselwirkung zwischen jeweils zwei gleichen Ladungen berücksichtigt. Das geschieht durch einen Term der folgenden Form:

(2.8) E_{elec} = \frac{1}{2} · ∑_{i∈EZ} ∑_{i′} \frac{q_{i} · q_{i′}}{r_{ii′}} ,
i, i′ = jeweils Kation-Kation- oder Anion-Anion-Paare.

Ich werde ihn im Folgenden als gleichnamigen Coulombterm bezeichnen. Die Berechnung dieser Summe gestaltet sich problematisch. Die äußere Summe durchläuft zwar nur die Ionen der Elementarzelle, die innere jedoch summiert über alle gleich geladenen Ionen des IR^{3}. Die Konvergenz der Summe (2.8) ist gesichert. Einen Beweis erhält man beispielsweise durch die später vorgestellte Methode der Ewald-Summen. Leider konvergiert die Summe nur sehr langsam. Das verdeutlicht das folgende Beispiel. Befinden sich zwei Ionen mit der formalen Ladung 1 in einem Abstand a von einander entfernt, so trägt das den Betrag \frac{1}{a} zur Summe bei. Befinden sich die gleichen zwei Ionen im Abstand 10·a voneinander entfernt, so beträgt deren Wechselwirkung \frac{1}{10} · \frac{1}{a} und beeinflußt damit deutlich die Gesamtsumme. Aus der Voraussetzung, daß sich ein Ion fest im Punkt (0, 0, 0) der Elementarzelle befindet, folgt direkt, daß dasselbe Ion durch Translation in allen Ecken aller Zellen erscheint (siehe Abbildung 2.6). Sollen alle Wechselwirkungen dieser Teilchen bis zum Anteil von 10% des Wechselwirkungsanteils in der Elementarzelle berechnet werden, so ist es notwendig, 20^{3} Elementarzellen abzuarbeiten. Das ist in akzeptabler Zeit nicht durchführbar. Da Kristalle eine regelmäßige Struktur besitzen, schlagen Pannetier et al.

[25] deshalb vor, die Summation nach der ”Minimum Image Convention for the Truncation of the Potential“ durchzuführen. Diese wird von Allem und Tildsley [3] vorgestellt. Ihre berlegung stellt sich etwa wie folgt dar: Der größte Beitrag zur Potentialsumme wird durch die nächsten Nachbarn zu einem Atom geliefert. Für Kräfte kurzer Reichweite verwendet man am besten eine Kugelhülle, in derem Inneren alle betrachteten Atome liegen. Der gestrichelte Kreis in Abbildung 2.7 stellt diese Kugelhülle im IR^{2} dar. Deshalb werden die Ionen 2, 4_{E} und 5_{C} in die berechnete Potentialsumme mit aufgenommen. Benutzt man anstatt der Kugelsphäre die Hülle eines affinen Quaders der Seitenlänge L, so unterscheidet sich die Anzahl der an der Potentialsumme

Abbildung 2.6: Translation im Gitter

beteiligten Ionen um etwa \frac{4πr^{3}}{3L^{3}} von der Anzahl der Ionen in der Kugelsphäre mit Radius r. Für L = r ergibt sich eine Differenz von etwa \frac{4π}{3} ≤ 5 Atomen. Der Vorteil eines affinen Quaders liegt auf der Hand. Verläßt ein Atom die-

Abbildung 2.7: Die ”Minimum Image Convention“ in einem zweidimensionalen System

sen Quader dadurch, daß eine seiner Begrenzungsflächen überschritten wird, so kann man erreichen, daß auf der gegenüberliegenden Seite ein translatiertes Atom vom gleichen Typ in die Elementarzelle eintritt. Diese Vorgehensweise hat den Vorteil, daß der elektrostatische Wechselwirkungsterm keine Sprungstellen hat. Ein ähnliches Vorgehen ist für den Fall der Kugelsphäre nicht möglich. Als Kantenlänge für den triklinen Quader bieten sich für die Kristallstruktur-Vorhersage Kantenlängen von einfachen Vielfachen der Kantenlängen der Elementarzelle an. In meiner Implementation der Pannetier-Zielfunktion habe ich zur Ermittelung des gleichnamigen Coulomb-Terms für jedes Atom einen Quader von der doppelten Größe der Elementarzelle gewählt. Es werden also insgesamt 27 Elementarzellen betrachtet, in deren Mitte sich die Ausgangs–Elementarzelle befindet (siehe Abbildung 2.7).

2.3.2 Die modifizierte Zielfunktion und der erweiterte Ansatz

Es zeigen sich zwei wesentliche Nachteile des Ansatzes von Pannetier et al. [25].

Die Verbesserung des Ansatzes von Pannetier et al. – Die modifizierte Zielfunktion

Das erste Problem besteht in der hohen Empfindlichkeit des gleichnamigen Coulomb-Terms. Der Grund für diese Empfindlichkeit ist in der Verwendung der ”Minimum Image Convention“ zu suchen. Aufgrund der sehr groben Abschätzung ergeben sich nur lokal Auswirkungen auf die Arrangements von Ionen um jedes Zentralatom. Man erhält eine eingeengte Beschreibung der tatsächlichen elektrostatischen Wechselwirkungen in der Elementarzelle. Das hat zur Folge, daß der gleichnamige Coulomb-Term lokal stark schwankt und daß die Wichtungen für den Term nur schwer anzugeben sind. Eine deutliche Verbesserung ergibt sich erst, wenn die Technik der Minimum Image Convention durch eine Methode ersetzt wird, die sich aus einer Berechnung des gesamten Coulomb-Terms ableitet. Eine Methode zur Berechnung des gesamten Coulomb-Terms beruht auf einer Fourierreihenentwicklung und wurde erstmals 1921 von Ewald [?] eingeführt. Die wesentlichen Ergebnisse sollen im Abschnitt 2.3.3 zusammengefaßt werden.

Ein alternativer Ansatz – Die erweiterte Zielfunktion

Der zweite Nachteil des Ansatzes von Pannetier et al. [25] besteht darin, daß die Zielfunktion nur in der Lage ist, Strukturen mit einfachen Konstruktionsprinzipien zu beschreiben: Der Valenzterm besitzt im Gegensatz zum gleichnamigen Coulomb-Term für die optimale Struktur immer ein Minimum mit Funktionswert 0. Der Funktionswert des elektrostatischen Wechselwirkungsterms dagegen ist in der Regel ungleich 0 und besitzt für die optimale Struktur nur dann ein Minimum, wenn die gleich geladenen Ionen in der Elementarzelle der optimalen Struktur soweit voneinander entfernt sind, daß sich ihre Wechselwirkungen gegenseitig die Waage halten. Nur in diesem Fall gelingt mit dem Ansatz von Pannetier et al. die Vorhersage der Struktur. So gelangt man zu Strukturen wie NaCl, CaF_{2}, ZnS, CsCl und mit leichten Anpassungen in der Wichtung sind auch Strukturen wie T iO_{2} und SrT iO_{3} vorhersagbar (siehe Kapitel 3). Wünscht man allerdings, kompliziertere Strukturen vorherzusagen, so reicht eine Beschreibung der Energien, wie sie von Pannetier et al. [25] vorgeschlagen wird, nicht aus. Dies ist z. B. bei Quarz (SiO_{2}) der Fall. Es wird notwendig, nach einer anderen Beschreibung für die Strukturen zu suchen. Man gelangt zu der folgenden Zielfunktion, die im wesentlichen drei Terme enthält:

• Der vollständige Coulomb-Term ergänzt den gleichnamigen Coulomb-Term um die ungleichnamigen Coulomb-Wechselwirkungen. Im Gegensatz zum gleichnamigen Coulomb-Term entsteht eine \underline{anziehende} Wechselwirkung.

• Eine abstoßende Wechselwirkung verhindert, daß Ionen mit gegensätzlichen Ladungen sich überlappen. Physikalisch betrachtet ist diese Kraft von nur sehr geringer Reichweite, da sie die Kräfte beschreibt, die auftreten, wenn sich voll besetzte Orbitale von Elektronen durchdringen. Anschaulich läßt sie sich am besten mit der Kompressibilität eines Festkörpers beschreiben. Ein Ansatz, diese Kraft zu beschreiben ist mit dem Gilbert-Term in 2.3 vorgestellt worden.

• Der Valenzterm kompensiert mögliche Fehler in den Parametern der abstoßenden Wechselwirkung und dient als Indiz für die Fehlerhaftigkeit der Lösung. Dieser ist identisch mit dem in der Zielfunktion von Pannetier et al. [25] verwendeten Term (siehe Formel (2.7)).

Letztendlich gelingt es mit den abstoßenden Wechselwirkungen, den Coulomb-Term so zu kompensieren, daß ein Minimum der Zielfunktion für die optimale Struktur entsteht. Die Parameter für die abstoßenden Wechselwirkungen

Abbildung 2.8: Wirkung von E_{coul} und E_{rep}

sind in der Literatur nur schwer zu finden. Da sie aus röntgenographischen Untersuchungen resultieren, sind sie zudem meist mit Fehlern behaftet. Numerische Werte der Konstanten A und B finden sich bei Busing ([7], Seite 126) (Na, Cl) und bei Burnham ([6], Seite 367) (Mg, Si, O). Sie sind jedoch strukturspezifisch und müssen daher unter Umständen angepaßt werden.

Um ein Mittel zu besitzen, mit welchem die Qualität des Ergebnisses beurteilt werden kann, wurde die Zielfunktion um den Valenzterm erweitert. Im Optimum sollte sich ein Wert nahe 0 für den Valenzterm einstellen.

2.3.3 Die konvergenzbeschleunigende Methode von Ewald, Bertaud und Williams

Beschreibt man einen Kristall als dreidimensional unendliche periodische Anordnung von Atomen und definiert man die Gitterenergie als Arbeit, die zu leisten ist, um die Atome eines Kristalls unendlich weit voneinander zu entfernen, dann läßt sich die elektrostatische Coulomb-Energie pro Zelle folgendermaßen beschreiben:

(2.9) E_{coul} = \frac{1}{2} ∑_{i∈EZ} ∑_{j≠i} \frac{q_{i} · q_{j}}{r_{ij}} ,
q_{i}, q_{j} : Punktladungen im IR^{3},
r_{ij} : Abstand zwischen den Punktladungen i und j.

Die formalen Ladungen q_{i} und q_{j} entsprechen numerisch den formalen Valenzen der Ionen i und j. Man beachte, daß der Index j in der Summe in \underline{allen} Elementarzellen des IR^{3} zu suchen ist. Der Faktor \frac{1}{2} wird deshalb benötigt, weil bei der Summation über das gesamte Gitter jedes Paar von Punktladungen doppelt vorkommt. Es ist bekannt, daß der Ausdruck (2.9) für die Coulomb-Summe nur sehr langsam konvergiert, es sei denn, es werden spezielle Summationstechniken verwandt. Die hier beschriebene Methode wurde erstmals 1921 von Ewald [?] vorgestellt und 1952 von Bertaud [?] weiterentwickelt. Die letzten Modifikationen erfuhr sie 1972 durch Williams [?]. Das Verfahren beruht auf einer Fourierreihen-Entwicklung und einer Thetafunktions–Transformation mit einigen sehr komplexen Umformungen. Detaillierte Hintergrundinformationen sind bei Nijboer und DeWette [23] oder Maradudin [20] zu finden. An dieser Stelle sollen lediglich die Resultate erläutert werden. Für den Umgang mit Punktladungen in Translationsgittern werden die folgenden Hilfsmittel benötigt.

Entfernungs- und Volumenberechnungen in affinen Räumen

Der metrische Tensor. Der metrische Tensor ist eine reellwertige 3 × 3-Matrix, die eine Abstandsberechnung zwischen zwei Punkten des IR^{3} zu einer nicht notwendigerweise orthogonalen Basis, wie zum Beispiel in einer triklinen Elementarzelle, ermöglicht. Gegeben seien die Basisvektoren a, b und c. Ferner sei α der Winkel zwischen den Basisvektoren b und c, β der Winkel zwischen den Basisvektoren a und c und γ der Winkel zwischen den Basisvektoren a und b. Der Abstand zweier beliebiger Punkte x = (x_{1}, x_{2}, x_{3}) und y = (y_{1}, y_{2}, y_{3}) in diesem Koordinatensystem ergibt sich wie folgt:

(2.10) d_{xy} = \sqrt{(x − y)^{T} · C · (x − y)}

mit

C = ( \begin{matrix}a^{2} & a · b · cos(γ) & a · c · cos(β) \\ a · b · cos(γ) & b^{2} & b · c · cos(α) \\ a · c · cos(β) & b · c · cos(α) & c^{2} \end{matrix} ) .

Das Volumen der Elementarzelle. Das Volumen einer Elementarzelle ergibt sich aus den Gitterkonstanten wie folgt:

(2.11) V = a · b · c\sqrt{1 − cos^{2}(α) − cos^{2}(β) cos^{2}(γ) + 2 cos(α) cos(β) cos(γ)}.

Das reziproke Gitter.

Grundlegend für die Methode der beschleunigten Konvergenz nach Ewald, Bertaud und Williams ist der Begriff des reziproken Gitters. Nach Bragg können die Röntgeninterferenzen an Kristallgittern als Reflektionen an verschiedenen Netzebenenscharen dieser Gitter aufgefaßt werden (siehe Borchart-Ott [4]). Um eine bessere bersicht über die Netzebenenscharen zu erhalten, weicht man auf eine Darstellung aus, die auf der Bestimmung der Netzebenennormalen beruht. Im Folgenden werde ich sie zur Verdeutlichung reziprokes Gitter nennen. Koordinaten im reziproken Gitter werden mit h̄ = (h, k, l) bezeichnet und Reflexe genannt. Im Gegensatz dazu spreche ich im Zusammenhang mit den ursprünglichen Koordinaten vom direkten Gitter. Das reziproke Gitter wird in der Mineralogie häufig für röntgenographischen Untersuchungen an Kristallen verwendet, weil es den Fourierraum des direkten Gitters beschreibt. Eine Darstellung von Punktladungen in ”direkten“ Koordinaten oder als ”Reflexe“ ist äquivalent. Die Gitterkonstanten des reziproken Gitters ergeben sich wie folgt aus den Gitterkonstanten des direkten Gitters.

Satz 2.4 (Die Gitterkonstanten des reziproken Gitters) Seien a_{0}, b_{0}, c_{0}, α, β und γ die sechs unabhängigen Gitterkonstanten des direkten Kristallgitters. Die Gitterkonstanten des reziproken Gitters ergeben sich wie folgt:

a_{0}^{∗} = \frac{b_{0} · c_{0} · sin(α)}{V},
b_{0}^{∗} = \frac{a_{0} · c_{0} · sin(β)}{V},
c_{0}^{∗} = \frac{a_{0} · b_{0} · sin(γ)}{V},

V ist das durch Ausdruck (2.11) berechnete Volumen der Elementarzelle.

Ebenso wie im direkten Kristallgitter existiert auch für das reziproke Gitter ein metrischer Tensor. Er läßt sich analog zu (2.10) aus den reziproken Gitterkonstanten erzeugen.

Formulierung der Ewald–Bertaud–Williams–Methode.

Nach diesen Vorbereitungen läßt sich nun die Ewald–Bertaud–Williams– Methode formulieren. Ich werde mich darauf beschränken, die Formel von Williams vorzustellen. Die Elemente der Formel werden im Anschluß daran erläutert. Eine umfassende Behandlung des Problems der unendlichen Gittersummen findet sich beispielsweise bei Maradudin [20].

Satz 2.5 (Williams) Williams zeigte, daß sich die Summe (2.9) folgendermaßen schreiben läßt:

(2.12) E_{coul} = \frac{1}{2} ∑_{i∈EZ} ∑_{j≠i} \frac{q_{i} · q_{j}}{r_{ij}} Erfc(a)

+ \frac{1}{2πV} ∑_{h̄≠0}^{reziprokes Gitter} (\frac{|F_{c}(h̄)|}{Q(h̄)})^{2} · e^{−b^{2}}
− K ∑_{i∈EZ} q_{i}^{2}.

Erfc ist hierbei das Komplement der sogenannten ”Gaußschen Fehlerfunktion“

Erfc(x) = 1 − \frac{2}{π} ∫_{0}^{x} e^{−t^{2}}dt.

Dieser Ausdruck läßt sich numerisch durch Integration über die Potenzreihenentwicklung mit folgendem Ausdruck approximieren:

Erfc(x) ≈ T (a_{1} + T (a_{2} + T (a_{3} + T (a4 + T · a_{5})))) · e^{−a^{2}},
T = \frac{1}{1 + 0.327511 · x},
a_{1} = 0.254829592,
a_{2} = −0.284496736,
a_{3} = 1.421413741,
a_{4} = −1.453152027,
a_{5} = 1.061405429.

K heißt ”Ewald–Konstante“. Sie wird problemabhängig gewählt und ist der einzige freie Parameter der Methode. ber die Ewald–Konstante wird die Gewichtung der Summen über das reziproke und das direkte Gitter vorgenommen. Die Wahl von K beeinflußt die Konvergenzgeschwindigkeit des Verfahrens. Gesicherte Erkenntnisse über eine geeignete Wahl von K liegen nicht vor und sind daher heuristisch zu ermitteln. Die Konvergenz der gesamten Summe bleibt jedoch von der Wahl von K unangetastet. Für K = 0 erhält man die Coulomb-Summe in ihrer ursprünglichen Form. Auf die Formelteile, die mit reziproken Gittervektoren arbeiten, werde ich erst später eingehen.

Legt man diese zugrunde, ergeben sich die übrigen Konstanten wie folgt:

a = \sqrt{πK^{2}r_{ij}^{2}},
b = \sqrt{\frac{πQ(h̄)^{2}}{K^{2}}} ,
V = Volumen der Elementarzelle.

Die erste Summe in (2.12) werde ich im Folgenden als direkte Gittersumme, die zweite als reziproke Gittersumme bezeichnen. Die Konvergenz des Verfahrens wird dadurch erreicht, daß jede der Summen einen Modifikationsfaktor besitzt. Dieser Modifikationsfaktor ist im Fall der direkten Gittersumme das Komplement der Gaußschen Fehlerfunktion. Im Fall der reziproken Gittersumme wird beschleunigte Konvergenz durch den Faktor e^{−b^{2}} erreicht. Beide Faktoren nähern sich für wachsendes r_{ij} schnell dem Wert 0. Der Betrag, der durch den Modifikationsfaktor in der direkten Gittersumme verloren geht, wird exakt durch die reziproke Gittersumme kompensiert. Die Summation über das reziproke Gitter erfolgt über die Vektoren h̄ = {h, k, l}. Im direkten Raum entspricht h̄ den Achsenabschnitten einer Netzebene. Da Röntgenstrahlen an den Netzebenen eines Kristalls ”reflektiert“ werden, spricht man auch von ”Röntgenreflexen“ (h, k, l). Sie ergeben sich aus den ”echten“ Koordinaten des Ions j nach Satz 2.4). Man erhält für F_{c}(h̄) den Ausdruck:

(2.13) F_{c}(h̄) = ∑_{i∈EZ} q_{i} · e^{2πih̄·x_{i}}

mit

x_{i} = Koordinatenvektor des Atoms i,
h̄ = Millersche Indizes (h,k,l) des aktuellen Reflexes.

Formt man diese Darstellung noch weiter um, so ergibt sich die im Programm verwendete Formel:

(2.14) |F_{c}(h̄)|^{2} = A_{c}^{2} + B_{c}^{2}

mit

A_{c} = ∑_{i∈EZ} q_{i} · cos(2 · π · h̄ · x_{i}),
B_{c} = ∑_{i∈EZ} q_{i} · sin(2 · π · h̄ · x_{i}).

Die Größe Q(h̄) = |h̄| stellt die Länge des aktuellen reziproken Gittervektors dar.

Vereinfachung der reziproken Gittersumme. Die reziproke Gittersumme muß nicht über das gesamte Gitter berechnet werden. Nach dem Friedelschen Gesetz gilt |F_{c}(h̄)| = |F_{c}(−h̄)|. Ebenso ist trivialerweise Q(h̄) = Q(−h̄). Man summiert daher nicht in alle Richtungen, sondern nur über einen Halbraum des reziproken Gitters. Dazu wird über zwei der drei reziproken Gitterrichtungen die vollständige Gittersumme errechnet. Die Summierung in die dritte Gitterrichtung dagegen erfolgt lediglich über eine Halbachse.

Allgemeine Bemerkungen. Die Implemetierung der Methode von Ewald, Bertaud und Williams beruht auf einer Beschreibung von Busing [7]. Busing geht von einem Sphärenmodell aus. Die Summation erfolgt über jeweils zwei Sphären des direkten und des reziproken Gitters. Unterscheiden sich die Ergebnisse der beiden Summen um eine festgesetzte kleine Zahl δ, so werden die Radien der Sphären vergrößert. Damit ist eine genügend gute Konvergenz der Coulomb-Summe gewährleistet. Im Programm ANNEAL wurden die Werte der Radien und die Ewald-Konstante konstant gewählt. Anfängliche Befürchtungen, die Werte könnten zu klein gewählt sein und müßten daher während des Programmablaufs geändert werden, haben sich für die berechneten Kristalle nicht bewahrheitet. Aus diesem Grund habe ich anfänglich programmierte Vergleichs- und Anpassungs-Routinen aus dem Programm genommen, um so Rechenzeit zu sparen. Sollten bei anderen Kristallstrukturen Unregelmäßigkeiten auftreten, sind die Werte der Sphärenradien anzupassen. Sie werden in der Unterroutine INIT belegt. Insgesamt hat sich das Verfahren von Ewald, Bertaud und Williams durchweg als wesentlich stabiler als die Minimum Image Convention erwiesen.

2.4 Zusammenfassung der Zielfunktionen

Für die modifizierte Zielfunktion, welche die Methode von Ewald, Bertaud und Williams verwendet, ergibt sich die folgende Darstellung:

(2.15) E_{ges} = w_{val} · E_{val} + w_{elec} · E_{elec} + w_{rep} · E_{rep}

Für die erweiterte Zielfunktion, welche ebenfalls die Methode von Ewald, Bertaud und Williams verwendet, ergibt sich die folgende Darstellung:

(2.16) E_{ges} = w_{val} · E_{val} + w_{coul} · E_{coul} + w_{rep} · E_{rep}

mit

w_{val} = Wichtung für den Valenzterm,
w_{elec} = Wichtung für den gleichnamigen Coulomb-Term,
w_{coul} = Wichtung für den Coulomb-Term,
w_{rep} = Wichtung für den Abstoßungsterm.
E_{val} = ∑_{i∈EZ} | q_{i} − ∑_{j=1}^{CN_{i}} S_{ij} |,
q_{i} =  formale Valenz des Ions i,
CN_{i} =  Koordinationszahl des Ions i,

S_{ij} = e^{\frac{r_{0}−rij}{B}} ,
B = 0.37.
E_{elec} = \frac{1}{2} ∑_{i∈EZ} ∑_{i≠i′} \frac{q_{i} · q_{i′}}{r_{ii′}} Erfc(a),

+ \frac{1}{2πV} ∑_{h≠0}^{reziprokes Gitter} (\frac{|F_{c}(h̄)|}{Q(h̄)})^{2} · e^{−b^{2}}
− K ∑_{i∈EZ} q_{i}^{2}, i,i’ gleichnamig geladen,
E_{coul} = \frac{1}{2} ∑_{i∈EZ} ∑_{j≠i} \frac{q_{i} · q_{j}}{r_{ij}} Erfc(a),

+ \frac{1}{2πV} ∑_{h≠0}^{reziprokes Gitter} (\frac{|F_{c}(h̄)|}{Q(h̄)})^{2} · e^{−b^{2}}
− K ∑_{i∈EZ} q_{i}^{2},
Erfc(x) ≈ T (a_{1} + T (a_{2} + T (a_{3} + T (a4 + T · a_{5})))) · e^{−a^{2}},
T = \frac{1}{1 + 0.327511 · x},
a_{1} = 0.254829592,
a_{2} = −0.284496736,
a_{3} = 1.421413741,
a_{4} = −1.453152027,
a_{5} = 1.061405429,
a = \sqrt{πK^{2}r_{ij}^{2}},
b = \sqrt{\frac{πQ(h̄)^{2}}{K^{2}}} ,
V = Volumen der Elementarzelle,
|F_{c}(h̄)|^{2} = A_{c}^{2} + B_{c}^{2},
A_{c} = ∑_{i∈EZ} q_{i} · cos(2 · π · h · x_{i}),
B_{c} = ∑_{i∈EZ} q_{i} · sin(2 · π · h · x_{i}),
x_{i} = Koordinatenvektor des Atoms i,
h̄ = Vektor (h,k,l) im reziproken Raum,
Q(h̄) = |h̄|,
K = 0.31.
E_{rep} = \frac{1}{2} ∑_{i∈ EZ} ∑_{j}(B_{i} + B_{j}) · e^{\frac{A_{i}+A_{j} −r_{ij}}{B_{i}+B_{j}}} ,

A_{i} =  Radiusparameter des Ions i,
B_{i} =  Kompressibilitätsparameter des Ions i.

Berechnete Kristalle

Die verschiedenen Zielfunktionen, die im letzten Abschnitt vorgestellt wurden, liefern die Grundlage für das Programm ANNEAL. Für die verwendeten Zielfunktionen ergibt sich folgende Tabelle:

• Der Ansatz von Pannetier besteht aus der von Pannetier et al. [25] vorgestellten Formulierung der Zielfunktion aus gleichnamigem Coulomb-Term und Valenzterm. Dieser Ansatz läßt sich noch weiter untergliedern:

1. Die Zielfunktion von Pannetier verwendet zur Berechnung der gleichnamigen Coulomb-Summe im Kristall die Minimum Image Convention.

2. Die Berechnung der gleichnamigen Coulomb-Summe beruht bei der modifizierte Zielfunktion auf der Methode von Ewald, Bartaud und Williams.

• Der erweiterte Ansatz besteht aus einer Formulierung der Zielfunktion mit drei Termen:

– die vollständige Coulomb-Summe umfaßt \underline{alle} elektrostatischen Coulomb-Wechselwirkungen im Kristall.

– die abstoßende Wechselwirkung geringer Reichweite beschreibt die berlappungsenergie von Elektronenschalen.

– der Valenzterm kompensiert Fehler in den Parametern der abstoßenden Wechselwirkung und dient als Indiz für die Fehlerhaftigkeit der Lösung.

ANNEAL stützt sich auf die Optimierungsroutine SIMULA. Die Verwendung beider Programme sowie die Auswahl der verschiedenen Zielfunktionen wird im Anhang beschrieben. Verschiedene Durchläufe von ANNEAL für Parametersätze von unterschiedlichen Kristallen lassen sich hinsichtlich der Anwendbarkeit der verschiedenen Simulated Annealing-Strategien und der Qualität der erzeugten Ergebnisse auswerten. Ziel der Diskussion ist es, die Leistungen und Grenzen der mathematischen Verfahren für konkrete Anwendungen in der Mineralogie herauszustellen. Des weiteren sollen die Ergebnisse der verschiedenen Zielfunktionsansätze gegenübergestellt werden. Der Präsentation liegen Anwendungen für die Kristallstrukturen von CsCl, NaCl, ZnS, CaF_{2}, T iO_{2}, SrT iO_{3} und SiO_{2} zugrunde. Die Graphiken wurden mit dem Programm ATOMS von E. Dowty^{1} erstellt. Eine Bewertung der Ergebnisse folgt im Anschluß. Numerische Werte für die Zellparameter findet man bei Megaw [22]. Im Programm ANNEAL wird, wie bei den meisten kristallographischen Anwendungen, mit normierten Koordinaten gerechnet. Das bedeutet, daß die Koordinaten des optimalen Vektors zwischen 0 und 1 angegeben werden. Die Berechnung der echten Abstände für die Zielfunktion erfolgt mit dem metrischen Tensor. Auf diese Art und Weise ist es leichter möglich, sich ein Bild über die Lage der Atome in der Elementarzelle zu machen und diese mit vorliegenden Tabellenwerken zu vergleichen. Mineralogisch interessant sind sogenannte metastabile Phasen. Sie geben Hinweis auf energetische Verwandschaften zwischen verschiedenen Strukturen. Mathematisch bedeutet dies, daß ein Nebenminimum der Zielfunktion existiert. Ich werde bei einigen Strukturen darauf eingehen.

3.1 Cäsiumchlorid

Cäsiumchlorid (CsCl) ist ein einfaches Beispiel für die Kristallstruktur- Vorhersage. Die Struktur läßt sich schon alleine mit dem Valenzterm sehr gut vorhersagen und steht als Beispiel für die Klasse von Kristallen, die sich aufgrund der geringen Problemdimension mit allen Simulated Annealing-Strategien bearbeiten läßt.

Beschreibung der Struktur

Das Gitter von CsCl ist kubisch innenzentriert. Die Kantenlänge der Elementarzelle beträgt 4,11 Å. Die Elementarzelle enthält eine Formeleinheit. Cs ist einfach positiv geladen und befindet sich im Punkt (0, 0, 0) der Elementarzelle, Cl ist einfach negativ geladen und befindet sich im Punkt (\frac{1}{2}, \frac{1}{2}, \frac{1}{2}). Jedes Chloridion ist von 8 Cs-Ionen umgeben, die sich in den Ecken der Elementarzelle befinden. Umgekehrt bilden jeweils 8 Chloridionen das Koordinationspolyeder für ein Cs-Ion. Der Abstand zwischen den unterschiedlichen Ionen beträgt \frac{\sqrt{3}}{2} · a ≈ 3, 56 Å. Andere Verbindungen desselben Strukturtyps lassen sich mit derselben Zielfunktion vorhersagen.

1Atoms – A Computer Program for Displaying Atomic Structures, Version 2.2. by Eric Dowty, 521 Hidden Valley Road, Kingsport TN 37663, USA

Bemerkungen zur Strukturvorhersage

Aufgrund der Tatsache, daß die Elementarzelle von nur zwei Atomen gebildet wird, ist es möglich, diese mit relativ geringen Mitteln vorherzusagen. Bezogen auf die benachbarten Elementarzellen ändert sich der Abstand der gleichnamigen Teilchen zueinander nicht. Die Zielfunktion ergibt sich somit aus einem Term, der die elektrostatischen Wechselwirkungen zwischen unterschiedlich geladenen Teilchen beschreibt. Hier reicht schon die Verwendung des Valenzterms. Aufgrund der geringen Problemdimension ist dies auch das einzige Beispiel der Kristallstruktur-Vorhersage, in welchem der Algorithmus von Krimphove zur Anwendung kam. In der folgenden Tabelle erkennt man, daß Kristallstruktur-Vorhersagen höherer Dimension nicht mehr in akzeptabler Zeit durchführbar sind.

konstant 3-dimension. Sphäre variabel^{2} variabel^{3}
Anz. Funktions- 5401 6201 10262 52942

auswertungen

Die nächste Tabelle gibt die Genauigkeit der berechneten Struktur wieder. Im Folgenden wurde für alle Tabellen ein Algorithmus mit Zustandswechseln auf 3-dimensionalen Kugelsphären konstanter Radien verwendet. Die Abweichung von echter und berechneter Struktur ist minimal, Winkel und Abstände werden fast vollständig eingehalten.

Algorithmus echte Struktur Abweichung
ø Abstand Cs − Cl  3,56Å  3,56Å  0Å
min. Abstand Cs − Cl  3,56Å  3,56Å  0Å
max Abstand Cs − Cl  3,56Å  3,56Å  0Å
min. Winkel Cl − Cs − Cl  71/109^{∘}  71/109^{∘}  < 1^{∘}
max. Winkel Cl − Cs − Cl  71/109^{∘}  71/109^{∘}  < 1^{∘}
min. Winkel Cs − Cl − Cs  71/109^{∘}  71/109^{∘}  < 1^{∘}
max. Winkel Cs − Cl − Cs 71/109^{∘}  71/109^{∘}  < 1^{∘}

25·5 Schrittweitenanpassungen

310·10 Schrittweitenanpassungen

Abbildung 3.1: Darstellung der berechneten Elementarzelle von Cäsiumchlorid

Abbildung 3.2: Darstellung der berechneten Elementarzelle von Cäsiumchlorid im Verbund mit umliegenden Zellen

3.2 Halit, NaCl

Die Struktur von Halit (NaCl) läßt sich mit allen Zielfunktionen vorhersagen und macht damit einen Vergleich der verschiedenen Ansätze möglich. Ferner hat sich NaCl als Versuchsstruktur für Verbesserungen am Algorithmus bewährt.

Beschreibung der Struktur

Das Gitter von Halit ist kubisch flächenzentriert. Die Kantenlänge der Elementarzelle beträgt 5,63 Å. Die Elementarzelle enthält vier Formeleinheiten. Jedes Na-Atom wird von jeweils 6 Cl-Atomen umgeben und umgekehrt bilden 6 Cl-Atome das Koordinationspolyeder für jeweils ein Na-Atom. Der Abstand zwischen Zentralatom und Ligand in den Koordinationspolyedern ist 2, 81 Å. Die Struktur von Halit ist sehr weit verbreitet und taucht bei nahezu allen Alkali-Halogeniden (mit Ausnahme von CsCl, CsBr und CsI) auf. Ferner findet man die NaCl-Struktur in allen Oxiden, Sulfiden, Seleniden und Telluriden von Mg, Ba, Ca und Sr. Es sollte keine Probleme bereiten, diese mit derselben Zielfunktion vorherzusagen.

Bemerkungen zur Strukturvorhersage

Halit läßt sich mit allen Arten von Zielfunktionen vorhersagen. Schon der Ansatz von Pannetier unter Verwendung der Minimum Image Convention für das elektrostatische Abstoßungspotential E_{elec} liefert für λ = 0, 95 gute Ergebnisse. Problematisch ist lediglich die Bestimmung von λ (siehe folgende Tabelle).

λ = 0.95  Algorithmus echte Struktur Abweichung
ø Abstand Na − Cl 2,82Å 2,81Å 0,01Å
λ = 0.5  Algorithmus echte Struktur Abweichung
ø Abstand Na − Cl  2,69Å  2,81Å  0,12Å

Pannetier et al. geben zwar eine Möglichkeit an, die Wichtung zu schätzen, doch führt dies vielfach nicht zu geeigneten Ergebnissen. Verwendet man hingegen die Methode der Ewaldsummen, um das elektrostatische Abstoßungspotential zu berechnen, so kommt man für Strukturen wie Halit generell mit einer Wichtung von λ = 1 zu denselben Ergebnissen. Erstmals wurde hier auch der neu vorgestellte erweiterte Ansatz, bestehend aus gesamter Coulombsumme und Gilbertschem Abstoßungspotential gewählt, um die Struktur zu erzeugen. Die folgende Tabelle zeigt die Qualität des erzielten Ergebnisses.

Algorithmus echte Struktur Abweichung
ø Abstand Na − Cl  2,81Å  2,81Å  0Å
min. Abstand Na − Cl  2,809Å  2,81Å  < 0.01Å
max Abstand Na − Cl  2,81Å  2,81Å  0Å
min. Winkel Na − Cl − Na  90^{∘}  90^{∘}  < 1^{∘}
max. Winkel Na − Cl − Na  90^{∘}  90^{∘}  < 1^{∘}

Metastabile Phasen

Bei der Berechnung der Halit-Struktur offenbart sich ein interessanter Bezug zu der Struktur von Zinkblende (ZnS). Berechnet man Halit mit der modifizierten Zielfunktion unter Verwendung der Terme von Pannetier et al. [25], so bildet sich im Falle eines Nebenminimums die Struktur von Zinkblende aus. Der Grund für dieses Verhalten ist in der Tatsache zu suchen, daß der gleichnamige Coulomb-Term für die Translationsgitter von Halit und Zinkblende gleich groß ist. Die folgende Tabelle zeigt, daß lediglich die Paulingsche Valenzsumme von 0 abweicht.

NaCl-Struktur ZnS-Struktur
Coulombsumme  29,74  29,75
Valenzsumme  0,02  3,12

Verschiedene Strukturen mit der modifizierten Zielfunktion.

Da sich jedoch auch Zinkblende mit diesem Ansatz beschreiben läßt (siehe Seite 3.3), entsteht ein Nebenminimum zur thermodynamisch stabilen Phase, eine ”metastabile Phase“.

Abbildung 3.3: Darstellung der berechneten Elementarzelle von Natriumchlorid

Abbildung 3.4: Darstellung der berechneten Elementarzelle von Natruimchlorid im Verbund mit umliegenden Zellen

3.3 Zinkblende, ZnS

Da die Zinkblendestruktur als metastabile Phase von Kochsalz auftaucht, bereitet der Versuch, die Struktur mit der modifizierten Zielfunktion und dem Ansatz von Pannetier et al. [25] vorherzusagen, keine Probleme.

Beschreibung der Struktur

Das Gitter von Zinkblende ist wie das Halit-Gitter kubisch flächenzentriert. Die Kantenlänge der Elementarzelle beträgt 5,41Å. Die Elementarzelle wird von vier Formeleinheiten ZnS gebildet. Jedes Zn-Atom wird im Gegensatz zum Kochsalz von nur 4 S-Atomen umgeben und umgekehrt. Das Koordinationspolyeder ist ein reguläres Tetraeder. Der Abstand zwischen Zentralatom und Ligand ist 2, 342 Å. Die Zinkblende-Struktur ist wie die Halit-Struktur weit verbreitet. Beispiele sind BeS, CuF, CdS, GaS und InSb.

Bemerkungen zur Strukturvorhersage

Wie Halit läßt sich auch Zinkblende mit der modifizierten Zielfunktion und dem Ansatz von Pannetier et al. [25] vorhersagen. Die Terme der Zielfunktion werden dabei gleich stark gewichtet. Der erweiterte Ansatz, bestehend aus gesamter Coulomb-Summe und Gilbertschem Abstoßungspotential, wurde aufgrund von fehlenden Parametern für das Gilbert-Potential nicht durchgerechnet. Daß ein solcher Ansatz zur Vorhersage der ZnS-Struktur nicht nötig ist, zeigt folgende Tabelle. Die Bindungswinkel und die Abstände in den Koordinations-Polyedern werden auch mit der modifizierten Zielfunktion und dem Ansatz von Pannetier et al. [25] genau genug beschrieben.

Algorithmus echte Struktur Abweichung
ø Abstand Zn − S  2,34Å  2,34Å  < 0, 01Å
min. Abstand Zn − S  2,33Å  2,34Å  0,01Å
max Abstand Zn − S  2,37Å  2,34Å  0,03Å
min. Winkel Zn − S − Zn  109^{∘}  109^{∘}  < 1^{∘}
max. Winkel Zn − S − Zn  110^{∘}  109^{∘}  1^{∘}

Metastabile Phasen

Aufgrund der engen Verwandtschaft zwischen der Zinkblende-Struktur und der Halit-Struktur liegen metastabile Phasen vom Kochsalztyp nahe. Diese treten jedoch erheblich seltener auf als umgekehrt der Typ der Zinkblende-Struktur bei der Vorhersage der Halit-Struktur.

Abbildung 3.5: Darstellung der berechneten Elementarzelle von Zinkblende, ZnS

Abbildung 3.6: Darstellung der berechneten Elementarzelle von Zinkblende, ZnS im Verbund mit umliegenden Zellen

3.4 Fluorit, CaF_{2}

Fluorit steht als erstes Beispiel für die Klasse von Kristallen, bei denen das Verhältnis von Kationen zu Anionen ungleich Eins ist. Die Struktur zeichnet sich nicht durch komplizierten Aufbau aus, ist jedoch aufgrund der höheren Anzahl von Variablen mathematisch von Interesse. Mit zwölf Atomen in der Elementarzelle erhält man einen Zustandsvektor mit 36 Variablen und kann sich einen Ausblick darüber verschaffen, wie sich die Laufzeit für größere Variablenzahlen bei der Kristallstruktur-Vorhersage entwickelt. Die nächste Tabelle gibt Laufzeiten für einige Strukturen wieder, die bei Rechnungen auf einer IBM RS-6000 erzielt wurden.

CsCl NaCl CaF_{2} Al_{2}O_{3}
Variablenzahl  6  24  36  90
Laufzeit (Min.) 2  53  201 3138

Sind Strukturen wie die von Halit und Zinkblende noch in akzeptabler Zeit auf einem Personalcomputer zu berechnen (486 DX-50: ungefähr 100 Minuten), so stößt man für mehr Variablen schnell an Grenzen. Die größte bisher mit dem Programm ANNEAL bearbeitete Struktur ist die Korund-Struktur (Al_{2}O_{3}) mit 90 Variablen. Ein Testlauf des Simulated Annealing mit der modifizierten Zielfunktion kostet für diese Struktur schon etwa 50 Stunden echte CPU-Zeit auf einer IBM RS-6000. Leider ist die modifizierte Zielfunktion nicht in der Lage, die Al_{2}O_{3}-Struktur vollständig zu beschreiben und Gilbert– Parameter für die erweiterte Zielfunktion liegen bis jetzt nicht vor. Trotzdem zeigt sich mit der Rechendauer des Testlaufs, wo die Grenzen der Kristallstrukturvorhersage mit dem Simulated Annealing-Algorithmus liegen.

Beschreibung der Calciumfluorit-Struktur

Wie die vorherigen beiden Gitter ist das Gitter des CaF_{2} kubisch flächenzentriert. Die Kantenlänge der Elementarzelle beträgt 5,45 Å. Die Elementarzelle enthält vier Formeleinheiten CaF_{2}, das heißt 12 Atome. Jedes Ca-Atom wird von 8 F -Atomen umgeben, jedes F -Atom verfügt jedoch nur über ein Koordinationspolyeder bestehend aus 4 Ca-Atomen. Dieses hat die Gestalt eines regulären Tetraeders. Der Abstand zwischen Zentralatom und Ligand in den Koordinationspolyedern ist konstant 2, 36 Å. Im gleichen Strukturtyp kristallisieren SrF_{2}, BaF_{2} und T hO_{2}.

Bemerkungen zur Strukturvorhersage

Wie schon in den vorherigen Beispielen läßt sich auch CaF_{2} mit der modifizierten Zielfunktion und dem Ansatz von Pannetier vorhersagen. Die Terme der Zielfunktion werden dabei immer noch gleich stark gewichtet. Das Ergebnis ist in der folgenden Tabelle dokumentiert.

Algorithmus echte Struktur Abweichung
ø Abstand Ca − F  2,35Å  2,36Å  0,01Å
min. Abstand Ca − F  2,39Å  2,36Å  0,03Å
max Abstand Ca − F  2,32Å  2,36Å  0,04Å
min. Winkel Ca − F − Ca  71^{∘}  71^{∘}  < 1^{∘}
max. Winkel Ca − F − Ca  109^{∘}  109^{∘}  < 1^{∘}
min. Winkel F − Ca − F  71^{∘}  71^{∘}  < 1^{∘}
max. Winkel F − Ca − F  109^{∘}  109^{∘}  < 1^{∘}

KAPITEL 3. BERECHNETE KRISTALLE 76 Abbildung 3.7: Darstellung der berechneten Elementarzelle von Fluorit, CaF_{2} Abbildung 3.8: Darstellung der berechneten Elementarzelle von Fluorit, CaF_{2} im Verbund mit umliegenden Zellen

3.5 Rutil, T iO_{2}

Ein Beispiel für eine Klasse von Strukturen, die nicht ohne nderung in der Wichtung mit dem Ansatz von Pannetier et al. [25] beschrieben werden kann, ist Rutil (T iO_{2}). Der Grund hierfür liegt in der Anordnung der Anionen. Die Rutil-Struktur ist ferner interessant, weil mit T iO_{2} erstmals ein Kristall vorhergesagt wird, dessen Elementarzelle nicht kubisch ist.

Beschreibung der Struktur

Das Gitter von Rutil (T iO_{2}) ist tetragonal primitiv. Die Kantenlängen der Elementarzelle sind a = 4, 594Å, b = 4, 594Å und c = 2, 958 Å. Die Elementarzelle enthält zwei Formeleinheiten T iO_{2}, das heißt 6 Atome. Jedes T i-Atom wird von 6 O-Atomen umgeben, jedes O-Atom verfügt jedoch nur über ein Koordinationspolyeder bestehend aus 3 T i-Atomen. Die Oktaeder um die T i-Atome sind nicht regulär. Es gibt zwei verschiedene T i-O-Bindungslängen von 1, 95 Å und 1, 98 Å. Unterschiedlicher sind jedoch die O-O-Abstände. Hiervon gibt es drei verschiedene Typen der Längen 2, 55, 2, 78 und 2, 98 Å. Im gleichen Strukturtyp kristallisieren Cassiterit (SnO_{2}), MgF_{2} und ZnF_{2}.

Bemerkungen zur Strukturvorhersage

Die Rutilstruktur läßt sich, trotz leichter Schwierigkeiten, immer noch mit der modifizierten Zielfunktion und dem Ansatz von Pannetier et al. [25] vorhersagen. Da die Struktur allerdings nicht mehr so ”regelmäßig“ wie die vorhergehenden Strukturen ist, ist eine stärkere Wichtung des Valenzterms nötig, damit sich die Bindungspolyeder um die T i-Atome besser herausbilden. Es zeigt sich, daß man sich an diesem Punkt an der Grenze des mit dem Ansatz von Pannetier Vorhersagbaren befindet. Es treten Kräfte im Kristall auf, die in dieser Form nur noch schwer zu beschreiben sind. Für die Rutilstruktur besitzt der Term der elektrostatischen Abstoßung, wie ihn Pannetier et al. vorschlagen haben, kein Minimum. Mit einer veränderten Wichtung des Valenzterms ist dies jedoch noch zu kompensieren. Gelingt dies bei anderen Strukturen nicht mehr, und dafür wird später noch ein Beispiel gebracht, so werden andere Ansätze notwendig. Die folgende Tabelle stellt die errechneten und die tatsächlichen Bindungslängen gegenüber. Es fällt auf, daß die Abweichungen in diesem Beispiel größer als in den vorangegangenen Strukturen sind.

Algorithmus echte Struktur Abweichung
ø 1.Abstand T i − O  1,81Å  1,95Å  0,14Å
min. 1.Abstand T i − O  1,75Å  1,95Å  0,20Å
max 1.Abstand T i − O  1,87Å  1,95Å  0,08Å
ø 2.Abstand T i − O  2,07Å  1,98Å  0,09Å
min. 2.Abstand T i − O  2,02Å  1,98Å  0,04Å
max 2.Abstand T i − O  2,12Å  1,98Å  0,14Å

KAPITEL 3. BERECHNETE KRISTALLE 78 Abbildung 3.9: Darstellung der berechneten Elementarzelle von Rutil, T iO_{2} Abbildung 3.10: Darstellung der berechneten Elementarzelle von Rutil, T iO_{2} im Verbund mit umliegenden Zellen

3.6 Perovskit, SrT iO_{3}

Perovskit (SrT iO_{3}) ist ein Beispiel für eine Struktur, in der drei verschiedene Sorten von Atomen vorkommen. Hierbei ist interessant, daß es gelingt, auch Strukturen mit unterschiedlichen Kationen und einer Sorte von Anionen richtig vorherzusagen, wenn man die modifizierte Zielfunktion anwendet.

Beschreibung der Struktur

Das Gitter von Perovskit (SrT iO_{3}) ist kubisch primitiv. Die Kantenlänge der Elementarzelle beträgt a = 3, 905 Å. Die Elementarzelle enthält eine Formeleinheit SrT iO_{3}. Jedes T i-Atom wird von 6 O-Atomen umgeben, jedes O-Atom hat genau zwei Nachbarn T i. Das zentrale Sr-Atom wird von 12 Atomen der beiden anderen Typen umgeben. Die Abstände zwischen T i und O betragen 1, 95 Å, die zwischen Sr und O betragen 2, 76 Å. Die Abstände zwischen den gleich geladenen Atomen T i und Sr betragen 3, 38 Å.

Bemerkungen zur Strukturvorhersage

Eine Vorhersage mit dem Ansatz von Pannetier et al. [25] unter Verwendung der Mini- mum-Image-Convention gelingt trotz Variation bei der Wahl der Wichtungen nicht mehr. Erfolgreicher verlief der Versuch, den Ansatz von Pannetier mit der modifizierten Zielfunktion zu berechnen. Auch hier deutet sich an, daß dieser Ansatz für bestimmte Strukturen schnell versagt. Da das Verhältnis von Anionen zu Kationen ungleich 1 ist, kompensieren sich der Valenzterm und der gleichnamige Coulomb-Term nicht vollständig. Darunter leiden vor allem die Bindungen der O-Atome, die nur kleine Koordinationszahlen besitzen. Die Bindungen und Winkel des Kristalls werden ”unscharf“ und der Kristall scheint leicht verzerrt. Eine Vorhersage ist erst möglich, wenn die Wichtung des Valenzterms groß genug gewählt ist. Da keine Parameter für die erweiterte Zielfunktion vorliegen, wurde dieses Beispiel damit nicht bearbeitet.

Algorithmus echte Struktur Abweichung
ø Abstand T i − O  1,93Å  1,95Å  0,02Å
min. Abstand T i − O  1,87Å  1,95Å  0,08Å
max Abstand T i − O  2,03Å  1,95Å  0,08Å
ø Abstand Sr − O  2,69Å  2,76Å  0,07Å
min. Abstand Sr − O  2,52Å  2,76Å  0,24Å
max Abstand Sr − O  2,71Å  2,76Å  0,05Å
ø Abstand Sr − T i  3,36Å  3,38Å  0,02Å

Abbildung 3.11: Darstellung der berechneten Elementarzelle von Perovskit, SrT iO_{3} im Verbund mit umliegenden Zellen

Abbildung 3.12: Darstellung der berechneten Elementarzelle von Perovskit, SrT iO_{3} im Verbund mit umliegenden Zellen

3.7 Quarz, SiO_{2}

Quarz (”Hochquarz“) ist die erste Struktur, die sich dem Ansatz von Pannetier et al. [25], bestehend aus Valenzterm und gleichnamigem Coulomb-Term, widersetzt. Sowohl die ursprüngliche Zielfunktion mit der Minimum Image Convention als auch die modifizierte Zielfunktion lieferten keine Quarzstruktur. Hier zeigt sich eine weitere Klasse von Kristallstrukturen, für die es notwendig ist, alternative Zielfunktionen zu entwickeln. Ein Grund für das Versagen des Ansatzes von Pannetier et al. [25] ist in der komplizierten Struktur der Elementarzelle des Quarz zu suchen. Die Atome der Elementarzelle von SiO_{2} bilden 3-dimensional verknüpfte Tetraeder, die schraubenförmig parallel zur kristallographischen \underline{c}-Achse angeordnet sind. Die Bindungswinkel an den Sauerstoffatomen sind nicht mehr, wie in den vorherigen Strukturen, regelmäßig, sondern schwanken um 140^{o}. Ein sinnvolles Bild der Quarzstruktur ergibt sich erst, wenn man sich einen berblick über mehrere Elementarzellen verschafft. Man erkennt ineinandergreifende Ringe, die abwechselnd aus Si- und O-Atomen bestehen. Interessant ist Quarz vom kristallographischen Standpunkt aus, weil es verschiedene Modifikationen von SiO_{2} gibt. Eng verwandt sind die Modifikationen des ”Hochquarz“ und des ”Tiefquarz“. Sie unterscheiden sich nur unwesentlich in der Ausrichtung der Bindungen. Eine weitere Modifikation ist der Tridymit. Interessanterweise taucht eine Struktur als metastabile Phase von Hochquarz auf, die der Tridymit-Struktur sehr ähnlich ist.

Beschreibung der Struktur

Das Gitter von Quarz (SiO_{2}) ist hexagonal primitiv. Die Kanten a und b der Elementarzelle stehen in einem Winkel von 120^{o} aufeinander. Der Winkel zwischen den Kante b und c sowie a und c is 90^{o}. Die Kantenlänge der Kanten a und b beträgt 5, 01 Å. Die Kante c hingegen ist 5, 47 Å lang. Die Elementarzelle enthält drei Formeleinheiten SiO_{2}, also 9 Atome. Jedes Si-Atom wird von jeweils 4 O-Atomen umgeben, jedes O-Atom hat genau zwei Nachbarn Si.

Bemerkungen zur Strukturvorhersage

Interessant für die Kristallstruktur-Vorhersage ist Quarz auf der einen Seite wegen der vielen möglichen Modifikationen. Auf der anderen Seite stellt die Struktur wegen ihrer Kompliziertheit gehobene Ansprüche an eine Zielfunktion. Der Ansatz nach Pannetier et al. [25] führt immer zu sehr regelmäßigen Elementarzellen. Eine regelmäßige ”Hochquarz“-Elementarzelle verletzt jedoch die Paulingsche Valenzregel erheblich. Es ist also nötig, einen anderen Ansatz zu verwenden, um sich der Quarzstruktur zu nähern. Dabei hat sich der erweiterte Ansatz, bestehend aus Valenzterm, vollständigem Coulomb-Term und Gilbertschem Abstoßungspotential bewährt. Die Parameter für den Gilbert-Term wurden der Veröffentlichung von Burnham [6] entnommen und sind angepaßt, da die angegebenen Parameter für die Struktur von Mg_{2}SiO_{4} ermittelt wurden. Die eigentliche Struktur von Hochquarz wurde schon mit den bei Burnham angegebenen Parametern sichtbar, jedoch waren Bindungslängen und Winkel nicht genau genug eingehalten. Eine abschließende Fitroutine, welche, ausgehend von der bis dahin erzeugten Struktur, den Valenzterm in Abhängigkeit von den Gilbert-Parametern mit einem Verfahren des steilsten Abstiegs optimiert, liefert verbesserte Werte. Die folgende Tabelle zeigt eine bersicht über die beiden erzeugten Strukturen. Für die metastabile Phase fällt erneut auf, daß der Valenzterm der errechneten Struktur schlechter ist, als der Valenzterm der echten Hochquarz-Struktur. besser istnicht 0 wird. Die Paulingsche Valenzregel wurde nicht eingehalten.

berechnete echte metastabile
Struktur Struktur  Phase
ø Abstand Si − O  1,59Å  1,61Å  1,58Å
min. Abstand Si − O  1,48Å  1,61Å  1,39Å
max Abstand Si − O  1,69Å  1,61Å  1,78Å
ø Winkel Si − O − Si  155^{∘}  149^{∘}  38^{∘}
min. Winkel Si − O − Si  149^{∘}  149^{∘}  10^{∘}
max. Winkel Si − O − Si  160^{∘}  149^{∘}  45^{∘}
errechnete Struktur  Valenzterm
Hochquarz  2,673
metastabile Modifikation  3,615

Metastabile Phasen

Tridymit ist eine weitere Modifikation des SiO_{2}. Die Struktur zeichnet sich durch geschlossene Ringe von 6 Si-Tetraedern aus. Jedes Si-Atom besitzt jeweils 4 benachbarte O-Atome. Die Möglichkeit, ähnliche Strukturen wie Tridymit als metastabile Phase von Hochquarz zu erhalten, liegt aufgrund der engen Verwandschaft zwischen Hochquarz und Tridymit, nahe. Daß ein Nebenminimum der Zielfunktion erreicht wurde, läßt sich leicht am Wert des Valenzterms im Optimum ablesen. Ist dieser ungleich 0, so wird ein Nebenminimum gefunden.

KAPITEL 3. BERECHNETE KRISTALLE 83 Abbildung 3.13: Darstellung der berechneten Elementarzelle von Quarz, SiO_{2} Abbildung 3.14: Darstellung der berechneten Elementarzelle von Quarz, SiO_{2} im Verbund mit umliegenden Zellen Abbildung 3.15: Darstellung der berechneten Elementarzelle der metastabilen Phase von Quarz

Abbildung 3.16: Darstellung der berechneten Elementarzelle der metastabilen Phase von Quarz im Verbund mit umliegenden Zellen

3.8 Diskussion

Die Vorhersage von Kristallstrukturen mit Mitteln der Optimierung ist ein relativ neues Gebiet. Lokale Optimierungsstrategien liefern erste Strukturen, die mit Hilfe von Symmetrieoperationen aus der asymmetrischen Einheit erzeugt wurden. Seit der Einführung von neuen, globalen Optimierungsstrategien und wesentlich leistungsfähigeren Rechnergenerationen ist es jedoch möglich, eine Vorhersage von Kristallstrukturen auf der Basis von ganzen Elementarzellen durchzuführen. Die in den vorangehenden Abschnitten vorgestellten Ergebnisse geben einen Einblick in die Möglichkeiten, die das Simulated Annealing bietet. Nach leichten Modifikationen hat sich das Verfahren durchweg als erfolgreich erwiesen. Die Tatsache, daß öfter auch lokale Minima gefunden werden, stört dabei den Gesamteindruck nicht, sondern paßt in das physikalische Modell, aus welchem das Simulated Annealing-Verfahren abgeleitet ist. Wird ein Kristall aus der Schmelze zu schnell abgekühlt, ordnen sich die Bausteine nicht im erwarteten Gitter an. Dieser Effekt ist mineralogisch nicht bedeutungslos, denn die gefundenen lokalen Minima signalisieren mögliche energetische Verwandschaften zwischen verschiedenen Kristallmodifikationen und bieten Ansatzpunkte für neue Versuche, Kristallstrukturen vorherzusagen. Leider, so zeigt sich im letzten Kapitel, ist eine Kristallstruktur-Vorhersage mit Mitteln der Optimierung immer nur so gut wie die verwendete Zielfunktion. Hier zeigt sich das Problem, auf das man bei der Arbeit in diesem Gebiet immer wieder stößt. Je detaillierter der Ansatz zur Beschreibung einer Struktur ist, um so schwieriger ist es, die notwendigen Parameter für die Zielfunktion zu sammeln. Sind einfache Strukturen wie die von CsCl oder Halit noch mit Energiefunktionen zu beschreiben, die nur lokale Zusammenhänge zwischen den Atomen im Kristall beschreiben, so ist es für etwas kompliziertere Strukturen wie Rutil nur noch möglich, mit genauen Approximationen der Kräfte langer Reichweite, wie dem gleichnamigen Coulombterm, ein sinnvolles Ergebnis zu erzielen. Für eine qualitativ ”gute“ Vorhersage ist für diese Strukturen sogar ein massives Eingreifen in die Wichtung nötig. Eine Vorhersage dieser Struktur ist mit der von Pannetier et al. [25] vorgestellten Zielfunktion und der Minimum Image Convention nicht gelungen. Ein Beispiel dafür, daß diese Zielfunktionen nur ein Ansatz zur Bestimmung von Kristallstrukturen sein können, findet man mit der Quarzstruktur. Alle Versuche, diese mit dem Ansatz von Pannetier et al. [25] vorherzusagen schlagen fehl. Erst ein verallgemeinerter Ansatz, der detaillierterer Vorinformationen bedarf, führt hier zum Ziel. Diese Zielfunktion beschreibt eine größere Klasse von Kristallen, doch stellt sich das Finden der nötigen Parameter als problematisch dar. Eine weitere Grenze, an die eine Kristallstruktur-Vorhersage mit dem Simulated Annealing stößt, ist der enorme Bedarf an Rechenzeit. Mit dem Programm ANNEAL sind heute Elementarzellen, die über mehr als 30 Atome verfügen, kaum noch bearbeitbar.

Trotzdem geht die Tendenz weiter in die Richtung, Energiefunktionen so zu modifizieren, daß möglichst wenig Hintergrundinformationen für eine Struktur-Vorhersage nötig sind.

3.9 Ausblick

Zielfunktionen, wie sie von Pannetier et al. vorgestellt wurden, haben sich, vor allem nach den Verbesserungen des Coulomb-Terms durch die Methode der Ewald-Summen, bewährt. Man stößt aber auch schon für relativ ”regelmäßige“ Strukturen schnell an die Grenzen der Methode. An dieser Stelle hat sich der erweiterte Ansatz bewährt. Mit ihm war es möglich, komplexere Strukturen wie Quarz vorherzusagen, ohne erheblich mehr Vorwissen in die Struktur zu investieren. Die Parameter sind zwar häufig nicht ohne weiteres zu finden, doch könnte die Nebenbedingung, daß das Paulingsche Valenzprinzip erfüllt sein muß, den Weg zu dem folgenden Optimierungsproblem mit Nebenbedingungen frei machen: Sei x ∈ IR^{n} der Zustandsvektor des folgenden Problems P und seien a ∈ IR^{p} und b ∈ IR^{q} die zu bestimmenden Gilbertparameter.

P : min_{x,a,b}E_{coul}(x) + E_{rep}(x, a, b)

unter der Nebenbedingung E_{val}(x) = 0

Auf diese Weise wäre man nicht mehr auf die Kenntnis der Parameter a und b angewiesen. Mit einigen weiteren Modifikationen des Simulated Annealing-Algorithmus sollte dies ein erfolgversprechender Ansatz für weitere berlegungen sein.

Anhang

A.1 Die Optimierungsroutine SIMULA

Der folgende Abschnitt beschreibt die Verwendung der in Pascal geschriebenen Optimierungsroutine SIMULA. Die Implementation erfolgte für IBM RS-6000/RS-9000 und umfaßt die Algorithmen von Jens Radzek [26] und Frank Krimphove [17]. Ein Crosscompiler für Turbopascal ab Version 6 liegt bei. Die Anpassung an andere Compiler und Rechnersysteme sollte leicht durchführbar sein. SIMULA ist eine Routine zur Lösung von globalen, nichtlinearen Optimierungsproblemen der Form:

P : Minimiere f(x), wobei x = (x_{1}, ..., x_{n})^{T}

unter den Nebenbedingungen a ≤ x_{i} ≤ b für i ∈ {1, …, m}

Das Optimierungsverfahren basiert auf dem Prinzip des Simulated Annealing. Eine Beschreibung der Strategie des Simulated Annealing findet man in Kapitel 1 (siehe auch Aarts und Korst [2]).

Um SIMULA zur Lösung des Problems P einzusetzen, ist es nötig, daß der Benutzer eine Routine schreibt, welche die Funktion f für alle Werte berechnet, die von den Variablen angenommen werden können. Desweiteren ist es notwendig, daß der Benutzer ein Programm schreibt, in welchem SIMULA aufgerufen wird. Die Details werden später erläutert.

A.1.1 Verwendung der Optimierungsroutine SIMU-LA

Der Aufruf von SIMULA erfolgt über den Befehl:

simula(create_output,add_info,mode,ausgabedat,anzkomp,
       schranke_u,schranke_o,eps,startkrit,chi_Null,t_start,
       absenkung,k,delta_a,lambda,u,lmk,stopkrit,Stoplaenge,
       Eps_Stop,Delta_Stop,algorithmus,minstep,
       maxstep,ch,testzahl,zyklen,konst_gitter,tauschmode);

Alle Parameter des Aufrufs sind in Kleinbuchstaben zu belegen. Mögliche Werte sind kursiv gesetzt. Ihre Bedeutung ergibt sich aus der folgenden Tabelle oder aus Kapitel 1.

Variable  Typ  Bedeutung

create output boolean Ausgabedatei wird neu angelegt (true) oder es wird an eine bereits bestehende Datei angehängt (false).

add info boolean Es wird zusätzlicher Output über die Anzahl der Zustandsübergänge und die Anzahl der Funktionsauswertungen generiert (true).

mode str20 Steuerung der Bildschirmausgabe während des Programmablaufs. Keine Ausgabe: quiet. Ausgabe auf dem Schirm: verbose.

ausgabedat str80 Name und Pfad der Ausgabedatei. anzkomp integer Anzahl der Variablen.

schranke o  real  Obere Schranke für die Variablen.
schranke u  real  Untere Schranke für die Variablen.
eps  real  Startschrittweite des Algorithmus.
startkrit  str20  Startkriterium. Mögliche Werte:

kirkpatrick, johnson, white, selbst.

chi null real Parameter für Startkriterien

kirkpatrik und johnson.

t start real Starttemperatur für das Startkriterium

selbst.

absenkung str20 Absenkungsbedingung. Mögliche Werte:

konstante,lundy&mees,aarts&van laarhoven.

k  real  Konstante für Absenkungsbedingung konstant.
delta a  real  Parameter für Absenkungsbedingungen

lundy&mees und aarts&van laarhoven.

lambda  real  Parameter für Absenkungsbedingung huang.
u  real  Parameter für Absenkungsbedingung

lundy&mees.

lmk integer Länge einer Markovkette. stopkrit str20 Stopkriterium. Mögliche Werte:

kirkpatrick, huang, lundy&mees,

aarts&van laarhoven.

stoplaenge integer Länge der Stopkette für das Stopkriterium

kirkpatrick.

eps stop real ϵ für Stopkriterien

aarts&vanlaarhoven, huang,lundy&mees

delta stop  real  δ für Stopkriterium lundy&mees
algorithmus  str20  Algorithmusart. Mögliche Werte:

konstant (Algorithmus mit konstanter Schrittweite), variabel (Algorithmus

mit veränderlicher Schrittweite).

minstep real Minimale Schrittweite bei der Schrittweiten-

steuerung im variablen Algorithmus.

maxstep real Maximale Schrittweite bei der Schrittweiten-

steuerung im variablen Algorithmus.

ch real Parameter in der Schrittweitensteuerung testzahl integer Anzahl der versuchten Zustandsübergänge

pro Richtung in der Schrittweitensteuerung

zyklen integer Anzahl der Durchläufe durch alle

Koordinatenrichtungen (siehe Diplomarbeit).

konstgitter str20 Art der Abschlußverfeinerung. Mögliche Werte: keine, random abstieg nach Radzek, konst abstieg, steepest descend, dfp verfahren.

tauschmode str20 Art der Zustandsübergänge im Algorithmus. Mögliche Werte: fix, sphaere, variabel.

Neben diesen Variablen müssen im Hauptprogramm noch einige Typen und Konstanten deklariert werden:

Konstante  Bedeutung
maxkomp  Maximale Zahl von

Variablen

Typ  Beschreibung  Bedeutung
str20  string(20)  String der Länge 20
str80  string(80)  String der Länge 80
vektor  array [1..maxkomp] of real Vektor für die Variablen

In der PC Version von SIMULA ist zu beachten, daß alle real Variablen durch double Variablen zu ersetzen sind. Dies ist aufgrund der Tatsache, daß der Pascal-Compiler der RS 6000/9000 real Variablen im 8-Byte-Format ablegt, notwendig. Dies entspricht dem Variablen Typ double unter Turbopascal.

A.1.2 Vom Benutzer zu erstellende Routinen

Die vom Benutzer zu schreibende Routine, welche den Funktionswert von f zum aktuellen Vektor x = (x_{1}, …, x_{n}) ermittelt, ist auf die folgende Art und Weise zu erzeugen:

 function kosten(x:vektor):real;

Steht eine solche Routine zur Verfügung, so kann eine übergeordnete Routine, in welcher die Parameter von SIMULA korrekt belegt sind und in die SIMULA eingebunden ist, die Optimierungsroutine aufrufen. Die folgende Routine stellt ein Beispiel für die Verwendung der Optimierungsroutine dar.

program testsimula;
const maxkomp     = 30;        (* maximale Anzahl der Variablen            *)
type  str20 = string(20);
 str80 = string(80);
 vektor = array[1..maxkomp] of real;
var
ausgabedat     : string;  (* Name der Ausgabedatei                    *)
mode           : str20;   (* verbose/quiet fuer Hintergrundaufruf     *)
anzkomp        : integer; (* Anzahl der Komponenten im Teilchenvek.   *)
schranke_o     : real;    (* obere Schranke fuer Optimierung          *)
schranke_u     : real;    (* untere Schranke fuer Optimierung         *)
eps            : real;    (* Parameter Startschrittweite              *)
t_start        : real;    (* Startwert t                              *)
startkrit      : str20;   (* kirkpatrick,johnson,white,selbst         *)
chi_null       : real;    (* fuer Startkrit.: kirkpatrick und Johnson *)
absenkung      : str20;   (* konstante,lundy&mees,aarts&van_laarhoven,*)
  (* oder huang je nach Absenkung von T       *)
k              : real;    (* fuer Absenkung mit Konstanter            *)
u              : real;    (* fuer Absenkung lundy&mees                *)
delta_a        : real;    (* fuer Absenkung lundy&mees und            *)
  (* aarts&van_laarhoven                      *)
lambda         : real;    (* fuer Absenkung huang                     *)
tauschmode     : str20;   (* fix,variabel,sphaere je nach Austausch   *)
  (* modus fix entlang einer Achse            *)
  (*       variabel entlang einer Achse       *)
  (*       auf fixer Kugelsphaere             *)
stopkrit       : str20;   (* kirkpatrick, huang, lundy&mees,          *)
  (* aarts&van_laarhoven je nach Stopkriterium*)
stoplaenge     : integer; (* fuer Stopkrit. kirkpatrick               *)
eps_stop       : real;    (* eps des jeweiligen Stopkriteriums        *)
delta_stop     : real;    (* delta des jeweiligen Stopkriteriums      *)
konst_gitter   : str20;   (* Abschlussverfeinerung                    *)
  (* vgl. Kapitel 5                           *)
lmk            : integer; (* Laenge der Mks : Anzahl der Nachbarn     *)
testzahl       : integer; (* Anzahl der Tests pro Koordinaten-        *)
  (* richtung (Krimphove)                     *)
ch             : real;    (* Parameter Schrittweitensteuerung         *)
zyklen         : integer; (* Anzahl der Durchlaeufe durch alle        *)
  (* Koordinatenrichtungen                    *)
algorithmus    : str20;   (* Art des Simulated Annealing-Algorithmus: *)
  (* radzek/krimphove                         *)
minstep,maxstep: real;    (* minimale,maximale Schrittweite bei       *)
  (* Veraenderung der Schrittweite            *)
create_output  : boolean; (* Output Datei neu anlegen ja/nein         *)
add_info       : boolean; (* Zusaetzliche Infos ausgeben ja/nein      *)

function kosten(x:vektor): real;
(* Berechnet den Funktionswert der Rastringin Funktion *)
begin (* kosten *)
kosten:=x[1]*x[1]+x[2]*x[2]-cos(18*x[1])-cos(18*x[2]);
end; (* kosten *)
%INCLUDE ./SIMULA.PAS
begin;
create_output:=true;
add_info:=true;
mode:='quiet';
ausgabedat:='test.out';
anzkomp:=2;
schranke_o:=1;
schranke_u:=-1;
eps:=0.1;
startkrit:='kirkpatrick';
chi_null:=0.95;
absenkung:='konstante';
k:=0.95;
lmk:=20;
stopkrit:='kirkpatrick';
stoplaenge:=10;
algorithmus:='konstant';
minstep:=0.00001;
maxstep:=0.5;
ch:=1.5;
testzahl:=10;
zyklen:=10;
konst_gitter:='random_abstieg';
tauschmode:='sphaere';
(* 1. Aufruf SIMULA mit neuer Datei *)
simula(create_output,add_info,mode,ausgabedat,anzkomp,schranke_u,schranke_o,
eps,startkrit,chi_Null,t_start,absenkung,k,delta_a,lambda,u,

  lmk,stopkrit,Stoplaenge,Eps_Stop,Delta_Stop,algorithmus,minstep,
  maxstep,ch,testzahl,zyklen,konst_gitter,tauschmode);
end.

A.2 Das Kristallstruktur-Vorhersage -Pro-gramm ANNEAL

Das Programm ANNEAL basiert auf der Optimierungsroutine SIMULA. Es dient der Vorhersage von Kristallstrukturen auf der Basis von Energiefunktionen. Die Eingabe der Parameter für die einzelnen Kristalle erfolgen über ein Datenfile. Die verschiedenen Zielfunktionen aus Kapitel 2 stehen im Programm zur Verfügung.

A.2.1 Verfügbare Versionen

Es gibt zwei Versionen:

1. Eine Grundversion ANNEAL (lauffähig auf IBM RS6000/IBM RS9000) unter XLP (Risc Pascal).

2. Eine crosscompilierte Version ANNEALPC für IBM-kompatible PC (lauffähig unter TURBOPASCAL ab Version 6).

Anweisungen zum Compilieren der beiden Versionen und ein Crosscompiler für MSDOS Rechner befinden sich in den entsprechenden Unterverzeichnissen auf der Diskette.

A.2.2 Maximale Anzahl von Atomen

Die maximale Anzahl von Atomen in der Elementarzelle liegt für die PC-- Version ANNEALPC.EXE bei 20. Sollen größere Probleme gerechnet werden, so ist man gezwungen, die UNIXVERSION von ANNEAL zu verwenden. Diese ist im Moment für bis zu 50 Atome pro Elementarzelle kompiliert. Die maximale Anzahl von Atomen pro Zelle wird von der Konstanten MAX-TEILCHEN im Programmkopf von ANNEAL.PAS bestimmt.

A.2.3 Syntax des Programmaufrufes

Ist aus dem Sourcecode ein ausführbarer Code erzeugt worden, so ist die Syntax für den Aufruf von ANNEAL und ANNEALPC die folgende:

ANNEAL EINGABEDATEI AUSGABEDATEI
ANNEALPC EINGABEDATEI AUSGABEDATEI

A.2.4 Aufbau eines Datenfiles für ANNEAL

Beispiele für Ein- und Ausgabedateien befinden sich in Unterverzeichnissen auf der beiliegenden Diskette. Generell sieht eine Zeile der Eingabedatei folgendermaßen aus:

 30 Zeichen Kommentar          Eingetragene Option

Zwischen den Eingabezeilen tauchen immer wieder einzelne Kommentarzeilen oder Blöcke von Kommentarzeilen auf. Sie haben das Format:

 60 Zeichen beliebiger Kommentar

Die Anzahl und Abfolge von Options- und Kommentarzeilen ist festgelegt. Als Beispiel für eine Eingabedatei sei das Datenfile zur Vorhersage von Zinkblende vorgestellt. Diese Datei läßt sich auch gut als Formular für eigene Datenfiles benutzen.

KOMMENTARZEILE
KOMMENTARZEILE     Datenfile f"ur ZnS 08.05.94
KOMMENTARZEILE
************** Steuerungsanweisungen ***********************
Mode:                        verbose
************************************************************
KOMMENTARZEILE
KOMMENTARZEILE     Parameter fuer das Simulated Annealing
KOMMENTARZEILE
************** A. Algorithmusart: **************************
Algorithmus:                 konstant
************** A.1. Schrittweitensteuerung *****************
Minimale Schrittweite:       -
Maximale Schrittweite:       -
ch:                          -
************** A.2. Anzahl der Tests fuer Schrittweite *****
Anzahl der Tests/Zyklus:     -
Anzahl der Zyklen:           -
************** B. Startkriterium und Aufheizakzeptanzrate :*
Startkriterium:              kirkpatrick
Chi_Null:                    0.98
Starttemperatur:             -
************** C. Art der Abkuehltabelle mit Parametern : **
Absenkung von T:             konstante
k:                           0.95
Delta_a:                     -
Lambda:                      -
u:                           -
************** D. Stopkriterium mit Parametern : ***********
Stopkriterium:               kirkpatrick
Stoplaenge:                  10
Eps_Stop:                    -
Delta_Stop:                  -
************** E. Laenge der Markovkette : *****************
Laenge der Markowkette:      20
************** F. Gitterart und Abschlussverfeinerung: *****
Gittererart:                 sphaere3d
Abschlussverfeinerung:       steepest_descend
************** G. Startschrittweite : **********************
Startschrittweite (Eps0):    0.10
************************************************************
KOMMENTARZEILE
KOMMENTARZEILE    Zielfunktionsart und Wichtungen :
KOMMENTARZEILE
************** i. Zielfunktion : ***************************
Zielfunktion=                modified
************** i.a. Parameter fuer Pannetier Zielfunktion **
Best. der Wichtung Lambda:   fix
Eingabe der Wichtung (fix):  0.50
Anzahl der Tests (variabel): -
Anzahl der Zeilen von i.b.): 8
************** i.b. Definition Zielfunktion : **************
Anzahl der Terme:            2
Term1:                       elektroterm=
Wichtung Term1:              1

Term2:                       valenzterm+
Wichtung Term2:              1
************** iii. Summationsart Coulombterm **************
Summationsart:               ewald_barteaud
************************************************************
KOMMENTARZEILE
KOMMENTARZEILE    Kristallographische Konstanten :
KOMMENTARZEILE
************** 1. Zellparameter : **************************
a:                           5.4100
b:                           5.4100
c:                           5.4100
Alpha:                       90.000
Beta:                        90.000
Gamma:                       90.000
************** 2. Problemumfang : **************************
Anzahl der Teilch pro Zelle: 8
************** 4. Koordinationen und Teilchenladungen : ****
Name Teilchen 1:             Zn
Formale Ladung Teilchen1:    2
Koordinationszahl Teilchen1: 4
Name Teilchen 2:             Zn
Formale Ladung Teilchen2:    2
Koordinationszahl Teilchen2: 4
Name Teilchen 3:             Zn
Formale Ladung Teilchen3:    2
Koordinationszahl Teilchen3: 4
Name Teilchen 4:             Zn
Formale Ladung Teilchen4:    2
Koordinationszahl Teilchen4: 4
Name Teilchen 5:             S
Formale Ladung Teilchen5:    -2
Koordinationszahl Teilchen5: 4
Name Teilchen 6:             S
Formale Ladung Teilchen6:    -2
Koordinationszahl Teilchen6: 4
Name Teilchen 7:             S
Formale Ladung Teilchen7:    -2
Koordinationszahl Teilchen7: 4
Name Teilchen 8:             S
Formale Ladung Teilchen8:    -2
Koordinationszahl Teilchen8: 4
************** 6. Matrix der optimalen Abstaende (R0) : ****
Anzahl:                      1
Name 1. Teilchen:            Zn
Name 2. Teilchen:            S
Optimaler Abstand (R0):      2.086
************** 7. Matrix der optimalen Abstaende (D0) : ****
Anzahl:                      2
Name:                        Zn
Atomradius (nach Gilbert):   -
Softness Parameter (Gilb.):  -
Erfolgt Fit der Parameter:   nein
Name:                        S
Atomradius (nach Gilbert):   -
Softness Parameter (Gilb.):  -
Erfolgt Fit der Parameter:   nein

A.2.5 Erläuterung der Parameter

Die Eingabe der Parameter erfolgt (wie oben angegeben) \underline{immer} ab Spalte 30. Alle Parameter müssen mit kleinen Buchstaben belegt werden. Für die Wahl der Parameter gibt es folgende Möglichkeiten:

0. Steuerungsanweisungen.

Variable  Mögliche Werte  Bedeutung

mode  quiet  keine Bildschrimausgabe

während des

Programmablaufs.

verbose Zwischeninformationen werden zur Ausführungszeit mitprotokolliert.

A. Algorithmusart.

Variable  Mögliche Werte  Bedeutung
algorithmus  konstant  Algorithmus mit konstanter

Schrittweite nach Radzek [26].

variabel Algorithmus mit variabler Schrittweite nach Krimphove [17] (nur für Probleme mit sehr wenigen Variablen empfehlenswert, da sehr rechenintensiv).

A.1. Schrittweitensteuerung. Wird im Falle des konstanten Algorithmus ignoriert!

Variable Bedeutung

maximale Schrittweite Maximale Schrittweite bei der Schrittweiten-

steuerung.

minimale Schrittweite Minimale Schrittweite bei der Schrittweiten-

steuerung.

ch Steuerungsparameter für Schrittweitenberechnung.

A.2. Anzahl der Tests für Schrittweite. Wird im Falle des konstanten Algorithmus ignoriert!

Variable Bedeutung

Anzahl Tests pro Zyklus Anzahl der zufällig versuchten Zustandsüber-

gänge pro Koordinatenrichtung.

Anzahl der Zyklen Anzahl der Durchläufe durch alle Koordinatenrichtungen.

B. Startkriterium. Zum testen ist die Eintragung kirkpatrick empfehlenswert. Ist für ein Datenfile dann eine ungefähre Starttemperatur bekannt, kann man diese mit der Eintragung selbst konkret vorgeben, das spart Zeit.

Variable  Mögliche Werte  Bedeutung
Startkriterium  kirkpatrick  Startkriterium nach Kirk-

-patrick [14]. white Startkriterium nach White [30]. johnson Startkriterium nach

Johnson [12].

selbst Fest vorgegebene Start-

temperatur.

k 0.85 … 0.99 Kontrollparameter für die Startwertbestimmung. Wird ein Prozentsatz von k · 100 zufälligen Zustandsübergängen akzeptiert, so ist eine geeignete Starttemperatur gefunden.

Starttemperatur Eigene Vorgabe der Starttemperatur. Wird außer für Startkriterium selbst ignoriert.

C. Art der Abkühltabelle. Die konstante Absenkung bringt für erste Versuche mit dem Programm schon brauchbare Ergebnisse.

Variable  Mögliche Werte  Bedeutung
Absenkung von T  konstante  Absenkung Multiplikation

mit einer Konstanten.

lundy&mees Absenkung nach Lundy und

Mees [19].

aarts&van laarhoven Absenkung nach Aarts und

van Laarhoven [1].

huang Absenkung nach Huang.

k 0.85 … 0.97 Bei konstanter Absenkung der Multiplikator für die Temperatur.

Delta a Parameter bei Lundy & Mees und Aarts & van Laarhoven.

Lambda Parameter bei Huang. u Parameter bei Lundy

& Mees.

D. Stopkriterium.

Variable  Mögliche Werte  Bedeutung
Stopkriterium  kirkpatrick  Stopkriterium nach

Kirkpatrick [14]. lundy&mees Stopkriterium nach Lundy

& Mees [19].

aarts&van laarhoven Stopkriterium nach Aarts

& van Laarhoven [1].

huang Stopkriterium nach Huang [11].

Stoplänge 5 … 15 Länge der Stopkette

bei Kirkpatrick.

Eps Stop Parameter bei Aarts & van Laarhven, Lundy & Mees sowie Huang.

Delta Stop Parameter bei Lundy & Mees.

E. Länge einer Markovkette.

Variable Bedeutung

Länge der Markovkette In der Regel sagt man: LMK ≤ Anzahl der Variablen. Diese Schätzung ist jedoch rein empirisch.

F. Gitterart und Abschlußverfeinerung.

Art der Nachbarnbestimmung und Art der Abschlußverfeinerung.

Variable  Mögliche Werte  Bedeutung
Gitterart  variabel  Variable Schrittweite (nur

variabler Algorithmus).

sphaere3d Konstante Schrittweite vom Mittelpunkt einer 3-dimen-sionalen Kugel auf einen Punkt auf ihrer Hülle.

fix Feste Schrittweite entlang

der Achsen.

Abschlußverfeinerung keine Keine Abschlußverfeinerung.

random abstieg Randomisiertes Abstiegsverfahren nach Radzek [26]. Dieses hat sich nicht bewährt, da zu viele Funktionsauswertungen nötig sind.

dfp verfahren DFP-Verfahren (nur dif-

ferenzierbare Probleme).

konst abstieg Einfaches sequentielles Ab-

stiegsverfahren.

steepest descend Verfahren des steilsten Abstiegs.

Die letzten beiden Verfahren der Abschlußverfeinerung haben sich für die bisherigen Zielfunktionen bei der Kristallstruktur-Vorhersage bewährt.

G. Startschrittweite.

Variable Bedeutung

Startschrittweite (EPS0) Startschrittweite für das Simulated Annealing. Sollte nicht zu klein gewählt werden, da sonst keine Konvergenz erreicht wird. Bei bis zu 30 Variablen (10 Teilen in der Elementarzelle) hat sich ein Wert von 0.1 bewährt.

i. Zielfunktion.

Variable  Mögliche Werte  Bedeutung
Art der Zielfunktion  pannetier  Zielfunktion nach Pannetier

et al [25]. Die Art der Wichtung wird in Punkt i.a.) festgelegt.

modified Selbst definierte Zielfunktion. Die Definition erfolgt in Punkt i.b.). Auch die erweitete Zielfunktion kann hier definiert werden.

i.a. Parameter für die Pannetier-Zielfunktion. Diese Eingaben werden im Falle der modifizierten (erweiterten) Zielfunktion übergangen!

Variable  Mögliche Werte  Bedeutung
Definition Wichtung  fix  Eigene Vorgabe des Kon-

trollparameters Lambda.

variabel Ermittelung von Lambda nach einer vorgegebenen Anzahl von Funktionsauswertungen mit der Formel von Pannetier (hat bis jetzt nie brauchbare Ergebnisse geliefert).

Eingabe der Wichtung Im Fall der festen Vorgabe von Lambda wird der Wert hier eingegeben.

Anzahl der Tests  100 … 200  Im Fall der variablen Er-

mittlung von Lambda wird hier die Anzahl der zu berechnenden Funktionsauswertungen vorgegeben.

Anz. der Zeilen v. i.b.) Anzahl der in i.b.) vorhandenen Zeilen mit berschrift. Ist als Zielfunktion Pannetier vorgegeben, müssen diese überlesen werden.

i.b. Definition der Zielfunktion. Wird im Fall der Pannetier-Zielfunktion überlesen!

Variable  Mögliche Werte  Bedeutung
Anzahl der Terme  1 … 3  Anzahl der folgenden Ziel-

funktionterme. Diese bestehen jeweils aus zwei Parametern:

Termart elektroterm<> Coulomb-Term für ungleich

geladene Atome.

elektroterm= Coulomb-Term für gleich ge-

ladene Atome.

elektroterm Gesamter Coulomb-Term für

alle Atome.

valenzterm Valenzterm ohne Vorgabe

der Koordinationszahlen.

valenzterm+ Valenzterm mit Vorgabe der

Koordinationszahlen.

ionenradienterm Gilbertscher Abstoßungs-

term.

Wichtung Wichtung der Terme in der Zielfunktion.

ii. Summationsart Coulombterm.

Variable  Mögliche Werte  Bedeutung
Summationsart  minimum image  Berechnung des Coulomb-

terms nach der ”Minimum Image Convention for the Truncation of the Potential“. ewald berteaud Berechnung des Coulomb-

terms nach der Ewald-Bertaud-Williams Methode.

1. Zellparameter. Eingabe der Zellkonstanten in Å. Die Eingabe der Winkel erfolgt in Grad.

2. Problemumfang. Eingabe der Anzahl der Atome in der Elementarzelle.

3. Koordination und Teilchenladung. Für alle Atome wird jeweils ein Block von drei Zeilen eingelesen.

Variable Bedeutung

Name des Atoms Dieser dient als Identifikationsmerkmal und darf

zwischen 1 und 5 Buchstaben lang sein!

Formale Ladung Formale Valenz des Atoms. Koordinationszahl Koordinationszahl des Atoms.

4. Matrix der optimalen Abstände R0. Eingabe der r0 Werte für die Paulingsche Valenzregel. Jeder r0–Block besteht aus drei Zeilen.

Variable Bedeutung

Anzahl Anzahl der folgenden r0–Blöcke. Name Atom 1 Name des ersten an der Bindung beteiligten

Atoms.

Name Atom 2 Name des zweiten, an der Bindung beteiligten

Atoms.

r0–Wert r0–Wert der Bindung (nach Busing [7]).

5. Matrix der optimalen Abstände D0. Eingabe der Parameter für den Gilbert–Term. Dieser Abschnitt wird außer bei der Anwahl des Ionenradiusterms für die Zielfunktionsart modified überlesen. Jeder D0–Block besteht aus fünf Zeilen.

Variable Bedeutung

Anzahl Anzahl der folgenden D0–Blöcke.

Name Atom Name des Atoms, für welches die Parameter

folgen.

Radius Atom Gilbert–Atomradius des Atoms (In der

Formel: A).

Softness Atom Gilbert–Softness–Parameter des Atoms

(In der Formel: B).

Fit Fit (ja/nein). Hier besteht die Möglichkeit, einen Abschließenden Fit der Gilbert-Parameter zu aktivieren. Leider ist diese Prozedur sehr zeitaufwendig, hat sich jedoch beispielsweise für Quarz bewährt.

A.3 Beschreibung einer Ausgabedatei von ANNEAL

In diesem Abschnitt soll das Aussehen einer Ausgabedatei der Kristallstruktur-Vorhersage anhand der Ausgabedatei für Zinkblende erläutert werden.

Das Protokoll der Eingabedatei

Bei der Arbeit mit der Kristallstruktur-Vorhersage hat es sich als nützlich erwiesen, daß eine Kopie der Eingabedatei mit in die Ausgabedatei aufgenommen wird, da kristallographische Parameter und Wichtungen häufig geändert wurden und nach einiger Zeit nicht mehr nachvollziehbar war, welche Ausgaben mit welchen Eingaben erzeugt wurden. Die Zeilen, die das Protokoll der Eingabedatei enthalten, sind mit der Kennung PROTO> gekennzeichnet.

Informationen zum Ablauf der Optimierung

Nach dem Protokoll der Eingabedatei werden die Informationen über den Ablauf der Optimierung in die Ausgabedatei geschrieben. Hierbei handelt es sich im Wesentlichen um die Standardausgaben der Optimierungsroutine SIMULA, die jedoch, um das Lesen der Informationen zu vereinfachen, etwas erweitert wurden. So werden die Werte der einzelnen Zielfunktionselemente einzeln ausgegeben und die Ausgabe des Zustandsvektors erfolgt zusammen mit den Kennungen für die einzelnen Atome gruppiert zu je drei Variablen.

Ausgabe der Bindungsinformationen

Den Abschluß der Ausgabedatei bilden die Bindungsinformationen der einzelnen Koordinationspolyeder. Die Ausgabe besteht aus Datenblöcken, die jeweils ein Koordinationspolyeder mit allen Abständen und Winkeln darstellen. Das zugehörige Zentralatom wird im Kopf eines jeden Blocks aufgeführt.

Die Ausgabedatei von Zinkblende

**********************************************************************
*   PROTOKOLL DER EINGABEDATEI DES ALGORITHMUS !!!                   *

**********************************************************************
PROTO>
PROTO>    Datenfile fuer ZnS 08.05.94
PROTO>
PROTO>    ************** Steuerungsanweisungen ***********************
PROTO>    Mode:                        verbose
PROTO>    ************************************************************
PROTO>
PROTO>        Parameter fuer das Simulated Annealing
PROTO>
PROTO>    ************** A. Algorithmusart: **************************
PROTO>    Algorithmus:                 konstant
PROTO>    ************** A.1. Schrittweitensteuerung *****************
PROTO>    Minimale Schrittweite:       -
PROTO>    Maximale Schrittweite:       -
PROTO>    ch:                          -
PROTO>    ************** A.2. Anzahl der Tests fuer Schrittweite *****
PROTO>    Anzahl der Tests/Zyklus:     -
PROTO>    Anzahl der Zyklen:           -
PROTO>    ************** B. Startkriterium und Aufheizakzeptanzrate :*
PROTO>    Startkriterium:              kirkpatrick
PROTO>    Chi_Null:                    0.98
PROTO>    Starttemperatur:             -
PROTO>    ************** C. Art der Abkuehltabelle mit Parametern : **
PROTO>    Absenkung von T:             konstante
PROTO>    k:                           0.95
PROTO>    Delta_a:                     -
PROTO>    Lambda:                      -
PROTO>    u:                           -
PROTO>    ************** D. Stopkriterium mit Parametern : ***********
PROTO>    Stopkriterium:               kirkpatrick
PROTO>    Stoplaenge:                  10
PROTO>    Eps_Stop:                    -
PROTO>    Delta_Stop:                  -
PROTO>    ************** E. Laenge der Markovkette : *****************
PROTO>    Laenge der Markowkette:      20
PROTO>    ************** F. Gitterart und Abschlussverfeinerung: *****
PROTO>    Gittererart:                 sphaere3d
PROTO>    Abschlussverfeinerung:       steepest_descend
PROTO>    ************** G. Startschrittweite : **********************
PROTO>    Startschrittweite (Eps0):    0.10
PROTO>    ************************************************************
PROTO>
PROTO>        Zielfunktionsart und Wichtungen :
PROTO>
PROTO>    ************** i. Zielfunktion : ***************************
PROTO>    Zielfunktion=                modified
PROTO>    ************** i.a. Parameter fuer Pannetier Zielfunktion **
PROTO>    Best. der Wichtung Lambda:   fix
PROTO>    Eingabe der Wichtung (fix):  0.50
PROTO>    Anzahl der Tests (variabel): -
PROTO>    Anzahl der Zeilen von i.b.): 8
PROTO>    ************** i.b. Definition Zielfunktion : **************
PROTO>    Anzahl der Terme:            2
PROTO>    Term1:                       elektroterm=
PROTO>    Wichtung Term1:              1
PROTO>    Term2:                       valenzterm+
PROTO>    Wichtung Term2:              1
PROTO>    ************** iii. Summationsart Coulombterm **************
PROTO>    Summationsart:               ewald_barteaud
PROTO>    ************************************************************
PROTO>
PROTO>         Kristallographische Konstanten :
PROTO>
PROTO>    ************** 1. Zellparameter : **************************
PROTO>    a:                           5.4100
PROTO>    b:                           5.4100
PROTO>    c:                           5.4100
PROTO>    Alpha:                       90.000
PROTO>    Beta:                        90.000
PROTO>    Gamma:                       90.000
PROTO>    ************** 2. Problemumfang : **************************
PROTO>    Anzahl der Teilch pro Zelle: 8
PROTO>    ************** 4. Koordinationen und Teilchenladungen : ****
PROTO>    Name Teilchen 1:             Zn
PROTO>    Formale Ladung Teilchen1:    2
PROTO>    Koordinationszahl Teilchen1: 4
PROTO>    Name Teilchen 2:             Zn
PROTO>    Formale Ladung Teilchen2:    2
PROTO>    Koordinationszahl Teilchen2: 4
PROTO>    Name Teilchen 3:             Zn
PROTO>    Formale Ladung Teilchen3:    2
PROTO>    Koordinationszahl Teilchen3: 4
PROTO>    Name Teilchen 4:             Zn
PROTO>    Formale Ladung Teilchen4:    2
PROTO>    Koordinationszahl Teilchen4: 4
PROTO>    Name Teilchen 5:             S
PROTO>    Formale Ladung Teilchen5:    -2
PROTO>    Koordinationszahl Teilchen5: 4
PROTO>    Name Teilchen 6:             S
PROTO>    Formale Ladung Teilchen6:    -2

PROTO>    Koordinationszahl Teilchen6: 4
PROTO>    Name Teilchen 7:             S
PROTO>    Formale Ladung Teilchen7:    -2
PROTO>    Koordinationszahl Teilchen7: 4
PROTO>    Name Teilchen 8:             S
PROTO>    Formale Ladung Teilchen8:    -2
PROTO>    Koordinationszahl Teilchen8: 4
PROTO>    ************** 6. Matrix der optimalen Abstaende (R0) : ****
PROTO>    Anzahl:                      1
PROTO>    Name 1. Teilchen:            Zn
PROTO>    Name 2. Teilchen:            S
PROTO>    Optimaler Abstand (R0):      2.086
PROTO>    ************** 7. Matrix der optimalen Abstaende (D0) : ****
PROTO>    Anzahl:                      2
PROTO>    Name:                        Zn
PROTO>    Atomradius (nach Gilbert):   -
PROTO>    Softness Parameter (Gilb.):  -
PROTO>    Erfolgt Fit der Parameter:   nein
PROTO>    Name:                        S
PROTO>    Atomradius (nach Gilbert):   -
PROTO>    Softness Parameter (Gilb.):  -
PROTO>    Erfolgt Fit der Parameter:   nein
**********************************************************************
*   HOMOGENER SIMULATED ANNEALING ALGORITHMUS !!!                    *
**********************************************************************
Eingabedatei                     : .\..\input\zns.ann
Anzahl der Variablen             :   24
Komponenten des Vektors zwischen  0.000 und  1.000
Die Startschrittweite war        : 0.1000
Bewegungen im Gitter             : auf 3 dimensionalen Kugelsphaeren
                                   mit konstanten Radien
**********************************************************************
Laenge einer Markovkette         :    20
Bestimmung des Startwerts der Temperatur T nach  Kirkpatrick
Startwert war                    : 0.00000000
Absenkung der Temperatur durch Multiplikation mit einer Konstanten
Konstante                        : 0.950
Bestimmung Stopkriteriums nach Kirkpatrick
**********************************************************************
Startvektor                          :
 Zn x[  1]=   0.00000 x[  2]=   0.00000 x[  3]=   0.00000
 Zn x[  4]=   0.60000 x[  5]=   0.30000 x[  6]=   0.60000
 Zn x[  7]=   0.30000 x[  8]=   0.80000 x[  9]=   0.60000
 Zn x[ 10]=   0.80000 x[ 11]=   0.10000 x[ 12]=   0.40000
 S x[ 13]=   0.60000 x[ 14]=   0.10000 x[ 15]=   0.50000
 S x[ 16]=   0.80000 x[ 17]=   0.40000 x[ 18]=   0.10000
 S x[ 19]=   0.00000 x[ 20]=   0.10000 x[ 21]=   0.70000
 S x[ 22]=   0.90000 x[ 23]=   0.70000 x[ 24]=   0.40000

**********************************************************************
Stopwert von t war               : 0.0000261626535585
**********************************************************************
Funktionswert Simulated Annealing: 133.45822328
Funktionsauswertungen S. A.      : 8501
Anzahl der Zustandsaenderungen   : 2344
**********************************************************************
Ergebnisvektor Sim. Annealing    :
 Zn x[  1]=   0.00000 x[  2]=   0.00000 x[  3]=   0.00000
 Zn x[  4]=   0.49995 x[  5]=   0.00006 x[  6]=   0.49996
 Zn x[  7]=   1.00003 x[  8]=   0.49990 x[  9]=   0.50023
 Zn x[ 10]=   0.49987 x[ 11]=   0.49997 x[ 12]=   0.00003
 S x[ 13]=   0.75384 x[ 14]=   0.74688 x[ 15]=   0.75369
 S x[ 16]=   0.75379 x[ 17]=   0.24706 x[ 18]=   0.25375
 S x[ 19]=   0.25368 x[ 20]=   0.74680 x[ 21]=   0.25378
 S x[ 22]=   0.25395 x[ 23]=   0.24655 x[ 24]=   0.75395

**********************************************************************
Anzahl der Funktionsauswertungen : 10968
Anzahl der Zustandsaenderungen   : 2344

Zusammensetzung der Zielfunktion :
Modifizierte Zielfunktion :
Term[1] = elektroterm=  Wichtung  :  1.00000
Term[2] = valenzterm+  Wichtung  :  1.00000
**********************************************************************
Der optimale Funktionswert ist   : 118.75973055
Valenzterm                       :   0.00865627
**********************************************************************
Der Ergebnisvektor lautet        :
Zn x[  1]=   0.00000 x[  2]=   0.00000 x[  3]=   0.00000
Zn x[  4]=   0.49995 x[  5]=   0.00006 x[  6]=   0.49996
Zn x[  7]=   1.00003 x[  8]=   0.49990 x[  9]=   0.50023
Zn x[ 10]=   0.49987 x[ 11]=   0.49997 x[ 12]=   0.00003
S x[ 13]=   0.75384 x[ 14]=   0.74688 x[ 15]=   0.75369
S x[ 16]=   0.75379 x[ 17]=   0.24706 x[ 18]=   0.25375
S x[ 19]=   0.25368 x[ 20]=   0.74680 x[ 21]=   0.25378
S x[ 22]=   0.25395 x[ 23]=   0.24655 x[ 24]=   0.75395
Koordinationspolyeder            :
CN [1]              = 4
CN [2]              = 4
CN [3]              = 4
CN [4]              = 4
CN [5]              = 4
CN [6]              = 4
CN [7]              = 4
CN [8]              = 4
**********************************************************************
*       ---     WINKEL UND ABSTANDSBERECHNUNG              ---       *
**********************************************************************
**********************************************************************
Analyse der Koordinationspolyeder Teilchen Nr. 1 (Zn):
1. Ligand: Nr.: 5 Name:   S (-0.246,-0.253,-0.246) Abstand : 2.32903
2. Ligand: Nr.: 8 Name:   S ( 0.254, 0.247,-0.246) Abstand : 2.33205
3. Ligand: Nr.: 6 Name:   S (-0.246, 0.247, 0.254) Abstand : 2.33349
4. Ligand: Nr.: 7 Name:   S ( 0.254,-0.253, 0.254) Abstand : 2.37589
Winkel: (1.Ligand,Zentralatom,2.Ligand) = (  S, Zn,  S)    : 110.278
Winkel: (1.Ligand,Zentralatom,3.Ligand) = (  S, Zn,  S)    : 110.303
Winkel: (2.Ligand,Zentralatom,3.Ligand) = (  S, Zn,  S)    : 110.149
Winkel: (1.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.779
Winkel: (2.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.620
Winkel: (3.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.664
**********************************************************************
Analyse der Koordinationspolyeder Teilchen Nr. 2 (Zn):
1. Ligand: Nr.: 7 Name:   S ( 0.254,-0.253, 0.254) Abstand : 2.32939
2. Ligand: Nr.: 8 Name:   S ( 0.254, 0.247, 0.754) Abstand : 2.33189
3. Ligand: Nr.: 6 Name:   S ( 0.754, 0.247, 0.254) Abstand : 2.33361
4. Ligand: Nr.: 5 Name:   S ( 0.754,-0.253, 0.754) Abstand : 2.37634
Winkel: (1.Ligand,Zentralatom,2.Ligand) = (  S, Zn,  S)    : 110.294
Winkel: (1.Ligand,Zentralatom,3.Ligand) = (  S, Zn,  S)    : 110.307
Winkel: (2.Ligand,Zentralatom,3.Ligand) = (  S, Zn,  S)    : 110.166
Winkel: (1.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.775
Winkel: (2.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.608
Winkel: (3.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.641
**********************************************************************
Analyse der Koordinationspolyeder Teilchen Nr. 3 (Zn):
1. Ligand: Nr.: 6 Name:   S ( 0.754, 0.247, 0.254) Abstand : 2.32890
2. Ligand: Nr.: 5 Name:   S ( 0.754, 0.747, 0.754) Abstand : 2.33225
3. Ligand: Nr.: 7 Name:   S ( 1.254, 0.747, 0.254) Abstand : 2.33347
4. Ligand: Nr.: 8 Name:   S ( 1.254, 0.247, 0.754) Abstand : 2.37696
Winkel: (1.Ligand,Zentralatom,2.Ligand) = (  S, Zn,  S)    : 110.272
Winkel: (1.Ligand,Zentralatom,3.Ligand) = (  S, Zn,  S)    : 110.208
Winkel: (2.Ligand,Zentralatom,3.Ligand) = (  S, Zn,  S)    : 110.110
Winkel: (1.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.817
Winkel: (2.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.716
Winkel: (3.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.673

Analyse der Koordinationspolyeder Teilchen Nr. 4 (Zn):
1. Ligand: Nr.: 8 Name:   S ( 0.254, 0.247,-0.246) Abstand : 2.32854
2. Ligand: Nr.: 7 Name:   S ( 0.254, 0.747, 0.254) Abstand : 2.33272
3. Ligand: Nr.: 5 Name:   S ( 0.754, 0.747,-0.246) Abstand : 2.33412
4. Ligand: Nr.: 6 Name:   S ( 0.754, 0.247, 0.254) Abstand : 2.37553
Winkel: (1.Ligand,Zentralatom,2.Ligand) = (  S, Zn,  S)    : 110.320
Winkel: (1.Ligand,Zentralatom,3.Ligand) = (  S, Zn,  S)    : 110.294
Winkel: (2.Ligand,Zentralatom,3.Ligand) = (  S, Zn,  S)    : 110.151
Winkel: (1.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.761
Winkel: (2.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.655
Winkel: (3.Ligand,Zentralatom,4.Ligand) = (  S, Zn,  S)    : 108.611
**********************************************************************
Analyse der Koordinationspolyeder Teilchen Nr. 5 (S):
1. Ligand: Nr.: 1 Name:  Zn ( 1.000, 1.000, 1.000) Abstand : 2.32903
2. Ligand: Nr.: 3 Name:  Zn ( 1.000, 0.500, 0.500) Abstand : 2.33225
3. Ligand: Nr.: 4 Name:  Zn ( 0.500, 0.500, 1.000) Abstand : 2.33412
4. Ligand: Nr.: 2 Name:  Zn ( 0.500, 1.000, 0.500) Abstand : 2.37634
Winkel: (1.Ligand,Zentralatom,2.Ligand) = ( Zn,  S, Zn)    : 110.285
Winkel: (1.Ligand,Zentralatom,3.Ligand) = ( Zn,  S, Zn)    : 110.266
Winkel: (2.Ligand,Zentralatom,3.Ligand) = ( Zn,  S, Zn)    : 110.120
Winkel: (1.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.790
Winkel: (2.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.703
Winkel: (3.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.630
**********************************************************************
Analyse der Koordinationspolyeder Teilchen Nr. 6 (S):
1. Ligand: Nr.: 3 Name:  Zn ( 1.000, 0.500, 0.500) Abstand : 2.32890
2. Ligand: Nr.: 1 Name:  Zn ( 1.000, 0.000, 0.000) Abstand : 2.33349
3. Ligand: Nr.: 2 Name:  Zn ( 0.500, 0.000, 0.500) Abstand : 2.33361
4. Ligand: Nr.: 4 Name:  Zn ( 0.500, 0.500, 0.000) Abstand : 2.37553
Winkel: (1.Ligand,Zentralatom,2.Ligand) = ( Zn,  S, Zn)    : 110.290
Winkel: (1.Ligand,Zentralatom,3.Ligand) = ( Zn,  S, Zn)    : 110.250
Winkel: (2.Ligand,Zentralatom,3.Ligand) = ( Zn,  S, Zn)    : 110.103
Winkel: (1.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.863
Winkel: (2.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.668
Winkel: (3.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.622
**********************************************************************
Analyse der Koordinationspolyeder Teilchen Nr. 7 (S):
1. Ligand: Nr.: 2 Name:  Zn ( 0.500, 1.000, 0.500) Abstand : 2.32939
2. Ligand: Nr.: 4 Name:  Zn ( 0.500, 0.500, 0.000) Abstand : 2.33272
3. Ligand: Nr.: 3 Name:  Zn ( 0.000, 0.500, 0.500) Abstand : 2.33347
4. Ligand: Nr.: 1 Name:  Zn ( 0.000, 1.000, 0.000) Abstand : 2.37589
Winkel: (1.Ligand,Zentralatom,2.Ligand) = ( Zn,  S, Zn)    : 110.279
Winkel: (1.Ligand,Zentralatom,3.Ligand) = ( Zn,  S, Zn)    : 110.265
Winkel: (2.Ligand,Zentralatom,3.Ligand) = ( Zn,  S, Zn)    : 110.142
Winkel: (1.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.764
Winkel: (2.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.651
Winkel: (3.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.694
**********************************************************************
Analyse der Koordinationspolyeder Teilchen Nr. 8 (S):
1. Ligand: Nr.: 4 Name:  Zn ( 0.500, 0.500, 1.000) Abstand : 2.32854
2. Ligand: Nr.: 2 Name:  Zn ( 0.500, 0.000, 0.500) Abstand : 2.33189
3. Ligand: Nr.: 1 Name:  Zn ( 0.000, 0.000, 1.000) Abstand : 2.33205
4. Ligand: Nr.: 3 Name:  Zn ( 0.000, 0.500, 0.500) Abstand : 2.37696
Winkel: (1.Ligand,Zentralatom,2.Ligand) = ( Zn,  S, Zn)    : 110.336
Winkel: (1.Ligand,Zentralatom,3.Ligand) = ( Zn,  S, Zn)    : 110.306
Winkel: (2.Ligand,Zentralatom,3.Ligand) = ( Zn,  S, Zn)    : 110.213
Winkel: (1.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.715
Winkel: (2.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.621
Winkel: (3.Ligand,Zentralatom,4.Ligand) = ( Zn,  S, Zn)    : 108.599

Literaturverzeichnis

[1] Aarts & van Laarhoven – Statistical Cooling: A General Approach to Combinatorial Optimization Problems. Phillips Journal of Research Nr. 40, 1985. Seite 193 – 226.

[2] Aarts & Korst – Simulated Annealing and Bolzmann Machines. John Wiley & Sons Ltd., 1989.

[3] Allem, Tildesley – Computer Simulation of Liquids. Oxfort Science Publications, 1987. Seite 28 – 31.

[4] Bartholomew-Biggs Recursive Quadratic Programming Based On Penalty Functions For Constrained Minimisation in Nonlinear Optimisation, Theory and Algorithms edited by Dixon, Spedicato, Szegö. Academic Press, 1980. Seite 295 ff.

[5] Bertaud The Electrostatic Energy of Ionic Lattices. J. Phys. Radium, 37, 1952. Seite 499 ff.

[6] Borchart-Ott Kristallographie. Springer Lehrbuch, 3. Auflage Herbst 1990.

[7] Brown, Altermatt – Bond-Valence Parameters Obtained from a Systematic Analysis of the Inorganic Structure Database. Acta Cryst. (1985). Seite 244 – 247.

[8] Burnham – Mineral Structure Energetics and Modelling Using the Ionic Approach. American Mineralogist. Seite 347 – 388.

[9] Busing – WMIN, A Computer Program to Model Molecules and Crystals in Terms of Potential Energy Functions. OAK Ridge National Labatory, Tennessy, 1981.

[10] Busing, Martin, Levy ORFFE, a Fortran Crystallographic Function and Error Program. U.S. Natature Technical Information Service. ORNL-TM-306.

[11] Corona, Marchesi, Martini, Ridella – Minimizing multimodal functions of continous variables with the ”Simulated Annealing“ Algorithm. ACM Trans. on Mathematical Software 13, 262 – 280.

[12] Dove – On the Computer Modelling of Diopside: Toward a Transferable Potential for Silicate Minerals. American Mineralogist, Volume 74, 1989. Seite 774 – 779.

[13] Ewald The Calculation of Optical and Electrostatic Lattice Potentials. Ann. Physik, 19, 1921. Seite 253 ff.

[14] Gilbert – Soft-phere Modell for closed-shell atoms and ions. J. Chem. Phys. 49. Seite 2640 – 2642.

[15] Huang, Romeo, Sangiovanni-Vincentelli – An Efficient General Cooling Shedule for Simulated Annealing. Proc. IEEE International Conference on Computer aided Design, Santa Clara 1986. Seite 381 – 384.

[16] IMSL Users Manual, Fortran Subroutines for Mathematical Applications. IMSL Inc. Volume 3, 1987. Seite 872 – 875.

[17] Johnson, Aragon, McGeoch, Schevon – Optimization by Simulated Annealing: An Experimental Evaluation. List of Abstract, Workshop on Statistical Physics in Engineering an Biology, Yorktown Hights, April 1984, revidierte Version 1986.

[18] Khachaturyan, Semenovskaya, Vainstein – A Statistical-Thermodynamic Approach to Determination of Structure Amplitude Phases. Kristallografiya 24, 1979. Seite 905 – 916. Sov. Phys. Crystallogr. 24, 1979. Seite 519 – 524.

[19] Kirkpatrick, Gelatt Jr., Vecchi – Optimization by Simulated Annealing Chains. IBM Reasearch Report RC 9355, 1982.

[20] Kittel – Einführung in die Festkörperphysik. Verlag Oldenbourg, 8. Auflage.

[21] Krabs – Einführung in die lineare und nichtlineare Optimierung für Ingeneure. Verlag B. G. Teubner, Stuttgart 1983.

[22] Krimphove – Simulated Annealing im Vergleich mit stochastischen Optimierungsstrategien. Diplomarbeit, Mai 1991.

[23] Kroll, Maurer, Stöckelmann, Becker, Fulst, Krüsemann, Stutenbäumer – Simulation of Crystal Structures by a Combined Distance Least Sqares and Valence Rule Method. Programmpaket. Institut für Mineralogie, WWU Münster.

[24] Lundy, Mees – Convergence of an Annealing Algorithm. Mathematic Programming 34, 1986. Seite 111 – 124.

[25] Maradudin, Montrall, Weiss, Ipatowa – Theory of Lattice Dynamics in the Harmonic Approximation. Solid State Physics, Suppl. 3, New York, 1971.

[26] Maurer – Vorlesung zur nichtlinearen Optimierung / Optimierung II. Vorlesung WS 1987, WWU Münster.

[27] Megaw – Crystal Structures: A Working Approach. Verlag Saunders, 1973.

[28] The NAG Fortran Library Introductory Guide – Mark 13. NAG Ltd., Wilkinson House – Jordan Hill Road, Oxford (UK). Chapter E04.

[29] Nijboer, DeWette – On the Calculation of Lattice Sums. Physica 23, 1957. Seite 309.

[30] O’Keeffe, Navrotzky – Structure and Bonding in Crystals, Vol. 2. Academic Press, 1981. Seite 19 – 28.

[31] Pannetier, Bassas-Alsina, Rodriguez-Carjaval, Calgnert – Prediction of Crystal Structures from Crystal Chemistry Rules by Simulated Annealing. Nature Vol. 346, July 26^{th} 1990.

[32] Pauling. The Principles Determining then Structure of Complex Ionic Crystals. Journal of the American Chemistry Society, 51, 1926. Seite 1010 – 1026.

[33] Radzeck – Optimierung durch Simulated Annealing. Diplomarbeit, Dezember 1990.

[34] Romeo, Sangiovanni-Vincentelli – Probabilistic Hill Climbing Algorithm: Properties and Applications. Proc. Chapel Hill Conference on VLSI, Mai 1985. Seite 393 – 417.

[35] Schittkowski – Computational Mathematical Programming NATO ASI Series F. Computer and System Sciences, Vol. 15, Springer Verlag, 1985. Seite 171 – 178.

[36] Törn & Zilinskas – Global Optimization. Springer Verlag, 1988.

[37] Williams Accelerated Convergence of Crystal-Lattice Potential Sums. Acta Cryst, A27, 1972. Seite 452 ff.

[38] White – Conceps of Scale in Simulated Annealing. IEEE International Conference on Computer Design, Port Chester, November 1984. Seite 646 – 651.

Hiermit erkläre ich, daß ich die vorgelegte Diplomarbeit selbstständig angefertigt und keine anderen als die im Literaturverzeichnis angegebenen Quellen benutzt habe.

Münster, im Juli 1994