New inversion coordinate for ammonia: Application to a CCSD(T) bidimensional potential energy surface

JOURNAL OF CHEMICAL PHYSICS

VOLUME 115, NUMBER 3

15 JULY 2001

New inversion coordinate for ammonia: Application to a CCSD„T… bidimensional potential energy surface Janne Pesonen and Andrea Miani Laboratory of Physical Chemistry, P.O. Box 55 (A.I. Virtasen aukio 1), FIN-00014 University of Helsinki, Finland

Lauri Halonena) JILA, University of Colorado, Boulder, Colorado 80309-0440

共Received 2 November 2000; accepted 26 April 2001兲 A new inversion coordinate is defined for ammonia as a function of the valence angles. Its square is similar to the often used totally symmetric bending displacement coordinate for the pyramidal XY 3 –type molecules. We have used this in a two-dimensional calculation including the totally symmetric stretching and the inversion mode. A conventional symmetrized internal coordinate is employed for the symmetric stretch. A two-dimensional potential energy surface is calculated using the ab initio CCSD共T兲 method together with the aug-cc-pVTZ, cc-pVQZ, and aug-cc-pVQZ basis sets. The corresponding eigenvalues are calculated variationally using a Morse oscillator basis set for the stretch and a harmonic oscillator basis set for the inversion. A good agreement is obtained between the calculated and 22 experimental inversion levels, 9 of 14NH3 and the others involving 4 other isotopomers ( 14ND3 , 15NH3 , 15ND3 , and 14NT3 ).With the aug-cc-pVTZ basis, a mean absolute error of 5.0 cm⫺1 is obtained whereas with the aug-cc-pVQZ basis set the error becomes 7.9 cm⫺1 . © 2001 American Institute of Physics. 关DOI: 10.1063/1.1379752兴 Plesset techniques.11 In the last 10 years, the coupled cluster CCSD共T兲 method 共singles, doubles, and perturbatively most important triples兲12 has provided results of good accuracy when employed to obtain force fields for semi-rigid molecules 共see, for example, Refs. 13 and 14兲. This method has been used in 1992 by Martin, Lee, and Taylor15 to determine an anharmonic force field for ammonia at the equilibrium configuration. There exist different kinds of theoretical approaches. In some of them, the main emphasis is on the inversion motion but high-frequency modes have also been included using normal coordinates. Van Vleck perturbation theory is employed to take nonresonance couplings into account. Strong resonances such as Fermi interactions between stretches and bends have also been explicitly included.9 In another type of solution, other motions explicitly excluding the inversion are treated. The effects of the inversion vibration have been taken into account by second-order Van Vleck perturbation theory in a model, which is based on the usage of curvilinear internal valence coordinates.16 In this way, it has been possible to vibrationally treat highly excited states. In a recent contribution, the vibrational Hamiltonian is expressed in terms of unsymmetrized internal coordinate type variables, including inversion,17 and eigenvalues are obtained variationally. Basically, the same coordinate system has been used in a contribution where the eigenvalues have been obtained with a new discrete variable method in a six-dimensional calculation.18 In another kind of work, the Hamiltonian is represented in a ‘‘vectorial’’ formalism with a combination of Jacobi and hyperspherical coordinates,19 and a Lanczosbased iterative diagonalization scheme20 is used to obtain the eigenvalues. In Ref. 21, the authors propose a system of symmetric hyperspherical coordinates. The inertia tensor is

I. INTRODUCTION

Ammonia is a prototype pyramidal molecule which provides a good test case for theoretical spectroscopic models. Unlike other similar molecules, phosphine, arsine, and stibine, ammonia is complicated due to a low-frequency large amplitude inversion motion, which makes it a challenge in understanding, globally, the vibrational problem. There is a large body of high-quality experimental spectroscopic data for the ground electronic state of ammonia 共see Refs. 1 and 2, and references therein兲. It extends up to the visible region, where mainly stretching vibrational overtone data are available. Ammonia is sufficiently small to allow the calculation of its electronic structure with accurate ab initio methods and large basis sets, but it is still too large to explore the whole potential energy surface 共PES兲 in detail. As a consequence, it is necessary to limit the study to regions which, for physical reasons, are significant for the problem under study. Thus, at the moment, there does not exist a six-dimensional PES, which treats accurately both the inversion and the other modes. Many of the calculations have particularly concentrated on the inversion potential 共see Refs. 3– 8兲. The value generally accepted for the inversion barrier is 1885 cm⫺1 as reported by Sˇpirko et al.,9 who fitted the inversional experimental data10 using their model Hamiltonian. This can be compared with the work by East and Radon who have recently published the value 1821⫾30 cm⫺1 with ab initio methods by investigating systematically the effect of diffuse functions on nitrogen and of core correlation in Møller– a兲

JILA Fellow 1999/2000. Present address: Laboratory of Physical Chemistry, P.O. Box 55 共A.I. Virtasen aukio 1兲, FIN-00014 University of Helsinki, Finland.

0021-9606/2001/115(3)/1243/8/$18.00

1243

© 2001 American Institute of Physics

Downloaded 21 Sep 2001 to 128.40.5.14. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

parametrized by three angular variables, which span the configurations of opposite chirality. The internal motions, which do not change the moments of inertia, are parametrized by the kinematic angles. We are presenting a new definition for the inversion coordinate in ammonia. This coordinate is unique in the sense that its square is similar to the usual definition for S 2 in the XY 3 -type molecules belonging to the C 3 v point group. We use it in a two-dimensional calculation, which includes the symmetric stretch and the inversion motion. For this purpose, we have calculated a new two-dimensional potential energy surface using the CCSD共T兲 method. We have also used the density functional B3LYP method22,23 and compared our results with those by Aquino, Campoy, and Yee-Madeira.24

