THE JOURNAL OF CHEMICAL PHYSICS 139, 194304 (2013)

Ab initio potential energy surface for the highly nonlinear dynamics of the KCN molecule H. Párraga,1 F. J. Arranz,1,a) R. M. Benito,1,b) and F. Borondo2,c) 1

Grupo de Sistemas Complejos, ETSI Agrónomos, Universidad Politécnica de Madrid, 28040 Madrid, Spain 2 Departamento de Química and Instituto de Ciencias Matemáticas (ICMAT), Universidad Autónoma de Madrid, Cantoblanco, 28049-Madrid, Spain

(Received 9 September 2013; accepted 30 October 2013; published online 18 November 2013) An accurate ab initio quantum chemistry study at level of quadratic configuration interaction method of the electronic ground state of the KCN molecule is presented. A fitting of the results to an analytical series expansion was performed to obtain a global potential energy surface suitable for the study of the associated vibrational dynamics. Additionally, classical Poincaré surfaces of section for different energies and quantum eigenstates were calculated, showing the highly nonlinear behavior of this system. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4830102] I. INTRODUCTION

There is an increasing interest in the study of chaos in realistic systems, namely, those susceptible of experimental verification. The vibrational dynamics of floppy molecules clearly fall into this category,1 being also relevant in many interesting chemical processes, such as unimolecular and isomerization reactions, photodissociation, or intramolecular energy relaxation.2 Laser technology and other new experimental techniques have made the experimental observation of these processes possible,3 even in real time.4 These studies are also interesting in transition state theory, which has recently been the subject of a renewed interest, after the improvement brought by the introduction of normally hyperbolic invariant manifolds in phase space.5–8 Here, the availability of benchmark realistic Hamiltonian models is important. Also, some interesting mathematical tools have been applied to the study of the quantum-classical correspondence9, 10 to this highly nonlinear realistic molecular systems. In this work, we present a global analytic potential energy surface (PES) for the KCN isomerizing system, based on accurate ab initio quantum chemistry calculations. The vibrational dynamics of KCN is a very interesting problem in nonlinear research for several reasons: First, it can be described realistically, second it is very highly nonlinear, and third it behaves as a two-dimensional system (due to adiabatic separation), thus providing mathematical simplicity without losing the essence of chaos. Two ab initio quantum chemistry calculations have been previously reported in the literature. The first one, published a long time ago,11 was carried out at a very low level of calculation, and incorrectly predicts the existence of only one minimum. The second one is accurate, but it is only restricted to the calculation of the different equilibrium points in the PES,12 not being then suitable for the dynamical studies that we pursue. Also, some preliminary classical and a) Electronic mail: [email protected] b) Electronic mail: [email protected] c) Electronic mail: [email protected]

0021-9606/2013/139(19)/194304/10/$30.00

quantum calculations on the dynamics of KCN are reported in this work, showing some interesting dynamical results and the way for further work. The organization of the paper is as follows. Section II is devoted to the PES calculations: in Subsection II A the coordinate system is defined, in Subsection II B details about ab initio calculations are given, in Subsection II C the fitting procedure is explained, and in Subsection II D the results are presented and discussed. Section III is devoted to the vibrational dynamics calculations: in Subsection III A the Hamiltonian model is introduced, and in Subsections III B and III C the classical and quantum calculations are presented and discussed. In Sec. IV, the present work is summarized and the possibilities for future work indicated. Finally, the relationship existing between the KCN normal modes and the coordinates used in our calculations is shown in the Appendix.

II. POTENTIAL ENERGY SURFACE A. Coordinate system

In order to obtain a PES suitable for vibrational calculations, Jacobi coordinates are used. As it is well known, these coordinates are the natural choice for studying the spectroscopically relevant stretching and bending motions taking place in triatomic molecules, especially if one of the bonds is much stronger than the others. For K–CN, the Jacobi coordinates (R, r, θ ) are defined as follows. Coordinate r directly corresponds to the bond length of the dimer C–N, coordinate R is the length between atom K and the center of mass of the C–N fragment, and finally, coordinate θ is given by the angle formed by the R and r directions. In this way, the linear configurations K–CN and CN–K correspond to angles θ = 0 and θ = π rad, respectively. In spectroscopic terms, r represents the stretching of the C–N bond, and R and θ the stretching and bending of the K–CN bond, respectively. Notice that the range for these coordinates are 0 ≤ R, r < ∞, and 0 ≤ θ < π .

139, 194304-1

© 2013 AIP Publishing LLC

194304-2

Párraga et al.

J. Chem. Phys. 139, 194304 (2013)

B. Ab initio calculations

7

R (a.u.)

