Lec Notes MEGHalle

See discussions, stats, and author profiles for this publication at: https://www.researchgate.
net/publication/251191601
10. Computational Materials Science with Materials Studio$^{\rlap{\LARGE\

(\circ\)}\,\hskip.5pt\raise3pt\hbox{\tiny\rm R}}\,$: Applications in Catalysis
Chapter · May 2004

DOI: 10.1007/978-3-540-39915-5_10
CITATIONS READS
6 915
5 authors, including:
Maria-Elena Grillo Niranjan Govind
41 PUBLICATIONS 364 CITATIONS
Pacific Northwest National Laboratory
158 PUBLICATIONS 5,924 CITATIONS
SEE PROFILE
SEE PROFILE
George Fitzgerald Klaus Stark

Universal Display Corporation Scienomics
80 PUBLICATIONS 3,653 CITATIONS 44 PUBLICATIONS 1,726 CITATIONS
SEE PROFILE SEE PROFILE
Some of the authors of this publication are also working on these related projects:
Natural gas conversion to aromatics over zeolite supported molybdenum catalysts View project
Advanced FMR and EPR studies on gold nanoclusters View project
All content following this page was uploaded by Maria-Elena Grillo on 06 September 2015.
The user has requested enhancement of the downloaded file.

Lecture Notes at the WE-Heraeus course, Martin-Luther-University Halle-Wittenberg, 2002
Preface
Computational modelling of novel materials is an increasingly powerful tool

being used in the development of advanced materials and their device appli-
cations. Computational materials science is a relatively new scientific field, in
which known concepts and recent advancements in physics, chemistry, math-
ematics and computer science are combined and applied numerically. The
unique advantage of such modelling lies in the possibility to predict macro-
scopic properties of materials based on calculations of microscopic quanti-
ties, i.e. at the atomic level. This has been made possible by the spectacular
increase in computational power over recent decades, allowing us to solve
numerically and with unprecedented accuracy, fundamental equations at the
atomic level. Today, based only on our knowledge of a single atom, we can
predict how the material formed by that atom type will look, what properties
that material will have and how it will behave under certain conditions. By
simply changing the arrangement of constituent atoms, or by adding atoms
of a different type, the macroscopic properties of all materials can be modi-
fied. It is in this way that one can learn how to improve mechanical, optical
and/or electronic properties of known materials, or one can predict proper-
ties of new materials, those which are not found in nature but are designed
and synthesized in the laboratory.
Supercomputers (both vector and parallel) and modern visualization tech-
niques are utilized to generate direct comparisons with experimental condi-
tions, and in some cases experiments may become redundant.
The authors of this book have endeavoured to give an overview of the
techniques, which operate at various levels of sophistication to describe mi-
croscopic and macroscopic properties of wide range of materials. The most
important methods used today in computational physics are addressed and,
in general, each topic is illustrated by a number of applications.
The book starts with basic aspects of density functional theory and the
discussion of modern methods to calculate the electronic structure of ma-
terials. A rapidly developing field of scientific interest over the last years is
nanophotonics. Two articles discuss how properties of photonic nanostruc-
tures can be computed.
The main part of the book contains contributions dealing with different
aspects of simulation methods. Ab initio calculations of free and supported
For personal use only

VI Preface
molecules and clusters are discussed. The application of molecular-dynamics

in biology, chemistry and physics is studied. The articles give a representative
cross section of different simulation methods on the one hand and of their
application to different materials on the other hand.
Essential for the field Computational Material Science is the availability
of effective algorithms and numerical methods. Therefore multigrid methods
and strategies for the implementation of sparse and irregular algorithms are
discussed as well.
The editors are grateful to the authors for their valuable contributions to
the book. The chapters are based, to some extend, on lectures given at the
WE-Heraeus course of the same name, held from 16th to 27th September
2002 in Halle. We gratefully acknowledge the support of the Wilhelm und
Else Heraeus Stiftung.
For the completion of the present book, the encouragement and patience
of Mr. M. Däne was essential.
Halle/Saale, Wolfram Hergert

October 2003 Arthur Ernst

Contents
1 Introduction
............................................................... 1
Part I Basic Description of Electrons and Photons in Crystals
Part II Simulation from Nanoscopic Systems to Macroscopic

Materials
2 Computational Materials Science with Materials Studio:

Applications in Catalysis
M. E. Grillo, J. W. Andzelm, N. Govind, G. Fitzgerald, K. B. Stark . . 9
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Geometry Optimization in Delocalised Internal Coordinates . . . . . . 10
2.3 Transition State Searching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Transition State Confirmation Algorithm . . . . . . . . . . . . . . . . . . . . . . . 17
2.5 Chemical Bonding and Elastic Properties of Corundum-Type
Oxides: The Rhodium Oxide Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3 Integration of Modelling
at Various Length and Time Scales
S. McGrother, G. Goldbeck-Wood, Yeng Ming Lam . . . . . . . . . . . . . . . . . . 25
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Structure-Activity and Structure-Property Approaches . . . . . . . . . . . 26
3.3 Atomistic and Mesoscale Simulations and their Parameterisation . . 27
3.3.1 Atomistic Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.2 Mesoscale Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3.3 Applications of Mesoscale Modeling . . . . . . . . . . . . . . . . . . . . . . . 31
3.4 Multiscale Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.4.1 From the Molecular to the Mesoscale . . . . . . . . . . . . . . . . . . . . . . 33
3.4.2 From Mesoscale to Finite Element Simulation . . . . . . . . . . . . . . 33
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

VIII Contents
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
Part III Modern Methods of Scientific Computing
Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

List of Contributors
Prof. Dr. H. Eschrig D. Hermann

Leibniz Institute for Solide State Institut für Theorie der Konden-
and Materials Research Dresden sierten Materie,
PF 27 01 16 Universität Karlsruhe
D-01171 Dresden D-76128 Karlsruhe
[email protected] [email protected]
Dr. A. Ernst Dr. S.F. Mingaleev

