Document [original]

Computational studies of hybrid

interface formation

Jan M. Knaup

Title image: APTES coadsorption by Jan M. Knaup

Computational studies of hybrid interface

formation

Dissertation zur Erlangung des akademischen Grades

Doktor der Naturwissenschaften (Dr. rer. nat.)

vorgelegt dem

Department Physik der Fakult¨

at f¨

ur Naturwissenschaften

Universit¨

at Paderborn

Dipl. Phys. Jan M. Knaup

Paderborn 2008

Abstract

Hybrid interfaces are interfaces between different classes of materials; among these, the hybrid interface

between native metal oxides and polymers is of great interest for the technologies of adhesive joining,

polymer coatings and metal/polymer hybrid materials. The properties of these metal oxide/polymer

interfaces are performance critical. A deep understanding of the interface properties and the chemical

processes from which they emerge is prerequisite for technological advancement.

As experimental observation of processes at buried interfaces at atomic resolution is highly challenging,

computer simulations are a valuable tool to understand them. However, preexisting computational meth-

ods that are able to describe chemical reactions at solid surfaces with the necessary accuracy, can barely

handle the system sizes necessary for a realistic modeling of the interface. Techniques to couple atomically

resolved simulation methods of different accuracy in order to limit the computational cost of simulations

are widely employed in protein and polymer chemistry. In these approaches, quantum mechanical elec-

tronic structure (QM) methods are only applied to those parts of the system, where their accuracy is

needed, while treating the rest of the system with less complicated molecular-mechanical (MM) methods.

Due to the higher density of bonds, these coupling schemes cannot be directly applied to solids, especially

polar ones. Using the density-functional based tight-binding method as the quantum mechanical part,

the existing QM/MM coupling schemes have been adapted to polar solids. Suitable methods to treat the

QM/MM boundary in directly coupled simulations have been transferred from biophysical applications.

New methods to treat cluster charge artifacts, resulting from selecting a QM zone in a polar material,

were developed and tested on different classes of solids. These methods are applicable to various oxide

surfaces, including metal oxides. The QM/MM approach allows considerable savings in computer time

and thus enables simulations of more complex systems, e.g. taking into account effects of the environment.

Aluminum is gaining importance in automotive industry and is the main structural material in aerospace

applications. Under ambient conditions, aluminum is covered by a native oxide layer. For corrosion

protection as well as for the construction of lighter parts, combining metal and fiber reinforced polymers,

the adhesion of polymers on this native oxide layer is extremely important. The formation of the hybrid

interface between native alumina and a model epoxide adhesive system was simulated. Fully quantum me-

chanical simulations were used to calculate the thermodynamics and barrier energies of the chemisorption

of the model adhesive’s components. The results are in good agreement with experimental observations

and offer a good explanation for the effect of silane adhesion promoters: their chemisorption is favorable

at more surface sites than that of the other components, leading to an increase in the density of covalent

bonds across the interface. The new possibility of QM/MM simulations of this hybrid interface allows to

include the effect of the environment on single adsorption reactions. Initial results indicate that water at

the surface, which is present natively and formed during the condensation of adhesive components, has a

crucial influence on the interface formation.

Kurzzusammenfassung

Hybridgrenzfl¨

achen sind Grenzfl¨

achen zwischen Materialien verschiedener Klassen. Unter diesen ist die

Hybridgrenzfl¨

ache zwischen nat¨

urlichen Metalloxyden und Polymeren von großem Interesse f¨

ur die Tech-

nologien von Klebeverbindungen, Polymerbeschichtungen und Metall/Polymer Hybridmaterialien. Die

Eigenschaften dieser Metall/Polymer Grenzfl¨

achen ist entscheidend f¨

ur die Leistung entsprechender Bau-

teile. Ein tiefes Verst¨

andnis der Grenzfl¨

acheneigenschaften sowie der chemischen Prozesse aus denen diese

hervorgehen, ist die Voraussetzung f¨

ur technologischen Fortschritt auf diesem Gebiet.

Da die experimentelle Beobachtung von Prozessen an verdeckten Grenzfl¨

achen mit atomarer Aufl¨

osung

sehr schwierig ist, stellen Computersimulationen ein wertvolles Mittel f¨

ur deren Verst¨

andnis dar. Aller-

dings sind existierende Simulationsmethoden, die chemische Reaktionen auf Festk¨

orperoberfl¨

achen mit

hinreichender Genauigkeit beschreiben k¨

onnen, kaum in der Lage, f¨

ur die realistische Modellierung dieser

Grenzfl¨

achen hinreichend großer Systeme zu beschreiben. Techniken, atomar aufgel¨

oste Simulationsme-

thoden verschiedener Genauigkeiten direkt miteinander zu koppeln sind in den Bereichen der Protein- und

Polymerchemie weit verbreitet. Bei diesen Ans¨

atzen werden quantenmechanische Elektronenstrukturme-

thoden (QM) nur auf diejenigen Teile des Systems angewandt, bei denen diese Genauigkeit erforderlich ist,

w¨

ahrend der Rest des Systems mit einfacheren molekularmechanischen (MM) Methoden simuliert wird.

Wegen der h¨

oheren Dichte an Bindungen, k¨

onnen diese Kopplungsschemata nicht direkt auf Festk¨

orper-

systeme angewandt werden, insbesondere nicht auf Polare Festk¨

orper. Unter Verwendung von dichtefunk-

tionalbasiertem Tight-Binding als QM Methode, wurden die existierenden QM/MM Kopplungsschemata,

von biophysikalischen Anwendungen kommend, auf Festk¨

orpersysteme transferiert. Geeignete Methoden

zur Behandlung der QM/MM Grenze wurden transferiert. Gleichzeitig wurden neue, geeignete Metho-

den zur Behandlung von Ladungsartefakten, die durch die Auswahl einer QM Zone aus einem polaren

Festk¨

orper entstehen, entwickelt und an verschiedenen Klassen von Festk¨

orpern getestet. Diese Metho-

den k¨

onnen auf verschiedene Oxydoberfl¨

achen, einschließlich Metalloxyden, angewandt werden. Dieser

QM/MM Ansatz erlaubt erhebliche Einsparungen an Rechenzeit, wodurch die Betrachtung komplizierter

Systeme, z.B. unter Ber¨

ucksichtigung von Umgebungseinfl¨

ussen, erm¨

oglicht wird.

Aluminium gewinnt an Bedeutung f¨

ur die Automobilindustrie und ist nach wie vor das vorherrschende

Strukturmatieral in der Luft- und Raumfahrt. Unter normalen Umgebungsbedingungen, ist Alumini-

um stets mit einer nat¨

urlichen Oxydschicht bedeckt. Zum Korrosionsschutz wie auch f¨

ur die Herstellung

leichterer Bauelemente, die Metall- und faserverst¨

arkte Kunststoffteile vereinen, ist die Haftung von Poly-

meren auf dieser nat¨

urlichen Oxydschicht von gr¨

oßter Bedeutung. Die Ausbildung der Hybridgrenzfl¨

ache

zwischen ist daf¨

ur von besonderer Bedeutung. Sie wurde f¨

ur die Grenzfl¨

ache zwischen nativem Aluminiu-

moxyd und einem Modellsystem f¨

ur Epoxydklebstoffe simuliert. Mittels voll quantenmechanischer Simula-

tionen wurden die Thermodynamik und die Energiebarrieren der Chemisorptionsreaktionen der einzelnen

Klebstoffkomponenten berechnet. Die Ergebnisse sind in guter ¨

Ubereinstimmung mit experimentellen Be-

obachtungen und bieten eine gute Erkl¨

arung f¨

ur den Wirkmechanismus von Silan-Haftvermittlern: deren

Chemisorption ist an mehr Stellen der Oberfl¨

ache energetisch g¨

unstig, als die anderer Komponenten. Auf

diese Weise f¨

uhrt der Haftvermittler zu einer Erh¨

ohung der Dichte an kovalenten Bindungen ¨

uber die

Substrat/Polymer Grenzfl¨

ache. Die neue M¨

oglichkeit, QM/MM Simulationen dieser Hybridgrenzfl¨

ache

durchzuf¨

uhren, erlaubt die Ber¨

ucksichtigung von Umgebungseinfl¨

ussen auf die einzelnen Adsorptionsre-

aktionen. Erste Ergebnisse zeigen, dass Wasser an der Oberfl¨

ache, welches sowohl von Natur aus dort

vorhanden ist, als auch im Laufe der Kondensation von Klebstoffkomponenten dort gebildet wird, ent-

scheidenden Einfluss auf die Grenzfl¨

achenausbildung aus¨

ubt.

Contents

1 Introduction 1

2 QM Treatment 3

2.1 Notation and units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.2 Condensed matter as a many body problem . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.3 The Born-Oppenheimer approximation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.4 Modeling the Solid State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.4.1 The perfect crystal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.4.2 Defective crystals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.5 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.5.1 Exact DFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.5.2 The local density approximation (LDA) . . . . . . . . . . . . . . . . . . . . . . . . 10

2.5.3 Treating electron spin (LSDA) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.5.4 Beyond LDA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.5.5 The DFT gap error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.5.6 Hybrid exchange functionals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.6 The pseudopotential method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.7 Density Functional based Tight Binding . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.7.1 Tight Binding schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.7.2 DFTB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.7.3 Extensions to SCC-DFTB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.8 Methods of QM/MM coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.8.1 The QM/MM boundary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3 Modeling Chemical Reactions 23

3.1 Reaction energetics at 0 K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.1.1 Finding barrier energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.1.2 Nudged elastic band . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.2 Beyond the minimum energy picture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.3 Modeling rare events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

vii

viii CONTENTS

4 DFTB QM/MM embedding 31

4.1 Challenges of QM/MM embedding in solids . . . . . . . . . . . . . . . . . . . . . . . . . . 32

4.2 Embedding approaches for solid state materials . . . . . . . . . . . . . . . . . . . . . . . . 33

4.2.1 Representation of the external field . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.2.2 Treatment of the dangling QM–MM bonds . . . . . . . . . . . . . . . . . . . . . . 34

4.2.3 Cluster neutralization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

4.3 Validation of the QM/MM embedding in different materials . . . . . . . . . . . . . . . . . 37

4.3.1 Evaluation Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

4.3.2 The fully QM treated surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

4.3.3 Evaluation of coupling schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

5 Epoxy adhesives on native Al2O359

5.1 The native Al2O3surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

5.2 A model adhesive system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

5.3 Surface model preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

5.3.1 Bulk γ-Al2O3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

5.3.2 The hydroxylated γ-Al2O3surface . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5.4 The H2O-assisted ring-opening of DGEBA . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5.5 Adsorption of single adhesive component molecules in vacuum . . . . . . . . . . . . . . . . 64

5.5.1 DETA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.5.2 t-DGEBA(open) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.5.3 AMEO(hyd). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.5.4 Comparison of adsorption mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . 66

5.6 The effect of neighboring adsorbates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

5.6.1 Adsorption energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

5.6.2 Adsorption paths . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

5.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

5.8 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

6 Conclusion 81

Bibliography 83

A QM/MM Evaluation Data 93

A.1 Simple Link Atoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

A.2 Neutralized Clusters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

Acknowledgments 117

CONTENTS ix

Colophon 119

List of Figures 121

List of Tables 125

xCONTENTS

Chapter 1

Introduction

Two of the most grave challenges faced by humankind today are the global warming and the expiration of

fossil fuel supplies. Both are closely related, since a large part of the greenhouse effect is caused by CO2

emissions from the burning of these fossil fuels [1]. As fossil fuels are currently the most important source

of energy for human societies, this will have deep ramifications on the future lifestyle of all humankind,

if no solutions are found. One necessary part of any solution will be energy conservation measures: in

the short term to slow down CO2emissions, in the middle run to facilitate the transition to sustainable

energy sources and on the long term to allow the worldwide standard of living to increase sustainably

to the level of industrialized countries. At the same time, the rising prices of fossil fuels, especially oil,

create a strong economic incentive to reduce consumption of these commodities.

One way to conserve energy is to reduce the use of energy-consuming technologies (e.g. all kinds of

electronics). This is in strong conflict with the natural human desire to improve the individual standard

of living and the humanitarian aim to reduce poverty. Therefore, fossil fuel conservation by abstinence

from technology can only be applied to a limited extent without raising severe social problems. Another

possibility is to increase the energy efficiency of current technologies. In this work, I will concentrate on

one aspect: the development of more capable adhesive joining technologies, which, among others, enable

lighter and thus more energy efficient transportation.

Adhesive joints work from the principle of inserting a layer of a third material between the surfaces of

the parts which are to be joined. This intermediate layer must at the same time form a strong bond

to the surfaces of both joined parts and have sufficient cohesion to be able to transmit the mechanical

forces between them. Additionally, the intermediate layer fills the roughness of the joined part’s surfaces,

thus filling the gap between the parts and sealing the joint against gases or liquids. Very closely related

to adhesives are coatings applied in liquid form, i.e. paints. In principle, a paint is nothing more than

a mixture of an adhesive and some type of dye or pigment. In both cases, the adhesion between the

polymer phase and the substrate is of paramount importance. In my work, I concentrate on the adhesion

at the polymer/substrate interface.

In most cases, the adhesive is applied to the surfaces in liquid form and then cured to form a solid layer.

Curing can be achieved in several ways: e.g. by evaporation of a solvent, by means of a phase change upon

cooling, or by polymerization of an adhesive mixture. For automotive and aerospace applications, the

latter class of adhesives are of great interest, because of their high mechanical and chemical endurance.

Therefore I will concentrate on this type of adhesives, especially on a model epoxide resin adhesive. Its

chief characteristics are being a thermoset type polymers, i.e. the polymerization is thermally activated,

and the fact that a strongly cross-linked network is formed from relatively low molecular weight precursors.

Although automotive industry has moved away from the formerly widely employed epoxy paints, for

reasons of work safety and expensive waste disposal, these paints still play a very important role in

aerospace applications. Results obtained on the interface properties between a substrate surface and an

epoxy polymer can be equally applied to joining and coating technologies.

2CHAPTER 1. INTRODUCTION

I focus my interest on the interface between the natively oxidized aluminum surface and a model epoxide

adhesive polymer, especially on the chemical processes from which the properties of this interface arise.

This choice of materials is governed by the tradition of using aluminum as a lightweight structural material

in aerospace applications and the trend in automotive industry to (partly) replace steel by aluminum to

conserve weight. A good understanding of the fundamental mechanisms is critical for being able to

improve any technology in a systematic and efficient manner.

When aiming to understand these underlying mechanisms of a chemical process, it is necessary to analyze

the behavior of the atoms involved. In many cases experimental observation is very difficult. The direct

observation of individual atoms during a reaction is only possible with very large effort and only for a few

specific reactions at solid surfaces in (nearly) vacuum. Here, the theoretical modeling of the reaction of

interest in a bottom-up approach, i.e. simulating the behavior of the individual atoms from first principles,

can give extremely valuable insights. However, since modeling reactions at hybrid interfaces requires the

simulation of comparatively large numbers of atoms1, usual ab initio methods for atomistic simulation

are too computationally involved. At the same time, simpler methods lack transferability, the ability to

describe bond formation or breaking (or both) and they are less accurate. Methods that model parts of

an examined system at different resolutions, e.g. taking into account the electronic structure in one part,

while using empirical interatomic potentials in another, depending on the accuracy required for each part,

have been successfully applied in biophysics and polymer science. However, these approaches are cannot

be readily transferred to solids.

In this work, I present my results on porting the approaches to couple the approximate density functional

theory method DFTB to classical force-fields from biophysical applications to solid surfaces. I apply both

DFTB and DFTB/MM to examine the formation of the interface between a model epoxide adhesive and

the native oxide surface on aluminum. I shall demonstrate that DFTB coupled to classical force fields is

capable of providing deep insight into the formation processes of hybrid interfaces.

1Compared to what is generally referred to as “quantum chemistry”, where mostly systems of <100 atoms are regarded.

Chapter 2

Quantum Mechanical Treatment of

Molecules and Solids

The description of the electronic properties of solids and molecules requires a quantum mechanical treat-

ment. This section shall describe the quantum mechanical methods employed in this work.

My aim is to simulate solid surfaces and interfaces at atomic resolution in a computationally efficient

yet accurate manner, in order to model chemical reactions at hybrid interfaces. I will employ a coupling

scheme between density-functional based tight-binding (DFTB) and molecular mechanics (MM) force-

fields, concentrating on the quantum mechanics side of the coupling. As a foundation for describing the

DFTB method and DFTB/MM coupling in solids, I shall first give an introduction into the quantum

mechanical description of solids in general and into density functional theory, to which DFTB is an

approximation.

In this chapter, I shall give a concise description of the theoretical foundations of my research presented

in this thesis. For an in-depth introduction, the reader should consider one of the excellent textbooks

available on this matter, e.g. Refs [2–5].

2.1 Notation and units

Since in this work I treat both molecules and solids, I shall use the term system unless a statement

explicitly refers to either an infinite crystal, or a finite molecule or section of a crystal.

The number of electrons in a system is Nand indices irun over the electrons. Likewise, Mis the number

of atoms, indexed by j. Indices k, l refer to electrons or nuclei, depending on whether they are used in

conjunction with ior j.

Vectors rwithout an index contain the spacial coordinates of the whole electronic system, while rigive

the coordinates of a single (quasi-) particle. In contrast, vectors ~r are three-dimensional vectors in real

space. Similarly vectors Rand Rirefer to the nuclei.

Φ(r,R, t) and Ψ(r, t) are the wave functions of the whole system and the electronic subsystem, respec-

tively. If present, αandβrefer to the two spin components of the electronic system. (Note that the

independent variable ~r is omitted by text-book convention, except for some special occasions.) If the

time variable tis absent, wave functions give the stationary ground-state. Similarly, n(~r) gives the total

electronic charge density of the system, while ni(~r) is the density contribution of a single quasi-particle.

In the following, equations will be given in atomic units[6,7], which are defined by assigning unity values

to the electron rest mass m, elementary charge e,1

4πǫ0and reduced Planck constant ~. From this definition

4CHAPTER 2. QM TREATMENT

follows, that the Bohr radius rB=~2

me2≈0.5291 ˚

A and Hartree energy EH=me4

~2≈27.2114 eV assume

unity values as well. Since this makes most constant prefactors vanish from the equations, these units are

very often applied in quantum mechanics. Additionally, atomic units also reduce numerical complexity.

The units of length and energy are also referred to as Bohr and Hartree. Using them, new dimensionless

variables rand Eare defined:

r=˜

, E =˜

.(2.1)

Here, ˜rand ˜

Edenote the corresponding variables with dimension. Using atomic units, e.g. the Schr¨

odinger

equation of a single electron in the potential of a nucleus with the charge Z, changes from

−~2

2me

∆−Ze2

4πε0|˜

r−˜

R|Ψi(˜

r) = ˜