All ab initio calculations have been performed using the G AUSSIAN13 software package. In order to account for the electron correlation, the quadratic configuration interaction method with single and double excitations, and also perturbative triple excitations [QCISD(T)] was selected. This ensures the size consistency14 needed for the PES calculation. In the geometry optimization calculations, the triple excitations were discarded to take advantage of the analytic gradients computation in the purely variational QCISD method. Besides, a 6-311+G(2d) G AUSSIAN standard split valence triple zeta basis set,15, 16 improved with diffuse functions,17 and also with polarization functions,18 is used in all calculations. The triple zeta basis set provides an accurate representation, and the diffuse and polarization functions add the needed flexibility to the basis set for nuclear configurations far from equilibrium. The default energy convergence value of 10−7 a.u. was used in all cases. Mass numbers corresponding to the most abundant isotope 39 K12 C14 N were selected; from them the actual masses, considering mass defect, are calculated by G AUSSIAN. The most stable equilibrium configuration, namely, the absolute minimum in the PES, was obtained using the Berny algorithm.19 This produced an equilibrium point at (R, r, θ ) = (4.903, 2.224, 0.547) (a.u., a.u., π rad), with normal mode frequencies of (ω1 , ω2 , ω3 ) = (159, 292, 2097) cm−1 . Furthermore, a second stable equilibrium point was found for the linear configuration K–CN at (R, r, θ ) = (6.141, 2.210, 0) (a.u., a.u., π rad), with frequencies (ω1 , ω2 , ω3 ) = (81, 273, 2162) cm−1 and also two activated complexes, i.e., saddle points, were found at the other linear configuration CN–K at (R, r, θ ) = (5.683, 2.223, 1) (a.u., a.u., π rad) with frequencies (ω1 , ω2 , ω3 ) = (i40, 304, 2113) cm−1 and for a triangular configuration at (R, r, θ ) = (5.893, 2.212, 0.175) (a.u., a.u., π rad) with frequencies (ω1 , ω2 , ω3 ) = (i88, 316, 2151) cm−1 . As it is shown in the Appendix, the normal modes (η1 , η2 , η3 ) correspond approximately to the Jacobi coordinates (θ , R, r) so that the normal modes frequencies (ω1 , ω2 , ω3 ) also correspond approximately to the Jacobi frequencies (ωθ , ωR , ωr ). A very important fact is that the frequency ωr is one order of magnitude greater than the other two frequencies, ωR and ωθ . Then, an adiabatic separation of the motion in the coordinate r from the rest can be expected, thus making it possible to study the KCN molecule as a two degrees of freedom system in (R, θ ) with r fixed at the equilibrium value req = 2.224 a.u. (absolute minimum) to a good approximation. (Notice that the equilibrium value of the coordinate r is very similar in all critical points found.) To obtain a representative region for the energy fitting to a series expansion, a grid in the configuration space (R, θ ) is defined as follows. First, the angular range 0 ≤ θ ≤ π is discretized by taking steps θ = π /18 rad (10 ◦ ). Then, for each angular value of θ i the coordinate R is discretized by taking steps of R = 0.1890 a.u. (0.1 Å), but keeping the electronic energy Eij ≤ 7000 cm−1 . In this way, the grid of 285 points shown in Fig. 1 was obtained, and then ab initio electronic energy calculations were performed at each point in the grid (with r = req fixed).

8

6

5

4

0

0.5

θ (π rad)

1

FIG. 1. Grid of 285 points (Ri , θ j ) used in the ab initio electronic energy Eij calculations. In all cases, Eij ≤ 7000 cm−1 . The color code corresponds to the distribution of the energy fitting absolute error in Fig. 2 (see text for details).

C. Series expansion fitting

An analytic expression for the PES, suitable for vibrational dynamics calculations, is then obtained by fitting the ab initio electronic energy values in the grid (Ri , θ j ) to the series expansion nmax m max   Cnm [1 − e−α(R−R0 ) ]n cos(mθ ), V (R, θ ) =

(1)

n=0 m=0

which contains (nmax + 1) × (mmax + 1) linear fitting parameters Cnm and two nonlinear ones α and R0 . This series expansion can be also expressed as follows: For a given θ , we have

Vθ (R) =

nmax 

Aθn [1 − e−α(R−R0 ) ]n ,

(2)

n=0

with the θ dependence being in the coefficients Aθn defined as Aθn =

m max 

Cnm cos(mθ ).

(3)

m=0

The expansion in the R coordinate in Eq. (2) corresponds to the perturbed Morse potential introduced by Huffaker,20 which is widely used for diatomic molecules, while the expansion in the coordinate θ in Eq. (3) is an even Fourier series. Let us remark that the Huffaker series has a radius of convergence given by (R0 − log 2/α, ∞),21 so that no fitting points (Ri , θ j ) should be with a value of Ri outside this interval. The actual fitting procedure is as follows: For each fixed value of the angular coordinate θ j , the Nj energy Eij points computed at (Ri , θ j ) were fitted to the series expansion of

Párraga et al.

194304-3

J. Chem. Phys. 139, 194304 (2013)

TABLE I. Numerical values of the fitting parameters in the KCN potential energy surface series expansion of Eq. (1). α = 0.3500000 a.u.−1

R0 = 4.8000000 a.u. Cmn (a.u.)

m\n 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

0

1

2

3

4

5

0.0300883 0.0218704 0.0365412 0.0072045 0.0077335 0.0011877 0.0008668 0.0001767 0.0001047 0.0000256 0.0000121 − 0.0000014 − 0.0000046 − 0.0000094 − 0.0000083 − 0.0000095 − 0.0000071 − 0.0000081 − 0.0000030

− 0.2134409 − 0.1004659 − 0.2334365 − 0.0310209 − 0.0400981 − 0.0068288 − 0.0060092 − 0.0015509 − 0.0010527 − 0.0002911 − 0.0001376 0.0000799 0.0000893 0.0001803 0.0001251 0.0001536 0.0000813 0.0000945 0.0000251

0.5896668 0.0945985 0.4282943 0.0427741 0.0820251 0.0173076 0.0184400 0.0059104 0.0045240 0.0011267 0.0006148 − 0.0007471 − 0.0005417 − 0.0012336 − 0.0006876 − 0.0009339 − 0.0003228 − 0.0003567 − 0.0000449