Max Planck Institute of Microstruc- Institut für Theorie der Konden-
ture Physics sierten Materie,
Weinberg 2 Universität Karlsruhe
D-06120 Halle D-76128 Karlsruhe
Dr. M. Lüders
M. Schillinger
Daresbury Laboratory,
Institut für Theorie der Konden-
Daresbury, Warrington,
sierten Materie,
Cheshire.
Universität Karlsruhe
WA4 4AD. UK
D-76128 Karlsruhe
[email protected]
[email protected]
Prof. Dr. K. Busch
Department of Physics and School Dr. A. Klaedtke
of Optics: CREOL & FPCE, Advanced Technology Institute
University of Central Florida, School of Electronics and Physical
Orlando, FL 32816, Sciences
U.S.A. University of Surrey
[email protected] GU2 7XH
United Kingdom
F. Hagmann [email protected]
Institut für Theorie der Konden-
sierten Materie, Dr. J. Hamm
Universität Karlsruhe Advanced Technology Institute
D-76128 Karlsruhe School of Electronics and Physical
Sciences
[email protected]

X List of Contributors
University of Surrey Dr. M. Celino

GU2 7XH Ente per le Nuove Tecnologie,
United Kingdom l’Energia e l’Ambiente,
[email protected] C.R. Casaccia, CP 2400, I-00100
Roma, Italy
Prof. Dr. O. Hess and Istituto Nazionale per la Fisica
Advanced Technology Institute della Materia,
School of Electronics and Physical Unità di Ricerca Roma1, Italy
Sciences [email protected]
University of Surrey
GU2 7XH
United Kingdom Dr. Y. Pouillon
[email protected] Ente per le Nuove Tecnologie,
l’Energia e l’Ambiente,
Dr. D. Ködderitzsch C.R. Casaccia, CP 2400, I-00100
Martin-Luther-Universität Halle- Roma, Italy
Wittenberg and Istituto Nazionale per la Fisica
Fachbereich Physik, Fachgruppe della Materia,
Theoretische Physik Unità di Ricerca Roma1, Italy
Von-Seckendorff-Platz 1 [email protected]
[email protected]
Dr. I.M.L. Billas
Prof. Dr. W. Hergert
Martin-Luther-Universität Halle- Département de Biologie et de
Wittenberg Génomique Structurales,
Fachbereich Physik, Fachgruppe Institut de Génetique et de Biologie
Theoretische Physik Moléculaire et Cellulaire,
Von-Seckendorff-Platz 1 1 rue Laurent Fries, BP 10142
D-06120 Halle F-67404 Illkirch, France
M. Däne Prof. Dr. V. Stepanyuk

Martin-Luther-Universität Halle- Max Planck Institute of Microstruc-
Wittenberg ture Physics
Fachbereich Physik, Fachgruppe Weinberg 2
Theoretische Physik D-06120 Halle
Von-Seckendorff-Platz 1 [email protected]
D-06120 Halle
[email protected]
Prof. Dr. P. Entel
Dr. C. Massobrio Universität Duisburg-Essen
Institut de Physique et Chimie des Institut für Physik, Theoretische
Matériaux de Strasbourg Tieftemperaturphysik
23 Rue du Loess BP 43 Duisburg Campus, Lotharstraße 1
F-67034 Strasbourg Cedex 2, France D-47048 Duisburg

List of Contributors XI
W. A. Adeagbo D-47048 Duisburg

Universität Duisburg-Essen [email protected]
Institut für Physik, Theoretische
Tieftemperaturphysik Dr. M. E. Grillo
Duisburg Campus, Lotharstraße 1 Accelrys, GmbH
D-47048 Duisburg Inselkammerstr. 1
[email protected] D-82008 Unterhaching
Germany
Dr. M.Sugihara [email protected]
Universität Duisburg-Essen
Institut für Physik, Theoretische Dr. J. W. Andzelm
Tieftemperaturphysik Accelrys Inc.
Duisburg Campus, Lotharstraße 1 9685 Scranton Rd
D-47048 Duisburg San Diego, CA 92121
[email protected] USA
[email protected]
G. Rollmann
Universität Duisburg-Essen
Dr. N. Govind
Institut für Physik, Theoretische
Accelrys Inc.
Tieftemperaturphysik
9685 Scranton Rd
Duisburg Campus, Lotharstraße 1
San Diego, CA 92121
D-47048 Duisburg
USA
[email protected]
[email protected]
A. T. Zayak
Universität Duisburg-Essen Dr. G. Fitzgerald
Institut für Physik, Theoretische Accelrys Inc.
Tieftemperaturphysik 9685 Scranton Rd
Duisburg Campus, Lotharstraße 1 San Diego, CA 92121
D-47048 Duisburg USA
M. Kreth Dr. K. B. Stark

Universität Duisburg-Essen Accelrys Inc.
Institut für Physik, Theoretische 9685 Scranton Rd
Tieftemperaturphysik San Diego, CA 92121
Duisburg Campus, Lotharstraße 1 USA
D-47048 Duisburg [email protected]
[email protected]
Dr. S. McGrother
K. Kadau Accelrys Inc.
Universität Duisburg-Essen 9685 Scranton Rd
Institut für Physik, Theoretische San Diego, CA 92121
Tieftemperaturphysik USA
Duisburg Campus, Lotharstraße 1 [email protected]

XII List of Contributors
Dr. G. Goldbeck-Wood [email protected]

Accelrys Inc
334 Science Park Prof. dr. G. Rünger
Cambridge CB4 0WN, UK Technische Universität Chemnitz
[email protected] Fakultät fr Informatik
Professur Praktische Informatik
Prof. Dr.-Ing. H. Altenbach Strae der Nationen 62
Martin-Luther-Universität Halle- D-09107 Chemnitz
Wittenberg [email protected]
Fachbereich Ingenieurwissenschaften
D-06099 Halle Prof. Dr. G. Wittum
[email protected] Interdisziplinäres Zentrum fr
Wissenschaftliches Rechnen (IWR)
Prof. Dr. Th. Rauber Simulation in Technology Center
Universität Bayreuth Im Neuenheimer Feld 368
Fakultt für Mathematik und Physik D-69120 Heidelberg
D-95440 Bayreuth [email protected]

