Global Analytical Potential Energy Surface for Large Amplitude

J. Phys. Chem. B 2005, 109, 8439-8451

8439

Global Analytical Potential Energy Surface for Large Amplitude Nuclear Motions in Ammonia† Roberto Marquardt* and Kenneth Sagui Laboratoire de Chimie The´ orique, UniVersite´ de Marne-la-Valle´ e, 5, Bd Descartes (Champs sur Marne), F-77454 Marne-la-Valle´ e Cedex 2, France

Wim Klopper* Lehrstuhl fu¨r Theoretische Chemie, Institut fu¨r Physikalische Chemie, UniVersita¨t Karlsruhe (TH), D-76128 Karlsruhe, Germany

Martin Quack* Physical Chemistry, ETH Zu¨rich, CH-8093 Zu¨rich, Switzerland ReceiVed: February 10, 2005; In Final Form: February 22, 2005

An analytical, full-dimensional, and global representation of the potential energy surface of NH3 in the lowest adiabatic electronic state is developed, and parameters are determined by adjustment to ab initio data and thermochemical data for several low-lying dissociation channels. The electronic structure is calculated at the CASPT2 level within an [8,7] active space. The representation is compared to other recently published potential energy surfaces for this molecule. The present representation is distinguished by giving a qualitatively correct description of the potential energy for very large amplitude displacements of the nuclei from equilibrium. Other characteristic features of the present surface are the equilibrium geometries req(NH3) ≈ 101.24 pm, req(NH2) ≈ 102.60 pm, Req(NH3) ≈ 106.67°, and the inversion barrier at rinv(NH3) ≈ 99.80 pm and 1781 cm-1 above the NH3 minimum. The barrier to linearity in NH2 is 11 914 cm-1 above the NH2(2B1) minimum. While the quartic force field for NH3 from the present representation is significantly different from that of the other potential energy surfaces, the vibrational structures obtained from perturbation theory are quite similar for all representations studied here.

1. Introduction Understanding the dynamics of some prototype molecules has been crucial for the development of molecular physical chemistry, such as the studies on intra- and intermolecular energy transfer in CO2 pioneered by George Flynn and co-workers.1,2 Ammonia is a prototype molecule for largeamplitude nuclear motion. Stereomutation in general and the inversion motion of ammonia in particular have been extensively studied for more than 75 years,3-6 with results relevant in quite different but related fields such as MASER action,7 microwave and infrared spectroscopy,8-22 chirality,3,23,24 and the explicit time-dependent multidimensional quantum wave packet dynamics under coherent laser excitation.19 Ammonia and its isotopomers have also been interesting prototypical systems for the study of photodissociation dynamics25-32 and for the study of atomic insertion and molecular abstraction reactions.33,34 Such reactions are challenging because of the complexity due to avoided crossings, conical intersections, and intersystem crossings taking place in the asymptotic regions. The lowest dissociation channel of ammonia is the elimination of molecular hydrogen, with formation of azanediyl (“triplet” imidogen, nitrogenmonohydride, NH) †

Part of the special issue “George W. Flynn Festschrift”. * Authors to whom correspondence should be addressed. Fax: 0(033)1 60 95 73 20 (R.M.). E-mail: [email protected]; klopper@ chem-bio.uni-karlsruhe.de; [email protected].

-1 35 ≠ NH3(1A1) f NH(3Σ-) + H2(1Σ+ g ) ∆rH 0 ≈ 414 kJ mol (1)

This reaction is generally expected to be slow due to intersystem crossing. The next reaction channel in energy is the single bond fission reaction for the formation of ground-state azanyl (amidogen, nitrogendihydride, NH2(2B1))

NH3(1A1) f NH2(2B1) + H(2S(1/2)) ∆rH ≠0 ≈ 448 kJ mol-1 35,36 (2) followed on the energy ladder by the dissociation into the first excited 2A1 state of NH2

NH3(1A1) f NH2(2A1) + H(2S(1/2)) ∆rH ≠0 ≈ 590 kJ mol-1 35 (3) The 1A1 state of NH3 correlates diabatically with the 2A1 state of NH2 + H in eq 3. For planar ammonia, both channels cross at a conical intersection. If additional interactions such as spinorbit couplings or vibronic couplings are considered, then all three channels described above will be mixed. In an extended picture of the Born-Oppenheimer approximation, the set of points of lowest electronic energy in configuration space yields

10.1021/jp0507243 CCC: $30.25 © 2005 American Chemical Society Published on Web 04/14/2005

8440 J. Phys. Chem. B, Vol. 109, No. 17, 2005

Marquardt et al.

a thoroughly adiabatic surface for the nuclear motion that might be represented by a smooth, analytical function. In the present work, we have looked at modeling qualitatively such a hypersurface, without seeking to include couplings explicitely. Within the Born-Oppenheimer approximation, the topography of such a potential energy landscape, which might be largely influenced by the occurrence of conical intersections or intersystem crossings, can only be described approximately. We believe, however, that the analytical representation should be appropriate, for instance, to describe the recombination reaction

NH2(2B1) + H(2S(1/2)) f NH3(1A1)

(4)

at sufficiently low collision energies of vibronically cold NH2 molecules. In our picture, “adiabatically” low collision energies should be smaller than possible couplings between the channels. In ammonia, spin-orbit couplings might be significant, because a triplet state is in the vicinity of the conical intersection, as discussed below. A quite different approach is adopted in ref 37, where different seams of conical intersections of up to three states are treated simultaneously to discuss the photodissociation dynamics of electronically excited ammonia within a picture that could be called “diabatic” with respect to the picture used in the present work. Both treatments remain approximations, of course, since the nuclear motion problem cannot be solved exactly for any initial condition on single-valued potential energy surfaces (PES) in such cases. Here, we restrict the discussion to single-valued PES, since analytical representations for multivalued PES are much more complex.38 The alternative channel

NH2(2B1) + H(2S(1/2)) f NH3(3A1)

(5)

is probably also relevant for the recombination of azanyl with atomic hydrogen, because the first excited triplet state of ammonia has an energy similar to that of the conical intersection for the process of eq 2.33,34 For similar reasons, elimination of molecular hydrogen following eq 1 may also be a relevant reaction channel, because the extra activation energy for the reverse recombination 3 NH(3Σ-) + H2(1Σ+ g ) f NH3( A1)

(6)

is probably not very high. We finally mention here also the ionization channel 2 NH3(1A1) f NH+ 3 ( A′′ 2) + e

PES representations are normally obtained from an adjustment of parameters to data from electronic structure calculations, and the accuracy of the representation is related to the accuracy of the ab initio calculations in the first place. Therefore, a second criterion for adequate analytical forms requires them to be flexible. Recent, quite accurate experiments and calculations on infrared transition line positions from ab initio PES of ammonia17-22 give state-of-the-art results in this respect. The PES representations used in most calculations are high-order Taylor expansions in the neighborhood of the equilibrium molecular structure and normally not global. We may say that PES representations obtained so far from ab initio calculations are either accurate, but not global, or they are global, but not sufficiently accurate to describe the complete set of experimental data related to a global PES. The strategy to improve a PES representation a posteriori within an “experimental refinement”, in which some or, ideally, only a few parameters are varied appropriately, implies that, for highly flexible PES representations, the task of choosing the right parameters for the experimental refinement may be very tedious. Therefore, a third criterion for adequate analytical forms is that they are compact, which means flexible enough but with few parameters. Compact forms usually make use of a physically correct asymptotic behavior of the PES. The fourth criterion for adequate forms is robustness; small variations of parameter values should lead to small qualitative variations of a robust PES representation. Robust analytical forms are also important to ensure a physically correct behavior in regions of configuration space that are not well-sampled by data from electronic structure calculations. In the present work we concentrate on the exploratory derivation of a global representation of the PES of ammonia within the first step of our approach, i.e., basically from an adjustment to data from ab initio calculations at the CASPT2 level. We found this level of calculation to be an adequate compromise between efficiency, to cover a large range of geometries in the 6D configuration space, and to have sufficient accuracy for our purposes. We also used scaling of the ab initio data to experimental thermochemical data to facilitate the task of experimental refinement of the representation in a second step in the future. Alternatively, we might also consider to adjust the parameters to very accurate ab initio data, such as calculated in ref 18, including geometries in the intermediate and asymptotic region. Further work along these lines is in progress.40 Figure 1 shows ammonia with the definition of some coordinates used here.

(7)

the threshold of which has recently been determined with high accuracy.39 A global, analytical representation of the lowest adiabatic electronic state of ammonia has not been determined, to our knowledge (but see ref 19). Our general strategy is to first derive analytical forms that give a qualitatively correct representation of the topography of the lowest adiabatic PES. In this sense, it is reasonable to model the surface with a set of data obtained from a sufficiently low-level calculation of the electronic structure on a large grid of points in configuration space. In a second step, parameters might be adjusted to experimental data, to improve quantitatively the PES. As pointed out earlier, the first criterion for such forms is that they are global, i.e., allowing for a representation of a potential energy hypersurface that is analytical in the complete configuration space, including all possible reaction channels in a given energy range.