− 0.4206166 0.0530252 − 0.3007120 − 0.0159284 − 0.0900044 − 0.0241703 − 0.0325037 − 0.0113799 − 0.0096980 − 0.0018664 − 0.0014319 0.0026259 0.0013756 0.0039793 0.0017723 0.0026254 0.0005810 0.0004236 − 0.0001024

0.0729087 − 0.1331254 0.0266408 − 0.0134304 0.0598689 0.0186429 0.0322706 0.0105288 0.0099876 0.0011834 0.0017249 − 0.0039172 − 0.0013991 − 0.0062265 − 0.0021972 − 0.0033585 − 0.0005290 0.0001918 0.0003809

0.0578923 0.0655379 0.0455366 0.0114334 − 0.0209769 − 0.0064511 − 0.0137776 − 0.0035557 − 0.0038594 − 0.0000828 − 0.0008255 0.0020593 0.0003409 0.0038301 0.0010722 0.0015222 0.0002432 − 0.0004699 − 0.0003003

Eq. (2), defining the quadratic error 2j

tion is iterated, increasing the expansion degree, until the root mean square σ C from the differences between ab initio data and fitted expansion values is σ C < 1 cm−1 .

Nj   2 Eij − Vθj (Ri ) = i=1

=

Nj 

 Eij −

i=1

nmax 

 n Aθj n 1 − e−α(Ri −R0 )

2 ,

(4)

n=0

and numerically solving the minimization equation ∂2j ∂Aθj n

= 0,

(5)

to obtain the optimal values of coefficients Aθj n for some initial guess of the nonlinear parameters α and R0 and the expansion degree nmax . Also, the root mean square σ A from the differences between ab initio energy data and the values given by the fitted expansion is computed. This calculation is iterated by changing the initial values of the nonlinear fitting parameters, guided by the numerical value of the gradient of σ A in the (α, R0 ) space, and also increasing the expansion degree, until the root mean square σ A < 1 cm−1 . Afterwards, the optimized coefficients Aθj n are introduced into the even Fourier series (3), and by using the quadratic error 2  m M max   2 n = Cnm cos(mθj ) , (6) Aθj n − j =1

m=0

with M = 19 (number of grid points in θ ) and solving the minimization equation ∂2n = 0, ∂Cnm

(7)

the linear fitting parameters Cnm are obtained, for an initial value of the expansion degree mmax . As before, the calcula-

D. Results and discussion

The numerical values of the optimized parameters in Eq. (1) for the KCN global PES obtained with the fitting procedure described in Sec. II C are given in Table I. Notice that the radius of convergence for the Huffaker series of (R0 − log 2/α, ∞) = (2.82, ∞) a.u. has been taken into account, and accordingly all values Ri in the grid of points used are within this value (see Fig. 1). A measure of the quality of our fitting procedure is obtained from the value of the absolute error fit , i.e., the difference between the ab initio energies at the grid points and the corresponding energies given by the series expansion. Because the convergence in the ab initio energies is 10−7 a.u., only differences up to 10−6 a.u. are considered. The corresponding distribution of point absolute errors is depicted in Fig. 2, where it can be seen that, at worst, only two points from a total of 285 have an error fit = 4 × 10−6 a.u. = 0.878 cm−1 . That is, the fitted series expansion reproduces the ab initio data to within 1 cm−1 in all cases. The distribution of errors fit in the configuration space, (Ri , θ j ) is essentially homogeneous, as shown by the color code in Fig. 1. The KCN PES obtained from the fitted series expansion (1) is shown in the bottom panel of Fig. 3. As can be seen, there are two minima: the deep triangular minimum at θ ≈ π /2 rad and the shallow linear minimum at θ = 0 and there are also two saddles: the curved triangular saddle at θ ≈ π /6 rad and the quasi-flat saddle at θ = π rad. This is a qualitative difference with the PES of Wormer and

194304-4

Párraga et al.

J. Chem. Phys. 139, 194304 (2013)

2

3

−1

Veq (10 cm )

200

N

150

137

1

0

100

8 67 58 50 7 21 2 0

1

2

Δ

fit

−6

(10

3

R (a.u.)

0

4

a.u.)

FIG. 2. Distribution of the energy fitting absolute error at the 285 calculated ab initio points.

5

Tennyson,11 which predicts only the triangular minimum, being saddles the two linear configurations in their calculation. The minimum energy path (MEP) connecting all equilibrium points (minima and saddles) was numerically calculated, and fitted to the (even) Fourier series Req (θ ) =

9 

an cos(nθ ),

(8)

n=0

where the optimal values of the coefficients an are given in Table II. The MEP is plot superimposed to the PES in Fig. 3 (bottom), and the corresponding energy profile along the MEP, Veq (θ ) = V [Req (θ ), θ ] is also given in the top panel. Here, the deep or shallow character of the minima, and the curved or flat form of the saddles is clearly visible. The molecular properties, i.e., bond lengths, electronic energy, and normal mode frequencies obtained from the ab initio QCISD optimized equilibrium points are reported in Table III for the minima, and in Table IV for the saddles, respectively. The corresponding values obtained from the fitted PES are also shown in the tables for comparison. As can be seen, there is a good agreement in all cases between these two calculations. Actually, the greatest disagreement corresponds to the electronic energy values (obviously except for the absolute minimum, which is taken zero by definition). This is not surprising, since it is due to the fact that the ab initio calculations for the fitting data include perturbative triple excitations. Moreover, we also include in the tables other values existing in the literature, i.e., the results from the ab initio coupled cluster with single, double, and non-iterative triples substitutions [CCSD(T)] of Lee et al.,12 and the experimental results of van Vaals et al.,22 Ismail et al.,23 and Leroi and Klemperer.24 As can be seen, there is a good agreement with our results, in all cases.

