Spin Density Functional Approach to the Chemistry of Transition Metal

Chapter 16

Spin Density Functional Approach to the Chemistry of Transition Metal Clusters Gaussian-Type Orbital Implementation

Downloaded by COLUMBIA UNIV on March 8, 2013 | http://pubs.acs.org Publication Date: June 8, 1989 | doi: 10.1021/bk-1989-0394.ch016

1

1

2

J. Andzelm , E. Wimmer , and Dennis R. Salahub 1

Cray Research, Inc., 1333 Northland Drive, Mendota Heights, M N 55120 Département de Chimie, Université de Montréal, Montréal, Québec H3C 3J7, Canada 2

In this contribution we review the accuracy and computational efficiency of the local spin density functional (LSDF) approach using linear combinations of Gaussian-type orbitals (LCGTO) for the calculation of the electronic structures, ground-state geometries, and vibrational properties of transition metal compounds and clusters. Specifically this is demonstrated for (1) bis(π-allyl) nickel where this approach gives an excellent qualitative and quantitative interpretation of the observed photoemission spectrum; (2) chemisorption of C atoms on a Ni(100) surface, where the present com putational approach determines the adsorption site of C as the four fold-hollow position above the surface with a calculated C - N i bond length of 1.79 Å(exp.: 1.75 ± 0.05 Å) and a vibrational frequency of 407 cm , (exp.: 410 cm ); (3) vibrational frequencies of C O on Pd: the calculations reveal that inclusion of surface/subsurface Pd-Pd motions couple significantly to the CO-Pd vibration leading to a reduction of the vibrational frequency by about 20% compared with a rigid substrate model. Inclusion of an external electrical field shows a stiffening of the C-O vibration with increased positive potential of the electrode; and (4) the electronic structure of Zn clusters where it is found that properties such as the first (s-electron) ionization potential converge rather slowly towards the value of the extended system requiring at least 20 transition metal atoms for an accurate description of the surface electronic structure. It is demonstrated that the computation of threeindex two-electron integrals can be achieved with a highly efficient vector/parallel algorithm based on recursive integral formulas recently published by Obara and Saika. Furthermore, we present the theoretical framework for L C G T O - L S D F gradient calculations. -1

-1

0097-6156/89/0394-0228$06.00/0 © 1989 American Chemical Society

In The Challenge of d and f Electrons; Salahub, D., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1989.


16. A N D Z E L M ET A L

Spin Density Functional Approach

229

There is growing evidence that local spin density functional (LSDF) theory (1) provides a unified theoretical framework for the study of electronic, geometric, and vibrational structures of solids, surfaces, interfaces, clusters and molecular systems (2,3) encompassing metallic as well as covalent bonding. Although widely used in solid state physics, including semiconductors, transition metals, lanthanides and actinides, density functional methods are still rather rarely applied to problems in chemistry. One of the reasons appears to be the lack of experience with this method i n addressing typical chemical questions such as molecular conformations, vibrations, and reactivity. Such investigations require the capability to calculate accurate analytical derivatives of the total energy with respect to displacements of the atomic nuclei. Evidently, the potential of L S D F gradient calculations has not yet been fully developed. However, i f such a goal could be achieved for molecules and large clusters including transition metals, one would not only have an additional theoretical/computational tool to investigate large molecules with first and second row atoms (which is the current domain of Hartree-Fock calculations) but one could also study organometallic compounds, investigate reactions on metallic surfaces, and simulate large and complex systems such as zeolites and enzyme catalysts at an unprecedented level of detail. In this paper, we present results obtained with a particular molecular orbital implementation of local spin density functional theory using a linear combination of Gaussian-type orbitals (LCGTO's). The results, derived from single-point total energy calculations, illustrate the applicability of the L C G T O - L S D F approach to questions of electronic structure, vibrational properties, and geometries of transition metal clusters and compounds and, in addition, shed light on the computational issues encountered in this kind of large-scale molecular simulation. Furthermore, we demonstrate that the L C G T O implementation allows a compact analytical formulation of energy gradients thus setting the stage for future exploitations of this chemically important feature. The paper is organized in the following way. First the key features of the L C G T O L S D F method are reviewed (3); for a more detailed description of L S D F calculations, the reader is referred to recent reviews (2,3). Four examples, discussed next, highlight the performance of the L C G T O - L S D F approach to predict (1) the electronic structure (photoelectron spectrum) of a transition metal complex, bis(TC-allyl) nickel; (2) chemisorption geometries of carbon atoms on the Ni(100) surface; (3) the vibrational properties o f C O on Pd including the influence of an external electric field; and (4) electronic properties of Z n clusters as a function of cluster size. The last example also illustrates the dependence of computational requirements on the size of the system.

The Linear Combination of Gaussian-Type Orbitals (LCGTO) Implementation Basis Sets. Since the suggestion of Boys (4), Gaussian-type basis functions have become the standard i n quantum chemical ab initio methods. In the Hartree-Fock theory the occurence of four-center two-electron integrals makes this choice a computational necessity. In local density functional theory, on the other hand, a conceptually simpler Hamiltonian gives greater freedom in selecting the variational basis set. For example, in many solid state calculations it has become common practice to use plane waves and augmented plane waves (5) with numerically generated radial functions i n a linearized form (6). A n elegant alternative consists in the use of numerically generated atomic orbitals as basis for molecular orbitals (7). Both approaches using numerical basis functions, the solid state linearized augmented plane wave (LAPW) method (6,8) and the molecular/cluster approach (7,9) have proven extremely useful in carrying out highprecision local spin density functional calculations for solids, surfaces, clusters, and


230

THE CHALLENGE OF d AND f ELECTRONS

molecules (2). On the other hand, the accurate evaluation of analytic energy derivatives within these implementations turns out to be a considerable challenge (10).


In the present studies, we use Gaussian-type basis functions for the molecular orbitals to solve the local density functional equations (11). This choice offers several advantages: (i) there is a wealth of experience from numerous ab initio Hartree-Fock calculations i n using G T O ' s for molecular calculations (12), (ii) computationally, this approach can be implemented in a highly efficient way, as will be discussed below, (iii) the analytic nature of the basis functions opens the possibility for accurate analytic calculations of energy gradients for geometry optimizations and density gradients for non-local corrections (13), and (iv) effective core potentials or model potentials can be readily incorporated (14). Besides the Gaussian basis set for the wavefunctions, there are two other sets of Gaussian expansions used in the present approach, one for the electron density and one for the exchange-correlation potential. The expansion of the electron density is used in the evaluation of Coulomb integrals. Hence the expansion coefficients of the electron density are chosen (lib) such as to minimize the error in the Coulomb energy arising from the difference between the "exact" electron density (i.e. the density originating directly from the wavefunctions) and the fitted electron density. A l l necessary steps to obtain the expansion coefficients of the electron density can be carried out analytically. On the other hand, the expansion coefficients for the exchange-correlation potential have to be obtained numerically by generating the values of the exchange-correlation potential on a grid, which are then used to fit a Gaussian expansioa After this numerical step, the matrix elements of the exchange-correlation potential operator are calculated analytically.

Integral Evaluation. In contrast to Hartree-Fock methods, the L C G T O - L S D F approach requires evaluation of only three-index integrals, thus representing an N algorithm (with Ν being the number of basis functions). A l l examples discussed below, except the Z n clusters, were calculated using the Hermite Gaussian basis originally implemented by Sambe and Felton (11a) and further developed by Dunlap, Connolly and Sabin (lib). In contrast, Cartesian, not Hermite, Gaussians are the most widely used choice i n ab initio quantum chemistry. The scheme of recursive computation of four-index cartesian Gaussian integrals, originally developed by Obara and Saika (15) for the Hartree Fock method, has now been reformulated for the three-index integrals needed i n the present method. A s shown below, a computationally highly efficient scheme results from this approach. 3

Two kinds of three-index integrals are needed, Coulomb integrals of the form (adopting the notation of Obara and Saika (15)) I = [a(l)b(l)llc(2)] c

(1)

where "II" refers to the l / r operator, and overlap-like integrals to calculate exchangecorrelation potentials and energies of the form 1 2

I

xc

= [a(l)b(l)c(l)]

(2)

Here a and b stand for orbital basis functions and c denotes Gaussian functions used in the fitting of the electron density or the exchange-correlation potential. Most of the time is spent in the computation of I . We can rewrite the original Obara and Saika formula (15) in a form suitable for computation of three-index integrals: c


16.

A N D Z E L M ET AL.

231


[àblICc+lj)]^ =