1 Introduction
The longstanding goal of computer simulations of materials is the accurate

knowledge of structural and electronic properties for realistic condensed-
matter systems. Since materials are complex in nature, in the past theoreti-
cal investigations were restricted only on simple models. Recent progress in
computer technology made possible to afford a realistic description of a wide
range of materials. Since then the computational material science is rapidly
developing and expanding onto new fields in science and industry and the
computer simulation becomes an important part in many sciences and mod-
ern technologies such as physics, chemistry, biology and nanotechnology. This
book focuses on the background and practical aspects of the computational
physics.
Most of the methods in the computational physics are based on the density
functional theory (DFT). The density functional theory deals with inhomo-
geneous systems of identical particles. This approach is formulated on the
basic theorem of Hohenberg and Kohn which ensures that the ground state
properties of a many-particle system can be exactly represented in terms of
the ground state density. This makes possible to eliminate the monstrous
many-particle wave function from the formulation of the theory and to ex-
press chosen quantities of the system directly in terms of the particle den-
sity or the particle current density. The desired ground state quantities can
be obtained by minimization of corresponding unique functionals which are
in practice usually approximated by some model density functionals. One
of the most popular functionals is the Kohn-Sham local density approxi-
mation (LDA), in which all many-body effects are included on the level of
the homogeneous electron gas in the local exchange-correlation potential. In
this case the variational problem can be reduced to an effective Kohn-Sham
equation for single-particle orbitals representing uniquely the ground state
density. Such approach enables to carry out so-called “first-principle” or “ab-
initio” calculations (direct calculations of material properties from fundamen-
tal quantum mechanical theory). Many fundamental properties, for example
bond strength and reaction energies, can be estimated from first principles.
By definition the density functional theory is designed for the ground state
and fails to describe excited state properties. A classical tool for this problem
is the Green’s function formalism which is an important tool for studying the

correlations in many-body systems. The Green’s function provides spectral

densities for occupied and unoccupied states and its poles can be interpreted
as single-particle excitations. To calculate the Green’s function is in general a
non-trivial task and it requires to solve a the Dyson equation which involves
a non-interacting Green’s function and a self-energy operator. In practice it is
necessary to make approximations for the self-energy which can be handled.
In most first-principle methods the self-energy is approximated by an local
exchange-correlation potential. In this case the Kohn-Sham eigenvalues are
interpreted as excitation energies. This simple approach works surprisingly
well in many applications and is widely used in a variety of methods on the
first-principle level. However, the electron correlations and many-body effects
are very difficult issues for the Kohn-Sham approximation which fails to de-
scribe, for example, band gaps in semiconductors or life time in spectroscopy.
On ab-initio level the self-energy can be implemented within the random-
phase approximation (RPA). This method treads the electron correlations
on the basis of many-particle theory and is much more accurate but much
more time-consuming as the Kohn-Sham approach.
A development of first-principle methods is a difficult and challenging
task. Firstly, an “ab-inito” method should be constructed so general as pos-
sible to be applicable to a wide class of problems. Usually, a first principle
code contains of many thousand lines. Secondly, a development of such pro-
grams requires a deep knowledge of numerical methods and programming
tools. Because the majority of physical systems exhibit intrinsic symmetries,
one should use a symmetry analysis and the group theory to optimize and
to speed up computational process. The group theory as a mathematical
tool plays an important role to classify the solutions within the context of
the underlying symmetries. Extensive use of group theory has been made to
simplify the study of electronic structure or vibrational modes of solids or
molecules. The group theory is as well especially useful in understanding the
degeneracy of electronic energy levels.
First-principle methods enable to describe fundamental processes in bi-
ology, chemistry and physics as accurately as possible with moderate com-
putational effort. Nevertheless, the application of first-principle methods is
limited to study the evolution of real system in time. This can be done by
using first-principle molecular-dynamics methods (MD), which can simulate
the evolution of atomic and electronic motions without assuming empirical
parameters. Molecular dynamics has proved to be an optimum numerical
recipe applicable to problems with many degrees of freedom from quite dif-
ferent fields of science. The knowledge of the energy or potential landscape
of interacting particles, like electrons and atoms, enables one to calculate
the forces acting on the particles and to study the evolution of the system
with time. Overall, first-principles molecular dynamics appears as a convinc-
ing method to corroborate experimental work and make reliable predictions
based on well-established electronic structure techniques.

1 Introduction 3
Because of similarities between the Schrödinger equation and the Maxwell’s

equations, the concept of electronic structure methods can be adopted for
the study of photonic nano-materials, which have important applications in
micro- and nano-electronics. Such computer simulations allow now to con-
ceive quantum dots or photonic crystals that act as new source of a coherent
radiation, cages or guided pathways, operating on smaller and smaller scales
at steadily increasing speed.

Part II
Simulation from Nanoscopic Systems to

Macroscopic Materials

2 Computational Materials Science with

Materials Studio: Applications in Catalysis
M. E. Grillo1 , J. W. Andzelm2 , N. Govind2 , G. Fitzgerald2 , and

K. B. Stark2
1
Accelrys GmbH, Inselkammerstr. 1, D-82001 Unterhaching, Germany
2
Accelrys Inc, 9685 Scranton Rd, San Diego, CA, 92121, USA
Summary. In this article, developments in the primary ab initio quantum me-

chanics codes (CASTEP and DMol3 ) to solve critical problems in materials design
and catalysis research are reviewed. A novel, general-purpose internal coordinate
optimization scheme in DMol3 for periodic systems will be presented. Applications
of this robust optimizer to a full range of solid-state systems will be discussed. Fur-
thermore, the implementation of a transition state confirmation algorithm based
on the nudged elastic band method to validate a transition state by connecting it
to the proper reactant and product is explained.
Understanding and controlling the chemical and physical processes behind