II. INVERSION COORDINATE

This paper deals with symmetric vibrational modes in ammonia, i.e., the symmetric stretch and the inversion. Ammonia belongs to the C 3 v point group at the equilibrium geometry, but we are also allowing the planar configuration due to the inversion motion. Thus, the point group D 3h or more precisely the isomorphic molecular symmetry group is used.25 The symmetric stretch is defined in the standard way26 S 1 共 A 1⬘ 兲 ⫽

1

共 ⌬r 1 ⫹⌬r 2 ⫹⌬r 3 兲 ,

冑3

共1兲

where ⌬r i ⫽r i ⫺R e (D 3h ), r i is the length of the bond ri ⫽xHi ⫺xN , vectors x␣ give positions of the nuclei relative to some chosen origin, R e (D 3h ) is the equilibrium value for the bond length in the planar D 3h reference configuration. The other symmetric displacement coordinates S 3a (E ⬘ ), S 3b (E ⬘ ), S 4a (E ⬘ ), and S 4b (E ⬘ ), if needed, are defined as linear combinations of the bond stretching and valence angle displacements in the usual way.26 An ideal inversion coordinate S 2 (A 2⬙ ) should fulfill some strict requirements. First of all, it must be a continuous monotonic function. It should be single valued with respect to the totally symmetric bend 共i.e., for each value of equal valence angles ␪ 12⫽ ␪ 13⫽ ␪ 23 there must be a unique value of S 2 ). Second, it must possess the correct symmetry (A 2⬙ ) and it should be easily related to the totally symmetric (A 1 ) bending coordinate for ammonia when the tunneling is not allowed for 共i.e., when ammonia is considered to belong to the C 3 v point group兲. Third, it should not include stretching character, that is, its value should be the same when the valence angles are constant and the bond lengths change. Finally, in order to have the possibility to derive an exact expression for the kinetic energy operator, the inverse coordinate relations should be easily found in closed form. We choose S 2 as S 2 共 A 2⬙ 兲 ⫽⫾ ⫽⫾

1 3

1/4

1 3 1/4

共 2 ␲ ⫺ ␪ 23⫺ ␪ 13⫺ ␪ 12兲 1/2 关 ⫺ 共 ⌬ ␪ 23⫹⌬ ␪ 13⫹⌬ ␪ 12兲兴 1/2,

共2兲

where the ⫹ sign is for the right-handed and the ⫺ sign for the left-handed configuration. The displacements ⌬ ␪ i j are measured from ␪ e (D 3h )⫽ 2 ␲ /3. We have inserted a proportionality factor 1/31/4

J. Chem. Phys., Vol. 115, No. 3, 15 July 2001

New inversion coordinate for ammonia

is the case, if one defines S 2 as the sum of the angles ␾ i between the bond ri and the plane r j ⵩rk , 30 i.e., if S 2 共 A 2⬙ 兲 ⫽

1

冑3

共 ⌬ ␾ 1 ⫹⌬ ␾ 2 ⫹⌬ ␾ 3 兲 ,

共7兲

because there exist three nonlinear redundancy conditions involving ␾ i . The same problem is encountered, if S 2 is defined as the sum of the cosines of the angles 肀 i between the bond ri and the axis r1 ⫻r2 ⫹r2 ⫻r3 ⫹r3 ⫻r1 n⫽ 兩 r1 ⫻r2 ⫹r2 ⫻r3 ⫹r3 ⫻r1 兩 perpendicular to the hydrogen plane, i.e., S 2 共 A 2⬙ 兲 ⫽

1

冑3

共 cos 肀 1 ⫹cos 肀 2 ⫹cos 肀 3 兲 ,

共8兲

III. TWO-DIMENSIONAL VIBRATIONAL ENERGY LEVEL CALCULATION A. Kinetic energy operator

By using some suitable internal scalar displacement coordinates ⌬q i ⫽q i ⫺q (e) i , the expectation value of the vibrational kinetic energy,31–33 ⳵ ប2 1 ⳵J d ␶ ⌿ * J 1/2w ⫺1/2 ⫹ 具 Tˆ 典 ⫽⫺ 2 ⳵ ⌬q i J ⳵ ⌬q i i, j

冕

J⫽

⫽

冏冏

冉

冊

⳵ J ⫺1/2w 1/2⌿, 共9兲 ⳵ ⌬q j where ⌿ is the eigenfunction of the vibrational Hamiltonian in question, the volume element of the integration is d ␶ ⫽wd⌬q 1 d⌬q 2 . . . , J⫽

兩 det g (q i q j ) 兩 ⫺1/2 兿 ␣ m ␣3/2

兩 sin ␤ 兩

共10兲