EiΨi(˜

into the dimensionless form −1

2∆−Z

|r−R|Ψi(r) = EiΨi(r).

2.2 Condensed matter as a many body problem

The theory behind the description of molecules and solids is in many aspects the same. In practically all

aspects, the description of the solid state contains all the necessary elements for the modeling of molecules

plus some additional features introduced by the symmetries of the solid state, covered in section 2.4.

In order to examine different geometrical or chemical configurations of a solid or molecule, one needs to

know its total energy depending on the geometry as well as the derivatives of the total energy in the

systems coordinates.

The quantum mechanical state ground state of a system of N electrons and M nuclei is described by the

time independent Schr¨

odinger equation

HΦ = EΦ (2.2)

with the (non-relativistic) Hamiltonian

H=−

j=1

2Mj

∆Rj

|{z }

j6=j′

ZjZj′

|Rj−Rj′|

|{z }

−1

l=1

∆rl

|{z }

Tel

l6=l′

|rl−rl′|

|{z }

Cel

−

M,N

j,l

|Rj−rl|

|{z }

Cel,N

.(2.3)

The corresponding wave function Φ describes the Mnuclei and Nelectrons. Since no general analytic

solution to this problem is known, approximations must be made, which will be described below.

2.3 The Born-Oppenheimer approximation.

The first approximation generally made in the simulation of solids and molecules is the Born-Oppenheimer

approximation, which is based on the mass ratio of nuclei versus electrons of about 103[8]. As the forces

acting on nuclei and electrons can be regarded as within the same order, the momenta of both will also

be similar. Comparing TNand Tel in equation 2.3 shows that the nuclear kinetic energy will then be

only about 10−3of kinetic energy of the electrons. This allows to separate the nuclear and electronic

wave functions according to the adiabatic theorem[9]. Therefore, the motion of electrons and atoms may

be treated separately, with the atomic coordinates at a given time as a (constant) input parameter to

the electronic equation and the time average of the electronic configuration as a constant in the atomic

2.4. MODELING THE SOLID STATE 5

equations of motion. This separation is exact if the potential acting on the nuclei is strictly harmonic in

the vicinity of its local minima1.

In the context of chemical reactions, the system of nuclei is usually treated classically, i.e. by solving

the Newton’s equations of motion for the nuclei. The electron problem remains to be treated quantum

mechanically with the electron Schr¨

odinger equation

Hel (Rj) Ψ = EelΨ (2.4)

and its Hamiltonian

Hel =−

i=1

∆ri+1

i6=i′

|ri−ri′|−

M,N

j,i

|Rj−ri|.(2.5)

Here Ψ is the electronic wave function. In the following, ˆ

Hwill be used as a shorthand to denote the

electronic Hamiltonian, where not specified otherwise. Still, even the electron many-body-wave function

cannot be solved analytically. This makes further approximations necessary to be able to solve the

problem numerically. Before describing one possible approaches to achieve this in section 2.5, it is

necessary to discuss some general concepts for the description of solids and the electronic wave functions.

2.4 Modeling the Solid State

When dealing with the defective solid, one cannot handle every atom of a system of technologically

relevant size explicitly. Therefore a model of adequate size is required. In the case of the perfect crystal,

the translational symmetry can be utilized, but that is broken by the introduction of a defect. In the

following section, starting from the treatment of the perfect crystal, I shall describe possibilities to handle

the defective solid.

2.4.1 The perfect crystal

The perfect crystal is by definition invariant under translation by the lattice vectors

Rl=

i=1

liai(2.6)

with the primitive unit cell vectors aiand arbitrary integers li. As a consequence, a reciprocal lattice in

momentum space can be constructed with the reciprocal lattice vectors

Gg=

j=1

gjbj(2.7)

with arbitrary integers gjand the primitive reciprocal lattice vectors bj

ai·bj= 2πδij.(2.8)

The symmetric unit cell of the reciprocal lattice is defined as the set of points closer to one lattice point

than to all others and is called the Brillouin zone (BZ), i.e. the Wigner-Seitz cell of the reciprocal lattice.

1Which, is often not the case, especially for point defects in solids. Still the harmonic approximation must be made,

because the full Schr¨

odinger equation is computationally too involved to solve for systems large enough to allow the modeling

of the relaxations around such point defects.

6CHAPTER 2. QM TREATMENT

Bloch’s theorem states that, due to the translational symmetry, the electronic wave functions of the solid

(called Bloch-electrons) can be expressed by the Bloch-wave functions

ϕnk(~r) = unk(~r)eik·~r.(2.9)

With the wave-vector krestricted to the BZ and the function unk(~r) being periodic on the lattice, the

Bloch-wave functions satisfy Bloch’s periodicity condition

ϕnk(~r +Rl) = ϕnk(~r)·e(ikRl)(2.10)

Since the translational invariance along arbitrary lattice vectors implies an infinite crystal, while the

real solid is finite (Niunit cells along the lattice vector ai,) the perfect translational symmetry must be

ensured by imposing artificial macroscopic periodicity, i.e. applying the cyclic boundary conditions of

Born and von K´arm´an:

ϕnk(~r +Niai) = ϕnk(~r) (2.11)

which makes the number of k-vectors in the BZ equal to N=N1N2N3and their density equal to V/8π3

with the volume of the solid V.

Due to this construction, an identical set of k-dependent equations applies for all unit cells. In principle,

these equations should be solved for all k-points in the BZ, since they contain an average over all k.

Instead, the so-called special k-point theorem [10] is often applied to reduce the number of equations to

be solved. It states, that the average of the function f(k)

f=1

VBZ ZVBZ

f(k)d3k(2.12)

can be approximated by a weighted sum over special k-points

f=X

ωqf(kq).(2.13)

If f(k) is invariant under translation by the reciprocal lattice vectors Ggand to the symmetry operation

ˆαof the point group Gof the lattice, then the Fourier expansion of f(k) can be split at a term M:

f(k) = f0+

m=1

fmAm(k) +

∞

m=M+1

fmAm(k) (2.14)

with

Am(k) = 1

NGX

ˆα∈G

eik(ˆ

αRm),(2.15)

where NGis the number of elements in the point group. Since for any m, the BZ-integral of Am(k)

vanishes, fequals f0. Substituting equation 2.14 into equation 2.13 then yields

f=X

ωqf0+

m=1

fmAm(k) +

∞

m=M+1

fmAm(k).(2.16)

For sufficiently high M, the Fourier coefficients fm, m > M become negligible, so that the last term can

be neglected. Provided, that the sum of the weighting factors ωqis unity, the error of the approximation

is determined by

f−f≈

m=1

fmX

ωqAm(kq).(2.17)

Therefore, the special k-point set and the corresponding weighting factors should be chosen so that

PqωqAm(kq) vanishes for m < M. The error of this approximation diminishes almost exponentially

with M.

2.4. MODELING THE SOLID STATE 7

Generally, the procedure proposed by Monkhorst and Pack [11] is used to generate special k-point sets

with Mas high as possible. For a given number Q, the special k-point set is generated from the primitive

reciprocal lattice vectors with the coefficients

qi=2pi−Q−1

2Q+q0,ipi, i = 1,2,3; pi= 1 . . .Q (2.18)

with the weighting factors ωq= 1/Q3and global shifts q0,i ∈[0..1], generally chosen as either 0 or

0.5 to specifically include or exclude k=~

0. The set obtained from this procedure is called a QxQxQ

MP set. For lattices with lower than cubic symmetry, Q may be different in the three directions of the

reciprocal lattice vectors. Under the symmetry operations of the group, some of the resulting k-points

may be equivalent, in such cases, only one k-point has to be evaluated with the respective weighting

factor multiplied by the degeneracy of the k-point. In this case, non-scalar properties are automatically

symmetrized.

2.4.2 Defective crystals

A single point defect introduced into a crystal destroys the translational symmetry completely. Thus, in

principle, equations for all atoms in the solid have to be solved. One could, in principle, treat a smaller

piece but in order to model the behavior of a defect in an otherwise undisturbed crystal, the model has

to be large enough, that a sufficient amount of atoms remain undisturbed by the defect to preserve the

crystalline properties of the solid. This means that at the edges of the model,

1. the amplitude of the localized wave functions,

2. the deviations of charge density and

3. the deviations of the host atom positions with respect to the perfect crystal

must vanish at least approximately. Since the perturbations from point defects can have ranges of 10 –

50˚

A, ∼2500 – 25000 atoms would be necessary to achieve this. Such numbers of atoms are orders of

magnitude larger than what can be handled routinely with ab initio methods on computers of the present

and the foreseeable future. Therefore, simplified models become necessary, which basically means, that

the fulfillment of the conditions mentioned above is assumed for a smaller group of atoms around the

defect usually referred to as the cluster. When constraining the calculations to manageable cluster sizes,

these criteria usually cannot be met, therefore some allowance for deviations must be made and some

kind of convergence test for all three of them becomes necessary in order to assess the deviations.

The first approach is to handle only the atoms of the cluster explicitly and embed them into an unper-

turbed crystalline background.

The simplest way to achieve this is the so-called molecular cluster model (MCM), which is based on

the assumption, that the electronic states of some solids are well described by localized bonds. This is

mostly the case with semiconductors. Ideally, the quantum mechanical interaction of the cluster with

the background is completely described by a set of localized bonds crossing the interface between the

two regions. If these bonds could be kept unchanged when cutting off the atoms of the background, the

embedding would still be perfect, except for long range Coulomb interactions. Since the dangling bonds

at the interface can never be saturated in a way that keeps them perfectly unchanged, in most cases

hydrogen atoms are used to saturate them. The hydrogen saturated cluster is a molecule in the chemical

sense, hence the name molecular cluster model. The MCM also forms the basis for quantum mechanics

to molecular mechanics (QM/MM) coupling schemes, which will be described in detail in section 2.8.

More complicated are the perturbed crystal and perturbed cluster approaches. In the case of the per-

turbed crystal approach, the Hamiltonian is separated into the Hamiltonian of the perfect crystal and a

perturbation ˆ

Urelative to that:

H=H+ˆ

U(2.19)

8CHAPTER 2. QM TREATMENT

which can be solved using a Green’s-function approach [12–14]. In the perturbed cluster approach, relying

on a localized basis, the Hamiltonian is partitioned in direct space according to basis functions in the

cluster (C) and the background (B)[15]

H=HCC HCB

HBC HBB .(2.20)

Another possibility is called the supercell model (SCM). In this approach, the translational symmetry

is restored artificially by constructing a periodic superlattice with the cluster as unit cell [16,17]. This

requires that the cluster is shaped in a way that allows periodic repetition without gaps or overlaps.

Clusters of this kind are generally referred to as supercells. In the supercell model, the momentum space

description of the perfect crystal outlined above can be applied directly without modification. One must

consider however, that the supercell usually contains a multitude of primitive unit cells in terms of the

perfect lattice. Therefore, the BZ of the supercell is reduced relative to the primitive Brillouin zone (PBZ)

and is consequently called the reduced BZ (RBZ). In a perfect supercell, i.e. one without any defect, the

vectors Kof the RBZ each represent a set of primitive k-vectors, so that, if the superlattice has the same

point group as the primitive one, the set {k}0corresponding to the center of the RBZ (K= 0) forms

a special k-point set [18]. Therefore, with a properly chosen SCM, a calculation limited to the Γ-point

(K= 0) of the RBZ may give a good description of the extend crystalline states.

If, however, the supercell contains a defect, the Γ-point is definitely no special k-point for the superlattice

of defects and, unless the supercell is sufficiently large, that K= 0 approximation might not be adequate.

The neutral vacancy in Si is a good test case for convergence. It has been shown [19] that the formation

energy computed in the Γ-point approximation and with a 2x2x2 MP set in 32, 64, 128 and 216 atom

supercells differs by 3.05,0.56,0.30 and 0.04eV respectively. This shows the importance of testing the

K-point set used in supercell calculations, especially in relatively small supercells.

In practice, the atoms at the boundary of the supercell are sometimes fixed during calculations. This

assumes, that the third criterion from page 7is met. (For deep levels, the defect wave function is only

considerably influenced by the relaxation in the first 1–2 neighbor shells of the defect.) This artificial

restriction of the host atom relaxation is desired in order to limit the effect of defect-defect interactions

introduced by long range relaxation beyond the SCM boundary.

Criterion 2 is usually met by SCM-s in the size range of 50–100 atoms in metallic or purely covalent solids,

but may require considerably larger cell sized for strongly ionic materials. Due to the long range of the

Coulomb interaction, the defect may cause far reaching polarization effects which have to be accounted

for in a self-consistent manner.

It should be noted, that the periodic boundary conditions can lead to severe problems in the coulomb

energy terms. If the SCM has a total charge, the interaction between this charge and its infinite number

of periodic images will lead to an infinite coulomb energy. Therefore a total charge of the supercell

must be neutralized, which is generally done by applying a uniform compensating charge density over

the whole volume of the supercell: ρcomp =−Qtot

Vcell . However, this approach can lead to severe artifacts

which depend, among others, on the supercell size. Therefore, the convergence of the calculation results

with supercell size should be checked. Similarly, an overall dipole moment of a supercell will also lead to

divergent coulomb energy. This is generally compensated by a shift of the charge density with respect to

the nuclei. This effect leads to an artificial energetic destabilization of defect configurations with dipole

moments. In covalent solids, this is generally not a serious problem, while in polar materials, the error

can become severe. Again, convergence checks with supercell size are necessary to check for this problem,

however, increasing the supercell size is often impractical.

2.5. DENSITY FUNCTIONAL THEORY 9

2.5 Density Functional Theory

2.5.1 Exact DFT

In materials science, the most frequently employed method to tackle the electronic problem is the density

functional theory (DFT). An in-depth treatment of this method can be found in the original articles

by Hohenberg and Kohn [20], Kohn and Sham [21] and the more mathematical treatment by Lieb [22].

DFT is based on the theorem by Hohenberg and Kohn [20] that the non-degenerate ground state of an

electronic system and its ground state energy are unique functionals of its particle density n(r)

E0=E[n(~r)].(2.21)

This allows to replace the wave function Ψ(r1,...,rN) with the particle density n(r) as the central

variable. The particle density2is an observable of the quantum mechanical system given by

n(~r) = hΨ|

i=1

δ(~r −ri)|Ψi.(2.22)

The total energy from equation 2.21 is given by the Kohn-Sham energy functional H[n]:

H[n] = [−T[n(~r)] + veff[n(~r)]] ,(2.23)

where the first term is the kinetic energy T[n(~r)] and veff [n(~r)] is the effective potential, which comprises

the terms

veff[n(~r)] = vext(r) + Z2n(~

r′)

|~r −~

r′|d~

r′+VXC[n(~r)].(2.24)

Here vext describes the external potential (the potential of the nuclei and external electromagnetic field),

the middle term is the Coulomb potential, vc, and VXC[n(r)] is the exchange-correlation potential, which

describes all many-particle effects.

With an expansion of the particle density into effective single particle wave functions Ψi(~r), the so-called

Kohn-Sham orbitals

n(r) =

i=1

|Ψi(~r)|2,(2.25)

the N-electron equation can be converted into a system of Ncoupled effective-single-electron equations,

the Kohn-Sham equations [21]

−1

2∆ + veff[n(~r)] − EiΨi(~r) = 0.(2.26)

Both vcand VXC are functionals of the particle density n(r) and therefore all other Kohn-Sham orbitals

Ψj(r) with j6=i. Thus the effective potential couples all effective single particle equations – the Ψiare

only formally independent.

Like the potential, the total energy can be split up into its known components:

Etot[n(~r)] = T[n(~r)] + ZZ n(~r)n(~

r′)

|~r −~

r′|d~r d~

r′+Zvext(~r)n(~r)d~r +EXC[n(~r)].(2.27)

In equations. (2.24) and (2.27), all terms but the last, follow directly from regarding the system as the

composition of Nindependent, non-interacting particles. Thereby all many-body effects, except the

2The particle density of the electronic system is also referred to as the charge density. Both terms are used synonymously

here, as in atomic units the electron charge, e, is unity.

10 CHAPTER 2. QM TREATMENT

coulomb repulsion, have been combined in the VXC and EXC terms. The exchange-correlation potential

and the exchange-correlation energy are connected by the functional derivative

VXC =∂EXC[n(~r)]

∂n(~r).(2.28)

Up to this point, DFT is a reformulation of the Schr¨

odinger equation without neglect of any part, as long

as the ground state can be expressed as a sum of single-determinant wave functions. However, neither the

exact form of EXC nor that of VXC is known. Therefore, the computational advantage of transforming

the total energy expression from being dependent on a 3Ndimensional function in space to just one, has

been bought with the need to introduce an approximation for EXC or VXC.

2.5.2 The local density approximation (LDA)

The first approximation to the exchange-correlation functional in DFT to be presented here is the local

density approximation. It is based on the assumption that VXC is only locally dependent on n(r). Thus,

it can be expressed in an integral form with an integrand ǫXC [n(r)], that only locally depends on the

density

EXC[n(~r)] ≈Zn(~r)ǫXC[(n(~r)]d~r (2.29)

for the exchange-correlation potential, this leads to:

vXC =∂EXC [n(~r)]

∂n(~r)≈∂

∂n(~r)Zn(~r)ǫXC[n(~r)]d~r (2.30)

=ǫXC[n(~r)] + n(~r)∂ǫXC[n(~r)]

∂n(~r)(2.31)

which then gives the Kohn-Sham equations in this local approximation

−1

2∆ + veff −EiΨi(~r) = 0 (2.32)

with veff(~r) = vc(~r) + vext(~r) + vXC(~r). The equations are still coupled via the effective potential veff, the

Kohn-Sham equations must be solved self-consistently.

It should be noted at this point, that all standard implementations of local exchange-correlation func-

tionals suffer from the problem, that the gap energies of solids are underestimated by 30-50%. This will

be discussed further in section 2.5.6. In some pathological cases, narrow gap semiconductors even appear

to be metallic.

2.5.3 Treating electron spin (LSDA)

For non closed-shell systems, it is necessary to include spin-polarization effects into the calculations. The

easiest way to describe a spin polarized system is, to split the charge density n(r) into two components

for the two possible spin orientations

n(~r) = n↑(~r) + n↓(~r) (2.33)

in this case the DFT becomes spin density functional theory (SDFT). The two components are then

treated in exactly the same way, as the particle density in the case of neglected spin effects, including

any approximations. (I.e. the local density approximation becomes the local spin density approximation

— LSDA.) With this, the one can define the spin- or magnetization density

m(~r) = n↑(~r)−n↓(~r) (2.34)

2.5. DENSITY FUNCTIONAL THEORY 11

and the corresponding spin-polarization

ξ=m(~r)

n(~r)=n↑(~r)−n↓(~r)

n↑(~r) + n↓(~r).(2.35)

This allows to treat spin-polarization effects of the system and perform open-shell calculations.3

2.5.4 Beyond LDA

Although LDA-DFT performs well for most applications, in some cases the purely local exchange-

correlation functionals do not produce sufficiently accurate results. To improve on the local functionals,

without having to apply a fully non-local formulation, which would be computationally far more expen-

sive, the generalized gradient approximation (GGA) was developed. It works by extending the LDA EXC

with a term dependent on the local gradient of the charge density

EGGA

XC [n(~r)] = ELDA

XC [n(~r)] + EXC[∇n(~r)].(2.36)

As they take into account the surrounding density around each point rwithout being fully non-local,

these functionals are also called semi local. For molecular systems, functionals combining the exchange

functional by Becke [23] with a correlation term proposed by Lee, Yang and Parr [24], called BLYP[25]

is very popular, in solid state physics, the exchange-correlation functional by Perdew, Burke and Ernz-

erhof [26] (PBE) is widely used. One must, however, be mindful that that the GGA does not always

improve on LDA, one prominent example being metals.

2.5.5 The DFT gap error

As mentioned earlier, the HF method, using exact exchange, is rather costly (∼N4) even without

considering correlation corrections, the calculation of which scales > N4. On the other hand, the very

effective local density or gradient corrected approximations of DFT fail to describe the excited electronic

states of most systems.

The LDA, commonly used in defect physics, is a very good approximation for systems where the electrons

are rather delocalized [27]. However, when the static electron-electron interaction becomes strong enough

that electrons get localized on atomic sites in a solid, the LDA, as well as GGA, fail to describe the correct

ground states. One of the most prominent examples are the 3dtransition metal oxides like NiO (and

probably also TiO2) which are Mott insulators characterized by localized d-electrons [28]. The LSDA in

contrast predicts somewhat too small magnetic moments and vanishing or at least very small gaps, also

observed in case of semiconductors [29]. Unfortunately, this error also influences the accuracy of calcu-

lating spectroscopic properties of defects therein. Wrong excited states and underestimated transition

levels between the different charge states are the consequence.

An important insight on the central origin of these problems can be obtained by the observation that

localized orbitals in the LDA give rise to an error due to a spurious self-interaction contained in the

effective potential

ELDA

SI [ni] = Zni(~r)ni(~

r′)

|~r −~

r′|d~rd~

r′+Zni(~r)ǫLDA

XC (ni(~r))d~r 6= 0.(2.37)

This error increases the more localized the orbitals become. For rather delocalized electrons, however,

the self-interaction is negligible or zero. In HF theory, where the exact exchange is calculated, the

self-interaction vanishes.

3In order not to overload the formalism, the spin indices will be neglected in the following. The approach for spin

polarized systems is completely analogous.

12 CHAPTER 2. QM TREATMENT

2.5.6 Hybrid exchange functionals

Since the exchange energies are typically an order of magnitude higher than the correlation energies, it

could be shown [30], that using exact exchange within DFT can largely eliminate the gap error. Still,

at present, no implementation of exact exchange DFT exists, which could be used efficiently on large

systems [31].

For the description of exited states a hybrid functional, e.g. the B3-LYP functional [25], is often used.

This is constructed by mixing exact exchange to the DFT exchange.

It has been shown [32] that with the mixing in the form of

EXC =ELSDA

XC +λ(Eexact

X−ELSDA

X),(2.38)

where λdenotes the empirical mixing factor, the gap error can be effectively eliminated.

Previous work[25,33] has shown, that a 20% admixture of exact exchange can reproduce the experimental

band gaps of 3C- (2.4 eV) and 4H-SiC (3.3 eV) as well as the positions of known defect levels very well.

However the empirical mixing factors, fitted to the fundamental gap, are not well transferable. The same

mixing factor only yields a band gap of 8.1 eV for α-quartz, as opposed to the experimental value of 9 eV,

while a mixing weight of 28% proves optimal for this material.

Still, within one material, consistently good results can be obtained for all properties. E.g., for silicon, a

mixing factor of 12% optimizes not only the gap, but the whole electronic structure (e.g. the valence band

width and all observed direct transitions between the VB and the CB,) as well as the elastic properties.

(The mixing has no significant effect on the binding energies and the lattice constant[34].)

The price for the improved precision is the need to evaluate four-center integrals, which entails a sig-

nificantly larger computational cost and larger memory requirements compared to pure DFT methods,

which make geometry optimizations of surface or interface models using hybrid exchange (HE) imprac-

tical. However, tests have shown, that an additional geometry optimization is only necessary, when the

system contains a shallow defect[34]. Therefore, it is possible to obtain an improvement of the electronic

structure derived from LDA calculations, by making a single HE calculation on LDA-optimized geometry.

2.6 The pseudopotential method

The bonding behavior of atoms in solids is mainly determined by their valence shell electrons. The

inner shell electrons basically act the same way as in an isolated, i.e. unbonded, atom and also screen

the valence electrons from the nuclear charge. It seems to be a promising approach to neglect the core

electrons in the self-consistent calculation of the charge density and instead include them into a new

external potential, called a pseudopotential vps, which now combines the potential of the nuclear charge

with the core electron potential [35–38]. This pseudopotential converts the wave function into a pseudo-

wave-function which is smooth and node-less close to the core but fits the all electron wave function

accurately in the outer region of the atom, and which yields the same valence energy eigenvalues as the

all electron wave function

Eps

v=EAE

v.(2.39)

Forming the pseudopotential entails the substitution of the charge density n(~r) by the valence charge

density nv(~r) that now only contains the valence shell electrons. The effective potential then takes the

form

veff[nv(~r)] = vps(~r) + vc(~r) + vXC[n(~r)] (2.40)

or in the local density approximation

veff(~r) = vps(~r) + vc(~r) + vXC(~r) (2.41)

2.6. THE PSEUDOPOTENTIAL METHOD 13

here vps(~r) is given by the superposition of the single core ion pseudo-potentials

vps =

j=1

vps

j(~r −Rj) (2.42)

where Rjare the positions of the respective nuclei.

The pseudopotential vps can be obtained by calculating the effective potential of a single, isolated atom

and then subtracting the valence electron contributions

vps

j(~r) = veff(~r)−Znv

j(~

r′)

|~r −~

r′|d~r −vXC[nj(~r)].(2.43)

As long as a frozen core is assumed, the same pseudopotential can be used for any spatial distribution of

ionic cores.

In practice, several different forms of pseudopotentials are used:

norm-conserving pseudopotentials [38,39] meet the conditions that the normalized radial electronic

pseudo-wave-function at the angular momentum lis equal to the corresponding all-electron radial

wave function beyond a chosen cutoff radius rcl

Rps

l(r) = RAE

lfor r > rcl (2.44)

or converges rapidly to that value, and that the charge enclosed within rcl for the two wave functions

is equal:

Zrcl

|Rps

l(r)|2r2dr =Zrcl

|RAE

l(r)|2r2dr.

soft or ultrasoft pseudopotentials [40–42] here the restraint of norm-conservation is dropped in order

to be able to create smoother pseudo-wave-functions and thus enable the use of smaller plane-wave

basis sets.

the projector augmented-wave method [43] indirectly uses pseudopotentials. It divides the wave

functions into an atom-centered sphere in which the wave function is represented by partial waves

and an outer (bond-) region, where it is represented by an envelope function, which in turn can be

expanded into an arbitrary basis, usually plane waves.

Here the pseudopotential approach is applied to the construction of the partial-wave basis for the

core region.

Beside the striking advantage of reducing the number of electrons to be considered quantum-mechanically,

and thus greatly reducing the computational cost, the pseudopotential method has another, not so obvious

point in its favor — When total energies are compared, the differences are relatively small compared to

the energy values; therefore small errors, that occur in the treatment of each single electron can sum up to

significant values, in total energy differences. By using pseudo-potentials, the number of treated electrons,

and therefore the number of possible sources of error is reduced. In addition to this, the numerical value

of the total energy is reduced, thus decreasing the effect of possible errors on total energy differences.

As long as only total energies and electronic properties of the valence shell are of interest, as is the case

in this work, the impact of neglecting the interaction of inner shell electrons on the valence shell is very

small. However, the ability to calculate some experimentally observable properties, such as electron spin

resonance (ESR) and nuclear magnetic resonance (NMR) is lost.

14 CHAPTER 2. QM TREATMENT

2.7 Density Functional based Tight Binding

2.7.1 Tight Binding schemes

Tight-binding (TB) approaches work on the principle of treating the electronic wave function of a system

as a superposition of atom-like wave functions. The electrons are tightly bound to the cores in the sense

that the electrons are not allowed to delocalize beyond the confines of a minimal LCAO basis.

In the popular semi empirical tight-binding (SETB) schemes, the Hamilton matrix elements are approx-

imated by analytical functions, the parameters of which are optimized to reproduce experimental data,

as are the atomic basis functions. This allows a much faster calculation of the total energy than DFT or

other ab initio methods, while retaining an ability to account some for quantum mechanical effects. In

such an approach, total energy can, e.g. be expressed as

ETB =

ǫi+1

j6=k

vj,k(|Rj−Rk|).(2.45)

Here vj,k(r) is an interatomic potential depending on the distance. The ǫiare taken to be eigenvalues of

a non self-consistent Schr¨

odinger-like equation

ǫiΨi(~r) = −1

2∆ + V(~r)Ψi(~r),(2.46)

which is solved variationally using the atom-like basis set leading to a secular equation

|H−ǫS|= 0 (2.47)

with the Hamilton and overlap matrices Hand S. In the most stringent application of the tight bind-

ing principle, only the diagonal matrix elements would be regarded, however, this would exclude any

interatomic interactions, so that they are instead assumed to vanish beyond the first or second shell of

neighbor atoms[44]4.

2.7.2 DFTB

In contrast to SETB, the Density Functional based Tight Binding (DFTB) derives its parametrization

from ab initio DFT calculations. In the following I shall sketch the formalism and the approximations

therein. A more detailed description of the method and its derivation can be obtained from the original

works of Seifert, Porezag et. al. [45–47], later reviews [48] and [49]. Ref [50] extends the original method

to include spin polarization effects.

Basic Formalism

Based on the DFT total energy expression after the introduction of the effective single particle wave

functions (equation 2.27),

E[n(~r)] =

occ

ni(~r)ZΨ∗

i(~r)





−∇2

2+1

2Zn(~

r′)

|~r −~

r′|d~

r′+vnuc(~r)

|{z }

H0[n0(~r)]







Ψi(~r)d~r +EXC[n(~r)],(2.48)

4Often, the basis functions are assumed to be orthonormal, so that Sbecomes the unit matrix.

2.7. DENSITY FUNCTIONAL BASED TIGHT BINDING 15

with the normalized occupation numbers nifollowing Janak

n(~r) =

occ

ni|Ψi(~r)|2;N=

occ

ni,(2.49)

the first step is to approximate equation (2.48) by Taylor-expanding it in a norm-conserving density

fluctuation ∆n(~r), around a reference density n0(~r). This expansion is then truncated after the second

order. With ∆n(~r) = n(~r)−n0(~r), ˆ

H[n(~r)] as defined in (2.23) and the internuclear energy ENNthe

Taylor-expansion leads to

E[n(~r)] =

occ

niZΨ∗

i(~r)ˆ

H0[n0(~r)]Ψi(~r)d~r

|{z }

EBS

2Z Z n0(~r)n0(~

r′)

|r−r′|d~

r′d~r +EXC [n0(~r)] −ZvXC[n0(~r)]no(~r)d~r +ENN

|{z }

Erep

2Z Z "1

|~r −~

r′|+∂2EXC [n(~

r)]

∂n(~

r)2#n0

∆n(~r)∆n(~

r′)d~

r′d~r

|{z }

E2nd

+O(∆n(~

r)3).(2.50)

Here the terms in the first two lines only dependent on n0and the third line contains the second order

contributions and truncation error. (Note that the first order terms cancel out.) More concisely, the total

energy is expressed as

E[n]≈ EDFTB =EBS +Erep(R) + E2nd.(2.51)

In systems in which interatomic charge transfer and long range Coulomb interactions play no significant

role, E2nd can be neglected. The energy expression then only depends on the reference density n0(~r),

leading to the non self-consistent version of DFTB. Otherwise, E2nd can be calculated self-consistently[48].

This addition to the original DFTB method is called self-consistent charge DFTB (SCC-DFTB).

Following the tight-binding philosophy, all interactions are assumed to be short-ranged and pairwise, no

three- or more center interactions are regarded.

Hamilton matrix elements

The elements of the zeroth order Hamilton matrix H0depend on the basis functions used to expand the

Ψi. DFTB uses a minimal localized valence basis set

|Ψi(~r)i=X

cνi φj

ν(~r −Rj),(2.52)

where the φj

νare taken to be linear combinations of Slater orbitals obtained from a self-consistent DFTB

reference calculation using a neutral, spherical symmetric pseudo atom. The reference density then

becomes the superposition of atomic reference densities

n0(~r) =

ν∈j

νφj

ν(~r)φj

ν(~r)=

0(~r).(2.53)

Within the two-center approximation and neglecting the crystal-field terms, the Hamilton matrix elements

thus become

µ,ν =









ǫfree atom

νµ=ν

Dφµ|ˆ

t+ˆ

V[nj(µ)

0+nk(ν)

0Eµ6=ν, j(µ)6=k(ν)

0µ6=ν, j(µ) = k(ν)

16 CHAPTER 2. QM TREATMENT

with j(µ) and k(ν) the corresponding atoms to the atomic orbitals µand ν. The diagonal elements of

this matrix would be the respective eigenvalues of the DFT reference calculation. However, in order

to obtain the correct dissociation limit, the eigenvalue of the free atom is used. The non-zero matrix

elements are calculated element-pairwise and tabulated over a range of interatomic distances. At the

same time, also the overlap matrix elements are calculated and tabulated. In practice, the wave function

and charge density of free atoms (or dimers) are not well suited to describe the more localized electrons

of condensed matter. Therefore, during the DFT reference calculations, the wave functions and densities

are artificially compressed to resemble the circumstances in the target system.

Erep(R)

In the same two-center approach as for the Hamilton matrix elements, the repulsive energy term is

expressed as a sum of pairwise, radial symmetric atomic contributions

Erep(R) = 1

j6=k

vj,k

rep(|Rj−Rk|).(2.54)

After the Hamilton matrix elements have been calculated, the repulsive potentials vj,k

rep(r) correct the

difference between the DFTB band structure energy EBS and the total energy from the reference DFT

method. They can be calculated from the energy differences over a range of interatomic separations

vj,k

rep(rj,k) = Etot

DFT(rj,k)−EBS(rj,k)−E2nd(rj,k).(2.55)

(E2nd ≡0 if second order terms are neglected.)

In practice, the repulsive potentials thus obtained often do not yield satisfactory accuracy. Therefore the

repulsive potentials are rather constructed by fitting a predefined function to the atom–atom distance

dependent energy differences in more extended systems, e.g. the C–C distances in HC≡CH, H2C=CH2

and H3C–CH3. By this approach, the error from neglecting multicenter and crystal-field terms is partly

corrected in the repulsive potentials

vj,k

rep(r)≈Etot

DFT(r)−EBS(r)−E2nd(r)fit systems .

Evolutionary algorithms can help to fit vj,k

rep(r) give an accurate reproduction of different reference DFT

methods [51].

Second order corrections

In systems which show significant charge transfer or long-range Coulomb interactions, the density-

fluctuation dependent terms summarized in E2nd must be taken into account. To this end, ∆nis

decomposed into atomic contributions, as it is done for n0in equation 2.53

∆n(~r) =

l,m

lmFj

lmγj

lm ≈

00γj

00∆qj,(2.56)

where truncation after the monopole term accounts for the most important contributions and the coeffi-

cients cj

lm are set to fluctuations of the atomic Mulliken charges defined as

∆qj=qj−qj

0(2.57)

with

∆qj=

niX

µ,ν∈j

cµ,icν,iSµ,ν.(2.58)

2.7. DENSITY FUNCTIONAL BASED TIGHT BINDING 17

Here, qj

orepresents the number of valence electrons of the neutral atom j. Substituting this into the

expression for E2nd given in Equation 2.50 leads to

E2nd =1

j,k

γjk(|Rj−Rk|)∆qj∆qk

where γjk =ZZ 1

|r−r′|+∂VXC[n]

∂n(r′)Fj

00(rj)Fk

00(r′

k)1

4πd~rd~

r′.

Here, the F00(~r) are fixed and radial orbital relaxation functions, thus the γjk depend only on the

interatomic distances. This assumes, that the electrostatics of the system can be described by atom-

centered monopoles. The second-order corrections to DFTB are discussed in detail in Refs [48].

Total energy and secular equation

With the approximations and contributions discussed so far, the DFTB total energy expression from

equation 2.48 becomes

Etot =

occ

µν

c∗

µicνih0

µν +1

γjk(Rjk)∆qj∆qk+Erep(R).(2.59)

The total energy is minimized by variation of the LCAO coefficients c. This leads to a eigenvalue problem

with the secular equation

cνi (hµν −ǫνSµν) = 0 ∀i, µ. (2.60)

As the cνi and the atomic Mulliken charges depend on each other, this equation must be solved self-

consistently.

2.7.3 Extensions to SCC-DFTB

As a prerequisite to the development of the QM/MM embedding schemes laid out in detail in chapter 4,

further extensions must be added to the SCC-DFTB method: the influence of an external electrostatic

field and the possibility to constrain atomic charges during the SCC iterations.

External electric field

Similar to the SCC-DFTB QM/MM coupling approach introduced in Ref [52] for biological molecules, the

DFTB Hamiltonian was extended by an additive term describing the coulomb interaction of the DFTB

treated atoms with external monopoles. The energetic contribution UXof external point charges in a

Coulomb energy of the form

Ux=

ext

atoms

rIA

qA·qI.(2.61)

Starting from the derivative of the Mulliken populations with respect to the ciλ

∂qA

∂ciλ

2ni

X

δB(λ)AciνSλν +X

µ(A)

ciµSλµ

(2.62)

18 CHAPTER 2. QM TREATMENT

we can formulate the derivative of UX

∂UX

∂ciλ

=∂

∂ciλ "X

rIA

qA·qI,#="X

rIA

∂qA

∂ciλ #(2.63)

rIA

qI·



X

δB(λ)AciνSλν +X

µ(A)

ciµSλµ





(2.64)

qI



X

rIA

δB(λ)AciνSλν +X

µ(A)

rIA

ciµSλµ





(2.65)

ν(X

qI·1

2Sλν 1

rIA(ν)

rIB(λ)·ciν ),(2.66)

which leads to the expression for the matrix element contribution of the external point charges:

µν =X

2Sµν 1

rIA(µ)

rIB(ν).(2.67)

The external point charges can be replaced by Gaussian charge distributions of variance σ2to better

represent charged atoms in the surroundings, which are not point-like but spatially extended. In that

case, the coulomb energy transforms to

gauss =qAqI

rIA

·erf rIA

σ(2.68)

and the respective matrix element contributions become

µν =X

2Sµν erf rIA(µ)

σ·1

rIA(µ)

+ erf rIB(ν)

σ·1

rIB(ν).(2.69)

Effectively, this scales the influence of external charges on close QM atoms down. The choice of a single

Gaussian to represent the spatial extent of the external charge distribution is arbitrary and based on

mathematical convenience.

Charge constraints

To achieve a bulk-like electronic structure of an embedded cluster, it can be helpful to constrain the

electronic structure of certain atoms of such a cluster, usually that of the boundary atoms. Since the

varied quantity in DFTB is the Mulliken population, this is where the electronic constraint must be

applied. This is done by adding a penalty potential Uconstr to the total energy expression of the form [53]

Uconstr =X

UA(qA−q0

|{z }

∆qA

)2.(2.70)

Here, UAis an atom-wise potential strength, qAis the Mulliken population of atom A,q0

Aits target

population and ∆qAits charge deviation. The quadratic form ensures a symmetric penalty for positive

and negative deviations. As the penalty potential affects the atomic charges, it leads to a perturbation

of the Hamiltonian which can be calculated by variation of Uconstr in the LCAO coefficients cλk of the

basis orbital λin eigenvector k:

∂Uconstr

∂c∗

λk

.(2.71)

2.7. DENSITY FUNCTIONAL BASED TIGHT BINDING 19

Starting from the derivative of the Mulliken population qA

∂qA

∂c∗

λk

∂

∂c∗

λk

occ.

ni

X

µ(A)X

c∗

µicνiSµν +X

µ(A)X

cµic∗

νiSνµ

(2.72)

in its symmetric form that ensures real population values at all k-points:

∂qA

∂c∗

λk

occ.

ni

δX

δkiδλµc∗

µicνiSµν +X

µ(A)X

cµic∗

νiSνµ

(2.73)

2nk

X

δB(λ)Ackν Sλν +X

µ(A)

ckµSλµ

(2.74)

the derivation in equation 2.71 can now be executed:

∂Uconstr

∂c∗

λk

UA2δqA

∂qA

∂c∗

λk

ΩA

∂qA

∂c∗

λk

(2.75)

2nkX

ΩA

X

δB(λ)Ackν Sλν +X

µ(A)

ckµSλµ

(2.76)

2nk

X

νX

ΩAδB(λ)Ackν Sλν +X

µ(A)

ΩAckµSλµ

(2.77)

2nk"X

ΩB(λ)ckν Sλν +X

ΩA(ν)ckµSλµ#(2.78)

2nkX

ckν Sλν ΩB(λ)+ ΩA(ν)(2.79)

Thus, the perturbation becomes

∆Hconstr

µν =Sµν UA(µ)∆qA(µ)+UB(ν)∆qB(ν).(2.80)

In addition to this, the constraint potential can also lead to a perturbation in the forces, which can be

derived in a similar manner to the Hamiltonian perturbation:

∆Fc=−∂Uconstr

=−X

UA2∆qA·∂qA

∂RC

with Eq 2.72 (2.81)

=−X

UA2∆qA

1

niX

µ(A)X

c∗

µicνi

∂Sµν

∂Rc

+cµic∗

νi

∂Sνµ

∂Rc

(2.82)

=−X

niX

µX

ν"X

c∗

µicνi

∂Sµν

∂Rc

+cµic∗

νi

∂Sνµ

∂Rc#UA(µ)qA(µ)(2.83)

=−X

ni(X

µX

c∗

µicνi

∂Sµν

∂Rc

UA(µ)∆qA(µ)(2.84)

νX

cνic∗

µi

∂Sµν

∂Rc

UA(ν)∆qA(ν))

niX

µX

c∗

µicνi

∂Sµν

∂RcUA(µ)∆qA(µ)+UA(ν)∆qA(ν).(2.85)

20 CHAPTER 2. QM TREATMENT

Equations 2.80,2.85 show that both the energy and force perturbation from the electronic constraint

potential vanish, if the constraints are well fulfilled after electronic convergence.

2.8 Methods of QM/MM coupling

My primary interest rests on the description of reactive events at hybrid interfaces, which imposes the ne-

cessity to model such interfaces with sufficient realism. The size of the models needed to properly describe

a surface or interface leads to extreme computational cost. However, at solid/liquid or solid/polymer hy-

brid interfaces, the large difference in rigidity between the solid and liquid or polymer sides causes large

parts of the solid surface to be subject to only small mechanical distortions by a single reactive event.

Therefore, this part of the surface can be regarded as an external bath for the QM calculation. Classical

force fields are in most cases sufficiently accurate to model this bath. In many cases, even fixing the

external bath, thus avoiding to model its internal interactions at all, can give satisfactory results.

To couple theoretical modeling at different levels of detail within a single system, two different approaches

are currently in use: additive and subtractive schemes. Important implementations of the latter are the

ONIOM[54] method and its predecessor IMOMO[55]. Here, in the first step, energy and gradients for the

whole system are calculated using the less complex low-level method. Subsequently, energy and gradients

of the part of the system to be treated with the more detailed high-level method, the high-level or QM

cluster, are calculated with both the low-level and high-level methods. Finally the low-level contributions

of the high-level cluster are subtracted from the full system result and the high-level energy and gradients

are added:

EONIOM

full =Elow

full −Elow

hlc +Ehigh

hlc (2.86)

FONIOM

full =~

Flow

full −~

Flow

hlc +~

Fhigh

hlc (2.87)

where the superscripts denote the low-level, high-level and ONIOM methods, the subscripts denote the

full system and the high-level cluster (HLC). A major drawback of this method is, that both parts of the

system must be described using the low-level method. This makes the subtractive scheme unsuitable for

the QM/MM examination of reactive events, unless a reactive force field is available. However, such force

fields are rather difficult to handle and their development is not as far advanced as for the non-reactive

case. Based on these considerations, I choose to focus on the additive coupling scheme:

In contrast to the subtractive scheme, the additive scheme directly couples the low- and high-level zones,

by partitioning the Hamiltonian:

H=ˆ

Hhigh +ˆ

Hlow +ˆ

Hhigh/low (2.88)

which avoids low-level calculations on the high-level zone, but on the other hand imposes the problem of

formulating a suitable high-/low-level coupling Hamiltonian.

The popular hybrid quantum mechanical/molecular mechanical (QM/MM) [56–63] approach combines

the advantages of quantum mechanical methods in describing chemically active regions with a description

of the remaining system by interatomic potential functions as force fields [64] or by ion-pair potential

functions [65].

Since the introduction of the QM/MM idea around the mid-1970s [56], numerous studies have been

conducted examining both performance of different schemes as well as interesting applications to various

solution and catalytic systems [61,62,66–80]. These previous studies have clearly demonstrated that

carefully applied QM/MM methods can provide useful insights into chemical mechanisms in complex

systems that are difficult to obtain otherwise.

In commonly used QM/MM schemes, the Hamiltonian operator of the entire system, ˆ

H, is written as

the sum of those for the QM partition, ˆ

HQM, the MM partition, ˆ

HMM and the interaction between

2.8. METHODS OF QM/MM COUPLING 21

the two, ˆ

HQM/MM, The precise expression of ˆ

HQM/MM varies but generally has contributions from

electrostatic, van der Waals, bonded interactions, and possibly additional constraints [57],

HQM/MM =ˆ

Helec

QM/MM +ˆ

HvdW

QM/MM +ˆ

Hbonded

QM/MM +ˆ

Hcons

QM/MM (2.89)

The bonded terms and constraints (e.g., fixed bond distance between boundary QM and MM atoms) are

used for keeping the proper connectivities and geometries when cutting through covalent bonds at the

QM/MM interface. The QM/MM van der Waals terms can be optimized to improve properties such as

distribution of MM groups around the QM group [81]. The electrostatic component, which is missing in

some early implementations for organometallic systems [60], is crucial for the investigation of reactions

in polar surroundings like water or many solids.

2.8.1 The QM/MM boundary

Figure 2.1: Illustration of a covalent

bond crossing the QM/MM bound-

ary, showing the MM host atom

(MMH), QM host atom (QMH) and

QM link atom (QML).

An issue that has been repeatedly raised, concerns the treatment

of the QM/MM boundary. The interaction between MM atoms

and nearby QM atoms should be carefully treated to reliably

describe the effect of the environment on chemical properties

of the QM region. Complications in the calculation of equa-

tion 2.88 arise in particular if the QM/MM partition involves

dividing the system across covalent bonds.

The large number of proposed hybrid methods differ in the

way this problem is dealt with and which contributions of

equation 2.89 are calculated quantum mechanically. The most

straightforward approach involves saturating dangling bonds by

link atoms, which are typically chosen as hydrogen atoms be-

tween the QM host atom (QMH, see Fig. 2.1) and the MM

host atom (MMH) [57]. The link atom (QML) is treated at

the QM level, and may be subject to an angular and distance

constraint to lie along the bond between QMH and MMH at

a fixed bond distance. The link atom typically interacts with

MM atoms through electrostatic terms but not through van der

Waals terms.

Instead of regular hydrogen atoms, hydrogen-like atoms or

pseudo halogens have been used to terminate the QM region [82,83]. In those approaches, the link

atom coincides with the MMH. The electronic nature of this atom is modified to mimic the behavior

of the MM host atom or MM host group (MMHG). Instead of using pseudo-atoms as QML, the bond

QMH–QML distance is often scaled so smaller values than their equilibrium bond length. This approach

is taken from the general practice of constructing 2D-slab models of solid state interfaces, where hydrogen

atoms are often used to terminate the dangling bonds on the backsides of the adjoined layers. As the

symmetric-antisymmetric splitting of two molecular orbitals is proportional to the overlap Sof the un-

derlying atomic orbitals, which in turn varies inversely with the bond distance, reducing the QMH–QML

distances can often shift spurious gap states originating from the QMH–QML bonds out of the band gap.

Other approaches to the QM capping problem include the use of hybrid orbitals [56] such as the localized

self consistent field (LSCF) [84–87] and the generalized hybrid orbital method (GHO) [88,89].

Here, I write the potential energy of the QM/MM system as (see equation 2.88):

Utot =DΨˆ

HQM +ˆ

Helec

QM/MM ΨE+ (2.90)

UvdW

QM/MM +Ubonded

QM/MM +UMM

22 CHAPTER 2. QM TREATMENT

where Ψ is the electronic wave function of the QM region using the SCC-DFTB method. The operator

describing QM/MM electrostatic interactions, ˆ

Helec

QM/MM , has the form of Coulomb interaction between

the MM point charge QAand the Mulliken charge ∆qBon the QM atom [52]

Helec,pointcharges

QM/MM =X

A∈MM X

B∈QM

QA∆qB

|rA−rB|(2.91)

(cf. equation 2.67).

In practice, the QML often fall very close to the MMH, which can cause severe problems in the quantum

mechanical description. Many different approaches have been developed to overcome this difficulty in

biochemical applications [61,90], however, most of these are not applicable to the solid state systems,

which are the primary interest of my work. One example is the charge deletion approach,relies on the

(formal) neutrality of groups within the examined amino acids or monomers. The charge of the MMH

atom is simply deleted and the neutrality of the group is restored by applying compensating charges to

the other atoms in the group.

In contrast, the approach to apply a Gaussian broadening to the point-charge distribution of the MM

atoms, however, can be transferred without difficulties. The idea of this approach is to reduce the short-

range influence of an atomic charge, while retaining the point-charge character at long distances. This is

achieved by applying a Gaussian broadening of standard deviation σto the point charge. This is justified

by the fact that atomic point charges actually represent a point-like nuclear charge which is surrounded

and (partly) shielded by a spatially extended electron distribution. The mathematical formulation by

means of a Gaussian function is chosen for the sake of mathematical simplicity: the coulomb interaction

is simply scaled by a an error-function prefactor, so that the coupling Hamiltonian becomes:

Helec,gaussharges

QM/MM =X

A∈MM X

B∈QM

QA∆qB

|rA−rB|·erf |rA−rB|

σ.(2.92)

(cf. equation 2.69) The property of the error-function to vanish at x→0 and quickly converge to unity

at x > σ, scales down the influence of a close-by point-charge while maintaining the influence of the

further parts of the charge distribution. Since the error-function term converges to unity for σ→0,

equation 2.91 can be regarded as the limit of equation 2.92 for vanishing size of the charge-distribution. σ

should be chosen appropriately, to alleviate the quantum mechanical problems between QML and MMH,

while keeping the coulomb interaction between QMH and MMH intact.

Chapter 3

Modeling Chemical Reactions

3.1 Reaction energetics at 0 K

At zero temperature, the transition between two stable states of a system follows a minimum energy

path (MEP) on the system’s potential energy hypersurface (PES).It is defined as the path connecting the

minima in a manner that the potential energy gradients vanish in all directions orthogonal to the path.

The MEP is the transition route which requires the least energy to follow. At finite temperatures the

system can deviate from the MEP with the extent of the deviations governed by the relation between the

depth of the MEP valley and the thermal energy. Maxima along the MEP and their associated geometries

are termed transition states1and are saddle points of first order on the PES. For a single step reaction,

the energy difference between the reactant and the transition state is called the barrier energy

Eb=Emax

tot −

reactants

tot.(3.1)

(Note that by definition Eb≥0.)

The energy difference between two local minima defines the reaction energy ∆Eof the corresponding

reaction. By convention, it is defined as

∆E=

products

tot −

reactants

tot.(3.2)

A reaction with negative reaction energy is called exothermic, as energy is released by the system. Analo-

gously, a reaction with positive energy is called endothermic. The terms exothermicity and endothermicity

are synonyms for the energies of exothermic and endothermic reactions, respectively. The calculation of

∆Eof any reaction is straightforward, it only requires the calculation of product and reactant energies

by means of geometry optimization of the respective compounds.

3.1.1 Finding barrier energies

Over the years, a large number of methods to find the saddle point(s) of a given transition on a specific

PES have been developed. The conventional approach is to search for a system configuration in which

1Some authors distinguish between the transition state and the transition structure [2], in which the former refers to a

free energy path and only the latter to an MEP. I choose not to follow this convention, as it is not followed consistently

throughout the literature and prone to confusion. Therefore I prefer to clarify in each case, whether I am referring to free

or potential energies.

24 CHAPTER 3. MODELING CHEMICAL REACTIONS

−2 −1.5 −1 −0.5 0 0.5 1 1.5

−1

−0.5

0.5

1.5

2.5

−2 −1.5 −1 −0.5 0 0.5 1 1.5

−1

−0.5

0.5

1.5

2.5

a) b)

Figure 3.1: Illustration of two adiabatic mapping path-searches in a M¨

uller-Brown Potential [91]

using differently defined reaction coordinates. The parallel solid lines mark the lines of minimization

perpendicular to the chosen reaction coordinate. After [92].

its matrix of energy second derivatives with respect to the atomic degrees of freedom rj,j∈1...3M







∂2Etot

∂r1∂r1

∂2Etot

∂r1∂r2··· ∂2Etot

∂r1∂r3M

∂2Etot

∂r2∂r1

....

.....

∂2Etot

∂r3M∂r1··· ··· ∂2Etot

∂r3M∂r3M







,(3.3)

called the Hessian matrix, has exactly one negative eigenvalue. This so-called surface walking is performed

by following the Hessian’s eigenvectors, which can be very computationally demanding if analytical second

derivatives are not available. Since the numerical evaluation of the Hessian demands 2(3M)2total energy

calculations, such algorithms are impractical for systems with a large number of degrees of freedom.

Interface reactions are such problems, since a realistic description of the surfaces generally demands a

large number of movable atoms.

One class of methods developed to reduce the computational costs are path sampling or chain-of-states

methods. They explore reaction paths and do not aim at locating the exact minima and maxima along a

reaction pathway, but rather on sampling the total energies along the path. The simplest of these methods

is the adiabatic mapping approach, but other methods have been developed to overcome the limitations

of this rather simplistic approach. These methods work on the principle of iteratively refining an initial

reaction coordinate guess, either during the construction of the reaction path, as in the growing string

method [93], or by refining an initial guess of the whole path, as in the nudged elastic band approach,

outlined in section 3.1.2 below.

Adiabatic mapping [94,95], also referred to as the drag method, scans the PES along one or more pre-

determined reaction coordinates, while minimizing the total energy in all orthogonal degrees of freedom.

This reduces the problem of finding the MEP to a number of constrained geometry optimizations equal

to the desired number of samples along the path, completely avoiding the calculation of second deriva-

tives for the path search. However, the necessity to define a suitable reaction coordinate in advance is

a serious drawback, as it can be very challenging to do so without excluding large portions of the PES

a priori. This problem is illustrated in Fig 3.1. Note that, although the choice seems obvious in this

simple, two-dimensional case, in practical systems with a large number of degrees of freedom, this can

be very different. Additionally, the saddle point is not obtained directly but only approximated by the

samples closest to it, so that the barrier energy and structure must be interpolated. This limitation is,

3.1. REACTION ENERGETICS AT 0 K 25

however, offset by the robustness of these methods in terms of convergence[94] and their inherent ability

to describe more complicated concerted reactions, which exhibit several saddle points along along the

MEP.

More sophisticated techniques of surface walking [96–101] rely on the educated guess that the eigenmode

of the smallest force constant leads to a saddle-point. Therefore they rely on optimizing the geometry

to a saddle point, starting from the reactant configuration, climbing uphill along the smallest frequency

eigenmode. But still, these methods require the costly calculation of the full Hessian.

To overcome this problem, a large number of approaches have been developed. They comprise, among

others, the method by Dewar, Healy and Steward[102], the ridge method by Ionova and Carter[103] or the

dimer method[104] by Henkelman and J´onsson, but also the conjugate peak refinement algorithm (CPR)

by Fischer[105]. The first three methods use a pair of images which is advanced towards the saddle point,

in an attempt to place the images close to each other on both sides of the transition state. The idea

behind these methods is to perform surface walking without the calculation of the whole Hessian. Instead,

the second derivative is only calculated in the direction of smallest curvature, using a finite-difference

approximation. An improved version of the dimer method by Heyden and Keil [106] uses a modified

representation of the dimer to avoid two total energy calculations per geometry iteration.

The CPR method relies upon line minimizations along linearly interpolated reaction path segments.

Hence it is able to describe reaction paths with several local minima and maxima, yielding the exact

transition state(s) and intermediary configurations. It can be regarded as a hybrid method between

surface walking and the chain-of-states methods described below.

3.1.2 Nudged elastic band

Within the category of the chain-of-states methods, the nudged elastic band[107] (NEB) approach is one

of the most prominent. As it is the main path-search method used in this work, I shall describe it in

detail. In NEB, the intermediate states are called images. They are are connected by virtual springs to

distribute the images evenly along the MEP. NEB differs from other chain-of-states methods[94] in that

along the path direction, only spring forces are regarded, the tangential components of the calculated

atomic forces are neglected. The chain of images is then optimized to minimize the NEB force for all

intermediate images along the path, usually under the constraint of fixed start- and end-geometries. The

even spacing of images along the sampled MEP helps to avoid the so-called valley-hopping problem: for

many MEP search methods it is not assured, that the path segments on the product and reactant sides

belong to the same PES valley. In that case, the maximum along the calculated MEP does not refer to

the reaction barrier, but to the (higher) barrier of changing between the two reaction mechanisms (cf.

figure 3.1, which illustrates the problem for adiabatic mapping).

As stated in section 3.1.1, the barrier energy must be interpolated when using a path-sampling method.

This interpolation is frequently performed by a cubic spline fit of the image energies, with the gradients

at the start- and end-images set to zero. Another approach is to fit a cubic polynomial to each segment

connecting two images, using the tangential real force (i.e. calculated force) components to determine

the energy gradient at each image[108].

A more elaborate version of NEB is the ANEB method[109], which uses adaptive spring forces proportional

to the image energy, in order to achieve an aggregation of the images around the transition state. This

improves the quality of the MEP estimate in the vicinity of the transition state, without imposing the

greater computational cost of using a larger number of images, by locally raising the sampling resolution.

The Climbing Image NEB method (CI-NEB)[110] uses the image closest to the transition state (identified

by its total energy), as a starting geometry for a saddle point refinement, using the rest of the MEP

estimate to obtain a local reaction coordinate at the image being moved to the transition state and

inverting the gradient component, parallel to this coordinate. In this context, it can be understood as a

synthesis of standard NEB and a surface walking method. Similarly, it is a frequently employed approach,

to use the highest energy NEB image as the input configuration for a surface walking algorithm.

26 CHAPTER 3. MODELING CHEMICAL REACTIONS

−2 −1.5 −1 −0.5 0 0.5 1 1.5

−1

−0.5

0.5

1.5

2.5

−2 −1.5 −1 −0.5 0 0.5 1 1.5

−1

−0.5

0.5

1.5

2.5

−2 −1.5 −1 −0.5 0 0.5 1 1.5

−1

−0.5

0.5

1.5

2.5

−2 −1.5 −1 −0.5 0 0.5 1 1.5

−1

−0.5

0.5

1.5

2.5

a) b)