6

4

0

0.5

1

θ (π rad)

FIG. 3. (Bottom) Potential energy surface for KCN obtained by fitting the analytical expression (1) to ab initio energies computed at the grid of points shown in Fig. 1. It is represented with contour plots spaced 700 cm−1 and using also color scale with darker colors representing smaller values of the energy. The minimum energy path connecting minima and saddles in the surface has been plotted superimposed in solid blue line. (Top) Potential energy profile along the minimum energy path.

With respect to our calculations in Tables III and IV there is one point worth mentioning. The two normal modes frequencies obtained from the series expansion (1) are calculated in a standard fashion,25 by using the vibrational Hamiltonian function in Jacobi coordinates defined below in Sec. III A. In this way, the normal modes frequencies ω1 and ω2 at the equilibrium point (Req , θeq ) are given by the roots of the fourth order equation Aω4 + Bω2 + C = 0,

(9)

TABLE II. Numerical values of the fitting coefficients for the minimum energy path series expansion in Eq. (8). n

an (a.u.)

n

an (a.u.)

0 1 2 3 4

5.4462444732 0.2334564415 0.4998689054 − 0.0022285587 − 0.0386511943

5 6 7 8 9

− 0.0010130958 0.0058189823 − 0.0003533387 − 0.0014070227 0.0000837943

194304-5

Párraga et al.

J. Chem. Phys. 139, 194304 (2013)

TABLE III. Comparison of bond lengths (in a.u.), electronic energy (in cm−1 ), and normal modes frequencies (in cm−1 ), obtained from the ab initio optimization, the series expansion fitting, and the literature, for the two minima in the KCN PES. Triangular minimum

Ab initio QCISD Fitted series expansion Lee et al.12 [CCSD(T)] van Vaals et al.22 (expt.) Ismail et al.23 (expt.) Leroi et al.24 (expt.) a

Linear minimum

KC , KN , CN

E

ω1 , ω2 , ω3

KC , KN , CN

E

ω1 , ω2 , ω3

5.22, 4.86, 2.22 5.19, 4.87, 2.22 5.16, 4.82, 2.23 5.14, 4.82, 2.21 ..., ..., ... ..., ..., ...

0 0 0 ... ... ...

159, 292, 2097 160, 291, . . . 166, 299, 2048 157, . . . , . . . 139, 288, 2050 207, . . . , 2158

4.95, 7.16, 2.21 4.95, 7.17, 2.22a 4.91, 7.12, 2.21 ..., ..., ... ..., ..., ... ..., ..., ...

1420 1376 1530 ... ... ...

81, 273, 2162 80, 275, . . . ..., ..., ... ..., ..., ... ..., ..., ... ..., ..., ...

Fixed r coordinate value corresponding to the absolute (triangular) minimum req = 2.22 a.u.

where 1 A = μ1 μ12eq , 4

momenta, and V (R, r, θ ) is the potential energy function. Then, making the two degrees of freedom approximation obtained by fixing r to its equilibrium value r = req the Hamiltonian function becomes  1 1 Pθ2 PR2 + V (R, req , θ ), (13) + + H = 2 2μ1 μ1 R 2 μ2 req 2

(10a)



2  ∂ 2V ∂ V 1 + μ12eq B=− , (10b) μ1 4 ∂θ 2 eq ∂R 2 eq C=

1 4



2

∂ V ∂R 2

and

eq

 μ12eq =

2

∂ V ∂θ 2



− eq

1 1 + 2 2 μ1 Req μ2 req

2

∂ V ∂R∂θ

where the two-dimensional potential V (R, req , θ ) is given by the fitted PES of Eq. (1).

2  , (10c) eq

B. Classical dynamics

−1

The Hamilton equations of motion corresponding to Eq. (13) are numerically integrated to obtain classical trajectories. In order to get a suitable graphical representation of the corresponding phase space, composite Poincaré surfaces of section (PSSs) taking the MEP as the sectioning plane have been calculated for different energies. For this purpose, we use the analytical expression for the MEP, Req (θ ) given in Eq. (8), and make the PSS an area-preserving map with the following canonical transformation:

(11)

.

Observe that the required derivatives are very easily obtained analytically from Eq. (1). The use of analytic derivatives is very important in order to ensure a good accuracy in both the classical and quantum vibrational calculations that are described below.

ρ = R − Req (θ ), III. VIBRATIONAL DYNAMICS

Pρ = PR ,

A. Hamiltonian model

ϑ = θ,   Pϑ = Pθ + PR dReq (θ )/dθ .

(14)

In this way, at a given energy E the PSS is defined as the (ϑ, Pϑ ) pairs making ρ = 0 and being in a predetermined branch (the negative one in the present calculations) in the second order equation for Pρ that arise from the Hamiltonian conservation H(ρ, ϑ, Pρ , Pϑ ) = E. Finally, all the PSS points are folded into the interval ϑ ∈ (0, π ) to take into account the symmetry of the molecular system.26 A representative selection of such composite PSS’s is shown in Fig. 4. The panels in the top row illustrate the onset

The rotationless (J = 0) vibrational Hamiltonian for KCN in Jacobi coordinates (R, r, θ ) is given by

2 Pθ PR2 Pr2 1 1 H = + + + + V (R, r, θ ), 2μ1 2μ2 μ1 R 2 μ2 r 2 2 (12) where μ1 = mK (mC + mN )/(mK + mC + mN ) and μ2 = mC mN /(mC + mN ) are reduced masses, mX the corresponding atomic masses, PR , Pr , and Pθ are the conjugate