2. Theory 2.1. Ab Initio Calculations. The ab initio data set used to establish the general shape of the ground-state potential surface was largely built up with results from CASPT2 calculations. For our purposes, a rather small active space was chosen that is appropriate to describe the lone pair at the nitrogen atom, three NH bonding orbitals, and three NH antibonding orbitals. This generates an [8,7] active space for eight valence electrons. The 1s orbital of the nitrogen atom is excluded from the active space. The complete active space (CAS) yields 490 configurations in C1 symmetry, which is the smallest CAS needed to describe the complete atomization of NH3 into N and three H atoms. Calculations were done with equal weights for the two lowest electronic states in the CAS self-consistent field (SCF) step at each point of the spatial grid. More than 5000 nonequivalent points were calculated. For the calculations, a contracted [4s3p2d/3s2p] atomic basis set41 was used within the MOLCAS42 program package. Certain

PES for Large-Amplitude Nuclear Motion in NH3

J. Phys. Chem. B, Vol. 109, No. 17, 2005 8441 when bonds are broken or when bond lengths are much smaller than the equilibrium bond length. 2.2.1. Stretching Potential. The stretching potential is given by the expression 3

∑ i)1

Vstretch )

fs ys(ri)2 (9)

2

where ri is the ith bond length and

ys(r) )

Figure 1. Ammonia with coordinate definitions used in this work. The analytical representation is set up as a multidimensional function of the bond lengths ri and bond angles Rij; the out-of-plane angle χ is used in some of the graphical representations below.

data were calculated with the MOLPRO43 program package using the s, p, and d functions of the 6-311+G(2df,2pd) basis set for the nitrogen atom as well as the s and p functions for the hydrogen atoms, conserving the same contraction scheme given above. They were found to agree well with the results obtained with the MOLCAS program package. For instance, for a test geometry where all bonds are 100 pm long and all angles equal to 115°, the nondiagonal CASPT2 result from the MOLCAS calculations yields -56.436312Eh (1 Eh ≈ 4.360 × 10-18 J), roughly 1 mEh lower than the RSPT2 result of -56.435344Eh from MOLPRO. The lowest calculated energy is -56.439255Eh at equivalent bond lengths of 101.5 pm and angles of 106.5°. Geometry optimization yields req(NH3) ≈ 101.28 pm and Req(NH3) ≈ 106.65°, which are close to the equilibrium values found in recent theoretical work.15,17-19,21 Mixed theoretical/experimental estimates of the equilibrium values are req(NH3) ≈ 101.16 pm and Req(NH3) ≈ 107.25°, taken from an analysis made by Helgaker et al.44 of the experimental data of Benedict and Plyler.45 In ref 45, req and Req values were obtained from the usual analysis of the vibrational state dependence of rotational constants. The currently most accurate ab initio optimization has yielded req(NH3) ≈ 101.09 pm and Req(NH3) ≈ 106.76°, at the CCSD(T)(FC)/cc-pV6Z level including corrections for core correlation and connected triple and quadruple excitations.46 2.2. Global Analytical PES for Ammonia. The analytical PES for ammonia is based on a general form developed for molecules of the type AX1X2X3...XN.47 An application of this general form to describe the PES of methane was published and discussed previously.48,49 This form is basically a sum of three, positive definite terms

exp(-as[r - req]) as

(

1+

∑n

( ( ) ))

n exp -

r

n

(10)

is a generalization of the original Morse coordinate50 that allows for a more flexible description of both the anharmonicity of the stretching potential (by the parameter as) and the energy related to a single bond rupture

De )

fs 2a2s

(1 +

∑n n)2

(11)

The parameters rn in eq 10 should satisfy the condition rn . req, to ensure the interpretation of req as the equilibrium bond length. Asymptotically, say for r1 f ∞, the stretching potential becomes

where V′stretch is the stretching potential of the product molecule NH2 and De is the dissociation energy given in eq 11. The potential energy thus obtains asymptotically radial contributions that are of the type of a multipole expansion. The lowest power in the expansion could, for instance, be n ) 3, if the dissociation products have permanent dipole moments. In the present case of neutral products and one nonpolar product molecule, the lowest power is n ) 6, and the corresponding expansion parameter 6, which can be interpreted as being proportional to the dispersion constant C6, should be positive. For practical reasons, we have limited our description to use one or two expansion coefficients, e.g., 6 and 8. But in principle, no additional limitation is imposed.51 In order that V′stretch describes correctly the stretching potential of NH2, the parameters fs, as, etc. must change values. The stretching force constants of NH3 and NH2 are not very different. The variation can therefore be described by assuming the parameters entering V′stretch to be smoothly varying functions of the bond lengths. For instance, the stretching “force constant” fs may be written as a function

fs ) fs(r1,r2,r3;rsb,nsb) where Vs(i) and Vb(ij) are bond-stretching and bond-bending potentials, respectively. The indices i and j designate any of the peripheral atoms or atomic groups Xi or Xj. While Vb(ij) can be viewed as three-body potentials related to the interaction in the group XiAXj, the pair potential Vp(ij) describes the twobody interaction between Xi and Xj. Such pair interactions are expected to be small or at least less important than the bondbending potential for strongly binding systems having a distinguished central atom. But they may become important, either

rn

(13)

where rsb and nsb are switching parameters. Variations of parameters as a function of bond lengths may generally be monitored via switching functions, which normally need only one or two switching parameters to connect different subsets of parameters. Such multidimensional switching functions can be built up of one-dimensional 01-switching functions of the form

( ())

Sq(r;rs,ns) ) 1 - exp -

rs r

ns

(14)


Marquardt et al.

They have been described in detail in ref 48. The parameter rs in eq 14 can loosely be interpreted to be a threshold bond length that separates different regions of configuration space corresponding to different degrees of atomization. For ammonia, we distinguish three different subsets of parameters, e.g., for fs eq eq eq eq eq f(I) s ) fs(r (I),r (I),r (I)) r (I) ) r (NH3)

(15)

fs(II) ) fs(∞,req(II),req(II)) req(II) ) req(NH2)

(16)

fs(III) ) fs(∞,∞,req(III)) req(III) ) req(NH)

(17)

where the NHn (n ) 1, 2, 3) are always in their lowest electronic states. These functions are totally symmetric with respect to the permutation of the bond lengths, and we normally have ∂m fs/Πi ∂ki ri ≈ 0 for up to quite high orders m ) ∑i ki. The same idea applies to all other parameters with the exception of req, for which a different switching function is chosen such that ∂req/∂ri * 0, which generally allows for a good description of the coupling between the stretching modes already in the harmonic force field. 2.2.2. Bond-Pair Coordinates and Damping Functions. In contrast to bond-stretching potentials, there is no simple, compact analytical form that can be used to describe the anharmonicity of bending potentials. There is also a much larger choice of possible coordinates. We prefer to use general bondpair coordinates of the type eq (k) x(k) ij ) (cos(Rij) - cos(Rij ))yd (ri,rj,Rij)

(18)

where Rij is the angle between bonds i and j, with Rij ) Req ij at equilibrium and ri are bond lengths. Expressions using the cosine of valence angles have been used before (e.g., in refs 52-54). The functions y(k) d are given by

y(k) d (ri,rj,Rij) )

( (

exp -exp(ad(k)l (Rij))

( )( ( )(

× exp -exp(ad(k)l (Rij))

ri

nd1

ri

r rj

nd1

r rj

-1 eq

r

-1 eq

- ad(k)x (Rij) eq

)) ))

- ad(k)x (Rij) req

ndx

ndx

(19)

where nd1 and ndx are integer numbers. These functions have the value 1 for ri ) rj ) req and the asymptotic behavior y(k) d (ri f ∞) f 0 to damp out contributions to the bending potential of bonds that are broken. An exponential damping behavior was discussed previously.55-57 A polynomial as argument for the exponential function allows introducing an intermediate “antidamping” behavior, which turned out to be quite useful to decribe the PES in the vicinity of conical intersections, avoided crossings, or intersystem crossings. Larger polynomial expansions are possible, in principle. The PES of ammonia has a nonmononotic behavior as a function of the NH bond length in the region where the conical intersection occurs. One possibility to describe such a behavior is to consider a polynomial as argument for the exponential function in eq 19, as mentioned above, because an appropriate choice of the coefficients allows the induction of an intermediate antidamping behavior. We found it also necessary to consider varying the coefficients of the polynomial as a function of the corresponding bond angle. Both ad(k)l (Rij) and a(k) dx (Rij) in eq 19 are functions of the following type