c) d)

Figure 3.2: Illustration of an NEB path optimization in a M¨

uller-Brown Potential [91]. a) initial

guess b) intermediate path estimate, c) final path estimate, d) path interpolation between images.

3.2. BEYOND THE MINIMUM ENERGY PICTURE 27

3.2 Beyond the minimum energy picture

For many applications, the reaction energies and MEPs provide sufficient information. Naturally, the

quasi-static picture, neglecting temperature effects, is most suitable to situations in which reaction en-

ergonicities and barriers are large, when compared to the thermal energies. If, however, that is not the

case, e.g. when solvent effects or van-der-Waals bonding are important, entropic effects must be taken

into account and therefore free energy differences and paths must be calculated.

The Gibbs free energy G, in the following just free energy, is defined as

G=U−T S +P V, (3.4)

where the internal energy U=Etot +Ekin is the sum of potential and kinetic energies, Tthe absolute

temperature, Sthe entropy and P V the volume work needed to “make room” for the system at constant

pressure P. Note, that the last term is only meaningful, if the simulated system interacts with an outside

system via an identifiable wall. The free energy difference ∆Gof a reaction gives the outside work that can

be extracted from or must be invested into a reaction, the free energy barrier of a reaction describes the

amount of outside work necessary activate it. Via the entropy term, the free energy can differ significantly

from the (non-free) reaction energy and barrier.

Thermodynamic simulation Macroscopic thermodynamic properties, are always averages of their

microscopic constituents. Therefore, either a time- or ensemble average must be performed to calculate

them. The latter approach is the basis of Monte-Carlo (MC) methods, first employed by Metropolis

et. al. [111] in which a statistical sampling of configurations according to a Boltzmann distribution is

performed. I shall not go into any more detail here, as I do not use atomistic Monte-Carlo methods in

this work.

The other possibility to obtain a thermodynamic average is by time-averaging over the trajectory of the

examined system. To obtain these, the equations of motion of the pertinent system must be solved. In

the case of a chemical reaction, these are the EOMs of the atomic nuclei. Most frequently, the atoms are

treated classically, using Newton’s EOM

R(t) = 1

dEtot

dRt

,(3.5)

regardless whether the interatomic forces are calculated quantum-mechanically or by means of classical

(empirical) force-fields. Thus a time-resolved trajectory of the molecular dynamics (MD) is obtained.

The EOM are generally solved by integrating them under the starting conditions of some initial atomic

positions, which are often chosen as a minimum energy configuration, and a set of Boltzmann-distributed,

random atomic velocities. The correct dynamics of the system can only be captured, if the integration

time-step is less than 1

2of the shortest atomic vibration period in the system2, which leads to time steps

of ∼10−15 sec, if hydrogen is present in the system. This severely limits the time-scales accessible by

molecular dynamics. Approaches to overcome this limitation are based on the assumption that certain

internal vibrations of a molecule are irrelevant to the dynamics on the timescale of interest. In this case,

they can be neglected, either by constraining the relative positions of the respective atoms [2], or by

coarse-graining approaches, in which whole groups of atoms are mapped to bead-like superatoms [115–

117]. Nevertheless, explicitly following the time evolution of a system severely limits the accessible time

scales, typically to a maximum of a few ps with ab initio energies and forces and <1µs with classical

force fields.

3.3 Modeling rare events

2This can be understood in the context of Nyquist-Shannon-Whittaker sampling theorem [112–114].

28 CHAPTER 3. MODELING CHEMICAL REACTIONS

Figure 3.3: Illustration of a rare reactive event.

The system performs a large number of vi-

brations, attempting to cross the barrier. Af-

ter [118].

When faced with rare reactive events, i.e. transi-

tions between minima separated by barriers high

in energy with respect to thermal energies at tem-

peratures of interest, the methods described in

section 3.2 generally fail to describe the time-

evolution of the system from state to state3. Fig-

ure 3.3 illustrates the problem: The system un-

dergoes a large number of vibrations, while most

of these attempts to cross the barrier remain un-

successful. It is not only uneconomical to sim-

ulate each of these, but in many cases the com-

putational cost is prohibitive. At the same time,

the sampling of the PES in the vicinity of the

transition state is coarse, so that the quality of

the thermodynamical averaging in this crucial re-

gion is relatively low. The latter problem can

be solved by biased MD techniques, like umbrella

sampling [2], however, they do not give access to

longer time scales.

A very successful ansatz to model the time development of such rare events is the Kinetic Monte Carlo

(KMC) method [118,119], which as been applied to a wide range of topics, such as radiation damage

annealing [120], surface adsorption, diffusion and growth [121–123], or semiconductor processing [124].

The basis of this method is to solve the time-scale problem, by operating on a dynamic time step. In

contrast to the MD technique, where the system is developed along the time axis by integrating the EOMs

over a pre-determined time step, here the time the system stays in its current state ibefore making a

transition, is determined. To this end, the rate constants kij for the transitions from state ito each

adjacent state jare calculated. The probability pij(t) that a system has not undergone the transition to

state jis assumed to be

pij(t) = kije−kij t.(3.6)

This assumption is true, as long as the system stays in each state ilong enough to “forget” from which

state it arrived, i.e. subsequent transitions are uncorrelated. This is true for rare events, as defined above.

If all rate constants for all transitions are known, the residence time of a system in state iand the final

state of the transition can be determined randomly. This way, a trajectory of the system trough the

minima of its PES can be determined. The great advantage of this method is, that not every potential

transition must be simulated, as in MD, but only actual transitions. The crucial point of this method,

however, is that all possible transitions with their associated rate constants must be known in advance.

The difference between KMC and the MC methods, mentioned in section 3.2, is that the latter aim at

sampling the PES of a system with thermodynamic weighting to calculate the free energies of certain

points of the PES. KMC, in contrast, aims at sampling the kinetics of the system between them.

Therefore, KMC does not only work on a (usually) larger time-scale than MD, but also on a larger length

scale — the degrees of freedom along the pathway(s) are only sampled at the minima and saddle points

while the degrees of freedom normal to the path are neglected.

The rate constants for each pathway can be calculated from the potential energy surface, using the

transition state theory4[125,126] (TST). The rate constant kij is taken to be the equilibrium flux

through a dividing surface, separating states iand jin phase space [118]. As an equilibrium property,

the TST kij can be determined without examining dynamical trajectories, it is simply proportional to the

Boltzmann probability of the system being at the dividing surface, relative to the probability of being in

state i. Often, the harmonic transition state theory (HTST) is used. Here it is assumed, that the PES is

3Here, a state is defined as a local minimum of the PES, separated by barriers from other local minima.

4Or, more accurately, chemical transition state theory, which should not be confused with Slater transition state theory.

The latter describes the transition of electrons between electronic eigenstates of a system.

3.3. MODELING RARE EVENTS 29

harmonic in the vicinity of the local minima iand j, as well as in all directions normal to the transition

at the transition state between them. In this case, the rate constant becomes

kHTST

ij =Q3N

iνmin

Q3N−1

iνsad

e−Estatic

kBT,(3.7)

where Estatic

bis the static barrier height of the transition, and νiare the non-imaginary mode frequen-

cies at the PES minimum and saddle points. As can be seen in equation 3.7, the temperature depen-

dence is strictly exponential. Experience shows, that the prefactors are frequently within the range of

1012 ...1013 sec−1, wherefore it is often arbitrarily set to a value from that range, to save the effort of

calculating all vibrational modes at the minima and saddle points[118].

Recent developments aim at coupling KMC simulations to continuum descriptions, e.g. using phase-field

models [127]. On the one hand, this enables the simulation of even larger length scales, as with KMC

alone, on the other hand, it allows to include processes like diffusion within the bulk phases on both sides

of an interface, without having to model these in the KMC framework, the computational cost of which

is often prohibitive.

30 CHAPTER 3. MODELING CHEMICAL REACTIONS

Chapter 4

QM/MM Embedding of DFTB

clusters at solid surfaces

The appeal of QM/MM embedding schemes lies in their capacity to describe an extended system with

nearly quantum-mechanical accuracy at much lower computational cost than a fully quantum mechanical

calculation. Figure 4.1 shows the mechanical and electronic changes of a γ-Al2O3surface along the MEP

of the adsorption of a small organic molecule. It demonstrates the local character of the mechanical

and electronic changes during the examined reactions and so proves that it is indeed justified to treat

different parts of this system at different resolutions. These representative results were obtained from

fully quantum mechanical calculations using SCC-DFTB (cf. chapter 5). The surface atoms are only

displaced significantly within a very short range around the adsorption site (cf. figure 4.1(a)). The

movement further away from the adsorbate is limited to re-orientations of polar OH groups, which are

driven by the Coulomb interaction and can be well described using molecular mechanics. At the same

time, significant charge transfer virtually only occurs within the adsorbate and up to its second-nearest

neighbors on the surface (cf. Figure 4.1(b)). They suggest, that full quantum mechanical precision is

only required within a limited range around the adsorption site, while the influence of the largest part

of the surface is electrostatic. These parts of the surface undergo only small movements and electronic

changes and can therefore be treated by classical force-fields or even kept fixed.

By determining, which part of the system requires which level of detail in its description, and only

applying methods of the required level to each part, QM/MM coupling allows to obtain the desired

(a) RMSD in ˚

A (b) max(dq) in e

Figure 4.1: RMS displacement and maximum atomic charge difference during the adsorption of

AMEO(hyd) on γ-Al2O3(cf. chapter 5).

32 CHAPTER 4. DFTB QM/MM EMBEDDING

results economically. The techniques to calculate total energies and forces in a QM/MM coupling scheme

have already been described in section 2.8. In this chapter, I shall focus on the treatment of the QM/MM

boundary on the QM side of such an additive coupling scheme. The MM part of the calculation is modeled

as a distribution of fixed charges, and the mechanical coupling by fixing the positions of the QMH and

QML atoms.

4.1 Challenges of QM/MM embedding in solids

(a) surface slab model with adsorbate (b) smal QM cluster

Figure 4.2: Illustration of the bonding situation of a fictitious QM cluster, cut from a surface slab

model of γ-Al2O3around the adsorption site of a small organic molecule.

QM/MM embedding approaches have been applied very successfully to biochemical and polymer systems

for several years (cf. section 2.8). Characteristically, the primary structure of these systems is a one

dimensional chain1. Additionally, polypeptides and polymers comprise only a limited number of monomer

types, many of which contain –H2C–CH2– motives. The, nearly unpolar central C–C bond of these is

a natural choice for a QM/MM boundary. All of this lets the number of covalent QM/MM bonds be

comparatively limited and the bonds easily be described with good precision.

In solids, however, the situation is different. One of the characteristics of solids is the high coordination

of each heavy atom by other heavy atoms. In the systems of interest for this work, i.e. SiO2and Al2O3,

(aside from interfaces or defects) each atom has at least as many nearest neighbors, as it has natural

valences. Additionally, every single bond in these systems is polar. This has important consequences for

the coupling between QM and MM zones within a single calculation.

While in biochemical systems, most problems can be avoided by appropriate choice of the QM cluster,

i.e. problematic QM–MM links can be avoided by extending the cluster to include them in the QM zone

1with some side chains

4.2. EMBEDDING APPROACHES FOR SOLID STATE MATERIALS 33

and place the boundary at a more favorable site, such possibilities are very limited in a solid. First of all,

polar bonds cannot be avoided in a polar material. Secondly, extending the QM cluster also increases

the number of QM–MM links and is therefore likely to cause more problems2than it solves. And lastly,

a QM cluster containing an interface in a polar solid is usually charged, if it does not have the bulk

stoichiometry and symmetry and no manipulation of the charges is applied. In a bulk polar solid, where

the polarity is caused by charge transfer between atoms of different elements, any cluster of the correct

stoichiometry and symmetry should be neutral. As the trials show, this can lead to serious distortion of

the electronic structure of the QM cluster compared to the fully QM treated system.

4.2 Embedding approaches for solid state materials

4.2.1 Representation of the external field

The first issue to be addressed in a QM/MM coupling scheme is, how to implement the coupling terms

in the QM zone Hamiltonian. The quantum mechanical issues in this, have already been discussed in

chapter 2.8, in the following I describe the tested approaches in their practical implementation.

Simple point charge embedding

The simplest scheme to implement is the representation of the electric field of the MM zone as a distri-

bution of point charges corresponding to the MM partial charges at the respective atomic positions. To

take into account the effect of the external point charges on the QM calculation, a coupling term is added

to the total Hamiltonian, as described in equations 2.67,2.91 (cf. chapter 2.8 on page 20). The point

charges can either be obtained from the force-field or, as it is done here, from a QM reference calculation

on a small periodic unit cell. In this work, the DFTB Mulliken populations obtained from these reference

calculation form the basis of all charge distributions.

One of the problems that can arise from using point charges, are unphysical populations or even conver-

gence failure, in case QML atoms are too close (<∼0.4˚

A) to the point charge of their respective MMH

atom. This can lead to severe problems in the calculations, as discussed in chapter 2.8.1 on page 21.

Uniformly Gaussian blurred point charges

To eliminate such problems at the QM/MM frontier, the charge of the MMH atom is often removed

and distributed among the remaining atoms of its neutrality group, e.g. a single amino acid. This

approach [61], which is successful for biochemical applications, it cannot be transferred to the solid state,

as no general definition of a neutrality group is possible there3.

In a more elaborate approach, the point charges are replaced by Gaussian-shaped charge distributions