[abllc]

(m+1)

(m+1

+ Vi η N^c) {[ablKc-lj)]^ - (ρ/η) [abll(c-l.)] >} + ^ ( ζ + η ) ) {N.(a) [(a-lj)bllc]< > + N.(b) [a(b-ipilc]^ >) m

n+1

(3)

Here, l is a short-hand notation for a p , D or p function and function c+lj has an angular momentum one order higher than c. The recursive nature of Equation 3 allows to build, for example, integrals with d-functions from integrals with only s- and p- functions. The superscript (m) refers to the order of the incomplete Γ function. W. and Q are related to geometries and values of exponents of Gaussian functions; η , ζ and ρ depend on these exponents and N (c) is a generalized Kronecker delta. For details the reader is referred to the original paper by Obara and Saika (15), Equation 39. i

x

z

{


{

Expansion of the integral with a single function c, instead of a (or b) already cuts the number of arithmetic operations by as much as 30% for d-type integrals. Substantial savings in the time of integral evaluation occur i f we calculate entire groups of integrals with shared exponents in various symmetries. Particularly, we can design basis sets for the electron density and exchange-correlation potential with shared exponents without loosing accuracy in the fitting process (16,17). Hegarty and van der Velde (18) analyzed the number of arithmetic operations necessary to calculate 4-index integrals. The best algorithm (18) requires about 22 operations per integral with d functions. In contrast, the present method requires about 8 operations per 3-index integral. The new formulas for the integral calculations can be efficiently programmed on a vector computer and the algorithm is amenable to parallelization as will be shown below. As for Hartree-Fock calculations, storage of integrals becomes the computational bottleneck for systems with a large number of basis functions. In this case, a "direct S C F ' (19) scheme can be adopted (20). For the L C G T O - L S D F method we deal with a smaller number of integrals (of the order N rather than N as i n Hartree-Fock calculations) and therefore we may use the standard approach for up to about N=1000, as discussed below for the case of large Zn clusters. There is an additional advantage in the direct scheme as we have to calculate the full Hamiltonian matrix only once and then, in each iteration, add the changes to the matrix which correspond to modifications i n the density and the exchange-correlation potential, but only for those matrix elements where this change is greater then a threshhold. During the course of iterations towards selfconsistency, a smaller and smaller number of matrix elements needs to be updated. 3

4

Gradients. The calculation of energy gradients within the L S D F method using localized basis sets has been investigated by a number of researchers (10,21). However, there has been no published formulation for the case of the L C G T O - L S D F method. We will now give an outline of the formulas that have recently been implemented (22). Full details will be published i n due course. In this method, both Coulomb and exchange-correlation energies are calculated analytically once the fitting coefficients for the electron density, p with ρ = e , and exchange-correlation energy-density, e , and potential, μ , are obtained (11). The total energy ( E ^ p ) can then be expressed as r

s

E

LSDF =

Σ

P Ν

h

pq ί p ,

+

*r Pr ^

+

Σ

£

, , ^

>"*

Σ

8

U

« Pr P,

+ »

4

Here Ρ denotes the density matrix, h contains kinetic energy and electron-nuclear attraction operators; [pqllr] and [rllt] are Coulomb repulsion integrals with 3 and 2 indices, respectively, [pqs] denote one-electron 3-index integrals and U is the nuclear-nuclear repulsion term. The form of Equation 4 ensures (lib) that the Coulomb energies are accurate up to second order in the difference between the fitted density and the "exact" density obtained directly from the wavefunction. n



232

In order to obtain the energy gradient of E ^ p by differentiation of Equation 4 with respect to a nuclear coordinate, x, one can, in the first step, closely follow the procedure used in Hartree-Fock theory (23). In addition, derivatives of the density fitting coefficients occur in the present approach. These terms can be eliminated using the equation for the density fitting together with the normalization condition of the total density. At this point, the intermediate gradient formula is given by a E

LSDp/

=ς ρ - vt ς

a x

μ

α

a


+

Z

ν

P r

Pq

{atyax + ς a[pqiir]/ax+ ς e aipqsi/ax j atriiti/ax + a u y a x - ς ^ w atpqi/ax γP r

8

P t

V

s

m

a x

{

Σ

(e

* s - V ) [pqs]} + Σ s

P

Ν

M

Σ,

Βεβχ

[pqs]

(5)

Here, W is an energy-weighted density matrix element as in the Hartree-Fock gradient formula {23). Equation 5 contains two "difficult" terms (the last two), the derivative of density matrix elements and the fitting coefficients e . It turns out that these terms can be eliminated by using the relationship s

ap/ax^

-e

x c

x c

) = p ae /ax

(6)

xc

which can be obtained by differentiation of the fundamental L S D F formula (1) =

d

£

a

£

< xc Ρ > / P = xc

+

Ρ

a

E

xc /

d

7

Ρ

after multiplying both sides by p. In practice, μ and e are obtained through a fitting procedure. The fitted quantities do not correspond strictly to the original density and eq. (6) no longer holds exactly. However, with the fitting basis sets currently used, this approximation appears to be reasonable (22). Assuming therefore the validity of Equation 6, we obtain the following L S D F energy gradient formula which is valid also for the spin-polarized case and in the presence of non-local corrections to the Hamiltonian (20)) xc

a E

LSDF^

X

=

F

HFB

+

F

8

D

()

with F

S

P

+

HFB = pq pq {

Σ

γ Pr U

F is the Hellman-Feynman force plus the correction for the orbital basis set dependence on the nuclear coordinate, x. The term F is specific to the present L C G T O L S D F implementation. Equation 10 is equivalent to R F B

D

F = D

f

^p [a(r)/axll(p-p )] r

f

(11)

with ρ being the "exact" density and p the fitted density. Qearly, in the case of a perfect fit F vanishes. It is important to realize that within the L C G T O implementation, evaluation of L S D F gradients boils down to computations of 2 and 3 index integrals D


16.

ANDZELMETAL.

233


which can be accomplished by the same efficient technique of integral calculation as described above.


Illustrative Examples Electronic Structure and Photoelectron Spectrum of bis(TC-allyl) nickel. The photoelectron (PE) spectrum of bis(TC-allyl) nickel has attracted considerable attention since it was first measured in 1972 by LLoyd and Lynaugh (24). Later experiments (25) established an accurate basis for comparisons with theory. The first theoretical investigation by Veillard et al. (26) using the Hartree Fock approach (cf. Table I) revealed that the interaction between d orbitals of N i and π* orbitals of the allyl group is responsible for most of the bonding. Koopmans* theorem was found to be invalid in this case and total energy A S C F calculations are required. However, even large-scale CI calculations (27) could not correctly identify the first band as ionization from the ligand π-type orbital (7a orbital assuming symmetry of the molecule). A semi-empirical Green function approach (28) (cf. Table I) provided an efficient calculation of all of the ionization potentials (IP's) and gave satisfactory assignment for most of the P E bands. Hancock et al. (29) performed scattered-wave (SW-Xot) calculations of IP's applying the transition state method. Compared with experiment, the calculated spectra are shifted considerably towards lower energies and especially ionizations out of π-type orbitals seem to be in error. u

The present calculations were performed with an all electron version of the L C G T O L S D F method. Triple zeta basis sets for nickel and carbon atoms including polarization functions were employed. Details of basis sets and the method of their optimization are given in Ref. 16. A recent neutron diffraction study (30) revealed a pronounced bending of the anti-hydrogen atoms (by 30°), as a result of their strong repulsion from the nickel atom. First we discuss results obtained assuming planar geometry of the allyl groups as this allows a direct comparisons with the other theoretical calculations of the P E spectra performed so far. The interaction diagram between the 3d and 4s orbitals of N i and the ( C H p fragment together with the resulting orbital levels of bis(TC-allyl) nickel are shown in Fig. 1. The basis set superposition error (BSSE) (17) is corrected for by using the same basis set in all calculations except for the isolated N i atom. This accounts for the splitting in the N i ( C ) case. 3

2

2h

Bonding is caused by the interaction between d-electrons of N i and the π* electrons of the ( C H ) fragment. By symmetry, bonding is allowed within the b and a manifolds. The donation of electrons from occupied 2b orbitals of the metal atom to the unoccupied 5b orbital of the ( C H ) fragment results in the bonding molecular orbital, 5b . Another bonding orbital, 10a is formed as a combination of 7a and 6a orbitals of Dis(TC-allyl) and nickel, respectively. This bonding effect is, however, largely cancelled by the antibonding partner, the 13a orbital. Close to the H O M O (6b ) there are levels of mainly 3

5

2

g

g

g

g

3

g

5

2

g

g

Mulliken charges reveals that the corresponding molecular orbitals (7a and l l b ) have some admixture of metal 4p states. The coupling with 4p states causes a transfer of the ligand*s charge to originally unoccupied 4p orbitals of the metal. This "back donation" mechanism was first discovered in SW-Χα calculations (29) and is confirmed by the present L C G T O - L S D F study. u


u



234

Table I. Experimental and theoretical results for the photoelectron spectrum of bis(TC-allyl) nickel (in eV) The principal orbital character of the ionized electron is added in parenthesis

band exp. (a)

ASCF (b)

1 2

7.8 (π) 8.2 (d)

3

8.6 (d)

4

9.4 (d, π)

5 6 7

10.4 (π) 11.6 (π) 12.7

8

14.2

5.5 5.6 6.8

11.0

ASCF-CI (c)

GF (d)

SW-Xcc (e)

6.4 6.6 6.7

8.7 8.9 9.2 9.2 9.5 10.0

2.5 4.5

10.8

10.9 12.2 13.013.2 15.4 15.6

5.0 5.1 5.5 5.6 6.6 7.9-8.2 9.09.2

LCGTO-LSDF(f) planar bent

7.8 8.0 8.2 8.4

8.1 8.1 8.4 8.4

(π, (d, (d, (d,

9.2 9.5 10.7 11.9 12.112.5 13.6 13.8

9.4 (d, l i a ) 9.9 (d, π, 10.3 (π, lib*) 11.2 (π A 10a ) 11.8- (σ) 13.3 14.3 (σ) 14.6

(a) experimental data : Batich, Réf. 25 (b) Veillard; Rohmer et al., Réf. 26 (c) Moncrieft et al.; Hillier, Réf. 27 (d) semi-empirical Green Function calculation: Bohm and Gleiter, Réf. 28 (e) scattered-wave Χ α calculation: Hancock et al., Réf. 29 (f) present calculations: allyl groups are planar or hydrogen atoms are bent


7a ) 13aJ 12ap u

6bJ

Aj


Figure 1. Orbital diagram of the L S D F one-particle energy levels of isolated nickel, nickeKC^), bis(TC-allyl) nickeKC^) and the ( C j H ^ C L ) fragment. • and Δ indicate spin up and spin down electrons. The three-dimensional structure of bis(7c-allyl)nickel is shown on the left-hand side.


236

THE CHALLENGE O F d A N D f ELECTRONS


Using a planar geometry for the H atoms, the calculated L S D F A S C F values for the P E agree rather well (within 0.2 eV) with experiment (cf. Table I). Using a bent geometry (30), the total energy of the transition metal complex is lowered by 0.9 eV. The corresponding calculated photoelectron spectrum agrees less convincingly with experiment. The main bonding mechanism, i.e. donation from N i 3d to allyl π*, and back donation from allyl π to Ni-4p, remains, however unchanged. The d orbitals of N i are barely affected by this geometry change, and the important shifts are found in the ionizations from ligand π and σ orbitals. It should be noted that the experimental P E was taken in the gas phase while the "bent" structure was deduced from measurements on the crystalline solid (30). Clearly, a full geometry optimization of the isolated complex would settle the question of crystal packing effects and clarify the details of the photoelectron spectrum.

Chemisorption of Carbon Atoms on the Ni(100) Surface. The next example demonstrates the capability of the L C G T O - L S D F method to predict the geometrical structures and vibrational frequencies of carbon atoms chemisorbed on the Ni(100) surface. The C / N i system is of fundamental importance in catalytic petrochemical processes and thus has been the subject of many experimental and theoretical studies. Despite these efforts, many aspects such as the equilibrium position of the C atoms (above or below the surface N i atoms) remained unsettled. Recently, a joint experimental and theoretical study of the chemisorption of carbon on Ni(100) led to a clearer understanding of this system (31). It was found that in the ground state, C is adsorbed in four-fold hollow sites above the surface with a N i - C distance of 1.79 Â, in agreement with the experimental value of 1.75 ± 0.05 Â, obtained from surface extended energy loss fine structure (SEELFS) measurements. Furthermore, the calculated vibrational frequency of the perpendicular mode of the adsorbed C atom, 407 cm is in excellent agreement with the experimental value of 410 c m . _ 1

_ 1

Vibrational Frequency of CO Adsorbed on Pd Clusters. In the next example, we present results for the vibrational properties of a C O molecule on a Pd surface, modelled by a cluster of 14 transition metal atoms (32,33). Furthermore, the influence of an external electric field is investigated using a smaller cluster, P d C O (33,34). In both cases, C O is assumed to be bonded in the bridge position. The chemically inactive Pd core electrons are described by a relativistic model potential as described in detail in Ref. 14. 2

Two coupling modes are considered: for the P d C O cluster the first mode (denoted as h) represents vibration of the rigid C O molecule with respect to the transition metal surface. The second mode is either the Pd-Pd vibration within the plane of Pd surface atoms (r) or out-of-plane stretch of the surface/sub-surface Pd-Pd bond (z). The total energy surfaces (h,r) and (h,z) are calculated for discrete points and then fitted to a fourth order polynomial. Variational and Quantum Monte Carlo (QMC) methods were subsequently applied to calculate the ground and first excited vibrational states of each two-dimensional potential surfaces. The results of the vibrational frequences ω using both the variational and Q M C approach are displayed in Table II. 14

If one assumes a rigid substrate which, at first, seems reasonable because of the high mass of Pd compared with the C O molecule, a frequency of almost 500 cm" is obtained for the vibration of the entire (rigid) C O molecule perpendicular to the surface. This value is in significant disagreement with the experimental value of 340 cm" (35). As can be seen from Table II, a substantial lowering of this "beating" mode occurs due to anharmonic effects in the coupling of the ζ and h modes. In other words, the vibration of the surface/sub-surface Pd-Pd bond stretching couples to the C O beating mode and lowers 1

1


16.

ANDZELMETAL.

237

Spin Density Functional Approach 1

Table II. Vibrational frequencies (cm' ) of the two coupled modes (h,r) and (h,z) obtained from a harmonic (Har), variational (Var) and Quantum Monte Carlo ( Q M Q approach

mode

Har

Var

QMC

240 498

192 521

192 521

196 498

65 402

66 397

(h,r): ω

Γ

Spin Density Functional Approach to the Chemistry of Transition Metal

Recommend Documents