J. Phys. Chem. 1995,99, 7731-1753
7731
Studies of Excess Electrons in Sodium Chloride Clusters and of Excess Protons in Water Clusters R. N. Barnett, H.-P. Cheng? H. Hakkinen,S and Uzi Landman” School of Physics, Georgia Institute of Technology, Atlanta, Georgia 30332 Received: October 5, 1994@
Electronic structure calculations using the local density functional method with nonlocal norm-conserving pseudopotentials, ab-initio molecular dynamics simulations, and a novel method of all-quantum simulations, combining a quantum path-integral description of the nuclear degrees of freedom with concurrent electronic structure calculations of the Born-Oppenheimer potential energy surface, were employed in investigations of the structure and dynamics of sodium chloride and protonated water clusters. For the sodium chloride clusters stoichiometric (NaCl),, 1 In I4,and nonstoichiometric metallization sequences NaCl,,,, 0 Im I 3, Na3Clm, 0 Im I2, Na2Clm, m = 1 and 2, and Na14C1,4-,,, 1 In I6 and n = 10, were investigated. Analysis of various electronic properties (Kohn-Sham one-electron eigenvalue spectra, ionization potentials, participation ratios of Kohn-Sham orbitals), as well as multiple F-center formation energies, suggests that Na14C114-, clusters, starting from the cuboid Na~C113nanocrystal, can be characterized as Nan(NaC1)~4-n having a “phase-separated” metallic part Nan. The Nal4C19 (or Na&19+) cluster exhibits a face (surface) segregated metallic Na5 (Nas+) over-layer, the stability of which is demonstrated by a molecular dynamics simulation at 660 K. In analogy with color centers in bulk ionic crystals, the excess electrons in the halogendeficient clusters in the various metallization sequences occupy energy levels which are split from the bottom of the unoccupied “conduction band” of the corresponding stoichiometric parent. Analysis of the electronic spatial distributions and participation ratios indicates that the excess electrons are of a more delocalized nature in comparison with the electrons occupying the p-like “valence band”. Using electronic structure, structural optimization, and all-quantum simulations, structures and energetics of H30+, (H20)2H+, (H*O),M+, n = 1 and 2 for M = Li, Na, and NI&+, (NH3)2H+, and the mixed (NH3)(H20)H+ cluster are described and analyzed. The quantum nature of the hydrogens in the protonated water clusters, as well as tunneling enhancement of the inversion isomerization in H30+ at 150 K. is demonstrated and discussed.
1. Introduction
One of the main themes motivating scientific endeavors in the emerging broad field of investigations concerned with the physics and chemistry of finite systems (clusters) is that of sizeevolutionary patterns of materials properties, from the atomic and molecular scale to the condensed-phase regime.Is2 Guiding these investigations are a large number of observations that energetics and structural, thermodynamic, spectroscopic, and chemical properties of clusters depend on their chemical constituents and on the degree (size) and state of aggregati~n.’-~ Furthermore, the manner in which cluster properties evolve as a function of size and the physical origins underlying the sizeevolutionary patterns often depend on the property whose size dependence is being interrogated; Le., for a given material different properties, such as cluster geometries, thermodynamic stability, abundance distributions, electronic and vibrational spectra, dissociation and ionization energies, and chemical reactivity, may exhibit different size-evolutionary patterns as the size of the system is ~ a r i e d . ~ In juxtaposition with such explorations of “global” evolutionary pattem~,~. clusters provide unique opportunities for designing and investigating finite-size analogues of bulk materials systems, whose properties may be “tuned” via variations of the degree of aggregati~n.~ Furthermore, investigations of certain properties and processes in isolated gas-phase clusters are of signifi-
’Present address: Department of Physics, University of Florida, Gainesville, Florida, 32611. Permanent address: Department of Physics, University of Jyvaskyla, 4035 1 Jyvaskyla, Finland. Abstract published in Advance ACS Abstracts, May 1, 1995. @
QQ22-3654/95/2O99-1131$O9 .OO/O
cance in explaining gas-phase ambient phenomena (such as chemical reactions), as well as in attempting to explore in a systematic manner the transition from gas-phase to liquid-phase (or solution-phase) chemical and physical processes. Guided by the above observations we present in this paper results of investigations of two cluster systems. First we discuss size-evolutionary patterns of energetics and structure in small stoichiometric sodium chloride clusters and excess electrons (“metallization”) in nonstoichiometric NanCl, (n > m) which form finite-size analogues6-” to color centers in bulk ionic crystals.I2 Second we present results for protonated small water clusters (H@+ and (H20)2Hc), which are basic molecular ions of importance in certain atmospheric ion-clustering reaction^'^-'^ are detected in flames: moist air,I7 and the D region of the ionosphere,’s-21 and have been the subject of a number of t h e ~ r e t i c a l ~and ~ -e~~~p e r i m e n t a l * investigations, ~-~~ providing insights into the nature of intermolecular hydrogen bonding and proton transfer p r o c e ~ s e s . ~These ~ - ~ ~investigations also allow us to focus on recent developments of theoretical approaches for studies of c l u ~ t e r s . ~ ~In’particular, ~ ~ ~ ~ - in ~ our investigations of alkali halide clusters we employ ab-initio molecular dynamics (MD) simulations (the BO-LDA-MD method39)where the ionic degrees of freedom evolve on the Bom-Oppenheimer (BO) electronic potential energy surface, calculated concurrently via the local density functional (LDA) method with the use of nonlocal pseudopotentials. Furthermore, our studies of the protonated water clusters are the first demonstration of a newly developed method>O where the nuclear degrees of freedom (as well as the electronic ones) are treated quantum mechanically 0 1995 American Chemical Society
Bamett et al.
7732 J. Phys. Chem., Vol. 99,No. 19, 1995
via adaptation of the path-integral technique, thus allowing for “all-quantum” simulations (AQMD). Following a discussion in section 2 of the calculational methods, our results for alkali halide clusters are presented in section 3. Our investigations of protonated water clusters are described in section 4, and a summary of our results is given in section 5. 2. Methods
On a most fundamental level, all properties of matter are of quantum mechanical origin. Nevertheless, physical and/or practical considerations often motivate and, in many cases, justify employment of certain limits of quantum behavior (often referred to as quasi-classical approximations) and use of hybrid methodsi0where some degrees of freedom are treated quantally and other classically (e.g., the Bom-Oppenheimer (BO) approximation based on separation of time scales); in many cases investigations wholly restricted to the classical domain are appropriate (sometimes complemented by estimates of quantum effects via perturbation theory and other methods). Classical molecular dynamics (MD) simulations were introduced over three decades ago as a technique for investigating properties and processes in condensed matter systems, allowing studies of physical systems with refined spatial and temporal r e s o l u t i ~ n . ~ Since ~ . ~ ~its- ~inception$’ ~ the most commonly used form of MD simulation involves numerical integration of the classical equations of motion of many-particle systems, with prescribed model interparticle interaction potentials determined (prior to the simulation) through various methodologies, such as fitting of parametrized functional forms to experimental data or to the results of quantum-chemical calculation^^^ and developing embedding functional scheme^.^^,^^ Such simulations proved to be very useful in a broad range of complex materials systems of various compositions, in different phases and degrees of aggregation, and under equilibrium as well as nonequilibrium condition^.^^^^^-^^ However, such simulations are restricted by the inherent limitations of the interaction potentials employed particularly under conditions which are outside the range of validity of the data-base used for parametrization or fitting of the potentials and in circumstances involving electronic rearrangements (e.g., chemical processes or changes in the nature of bonding). The recent introduction of ab-initio MD scheme^,^^^^^,^' where the classical dynamics of the ions evolves on the concurrently quantum-mechanically calculated ground-state Bom-Oppenheimer potential energy surface, alleviate certain of the aforementioned limitations and open new avenues for explorations of a large number of chemical and physical problems in bulk systems5I and c1usters.8-1’~38-40~52-58 In our calculations we have used a Bom-Oppenheimer local density functional molecular dynamics (BO-LDA-MD) method which we developed recently, specifically constructed for investigations of finite-size systems.39 Since the method has been described in detail elsewhere, we limit ourself here to a brief description (for details see ref 39), followed by presentation of a new all-quantum simulation method (AQMD),@ which unifies the BO-LDA-MD methodology with a quantum pathintegral description of the nuclear degrees of freedom. 2.1. The Born-Oppenheimer Local Density Functional Molecular Dynamics (BO-LDA-MD)Method. In simulations of the classical dynamics of ions on the ground-state BO electronic potential energy surface, the dynamical evolution is generated via integration of the Newtonian equations of motion
where 71denotes the vector position of the Zth ion of charge Z1, the second term on the right is the force on the ion due to interionic Coulomb repulsion, and the first term is the force due to its interaction with the ground-state electronic density denotes collectively the coordinates of all the ions in the system). Evaluation of the electronic force, via the HellmanFeynman theorem,59requires an efficient and accurate method for calculating the electronic energy Eelecfor any ionic configuration, T;}. In our calculations
(m
where Teis the kinetic energy, Eel is the electron-ion interaction energy, and Eee is the electron-electron interaction energy, is evaluated after each classical molecular dynamics integration step (Le., solution of eq l), with the use of the local density functional (LDA) or the local spin density functional (LSD) method. In the Kohn-Sham (KS) formulation of the LSD method,6016’ one solves for the wave functions VjU,which are eigenfunctions of a self-consistent single-particle Schrodinger equation (the KS equation)
(3) where a is the spin index (a = f l ) . The corresponding density for the spin manifold LT is given by
(4) where Gu}are the occupation numbers, with Ej,& = Nelec, and NelWis the number of valence electrons. For a finite system, provided that the ground state is not degenerate, the lowest Nelec orbitals should be occupied; however, if the gap between occupied and unoccupied eigenvalues is very small it may be necessary to have fractional occupation numbers to achieve selfconsistency, and we use a Fermi with kBT I 0.001 au in such instances. The KS potential operator, id of eq 3, is given by the functional derivative of Eelec
and since it depends on the electronic density it must be obtained self-consistently. This is achieved by employing an iterative density mixing scheme in the solution of eqs 3-5. The lowest 2, (n, 2 N E ~ eigenfunctions ~J of the KS equation for a given Vu are obtained by an iterative diagonalization method similar to the block-Davidson methoda which we have developed.39 We tum next to a description of the contributions to the ground-state energy of the valence electrons, eq 2, and the KS potential, eq 5. Te is the kinetic energy of the independent particle KS wavefunctions, given by
where the kinetic energy operator ? = -VI2 is diagonal in momentum space. In our dual-space plane-wave method, this term is evaluated in momentum space (i.e., using a plane wave).
J. Phys. Chem., Vol. 99, No. 19, 1995 7733
Studies of Excess Electrons in Clusters The electron-ion interaction is treated via separable nonlocal
pseudo potential^,^^ (7) where, in real space,
vx3,;3,3‘)= *(I3 - 3,l) d(3 - 3’) +
The nonlocal term (for an ion at the origin) is obtained from a semilocal pseudopotential via the Kleinman-Bylander prescription:66
6Gm(3) Gm(3’)
vy(3,3’)=
6= {hm dr
[ I #(r)]’
(9)
A#(r))-l
where A e ( r ) and &(r) are the semilocal pseudopotential and radial pseudo-wave function, respectively. For the purpose of determining the Kohn-Sham potential operator, eq 5, it is convenient to separate Eel into local and nonlocal parts:
and
The local contribution to
Pois given by
and can be combined with the LSD potential.39 The nonlocal KS potential is
where q“ is given by eqs 9- 11. The electron-electron many-body interaction is given by Eee
= Jd3r [@+(r)-I- @-(r)I[€HF)+ € x c ( 3 ) 1
(14)
where
1 +(3) = ?Jd3J
[e+(?‘’)+ ~-(3’)]/13- 3’1
(15)
is the Hartree energy functional and
E x c ( 3 ) = ~XC[@+(~),@-(~);VQ+(~),V@-(~)I (16) is the exchange-correlation energy functional. In LSD, exc depends only on the densities and is simply the exchangecorrelation energy per electron of a uniform electron gas with
densities Q+O and Q-0; in our calculations we use the Vosko and W i l k parametrization ~~~ of the Ceperly and AldeP8 result. If a gradient correction is added, then exc also depends on the density gradients, V Q + O and V Q - 0 . In our calculations of alkali halide clusters we have used the exchange gradient correction of B e ~ k and e ~ the ~ ~correlation gradient correction of Perdew70in the post-LSD approximation; that is, we calculate the gradient correction (xcg) to the energy non-self-consistently using the (self-consistent) KS-LSD generated density. This approach has been shown to be justified via calculations for a number of molecule^.^' In our optimization of the structures of protonated water clusters we used the generalized gradient correction (GGC) method, where the xcg corrections are incorporated into the electronic structure calculations (and evaluation of forces on the ions) in a self-consistent manner. Finally, the contribution of the electron-electron interaction to the KS potential at the LSD level is given by
where EHO is given by eq 15 and exc and its derivatives are given in ref 67. The gradient correction contributions to the potential have also been p~blished.~~-’l Our method39for solving the electronic structure problem is a dual-space method (both real- and momentum-space representations of the wave functions are used), similar to that developed for band-structure calculationsand ab-initio molecular dynamics with periodic boundary condition^.^^,^^ We emphasize, however, that while we use a plane-wave basis set for finite systems, we do not employ a supercell procedure, Le., there are no periodically repeated images of the system.39 The advantage of this approach is that it enables us to treat systems which have large multipole moments, avoiding the need for the large length scale which would be necessary in a supercell method to make the interaction with images negligible. Since in our method the electronic energy and forces on the nuclei are calculated during a dynamical simulation for each nuclear configuration, we are assured of remaining on the ground-state BO surface throughout, thus allowing a relatively large integration time step At of the nuclear equations of motion (determined by the highest characteristic vibrational frequencies in the system; thus, in simulations of water systems we used At = 0.8 fs, and simulations of alkali metal clusters were performed with At 2-5 fs, using the Gear fifth-order predictor-corrector algorithm4I and conserving energy to of the total energy). We note that since in our method the solutions for the Kohn-Sham (KS) eigenstates at time t “precondition” the self-consistentiterations at a successive time, t At, convergence is achieved faster when the ionic configuration at the later time is close to that at the earlier one. In this sense in our MD simulations the number of Hamiltonian operations in the self-consistent solution of the KS-LDA (or LSD) equations (see discussion in ref 39) depends on the integration time step At of the ionic equations of motion. Consequently, the optimal At is determined by requiring energy conservation in the integration of the ionic equations of motion, thus also determining the time interval between electronic structure evaluations, in juxtaposition with minimization of the number of Hamiltonian operations required for convergence of the solutions to the KS equations. In all our calculations we used nonlocal norm-conserving pseudo potential^^^ for the valence electrons (7 for chlorine, 6 for oxygen, and 1 for sodium and hydrogen).
-
+
Barnett et al.
7734 J. Phys. Chem., Vol. 99, No. 19, 1995 Our calculations for Na,Cl, (m 5 n ) clusters for n 5 4 were performed on three-dimensional grids of n1 x n2 x n3 points, with 25 5 nl,n2,n3 5 45 depending on the cluster being investigated.IIb We used a grid spacing of 0.7&, corresponding to a kinetic energy cutoff of the plane wave basis of E, = 20.1 Ry. In our calculations for larger clusterslla Na14C1, (m < 14) we used a grid of 363 points, with a grid spacing of lao, corresponding to E, = 9.9 Ry. By comparing with calculations with larger cutoffs up to 61.7 Ry we find that the total energy of the NaCl monomer was converged to better than 0.5% with E, = 9.9 Ry and that the optimized bond length (2.36 8, with Ec = 9.9 Ry, 2.3 8, with E, > 20 Ry) and vibrational frequency (362 cm-’) are in good agreement with the experimental (2.36 8, and 366 cm-I). Examining the binding energies (with respect to atoms and ions) we find that these energies are -0.4 eV higher than the experimental values74b (4.22 eV for the atomization and 5.57 eV for the ionic binding energy, respectively), which we attribute to the use of the local density approximation. 2.2. All-Quantum Simulations: The AQMD Method. In all the above-mentioned simulations the nuclei were restricted to behave classically. However, in certain cases, such as systems containing light nuclei (particularly protons) and/or small potential barriers, considerations of the quantum-mechanical nature of the nuclear degrees of freedom are of importance. To this end simulation methods based on the path-integral (PI) formulation of quantum statistical mechanics75 have been de~eloped.’~ According to the Feynman path-integral formulation of quantum statistical mechanics, the partition function
Z = TI(-’”)
(18)
(where /3 = l/k& and Tr denotes a trace) for a system of N interacting distinguishable particles with the Hamiltonian H given by
mass mi each of the factors E(7:fa),Ty1)) approaches &?:fa’ ?la+’)) and the corresponding necklace “collapses” to a single point in space, 7i.This is the classical limit which can similarly be achieved for high temperature, i.e., small /3, with necklaces corresponding to particles of larger mass collapsing first. Since ,I2 = 48 (A2 amu K)/T, we note that for hydrogen A2(H;r) = 48 (A2 K)/T while for oxygen A*(O;r) = 3 (A2 K)/T. Consequently, in a system consisting of oxygens and hydrogens, the oxygens may be assumed classical (i.e., single-bead necklaces, although as we will see below this does not bring a significant saving in computational time). Specializing for a system consisting of NH hydrogen atoms (protons) and NO oxygens, eq 21 can be rewritten as
nn n NH
P
NO
d7-Y)
i=l
dTj e-Bvp(23)
j= 1
a=l
where
1
P
c P
V(7?’, ...,Fc;F1,...,FN0) (24)
a=l
where mH is the proton mass and7j, 1 5 j 5 NO, are the positions of the oxygen atoms (taken as classical particles). The “classical form” of the expression for the partition function in eq 23 is evident. However we iterate that in the limit of P = it yields (rigorously) the value for the quantum partition function. The quantum expectation value for any quantity A m , where T denotes collectively the coordinates of all particles in the “real” physical system under consideration (hydrogens and oxygens in our case) is given in terms of averages over the canonical ensemble of a classical system with a potential Vp (eq 24),
-
The partition function may be written as
Z = lim Zp P-W
1
No
where N
P
N
P
The Boltzmann-weighted canonical averages in eq 25 may be evaluated via a Monte Carlo sampling of the configurational space of the isomorphic system, or (as in our simulations) through averaging over the phase-space trajectories generated in a MD simulation77via integration of the equations of motion obtained from the Hamiltonian,
P
a= 1
In eq 21, E = BfP, and
1
(22) where the characteristic length 2: = h2/3/mi is the thermal wavelength (or uncertainty spreading) of a free particle of mass mi. In this description each quantum particle maps onto an isomorphic “ring polymer” (or cyclic “necklace”), consisting of P pseudoparticles (or “beads”; in general we could have assigned to each particle a different number of beads), held together by harmonic springs. We note that in the limit of large
No
where mH* and mo* are auxiliary (fictitious) masses of the pseudoparticles (beads), chosen so as to enhance proper sampling of the phase space of the isomorphic system. The potential V in eq 24 is general. All PI simulations to date employed prescribed model potentials and have been applied with significant success to various s y s t e m ~including ~~,~~ excess electrons in fluids?8 excess electrons in clusters,79 hydrogen diffusion in metals,80 hydrogen adsorptionsia and diffusion81bon metal surfaces, electron tunneling in biophysical systems,82and quantum However, these simulations
Studies of Excess Electrons in Clusters share siinilar limitations, pertaining to the interaction potentials employed, as those mentioned above in the context of classical simulations. Our new method40 unifies the ab-initio BornOppenheimer MD methodology, where the ground-state electronic potential energy surface is evaluated with the use of density functional theory, with the quantum path-integral description of the nuclear degrees of freedom. Consequently, the method which we developed, where the quantum nature of the nuclei is included (via the path-integral formulation) and the energetics (of the isomorphic system) on the self-consistent BO electronic potential energy surface are concurrently evaluated (via LDA or LSD), may be termed as an all-quantum MD simulation method (AQMD). We remark in this context that while the nuclei are restricted to sample the (evolving) BO ground-state electronic potential energy surface, they are free to thermally populate (with canonical ensemble statistical weights) excited (vibrational-like) states of the nuclear degrees of freedom. The new method, which we illustrate in this paper via application to protonated H30+ and (H20)2H+ clusters, offers new ways for investigations of physical and chemical issues in molecular and condensed matter systems, particularly under circumstances involving light species, such as proton transfer processes, hydrogen bonding, tunneling, molecular liquids (e.g., water and ammonia), ice, compressed hydrogen, quantum liquids, and clusters. Finally, we note from eq 24 that while the harmonic term couples consecutive beads (a and a 1) of particle i, the interaction potential V(7~’,...,7~?1,...~~~) couples only beads with the same index a on different necklaces (corresponding to different particles). Thus in performing an ab-initio BOLSD-MD simulation for an isomorphic system, one needs to evaluate, at each integration step of the (Newtonian) equations of motion for the beads (derived from the Hamiltonian in eq 26), the Bom-Oppenheimer ground-state energies for P sets of coordinates; (7;’),...,+;:?I,... 7 ~ ,...,(7:‘p’, ~ )...,7rj?l, ...7 ~ ~ Consequently, in such all-quantum simulations, where both the electrons and nuclei are treated quantum mechanically, the computational effort is multiplied P-fold in comparison to a simulation where the nuclei are treated classically (since the dominant portion of time in the simulation is spent in solving the electronic structure for given positions of the hydrogen beads (and oxygen nuclei). We additionally remark that the structure of the equations allows for significant gains in computational efficiency via the use of parallel computational architectures. In our simulations of protonated water and ammonia clusters and metal ion hydrated clusters we have employed normconserving nonlocal pseudo potential^^^ for the oxygen, nitrogen, and metal atoms (s- and p-components, with s-nonlocality) and a local pseudopotential for the hydrogens. In structural optimization, plane-wave kinetic energy cutoffs, E,, of 96 and 158 Ry were used. Both LDA calculations and generalized gradient correction (GGC) calculations were performed. In the BO-LDA-MD and AQMD simulations for H30f and (H20)2Hf at 150 K, the BO potential energy surfaces were calculated via LDA with Ec = 98 Ry. In the path-integral MD the number of pseudoparticles P = 8, and an integration time step of At = 0.8 fs was used. Statistical averages were obtained from simulations of up to (1.2 x lo4)&
J. Phys. Chem., Vol. 99, No. 19, 1995 7735
+
3. Stoichiometric and Nonstoichiometric Sodium Chloride Clusters
Investigations of nonstoichiometric ionic (e.g., alkali halide) clusters containing single or multiple excess electrons (e.g., those electrons which substitute C1- anions in NanCl,, 0 5 m 5 n) open a new dimension for studies of excess-electron localization
)
and bonding in cluster^.^^^-"^^^-^^ Particularly interesting is the “metallization” sequence (MS) of an initially stoichiometric ionic cluster (i.e., starting from a neutral stoichiometric ionic cluster with m = n and successively substituting anions by electrons, ultimately resulting in a neutral metallic cluster Nan) which may portray a transition from an insulating to a metallic state in a finite s y ~ t e m . ~ s Furthermore, ~-~~ such nonstoichiometric ionic clusters provide finite-size analogues of bulk color centers,I2 and consequently, systematic investigations of such clusters allow investigations of size-evolutionary pattems of defects in condensed matter s y ~ t e m s . ~ ~ ~In- ’this l context we note that while color centers have been studied extensively in bulk crystalline alkali halidesi2 and transitions from F-center to metallic behavior have been observed in bulk molten alkali halides at high excess-metal concentrations,’0’$’02 only limited theoretica18-11s100and experimental informations9,90-93b,98 on metal-rich alkali halide clusters is available. R e ~ e n t l ywe ~ ~have ~ . ~ investigated NanF, clusters (1 5 n 5 4 and, for a given n, 0 5 m < n; also n = 14, 0 5 m 5 14), revealing certain trends of the metallization process in small clusters. These trends include: (i) the structural and energetic substitutional nature of excess electrons in halogen-deficient alkali halide clusters, (ii) cluster structural transitions (which can involve a change in dimensionality) upon reaching the metallization limit, (iii) electronic structure odd-even and shellclosing effects portrayed in formation energies and ionization potentials, indicating that halide-deficient clusters NanFm(m -= n) may be regarded as composed of a metallic and an ionic (alkali halide) components, i.e., Nan-,(NaF),, and (iv) layer metallization (i.e., segregation of the excess metallic component on a face of the cluster) in an intermediate-size cluster (Na14F9). In our earlier investigation^^-*,^ we have used electronic structure local spin density functional (LSD) calculations in conjunction with classical molecular dynamics (MD) of the ions . the ground-state Bom-Oppenheimer (BO) potential energy on surface (a LSD-BO-MD method39). However, in those earlier ~ t u d i e s ~only . ~ . the ~ excess electrons were considered (i.e., n m electrons in NanF,, n > m ) and their interaction with the halide anions and alkali cations were modeled by a Coulomb repulsive potential and by a local pseudopotential, respectively. The interionic interactions were described via semiempirical parametrized potentials. In the present calculations of sodium chloride clusters,” using the LSD-BO-MD method,39 we employ an accurate ab-initio description of the system, by considering all the valence electrons in a given cluster (i.e., 7 electrons for each chlorine atom and 1 electron for each sodium), interacting with the ions via nonlocal norm-conserving pseudopotentials. Moreover, comparative studies between spin density functional calculations using a local approximation for the exchange-correlation functional (xc) and calculations including exchange-correlation gradient (xcg) corrections have been performed. To investigate the structure, energetics, and bonding in alkali halide clusters and the evolution of these properties with the size of the clusters, we have performed calculations for stoichiometric and nonstoichiometric (halogen-deficient) sodium chloride clusters in a range of sizes. We focus here on stoichiometric (NaCl),, 1 5 n 5 4 clusters and on the metallization sequences NaCl,, 0 5 m 5 4, Na3Cl,, 0 5 m 5 3, Na2Clm,0 5 m 5 2, and Na14C1,, 8 5 m 5 13 and m = 4. The presentation of our results is organized as follows: In section 3.1 we discuss the energetics and structure of stoichiometric clusters (NaCl),, 1 5 n 5 4; the smallest nonstoichiometric cluster, i.e., the NazCl molecule, and nonstoichiometric clusters in the sequence Na3Clm,1 5 m 5 2 , are
7736 J: Phys. Chem., Vol. 99, No. 19, 1995
Barnett et al.
TABLE 1: Energetics of Stoichiometric (NaCI), Clusters: Total Energy, ET, Core-Core Coulomb Energy E,, NaCl Removal Energy A(")(NaCI)= ET[(NaCI),-~]4- ET(NaC1) - E~[(Nacl),], c1 Removal Energy A("'(C1)= E~[(Na,c1,-1)] -k E(C1) ET[(N~,CI,)], Cluster Dissociation Energy to Ions A:.) = nE(CI-) 4- nE(Na+) - ET[(NaCI),], and Cluster Atomization Energy A:) = nE(CI) 4- nE(Na) - ET[(Na,Cl,)~ (NaCI)
A(")(NaCI) (LSD) A(")(NaCI)(xcg) A("'(C1) (LSD) A(n)(CI)(xcg) A?) (LSD) A?) (xcg) A): (LSD) A): (xcg)
5 -2 109.26 -2 120.2 I 24 14.78 I .40 I .34 36. I5 35.89 30.59 29.25
4 - 1687.54 - 1696.28 1639.33 2.5 1 2.22 6.15 5.94 29.06 28.82 24.61 23.52
3 - 1264.72 - 127 1.47 763.80 1.67 1.67 5.76 5.62 20.86 20.87 17.52 16.90
3 - 1264.69 - I27 1.37 809.22 1.64 1.57 5.70 5.53 20.83 20.78 17.49 16.80
2 -842.73 -847.21 34 1.44 2.09 2.03 5.73 5.54 13.49 13.48 11.26 10.82
1 -420.32 -422.59 42.25 4.58 4.40 5.7 5.73 4.58 4.40
For (NaCI)3, (R) and (L) denote ring and ladder (double-chain) isomers, respectively. LSD denotes results obtained using local spin density functional calculations, and xcg denotes results obtained using LSD with exchange correlation gradient corrections. In calculations of A(")(CI), A?', and A): the following electronic energies of the atoms and their ions are used (values in parentheses refer to xcg calculations): E(CI) = -410.53 eV (-412.78 eV); E(CI-) = -414.62 eV (-416.865 eV); E(Na) = -5.205 eV (-5.41 eV), and E(Na+) 0. Energies in eV.
discussed in section 3.2; results for the sequence NaCl,, 1 5 m 5 3, are presented in section 3.3, and in section 3.4 the metallization sequence Nal4Cl,, 8 5 m 5 13 is discussed. 3.1. Stoichiometric Neutral and Ionized Clusters: (NaCI), (1 In i 4). Experimental and theoretical studies of stoichiometric (NaCl), alkali halide clusters indicate that for small n the low-energy isomeric structures of these clusters consist of "ring", "ladder" (or "double chain"), and "chain" motifs.Io3As the number of molecules in the cluster increases, energetically optimal configurations based on the rock salt (cubic) structural motif become prevalent. Results for the energetics of (NaCl),, 1 i n 5 4, clusters are given in Table 1, and the corresponding minimum energy structures are shown in Figure 1. Also included in Table 1 are results for (NaC1)5, whose structure, of D2h symmetry, is shown in Figure 2. Of the minimum energy structures shown in Figure 1 those for (NaC1)S and (NaC1)4 are three-dimensional (3D), while the minimum energy structures of (NaCl), for n i 3 are twodimensional (2D), with D3h and C2" symmetries for the ring (R) and ladder (L) isomers of (NaCl)3 and a D2h symmetry for the (NaC1)2 cluster. Our results, which are in good correspondence with experimental (when available) and other ab-initio calculated96Jw values, show (see Figure 1) that in (NaCl),, n 3 2, clusters the distances between neighboring alkali and halogen atoms, dNa-CI, are increased compared to that in the monomer molecule. Thus while in the monomer we find dNa-CI = 2.386 8, (compared to an experimental value of 2.361 AIos and 2.388 f 0.008 AIM), dNa-C1 = 2.542 8, in (NaC1)2 (compared to an experimental valueIMof 2.584 f 0.034 A), dNa-CI = 2.524 8, in the R-isomer of (NaCl)3 (averaged over all the bonds), and dNa-C1 = 2.640 8, in (NaC1)4. We note that these values are smaller than the measured interionic distance in the bulk cubic crystalIo7where dNa-CI = 2.814 A. Inspection of the energetics of the (NaCl),, 1 i n i 4, clusters displayed in Table 1 reveals several trends: (i) While the permolecule Coulomb core-core repulsion energy becomes larger as n increases, the total energy (electronic plus core-core Coulomb repulsion) per molecule (ET/n, from Table 1) is approximately the same for 1 5 n 5 4. (ii) The energy for removal of a halogen atom from a given cluster (A(,)(Cl) in Table 1) decreases with n, while the energy required to remove an NaCl molecule from the cluster (A(,)(NaCl) = E[(Nacl),-~] = E[(NaCl),]) is smallest for n = 3. The nonmonotonic
1
2
Figure 1. Optimal geometries for (NaCI),, 2 5 n 5 4, clusters, with the optimization performed at the LSD level (Le., local exchangecorrelation functional). Large and small spheres correspond to chlorine and sodium atoms, respectively. For the T d optimal (NaC1)4 cluster, dNa-CI = 2.64 A, dc1-a = 3.82 A, L(123) = 82.1', and L(234) = 97.4' ( L ( i j k ) denotes the angle subtended by atoms i, j, and k, withj at the apex). For the lowest energy 2D (NaC1)3 ring isomer of D3h symmetry, dNaCl= 2.52 A, L(123) = 104.4', and L(234) = 135.6'. For the other 2,D (NaC1)3 lad$er isomer ofoC2, symmetry, d12= 2.58 A, d23 = 2.49 A, d34 = 2.56 A, d14 = 2.93 A, d24 = 4.00 A, L(123) = 85.5', L(234) = 104.4', L(341) = 77.3", and L(412) = 92.8'. For the 2D (NaC1)2 cluster of Da symmetry,d12= 2.54 A, L(123) = 102.2', and L(234) = 77.8'. In the monomer NaCl molecule dNaCl = 2.39 A.
variation of A(")(NaCl)with n may be correlated with a change in the dimensionality of the ground-state structures; i.e., removing an NaCl molecule from (NaC1)4 to form (NaC1)3 (that is, A(4)(NaCl)) involves a transformation from a threedimensional (3D) structure to a two-dimensional (2D) one, with the alkali and halide ions better coordinated in the 3D structure. The small value of A(3)(NaCl)is associated with a transformation between two 2D clusters, with the ions in the final one, (NaC1)2, somewhat better coordinated than in the initial (NaCl)3 cluster. (iii) The binding energy per ion pair (Aj"'/n, from Table 1) increases monotonically, starting from 5.73 eV for the NaCl molecule (compared to an experimental valueIo3 of 5.57 eV)
Studies of Excess Electrons in Clusters
J. Phys. Chem., Vol. 99, No. 19, 1995 7737
1 4
Figure 2. Optimal structure of (NaC1)S. Large and small spheres correspond to chlorine and sodium atoms. Comparison of the structure to that of (NaC1)4, shown in the upper left of Figure 1, indicates that the addition of the NaCl molecule (atoms 2 and 3 in this figure) distorts the T d cuboid structure of (NaC1)4. S o y e of the inter$tomic distances and angles in the cluster are d12= 2.62 A, d23 = 2.47 A, d34 = 2.59 A, d14= 3.46 A, and L(123) = 112.6", L(234) = 89.5", L(341) = 87.6", and L(412) = 70.3".
2
7+-4@3
1
1
4
Figure 3. Optimal structures of (NaCl),+ (1 In I4) clusters. The structural pa!ameters of the distorted cuboid structure of (NaC1)4+ are dl7 = 2.73 A, d48 = 2.76 A, and all the other noted bonds between sodium and chlorine atoms are of leng!h 2.65 f O.Olo A. Some of the interchlorine distances are d47 = 3.19 A, d26 = 3.42 A, and d46 = 3.65 A. The (NaC1)3+cluster was obtained via energy minimization starting from that of the ionized L isomer of (NaC1)3 (see bottom right of Figure 1). The structural parameters for the C2" structure of (NaC1)3+are d12 = 2.67 A, d23 = 2.57 A, d34 = 2.70 A, d14 = 2.70 A, d24 = 2.97 A, L(123) = 114.5", L(234) = 68.7", L(341) = 109.4", and L(412) = 67.4". The structural parameters of the rhombus geometry of (NaC1)2+ are d12 = 2.67 A, dl3 = 2.59 A, and L(123) = 122". In the ionized NaCl molecuk d ~ = ~2.82~ A.l
and achieving a value of 7.205 eV for (NaC1)4 which is 90% of the "lattice energy" in a bulki0*sodium chloride crystal (7.98 eV). (iv) The (NaC1)4 cluster is the basic building block for larger cuboid sodium chloride clusters. As seen from Figure 2 addition of a NaCl molecule to (NaC1)4 to form the (NaC1)s cluster, while distorting the structure of the (NaC1)4 unit, maintains its integrity. We note that the energy for removal of a NaCl from the (NaC1)4 cluster is significantly larger than the corresponding energy for (NaC1)s (see A(")(NaCl)in Table 1). In addition the binding energy per ion pair for (NaC1)5 is somewhat lower (7.18 eV) than that of the (NaC1)4 cluster. The structure of the ionized (NaC1)4+ cluster bears a certain similarity to that of the neutral parent; i.e., the structure of (NaC1)4+ is a distorted cuboid (see Figure 3) with variable dNaCl bond distances. We note in particular that in the optimized structure some of the inter-chlorine dcl-cl distances shorten significantly, the shortest being 3.19 A, in comparison to their value in the neutral (NaC1)4 parent where ~ C I - C I = 3.82 A. Ionization of smaller clusters may induce even larger structural variations. As an example we show in Figure 3 the minimum energy structure of (NaC1)3+ obtained via structural optimization of the ionized ladder neutral isomer (see bottom right in Figure 1). Note in particular the significant reduction of one of the inter-chlorine distances (d24 = 2.97 A), in comparison to its value (4.00 A) in the neutral (NaC1)3 cluster. For the minimum energy structure of (NaC1)2+ shown in Figure 3 we find a 2D rhombus (D2h symmetry) with dNaCl =
Figure 4. Optimal structure of Na2CI (C2" symmetry) and contours of the electron density distribution of the highest occupied Kohn-Sham orbital (populated by the excess electron), plotted in the plane of the molecule. Empty circles correspond to sodium atoms, and the filled circle represents the chlorine atom. dNaCl = 2.54 A, and L(NaCINa) = 80.7" (compare to Na2C12 shown in bottom left of Figure 1). Upon ionization the resulting cationic Na2CI+ cluster has a linear ground state structure (Na-Cl-Nt)+ of D,h symmetry, with ~ N ~ = C I 2.50 8, (in NaCl+, ~ N ~ = C I2.82 A).
2.67 8, and an angle L(NaC1Na) = 122.0" (compared to = 2.54 8, and L(NaC1Na) = 77.8" in the neutral (NaC1)2 cluster, see Figure 1). Interestingly the distance dcl-cl = 2.59 %, in the ionized (NaC1)2+ cluster (compared to 3.96 8, in the neutral (NaC1)2 cluster) is close to our calculated equilibrium distance of 2.56 A in C4-. This result is at variance with that obtained by a Hartree-Fock calc~lation*~ where a minimum energy structure of C2" symmetry was found (which may be viewed as a linear [Na-Cl-Na]+ cluster with a chlorine atom weakly bound to it, forming a 2D triangular structure). We have verified our result for (NaC1)2+ via minimization of the energy of the cluster starting from the (C2", 'AI,) structure suggested in ref 85, resulting in the aforementioned D2h rhombus geometry. Furthermore, we performed a constrained energy minimization of (NaC1)2+, where the dNaC1 distance in a linear arrangement of NaClNa was kept constant at 2.50 8, (which is the distance we determined for the optimal linear structure of (Na2C1)+), and the chlorine was placed initially at the apex of a equilateral triangle at a distance of 3.45 8, (as in the C2", 'A,, structure of ref 85). Upon (constrained) minimization of the energy with respect to the position of the extra chlorine, a distance ~ C I - C I= 2.52 A was obtained, again close to the internuclear distance in C12-), and the total energy of the so optimized triangular C2" configuration is higher by -0.5 eV than our optimized D2h rhombus configuration. This leads us to suggest that the groundstate (NaC1)2+ cluster may be described symbolically as Na+(C12-)Na+, in a planar D2h rhombus geometry. In this context it is interesting to note that the above tendency, which we have observed in (Na,Cl,)+ (n 5 4) for shortening of the bond distance between a pair of chlorine atoms, is reminiscent of the structural relaxation associated with self-trapped excitons in alkali halide crystals,IW involving formation of a C12molecular ion. 3.2. Na2C1, and Na3C1, (m = 1 and 2). (i) NaZCZ. The smallest halide-deficient alkali halide specie is the Na2Cl molecule. In its ground electronic state the Na2Cl molecule is bent, with an angle L(NaC1Na) = 81.3" (see Figure 4), compared to 77.8" in (NaC1)2. The electronic charge density contours of the highest occupied KS orbital (corresponding to the excess electron, see Figure 5) shown in Figure 4 illustrate a somewhat delocalized distribution with a concentration of charge density opposite the chlorine along the bisector of the NaClNa angle. The distance dNaCl in the molecule (2.506 A)
Barnett et al.
7738 J. Phys. Chem., Vol. 99, No. 19, 1995
w
-20
Na,Cl,
-
Na,C1
-25 Figure 5. LSD Kohn-Sham energy levels for Na2C1, Na3Cl, Na3Cl+, and the L and R isomers of Na3C12. For each cluster, levels for both spin orientations are shown, with up-spin on the left and down-spin on the right (when these are degenerate a long horizontal line is shown). The shortest lines correspond to unoccupied levels and the longer ones to occupied levels.
..
I'
h
-Na-Cl+
TABLE 2: Energetics of Na3Clm,m = 1 and 2, and Na2Cl Clusters" Na3C12(L) m ET(LSD) ET(xcg) E, A("l(NaC1) (LSD) A("'(NaC1) (xcg) A(,)(Cl) (LSD) A(")(Cl) (xcg) A("(Na) (LSD) A(")(Na) (xcg) A?) (LSD) AT) (xcg) A ): (LSD) AYJ (xcg)
2 -848.46 -853.06 400.72 1.68 1.59 5.09 4.94 0.53 0.44 14.02 ' 13.92 11.77 11.27
Na3C12(R) 2 -848.41 -853.06 364.62 1.62 1.59 5.04 4.94 0.47 0.44 13.96 13.92 11.72 11.27
Na3C1
Na2Cl
1 -432.84 -435.34 113.29 1.30 1.21 5.49 5.26 1.17 1.05 7.80 7.65 6.68 6.33
1 -426.47 -428.88 83.90 0.95 0.88 4.72 4.55 0.95 0.88 6.64 6.60 5.52 5.28
200
-+
TABLE 3: Vertical (VIP)and Adiabatic (aIP) Ionization Potentials and Reorganization Energies ER = VIP - aIP, for Na3Clm,m = 1 and 2 and NazCl Clusters"
~~
Na3C12 (L) Na3C12 (R) Na3C1 NazCl
4.89 5.58 5.49 4.90
6.00 5.71 5.51 5.02
4.24 4.18 5.04 4.12
4.35 4.27 5.09 4.24
~
0.65 1.40 0.45 0.78
~~
~
0.65 1.44 0.42 0.78
a LSD denotes results obtained using the local spin density functional calculations, and xcg denotes results obtained using LSD with exchangecorrelation gradient corrections. Energies in eV.
is shorter than in the (NaC1)z dimer (2.54 A). From Table 2 we observe that the energy for removal of a Na atom from NazC1 (0.88 eV, calculated with xcg corrections) is in adequate agreement with the experimentally determineds7value (0.78 & 0.03 eV), and other calculation^.^^ Additionally our calculated energies of the process NazC1- Naz+ C1- are 5.68 eV (using LSD) and 5.94 eV (including xcg corrections) compared to an experimental value of 5.51 f 0.03 eV.87 Upon ionization the quadrilateral NazCl cluster transforms into a linear structure. This structural change is reflected in a large reorganization energy E R (0.78 eV, see Table 3), with the calculated value of the adiabatic ionization potential, aP = 4.12
+
300 w (cm-')
350
400
Figure 6. Vibrational densities of states obtained via Fourier transformation of the velocity correlation functions recorded in BomOppenheimer molecular dynamics simulations at -50 K. Results are shown for (NaC1)2 (top) and for NaZCl and Na2W (bottom). Also indicated is the calculated vibrational frequency of the NaCl molecule (362 cm-I). For Na2Cl and Na2Cl+ only frequencies in the stretch range are shown. For Na2C1 the symmetric and antisymmetric stretch frequencies are 256 and 350 cm-', respectively, and the corresponding frequencies for Na@ are shifted downward to 242 and 294 cm-I.
'For the Na3C12 cluster, energies for the ladder (L) and ring (R) isomers are given. The definitions of all entries are as given in Table 1. A(")(Na) is the sodium removal energy which for a cluster Na,Cl, (m < n) was calculated as A")(Na) = ,??~[Na,-~cl,]+ E(Na) &[Na,Cl,]. For the cluster dissociation energy into ions we calculated for Na,Cl, ( m < n), APJ = mE(C1-) mE(Nat) (n - m)E(Na) &[(Na,&)]. Energies in eV.
+
250
eV (4.24 eV), determined without (with) xcg corrections, close to the experimentally measureds7 ionization potential (4.1 & 0.1 eV). Further characterization of the NazCl cluster and insight into its nature may be obtained via molecular dynamics simulations of the nuclear motions on the concurrently calculated electronic BO potential energy surface. The results of such simulation for (NaCl)z, NazC1, and Na@, performed at a kinetic ionic temperature of -50 K, are shown in Figure 6, where we display the vibrational densities of states, evaluated via Fourier transformation of the ionic velocity correlation functions obtained from 10 ps simulations (a time step of 3.1 fs was used in integration of the equations of motion of the ions). First we note that our calculated vibrational frequency for the NaCl molecule (we= 362 cm-') agrees well with the measured74avalue, 364.6 cm-I. The (NaC1)z dimer cluster has 6 normal modes. The low-frequency peak in the density of states shown in Figure 6 corresponds to rotations of the molecule. The location of the other peaks are in good correspondence with recent theoretical vibrational analysis of the molecule, based on quantum-chemical calculations. In order of increasing frequency they have been assigned as96B1, (89 cm-I), A, (128 cm-'1, B3u and B I , (or ~ ~B2u and B1,"O) under the peak centered at w = 224 cm-I, A, (252 cm-I), and B2u96 (or B3,'I0) for the peak at w = 274 cm-I. The vibrational characteristics of NaZCl and NaZCl+ in the stretch region, shown in the bottom panel of Figure 6, correspond to the symmetric (lower frequency) and antisymmetric stretch modes, exhibiting shifts to higher frequencies upon ionization of the Na2C1 molecule. A similar trend was obtained by us previouslyI0 in simulations were the interionic interactions were modeled by semiempirical Born-Huang interaction potentials, where we have also found that the lower frequency bending mode of NazCl (w x 105 cm-I) shifts to lower frequency (wsz 85 cm-I) upon ionization.I0 The pattern
Studies of Excess Electrons in Clusters
.2
.2
0 3
J. Phys. Chem., Vol. 99,No. 19,1995 7739
0 3
Figure 7. Optimal structures of the C2u ring (R) and C, ladder (L) isomers of Na3C12 (top and bottom left, respectively) and of the corresponding (Na3C12)+ ions (top and bottom right). Contours of the excess electron distribution plotted in the plane of the structures are shown for the neutral clusters, exhibiting localization in the vicinity of the missing C1 atom (compare with the corresponding structures of the R and L isomers of the stoichiometric parent (NaC1)3 cluster, shown in Figure 1). Empty and filled circles represent sodium and chlori!e atoms, respective1 Structural data: Na3C12 (R), d12 = d45 = 2.54 A, d23 = d34 = 2.52 L(123) = 4 3 4 5 ) = 98.6", and L(234) = 134.9'; linear Na3C12+, d12 = 2.45 A, and d23 = 2.53 A; Na3C12 (L), d12 = 2.60 A, d23 = 2.51 A, d34 = 2.59 A, d45 = 2.62 A, d14 = 2.12 A, L(123) = 81.9", L(234) = 102.9', L(341) = 78.1", L(412) = 97.0", and L(145) = 17.9'; C2" Na3C12+, d12 = 2.49 A, d34 = 2.14 A, d45 = 2.51 A, L(123) = 90.5", L(234) = 94.4", and L(341) = 80.7".
1;
of these vibrational frequency shifts correlates with partial screening of the sodium atoms by the excess electron, causing a softening of the stretch vibrations in the neutral NazC1 molecule compared to the ionized molecule, while localization of the excess electron opposite the chlorine in the 2D quadrilateral structure stiffens the bending vibration of the neutral, as compared to that in the ionized molecule. (ii) Na3C1, (m = 1 and 2). The two lowest energy isomers of Na3C12, shown in Figure 7, are two-dimensional structures which, in analogy with the isomers of the (NaC1)3 stoichiometric cluster, we name as ring (R)and ladder (L) isomers (of CzV and C, symmetry, respectively). As seen from Figure 7 the excess electron in these nonstoichiometric clusters is localized in the region of the missing chlorine (compare with the structures of (NaC1)3 in Figure 1). Starting from the corresponding isomer of (NaC1)3 the energies required to form the Na3C12 R and L isomers by removal of a C1 atom are 5.62 and 5.52 eV, respectively (calculated with xcg corrections). While the overall energetics of the two isomers of Na3C1~is rather close (note, however, the larger Coulomb repulsion for the L isomer; see Ec in Table 2; see also the KS energies for the two isomers in Figure 5), they may be distinguished via measurement of their vertical ionization potentials (see VIP in Table 3). Upon ionization, the R isomer transforms to an optimal linear D,h structure of Na3C1zC (see Figure 7), with the transformation accompanied by a large reorganization energy (see E R = 1.44 eV in Table 3). A smaller reorganization energy is associated with the structural transformation of the L isomer (see Figure 7) after ionization ( E R = 0.65 eV). Removal of a chlorine atom from Na3C12 results in a Na3C1 cluster whose optimal structure is of CzVsymmetry, resembling that of the R isomer of Na2C13, with the density distribution of the two excess electrons delocalized over the three sodiums, as well as populating a p-like orbital of the chlorine (see Figure 8). While the energy for removal of a NaCl molecule from the cluster (Le., Na3C1- Na2 NaC1) is energetically less costly
+
Figure 8. Optimal C2, structures of the Na3Cl (left) and (Na3Cl)' (right) clusters. Empty and filled circles correspond to sodium and chlorine atoms, respectively. Superimposed are contours of the distribution of the two excess electrons for Na3C1and that of the single excess electron for (Na3Cl)+, plott5d in the plane of the cluster. For Na3C1, d12 = 2.54 A, d23 = 3.29 A, d24 = 3.16 p\, L(214) = 77.0°, L(234) = 51.4", and L(324) = 61.3". For (Na3Cl)+, d12 = 2.47 A, d23 = 2.54 A, d24 = 4.36 A, L(214) = 123.9", L(234) = 57.4". and L(324) = 61.3'.
-
(1.21 eV, calculated with xcg corrections, see Table 2) than that for the Na3C12 cluster (Le., Na3C12 Na2Cl 4- NaC1, requiring an energy of 1.59 eV), partially due to the stability of the closed-shell Naz product, the reverse holds for removal of a chlorine atom (Le., Na3Cl Na3 C1 requiring an energy of 5.26 eV, as compared to Na3C1~ Na3C1+ C1 at a cost of 4.94eV). Similarly, the energy for the dissociation reaction Na3Cl -* Na2Cl Na is 1.05 eV, as compared to 0.44 eV required in NasCh (NaC1)z Na (both calculated with xcg corrections). Upon ionization of the Na3Cl cluster the KS energies shift downward (see Figure 5 ) with the KS level occupied by the single excess electron in Na3C1+ shifted to an energy close to that of the p-like levels in the neutral Na3Cl cluster. We note the larger adiabatic ionization potential for Na3Cl in comparison to NazC1, while the vIPs for the two clusters are closer in value, due to the similarity of the structure of the Na@ ion to that of the neutral parent (see Figure 8) and consequently relatively smaller reorganization energy upon ionization (see Table 3). 3.3. Metallization Sequence: NaCl,,, (0 Im < 4). As discussed in the Introduction the main motivations for investigations of nonstoichiometric alkali halide clusters pertain to the nature of F-centers and aggregate excess electron centers in finite systems and the systematics of insulator-to-metal transitions upon metallization in a sequence of cluster^.^,*-^^ Such a metallization sequence, MS, (Na,Cl,, 0 Im In ) starts from a stoichiometric ionic salt, Na,Cl, (i.e., m = n), and ends as a metal particle, Nan (Le., m = 0). As mentioned earlier these issues were discussed by us p r e v i o u ~ l yusing ~ ~ ~an~approximate ~ description of the electronic structure, in a study of sodium fluoride clusters. In the following we readdress them using our current calculations for small sodium chloride clusters. Results pertaining to the energetics of the clusters in the MS NaCl, (1 5 m < 4)are given in Tables 4 and 5 and in Figure 9, and structural information is shown in Figures 10 and 11 (in Figure 11 excess electron density contours are superimposed on the ionic structures). Excess metal nonstoichiometric alkali halide clusters may be regarded as being composed of a "metallic" component and a molecular ionic one, represented symbolically by Na,Cl, Nan-,(NaCl),. We note first the significant reduction in core-core repulsion upon successive removal of C1 atoms (i.e., decreasing m). Also the total energy per NaCl molecule in the cluster (i.e., ETh), as well as the energy for removal of a NaCl molecule from the cluster, A(")(NaCl) = &(Na3C1,-1) f ET(NaC1) - ET(N~CI,), decrease as m varies from m = 4 to m = l . The excess-electron defect formation energy A(")(Cl) decreases in the sequence in a nonmonotonic manner, being largest for NaC12. A similar behavior is exhibited by the ionization potentials (particularly the adiabatic ionization potential aIP, see Table 5). Such oddeven oscillations with (n - m) of Acm)(C1)and the ionization
- -+
+
-
+
Barnett et al.
7140 J. Phys. Chem., Vol. 99, No. 19, 1995
TABLE 4: Energetics of Nonstoichiometric NqCl,, n 5 3. Total Energy ET, Core-Core Coulomb Energy E,, NaCl Removal Energy A(")(NaCl),C1 Removal Energy A(")(CI),Na Removal Energy A("(Na), Dissociation Energy to Ions A?', and Atomization Energy AT'a m
ET(LSD) ET (xcg) EC A("(NaC1) (LSD) A("l(NaC1) (xcg) A("'(C1) (LSD) A(")(Cl) (xcg) A(")(Na)(LSD) A(")(Na) (xcg) A?) (LSD) A?) (xcg) A?) (LSD) A?) (xcg) 'A: (LSD) A(m) xc 4 % . cI ( g,
Na4C13
NaC12
NaCl (A2d)
NqCl (B2d)
NaCl(3d)
NaCl (C2d)
N%C1 (D2d)
3 -1270.87 -1277.55 967.09 2.10 1.90 5.45 5.32 0.94 0.68 21.80 21.55 18.46 11.57 0.008 0.006 1.76
2 -854.89 -859.45 468.22 1.73 1.52 5.77 5.40 1.22 0.98 15.24 14.90 13.01 12.25 0.03 0.05 7.28
1 -438.59 -441.27 142.70 1.46 1.37 5.17 5.05 0.55 0.52 8.36 8.18 7.25 6.85 0.08 0.17 18.57
1 -438.59 -441.21 158.05 1.46 1.32 5.17 4.99 0.55 0.47 8.36 8.12 7.25 6.79 0.08 0.11 3.22
1 -438.56 -441.16 161.02 1.43 1.26 5.13 4.93 0.52 0.41 8.33 8.07 7.22 6.74 0.05 0.05 0.25
1 -437.62 -440.34 109.60 0.48 0.45 4.19 4.12 0.43 0.41 7.38 7.25 6.27 5.92
1 -437.69 -440.40 122.79 0.55 0.50 4.26 4.17 0.36 0.35 7.46 7.30 6.34 5.98
For definitions of the energies see Tables 1 and 2. 'A: is the total reorganization energy of the cluster (see text), and AE: is the part of the reorganization energy associated with only the interionic Coulomb interactions. For NQCl the energies corresponding to four two-dimensional isomers (A2d, B2d, C2d, and D2d) and for one three-dimensional isomer (3d) are given. LSD denotes results obtained using local spin density functional calculations, and xcg denotes results obtained using LSD with exchange-correlation gradient corrections. Energies in eV.
TABLE 5: Vertical (VIP)and Adiabatic (aIP) Ionization Potentials and Reorganization Energies ER = VIP - aIP, for NqCl,, 1 5 m 5 4, Clusters" VIP N&l3 N&12 NaCl A2d B2d C2d D2d 3d
aIP
ER
LSD 4.78 4.91
xcg 4.85 4.85
LSD 3.89 4.48
xcg 3.95 4.46
LSD 0.89 0.43
xcg 0.90 0.39
4.66 4.36 5.26 5.25 4.09
4.80 4.49 5.39 5.39 4.21
3.92 3.92 5.17 c3.6 3.92
4.10 4.05 5.31
0.74 0.44 0.09
0.70 0.44 0.08
4.05
0.17
0.16
'LSD denotes results obtained using local spin density functional calculations, and xcg denotes results obtained using LSD with exchangecorrelation gradient corrections. Energies in eV.
-s v
-5
-----
--I
--
--_-
c
- - - - -- -
1
Figure 9. LSD Kohn-Sham energy levels for NaCl,, 1 5 m 5 4, clusters. For the NaC1 cluster, energy levels for four 2D isomers (with A2d and B2d being the optimal ones) and for a 3D isomer are shown. Also included are the KS levels for the C1- ion (shifted so that the s-state eigenvalue coincides with that in N%C4). For each cluster levels for both spin orientations are shown, with up-spin on the left and downspin on the right (when these are degenerate a long horizontal line is drawn). The shortest lines correspond to unoccupied levels and the longer ones to occupied levels. potentials were predicted by us p r e v i o u ~ l y in ~ , our ~ ~ ~investigations of Na,-,(NaF), metallization sequences (for n = 4 and n = 14) and were correlated with characteristic odd-even oscillations of energetic properties of the metallic component,"' Nan-,, of the nonstoichiometric halogen-deficient clusters (in this context we provide in Table 6 ionization potentials and
energetics for Nan, 1 I n I 4, clusters). While such a correlation is supported by our results for larger cluster^^^^-'^ (see section 3.4), we remark that for the Na4C1, (0 Im 5 4 ) short metallization sequence the interpretation is complicated due to the change in the dimensionality of the optimal cluster geometries in going from NaC12 to NQC1 (see below). W e also note from Table 4 that the total reorganization energy of the clusters, associated with formation of a NaC1, cluster by removal of a chlorine atom from Na&l,+,, is rather small; A?' = Z$(Na&l,) - ET(Na&l,), where and ET are the total energies of the unrelaxed (that is NaC1, in the geometry of the parent NaCl,+I cluster) and relaxed cluster, respectively. On the other hand the interionic Coulomb part of the cluster reorganization energy AC: (defined in analogy with A$'", but using only the Coulomb energies of the clusters) is rather large. Effective cancellation between AF: and the change in the electronic energy accompanying the relaxation of the structure results in the small calculated values for the total reorganization energies, A r ) . In this context we remark that the reorganization energies associated with ionization, Le., E R = VIP - aIP (where VIP and aIP are the vertical and adiabatic ionization potentials, respectively; see Table 5 ) , are quite significant, reflecting structural variations between the neutral and charged states of the clusters. The electronic structure and bonding in bulk alkali halide crystals are commonly classified as ionic in nature.12 In sodium chloride this results in cubic (rock salt structure) crystals with a 3p-like occupied valence band (VB) separated by a wide band gap from an s-like unoccupied conduction band. Experimentally the band gap in sodium chloride was determined to be Eg = 8.6 eV,Il2 and values for the V B width vary from 1.7 to 4.5 e V (the large degree of uncertainty is attributed to numerous reasons such as instrumental broadening, thermal broadening, surface charging effects, etc.' 13). The VB electronic binding energy, determined from ultraviolet photoelectron spectrometry,Il4 is 10.5 eV. Inspection of the Kohn-Sham energy levels (see Figure 9) for the metallization sequence NWC1, (1 5 m I4)shows that for the stoichiometric limit, (NaC1)4, a p-like valence band of width w = 1.03 e V developed with Eg = 5.31 eV. The VB
E",
Studies of Excess Electrons in Clusters
J. Phys. Chem., Vol. 99, No. 19, 1995 7741 TABLE 6: Energetics of Optimal Ground-State Nan, 1 In 5 4, ClustersD
2
N&
Naq
Na2
Na
-22.90 -23.44 24.01 0.88 0.73 4.35 4.37 4.33 4.36 0.02 0.01
- 16.82
- 1 1.22 -1 1.54 4.86 0.8 1 0.72 5.18 5.25 5.06 5.15 0.12 0.10
-5.205 -5.41
~
ET (LSD) ET (xcg) E, E")(Na) (LSD) E'")(Na)(xcg) VIP(LSD) VIP (xcg) aIP (LSD) aIP (xcg) ER (LSD) ER (xcg)
- 17.30 12.56 0.39 0.35 4.24 4.39 4.07 4.2 1 0.17 0.18
5.21 5.41
a ET, E,, and E")(Na) are the total, Coulomb interionic repulsion, and sodium removal @)(Nan) = ET(Nan-I) E(Na) - &(Nan), energies, respectively. VIP and aIP are the vertical and adiabatic ionization potentials, and E R = VIP - aIP is the reorganization energy upon ionization. Energies in eV.
+
mz
651,
(nr
CI,
I)
Figure 10. (A) Optimal structures for NaCI,,, (1 5 m 5 4) clusters, with LSD optimization. The sequence starts with NaC14 at the upper left comer and, going clockwise, ends with the 2D NaCI (A2d) isomer shown at the lower left comer. Large and small spheres correspond to sodium and chlorine atoms, respectively. For the 3D NaCI:, cluster dNaCl = 2.65 A, L(123) = 97.9", L(234) = 83.4", and 4 3 5 ) = 8 1.9". For the 3D NaC12 cluster d12 = 2.73 A, d23 = 2.63 A, L(123) = 79.6", L(234) = 193.3", and 4 3 4 5 ) = 76.2". For the NaCI (A2d) cluster d12= 3.34 A, d23 = 3.41 A, d34 = 3.40 A, d14 = 3.46 A, dls = 2.51 A, d45 = 2.50 A, L(123) = 61.1", L(234) = 119.2", L(341) = 60.1",L(412) = 1 19.6', and 4 4 5 1) = 87.5". (B) Optimal geometries for the 3D isomer of NaCI (at the upper left corner) and, going clockwise, for the B2d, C2d, and D2d isomers.o For the NaCI (B2d) isomer (upper right), d12= 3.20 A, d23 = 3.47 A, d34 = 3.34 A, d45 = 2.60 A, dls = 2.57 A, d,s = 2.69 A, L(354) = 78.3", L(351) = 75.8", L(123) = 57.9", L(231) = 57.0°, and L(213) = 65.1".
narrows progressively upon depletion of the halide component in the MS (w = 0.9,0.62, and 0.18 eV for Na4Cl3, Na4C12, and NaCl (3d), respectively. We note that the 3D isomer NuCl (3d) is not the lowest energy configuration for this composition. Rather, the minimum energy structures for this cluster are predicted to be two-dimensional (although the differences in total energies between the optimal 2D structures and the 3D one are small; see A2d and B2d in Table 4; see also the ionic structures shown in Figures 10 and 11). The VB width for the A2d and B2d isomers of NaCl are 0.43 and 0.67 eV, respectively.
The excess electrons in halide-deficient sodium chloride clusters occupy energy levels which are split from the bottom of the unoccupied conduction band (CB) of the stoichiometric (NaC1)4 cluster and are located in the gap between the VB and the CB. This behavior is reminiscent of that characteristic of color centers in bulk alkali halides.12 The calculated gap in energy separating the top of the occupied p-like valence band from the lowest energy level associated with the excess electrons decreases gradually; i.e., 3.97, 3.77, and 3.76 eV for Na4Cl3, Na4C12, and N&Cl(3d), respectively, and 2.31 and 3.35 eV for the lowest energy isomers A2d and B2d of NQCl (see Figure 9). Further characterization of the excess electrons in the halidedeficient clusters, whose spatial density distributions are shown in Figure 11, can be given in terms of the inverse participation ratio,'I5 defined as
vi
where is the ith occupied KS wave function, such that pi = 1 for a localized state and 1/Q, where 52 is the volume of the system (in our case, the volume of the wave functions' calculational grid), for a completely delocalized state. Comparison of the values of pi calculated for the excess electrons (Le., highest occupied KS levels) given in Table 7 with those corresponding to states in the VB of the excess-metal clusters indicates that the excess electrons are of a more delocalized nature. The angular momentum character of the excess electrons, obtained via a spherical-harmonics decomposition' l6 of the corresponding KS eigenfunctions is so.*p0.07d0.04f0.03 for the single excess electron in Na4Cl3 and s0~82po~05do~04fo~04 for the two spinpaired excess electrons in Na4C12, portraying their s-like character. The structures of the clusters in the MS are shown in Figure 10 (for Na4Cl4, Na4Cl3, and Na4C12 only the lowest energy structures are shown, while for NQCl the lowest energy structures, denoted as A2d and B2d, as well as those of several other isomers, including a three-dimensional one, denoted as NaCl (3d), are shown). In Figure 11 we exhibit the ionic structures with superimposed contours of the excess electron(s) densities. As seen from these figures, the minimum energy structures for Na4Cl3 and Na4C12 are three-dimensional, deriving from the cubic structure of the stoichiometric NaC14 parent cluster. The excess electron@)in these clusters are "substitutional" in nature occupying mainly the regions of the missing halide anions (see
Bamett et al.
7142 J. Phys. Chem., Vol. 99, No. 19, 1995
(A) 4-
0
2 0 s-
-
"
xs-
?-
9-
-5.0 -3.0 -LO LO x (a,)
3.0 5.0 7.0
I
I ,
I
-5.0 -3.0
-io io
3.0 5.0 7.0
x (a,)
1 1
-5.0
0.0
5.0
-5.0
.-.J
0.0
5.0
1 . 0
5.0
10.0
x (a,)
5.0 -3.0
-LO LO 3.0 5.0 7.0 x (a,)
2
5.0
0.0 x (a,)
7
4
A
n
3-9 x 0
r'
,@, -5.0
0.0
5.0
LJ
4.0
0.0
i.0
0
x (a,)
Figure 11. (A) Contours of the excess electron(s) density (highest occupied KS orbitals, separated from the VB orbitals) for the optimal structures of Na4C1, (0 5 m 5 3). The figure at the upper left comer is for NqCl3 and proceeding clockwise additional chlorine atoms are removed successively, reaching the 2D Na4 cluster shown at the bottom left. The sodiums are denoted by open spheres and the chlorines by filled ones. The contours for Na4Cl3 are plotted in the plane containing atoms 1, 2, and 3: those for Na& are plotted in the plane containing atoms 1 and 2 and bisecting
the line connecting the two chlorine atoms. For the 2D NaC1 (A2d) and NQ clusters the contours are in the plane of the atoms. (B) Same as Figure 1 l A for the 3D NaC1 isomer (top left comer: contours in the plane containing atoms 1, 2, and 3) and, proceeding clockwise, for the B2d, C2d, and D2d isomers of the cluster.
Figure 10A). For NQCl the lowest energy structure is twodimensional (A2d or the close-lying isomer B2d see also Table 2), with the single halide atom bridging an edge of the rhombus made by the four alkali atoms. The density of the excess electrons in these 2D structures of NQC1 is delocalized over part of the metal atom component of the cluster (for comparison
we included in Figure 10A contours of the electron density in the NQ cluster). As previously mentioned, two of the 2D isomers of NqC1 (A2d and B2d) are essentially energetically degenerate. To estimate the barrier for transition between these two isomers we have performed constrained optimization of the energy of
Studies of Excess Electrons in Clusters
J. Phys. Chem., Vol. 99, No. 19, 1995 7743
TABLE 7: Average Participation Ratios for the Excess-Electron(s) Wave Function(s) (Highest Occupied KS Levels), p , and for the Other Occupied States (Valence Band), p', for NaCI,,,, m 5 3, Cluster9 Na4C13 NaC12 NaC1 A2d B2d 3d Na4
Q-I
P'
1.4 x 10-3 1.2 10-3
1.0 x 10-2 1.2 x 10-2
4.6 x 10-5 4.6 10-5
1.2 x 1.2 x 1.1 1.0 x
2.4 x 2.4 x lo-* 2.4 x
1.9 1.9 x 4.6 x 1.6 x
5
4
L1
P
10-3 10-3 10-3 10-3
0
10-5 10-5 10-5 10-5
a R-' is the inverse of the volume of the wave functions' calculational grid.
100
-k? 1 ?ij
40
'
i
I
I
I
8
B2d
I
I
A2d
I 1 I
204
Figure 12. Energy barrier for transformation between the A2d (right) and B2d (left) isomers of NaCl, plotted as a function of the distance d between the C1 atom (filled circle) and the opposing sodium, as indicated. The data was obtained via constrained minimization. Energy values are in eV.
the cluster, starting from the A2d isomer and optimizing the energy as a function of the distance d between one of the Na atoms and the C1 atom, while relaxing the positions of all the other ions in the cluster. From these calculations, the barrier height is determined to be -80 meV (see Figure 12), indicating that at high temperature both isomers may coexist. Finally, we note that our calculations show sensitivity of the vertical ionization potentials to the isomeric structure of the cluster, suggesting that their measurement may allow determination of the lowest energy one. The adiabatic ionization potentials on the other hand are similar for the A2d and B2d isomers, since ionization and subsequent relaxation of both clusters results in the Na&1+ cluster whose structure is shown in Figure 13. Usually, one observes in mass spectrometry nonstoichiometric alkali halide clusters; Le., one detects a large abundance of halogen-deficient (MnXn-l)+clusters (where M is an alkali atom and X a halogen) relative to the stoichiometric (M,X,)+ clusters. Commonly, this finding is described in terms of an ionic model, where it is assumed that the ionization of (MX), results in a weakly bound neutral X atom, interacting with the rest of the cluster via a relatively weak monopole-induced dipole interaction potential.Io3 Consequently, dissociation of the X atom from the charged cluster results in the formation of (MnXn-~)+, Le., (MnXn)+ M,X,-lf X. While to our knowledge no direct measurements of the energetics of the above dissociation reactions are available, several calculations based on ionic models yielded small binding energies of a C1 atom to Na,Cl,-I+ clusters (for n I3, see ref 117; for the case of n = 2, a HartreeFock c a l c ~ l a t i o nyielded ~~ a binding energy of 0.07 eV). In contrast to these results, our calculations yield larger binding energies (see left column of Table 8). Our calculated values, show a strong dependence on the state of the parent Na,Cl,+
-
+
2 Figure 13. Optimal C2, structure of the (NaCl)' ion. Empty and filled circles correspond to sodium and chlorine atoms. Contours of the distribution of the two excess electrons, plotted in the plane of the cluster, are superimposed. d12 = 3.3 A, dl3 = 3.13 A, di4 = 2.70 A, d45 = 2.48 A, L(123) = 56.2", L(231) = 61.9', L(134) = 54.6", and L(341) = 70.8'.
-
cluster; thus smaller values are calculated for the energy required in the dissociation reaction Na,Cl,+ Na,Cl,-I+ C1 when the parent ionized cluster (Na$l,+) is taken in the ground-state configuration of the corrCsponding neutral one (Le., resulting from a vertical ionization of Na,Cl,). Having formed Na,Cl,-l+ clusters other modes of decay are possible. Calculated energetics of two such decay modes, Na,Cl,-l+ NaCl Na,-1Cl,-2 and Na,Cl,-l+ Na+ (NaCl),-l are also given in Table 8). These results show that, energetically, the above decay modes are highly competitive (however, the rates for the above reactions would depend, of course, on the reaction pathways and barriers). 3.4. Metallization of the Na14Cl13Cluster: Na14CIm,8 I m 5 13. Motivated by previous s t u d i e ~ ~of > *Na14F14-, ,~ (1 I n 5 14) we minimize the total energy of Na14C114-, (1 5 n I 6) formed by successive removal of neighboring neutral C1 atoms from one face of the cluster (see Figure 14 for a schematic description of the C1 removal sequence), starting from Nal4Cl13 whose ionic counterpart is a stable "magic" 3 x 3 x 3 filled cuboid with ions of alternating charges located on the edges (other geometrical halide atom removal sequences yield clusters of higher energy). The energetics of the metallization sequence is reflected in the calculated energy of removal of successive C1 atoms (Ef)and ionization potentials (IP) (both vertical and adiabatic), shown in Figure 15 (see also Tables 9-11). We find that the single excess electron in Nal4C113 occupies a diffuse surface state, characterized by relatively low vertical (2.86 eV) and adiabatic (2.83 eV) ionization potentials and, consequently, a very low reorganization energy (0.03 eV) associated with the formation of the Nal&l13+ cluster from the neutral one. In the subsequent halogen-deficient clusters the excess electrons are localized substitutionally in the vicinity of the halogen v a c a n ~ i e s . ~Pertaining .~~~ to the metallization process, it is interesting to note that both the Ef and the IPS show clear oscillations as a function of the number n of excess electrons (Le., Na,4C114-, may be regarded as Na,(NaC1)14-,), resembling the shell-closing effects of free Nan clusters (see ref 111; see also Table 6). At the same time, the difference AIP = IP(Na,) - IP(Nal4C114-,) rapidly decreases, having a local minimum for the Nal4Cllo (Le., Na(NaC1)lo) cluster, which contains an even number (four) of excess electrons. We note that incorporating gradient corrections for exchange69 and
-
+
+
-
+
7744 J. Phys. Chem., Vol. 99, No. 19, 1995
Barnett et al.
TABLE 8: Energetics of Reactions of Ionized Sodium Chloride Clusters (Difference between the Total Energies of the Products and the Parent Clusters), Calculated Using LSD with xcg Correction# Na,CI,+ n
1
2 3 4
-
Na,CI,-I+
LSD 0.29 (0.07) 1.59 (0.42) 1.446(0.96)c 1 (0.97)e 1.34 (0.87)
+ CI
Na,C1,-1+
xcg
n
0.33 (0.12) 1.36 (0.40) 1.096(0.88)c 1.17d (0.96)e 1.06 (0.75)
2 3 4
-NaCl + Nafl-1C1,-2
Na,CI,-I+
LSD 2.03 1.56r 1.5@
2.05 1.59 1.488
2 3
2.43h 2.44'
2.22h 2.30'
4
xcg
-
Na+
+ (NaCI),-I
LSD 2.03
n
xcg
2.05
1.5W
1.58f
1.5W
1.5W
2.29h 2.25'
2.27h 2.14'
a For the dissociation reactions of Na,Cl,+, 1 In I4, results are given for the relaxed parent cluster and (in parentheses) for the unrelaxed one (Le., in the geometry of the corresponding neutral) obtained via vertical ionization of the corresponding Na,Cl, cluster. Energies in eV. Parent Na3C13' in L structure (see Figure 3) and product Na3C12+ in linear configuration (see Figure 7). Parent Na3C13+ in L structure of Na3CI3 (see Figure 1) and product Na3C12+ in linear configuration (Figure 7). dParent Na3C13+ in L structure (see Figure 3) and product Na3C12+ in 2D C2" structure (see Figure 7). e Parent Na3Cl3+ cluster in L structure of Na3C13 (see Figure 1) and product Na3C12+ in 2D C2" structure (see Figure 7). 'Parent Na3C12+ cluster in linear configuration (see Figure 7). 8 Parent Na3C12+ cluster in 2D C2" structure (see Figure 7). Product Na3C12+ in linear configuration (see Figure 7). Product Na3C12+ cluster in 2D C2" structure (see Figure 7).
'
4n
% 3-
v
%2Figure 14. The ground-state structure of Nal4Cl13and designation of the optimal metallization sequence (numbered CI atoms correspond to the optimal chlorine removal sequence). Starred numbered atoms indicate alternative removal sequences, resulting in higher energy isomers (see Tables 9 and 11). In the optimal sequence one successively removes the atoms marked 1-6, that is, starting from Nal4C113 removing the C1 atom marked 1 first yields the lowest energy Na14C112 cluster, continuing by removing the CI atom marked 2 from Nal4C112 yields the lowest energy Nal4Clll cluster, and so on. If instead one removes the internal C1 (marked 1") from Nal4Cl13a higher energy isomer of Na14C112is obtained (see Table 9). Starting from the optimal Nal4C112 cluster (Le., atom 1 removed) and removing next the CI atom marked 2*, instead of that marked 2, yields a higher energy isomer of Nal4C111. Similarly, starting from the optimal Nal4C19 cluster (Le., CI atoms marked 1-5 removed) and removing next the internal atom marked 6*, instead of the edge atom 6 as in the optimal sequence, leads to a high-energy isomer of NalsCls.
c ~ r r e l a t i o nenergies ~~ in a post-LSD approximation (see section 2) does not change any of the trends shown in Figure 15 but yields Efand IP values which are 0.1-0.3 eV smaller. Next we turn to a more detailed characterization of the electronic structure throughout the metallization sequence (see Figures 16 and 17). Examination of the Kohn-Sham oneelectron energy eigenvalues (Figure 16) reveals two main features. First, each cluster shows two perturbed bands, easily interpreted as being formed by s and p electrons of the C1ions (see the level spacing of s and p states of a single C1- ion in Figure 16). These states contribute to the strong ionic bonding characteristic of the clusters. Second, the excess electron states are split off from the bottom of the unfilled conduction band and are well separated from the valence band (the energy gap between the highest occupied valence level and the lowest occupied excess electron level varies from 5.5 eV in Na14C113to 3.3 eV in Nal4C18, Le., n = 1 and n = 6 in Figure 16), in close analogy with the F-center states in ionic solids.l2 The separation of the excess electron states from other states of the cluster (Le., weak mixing) justifies analysis of the spatial
1-
ol
W I
I
I
"
'
1 2 3 4 5 6 n
I
Figure 15. (a) The formation energy Er(Nal4C114-,) = E(Nal4C114-,)
+
-
E(CI) - E(Na14C114-,+1) for successive removal of C1 atoms along the metallization sequence Na14C113 Nal4Cls. (b) Ionization potentials (IP) for Nal4C114-, (1 5 n 5 6), upper curves, and AIP = IP(Na,) IP(Nal4Cl14-,), lower curves. Circles and triangles refer to vertical and adiabatic IPS, respectively.
properties of the n excess electrons by considering only those Kohn-Sham wave functions associated with the n highest energy eigenvalues. To quantify the nature of the spatial distribution of electrons in our systems we use the inverse participation ratio (see eq 27). For each cluster in the metallization sequence, we calculated the average of the pi's corresponding to the n excess electron states as well as the average participation ratio for all the other occupied states. The result is shown in Figure 17 (for Nal4C114-,, 1 I n 5 6) together with the average participation ratios for electrons in free sodium clusters Nan, 1 I n 5 6. Again, we observe remarkably similar behavior between the excess electron states in the halogen-deficient clusters and the electronic states of the corresponding free metal clusters, the only exception being Nal4C113 (i.e., n = 1 in Figure 17), where the single excess electron occupies a rather delocalized surface state as discussed above, covering a much larger spatial volume than the 3s' electron of a free Na atom. We also note that the excess electrons in the Na&l14-, clusters are of a more delocalized nature (i.e., smaller participation ratio) than that of the other (valence) electrons in the cluster. Removing four neighboring Cl atoms leaves one face of the initial Na14C113cuboid (rock salt) cluster covered solely by Na atoms (Figure 18a), portraying face (surface) segregation of a
Studies of Excess Electrons in Clusters TABLE 9: Energetics of Na14Cb, 8 5 m Na14Cl1 3 -5491.44 -5519.20 15602.54 5.54 5.23 96.18 94.55 81.68 77.32 0.03 -0.03 -86.58
ET(LSD) ET (xcg) Ec A(C1) (LSD) A(C1) (xcg) A?) (LSD) A?) (xcg) A ): (LSD) A:) (xcg) AE’ (LSD) AE’ (xcg) A(m) R,c
J. Phys. Chem., Vol. 99, No. 19, 1995 7745 5
13, Clusters, in Their Optimal Structures?
Na14Ch
Na14ChI -4659.11 -4682.91 11771.58 5.29 5.11 82.68 81.16 70.41 66.59 0.30 0.23 - 120.49
-5075.37 -5 I 0 1.20 13517.32 5.73 5.51 89.52 88.00 76.14 72.10 0.28 0.16 - 174.57
Na14Cl1o -4243.29 -4265.02 10131.19 5.76 5.54 76.27 74.73 65.12 61.48 0.12 0.01 -78.09
Na14C19
Na14Cla
-3827.OO -3846.70 8630.87 5.51 5.39 69.40 67.86 59.36 55.94 0.29 0.26 -41.02
-34 10.97 -3428.53 721 1.93
62.78 61.15 53.86 50.55 0.21 0.22 -48.52
For definitions of all entries see Tables 1 and 2. AE’ is the total reorganization energy of the cluster (see text), and Ak”(c‘is the part of the reorganization energy associated with only the interionic Coulomb interactions. Energies in eV. a
TABLE 10: Energetics of Higher Energy Isomers of Na14Cl12, Na14Cl11,and Nal4C18 Clusters, Obtained via Nonoptimal C1 Removal Sequences (See Figure 14)” ET(LSD) ET
(xcg)
E, A(C1) (LSD) W l ) (xcg) AT) (LSD) AT) (xcg) A?’ (LSD) A?’ (xcg) AE’ (LSD) A(m) xc
A
Na14C112*
N~I~CIII*
Na14C18*
-5075.04 -5101.02 12736.41 5.40 5.33 89.19 87.82 75.81 71.92 0.62 0.35 -389.44
-4658.95 -4682.71 11712.72 5.71 4.91 82.51 80.96 70.25 66.39 0.37 0.25 -107.33
-34 10.60 -3428.24 6823.89
62.41 60.86 53.49 50.26 0.35 0.35 -12.62
km)( g, ARS All entries are as in Table 9 (see captions to Tables 1 and 2). Energies in eV.
01
I
I
1
w -21 I
1
2
3
4
5
6
c
n Figure 16. Local spin density functional Kohn-Sham one-electron energy eigenvalues for the C1- ion (shifted so that the s-state eigenvalue equals the lowest eigenvalue in Nal4C113) and for the clusters Nal4ClI4-, (1 5 n 5 6) in the metallization sequence. For each cluster, levels for both spin orientations are shown, with up-spin on the left and downspin on the right (when these are degenerate a long horizontal line is drawn). The shortest lines correspond to unoccupied levels and the longer ones to occupied levels.
TABLE 11: Vertical and Adiabatic Ionization Potentials (VIPand aIP) and Cluster Reorganization Energies CR = VIP - aIP, Calculated Using LSD, and with xcg Corrections, for Na14Cl,, 8 5 m 5 13, Clusters‘ VIP
aIP
ER
LSD
xcg
LSD
xcg
LSD
xcg
2.86 3.68 3.97 3.57 4.07 4.06 3.58 3.93 3.80
2.62 3.49 3.68 3.56 4.07 3.98 3.58 3.85 3.72
2.83 3.40 3.35 3.27 3.70 3.93 3.29 3.66 3.45
2.65 3.33 3.33 3.34 3.81 3.97 3.32 3.63 3.37
0.03 0.28 0.62 0.30 0.37 0.13 0.29 0.27 0.35
-0.03 0.16 0.35 0.22 0.26 0.01 0.26 0.22 0.35
Unstarred clusters correspond to the most stable isomers (see the optimal chlorine removal sequence shown in Figure 14). Starred clusters correspond to higher energy isomers (see Figure 14). Energies in eV.
metallic o~erlayer.~.~-” It is intriguing to note that the sodium overlayer of the ionized Na14C19+(with four excess electrons, Le., Nas+(NaCl)g) exhibits the same geometry as the free Na5+ cluster (a planar structure made of two isosceles triangles sharing a common atom). In fact, we find that starting from Na14C19+ in the geometry of the neutral cluster, the ratio of the length of the edges of the isosceles triangle changes during structural optimization toward the ratio which we find for the free Nas+ cluster. For Nas+ we find the ratio of the shorter edge to the longer one to be 3.1Y3.49 (lengths in A), and in the optimal structure of Nal4C19+ we find it to be 3.42/3.58; in neutral Nal4Cl9 the corresponding ratio is 3.6U3.34. Figure 18b also shows
o!
? 1
I 1
2
3
4
5
6
n Figure 17. The inverse participation ratio p along the metallization sequence Na14C114-,, 1 5 n 5 6 (triangles and circles) and for free Nan clusters (asterisks). Circles correspond to the participation ratios averaged over the Kohn-Sham orbitals of the excess electrons, and triangles denote the average participation ratios of the valence electron orbitals.
that the spatial distribution of the excess electron density in the overlayer is almost indistinguishable from the electron density in the free Na5+ cluster. The stability and dynamics of the metallic overlayer were studied in a relatively long molecular dynamics run (-7.4 ps) starting from the optimized structure and adding kinetic energy by randomizing the velocities at random intervals during the first 0.5 ps. The vibrational temperature during the subsequent constant-energy period was 660 f 90 K, and the total energy was conserved to better than 0.7 meV/ion ( < 5 x of the total energy). No structural transformation was seen during the run. The dynamics of the stable overlayer is portrayed in Figure 19, which shows the vibrational density of states, calculated by
Barnett et al.
7746 J. Phys. Chem., Vol. 99, No. 19, I995
"9
/ \ A '
i Figure 18. (a) The optimized structure of the neutral Na14C19cluster showing a segregated sodium overlayer. Small and large spheres depict sodium and chlorine ions, respectively. (b) Excess electron density in the plane of the sodium overlayer of Nal4C19+ (top), and electron density of the free Nas+ cluster (bottom). The spacing between contours in b is 0.001 au.
I\/
V
/
c)
3o
0
1
I
100
200
"\
I
300
o (cm-') Figure 19. Vibrational density of states (VDOS) in the Na14C19+cluster calculated separately for the sodium overlayer (solid curve) and for all the other sodium ions (dashed curve) from the velocity autocorrelation function obtained in BO-LDA-MD simulations at 660 K. The peaks at o = 0 correspond to the rotational motion of the cluster.
Fourier transforming the velocity-velocity autocorrelation function, averaged separately for the five Na ions belonging to the overlayer and for the rest of the Na ions in the ionic part of the cluster. The shift to lower frequencies of the vibrational spectrum associated with the sodiums in the metallic overlayer is correlated with a softening of the interatomic forces in that layer, as compared to those in the ionic part of the cluster. Finally we comment on clusters at the other end of the metallization sequence, Le., Na,CI, with n - m > m, for example Nal4C14. Our previous simulations of sodium fluoride c l ~ s t e r s ~indicate . ~ - ~ that for clusters of this size (Le., n = 14) and in this composition range the alkali halide component (i.e., (NaF), in Na,-,(NaF),) tends not to be uniformly surrounded by the excess sodiums. Rather, in the optimal structures of these clusters, only partial wetting of the ionic component by the metallic one is found. As a demonstration of this effect we have performed optimization calculations for Na14C14,starting from three different atomic arrangements which may be obtained from the Nal4C113 cuboid structure shown in Figure 14 via alternative sequences of chlorine atom removal (we have not attempted here to determine the global optimized structure for Na14C14). The local structural minima corresponding to these alternative starting configurations are shown in Figure 20 (isomer A was obtained from a starting configuration consisting of a Na5C14 planar facet with the other nine sodium atoms on one side of it; isomer B corresponds to a starting configuration obtained via removal of C1 atoms except for a corner (NaC1)4 cube; and isomer C was obtained from a starting configuration where C1 atoms were removed above and below a N a C h plane). While the structures shown in Figure 20 do not represent the global minimum for the Nal4C14 cluster, our results for these
Figure 20. Three isomers (local potential energy minima) of the Na14C14 cluster. Isomer 1 is the energetically favorable one. Large spheres
correspond to chlorines and smaller ones to sodiums.
structures indicate that isomer A is energetically the most favorable one, followed by isomer B, and isomer C is least favorable among the three (in this context it is interesting to note that in the short-time finite-temperature MD simulations which we performed we observed a tendency for structure C to transform to configurations with a more segregated metal component (i.e., similar to isomers B and A). 4. Protonated Water Clusters:
H30+
and
(H20)2H+
Protonated water clusters (the hydronium ion, H30+, and the hydrated hydronium, H30+*nH20, or (H2O),H+) are systems of importance in a number of areas of chemical and physical research. Such systems, have been detected in the D region of the ionosphere, as well as being the dominant naturally occurring in flames,16 and ions in the troposphere and in moist air.I7 Additionally, they are produced in chemical ionization mass spectrometers where they act as protonating and solvating reagents. I I Furthermore, investigations of such species provide insights into the nature of intermolecular hydrogen b ~ n d i n g , ~proton ~ - ~ ~transfer p r o c e ~ s e s , ~ atmo~,~~ spheric ion-clustering reactions, and the transition from gasphase to solution-phase chemistry via ion solvation processes. Motivated by the above, several experimental studies and quantum-chemical theoretical investigations of protonated water clusters have been r e ~ o r t e d . ~ ~ - ~ ~ In this study we focus on the H3O+ and H5O2+ ions, using the AQMD method, which, as described in section 2, combines
Studies of Excess Electrons in Clusters
J. Phys. Chem., Vol. 99, No. 19, 1995 7747
TABLE 12: Energetics of Protonated Clusters A,X+, n = 1 and 2, for A = H20, NHJ, and (NHJ)(H~O), and X = H, Li, and Naa H20 H20 H20
H H H Li Li Li Na Na Na H H H H H H Li Li Li Na Na Na H H H H H H
1 LDA 7.361 1 GGC 7.5 19 7.229,h 1 exp 7.19," 7.02d 1 LDA 1.603 1 GGC 1.366 1.47' 1 exP 1 LDA 1.097 1 GGC 0.81 1 1 exP 1.04' I LDA 9.065 I GGC 9.30 1 exP 8.871f 2 LDA 8.90 2 GGC 8.955 2 exP 2 LDA 2.577 2 GGC 2.354 2 exP 2.36' 2 LDA 1.679 2 GGC 1.412 2 exP 1.67' 2 GGC 10.452 2 exP 1 GGC 9.936 1 GGC I exP I exP
J
/
J
/
1.934 1.549 1.370: 1.431,R 1.56Ih 1.366 1.151 1.12' 0.974 0.833 0.86' 1.233 1.075,' 1.17! l . l @ 0.877" 2.646" 0.863'*" 2.37"~"
a The binding energy of X+ to An is AI = E[An] - E[AnX+] and the hydration energy of AX+ is given by A2 = E[A~-IX+]+ E[Al E[AnX+]. In the case of (NHs)(H20)H+ results for A2 are given for (NH3)(HzO)H+ NHa+ H20 and (NHs)(H20)H+ HsO+ NHa. Results are given for LDA calculations and for LDA calculations including exchange-correlation gradient corrections self-consisteatly (GGC), with a plane-wave energy cutoff Ec = 96 Ry. Energies in eV. Reference 12 1. Reference 122. Reference 123. e Reference 124. f Reference 125. * Reference 126. Reference 127. Reference 128. j Reference 129. Reference 130. Reference 131. " Reference 132. 'I (NHs)(H20)H+ NH4+ H20. * (NHs)(H20)H+ NHs H3O+.
-
-
+
+
-
-
+
+
a quantum path-integral treatment of the ionic degrees of freedom (in our calculations the protons are treated quantum mechanically and the oxygen ions classically) with concurrent electronic structure calculations, based on the local density functional method, of the Born-Oppenheimer potential energy surface. Energetics and structural information pertaining to protonated water clusters (H20),H+ (n = 1 and 2), (NH3),H+ (n = 1 and 2), (NH3)(H20)H+,and (H20),M+ (n = 1 and 2) for M = Li and Na, are given in Table 12 and Figures 21-23. As seen from Figure 21 the optimal structure of the hydronium ion H3O+ is of C3" symmetry, with the three protons equidistant from the oxygen, do^ = 0.984 A. The optimal structures of (H20)Li+ and (H20)Na+ are of C2" symmetry with dLi-0 = 1.93 A,and ~ N ~ -=o 2.48 A, and the structure of the ion is of tetrahedral symmetry, with ~ N - H = 1.02 A. For comparison, the calculated distances in H20 and NH3 are &H = 0.964 8, and &H = 1.01 A,respectively, and L(H0H) = 104.7" and L(HNH) = 107.2". For the (H20)2H+ ion the optimal structure shown in Figure 22 is of C2 symmetry (compare to the optimal structure of the water dimer, (H20)2, shown at the top of Figure 23), with the distance between the oxygens to the hydrogen-bonding proton do++ = 1.196 A. For the (H20)2Li+ and (H20)2Na+ ions the optimal structures are of C2" symmetry with do-Li+ = 1.94 A, and = 2.54 A,and in the optimal structure of (NH3)2H+ &-H+ = 1.307 8, (see Figure 22). In the mixed (H2O)(NH3)-
m+
\
ks)
Figure 21. Optimized structures (at the GGC level, with a plane-wave cutoff energy of 96 Ry) for H3O+, (H20)Li+, (H20)Na+, and NH4+ (from top to bottom). The smallest spheres correspond to hydrogens and the largest spheres to an oxygen or a nitrogen; the darker intermediate size sphere corresponds to a metal atom. The geometrical parameters of the molecular ions are, for H30+, OH = 0.984 A L(H0H) = 114.2', thc distance of the oxygen from the plane of hydrogens ho = 0.244 A, and the angle betwten ho and the OH bond 8 = 75.8'; for (H20)Li+, &i+ = 1.930 A, &H = 0.96 A, and L(H0H) = 105.7'; for (H20)Na+, = 2.48 A, do^ = 0.96 A, and L(H0H) = 105.2'; for N&+, ~ N = H 1.019 A, and L(HNH) = 109.47'.
H+ molecular ion (not shown) the proton bond is nearly linear (Le., L(NH+O) = 177.9') with &H+ = 1.56 A,and ~ N H += 1.07 A. The water molecule in the cluster is in the same plane as the oxygen, H+, and nitrogen, with &H = 0.961 8, and L(H0H) = 107.2". For the ammonia molecule in the mixed cluster ~ N = H 1.016 8, and L(HNH) = 108.7" and one of the NH bonds of the ammonia molecule is in the NH+O plane. Included in Figure 23 are isomeric structures of higher energy for the metal-water dimer cations and for (NH&H+ (the energy increments A of these isomeric structures, above the corresponding optimal ones, are A[(H20)2Li+] = 0.054 eV, A[(H20)2Na+] = 0.144 eV, and A[(NH3)2H+] = 0.012 eV). From inspection of the energetics data displayed in Table 12 we observe the following (based on GGC results): (i) The calculated proton affinity of the ammonia molecule is 22.7% larger than that of H20 (in excellent agreement with the experimental value of 23.1%). (ii) The proton binding energy to H20 is significantly larger than those for alkali cations, with the values for the latter decreasing for the heavier metal. (iii) The proton affinity of (H20)2 is 19% larger than for the H20 monomer, and that of (NH3)2 is 12% larger than that of NH3 (see A, in Table 12). (iv) The proton affinity of (NH3)2 is larger than the proton affinity of (H20)2 (by 16.7%), and the value for (NH3)(H20) is intermediate between the two (larger by -1 1% than that of (H20)2). (v) The hydration energy (Le., the binding energy of H20 to H30+) of the hydronium ion is 25% larger than the binding energy of NH3 to NH4+ (see A2 in Table 12). (vi) The hydration energy of H3O+ is 76% larger than that of N&+ (see A2 in Table 12), and the dissociation energy of (NH3)(H*O)H+ into H3O+ NH3 is about 3 times larger than that for dissociation into NH4+ H20. (vii) The binding energies of Li+ and Na+ to (H20)2 are higher than the corresponding binding energies of the ions to H20, exhibiting
+
+
Barnett et al.
7748 J. Phys. Chem., Vol. 99, No. 19, 1995
i
k
Figure 22. Optimized structure for (H20)2H+,(H20)2Li+,(H20)2Na+, and (NH3)2H+(from top to bottom), calculated via LDA-GGC with E, = 96 Ry. The geometrical parameters are, for the C2 (H20)2H+cluster, do^+ = 1.197 A, do^ = 0.964 A, L(OH+O) = 174.5', L(H0H) = 110.4', L(H1OH') = 121.1', and L(H20H+) = 119.7', (calplations using LDA-GGC, with E, = 158 Ry yielded do^+ = 1.196 A, d m = 2.389 A, do^ = 0.973 A, L(OH+O) = 174.5", L(H0H) = 110.9', L(H10H+) = 120.9', and L(H20H+) = 119.3"). Also shown is a view along the 0 1 0 2 axis, the angles 61 and 6 2 between the projections of the OIHI and OlH2 bonds on a plane normal to the 0 1 0 2 axis and the vertical line indicated in the figure are 22.4' and 56.3', respectively (61=;2.8" and 6 2 = 56.3" for Ec = 168 Ry); for (H20)Na+, = 2.54 A, d m = 5.08 A, OH = 0.96 A, L(H0H) = 104.5', and the angle between the normals to the two H20 molecular planes and between the normal to each of the moleculoar planes and d m is 90"; for the (NH3)2H+ molecule, ~ N H += 1.307 A, ~ N = N 2.614 A, ~ N = H 1.014 A, L(HNH) = 107.5', and the two proton-bonded NH3 molecules are rotated by 60" with respect to each other about the NN axis.
a much larger dependence on the size of the water cluster than in the case of proton binding (compare AI values in Table 12; see also ref 53 for a discussion of Na(H20),, 1 5 n 5 6, clusters). We turn next to a discussion of results obtained via AQMD simulations of H30+ and (H20)2H+ at 150 K. In these simulations the electronic structure was calculated within LDA (Le., no exchange-correlation gradient corrections; a plane-wave energy cutoff Ec = 96 Ry was used in the simulations). First we show in Figure 24a for H30+ the distribution of distances of the protons from the oxygen, obtained via the AQMD method (solid line), and compare it with the results obtained via BOLDA-MD simulations (dashed line where the protons were treated classically) at the same temperature. As seen from Figure 24a the OH distances are distributed over a broader range when the protons are treated quantum mechanically. The average value at 150 K for do^ in calculated from the BO-LDA-MD simulations is &H = 1.852~0(a0 = 0.52918 A), and the expectation value from the AQMD simulation is @OH) = 1.866~0compared to a value of 1.844~0calculated for the optimal (zero temperature) structure of H30+ with no xcg corrections.
d Figure 23. Optimized structure of the water dimer (H20)2 (top) and isomeric structures for (H20)2Li+, (H20)2Na+, and (NH&H+ (where the two proton-bonded NHj molecules are in an eclipsed relative orientation), whose corresponding lowest energy structures are shown in Figure 22.
For (H20)2H+, the distance distributions for do^ in the H20 molecules and for the distance &H+ between the oxygens and the "intermolecule" proton, obtained from the BO-LDA-MD (dashed) and AQMD (solid line) simulations, are shown in Figure 24b. Again the distance distributions obtained via the AQMD simulations are broader, and the &H+ distributions are broader than those for do^. Th~averagevalues obtaked via the BO-LDA-MD method are &H = 1.833~0and &H+ = 2.256~0,and the expectation values from the AQMD simulation are (&H) = 1.849uo and @OH+) = 2.281~0,compared to do^ = 1.827~0and do"+ = 2.242~0in the optimal (zero temperature) LDA (with no xcg) structure. Further characterization of the intermolecular proton in (H20)2H+ is given in Figure 25, where we show the distributions (obtained via BO-LDA-MD and AQMD simulations) of the magnitude of the projection of the position vector of the proton on the normal plane bisecting the interoxygen vector (De in Figure 25a), the distribution of the projection on the interoxygen vector (D,in Figure 25b), and the distribution of the magnitude of the distance between the proton and the midpoint of the interoxygen distance. We obseve that the effect of the quantum treatment of the proton is to broaden these distributions and that the width of the distributions reflects a rather shallow potential energy surface for motion of the "bonding" proton. Further analysis of the dihedral angle distribution and that of the angle between the two H20 molecular planes indicates that (H20)2H+ is a rather floppy molecular ion, with low-frequency librations of the H20 molecules. Another interesting finite-
J. Phys. Chem., Vol. 99, No. 19, 1995 7749
Studies of Excess Electrons in Clusters
8.0 1. 6.0 n
I
1:
TaO
Y
4.0
ni 2.0 0.0 12.0
I
I
1
8.0
1
0.0
'
/ I
> .' '
1.6
0.1
0.3
0.5
0.7
0.9
j,
,
. ,, 2.0
2.4
2.8
rOH
-.-
-0.5 -0.3 -0.1 0.1
Figure 24. Normalized distributions of distances between the protons and the oxygen in (H*O)H+ (shown in a) and in (H20)2H+ (shown in b). Results are shown for simulations at 150 K where the protons were treated classically (BO-LDA-MD, denoted by dashed lines) and for AQMD simulations where the protons were treated quantum mechanically (solid lines). In b the distributions about - 1 . 8 ~correspond to the TOH in the water molecules and those about 2 . 2 ~ correspond 0 to the distance between the intermolecular proton and the molecular oxygens. Distance in units of the Bohr radius a0 = 0.52918 A.
temperature feature of the (HzO)zH+ molecular ion revealed by the simulations is exhibited in Figure 25d where the calculated distributions of the difference between the distances of the intermolecular proton from the two oxygens of the water molecules are shown. We note that values of the difference Ad = [do,H+- d o 2 ~ +(where ] 01and 0 2 are the oxygens of the proton-bonded water molecules) obtained in the AQMD simulations (at 150 K) can achieve values which are up to twice those obtained in BO-LDA-MD simulations where the protons are treated classically (dashed line). Moreover we observed some correlation between Ad and the magnitude of the interoxygen distance, indicating that larger values of Ad correspond to larger intermolecular distances (Le., larger interoxygen separations). We expect that at higher temperature one may observe even larger variations in Ad, correlated with larger values of the interoxygen distance, resulting in the occurrence of configurations of ionic character (H30+.H20 and H2OvH30+) in the equilibrium ensemble (in this context we note for a certain interoxygen distance the transfer of a proton from H30+ to H2O involves a potential barrier, which collapses for smaller intermolecular distances). Further characteristics of the quantum nature of the protons in H30+ and (H20)2H+ are provided by
where is the position of the ith bead (pseudoparticle) on the isomorphic necklace (see section 2 ) , and the imaginary time correlation
0.3 0.5
n
Tg U L
n
0.1
0.3
0.5
0.7
0.9
r (ad
Figure 25. (a,b,c) Distributions of the projections of the vector position of the intermolecular proton in (H20)2H+ on the normal plane bisecting the interoxygen vector (De in a) and on the inter-oxygen vector (0:in b). Distributions of the distance between the intermolecular proton and the midpoint of the inter-oxygen vector are shown in c. Shown in d are distributions of the difference between the distances of the intermolecular proton and the two molecular oxygens, Ad = I ~ o , H & do2"+I. Results obtained via AQMD simulations are given by solid lines, and those obtained from BO-LDA-MD simulations are given by dashed lines. Both simulations were performed at 120 K. Distances are given in units of the Bohr radius a0 = 0.52918 A.
The above quantities, which can be calculated for each of the protons, reduce for a free particle to R T = ~ 3'I2A where 1 is the thermal wavelength ( R ~ ~ ( 1 5K,0 H+) = 1 . 8 5 5 ~and ) @(t t') = 3A2(t - t')@i - ( t - t'))/@~)*, which for t - t' = @h/2
7750 J. Phys. Chem., Vol. 99, No. 19, 1995
Bamett et al.
0.8
0.0
1
-0.4
0.0
0.4
U 0.6 Figure 27. Distributions of u = cos 0, where 8 is the angle in HsO+ between an OH bond and the vector normal to the plane of the hydrogens to the oxygen atom. Results obtained from simulations at 0.4 150 K, using the BO-LDA-MD method (classical protons, dashed line) and the AQMD method (quantum protons, solid line), are shown. u = 0.2 0 corresponds to the inversion planar transition state configuration.Note the quantum enhancement of the inversion probability indicated by an increase in the probability to find the system in the vicinity of the 0.0 transition state. 1 2 3 4 5 6 7 8 LDA calculations with no xcg corrections and with a lower k plane-wave energy cutoff Ec = 96 Ry, a lower barrier of 0.033 Figure 26. Imaginary time correlation functions G?, plotted versus k eV was obtained). In the optimal pyramidal geometry of H30+ (where 0 5 k @hP)5 @),I for the protons in H30+ (in a) and for the the calculated (Ec = 96 Ry with GGC) OH distance is doH = molecular hydrogens and intermolecular proton in (H20)2H+ (in b). 1.859~0(1.851~0for Ec = 158 Ry) with an angle L(H0H) = The results were obtained from AQMD $mulations at 150 K, R is in 114.2' (the height of the oxygen above the three hydrogen plane units of the Bohr radius a0 = 0.52918 A. is ho = 0.461m (0.452~0for Ec = 158 Ry) and the angle between the normal to the hydrogen plane to the oxygen and (diameter of the necklace) takes the value .@&/3W2)= (3/4)A2an OH bond is 8 = 75.85') (the L(H0H) and 6 angles are the ( a ( p W 2 ) ; 150 K, H+) = 0.927~0). same for both Ec = 96 and 158 Ry calculations). In the The value for RT (averaged over the three protons) for H30+ transition state for inversion the geometry of H30+ is planar obtained in our simulations is 1.588~~0; RT for the molecular with doH = 1.849~0 (1.837~0 for E, = 158 Ry) and L(H0H) = hydrogens in (H20)2Hf was obtained as 1.523~0, while the value 120". (In calculations with Ec = 96 Ry and no xcg corrections for the intermolecular proton is 1.243~0, all of which are smaller doH = 1.849~0,and ho = 0.417~~0, L(H0H) = 115' and 8 = than R T at ~ 150 K. 76.95' for the optimal geometry, and do" = 1.849~0for the The imaginary time correlation functions G?(t - t') for planar transition state configurations). These values are in and (H20)2H+ are shown in Figure 26a and b, respectively (in agreement with those used in the analysis of spectroscopical Figure 26b we display @,for both the molecular hydrogens and data for H30+.28*29 for the intermolecular proton). First we observe that for both While a determination of the inversion rate using the paththe protonated clusters GZ @W2) is smaller than its value for a integral-based AQMD simulations is feasible using the centroid free proton at 150 K. We also note that the character of the method,76css0we show in Figure 27 results for the distribution variation of G? with k (where k is an index of the pseudoparticle of u = cos 8 (8 is the angle between the normal to the hydrogen on the classical isomorphic closed necklace, Le., 0 5 t - t' = plane to the oxygen and an OH bond vector) obtained at 150 K k(pWP) 5 tip) is different for H30+ and for (H20)2H+ (see in via BO-LDA-MD simulations (with classical protons, dashed particular Figure 26b for H+). Thus while for &O+ 5% varies line) and AQMD simulations (with the protons treated quantum in a manner similar to that predicted for a free particle (see mechanically via the PI method, solid line). As noted before expression for G7dt - t') following eq 29), those for (H20)2Hf in these simulations, the electronic Born-Oppenheimer potential (and in particular the one corresponding to the intermolecular surfaces were calculated using LDA, with no xcg correction, proton) exhibit a region in t - t' where B saturates. This characteristic behavior is termed ground-state d ~ m i n a n c e , ~ ~ ~ -with ~ * ~Ec = 96 Ry; the corresponding calculated inversion barrier is 0.033 eV, or 373 K, which is over twice the temperature of occurring when PAE >> 1, where AE is the energy gap between the simulations (150 K). A value of u = 0, in Figure 27, the ground state and the first excited state of the particle. From corresponds to the inversion saddle (transition state) planar the above we conclude that while for the protons in H30+ and configuration (see a snapshot, taken from the AQMD simulation, (to a somewhat lesser extent) the molecular hydrogens in shown in Figure 28), and other values of u (positive and negative (H20)2H+, low lying excitations of the quantized nuclear on the two sides of the saddle) correspond to pyramidal reactions are populated at 150 K, the ground-state dominates configurations, see e.g., top of Figure 28, (Le., a change of sign for the intermolecular proton in (H20)2H+. of u corresponds to an inversion). Comparison between the The structure and the nature of the inversion tunneling in results obtained via the two simulations provides a measure of the hydronium ion, H@+, have been the subject of a number the quantum enhancement of the inversion process. Further of experimental and theoretical s t ~ d i e s . ~The ~ - structure ~ * ~ ~ ~of simulations (including GGC) to determine the inversion rate H30+ has a pyramidal geometry, similar to its isoelectronic are in progress. counterpart NH3. Generalized gradient correction (GGC) LDA calculations with a plane-wave cutoff E, = 158 Ry (a 543 grid, 5. Summary with a grid spacing of 0.25~0) yield an inversion barrier of 0.064 eV (5 16 cm-I) compared to a barrier of 672 cm-I obtained via This paper has two aims, a methodological one and a physical analysis of high-resolution spectroscopic measurement^^^,^^ (in one. Pertaining to methodology, we have discussed and
Studies of Excess Electrons in Clusters
Figure 28. Snap shots of configuration taken from AQMD simulations of HaO+ at 150 K. The large sphere corresponds to the oxygen atom and the small spheres represent the pseudoparticles of the classical isomorphism (P = 8 pseudoparticles representing each proton). The configuration at the top corresponds to a pyramidal structure of HxO+, and the ones at the middle (side view) and bottom (top view) of the figure show a representative transition state structure.
employed electronic structure calculations based on the local density functional (LDA) and local spin density functional (LSD) theory, with the use of nonlocal norm-conserving pseudopotential and plane-wave basis sets, for evaluation and optimization of the structures of sodium chloride and protonated water (and ammonia) clusters and investigations of sizeevolutionary patterns of their energetics and properties (energy levels, binding energies, ionization potentials, and the energetics of certain reactions). These calculations used a newly developed method for explorations of optimal structures and investigations of the dynamics of finite systems, via ab-initio molecular dynamics (MD) simulations of the evolution of the ionic degrees of freedom on concurrently calculated electronic Born-Oppenheimer (BO) potential energy surfaces (the BO-LDA-MD method, see ref 39, and section 2). Furthermore, we have presented a novel method40 unifying the quantum path-integral description of the ionic degrees of freedom with the BO-LDAMD methodology, thus allowing investigations of physical systems where both the electronic and ionic degrees of freedom are treated quantum mechanically (the all-quantum MD simulation method, AQMD, described in section 2; see also ref 40). In section 3 we presented an extensive study of stoichiometric and nonstoichiometric sodium chloride clusters [(NaCI),, 1 5 n 5 4, NaCI,, 0 5 m 5 3, Na3CI,, 0 5 m 5 2, Na2Cl,, m = 1 and 2, and Na&114-, for 1 5 n 5 6 and n = lo]. Our findings include (i) The optimized structures of (NaCI), clusters for n 2 4 are three-dimensional (3D) cuboids (rock salt), while for n = 2 and 3 they are two-dimensional (2D). (ii) In ionized (NaCl),+, 2 5 n 5 4 clusters, we observed a tendency for contraction of an inter-chlorine distance, which in (NaC1)2+ reduces to a value close to that calculated for the C12- molecule,
J. Phys. Chem., Vol. 99, No. 19, 1995 7751 suggesting that the planar D2h rhombus (NaC1)2+ molecule may be described as Naf(C12-)Na+ (see section 3.1). This pattern is reminiscent of the structural relaxation associated with selftrapped excitons in bulk alkali halide crystals, involving (iii) Of the halideformation of a C12- molecular ion." deficient clusters the optimal structures are 3D for NaCl3 and NaC12, with the ionic arrangements similar to that in the parent cubic NaC14 cluster, and the excess electron(s) distribution is localized in the regions of the missing halide atom(s). For N a C1 the lowest energy isomeric structures are 2D, with the metal ions forming an approximate rhombus (resembling the structure of N a ) and the single chlorine ion capping one of the edges of the rhombus (the two lowest isomers are practically degenerate in energy, and the transformation between them involves a barrier of -80 meV). The electronic distribution of the excess electrons in these 2D clusters is of a delocalized character. Analysis of the electronic distributions and inverse participation ratios indicates that the excess electrons are of a more delocalized character than those involved in ionic bonding (the valenceband electrons). (iv) In analogy with F-centers in bulk alkali halide crystals and their defect aggregates (R and M centers12), the energy levels of the excess electrons in nonstoichiometric alkali halide clusters are split from the bottom of the unoccupied manifold (conduction band) lying in the gap between the occupied (valence-band) and unoccupied levels. (v) Our investigations of nonstoichiometric clusters, particularly for larger clusters, such as the metallization sequence Nal4CI,, 0 5 m 5 14, suggest a description of nonstoichiometric alkali halide clusters Na,CI, (m < n) as composed of phase segregated metallic and ionic components, Le., Na,-,(NaCI),, which in sufficiently large clusters leads to structures exhibiting a face (or surface) segregated metallic layer (such as in Nas(NaC1)9 (see Figure 18) Nas(NaF)s, which are composed of a segregated Nas layer adsorbed on a face of a 3 x 3 x 2 ionic crystallite). The stability and the metallic character of the Nas(Nas+) overlayer in Nal4C19 (Na,4C19+) was demonstrated via BO-LDA-MD simulation at 660 K. This picture is also supported by odd-even oscillations (with n - m) in defect formation energies (i.e., the energies for removal of successive halide atoms along the metallization sequence) and in the ionization potentials. These results suggest further structural and electronic studies (both experimental and theoretical) of larger clusters and of alkali metal overlayer(s) on alkali halide substrates. 133 In section 4 we discussed optimal structures and energetics, obtained via the BO-LDA-MD method with self-consistent exchange-correlation gradient corrections (GGC), for protonated AH20),H+ and (NH3),H+ for n = 1 and 2, as well as for (H20)(NH3)H+, and (H20),M+, for n = 1 and 2 with M = Li and Na. Furthermore, we presented the first results from our newly developed all-quantum molecular dynamics (AQMD) for H30+ and (H20)2H+ clusters, simulated at 150 K. Our investigations demonstrate broadening of the distributions of internuclear distances involving the hydrogens due to their quantum nature. In the proton-bonded (H20)2H+, [(H20)H+(HzO)], clusters, the state of the intermolecular proton is characterized by ground-state dominance (see Figure 26) and our results for the spatial distribution of the proton suggest that at higher temperatures configurations of ionic character (i.e., H30+*H20 and H2O*H3O+) may develop. Finally, quantum enhancement of the inversion isomerization in H3O+, at 150 K, was demonstrated. Acknowledgment. This research is supported by the U.S. Department of Energy (Grant No. DE-FGO5-86ER-45234). H.H. acknowledges partial support from the Academy of Finland.
7752 J. Phys. Chem., Vol. 99, No. 19, 1995
Calculations were performed on CRAY computers at the Florida State University Computing Center, the National Energy Research Supercomputer Center, Livermore, California, and the GIT Center for Computational Materials Science. Useful conversations with R. L. Whetten and the assistance of C. L. Cleveland in generating the figures containing molecular configurations are gratefully acknowledged. References and Notes (1) See articles in Physics and Chemistry of Finite Systems: From Clusters to Crystals; Jena, P., Khanna, S. N., Rao, B. K., Eds.; Kluwer:
Dordrecht, 1992. (2) See articles in Clusters of Atoms and Molecules; Haberland, H., Ed.; Springer: Berlin, 1994. (3) See Landman, U.; Bamett, R. N.; Cleveland, C. L.; Rajagopal, G. in ref 1, Vol. 1, p 165. (4) Jortner, J. Z. Phys. D 1992, 24, 247. (5) Steigerwald, M. L.; Brus, L. E. Annu. Rev. Mater. Sci. 1989, 19, 471. Siegel, R. W. Mater. Sci. Eng. B 1993, 19, 37. (6) Landman, U.; Scharf, D.; Jortner, J. Phys. Rev. Lett. 1985,54, 1860. (7) Bamett, R. N.; Landman, U.; Scharf, D.; Jortner, J. Acc. Chem. Res. 1988, 22, 350. (8) Rajagopal, G.; Bamett, R. N.; Landman, U. Phys. Rev. Lett. 1991, 67, 727. See also ref 3. (9) Rajagopal, G. Ph.D. Thesis, Georgia Institute of Technology, Atlanta, 1992. (10) Bamett, R. N.; Landman, U.; Rajagopal, G.; Nitzan, A. Isr. J. Chem. 1990, 30, 85. (1 1) (a) Hakkinen, H.; Bamett, R. N.; Landman, U. Europhys. Lett. 1994, 28, 263. (b) Hakkinen, H.; Bamett, R. N.; Landman, U. Chem. Phys. Lett. 1995, 232, 79. (12) (a) Physics of Color Centers; Fowler, W. B., Ed.; Academic: New York, 1968. (b) Markham, J. J. F-Centers in Alkali Halides; Academic: New York, 1966; Solid State Physics, Supplement 8. (c) Mott, N. F.; Gurney, R. W. Electronic Processes in Ionic Crystals; Dover: New York, 1964. (d) Stoneham, A. M. Theory of Defects in Solids: Electronic Structure of Defects in Insulators and Semiconductors; Clarendon: Oxford, 1975. (13) Arijs, E.; Nevejans, D.; Ingels, J.; Frederick, P. J. Geophys. Res. 1985, 90, 5891. (14) Viggiano, A. A.; Dale, F.; Paulson, J. F. J. Chem. Phys. 1988, 88, 2469. (15) Ferguson, E. E. NATO Adv. Study Inst. Ser., Ser. B 1979, VB40, 377. (16) Knewstubb, P. F.; Sugden, T. M. Proc. R. Soc. London 1960, A255, 520. (17) Kebarle, P.; Godbole, E. W. J. Chem. Phys. 1963, 39, 1131. (18) Narcisi, R. S . ; Barley, A. D. J. Geophys. Res. 1965, 70, 3687. (19) Fehsenfeld, F. C.; Dotan, 1.; Albritton, D. L.; Howard, C. J.; Ferguson, E. E. J. Geophys. Res. 1978, 90, 1333. (20) Amold, F.; Krankowsky, 0.;Marien, K. M. Nature 1977,267, 30. (21) Heurtas, M. L.; Fontan, J. Atmos. Environ. 1982, 16, 2521. (22) Frisch, M. J.; Del Bene, J. E.; Binkley, J. S.; Schaefer, H. F., 111. J. Chem. Phys. 1986, 84, 2279. (23) Lee, E. P. F.; Dyke, J. M. Mol. Phys. 1991, 73, 375 and references to earlier studies cited therein. (24) Xie, Y.; Remington, R. B.; Schaefer, H. F., 111. J. Chem. Phys. 1994, 101, 4878 and references therein. (25) See reviews in Saykally, R. J. Science 1988,239, 157. Coe, J. V.; Saykally, R. J. In Ion and Cluster Ion Spectroscopy and Structure; Maier, J. P., Ed.; Elsevier: Amsterdam, 1989; p 131. Guteman, C. S.; Saykally, R. J. Annu. Res. Phys. Chem. 1984, 35, 387. (26) Good, A.; Durden, D. A.; Kebarle, P. J . Chem. Phys. 1970, 52, 212, 222. (27) Olovsson, I. J. Chem. Phys. 1968, 49, 1063. (28) Sears, T. J.; Bunker, P. R.; Davies, P. B.: Johnson, S. A,; Spirko, V. J. Chem. Phys. 1985, 83, 2676. (29) Liu, D.-J.; Oka, T.; Sears, T. J. J . Chem. Phys. 1986, 84, 1312. (30) Yeh, L. I.; Okumura, M.; Myers, J. D.; Price, J. M.; Lee, Y. T. J . Chem. Phys. 1989, 91, 7320. (31) Okumura, M.; Yeh, L. I.; Myers, J. D.; Lee, Y. T. J. Phys. Chem. 1990, 94, 3416. Yeh, Y . I.; Lee, Y. T.; Hougen, J. T. J . Mol. Spectrosc. 1994, 164, 473. (32) Shi, Z.; Ford, J. V.; Wei, S.; Castleman, A. W., Jr. J . Chem. Phys. 1993, 99, 8009. (33) Intermolecular Forces; Huyskens, P. L., Luck, W. A. P., ZeegersHuyskens, T., Eds.; Springer: Berlin, 1991. (34) The Hydrogen Bond Schuster, P., Zundel, G., Sandorfy, C., Eds.; North-Holland: Amsterdam, 1976; Vols. 1-3.
Bamett et al. (35) Hobza, P.; Zahradnik, R. Intermolecular Complexes: The Role of Van der Waals Systems in Physical Chemistry and in the Biodisciplines; Studies in Physical and Theoretical Chemistry, Vol. 52; Elsevier: Amsterdam, 1988. (36) Bell, R. P. The Proton in Chemistry, 2nd ed.; Come11 University Press: Ithaca, 1973. (37) Proton-Transfer Reactions; Caldin, E., Gold, V., Eds.; Chapman and Hall: London, 1975. Proton Transfer in Hydrogen-Bonded Systems; Bountis, T., Ed.; Plenum: New York, 1992. (38) See review by Landman, U.; Bamett, R. N.; Cheng, H.-P.; Cleveland, C. L.; Luedtke, W. D. In Computations for the Nano-Scale; Blochl, P. E., Joachim, C., Fisher, A. J., Eds.; Kluwer: Dordrecht, 1993; p 75. (39) Bamett, R. N.; Landman, U. Phys. Rev. B 1993,48, 2081. (40) Cheng, H.-P.; Bamett, R. N.; Landman, U. Chem. Phys. Lett., in press. (41) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University Press: Oxford, 1987. (42) Atomistic Simulations of Materials: Beyond Pair Potentials; Vitek, V., Srolovitz, D. J., Eds.; Plenum: New York, 1989. (43) Computer Simulations in Materials Science; Meyer, M., Pontikis, V., Eds.; Kluwer: Dordrecht, 1991. (44) Many-Atom Interactions in Solids; Nieminen, R. N., Puska, M. J., Manninen, M. J., Eds.; Springer: Berlin, 1990. (45) Simulations of Liquids and Solids; Ciccotti, G., Frenkel, D., McDonald, I. R., Eds.; North-Holland: Amsterdam, 1987. (46) Landman, U.; Luedtke, W. D.; Ouyang, J.; Xia, T. K. Jpn. J . Appl. Phys. 1993, 32, 1444. (47) See the collection of articles in ref 45. (48) Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, W. A. Intermolecular Forces; Clarendon Press: Oxford, 1987. (49) Foiles, S . M.; Baskes, M. J.; Daw, M. S. Phys. Rev. B 1986, 33, 7983. (50) Jacobsen, K. W.; Norskov, J. K.; Puska, M. J. Phys. Rev. B 1987, 35, 7423. Chetty, N.; Stokbro, K.; Jacobsen, K. W.; Norskov, J. Phys. Rev. B 1992, 46, 3798. (51) Car, R.; Paninello, M. Phys. Rev. Lett. 1985,55, 2471. See review by: Galli, G.; Paninello, M. In ref 43, p 283. (52) Bamett, R. N.; Landman, U.; Rajagopal, G. Phys. Rev. Lett. 1991, 67, 3058. (53) Barnett, R. N.; Landman, U. Phys. Rev. Lett. 1993, 70, 1775. (54) Cheng, H.-P.; Bamett, R. N.; Landman, U.Phys. Rev. B 1993.48, 1820. (55) Brechignac, C.; Cahuzac, Ph.; Carlier, F.; de Frutos, M.: Bamett, R. N.; Landman, U. Phys. Rev. Lett. 1994, 72, 1636. (56) Kaukonen, H.-P.; Bamett, R. N.; Landman, U. J. Chem. Phys. 1992, 97, 1365. (57) Rothlisberger, U.; Andreoni, W. J. Chem. Phys. 1991, 94, 8129. Car, R.; Paninello, M.; Andreoni, W. In Microclusters; Sugano, S., Nishina, Y., Ohnishi, S . , Eds.; Springer: Berlin, 1987. (58) Laasonen, K.; Andreoni, W.; Parrinello, M. Science 1992, 258, 1916. Rothlisberger, U.; Andreoni, W.; Paninello, M. Phys. Rev. Lett. 1994, 72, 665. (59) Feynman, R. P. Phys. Rev. 1939, 56, 340. (60) Kohn, W.; Sham, L. Phys. Rev. 1965, 140, A1133. (61) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (62) Bamett, R. N.; Landman, U.; Nitzan, A.; Rajagopal, G. J . Chem. Phys. 1991, 94, 608. (63) Saito, S . ; Cohen, M. L. Phys. Rev. B 1988, 38, 1123. Fernando, G. W.; Qian, G.-X.: Weinert, M.; Davenport, J. W. Phys. Rev. B 1989,40, 7985. (64) (a) Wood, D. M.; Zunger, A. J . Phys. A 1985, 18, 1343. (b) Davidson, E. R. J . Comput. Phys. 1M5, 17, 87. (65) Troullier, N.; Martins, J. L. Phys. Rev. B 1991, 43, 1993. (66) Kleinman, L.; Bylander, D. M. Phvs. Rev. Lett. 1982. 48. 1425. (67) Vosko, S . H.; Wilks, L.; Nusair, M.'Can. J. Phys. 1980,58, 1200. Vosko, S . H.; Wilks, L. J. Phys. C 1982, 15, 2139. (68) Ceperley, D. M.; Alder, B. 1. Phys. Rev. Lett. 1980, 45, 566. (69) (a) Becke, A. D. Phys. Rev. A 1988,38, 3098. (b) Becke, A. D. J. Chem. Phys. 1992, 96, 2155. (70) Perdew, J. P. Phys. Rev. B 1986, 33, 8822; 1986, 34, 7046. (71) Fan, L.; Ziegler, T. J . Chem. Phys. 1991, 94, 6057. Ortiz, G.; Ballone, P. Phys. Rev. B 1991, 43, 6376. (72) Martins, J. L.; Cohen, M. L. Phys. Rev. B 1988,37,6134. Troullier, N.; Martins, J. S. Phys. Rev. B 1991, 43, 8861. (73) Wentzcovitch, R. M.; Martins, J. L. Solid State Commun. 1991, 78, 831. (74) (a) Huber, K. P.; Herzberg, G. Molecular Spectra and Molecular Structure IV. Constants of Diatomic Molecules; Van Nostrand-Reinhold: New York, 1979. (b) Handbook of Chemistry and Physics: CRC Press: Cleveland, 1974. (75) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965.
Studies of Excess Electrons in Clusters (76) See the following reviews: (a) Chandler, D. In Liquids, Freezing and Glass Transition; Hansen, J. P., Levesque, D., Zinn-Justin, J., Eds.; Elsevier: Amsterdam, 1991. (b) Sprik, M. Quantum Simulations Using Path Integrals. In ref 43, p 305. (c) Gillan, M. J. In Computer Modeling of Fluids, Polymers and Solids; Catlow, C. R. A,, et al., Eds.; Kluwer: Dordrecht, 1990; p 155. (d) Beme, B. J.; Thirumali, D. Annu. Rev. Phys. Chem. 1986, 37, 401. (77) (a) Parrinello, M.; Rahman, A. J . Chem. Phys. 1984, 80, 860. (b) For applications to clusters see refs 6 and 7. Landman, U.; Barnett, R. N.; LUO,J.; Scharf, D.; Jortner, J. In Few-Body Systems and Multiparticle Dynamics; Micha, D. A., Ed.; AIP Conf. Proc. No. 162; AIP: New York, 1987; p 200. Cleveland, C. L.; Landman, U.; Bamett, R. N. Phys. Rev. E 1989, 39, 117. (78) See review by (a) Chandler, D.; Leung, K. Annu. Rev. Phys. Chem. 1994, 45. (b) Rossky, P. J.; Schnitker, J. J . Phys. Chem. 1988, 92, 4277. (c) Coker, D. F.; Beme, B. J. In Excess Electrons in Dielectric Media; Verrandini, C., Jay-Gerin, J. P., Eds.; CRC Press: Boca Raton, 1991; p 211. (79) (a) See reviews in refs 7, 10, and 79b. (b) Cleveland, C. L.; Landman, U.; Bamett, R. N. Phys. Rev. E 1989, 39, 117. (80) Gillan, M. J. Phys. Rev. Lett. 1987, 58, 563; Philos. Mag. A 1988, 58, 257. (81) (a) Landman, U.; Barnett, R. N.; Cleveland, C. L.; Norlander, P. In Tunneling; Jortner, J., Pullman, B., Eds.; Reidel: Dordrecht, 1986; p 269. (b) Mattson, T. R.; Engberg, U.; Wahnstrom, G. Phys. Rev. Lett. 1993, 71, 2615. (82) See review by Kuki, A. In Long-Range Electron Transfer in Biology; Springer: Berlin, 1991; p 49. Kuki, A.; Wolynes, P. G. Science 1987, 236, 1647. Wong, C. F.; Zheng, C.; Shen, J.; McCammon, J. A.; Wolynes, P. G. J. Phys. Chem. 1993, 97, 3100. (83) Ceperley, D. M.; Pollock, E. L. Phys. Rev. Lett. 1986, 56, 351. Ceperley, D. M.; Jacucci, G. Phys. Rev. Left. 1987,58,1648. See ref 79b. (84) Scharf, D.; Landman, U.; Jortner, J. J . Chem. Phys. 1987,87,2716. (85) Galli, G.; Andreoni, W.; Tosi, M. P. Phys. Rev. A 1986, 34, 3580. (86) Peterson, K. I.; Rao, P. D.; Castleman, A. W., Jr. J . Chem. Phys. 1983, 79, 777. (87) Kappes, M.; Radi, P.; Schar, M.; Schumacher, E. Chem. Phys. Lett. 1985, 113, 243. (88) Miller, T. M.; Leopold, D. G.; Murray, K. K.; Lineberger, W. C. J . Chem. Phys. 1986, 85, 2368. (89) Bergman, T.; Limberger, H.; Martin, T. P. Phys. Rev. Lett. 1988, 60, 1767. (90) Honea, E. C. Ph.D. Thesis, University of California, Los Angeles, 1990. (91) Pandey, R.; Sed, M.; Kunz, A. B. Phys. Rev. E 1990, 41, 7955. (92) Honea, E. C.; Homer, M. L.; Labastie, P.; Whetten, R. L. Phys. Rev. Lett. 1989, 63, 394. Honea, E. C.; Homer, M. L.; Whetten, R. L. Phys. Rev. B 1993, 47, 7480 and references therein. (93) (a) Pollack, S.; Wang, C. R. C.; Kappes, M. M. Chem. Phys. Lett. 1990, 175, 209; (b) Z . Phys. D 1989, 12, 241. (94) Yang, Y. A,; Conover, C. W.; Bloomfield, L. A. Chem. Phys. Lett. 1989, 158, 279. (95) Rajagopal, G.; Bamett, R. N.; Nitzan, A,; Landman, U.; Honea, E.; Labastie, P.; Homer, M. L.; Whetten, R. L. Phys. Rev. Lett. 1990, 64, 2933. (96) Weiss, P. W.; Ochsenfeld, C.; Ahlrichs, R.; Kappes, M. M. J. Chem. Phys. 1992, 97, 2553. (97) Xia, P.; Yu, N.; Bloomfield, L. A. Phys. Rev. E 1993, 47, 10040. (98) Xia, P.; Bloomfield, L. A. Phys. Rev. Lett. 1993, 70, 1779. Xia, P.; Cox, A. J.; Bloomfield, L. A. Z. Phys. D 1993, 26, 1841. (99) Sarkas, H. W.; Kidder, L.; Bowen, K. H. Preprint, 1994. 1100) Bonacic-Kouteckv, V.; Fuchs, C.; Gaus, J.; Pittner. J.; Kouteckv, J. Z . Phys. D 1993, 26, 192. (101) See review by Warren, W. W., Jr. In The Metallic and NonMetallic States of Matter; Edwards, P. P., Rao, C. N., Eds.; Taylor and
J. Phys. Chem., Vol. 99, No. 19, 1995 7753 Francis: London, 1985. (102) Selloni, A,; Fois, E. S.; Paninello, M.; Car, M. Phys. Scr. 1989, R5,261. (103) Martin, T. P. Phys. Rept. 1983, 95, 167. (104) Ochsenfeld, C.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 3487. (105) Varshui, Y. P.; Shukla, R. C. J. Mol. Spectrosc. 1965, 16, 63. (106) Mawhorter, R. J.; Fink, M.; Hartley, J. G. J . Chem. Phys. 1985, 83, 4418. (107) Wells, A. F. Structural Inorganic Chemistry, 5th ed.; Clarendon: Oxford, 1984; p 444. (108) Huheey, J. E. lnorganic Chemistry, 3rd ed.; Harper International, SI Edition: New York, 1983; p 67. (109) Williams, R. T.; Song, K. S.; Faust, W. L.; Leung, C. H. Phys. Rev. E 1986, 33, 7232. (110) Dickey, R. P.; Maurice, D.; Cave, R. J.; Mawhorter, R. J . Chem. Phys. 1993, 98, 2182. (111) See review on metallic clusters by de Heer, W. A. Rev. Mod. Phys. 1993, 65, 611. (112) Nakai, S.; Sagawa, T. J . Phys. SOC. Jpn. 1969, 26, 1427. (113) For a compilation of experimental values and discussion see Erwin, S. E.; Lin, C. C. J . Phys. C 1988, 21, 4285 (Table 1). (114) Poole, R. T.; Jenkin, J. G.; Liesegang, J.; Leckey, R. C. G. Phys. Rev. E 1975, 11, 5179. (1 15) Elliot, S. R. Physics ofAmorphous Materials; Longman: London, 1984; p 194. (116) The angular momentum character is calculated via analysis of the KS wave functions of the ith occupied state in spherical harmonics, is., Vp)= & , @ i m ( rY!m(Q), ) with the weight of the (1,m)component given by = J[&(r)l2? dr. In our analysis the expansion was centered about the electronic charge center of the wave function. (117) Berkowitz, J.; Batson, C. H.; Goodman, G. L. J . Chem. Phys. (Paris) 1980, 77, 631. (118) Nicol, G.; Sunner, J.; Kebarle, P. Int. J . Mass Specrrom. Ion Processes 1988, 84, 135. (119) Gillan, M. J. J. Phys. C 1987, 20, 3621. (120) Voth, G. A.; Chandler, D.; Miller, W. H. J. Chem. Phys. 1989, 91, 7749. See also Bader, J. S.; Kuharski, R. A.; Chandler, D. J. Chem. Phys. 1990, 93, 230. (121) Cunningham, A. J.; Payzant, J. D.; Kebarle, P. J. Am. Chem. SOC. 1972, 94, 7627. (122) Ng, C. Y.; Trevor, D. J.; Tiedemann, P. W.; Ceyer, S. T.; Kronebusch, P. L.; Mahan, B. H.; Lee, T. T. J. Chem. Phys. 1977, 67, 4235. (123) Collyer, A.; McMahon, B. J. Phys. Chem. 1983, 87, 909. (124) Dzidik, I.; Kebarle, P. J . Phys. Chem. 1970, 74, 1466. (125) See compilation by Th. Zeegers-Huyskens and P. Huyskens, in ref 33. p 23. (126) Meot-Ner, M.; Field, F. H. J. Am. Chem. SOC. 1977, 99, 998. (127) Kebarle, P.; Searles, S. K.; Zolla, A,; Scarbrough, J.; Arshadi, M. J. Am. Chem. SOC.1967, 89, 6393. (128) Payzant, J. D.; Cunningham, A. J.; Kebarle, P. Can. J. Chem. 1973, 51, 3242. (129) Searles, S. K.; Kebarle, P. J . Phys. Chem. 1968, 72, 742. (130) Tang, I. N.; Castleman, A. W., Jr. J. Chem. Phys. 1975,62,4576. (131) Meot-Ner, M. J . Am. Chem. SOC. 1984, 106, 1265. (132) Keesee, R. G.; Castleman, A. W., Jr. J. Phys. Chem. Ref. Data 1986, 15, 1011, 1040. (133) Kisiel, W.; Stankiewicz,B. Sur$ Sci. 1991,247,120 and references therein. Stoffel, N. G.; Riedel, R.; Colvita, E.; Margaritondo, G.; Haglund, R. F., Jr.; Taglauer, E.; Tolk, N. H. Phys. Rev. E 1985, 10, 6805. Singh, G.; Gallon, T. E. Solid State Commun. 1984, 51, 281. JB42706J