of variance σ. This effectively reduces the close-range intensity of the electric field, while maintaining

the overall charge of the full charge distribution, introducing only a small perturbation of the long-range

field (cf. equations 2.69,2.92 in section 2.8.1). Physically, replacing the point charge by a spatially

extended charge distribution is justified, since atoms are not points but formed by the combination of

nuclear charge and the charge distribution of its electrons. The choice of a Gaussian charge distribution

is arbitrary and motivated by mathematical convenience. σshould be chosen so that the embedded QM

cluster reproduces the properties of a fully QM calculation as closely as possible.

For a heteronuclear material element-specific, i.e. non-uniform, Gaussian blurring of the atomic partial

charges would be possible. I choose not to do this, as it would increase the number of parameters in the

embedding too much.

2In the treatment of the QM–MM boundary. Of course, the primary criterion for the cluster size is determined by the

extent of the effects which have to be treated quantum mechanically.

3Charge neutrality groups are often defined for ionic crystals, but in the case of disordered systems, such as γ-alumina,

or interfaces between solid phases, the description of which is the goal of this work, these approaches cannot be applied.

34 CHAPTER 4. DFTB QM/MM EMBEDDING

4.2.2 Treatment of the dangling QM–MM bonds

The second issue of QM/MM coupling is the treatment of the dangling bonds which occur on the QM

cluster surface. The following approaches can be combined in any manner with the electrostatic coupling

methods described above.

Simple link atoms embedding (SLA)

In the simplest and most popular approach, the dangling QM–MM bonds are completely ignored on the

MM side and the QMH–MMH bond distance is constrained to a fixed value. The bond on the QM side

is then saturated by adding a hydrogen atom at its equilibrium bond distance to to the QMH.

Link-atom distance scaling

However, additional gap states are often induced by the QMH–QML bonds. If these states are unoccupied,

they can be shifted out of the gap, by increasing the overlap between QMH and QML orbitals involved in

the bond. This is achieved by shortening the QMH–QML distance below the equilibrium bond distance.

Additionally, sometimes erroneously high charge transfer into the QMH–QML bonding state(s) can be

alleviated by modulating the bond distance. When scaling the QMH-QML distance to higher values,

it should be kept in mind, that the QML should not come too close to the MMH point charge, since

convergence problems and other difficulties may arise (cf. section 4.2.1).

The link atom distance scaling was implemented by multiplying the QMH–QML equilibrium distance,

determined from their tabulated covalent radii, with a link-atom distance factor (LADF).

Charge constraints

In addition to placing hydrogen saturators at the QM/MM boundary, it is also possible to constrain the

Mulliken populations of the QMH atoms to the charges they (or their symmetry equivalents in a smaller

cell) show in a full DFTB single-point reference calculation. The charge constraints are implemented as

described in section 2.7.3.

4.2.3 Neutralization of the spurious cluster charge

The last, important issue is the neutralization of the artificial cluster charge, to avoid erroneous filling or

emptying of the band edges. The schemes developed here aim at removing this artificial total charge of a

cluster that arises as a result of cutting the cluster from the polar bulk (cf. figure 4.3). Desired charges,

e.g. when examining the interaction of a charged amino acid with a surface, must not be neutralized by

any of these methods. The cluster neutralization is independent on the electrostatic embedding, yet the

approaches discussed here are tailored to link-atom termination of the dangling bonds.

Homogeneous Charge Subtraction (HCS) embedding scheme

As discussed in section 4.1, in polar systems, only special QM clusters are neutral. To cure this problem I

artificially enforce the neutrality of the QM and MM zones by shifting part of the MMH’s partial charge

4.2. EMBEDDING APPROACHES FOR SOLID STATE MATERIALS 35

QMH QML MMH

- -

-++

QMZ

QMH QML MMH

- -

QMZ

Figure 4.3: Illustration of the QM zone neutralization. In a neutral polar material, each neighbor

of an atom compensates a portion of its partial charge. (In most cases, this is formally achieved

by the nearest neighbors only, and I will adhere to this convention.) The QM/MM boundary leads

to the effect that the partial charges of the QMH atoms are not fully compensated within the QM

zone(left). Overall neutrality of the QM cluster is restored by moving the missing compensating

charge from the MMH atom onto the QML atom and thus into the QM zone (right).

to its associated QM link atom4( cf. figure 4.3):

Q′

i=Qi−Qtot

NQML

∀i∈QML, j ∈MMH, i terminates bond to j(4.1)

Q′

j=Qj+Qtot

NQML

(4.2)

where Q′

i,Q′

jare the modified QML and MMH charges; Qi,Qjare the respective original charges,

Qtot

QM is the total QM cluster charge and NQML is the number of QM–MM links (≡number of QML

atoms, assuming that no higher-order or aromatic bonds cross the QM/MM boundary). After this re-

arrangement, both the QM cluster and MM zone are neutral, provided the whole system was neutral

before. This eliminates the imperfect filling of valence or conduction band in the charged clusters.

However, a spurious dipole moments at the QM/MM border is introduced:

derr =Qtot

NQML

(~rMMHA −~rQML).(4.3)

Since the MMH–QML distance is often very small (<0.1 ˚

A) and the dipoles are normal to the QM/MM

boundary, their influence is comparatively small. The error caused by these dipoles is smaller than the

error resulting from the charged-cluster problem, as the results presented in section 4.3.3 show.

It should be noted, that this scheme depends on the QM cluster to be constructed in a way which ensures

that all QMH-s and all MMH-s are of the same type, respectively.

Bond Charge Transfer Compensation (BCTC) embedding

The HCS scheme introduced above relies on the ability to construct a cluster where only one type of bonds

cross the QM/MM boundary. This imparts severe limitations not only on the choice of QM clusters, but

4The charges of the QM zone atoms are not necessarily constrained during the embedded calculations, in fact the atomic

charges of the QM zone core atoms must be free to be able to describe chemical reactions. Here the atomic partial charges

are only used to determine the cluster total charge by summing them up.

36 CHAPTER 4. DFTB QM/MM EMBEDDING

also on the choice of examined systems. In compound materials, where a portion of the cluster surface

corresponds to a non-polar crystal surface, the need for only one certain bond type crossing the border,

will lead to a very rough cluster surface and an increased number of interrupted bonds which have to be

saturated (cf. figure 4.2 on page 32). Additionally, this situation will tend to enforce increasing sizes,

with unfavorable consequences on the computational cost. Furthermore, solid interfaces like SiC/SiO2

cannot be treated in the QM zone, as the two phases make it impossible to limit the QM–MM links to

just one bond type.

To overcome this limitation, I extended the HCS scheme to a more general formulation: The assumption

of a homogeneous charge transfer for a specific combination of bond partners is retained. However, a

matrix of per bond-type charge transfer parameters dqa,b is constructed. For a system with Nelements,

the N×Nmatrix Qhas the form:

Q=





dq1,1dq1,2··· dq1,N

dq2,1dq2,2··· dq2,N

.....

dqN,1dqN,2··· dqN,N







=





0dq1,2··· dq1,N

−dq1,20··· dq2,N

.....

−dq1,N −dq2,N ··· 0







(4.4)

where the physical assumption that the charge transfer is caused by differences between the chemical

elements leads to the zero diagonal matrix elements and the preservation of charges leads to the dqa,b =

−dqb,a antisymmetry. This approach is inspired by the charge increments formalism used in many force

fields to set the atomic partial charges. To determine the dqa,b for a specific system,a set of linear

equations, expressing the Mulliken charge Qiof each atom of element ain terms of charge transfers with

its bond partners was constructed:

Qi=

b=1 na,b ·dqa,b a < b

−na,b ·dqa,b a≥b,(4.5)

where na,b is the number of bond partners of element b. In a model with Matoms, this procedure

yields a set of Mequations for the upper triangle of Qfrom ( equation 4.4), including the diagonal

elements. The diagonal elements are included for two reasons: first, it implementation is easier and

secondly, the deviations of the diagonal elements from zero allow to assess the quality of the solution

found for equations 4.5.

The equations are constructed from the results on the fully QM treated system. To solve equation. 4.5, I

use the linear least squares method, since it is well known to be an efficient and precise method to solve

an overdetermined system of linear equations[128]. It is robust against degeneracy introduced by periodic

expansion of the underlying molecular model and scales favorably with the number of equations.

The modification of the formal atomic charges of MMH and QML is performed equivalently to the HCS

method. Equations 4.1,4.2 become:

Q′

i=Qi−dqel(i),el(j)−qhom ∀i∈QML, j ∈MMH, i terminates bond to j(4.6)

Q′

j=Qj+dqel(i),el(j)+qhom (4.7)

where el(i) is the element of the QMH atom saturated by the QML atom i, el(j) is the element of

MMH jand qhom is a small homogeneous compensating charge. It becomes clear that HCS is a special

case of the BCTC approach where dqa,b is approximated by Qtot

NQML . The compensating charge qhom is

necessary since, unlike with HCS, the sum of dqa,b-s applied to the cluster will only neutralize it within

the variance of the dqa,b.qhom compensates for this residual charge error by an additional homogeneous

charge transfer between QM and MM zone, calculated as:

qhom =1

NQML 

Qtot

QM −

QML

MMHA

dqel(i),el(j)

.(4.8)

4.3. VALIDATION OF THE QM/MM EMBEDDING IN DIFFERENT MATERIALS 37

4.3 Validation of the QM/MM embedding in different materials

4.3.1 Evaluation Criteria

Electronic properties

In the framework of SCC-DFTB the atomic Mulliken population of the valence shell is varied self-

consistently. Therefore, the converged Mulliken charge of each atom should be the principal property

for the evaluation of a QM/MM embedding scheme. From this follows, that the first evaluation criterion

should be the sum of Mulliken charge differences on the atoms of the QM cluster between full SCC-DFTB

and the QM/MM embedding scheme. This indicates an erroneous charge transfer between QMH and

QML atoms:

∆Qtot =X

i∈QM

QQM/MM

i−QDFTB

i(4.9)

where QM is the set of atoms in the QM cluster. It should, however, be considered, that the bond-

termination might induce a tolerable amount of disturbance on the QM border atoms.

Besides erroneous charge transfer to or from the QML-s, charge oscillations within the cluster can be

induced by the QM zone termination. Additionally, the choice of a suitable QM cluster, regarding its

size, stoichiometry and the minimization of cross-zone bonds, can influence on the achievable quality of

the QM/MM coupling. To detect intra-cluster charge oscillations and to facilitate comparison at different

cluster sizes, the root of mean squares of QM atom charge deviations will be the second criterion:

RMSQ = v

i∈QM QQM/MM

i−QDFTB

i2

NQM (4.10)

The second, equally important property of the embedded QM cluster, is the reproduction of the density

of states (DOS) and the Fermi energy Ef. Ideally, the QM cluster should exhibit no additional gap states,

although some very shallow ones might be acceptable, as long as their nature is similar to those of the

band edges of the reference. Any other gap states would lead to severe distortion of the chemical behavior

of the cluster, rendering the examination of chemical reactions impossible. Deviations in the Fermi level

can be accounted for in reaction energetics but should, obviously, be avoided. However, comparison of the

Fermi level between periodic and cluster calculations is difficult, since the periodic boundary conditions

may introduce a global shift of the band structure, which does not influence energy differences between

structures. Since my primary focus rests on molecules and clusters, I only regard the DOS at the γ-point,

calculated in the conventional way by convolution of the eigenspectrum with a Lorentz function (which

is simplified by the discrete eigenspectrum obtained from DFT(B)):

DOS(E) =

i=1

s2+ (E−ǫi)2(4.11)

where Nis the number of effective single-particle eigenstates, ǫithe i-th eigenvalue and sis the (arbitrarily

chosen) half-width at half-maximum (HWHM) of the Lorentzians. I then examine the DOS of each cluster

to check for spurious states within the gap of the periodic reference system. I write Efwith reference to

the highest occupied molecular orbital HOMO energy of the respective reference system.

I check the DOS replication by two criteria: first I determine the overall shift of the DOS by finding

the maximum in the cross correlation between reference DOS and the DOS of the examined embedded

cluster:

R(∆E) = ZE

DOS(E−∆E)·DOSref(E)dE. (4.12)

38 CHAPTER 4. DFTB QM/MM EMBEDDING

As long as the densities of states are of reasonably similar shape, the cross-correlation is a good means to

find the overall shift of the band-structure. Secondly, I integrate the number of states within the shifted

energy range, corresponding to the gap of the reference DOS

G = ZEref

C+∆E

Eref

V+∆E

DOS(E)dE, (4.13)

where Eref

Vand Eref

Ccorrespond to the valence and conduction band edges of the reference DOS.

-10

-5

DOS (arb. units)

cross correlation

Energy [eV]

reference

cluster correlation

(a) DOS cross correlation

States in shifted gap

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

(b) states in shifted band gap

Figure 4.4: The cross correlation between a reference density of states and the DOS of an embed-

ded cluster (a) and sum of states within that shifted range (b). Here, the full model of the 2x1

reconstructed Si surface and a cluster of 29 Si atoms and 28 saturator atoms, embedded in 5.0 a.u.

Gaussian-broadened charges are shown.

Figure 4.4(a) illustrates, how the band shift can be determined using the cross-correlation function. The

cross correlation (blue line) has a distinct absolute maximum at +0.16 eV. This means that the cluster

DOS must be shifted up by 0.16 eV to coincide best with the reference DOS. The two side maxima

around ±6.5 eV correspond to the superposition of the reference CB with the cluster VB and vice-versa.

Figure 4.4(b) illustrates, how the number of gap states varies with the parameters of embedding. The

Gaussian blur width has very little effect on the number of gap states in contrast to the link-atom distance

scale. Here, downscaling the link-atom distance obviously introduces a large number of new states within

the gap (cf. section 4.2.2).

Due to the product in the cross-correlation function (cf. equation 4.12), higher features of the DOS (i.e.

the bulk bands) have a greater weight in determining the shift than, e.g. single gap states. Care must be

taken in examining the cross correlation, if one of the densities of states has a high peak not present in

the other, in this case the cross-correlation maximum can coincide with the superposition of this singular

peak with one of the bulk bands.

Geometry reproduction

To check the quality of geometry reproduction, the geometry of the embedded cluster is optimized,

keeping the QML and QMH atoms fixed and employing the conjugate gradients method. The root of

mean squares displacement (RMSD) of the QM zone atoms, excluding the QML atoms is then calculated

as the criterion for geometry reproduction:

RMSD = sPNQM /∈QML

i=1 |r′

i−ri|2

NQM/∈QML

4.3. VALIDATION OF THE QM/MM EMBEDDING IN DIFFERENT MATERIALS 39

where r′are the optimized and rthe original atomic coordinates.

In addition to the RMSD, I also visually check the optimized geometries for qualitative changes, e.g.

bond-breaking or -formation.

Since the geometry optimization is computationally highly involved, I only perform this check for the

embedding parameters chosen from the electronic evaluation results.

Reaction Energies

Since reactions are the focus of this work, the ability of an embedding scheme to reproduce reaction

energetics is the ultimate requirement: The enthalpy of a simple interface reaction is calculated, if an

embedding scheme appears promising for a specific substrate, and compared to the results from a fully

SCC-DFTB treated, periodic interface model.

4.3.2 The fully QM treated surface

The gain in computational speed, achieved by employing a QM/MM coupling scheme can only be utilized

to gain efficiency, when the error imparted by this simplification remains within the bounds of tolerance

imposed by the application. To asses the error introduced by partitioning a surface or interface model into

parts which are treated classically or quantum mechanically, respectively, I start from the assumption, that

any deviation from a fully quantum mechanical treatment, exhibited by a QM/MM coupled simulation,

is erroneous5. The first step towards evaluating different QM/MM coupling schemes therefore lies in

the examination of the fully quantum mechanically described surface to obtain references for the criteria

introduced in section 4.3.1.

Computational details

All reference and test calculations are performed using SCC-DFTB, using the pbc-1 set of Slater-Koster

files for Si-containing systems and the alo-1 set for the two alumina examples (cf. Chapter 5). I iterate

the SCC calculation, until the total charge difference between two subsequent cycles falls below 10−5a.u.

Unless stated otherwise, I perform geometry optimization using a CG algorithm until the maximum

force component is below 10−4a.u. To avoid possible convergence problems due to nearly degenerate

states around the gap, I apply a Fermi-Dirac distribution function with a temperature of 300 K to the

occupation of the eigenstates. For the calculation of the DOS, I employ a half-width of 0.1 eV.

Classes of materials

Based on the intended application of the QM/MM coupling scheme, i.e. reactions at the γ-alumina sur-

face, which has rather complicated properties, I will test the different embedding approaches on materials

from different classes. I will start with a simple, ideally coordinated, homonuclear material – silicon.

From there I shall advance to silicon oxides, which introduce polar, heteronuclear bonding but remain

ideally coordinated. The third class of materials will be different phases of alumina, i.e. α-Al2O3which is

over-coordinated6but still crystalline, and native or γ-Al2O3, which is over-coordinated and amorphous.

5The underlying assumption, that the fully periodic quantum mechanical simulation is free of errors, is somewhat

problematic: It beats the idea that QM/MM coupling should make larger and hence more realistic models accessible.

However, within the set of test-systems presented here, which are mostly undisturbed surfaces and interfaces, it can be

assumed to hold.

6More precisely: each atom has more direct neighbors than formal valences.

40 CHAPTER 4. DFTB QM/MM EMBEDDING

Example surfaces

Since this work aims at surface and interface studies, all test- and reference systems chosen for the

evaluation of QM/MM embedding schemes, contain some kind of surface. I perform the SCC-DFTB

reference calculation for a small periodic model of each test system. In the following paragraphs, I

describe the characteristics of of each test system as found in the reference calculations.

Si (Si-2x1)

(a)

-7

-6

-5

-4

-3

-2

-1

DOS (arb. units)

Energy [eV]

surface DOS

Evref Efref

Ecref

(b)

-7

-6

-5

-4

-3

-2

-1

DOS (arb. units)

Energy [eV]

Γ DOS

MP DOS EvΓ

EfΓEcΓ

Figure 4.5: SCC-DFTB total densities of states around the band gaps of the Si (2x1) (0001) surface

(a) and bulk Si (b). For the bulk, the Γ-point DOS and a DOS from an 8x8x8 MP calculation are

given, for comparison.

(a) (b)

-0.25

-0.2

-0.15

-0.1

-0.05

0.05

0.1

0.15

0.2

atomic Mulliken charge in e

Figure 4.6: Structure and charge distribution (a) and atomic charge histogram (b) of the 2x1 Si(001)

surface, determined using SCC-DFTB.

4.3. VALIDATION OF THE QM/MM EMBEDDING IN DIFFERENT MATERIALS 41

Figure 4.7: Isosurface (in arbitrary units) of

the localized density of states in a 8x8 ˚

Alat-

eral cutout of the Si (001) 2x1 reconstructed

surface. Energy range -5.26 – -3.96 eV.

The first reference surface is the 2x1 reconstructed

Si (001) surface. I model the surface as a 2D

slab showing the 2x1 reconstruction on one side

and hydrogen termination of the dangling bonds

at the opposite face. To accommodate a QM zone

later on, the model is periodically expanded to in-

clude 2x2 reconstruction units, as it is shown in

figure 4.6 (a)). The surface model is character-

ized by a non-polar bulk (a small perturbation in

the form of a ∼+0.1 e charge of the lowermost Si

layer is introduced by the H termination of dan-

gling bonds) and charges of ∼ ±0.2 e on the recon-

structed atoms. In this homonuclear system (the

saturator hydrogens are not part of the system be-

ing modeled), the calculation of BCT coefficients

makes no sense, although the system shows a dis-

tinct charge transfer between neighboring silicon

atoms in the reconstruction dimers. The BCTC

approach, however, only uses the core charge to

identify an atom type, so that it cannot describe

this particular situation. Since the bulk is apolar,

it makes much more sense to just avoid cutting

through a reconstructed dimer than to further ex-

tend the cluster neutralization formalism.

Figure 4.5(b) shows the DOS of bulk Si at the

gamma-point and using a 8x8x8 Monkhorst-Pack k-point sampling [11]. The DOS shows the typical

behavior of SCC-DFTB, which gives a good reproduction of the direct band gap, but (due to the omission

of d-orbitals) fails to reproduce the indirect gap of Si. The SCC-DFTB direct HOMO-LUMO gap of the

surface model is ∼0.89 eV. The DOS plot in Fig 4.5 (a) shows a narrowing of the HOMO-LUMO gap

compared to bulk Si, which is caused by a peak of unoccupied states around mid-gap. In comparison to

the bulk data, (taking into account an overall shift of the surface eigenspectrum towards lower energies),

it appears that the VB–CB gap actually widens from ∼1.4 eV to nearly 2 eV, while a peak of unoccupied

states is introduced by the surface around mid-gap. Figure 4.7 shows the localized density of states,

defined as the sum of the density contributions of the effective single particle eigenstates within a certain

energy range:

LDOS = Pi, Elo<ǫi<Ehi |Ψ|2

Ehi −Elo

(4.14)

where Elo and Ehi represent the high and low boundaries of the energy interval in question. It shows

clearly, that the eigenstates corresponding to the DOS peak just below -4 eV are strongly localized to the

reconstructed surface dimers. These results are in good agreement with earlier theoretical work, using

DFT to examine Si surface reconstructions (cf. e.g. [129]).

42 CHAPTER 4. DFTB QM/MM EMBEDDING

(a)

-15

-10

-5

DOS (arb. units)

Energy [eV]

DOS Evref Efref Ecref

(b)

-6

-4

-2

DOS (arb. units)

Energy [eV]

DOS Evref Efref Ecref

Figure 4.8: Densities of states around the band gaps of the α-SiO2(a) and SSZ60 (b) reference

structures.

α-Quartz

(a) (b)

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

1.2

atomic Mulliken charge in e

Si H O

Figure 4.9: Structure of the 1x1 reconstructed α-SiO2surface model (a) and atomic charge his-

tograms (b) from the SCC-DFTB reference calculation

The next evaluation system is one of the (0001) surfaces of α-quartz. It shows a 1x1 surface reconstruction

by re-arrangement of siloxane bridges, which eliminates all dangling bonds. As in the Si surface model,

one face of the 2D slab model is reconstructed, while on the other face dangling bonds are saturated by

H atoms. The silica class materials are characterized by stronger polarity of the bulk, while the recon-

structed surface shows no significant difference in polarity, compared to the underlying bulk (figure 4.9).

The heteroatomic BCT coefficients in this model, listed in table 4.1, conform to the expectations from

4.3. VALIDATION OF THE QM/MM EMBEDDING IN DIFFERENT MATERIALS 43

Elements H–H H–O H–Si O–O O–Si Si–Si

dqa,b 0.000 0.000 -0.159 0.000 -0.262 0.000

Table 4.1: Bond charge transfer coefficients of the α-quartz surface model. The error is σ=0.019.

experience. The homonuclear coefficients are <10−8and their variance is reasonable, showing that the

coefficients are of suitable quality.

The reference direct band gap of this surface is ∼7.94 eV. As can be seen in Fig 4.8(a), the surface does

not cause any deep gap states in this case, however, small DOS peaks close to the band edges appear,

leading to a narrowing of the gap of ∼1 eV.

SSZ60-zeolite

(a) (b)

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

1.2

atomic Mulliken charge in e

Si O

Figure 4.10: Geometry and Mulliken charge distribution (a) and atomic charge histogram (b) of

SSZ60 Zeolite, calculated using SCC-DFTB.

The second silica class test system is a model bulk SSZ60 zeolite. It stands out from other references

examined here, by not having an outer surface, but inner surfaces on the walls of its pores. In its

electronic properties, it differs only slightly from the α-quartz phase, it is included to be able to assess

the transferability of coupling schemes within a class of materials. Figure 4.10 shows a 2x2x2 periodically

expanded version of the basic model for clarity. As can be seen from the charge histograms, the partial

charges are very similar to the α-quartz phase. The O–Si BCT coefficient is equal to that in the α-quartz

(cf. tables 4.1,4.2) with a very reasonable error.

Elements O–O O–Si Si–Si

dqa,b 0.000 -0.261 0.000

Table 4.2: Bond charge transfer coefficients of the bulk SSZ60 model. The error is σ=0.015.

44 CHAPTER 4. DFTB QM/MM EMBEDDING

The HOMO-LUMO gap found in the reference calculation is ∼9.33 eV. In this system, the DOS, plotted

in Fig 4.8(b), does not show any striking features, which is to be expected, as a perfect bulk system is

examined here.

α-Al2O3(corundum)

(a)

-10

-5

DOS (arb. units)

Energy [eV]

DOS

bulk Evref

Efref Ecref

(b)

-10

-5

DOS (arb. units)

Energy [eV]

DOS Evref Efref Ecref

Figure 4.11: SCC-DFTB densities of states around the band gaps of the α- (a) and γ- (b) Al2O3

reference structures, calculated using SCC-DFTB.

(a) (b)

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

atomic Mulliken charge in e

Al H O

Figure 4.12: Geometry and Mulliken charges (a) and atomic mulliken charge histograms (b) of the

α-Al2O3(0001) surface from the SCC-DFTB reference calculation.

4.3. VALIDATION OF THE QM/MM EMBEDDING IN DIFFERENT MATERIALS 45

The last class of test materials are different phases of alumina. The first example is the highly ordered

α-Al2O3phase, also known as corundum or sapphire7. It has high symmetry and evenly distributed

polarity of the bonds. The high degree of coordination (all atoms are over-coordinated) hints, that the

bonding is much less covalent in nature than in the preceding examples. Yet, the polarity is comparable

to silica and by no means ionic. As shown in table 4.3, the O–Al bond charge transfer is much smaller

Elements H–H H–O H–Al O–O O–Al Al–Al

dqa,b 0.000 0.278 0.000 0.000 -0.100 0.000

Table 4.3: Bond charge transfer coefficients of the bulk α-Al2O2model. The error is σ=0.049.

than the O–Si value in the silica systems, with a large relative error. This hints that the BCT model

does not describe the bonding situation perfectly.

I employ a 2D slab model, with the surface in (0001) direction. The Al terminated surface is 2x2

reconstructed, while on the O terminated side, the dangling bonds are terminated using hydrogen atoms.

Figure 4.12 shows, that the distribution of partial charges is uniform throughout the bulk, with only a

small distortion introduced by the hydrogen termination of the (000¯

1) surface. The surface Al-Atoms

however differ by having much smaller partial charges than in the bulk. Curiously, this is the opposite

behavior from the Si (001) surface (cf. figure 4.6). The comparison of the surface model DOS and bulk

DOS in figure 4.11shows that the surface leads to a downwards shift of the whole eigenspectrum by about

2 eV and to the introduction of a whole spectrum of, mostly empty, gap states.

Hydroxylated γ-Al2O3

(a) (b)

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

atomic Mulliken charge in e

Al H O

Figure 4.13: Results of the SCC-DFTB reference calculation on the hydroxylated γ-Al2O3surface

model used for the adhesion studies in chapter 5(a): structure and charge distribution, (b): atomic

charge histograms

7In geology, sapphire denotes all natural corundum precious stones, except the red varieties, which are generally known

as ruby. The colors of these materials are caused by (usually metal ion) impurities, e.g. Cr3+ ions in rubies.

46 CHAPTER 4. DFTB QM/MM EMBEDDING

The last, and definitively most complicated evaluation system is a model of native alumina. Its most

striking difference to corundum is its pseudo-amorphous structure, a more in-depth discussion of which

is given in chapter 58. I employ the surface model described in chapter 5as a reference for the desired

application scenario. The results for the stoichiometric unit cell are illustrated in figure 4.13. The

difficulties in the description of this material, in addition to those described above for corundum, stem

from the amorphous structure: Al appears in tetrahedral and octahedral configurations imposed by

the spinel structure and the defects (i.e. Al vacancies) lead to fluctuating coordinations for each atom

leading to distinctly fluctuating atomic partial charges, as can be clearly seen when comparing the atomic

Mulliken charge histograms (cf. figures 4.6,4.9,4.10,4.12, and 4.13). The average BCT coefficients (cf.

Elements H–H H–O H–Al O–O O–Al Al–Al

dqa,b 0.000 0.285 0.000 0.000 -0.119 0.000

Table 4.4: Bond charge transfer coefficients of the γ-Al2O3surface model. The error of σ=0.14733

hints, that the definition of dqa,b has difficulties in describing this system.

table 4.4) match those found for corundum, but with a relative error >1 (for O–Al), considerable doubt

is cast on the BCT model. However it should be considered, that the –OH groups attached to the surface

Al atoms will disturb the BCT coefficient calculation.

The HOMO–LUMO gap of this surface model is only ∼0.019 eV, but as can be seen in figure 4.11(b),

this is not the band gap of the bulk material but caused by surface states.

4.3.3 Evaluation of coupling schemes

When chemical reactions are examined, the major contributions9originate from changes in the chemical

bonding, while changes within the more distant environment are small (cf. figure 4.1). Therefore, the first

step in testing a QM/MM coupling scheme is to ensure that the description of the embedded QM zone

meets the requirements of the intended application. This means that the representation of the QM zone

should match fully QM calculated references as closely as possible with respect to the criteria defined

in section 4.3.1. The internal dynamics of the MM zone should be neglected during these tests, since a

molecular-mechanical description of the MM zone may inherently introduce systematic deviations from

the QM description, e.g. by slightly shifting the positions of the external charges. The effects of such

changes could not be easily distinguished from errors introduced by the transition from a fully QM to an

embedded modeling. Additionally, a force-field description of the MM zone adds to the overall complexity

of the system and therefore introduces further potential errors. I therefore use the atomic equilibrium

positions and partial charges as a fixed model of the MM part in my evaluation of the different embedding

schemes.

Since the QM calculations are limited to small periodic supercells, or even the unit cell, I periodically

expand the distribution of “MM atoms” in order to obtain a sufficiently large MM zone10 (cf. figure 4.14).

In addition to the computational details described in Section 4.3.2, I keep the link-atoms and QM border

atoms fixed during geometry optimization, unless stated otherwise. I keep the whole QM cluster geometry

fixed while searching for embedding parameters, in order to limit the computational cost. After that, I

select a set of embedding parameters and relax the geometry using these parameters.

The complete set of plotted evaluation results, can be found in appendix A.

8Suffice to say here, that the native Al2O3phase, labeled γ-Al2O3, is best described as a defective spinel structure,

with the oxygen part forming a perfect fcc lattice and the Al part on tetrahedral and octahedral positions. To ensure

stoichiometry, Al atoms must be randomly removed from the perfect spinel structure. Additionally, the native oxide surface

is generally hydroxylated.

9energetic and otherwise

10The original, periodic supercell is split into the QM- and MM zones. The external charges are then the atomic Mulliken

charges of the MM zone. Additionally, one or more shells of charge distributions of the full supercell (containing the charges

of QM- and MM zone atoms) may be added around the original cell.

4.3. VALIDATION OF THE QM/MM EMBEDDING IN DIFFERENT MATERIALS 47

Figure 4.14: Illustration of the embedding of a QM zone (green) within the charge distribution of

the MM zone, which comprises the MM atoms of the original supercell (red) and further replicas of

the whole original charge distribution (blue).

Simple link atoms (SLA)

The simple link atoms approach (cf. section 4.2.2) with point charges performs well for the (2x1) recon-

structed Si surface. Due to the nearly vanishing point charges, the cluster charge is very small (0.081 e),

so that only small errors are expected to arise from this; the Fermi level will be shifted to the VB edge

but the difference in occupation of <1 e will not influence reaction energetics significantly.

(a) (b)

-0.3

-0.25

-0.2

-0.15

-0.1

-0.05

0.05

0.1

0.15

0.2

0.25

atomic Mulliken charge in e

ref. QM/MM

Figure 4.15: Charge distribution (a) and atomic charge histogram (b) of the Si13 2x1 Si(001) surface

using point-charge modeling of the MM zone.

48 CHAPTER 4. DFTB QM/MM EMBEDDING

Using a QM core of 29+28 (29 QM + 28 QML) atoms in a 3x3 expanded surface supercell to accommodate

the larger cluster, the charge distribution of the inner QM atoms matches that of the reference quite

closely, as can be seen from figures 4.15(a) and 4.6. The calculated HOMO-LUMO gap of 0.97 eV is

slightly larger than the reference of 0.84 eV, as would be expected as a result of the electronic confinement.

1.838

1.84

1.842

1.844

1.846

1.848

1.85

0.177

0.178

0.179

0.18

0.181

0.182

0.183

0.184

0.185

dQtot [e]

RMSQ [e]

gaussian blur width [a.u.]

dQtot

RMSQ

Figure 4.16: QMH atom total and RMS charge errors of the 29+28 atom cluster embedded into

Gaussian blurred charges, dependent on the Gaussian blur width. Crosses mark calculated values,

lines serve to guide the eye.

From figure 4.15(b), it can be seen that simple hydrogen termination causes some charge distortions on

the QMH atoms (cf. the histogram peaks around +0.02 e in the reference and -0.1 e in the embedded

calculation). To investigate chemical reactions at one of the reconstructed dimers, a larger QM cluster

which contains at least three of these dimers would be necessary, to isolate the central dimer from

termination-induced polarization effects by at least one shell of non-border neighbor atoms. Gaussian

broadening has a detectable effect on the QM atom charge deviations with distinct minima, however,

as can be seen from figure 4.16, this effect is very small. (The maximal difference is∼0.01 e in dQ and

RMSQ.) This is to be expected, since the point charges themselves are very small.

Changing the link-atom distance has no beneficial effect on this embedded cluster, as can be seen in

figure 4.17. Outside an unfavorable region of short QMH–QML distances at small charge broadening,

the RMSQ is virtually independent on the distance scale. The band shift has its minimum for a QML

distance factor of 1.0, independent of the Gaussian broadening. Figure 4.4(b) on page 38 shows a very

similar picture for the number of gap states as for the band shifts.

Overall, it appears that using an SLA cluster combined with a Gaussian broadening of the external

charges by 5.0 Bohr radii and no link-atom distance scaling gives a very good representation of the 2x1

reconstructed Si surface. When simulating chemical reactions, the cluster size must me checked, to ensure

that the QMH atoms are not the nearest neighbors of atoms involved in the reaction.

The picture changes, however, when applying the SLA embedding approach to silicon and aluminum

oxide systems. First of all a considerable number of gap states appear (cf. figure 4.18(a)) in SiO2, when

placing the QML atoms at their equilibrium bond distance to the QMH atoms. Depending on the cluster

size and QMH atom type, the range of useful link-atom distances and Gaussian blur widths can be very

limited(cf. figures A.2–A.6 or figure 4.18 as an example for α-quartz) . E.g. for the Si7O8α-quartz cluster,

the LADF range is ∼0.8–0.9 and gaussian blur widths should be >∼3.75 a.u.(cf. figure 4.18(a)).Within

these constraints, the choice of embedding parameters is then governed by the minimization of the RMSQ.

4.3. VALIDATION OF THE QM/MM EMBEDDING IN DIFFERENT MATERIALS 49

System CC NQM NQML NMM Qtot

clst NSCC RMSQ RMSD ∆Eb[eV] σgb LADF QMH QQMH Figure

Si 2x1 Si29 29 28 4371 +0.08 17 0.18 0.01 -0.16 5.0 1.0 Si free A.1

α-quartz Si7O815 12 885 +3.11 1001∗0.10 0.27N-1.11 4.0 0.85 Si free A.2

Si7O20 27 12 873 -3.15 89 0.20 1.71N+0.74 1.8 0.85 O free A.3

SSZ60 Si3O10 13 8 10739 -2.07 121 0.16 0.41N-1.65 2.8 1.0 O free A.4

Si11O10 21 24 10731 6.26 102 0.11 0.15N-2.80 4.0 0.80 Si free A.5

Si11O34 45 24 10707 -6.21 55 0.28 0.08 -0.09 4.0 0.95 O free A.6

α-Al2O3Al20O30 50 60 1138 -1.63 24 0.18 0.11 -0.67 1.6 1.075 N/A free A.7

40 0.06 0.06 -0.30 1.7 1.175 constr A.8

Al20O51 71 102 1117 -9.71 27 0.21 0.68N-1.35 1.7 1.0 O free A.9

38 0.11 0.46N-2.10 1.0 1.05 constr A.10

Table 4.5: Summary of the test results for SLA embedding. The column CC gives the cluster composition without hydrogen link atoms; NQM,

NQML and NMM give the numbers of QM atoms, saturators and and external charge centers, respectively. Qtot

clst is cluster total charge in e,

NSCC, RMSQ, RMSD and ∆Ebare the number of SCC cycles until convergence, RMS charge deviation of the QM atoms in e, RMS movement

of the core QM atoms and band shift at the selected embedding parameters σgb, the Gaussian blur width in a.u. , and LADF, the link-atom

distance factor. QMH gives the element of the QMH atoms of the cluster, if all QMH are of the same element. ∗SCC did not converge for this

calculation. NThe geometry underwent qualitative reconfiguration during the optimization run.

50 CHAPTER 4. DFTB QM/MM EMBEDDING

(a)

0.16

0.17

0.18

0.19

0.2

0.21

0.22

0.23

0.24

0.25

RMSQ

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

(b)

-1

-0.8

-0.6

-0.4

-0.2

0.2

band shift

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure 4.17: RMS charge deviations of the QM atoms (excluding QML)(a) and overall band struc-

ture shift(b) of 29+28 atom the 2x1 Si(001) surface cluster, dependent on QML distance scale and

gaussian broadening of the external charges.

(a)

States in shifted gap

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

(b)

States in shifted gap

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure 4.18: States in the shifted gaps of the 15+12 atom (a) and 27+12 atom (b) α-quartz clusters

dependent on gaussian blur width and QML distance factor.

Table 4.5 summarizes the evaluation results obtained with the SLA scheme. It turned out that any cluster

of reasonable shape11 cut from polar materials was severely charged, even if it had the correct stoichiom-

etry, like the Al20O30(+H60). Additionally, in all cases, the number of SCC cycles until convergence was

high in the regions of good charge and band-structure reproduction. This is understandable, since in the

bulk, the band edges are generally degenerate at the Γ-point, which makes convergence of charged bulk

structures extremely difficult. More series however, are the large geometry distortions present in all SLA

embedded systems, except Si and the 45-atom SSZ60 cluster, when relaxing the non-host QM atoms.

The large RMSD values for these structures (table 4.5) indicate severe problems. At visual inspection, it

turns out that in all cases with RMSD >0.2 bonds inside the QM cluster were broken.

11defined as a mostly convex surface with only one QM–MM link per MMH atom. The latter cannot be fully achieved

for alumina materials, because of the high atomic coordination numbers.

4.3. VALIDATION OF THE QM/MM EMBEDDING IN DIFFERENT MATERIALS 51

In the Al2O3test cases, the stoichiometric Al20O30 cluster showed reasonable geometry (RMSD 0.11). In

this cluster, the bond breaking occurred in the vicinity of the QMH atoms. This suggests that this problem

could be caused by erroneous charge transfers to or from these atoms. Therefore, test calculations with

the charges of the QMH atoms constrained to their reference values were performed for the two Al2O3

clusters. The results show that, although the charge constraints lead to an overall improvement in all

benchmark values (e.g. RMSD 0.11 →0.06 and 0.68 →0.48 for Al20O30 and Al20O51, respectively), the

qualitative failure in geometry representation in the Al20O51 cluster could not be cured. Additionally,

the charge constraints make the SCC convergence more difficult, resulting in a increase of about 50% in

the needed SCC iterations and therefore the calculation time.

Still, the Fermi level was pinned to the VB or CB edge in each polar cluster, which would lead to severe

errors in the description of chemical reactions. Therefore, embedding schemes which leave a residual

cluster charge are unsuitable for the simulation of polar materials.

Neutralized cluster schemes

To solve the problems arising from the cluster charge introduced by the SLA embedding scheme in polar

materials, HCS and BCTC embedding (see section 4.2.3f.) were employed. The comparison between the

results from SLA calculations in table 4.5 and the neutralized cluster schemes (table 4.6) shows that they

are successful in this. Except for the very small clusters12, the RMSQ is significantly smaller for the HCS

and BCTC schemes, and in each case, the number of SCC cycles until convergence is much smaller.

The biggest improvement, however, is found in the atomic RMSD of the test clusters. Bond reconfigura-

tion only occurred in the case of the oxygen terminated Al20O51 cluster.

A detail of note is the fact that the band offsets are somewhat larger than with SLA embedding. This

is mostly caused by the additional dipole moments introduced by the charge relocation, as discussed in

sections 4.2.3 f.. However, as this external field acts uniformly on all eigenstates, it should only result in

a small global shift of total energies, while not significantly influencing relative energetics, i.e. reaction

energies, -barriers etc.

For the very small quartz and SSZ60 clusters, some rotation of most Si–O–Si bridges around the O–O

axes is observed and causes the slightly above average RMSD values.

The results for Al2O3indicate, that oxygen QMH clusters do not provide adequate results, with or

without charge neutralization (cf. tables 4.5 and 4.6). Increasing the MM zone size does not cure

these problems. An Al-QMH cluster provides reasonable results with charge constraints, however, the

neutralized stoichiometric clusters clearly give the best reproduction of the results from a fully periodic

model. In all cases, charge constraints lead to an improvement of the RMSQ by about 50% at the cost

of an increase in calculation time by >50% (cf. tables 4.5 and 4.6). Embedding parameters in the

γ-alumina clusters are quite close to each other and comparable to those in the stoichiometric Al20O30

α-Al2O3cluster. This observed transferability between clusters as well as different phases of one material

shows, that the embedding approach is stable and not over-parametrized. Therefore it can be assumed

with reasonable conviction that the clusters can also be used to examine chemical reactions at the surface.

Based on the α-alumina results, the hydroxylated γ-Al2O3surface is modeled using neutralized, stoi-

chiometric clusters (i.e. keeping the Al2O3stoichiometry of the alumina phase — the surface hydroxyl

groups distort the overall stoichiometry). The Al24O40H6(+saturators) cluster is centered around the

-OH group called site 1 in the surface adsorption studies (cf. chapter 5). It contains three binary and

three ternary OH groups, formed from the dissociation of three water molecules. The stoichiometry of

the Al2O3part is not ideal — it has one excess oxygen, which cannot be avoided with reasonable effort

in this pseudoamorphous phase. The results show, that the embedding is quite successful in reproducing

a fully periodic model. One surface Al atom, which forms an O-Al-O bridge between two QMH oxygens,

rotates around the O–O axis during relaxation. Since this case is somewhat pathological and the Al is

located at the edge of the QMZ, I extend the QMH geometry constraints to include this atom.

12Which are anyway problematic for the description of chemical processes (cf. section 4.3.3).

52 CHAPTER 4. DFTB QM/MM EMBEDDING

System CC NQM NQML NMM NSCC RMSQ RMSD ∆Eb[eV] σgb LADF QMH QQMH Figure

α-quartz HCS Si7O815 12 885 13 0.15 0.12 0.82 2.2 0.80 Si free A.11

Si7O20 27 12 873 4 0.05 0.08 -2.06 2.0 0.80 O free A.12

SSZ60 HCS Si11O10 21 24 10731 13 0.32 0.10 -0.01 1.6 0.80 Si free A.13

Si11O34 45 24 10707 12 0.04 0.07 -3.77 1.4 0.65 O free A.14

α-Al3O3BCTC Al20O51 71 102 1117 23 0.19 0.33N0.85 1.5 1.0 O free A.15

39 0.22 0.48N7.04 1.0 1.0 constr A.16

10621 22 0.15 0.34N-21.80 3.0 1.0 free -

BCTC Al20O30 50 60 1138 23 0.11 0.07 0.25 1.5 1.2 N/A free A.17

38 0.05 0.07 0.05 1.5 1.2 constr A.18

BCTC Al32O18 50 81 1138 33 0.23 0.10 -1.60 2.5 0.75 Al free A.19

68 0.04 0.08 -1.20 3.0 0.75 constr A.20

γ-Al3O3BCTC Al24O40H670 52 13955 37 0.17 0.14 -2.625 1.5 1.0 N/A free A.21

63 0.03 0.13 -2.85 2.2 1.0 constr A.22

BCTC Al36O56H12 104 67 13921 32 0.11 0.18 -2.425 2.0 1.05 N/A free A.23

57 0.02 0.22 -0.275 2.0 0.80 constr -

69 0.03 0.19 -1.625 2.0 1.0 constr A.24

Table 4.6: Summary of the test results for neutralized cluster embedding. The column CC gives the cluster composition without hydrogen link

atoms; NQM,NQML,NMM and give the numbers of QM atoms, saturators and external charge centers, respectively; NSCC, RMSQ, RMSD and

∆Ebare the number of SCC cycles until convergence, RMS charge deviation of the QM atoms in e, the RMS displacement of the core QM atoms

in ˚

Aand band shift in eV at the selected embedding parameters σgb, the Gaussian blur width, and LADF, the link-atom distance factor in a.u.

QMH gives the element of the QMH atoms of the cluster. NThe geometry underwent qualitative reconfiguration during the optimization run.

4.3. VALIDATION OF THE QM/MM EMBEDDING IN DIFFERENT MATERIALS 53

(a) site 1 (b) both sites

Figure 4.19: The tested γ-Al2O3clusters. The cluster containing both sites is an extension of the

site 1 cluster.

For the calculations in shown in table 4.6, the γ-alumina clusters are surrounded with a distribution

of external charges containing 5x5 original supercells (a total of 13955 external charge centers for the

site 1 cluster). To check for convergence with the size of the external charge distribution, I repeat the

geometry check of the non charge-constrained BCTC embedded Al24O40H6cluster with external charge

distribution sizes ranging from 1x1 to 11x11 original supercells. The results are listed in table 4.713. They