1 ad(k)κ (Rij) ) (ad(k)κc(1 + cos(Rij))ndcos + ad(k)κo(1 - cos(Rij))ndcos) 2 κ ) l, x (20) (k) (k) The parameters adl(k) c and adlo are the values of adl (Rij) at Rij ) (k) (k) 0° and 180°, respectively; adxc and adxo are defined similarly. The index (k) indicates different sets of bond-pair coordinates that may be needed to increase the flexibility of the analytical representation. For instance, A- and E-type coordinates may have different damping behaviors or the damping may be different for different powers of the bending potential expansion to be discussed below. 2.2.3. Bending Potential. A flexible, yet robust, always positive definite form for a bending potential between bonds i and j is given by

Vb(ij) ) fbij p2n(xij)

(21)

where pn(x) is a homogeneous polynomial of order n

pn(x) ) x + ab2x2 + ‚‚‚ + abnxn

(22)

and xij are bond-pair coordinates such as described in eq 18. The coefficients abk are adjustable parameters that play the role of bending anharmonicity parameters. A major gain of flexibility is achieved if we allow the damping parameters in eq 18 to have different values depending on the power of the bond-pair coordinate in the bending coordinate expansion eq 22. The superindex (k) might be used to indicate this order dependence. In molecules with C1 symmetry, an independent set of bending anharmonicity parameters may be adjusted for each bond pair ij. However, for more highly symmetric molecules, such as NH3, these parameters are not independent. Symmetry may, indeed, lead to a strong a priori reduction of the number of independently adjustable parameters. Symmetry constraints may be considered by choosing symmetry-adapted analytical expressions. For methane, for instance, the space spanned by the six bond angle coordinates may be written in terms of one one-dimensional, one two-dimensional, and one three-dimensional irreducible representation of the Td point group.48 This reduces the number of quadratic bending force constants fbij in methane first from six to three and then, because of an additionally geometrical redundancy between the angular coordinates, to two independently adjustable parameters. The reader is referred to a recent discussion on geometrical redundancies and their implication on the use of dynamical variables in polyatomic molecules.58 In the ammonia molecule, the three bond-pair coordinates may be expressed as linear combinations of a one-dimensional and a two-dimensional symmetry-adapted coordinate, sb1 and sb2 ) (sb2a,sb2b), respectively, or vice versa

s b1 )

{

1 (x12 + x13 + x23) x3

1 (2x12 - x13 - x23) x6 s b2 ) 1 (x13 - x23) sb2b ) x2

(23)

sb2a )

(24)

We point out here that these symmetry-adapted coordinates are highly useful for the definition of the potential. All of the


J. Phys. Chem. B, Vol. 109, No. 17, 2005 8443

possible angular displacements can be modeled herewith, in particular, as will be shown below, the inversion motion. However, none of these coordinates can be directly interpreted as the inversion coordinate. The latter is rather given by the out-of-plane angle χ in Figure 1. The mathematical relation between the out-of-plane angle χ and the bond angles Rij is readily obtained by first transforming a molecular structure given in terms of χ into a set of Cartesian coordinates for the atomic positions and then transforming these coordinates into a set of bond lengths and bond angle coordinates. In the special case when C3V symmetry is conserved during the inversion (sb2a ) sb2b ) 0; R12 ) R13 ) R23 ≡ R), the relation can be cast into the simple mathematical form cos2(χ) ) (2 cos(R) + 1)/3. There are finally two distinguishable symmetry species of homogeneous polynomials to be used for the construction of a positive definite form of the kind shown in eq 21, one onedimensional of A symmetry and one two-dimensional of E symmetry. The general form adopted in this work is

1 1 Vbend ) fb1S2b1 + fb2(Sb22a + S2b2b) 2 2

(25)

with the following symmetry-adapted homogeneous polynomials 2 (S1) 2 2 Sb1 ) x3sb1 + a(S1) b21 sb1 + ab22 (sb2a + sb2b) + ‚‚‚

(26)

(S2)1 2 2 Sb2a ) x2sb2a + a(S2) b21 sb1sb2a + ab22 (sb2a - sb2b) + ‚‚‚ 2

(27)

(S2) Sb2b ) x2sb2b + a(S2) b21 sb1sb2b + ab22 (-sb2asb2b) + ‚‚‚ (28)

We have developed a program code that can be used in the computer algebra program MAPLE59 to calculate analytically the symmetry reduction of products of symmetry-adapted coordinates at any power. More details of this procedure will be published separately.60 The present model for the potential energy surface of ammonia includes symmetrized polynomial expressions to up to the fourth order, which are summarized in the Appendix. For the methane molecule, different symmetryadapted polynomial expressions have been published earlier.47,48 Upon rupture of a given bond, say for r1 f ∞, the remaining bending potential will contain terms involving x23 only, because of the damping functions used in x12 and x13. It may therefore be used to describe the bending potential of the product molecule NH2, if the parameters are varied accordingly. Therefore, as for the stretching potential, all parameters are considered to be smoothly varying functions of the bond lengths. The special choice ∂Req/∂ri * 0 allows for a good description of the coupling between the stretching and the bending modes already in the harmonic force field. Higher-order couplings can be described by the damping parameters adl introduced in eq 20, and in the pair potentials described below. All switching functions used in the present work have been described in ref 48. 2.2.4. Pair Potential. The direct (two-body) interaction potential between two peripheral atoms has been studied before (e.g., ref 52). For the purpose of the present work, it is sufficient to have a totally symmetric sum 3

Vpair )

Vp(ij) ∑ j>i

(29)

of Morse-type pair potentials of the form

Vp(ij) ) VMorse(ij) - D∞ij Sp(ri)Sp(rj)

∏ Sp(rki)Sp(rkj),

k*i,j

(30)

with 2 VMorse(ij) ) Dij(ri,rj,Rij)(1 - exp(-aij[rij - req ij ])) (31)

In eqs 30 and 31, rij is the interatomic distance between the atoms Hi and Hj, with rij ) req ij at equilibrium and

Dij(ri,rj,Rij) ) ∞ D0ij y(p) d (ri,rj,Rij) + Dij Sp(ri)Sp(rj)

∏ Sp(rki)Sp(rkj)

(32)

k*i,j

with 01-switching functions Sp(r) ≡ 1 - Sq(r;rdb,ndb) (the function Sq(r;rs,ns) is defined in eq 14); rdb and ndb are switching parameters for a “double bond enlongation”, which might differ from the parameters rsb and nsb for “single bond enlongation” that were introduced in eq 13. These switching functions have the properties Sp(r ≈ req) ) 0 and Sp(r f ∞) f 1. In eq 32, D∞ij is the dissociation energy of the isolated H2 molecule. The damping function yd(p) used in eq 32 has the same form as that given in eq 19; in this function, the special parameters (p) ad(p)lc, ad(p)lo, a(p) dxc, and adxo are used (see eq 20). The forms of eqs 30, 31, and 32 allow us to describe the formation of H2 from NH3 in the limit of simultaneous bond ruptures ri f ∞ and rj f ∞. However, to avoid the subsequent subtraction of binding energies in the asymptotic limit of three and more bond ruptures, which could lead to deep unphysical energy minima on the global potential surface, the H2 potential pertaining to the pair ij is evoked only if all interatomic distances rkj and rki, with k * i and k * j, are sufficiently large, which is guaranteed by the factor Πk*i,j Sp(rki)Sp(rkj) in eq 30. 3. Results and Discussion 3.1. Adjustment Procedure. The analytical forms described in section 2.2 above were adjusted to the ab initio data calculated as described in section 2.1 within a modified version of the Levenberg-Marquardt algorithm61 that allows us to consider further analytical relations between adjustable parameters during the fitting procedure by introduction of Lagrange multipliers as described in refs 47 and 48. The ab initio data might be weighted, the weights being reciprocal to the “uncertainties”

σ(E)/hc ) 100.048(lg((E-Emin)/hccm

-1))2.74

cm-1

(33)

which were estimated heuristically from the discussion of the PES of methane in ref 48. The ab initio data reported in the present work are expected to give a qualitatively correct picture of the global PES of ammonia. We do not expect, however, to obtain quantitatively accurate values for the equilibrium geometries and the dissociation energies at the several dissociation limits or for special energies such as the inversion barrier. Inclusion of electronic correlation decreases as a function of the degree of atomization in the present case. In the limit N + 3H, our CAS-SCF results correspond to the Hartree-Fock level of theory. To obtain a more realistic model potential, we have considered to scale the ab initio data prior to the adjustment procedure, both NH bond lengths and energies. Such a strategy was already adopted for the derivation of the PES for methane.48 In that case, scaling factors were determined by comparison to experimental data. Here, the NH bond distances are multiplied with the factor 0.999615, which is the ratio of the equilibrium NH bond length 101.24 pm calculated at the CCSD(T)/cc-pVQZ level in ref 15, to the


Marquardt et al.

TABLE 1: Parameters Used in the Present Global Analytical Representation and Adjusted Parameter Values AMMPOT1 NH3 Vstretch Parameters ar /Å-1 re/Å Fs/aJ Å-2 as /Å-1 6 8 rs/Å

2.4560360 1.0124000 6.3867210 2.5439580 0.2305466 0.0000000 2.1865750

1.9888940 1.0240000 6.1978500 2.2254780 0.0488106 0.0000000 2.0000000

Vpair Parameters D0ij /aJ

0.3219573

0.0647524

a0ij /Å-1 r∞ij /Å a∞ij /Å-1 D∞ij /aJ

AMMPOT2

NH2

NH 1.9888940 1.0380000 6.5478080 1.9231320 0.0488106 -0.3195010 2.0000000 a

NH3

NH2

2.0948560 1.0124000 6.3867220 2.1908770 0.1071329 0.0000000 2.0000000

1.6844120 1.0260000 6.3253380 2.1855660 0.0349995 0.0000000 2.0000000

0.4458559

1.0000000

0.9000000

2.0000000

0.7647939

0.5064664

0.7414000

0.7414000

0.7414000

0.7414000

1.9453400

1.9453400

1.9453400

1.9453400

0.7605968

0.7605968

0.7605968

0.7605968

Vbend Parameters ac /Å-1 1.3388740 ce -0.2868590 Fb1/aJ 0.0639957 Fb2/aJ 0.1374017 -0.2992830 ab(S1)

1.3388740 -0.2249511 0.2611156 0.0000000 -0.6000000

1.6570470 -0.2868590 0.0639957 0.1374017 -0.0517006

1.6570470 -0.2368381 0.1316461 0.1000000 -0.2961002

ab(S221)

-0.7163208

0.0000000

-1.6316170

-0.2961002

ab(S311)

-0.4554199

0.1631592

-0.2171079

0.0333953

ab(S321)

0.0000000

0.0000000

-0.3605843

0.0333953

ab(S331)

-0.2812387

0.0000000

0.6491018

0.0333953

ab(S212)

0.0000000

0.0000000

0.3746296

0.0000000

ab(S222)

-0.7088556

0.0000000

0.0657369

0.0000000

ab(S312)

0.0000000

0.0000000

0.0498214

0.0000000

ab(S322)

0.0000000

0.0000000

-1.0329350

0.0000000

ab(S332)

21

-0.1434478

0.0000000

-0.2760892

0.0000000

Damping Parametersb 1.8932800 ad(a1) xc

0.9436250

1.3918580

1.0000000

ad(a2) xc

3.8415090

3.1621130

1.3918580

1.0000000

ad(a3) xc

1.0000000

1.0000000

1.3918580

1.0000000

ad(a1) xo

1.8932800

0.9436250

2.2948900

1.0000000

ad(a2) xo

3.8415090

3.1621130

2.2948900

1.0000000

ad(a3) xo

1.0000000

1.0000000

2.2948900

1.0000000

ad(a1) lc

-0.6685000

-0.6457140

1.4491060

2.0000000

ad(a2) lc

-0.6685000

-0.6457140

0.0000000

2.0000000

ad(a3) lc

-0.6685000

-0.6457140

0.0000000

2.0000000

ad(a1) lo

-0.6685000

-0.6457140

0.0000000

-1.0000000

ad(a2) lo

-0.6685000

-0.6457140

0.0000000

-1.0000000

ad(a3) lo

-0.6685000

-0.6457140

0.0000000

-1.0000000

ad(e1) xc

1.8932800

0.9436250

2.9111283

1.0000000

ad(e2) xc

3.8415090

3.1621130

2.9111280

1.0000000

ad(e3) xc

1.0000000

1.0000000

2.9111280

1.0000000

ad(e1) xo

1.8932800

0.9436250

-0.8862324

1.0000000

ad(e2) xo

3.8415090

3.1621130

-0.8862324

1.0000000

ad(e3) xo

1.0000000

1.0000000

-0.8862324

1.0000000

ad(e1) lc

-0.6685000

-0.6457140

-0.5275110

2.0000000

NH 1.6844120 1.0363500 6.3928490 2.1560860 0.0023672 -0.1100000 2.1500000 a



TABLE 1: (Continued) AMMPOT1 NH3

AMMPOT2

NH2

NH

NH3

NH2

Damping Parametersb -0.6685000 ad(e2) lc

-0.6457140

0.0000000

2.0000000

ad(e3) lc

-0.6685000

-0.6457140

0.0000000

2.0000000

ad(e1) lo

-0.6685000

-0.6457140

0.0000000

-1.0000000

ad(e2) lo

-0.6685000

-0.6457140

0.0000000

-1.0000000

ad(e3) lo adpxc adpxo adplc adplo nd1

-0.6685000

-0.6457140

0.0000000

-1.0000000

c c c c 1.0000000 1.0000000 1.0000000

c c c c 1.0000000 1.0000000 1.0000000

1.0000000 1.0000000 -0.0563178 -1.0994810 2.0000000 1.0000000 1.0000000

1.0000000 1.0000000 0.0000000 -1.0000000 2.0000000 1.0000000 1.0000000

6.0000000 2.0000000 6.0000000 2.0000000

6.0000000 2.0000000 6.0000000 2.0000000

ndx ndcos

Switching Parameters nsb 8.0000000 rsb/Å 1.9209602 ndb 6.0000000d rdb/Å 1.7000000d

8.0000000 1.9209602 6.0000000d 1.7000000d

8.0000000 1.9209602

NH

6.0000000 2.0000000

a Undefined parameters were either found to be effectively irrelevant for the fit or cannot be defined within certain subsets; where appropriate, the value 0.0 can be set. b The index p on adx parameters indicates parameters used in Vpair and indices a or e indicate parameters used in A- and E-type expressions of power n in Vbend, respectively. c These parameters are not used in AMMPOT1, for which the expression Sq(ri)Sq(rj) was used (p) instead of the expression yd in eq 32. d In AMMPOT1, different switching parameters were used to vary the equilibrium distance of a pair of H atoms and the Morse parameter between a0ij and a∞ij ; these parameters are rrij ) 1.35 Å, nrij ) 12 and raij ) 1.80 Å, naij ) 4.

equilibrium value 101.28 pm obtained in section 2.1. Energies relative to the minimum of the PES were multiplied by the factor the exp Dexp ≈ 41 051 cm-1 is an estimae /De ) 1.0899, where De tion of the experimental value for the dissociation energy of reaction eq 2 from the thermochemical data of ref 35, including differences of the zero point energy in the harmonic approxima-1 is the value for this energy from tion, and Dthe e ≈ 37 665 cm the present ab initio calculations. The values for the equilibrium bond angles from the present data, Req(NH3) ≈ 106.65° are very close to those reported in ref 15; thus bond angles were not further scaled here. In addition to scaling the geometries, we considered the following set of constraints. The harmonic force field should be identical to that given in ref 15; the dissociation energy should be equal to the value for Dexp e determined above; the potential energy at the inversion barrier geometry found in ref 17 should coincide with the inversion barrier determined there. For this purpose, eight analytical conditions involving several parameters introduced in section 2.2 were developed as described in ref 48 and used in connection with Lagrange multipliers during the fitting. A grand set of parameters was finally obtained by first finding an optimal subset of parameters p(III) suitable to describe the potential for (3Σ-)NH in the dissociation limit of eq 1. Instead of using the ab initio data from the present work for this purpose, we used the potential function obtained in ref 62. We next adjusted the subset of parameters p(II) that describe the potential for (2B1)NH2 in the dissociation limit eq 2. Here too, instead of using the present ab initio data, we adjusted these parameters to a potential for NH2 reported in ref 63. In the final step, we adjust only the subset of parameters p(I) pertaining to the parent molecule (1A1)NH3. The scaled data described above were used for the adjustment, in this case. In this step, the switching parameters needed to connect the different subsets were also adjusted.

Two resulting parameter sets, AMMPOT1 and AMMPOT2, are collected in Table 1. Parameter set AMMPOT1 was used in our previous work19,64 and was obtained from an adjustment to weighted data with weights as given in eq 33. The square root of the sum of the squares of weighted residuals (the “weighted root mean square (rms)”) is on the order of 64 cm-1 for more than 5000 adjusted energies falling in the range below correspondingly 60 000 cm-1. The corresponding sum of the squares of unweighted residuals is 1590 cm-1. Parameter set AMMPOT2 was obtained from an adjustment to ab initio data with equal weights and is the more refined development of the global PES for ammonia in the present work. The sum of the squares of unweighted residuals is 1275 cm-1. Although symmetrized bending expressions were developed in the present work up to the fourth order, we found it sufficient to adjust the bending “anharmonicity” parameters ab up to the third order only, in both representations given here. 3.2. Some Characteristic Features of the PES. Figures 2-5 show one-dimensional sections of recently published, analytical, six-dimensional representations of the PES for NH3. Figure 2 shows functions of one NH bond length with bond angles Rij ) 80° (Figure 2a; the out-of-plane inversion angle is χ ) 132.1°), 106.5° (Figure 2b; χ ) 112.3°), and 120° (Figure 2c; χ ) 90°). The scaled CASPT2 data are shown with marks, where a ] indicates the lower and a * the higher of the two lowest states from the CAS calculation. The representations from refs 21 and 22 are of the form of a polynomial expansion in bond length and bond angle coordinates that have been developed to give highly accurate descriptions of the infrared spectrum in the mid- to near-infrared spectral range. However, they have a limited range of validity, typically for bond length displacements from equilibrium smaller than 50 pm. In contrast to these representations, the analytical representation developed in the present work (AMMPOT2), shown as a continuous line in Figure 2 as well as that from ref


Marquardt et al.

Figure 3. Decomposition of the potential curve shown in Figure 2c into the three terms of eq 8. The different symbols and curves correspond to CASPT2 data from the present work (AMMPOT2, as in Figure 2c): s present work (from eq 8, as in Figure 2c); - - Vstretch from eq 9; - ‚ - Vbend from eq 25; - - - Vpair from eq 29.

Figure 2. One-dimensional sections of the potential energy surface of ammonia as a function of one NH bond length. (a) Pyramidal ammonia. Values for the remaining coordinates are bond lengths at roughly 101.2 pm and bond angles at 80.0°; the inversion angle out of the planar structure is χ ) 132.1°. (b) Pyramidal ammonia. Values for the remaining coordinates are bond lengths at roughly 101.2 pm and bond angles at 106.5°; this is almost the equilibrium bond angle. The inversion angle out of the planar structure is χ ) 112.3°. (c) Planar ammonia. Values for the remaining coordinates are bond lengths at roughly 101.2 pm and bond angles at 120°; the inversion angle is χ ) 90°. The different symbols correspond to CASPT2 data from the present work, scaled as described in section 3.1; ] designates the lowest state and * the first excited state: s present work (AMMPOT2); - Rajama¨ki et al.;22 - - - Leonard et al.;21 - ‚ - Lin et al.18

18 are suitable to describe bond fisson. The representation from ref 18 uses Morse potential functions for the stretching modes and was developed from adjustments to data from high-quality ab initio calculations at the coupled cluster level of theory with extrapolation to the basis set limit in the neighborhood of the equilibrium structure. This representation yields also a very

accurate description of the infrared spectrum up to the first and second overtones, in addition to having a global character. The representation from the present work describes the lowest adiabatic electronic state. For planar ammonia, the corresponding PES has a transition structure in the region where the channels given by eqs 2 and 3 have a conical intersection. The conical intersection is indicated in Figure 2c to occur at a bond extension of approximately 200 pm. The representation from ref 18 seems to correlate with the higher-lying (2A1) NH2 state. Not unexpectedly, when data in the vicinity of the equilibrium structure are used to adjust the parameters of the stretching potential, the diabatic potential surface can be described correctly with simple Morse potential functions, at least qualitatively. For the description of the global, lowest adiabatic potential surface, a more flexible form such as that described in eq 10 is needed. This form guarantees the correct asymptotic limit, at least qualitatively with a -1/r6 potential. However, a detailed analysis showing the decomposition into all three terms of eq 8 reveals, in Figure 3, that the transition structure in the planar configuration is mainly given by contributions from the pair potential. The present analytical representation was developed such that these contributions are decreased when the molecule bends. As indicated by the shape of the continuous line in Figure 2b, in the presently best representation AMMPOT2, these contributions are still apparent at the equilibrium bond angle for ammonia. Although the several potential representations shown in Figure 2 are quite different, the corresponding force fields at the equilibrium geometry are rather similar to up to the third order, as shown in Table 2. These force constants were defined, and their values were obtained numerically, for all representations discussed here, with a finite difference method as described in ref 49. Only the F122 and F222 force field coefficients from the PES developed in this work differ significantly from those of other representations. Some coefficients of fourth order change quite substantially from one representation to the other, the values for F1111 and F1112 from ref 22 being particularly different. Quite different quartic force fields may thus yield quite similar graphical representations in the vicinity of the equilibrium structure, which underlines the fact that force constants do not necessarily always have a physically unique meaning. However, this discussion must be extended to also include the geometry and vibrational structure. As described above, the PES developed in this work was constrained to have the equilibrium geometry and a harmonic force field as that of ref 15 (CCSD(T)/ccpVQZ).



TABLE 2: Quartic Force Field of Ammonia in the Electronic Ground Statea

req/Å Req b F11/aJ Å-2 F12/aJ Å-1 F22/aJ F3a3a/aJ Å-2 F3a4a/aJ Å-1 F4a4a/aJ F111/aJ Å-3 F112/aJ Å-2 F122/aJ Å-1 F13a3a/aJ Å-3 F13a4a/aJ Å-2 F14a4a/aJ Å-1 F222/aJ F23a3a/aJ Å-2 F23a4a/aJ Å-1 F24a4a/aJ F3a3a3a/aJ Å-3 F3a3a4a/aJ Å-2 F3a4a4a/aJ Å-1 F4a4a4a/aJ F1111/aJ Å-4 F1112/aJ Å-3 F1122/aJ Å-2 F113a3a/aJ Å-4 F113a4a/aJ Å-3 F114a4a/aJ Å-2 F1222/aJ Å-1 F123a3a/aJ Å-3 F123a4a/aJ Å-2 F124a4a/aJ Å-1 F13a3a3a/aJ Å-4 F13a3a4a/aJ Å-3 F13a4a4a/aJ Å-2 F14a4a4a/aJ Å-1 F2222/aJ F223a3a/aJ Å-2 F223a4a/aJ Å-1 F224a4a/aJ F23a3a3a/aJ Å-3 F23a3a4a/aJ Å-2 F23a4a4a/aJ Å-1 F24a4a4a/aJ F3a3a3a3a/aJ Å-4 F3a3a3a4a/aJ Å-3 F3a3a4a4a/aJ Å-2 F3a3a4b4b/aJ Å-2 F3a3b4a4b/aJ Å-2 F3a4a4a4a/aJ Å-1 F4a4a4a4a/aJ

from ref 15

from ref18

from ref21

from ref 22

this workc

1.01240 1.85320 6.95959 0.52373 0.54312 7.07046 -0.20045 0.69467 -26.15634 -0.17502 0.07513 -25.99714 0.24514 -0.34863 -0.15694 0.10298 0.18061 -0.54058 -18.45040 -0.37483 0.06966 -0.04432 83.02676 -0.98588 0.88696 85.64277 0.42450 -0.13903 0.72614 -0.58988 -0.13530 0.20123 61.10549 0.34261 0.04796 0.09635 -0.21718 -0.01396 -0.23835 -0.20015 -0.57754 0.21735 -0.01690 -0.48039 130.62156 0.65926 -0.31791 -0.79120 0.24058 -0.29254 0.19583

1.01029 1.86319 7.03967 0.49382 0.51320 7.15609 -0.18814 0.66667 -26.34082 -0.10358 0.08344 -26.12849 0.20492 -0.34865 -0.19755 0.16044 0.17435 -0.53073 -18.48626 -0.41673 0.07662 -0.06312 84.73963 -0.36294 1.09496 85.80742 0.27907 0.03574 0.84866 -0.25333 -0.10108 0.27931 60.98271 0.24535 0.11313 0.11566 -0.22920 0.11166 -0.34557 -0.06542 -0.47444 0.32046 -0.04583 -0.43282 130.82443 0.26165 -0.12068 -0.87447 0.37689 -0.38976 0.34469

1.01090 1.86323 7.02511 0.55092 0.51848 7.12481 -0.18209 0.66748 -26.92136 -0.03062 0.00200 -26.25481 0.18195 -0.33905 -0.25955 0.25692 0.14361 -0.49881 -19.39907 -0.39065 0.05182 -0.04444 84.44085 0.56051 1.24976 88.63932 0.39690 0.19555 -0.20716 -1.69701 0.06292 0.28945 56.19193 1.06002 -0.17812 0.31268 -0.49316 -0.44531 -0.24105 -0.22177 1.07745 0.10939 -0.20620 -0.48238 136.16906 -2.83254 0.28558 -0.41627 0.30203 -0.50049 0.34121

1.01280 1.85942 6.79623 0.48532 0.50068 7.10125 -0.19197 0.67186 -24.24539 -0.40350 0.04211 -26.05830 0.22629 -0.33845 -0.26392 0.12507 0.17580 -0.53018 -18.49670 -0.37178 0.06818 -0.05932 114.77737 -12.06942 0.71897 85.37994 0.24254 0.03201 -0.18908 -0.47000 -0.18557 0.27294 60.82228 0.39667 -0.00756 0.14924 0.23838 0.03816 -0.25830 -0.14000 -0.56901 0.21693 -0.03301 -0.44869 130.27444 -0.00017 -0.24021 -0.76536 0.26257 0.00000 0.22292

1.01240 1.85878 6.96000 0.52401 0.54300 7.07101 -0.20000 0.69500 -22.60064 -0.45018 -0.36922 -24.68577 0.23138 -0.34157 -0.49671 0.00862 0.20621 -0.22973 -17.85480 -0.16584 0.11674 -0.51710 63.01759 0.19814 1.34921 69.10127 -0.12095 -0.48869 0.05397 -0.47739 -0.24331 0.11940 49.72865 0.33658 0.18207 0.06397 -0.99764 1.21149 -0.11164 -0.49945 -0.62310 0.10297 -0.08597 -0.44790 116.80312 0.17054 -0.62356 -1.32689 0.35166 -0.03291 0.68758

a For a definition of the force field, see ref 69. The equilibrium geometry is from a present numerical evaluation of published potential energy surfaces. 1 Å ) 100 pm. b Angles are given in radians. c AMMPOT2.

The quartic force fields of the PES representations given in Table 2 being quite different one from the other, one could expect the vibrational spectrum to differ already in the fundamental transitions. Table 3 shows, however, that all PES

representations discussed here yield rather similar fundamental transitions for NHD2, when perturbation theory is used. In this table, we give, for brevity of the comparison, only the lowest tunneling levels (see also ref 20). We explain remaining differences to the experimental data on the order of 2-30 cm-1 with the lack of precision of perturbation theory. Only ν˜ 3a and ν˜ 4a from the evaluation of ref 21 differs substantially from the general trends. The corresponding harmonic wavenumbers are 2552 and 1270 cm-1, respectively; their ratio is close to 2:1, and perturbation theory clearly breaks down in this particular case as also discussed in refs 19 and 64. Where available, sixdimensional variational calculations are also reported. One might expect, on the basis of a comparison between results from variational calculations and perturbation theory, that all representations have a similar accuracy, at least in regions of low vibrational excitations. A concluding analysis of the PES developed here requires a six-dimensional variational treatment of the vibrational structure, which is the subject of separate publications.19,40,64,65 In some of our previous work, we have compared approximate and full-dimensional dynamical results for the isotopomers NHD2 and NH3 using the potential AMMPOT1. For instance, tunneling splittings in the fundamental transition of the inversion motion are 0.12 cm-1 (0.17 cm-1) for NHD2 19 and 0.55 cm-1 (0.79 cm-1) for NH3 64 without any specific adjustment of the surface to the experimental values, indicated in parentheses. These results have shown that given the global structure of our surface the success is quite good. Further work using AMMPOT2 is in progress. Figure 4 shows sections of PES representations along the inversion coordinate (the out-of-plane angle χ; planar ammonia is obtained with χ ) 90°; see also the discussion following eq 23 above). All bond lengths are kept at 100 pm, close to their equilibrium values (of roughly 101.2 pm). In Figure 4a, C3V symmetry is conserved, leading to regular pyramidal structures upon displacement of χ from 90°. The inset shows the projections of the NH bonds on the inversion plane denoted by the coordinates x and y. They form the angles 180°, 60°, and -60° with the x axis. The inversion barrier of the present analytical PES representation is at 1781 cm-1 and bond lengths 99.8 pm, which are close to the values obtained in ref 17: 1788 cm-1 at 99.4 pm. As mentioned above, the coincidence of the potential energies at the barrier geometry from ref 17 was used as an additional constraint during the adjustment procedure. For the other representations shown in this figure, the barriers are at 1769 cm-1 (99.4 pm),21 1791 cm-1 (99.6 pm),22 and 1790 cm-1.18 In Figure 4b, the symmetry is lowered to C1 with highly distorted pyramids, in which the projections of the NH bonds on the inversion plane form the angles 180°, 90°, and -30°, respectively. Figure 4c has symmetry Cs with projections of the NH bonds on the inversion plane that form the angles 180°,

TABLE 3: Fundamental Vibrational Transitions in NHD2 theory from Lin et al.18

from Le´onard et al.21

this workb

PTc

PTd

PTc

vare

PTc

PTc

PTc

3373.8 851.9 2426.5 2537.5 1248.8 1472.8

3377.6 871.5 2419.2 2531.2 1246.9 1466.6

3378.9 813.5 2416.7 2544.0 1229.8 1454.5

3407.1 810.2 2430.1 2558.0 1233.1 1461.6

3398.6 815.4 2584.6 2556.7 1232.0 1457.8

3380.1 822.7 2809.5 2547.6 1375.8 1458.5

3387.6 836.6 2436.5 2544.9 1251.8 1475.7

from Martin et -1

ν˜ 1/cm ν˜ 2/cm-1 ν˜ 3a/cm-1 ν˜ 3b/cm-1 ν˜ 4a/cm-1 ν˜ 4b/cm-1

Γ (Cs)

experiment

A′ A′ A′ A′′ A′ A′′

3404.220 810.270 2430.820 2559.820 1233.420 1461.820

a

al.15

from Rajama¨ki et

al.22

a The fundamental wavenumbers of the lower tunneling level are given; see also ref 20 for data for the higher level. b AMMPOT2 c Perturbation theory following ref 71, p 160, with corrections as discussed in ref 72. d Perturbation theory following ref 15. e Variational calculation following ref 22.


Figure 4. One-dimensional sections of the potential energy surface of ammonia as a function of the inversion angle. The insets show the projections of the NH bonds on the inversion plane having the axes x and y. (a) Bond lengths are at roughly 101.2 pm, and bond angles change upon conservation of C3V symmetry (regular pyramids). (b) As in part a, but the structure has C1 symmetry, and pyramids are highly distorted. (c) As in part a, but bond angles change upon conservation of Cs symmetry (planar ammonia forms a “T” and pyramids are distorted). See also text and Figure 2.

90°, and -90°. All representations describe the ab initio data well; the representation from the present work slightly overestimates the energies at planar but distorted structures in Figure 4b. All representations describe the ab initio data in a qualitatively satisfactory way. In comparison to the global representation from the present work, which needs only on the order of 30 adjustable parameters and is likely to be more robust, the representation from ref 18 needs between 73 and 78 parameters. The representation from ref 22 is highly complex and requires on the order of 850 parameters. Due to a misinterpretation of these

Marquardt et al.

Figure 5. One-dimensional sections of the potential energy surface of planar ammonia as a function of the bond angle R23 ) R. Bond lengths are such that r2 ) r3 ) 100 pm, and bond angles are R12 ) R13 ) (180° - R)/2. (a) Bond length r1 ) 150 pm. (b) Bond length r1 ) 300 pm. (c) Bond length r1 ) 500 pm. See also text and Figure 2, where references to the potentials described by the various symbols are given.

parameters, the graphical representations of the PES from ref 22 discussed in ref 66 are incorrect. The PES from ref 22 yields accurate values for fundamental and low overtone transitions and combination bands of ammonia and isotopomers, but its complexity and the large number of parameters decreases its robustness and renders it less appropriate to ensure even a qualitatively correct description of the PES in regions of configuration space that have not been sufficiently well-sampled, as shown in Figure 5a. Figure 5 shows sections of the PES for planar ammonia, in which the angle R23 between the bonds 2 and 3 is varied at different values of the opposite NH bond r1. At a bond length r1 ≈ 150 pm (Figure 5a), which is close to the switching

PES for Large-Amplitude Nuclear Motion in NH3 parameter rsb (see eq 13 and Table 1), parameters start to switch from their values in set I (NH3) to their values in set II (NH2). In this region, the global PES is indeed a mixture representation with parameter sets that describe simultaneously the NH3 and the NH2 moieties and approximate mixture coefficients of 96% and 4%, respectively. These coefficients are given by Sq and 1 - Sq from eq 14 with r ) 150 pm, rs ≡ rsb ≈ 200 pm, and ns ≡ nsb ) 6. At r1 ) 300 pm (Figure 5b), the mixture of parameter sets is 8% for the NH3 moiety and 92% for NH2; the representaion is thus almost that of the NH2(2B1) PES. In Figure 5c, for r1 ) 500 pm, the representation shows the double minimum PES of NH2 with a barrier to linearity of roughly 12 000 cm-1. This barrier corresponds to previous findings.25 This value was indeed considered implicitly during the adjustment procedure by constraining the parameters of subset II to obey the appropriate analytical relation for the barrier to linearity within our modified version of the Marquardt algorithm as described in section 3.1 above. We point out though that the Renner-Teller effect of linear NH2 was not explicitely considered in the CASPT2 method used in the present work and shown in Figure 5. At r ) 500 pm, the present PES representation therefore already describes the dissociation limit in eq 2. Clearly, none of the other representations shown in Figure 5 is capable to give even qualitatively the correct behavior of the PES in this region of configuration space, which is likely to be important for the description of single bond fission and formation reactions of ammonia. Further characteristic features of the present PES representation of ammonia will be discussed in a future publication.60 Work is also in progress to extend the description of the potential by means of a much improved ab initio data set.40 4. Conclusions Potentials for multidimensional dynamics in polyatomic molecules open new dimensions to our theoretical understanding of spectroscopy and kinetics.67,68 In particular, analytical representations of global potential energy hypersurfaces are important tools for the calculation of spectroscopic, structural, and dynamical properties of polyatomic molecules.48,49 The derivation of global and compact analytical representations of PES may occur in two or more steps, the first of which is the adjustment to data from ab initio calculations. In a second step, parameter values may be “experimentally refined” either by new, independent adjustments, with the previously determined set of parameter values as a highly valuable starting point, or by simultaneous consideration of experimental conditions during the fit and analytical relations between the adjusted parameters. Within this approach, it is interesting to perform ab initio calculations at a moderate level of theory in a large region of configuration space. Alternatively, one may try to develop global analytical representations of the PES from adjustments to highlevel ab initio data right away.18,40 In any case, it is important to have robust, compact analytical forms to represent the PES. The present work reports results for a global analytical representation of the PES for ammonia that are based on ab initio calculations at the CASPT2 level of theory. The comparison with other recently published representations of this PES18,21,22 show that a global representation with a rather small number of adjustable parameters is feasible. Some applications of the AMMPOT1 PES of ammonia have been reported previously.19,47,64 These dynamical calculations on AMMPOT119,47,64 and their comparison with experiment20,64 have already shown the success of the present approach to formulating

J. Phys. Chem. B, Vol. 109, No. 17, 2005 8449 the global analytical hypersurface. Further work on dynamical calculations using AMMPOT2 is underway.65 For AMMPOT2, we may conclude, in particular: 1. In contrast to the PES from refs 18, 21, and 22, the present representation gives a qualitatively correct description of the potential energy in several low-lying dissociation channels. 2. It describes both the equilibrium geometries with req(NH3) ≈ 101.24 pm, req(NH2) ≈ 102.60 pm, Req(NH3) ≈ 106.67°, and the inversion barrier at rinv(NH3) ≈ 99.80 pm with an energy coresponding to 1781 cm-1 above the NH3 minimum. The barrier to linearity in NH2 is 11 914 cm-1 above the NH2(2B1) minimum. These results agree well with recent theoretical work.17,44,63 3. Although the quartic force field for NH3 from the present representation is significantly different from the corresponding quartic force field of the other potential energy surfaces studied here, the vibrational structure of NHD2 obtained from perturbation theory is quite similar to the corresponding spectra of the other representations. Nevertheless, we think that this representation can be significantly improved by adjustment to high-level electronic structure data without extension of the parameter set soon.40 Such global potentials will be of great value in the description of the dynamics of ammonia isotopomers including all largeamplitude motions and bond fission channels. Acknowledgment. We thank I. Thanopulos and D. Luckhaus for important discussions in the initial phase of this work. We also thank P. Rosmus, H. Hollenstein, W. Thiel, F. Mariotti, and J. Zheng for fruitful discussions. This work has received financial support from the French Ministry for Research and the Schweizerischer Nationalfonds as well as ETH Zu¨rich (including C4, AGS, and CSCS). Supporting Information Available: FORTRAN code. This material is available free of charge via the Internet at http:// pubs.acs.org. Appendix Definitions

Sb1 ) x3 sb1 + Sb2a ) x2 sb2a + Sb2b ) x2 sb2b +

4 jmax(k)

∑∑ k)2 j)1

ab(Skj1) σk(Sj 1)

(A1)

ab(Skj2) σk(Sj 2) (a)

(A2)

ab(Skj2) σk(Sj 2) (b)

(A3)

4 jmax(k)

∑∑ k)2 j)1

4 jmax(k)

∑∑ k)2 j)1

The expressions sbi are defined in eqs 23 and 24, the abkj are parameters, and the expressions σkj are given below. Number of Independent Symmetrized Polynomial Expressions S1 S2

order k

2

3

4

jmax(k) jmax(k)

2 2

3 3

4 5

List of Symmetrized Polynomial Expressions 2 1) σ(S 21 ) s b1

(A4)

2 2 1) σ(S 22 ) sb2a + sb2b

(A5)

8450 J. Phys. Chem. B, Vol. 109, No. 17, 2005 3 1) σ(S 31 ) sb1

(A6)

2 2 1) σ(S 32 ) (sb2a + sb2b)sb1

(A7)

3 2 1) σ(S 33 ) sb2a - 3sb2bsb2a

(A8)

4 1) σ(S 41 ) sb1

(A9)

2 2 2 1) σ(S 42 ) (sb2a + sb2b)sb1

(A10)

2 2 1) σ(S 43 ) (sb2a - 3sb2b)sb2asb1

(A11)

2 2 2 1) σ(S 44 ) (sb2a + sb2b)

(A12)

2) σ(S 21 )

2) σ(S 22

{

2) σ(S 32 )

)

2) σ(S 43

2) σ(S 44

)

{

2) σ(S 45

)

{

2 2) σ(S 31 (a) ) sb2asb1 2 2) σ(S 31 (b) ) sb2bsb1

{

2 2 2) σ(S 33 (a) ) (sb2a + sb2b)sb2a 2 2 2) σ(S 33 (b) ) (sb2a + sb2b)sb2b

{

{

3 2) σ(S 41 (a) ) sb2asb1 3 2) σ(S 41 (b) ) sb2bsb1

1 2 2 2 2) σ(S 42 (a) ) (sb2a - sb2b)sb1 2 ) 2 2) σ(S 42 (b) ) -sb2asb2bsb1

)

)

2) σ(S 21 (b) ) sb2bsb1

1 2 2 2) σ(S 32 (a) ) (sb2a - sb2b)sb1 2 2) σ(S 32 (b) ) -sb2asb2bsb1

2) σ(S 41 )

2) σ(S 42

{

2) σ(S 21 (a) ) sb2asb1

1 2 2 2) σ(S 22 (a) ) (sb2a - sb2b) 2 ) 2) σ(S 22 (b) ) -sb2asb2b

2) σ(S 31

2) σ(S 33

Marquardt et al.

{ {

{

2 2 2) σ(S 43 (a) ) (sb2a + sb2b)sb2asb1 2 2 2) σ(S 43 (b) ) (sb2a + sb2b)sb2bsb1

4 4 2) σ(S 44 (a) ) sb2a - sb2b 3 3 2) σ(S 44 (b) ) -2sb2asb2b - 2sb2asb2b

4 2 2 4 2) σ(S 45 (a) ) sb2a - 6sb2a sb2b + sb2b 3 3 2) σ(S 45 (b) ) 4sb2asb2b - 4sb2asb2b

(A13)

(A14)

(A15)

(A16)

(A17)

(A18)

(A19)

(A20)

(A21)

(A22)

References and Notes (1) Hewitt, S. A.; Hershberger, J. F.; Flynn, G. W.; Weston, R. E. J. Chem. Phys. 1987, 87, 1894-1895. (2) Flynn, G. W.; Weston, R. E. Annu. ReV. Phys. Chem. 1986, 37, 551-585. Weston, R. E.; Flynn, G. W. Annu. ReV. Phys. Chem. 1992, 43, 559-589. (3) Hund, F. Z. Phys. 1927, 43, 805-826. (4) Fermi, E. NuoVo Cimento 1932, 9, 177. (5) Herzberg, G. Molecular Spectra and Molecular Structure II. Infrared and Raman Spectra of Polyatomic Molecules; Van Nostrand Reinhold: New York, 1991 (reprint). (6) Longuet Higgins, H. C. Mol. Phys. 1963, 6, 445.

(7) Cheung, A. C.; Rauh, D. M.; Townes, H. C.; Thornton, D. D.; Welch, W. J. Phys. ReV. Lett. 1968, 21, 1701. (8) Takami, M.; Jones, H.; Oka, T. J. Chem. Phys. 1979, 70 (7), 35573558. (9) Sasada, H.; Hasegawa, Y.; Amano, T.; Shimuzu, T. J. Mol. Spectrosc. 1982, 96, 106-130. (10) Sˇ pirko, V. J. Mol. Spectrosc. 1983, 101, 30-47. (11) Papousˇek, D. J. Mol. Struct. 1983, 100, 179-198. (12) Tanaka, K.; Ito, H.; Tanaka, T. J. Chem. Phys. 1987, 87 (3), 15571567. (13) Sˇ pirko, V.; Kraemer, W. P. J. Mol. Spectrosc. 1989, 133, 331344. (14) Snels, M.; Fusina, L.; Hollenstein, H.; Quack, M. Mol. Phys. 2000, 98, 837-854. (15) Martin, J. M. L.; Lee, T. J.; Taylor, P. R. J. Chem. Phys. 1992, 97 (11), 8361-8371. (16) Marshall, M. D.; Izgi, K. C.; Muenter, J. S. J. Chem. Phys. 1997, 107 (4), 1037-1044. (17) Klopper, W.; Samson, C. C. M.; Tarczay, G.; Csa´sza´r, A. G. J. Comput. Chem. 2001, 22(13), 1306-1314. (18) Lin, H.; Thiel, W.; Yurchenko, S. N.; Carvajal, M.; Jensen, P. J. Chem. Phys. 2002, 117, 11265-11276. Thiel, W. Chimia 2004, 58 (5), 276-280. Yurchenko, S. N.; Carvajal, M.; Jensen, P.; Lin, H.; Zheng, J.; Thiel, W. Mol. Phys. 2005, 103, 359-378. (19) Marquardt, R.; Quack, M.; Thanopulos, I.; Luckhaus, D. J. Chem. Phys. 2003, 118 (2), 643-658. (20) Snels, M.; Hollenstein, H.; Quack, M. J. Chem. Phys. 2003, 119, 7993. (21) Leonard, C.; Handy, N. C.; Carter, S. Chem. Phys. Lett. 2003, 370, 360-365. (22) Rajama¨ki, T.; Miani, A.; Halonen, L. J. Chem. Phys. 2003, 118 (14), 6358-6369. (23) Quack, M. Angew. Chem., Int. Ed. Engl. 1989, 28, 571-586. (24) Quack, M. Angew. Chem., Int. Ed. 2002, 41, 4618-4630. (25) McCarthy, M. I.; Rosmus, P.; Werner, H.-J.; Botschwina, P.; Vaida, V. J. Chem. Phys. 1987, 86 (12), 6693-6700. (26) Biesner, J.; Schnieder, L.; Ahlers, G.; Schmeer, J.; Xie, X.; Welge, K. H.; Ashfold, M. N. R.; Dixon, R. N. J. Chem. Phys. 1988, 88 (6), 36073616. (27) Biesner, J.; Schnieder, L.; Ahlers, G.; Xie, X.; Welge, K. H.; Ashfold, M. N. R.; Dixon, R. N. J. Chem. Phys. 1989, 91 (5), 2901-2911. (28) Ma¨nz, U.; Reinsch, E.; Rosmus, P.; Werner, H.; Neil, S. O. J. Chem. Soc., Faraday Trans. 1991, 87 (12), 1809-1814. (29) R., P.; Paidarova´, I.; Sˇpirko, V.; Kuntz, P. J. Int. J. Quantum Chem. 1996, 57, 429-440. (30) Reid, J. P.; Loomis, R. A.; Leone, S. R. J. Phys. Chem. A 2000, 104, 10139-10149. (31) Bach, A.; Hutchinson, J. M.; Holiday, R. J.; Crim, F. F. J. Chem. Phys. 2002, 116 (21), 9315-9325. (32) Bach, A.; Hutchinson, J. M.; Holiday, R. J.; Crim, F. F. J. Phys. Chem. A 2003, 107, 10490-10496. (33) Rosmus, P.; Botschwina, P.; Werner, H.-J.; Vaida, V.; Engelkind, P. C.; McCarthy, M. I. J. Chem. Phys. 1987, 86 (12), 6677-6692. (34) Ma¨nz, U.; Rosmus, P.; Werner, H.-J.; Botschwina, P. Chem. Phys. 1988, 122, 387-393. (35) Chase, M. W., Jr.; Davies, C. A.; Downey, J. R., Jr.; Frurip, D. J.; McDonald, R. A.; Syverud, A. N. JANAF Thermochemical Tables, 3rd ed. J. Phys. Chem. Ref. Data 1985, 14 (Suppl. 1). (36) Berkowitz, J.; Ellison, G. B.; Gutman, D. J. Phys. Chem. 1994, 98 (11), 2744-2765. (37) Yarkony, D. R. J. Chem. Phys. 2004, 121 (2), 628-631. (38) Varandas, A. J. C. In Lecture Notes in Chemistry; Lagana`, A., Riganelli, A., Eds.; Springer: Berlin, 2000; Vol. 75, pp 33-56. (39) Seiler, R.; Hollenstein, U.; Softley, T. P.; Merkt, F. J. Chem. Phys. 2003, 118, 10024-10033. (40) Marquardt, R.; Sagui, K.; Zheng, J.; Thiel, W.; Mariotti, F.; Quack, M., to be submitted for publication. (41) Widmark, P. O.; Malmqvist, P.-A.; Roos, B. Theor. Chim. Acta 1990, 77, 291. (42) Andersson, K.; Fu¨lscher, M.; Karlstro¨m, G.; Lindh, R.; Malmqvist, P.-A.; Olsen, J.; Roos, B.; Sadlej, A.; Blomberg, M.; Siegbahn, P.; Kello¨, V.; Noga, J.; Urban, M.; Widmark, P.-O. MOLCAS, version 3; Lund University: Lund, Sweden, 1994. (43) H.-J. Werner and P. J. Knowles, version 2002.1. Amos, R. D.; Bernhardsson, A.; Berning, A.; Celani, P.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Hampel, C.; Hetzer, G.; Knowles, P. J.; Korona, T.; Lindh, R.; Lloyd, A. W.; McNicholas, S. J.; Manby, F. R.; Meyer, W.; Mura, M. E.; Nicklass, A.; Palmieri, P.; Pitzer, R.; Rauhut, G.; Schu¨tz, M.; Schumann, U.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Werner, H.-J. MOLPRO; University of Birmingham, U. K., 2002. (44) Helgaker, T.; Gauss, J.; Jørgensen, P.; Olsen, J. J. Chem. Phys. 1997, 106 (15), 6430-6440.

PES for Large-Amplitude Nuclear Motion in NH3 (45) Benedict, W. S.; Plyler, E. K. Can. J. Phys. 1957, 1235-1241. (46) Heckert, M.; Ka´llay, M.; Gauss, J. Mol. Phys., in press 2005. (47) Marquardt, R. Analytische Darstellungen Von Potentialhyperfla¨chen zur Beschreibung Von Bewegungen grosser Amplitude in XYn Moleku¨len Habilitationsschrift; Eidgeno¨ssische Technische Hochschule: Zu¨rich, Switzerland, 1997. (48) Marquardt, R.; Quack, M. J. Chem. Phys. 1998, 109 (24), 1062810643. (49) Marquardt, R.; Quack, M. J. Phys. Chem. A 2004, 108 (15), 31663181. (50) Morse, P. M. Phys. ReV. 1929, 34, 57-64. (51) Hobza, P.; R., Z. Weak Intermolecular Interactions in Chemistry and Biology; Studies in Physical and Theoretical Chemistry 3; Elsevier: New York, 1980. (52) Furue, H.; LeBlanc, J. F.; Pacey, P. D.; Whalen, J. M. Chem. Phys. 1991, 154, 425-435. (53) Jensen, P. J. Mol. Spectrosc. 1989, 133, 438-460. (54) Sˇ pirko, V.; Bunker, P. R. J. Mol. Spectrosc. 1982, 95, 381-390. (55) Quack, M.; Troe, J. Ber. Bunsenges. Phys. Chem. 1974, 78 (3), 240-252. (56) Wolf, R. J.; Bhatia, D. S.; Hase, W. L. Chem. Phys. Lett. 1986, 132 (6), 493-497. (57) Isaacson, A. D. J. Phys. Chem. 1992, 96 (2), 531-537. (58) Mladenovic´, M. J. Chem. Phys. 2003, 119 (22), 11513-11525. (59) Maple V, release 3; Waterloo Maple Software, University of Waterloo: Waterloo, ON, Canada, 1994.

J. Phys. Chem. B, Vol. 109, No. 17, 2005 8451 (60) Marquardt, R.; Sagui, K., to be submitted. (61) Marquardt, D. W. J. Soc. Ind. Appl. Math. 1963, 11(2), 431-441. (62) Martin, J. M. L. Chem. Phys. Lett. 1998, 292, 411-420. (63) Gabriel, W.; Chambaud, G.; Rosmus, P.; Carter, S.; Handy, N. C. Mol. Phys. 1994, 81 (6), 1445-1461. (64) Marquardt, R.; Quack, M.; Thanopulos, I.; Luckhaus, D. J. Chem. Phys. 2003, 119 (20), 10724-10732. (65) Luckhaus, D.; Marquardt, R.; Quack, M., to be submitted. (66) Cuvelier, F.; Herve´, S.; Marquardt, R.; Sagui, K. Chimia 2004, 58 (5), 296-305. (67) Yamaguchi, Y.; Osamura, Y.; Goddard, J. D.; Schaefer, H. F., III A New Dimension to Quantum Chemistry. Analytic DeriVatiVe Methods in Ab Initio Molecular Electronic Structure Theory; International Series of Monographs in Chemistry 29; Oxford University Press: New York, 1994. (68) Schatz, G. C. In AdVances in Molecular Electronic Structure Theory; Dunning, T. H., Jr., Ed.; JAI Press: Greenwich, CT, 1990; Vol. 1, pp 85-127. (69) Duncan, J. L.; Mills, I. M. Spectrochim. Acta 1964, 20, 523-546. (70) Kartha, S. B.; Singh, K.; Job, V. A.; Kartha, V. B. J. Mol. Spectrosc. 1988, 129, 86-98. (71) Papousek, D.; Aliev, M. R. Molecular Vibrational-Rotational Spectra; Elsevier Scientific Publishing: Amsterdam, 1982. (72) Hollenstein, H.; Marquardt, R.; Quack, M.; Suhm, M. A. J. Chem. Phys. 1994, 101 (5), 3588-3602.

Global Analytical Potential Energy Surface for Large Amplitude

Recommend Documents