is the absolute value of the Jacobian of the coordinate transformation including the rotational part ( ␤ is an Euler angle between the molecule and space-fixed z axes兲, and N

g (q i q j ) ⫽

1

兺␣ m ␣ 共 ⵜ ␣ q i 兲 • 共 ⵜ ␣ q j 兲

共11兲

is an element of the mass weighted reciprocal metric tensor. The quantity ␣ is an index over the nuclei, m ␣ is the mass of the ␣ th nucleus, and N is the number of the nuclei. We use symmetry coordinates but this is not a problem because it is possible to express their gradients in terms of the already known internal coordinate gradients (ⵜ Nr i , ⵜ N␪ i j , etc.兲.29,30 The gradients for all internal and symmetric displacement coordinates are tabulated as supplementary material34 and they can be expressed as the functions of the symmetry coordinates S r by using inverse relations between the internal and the symmetry coordinates. Thus, the g–tensor elements g (S r S s ) are formed from Eq. 共11兲 by replacing q i coordinates by S r coordinates. By using the chain rule,35 the absolute value of the new Jacobian is36

冏

⳵ 共 x H1 ,y H1 ,z H1 , . . . ,x N ,y N ,z N兲 ⳵ 共 r 1 ,r 2 ,r 3 , ␪ 12 , ␪ 13 , ␪ 23 , ␤ 兲 ⳵ 共 r 1 ,r 2 ,r 3 , ␪ 12 , ␪ 13 , ␪ 23 , ␤ 兲 ⳵ 共 S 1 ,S 2 ,S 3a ,S 3b ,S 4a ,S 4b 兲

冏

r 21 r 22 r 23 sin ␪ 12 sin ␪ 13 sin ␪ 23 sin ␤

冏

冏

共 1⫺cos2 ␪ 12⫺cos2 ␪ 13⫺cos2 ␪ 23⫹2 cos ␪ 12 cos ␪ 13 cos ␪ 23兲 1/2

The factor sin ␤ can be neglected because it can be treated as a multiplicative constant in the pure vibrational problem. We obtain our final two-dimensional kinetic energy operator by rearranging terms and by using the volume element

兺

⫻g (q i q j )

where ri •n ri •r j ⫻rk ⫽ cos 肀 i ⫽ ri r i 兩 r1 ⫻r2 ⫹r2 ⫻r3 ⫹r3 ⫻r1 兩 共indices i, j, and k are in cyclic order兲. In a recent contribution, the inversion problem is formulated in a somewhat different manner.17 The inversion coordinate is defined as an angle which each bond makes with the trisector of them. The other two angle variables are two valence angles projected to a plane perpendicular to the bisector plane. In this approach, the kinetic energy operator is obtained in closed form. A disadvantage, when compared with our approach, might be the possibility that the bending coordinates do not necessarily resemble much of the physical intuition of the form of coordinates which describe in the best way the modes involved. This comment also applies to the work in Ref. 20. In any case, these approaches have been successful in modeling the high frequency modes. Another point worth mentioning is the possibility that when compared with previous work,17,20 the potential energy function in our coordinate system might provide the best and physically the most meaningful representation in the regions of interest, i.e., both in the inversion and stretching vibrational region.

1245

⫻ 兩 ⫺2S 2 兩 .

共12兲

d ␶ ⫽dS 1 dS 2 as

Tˆ ⫽⫺

ប2 2

2

兺

r,s⫽1

冋

(0) (1) f rs ⫹ f rs

册

⳵ ⳵2 (2) ⫹ f rs , ⳵Ss ⳵ S r⳵ S s

共13兲

Downloaded 21 Sep 2001 to 128.40.5.14. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

1246

J. Chem. Phys., Vol. 115, No. 3, 15 July 2001

Pesonen, Miani, and Halonen

TABLE I. Ammonia geometry and inversion energy calculated using many basis sets and the CCSD共T兲 method and the B3LYP density functional.a Number basis

r eq 共Å兲

␪ eq 共Deg.兲

(eq) 共au兲 E tot

r st 共Å兲

(st) 共au兲 E tot

⌬ inv 共cm⫺1 )

20 72 115 29 50 72 115 145 218

1.018 868 1.013 837 1.013 049 1.027 316 1.023 670 1.014 098 1.014 950 1.012 407 1.012 769 1.011 6b

106.6311 106.4766 107.2293 103.5536 105.9298 105.6518 106.3911 106.1759 106.5304 106.68b

⫺56.564 040 5 ⫺56.584 725 3 ⫺56.588 857 5 ⫺56.402 802 2 ⫺56.425 519 9 ⫺56.473 197 3 ⫺56.480 562 6 ⫺56.493 053 0 ⫺56.495 732 4

1.004 190 0.996 617 0.997 438 1.005 140 1.005 433 0.995 127 0.997 501 0.994 873 0.995 978

⫺56.554 867 4 ⫺56.576 286 4 ⫺56.581 536 7 ⫺56.388 657 2 ⫺56.416 255 ⫺56.463 000 8 ⫺56.471 739 2 ⫺56.483 823 6 ⫺56.487 226 1

2103.3 1852.1 1606.7 3104.5 2033.4 2237.9 1936.5 2025.6 1866.9 1884.7c

B3LYP/DZVP B3LYP/cc-pVTZ B3LYP/aug-cc-pVTZ CCSD共T兲/cc-pVDZ CCSD共T兲/aug-cc-pVDZ CCSD共T兲/cc-pVTZ CCSD共T兲/aug-cc-pVTZ CCSD共T兲/cc-pVQZ CCSD共T兲/aug-cc-pVQZ Expt.

Number basis is the number of Gaussian basis functions in the ab initio calculation; r eq and ␪ eq are the bond length and valence angle for ammonia in its (eq) (st) equilibrium pyramidal configuration; r st is the bond length for ammonia in its planar configuration; E tot and E tot are the total energy 共electronic plus nuclear repulsion兲 for ammonia in the equilibrium and planar configurations, respectively. ⌬ inv is the difference between the energy of the minimum of the PES and that of the saddle point. b Taken from Ref. 37. c Taken from Ref. 9.

a

B. Potential energy surface

where

⳵ J ⳵ J 1 ⫺1 ⳵ g (S r S s ) ⳵ J 1 (0) f rs ⫽ J ⫺2 g (S r S s ) ⫺ J 4 ⳵Sr ⳵Ss 2 ⳵Sr ⳵Ss ⳵ 2J 1 ⫺ J ⫺1 g (S r S s ) , 2 ⳵ S r⳵ S s (1) ⫽ f rs

⳵ g (S r S s ) , ⳵Sr

共14兲

(2) ⫽g (S r S s ) . f rs

The Jacobian for the two-dimensional problem is given by sin3 ␪

J⫽2r 6 兩 S 2 兩

,

冑1⫺3 cos2 ␪ ⫹2 cos3 ␪ ␪ ⫽ 2 ␲ /3 ⫺ S 22 / 冑3 and r⫽R e ⫹ S 1 / 冑3.

where cal metric tensor elements are g (S 1 S 1 ) ⫽ g (S 2 S 2 ) ⫽

The recipro-

1 2 1 ⫹ ⫹ cos ␪ , mN mH mN 1 r 2 S 22 ⫹

g (S 1 S 2 ) ⫽

共15兲

冉

1 1⫹2 cos ␪ 2m H 1⫹cos ␪

冊

1 1⫹cos ␪ ⫺2 cos2 ␪ , mN 1⫹cos ␪

共16兲

1 cos ␪ ⫺cos 2 ␪ . m NrS 2 sin ␪

The first and the second derivatives of the Jacobian with respect to S 1 and S 2 , and the first derivatives of the reciprocal metric tensor elements are given as supplementary material.34 By applying l’Hospital’s rule, it is seen that there are no singularities in the equations expressing the Jacobian, the g–tensor matrix elements and their derivatives, even though many of these terms possess the indeterminate form 0/0. This can create some numerical instabilities that will be dealt with later.

Before calculating the bidimensional PES, we tested some basis sets to be used both with the CCSD共T兲 ab initio method and the B3LYP density functional approach. The results together with a comparison with experimental data37 are reported in Table I. As it can be seen, there is a large improvement in the calculation of the inversion barrier when diffuse functions are added to the basis set. The basis sets, which were selected for the calculations, are the aug-ccpVTZ, cc-pVQZ, and aug-cc-pVQZ bases of Dunning and coworkers.38,39 The acronyms cc-pVXZ indicate correlation consistent polarized valence X zeta basis set, where X reads T 共triple兲 or Q 共quadruple兲 in our case. The prefix ‘‘aug’’ means augmented and refers to the presence of diffuse functions necessary to describe accurately molecular anions, Rydberg states or, as in the present case, the lone pair of electrons present in ammonia. The time required for the computation of the energy at one point of the PES is about 15 min with the aug-cc-pVTZ basis set, 1 h with cc-pVQZ, and about 8 h with the largest basis set we tried, aug-cc-pVQZ. All calculations have been performed on an SGI ORIGIN 2000 computer at the CSC computer center in Helsinki using the GAUSSIAN 98 program.40 With the smallest basis set 共aug-cc-pVTZ兲, we calculated a grid of 630 points 共45 for the S 1 coordinate and 14 for S 2 ). With the cc-pVQZ basis, we calculated 208 single points 共16 ⫻13兲, whereas 162 points have been calculated for the augcc-pVQZ basis set 共18⫻9兲. The most important regions of the PES are the two minima and the region around the saddle point, but because we have to calculate integrals involving highly exited states, it is also important to have a good description of the PES quite far from these regions. We could determine how large the grid should have been after some preliminary tests were made with the smallest basis sets. We decided to calculate a grid within the following boundaries: 0.6 Å⭐r⭐2.5 Å,