cells ext. charges Etot [H] dEtot [H] RMSQ [e] RMSD [˚

1x1 491 -176.951781 -0.0602 0.14024 0.157802

3x3 4979 -176.913209 -0.0572 0.14050 0.157484

5x5 13955 -176.910081 -0.0567 0.14053 0.156064

7x7 27419 -176.908993 -0.0565 0.14056 0.155394

9x9 45371 -176.908452 -0.0564 0.14058 0.155023

11x11 67811 -176.908130 -0.0563 0.14060 0.154787

Table 4.7: Charge distribution convergence check results for the Al24O40H6γ-alumina cluster.

RMSQ here is after geometry optimization and therefore not directly comparable to table 4.6.

dEtot gives the energy gain by geometry optimization (see footnote 13 on p 53).

show, that the originally chosen 5x5 point charge distribution is reasonably sized, however, a 7x7 charge

distribution would provide somewhat better results. In all examined criteria, the results can be regarded

as converged at the 7x7 external charge distribution. Note that the RMSQ differs only within the same

order of magnitude as the SCC convergence criterion (10−5electrons).

The second γ-alumina cluster examined here is an extended version of the site 1 cluster, constructed so

as to also include the surroundings of site 2 (cf. chapter 5). The results show, that the larger γ-alumina

cluster exhibits a slightly better behavior than the site 1 cluster, in terms of SCC convergence and charge

deviation. This can be expected from the smaller ratio of QM core to QMH and QML atoms, leading to

more bulk-like behavior of the cluster. (Note that the RMSD value entails a sight rotation of the O-Al-O

bridge, which was kept fixed in the site 1 cluster. This accounts for 0.2 in the RMSQ in each case.)

13It should be noted, that the final geometric configuration of the surface –OH groups strongly depends on the initial

geometry and the resulting charge distribution — starting a new calculation from the final geometry of a previous one,

without reading the charges, leads to relaxation to a different (albeit very shallow) local minimum with visible displacement

of the hydroxyl groups. This behavior can also be observed in fully QM calculations and is not inherent to the embedding.

54 CHAPTER 4. DFTB QM/MM EMBEDDING

Differences between SiO2and Al2O3

The embedding evaluations so far show one interesting and unexpected result: Although SiO2and Al2O3

appear to be quite similar materials at first glance: both materials are insulators with direct band gaps of

∼9 eV for SiO2and ∼6 eV for alumina. As can be seen in section 4.3.2, their polarities, in terms of atomic

charges quite comparable and the electronegativity differences between oxygen and the heavy atom are

similar, as are those between the material’s atoms and the hydrogen saturators. All these similarities

lead to the expectation that both materials should behave similarly in a QM/MM framework.

Yet the cluster compositions which work very well for the SiO2systems, i.e. clusters with homogeneous

QMH elements, fail for alumina. Instead, in this system the correct cluster stoichiometry plays a central

role. The reason for this difference is not clear, but the observed bond-breaking in the non-stoichiometric

Al2O3clusters suggests a connection to the structure of the valence band. Figure 4.20 shows the wave

function coefficients of the HOMOs and LUMOs of one unit cell of bulk quartz and sapphire. Comparing

the HOMOs of both materials, one finds that the VB edge of quartz is almost exclusively composed of

oxygen orbitals, as it is common text-book knowledge (note that in the SK parametrization employed, the

Si-d orbitals act as diffuse functions for the oxygens[130]). In comparison, the contribution of Al orbitals

to the VB edge in α-Al2O is much larger than that of Si in quartz (the ratio of oxygen to aluminum

contributions is considerably smaller in α-Al2O). This may explain, why the valence band of alumina is

much more sensitive to deviations from stoichiometry, which disturb the composition of the VBE states,

than that of silica.

0.1

0.2

0.3

0.4

0.5

abs(coeff.)

index

Si O

(a) α-quartz LUMO

0.1

0.2

0.3

0.4

0.5

abs(coeff.)

index

Al O

(b) α-Al2O3LUMO

0.1

0.2

0.3

0.4

0.5

abs(coeff.)

index

Si O

0.1

0.2

0.3

0.4

0.5

abs(coeff.)

index

Al O

(d) α-Al2O3HOMO

Figure 4.20: Atomic orbital coefficients of the HOMO and LUMO of bulk α-quartz and α-Al2O3,

ordered by Atom and, shell nand magnetic quantum number m. Each consecutive line represents

one atom.

4.3. VALIDATION OF THE QM/MM EMBEDDING IN DIFFERENT MATERIALS 55

Reaction energies

After having established suitable embedding parameters for surface and interface clusters in the test

systems, the final test of the developed embedding scheme will be the reproduction of reaction energies.

In order to limit the computational cost, I concentrate on the α-quartz surface and the γ-alumina surface,

examining the reactions:

α−quartz + H2O−→ [α−quartz + H + OH] ,(4.15)

γ−Al2O3−OH + AMEO −→ [γ−Al2O3+ AMEO] + H2O.(4.16)

Reaction 4.15 is the dissociative water adsorption on α-quartz. Reaction 4.16 describes the condensation

of a 1,3-trihydroxiaminopropylsilane at a surface hydroxyl group on native alumina, as described in

chapter 5.

I examine each reaction using BCTC embedding with the cluster which has shown the best results and

the appropriate parameters determined in the previous section. The last parameter that remains to be

studied in a systematic manner, is the size of the external charge distributions.

α−quartz + H2O−→ [α−quartz + H + OH]

Figure 4.21: Si-(OH)-Si bridges formed by adsorption of an H2O molecule on α-quartz. Atoms from

the dissociatively adsorbed water molecule are marked blue.

I use the Si7O20 α-quartz HCS neutralized surface cluster with the embedding parameters from table 4.6.

The geometry of the dissociatively adsorbed water molecule is shown in figure 4.21. To obtain optimal

transferability, I treat the reactants in isolated calculations, i.e. the water is not in vicinity of the surface

before the adsorption. As a reference, I calculate the reaction energy in the full periodic 6x6 surface cell

(cf. figure 4.9).

The reaction is endothermic by 1.02 eV in the periodic reference calculation, while the endothermicity

converges to 0.87 eV in the (non-periodic) embedded calculations with up to 42x42 lateral cells in the

external charge distribution ( figure 4.22). It can be seen, that the reaction energy converges to a value

about 0.2 eV less endothermic than in the fully periodic reference.

56 CHAPTER 4. DFTB QM/MM EMBEDDING

0.7

0.75

0.8

0.85

0.9

0.95

1.05

1000

2000

3000

4000

5000

6000

7000

8000

Er [eV]

surface area [Å2 ]

ext. charges

6x6

18x18 30x30 42x42

periodic ref.

Figure 4.22: Calculated reaction energies for the dissociative H2O adsorption on α-quartz plotted

over the charge distribution surface area. Text labels give the number of lateral surface cells.

γ−Al2O3−OH + AMEO −→ [γ−Al2O3+ AMEO] + H2O

−2.2

−2.15

−2.1

−2.05

−2

−1.95

−1.9

−1.85

−1.8

10000

20000

30000

40000

50000

60000

70000

Er [eV]

surface area [Å2 ]

ext. charges

3x1

9x3 15x5 21x7 27x9 33x11

periodic ref.

Figure 4.23: Calculated reaction energies for the condensation of AMEO on γ-alumina plotted over

the charge distribution surface area. Text labels give the number of lateral surface cells.

(1,3)-trimethoxyaminopropylsilane, also known as APTES or AMEO, plays an important in the alu-

mina/adhesive interface formation examined in chapter 5. Reaction 4.16 is the acid-base condensation of

hydroxylated AMEO with the –OH group at site 1 of the Al36O56H12 cluster containing both adsorption

sites used in chapter 5. Using BCTC embedding with free charges and the parameters listed in table 4.6, I

calculate the condensation energy of AMEO, following the same procedure as for reaction 4.15. As shown

in figure 4.23, the reaction energies calculated with QM/MM converge to -2.10 eV, while the periodic

full-DFTB calculation gives a reaction energy of -1.84 eV (cf. table 5.2), i.e. the QM/MM calculation

finds the adsorbate to be ∼0.25 eV more stable.

General considerations

In all examined reactions, the adsorbate appears ∼0.2–0.3 eV more stable in the non-periodic embedded

calculations than in the periodic full-DFTB calculation. This raises the question, whether this deviation

indicates a systematic error introduced by the QM/MM embedding. It must be noted here, that the

embedded calculations and the full DFTB calculations differ in one very important respect: to avoid the

problems of bond termination, the full DFTB calculations employ periodic boundary conditions. This

4.4. SUMMARY 57

means, that the full DFTB calculations model infinite surfaces with an infinite number of defects, one

per supercell. In contrast, the embedded calculations model one single defect on a finite crystallite.

All examined reactions lead to changes in the dipole moments of the surfaces, especially the added

hydroxyl groups and the AMEO molecule introduce additional dipole moments. The self-interaction of

these dipoles across the periodic boundary may explain the observed reaction energy differences. To

test this hypothesis directly, it would be necessary to re-calculate the reaction energies in a series of

supercells of different sizes, to check whether the reaction energy difference scales with 1

r3, like the dipole-

dipole interaction does. At the current state of the art this is impossible, since the employed supercells

are already at the boundary of what can be handled with reasonable computational effort. Shrinking

the supercells is not a sensible option, as this incurs the danger of introducing further self-interaction

problems, beyond the electrostatic self-interaction. Since the number of atoms Nscales as a2with the

lateral length of the supercell and the computational cost scales proportional to N3, doubling the supercell

size leads to four times the memory requirements and (at least) 64 times the CPU requirements. Due to

these constraints, a systematic study of this problem is impractical.

4.4 Summary

The evaluation of the tested embedding schemes leads to two results. Firstly, it could be shown that

QM/MM embedding with DFTB as the QM method can give reliable results, so that the method is suit-

able to simulate solid surfaces and interfaces which would would be inaccessible to pure QM calculations,

due to their size. Secondly, the tests show that the necessary embedding scheme strongly depends on the

material system at hand. While in SiO2, non-stoichiometric QM clusters with homogeneous QMH atoms

can be successfully employed, this cluster construction leads to severe problems in Al2O3. It becomes

also clear, that in polar materials, the clusters must be carefully neutralized; here the BCTC method,

developed in the course of this work proves well suited for solid surfaces.

58 CHAPTER 4. DFTB QM/MM EMBEDDING

Chapter 5

Epoxy adhesives on native Al2O3

During the past two decades aluminum has constantly gained importance in technical applications, with

its use spreading from aerospace to automotive applications and now covering nearly every area of indus-

trial and consumer appliances. Parallel to this development, adhesive technology has achieved advances

which enable adhesive bonding of metal parts to supplement or even replace traditional metal joining

technologies such as welding, bolting or riveting. Related to this is the development of metal-resin or

metal-resin-fiber compound materials, which have recently advanced to the point of technical application.

The demands for lighter and at the same time stronger materials, in order to build more fuel-efficient, i.e.

lighter, vehicles and aircraft without sacrificing structural strength leads to a strong interest in the devel-

opment and improvement of fiber reinforced aluminum-polymer hybrid materials and adhesive bonding

technology for aluminum. For both of these technologies, the aluminum-polymer adhesion is of crucial

importance.

Since, outside high vacuum environments, aluminum instantly develops a surface oxide layer[131,132],

the problem of adhesion on aluminum translates into the problem of adhesion on native aluminum oxide.

The improvement of adhesion technology requires an understanding of the underlying chemical processes

of the bonding of organic adhesives to the alumina surface. In this study I aim at finding a suitable

methodology for gaining insight into the initial bonding of adhesive component molecules as well as the

bonding competition between different organic species at the native Al2O3-surface. To achieve this, I

concentrate on a model adhesive system, comprising diglycidylesterbisphenol-A (DGEBA) as the resin,

diethyltriamine (DETA) as the curing agent and 3-aminopropyltrimethoxysilane (sold as: Dynasilan

AMEO, here referred to as AMEO, also known as APTES) as an adhesion promoter component. Although

our results are focused on an adhesive, some of the examined components are also used in compound

materials or paints, and the method is sufficiently general, to be transferable to the adsorption of any

organic molecule of similar size on inorganic surfaces.

5.1 The native Al2O3surface

Native alumina, created by corrosion of Al in air, is usually described as γ-Al2O3, which is formed by a

pseudo-amorphous defective spinel structure in which Al atoms and vacancies are randomly distributed

among the tetrahedral and octahedral sites within a perfect fcc sub-lattice of oxygen atoms[133–138].

However there is still controversy in the literature concerning the crystallography of γ-alumina. An

extensive study, combining experimental and theoretical approaches, suggests, that the description of γ-

alumina requires a significant number of Al atoms on non-spinel positions to be accurate[139]. Extensive

review of the literature concerned with the bulk structure of γ-alumina is given in Refs. [140,141]. As

the underlying bulk structure is amorphous, the surface does not exhibit any regular reconstruction. The

main characteristics of the surface are the distribution statistics and density of surface Al on tetrahedral

60 CHAPTER 5. EPOXY ADHESIVES ON NATIVE AL2O3

and octahedral sites as well as on the distribution of hydroxyl groups at the surface. A mesoscopic model

of the (100) surface of spinel-like γ-alumina and its hydroxylation has already been presented in Ref. [142].

It is well known that in a moist environment, such as air, the Al2O3surface quickly absorbs water,

both molecularly and dissociatively[142–145]. It is therefore essential to include the hydroxylation of the

surface in any study of its chemistry. Since alumina in general and especially γ-alumina is of high interest

as a catalyst for the chemical industry, a large number of studies on the surface reactivity have been

performed, most of which are concerned with the surface acidity. An extensive overview can be found in

[138,146] and references therein.

5.2 A model adhesive system

Investigations on the adsorption of methanol at the hydroxylated surface show that these react with hy-

droxyl groups attached to surface aluminum atoms by condensation.[147,148]. A study of the adsorption

of maleic acid anhydride[149], and theoretical investigations[150] on the adsorption of adhesive compo-

nents on amorphous SiO2and hydroxylated α-Al2O3also find condensation reactions on the alumina

surface, except for DGEBA, which is found to react with the surface by additive ring-opening. I therefore

assume that similar reactions will occur in the system studied here, however, we assume the ring-opening

of DGEBA to occur in the liquid phase of the adhesive mixture, independent of the surface. Therefore

we examine the ring-opening of DGEBA by water and the condensation reaction of the opened DGEBA

separately.

OO OO

H3C CH3

H2NSi(OH)3

H2N

NNH2

H3C

H2NSi(OC2H5)3

H3C

(a)

(b) (c)

(d)

(e) (f)

Figure 5.1: The molecules of the model adhesive system: (a) DGEBA, (b) t-DGEBA, (c) t-

DGEBA(open), (d) DETA, (e) AMEO, (f) AMEO(hyd). From [151].

In a recent study of the DETA absorption energetics[152] reaction energies of ∼35kcal mol−1have been

observed for a chemisorption at an octahedral surface aluminum site. In this study, I concentrate on the

interaction of the adhesive components with surface hydroxyl groups, which results in a condensation

reaction for the adsorption of DETA.

As DGEBA has a mirror plane, I consider a truncated molecule which only contains one phenol-ring and

one epoxy group, referred to as t-DGEBA (c.f. Fig. 5.1(a),(b)). Since AMEO can be expected to rapidly

lose its ethyl groups and hydrolyze in an aqueous solution, I only study the interaction of the hydrolyzed

form of AMEO which I refer to as AMEO(hyd) (c.f. Fig. 5.1(e),(f)). In my model, the adsorption of

DGEBA depends on the prior opening and hydroxylation of at least one oxirane ring. This catalyzed

5.3. SURFACE MODEL PREPARATION 61

reaction is not limited to the alumina surface but must happen throughout the adhesive mixture, the ring-

opened, truncated DGEBA molecule will henceforth be referred to as t-DGEBA(open) (c.f. Fig. 5.1(c)).

After these considerations, I examine the following reactions at the surface:

t-DGEBA(open) + surf. ∆E

−→ (t-DGEBA + surf.) + H2O (5.1)

DETA + surf. ∆E

−→ (DETA + surf.) + H2O (5.2)

AMEO(hyd) + surf. ∆E

−→ AMEO(hyd) + surf.+ H2O (5.3)

where surf. denotes the hydroxylated surface model and (species+surf.) denotes the structure of the

organic species adsorbed at the surface. Additionally, I also consider the reaction

t-DGEBA + 2H2O∆E

−→ t-DGEBA(open) + H2O.(5.4)

Here, the second water molecule stabilizes the ring-opening.

5.3 Surface model preparation

5.3.1 Bulk γ-Al2O3

The results presented in this subsection were published in a joint publication with Blumenau and

Amkreutz [151]. They are presented here, as the foundation of the hydroxylated model construction

described below. As stated in section 5.1, the structure of γ-alumina is still under discussion in the liter-

ature. The consensus is however, that it resembles a close packed oxygen (fcc-) lattice with Al atoms and

vacancies statistically distributed among the tetrahedral and octahedral sites. It can be generated either

from a strained hausmannite or a spinel structure, by removing Al atoms to achieve Al2O3stoichiometry.

Due to the statistical distribution of vacancies, the structure can be regarded as pseudo-amorphous.

structure Ncell EDFTB

tot [eV/Al2O3]ELDA

tot [eV/Al2O3]DFTB [g cm−3] Al–ODFTB [˚

hausmannite 80 −334.98 −1425.70 3.65 1.88 (0.07)

spinel cub. 1 160 −335.35 −1426.30 3.66 1.88 (0.07)

spinel cub. 2 160 −335.41 −1426.36 3.66 1.88 (0.07)

spinel cub. 3 160 −335.39 −1426.37 3.64 1.88 (0.06)

spinel rhomb. 80 −335.05 −1425.74 3.66 1.88 (0.06)

Table 5.1: SCC-DFTB and LDA results on γ-Al2O3.Ncell is the number of atoms per unit cell,

Al-O is the mean Al-O bond length and its statistical error. From [151].

Five different models of γ-Al2O3were constructed, one based on a hausmannite cell, one based on a

rhombohedral representation of spinel and three based on a cubic representation of the spinel structure.

For the cubic spinel case, three different starting structures in which different Al atoms have been ran-

domly removed were examined. The atomic positions were relaxed in all cells, optimizing the cell volume.

Table 5.1 lists the characteristic data of the relaxed cells. We find that the differences in terms of density

of the optimized unit cell and Al-O bond lengths are negligible and that the densities correspond well to

the experimentally reported value of 3.67 g cm−3[153]. The energy differences between the three different

basic structures (i.e. hausmannite, cubic spinel and rhombohedral spinel) deviate slightly, indicating that

the cubic representation is slightly favorable. Notably, the differences between the three different cubic

spinel based models are negligible in all respects.

62 CHAPTER 5. EPOXY ADHESIVES ON NATIVE AL2O3

Figure 5.2: Views along the (001) (left) and (010) (right) directions of the hydroxylated γ-Al2O3

surface model. The selected adsorption site 1 is marked blue, site 2 is marked green. (Note that,

due to the periodic repetition, each adsorption site appears twice.)

5.3.2 The hydroxylated γ-Al2O3surface

To generate a model of the γ-Al2O3surface, the first of our 4 oxygen-layer thick, 160-atom cubic-spinel

based models was chosen (second line in table 5.1). The bulk alumina was cut by increasing the length

of the cell vector along the (001) direction to 60 ˚

Awithout changing any of the atomic positions. The

resulting slab model was then relaxed using conjugate gradients (CG).

To hydroxylate the surface, a 10 ˚

Athick layer of H2O molecules, cut from a larger, geometry-optimized

supercell of liquid water, was placed on top of the surface model and then performing a 1 ps SCC-DFTB

MD run at 300 K, followed by a second 1 ps MD run at 600 K to remove the non-adsorbed water

molecules. This procedure does not lead to chemisorption of water molecules at the surface, but since

sufficiently long MD runs to achieve this are computationally far too expensive, I then manually split

the physisorbed water molecules, forming an OH-group at the adsorption site and using the split-off

hydrogens to protonate nearby Al-O-Al bridges. After this procedure, I relax the geometry again using

CG. During the whole procedure, the lowermost oxygen and aluminum layers were held fixed.

The surface supercell generated this way has a surface area of ∼24 ˚

A x 8 ˚

A, which is too small in (0,1,0)

direction, wherefore it was duplicated model along this axis. The resulting model contains 374 atoms in

a 23.7 x 15.8 x 60 ˚

A3periodic supercell (figure 5.2). To limit the computational cost of the study, the

two sites marked in in the figure were chosen out of the 13 inequivalent –OH groups for the studies of

adsorption reactions, based on the different distribution of neighboring hydroxyl groups and tetrahedrally

coordinated aluminum atoms.

5.4 The H2O-assisted ring-opening of DGEBA

In order to gain an understanding of the chemical reactions from which the alumina/adhesive interface

arises, I calculate the reaction enthalpies and MEPs of the condensation reactions between the model

adhesive system’s components and the surface hydroxyl groups at the two selected reaction sites.

As described in section 5.2, the nucleophilic ring-opening of the oxirane rings in DGEBA by water

molecules occurs throughout the adhesive mixture and is independent of the surface. Therefore this

reaction is modeled independently of the surface adsorption reactions. To conserve calculation time,

5.4. THE H2O-ASSISTED RING-OPENING OF DGEBA 63

-2

-1.5

-1

-0.5

0.5

1.5

0.2

0.4

0.6

0.8

images

interpol.

(a) (c)(b)

Figure 5.3: Minimum Energy path of the ring-opening of t-DGEBA by H2O. Atoms of the phenyl

group are rendered transparent for clarity, circles mark energies of the configurations shown above.

I regard only one of the oxirane rings in the molecule and therefore simulate a t-DGEBA, instead of

DGEBA (cf. figures 5.1,5.3). In the mixture, the DGEBA molecules are surrounded by other components

with polar groups as well as some water molecules. As a simple model of the stabilizing effect of such

surrounding polar groups, a water molecule is included in the ring-opening path-search calculations[154].

I use the SCC-DFTB method, applying the alo-1 Slater-Koster data sets. The SCC is iterated until the

energy difference between two subsequent iterations falls below 10−6H and determine the minimum energy

configurations of reactants and products by CG optimization until the forces fall below 10−3a.u. The

MEP is then searched using NEB, where the path optimization was started by an adaptive displacement

steepest-descent algorithm for the first few hundred iterations, then switching to a projected velocity-

Verlet scheme, which is iterated until the maximum forces at the images fall below 10−3a.u. I construct

the initial path by linear interpolation of the Cartesian coordinates between the reactant, a guessed

transition state and product configurations. The transition sate guess is necessary, since the reaction

coordinate has dihedral rotation components, which cannot be reproduced by linear interpolation in

Cartesian coordinates. A reaction coordinate is segment-wise as the linear displacement between images:1

ci=~

R−~

Ri

with

li=~

Ri+1 −~

Ri

1Note that this differs from the reaction coordinate used in Ref [151].

64 CHAPTER 5. EPOXY ADHESIVES ON NATIVE AL2O3

for a position ~

Ralong the path in the segment between images iand i+ 1, leading to the normalized

position along the reaction path

x=Pi

j=1 lj+ci

Pjlj

The total energies are interpolated along this reaction coordinate, following the interpolation procedure

described in Ref [108]2:

E(xi) = −Ei+xi· −Ft

i(5.5)

+xi2 3 (Ei+1 −Ei)

+−2Ft

i−Ft

i+1

li!

+x3

i−2 (Ei+1 −Ei)

−Ft

i−Ft

i+1

i.

The climbing image NEB method is not employed in this study for two reasons. First, the interest rests

not only on finding the barrier energy, but on obtaining an overview of the whole reaction mechanism.

Second, not all reaction paths examined here have well defined barrier, e.g. the AMEO adsorption

reactions (cf. section 5.5.3).

Fig 5.3 shows the calculated image energies and structures (crosses), as well as the cubic interpolation of

the total energy (red line). I find a barrier energy of ∼1.91 eV and an exothermicity of ∼1.88 eV. No

direct comparison data is available, however for the similar reaction

C2H4O

|{z }

oxirane

+H2O−→ H2HOC −CH2OH

|{z }

1,2 Ethanediol

reference data using PBE and an 6-31G basis [155] gives a reaction enthalpy of -1.31 eV, which is in good

agreement with my calculation results, taking into account the formation of a hydrogen bond between

one of the formed OH-groups and the catalyzing H2O. Before the formation of this bond, the reaction

energy is ∼ −1.2 eV (cf. figure 5.3).

5.5 Adsorption of single adhesive component molecules in vac-

uum

To examine the energetics of reactions 5.1–5.3, I follow the procedure verified above for the t-DGEBA

ring-opening. Each of these reactions is simulated at both adsorption sites described in section 5.3.2.

The results of these calculations are summarized in table 5.2. Two reaction energies for the condensation

of each species are given: Ecomb.

ads and Esep.

ads . The former is extracted from the adsorption paths, as the

difference between final and first image. It includes an additional contribution from the physisorption of

the water molecules generated during the reaction. The latter was calculated from separate simulation of

reactants and products and is therefore free of such influences. It can be seen, that the water adsorption

results in an exothermic shift of the reaction energies of ∼0.5–∼2 eV. These shifts are unrealistic

under real conditions, since the oxide surface is known to be covered by at least one monolayer of water

molecules [142–145] — a site at which water may be adsorbed with an energy gain of 2 eV will surely

be already occupied. On the other hand, Esep.

ads is unrealistically low, as the solvation of the condensated

water molecule will always release some energy. However, as table 5.2 shows, these uncertainties do not

change the relative energetics qualitatively. In the following I shall give a detailed description of the

reaction mechanisms found for each compound. (Figures showing the adsorption MEPs are given on

pages 67–72.)

2Note, that the equations given in reference [108] do not match the derivation described in the text, due to several errors.

Here, I present the corrected energy interpolation formula.

5.5. ADSORPTION OF SINGLE ADHESIVE COMPONENT MOLECULES IN VACUUM 65

component site Ecomb.

ads [eV] Esep.

ads [eV] Ebarr [eV]

t-DGEBA ringopen -1.8780 - 1.9081

t-DGEBA(open) 1 -0.42741 -0.11970 2.0726

2 -2.5298 -1.6108 0.1282

AMEO(hyd) 1 -4.0815 -2.1364 0.0∗

2 -4.4355 -2.0329 ∼0.35∗

DETA 1 0.20411 2.3602 1.1583

2 -0.47487 2.0022 ∼0.8∗

Table 5.2: Calculated reaction energies and barriers of the adsorption and ring-opening reactions of

the model adhesive system on Al2O3.Ecomb.

ads gives the adsorption energy with physisorption of the

condensed H2O, Esep.

ads with the condensed water isolated in vacuum. ∗Cf. sections 5.5.2 and 5.5.3

for a discussion of these values.

5.5.1 DETA

Table 5.2 clearly shows, that the condensation of the curing agent DETA is very unfavorable. The

reaction energies without the contribution of the physisorption of the condensed water molecule Esep.

ads are

endothermic by more than 2 eV. At site 2 the adsorption path has a local minimum before the actual

condensation, (cf. figure 5.6) which is in agreement with findings, that DETA can bond to the surface

non-covalently via hydrogen bonding [152]. The Barrier given for this reaction is calculated between the

local minimum and the following maximum.

5.5.2 t-DGEBA(open)

The picture for t-DGEBA(open) shows significant differences between both adsorption sites. While at

site 2 a second, ternary hydroxyl bond to the surface is formed, between the second hydroxyl group of

the opened epoxy ring and a clean Al atom at the surface(cf. figure 5.8), at site 1 only the C-O-Al bridge

from the condensation reaction appears (cf. Fig 5.7). This leads to large differences in the reaction

energetics at both sites — at site 1 the adsorption is only slightly exothermic, while being hindered by

a 2 eV barrier3; at site 2 only a negligible barrier is found, while the condensation reaction is highly

favorable. At both sites, the transition state is characterized by a strong repulsion between the surface-

and molecules OH groups — the Al-O-C bridge can only form after one of the hydroxyl groups has been

release by proton transfer from the other. At site 2 this transition state is stabilized by the formation of

a the ternary OH group with the free surface Al atom.

5.5.3 AMEO(hyd)

I only examine the adsorption of the silanol function of AMEO, since its amine group will behave very

similar to the amine groups of DETA. Their interaction with the surface is so unfavorable, that a

qualitatively different behavior of the amine function in AMEO cannot be expected.

For the condensation of AMEO, no barriers appear between the starting configurations ∼5˚

Aabove the

surface. The barrier given for site 2 is most likely caused by a local valley of the electrostatic potential,

which the AMEO molecule slides into during the approach to the surface, but has to leave before the

condensation reaction. At both sites, the AMEO molecule forms two covalent Si-O-Al bridges to the

surfaces, by reacting with a second, unhydroxylated Al atom nearby (cf. figures 5.9,5.10). Also, at

3which is large in comparison with thermal energies at ambient temperatures

66 CHAPTER 5. EPOXY ADHESIVES ON NATIVE AL2O3

-6

-5

-4

-3

-2

-1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

E [eV]

reaction coordinate

t-DGEBA AMEO DETA

Figure 5.4: Comparison of the calculated adsorption MEPs on the Al2O3surface.

both sites the adsorption of the condensated water molecule is very exothermic. Fig 5.10 shows, that at

site 2 the transition state is marked by a five-fold coordination of the Si atom by four hydroxlys (three

binary and one ternary) and one carbon. Subsequently one of these binary OH groups is released by a

short-range proton transfer from the ternary hydroxyl bridge. Neither the OH group formation, nor the

proton transfer show significant barriers.

5.5.4 Comparison of adsorption mechanisms

Figure 5.4 illustrates the energetic differences between the adsorption mechanisms of the model adhesive

system. It shows very clearly, that AMEO is the most favorable adsorbate at both surface sites (green

lines). It is also the the only adsorbate, which shows no significant adsorption barriers at both sites.

In contrast, the chemisorption of DETA (blue lines) is unfavorable. It can rather form a hydrogen

bond, which is, however, much weaker than the favorable covalent bonds possible in this system. The

resin component, DGEBA (red lines) shows diverse behaviors in the sense, that it can compete with

the AMEO adsorption at site 2, while its chemisorption at site 1 is thermodynamically on par with the

physisorption of DETA, while being hindered by a high barrier.

The differences between DGEBA and AMEO can be understood, when scrutinizing the transition states.

The silicon atom of AMEO allows for a five-fold coordinated configuration, which makes it possible to

form the Al-O-Si bridge, before the proton transfer and H2O condensation. In contrast to this, the

carbons of the opened oxirane rings in DGEBA do not have this possibility. In addition to this, DGEBA

with its bulky phenyl ring exhibits a much higher steric hindering, so that it is not always possible to form

additional, ternary OH-group links to surface Al atoms. The much more linear and smaller AMEO(hyd),

in contrast, can adapt much easier to the local conditions at the adsorption site.

These findings allow an explanation of the working mechanism AMEO as an adhesion promoter. In an

adhesive mixture without adhesion promoter, only the resin component can form strong, covalent bonds

to the native alumina surface. However, only a part of the available surface hydroxyl groups is available for

adsorption. At other sites, the formation of covalent bonds between DGEBA and surface is not distinctly

favorable and must compete with the physisorption of curing agent molecules. The silane-based adhesion

promoter AMEO can improve the linkage between surface and polymer, mostly because its adsorption

is favorable at more surface sites, especially such sites, where the formation of a covalent bond between

DGEBA and the surface is not. Therefore, it can be concluded, that the chief effect of the adhesion

promoter is to increase the number and density of covalent links between polymer and surface.

Since the favorable effect of AMEO is closely related to the chemistry of its Si atom, it can be concluded

that similar silanes should show a similar effect, as long as no too strong steric hindering from its backbone

or other functional groups occurs.

5.5. ADSORPTION OF SINGLE ADHESIVE COMPONENT MOLECULES IN VACUUM 67

Adsorption path figures

-0.2

0.2

0.4

0.6

0.8

1.2

0.2

0.4

0.6

0.8

images

interpol.

(a) (c)

(d) (e) (f)

(b)

Figure 5.5: Minimum Energy path of the condensation of DETA at site 1. Insets show a cutout of

the surface model in the vicinity of the adsorption site.

68 CHAPTER 5. EPOXY ADHESIVES ON NATIVE AL2O3

-1.2

-1

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

images

interpol.

(a) (c)

(d) (e) (f)

(b)

Figure 5.6: Minimum Energy path of the condensation of DETA at site 2. Insets show a cutout of

the surface model in the vicinity of the adsorption site.

5.5. ADSORPTION OF SINGLE ADHESIVE COMPONENT MOLECULES IN VACUUM 69

-1

-0.5

0.5

1.5

2.5

0.2

0.4

0.6

0.8

images

interpol.

(a) (c)

(d) (e) (f)

(b)

Figure 5.7: Minimum Energy path of the condensation of t-DGEBA at site 1. Insets show a cutout

of the surface model in the vicinity of the adsorption site.

70 CHAPTER 5. EPOXY ADHESIVES ON NATIVE AL2O3

-3

-2.5

-2

-1.5

-1

-0.5

0.5

0.2

0.4

0.6

0.8

images

interpol.

(a) (c)

(d) (e) (f)

(b)

Figure 5.8: Minimum Energy path of the condensation of t-DGEBA at site 2. Insets show a cutout

of the surface model in the vicinity of the adsorption site.

5.5. ADSORPTION OF SINGLE ADHESIVE COMPONENT MOLECULES IN VACUUM 71

-4.5

-4

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0.2

0.4

0.6

0.8

images

Interpol.

(a) (c)

(d) (e) (f)

(b)

Figure 5.9: Minimum Energy path of the condensation of AMEO(hyd) at site 1. Insets show a cutout

of the surface model in the vicinity of the adsorption site.

72 CHAPTER 5. EPOXY ADHESIVES ON NATIVE AL2O3

-6

-5

-4

-3

-2

-1

0.2

0.4

0.6

0.8

images

interpol.

(a) (c)

(d) (e)

(b)

Figure 5.10: Minimum Energy path of the condensation of AMEO(hyd) at site 2. Insets show a

cutout of the surface model in the vicinity of the adsorption site.

5.5. ADSORPTION OF SINGLE ADHESIVE COMPONENT MOLECULES IN VACUUM 73

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0.2

0.4

0.6

0.8

images

interpolation

(a) (b)

Figure 5.11: Minimum Energy path of the coadsorption of AMEO at site 1 in the vicinity of AMEO

at site 2.

74 CHAPTER 5. EPOXY ADHESIVES ON NATIVE AL2O3

-8

-7

-6

-5

-4

-3

-2

-1

0.2

0.4

0.6

0.8

images

interpolation

(a) (b)

Figure 5.12: Minimum Energy path of the coadsorption of DGEBA at site 1 in the vicinity of AMEO

at site 2.

5.6. THE EFFECT OF NEIGHBORING ADSORBATES 75

Site 1 Site 2 Etot

ads ∆Eads

AMEO AMEO -4.38 -0,50

DGEBA -3.22 0,23

DETA 0.26 0,08

DGEBA AMEO -2.17 -0,01

DGEBA -1.89 -0,17

DETA 1.90 -0,01

DETA AMEO 0.25 -0,06

DGEBA -7.01 -7,76

DETA 4.40 0,03

Table 5.3: Influence of nearest neighbor adsorbates on adsorption energies. Energies in eV, both

adsorbates in the same half of the surface model (cf. section 5.3). ∆E=Etot

ads −E1

ads + E2

ads.

Ads. 1 Ads. 2 Site Etot

ads ∆Eads

AMEO AMEO 1 -3,84 -0,16

2 -4,09 0,00

DGEBA 1 -2,27 -0,32

2 -3,66 -0,01

DETA 1 0,49 -0,02

2 2,62 2,64

DGEBA DGEBA 1 -0,22 0,00

2 -3,20 0,01

DETA 1 2,24 0,00

2 0,42 0,00

DETA DETA 1 4,71 0,00

2 4,05 0,01

Table 5.4: Influence of second neighbor adsorbates on adsorption energies. Energies in eV, both

adsorbates at the same site in different halves of the surface model (cf. section 5.3). ∆E=Etot

ads −

E1

ads + E2

ads.

5.6 The effect of neighboring adsorbates

Beyond the adsorption of isolated adhesive components, the influence of neighboring adsorbates on the

condensation reactions can be used to gain further insight into the alumina–adhesive interface formation.

5.6.1 Adsorption energies

To examine, how the presence of an adsorbate influences the condensation reaction of a second molecule

from the adhesive mixture, I calculate the reaction energies for these reaction in the presence of a second

adsorbate in first or second nearest neighbor position. Tables 5.3 and 5.4 list the results of these

calculations.

For the interaction between nearest-neighbor adsorbates, the adsorption of two AMEO(hyd) molecules in

close vicinity is very favorable, while putting a t-DGEBA(open) beside the AMEO is unfavorable. A slight

advantage of condensating two t-DGEBA molecules in close vicinity can be inferred from the data. The

co-adsorption of DETA at site 1 and DGEBA at site 2 appears to be extremely favorable. As can be

76 CHAPTER 5. EPOXY ADHESIVES ON NATIVE AL2O3

Figure 5.13: Co-adsorbed structure of t-DGEBA(open) at site 2 and DETA at site 1. Dashed lines

indicate hydrogen bonds.

seen in figure 5.13, the molecules are oriented nearly parallel in this case, leading to several hydrogen

bonds and a strong electrostatic interaction between both molecules. It should be kept in mind, that here

the simplified model of the adhesive components and surface can lead to unrealistically high adsorption

energy estimates. In this case, it is to be expected that in a liquid polymer mixture, the bond sites

interacting here, will already be occupied by hydrogen bond partners before the adsorption reactions. In

all other cases, the adsorption energy differences are <0.1 eV, which is insignificant, taking into account

general limitations of the employed model (cf. 5.3.2).

As expected, table 5.4 shows that the presence adsorbates on the second dearest neighbor sites has

much smaller influence on the condensation reaction energetics. In comparison to the nearest neighbors’

influences, it can be seen that co-adsorption of AMEO is still favorable at longer distances, at least at site 1.

Interestingly, the co-adsorption of AMEO and DGEBA is significantly more favorable at second nearest

neighbor sites, indicating that the pronounced shift towards endothermicity in the nearest neighbor site

case is caused by steric hindering. Similarly, steric hindering is culpable for the clear unfavorability of

co-adsorbing AMEO and DETA. A different orientation of the DETA, which adsorbs nearly parallel to

the surface at this site (cf. figure 5.6) may alleviate this, however it should be kept in mind that the

adsorption of DETA is highly unfavorable in any case. In the vast majority of examined reactions, the

presence of a second neighbor adsorbate did not influence the adsorption energetics significantly.

5.6.2 Adsorption paths

The computational effort needed to calculate the adsorption paths presented in section 5.5 makes a

systematic study of the influence of the surroundings on the adsorption reactions impractical. Therefore I

switch to the DFTB based QM/MM scheme presented in chapter 4. The Al36O56H12 QM Zone containing

both adsorption sites verified in that chapter is employed. Here I will only examine two coadsorption

reactions, to demonstrate the suitability of my QM/MM approach to simulate complex reactions at hybrid

interfaces. The adsorptions of an AMEO(hyd) molecule at site 1 in the presence of a second AMEO at

site 2 (with the water molecule created during the condensation of AMEO at site 2 still physisorbed at the

surface), and the adsorption of a t-DGEBA(open), also at site 1 in the same surroundings were simulated.

5.6. THE EFFECT OF NEIGHBORING ADSORBATES 77

-4.5

-4

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0.5

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

E [eV]

reaction coordinate

isolated coadsorption

(a) AMEO

-8

-7

-6

-5

-4

-3

-2

-1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

E [eV]

reaction coordinate

isolated coadsorption

(b) DGEBA

Figure 5.14: Comparison of site 1 adsorption energy profiles for isolated adsorption and in the

vicinity of a second adsorbed AMEO.

AMEO + AMEO Figure 5.14(a) shows the energy profiles for the condensation of AMEO at site 1

with and without a neighboring adsorbed AMEO. The first notable difference lies in the total energy of

the reaction. This difference is caused by the fact, that the water molecule created during the course of

the condensation reactions at both sites physisorbs at the same surface Al atom. Therefore, another place

had to be chosen for the second water molecule in the coadsorption case. As can be seen, the adsorption

of the water molecule is very favorable at the original surface atom. The second, stark difference is in

the existence of a significant secondary barrier of ∼0.95 eV for the physisorption of the second water

molecule. This difference is also related to the choice of a different site for the second water molecule (cf.

figure 5.11). However, the barrier is still small enough that it can be overcome by the energy gained from

the condensation reaction itself.

Besides the differences stemming from the different disposition of the water molecule, the adsorption

reactions occur in a very similar manner (cf. figures 5.9 and 5.11). It can be concluded, that the conden-

sation reaction of AMEO is not strongly influenced by the presence of neighboring AMEO adsorbates.

This is in good agreement with the results on the thermodynamics, presented in the preceding section.

AMEO+DGEBA In contrast to AMEO, t-DGEBA(open) reacts sensitively to the presence of a neigh-

boring AMEO adsorbate or, more precisely, to the presence of the water molecule formed during the

condensation of the latter. The presence of a neighboring AMEO does not significantly influence the

adsorption enthalpy. However, the deep intermediate minimum in the path shows that another reaction

that does not lead to condensation of the DGEBA(open) is actually much more favorable. When the

geometry is optimized, starting from one of the intermediate geometries between reaction coordinates

0.3–0.4 of the coadsorption path (cf. figure 5.12(c)), the resulting minimum energy configuration is very

different from the result of the condensation reaction (figure 5.15 4). The physisorbed water molecule

from the preceding AMEO condensation is split up to form a surface hydroxyl group. The proton released

by this reaction then attacks the the opened oxirane ring and passivates one of it’s –OH groups. Finally,

the dangling carbon bond resulting from the preceding step is passivated by another proton transferred

from a nearby Al–(OH)–Al bridge at the surface.

Even taking into account the inevitable error of this calculation due to the too small surface cluster, this

result is gives extremely interesting insight into the interface formation process at hand. It clearly shows

4Note that the final geometry from this calculation exhibits the problem, that the free products of the reaction, i.e.

the passivated DGEBA and the water molecule marked orange, have moved beyond the surface cluster. Therefore this

equilibrium geometry is not realistic and its energy is useless. Nevertheless it shows the main features of the alternate

reaction path found here.

78 CHAPTER 5. EPOXY ADHESIVES ON NATIVE AL2O3

Figure 5.15: DGEBA(open) passivated by a physisorbed water molecule stemming from the nearby

condensation of AMEO. The formerly physisorbed water forms an –OH group (marked blue),

transferring one of its protons to one of the –OH groups of the opened oxirane ring which then

forms a water molecule (orange). The dangling C- bond is then saturated by a hydrogen atom

taken from a formerly protonated Al–O–Al surface bridge (green).

that water at the surface plays a very important role. Further, systematic studies are necessary, which

take into account the presence of a physisorbed water layer between the surface hydroxyl groups. This

singular result also points at a mechanism of the effect of water at the surface that is likely to be highly

important. Depending on the pH value and surface hydroxylation density5, surface water can strongly

react with the surfactants by proton transfer processes. These can lead to the passivation of functional

groups of the adhesive components, which may lead to a degradation of the adhesion between polymer

and substrate.

Comparison of these first results on AMEO and DGEBA suggests, that the epoxy resin component is

more susceptible to this passivation process than the silane adhesion promoter. This is in good agreement

with the experimental observation that silane components improve the adhesion significantly. It is also

consistent with the conclusion from section 5.5.4, that the effect of the adhesion promoter appears to be

based on an increase of the density of covalent polymer–substrate bonds.

5.7 Summary

The results presented in this chapter, demonstrate, that minimum energy paths, obtained from NEB

calculations and the DFTB method can provide very useful insight into the basic mechanisms of adhe-

sive molecule adsorption on inorganic surfaces. The different components of the adhesive mixture show

qualitatively different adsorption energetics, so that, consistent with experimental observation, a clear

preference for the adsorption of the silane adhesion promoter could be established. At the same time,

it could be demonstrated that condensation of the amine curing agent is unfavorable. Furthermore, the

calculated MEPs also provide insight into the mechanisms leading to these differences. The condensation

of amine functionalities is unfavorable, because necessary the removal of the surface hydroxyl group is

extremely unfavorable. The silane adhesion promoter has more possibilities to form bidentate links to

the surface, which is more difficult for the DGEBA resin, chiefly because of steric hindering at part of

the available adsorption sites.

5Which had to be neglected so far, due to of the prohibitive computational cost of modeling them.

5.8. OUTLOOK 79

The thermodynamics of the adsorption reaction show some influence from the presence of neighboring

adsorbates, however these differences are small compared to the adsorption energies. First results on

the adsorption paths in the presence of neighboring adsorbates give two important results. Firstly, the

adsorption of AMEO is only weakly influenced by the presence of adsorbates on neighbor sites, which

supports the conclusion that AMEO can serve to increase the number of covalent links between substrate

and polymer. Secondly, the presence of water has a much larger influence on the behavior of DGEBA

at the interface, that neighboring adsorbates do. The results propose that, depending on the pH value,

physisorbed water can react with the ring-opened DGEBA in a manner which partly passivates the –OH

groups of the opened oxirane ring and leads to further hydroxylation of the surface.

Finally, it must be stated that the results obtained so far are not fully satisfactory in three respects:

first, only a small fraction of the adsorption sites present in the surface model could be examined, better

statistics would be highly desirable. Secondly, the influence of neighboring adsorbates could only be

analyzed in the reaction energies but not the barriers. The third and most severe limitation lies in the

examined model itself. The influence of physisorbed surface water and the surrounding polymer mixture

had to be neglected completely. All of these limitations are based on the prohibitive computational cost

that would otherwise arise. QM/MM coupling promises to alleviate these problems, and first results on

the coadsorption of different components support this expectation. They show that QM/MM is indeed

capable to give further insight into interface chemistry at greatly reduced computational cost, compared

to pure QM simulations. It has become clear, that further studies of the interface formation reactions

are needed, especially concerning the influence of water, and that these studies can now be conducted in

practice.

5.8 Outlook — Kinetics simulation of the adhesive-alumina in-

terface formation

Even considering the possibilities (and sometimes necessity) for improving on the findings so far obtained

on the mechanisms of polymer/solid interface formation on an atomic scale, the limits of the atomistic

perspective become apparent. Even with big improvements in computational efficiency, the accessible

length and time scales cannot be extended qualitatively. In so far, simulating the whole process of

interface formation on an atomic scale appears unrealistic. However, qualitative and quantitative insight

into the chemistry of interface formation at the high level of detail attained in this study, can serve as the

basis for examining larger time and length scales. One very promising way to achieve this lies in using the

available results as input data to kinetic Monte Carlo simulation of the the interface formation. Reaction

energies and barriers, which are the basic input data for this approach (cf. section 3.3) are available,

even allowing to (partly) take into account the influence of neighboring sites. Studies to employ the

data presented here as input to the KMC part of a coupled KMC/continuum model [127] are underway

in cooperation with the developers of this scheme. Valuable insights into the larger scale mechanism of

interface formation and resulting interface properties, e.g. the distribution of different adsorbates along

the surface, are expected from this ongoing research.

80 CHAPTER 5. EPOXY ADHESIVES ON NATIVE AL2O3

Chapter 6

Conclusion

During the work presented here, advancements in both computational methodology and the understanding

of the alumina/epoxy polymer interface were achieved.

Concerning methodology, the applicability of the DFTB/MM coupling method was extended from bio-

molecules to solid state systems and especially solid surfaces. To this end, a bond-termination scheme

employing ordinary hydrogen atoms at fixed non-equilibrium distances to their QM host atoms is used.

The bond distances are optimized to remove bonding states between saturators and host atoms from

the band gap of the solid model. The electrostatic embedding of the QM zone into the MM zone is

implemented by a Hamiltonian which allows to apply spherical Gaussian smearing to the point charges

of the MM zone. This allows to circumvent problems originating from the unphysical proximity of the

saturator atoms to the first layer of MM point charges. The two free parameters of this embedding

scheme, i.e. the link-atom distance scaling and the Gaussian smearing width, must be optimized for each

materials system to minimize charge distribution errors in the QM zone. The results show, however, that

the optimized parameters are, within their error bars, independent on the chosen QM zone. They are also

very similar between closely related materials, e.g. different phases of SiO2or Al2O3, clearly indicating

that the embedding scheme is not over-parametrized.

The choice of a proper QM cluster depends on the individual materials system. In silica systems, char-

acterized by ideal nearest neighbor coordination, a cluster in which all QM host atoms are oxygens,

regardless of the cluster stoichiometry, is favorable. In contrast, alumina systems, characterized by an

overcoordinated bonding situation, only small deviations from Al2O3stoichiometry in the solid part of

the surface are tolerable. Under this constraint it is not possible to construct a cluster with only one ele-

ment type as QM host atoms and a reasonable surface to volume ratio. Constraints on individual atom’s

Mulliken populations, especially those of the QM host atoms, can help to reduce charge oscillations intro-

duced by the QM/MM boundary, but they do not solve the problems introduced by non-stoichiometric

alumina cluster composition.

Regarding formation between natively oxidized aluminum and an epoxy-type adhesive system, an insight

into the working mechanism of silane-based adhesion promoters could be achieved from examining the

minimum energy paths of the condensation reactions of the model adhesive’s components in vacuum. One

important effect of the adhesion promoter is an increase in the number of covalent substrate–polymer

links, due to the fact that condensation of such adhesion promoters is favorable at more surface sites than

for the other components.

Furthermore, first results on the influence of neighboring adsorbates during the interface formation,

indicate that the direct influence of other adhesive components in the vicinity of the adsorption site is

small, both in absolute terms and in comparison to the effect of water. Even most of the effect of the

neighboring adsorbates is caused by the addition of water molecules formed during the condensation

reactions leading to covalent surface–polymer bonds. These initial results, which only take into account

the influence of water formed during the interface formation suggest, that proton-transfer processes at the

82 CHAPTER 6. CONCLUSION

surface, during the adsorption of adhesive components can strongly influence the surface hydroxylation

and may lead to the passivation of functional groups of epoxy resins.

From this, the course for further work to understand the formation of the alumina–adhesive interface

becomes clear. The influence of water is of paramount importance and must be studied in a system-

atic manner. The extensive calculations necessary for this have been enabled by the development of a

computationally efficient QM/MM coupling scheme applicable the surfaces of polar solids. Parallel to

this, the knowledge about interface reaction energetics obtained during these studies at atomic resolution

can serve as input parameters for Kinetic Monte-Carlo simulations aimed at understanding the interface

formation process on larger length- and time scales.

Bibliography

[1] Lenny Bernstein, Peter Bosch, Osvaldo Canziani, Zhenlin Chen, Renate Christ, Ogunlade Davidson,

William Hare, Saleemul Huq, David Karoly, Vladimir Kattsov, Zbigniew Kundzewicz, Jian Liu,

Ulrike Lohmann, Martin Manning, Taroh Matsuno, Bettina Menne, Bert Metz, Monirul Mirza,

Neville Nicholls, Leonard Nurse, Rajendra Pachauri, Jean Palutikof, Martin Parry, Dahe Qin,

Nijavalli Ravindranath, Andy Reisinger, Jiawen Ren, Keywan Riahi, Cynthia Rosenzweig, Matilde

Rusticucci, Stephen Schneider, Youba Sokona, Susan Solomon, Peter Stott, Ronald Stouffer, Taishi

Sugiyama, Rob Swart, Dennis Tirpak, Coleen Vogel, Gary Yohe, Terry Barker, Abdelkader Allali,

Roxana Bojariu, Sandra Diaz, Ismail Elgizouli, Dave Griggs, David Hawkins, Olav Hohmeyer,

Bubu Pateh Jallow, Lucka Kajfez-Bogataj, Neil Leary, Hoesung Lee, and David Wratt. Climate

Change 2007: Synthesis Report. An Assessment of the Intergovernmental Panel on Climate Change.

Technical report, Intergovernmental Panel on Climate Change, 2007. 1

[2] Andrew R. Leach. Molecular Modelling. Pearson, Harlow, UK, 2nd edition, 2001. 2,1,3.2,3.3

[3] Robert G. Parr and Weitao Yang. Density-Functional Theory of Atoms and Molecules. Oxford

University Press, 1989.

[4] A Szabo and S. Ostlund. Modern Quantum Chemistry. Dover Publications, Inc., Mineola, New

York, 1989.

[5] Efthimios Kaxiras. Atomic and Electronic Structure of Solids. Cambridge University Press, Cam-

bridge, 2003. ISBN-10: 0521523397. 2

[6] D. R. Hartree. Proc. Camb. Phil. Soc., 24:89, 1926. 2.1

[7] H. Shull and G. G. Hall. Atomic units. Nature, 184:1559–1560, 1959. 2.1

[8] K. Ohno, K. Esfarjani, and Y. Kawazoe. Computational Materials Science. Springer series in solid

state sciences. Springer, 1999. 2.3

[9] M. Born and V. Fock. Beweis des Adiabatensatzes. Zeitschrift f¨

ur Physik A Hadrons and Nuclei,

51:165, 1928. 2.3

[10] D. J. Chadi and L. Cohen. Special Points in the Brillouin Zone. Phys. Rev. B, 8(12):5747–5753,

Dec 1973. 2.4.1

[11] Hendrik J. Monkhorst and James D. Pack. Special points for brillouin-zone integrations. Phys.

Rev. B, 13(12):5188–5192, Jun 1976. 2.4.1,4.3.2

[12] J.-M. Spaeth and H. Overhof. Point Defects in Semiconductors and Insulators – Determination of

Atomic and Electronic Structure from Paramagnetic Hyperfine Interactions, volume 51 of Springer

Sieries in Material Sciences. Springer Verlag, Heidelberg, 2003. 2.4.2

[13] H. Overhof, H. Weihrich, and G. Corradi. Electronic structure of isolated aluminum point defects

and associated trigonal pairs and clusters in Si. Phys. Rev. B, 45(16):9032–9041, Apr 1992.

84 BIBLIOGRAPHY

[14] H. Weihrich and H. Overhof. Ground-state properties of isolated interstitial iron in silicon: Elec-

tronic structure and hyperfine interactions. Phys. Rev. B, 54(7):4680–4695, Aug 1996. 2.4.2

[15] G. A. Baraff and M. Schl¨

uter. The LCAO approach to the embedding problem. J. Phys. C: Solid

State Physics, 19:4384–4391, 1986. 2.4.2

[16] R. P. Messmer and G. D. Watkins. Inst. Phys. Conf. Ser., 16:255, 1973. 2.4.2

[17] S. G. Louie, M. Schl¨

uter, J. R. Chelikowsky, and M. L. Cohen. Self-consistent electronic states for

reconstructed Si vacancy models. Phys. Rev. B, 13(4):1654..1663, Feb 1976. 2.4.2

[18] R. A. Evarestov and V. P. Smirnov. Use of the Large Unit Cell Approach for Generating Special

Points of the Brillouin Zone. Phys. Stat. Sol. (b), 99:463, 1980. 2.4.2

[19] M. J. Puska, S. P¨

oykk¨

o, M. Pesola, and R. N. Nieminen. Convergence of supercell calculations for

point defects in semiconductors: Vacancy in silicon. Phys. Rev. B, 58(3):1318–1325, 1998. 2.4.2

[20] P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Phys. Rev., 136:12520–12536, 1964. 2.5.1

[21] W. Kohn and L. J. Sham. Self-consistent equations including exchange and correlation effects.

Phys. Rev., 140:A1133, 1965. 2.5.1,2.5.1

[22] E. H. Lieb. Densitiy functionals for Coulomb systems. International Journal of Quantum Chemistry,

24:243, 1983. 2.5.1

[23] A. D. Becke. Density-functional exchange-energy approximation with correct asymptotic behavior.

Phys. Rev. A, 38(6):3098–3100, Sep 1988. 2.5.4

[24] Chengteh Lee, Weitao Yang, and Robert G. Parr. Development of the colle-salvetti correlation-

energy formula into a functional of the electron density. Phys. Rev. B, 37(2):785–789, Jan 1988.

2.5.4

[25] Axel D. Becke. Density-Functional thermochemistry. III. The role of exact exchange. J. Chem.

Phys., 98(7):5648–5652, Apr 1993. 2.5.4,2.5.6,2.5.6

[26] John P. Perdew, Kieron Burke, and Matthias Ernzerhof. Generalized gradient approximation made

simple. Phys. Rev. Lett., 77(18):3865–3868, Oct 1996. 2.5.4

[27] W. Kohn. Nobel Lecture: Electronic structure of matter — wave functions and density functionals.

Rev. Mod. Phys, 71:1253, 1999. 2.5.5

[28] N. F. Mott. Metal-insulator transitions. Taylor and Francis, London, 1974. 2.5.5

[29] K. Terakura, A.R. Williams, T. Oguchi, and J. K¨

ubler. Transition-Metal Monoxides: Band or Mott

Insulators. Phys. Rev. Lett., 52:1830, 1984. 2.5.5

[30] D. Vogel, P. Kr¨

uger, and J. Pollmann. Structual and electronic properties of group-III nitrides.

Phys. Rev. B, 55(19):12836–12839, May 1997. 2.5.6

[31] J. Neugebauer. private communication. 2004. 2.5.6

[32] R. Dovesi, R. Orlando, C. Roetti, C. Pisani, and V.R. Saunders. The Periodic Hartree-Fock Method

and Its Implementation in the CRYSTAL Code. phys. stat. sol. (b), 217:63–88, 2000. 2.5.6

[33] A. Gali, P. De´ak, E. Rauls, N. T. Son, I. G. Ivanov, F. H. C. Carlsson, E. Janz´en, and W. J. Choyke.

The correlation between the anti-site pair and the D1center in SiC. Phys. Rev. B, 67:155203/1–5,

2003. 2.5.6

[34] Peter De´ak, Adam Gali, Andras Soly´om, Adam Buruzs, and Frauenheim Thomas. Electronic struc-

ture of boron-interstitial clusters in silicon. Journal of Physics: Condensed Matter, 17(22):S2141–

S2153, 2005. 2.5.6

BIBLIOGRAPHY 85

[35] Walter A. Harrison. Pseudopotentials in the Theory of Metals. Frontiers in physics. Benjamin, 1966.

2.6

[36] Arthur M. Stoneham. Theroy of defects in solids. Monographs on the physics and chemistry of

materials. Clarendon Pr., 1985.

[37] G. B. Bachelet, D. R. Hamann, and M. Schl¨

uter. Pseudopotentials that work from H to Pu. Phys.

Rev. B, 26:4199, 1982.

[38] N. Troullier and J. L. Martins. Efficient pseudopotentials for plane-wave calculations. Phys. Rev.

B, 43(43):1993, Jan 1991. 2.6,2.6

[39] J. C. Barthelat and Serafini A. Durand, P. Non-empirical pseudopotentials for molecular calcula-

tions I. The PSIBMOL algorithm and test calculations. Molec. Phys., 33:159, 1977. 2.6

[40] D. Vanderbilt. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys.

Rev. B, 41(11):7892–7895, Apr 1990. 2.6

[41] K. Laasonen, A. Pasquarello, R. Car, C. Lee, and D. Vanderbilt. Car-Parrinello molecular dynamics

with Vanderbilt ultrasoft pseudopotentials. Phys. Rev. B, 47(16):10142–10153, Apr 1993.

[42] K. Laasonen, R. Car, C. Lee, and D. Vanderbilt. Implementation of ultrasoft pseudopotentials in

ab initio molecular dynamics. Phys. Rev. B, 43(8):6797–6799, Mar 1991. 2.6

[43] P. E. Bl¨

ochl. Projector augmented-wave method. Phys. Rev. B, 50(24):17953–17979, Dec 1994. 2.6

[44] W. Matthew C. Foulkes and Roger Haydock. Tight-binding models and density-functional theory.

Phys. Rev. B, 39(17):12520–12536, Jun 1989. 2.7.1

[45] G. Seifert, H. Eschrig, and W. Bieger. Eine approximative Variante des LCAO-XαVerfahrens. Z.

Phys. Chem., 267:529, 1986. 2.7.2

[46] D. Porezag. Development of ab-initio and approximate density functional methods and their ap-

plication to complex fullerene systems. PhD thesis, Fakult¨

at f¨

ur Naturwissenschaften, Technische

Universit¨

at Chemnitz-Zwickau, 1997.

[47] D. Porezag, Th. Frauenheim, T. K¨

ohler, G. Seifert, and R. Kaschner. Construction of tight-

binding-like potentials on the basis of density-functional theory: Application to carbon. Phys. Rev.

B, 51:12947, 1995. 2.7.2

[48] Th. Frauenheim, G. Seifert, M. Elstner, Z. Hajnal, G. Jungnickel, D. Porezag, S. Suhai, and

R. Scholz. A self-consistent charge density-functional based tight-binding method for predictive

materials simulations in physics, chemistry and biology. phys. stat. sol. (b), 217:41, 2000. 2.7.2,

2.7.2,2.7.2

[49] Michael Sternberg. The atomic structure of diamond surfaces and interfaces. PhD thesis, Fach-

bereich Physik, Universit¨

at Paderborn, 2001. 2.7.2

[50] Christof K¨

ohler. Ber¨

ucksichtigung von Spinpolarisationseffekten in einem dichtefunktionalbasierten

Ansatz. PhD thesis, Univerasit”at Paderborn, 2004. 2.7.2

[51] Jan M. Knaup, Ben Hourahine, and Thomas Frauenheim. Initial Steps Towards Automating the

Fitting of DFTB EREP .J. Phys. Chem A, 111:5637, 2007. 2.7.2

[52] Q. Cui, M. Elstner, E. Kaxiras, T. Frauenheim, and M. Karplus. A QM/MM implementation of

the self-consistent charge density functional tight binding (SCC-DFTB) method. J. Phys. Chem.

B, 105:569, 2001. 2.7.3,2.8.1

[53] B. Aradi and B. Hourahine. Atomic charge constrains in dftb. private communication, 2008. 2.7.3

86 BIBLIOGRAPHY

[54] S. Dapprich, I. Komaromi, K. S. Byun, K. Morokuma, and M. J. Frisch. A new oniom implementa-

tion in gaussian98. part i. the calculation of energies, gradients, vibrational frequencies and electric

field derivatives. J. Mol. Struct. (Theochem), 461:1, 1999. 2.8

[55] S. Humbel, S. Sieber, and K. Morokuma. The imomo method: Integration of different levels of

molecular orbital approximations for geometry optimization of large systems: Test for n-butane

conformation and sn2 reaction: Rcl+cl−.J. Chem. Phys., 105:1959, 1996. 2.8

[56] A. Warshel and M. Levitt. Theoretical studies of enzymic reactions: dielectric, electrostatic and

steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol., 103:227, 1976.

2.8,2.8.1

[57] M. J. Field, P. A. Bash, and M. Karplus. A combined quantum mechanical and molecular mechanical

potential for molecular dynamics simulations. J. Comput. Chem., 11:700, 1990. 2.8,2.8.1

[58] J. Aaqvist and A. Warshel. Simulation of enzyme reactions using valence bond force fields and

other hybrid quantum/classical approaches. Chem. Rev., 93:2523, 1993.

[59] J. Gao. Hybrid quantum mechanical/molecular mechanical simulations: An alternative avenue to

solvent effects in organic chemistry. Acc. Chem. Res., 29:298, 1996.

[60] M. Svensson, S. Humbel, R. D. J. Froese, T. Matsubara, S. Sieber, and K. Morokuma. ONIOM:

A multi-layered integrated MO + MM method for geometry optimizations and single point energy

predictions. a test for diels-alder reactions and Pt(P(t-Bu)3)2+ H2oxidative addition. J. Phys.

Chem., 100:19357, 1996. 2.8

[61] P.H. K¨

onig, M. Hoffmann, Th. Frauenheim, and Q. Cui. A critical evaluation of different qm/mm

frontier treatments using scc-dftb as the qm method. J. Phys. Chem. B, 109:9082, 2005. 2.8,2.8.1,

4.2.1

[62] D. Riccardi, P. Schaefer, Y. Yang, H. Yu, N. Ghosh, X. Prat-Resina, P. K¨

onig, G. Li, D. Xu,

H. Guo, M. Elstner, and Q. Cui. Development of effective quantum mechanical/molecular mechan-

ical (qm/mm) methods for complex biological processes. J. Phys. Chem. B, 110:6458–69., 2006.

2.8

[63] D. Das, K. Eurenius, E. Billings, P. Sherwood, D. Chatfield, M. Hodoscek, and B. Brooks. Optimiza-

tion of quantum mechanical molecular mechanical partitioning schemes: Gaussian delocalization of

molecular mechanical charges and the double link atom method. J.Chem.Phys., 117:10534, 2002.

2.8

[64] J. W. Ponder and D. A. Case. Force fields for protein simulations. Advances in protein chemistry,

66:27–85, 2003. 2.8

[65] C.R.A. Catlow, M. Dixon, and W. C. Mackrodt. In Computer Simulations of Solids, Lecture Notes

in Physics, volume 166, page 130. Springer Verlag, Berlin, 1982. 2.8

[66] J. Gao and D. G. Truhlar. Quantum mechanical methods for enzyme kinetics. Annu. Rev. Phys.

Chem., 53:467, 2002. 2.8

[67] A. Warshel. Computer simulations of enzyme catalysis: methods, progress, and insights. Annu.

Rev. Biophys. Biomol. Struct., 32:425, 2002.

[68] Q. Cui and M. Karplus. Catalysis and specificity in enzymes: a study of triosephosphate isomerase

and comparison with methyl glyoxal synthase. Adv. Prot. Chem., 66:315, 2003.

[69] C. Lennartz, A. Schaefer, F. Terstegen, and W. Thiel. Enzymatic reactions of triosephosphate

isomerase: A theoretical calibration study. J. Phys. Chem. B, 106:1758, 2002.

[70] V. Guallar and R. A. Friesner. Cytochrome p450cam. J. Am. Chem. Soc., 126:8501–8508, 2004.

BIBLIOGRAPHY 87

[71] G. A. Cisneros, H. Y. Liu, Y. K. Zhang, and W. T. Yang. Ab into qm/mm study. J. Am. Chem.

Soc., 125:10384, 2003.

[72] T.K. Woo, P.M. Margl, L. Deng, L. Cavallo, and T. Ziegler. Toward more realistic computational

modeling of homogenous catalysis by density functional theory: combined qm/mm and ab initio

molecular dynamics. Catalysis Today, 50:479, 1999.

[73] A. Maiti, M. Sierka, J. Andzelm, J. Golab, and J. Sauer. J. Phys. Chem. A, 104:10932, 2000.

[74] C.H. Choi and M.S. Gordon. J. Am. Chem. Soc., 121:11311, 1999.

[75] P.E. Sinclair, A. de Vries, P. Sherwood, C.R.A. Catlow, and R.A. van Santen. Quantum-chemical

studies of alkene chemisorption in chabazite: A comparison of cluster and embedded-cluster models.

J. Chem. Soc., Faraday Trans., 94:3401, 1998.

[76] P. Sherwood, A. H. de Vries, M. F. Guest, G. Schreckenbach, C. R. A. Catlow, S. A. French, A. A.

Sokol, S. T. Bromley, W. Thiel, A. J. Turner, S. Billeter, F. Terstegen, S. Thiel, J. Kendrick, S. C.

Rogers, J. Casci, M. Watson, F. King, E. Karlsen, M. Sjovoll, Fahmi A., A. Schafer, and C. Lennartz.

QUASI: A general purpose implementation of the QM/MM approach and its application to problems

in catalysis. THEOCHEM, 632:1, 2003.

[77] I.H. Hillier. J. Mol. Struct. (Theochem), 463:45, 1999.

[78] J. Sauer and M. Sierka. J. Comput. Chem., 21:1470, 2000.

[79] J.O. Noell and K. Morokuma. J. Phys. Chem., 80:2675, 1976.

[80] N. Lopez, G. Pacchioni, F. Maseras, and F. Illas. Chem. Phys., 294:611, 1998. 2.8

[81] D. Riccardi, G. Li, and Q. Cui. Importance of van der waals interactions in qm/mm simulations.

J. Phys. Chem. B, 108:6467, 2004. 2.8

[82] Y. Zhang, T. Lee, and W. Yang. A pseudobond approach to combining quantum mechanical and

molecular mechanical methods. J. Chem. Phys. B, 110:46, 1999. 2.8.1

[83] I. Antes and W. Thiel. Adjusted connection atoms for combined quantum mechanical and molecular

mechanical methods. J. Phys. Chem. A, 103:9290, 1999. 2.8.1

[84] G. G. Ferenczy, J. L. Rivail, P. R. Surjan, and G. Naray-Szabo. NDDO fragment self-consistent

field approximation for large electronic systems. J. Comput. Chem., 13:830, 1992. 2.8.1

[85] V. Thery, D. Rinaldi, J. L. Rivail, B. Maigret, and G. G. Ferenczy. Quantum-mechanical compu-

tations on very large molecular systems: the local self-consistent field method. J. Comput. Chem.,

15:269, 1994.

[86] D. M. Philipp and R. A. Friesner. Mixed ab initio qm/mm modeling using frozen orbitals and tests

with alanine dipeptide and tetrapeptide. J. Comput. Chem., 20:1468, 1999.

[87] R. B. Murphy, D. M. Philipp, and R. A. Friesner. A mixed quantum mechanics/molecular mechanics

(qm/mm) method for large-scale modeling of chemistry in protein environments. J. Comput. Chem.,

21:1442, 2000. 2.8.1

[88] J. Gao, P. Amara, C. Alhambra, and M. J. Field. A generalized hybrid orbital (GHO) method for

the treatment of boundary atoms in combined QM/MM calculations. J. Phys. Chem. A, 102:4714,

1998. 2.8.1

[89] J. Pu, J. Gao, and D. G. Truhlar. Combining self-consistent-charge density-functional tight-binding

(SCC-DFTB) with molecular mechanics by the generalized hybrid orbital (GHO) method. J. Phys.

Chem. A, 108:5454, 2004. 2.8.1

88 BIBLIOGRAPHY

[90] Jan M. Knaup, Christof K¨

ohler, Michael Hoffmann, Peter H. K¨

onig, and Th. Frauenheim. <em>Ab

initio</em>simulation of interface reactions as a foundation of understanding polymorphism. Eur.

Phys. J. Special Topics, 149:127, 2007. 2.8.1

[91] Klaus M¨

uller and Leo D. Brown. Location of saddle points and minimum energy paths by a

constrained simplex optimization procedure. Theor. Chim. Act., 53:75, 1979. 3.1,3.2

[92] Klaus M¨

uller. Reaction paths on multidimensional energy hypersurfaces. Angewandte Chemie

International Edition in English, 19:1, 1980. 3.1

[93] Baron Peters, Andreas Heyden, Alexis T. Bell, and Arup Chakraborty. A growing string method for

determining transition states: Comparison to the nudged elastic band and string methods. Journal

of Chemical Physics, 120(17):7877–7886, 2004. 3.1.1

[94] Graeme Henkelman, G´ısli J´ohannesson, and Hannes J´onsson. Theoretical Methods in Condensed

Phase Chemistry, volume 5 of Progress in Theoretical Chemistry and Physics, chapter Methods

for finding saddlepoints and minimum energy paths, pages 269–300. Kluwer Academic Publishers,

Dodrecht, the Netherlands, 2002. 3.1.1,3.1.2

[95] B. R. Gelin and M. Karplus. Sidechain torsional potentials and motion of amino acids in proteins.

bovine pancreatic trypsin inhibitor. Proc. Natl. Acad. Sci., 72:2002, 1975. 3.1.1

[96] J Baker. An algorithm for the location of transition states. J. Comput. Chem., 7:385, 1986. 3.1.1

[97] C. J. Cerjan and W. H. Miller. On finding transition states. J. Chem. Phys., 75:2800, 1981.

[98] D. T. Nguyen and D. A. Case. On finding stationary states on large-molecule potential energy

surfaces. J. Phys. Chem., 89:4020, 1985.

[99] W. Quapp. A gradient-only algorithm for tracing a reaction path uphill to the saddle of a potential

energy surface. Chem. Phys. Lett., 253:286, 1996.

[100] H. Taylor and J. Simmons. Imposition of geometrical constraints on potential energy surface walking

procedures. J. Phys. Chem., 89:684, 1985.

[101] D.J. Wales. Finding saddle points for clusters. J. Chem. Phys., 91:7002, 1989. 3.1.1

[102] M. J. S. Dewar, E. F. Healy, and J. J. P. Stewart. Location of transition states in reaction mecha-

nisms. J. Chem. Soc., Faraday Trans. 2, 80:227, 1984. 3.1.1

[103] V. Ionova, Irina and Emily A. Carter. Ridge method for finding saddle points on potential energy

surfaces. J. Chem. Phys., 98(8):6377–6386, April 1993. 3.1.1

[104] Graeme Henkelman and Hannes J´onsson. A dimer method for finding saddle points on high dimen-

sional potential surfaces using only first derivatives. J. Chem. Phys., 111:7010, 1999. 3.1.1

[105] S. Fischer and M. Karplus. Conjugate peak refinement: an algorithm for finding reaction paths and

accurate transition states in systems with many degrees of freedom. Chem. Phys. Lett., 194(3):252–

261, June 1992. 3.1.1

[106] Andreas Heyden, Alexis T. Bell, and Frerich J. Keil. Efficient methods for finding transition states

in chemical reactions: Comparison of improved dimer method and partitioned rational function

optimization method. The Journal of Chemical Physics, 123(22):224101, 2005. 3.1.1

[107] H. J´onsson, G. Mills, and K. W. Jacobsen. Classical and Quantum Dynamics in Condensed Phase

Simulations, chapter Nudged elastic band method for finding minimum energy paths of transitions,

pages 387–405. World Scientific, Singapur, 1998. 3.1.2

[108] Graeme Henkelman and Hannes J´onsson. Improved tangent estimate in the nudged elastic band

method for finding minimum energy paths and saddle points. J. Chem. Phys., 113(22):9978–9985,

December 2000. 3.1.2,5.4,2

BIBLIOGRAPHY 89

[109] P. Maragakis, Stefan A. Adreev, Yisroel Brumer, Davird R. Reichmann, and Efthimios Kaxiras.

Adaptive nudged elastic band approach for transition state calculation. Journal of Chemical Physics,

117(10):4651–4658, September 2002. 3.1.2

[110] Graeme Henkelman, Blas P. Uberuaga, and Hannes J´onsson. A climbing image nudged elastic

band method for finding saddle points and minimum energy paths. Journal of Chemical Physics,

113(22):9901–9904, December 2000. 3.1.2

[111] Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and

Edward Teller. Equation of state calculations by fast computing machines. The Journal of Chemical

Physics, 21(6):1087–1092, 1953. 3.2

[112] H. Nyquist. Certain topics in telegraph transmission theory. Trans. AIEE, 47:617, 1928. 2

[113] C. E. Shannon. Communication in the presence of noise. Proc. Institute of Radio Engineers, 37:10,

1949.

[114] E. T. Whittaker. On the functions which are represented by the expansions of the interpolation

theory. Proc. Roy. Soc. Edinburgh, 35:181, 1915. 2

[115] G. Milano and F. M¨

uller-Plathe. Mapping atomistic simulations to mesoscopic models: A sys-

tematic coarse-graining procedure for vinyl polymer chains. Journal of Physical Chemistry B,

109(39):18609–18619, 2005. 3.2

[116] V.A. Harmandaris, N.P. Adhikari, N.F.A. vanderVegt, and K. Kremer. Hierarchical modeling of

polystyrene: From atomistic to coarse-grained simulations. Macromolecules, 39(19):6708–6719,

2006.

[117] Berk Hess, Le´on Salvador, Nico van der Vegt, and Kurt Kremer. Long time atomistic polymer

trajectories from coarse grained simulations: bisphenol-a polycarbonate. Soft Matter, 2:409, 2006.

3.2

[118] Arthur F. Voter. Radiantion Effects in Solids, chapter Introduction to the Kinetic Monte Carlo

Method, page 1. NATO Science Series. Springer Netherlands, 2007. 3.3,3.3,3.3

[119] A. P. J. Jansen. An introduction to monte carlo simulations of surface reactions.

http://arxiv.org/abs/cond-mat/0303028v1, 2003. Lecture Notes. 3.3

[120] J. R. Beeler. Displacement spikes in cubic metals. i. α-iron, copper, and tungsten. Phys. Rev.,

150(2):470–487, Oct 1966. 3.3

[121] Richard Gordon. Adsorption isotherms of lattice gases by computer simulation. The Journal of

Chemical Physics, 48(3):1408–1409, 1968. 3.3

[122] Michael Bowker and David A. King. Adsorbate diffusion on single crystal surfaces : I.the influence

of lateral interactions. Surface Science, 71:583–598, 1978.

[123] David A. Reed and Gert Ehrlich. Surface diffusivity and the time correlation of concentration

fluctuations. Surface Science, 105:603–628, 1981. 3.3

[124] I. Martin-Bragado, P. Castrillo, M. Jaraiz, R. Pinacho, J. E. Rubio, J. Barbolla, and V. Moroz.

Fermi-level effects in semiconductor processing: A modeling scheme for atomistic kinetic monte

carlo simulators. Journal of Applied Physics, 98(5):053709, 2005. 3.3

[125] R. Marcelin. Annales de Physique, 3:120–231, 1915. 3.3

[126] Henry Eyring. The activated complex in chemical reactions. The Journal of Chemical Physics,

3(2):107–115, 1935. 3.3

90 BIBLIOGRAPHY

[127] M.H. Radke de Cuba, H. Emmerich, and S. Gemming. Finding polymorphic structures during

vicinal surface growth. Eur. Phys. J. Special Topics, 149:43–56, 2007. 3.3,5.8

[128] Volker Blobel and Erich Lohrmann. Statistische und numerische Methoden der Datenanalyse. B.

G. Teubner Stuttgart, 1998. 4.2.3

[129] A. Ramsat, G Brocks, and P. J. Kelly. Theoretical study of the si(001) surface reconstruction.

Phys. Rev. B, 51:14504, 1995. 4.3.2

[130] Christof K¨

ohler, Zolt´an Hajnal, P´eter De´ak, Thomas Frauenheim, and S´andor Suhai. Theoretical

investigation of carbon defects and diffusion in α-quartz. Phys. Rev. B, 64(8):085333, Aug 2001.

4.3.3

[131] Karl Wefers and Chanaka Misra. Oxides and hydroxides of aluminum. Technical report, Alcoa

Laboratories, 1987. 5

[132] A. F. Holleman and E. Wiberg. Lehrbuch der Anorganischen Chemie, chapter Das Aluminium,

pages 1081–1083. de Gruyter, Berlin, 1995. 5

[133] M.-H. Lee, Chi-Feng Cheng, Volker Heine, and Jacek Klinowski. Distribution of tereahedral and

octahedral Al sites in gamma alumina. Chem. Phys. Lett., 265:673–676, 1997. 5.1

[134] C. Wolverton and K. C. Hass. Phase stability and structure of spinel-based transition aluminas.

Phys. Rev. B, 63:024102, December 2000.

[135] H. Saalfeld and B. Mehrotra. Elektronenbeugungsuntersuchungen an aluminiumoxiden. Ber. Dtsch.

Keram. Ges., 42:161–166, 1965.

[136] T. C. Chou and T. G. Nieh. Nucleation and concurrent anomalous grain growth if α-al2o3during

γ-αphase transformation. J. Am. Ceram. Soc., 74:2270, 1991.

[137] B. C. Lippens and J. H. de Boer. Study of phase transformations during calcination of aluminum

hydroxides by selected area electron diffraction. Acta. Cryst., 17:1312–1321, 1964.

[138] O. Maresca, A. Ionescu, A. Allouche, J.P. Aycard, M. Rajzmann, and F. Hutschka. Quantum study

of the active sites of the γ-alumina surface (II): QM/MM (LSCF) approach to water, hydrogen

disulfide and carbon monoxide adsorption. J. Mol. Struc.-Theochem, 620:119–128, 2003. 5.1

[139] G. G. Paglia, A. L. Rohl, C. E. Buckley, and J. D. Gale. Determination of the structure of γ-alumina

from interatomic potential and first-principles calculations: The requirement of significnt numbers

of nonspinel positions to achieve an accurate structural model. Phys. Rev. B, 71:224115, 2005. 5.1

[140] G. Paglia, C. E. Buckley, A. L. Rohl, B. A. Hunter, R. D. Hart, J. V. Hanna, and L. T. Burne.

Tetragonal structure model for boehmite-derived γ-alumina. Phys. Rev. B, 68:144110, 2003. 5.1

[141] G. G. Paglia. Determination of the Structure of γ-Alumina using Empirical and First Principles

Calculations combined with Supporting Experiments. PhD thesis, Curtin University of Technology,

Perth, 2004. 5.1

[142] J. B. Peri. A Model of the Surface of γ-Alumina. J. Phys. Chem., 69:220–230, 1965. 5.1,5.5

[143] J. B. Peri. Infrared and Gravimetric Study of the Surface Hydration of γ-Alumina. J. Phys. Chem.,

69:211–219, 1965.

[144] J. H. De Boer, J. M. Fortuin, B. C. Lippens, and W. H. Meijs. Study of the nature of surfaces with

polar molecules ii. the adsorption of water on aluminas. J. Catal., 2:1–7, 1963.

[145] D. S. Maciver, H. H. Tobin, and R. T. Barth. Catalytic aluminas i. surface chemistry of eta and

gamma alumina. J. Catal., 2:485–497, 1963. 5.1,5.5

BIBLIOGRAPHY 91

[146] O. Maresca, A. Allouche, J. P. Aycard, M. Rajzmann, S. Clemendot, and F. Hutschka. Quantum

study of the active sites of the γ-alumina surface: chemisorption and adsorption of water, hydrogen

sulfide and carbon monoxide on aluminum and oxygen sites. J. Mol. Struc.-Theochem, 505:81 –94,

2000. 5.1

[147] David A. De Vito, Fran¸cois Gilardoni, Lioubov Kiwi-Minsker, Pierre-Yves Morgantini, St´ephane

Porchet, Albert Renken, and Jacques Jacques Weber. Theoretical investigation of the adsorption

of methanol on the (110) surface of γ-alumina. J. Mol. Struc.-Theochem, 469:7–14, 1999. 5.2

[148] V. Movarek, M. Kraus, L. V. Malysheva, E. A. Paukshtis, and E. N. Yurchenko. Czech Chem.

Commun., 53:459, 1988. 5.2

[149] B. Schneider, O.-D. Hennemann, and W. Possart. J. Adhes., 78:779, 2002. 5.2

[150] T. Kr¨

uger, M. Amkreutz, P. Schiffels, B. Schneider, O.-D. Hennemann, and Th. Th. Frauenheim.

Theoretical Study of the Interaction between Selected Adhesives and Oxide Surfaces. J. Phys.

Chem. B, 109:5060–5066, 2005. C.U.B. 5.2

[151] Jan M. Knaup, Christof K¨

ohler, Thomas Frauenheim, Alexander T. Blumenau, Marc Amkreutz,

Peter Schiffels, Bernhard Schneider, and Otto-Diedrich Hennemann. Computational studies on

polymer adhesion at the surface of γ-al2o3. i. the adsorption of adhesive component molecules from

the gas phase. J. Phys. Chem. B, 110:20460, 2006. 5.1,5.3.1,5.1,1

[152] Schneider B., P. Schiffels, and R. Wilken. Bonding of amine curing agents to native aluminum oxide

surfaces. In Proceedings of the 28th Annual Meeting of the Adhesion Society, page 431, 2005. 5.2,

5.5.1

[153] R. W. G. Wyckoff. Crystal Structures. Interscience, New York, 1965. 5.3.1

[154] Peter Schiffles. private communication. Fraunhofer IFAM, 2007. 5.4

[155] Russell D. Johnson III. NIST Computational Chemistry Comparison and Benchmark Database,

NIST Standard Reference Database number 101. http://srdata.nist.gov/cccbdb, Sept 2006. Release

14. 5.4

[156] W. Humphrey, A. Dalke, and K. Schulten. Vmd -visual molecular dynamics. J. Molecular Graphics,

14:33–38, 1996. A.2

[157] B. Aradi, B. Hourahine, and Th. Frauenheim. DFTB+, a sparse matrix-based implementation of

the DFTB method. J. Phys. Chem. A, 111:5678, 2007. A.2

92 BIBLIOGRAPHY

Appendix A

QM/MM Evaluation Data

A.1 Simple Link Atoms

-1

-0.8

-0.6

-0.4

-0.2

0.2

band shift

0 1 2 3 4 5 6

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 1 2 3 4 5 6

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.16

0.17

0.18

0.19

0.2

0.21

0.22

0.23

0.24

0.25

RMSQ

0 1 2 3 4 5 6

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

1.4

1.6

1.8

2.2

2.4

2.6

2.8

0 1 2 3 4 5 6

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.1: SLA evalution data for the 29+28 Si cluster.

94 APPENDIX A. QM/MM EVALUATION DATA

-3

-2

-1

band shift

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.2

0.4

0.6

0.8

1.2

1.4

1.6

1.8

RMSQ

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-14

-12

-10

-8

-6

-4

-2

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.2: SLA evalution data for the 15+12 α-quartz cluster. The circle marks the selected

embedding parameters for this system.

A.1. SIMPLE LINK ATOMS 95

-7

-6

-5

-4

-3

-2

-1

band shift

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

RMSQ

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-4

-2

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.3: SLA evalution data for the 27+12 α-quartz cluster. The circle marks the selected

embedding parameters for this system.

96 APPENDIX A. QM/MM EVALUATION DATA

-1

-0.8

-0.6

-0.4

-0.2

0.2

band shift

0 1 2 3 4 5 6

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 1 2 3 4 5 6

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.16

0.17

0.18

0.19

0.2

0.21

0.22

0.23

0.24

0.25

RMSQ

0 1 2 3 4 5 6

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

1.4

1.6

1.8

2.2

2.4

2.6

2.8

0 1 2 3 4 5 6

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.4: SLA evalution data for the 13+8 SSZ60 cluster. The circle marks the selected embedding

parameters for this system.

A.1. SIMPLE LINK ATOMS 97

-10

-5

band shift

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.5

1.5

2.5

RMSQ

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-30

-25

-20

-15

-10

-5

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.5: SLA evalution data for the 21+24 SSZ60 cluster. The circle marks the selected embed-

ding parameters for this system.

98 APPENDIX A. QM/MM EVALUATION DATA

-10

-8

-6

-4

-2

band shift

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

RMSQ

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-10

-5

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.6: SLA evalution data for the 45+24 SSZ60 cluster. The circle marks the selected embed-

ding parameters for this system.

A.1. SIMPLE LINK ATOMS 99

-3

-2

-1

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

-10

-9

-8

-7

-6

-5

-4

-3

-2

-1

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

Figure A.7: SLA evalution data for the 50+60 atom α-Al2O3cluster. The circle marks the selected

embedding parameters for this system.

100 APPENDIX A. QM/MM EVALUATION DATA

-4

-2

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

-6

-5

-4

-3

-2

-1

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

Figure A.8: SLA with charge constraints evalution data for the 50+60 atom α-Al2O3cluster. The

circle marks the selected embedding parameters for this system.

A.1. SIMPLE LINK ATOMS 101

-12

-10

-8

-6

-4

-2

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

100

110

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.2

0.4

0.6

0.8

1.2

1.4

1.6

1.8

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-35

-30

-25

-20

-15

-10

-5

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.9: SLA evalution data for the 71+102 atom α-Al2O3cluster. The circle marks the selected

embedding parameters for this system.

102 APPENDIX A. QM/MM EVALUATION DATA

-3

-2

-1

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

-10

-9

-8

-7

-6

-5

-4

-3

-2

-1

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

Figure A.10: SLA with charge constraints evalution data for the 71+102 atom α-Al2O3cluster.

The circle marks the selected embedding parameters for this system.

A.2. NEUTRALIZED CLUSTERS 103

A.2 Neutralized Clusters

-2

-1

band shift

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.2

0.4

0.6

0.8

1.2

1.4

1.6

1.8

2.2

RMSQ

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-18

-16

-14

-12

-10

-8

-6

-4

-2

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.11: HCS evalution data for the 15+12 α-quartz cluster. The circle marks the selected

embedding parameters for this system.

104 APPENDIX A. QM/MM EVALUATION DATA

-7

-6

-5

-4

-3

-2

-1

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-2

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.12: HCS evalution data for the 27+12 α-quartz cluster.

A.2. NEUTRALIZED CLUSTERS 105

-3

-2

-1

band shift

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.5

1.5

2.5

RMSQ

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-35

-30

-25

-20

-15

-10

-5

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.13: HCS evalution data for the 21+24 SSZ60 cluster. The circle marks the selected

embedding parameters for this system.

106 APPENDIX A. QM/MM EVALUATION DATA

-10

-9

-8

-7

-6

-5

-4

-3

-2

-1

band shift

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

RMSQ

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-2

0 1 2 3 4 5

gaussian blur width [Bohr]

0.4

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.14: HCS evalution data for the 45+24 SSZ60 cluster. The circle marks the selected

embedding parameters for this system.

A.2. NEUTRALIZED CLUSTERS 107

-12

-10

-8

-6

-4

-2

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

100

110

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

0.2

0.4

0.6

0.8

1.2

1.4

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

-25

-20

-15

-10

-5

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

Figure A.15: BCTC evalution data for the Al20O51+H102 α-Al2O3cluster. The circle marks the

selected embedding parameters for this system.

108 APPENDIX A. QM/MM EVALUATION DATA

-4

-3

-2

-1

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

0.05

0.1

0.15

0.2

0.25

0.3

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

-10

-9

-8

-7

-6

-5

-4

-3

-2

-1

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

Figure A.16: BCTC with charge constraints evalution data for the Al20O51+H102 α-Al2O3cluster.

The circle marks the selected embedding parameters for this system.

A.2. NEUTRALIZED CLUSTERS 109

-2

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

-12

-10

-8

-6

-4

-2

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

Figure A.17: BCTC evalution data for the Al20O30+H60 α-Al2O3cluster. The circle marks the

selected embedding parameters for this system.

110 APPENDIX A. QM/MM EVALUATION DATA

-4

-2

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

0.05

0.1

0.15

0.2

0.25

0.3

0.35

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

-6

-5

-4

-3

-2

-1

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.4

0.6

0.8

1.2

1.4

1.6

link-atom distance factor

Figure A.18: BCTC with charge constraints evalution data for the Al20O30+H60 α-Al2O3cluster.

The circle marks the selected embedding parameters for this system.

A.2. NEUTRALIZED CLUSTERS 111

-5

-4

-3

-2

-1

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.1

1.2

1.3

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

-35

-30

-25

-20

-15

-10

-5

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

Figure A.19: BCTC with data for the Al32O18+H81 α-Al2O3cluster. The circle marks the selected

embedding parameters for this system.

112 APPENDIX A. QM/MM EVALUATION DATA

-4

-3

-2

-1

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

0.05

0.1

0.15

0.2

0.25

0.3

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

-10

-9

-8

-7

-6

-5

-4

-3

-2

-1

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1.05

link-atom distance factor

Figure A.20: BCTC with charge constraints evalution data for the Al32O18+H81 α-Al2O3cluster.

The circle marks the selected embedding parameters for this system.

A.2. NEUTRALIZED CLUSTERS 113

-4

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.55

0.6

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-8

-7

-6

-5

-4

-3

-2

-1

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.21: BCTC evalution data for the Al24O40H6+H52 (site 1) γ-Al2O3cluster. The circle

marks the selected embedding parameters for this system.

114 APPENDIX A. QM/MM EVALUATION DATA

-2

-1.5

-1

-0.5

0.5

1.5

2.5

3.5

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0.5

1.5

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.22: BCTC with charge constraints evalution data for the Al24O40H6+H52 (site 1) γ-Al2O3

cluster. The circle marks the selected embedding parameters for this system.

A.2. NEUTRALIZED CLUSTERS 115

-3

-2

-1

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-20

-18

-16

-14

-12

-10

-8

-6

-4

-2

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.23: BCTC evalution data for the Al36O56H12+H67 (both sites) γ-Al2O3cluster. The

circle marks the selected embedding parameters for this system.

116 APPENDIX A. QM/MM EVALUATION DATA

-3

-2

-1

band shift

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

shifted range states

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

0.02

0.04

0.06

0.08

0.1

0.12

0.14

RMSQ

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

-8

-7

-6

-5

-4

-3

-2

-1

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

gaussian blur width [Bohr]

0.5

0.6

0.7

0.8

0.9

1.1

link-atom distance factor

Figure A.24: BCTC with charge constraints evalution data for the Al36O56H12+H67 (both sites)

γ-Al2O3cluster. The circle marks the selected embedding parameters for this system.

Acknowledgments

I am grateful to Professor Thomas Frauenheim for providing me with a challenging theme for my disserta-

tion. I also thank him for his supervision and advice and for providing a very good working environment

and facilities to conduct my research.

Funding of this work from the DFG priority program 1155 — Molecular Modeling and Simulation in

Process Engineering is greatly appreciated.

The non-QM/MM calculations in chapter 5 were carried out on the ARMINIUS cluster of the Paderborn

Center for Parallel Computing (PC2). I want to extend my gratitude for providing the computer time.

I thank my colleagues at the theoretical physics department in Paderborn and at the Bremen Center for

Computational Materials Science for creating a pleasant working atmosphere and for always being willing

to answer my questions.

I want to thank Dr. Peter De´ak for being a mentor and good friend. This work might never have been

finished without your advice and moral support.

My gratitude goes towards Drs. Michael Hoffmann and Peter K¨

onig for sharing their experience in

QM/MM coupling as well as in reaction path search.

I thank Drs. C. K¨

ohler, B. Aradi, A. Gagliardi, S. Sanna, B. Hourahine and M. Hoffmann for their support

and many fruitful discussions during the course of my research. I am also grateful to C. K¨

ohler for keeping

the computers and network alive, even under difficult conditions, and to B. Aradi for implementing my

requests for extensions to DFTB+.

Thanks to Brigitte Farchmin, Juliane Schoppe, Christine Frauenheim, Karin Sch¨

utte and Simone Lange

for their support, especially for keeping a lot of administrative work off my back and for dealing with

complicated travel expenses.

I am grateful to Professor Keiji Morokuma, Professor Stephan Irle and Professor Yoshiyuki Kawazoe for

their hospitality on my very educating and inspiring research visits to Kyoto, Nagoya and Sendai.

I would like to thank Dr. Adam Gali for close cooperation, great hospitality and for showing me the

vineyards and woods around lake Balaton.

The very interesting and educating cooperation with Drs. Sanaz Mostaghim, B´alint Aradi, Michael Hoff-

mann and Benjamin Hourahine on the application of evolutionary algorithms to DFTB parametrization

is greatly appreciated.

My gratitude goes towards Drs. Marc Amkreutz, Peter Schiffels and Bernd Schneider from the Fraunhofer

IFAM, for the cooperation and advice on the adsorption of adhesive components on alumina, even under

difficult conditions.

I am blessed by the friendship of Michael Hoffmann, who is a great companion, and a dependable friend.

Thank you for loads of advice and for listening to my elegies of graduate student misery. Without your

help, I could not have succeeded.

Thanks to Marius, B´alint, Alessio, Pia and Boris for great trip, full of stunning sights, art and adventure.

117

118 APPENDIX A. QM/MM EVALUATION DATA

I greatly appreciate Peule, Game Over, Sinoed, Spinor and Scofield for cover fire, tanking and amazing

amounts of dps.

I am indebted to Arno, B´alint, Christoph, Jens, Katrin and Simone for transportation services.

Special thanks go to Drs. B. Aradi, M. Hoffmann and B. Hourahine for proofreading the manuscript and

a lot of helpful advice; also thanks to Simone Lange for logistics support in printing and submitting my

thesis.

Finally, I want to express my gratitude to my parents, family and friends for their moral and practical

support during the course this work. Far too many people to list helped me during the past four years.

Your contribution is bigger than you probably think.

Colophon

Most graphical representations of molecular and crystal structures were created using VMD [156].

DFTB calculations were performed using DFTB+[157] version 1.0.1 and development versions.

119

120 COLOPHON

List of Figures

2.1 Covalent bond crossing the QM/MM border . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.1 Illustration of adiabatic mapping path-searches in a M¨

uller-Brown Potential . . . . . . . . 24

3.2 Illustration of an NEB path optimization in a M¨

uller-Brown Potential . . . . . . . . . . . 26

3.3 Illustration of a rare reactive event . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

4.1 Ranges of electronic and mechanical influences of small molecule adsorption . . . . . . . . 31

4.2 Bonding situation of a QM cluster in γ-Al2O3. . . . . . . . . . . . . . . . . . . . . . . . . 32

4.3 Illustration of the QM zone neutralization . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

4.4 Cross correlation between two densities of states and states in shifted gap . . . . . . . . . 38

4.5 DOS around the band gaps of the Si reference structures . . . . . . . . . . . . . . . . . . . 40

4.6 Structure and charge distribution of the 2x1 Si(001) surface . . . . . . . . . . . . . . . . . 40

4.7 Isosurface of the gap LDOS at the Si (001) 2x1 reconstructed surface . . . . . . . . . . . . 41

4.8 DOS around the band gaps of the SiO2reference structures . . . . . . . . . . . . . . . . . 42

4.9 Structure of the 1x1 reconstructed α-SiO2surface model and its electronic characteristics 42

4.10 Geometry and Mulliken charge distribution of SSZ60 Zeolite . . . . . . . . . . . . . . . . . 43

4.11 DOS around the band gaps of the Al2O3reference structures . . . . . . . . . . . . . . . . 44

4.12 Geometry and Mulliken charges of the α-Al2O3(0001) surface . . . . . . . . . . . . . . . . 44

4.13 Reference geometry and charges of the hyrdoxylated γ-Al2O3surface model . . . . . . . . 45

4.14 Embedding of the QM supercell into en extended distribution of external charges . . . . . 47

4.15 Charge distribution of the Si13 2x1 Si(001) surface . . . . . . . . . . . . . . . . . . . . . . 47

4.16 Charge errors in the 2x1 Si surface cluster over Gaussian blurring . . . . . . . . . . . . . . 48

4.17 Characteristics of the 29 atom Si surface cluster dependent on link-atom distance and

gaussian broadening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

4.18 States in shifted band gap of the tested α-quartz clusters . . . . . . . . . . . . . . . . . . 50

4.19 The tested γ-Al2O3clusters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4.20 Atomic orbital coefficients of the HOMO and LUMO of bulk α-quartz and α-Al2O3. . . 54

4.21 Si-(OH)-Si bridges formed by adsorption of an H2O molecule on α-quartz . . . . . . . . . 55

4.22 Calculated reaction energies for the dissociative H2O adsorption on α-quartz . . . . . . . 56

121

122 LIST OF FIGURES

4.23 Calculated reaction energies for the condensation of AMEO on γ-alumina . . . . . . . . . 56

5.1 Molecules of the model adhesive system . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

5.2 Views along the (001) and (010) directions of the hydroxylated γ-Al2O3surface model . . 62

5.3 Minimum Energy path of the ring-opening of t-DGEBA by H2O . . . . . . . . . . . . . . 63

5.4 Comparison of the calculated adsorption MEPs on the Al2O3surface. . . . . . . . . . . . 66

5.5 Minimum Energy path of the condensation of DETA at site 1 . . . . . . . . . . . . . . . . 67

5.6 Minimum Energy path of the condensation of DETA at site 2 . . . . . . . . . . . . . . . . 68

5.7 Minimum Energy path of the condensation of t-DGEBA at site 1 . . . . . . . . . . . . . . 69

5.8 Minimum Energy path of the condensation of t-DGEBA at site 2 . . . . . . . . . . . . . . 70

5.9 Minimum Energy path of the condensation of AMEO(hyd) at site 1 . . . . . . . . . . . . . 71

5.10 Minimum Energy path of the condensation of AMEO(hyd) at site 2 . . . . . . . . . . . . . 72

5.11 Minimum Energy path of the AMEO–AMEO co-adsorption . . . . . . . . . . . . . . . . . 73

5.12 Minimum Energy path of the DGEBA–AMEO co-adsorption . . . . . . . . . . . . . . . . 74

5.13 Co-adsorbed structure of t-DGEBA(open) at site 1 and DETA at site 2. . . . . . . . . . . . 76

5.14 Comparison of site 1 adsorption energy profiles in the vicinity of AMEO . . . . . . . . . . 77

5.15 DGEBA(open) passivated by a physisorbed water molecule . . . . . . . . . . . . . . . . . . 78

A.1 SLA evalution data for the 29+28 Si cluster . . . . . . . . . . . . . . . . . . . . . . . . . . 93

A.2 SLA evalution data for the 15+12 α-quartz cluster . . . . . . . . . . . . . . . . . . . . . . 94

A.3 SLA evalution data for the 27+12 α-quartz cluster . . . . . . . . . . . . . . . . . . . . . . 95

A.4 SLA evalution data for the 13+8 SSZ60 cluster . . . . . . . . . . . . . . . . . . . . . . . . 96

A.5 SLA evalution data for the 21+24 SSZ60 cluster . . . . . . . . . . . . . . . . . . . . . . . 97

A.6 SLA evalution data for the 45+24 SSZ60 cluster . . . . . . . . . . . . . . . . . . . . . . . 98

A.7 SLA evalution data for the 50+60 atom α-Al2O3cluster . . . . . . . . . . . . . . . . . . . 99

A.8 SLA with charge constraints evalution data for the 50+60 atom α-Al2O3cluster . . . . . 100

A.9 SLA evalution data for the 71+102 atom α-Al2O3cluster . . . . . . . . . . . . . . . . . . 101

A.10 SLA with charge constraints evalution data for the 71+102 atom α-Al2O3cluster . . . . . 102

A.11 HCS evalution data for the 15+12 α-quartz cluster . . . . . . . . . . . . . . . . . . . . . . 103

A.12 HCS evalution data for the 27+12 α-quartz cluster . . . . . . . . . . . . . . . . . . . . . . 104

A.13 HCS evalution data for the 21+24 SSZ60 cluster . . . . . . . . . . . . . . . . . . . . . . . 105

A.14 HCS evalution data for the 45+24 SSZ60 cluster . . . . . . . . . . . . . . . . . . . . . . . 106

A.15 BCTC evalution data for the Al20O51+H102 α-Al2O3cluster . . . . . . . . . . . . . . . . 107

A.16 BCTC with charge constraints evalution data for the Al20O51+H102 α-Al2O3cluster . . 108

A.17 BCTC evalution data for the Al20O30+H60 α-Al2O3cluster . . . . . . . . . . . . . . . . . 109

A.18 BCTC with charge constraints evalution data for the Al20O30+H60 α-Al2O3cluster . . . 110

A.19 BCTC evalution data for the Al32O18+H81 α-Al2O3cluster . . . . . . . . . . . . . . . . . 111

LIST OF FIGURES 123

A.20 BCTC with charge constraints evalution data for the Al32O18+H81 α-Al2O3cluster . . . 112

A.21 BCTC evalution data for the Al24O40H6+H52 (site 1) γ-Al2O3cluster . . . . . . . . . . . 113

A.22 BCTC with charge constraints evalution data for the Al24O40H6+H52 (site 1) γ-Al2O3

cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

A.23 BCTC evalution data for the Al36O56H12+H67 (both sites) γ-Al2O3cluster . . . . . . . . 115

A.24 BCTC with charge constraints evalution data for the Al36O56H12+H67 (both sites) γ-

Al2O3cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

124 LIST OF FIGURES

List of Tables

4.1 BCTC coefficients of the α-quartz model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.2 BCT coefficients of the SSZ60 model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.3 BCT coefficients of the α-Al2O2model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.4 BCTC coefficients of the γ-alumina model . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.5 Summary of results for SLA embedding . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4.6 Summary of results for neutralized cluster embedding . . . . . . . . . . . . . . . . . . . . 52

4.7 Charge distribution convergence check results for the Al24O40H6γ-alumina cluster . . . . 53

5.1 SCC-DFTB and LDA results on γ-Al2O3. . . . . . . . . . . . . . . . . . . . . . . . . . . 61

5.2 Reaction energies and barriers of the examined reactions in the model adhesive system on

Al2O3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.3 Influence of nearest neighbor adsorbates on adsorption energies. . . . . . . . . . . . . . . . 75

5.4 Influence of second neighbor adsorbates on adsorption energies. . . . . . . . . . . . . . . . 75

125