TABLE IV. Comparison of bond lengths (in a.u.), electronic energy (in cm−1 ), and normal modes frequencies (in cm−1 ), obtained from the ab initio optimization, the series expansion fitting, and the literature, for the two saddles in the KCN PES. Triangular saddle

Ab initio QCISD Fitted series expansion Lee et al.12 [CCSD(T)] a

Linear saddle

KC , KN , CN

E

ω1 , ω2 , ω3

KC , KN , CN

E

ω1 , ω2 , ω3

4.92, 6.78, 2.21 4.91, 6.78, 2.22a ..., ..., ...

1539 1492 1613

i88, 316, 2151 i87, 318, . . . ..., ..., ...

6.88, 4.66, 2.22 6.88, 4.66, 2.22a ..., ..., ...

834 940 1087

i40, 304, 2113 i37, 305, . . . ..., ..., ...

Fixed r coordinate value corresponding to the absolute (triangular) minimum req = 2.22 a.u.

194304-6

Párraga et al.

E = 45

Pϑ (a.u.)

9

E = 65

9

0

0.55

0.65

E = 1300

40

0

−9 0.45

0.55

0.65

E = 1500

40

−9 0.45

0.65

E = 2400

40

0

0.55

0

ϑ

0

E = 145

9

0

−9 0.45

P (a.u.)

J. Chem. Phys. 139, 194304 (2013)

−40

0

0.5

ϑ (π rad)

1

−40

0

0.5

ϑ (π rad)

1

−40

0

0.5

ϑ (π rad)

1

FIG. 4. Composite Poincaré surfaces of section using the minimum energy path as the sectioning plane for increasing values of the vibrational energies.

of chaos as energy increases, following the Poincaré-Birkhoff and Kolmogorov-Arnold-Moser theorems.27 Thus, for E = 45 cm−1 the behavior of the system is completely regular, and all visible structures in phase space are invariant tori, or correspond to chains of islands. For E = 65 cm−1 , the principal (central) elliptic point is seen to undergo a pitchfork bifurcation, at which the original fixed point becomes hyperbolic (unstable), and two new elliptic (stable) points appear. As energy increases, the separatrix of this new central hyperbolic point breaks, giving rise to a noticeable stochastic band, thus marking the onset of widespread chaos (Chirikov’s mechanism27 ); this band is clearly seen in the third top panel corresponding to E = 145 cm−1 . Notice that in this system the phase space is already mixed at very low values of the energy, with large regions of chaos in between the regular tori. As it will be shown below in Sec. III C, this low energy threshold is relevant for the (very irregular) nodal structure of the corresponding quantum eigenstates. Panels in the bottom row show the phase space evolution at higher energies. Notice the very different axes range used here, as compared to those in the plots of the top row. As energy passes 145 cm−1 , all visible tori in the regular region are destroyed, and the chaotic region spreads throughout all the available phase space, as can be seen in the leftmost bottom panel corresponding to E = 1300 cm−1 . For E = 1500 cm−1 shown in the central bottom plot, the triangular and linear minima regions are connected, the motion around the linear minimum is mostly regular, while the dynamics around the triangular minimum is chaotic, and more remarkably, a small regular structure (tori) emerges at the linear saddle at (ϑ, Pϑ ) = (π , 0) in a region that was fully chaotic at lower energies. For E = 2400 cm−1 , chaos is (quasi-) fully developed over the available phase space since the region around the linear

minimum is now chaotic, although the little regular structure at the linear saddle still persists, having become larger and even showing now a clear (new) chain of islands structure. We think that the emergence of this regular structure at the linear saddle, with a periodic orbits bifurcation skeleton, is related to the normally hyperbolic invariant manifolds and its bifurcations as energy increases above the saddle point. This is a topic of interest that has been recently studied in the literature of open systems,28, 29 and also observed in closed systems due to saddle-node bifurcations,30 and it certainly deserves further consideration in KCN. C. Quantum eigenstates

The eigenfunctions and eigenenergies of the Hamiltonian operator associated to Eq. (13) were obtained using the discrete variable representation-distributed Gaussian basis method of Baˇci´c and Light31 and a final basis set of 2417 elements in 50 angular rays. In this way, we obtained the 400 low lying eigenfunctions, R θ |n n = 1, . . . , 400 with its eigenenergies converged to within 0.5 cm−1 . The numerical values of the eigenenergies corresponding to the first 100 eigenstates are given in Table V. By examining the nodal patterns32 of the lowest lying wavefunctions, it can be observed that they look very irregular, in the sense that nodal lines parallel and/or perpendicular to the MEP cannot be found, with the possible exception of the first excited state. For example, we show in Fig. 5 the first 9 eigenstates, which are all of them localized at the triangular minimum. When compared with the results of Tennyson and van der Avoird33 and Henderson et al.,34 both using the PES of Wormer and Tennyson,11 we see that the three calculations show a very early onset of irregularity in the wave functions,

194304-7

Párraga et al.

J. Chem. Phys. 139, 194304 (2013)

TABLE V. Rotationless vibrational eigenenergies (in cm−1 ) of the first 100 eigenstates. n

En

n

En

n

En

n

En

n

En

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

223 377 491 541 623 694 743 787 836 875 910 971 984 1030 1053 1076 1098 1123 1161 1167

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

1202 1220 1260 1278 1305 1318 1342 1368 1384 1396 1417 1437 1461 1487 1489 1512 1531 1550 1576 1587

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1589 1615 1625 1630 1654 1668 1684 1693 1707 1714 1739 1745 1755 1782 1796 1805 1817 1831 1854 1865

61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80

1872 1873 1883 1910 1918 1925 1948 1965 1969 1982 1992 2005 2008 2026 2031 2036 2063 2077 2084 2093

81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100

2107 2113 2135 2143 2153 2161 2167 2173 2195 2203 2219 2223 2231 2255 2259 2271 2291 2295 2301 2313

and ours are in general more irregular. The irregular nodal pattern is in general a quantum manifestation of the classical chaotic behavior, i.e., a signature of quantum chaos in the system. However, in KCN we have also the presence of a Fermi resonance between ω1 (mainly the R stretch) and ω2 (mainly θ bending). As it was shown in Sec. III B, the principal elliptic fixed point in the KCN classical dynamics bifurcates at

n=1

n=2

n=3

n=4

n=5

n=6

n=7

n=8

n = 40

n = 61

n = 85

n = 112

n = 142

n = 172

FIG. 6. Same as Fig. 5 for states n = 40, 61, 85, 112, 142, 172, with eigenenergies En = 1587, 1872, 2153, 2431, 2712, 2973 cm−1 , localized at the linear minimum.

very low energy, thus leading to stochastic bands and chains of islands. At the eigenenergy of the ground state, the classical phase space of KCN is mixed, containing large regions of chaos. Then, the eigenstates are not quantized over tori around the principal elliptic fixed point, and some do in the bifurcated chains of islands, thus giving rise to the observed irregular nodal pattern. Notice also that the potential around the triangular minimum is highly anharmonic and presents high mode couplings, so that normal modes (harmonic approximation) cannot properly account for the corresponding vibrational spectrum. Therefore, the anharmonic stretch and bending frequencies (ω1 , ω2 ) are more adequately obtained from the energy differences E2 − E1 and (E3 + E4 )/2 − E1 (due to the 1:2 Fermi resonance), respectively. This results in a value of (ω1 , ω2 ) = (154, 293) cm−1 , in agreement with the most recent experimental data22, 23 in Table III. Other eigenstates, such as those shown in Fig. 6, are localized on the linear minimum at θ = 0. In the figure, we have chosen states n = 40, 61, 85, 112, 142, and 172, which correspond to the quantum numbers (nR , nθ ) = (0, 0), (1, 0), (2, 0), (3, 0), (4, 0), (5, 0), respectively. Notice that all states of this class are excited in the coordinate R since any excitation in θ leads to the delocalization of the state over all the available

n = 101

n = 130

n = 164

n = 200

n = 238

n = 277

n=9

FIG. 5. Wavefunctions of states n = 1 − 9 of KCN, with eigenenergies En = 223, 377, 491, 541, 623, 694, 743, 787, 836 cm−1 , localized at the triangular minimum. The color scale indicates their values, with the darker colors corresponding to the higher values and the warm/cold tones indicating opposite signs. The minimum energy path and the corresponding isoeigenenergy contour have been represented in blue and black solid lines, respectively. The axes are the same as in Figs. 1 and 3 (bottom).

FIG. 7. Same as Fig. 5 for states n = 101, 130, 164, 200, 238, 277, with eigenenergies En = 2320, 2603, 2907, 3200, 3490, 3772 cm−1 , localized at the linear saddle.

194304-8

Párraga et al.

J. Chem. Phys. 139, 194304 (2013)

n = 242

n = 267

n = 324

n = 337

n = 371

n = 383

FIG. 8. Same as Fig. 5 for some hinge (top panels) and rotor-like (bottom panels) states of KCN, corresponding to n = 242, 267, 324, 337, 371, 383, with eigenenergies En = 3517, 3702, 4087, 4171, 4385, 4465 cm−1 .

θ domain. These eigenstates are initially quantized over the elliptic classical structure that can be seen in the middle bottom panel of Fig. 4. Eigenstates similar to these were found by Henderson et al.,34 although the linear configuration K– CN corresponds to a saddle point in the PES of Wormer and Tennyson11 that they used. Also, eigenstates localized over the linear saddle do appear in our calculations. As an illustration, we show in Fig. 7 states n = 101 130, 164, 200, 238, 277, which correspond to the quantum numbers (nR , nθ ) = (4, 0), (5, 0), (6, 0), (7, 0), (8, 0), (9, 0), respectively. Notice that these states are initially quantized over the emerging elliptic classical structure seen in the right bottom panel of Fig. 4. As in the K–CN linear configuration case, similar eigenstates were also found by Henderson et al.,34 although in that case the linear configuration CN–K was also a saddle point in the PES11 that they used. Finally, we show in Fig. 8, some eigenstates which are localized over classical periodic orbits, i.e., scarred states,35 similar to those found by Henderson et al.34 In the top panels, some hinge states are depicted, and in the bottom ones some rotor-like states are presented. As discussed in previous works,9, 10 the correlation diagram of vibrational eigenlevels versus, for example, the value of the Planck’s constant (artificially varied) and the distribution of zeros of the Husimi function are very powerful tools for the study and understanding of quantum chaos in molecular systems like KCN. Such study will be presented in a future publication. IV. SUMMARY AND FURTHER WORK

The main goal of this work was the construction of an accurate analytical expression for the PES of KCN that is suitable for detailed dynamical calculations. KCN is a highly nonlinear system, whose vibrational dynamics can be realistically described with a simple two-dimensional model. These characteristics make it a very interesting benchmark system in the field of nonlinear research and quantum chaos. For this purpose, we have performed an ab initio accurate calculation using the G AUSSIAN package at the QCISD(T)/6-