共17兲

70°⭐ ␪ ⭐120°.

共18兲

Downloaded 21 Sep 2001 to 128.40.5.14. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

The grid of ab initio points obtained with all the methods employed in this work was chosen to be more dense in the region close to the minima of the PES and to the saddle points. To fit the ab initio points, we used a series expansion in orthogonal Chebyshev polynomials41 T n (x) of degree n, which span the domain 关 ⫺1,1兴 . To employ a Chebyshev polynomial between 关 S min ,Smax兴, it is sufficient to use the following equation ¯S ⫽

2S⫺ 共 S min⫹S max兲 , S max⫺S min

共19兲

where S is S 1 or S 2 . We expressed our PES as the direct product of monodimensional Chebyshev polynomials, that is as V⫽

兺i j a i j T i共 S 1 兲 T j 共 S 2 兲 ,

共20兲

where S 1 and S 2 were calculated using Eq. 共19兲. Chebyshev polynomials are widely used in interpolation problems because they, among all polynomial approximations of the same degree, are similar to the minimax polynomial, that is a polynomial function which possesses the smallest maximum deviation from the true function f (x) 共which the Chebyshev series approximates兲.42 The a i j coefficients of Eq. 共20兲 were calculated using the NAG subroutine E02CAF43 and are available as supplementary material34 for the CCSD共T兲/aug-cc-pVTZ, CCSD共T兲/ccpVQZ, and CCSD共T兲/aug-cc-pVQZ basis sets. The analytical expressions obtained after fitting the points show some oscillations, but these are observed in regions which are sufficiently far away from the minima and the saddle point, not to influence the lowest eigenvalues obtained in our calculations. Finally, we have also tested the significance of oscillations between calculated and fitted points using splines. They make the calculation of integrals involving the PES more complicated but they have the advantage that the ab initio points can be fitted exactly without oscillations. We found a difference of less than 0.1 cm⫺1 on average when we compared the vibrational levels obtained with splines and with the Chebyshev polynomial fit. C. Calculation of the energy levels

