Vibrational spectroscopy using ab initio density-functional techniques

Journal of Molecular Structure 651–653 (2003) 3–17 www.elsevier.com/locate/molstruc

Vibrational spectroscopy using ab initio density-functional techniques Ju¨rgen Hafner* Institut fu¨r Materialphysik and Center for Computational Materials Science, Universita¨t Wien, Sensengasse 8/12, A-1090 Wien, Austria Received 2 October 2002; accepted 23 October 2002

Abstract The basic concepts for ab initio density-functional calculations of vibrational spectra for molecules and solids are reviewed. Applications to a wide range of problems, from the characterization of hydrogenated diamond surfaces to the spectroscopy of molecular adsorbates on metallic surfaces and of molecular crystals to the investigation of molecular processes in zeolites will be reviewed. q 2003 Elsevier Science B.V. All rights reserved. Keywords: Vibrational spectroscopy; Hydrogenated diamond surfaces; Zeolites

1. Introduction The spectroscopy of atomic vibrations in molecules and solids, as well of atomic and molecular adsorbates at solid surfaces and in the cavities of mesoporous materials is now firmly established as a very important tool in surface chemistry and for studying the dynamics of catalytic reactions. While for experimental spectroscopy a wide variety of accurate tools is available (electron energy loss, infrared adsorption, inelastic atom scattering to name only a few), theoretical spectroscopy has remained for a rather long time the domain of empirical and semiempirical valence-force-field techniques. Only very recently first-principles electronic structure calculations based on density-functional * Tel.: þ43-1-4277-51400; fax: þ43-1-4277-9514. E-mail address: [email protected] (J. Hafner).

theory [1,2] have been developed to a point that accurate Hellmann – Feynman forces acting on the atoms may be calculated even for very complex molecular and crystalline system. This has been made possible by the development of new iterative techniques for solving the Schro¨dinger-equation [3 –5] and of novel plane-wave representations of the Hamiltonian [6,7]. The possibility to calculate the interatomic forces with high accuracy and efficiency has opened new routes to a theoretical spectroscopy of atomic vibrations: (i) Within a harmonic approximation, the reaction forces induced by small atomic displacements may be calculated. This allows the determination of the interatomic force constants and furthers the calculation of the vibrational eigenstates via classical Born-van Karman theory. (ii) Beyond a harmonic approximation, ab initio molecular dynamics simulations may be used to calculate vibrational spectra via a Fourier transform of

0022-2860/03/$ - see front matter q 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0022-2860(02)00624-5

4

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

the velocity-autocorrelation functions. This opens even the access to the dynamical properties of liquids. In this paper we shall review very briefly the basic physics underlying the calculation of vibrational spectra in molecules and solids and illustrate state-of the-art application at a series of examples, ranging from carbon-based systems over molecular crystals and molecules adsorbed at surfaces to zeolites.

2. Structure and dynamics using density-functional techniques The basic task in the optimization of the structure of a molecule or solid or in calculation of the vibrational spectrum is the accurate mapping of the potential-energy surface (PES) of the system. Within self-consistent field techniques, the forces acting on the atoms may be determined via the Hellmann– Feynman theorem in terms of the expectation value of the derivative of the Hamiltonian with respect to the atomic coordinates, calculated for the electronic ground state corresponding to the given atomic configuration. For crystals, the stresses on the unitcell may be calculated in a similar way by taking the expectation value of the Hamiltonian with respect to the cell parameters. Using the Hellmann –Feynman forces and stresses, the molecular or crystalline structure may be optimized using static (conjugate gradients, quasi-Newton) or dynamic (molecular dynamics, Monte Carlo) techniques. Exploiting the Hellmann –Feynman forces, the vibrational spectrum may be calculated within the harmonic approximation using the ‘direct method’. There are two variants of this technique. In the ‘frozen phonon’ approach [8], the vibrational energy is calculated as a function of the displacement amplitude in terms of the difference in the energies of the distorted and ideal structure. Evidently this requires the knowledge of the polarization-vectors of each mode from symmetry alone, for crystals the wavelength of the phonons must in addition be compatible with the boundary conditions applied to the supercell used in the calculations. The alternative approach uses the forces on the atoms induced by the displacements of other atoms in the supercell. From these forces calculated via the Hellman – Feynman theorem, certain elements of

the Born –van Karman force-constant matrix may be calculated. Exploiting the point-group symmetry of the molecule or the space-group symmetry of the crystal, the full force-constant matrix may be determined from a minimal number of individual displacements [9 – 11]. For an isolated molecule or an adsorbate –substrate complex with a finite range of interatomic interactions, the vibrational spectrum and the eigenvectors of all modes are directly accessible from the force-constant matrix. For periodic crystals, the dynamical matrix has to be calculated for all wave-vectors within the Brillouin zone—the results for the vibrational eigenmodes are exact (in the sense of a frozen phonon calculation) if the wave-vector is compatible with the periodic boundary conditions applied to the computational supercell. For other wave-vectors, the result must be considered as a sophisticated interpolation. For polar phonons the use of periodic boundary conditions leads to an additional complication: due to incomplete screening, the displacement of charged ions creates an electrostatic multi-pole field that is periodically repeated due to the boundary conditions. This periodic repetition creates an unphysical additional field whose effect is to suppress the splitting of longitudinal and transverse optical modes in the Brillouin-zone center. Methods for constructing compensating fields for polar crystals, solid surfaces and molecules require the calculation of the tensor of the Born effective charges, details may be found in the literature [12 – 14]. Beyond the harmonic approximation, ab initio molecular dynamics may be used to determine the vibrational spectrum as the Fourier-transform of the velocity autocorrelation function calculated in long molecular dynamics run. The drawback of these techniques is that rather long simulation runs are required to minimize the error introduced by the unavoidable truncation of the Fourier integral and that eigenvectors are not available. However, information on the character of individual eigenmodes may be achieved by projecting the velocity-autocorrelation function on selected groups of atoms. Quantumeffects may become important for a correct description of the dynamics of loosely bound protons in hydroxyl groups, hydrogen-bonded systems or soft molecules can be incorporated in different ways: (1) ab initio DFT calculations may be used to construct a low dimensional PES for the protons. On this PES,

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

the Schro¨dinger equation describing the proton dynamics may be integrated numerically using different techniques [15,16]. (2) As a computationally extremely demanding alternative, a quantum-MonteCarlo treatment of the atomic dynamics may be grafted upon the ab initio DFT calculations of configurational total energies [17]. The challenge evidently is to perform the basic operations—i.e. the determination of the electronic ground-state and the calculation of the Hellmann– Feynman forces—for complex systems containing many inequivalent atoms in a large number of different configurations, and with very high accuracy. 2.1. Ab initio DFT calculations The seminal paper of Car and Parrinello demonstrating that the ionic and electronic degrees of freedom of a solid or molecules can be optimized simultaneously was published in 1985 [3]. It has triggered the development of new electronic structure codes opening the possibility to treat systems with ever increasing complexity. Essentially, two different routes have been followed. The construction of the most localized local orbitals, and the use of this basis for a down-folding of the Hamiltonian opens—at least in principle—the way towards a technique where the computational effort scales linearly with the number of the degrees of freedom (the so-called O(N)techniques)—the most advanced of these codes is probably the SIESTA program [18]. The drawback is that the transformation to the most localized basis necessary for the construction of a short-ranged Hamiltonian and to achieve O(N)-scaling unavoidable leads to come loss of accuracy. The alternative is to use a plane-wave basis set and to describe the electron-ion interaction either by pseudopotentials or projector-augmented-wave (PAW) potentials. If pseudopotentials are used for systems containing first-row or transition-elements, so-called ultrasoft pseudopotentials [19,20] have to be used to cut the basis set to a manageable size. The most advanced techniques, however, are based on the projector-augmented waves introduced by Blo¨chl [6] and adapted by Kresse and Joubert [7]. They allow to combine the computational efficiency of pseudopotentials with the accuracy of the most advanced—but

5

computationally much more demanding—all-electron codes. The cornerstone for the extension to very large systems is the efficient diagonalization of the Hamiltonian using a doubly iterative approach. The RMM-DIIS (residual minimization method—direct inversion in the iterative subspace) [4,5,21] performs an iterative optimization of the occupied and the lowest empty eigenvalues, requiring only a minimal number of orthogonalization operations. This greatly improves the scaling of the computational effort with the number of valence electrons in the system. The iterative optimization of eigenvalues has to be complemented by the iterative update of charge- and spin-densities and potentials, using highly optimized mixing-routines [4,5] to build an extremely efficient and highly accurate tool for determining the electronic ground-state of large and complex systems. A further advantage of the plane-wave based techniques is that the Hellmann –Feynman (HF) forces acting on the atoms and the HF-stresses acting on the unit cell may be calculated analytically and very efficiently. With plane-wave techniques, the test of basis-set convergence is also almost trivial—in evident contrast to local basis set methods. With PAW-potentials, absolute convergence can be achieved even for very complex systems. The flagship among the plane-wave techniques currently available is the Vienna-ab initio-simulationpackage VASP developed by Georg Kresse1 [4,5]. Even for the largest systems that can be treated at affordable computational effort (i.e. , 400 atoms for static optimization of the atomic geometry and , 200 atoms for ab initio molecular dynamics simulations over time-spans sufficiently extended to allow to extract high-resolution vibrational spectra—for illustrative examples see [22 –25]) the scaling behavior is still /N a with a , 1:5 – 2; i.e. far from being dominated by the orthogonalization operations. The large system-size accessible is also very close to the break-even point (for comparable accuracy) with the O(N)-techniques whose advantages begin to be felt only for much larger systems. As a final remark: it is a bit ironical that the tremendous development triggered by Car and Parrinello [3] did not follow their 1

Detailed informations on the VASP code may be found under http://cms.mpi.univie.ac.at/CMSPage/main/

6

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

strategy of a pseudo-Newtonian dynamics of the coupled electron-ion system on a near-adiabatic PES, but returned to the classical strategy of a determination of the electronic ground-state using iterative diagonalization and mixing, performing the atomic dynamics on the exact Born – Oppenheimer PES. 2.1.1. Exchange-correlation functionals For DFT calculations, the choice of the exchangecorrelation functional is always a crucial step. It is well known, that the local density approximation (LDA) shows a pronounced tendency towards overbinding (too high cohesive energy, too small lattice constants) even if based on very accurate functionals fitted to quantum-Monte Carlo simulations of the homogeneous electron gas [26,27]. The deficiencies of the LDA may be cured to a large extent by adding nonlocal correction in the form of a generalized gradient approximation (GGA). For the construction of GGA functionals, two distinctly different strategies have been adopted: (i) fitting a parameterized functional to a large data-set of molecular and crystalline structures and binding energies [28,29] and (ii) adapting the functional to as many of the known exact sum rules as possible [30,31]. Not only from the point of view of a purist theoretician, the latter approach is clearly to be preferred as it avoids the uncertainties arising if one moves outside the chosen data base. Recently, attempts have been made to improve GGA functionals by adding terms depending on the kinetic energy density [32 – 34], leading to the so-called meta-GGA’s. However, it is not yet clear whether this represents a systematic improvement. Hybrid-functionals attempt to correct the underestimation of the band-gap inherent in DFT by including a certain fraction of ‘exact’, i.e. Hartree – Fock exchange [35,36] (which notoriously overestimates the gap-width if no correlation is included). Sofar, the hybrid functionals have not been much used in the solid state community. This is due to the difficulty of implementing exact exchange in planewave-based codes, the increased computational effort and the uncertainty whether this will also work for metallic systems. While the gradient-corrected functionals lead to a very good description of systems with metallic, covalent or ionic bonding, their performance for

molecular systems held together by van der Waals forces or hydrogen bonds has been questioned. While we will show below that hydrogen-bonded networks are well described within the GGA, dispersion forces are definitely beyond the DFT. For some systems, the overbinding characteristic of the LDA accidentally leads to binding energies and geometries in reasonable agreement with experiment—but this does not imply that the van der Waals bond is correctly described. Attempts are being made, following ideas put forward by Kohn et al. [37], to use dynamic local polarizabilities calculated within DFT to derive approximate expressions for dispersion forces to be added to the GGA, but it is too early to decide whether this will be successful. [38] 2.2. Vibrational calculations using plane-waves and periodic boundary conditions The computational efficiency of plane-wave codes is largely due to the ease of passing from real- to reciprocal-space representation using fast Fouriertransforms—this requires the use of periodic boundary conditions. Even single-molecules must be put into a periodically repeated supercell, which has to be large enough to avoid intermolecular interactions. Within the harmonic approximation, the normal modes of vibration are calculated from a force constant matrix. Each row of this matrix is calculated from a set of single-point total energy calculations for one atom in the supercell being displaced from its equilibrium position in positive and negative senses by an amount d along one of the Cartesian directions. The force constant fij is given by the central ) difference of the reactions forces Fþ(2 induced by ij the displacement,

Fij ¼

Fijþ 2 Fij2 ›2 E ¼ : ›r i ›r j 2d

ð1Þ

The displacement d should be taken as large as possible, but small enough to remain within the regime of linear restoring forces. For n atoms per molecule or unit cell, in principle 3n total energy calculations are required to complete the dynamical matrix, but exploiting the symmetry allows to reduce the number of calculations. For cubic symmetry in principle a single pair of calculations for an atom

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

displaced in an off-symmetry direction is sufficient [9]. For isolated molecules and for zone-center modes of crystals, the force constant gives the dynamical matrix, divided by the geometric mean of the masses of the interacting atoms, 1 Dij ¼ pffiffiffiffiffiffi Fij ; mi mj

ð2Þ

to calculate the dispersion of off-center phonons, one has to multiply with the phase-factors for a given wavevector q~ . For q~-vectors compatible with the periodic boundary conditions applied to a supercell consisting of several elementary cells, the results of the ab initio force-constant approach are identical with those of a frozen-phonon calculation [8]. The range of the interatomic forces determines the necessary dimensions of the supercell. For covalently bonded systems such as the elemental semiconductors with the diamond structure, a (2 £ 2 £ 2) supercell with altogether 64 atoms is largely sufficient. Problems arise mostly for metallic systems with long-range forces arising from Friedel-oscillations. Here up to ten elementary cells have to be aligned— still results with tractable computational effort may be realized by a combination of supercells elongated along different symmetry directions [39]. For systems built of large molecules, the computational cell may be chosen to be identical with a unit cell. For polar phonons creating an electrical multipole moment, the periodic repetition creates an artificial electrostatic field suppressing the splitting of the longitudinal and transverse optical modes. This fake electrical field must be compensated—the compensating field is expressible in terms of the tensor of the Born effective charges—this requires in principle the calculations of the susceptibility matrix. A Berry-phase approach to the susceptibility has been implemented within VASP [41], but for low-symmetry systems to complete the Born charge-tensor remains a formidable task. For high symmetry where the tensor has only diagonal elements, simplified techniques for the construction of the compensating field have been developed [12 – 14]. The force-constant approach may also be used for the calculation of the vibrational spectra of adsorbate – substrate complexes at solid surfaces [40] or within mesoporous cavities [43]. In this

7

case it is usually sufficient to construct the forceconstant matrix for a subsystem consisting of the adsorbate and those atoms of the substrate interacting with the adsorbate. For aperiodic systems (liquids, amorphous materials) and molecules interacting loosely with their surrounding (extraframework particles in zeolites, weakly chemisorbed or physisorbed molecules), the full anharmonic vibrational spectrum may be calculated via ab initio molecular dynamics [24,25,42 – 44]. The total or partial velocity-autocorrelation function may be determined for a long simulation run—sufficiently long to reduce the cutoff errors in the Fourier transform to determine the vibrational spectrum. The disadvantage is that eigenvectors of the individual modes are not accessible—however an approximate assignment of modes may be achieved by a projection on individual atoms or groups of atoms. In the following, we shall discuss the application of the ‘direct’ approach to a range of examples—if not otherwise stated, all calculations have been performed VASP, using gradient-corrected functionals.

3. Carbon-based systems: from diamond to nanotubes 3.1. Diamond and graphite Carbon-based systems are a very interesting testing ground for ab initio DFT-based vibrational spectroscopy because of the wide variety of structures and bonding characteristics. Fig. 1 shows the phonon dispersion relations for diamond and graphite, compared with the experimental data from inelastic neutron scattering and HREELS experiments. The calculations have been performed for rather modest supercells: a cubic (2 £ 2 £ 2) cell containing 64 atoms for diamond, a hexagonal (3 £ 3 £ 2) cell with 72 atoms for graphite [9]. Nevertheless, the calculations correctly reproduces all the characteristic features that differentiate the phonon spectrum of diamond from that of the other elemental superconductors: (i) the steep increase of the TA modes away from the G-point, (ii) an X1 ðLA; LOÞ mode above the X4 (TO) mode, (iii)

8

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

Fig. 1. Phonon frequencies in diamond (left panel) and graphite (right panel). The full lines show the results of the ab initio calculations, the symbols the frequencies determined by inelastic neutron scattering and high-resolution EELS. After Kresse et al. [9].

the near degeneracy of the S3 modes near the Kpoint, and (iv) the ‘overbinding‘ of the optical dispersion relations, leading to a maximum frequency at q~ – 0: The calculation of the vibrational spectrum of graphite represents a particularly stringent test of the ab initio calculations because of the difficulty to simultaneously reproduce the high optical frequencies associated with deformations of the graphene sheets and the low-lying acoustic and optic modes along the c-direction and the negative, almost quadratic dispersion of the transverse acoustic modes along GM associated with shearing the sheets against each other—still reasonable agreement is achieved. The choice of the exchange-correlation functional deserves a comment: while for diamond LDA and GGA lead to almost equivalent results (lattice ˚ (LDA, GGA, constant a ¼ 3.528, (3.574, 3.567) A experiment, bulk modulus B ¼ 4.60 (4.24, 4.43) Mbar, optical frequency at the zone-center LO(G) ¼ 1310(1290,1315) cm21), for graphite only the LDA leads to an accurate description of the static ˚ (LDA, experlattice properties (a ¼ 2.443, 2.46 A iment), c/a ¼ 1.367(1.35 2 1.365), B ¼ 2.88 (2.86 2 3.19) Mbar) while the GGA leads to completely unrealistic values for the axial ratio as the cohesive forces between the graphitic layers are seriously underestimated—for detailed references see Ref. [9]. However, one should bear in mind that although the overbinding of the LDA compensates for the missing van-der-Waals forces, from the point of view of the physics this is not an entirely satisfactory description.

3.2. Hydrogen-covered diamond surfaces The fabrication of well defined diamond films by chemical-vapor deposition or from a carbon-hydrogen plasma is a prerequisite for the development of hightemperature diamond electronics. During growth, the diamond surface is in contact with an atmosphere containing a high concentration of hydrogen, and the hydrogen-induced surface-modification is known to be important for avoiding graphitization: hydrogen adsorbed on the surface saturates dangling bonds, suppresses surface-reconstruction and helps to preserve the sp3 hybridization of surface-C atoms. On the C(111) and C(110) surfaces, each surface-atom has only one dangling bond. Hence adsorption of a monolayer of hydrogen and the formation of CH—groups is sufficient to restore a bulk-like geometry around the surface atoms. On the C(100) surface, two bonds per surface-atom are cut and a monolayer of hydrogen is not sufficient to restore the (1 £ 1) bulk-terminated geometry. Whether a higher degree of hydrogenation, involving the formation of CH2—dihydride groups is possible has been the subject of an extended debate in the literature. Very recently, extended DFT-investigations by Steckel et al. [48] have investigated the structure and energetics of C(100) surfaces under heavier hydrogen loads. It was shown that when the coverage exceeds one monolayer, it is only marginally stable against the desorption of small hydrocarbon molecules or radicals. However, since desorption is an activated process, locally dihydride-units might exist. It was demonstrated that the dihydride H – C –H bending modes leave a unique signal in the IRAS spectrum at a

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

Fig. 2. Intensities of calculated vibrational frequencies of hydrogenated C(100) surfaces. On the (2 £ 1):H surface, only monohydride CH—units are present, on the (5 £ 1):1.2H and (3 £ 1):1.33H surfaces, mono- and dihydride CH2—units coexist. After Steckel et al. [48].

frequency of about 1475 cm21 which is absent on the clean and monohydrogenated C(100) surfaces (see Fig. 2), but has been reported by Lee and Apai [49] and Smentkowski et al. [50] in HREELS experiments on natural diamond surfaces. Hence the combination of theoretical and experimental spectroscopy allows a unique characterization of the surface which was impossible using diffraction techniques alone. 3.3. Carbon nanotubes During the recent years the discovery that graphene sheets may be rolled up to form ‘nanotubes‘ has attracted much attention. Nanotubes could be important elements to nanomechanics operating with single atoms or molecules. Depending on their diameter and chirality, nanotubes may show metallic or semiconducting behavior. The vibrational properties of carbon nanotubes are interesting from a fundamental point of view, as well as for the characterization of nanotubes samples. While the radial breathing mode can be used to determine the tube-radius, the high-frequency part of the Raman-spectrum (the G-band) is useful for distinguishing metallic and semiconducting nanotubes. Surprisingly, the Raman signal of metallic tubes exhibits a Breit – Wigner – Fano line with a significantly lower average frequency than the corresponding modes in semiconducting tubes [45,46]—see Fig. 3. Dubay et al. [47] have performed ab initio

9

Fig. 3. Phonon frequencies in the G-band of nanotubes as calculated via ab initio DFT (discrete symbols) and zone-folding (lines). Phonons are characterized by their symmetry and the polarization vector: L-parallel to the tube axis, T-perpendicular to the tube axis. Experimental results for the BWF line [46] are shown by diamonds. The dashed line interpolating between the stars highlights the softening of the A1 (L) mode After Kresse et al. [9].

calculations of the vibrational frequencies of nanotubes and demonstrated that the mode-softening is caused by a mechanism resembling a Peierls-distortion: for metallic tubes the atomic displacement characteristic for the totally symmetric longitudinal A1(L) mode leads to an opening of a gap at the Fermi level. The energy gain associated with the opening of the gap is not large enough to overcome the restoring forces and to induce a phase transition, but the oscillatory gap lowers the frequency of this mode compared to graphene and to insulating nanotubes. Since for metallic tubes the density of states at EF is inversely proportional to the diameter of the tubes, the softening also increases with decreasing tube-diameter.

4. Molecules, molecular adsorbates and molecular crystals Ab initio calculations of the vibrational spectra of gas-phase molecules are well-established in quantumchemistry. To explore the properties of molecules adsorbed on the inner or outer surfaces of solids or to study molecular crystals, calculations on extended

10

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

Fig. 4. Optimized geometries of benzene adsorbed on Ni(111) in the bridge and hollow positions. After Mittendorfer and Hafner [52].

systems using periodic boundary conditions are required. Here we will discuss the state of the art at the example of aromatic molecules adsorbed on metallic surfaces and of crystalline bases of RNA and DNA. 4.1. Aromatic molecules adsorbed on metallic surfaces The chemisorption of aromatic molecules such as benzene, pyridine or thiophene on metallic surfaces is of growing interest because of the importance of understanding the interaction and bonding mechanism at the interface between organic films and metals—the motivation reaches from the catalysis of hydrogenation and ring-cracking reactions to molecular electronics. The adsorption of benzene on transitionmetal surfaces has been studied most extensively as a prototype for the interaction between the p-bonded aromatic ring and d-orbitals. However, even in a simple case such as C6H6:Ni(111), experiment could not uniquely determine the stable adsorption geometry: the analysis of the high-resolution LEED spectra by Held et al. [51] found the lowest R-factor for an adsorption at a bridge-site, but because of the strong distortion of the molecular geometry predicted by these simulations, adsorption in the hcp hollow producing a somewhat higher R-factor, but a more

reasonable molecular geometry was adopted as the most likely solution. Mittendorfer and Hafner [52] have studied the geometry, energetics, electronic and vibrational spectra of benzene on all three low-index surfaces of Ni. On the Ni(111) surface, the bridgeposition was found to be slightly more favorable in energy (Ead ¼ 1.0 eV) than either the fcc and hcp hollows (Ead ¼ 0.94(0.91) eV for fcc (hcp), see Fig. 4. Only weak molecular distortions were calculated for both geometries. To calculate the vibrational spectrum of benzene in the gas-phase and adsorbed on Ni(111), ab initio MD

Fig. 5. Vibrational spectrum of gas-phase benzene, calculated by ab initio MD. After Mittendorfer and Hafner [52].

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

11

Table 1 Calculated and experimental values of selected vibrational modes of benzene (in meV)—gas-phase and adsorbed. After Mittendorfer and Hafner [52] Mode

Symm.

Calc. (gas)

Exp. (gas)

n9 n3 n13b n16b n15c n5c n12c n1c

B2u A2g E1u E2g E2g B1u E1u A1g

164 169 187 197 384 384 384 390

164 164 186 198 378 378 378 379

a b c

Calc. Follow

Calc.Bridge

Exp.a

160 167 378 378 378

168 175 374 374 374

175,177 373,372 373,372 373,372

After Ref. [53] and [54]. C– C ring stretching. C– H stretching.

simulations in a canonical ensemble have been performed. The simulation temperature was 300 K, the time-step 0.5 fs, the run was extended over 2.5 ps. To illustrate the resolution achieved with this set-up, the gas-phase spectrum is shown in Fig. 5. The interesting point, however, is the adsorption-induced shift of the C – H and C – C stretching modes which can be compared with EELS investigations of Bertolini et al. [53] and Lehwald et al. [54]. The comparison (Table 1) shows that while for the hollow site a large shift of the C –H modes and no change for the C – C ring-stretching modes is predicted, the EELS experiments and the calculations for the bridgeadsorbed C6H6 agree on a smaller shift of the C –H modes and a corresponding down-shift also of the C – C modes. Hence the vibrational analysis supports the predictions of the total-energy calculations. On a transition-metal surface, the interaction between adsorbate and substrate is of a weak covalent character. On a s,p-bonded metal such as Al, the adsorbate-substrate bonding is much weaker and dominated by the mutual polarization of adsorbate and substrate. Again this is confirmed by vibrational spectroscopy showing that there are only very small shifts of the molecular eigenmodes [55]. 4.2. Molecular vibrations in the crystalline bases of RNA and DNA The crystalline bases of RNA and DNA form a set of organic compounds bonded by extended hydrogen

networks in up to three dimensions. Vibrational spectroscopy of the crystalline bases can be regarded as a means to determine force constants of vibrating molecular units, which could then be transferred to related systems. Computational spectroscopy on these systems is a real tour de force, given the large number of degrees of freedom. In addition, the description of the relatively weak hydrogen bonds is a challenge to DFT calculations. Very recently, Johnson et al. [56 –58] have performed a series of pioneering studies of the vibrational spectra of molecular crystals using periodic DFT calculations. The calculations use standard GGA-functionals and perform a full optimization of the atomic coordinates and of the cell geometry prior to the calculation of the vibrational eigenmodes using the direct force-constant approach. Inelastic neutron scattering spectra were calculated from the normal modes using the CLIMAX program of Kearley [59,60]. The possibility to optimize the structure data turns out to be of great importance for the calculation of the vibrational spectra. Generally the structure of the molecular crystals has been determined by X-ray methods at room temperature and are therefore not consistent with the vibrational spectra measured at low temperature. In addition, proton positions cannot be determined very accurately by X-ray techniques, but they are needed to define the geometry of the hydrogen-bonded network. Although the changes in the internal coordinates and in the cell geometry predicted by the DFT

12

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

 and refinement are quite small ðldrlave # 0:1 AÞ energetically stable structures were always found close to the initial measured structures, they turn out to be quite important for the calculations of the eigenvalues. This is illustrated in Fig. 7 at the example of thymine whose crystal structure with the hydrogen-bonded network is shown in Fig. 6. The measured spectrum is compared with the spectra calculated with and without structure optimization. It is evident that the optimization of both coordinates and cell parameters is necessary to achieve a full agreement with experiment, especially in the low-frequency region (Fig. 7). The agreement between the theoretical and experimental spectra is good enough to establish a one-toone relation between the eigenmodes, such that the

Fig. 7. Measured and calculated INS spectra for thymine: (a) measured spectrum, (b) calculated after optimization of internal coordinates and cell geometry for the low-temperature structure, (c) calculated after optimizing the coordinates only, (d) calculated for the high-temperature structure after coordinate-optimization only. After Plazanet et al. [57].

calculated eigenvectors can be used to assign the character of the eigenmodes. Fig. 8 illustrates the mode assignment for thymine, from the lowest frequencies associated with CH3 torsion and rocking

Fig. 6. Crystal structure of thymine. The dashed lines represent intermolecular O –H hydrogen bonds. Courtesy of M.R. Johnson.

Fig. 8. Vibrational mode assignment for thymine, showing the eigenvectors for a characteristic mode from each band. After Plazanet et al. [57].

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

13

Fig. 9. Reaction of a hydrated Na-gmelinite (containing one water molecule per cavity) with NH4Cl: initial (left) and final configuration (middle), NH4-zeolite after elimination of the hydrated NaCl-complex (right). (a) and (b) show top and side view into the main cavity, respectively. After Benco et al. [62].

modes over ring deformations at medium frequencies to the ‘umbrella’-like deformation of the methylgroup in the highest frequency mode. Modes involving the O-atoms and therefore sensitive to a correct description of the hydrogen bonds occur mostly at intermediate frequencies, but they are also in reasonable agreement with experiment demonstrating that the GGA leads to a reasonable description of hydrogen-bonds. Altogether agreement is good, but not perfect. One of the reasons is that the spectra have been calculated for the G-point only. Future work will have to extend these studies to the full Brillouin-zone. There are indications that in particular the modes influenced by hydrogen-bonds will show appreciable dispersion, leading to a more structured vibrational spectrum.

5. Vibrational spectroscopy of chemical activity and molecular processes in zeolites Zeolites are mesoporous aluminosilicates with multiple applications as molecular sieves and as acid-based catalysts [61]. Ideal purely siliceous

Fig. 10. Vibrational spectra of a Na-zeolite (a), a hydrated Nazeolite with different water-concentrations (b,c), a zeolite contain2 þ ing a Naþ – water–NHþ 4 –Cl complex (c) and a NH4 -zeolite. After Benco et al. [62].

14

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

zeolites consist of a loosely connected network of SiO4-tetrahedra. In natural zeolites, a fraction of the tetrahedrally coordinated Si atoms is replaced by Al and the charged-balance is restored by introducing ‘counter-ions’ (alkali or alkaline earth atoms into the cavities whose valence electrons are transferred to the zeolitic framework). Brønsted acid-sites are produced by exchanging the counterions against H-atoms bound to an oxygen-atom coordinated to Al. Vibrational spectroscopy is the most widely used technique to characterize the chemical activity of the acid-sites. 5.1. Spectroscopy of ion-exchange reactions Acid zeolites can be prepared by treating the natural Na-zeolite with ammonium chloride,

producing an NH4-zeolite and solvated NaCl. Heating in vacuo is used to finally prepare the active Hform [62]. Benco et al. [63] have performed an ab initio MD simulation of this ion-exchange process in gmelinite, a zeolite of moderate complexity and a pore geometry identical to industrially widely used zeolites such as faujasite. Fig. 9 shows snapshots of the most important stages of this process— the initial phase where ammonium chloride is introduced into the hydrated Na-zeolite, the intermediate structure produced by the reaction, and finally the NH4-zeolite after elimination of the solvated NaCl complex. The whole process can be monitored by vibrational spectroscopy—for this system of loosely interacting molecules, ions and the zeolite framework it is very essential to go

Fig. 11. Top view of the (001) surface of mordenite, displaying surface hydroxyl groups saturating dangling bonds. Small dark balls-oxygen atoms, larger balls-tetrahedral sites that may be occupied by either Si or Al. Note that some of these surface hydroxyls form hydrogen bonds. After Ref. [42].

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

beyond the harmonic approximation and to perform a fully dynamic simulation. Fig. 10 shows the simulated vibrational spectra at various stages of the process. At high concentrations of water in the cavities, the O – H stretching (, 3700 – 3900 cm21) and bending modes (, 1550 cm21) are well separated from the eigenmodes of the framework (c). At lower water concentrations, the bending mode is down-shifted due to the strong interaction with the Naþ cation (c). In the intermediate configuration with the Naþ –water – Cl2 – NH4 complex a broad band ranging from , 1900 –2400 cm21 originates mainly from the N – H stretching modes (d). Finally in the NH4-zeolite a characteristic band of H –N –H bending (b) and N –H stretching (s) modes are identified. The shifts compared to the eigenfrequencies of gasphase NHþ 4 are asymmetric, indicating that not all N – H bonds are equivalent in the adsorbed ion—see also the final configuration shown in Fig. 9. Hence vibrational spectroscopy allows to characterize the complex ion-exchange process. 5.2. Spectroscopy of chemically active sites at zeolite surfaces Although there is ample evidence that certain chemical reactions occur mainly on the external surfaces of zeolites, and in particular at the poremouths of the large channels and cavities [64,65], little attention has been paid sofar to the structure of zeolite surfaces and to the presence of chemically active sites different from those present in the bulk. Very recently, Bucko et al. [42] have presented a detailed investigation of the (001) surface of mordenite—a widely used catalyst in industry. It was shown that if the coordinatively unsaturated oxygen atoms at the surface are saturated with hydrogen, the external surface undergoes only a modest relaxation—a top-view of the surface is shown in Fig. 11. Still, the chemical activity is profoundly influenced by the presence of the external surface: (1) The stability of the Brønsted acid sites is influenced by the modified environment of the hydroxyl groups at the free surface. (2) Terminal OH groups form a new class of chemically active sites. (3) Terminal OH groups may form hydrogen bonds with other surface oxygen atoms, leading to stretched O –H

15

Fig. 12. OH-stretching frequencies vs. O– H bond-length Triangles: Brønsted acid sites (open—bulk, full—in positions close to the surface), dots: terminal OH groups. The OH-groups with greater bond-length and lower stretching frequencies form hydrogen bonds with another O atom, cf. Fig. 11. After Ref. [42].

distances and hence modified chemical activity. In Fig. 11 these hydrogen bonds are indicated by dashed lines. Again vibrational spectroscopy is found to be very helpful for the characterization of the chemical activity of these sites. Fig. 12 shows the correlation between the OH-stretching frequencies of the Brønsted OH-groups and of the terminal OH-groups with the length of the O – H bond. It is shown that while the bond-lengths and frequencies of the Brønsted hydroxyls are distributed over a relatively modest interval, those of the terminal hydroxyl groups cover a wide range. In particular, the formation of secondary hydrogen bonds leads to a stretching of the O –H bondlength, a lowering of the eigenfrequency and hence a higher acidity of these sites. 6. Summary In this paper I have attempted to present a broad view on the foundations of vibrational spectroscopy of solids, surfaces, and molecules using ab initio densityfunctional techniques and on a wide range of applications. Although the technique has now reached a rather mature stage and offers many opportunities for fruitful applications, there is ample room for further developments. This concerns the fundaments

16

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17

of DFT (development of improved functionals, solving the notorious gap-problem of DFT and achieving and improved description of dispersion forces) as well as the efficiency of the codes with the aim of accessing even more complex systems.

[24] [25] [26] [27] [28] [29] [30]

Acknowledgements [31]

I would like to acknowledge the important contributions of the members of my research group, and in particular of Georg Kresse and Lubomir Benco, to the projects reviewed in this paper. I would like to thank G. Kresse and M.R. Johnson for the permission to reproduce figures from their papers.

[32] [33] [34] [35] [36]

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]

[19] [20] [21] [22] [23]

W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) 1133. R.O. Jones, O. Gunnarsson, Rev. Mod. Phys. 61 (1989) 5048. R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471. G. Kresse, J. Hafner, Phys. Rev. B 47 (1993) 588. G. Kresse, J. Furthmu¨ller, Phys. Rev. B 54 (1996) 11196. P. Blo¨chl, Phys. Rev. B 50 (1994) 17953. G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758. K. Kunc, R.M. Martin, Phys. Rev. Lett. 48 (1982) 406. G. Kresse, J. Furthmu¨ller, J. Hafner, Europhys. Lett. 32 (1995) 729. W. Frank, C. Elsa¨sser, M. Fa¨hnle, Phys. Rev. Lett. 74 (1995) 1791. K. Parlinski, Z.Q. Li, Y. Kawazoe, Phys. Rev. Lett. 78 (1997) 4063. G. Kern, G. Kresse, J. Hafner, Phys. Rev. B 59 (1999) 8551. J. Neugebauer, M. Scheffler, Phys. Rev. B 46 (1992) 16069. G. Makov, M.C. Payne, Phys. Rev. 51 (1995) 4014. P. Ugliengo, ANHARM —A programm to solve numerically the mononuclear Schro¨dinger equation. A. Eichler, J. Hafner, A. Groß, M. Scheffler, Phys. Rev. B 59 (1999) 13297. M. Benoˆit, D. Marx, M. Parrinello, Nature 392 (1998) 258. Among the O(N) code, the SIESTA code is probably the best suited for actual applications. See.E. Artacho, D. SanchezPortal, P. Ordejon, A. Garcia, J.H. Soler, Phys. Sat. Solidi (b) 215 (1999) 809 and the website http://www.vam.es/ departamentos/ciencias/fismaterial/siesta. D. Vanderbilt, Phys. Rev. B 41 (1990) 7892. G. Kresse, J. Hafner, J. Phys. Condens. Matter 6 (1994) 8245. A.D. Wood, A. Zunger, J. Phys. A: Math. General Phys. 18 (1985) 1343. A. Bogicevic, C. Wolverton, Phys. Rev. 64 (2001) 0141061. G. Kresse, Proc. XIth Intern. Conference on Liquid and Amorphous Metals (LAM11), J. Non-Cryst. Solids, 312-314 (2002) 52.

[37] [38] [39] [40] [41] [42] [43] [44] [45]

[46]

[47] [48] [49] [50] [51] [52] [53] [54] [55] [56]

O. Genser, J. Hafner, J.Phys.: Condensed Matter 13 (2001) 959. O. Genser, J. Hafner, J.Phys.: Condensed Matter 13 (2001) 981. J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048. D. Ceperley, B.J. Alder, Phys. Rev. Lett. 45 (1980) 566. F.A. Hamprecht, A.J. Cohen, D.J. Tozer, N.C. Handy, J. Chem. Phys. 109 (1998) 6264. T. van Voorhis, G.E. Scuseria, J. Chem. Phys. 109 (1998) 400. J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Pedersen, D.J. Singh, C. Fiolhais, Phys. Rev. B 46 (1992) 6671. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. S. Kurth, J.P. Perdew, P. Blaha, Int. J. Quantum Chem. 75 (1999) 889. J.P. Perdew, S. Kurth, A. Zupan, P. Blaha, Phys. Rev. Lett. 82 (1999) 2544. J.P. Perdew, S. Kurth, A. Zupan, P. Blaha, Phys. Rev. Lett. 82 (1999) 5179 Erratum. A.D. Becke, J. Chem. Phys. 98 (1993) 5648. P.J. Stevens, F.J. Devlin, C.F. Chabalowski, M.J. Frisch, J. Chem. Phys. 80 (1994) 11623. W. Kohn, Y. Meir, D.E. Makarov, Phys. Rev. Lett. 80 (1998) 4153. J. Angyan, M. Marsman, private communication. A. Eichler, K.P. Bohnen, W. Reichardt, J. Hafner, Phys. Rev. B 57 (1998) 324. A. Eichler, J. Hafner, Chem. Phys. Lett. 343 (2001) 383. M.Marsman, PhD Thesis, Techn. Univ. Delft 2001, unpublished. T. Bucˇko, L. Benco, Th. Demuth, J. Hafner, J. Chem. Phys. 117 (2002) 7295. L. Benco, Th. Demuth, J. Hafner, F. Hutschka, H. Toulhoat, J. Catalysis 209 (2002) 480. M.P. Allen, D.J. Tildesley, Computer Simulations of Liquids, Oxford Science Press, London, 1987. S.D.M. Brown, A. Jorio, P. Corio, M.S. Dresselhaus, G. Dresselhaus, R. Saito, K. Kneipp, Phys. Rev. B 63 (2001) 155414. A. Jorio, A.G. Souza-Filho, G. Dresselhaus, M.S. Dresselhaus, ¨ nlu¨, B.B. Goldberg, M.A. Pimenta, J.H. A.K. Swan, M.S. U Hafner, C. Lieber, R. Saito, Phys. Rev. B 65 (2002) 155412. O. Dubay, G. Kresse, H. Kuzmany, Phys. Rev. Lett. 88 (2002) 235506. J. Steckel, G. Kresse, J. Hafner, Phys. Rev. B, 66 (2002) 155406. S.T. Lee, G. Apai, Phys. Rev. B 48 (1993) 2684. V.S. Smentkowski, H. Ja¨nsch, M.A. Henderson, J.T. Yates, Surf. Sci. 330 (1995) 207. G. Held, M.P. Bessent, S. Titmuss, D.A. King, J. Chem. Phys. 105 (1996) 11–305. F. Mittendorfer, J. Hafner, Surf. Sci. 472 (2001) 133. J.C. Bertolini, G. Dalmai-Imelik, J. Rousseau, Surf. Sci. 67 (1977) 478. S. Lehwald, H. Ibach, J.E. Demuth, Surf. Sci. 78 (1978) 577. R. Duschek, F. Mittendorfer, R.I.R. Blyth, F.P. Netzer, J. Hafner, M.G. Ramsey, Chem. Phys. Lett. 318 (2000) 43. M. Plazanet, N. Fukushima, M.R. Johnson, A.J. Horsewill, H.P. Trommsdorff, J. Chem. Phys. 115 (2001) 3241.

J. Hafner / Journal of Molecular Structure 651–653 (2003) 3–17 [57] M. Plazanet, N. Fukushima, M.R. Johnson, Chem. Phys. 280 (2002) 53. [58] M.R. Johnson, H.P. Trommsdorff, Chem. Phys. Lett., 364 (2002) 34. [59] G.J. Kearley, J. Chem. Soc.-Faraday Trans. 82 (1986) 41. [60] G.J. Kearley, J. Chem. Soc. Nucl. Instr. Meth. Phys. Res. A 354 (1995) 53. [61] D.W. Breck, Zeolite Molecular Sieves, Wiley, New York, 1974.

17

[62] P. Chu, F.G. Dwyer, in: G.D. Stucky, F.G. Dwyer (Eds.), Intrazeolite Chemistry, ACS Symposium Series, vol. 218, ACS, Washington, DC, 1983, p. 59. [63] L. Benco, Th. Demuth, J. Hafner, F. Hutschka, Micropor. Mesopor. Mater. 42 (2001) 1. [64] S.H. Park, H.K. Ree, Catal. Today 63 (2000) 267. [65] J. Aguilar, A. Corma, F.V. Melo, E. Sastre, Catal. Today 65 (2000) 225.