materials design (e.g., chemical vapor deposition (CVD), corrosion, band-gap
engineering in semiconductors) and reactions in heterogeneous catalysis is the
principal goal of atomistic modeling and simulation technologies. Due to the
impressive methodological developments, density-functional theory (DFT) is
accepted in corporate, government and academic research labs as a powerful
tool to calculate free energies, as well as the electronic and atomistic structure
of medium-sized systems with predictive accuracy. Accelrys has focused a
great deal of R&D resources over the past few years to provide fast, robust
and accurate DFT codes for handling these systems and their processes.
In this article, developments in the primary ab initio quantum mechanics
codes (CASTEP and DMol3 ) to solve critical problems in materials design
and catalysis research are reviewed. Atomic-level knowledge of structures is
crucial to perform predictive simulations, be it of a carbon nanotube tip in a
field emission based flat-panel display, or of surface defects in a partial oxida-
tion catalyst. Therefore, geometry optimization of the structure of molecules
and crystals is a crucial task in simulations using quantum mechanics codes. A
novel, general-purpose internal coordinate optimization scheme in DMol3 for
periodic systems will be presented. Applications of this robust optimizer to a
full range of solid-state systems will be discussed. The energies of intermedi-
ates and transition states are also critical for catalytic reactions engineering.
For reactions at surfaces the dimension of phase space is so high that there
may exist more than one transition state. The transition state (TS) search
technique implemented in DMol3 , which does not require explicit calculation
of the Hessian matrix, is presented and validation examples are discussed.
Furthermore, the implementation of a transition state confirmation algorithm

10 M. E. Grillo et al.
based on the nudged elastic band method to validate a transition state by

connecting it to the proper reactant and product is explained. With input of
a reactant, a product and TS, it returns a trajectory containing any alterna-
tive minima in the pathway containing these three structures. This approach
can be extended to investigate potential energy surfaces containing several
minima and maxima for both molecules and solids.
2.1 Introduction
A catalyst is defined as a substrate capable to change the kinetics of a chem-
ical reaction increasing its rate by lowering the energy barrier to activate
the reactant molecules [1]. Catalysts are not consumed during a chemical
reaction, and their activity is defined as the reaction rate for conversion of
reactants into products. Understanding at a molecular level the relation be-
tween the catalyst surface composition and its activity, assists in the catalyst
design in industrial chemical processes.
Fig. 2.1. Schematic representation of energy profile for a catalyzed and uncat-
alyzed reaction.
2.2 Geometry Optimization

in Delocalised Internal Coordinates
Structural optimization of materials is a pre-requisite in most studies using

atomistic simulation. It is well established by now that internal coordinates
are more efficient for optimization of molecular systems than are Cartesian
coordinates. In the case of periodic systems the geometry optimization is

2 Materials Studio: Applications in Catalysis 11
typically done using Cartesian coordinates. Only recently, the first use of de-
localized internal coordinates for periodic systems was reported [2, 3]. This
is because use of internal coordinates for solids is more challenging than for
molecules. First, one needs to define a set of necessary primitive internal coor-
dinates such as bonds, angles, torsions that span the entire space of the solid.
For an infinite crystal, the number of internal coordinates is, in principle,
infinite; however, the number of degrees of freedom required for optimization
is just 3N − 3, where N is the number of the atoms in a single unit cell.
Therefore, only a subset of all internal coordinates of the crystal can be cho-
sen for geometry optimization.
According to Baker et al. [4], the best combination of internal coordi-

nates that span optimization space for molecules can be constructed by using
so-called delocalized coordinate method. We have adopted that method for
solid-state calculations. However, the number of primitive internal coordi-
nates is significantly larger for solids, particularly for dense packed systems
such as face centered cubic or hexagonal close packed metals. That requires
an efficient method of selecting the minimal number of internal coordinates
that span optimization space for solids. Such a ”pruning algorithm” [5] was
implemented recently in DMol3 [6]. We will briefly present our method here,
referring the interested reader to Ref. [2, 5] for a complete description.
The connection between internal and Cartesian displacements is defined

with the usual B matrix [7]. For a molecular systems B is defined through,
∂qα
Bαi = (2.1)
∂xi
where qα are the primitive internal coordinates and xi are the 3N Cartesian
coordinates of the atoms. For an infinite periodic crystal there are a finite
number of translationally unique primitive internals coordinates, 1 ≤ α ≤ M ,
which form a suitable set of coordinates to describe the system geometry. The
generalization of B for a periodic system is then defined [2] as
X ∂qα
Bαi = (2.2)
∂xiR
R
where the index α runs over a translationally unique set of primitive internal
coordinates, i runs over all 3N Cartesian coordinates of the atoms in one unit
cell and the sum over R runs over all unit cells in the crystal. B is M × 3N
dimensional giving the rate of change of primitive internal coordinate qα on
periodic displacement of the i-th atom in the unit cell. The infinite sum over
atoms in Equation 2 does not pose any problems in practice, because each
primitive internal depends at most on the position of four atoms in the crys-
tal. Hence, the terms in sum over R vanish after only a few cells.

Bakers scheme of construction of the best delocalized internal coordinates

requires diagonalization of the M × M dimensional matrix G=BBT to ob-
tain the eigenvectors U with non-zero eigenvalues. However, it is convenient
to introduce a closely related 3N × 3N matrix F=BT B [2, 5]. U can be
obtained directly from the eigenvectors of F. If FT=TΛ, with Λ being the
diagonal matrix of eigenvalues and T the matrix of eigenvectors with non-zero
eigenvalues, then the eigenvectors U of the matrix G are given by
1
U = BTΛ− 2 (2.3)
We can therefore obtain the matrix U that is needed to construct the de-
localized coordinates by diagonalizing the matrix F, which is much smaller
than G, since for solid-state systems M can be much greater than 3N . This
procedure will yield 3N − 3 vectors in U, which define our non-redundant
delocalized coordinates. The non-redundant delocalized internal coordinates
cnp of the system can be written as a linear combination of primitive internals
through X
cnp = qα Uαp (2.4)
α
An optimization procedure carried out in the space of non-redundant inter-

nal coordinates converges much faster than optimization in Cartesian coor-
dinates.
The natural bonding topology of some crystals yields a set of coordinates

that do not span the space of F, and optimization in such an internal coor-
dinate space would potentially fail to relax all degrees of freedom of the sys-
tem. Examples where this often occurs include graphite, molecular crystals,
molecules weakly bound to a surface and crystalline polymers. One solution
to this problem is to add additional bonds that connect the disconnected
fragments in the unit cell. However, this process is arbitrary and requires
user intervention.
We proposed an alternative, completely automatic solution to this prob-
lem. It can be proven that the primitive internal coordinates span the opti-
mization space if and only if F has precisely three zero eigenvalues [6]. If there
are more than three zero eigenvalues, supplementary coordinates have to be
added. These supplementary coordinates are constructed from the eigenvec-
tor space of F corresponding to the zero eigenvalues after three translational
vectors were removed from the space [5].
The performance of the new delocalized coordinate optimizer versus the

Cartesian based optimizer in DMol3 was evaluated for many solids [2, 5].
The computational overhead associated with the optimizer is very small, so
the speedup provided by the method can be measured as the number of
Cartesian optimization cycles divided by the number of internal coordinate
optimization cycles. Typically we find that delocalized internal coordinates

Fig. 2.2. Crystal of graphite. The internal coordinates of this crystal do not span
the optimization space
are 3-5 times more efficient than Cartesian coordinates. The efficiency of
internal coordinates increases with system size, and they are particularly
beneficial in the study of systems that can undergo significant rearrangement
of geometry. Periodic systems that belong to that category include surface
reactions, molecular crystals and zeolites. For example, our recent study on
Al-substituted chabazite revealed that optimization in delocalized internals
converges in 21 cycles, whereas Cartesian coordinates require 96 cycles [2].
The need for efficient geometry optimization in studying catalytic reac-

tions is illustrated in the DFT-study of the metal-catalyzed ring opening
of maleic anhydride (e.g., selective hydrogenation) [8]. This reaction is im-
portant for the commercial DuPont process of spandex fibers production. A
detailed understanding of the reaction pathway at a molecular level is ex-
pected to provide valuable insights into factors controlling catalytic activity
and selectivity.
Two adsorption modes for maleic-anhydride adsorption on Pd(111) were

studied: the atop (η 1 ) configuration, where the molecule is perpendicular to
the surface (Figure 3a); and the di-σ configuration (Figure 3b), where the
molecule is nearly parallel to the Pd surface.
DFT periodic slab calculations show that the di-σ adsorption configu-
ration (with the ring bound through the olefinic C=C group) is preferred
over the atop (η 1 ) adsorption mode. The calculated adsorption energy of 76
kJ/mol is in good agreement to the measured energy of 90 kJ/mol determined
from ultra-high vacuum temperature programmed desorption (UHV-TPD)
[9]. The di-σ adsorption configuration was found starting from the slightly
distorted η 1 mode. Optimization in delocalized internals took 34 cycles to
obtain the di-σ structure, whereas optimization using Cartesian coordinates

Fig. 2.3. η 1 (a) and di-σ (b) adsorption modes of maleic anhydride on Pd(111).
requires 141 cycles to reach the same structure.
Fig. 2.4. Reactant di-σ (a), transition state (b) and the stable intermediate (c) for
the C-O bond breaking in maleic anhydride on Pd(111) surface.
Figure 2.4 shows the reactant di-σ (a), transition state (b), and stable in-
termediate (c) structures for this reaction. The barrier of 59 kJ/mol is much
smaller than for the ring opening of the (η 1 ) structure. The ring opening of
di-σ structure is also an endothermic reaction, and leads to a stable interme-

diate structure of adsorption energy of about -39 kJ/mol.
2.3 Transition State Searching
The transition state (TS) search algorithm [10] implemented in our DFT
codes CASTEP and DMol3 is a generalized scheme based on the traditional
linear and quadratic synchronous transit (LST/QST) method [10] coupled
with previously-proposed conjugate gradient refinement ideas [11]. We refer
the interested reader to Ref. 10 and references therein for a complete descrip-
tion of the method.
The method can be summarized with a flow diagram (Figure 2.5). The
energies of suitable reactant and product structures are first calculated. This
is followed by a search for the LST maximum using essentially the method
of Halgren and Lipscpomb [10]. This maximum energy structure provides an
upper bound for the transition state. A conjugate gradient (CG) refinement
is initiated for further refinement of the estimated transition state, searching
in directions conjugate to the vector connecting reactants to products. CG
methods make intelligent use of gradient information, thereby requiring no
explicit Hessian to be calculated. As in a conventional optimization, the cal-
culation is considered converged if all the residual forces on the structure fall
below a certain threshold.
Fig. 2.5. Flow diagram illustrating the transition state search scheme.

If this is not achieved before the number of conjugate directions has been
explored, a QST maximum search is performed using the reactant, product,
and latest CG geometry. A new CG refinement cycle is started following the
QST maximum search. This is continued until the required convergence is
attained. A frequency analysis on this converged structure should yield ex-
actly one negative eigenmode, corresponding to the direction by which the
system would evolve away from the saddle point.
A number of calculations [12, 13, 14, 15] have been performed using the
new transition state search algorithm. For instance, it was used to investi-
gate reaction pathways in the methanol to gasoline conversion process over
zeolite catalysts [12, 13]. Reaction paths and energy barriers for the carbon-
oxygen (C-O) bond cleavage and first carbon-carbon (C-C) bond formation
were explored using all-electron periodic supercell calculations with DMol3 .
Calculations along these lines can help understand mechanism of a reaction
and ultimately improve existing catalysts.
The cleavage of the C-O bond of ethanol was studied by calculating the
reaction path for the methylation of a zeolitic oxygen at the aluminosilicate
Brønsted acid site. A reaction barrier of ∼ 54 kcal/mol was calculated for
such reaction. This is a concerted SN2-type reaction involving breaking of the
C-O bond in methanol and bond formation between the C atom and surface
oxygen, with the formation of a water molecule resulting from the transfer
of a zeolite proton to the hydroxyl group. A lower activation barrier of 44
kcal/mol was found in the presence of a second methanol molecule. The pres-
ence of a second methanol molecule in the zeolite cage causes the spontaneous
deprotonation of the zeolite, leading to the formation of the methoxonium
ion [12, 13].
Surface methoxyl species can undergo a reaction with incoming methanol

molecule to form the ethanol molecule. The barrier for such reaction is 25
kcal/mol, if the water molecule is present. If the water is not involved in this
reaction, the activation energy increases to about 50 kcal/mol. Therefore,
water seems to be an important catalytic agent for the initial formation of
the first C-C bond [12, 13].
A recent hypothesis concerning the formation of adsorbed ylide (CH2)

species as an initial precursor for the first C-C bond formation was also ex-
plored. A very high calculated barrier of 78.5 kcal/mol ruled out the insertion
of a ylide group from methoxyl into the Al-O bond. Such reaction would in-
volve a significant rearrangement of the zeolite framework and it is clearly not
competitive as compared with reactions involving adsorbed methoxyl species
[12, 13].

2.4 Transition State Confirmation Algorithm
The transition state found by the LST/QST/CG method in Section 3 may

not be the transition state connecting reactants and products or may be the
transition state for a different reaction or reaction step. It is also possible that
more than one transition state or other minima exist on a complex reaction
path. These additional stationary points may be missed by the LST/QST
procedure. Therefore, mapping of a reaction path is an important part of
studying reactivity. The simplest way to calculate a reaction path is to start
at a saddle point and take successive steps in the direction of the negative
gradient. This steepest descent approach leads to a minimum energy path
(MEP). If the coordinate system is mass-weighted, this is called an intrinsic
reaction coordinate (IRC). The MEP (or IRC) path may be quite complicated
and may have several minima. The highest saddle point is of most interest,
as the overall reaction rate depends on the height of this reaction barrier.
Following the reaction path allows to discover intermediate structures and
to connect the reaction barrier to the correct reactant and product. DMol 3
uses the nudged elastic band (NEB) method [16] to validate a transition state
by connecting it to the proper reactant and product [18]. NEB introduces a
fictitious spring force that connects neighboring points on the path to ensure
continuity of the path and projection of the force so that the system con-
verges to the MEP.
Fig. 2.6. Schematic representation of the minimum energy pathway. R: reactant,

P: product, and TS: transition state. M1 and M2 are local minima.
Consider the reaction pathway in Figure 2.6. The LST/QST algorithm

takes the endpoints R and P as input and locates one of the local maxima
on the reaction coordinate. Assume, in this example it locates the highest
energy barrier, TS. The NEB algorithm will take these three points as input

and return a trajectory that contains at least one point in the vicinity of
the new minima, in this case M1 or M2. Points are shown in Figure 2.6 as
blue (minima) and red (maxima). In general, the path will also contain a few
points that do not correspond to stationary points (green). A conventional
NEB is intended to start from the end points and yield the TS and the entire
reaction pathway, i.e., many green points. In contrast, the algorithm imple-
mented in DMol3 is intended to start at the TS, and to locate the alternative
stationary points in the direction of reactants and products.
The optimization consists of 2 phases termed macroiteration and microi-

teration. It is assumed that the points are held loosely in place by springs.
Macroiterations consist of an optimization over all the images; microiterations
consist of the constraint optimization of the molecular or crystal in a direc-
tion orthogonal to the reaction pathway. The calculation is complete when
both the macroiterations and microiterations have converged. The result is
an approximation to the MEP based on the number of structures used in
the defined trajectory between the reactants and products. The newly found
stationary points on the path will not be fully optimized because they have
been subjected to constraints. However, they will be sufficiently close to their
minimum structures, so that fu rther full optimization only takes very few
steps and the energetics will not change substantially.
To illustrate the usefulness of this new NEB algorithm we consider the

Diels-Alder reaction [19] of cis-butadiene and ethylene reacting to cyclohex-
ene. The reactants and the product were optimized using DMol3 with a DNP
(double numeric + polarization) basis set and the BLYP functional [20]. A
subsequent LST/QST calculation found the transition state structure de-
picted in Figure 2.7.
Fig. 2.7. Calculated transition state structure for the Diels-Alder reaction of bu-
tadiene and ethylene to cyclohexene.

A normal mode analysis performed on this activated complex resulted

in one imaginary frequency of 508i conforming that the LST/QST algorithm
successfully found the transition state of this reaction. However, the predicted
activation energy was only 17.6 kcal/mol being about 30 % smaller than the
experimental value of 25 kcal/mol [17]. Such a substantial deviation cannot be
attributed to approximations made by the DFT method. Therefore, a NEB
calculation was performed in order to check if any minima on the reaction
path were missed. The result is displayed in Fig. 2.8.
Fig. 2.8. Resulting graph of the NEB calculation, and geometry of this new min-
imum structure. A new minimum is found at path coordinate 0.23. Note that the
cis-butadiene exhibits a bent structure. The green curve is the energy vs QST path
1. The yellow, blue, pink, and magenta curves denote the energy vs MEP paths.
The graph in Fig. 2.8 shows that the NEB procedure has found a new
minimum on the reaction path, at path coordinate 0.23, where the path co-
ordinate is defined between the reactants (0.0) and the products (1.0). The
new minimum structure shown in Fig. 2.8 exhibits a reactant complex in
which the cis-butadiene is bent rather than planar as in the reactant input
structure. This new minimum corrects the activation energy to 23.9 kcal/mol

in good agreement with experiment. A subsequent full optimization of all

stationary points on the reaction path refined the activation energy to 23.3
kcal/mol, still in satisfying agreement with experiment.
2.5 Chemical Bonding and Elastic Properties of

Corundum-Type Oxides: The Rhodium Oxide Case
Commercial needs for improving durability and performance of supported
metal catalysts drive attention in investigating structural changes occurring
at high oxidizing environments. In Rh-based three-way emission control cat-
alysts, formation of reducible rhodium sesquioxide (Rh2 O3 ) has been associ-
ated to the catalyst deactivation [18]. Moreover, understanding the nature of
the bonding in transition metal oxides posses a challenge to ab initio meth-
ods, mainly due to the treatment of the exchange and correlation in partially
occupied d bands when the material is reduced.
The band structure and elastic properties of the rhodium sesquioxide

(Rh2 O3 ) in the corundum structure were calculated using density-functional
theory (DFT) and the generalized gradient approximation (GGA) of Perdew,
Burke, and Ernzerhof (PBE) [19] for the exchange-correlation functional. The
calculations were performed using the plane-wave pseudopotential-based code
CASTEP [20]. The Bulk properties of Rh2 O3 such as lattice constants (a0 ,
c0 ), bulk modulus (B0 ) and its pressure derivative (B0 ) were obtained by cal-
culating the total energy (E(V)) as a function of the unit-cell volume (V) and
by fitting the calculated energy values to the Birch-Murnaghan equation of
state (EOS) [21]. The minimum of the Murnaghan curve yields the optimum
value of the energy, volume (V0 ), lattice parameters (a0 and c0 ), and bulk
modulus (B0 ) at zero pressure, as reported in Table 2.1.
Table 2.1. Calculated unit-cell edges (a0 , c0 ), volume per oxide formula-unit (V0 ),
internal fractional coordinates (z(Rh), x(O)). Percentage errors are reported in the
second line for each quantity.
a0 (c0 /a0 ) V0 z(Rh) x(O)
Rh2 O3 5.135 Å 2.710 52.96 Å3 0.3488 0.5498
[22] +0.2 % +0.3 % +1.1 % +1.4 % +0.9 %

The bulk modulus as extracted from the Birch-Murnaghan equation of

state, B0 (∂ 2 E/∂ 2 V)V0 , is calculated to be 2.30 Mbars, and its pressure deriva-
tive 5.3. B0 is related to the form of the function E(V) around V0 . The
predicted B0 value is similar to that measured for other transition metal
sesquioxides in the corundum structure as Fe2 O3 (2.25 Mbars) [23] and Cr2 O3
(2.38 Mbars) [23]. Hence, the rhodium sesquioxide is predicted to be as stiff
as the corresponding iron and chromium structures under isotropic compres-
sion. Sofar, this value is not known experimentally for Rh2 O3 .
CASTEP was also used to calculate the full elastic constant tensor. The
calculated elastic constants and bulk modulus obtained as a linear combi-
nation of elastic constants are summarized in Table 2.2. CASTEP uses the
finite strain technique to predict elastic constants, which applies a given ho-
mogeneous deformation (strain) and calculates the resulting stress. The bulk
modulus obtained by this technique of 2.29 Mbar agrees well with that ob-
tained by fitting the energy data to the Birch-Murnaghan EOS of 2.11 Mbar.
Table 2.2. Calculated elastic constants for corundum Rh2 O3 .
C11 C33 C44 C12 C13 C15
Rh2 O3 406 287 181 150 177 -18

GPa
The band structure and density of states reveal that corundum rhodium
is a Mott-Hubbard insulator (see Fig. 2.9) with a metal-to-metal charge tran-
sition dominated by the dt2g (valence band edge) and deg metal (conduction)
at an energy of 1.4 eV. This is consistent with the observed semiconducting
behaviour of Rh2 O3 [24]. The band structure presents three groups of bands:
the low lying bands at about 20 eV correspond to oxygen 2p states, the group
of bands between 8 and 4 eV originates mainly from oxygen 2p states, the
bands from about 3 eV up to the Fermi level are of mainly Rh-4dt2g charac-
ter, and the conduction band consists of mainly oxygen 2p states hybridized
with Rh-4deg levels.

Fig. 2.9. Band structure and density of states (DOS). Red, green and blue curves
represent the p-, d- and s-DOS, respectively. The energies are relative to the Fermi
level in eV.
2.6 Summary
Industrial research in new materials development or improvement of existing
ones for catalysis applications, crystal growth (e.g. chemical vapour deposi-
tion), and microelectronics are commonly based upon phenomenological the-
oretical models, and high throughput experimentation. Recent advances in
quantum technology provide the materials researcher with tools to thought-
fully develop intelligent descriptors based upon the molecular understanding
of the underlying processes. A novel optimization technique for solid-state
calculations has been presented, which shows a superior performance for com-
plex problems, e.g. surface reactions. A powerful algorithm to optimize the
reaction path of a process by DFT within the transition state theory has

been described. This method allows not only the evaluation of the energy
difference between the particle at a minimum and the saddle point of the po-
tential energy surface (energy barrier), but also, it validates a transition state
by connecting it to the actual reactant and product. The constant pressure
optimization algorithm in CASTEP allows calculating the equation of state,
and the finite strain technique to predict the elastic properties for any kind
of existing or novel materials.
Acknowledgment
The authors are grateful to Professor Mathew Neurock for his advice in the
project on dissociation of maleic anhydride on Pd(111).
References
1. B.C. Gates. In: Catalytic Chemistry, (John Wiley & Sons 1992).
2. J. Andzelm, R.D. King-Smith, G. Fitzgerald: Chem. Phys. Lett. 335, 321
(2001).
3. K. Kudin, G.E. Scuseria, H.B. Schlegel: J.Chem.Phys. 114, 2919 (2001).
4. J. Baker, A. Kessi, B. Delley: J. Chem. Phys. 105, 192 (1996) .
5. J. Andzelm, G.Fitzgerald, N.Govind, D. King-Smith: J. Chem. Phys., submit-
ted (2003).
6. B. Delley,: J. Chem. Phys. 92, 508 (1990); J. Phys. Chem. 100, 6107 (1996) .
7. E. B. Wilson, J.C. Decius and P.C. Cross. In: Molecular Vibrations, (McGraw-
Hill, New York 1955).
8. V. Pallassana, M. Neurock, G. Coulston: J. Phys. Chem. 103, 8973 (1999).
9. Xu, Goodman, Langmuir 12, 1807 (1996).
10. T.A. Halgren, W.N. Lipscomb: Chem. Phys. Lett. 49, 225 (1977).
11. S. Bell, J.S. Crighton: J. Chem. Phys. 80, 2464 (1984); S. Fischer, M. Karplus:
Chem. Phys. Lett. 194, 252 (1992); J.E. Sinclair, R. Fletcher: J. Phys. C 7,
864 (1974).
12. N. Govind, J.W. Andzelm, K. Reindel, G. Fitzgerald: Int. J. Mol. Sci. 3, 423
(2002).
13. J. Andzelm, N. Govind, G. Fitzgerald, A. Maiti: Int. J. Quant. Chem. 91, 467
(2003).
14. A. Maiti, N. Govind, P. Kung, D. King-Smith, J. Miller: C. Zhang, G. Whitwell,
J. Chem. Phys. 117, 8080 (2002).
15. M. Petersen: Comp. Mater. Sci (submitted for publication).
16. G. Henkelman, G. Johannesson, H. Jonsson: Progress in Theoretical Chemistry
and Physics, ed. by S.D. Schwartz, (Kluwer Academic Publishers 2000).
17. K.N. Houk, R.J. Loncharich, J.F. Blake, W.L Jorgensen: J. Am. Chem. Soc.
111, 9172 (1989).
18. G.L. Kellog: J. Catal. 92, 167 (1985); G.L. Kellog: Surf. Sci. 171, 359 (1986).
19. J.P. Perdew, K. Burke, M. Ernzerhof: Phys. Rev. Lett. 77, 3865 (1996).
20. M.D. Segall, P.L.D. L indan, M.J. Probert, C.J. Pickard, P.J. Hasnip, S.J.
Clark, M.C. Payne: J. Phys.: Cond. Matt. 14, 2717 (2002).

21. F. Birch: J. Geophys. Res. 57, 227 (1952); F.D. Murnaghan: Proc. Natl. Acad.
Sci. U.S.A. 30, 224 (1944).
22. C.E. Boman: Acta Chem. Scannd. 24, 116 (1970); K.R. Poeppelmeier, J.M.
Newsam, J.M. Brown: J. Solid State Chem. 60, 68 (1985).
23. L.W. Finger and R.M. Hazen, J. Appl. Phys. 51, 5362 (1980).
24. G. Bayer and H.G. Wiedemann, Thermochimica Acta 15, 213 (1976).
View publication stats

Lec Notes MEGHalle

Uploaded by

Copyright:

Available Formats

Lec Notes MEGHalle

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Lec Notes MEGHalle

Uploaded by

Copyright:

Available Formats

See discussions, stats, and author profiles for this publication at: https://www.researchgate.

10. Computational Materials Science with Materials Studio$^{\rlap{\LARGE\

Chapter · May 2004

Maria-Elena Grillo Niranjan Govind

George Fitzgerald Klaus Stark

SEE PROFILE SEE PROFILE

Advanced FMR and EPR studies on gold nanoclusters View project

The user has requested enhancement of the downloaded file.

Computational modelling of novel materials is an increasingly powerful tool

For personal use only

molecules and clusters are discussed. The application of molecular-dynamics

Halle/Saale, Wolfram Hergert

For personal use only

Part I Basic Description of Electrons and Photons in Crystals

Part II Simulation from Nanoscopic Systems to Macroscopic

2 Computational Materials Science with Materials Studio:

For personal use only

Part III Modern Methods of Scientific Computing

For personal use only

Prof. Dr. H. Eschrig D. Hermann

Dr. A. Ernst Dr. S.F. Mingaleev

For personal use only

University of Surrey Dr. M. Celino

M. Däne Prof. Dr. V. Stepanyuk

For personal use only

W. A. Adeagbo D-47048 Duisburg

M. Kreth Dr. K. B. Stark

For personal use only

XII List of Contributors

Dr. G. Goldbeck-Wood [email protected]

For personal use only

The longstanding goal of computer simulations of materials is the accurate

For personal use only

correlations in many-body systems. The Green’s function provides spectral

For personal use only

Because of similarities between the Schrödinger equation and the Maxwell’s

For personal use only

Simulation from Nanoscopic Systems to

For personal use only

2 Computational Materials Science with

M. E. Grillo1 , J. W. Andzelm2 , N. Govind2 , G. Fitzgerald2 , and

Summary. In this article, developments in the primary ab initio quantum me-

Understanding and controlling the chemical and physical processes behind

For personal use only

based on the nudged elastic band method to validate a transition state by

2.2 Geometry Optimization

Structural optimization of materials is a pre-requisite in most studies using

For personal use only

2 Materials Studio: Applications in Catalysis 11

According to Baker et al. [4], the best combination of internal coordi-

The connection between internal and Cartesian displacements is defined

For personal use only

Bakers scheme of construction of the best delocalized internal coordinates

An optimization procedure carried out in the space of non-redundant inter-

The natural bonding topology of some crystals yields a set of coordinates

The performance of the new delocalized coordinate optimizer versus the

For personal use only

2 Materials Studio: Applications in Catalysis 13

The need for efficient geometry optimization in studying catalytic reac-