311+G(2d) level in a grid of points for the ground state of the K–CN isomerizing system, keeping the C–N length fixed at its equilibrium value. The results derived from this calculation are in good agreement with previous theoretical12 and experimental22–24 data reported in the literature. Furthermore, this collection of energy points has been fitted to an analytical expression in the vibrational coordinates consisting of a Huffaker20 series in the R-coordinate and an even Fourier series in θ -coordinate. The fitted series expansion reproduces the ab initio data to within 1 cm−1 . Also, the MEP connecting the equilibrium points was calculated and fitted to an even Fourier series. Both expressions are useful for our dynamical calculations. Additionally, some preliminary classical and quantum dynamical calculations have been performed, showing the highly nonlinear behavior of the KCN system and its effects. At very low energy, the principal stable classical periodic orbit bifurcates becoming unstable, so that only some chains of islands remain as regular structures for quantization of the corresponding vibrational states. Accordingly, the low lying quantum eigenstates have a very irregular nodal structure, except for ground and first excited state n = 1, 2. As energy increases, some interesting regular structures appear both in the classical and quantum dynamics of the molecule. Another relevant result in our study is the emergence at high energies of an elliptic point at the linear saddle, with the corresponding chains of islands structure around it. This is probably related to the associated normally hyperbolic invariant manifolds and its bifurcations, as the energy increases above the saddle point. This is a topic of recent interest,28, 29 and certainly deserves further study. On the other hand, regular R-excited eigenstates appear at both linear configurations (minimum and saddle). Also, eigenstates localized over periodic orbits (scars) emerge as energy increases. As it has been previously shown by us,9, 10 the different vibrational states of floppy molecules can be well characterized by the distribution of zeros of the Husimi function and the correlation diagram of eigenlevels versus Planck’s constant. This will be also the subject of a future publication. ACKNOWLEDGMENTS

This research was supported by the Ministry of Economy and Competitiveness-Spain under Grant Nos. MTM201239101 (F.J.A. and R.M.B.) and ICMAT Severo Ochoa SEV2011-0087 (F.B.). APPENDIX: NORMAL MODES AND JACOBI COORDINATES

In this appendix, we show graphically the relationship existing between Jacobi coordinates (R, r, θ ) and normal modes (η1 , η2 , η3 ) at each equilibrium point of the KCN molecular system, where the normal modes indices are taken in an increasing order of the frequency values, i.e., ω1 < ω2 < ω3 . For this purpose, the relationship between Jacobi coordinates and normal modes is obtained in a standard way by diagonalization of the force constant matrix in mass weighted

194304-9

q∼

1

Párraga et al.

J. Chem. Phys. 139, 194304 (2013)

coordinates R and θ . Moreover, at the two linear equilibrium configurations the normal mode η2 is essentially described the Jacobi coordinate R being the normal mode η1 approximately −θ although there is a little contribution from r. As can be observed, this latter (r) behavior is nearly quadratic, and then its relative contribution to the normal mode tends to zero in the limit of very small oscillations. Finally, at the two triangular equilibrium configurations, the approximate relation (η1 , η2 ) ≈ (θ , R) is worse. Although at the minimum this approximation could be acceptable (with a contamination of ∼20% from the other Jacobi coordinate), at the saddle the approximation fails (being here the contamination ∼40% in the mode η1 and ∼80% in η2 ). Clearly, a better approximation for the modes η1 and η2 has to be obtained for the two triangular configurations. This is easily done by using the coordinates along and perpendicular to the MEP defined in Eq. (14), i.e., (η1 , η2 ) ≈ (ϑ, ρ).

Triangular minimum (C ) s

A′

A′

A′

0 −1

q∼

1

Π

Linear minimum (C∞v ) Σ

Σ

0 −1

q∼

1

Triangular saddle (Cs ) A′

A′

A′

0

1 See,

−1

q∼

1

Π

Linear saddle (C∞v ) Σ

Σ

0 −1 −1

0

∼ η1

1 −1

0

∼ η2

1 −1

0

∼ η3

1

FIG. 9. Correspondence between normalized normal modes ηi (i = 1 2, 3)

and normalized Jacobi coordinates q at each equilibrium point, where q=R (green line), q = r (blue line), q = R eq θ (red line), for a cartesian elongation of [−1/4, 1/4] Å around the equilibrium. The symmetry point group of the molecule at the equilibrium point, and the corresponding irreducible representations of the normal modes are also indicated. The normal modes indices are taken in an increasing frequency order, i.e., ω1 < ω2 < ω3 .

Jacobi coordinates. Then, a cartesian elongation of [−1/4, 1/4] Å around each equilibrium point is introduced, and the corresponding values of the normal modes and Jacobi coordinates calculated. To obtain a better comparison, normalized versions of both normal modes and Jacobi coordinates are introduced as follows. The normalized ith normal mode is taken as ηi = (ηi − ηieq )/(ηi − ηieq )max , where the new coordinate is referred to the normal mode equilibrium value ηieq , and it is normalized by dividing by the maximum value (ηi − ηieq )max . In the same way, normalized Jacobi coordinates are defined qmax (with q = R, r, Req θ ), where qmax as q = (q − qeq )/ = max((R − Req )max , (r − req )max , Req (θ − θeq )max ), that is, the new coordinate is normalized with respect to the maximum value of the three maximum differences. Notice that, to define a dimensionally acceptable comparison for the angular coordinate, the arc length from the equilibrium point Req θ has to be used instead of the angle θ . In Fig. 9, we illustrate the relationship existing beη2 , η3 ) and normalized tween normalized normal modes, ( η1 ,

r, R θ ) at each equilibrium point of Jacobi coordinates, (R, eq the KCN PES. Notice that in all of four cases the normal mode η3 is essentially given by the Jacobi coordinate r (or −r at the linear saddle), being negligible the contribution from the other

for example: F. Borondo and R. M. Benito, in Frontiers of Chemical Dynamics, edited by E. Yurtsever, NATO ASI Series C Vol. 470 (Kluwer, Dordrecht, 1995), and references therein. 2 R. Schinke, Photodissociation Dynamics (Cambridge University Press, Cambridge, 1992). 3 V. S. Letokhov, Laser Spectroscopy of Highly Vibrationally Excited Molecules (Adam Hilger, Bristol, 1989); Molecular Dynamics and Spectroscopy by SEP, edited by H. L. Dai and R. W. Field (World Scientific, Singapore, 1995). 4 The Chemical Bond: Structure and Dynamics, edited by A. Zewail (Academic Press, San Diego, 1992). 5 T. Uzer, C. Jaffé, J. Palacián, P. Yanguas, and S. Wiggins, Nonlinearity 15, 957 (2002). 6 T. Komatsuzaki and R. S. Berry, Adv. Chem. Phys. 123, 79 (2002). 7 T. Bartsch, R. Hernandez, and T. Uzer, Phys. Rev. Lett. 95, 058301 (2005). 8 T. Bartsch, F. Revuelta, R. M. Benito, and F. Borondo, J. Chem. Phys. 136, 224510 (2012). 9 F. J. Arranz, L. Seidel, C. G. Giralda, R. M. Benito, and F. Borondo, Phys. Rev. E 82, 026201 (2010). 10 F. J. Arranz, L. Seidel, C. G. Giralda, R. M. Benito, and F. Borondo, Phys. Rev. E 87, 062901 (2013). 11 P. E. S. Wormer and J. Tennyson, J. Chem. Phys. 75, 1245 (1981). 12 D. Lee, I. S. Lim, Y. S. Lee, D. Hagebaum-Reignier, and G. Jeung, J. Chem. Phys. 126, 244313 (2007). 13 M. J. Frisch, G. W. Trucks, H. B. Schlegal et al., Gaussian 03, Revision C.01, Gaussian, Inc., Wallingford, CT, 2004. 14 J. A. Pople, M. Head-Gordon, and K. Raghavachari, J. Chem. Phys. 87, 5968 (1987). 15 R. Krishnan, J. S. Binkley, R. Seeger, and J. A. Pople, J. Chem. Phys. 72, 650 (1980). 16 J. P. Blaudeau, M. P. McGrath, L. A. Curtiss, and L. Radom, J. Chem. Phys. 107, 5016 (1997). 17 T. Clark, J. Chandrasekhar, G. W. Spitznagel, and P. v. R. Schleyer, J. Comput. Chem. 4, 294 (1983). 18 M. J. Frisch, J. A. Pople, and J. S. Binkley, J. Chem. Phys. 80, 3265 (1984). 19 C. Peng, P. Y. Ayala, H. B. Schlegel, and M. J. Frisch, J. Comput. Chem. 17, 49 (1996). 20 J. N. Huffaker, J. Chem. Phys. 64, 3175 (1976). 21 J. J. Camacho, A. Pardo, and J. M. L. Poyato, J. Chem. Soc. Faraday Trans. 90, 23 (1994). 22 J. J. van Vaals, W. L. Meerts, and A. Dymanus, Chem. Phys. 86, 147 (1984). 23 Z. K. Ismail, R. H. Hauge, and J. L. Margrave, J. Mol. Spectrosc. 54, 402 (1975). 24 G. E. Leroi and W. Klemperer, J. Chem. Phys. 35, 774 (1961). 25 H. Goldstein, C. Poole, and J. Safko, Classical Mechanics (AddisonWesley, 2001). 26 R. M. Benito, F. Borondo, J. H. Kim, B. G. Sumpter, and G. S. Ezra, Chem. Phys. Lett. 161, 60 (1989). 27 A. J. Lichtenberg and M. A. Lieberman, Regular and Chaotic Dynamics (Springer-Verlag, 1992).

194304-10 28 M.

Párraga et al.

Iñarrea, J. F. Palacián, A. I. Pascual, and J. P. Salas, J. Chem. Phys. 135, 014110 (2011). 29 H. Teramoto, M. Toda, and T. Komatsuzaki, Phys. Rev. Lett. 106, 054101 (2011). 30 F. Borondo, A. A. Zembekov, and R. M. Benito, Chem. Phys. Lett. 246, 421 (1995); F. Borondo, A. A. Zembekov, and R. M. Benito, J. Chem. Phys. 105, 5068 (1996).

J. Chem. Phys. 139, 194304 (2013) 31 Z.

Baˇci´c and J. C. Light, J. Chem. Phys. 85, 4594 (1986). M. Stratt, N. C. Handy, and W. H. Miller, J. Chem. Phys. 71, 3311 (1979). 33 J. Tennyson and A. van der Avoird, J. Chem. Phys. 76, 5710 (1982). 34 J. R. Henderson, H. A. Lam, and J. Tennyson, J. Chem. Soc. Faraday Trans. 88, 3287 (1992). 35 E. J. Heller, Phys. Rev. Lett. 53, 1515 (1984). 32 R.

The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp