Atomic-scale computational materials science

Acta mater. 48 (2000) 71±92 www.elsevier.com/locate/actamat

ATOMIC-SCALE COMPUTATIONAL MATERIALS SCIENCE p

J. HAFNER 1, 1

2

Institut fuÈr Materialphysik, UniversitaÈt Wien, Sensengasse 8, A-1090 Wien, Austria and 2Centre for Computational Materials Science, UniversitaÈt Wien, Sensengasse 8, A-1090 Wien, Austria (Received 1 June 1999; accepted 15 July 1999)

AbstractÐAn overview on methodology and leading-edge applications of atomistic computational materials science based on quantum mechanical concepts is presented. # 2000 Acta Metallurgica Inc. Published by Elsevier Science Ltd. All rights reserved. Keywords: Computational materials science; Theory and modelling; Electronic and magnetic structure

1. INTRODUCTION: THE ROLE OF COMPUTATIONS IN MATERIALS SCIENCE

With the rapid development of advanced technologies in all areas, from light-weight constructions in the automotive and aircraft industries over very large-scale integration in microelectronics (and the upcoming nanoelectronics) to the chemical and biomedical industries, advanced materials science is confronted with the need for an increasingly precise control of ®nal product properties and the creation and design of new materials to meet speci®c demands. In parallel with this need it has come to recognition that it can be met eciently only by achieving a thorough understanding of how atomic, molecular, and mesoscopic features in¯uence macroscopic behaviour. Quantum mechanical calculations provide a means to understand and to predict the interactions between atoms and molecules and to model chemical reactions on a microscopic scale. Methods based on statistical mechanics and on continuum mechanics can use this information as input to investigate the mesoscopic and macroscopic behaviour. The last decades have witnessed important theoretical and algorithmic advances in the theory of condensed matter and in quantum chemistry. Together with the revolution that has occurred in computer technology, these advances enable us to

p

The Millennium Special Issue Ð A Selection of Major Topics in Materials Science and Engineering: Current status and future directions, edited by S. Suresh.

tackle problems of practical importance. As a result, atomistic computational technologies (computational quantum mechanics, molecular simulations and molecular mechanics) begin to bridge the gap between fundamental materials science and materials engineering. Computational materials science, however, raises some speci®c questions. In principle, the laws governing the behaviour of electrons and ions in a solid, a ¯uid or in a molecule are well known: it is sucient to solve the SchroÈdinger equation (or its relativistic analogue, the Dirac equation, in the case where the material contains heavy elements or if we want to calculate properties that depend on the coupling of the spin of the electron to its spatial degrees of freedom). The challenge is to ®nd this solution not for the simple case of the hydrogen atom (where the analytical answer can be found in any elementary textbook), but for the complex problem of the almost in®nitely many electrons and ions forming a solid and to achieve at the same time the exceedingly high accuracy necessary to energetically di€erentiate ferromagnetic a-Fe from paramagnetic (or antiferromagnetic) g-Fe or to estimate the heat of reaction and the reaction barrier in a catalytic process. The complexity of the problem, combined with the need to achieve high accuracy has two important consequences: (a) the calculation always has the character of a computer experiment, and (b) the stability, accuracy and eciency of the numerical algorithms are of utmost importance. The role of computer experiments in materials science is two-

1359-6454/00/$20.00 # 2000 Acta Metallurgica Inc. Published by Elsevier Science Ltd. All rights reserved. PII: S 1 3 5 9 - 6 4 5 4 ( 9 9 ) 0 0 2 8 8 - 8

72

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

fold: ®rst, computation allows theory to be applied to a complex problem: the laboratory experiment tells us what happens, theory, combined with computation, allows us to ®nd out why it happens and hopefully also to predict what will happen under similar circumstances. However, as the actual computer experiment is often quite demanding, the predictive capacity is often rather restricted. It can be greatly enhanced if the computer experiment is used to build a model on which simple, yet reliable predictions can be based. The second important role of

Hion ˆ

N X Iˆ1

hÿ 2 @ 2 ÿ *2 2M @ R I

! ‡

* * * 1 X Vion±ion … RI ÿ RJ † ‡ Etot …fRI g† 2 I;Jˆ1

a computer experiment is to supplement the experiments performed in the laboratory. Although impressive progress has been achieved in manipulating single atoms (e.g. at a solid surface using the atomic force microscope) and even single electrons in quantum devices, atomically resolved experimental information can often not be drawn from laboratory experiments, but is readily accessible by computer. The important point is a proper design of the computer experiment: hypothesize a mechanism at the atomic level, decide which steps can be performed in the laboratory (and use these results to validate those of the computation) and which are better analysed via computer. The explanatory and predictive capability of computer-assisted theory, as well as the reliability of the computer experiment, can explain the ever-increasing impact of computation on modern materials science. The intention of this paper is to analyse the stateof-the art of computational materials science at the atomistic level. After a brief review of some of the

Helec ˆ

" n X iˆ1

ÿ

terials is the Hamiltonian of the many-ion-manyelectron system. The ®rst step towards a reduction of its complexity is the adiabatic or Born± Oppenheimer approximation. Since the dynamics of the electrons occurs at a time-scale that is much faster than that of the much heavier ions, it is legitimate to assume that irrespective of the instantaneous con®guration of the ions, the electrons are at every moment in their ground state. The ionic system is described by an e€ective Hamiltonian

I6ˆJ

with an e€ective potential composed of the direct ion±ion interaction and the electronic total energy * as a function of the ionic coordinates fRI g: In most cases (in fact for all but the lighter atoms) the quantum nature of the ionic dynamics may be neglected and the SchroÈdinger equation based on the Hamiltonian (1) may be replaced by the classical Newtonian equation of motion. The adiabatic approximation rarely breaks down, although for a description of certain phenomena (electronic transport, supraconductivity, . . . ) the coupling of the ionic and electronic dynamics (electron±phonon coupling ) is fundamental. 2.2. Reducing the many-electron Hamiltonian In the adiabatic approximation, the Hamiltonian of the many-electron system includes the kinetic energy of n electrons, their interaction with the * external potential of N ions located at the sites RI and the interaction of the electrons with each other

# N n X N X * * 1X e2 h2 @2 * * : ‡ V …r , † ‡ el±ion i RI *2 2 i j6ˆi ri ÿ rj 2m @ r Iˆ1 ÿ

…1†

…2†

i

fundamental aspects of a quantum theory of materials, the capabilities and limitations of the most commonly used program packagesÐthe tools in our computer laboratoryÐare discussed. Exemplary applications illustrating the present status of computational materials and possible future developments are outlined, concentrating on the use of parameterfree ab initio methods.

2. FUNDAMENTALS OF A QUANTUM THEORY OF MATERIALS

2.1. Adiabatic approximationÐdecoupling electrons and ions The starting point of a quantum theory of ma-

The solution of the time-independent SchroÈdinger equation may be attempted at di€erent levels of theory. Monte Carlo methods (MC), based either on variational quantum MC (VMC), di€usion quantum MC (DMC), or path-integral MC (PIMC) may be used to tackle the quantum-many-body problem [1, 2]. While the application of quantum MC techniques to the homogeneous electron gas [3] provides a very valuable input for the construction of approximate one-electron approaches, their application to real solids is still at a pioneering stage [4]. For applications in the ®eld of materials two approximations reducing the many-electron problem to an e€ective one-electron form are widely used: Hartree±Fock (HF) theory and density-functional theory (DFT).

HAFNER: COMPUTATIONAL MATERIALS SCIENCE *

2.3. Hartree±Fock theory The basic approximation of HF theory is to approximate the many-electron wavefunction by an antisymmetrized product of one-electron wavefunctions and to determine these by a variational condition applied to the expected value of the Hamiltonian. In the resulting one-electron equations, 2 6 ÿh @ 4ÿ ‡ Vel±ion ‡ 2m @*r2 2

2

ition r is equal to the exchange-correlation energy of an electron in an homogeneous electron gas of a * density n0 equal to the local electron density n…r†: A HF description of the electron gas leads to a simple form of the exchange-only energy functional * * Ex ‰n…r†ŠAn4=3 …r†: A much more accurate exchangecorrelation energy for the homogeneous electron gas as a function of density may be derived from

3 * 0 2 *0 … …  *0 X cj …r † * 0 7 * X cj …r †ci …r † * 0 * * 2 5 * * 0 dr cj …r† ˆ ei ci …r† * * 0 dr ci …r† ÿ e dsi sj r ÿ r r ÿ r j j

the pair-wise electron±electron repulsion is replaced by the interaction of the electron i with the average electrostatic ®eld created by the charge distribution of all other electrons and an additional exchange term keeping electrons of like spin …si ˆ sj † away from each other to account for the Pauli principle. While exchange is treated exactly in HF theory, correlations arising from short-distance Coulomb interactions are neglected. While in calculations for molecules and clusters, beyond-HF correlation can be incorporated by con®guration-interaction (HFCI) schemes, this is dicult in HF calculations for periodic systems [5]. For this reason HF and HF-CI calculations are most widely used in molecular quantum chemistry. 2.4. Density-functional theory (DFT) Density-functional theory is based on the famous theorem by Hohenberg and Kohn [6] who demonstrated that the total energy of a many-electron system in an external potential is a unique functional of the electron density and that this functional has its minimum at the ground-state density. Expressing * the electron density n…r† as a sum over one-electron * densities jci …r†j2 and using the one-electron wave functions as the variational parameters leads to the Kohn±Sham one-electron equations [7] "

73

h2 @2 ‡ Vel±ion ‡ ÿ 2m @*r2 ÿ

*

n…r† ˆ

2

e

…

i

quantum MC simulations [3] and used to construct exchange-correlation functionals within the framework of a local density approximation (LDA), the most widely used such functionals are due to von Barth and Hedin [8], Perdew and Zunger [9], and Vosko et al. [10]. The relaxation of the constraint of an equal occupation of spin-up and spin-down states and the minimization of the total energy with respect to spin occupancy leads to the local spin-density approximation (LSDA) which allows magnetic properties of atoms, molecules and solids to be treated [8]. However, within a non-relativistic or scalar-relativistic description the spin direction is not coupled to the spatial degrees of freedom. Spin-orbit coupling is introduced only through a relativistic description on the basis of the Dirac equation. Alternatively, it may be added to a scalar-relativistic approach as a ®nal perturbation [11]. The LDA is generally very successful in predicting structures and macroscopic properties, but some shortcomings are notorious as well. These concern in particular: (i) energies of excited states, and in particular band gaps in semiconductors and insulators are systematically underestimated. This is after all not so surprising since DFT is based on a theorem referring to the ground state only. (ii) There is a general tendency to overbinding, i.e. cohesive

# *0 n…r † * 0 * * * * * 0 dr ‡ Vxc …r† ci …r† ˆ ei ci …r† r ÿ r

n X * 2 c …r † , i

where the exchange-correlation potential has been expressed as the functional derivative of the exchange-correlation energy. While equation (4) is formally exact within DFT, the problem is that the * functional Exc ‰n…r†Š is unknown. 2.4.1. Local density approximation. The ®rst practical approximation is to assume that the exchangecorrelation energy of a single electron at the pos-

…3†

…4†

*

*

Vxc …r† ˆ

dExc ‰n…r†Š * dn…r†

…5†

energies are signi®cantly overestimated and lattice parameters are underestimated by up to 3%. (iii) The wrong ground state is predicted for some magnetic systems (the most notable example is Fe which is predicted to be hexagonal close packed and non-magnetic instead of body-centred cubic and ferromagnetic) and for strongly correlated systems (e.g. the Mott insulators NiO and La2CuO4 are predicted to be metallic in the LDA). (iv) Van der Waals

74

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

interactions are not appropriately described in the LDA, although there are some recent suggestions for overcoming this problem [12±14]. 2.4.2. Generalized gradient corrections. Attributing the limitations of the local description to the neglect of the dependence of the exchange-correlation functional on the local variations of the electron density, generalized-gradient approximations (GGA) have been introduced. In the GGA, there is an explicit dependence of the exchange-correlation functional on the gradient of the electron density. As a straightforward expansion in terms of the gradient violates the sum-rules for the exchange hole, generalized gradient expansions corrected for the sumrules have been proposed by a number of authors. However, at the moment there is as yet no consensus on the best GGA. For solid-state applications, the GGAs proposed by Perdew and co-workers [15±19] have been widely used and have proved to be quite successful in correcting some of the de®ciencies of the LDA: the overbinding is largely corrected (the GGAs lead to larger lattice constants and lower cohesive energies) [20, 21] and the correct magnetic ground state is predicted for ferromagnetic Fe [22] and antiferromagnetic Cr and Mn [23± 25] (this also includes a greatly improved prediction of magnetovolume and magnetostructural e€ects). However, there are also cases where the GGA overcorrects the de®ciencies of the LDA and leads to an underbinding. The most notable examples are noblegas dimers and N2 molecular crystals which according to the GGA would not be bound [26, 27]. A further limitation of some of the GGAs is to incorrectly estimate the strength of hydrogen bonds, leading, e.g. to a wrong estimate of the freezing point of water. Hydrogen bonds seem to be better described on the basis of hybrid functionals including an exact description of exchange. 2.4.3. Further improvements of DFT. For an improved prediction of orbital energies, various strategies have been developed: self-interaction corrections [28], optimized e€ective potentials [29, 30], or adding Hubbard corrections …LDA ‡ U [31]). E€orts to improve the description of excited states are based on time-dependent DFT [32, 33]. It is expected that some of these extensions of the LDA will be built into the standard tools for materialsoriented applications in the near future. 2.5. Basis functions *

In practice, the wavefunctions ci …r† must be represented as a linear combination of a ®nite number of basis functions. The choice of the basis determines the achievable accuracy and computational eciency. In addition, di€erent basis sets may be more or less convenient for calculating given materials properties.

Methods applicable to complex systems employ mostly one of three types of basis sets, namely: (i) linear combinations of atomic orbitals (LCAOs); (ii) linearized augmented plane waves (LAPWs); or (iii) plane waves (PWs) in combination with pseudopotentials for describing the electron±ion interaction. The choice of a basis is problematic because near to an atom wavefunctions and potential are atomic-like (i.e. almost spherical-symmetric and strongly varying with radial distance), while in the interstitial regions potential and wavefunctions are quite smooth. In LCAO calculations, atomic orbitals are expressed as products of angular momentum eigenfunctions and radial orbitals, in calculations with periodic boundary conditions appropriate phase factors are added to form correct Bloch functions. The radial functions are either represented numerically or, to facilitate the calculation of multi-centre integrals, in terms of linear combinations of Slatertype or Gaussian-type orbitals (STOs and GTOs). In a minimal basis only LCAOs corresponding to the quantum numbers of the occupied states are used. If orbitals are added which are unoccupied in the electronic ground state, these are called polarization functions. LCAOs are ecient for systems with localized electrons, they are the basis of choice for performing Hartree±Fock calculations (see, e.g. Pisani et al. [34]) and have recently enjoyed renewed popularity in attempts to construct algorithms for which the computational e€ort scales linearly with the number of inequivalent atomic sites (linear scaling or O(N ) methods, see Refs [35, 36]). In the LAPW approach space is divided into spheres centred around the atomic sites and the interstitial region. Within the spheres the basis functions are spherical waves, in the interstitial region plane waves are used. The coecients of the atomic-like expansion within the spheres and of its energy derivative are chosen such that they match the plane-wave expansions at the sphere boundaries in value and slope. The intrinsic energy dependence of the eigenfunctions is linearized by expanding the orbitals to ®rst order with respect to the deviation from an appropriately chosen reference energy [37], thus transforming the non-linear into a linear eigenvalue problem. Modern versions of the LAPW method do not use any shape approximation to the e€ective one-electron potential and are therefore referred to as full-potential-LAPW or short FLAPW [38±40]. Similar basis sets are used also in the linear-mun-tin-orbital (LMTO) [41, 42] and Kohn±Korringa±Rostocker [43] methods. The LMTO method using the so-called atomic-sphere approximation (ASA) replacing the Wigner±Seitz sphere of equal volume is attractive because of its high computational eciency. A canonical transformation to the most localized LMTO basis allows contact to be made with LCAO-based techniques in the form of the tight-binding LMTO technique [44].

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

The TB-LMTO method is particularly well suited for performing calculations in real space without using Bloch's theorem and is hence particularly well suited for aperiodic systems. Plane-wave basis sets o€er several advantages: (i) convergence with respect to the completeness of the basis set is easily checked by extending the cut-o€ energy (i.e. the highest kinetic energy in the PW basis); (ii) fast-Fourier-transforms facilitate the solution of the Poisson equation; and (iii) forces on the atoms and stresses on the unit cell may be calculated directly via the Hellmann±Feynman theorem, without applying Pulay corrections for the site-dependence of the basis set [45]. The drawback is that in order to achieve convergence with a manageable size of the basis set, the strong electron±ion interaction must be replaced by a suciently weak pseudopotential. The basic idea of the pseudopotential approach is to project the SchroÈdinger equation for the valence electrons onto the subspace orthogonal to the core orbitals [46, 47], eliminating the nodal structure of the valence orbitals close to the core without modifying them in the region where chemical bonding occurs. Much e€ort has been spent in constructing pseudopotentials designed to satisfy the con¯icting requirements of accuracy, transferability and computational eciency. For the s,p-bonded main-group elements a variety of pseudopotentials that work well have been proposed [48±51], but the dicult cases are the ®rst-row elements and the transition and rare-earth elements with nodeless 2p, 3d and 4f orbitals where the usual concepts fail. For these elements the so-called ultrasoft pseudopotentials (US-PP) supplementing a PW representation of the smooth part of wavefunctions, charge densities and potentials with localized augmentation functions are by far the most ecient solution [52, 53]. Recently, the close correlation between US-PPs and the FLAPW method has been established [54], a one-to-one correspondence between US-PPs and the projector-augmented-wave (PAW) method proposed by BloÈchl [55] has been demonstrated by Kresse and Joubert [56]. To date modern FLAPW, PAW and US-PP techniques provide the most accurate LDF calculations, with the computational eciency increasing in that order. Di€erences between the US-PP approach and the all-electron PAW or FLAPW approaches exist only in those cases where the linearization of the valence±core exchange interaction becomes problematic (some ferro- and antiferromagnetic systems, early 3d elements, . . . ) [56, 57]. 2.6. Solving the eigenvalue problem: diagonalization or O(N ) methods The central task in any electronic structure calculation is the determination of the eigenvalues and the eigenvectors of the Hamiltonian, represented in matrix form in an appropriately chosen basis. This

75

leads to a diagonalization problem. Direct diagonalization using standard numerical techniques [58] leads to a scaling of the computational e€ort with the third power of the number Mb of basis functions ‰O…M3b † scaling], hence the application of these techniques to larger systems is severely restricted. Replacing the direct diagonalization of the Hamiltonian by a minimization of the total energy with respect to the electronic wavefunctions, performed either in the spirit of a dynamical-simulated-annealing approach as proposed by Car and Parrinello [59] (see below for a more detailed discussion) or via preconditioned conjugate gradient minimizations [60, 61] represented a substantial algorithmic progress as the scaling with the number of basis functions is now Mb log Mb if a plane-wave basis is used (the scaling being determined by fastFourier-transform between direct and momentumspace representations). However, parts of the calculations, associated with the orthogonality requirement of the eigenstates, scale with the third power of the number of atoms Nat, so that very large systems are again inaccessible. Improvements (leading also to a better stability of the algorithm, especially for metals) can be achieved if total-energy minimization is replaced by a band-by-band minimization of the eigenvalues [62, 62a], combined with ecient charge±density mixing routines. However, this does not eliminate the orthogonalization constraint. On the other hand, the minimization of the norm of app the residual vector R‰ci Š ˆ k…H ÿ eapp i †ci k to each app app eigenstate (where ei and ci are the current approximations to eigenvalue and eigenvector of H ) is a constraint-free-variational condition (subject only to a subspace-diagonalization after running over all bands) [63, 63a, 63b], and allows the number of orthogonalization operations to be eciently minimized so that scaling is e€ectively O…N2at † even for the largest systems that can be treated on modern computers. Strategies for achieving O(N ) scaling are based on expressions of all quantities (energies, forces, . . . ) in * *0 terms of the density matrix n~ …r, r † and exploit the fast decay of the o€-diagonal matrix elements of * *0 n~ …r, r †: To obtain linear scaling one has to cut o€ the exponentially decaying quantities where they become suciently small, introducing the concept of a localization region. The existing approaches to achieve linear scaling have very recently been reviewed by Goedecker [64]. They fall essentially into two classes. The ®rst exploits the fact that the decay of the density matrix is exponential at ®nite temperature [65] and uses an expansion of the Fermi occupation function f(E ) in terms of Chebyshev polynomials [66] or rational functions to construct a computable approximation to the matrix function F ˆ f …H † (where H is the one-electron Hamiltonian of the system). The second class of methods determines the density matrix by minimizing a functional representing the grand potential either directly with

76

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

respect to n~ [67] or with respect to the coecients of a representation of the density matrix in terms of Wannier orbitals [35, 36, 68]. The crucial point for assessing the usefulness of O(N ) methods is the cross-over point, i.e. the system size where the O(N ) methods become faster than diagonalization methods. Compared with standard diagonalization methods, a cross-over point of about 100 atoms has been estimated [67]; compared with iterative diagonalization techniques it is considerably higher. Goedecker [64] also gives a critical evaluation of the state-of-the-art in the development and application of the various O(N ) techniques whose full potential will probably be accessible only after a further substantial increase in computer performance.

2.7. Beyond one-electron theories In principle, the exact many-body ground state and the excited states may be found using (multi-) con®guration-interaction (CI) techniques. However, CI techniques require computer times scaling as N! where N is the number of particles. In contrast, Monte Carlo simulations for the ground state of many-body systems scale as N 3 so that they can be applied to suciently large systems to allow extrapolation to the bulk limit. Various types of quantum Monte Carlo (QMC) techniques have been developed. Variational Monte Carlo (VMC) uses statistical sampling methods to evaluate expectation values of operators with many-body trial functions including correlation e€ects [1, 2, 69, 70]. By optimizing the parameters of the wavefunctions one can ®nd an upper bound to the energy and determine which types of electron±electron correlations are most important. Di€usion Monte Carlo [1±3, 71] starts from a variational wavefunction and iterates the SchroÈdinger equation in imaginary time to ®nd the best approximation to the exact ground-state energy and wavefunction. For bosonic systems with nodeless ground-state wavefunctions, QMC techniques can always ®nd the exact ground state. For Fermions due to the ``sign-problem'' arising from the antisymmetry of the wavefunctions, no exact method applicable to more than a few particles has been developed [1, 2], but the accuracy of techniques based on the ``®xed node'' approximation is well established. Core states with large eigenvalues and short characteristic time-scales are dicult to handle in QMC. The way to circumvent this diculty is to perform many-body calculations for the valence electrons only and to describe the core± valence interaction by pseudopotentials [69±71]. The most extensive applications of QMC methods have been made on the homogeneous electron gas [3] and metallic hydrogen at very high pressures [72, 73]. Applications to carbon- and silicon-based semiconductors show consistently improved predictions of the cohesive and structural

properties compared with LDA and HF calculations [69, 70]. These calculations also yield accurate results for the pair-correlation functions of electrons and hence bring valuable information on the form of the exchange-correlation hole. This information is very valuable for constructing improved DFT functionals, but QMC techniques will soon also become tools of computational materials science in their own right. 2.8. Moving the atoms Beyond the determination of the ground state and of the excited states of the electrons, the challenge is to determine the equilibrium atomic structure of the system, to compute its dynamic and thermodynamic properties and to model phase transformations and chemical reactions by solving the equations of motion based on the adiabatic Hamiltonian of the ionic subsystem, see equation (1). Except for systems involving very light atoms (mostly hydrogen), a classical treatment of the atomic dynamics is usually sucient. 2.8.1. Car±Parrinello, Hellmann±Feynman and tight-binding molecular dynamics. The variational property of the ground-state energy has a very im* portant consequence: for a given con®guration fRI g, * * the force FI acting on the atom at site RI is given by the expected value of the derivative of the * Hamiltonian H with respect to RI in the electronic ground state jF0 i   @H * FI …fRI g† ˆ ÿ F0 * F0 : @ RI

*

…6†

Equation (5) is the famous Hellmann±Feynman theorem which allows the full set of quantum-many* body forces FI to be calculated which can then be used to optimize the atomic geometry or to study the dynamics of the atoms by integrating the Newtonian equations of motion, * *  …t† ˆ * MI R FI …fRI …t†g†, I

…7†

the challenge being to perform a very accurate determination of the electronic ground state after each move of the ions. In a seminal paper published in 1985, Car and Parrinello [59] pointed out that the optimization with respect to the ionic and electronic degrees of freedom need not be done separately. If a dynamical-simulated-annealing strategy is used to determine the ground state described by the orbitals {ci}, one arrives at a coupled set of pseudoNewtonian equations of motion for the ionic coor* dinates fRI g and the electronic orbitals {ci}  …t† ˆ ÿ MI R I *

*

@ E…fRI g, fci g† * @ RI …t†

…8†

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

and me c i …t† ˆ ÿ

* @ E…fRI g, fci g† X Lij …hcj ci i ÿ dji † ‡ @ ci …t† i,j

…9† where the Lagrange multipliers Lij are used to ensure conservation of orthonormality of the {ci } during the dynamical evolution of the system. In this approach, the system is not exactly in its adiabatic ground state, but moves at a ``distance'' from the Born±Oppenheimer surface given by the sum of 2 the ®ctitious kinetic energies …1=2†me c_ i …t† of the electronic orbitals ci where me is the mass assigned to the electronic degrees of freedom. The problem is to prevent the system from drifting away from the adiabatic ground state. For insulating or semiconducting systems this is unproblematic: due to the existence of an energy gap the characteristic frequencies for electronic excitations are orders of magnitude higher than the ionic eigenfrequencies and electrons and ions are essentially in a metastable decoupled equilibrium. For metals on the other hand, electrons and ions readily exchange energy, and the system quickly drifts away from the Born±Oppenheimer surface. An ad hoc solution is to drain the excess kinetic energy through the coupling of the electrons to a Nose thermostat [73a]. Dynamical Car±Parrinello calculations have been performed for systems with up to a few hundred atoms and extending over a few picoseconds. The problem of non-adiabaticity does not exist if the dynamics is based on the exact Hellmann±Feynman forces [equation (6)], even for transition metals with a high electronic density of states at the Fermi level. Indeed with the improved eciency of ab initio LDF calculations (cf. Section 2.6), a straightforward Hellmann±Feynman molecular dynamics is now feasible. One has to remember that in a Car± Parrinello molecular dynamics the time-step for the integration of motion has to be adjusted to the fast dynamics of the electronic degrees. Even if the mass me of the electronic orbitals is chosen to be several hundred times the electron mass, the time-step in Car±Parrinello simulations is typically at least one order of magnitude shorter than for comparable Hellmann±Feynman simulations. Simulations extending to about 100 ps are now feasible and open an access to the calculation of transport properties (di€usion, viscosity, . . . ) even in systems with a relatively slow single-particle dynamics [74]. An attempt to extend the applicability of MD simulations based on quantum-many-body forces to larger systems exploits the higher computational eciency of semiempirical TB techniques or of TBbased O(N ) techniques (cf. Section 2.6). The TBMD techniques based on the diagonalization of the Hamiltonian (see, e.g. Refs [75, 76]) allow larger systems to be studied over larger time-scales,

77

although systems containing more than a few hundred atoms remain inaccessible due to the O(N 3)scaling problem. TB-based O(N ) methods have been used in a few pioneering studies: the orbital minimization method was used to perform crash tests of fullerene molecules impinging on a diamond surface [77], Ajayan et al. [78] studied the irradiation damage on carbon nanotubes. Applications of DFT-based O(N ) techniques using basis sets large enough so that convergence can be claimed have not been reported yet. However, feasibility studies of such calculations in cells containing up to 6000 silicon atoms have been performed [79]. TB techniques have been also used in combination with Monte Carlo techniques. Bichara et al. [80] used a semiempirical TB-Hamiltonian and classical MC techniques to study liquid and amorphous Se and other disordered systems close to a metal/ semiconductor transition. KrajcÏõ and Hafner [81] have developed a fuzzy TB-Monte Carlo technique with O(N ) scaling and a cross-over point compared with standard TB-MC calculations of the order of 100 atoms. Other attempts to bridge the gap between the microsystems accessible to ab initio calculations and mesoscopic systems are based on the mapping of the energetics (heats of reactions, energy barriers, . . . ) onto classical Hamiltonians that allow for the application of well-elaborated tools of statistical mechanics. Such techniques have been used to study the growth of thin ®lm on solid surfaces (layer-by-layer vs insular growth) or the temperature-programmed desorption of atoms and molecules from surfaces [82]. 2.8.2. Classical molecular dynamics with e€ective interatomic forcesÐbridging the gap in the length scales. The main limitation of the modern mixed quantum/classical molecular dynamics is the restricted system size. On the other hand, classical molecular dynamics simulations may now be performed for systems with a million and more atomsÐhere the limitations arise from the restricted accuracy of the empirical forces. An alternative is to attempt to derive the interatomic forces from a well-founded quantum mechanical concept. For the s,p-bonded or simple metals such a concept is provided by linear-response theory in combination with a pseudopotential perturbation approach [46, 47]. Electron-density-dependent e€ective pair potentials indeed permit a successful calculation of the structure and dynamical properties of crystalline metals and of the structure, dynamic and thermodynamic properties of liquid metals and alloys (see, e.g. Refs [83, 84]). Extensions to many-body forces have been proposed [85]. Linear-response theory based on free-electron response functions fails for covalent semiconductors and for transition metals. Here a tight-binding approach o€ers a sounder basis. Earlier approaches to e€ective interatomic forces

78

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

based on low-order moment expansions have been put on a sounder basis by the tight-binding-bond (TBB) theory developed by Pettifor [86, 87]. TBB potentials have been developed on various levels of sophistication, ranging from renormalized perturbation theory on a Bethe-reference lattice (resulting in e€ective pair forces [88] well suited for describing liquid and amorphous systems) to many-body cluster expansions [89±91] for the construction of many-body forces. The calculation of the forces is based on an extension of the recursion method [92] [which is itself a well-established O(N ) technique for the calculation of the density of statesÐi.e. of diagonal elements of the density matrix] to o€-diagonal elementsÐthis establishes a link to other O(N ) techniques. However, it is necessary to extract some generalized moments H n from the recursion coecients [86, 87, 91], and both the stability of this inversion procedure and the convergence of the cluster expansion are still a subject of debate (see, e.g. Ref. [64])). The embedded-atom [93], e€ective-medium [94, 95], Finnis±Sinclair [96] methods and the glue model [97] are all based on the decomposition of the total energy in the generic form Etot ˆ

X * 1X * F…ni † ‡ F…Ri ÿ Rj † 2 i I6ˆJ

…10†

where ni is an e€ective electron density at the site of the electron i, F is a nonlinear function describing the interaction of the atom with the surrounding electrons and F is a pair-potential (but not the same as that resulting from linear-response theory). The electron density ni in equation (10) is an appropriately averaged electron density coming from the surrounding atoms ni ˆ

X * * na …Ri ÿ Rj †

…11†

I6ˆJ

where na is not necessarily a free-atom density. The functions na, F and F may be determined empirically (as in the embedded-atom, Finnis±Sinclair, and glue models) or derived theoretically at various levels of sophistication, as in the di€erent variants of the e€ective-medium theory. In its most developed form [95], the e€ective medium theory is ab initio in the sense that all input parameters may be calculated within the LDA, without ®tting to experiment, but approximate in the sense that it uses simpli®cations in the total-energy expressions and numerical interpolations to make the approach computationally ecient. The particular advantage is that the ®rst term in equation (10) accounts to some extent for the many-body character of the interatomic forces and for their local ¯uctuations in inhomogeneous systems (such as clusters, surfaces, interfaces, crystals with vacancies . . . ). Classical force ®eld simulations have now been performed for systems comprising several millions

of atoms. The necessity to use such large systems arises not only from the long-range nature of elastic strain ®elds even when the interatomic forces themselves are short-ranged, but is also intimately connected with the propagation of collective excitations. During the collision of clusters or molecules with solid surfaces or during crack propagation in a crystal, phonons are emitted. It is important that these phonons are not re¯ected from the boundaries of the simulation box during the time-scale of the collision or crack-formation event [98]. This example also illustrates that not all parts of the problem need necessarily be treated at the same level of sophistication. Whereas in the immediate region of impact of the cluster on the surface or of the crack tip, the breaking and reformation of bonds should be described on a quantum mechanical level, phonons are usually quite well described by classical force ®elds, and the more distant regions of the strain ®eld might even be suf®ciently well described by continuum mechanics. The possibility to treat the entire problem even with the most advanced and fastest O(N ) methods seems to be extremely remoteÐand even if this could be done it would constitute an immense overkill in the sense of treating a very large number of essentially inactive atoms with highly accurate techniques. The strategy to be followed seems to be evident: the expensive DFT or TB methods are applied in the core region, with boundary conditions determined by the results of the force-®eld calculations in the outer region. Eventually a second embedding shell based on ®nite-element calculations might be used. In principle, the input data necessary for the force®eld and ®nite-element calculations could be derived from the full-¯edged quantum mechanical calculations. The combination of molecular methods of di€erent accuracy and speed is not an easy task, but this will certainly be one of the main lines of further development. 2.8.3. Quantum dynamics. In the preceding discussion it has been assumed that the positions and momenta of the atoms or ions may be treated as classical variables. This approximation breaks down if light particles are involved: a dramatic manifestation of the quantum nature of the atomic nuclei takes place in liquid helium where the l-transition to the supra¯uid 4He phase cannot be understood if the quantum character of the bosonic 4He nuclei is ignored. Most problems in condensed-matter physics where the quantum nature of the atomic dynamics is relevant arise from the presence of hydrogen. The most extreme case is solid hydrogen which undergoes several phase transformations which cannot be understood without accounting for the quantum dynamics of the nuclei [99, 100]. Other examples are phase transitions in ice, or the behaviour of hydrogen adsorbed in semiconductors or metals, or at solid surfaces. Here it has been

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

demonstrated that the dissociation dynamics of molecular hydrogen is strongly in¯uenced by the quantized nature of the rotational and vibrational excitations of the H2 molecule [101]. Small molecules may show such an extreme anharmonicity of their internal vibrations that the structure of the molecule is determined primarily by the quantum ¯uctuations of the nuclei and less by the electronic e€ects on the molecular bonding [102]. Most attempts to describe the quantum dynamics of atoms and nuclei are based on the path-integral formulation of the canonical partition function of a system with the standard Hamiltonian Hion ˆ

N *2 X * PI ‡ V…fRI g†, 2M I Iˆ1

…12†

*

where V…fRI g† is the potential-energy function for the nuclei [103]. The numerical techniques for evaluating the path integrals are described in many monographs (see, e.g. Refs [1, 2, 104, 105]). Diculties arise from the fact that the kinetic and potential energy terms of the Hamiltonian Hion cannot be diagonalized in a common basis. A possible solution is to re-formulate the operator function exp(ÿbHion) as a product of exponential functions of the non-commuting operators for the kinetic and potential energies using the Trotter theorem [106] 8 2 3 < N *2 X b P I 5 exp…ÿbHion † ˆ lim exp4 ÿ p 4 1: P Iˆ1 2MI   exp ÿ

9  =P

b * VfRI g : ; P

…13†

For practical calculations the minimal value of P leading to acceptable quantum-statistical averages must be used. Most applications are based on model-potential-energy functions, using, e.g. Lennard±Jones potentials to describe the interaction between helium atoms. However, this approximation breaks down if the interatomic forces are very complex, in particular in situations where chemical bonds are broken and reformed. This occurs for example at the dissociative adsorption of H2 on metallic surfaces or at the dissociation of H2O into H+ and OHÿ, a process which is responsible for the pH-value of water. To deal with such situations, an ab initio pathintegral method combining the advantages of a classical molecular dynamics with ab initio quantummany-body forces with a quantization of the nuclear degrees of freedom has been developed [107, 108]. Like the Car±Parrinello MD, the technique is based on the adiabatic decoupling of the electronic and ionic degrees of freedom. A discretization of the partition function on the basis of the Trotter theorem generates a doubly generalized

79

Lagrangian depending on N  P nuclear and n  P electronic degrees of freedom. The simultaneous propagation of the corresponding equations of motion generates a quantum mechanical trajectory for the nuclei with the forces derived from an LDA treatment of the electrons. This approach represents the most advanced quantum mechanical treatment of both electrons and ions. In systems where the quantum dynamics of the lighter ions is e€ectively decoupled from the much slower dynamics of the system, a less demanding approach may be adopted: ab initio calculations may be used to determine an e€ective low-dimensional potential-energy-surface (PES) for the light particles, e.g. for the dissociation of a diatomic molecule over a surface for example, a six-dimensional PES must be calculated. A quantum description of the dissociation dynamics proceeds by solving the time-dependent SchroÈdinger equation for the molecule moving in the 6D-PES, using wave-packet dynamics [108a] or a transmission/re¯ection matrix approach based on a discretization of the PES along the reaction path coordinate [109, 110]. However, it must be pointed out that neither the mapping out of a high-dimensional PES nor the integration of a multidimensional SchroÈdinger equation are trivial tasks. 2.8.4. Rare-event dynamicsÐbridging the gap in the time-scales. Microscopic structural phenomena or chemical reactions often proceed on time-scales that are very long compared to those for atomic oscillations. This is the case, for example, for the structural relaxation in a glassy material, for the di€usion of an atom in a solid, or for the reaction between molecules leading to the formation of a new molecular species. In each case, the system is initially in a deep minimum of the con®gurational energy, surrounded by energy barriers many times higher than its temperature. Only rather rare ¯uctuations will allow the system to overcome the barrier and to ®nd a new potential energy minimum. As the rate with which such events occur decreases exponentially with the barrier height, the time-scale of such processes may reach macroscopic values. This renders atomistic simulations, whose time-step is set by the period of the atomic oscillations, very dicult. The problem is particularly acute for molecular dynamics simulations where the entire duration of a computer experiment can hardly be extended over more than a few nanoseconds. Very recently it has been proposed to combine transition-state theory and MD [111]Ðfor a simple model system a signi®cant speed-up of the system has been demonstrated, but it remains to be seen whether this can also be realized for systems of technological relevance. In Monte Carlo simulations, the diculty in simulating rare events is not related to any inherent time-scale, but depends rather on the nature of the

80

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

attempted moves. In liquids, single-atom moves are rather ecient [112], but they are not well suited for describing processes depending on the collective motion of groups of atoms. In lattice models such as the Ising model one can turn from single-spin ¯ips to ¯ips of clusters of spins to describe collective events determining the behaviour of longer times. For continuum-based models an activation± relaxation mechanism based on a two-step mechanism has been proposed [113]: the ®rst step consists of an activation of the system pushing it from a local minimum to a nearby saddle-point, the second of a relaxation into a new minimum. The technique has been applied to study relaxation and di€usion mechanisms in disordered materials [114, 115] and should be generalizable to other types of problems. Another attempt to achieve, simultaneously, chemical accuracy (necessitating the most advanced quantum mechanical calculations) and the timescales necessary to describe non-equilibrium thermodynamics is based on the mapping of the energetics (heats of reactions, energy barriers, . . . ) onto model Hamiltonians used in classical statistical physics. For example, Stamp¯ et al. [82] have used a lattice-gas Hamiltonian with parameters derived from ab initio density-functional calculations to simulate the temperature-programmed desorption of oxygen from ruthenium surfaces.

2.9. Tools in the computer laboratory A prerequisite to the successful application of the quantum mechanical concepts sketched above to real materials science problems is the availability of highly ecient, accurate, and stable computer codes. Essentially two sets of codes are required: one allowing to solve the SchroÈdinger equation at the required level of sophistication, and the second to use the information on the structures, wavefunctions, charge densities, etc. to calculate the required materials properties or to model a process in a material. Of the existing quantum mechanical techniques, only methods based on density-functional or Hartree±Fock theory have been implemented in general purpose codes executable on a variety of computer architectures and capable of treating a wide range of systems di€ering in size, chemical bonding, and structure. On the Hartree±Fock level the CRYSTAL program package of Pisani et al. [5, 34] working in an LCAO basis is essentially the only widely accessible tool for performing calculations on periodic systems. In practice, the accessible system size is restricted to a few dozen atoms. GAUSSIAN o€ers the possibility to perform cluster calculations using both HF and DFT methods, using a variety of localized basis sets, DFT exchange±correlation functionals and approximations to correlations within HF. Within a large number of utilities for structure optimization, bond-

ing analysis, etc., GAUSSIAN is a very professional package. Within DFT, a wider choice of highly ecient tools is available. The full-potential linearized augmented-plane-wave technique has been developed to a high level of eciency by several groups. The existing codes, for example the WIEN code distributed by Schwarz et al. [40] and the FLEUR code being developed by a consortium of European research groups [41a], o€er a variety of interesting features, including Hellmann±Feynman forces and stresses for geometry optimizations. However, applications are restricted to systems with a maximum of a few dozen atoms, and the eciency does not seem sucient for extended dynamical simulations. Ab initio dynamics is at the moment the realm of plane-wave based techniques. Among the existing plane-wave codes CASTEP (developed by the TCM group at Cambridge University and using ultrasoft pseudopotentials in the latest version) [60] and VASP (developed by G. Kresse in Wien, allowing both ultrasoft pseudopotentials and all-electron projector-augmented wave techniques to be used) [56, 62, 63a, 63b] are probably at the leading edge. O(N ) techniques have been implemented in the SIESTA code distributed by the Spanish group around Ordejon, Artacho et al. [35, 36]. All these tools begin to develop a serious impact on industrial research, but academic groups are undoubtedly at the fore-front of the methodological development. 3. APPLICATIONS TO SPECIFIC MATERIALS AND PROCESSES

3.1. Structural materials Historically, modern materials science has evolved from the ®elds of metallurgy and metal physics. It is therefore perhaps not inappropriate to start our brief overview of selected current applications of quantum mechanically based computational techniques in this area. Some time ago, Pettifor published in the New Scientist an article entitled `` New alloys from the quantum engineer'' describing how structural maps constructed using chemical coordinates introduced on the basis of quantum mechanical reasoning could possibly be used to ®nd new technologically useful alloys. The examples chosen to illustrate this approach were largely drawn from the ®eld of light-weight alloys. Here we want to demonstrate that the potential of computational materials science has developed to a point were it begins to rival with such concepts in explanatory and predictive ability. 3.1.1. Light-weight alloys. Al±Li alloys have attracted considerable attention because they represent promising materials for aerospace applications due to the combination of light weight and high strength. It is generally believed that the

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

81

Fig. 1. Side view of a stepped diamond (111) surface. (a) Starting con®guration after static relaxation. (b), (c) Side and top view of the reconstructed terraces at T ÿ 500 K, the Pandey chains formed by the reconstructed surface layer are shown in white. (d)±(h) Heating-induced graphitization proceeding through the slab. After Kern and Hafner [123].

strengthening e€ect is mainly due to the precipitation of the metastable d-phase (ordered Al3Li in the L12 structure). Understanding (and control!) of the strengthening e€ect requires the knowledge of the metastable phase diagram and of the strengthening mechanism. Sluiter et al. [116] have used a combination of ab initio FLAPW calculations and

the cluster-variation method to determine the phase diagram. Guo et al. used the FLAPW approach to determine the elastic constants of metastable Al3Li [117] (which could hardly be measured in the laboratory). The strong increase of Young's modulus of Al3Li compared to pure Al found in the calculations supports the view that the strengthening

82

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

e€ect is primarily due to the modulus di€erence between the precipitates and the matrix and not related to the creation of crystalline defects. Intermetallic light-weight high-strength alloys based on g-TiAl (L10 structure) are promising candidates for high-temperature applications. However, the alloy is brittle below 7008C. A large enhancement of the tensile fracture strain can be achieved in two-phase alloys of g-TiAl and a2-Ti3Al (DO3 structure) with small additions of a few atomic percent Nb, V, Cr, or Mn or combinations of these elements. Self-consistent total-energy calculations on ternary TiAlX performed on supercells of the gand a2-phases with up to 32 atoms allows the determination of the lattice distortions due to the alloying with the ternary element and the site preference [118, 119]. This allows insight to be gained into the conditions at the g/a2 interfaces which are held to be responsible for the ductility improvements. 3.1.2. Refractory compounds. Refractory transition-metal carbides and nitrides exhibit a combination of properties which has attracted attention with respect to both technological applications and fundamental materials physics: high melting points and ultrahardness typical for covalent materials are combined with metallic electrical and thermal conductivities. High-speed steels may be considered as a two-phase composite consisting of carbides embedded in a martensitic matrix. Predicting the overall strength of the compound using ®nite-element methods [120] requires an accurate knowledge of the mechanical properties of the carbide inclusions which are extremely dicult to measure. Wolf et al. [121] have performed ®rst principles calculations of the elastic and thermal properties of a series of carbides as a function of volume, pressure and vacancy concentration. The results demonstrate that the predictive ability of DFT can be rated very high for this class of materials. 3.1.3. Ultrahard materials. The fascinating properties of diamond and cubic boron nitride such as extreme hardness, high melting point, high thermal conductivity, and wideband-gap have motivated widespread applications ranging from high-temperature microelectronics to the coating of high-duty tools. Problems arise from the fact that both materials exist not only in a hard cubic form, but also as soft polymorphs with a layered structure (graphite and hexagonal BN). Very recently the microscopic mechanism of surface-induced graphitization of diamond was studied by de Vita et al. [122] and Kern and Hafner [123] using ab initio molecular dynamics. It was demonstrated that as in the bulk, the diamond 4 graphite transition is strongly ®rst order (and hence has a pronounced tendency to overheating) for a perfectly ¯at surface, but that atomic defects such as steps can lead to a substantial lowering of the transition

temperature. Figure 1 shows the transformation history of a stepped C (111) surface, starting from the relaxed …1  1† surface and the spontaneous …2  1† reconstruction of the terraces and proceeding via a progressive breaking of the interlayer bonds. The production of cubic BN by chemical or physical vapour deposition is considerably more dif®cult than that of diamond ®lms. A qualitative explanation can be given by the Ostwald rule stating that a system does not reach its ground state directly, but passes ®rst through the less-stable stateÐhence the formation of h-BN is expected at temperatures where c-BN is stable. However, until recently, there was no consensus as to the c-BN/hBN coexistence line. Kern et al. [124] have used ab initio total-energy and phonon-dispersion calculations to construct the phase diagram of BN: in contrast to carbon where graphite and diamond are energetically almost degenerate at low temperature, c-BN is lower in energy than h-BN by about 0.06 eV/atom. The transformation to h-BN at higher temperature is driven by the low-lying shear and optical modes corresponding to the weak interlayer coupling, its pressure dependence is strongly in¯uenced by the strong mode-dependence of the GruÈneisen parameters gi …gi < 0 for TO modes leading to a buckling of the sheets, gi ' 5 for some of the low-frequency modes). Here the computation leads to much more complete information than experiment which was so far restricted to Raman spectroscopy on the high-frequency optical modes (mostly for preparative reasons). 3.2. Alloys Research on alloys is one of the oldest and widest ®elds of materials science. Here we shall content ourselves with a rapid look at new developments in a traditionally established ®eldÐmartensitic transformationsÐand a brief discussion of the contribution of computations to the investigation of two new classes of alloys: quasicrystals and metallic glasses. 3.2.1. Martensitic transformations. Martensitic transformations are di€usionless structural transformations of ®rst order. Typically, a close-packed low-temperature phase transforms with increasing temperature into a less dense body-centred phase, the transformation being driven by the rapid increase of the vibrational entropy caused by lowenergy phonons in the less-dense phase. The concentration-dependence of the transformation is correlated to electronic band-®lling e€ects, which are most pronounced in Hume±Rothery alloys. Martensitic transformations and phase stability in s,p-bonded metals and alloys were successfully described in the early stages of modern computational materials science, see, e.g. Ref. [47]. Structural transformations in magnetic transition-

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

metal alloys are still more complex due to the strong temperature dependence of magnetic longrange order coupled to the structural degrees of freedom, leading to Invar and Anti-Invar e€ects. First-principles calculations performed in the local spin-density approximation and including generalized gradient corrections have greatly contributed to clari®cation of this complex scenario (for recent reviews see, e.g. Refs [125±127]). It has been shown that f.c.c. g-Fe shows a unique distance dependence of the contributions of the t2g and eg states to the magnetic exchange coupling, leading to a competition between high-moment states stabilized by the occupation of antibonding t2g states (leading to lattice expansion) and energetically almost degenerate ferro- and antiferromagnetic low-moment states arising for the depletion of these antibonding states. This competition is at the origin of the magnetovolume instabilities in Fe-based transition metal alloys. Phonon-induced electronic transitions from antibonding majority t2g to nonbonding eg states can stabilize the Invar e€ect over considerable ranges of temperatures and composition. Anti-Invar behaviour has been shown to arise from the occupation of high-moment antibonding states caused by spin ¯uctuations of large amplitude at high temperature. 3.2.2. Quasicrystals. The discovery by Shechtman et al. [128] that the X-ray di€raction pattern of rapidly quenched Al±Mn alloys has icosahedral (and hence crystallographically forbidden) symmetry represented not only a revolution in crystallography, but raised also new questions about structure, vibrational lattice modes, and electronic properties. It was recognized very soon that the existence of non-crystallographic point-group symmetry leads to quasiperiodic instead of periodic translational order [129]. For one-dimensional models of quasicrystals it has been demonstrated using renormalization group techniques that the eigenstates of the Hamiltonian are neither extended (as in periodic crystals) nor localized (as in amorphous materials), but critical, with a power-law decay of the amplitudes. So far, the analytical studies could not be extended to realistic three-dimensional systems. Numerical studies of quasicrystalline (icosahedral or decagonal) alloys are based on periodic approximants, i.e. a hierarchy of complex crystalline lattices with the same local order as the quasicrystal and converging to the quasiperiodic limit as the lattice constant approaches in®nity. The computational studies are extremely demanding (systems with up to about 55,000 atoms have been studied using TB-based techniques), they have largely contributed to the development of a generalized Brillouin-zone concept for vibrational and electronic eigenstates and to an understanding of the role of electrons in stabilizing a quasiperiodic ground state (for a recent review, see the papers by Fujiwara and Hafner and KrajcÏõ in Ref. [130]). On

83

the basis of these results, the main lines of an understanding of the highly unusual physical properties of quasicrystals, ranging from high electrical resistivities approaching a metal±insulator transition to giant magnetic moments on isolated sites, begin to emerge. 3.2.3. Glassy alloys. Amorphous alloys produced by rapid quenching are of interest because of the change in their physical properties compared to their crystalline analogues caused by the absence of grain boundaries and extended lattice defects. Here the challenges are, on the one hand, that experiment only delivers a one-dimensional projection of the three-dimensional atomic arrangement and, on the other hand, that a calculation of the atomic structure must be based on realistic interatomic forces (which should account at least for the main features of the electronic spectrum). For amorphous transition-metal and transition-metal±metalloid alloys an at least approximately consistent description of their atomic and electronic structure has been given on the basis of tight-binding methods and simulated MD quenches [131±134]. These studies elucidate the important role of d±d and d±p hybridization, respectively, on the short- and medium-range order, stability, electronic and magnetic properties. 3.3. Magnetic materials Research on magnetic materials remains an extremely wide ®eld, with new developments in many areas ranging from crystalline and nanocrystalline permanent magnets over soft-magnetic metallic glasses to two-dimensional magnetism and technological impact in many areas, ranging from the electric power industry to information storage in computers. Here we shall concentrate on recent results on thin magnetic ®lms and magnetic multilayers. 3.3.1. Thin magnetic ®lms. The behaviour of thin magnetic ®lms di€ers from that of bulk materials because of surface and interface e€ects and because the epitaxial relationship of the magnetic ®lm and the nonmagnetic metallic or semiconducting substrate can lead to non-negligible variations in the crystalline structure. A very widely investigated system is thin Fe ®lms on Cu (100) substrates where the f.c.c. structure of g-Fe is stabilized up to a thickness of 11±12 monolayers. The interplay of an enhanced magnetic moment at the free surface and at the interface with the substrate with the unusually strong magneto-volume e€ects in f.c.c. Fe leads to a very complex structural and magnetic phase diagram which could be completely elucidated only with the help of ab initio LSDA ‡ GGA calculations [135, 136]: in the region of up to three monolayers surface e€ects stabilize a ferromagnetic high-spin state and the magnetovolume e€ect leads

84

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

Fig. 2. Interlayer-exchange coupling (IEC) in Fe5/Crn multilayers as a function of the thickness of the Cr spacer layers. Positive IEC indicates antiferromagnetic, negative IEC ferromagnetic coupling of the Fe layers (left). The right side shows the pro®le of the magnetic moments in a Fe5/Cr30 multi-layer for ferro- and antiferromagnetic coupling. After Hafner et al. [140].

to a tetragonal distortion of the iron ®lm. For four and more monolayers the surface bilayer is ferromagnetically coupled, in the deeper layers bilayer antiferromagnetic structures develop in ®lms with an even number of layers, while for odd numbers of monolayers a variety of energetically nearly degenerate complex magnetic structures coexist. Expansion of the interlayer distance between ferromagnetically coupled layers and contraction between antiferromagnetically coupled planes lead to a structure that is isotropically f.c.c. on average in this thickness range. The lattice distortions also in¯uence the magnetic anisotropyÐa striking example has been given in a recent study of Ni/ Cu(100) ®lms [137] demonstrating that the interplay of surface e€ects and substrate-induced tetragonal distortion leads to a highly unusual sequence of spin-reorientation transitions with increasing ®lm thickness. 3.3.2. Magnetic multilayers. The discovery of giant or colossal magnetoresistance in magnetic multilayers [138] and the associated potential of applications in magnetic information storage technology has certainly spurred the interest in this fascinating class of materials. Fe/Cr multilayers are a system where a particularly rich variety of physical phenomena have been observed (oscillatory magnetic interlayer coupling, giant magnetoresistance, noncollinear magnetic structures). Very

recently Niklasson et al. [139] have used an LMTOGreens-function technique to study the existence of spin-density-wave states in Fe/Crn/Fe sandwiches with up to n ˆ 51 Cr monolayers. Hafner et al. [140] have used a TB-LMTO technique in real space to determine the interlayer-exchange coupling in . . . Fe5/Crn/Fe5/Crn . . . multilayers. The results summarized in Fig. 2 demonstrate a short oscillation period of two monolayers, and a superposition with a longer period leading to the variation of the amplitude of the coupling, in accordance with the experiments of Unguris et al. [141]. 3.4. Minerals The application of ab initio DFT techniques is also expanding rapidly into areas such as crystallography (a recent review of theoretical crystallography is given by Winkler [142]) and mineralogy. Here we discuss very brie¯y applications to silica and aluminosilicates. 3.4.1. Polymorphism in silica. Silica is a fundamental building block of many rock-forming minerals, it is widely used in the ceramic and glassforming industries and has an enormous potential of applications in optical ®bres, microelectronics and catalysis. On the other hand, because it occurs in so many di€erent structures, it is one of the most dicult materials to study. Recent studies [143] of a

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

Fig. 3. Energy±volume curves calculated for a-cristobalite and three structural variants of b-cristobalite. Note the energetic degeneracy extending over a substantial range of volumes and pressures. This explains the coexistence of so many structural forms of this mineral. After Demuth et al. [143].

85

86

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

large number of tetrahedrally coordinated SiO2 polymorphs (a,b-quartz, a,b-cristobalite, tridymite, keatite, coesite) and of octahedrally coordinated stishovite have demonstrated that the LDA leads to an excellent structural prediction, but gradient corrections are necessary for a quantitatively correct description of the pressure-induced transition from quartz to stishovite. These studies also lead to new insights on the relation between di€erent structural variants proposed for b-cristobalite. Figure 3 shows the unit cells of a-cristobalite and of three structural variants proposed for b-cristobalite, as well as the energy vs volume (pressure) relations for all four phases. The near degeneracy of the b-phase variants over a substantial range of pressures explains their occurrence in minerals with di€erent thermophysical histories. Ab initio calculations have also been used to explore possible structural transformations at very high pressures [144]. On the basis of an exploration of possible octahedrally coordinated structures, a post-stishovite sequence of stable phases stishovite 4 CaCl2-type 4 a-PbO2-type 4 Pa3 has been proposed and their structural and cohesive properties have been calculated. Such studies are an important contribution to the geophysics of the earth's crust. 3.4.2. Zeolites. Zeolites are microporous aluminosilicates with important applications as molecular sieves and as catalysts. The basic framework of zeolites consists of a loosely connected network of SiO4 tetrahedra, with a certain percentage of the Si atoms substituted by Al. The charge-imbalance created by the Al/Si substitution has to be compensated by the introduction of counterions (mostly alkalis) or by attaching H atoms to an oxygen next to an Al site. It is the acidity of these OH groups, together with the con®nement of the reactant within the micropores that leads to the unique catalytic properties of zeolites. The problem is that the reactivity of the protonated sites depends on their local environment, which can be characterized only insuf®ciently by experimental techniques, due to the at least partially random Al/Si distribution. Despite the complexity of zeolites, DFT calculations have been used to optimize their crystal structures for cells as large as that of mordenite (144 atoms in the purely siliceous form), to characterize the geometrical, vibrational and chemical properties of the acid sites and to model chemical reactions [145±148]. 3.5. Liquids Traditionally research on ¯uids has concentrated on systems where the interatomic potentials can be modelled by simple pair interactions (plus eventually three-body terms), independent of the thermodynamic state of the system. Situations where the potentials are inherently state-dependent as, e.g. in materials close to a metal±insulator transition, or

mediated by d-electron bonds were considered as lying outside the realm of conventional statistical liquid-state physics. This situation has been completely reversed with the advent of ab initio MD methods. Below we summarize some recent applications opening new views into old problems. 3.5.1. Near-critical metallic melts. It is now well known that expanded liquid metals such as the alkalis or Hg undergo a metal±insulator (M/I) transition in the vicinity of the ¯uid/gas critical point. For liquid Hg it has been shown that the M/I transition occurs at the liquid side of the coexistence line, at about twice the critical density. Ab initio MD simulations on this material have demonstrated that it is possible to achieve full agreement with the available di€raction data in the entire range between the triple and critical points and that the M/I transition is a simple band-crossing transition between the 6s and 6p bands [149]. This discards much discussed many-electron scenarios for this transition. 3.5.2. Metal/semiconductor transitions in ¯uid materials. The semiconductor to metal (SC/M) transition occurring in the molten chalcogenides at temperatures not too far above the melting point, and the associated changes in the structural and thermophysical properties (thermal expansion, viscosity, . . . ) have been investigated over many decades. Only very recently ab initio MD simulations for liquid Se [56, 150, 150a] have led to a clear picture of how the progressive breaking and entanglement of the polymeric Se chains leads to the formation of mid-gap states associated with one- and three-fold coordinated defects driving the SC/M transitions. Composition-dependent M/SC transitions occur also in many alloys of alkali- or alkaline-earth alloys with polyvalent elements from groups III to VI. The very rich structural chemistry of these so-called Zint alloys that are at the borderline between metals, semiconductors and salts has been experimentally investigated for many decades [150b]. However, only recently ab initio electronic structure methods have been suciently developed to allow an investigation of these often extremely complex crystal structures characterized by the existence of polyanionic clusters, and to test the conjecture that these units might also continue to exist in the molten phase. Recent work on crystalline and molten K±Sb alloys illustrates the potential of computational techniques [151, 152]. 3.5.3. Extreme temperatures and pressures. The properties of materials under extreme thermodynamic conditions such as they exist, e.g. in the earth's core are certainly outside the range of laboratory experiments, but can certainly be studied by advanced computations. Very recently, de Wijs et al. [74] have performed ab initio MD calculations of the viscosity of liquid Fe and liquid Fe/Si alloys under such extreme temperatures and pressures.

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

87

Fig. 4. The dye molecule monomethine adsorbed on an NaCl (100) surface. (a) Initial con®guration. (b) Stable con®guration obtained by dynamical simulations. After Sugihara et al. [155].

This information will certainly have a strong impact on geophysical modelling. 3.6. Molecular materials The investigation of the structural and electronic properties of molecules has certainly a long tradition in quantum chemistry. However, even in this area the possibility to perform dynamical simulations (classical or quantum) based on ab initio total energies and forces has opened new issues. Here we illustrate the potential of the new techniques at two examples. 3.6.1. Floppy molecules. Floppy (or ¯uxional) molecules are characterized by the existence of many quasidegenerate minima on the potential-

energy hypersurface that are separated only by small barriers. This leads to large amplitudes of the intramolecular motions such that the molecule is no longer adequately described by a single reference structure. Examples are protonized acetylene C2 H‡ 3 or protonized water clusters like H5 O‡ 2 Ðin these cases the analysis is even more dicult because of the quantum character of the H motions. Marx and Parrinello [153] have presented ab initio path-integral QMC studies of C2 H‡ 3 demonstrating how the conventional structural models resulting from electronic structure calculations (with the added proton in a bridging position between the C atoms) are modi®ed by the quantum mechanical tunnelling of all protons between the available sites. Parrinello's group has also used the same technique to study

88

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

the in¯uence of quantum e€ects on pressureinduced phase transitions in ice [154]. 3.6.2. Dye molecules. Dye molecules are of a very broad interdisciplinary interest ranging from biology to physics. One of the most discussed fundamental problems is the adsorption of the dye molecules on surfaces and the subsequent change of

conformation leading to a change in its physical properties. Very recently Sugihara et al. [155] have used ab initio MD to study the adsorption of charged dye molecules (trimethine [C19H17N2O2]+ and monomethine [C21H23N2]+) on an NaCl (100) surface. Figure 4 shows the change of the initially ¯at structure of the molecule to a twisted con®gur-

Fig. 5. Electron localization function calculated along a reaction path for the catalytic oxidation of CO on Pt (111). Note the localization of electrons in the bond charge illustrating the degree of covalent bonding. After Eichler and Hafner (unpublished work).

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

ation obtained by a dynamical relaxation, in agreement with experimental observation. 3.7. Interfaces Interfaces between di€erent types of materials are being studied in many di€erent contexts, ranging from metal/oxide interfaces (which are important in such widely di€erent areas as composite structural materials and catalysis), over microelectronics to electrochemistry (solid/liquid interfaces). Here we can mention only very brie¯y one example of stateof-the-art simulations. 3.7.1. Epitaxy of transition-metal silicides. The epitaxial growth of silicides on silicon or metallic substrates is of substantial interest in microelectronics. CoSi2 and NiSi2 are known to grow with excellent structural quality on Si (111) substrates, but the growth of Fe silicides is much more complex, leading eventually to the formation of metallic phases di€erent from the stable semiconducting bulk phases. Very recently Moroni et al. [156] have demonstrated that ab initio DFT calculations may be used to search for epitaxially stabilized structures of Fe, Co, and Ni silicides on Si or metallic substrates. It was shown that pseudomorphic phases of FeSi2 in the C1 (CaF2-type) structure and of CoSi and NiSi in the B2 (CsCl-type) structure soften dramatically under compressive strain induced by epitaxy on a (100) substrate. This supersoft e€ect is characterized by zero strain energy, constant volume and constant bond energies. This is a very good example for the predictiveÐand not merely explanatoryÐabilities of modern DFT calculations. Detailed calculations of CoSi2/Si interfaces, including the determination of Schottky-barrier height have recently been performed by Stadler et al. [157]. 3.8. Surfaces and gas±surface interactions The use of ab initio DFT techniques for calculating the structural and electronic properties of surfaces and the interaction of atoms and molecules with surfaces have now reached a quite mature stage. A recent review of ab initio studies of the adsorption of diatomic molecules on metallic and semiconducting surfaces, including classical and quantum simulations of dissociative adsorption/associative desorption has been given recently by Gross [158, 159]. Adsorbate-induced reconstructions have been discussed by Liem et al. [160] for O/Cu (110) and by Dong et al. [161] for H/Pd (110). A more complex situation has been considered by Raybaud et al. [162] who studied the adsorption of thiophene (C4H4S) on the catalytically active edge surfaces of MoS2, allowing for a complete relaxation of the adsorbate/substrate complex. The resulting equilibrium geometry shows that the molecule is bound by a direct Mo±S bond and a so-

89

called Z 5-bond (in the terminology used in organometallic chemistry) between the aromatic ring of thiophene and a metal atom of the substrate. The important point is that the molecule is activated with respect to both C±S bond-breaking and reduction of the aromatic character. An important development is also the ab initio simulation of scanning-tunnelling-microscopy images [163, 164] via a calculation of the current between the tip and the adsorbate±substrate complex. The particularly interesting application to the STM imaging of acetylene adsorbed on Pd (111) [163] shows that without a theoretical analysis the straightforward interpretation of the STM images is very likely to be misleading: the STM reproduces the p-orbitals of the C2H2 molecule adsorbed in a slightly tilted bridge con®guration by dark and bright spots, depending on the distance of the opposite p-lobes from the substrate plane (and depending to a certain extent on the STM tip itself). We expect that the combined use of STM imaging and DFT calculations will vastly improve the insight to be gained from such studies. We refer also to the recent review of Briggs and Fisher [164] on the combined use of STM experiments and atomistic modelling of molecules on semiconductor surfaces.

3.9. Catalysis Catalysis is one of the most complex ®elds of chemistry, and also one in which computational studies require the highest computational e€ortÐ the aim being not only to determine a ground-state con®guration, but to explore a reaction scenario and to determine the arrangement leading to the lowest reaction barrier. Nevertheless, it is precisely in this area where some of the most exciting studies in modern computational chemistry have been achieved. Nielsen et al. [165, 166] have demonstrated that Au, normally immiscible with Ni, forms a surface alloy, and predicted that the surface alloy has improved activity as a catalyst in steam-reforming reactions (i.e. the production of hydrogen from methane or other alkanes). This prediction has received striking con®rmation by laboratory experiments so that this can be considered as the ®rst example of the design of a new catalyst on the basis of DFT calculations. The oxidation of CO on Pt (111) surfaces (important for post-combustion catalysis in car exhausts) was studied by Alavi et al. [167] and Eichler and Hafner [168], providing complete reaction scenarios for CO ‡ O and CO ‡ O2 reactions at di€erent surface coverages. Figure 5 shows the variation of the electron-localization functions along the reaction path, illustrating the breaking and formation of molecular bonds.

90

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

Raybaud et al. [169] have calculated energy pro®les for hydrodesulfurization reactions of thiophene on MoS2 catalysts under di€erent sul®ding conditions and also studied the e€ect of Co atoms acting as catalytic promotors. These studies were extended by Cristol et al. [170] to more complex species such as benzothiophene and methylbenzothiophene. Ab initio techniques are also increasingly important in the investigation of acid-based catalytic reactions inside zeolites [145±148]. 4. SUMMARY

In the ®rst part of this paper I have attempted to present an overview on the methodology on which modern atomistic computational materials science is based, and to characterize the tools that are now available in the computer laboratory. In the second part I have discussed exemplary applications in a variety of ®elds, ranging from structural materials over magnetism to catalysis. Of course the choice of subject in both parts is largely in¯uenced by personal experience, expertise and interest. Nevertheless I hope to have succeeded in convincing the reader that atomistic computational materials science is, at this turn of the millennium, a ¯ourishing subject all set to explore new frontiers. REFERENCES 1. Ceperley, D. M. and Kalos, M. H., in Monte Carlo Methods in Statistical Physics, ed. K. Binder. Springer, Berlin, 1986, p. 145. 2. Nithingale, M. P. and Umrigar, M. C., Quantum Monte Carlo Methods in Physics and Chemistry. Kluwer Academic, New York, 1998. 3. Ceperley, D. M. and Alder, B. J., Phys. Rev. Lett., 1980, 45, 566. 4. Foulkes, W. M., Nekovee, M., Gaudoin, R. L., Stedman, M. L., Needs, R. J., Hood, R. Q., Rajagopal, G., Towler, M. D., Kent, P. R. C., Lee, Y., Leung, W. K., Porter, A. R. and Breuer, S. J., in High Performance Computing, ed. R. J. Allan, M. F. Guest, A. D. Simpson, D. S. Hentys and A. D. Nicole. Plenum Press, New York, 1988. 5. Dovesi, R., Roetti, C., Freyria-Fava, C., Apra, E., Saunders, V. R. and Harrison, N. M., Phil. Trans. R. Soc. Lond., 1992, A341, 203. 6. Hohenberg, P. and Kohn, W., Phys. Rev., 1964, B136, 864. 7. Kohn, W. and Sham, L. J., Phys. Rev., 1965, 140, A1133. 8. von Barth, U. and Hedin, L., J. Phys. C: Solid St. Phys., 1972, 5, 1629. 9. Perdew, J. P. and Zunger, A., Phys. Rev., 1981, B23, 5075. 10. Vosko, S. H., Wilks, L. and Nusair, M., Can. J. Phys., 1980, 58, 1200. 11. White, R. M., Quantum Theory of Magnetism. Springer, Berlin, 1983. 12. Hult, E., Andersson, Y. and Lundqvist, B. I., Phys. Rev. Lett., 1996, 77, 2029. 13. Andersson, Y., Langreth, D. C. and Lundqvist, B. I., Phys. Rev. Lett., 1996, 76, 102. 14. Kohn, W., Meir, Y. and Makarov, D. E., Phys. Rev. Lett., 1998, 80, 4153.

15. Perdew, J. P., Phys. Rev., 1986, B33, 8822. 16. Becke, A. D., Phys. Rev., 1988, A38, 3098. 17. Perdew, J. P., Chevary, A., Vosko, S. H., Jackson, K. A., Pedersen, M. R., Singh, D. J. and Fiolhais, C., Phys. Rev., 1992, B46, 6671. 18. Perdew, J. P. and Wang, Y., Phys. Rev., 1992, B45, 13244. 19. Perdew, J. P., Burke, K. and Ernzerhof, M., Phys. Rev. Lett., 1996, 77, 3865. 20. Kresse, G., FurthmuÈller, J. and Hafner, J., Phys. Rev., 1994, B50, 13181. 21. Seifert, K., Hafner, J. and Kresse, G., J. Phys.: Condens. Matter, 1995, 7, 3683. 22. Leung, T. C., Chan, T. C. and Harmon, B. N., Phys. Rev., 1991, B44, 2923. 23. Asada, T. and Terakura, K., Phys. Rev., 1993, B47, 15992. 24. Eder, M., Moroni, E. G. and Hafner, J., Surf. Sci. Lett., 1999, 243, 244. 25. Eder, M., Moroni, E. G. and Hafner, J., Phys. Rev. B. in press. 26. Perez-Jorda, J. M. and Becke, A. D., Chem. Phys. Lett., 1995, 233, 134. 27. Kern, G. and Hafner, J., Unpublished. 28. Perdew, J. P. and Zunger, A., Phys. Rev., 1991, B23, 5048. 29. Bylander, D. M. and Kleinman, L., Phys. Rev. Lett., 1995, 74, 3660. 30. Grabo, T. and Gross, E. K. U., Chem. Phys. Lett., 1995, 240, 141. 31. Anisimov, V. I., Aryasetiawan, F. and Lichtenstein, A. I., J. Phys.: Condens. Matter, 1997, 9, 767. 32. Gross, E. K. U., Dobson, J. F. and Petersilka, M., in Density Functional Theory, ed. R. F. Nalejawski, Topics in Current Chemistry, Vol. 181. Springer, Berlin, 1996, p. 81. 33. Petersilka, M., Gossmann, U. J. and Gross, E. K. U., Phys. Rev. Lett., 1996, 76, 1212. 34. Pisani, C., Dovesi, R. and Roetti, C., Hartree±Fock Ab-initio Treatment of Crystalline Systems, in Lecture Notes in Chemistry, Vol. 48. Springer, Berlin, 1988. 35. Ordejon, P., Artacho, E. and Soler, J. M., Phys. Rev., 1996, B53, 10441. 36. Sanchez-Portal, D., Ordejon, P., Artacho, E. and Soler, J. M., Int. J. Quant. Chem., 1997, 65, 453. 37. Andersen, O. K., Jepsen, O. and GloÈtzel, D., in Highlights of Condensed Matter Theory, ed. F. Bassani, F. Fumi and M. Tosi. North-Holland, Amsterdam, 1985, p. 59. 38. Wimmer, E., Krakauer, H., Weinert, M. and Freeman, A. J., Phys. Rev., 1981, B24, 864. 39. Singh, D. J., Pseudopotentials and the LAPW Method. Kluwer Academic, Norwell, 1994. 40. Schwarz, K. and Blaha, P., in Quantum Mechanical Calculations of Properties of Crystalline Materials, ed. C. Pisani, Lecture Notes in Chemistry, Vol. 67. Springer, Berlin, 1996, p. 139. 41. Skriver, H. L., The LMTO Method. Springer, Berlin, 1984. 41a. Di Pomponio, A., Continenza, A., Podloucky, R. and Vackar, J., Phys. Rev., 1996, B53, 9505. 42. Andersen, O. K., Phys. Rev., 1975, B12, 3060. 43. Turek, I., SÏob, M., Kudrnovsky, J., Drchal, V. and Weinberger, P., Electronic Structure of Disordered Alloys, Surfaces and Interfaces. Kluwer Academic, Boston, 1997. 44. Andersen, O. K., Jepsen, O. and SÏob, M., in Electronic Band Structure and its Applications, ed. M. Yussouf. Springer, Berlin, 1987, p. 1. 45. Pulay, P., Molec. Phys., 1969, 17, 197. 46. Heine, V., Solid State Physics, Vol. 24. Academic Press, New York, 1970, pp. 1 and 247.

HAFNER: COMPUTATIONAL MATERIALS SCIENCE 47. Hafner, J., From Hamiltonians to PhaseDiagrams. Springer, Berlin, 1987. 48. Bachelet, G. B., Hamann, D. R. and SchluÈter, M., Phys. Rev., 1982, B26, 4199. 49. Rappe, A. M. and Joannopoulos, J. D., in Computer Simulation in Materials Sciences, ed. M. Meyer and V. Pontikis. Kluwer Academic, Boston, 1999, p. 409. 50. Lin, J. S., Qteish, A., Payne, M. C. and Heine, V., Phys. Rev., 1993, B47, 4174. 51. Vanderbilt, D., Phys. Rev., 1985, B32, 8412. 52. Vanderbilt, D., Phys. Rev., 1990, B41, 7892. 53. Kresse, G. and Hafner, J., J. Phys.: Condens. Matter, 1994, 6, 8245. 54. Holzwarth, N. A. W., Matthews, G. E., Dumring, R. B., Tackett, A. R. and Zeng, Y., Phys. Rev., 1997, B55, 2005. 55. BloÈchl, P. E., Phys. Rev., 1994, B50, 17953. 56. Kresse, G. and Joubert, D., Phys. Rev., 1999, B59, 1758. 57. Moroni, E. G., Kresse, G., FurthmuÈller, J. and Hafner, J., Phys. Rev., 1997, B56, 15629. 58. Press, W., Flannery, B. P., Teukorsky, S. A. and Vetterling, W. T., Numerical Recipes: The Art of Scienti®c Computing. Cambridge University Press, Cambridge, 1986. 59. Car, R. and Parrinello, M., Phys. Rev. Lett., 1985, 55, 2471. 60. Payne, M., Teter, D., Allan, D., Arias, T. and Joannopoulos, J., Rev. mod. Phys., 1992, 64, 1045. 61. Stich, J., Car, R., Parrinello, M. and Baroni, S., Phys. Rev., 1989, B39, 4997. 62. Kresse, G. and Hafner, J., Phys. Rev., 1994, B49, 14251. 62a. Kresse, G. and Hafner, J., Phys. Rev., 1993a, B47, 588. 63. Wood, D. M. and Zunger, A., J. Phys. A: Math. Gen., 1985, 18, 1343. 63a. Kresse, G. and FurthmuÈller, J., Phys. Rev., 1996, B54, 11169. 63b. Kresse, G. and FurthmuÈller, J., Comp. Mat. Sci., 1996, 6, 15. 64. Goedecker, S., Rev. mod. Phys., 1999, 71, 1085. 65. Ismail-Beigi, S. and Arias, T., Phys. Rev., 1998, B57, 11923. 66. Goedecker, S. and Teter, M., Phys. Rev., 1995, B51, 9455. 67. Li, X. P., Nunes, R. and Vanderbilt, D., Phys. Rev., 1993, B47, 10891. 68. Ordejon, P., Drabold, D., Grumbach, M. and Martin, R. M., Phys. Rev., 1993, B48, 14646. 69. Fahy, S., Wang, X. W. and Louie, S. G., Phys. Rev. Lett., 1988, 61, 1631. 70. Fahy, S., Wang, X. W. and Louie, S. G., Phys. Rev., 1990, B42, 3503. 71. Li, X. P., Ceperley, D. M. and Martin, R., Phys. Rev., 1991, B44, 10929. 72. Ceperley, D. M. and Alder, B. J., Phys. Rev., 1987, B36, 2092. 73. Wang, X. W., Fahy, S. and Louie, S. G., Phys. Rev. Lett., 1990, 65, 2414. 73a. BloÈchl, P. E. and Parrinello, M., Phys. Rev., 1992a, B45, 9413. 74. de Wijs, G. A., Kresse, G., VocÏadlo, L., Dobson, D., AlfeÁ, D., Gillan, M. J. and Price, G. D., Nature, 1998, 392, 805. 75. Wang, C. Z., Chan, C. T. and Ho, K. M., Phys. Rev., 1989, B39, 8586. 76. Winn, M. D., Rassinger, M. and Hafner, J., Phys. Rev., 1997, B55, 5364. 77. Galli, G. and Mauri, F., Phys. Rev. Lett., 1994, 73, 1994.

91

78. Ajayan, P., Ravikumar, V. and Charlier, C., Phys. Rev. Lett., 1998, 81, 1437. 79. Goringe, C., Hernandez, E., Gillan, M. and Bush, I., Comp. Phys. Commun., 1997, 102, 1. 80. Bichara, C., Pellegatti, A. and Gaspard, J. R., Phys. Rev., 1994, B49, 6581. 81. KrajcÏõ , M. and Hafner, J., Phys. Rev. Lett., 1995, 74, 5100. 82. Stamp¯, C., Kreuzer, H. J., Payne, S. H., PfnuÈr, H. and Sche‚er, M., Phys. Rev. Lett., 1999, 83, 2993. 83. Hafner, J., Mater. Sci. Engng A, 1994, 178, 1. 84. Hafner, J., in Electron Theory of Alloy Design, ed. A. Cotrell and D. G. Pettifor. The Institute of Materials, London, 1992, p. 47. 85. Moriarty, J. A., in Many-Atom Interactions in Solids, ed. R. M. Nieminen, M. J. Puska and M. J. Manninen. Springer, Berlin, 1990, p. 158. 86. Pettifor, D. G., Phys. Rev. Lett., 1989, 63, 2480. 87. Pettifor, D. G., Bonding and Structure of Molecules and Solids. Clarendon Press, Oxford, 1995. 88. Hausleitner, Ch. and Hafner, J., Phys. Rev., 1992, B45, 115. 89. Pettifor, D. G. and Aoki, M., Phil. Trans. R. Soc. Lond., 1991, A334, 439. 90. Aoki, M., Phys. Rev. Lett., 1993, 71, 3842. 90a. Hors®eld, A., Bratkovsky, A., Pettifor, D. G. and Aoki, M., Phys. Rev., 1996a, B53, 1656. 91. Kreuch, G. and Hafner, J., Phys. Rev., 1997, B55, 8304. 92. Heine, V., Haydock, R. and Kelly, M. J., Solid State Phys., Vol. 35. Academic Press, New York, 1981. 93. Daw, M. S. and Baskes, M. I., Phys. Rev. Lett., 1983, 50, 1285. 94. Jacobsen, K. W., Nùrskov, J. K. and Puska, M. J., Phys. Rev., 1987, B35, 7423. 95. Chetty, N., Stokbro, K., Jacobsen, D. W. and Nùrskov, J. K., Phys. Rev., 1992, B46, 3798. 96. Finnis, M. W. and Sinclair, J. C., Phil. Mag. A, 1984, 50, 45. 97. Ercolessi, E., Tosatti, E. and Parrinello, M., Phys. Rev. Lett., 1986, 57, 719. 98. Zhou, S., Beazley, D., Lomdahl, P. and Holian, B., Phys. Rev. Lett., 1997, 78, 479. 99. Mao, H. K. and Hemley, R. J., Rev. mod. Phys., 1994, 66, 671. 100. Mazin, I. I., Hemley, R. J., Goncharov, A. F., Han¯and, M. and Mao, H. K., Phys. Rev. Lett., 1997, 78, 1066. 101. Gross, A. and Sche‚er, M., Phys. Rev., 1998, B57, 2493. 102. Benderskii, V. A., Makarov, D. E. and Wight, C. A., Adv. chem. Phys., 1994, 88, 1. 103. Feynman, R. P. and Hibbs, A. R., Quantum Mechanics and Path Integrals. McGraw-Hill, New York, 1965. 104. Ceperley, C. M., Rev. mod. Phys., 1995, 67, 279. 105. Chandler, D., in Liquids, Freezing, and Glass Transition, ed. J. P. Hansen, D. Levesque and J. Zinn-Justin. Elsevier, Amsterdam, 1991. 106. Schulman, L. S., Techniques and Applications of Path Integrations. Wiley, New York, 1981. 107. Marx, D. and Parrinello, M., Z. Phys., 1994, B95, 143. 108. Marx, D. and Parrinello, M., J. chem. Phys., 1996, 104, 4077. 108a. Kroes, G. J., Baerends, E. J. and Mowrey, R. C., Phys. Rev. Lett., 1997a, 78, 3583. 109. Brenig, W., Brunner, T., Gross, A. and Russ, R., Z. Phys., 1993, B93, 91. 110. Gross, A., Wilke, S. and Sche‚er, M., Phys. Rev. Lett., 1995, 75, 2718. 111. Voter, A. F., Phys. Rev. Lett., 1997, 78, 3908.

92

HAFNER: COMPUTATIONAL MATERIALS SCIENCE

112. Hansen, J. P., in Monte Carlo Methods in Statistical Physics and Chemistry, ed. K. Binder. Springer, Berlin, 1986. 113. Barkema, G. T. and Mousseau, N., Phys. Rev. Lett., 1996, 77, 4358. 114. Barkema, G. T. and Mousseau, N., Phys. Rev., 1998, B57, 2419. 115. Barkema, G. T. and Mousseau, N., Phys. Rev. Lett., 1998, 81, 1865. 116. Sluiter, M., de Fontaine, D., Guo, X. Q., Podloucky, R. and Freeman, A. J., Phys. Rev., 1990, B42, 10460. 117. Guo, X. Q., Podloucky, R. and Freeman, A. J., J. Mater. Res., 1991, 6, 324. 118. Wolf, W., Podloucky, R. and Rogl, P., Mater. Res. Soc. Symp., 1995, 364, 1005. 119. Ding, J. J., Schweiger, H., Wolf, W., Rogl, P., Vogtenhuber, D. and Podloucky, R., Proc. TMS Conf., San Diego, in press. 120. Antretter, T. and Fischer, F. D., Comput. Mater. Sci., 1996, 7, 247. 121. Wolf, W., Podloucky, R., Antretter, T. and Fischer, F. D., Phil. Mag. B, 1999, B79, 839. 122. de Vita, A., Galli, G., Canning, A. and Car, R., Nature, 1996, 379, 523. 123. Kern, G. and Hafner, J., Phys. Rev., 1998, B58, 13167. 124. Kern, G., Kresse, G. and Hafner, J., Phys. Rev., 1999, B59, 8551. 125. Entel, P., Ho€mann, E., Clossen, M., Kadau, K., SchroÈter, M., Meyer, R., Herper, H. C. and Yang, M. S., in The Invar E€ect, ed. J. Wittenauer. The Minerals, Metals and Materials Society, London, 1997, p. 87. 126. Entel, P., Herper, H. C., SchroÈter, M., Ho€mann, E., Kadau, K. and Meyer, R., in Displacive Phase Transformations and Their Applications in Materials Engineering, ed. K. Inoue, K. Mukherjee, K. Otsuka and H. Chen. The Minerals, Metals and Materials Society, London, 1998, p. 319. 127. Herper, H. C., Ho€mann, D. and Entel, P., Phys. Rev. B, 1999, B60, 3839. 128. Shechtman, D., Blech, I., Gratias, D. and Cahn, J. W., Phys. Rev. Lett., 1984, 53, 1951. 129. Levine, D. and Steinhardt, P., Phys. Rev. Lett., 1984, 53, 2477. 130. Stadnik, Z. M. (ed.), Physical Properties of Quasicrystals. Springer, Berlin, 1998. 131. Hausleitner, Ch. and Hafner, J., Phys. Rev., 1992, B45, 115. 132. Hausleitner, Ch. and Hafner, J., Phys. Rev., 1992, B45, 128. 133. Hausleitner, Ch. and Hafner, J., Phys. Rev., 1993, B47, 5689. 134. Hafner, J., Tegze, M. and Becker, Ch., Phys. Rev., 1994, B49, 285. 135. Asada, T. and BluÈgel, S., Phys. Rev. Lett., 1997, 79, 507. 136. Moroni, E. G., Kresse, G. and Hafner, J., J. Phys.: Condens. Matter, 1999, 11, L35. 137. Uiberacker, C., Zabloudil, J., Weinberger, P., Szunyogh, L. and Sommers, C., Phys. Rev. Lett., 1999, 82, 1289. 138. Baibich, M. N., et al., Phys. Rev. Lett., 1988, 61, 2472. 139. Niklasson, A. M. N., Johansson, B. and NordstroÈm, L., and , Phys. Rev. Lett., 1999, 82, 4544. 140. Hafner, R., SpisÏ aÂk, D., Lorenz, R. and Hafner, J., to be published. 141. Unguris, J., Celotta, R. J. and Pierce, D. T., Phys. Rev. Lett., 1992, 69, 1125.

142. Winkler, B., Z. Kristallogr., in press. 143. Demuth, T., Jeanvoine, Y., Hafner, J. and AÂngyaÂn, J., J. Phys.: Condens. Matter, 1999, 11, 3833. 144. Teter, D. M., Hemley, R. J., Kresse, G. and Hafner, J., Phys. Rev. Lett., 1998, 80, 2145. 145. Nusterer, E., BloÈchl, P. E. and Schwarz, K., Chem. Phys. Lett., 1996, 253, 448. 146. Jeanvoine, Y., AÂngyaÂn, J., Kresse, G. and Hafner, J., J. Phys. Chem. B, 1998, 102, 5573. 147. Demuth, T., Benco, L., Hafner, J. and Toulhoat, J. J. Phys. Chem. B, submitted. 148. Benco, L., Demuth, T., Hafner, J. and Hutschka, F., J. chem. Phys., submitted. 149. Kresse, G. and Hafner, J., Phys. Rev., 1997, B55, 7539. 150. Kirchho€, F., Kresse, G. and Gillan, M. J., Phys. Rev., 1998, B57, 10482. 150a. Kresse, G., Kirchho€, F. and Gillan, M. J., Phys. Rev., 1999a, B59, 3501. 150b. Winter, R., in Thermodynamics of Alloy Formation, ed. Y. Chang and F. Sommer. The Minerals, Metals and Materials Society, London, 1997, p. 143. 151. Seifert-Lorenz, K. and Hafner, J., Phys. Rev., 1999, B59, 829. 152. Seifert-Lorenz, K. and Hafner, J., Phys. Rev., 1999, B59, 843. 153. Marx, D. and Parrinello, M., Science, 1996, 271, 179. 154. Benoit, M., Bernasconi, M., Focher, P. and Parrinello, M., Phys. Rev. Lett., 1996, 76, 2934. 155. Sugihara, M., Meyer, H., Entel, P., Sakamoto, Y., Hafner, J. and Buss, V., Proceedings SHDS Conference, Duisburg, 1999. World Scienti®c, Singapore, in press. 156. Moroni, E. G., Podloucky, R. and Hafner, J., Phys. Rev. Lett., 1998, 81, 1969. 157. Stadler, R., Wolf, W., Vogtenhuber, D. and Podloucky, R., Phys. Rev., 1999, B60. 158. Gross, A., Surf. Sci. Rep., 1998, 32, 291. 159. Eichler, A., Gross, A., Hafner, J. and Sche‚er, M., Phys. Rev., 1999, B59, 13297. 160. Liem, S. Y., Kresse, G. and Clarke, J. H. R., Surf. Sci., 1998, 415, 194. 161. Dong, W., Letendu, V., Sautet, Ph., Kresse, G. and Hafner, J., Surf. Sci., 1997, 877, 56. 162. Raybaud, P., Hafner, J., Kresse, G. and Toulhoat, H., Phys. Rev. Lett., 1998, 80, 1481. 163. Dunphy, J. C., Rose, M., Behler, S., Ogletree, D. F., Salmeron, M. and Sautet, Ph., Phys. Rev., 1988, B58, R12705. 164. Briggs, G. A. D. and Fisher, A. J., Surf. Sci. Rep., 1999, 33, 1. 165. Nielsen, L. P., Besenbacher, F., Stensgaard, I., Laegsgaard, E., Engdahl, C., Stoltze, P., Jacobsen, K. W. and Nùrskov, J. K., Phys. Rev. Lett., 1993, 71, 754. 166. Nielsen, L. P., Besenbacher, F., Stensgaard, I., Laegsgaard, E., Engdahl, C., Stoltze, P., Jacobsen, K. W. and Nùrskov, J. K., Science, 1996, 279, 1913. 167. Alavi, A., Hu, P., Deutsch, Th., Silvestrelli, P. L. and Hutter, J., Phys. Rev. Lett., 1998, 80, 3650. 168. Eichler, A. and Hafner, J., Phys. Rev., 1999, B59, 5960. 169. Raybaud, P., Hafner, J., Kresse, G. and Toulhoat, H., Proc. Int. Symp. on Hydrotreatment and Hydrocracking of Oil Fragments, in press. 170. Cristol, S., Paul, J. F., Payen, E., Bougeard, D., Hutschka, F. and Hafner, J., Proc. Int. Symp. on Hydrotreatment and Hydrocracking of Oil Fragments, in press.