In all calculations in this work, as basis sets we have used a set of Morse oscillator eigenfunctions for the symmetric stretch and a set of harmonic oscillator eigenfunctions for the inversion. The definition of the Morse functions can be obtained from Refs. 44 – 46. The analytic expressions for the integrals of the first and the second derivative operators, which are used in this work, can be found either in Ref. 46 or 47. The analytical matrix elements for the first- and the second-order derivatives of the harmonic oscillator functions can be found in books of molecular spectroscopy 共e.g., see Ref. 30兲. There are three parameters present in the basis functions: The Morse dissociation energy (D 1 ) and the Morse steepness parameter ( ␣ 1 ), and the harmonic oscillator scaling parameter ␣ 2 ⫽2 ␲ c 0 ␮ ␻ /ប, where ␻ is the harmonic wave number and ␮ is the reduced mass of the appropriate har-

monic oscillator. They do not appear in the Hamiltonian, and their presence is useful for two different purposes. If optimized, they allow to increase the convergence rate in the variational calculation and they make it possible to check the correctness of the program, considering that the eigenvalues must not depend on their values. In this bi-dimensional calculation, we did not take the first aspect fully into account 共i.e., we did not try to find fully optimized parameters兲, considering the relatively small matrices to be diagonalized. We found approximately an optimal set of parameters from preliminary monodimensional calculations. In fitting a Morse potential to an ab initio aug-cc-pVTZ monodimensional PES at a fixed pyramidal configuration, we obtained the values for D 1 and ␣ 1 to be 112 301.16 cm⫺1 and 1.250 95 Å ⫺1 , respectively. The value for ␣ 2 was fixed at 36.0 from considerations on the behavior of the monodimensional basis set involving S 2 . The value adopted corresponds to a harmonic surface which roughly coincides with the outer part of the inversion surface. A much larger value for ␣ 2 would contract the basis functions too much and make the convergence of the eigenvalues slower, but a smaller value would let the wave functions for high quantum numbers not vanish before S 2 ⫽ 冑2 ␲ / 冑3 共or ␪ ⫽0兲toaa4 Tf91-[(tnsis)-60339

1248

J. Chem. Phys., Vol. 115, No. 3, 15 July 2001

Pesonen, Miani, and Halonen

TABLE II. Calculated and experimental inversion levels for ammonia. GS is the ground state. All values are expressed in cm⫺1 . This work States

B3LYP24

CISD5

MP24

CEPA8

CBS-QCI3

6D18

GS⫹ GS⫺ ␯⫹ 2 ␯⫺ 2 2␯⫹ 2 2␯⫺ 2 3␯⫹ 2 3␯⫺ 2 4␯⫹ 2 4␯⫺ 2

0 0.84 931.72 968.67 1596.76 1885.33 2389.14 2902.99

0 0.5 1031.2 1071.4 1627.4 1974.9

0 1.02 986.7 1032.3 1682.0 2011.6

0 0.95 958.0 998.7 1640.7 1945.6 2466.0 2998.1 3585.6 4205.9

0 0.99 905.5 946.84 1553.20 1852.91 2353.99 2869.27 3439.37 4046.53

0 0.22 1018.33 1030.59 1805.46 1975.92 2500.69 2958.38 3504.40 4079.20

1D/augTZa 2D/B3LYPb 2D/augTZc 2D/QZd 2D/augQZe 0 0.62 969.62 1002.26 1664.71 1947.32 2467.66 2995.66 3584.86 4212.21

0 0.74 959.72 994.39 1650.00 1941.24 2455.60 2987.68 3573.99 4196.56

0 0.79 936.01 971.23 1606.52 1889.58 2394.76 2907.65 3477.95 4080.17

0 0.62 960.66 990.37 1652.58 1918.60 2427.48 2937.04 3507.42 4110.81

Experimentf

0 0.96 922.92 964.74 1577.97 1882.32 2387.96 2909.76 3485.55 4093.93

0 0.793 932.43 968.12 1598.47 1882.18 2384.17 2895.61

a

Monodimensional CCSD共T兲/aug-cc-pVTZ PES calculated at r⫽1.014 95 Å. Bidimensional PES calculated with the B3LYP/DZVP basis set. c Bidimensional PES calculated with the CCSD共T兲/aug-cc-pVTZ basis set. d Bidimensional PES calculated with the CCSD共T兲/cc-pVQZ basis set. e Bidimensional PES calculated with the CCSD共T兲/aug-cc-pVQZ basis set. f Experimental data taken from Ref. 10. b

dimensional results are clearly not good for inversion states because the surface has been constructed near the equilibrium and, consequently, it does not provide an accurate representation near the saddle point. It can also be seen that the coupling between the ␯ 1 (A 1⬘ ) and the ␯ 2 (A 2⬙ ) modes causes a relaxation of all the inversion levels of several cm⫺1 , as it is logical to expect because the PES we used in one-dimensional calculations was not relaxed, but it was a section of the bidimensional PES corresponding to the pyramidal configuration, giving a barrier of about 2110 cm⫺1 . From Table I, it is evident that, despite the total electronic and repulsive energy decreases for both the minima and the saddle point, the variation is not the same for the two states, causing the inversion barrier to change irregularly. Going from the aug-cc-pVTZ basis set to the cc-pVQZ and aug-cc-pVQZ basis sets, the diffuse functions in the triple basis set allow to better describe the inversion barrier but not the geometry which is better described by the quadruple basis sets. The values of both the pyramidal geometry and of the inversion barrier are good and compare well to what are considered to be the best experimental values so far. With the aug-cc-pVTZ basis set, the bond length is a little

basis sets. In Table III, we publish the inversion levels for the isotopomers 14ND3 , 15NH3 , 15ND3 , and 14NT3 obtained with the CCSD共T兲/aug-cc-pVTZ and CCSD共T兲/aug-ccpVQZ ab initio methods and the experimental values obtained by Sˇpirko et al.10 Finally, in Table IV, we give the calculated and experimental splittings between the antisymmetric and symmetric inversion levels for 14NH3 and the other isotopomers. From the results in Table II, it is possible to see that the inversion levels for ammonia are reproduced much better than previously reported. Probably, the best ab initio PES so far that of Rush and Wiberg3 underestimates all vibrational levels about 30 cm⫺1 . We also give in Table II the levels calculated using the DZVP basis set50 together with the B3LYP density functional to allow a comparison with the result of Aquino et al.,24 who obtained a remarkably good agreement with their monodimensional potential. Our results indicate that their numbers must be considered fortuitous. Inversion states from a six-dimensional calculation 18 using the potential energy surface of Martin et al.15 are also given in Table II. The other six-dimensional calculations17,20 give essential identical results with the same surface because these models use exact kinetic energy operators. The six-

TABLE III. Experimental and calculated inversion levels for isotopomers of ammonia. All values are expressed in cm⫺1 . 14

15

ND3

a

15

NH3

b

c

a

States

augTZ

augQZ

Expt.

augTZ

GS⫹ GS⫺ ␯⫹ 2 ␯⫺ 2 2␯⫹ 2 2␯⫺ 2 3␯⫹ 2 3␯⫺ 2

0.00 0.05 746.85 750.34 1364.20 1433.73 1834.39 2113.17

0.00 0.07 740.41 744.89 1339.62 1422.56 1810.29 2104.66

0.00 0.05 745.60 749.15 1359.0 1429.0 1830.0 2106.6

0.00 0.74 932.69 966.38 1601.64 1878.67 2380.30 2888.18

14

ND3

b

c

a

NT3

b

c

a

augQZ

Expt.

augTZ

augQZ

Expt.

augTZ

augQZb

Expt.c

0.00 0.91 919.81 959.89 1572.91 1871.24 2373.15 2890.01

0.00 0.76 928.46 962.89 1591.19 1870.86 2369.32 2876.13

0.00 0.05 741.22 744.36 1357.44 1422.07 1822.35 2092.95

0.00 0.06 734.93 738.97 1333.36 1410.85 1797.55 2084.10

0.00 0.05 739.53 742.78

0.00 0.01 656.80 657.59 1236.94 1258.91 1663.97 1826.74

0.00 0.01 652.00 653.06 1220.03 1248.16 1633.18 1814.64

0.00 0.01 656.37 657.19

a

Basis set aug-cc-pVTZ. Basis set aug-cc-pVQZ. c Experimental data taken from Ref. 10. b

Downloaded 21 Sep 2001 to 128.40.5.14. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 115, No. 3, 15 July 2001

New inversion coordinate for ammonia

1249

TABLE IV. Experimental and calculated inversion splittings for ammonia and its isotopomers. All values are expressed in cm⫺1 . 14

14

NH3

States

a

augTZ

GS ␯2 2␯2 3␯2

0.79 35.22 283.06 512.89

15

ND3

b

augQZ

c

Expt.

augTZ

a

0.96 41.82 304.35 521.8

0.79 35.69 283.71 511.44

0.05 3.49 69.53 278.70

15

NH3

b

augQZ

c

Expt.

a

augTZ

0.07 4.48 82.94 294.37

0.05 3.55 70.0 276.6

0.74 33.69 277.03 507.88

14

ND3

b

NT3

augQZ

c

Expt.

a

augTZ

b

augQZ

0.91 40.08 298.33 516.86

0.76 34.43 279.67 506.82

0.05 3.14 64.63 270.60

0.06 4.04 77.49 286.55

c

Expt.

a

augTZ

augQZb

Expt.c

0.05 3.25

0.01 0.79 21.97 162.77

0.01 1.06 28.13 181.46

0.01 0.82

a

Basis set aug-cc-pVTZ. Basis set aug-cc-pVQZ. c Experimental data taken from Ref. 10. b

overestimated and the valence angle and the inversional barrier are well-reproduced by the calculation. From Table II, it is clear that the above-mentioned nonmonotonic behavior with respect to the basis set of the structural properties for ammonia is also present in the calculated spectrum. With the aug-cc-pVTZ and cc-pVQZ basis sets, all of the inversional levels are calculated to be above the experimental values, behavior which is confirmed for all the isotopomers 共see Table III兲. The aug-cc-pVTZ basis reproduces better the observed spectrum, as a consequence of the lower calculated inversion barrier giving a mean absolute error of 5.0 cm⫺1 on all the 22 vibrational levels of the studied ammonia isotopomers 共excluding the precisely reproduced splittings of the ground state for the various species兲. The calculation made with the aug-cc-pVQZ basis set also reproduces well the experimental levels 共mean absolute error of 7.9 cm⫺1 , calculated as before兲 even though their behavior is less regular than that observed with the aug-cc-pVTZ basis set. The lowest levels are below the experimental values, while the higher ones are above both the experimental and of the aug-cc-pVTZ basis set values. This behavior can be understood from Table IV, where the calculated and the experimental splittings between the symmetric and the antisymmetric states for ammonia and its isotopomers are reported. The splittings should be less influenced by the coupling with the other modes than the levels themselves, assuming that the harmonic and anharmonic coupling, which here is neglected, is similar for both states. They should reflect more the influence of the inversion PES and, as a consequence, should compare better to the experimental data. In Table IV, the agreement between the aug-cc-pVTZ results and the experimental data is good for all studied ammonia isotopomers. It can be seen that, as far as the aug-cc-pVQZ basis set is concerned, all the calculated splittings are larger than the experimental ones. It seems reasonable to deduce from this result that the aug-cc-pVQZ calculation slightly underestimates the inversion barrier for ammonia, which could be closer to the value predicted by the aug-cc-pVTZ basis set 共1936.5 cm⫺1 ) than to that calculated with the augcc-pVQZ basis 共1866.9 cm⫺1 ). A definitive conclusion on this issue will be possible after a six-dimensional calculation will be completed at least in the case of the aug-cc-pVTZ basis set. It might well be that the basis set aug-cc-pVQZ is too expensive to be used at the present moment. We may also compare our results for one- and twodimensional calculations 共see Table II兲. It is clear that the one-dimensional inversion model is not as good in explain-

ing experimental observations. One problem is the value of the NH bond length, which is needed for the reduced mass of the one-dimensional model and, which according to our electronic structure calculations, change a bit during inversion 共see Table I兲. This change is probably the physical origin of the coupling between the symmetric modes and is the reason for the success of the two-dimensional Hamiltonian. Further support of the two-dimensional approach is the inversion manifold of the symmetric stretching vibration ␯ 1 in 14NH3 . Due to neglected anharmonic interactions between ␯ 1 and the other high-frequency modes 共both nonresonance and resonance interactions, such as Fermi resonances between ␯ 1 and 2␯ 4 ), we are unable to calculate accurately the position of ␯ 1 or ␯ 1 ⫹ ␯ 2 共our calculated numbers are about 80 cm⫺1 too high when compared with experimental results兲 but the inversion splittings are well reproduced with the CCSD共T兲/aug-cc-pVTZ method: the observed splitting9 0.99 cm⫺1 in ␯ 1 and 26.42 cm⫺1 in ␯ 1 ⫹ ␯ 2 should be compared with our calculated values 1.57 cm⫺1 and 23.6 cm⫺1 , respectively. This is pleasing in the light of the fact that the inversion states in ␯ 1 are less isolated than the corresponding states in the ground vibrational state. The six-dimensional calculations17,18,20 with the potential energy surface of Martin et al.,15 do not produce these splittings any better: 0.74 cm⫺1 in ␯ 1 and 10.15 cm⫺1 in ␯ 1 ⫹ ␯ 2 . In fact, our results are clearly better for ␯ 1 ⫹ ␯ 2 . This is probably an indication of the superiority of our surface in describing inversion. The absolute values for the ␯ 1 and ␯ 1 ⫹ ␯ 2 states, about 30 and 110 cm⫺1 too high, are not produced well either. This again reflects the quality of the potential energy surface used in the six-dimensional calculations. It seems that a new six-dimensional surface is needed that describes more accurately the inversion part of the space. Our experience with the two-dimensional approach should help in this. Finally, we briefly compare our kinetic energy operator with others, particularly those presented in Refs. 17 and 20. Neither the Hamiltonian nor the coordinates used here contain singularities in important regions, for example, in the symmetrical planar configuration. Models based on different coordinate systems give different insight to the problem under study. The coordinate system we present, curvilinear internal symmetry coordinates, together with our new inversion coordinate, might be a more natural choice when these models are extended to larger molecular systems.

Downloaded 21 Sep 2001 to 128.40.5.14. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

1250

J. Chem. Phys., Vol. 115, No. 3, 15 July 2001

V. CONCLUSION

We have presented a new symmetry coordinate for ammonia defined as a function of the curvilinear internal valence angles. We have obtained the gradients of the symmetry coordinates with which it is easy to obtain the sixdimensional g metric tensor elements if needed. We have tested our kinetic energy operator with a preliminary bidimensional potential energy surface that has been calculated using highly accurate ab initio methods and basis sets. From the variational calculations, we could obtain the best agreement obtained so far between ab initio and experimental data. A noticeable outcome of our results is the importance of the coupling of the symmetric stretch with the inversion motion in calculating the inversion energy levels. This finding might be of importance in modelling chemical reactions in a way that the reaction mode 共large amplitude motion corresponding to inversion in our problem兲 is assumed to be decoupled from the vibrational degrees of freedom. To proceed further in modelling the vibrational motion in ammonia, it is necessary to perform full six-dimensional calculations. This can be done with the kinetic energy operator presented in this paper. A new full six-dimensional PES is necessary, and we believe that the bidimensional PES presented here is a good starting point on which an accurate six-dimensional PES can be built. The results from this kind of analysis will show how good our present two-dimensional approach is. ACKNOWLEDGMENTS

The authors wish to thank both The Academy of Finland and The European Commission 共membership agreement HPRN-CT-1999-00005兲 for financial support. The authors thank CSC Scientific Computing for computer time. 1

K. K. Lehmann and S. L. Coy, J. Chem. Soc., Faraday Trans. 2 84, 1389 共1988兲. 2 S. L. Coy and K. K. Lehmann, Spectrochim. Acta, Part A 45, 47 共1989兲. 3 D. Rush and K. Wiberg, J. Phys. Chem. A 101, 3143 共1997兲. 4 P. Wormer, E. Olthof, R. Engeln, and J. Reuss, Chem. Phys. 178, 189 共1993兲. 5 P. Bunker, W. Kraemer, and V. Sˇpirko, Can. J. Phys. 62, 1801 共1984兲. 6 G. Campoy, A. Palma, and L. Sandoval, Int. J. Quantum Chem., Quantum Chem. Symp. 23, 355 共1989兲. 7 A. M. Halpern and E. D. Glendening, Chem. Phys. Lett. 333, 391 共2001兲. 8 P. Rosmus, P. Botschwina, H.-J. Werner, V. Vaida, P. C. Engelking, and M. I. McCarthy, J. Chem. Phys. 86, 6677 共1987兲. 9 V. Sˇpirko and W. P. Kraemer, J. Mol. Spectrosc. 133, 331 共1989兲. 10 V. Sˇpirko, J. Mol. Spectrosc. 101, 30 共1983兲. 11 A. East and L. Radom, J. Mol. Struct. 376, 437 共1996兲. 12 K. Raghavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon, Chem. Phys. Lett. 157, 479 共1989兲. 13 J. M. L. Martin, T. J. Lee, P. R. Taylor, and J. Franc¸ois, J. Chem. Phys. 103, 2589 共1996兲.

Pesonen, Miani, and Halonen 14

T. J. Lee, J. M. L. Martin, and P. R. Taylor, J. Chem. Phys. 102, 254 共1995兲. 15 J. Martin, T. Lee, and P. Taylor, J. Chem. Phys. 97, 8361 共1992兲. 16 E. Kauppi and L. Halonen, J. Chem. Phys. 103, 6861 共1995兲. 17 N. C. Handy, S. Carter, and S. M. Colwell, Mol. Phys. 96, 477 共1999兲. 18 D. Luckhaus, J. Chem. Phys. 113, 1329 共2000兲. 19 F. Gatti, J. Chem. Phys. 111, 7225 共1999兲. 20 F. Gatti, C. Iung, and C. Leforestier, J. Chem. Phys. 111, 7236 共1999兲. 21 V. Aquilanti and S. Cavalli, J. Chem. Soc., Faraday Trans. 93, 801 共1997兲. 22 C. Lee, W. Yang, and R. Parr, Phys. Rev. B 37, 785 共1988兲. 23 A. Becke, J. Chem. Phys. 98, 5648 共1993兲. 24 N. Aquino, G. Campoy, and H. Yee-Madeira, Chem. Phys. Lett. 296, 111 共1998兲. 25 P. R. Bunker and P. Jensen, Molecular Symmetry and Spectroscopy, 2nd ed. 共NRC Research, Ottawa, 1998兲. 26 F. A. Cotton, Chemical Applications of Group Theory, 3rd ed. 共Wiley, New York, 1990兲. 27 A. R. Hoy, G. Strey, and I. M. Mills, Mol. Phys. 24, 1265 共1972兲. 28 Y. Pak, E. L. Sibert, and R. C. Woods, J. Chem. Phys. 107, 1717 共1997兲. 29 J. Pesonen, J. Chem. Phys. 112, 3121 共2000兲. 30 E. B. Wilson, J. C. Decius, and P. C. Cross, Molecular Vibrations 共McGraw-Hill, New York, 1955兲. 31 B. Podolsky, Phys. Rev. 32, 812 共1928兲. 32 G. Gruber, Int. J. Theor. Phys. 6, 31 共1972兲. 33 A. Nauts and X. Chapuisat, Mol. Phys. 55, 1287 共1985兲. 34 See EPAPS Document No. E-JCPSA6-115-008127 for latex files containing the tables with internal coordinate gradients, first and second derivatives of the bidimensional Jacobian and first derivatives of the bidimensional g matrix elements, and the coefficients of the Chebyshev polynomials obtained fitting the CCSD共T兲/aug-cc-pVTZ, CCSD共T兲/ccpVQZ and CCSD共T兲/aug-cc-pVQZ ab initio points. This document may be retrieved via the EPAPS homepage 共http://www.aip.org/pubservs/ epaps.html兲 or from ftp.aip.org in the directory /epaps/. See the EPAPS homepage for more information. 35 H. Margenau and G. M. Murphy, The Mathematics of Physics and Chemistry 共Van Nostrand, Princeton, 1956兲. 36 T. Lukka, J. Chem. Phys. 102, 3945 共1995兲. 37 Landolt-Bo¨rnstein, Numerical Data and Functional Relationships in Science and Technology II共7,2兲 共Springer, New York, 1976兲. 38 T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 共1989兲. 39 R. A. Kendall, T. H. Dunning, Jr., and R. J. Harrison, J. Chem. Phys. 96, 6796 共1992兲. 40 M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., GAUSSIAN 98, Revision A.7, Gaussian, Inc., Pittsburgh, PA, 1998. 41 S. Wolfram, The Mathematica Book, 4th ed. 共Wolfram Media/Cambridge University Press, Cambridge, MA, 1999兲. 42 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran77: The Art of Scientific Computing 共Cambridge University Press, Cambridge, MA, 1992兲. 43 See http://www.nag.com/numeric/flolch/mk17/e/e02caf.html for details on the numerical method. 44 P. M. Morse, Phys. Rev. 34, 57 共1929兲. 45 A. Matsumoto, J. Phys. B 21, 2863 共1989兲. 46 M. Bancewicz, J. Phys. A 31, 3461 共1998兲. 47 E. Kauppi and L. Halonen, J. Phys. Chem. 94, 5779 共1990兲. 48 M. Mori, J. Comput. Appl. Math. 12–13, 119 共1985兲. 49 Subroutines DEFINT and DEHINT available from the Gams software repository at http://gams.nist.gov 50 N. Godbout, D. R. Salahub, J. Andzelm, and E. Wimmer, Can. J. Chem. 70, 560 共1992兲.

Downloaded 21 Sep 2001 to 128.40.5.14. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp