Global Potential Energy Surface for HO2+ Using the CHIPR Method

Subscriber access provided by University of Winnipeg Library

A: Spectroscopy, Molecular Structure, and Quantum Chemistry 2+

Global Potential Energy Surface for HO Using the CHIPR Method George Densingh Xavier, Marco Martíñez, and António J.C. Varandas J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b12005 • Publication Date (Web): 01 Feb 2019 Downloaded from http://pubs.acs.org on February 3, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Global Potential Energy Surface for HO2+ Using the CHIPR Method F George D. Xavier,† M. Martínez González,†,‡ and A. J. C. Varandas∗,¶,† †Department of Chemistry, University of Coimbra, 3004-525 Coimbra, Portugal ‡Universidad de La Habana, Facultad de Química, calle San Lázaro sn., 10400 La Habana, Cuba ¶School of Physics and Physical Engineering, Qufu Normal University, 273165 Qufu, China E-mail: [email protected]

Abstract

lifetime and thermodynamic stability. An early extensive theoretical calculation 3 provided negative results concerning its formation and existence owing to the thermoneutral nature of the forward and reverse reactions. On the other hand, the investigation of its vibrational frequencies either in its free state or in tagged form gained momentum in recent years. The latest accurate work appeared in 2018 with the estimation of three modes of vibration by producing the target molecules from its precursor H2 O2 in laboratory conditions. Previously, there were a few attempts 5–8 to estimate frequencies for this ion, with the latest 4 being the most accurate as well as up to date. From the computational point of view, this ion has been extensively studied 3,9–15 starting from the eighties to understand the nature of bonding and interaction as well as to estimate its structural parameters. However, they were all restricted to compute a few selected cuts of the full dimensional PES and hence remained insufficient to shed light on the finer details of topograhical features occuring in the low and high energy landscapes. A good amount of progress has been made in this aspsect for its neutral counterpart O2H using double manybody expansion (DMBE) theory as well as the present CHIPR method. 16–19 But, the same is not true of the ionic case, due to lack of a full dimensional accurate global potential energy surface (PES) in the literature. A full dimensional

An analytical potential energy function for the title ion based on the Combined-HyperbolicInverse-Power-Representation (CHIPR) method and its characteristics are discussed at length in the present work. The curves of two diatomic ions, O2+ and OH+, are also obtained within the same approach. The model PES so obtained exhibits extraordinary flexibility in describing with subchemical accuracy even the weak topological features near the higher energy regions. Thus, structural properties predicted by the model may help spectroscopists who want to compare their experimental values with the ones from theory. The relaxed PESs in various coordinates have been calculated by relaxing the O2 bond distance using the present model, thus throwing light on all the possible isomers and their interconversions. The latest estimates of IR frequencies for three vibrational modes have been compared with the computed frequencies using the present model and the agreement seems encouraging.

1

Introduction

The hydroperoxyl ion HO2+, was thought to be present in the atmosphere after the presence of its neutral partner was confirmed. 1,2 Yet, preferably now, the inablity to prove its existence beyond doubt raises concerns as to its

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

global analytical surface describing its possible dissociation channels via four dimensional interpolations (energy included) were reported by us. 15 In this report, we revisit the problem using an analytical function based on a high degree polynomial of three independent variables written in terms of a contracted basis set of hyperbolic functions as proposed in CHIPR. 16–18 The accumulated and stratified errors in the fitting process are ensured to be as small as possible, preferably below 1 kcal mol−1 . The paper is organized as follows. In section 2, the methodology adopted in CHIPR is provided from a mathematical perspective. The results and discussion are in section 3. Specifically, section 3.1 is devoted to the fitting details of the diatomic molecular ions, then followed by section 3.2 on the triatomic fit and the properties as predicted from fitted potential form. This will be followed by the summary and conclusions in section 4.

Page 2 of 15

tial which is obtained by subtracting the sum of individual diatomic potentials from the total potential. In Varandas’ CHIPR method, the assumption is that one obtains first contracted functions, expressed as the linear combination of secant functions with appropriate weighting factors (to be known from a nonlinear least-squares fit) and the expression then used to fit the data points using a high-degree polynomial in the basis of contractions. For a triatomic molecule, three such contracted basis functions representing each mode of nuclear motion are required. Therefore, the model function will be a mixed polynomial with each of its terms involving the product among these three contracted functions. In general, the number of such functions required will be equal to degrees of freedom of nuclear motion of molecule. Therefore, we write under the CHIPR approach: 16–19 yi =

2

Method

M X

cij φij

(2)

j=1

where yi is the contracted basis function meant for a given degree of freedom of nuclear motion (channel) as specifed by the index i, cij are the unknown coefficients to be determined, φij is the j th hyperbolic secant function with the coefficient cij for ith channel and M is the total number of terms chosen in the basis, an integer. In general, the hyperbolic basis can be written as follows 16–19

According to theory of many-body expansion, 20,21 the total interaction potential of the triatomic system is decomposed into sum of all possible two-body and three-body interaction. Thence, we can write V (R1 , R2 , R3 ) =V2 (R1 ) + V2 (R2 ) + V2 (R3 ) + V3 (R1 , R2 , R3 ) (1) where Ri (i = 1 − 3) are interatomic coordinates of the triatomic system. In the present work, R1 is treated as the distance between two oxygen atoms and R2 , R3 refer to the distance between an oxygen and hydrogen atom. There are two such equivalent OH distances based on permutational symmetry. V refers to the total interaction potential and V2 and V3 are the twoand three-body potentials, respectively. Specifically, the first term on the right-hand-side of Eq. (1) refers to the O2+ interatomic potential and the next two terms refer to the permutationally equivalent OH potential due to the presence of two identical nuclei and the final term denotes the three-body interaction poten-

φij = sechηij (γij ρij )

(3)

ref with ρij = Ri − Rij where Ri is distance along ref th the i channel and Rij is the reference geomth etry for i channel and j th term in the expansion given by Eq. (2). To make the reference geometry independent of linear expansion and to achieve the distributed origin of secant funcref tions along the distance Ri , Rij is further expressed as 16–19 ref Rij = ζi (Riref )j−1

(4)

where ζi and Riref are the two channel dependent parameters. Combining all above details


2

Page 3 of 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for i = 1, j = 2, Eq. (2) can be explicitly given as h i ref 0 η11 γ11 (R1 − ζ1 (R1 ) y1 =c11 sech h i (5) + c12 sechη12 γ12 (R1 − ζ1 (R1ref )1

where τ is the total number of degrees of permutational of the molecule (τ = 3 for the ti˜ is highest degree of polynomial, tle ion), M ˜ . For the specific and i1 + i2 + · · · + iτ ≤ M + case of HO2 , where permutational symmetry exists, on exchanging H+ from one oxygen to the other, the number of terms in each power of polynomial reduces relative to the case where there is no such symmetry. For example, terms with power equal to 2, there are 3 terms without permutational symmetry but only 2 when there is symmetry and also for power equal to 3, there are 7 terms with no symmetry but only 4 with symmetry. Thence, the reduction in terms count leads us to write the above expression, Eq. (9), in slightly different form reflecting its symmetry behaviour 16–19

where c11 , c12 , η11 , η12 , ζ1 , R1ref are unknowns to be estimated from the fitting procedure. Equation (5) can be extended for any i and j values. Using a similar contracted basis, the proposed diatomic model potential assumes the form m Na Nb X V2 = Ck y1k R1 k=1

(6)

where NA , NB are the atomic number of A and B of AB diatomic molecule and m is the highest power of the polynomial. Setting m = 3, Eq. (6) becomes

V3 =

˜ M X

Ci1 i2 i3 y1i1 y2i2 y3i3 + y2i3 y3i2

(10)

i1 ,i2 ,i3 =0

V2 =

Na Nb C1 y11 + C2 y12 + C3 y13 . R1

(7)

where yi (i = 1 − 3) are the contracted basis functions for the three modes of motion. This form is used later to fit the three-body ab initio ˜ , where M ˜ = 15 points keeping i1 + i2 + i3 ≤ M in this work. In fact, the only difference is that the secant function in Eq. (3) is accompanied by an extra term which accounts for the shortrange behaviour, Ri → 0. So, Eq. (3) becomes i h  tanh βi (Ri − ζi (Riref )j−1  φij =  (Riref )j−1 h i × sechηij γij Ri − ζi (Riref )j−1 . (11)

where C1 , C2 , and C3 are the parameters to be evaluated from the fit, and y1 is defined in Eq. (5). For some potential energy curves, there may exist long range minima due to charge quadrupole-quadrupole or charge-induced (e.g., dipolar) moment interactions, Eq. (2) has to be modified as follows by adding an extra polynomial term, M X yi = cij qij φij (8) j=1

where qij is a pre-sech (i, j)-dependent polynomial that has previously 16–18 been assumed for simplicity as unit; thence yi refers to atompair (i, j). This type of contracted function is discussed when the fitting of O2 ab initio data points is taken up in the next section. In the case of the triatomic system here examined, to fit three-body potentials three contracted functions have to be used to write the global model function. It assumes the general form 20,21 V3 =

˜ M X i1 =0,i2 =0,··· ,iτ =0

Ci1 ,i2 ,··· ,iτ

τ Y

ypip

With the help of Eq. (11), the expression for the various yi can be explicitly given as h i  ref 0 tanh β1 (R1 − ζ1 (R1 )  y1 =c11  (R1ref )0 h i × sechη11 γ11 R1 − ζ1 (R1ref )0 h i  tanh β1 (R1 − ζ1 (R1ref )1  + c12  (R1ref )1 h i × sechη12 γ12 R1 − ζ1 (R1ref )1

(9)

p=1


3


i h tanh β1 (R1 − ζ1 (R1ref )2  + c13  ref 2 (R1 ) i h (12) × sechη13 γ13 R1 − ζ1 (R1ref )2 

Page 4 of 15

finitely seperated. One obtains the following fitting expression: V3 =C110 y1 (y2 + y3 ) + 2C011 y2 y3 + 2C111 y1 y2 y3 + C120 y1 (y22 + y32 ) + C210 y12 (y2 + y3 ) + C012 (y2 y32 + y22 y3 ) (15)

h i ref 0 tanh β2 (R2 − ζ2 (R2 )  y2 =c21  (R2ref )0 h i × sechη21 γ21 R2 − ζ2 (R2ref )0 i h  tanh β2 (R2 − ζ2 (R2ref )1  + c22  (R2ref )1 h i × sechη22 γ22 R2 − ζ2 (R2ref )1 i h  tanh β2 (R2 − ζ2 (R2ref )2  + c23  (R2ref )2 h i × sechη23 γ23 R2 − ζ2 (R2ref )2 (13) 

 y3 =c31 

h

tanh β3 (R3 −

ζ3 (R3ref )0

There are 3 coefficients for power two, but because of permutational symmetry, C110 = C101 , only two coefficients appear in Eq. (15). Similar arguments hold for power 3 or any other higher power. Thence, number of terms is reduced when symmetry arguments are taken into consideration. In this work, we have gone parame˜ = 15 and the total number of coeffients ters M in the polynomial expansion alone is 425 in additon to the parameters from the contracted basis defined in Eq. (12), Eq. (13), Eq. (14) where ηij is fixed at 1 and not subject to fitting.

3 3.1

i

Results and discussion O2+ and OH+ curves

The data points utilized for fitting are obtained via a double-level extrapolation 22,23 of ab initio points calculated at MRCI level of accuracy with cc-pVT Z and cc-pVQZ basis set in split form as discussed in a previous paper. 15 In the fitting expression, the sum in the contracted basis in both cases has gone up to M = 4 as defined in Eq. (2). In the special case of O2+, a linear function is multiplied with hyperbolic function to accurately fit long-range forces due to the leading charge-quadrupole interaction. The required basis functions are given below while their optimal parameters are in supplementary information (SI). For OH+ it reads h i y1 = c11 sechη11 γ11 R1 − ζ1 (R1ref )0 h i + c12 sechη12 γ12 R1 − ζ1 (R1ref )1 h i + c13 sechη13 γ13 R1 − ζ1 (R1ref )2 h i + c14 sechη14 γ14 R1 − ζ1 (R1ref )3 (16)

 (R3ref )0 i h ref 0 η31 γ31 R3 − ζ3 (R3 ) × sech i h  tanh β3 (R3 − ζ3 (R3ref )1  + c32  (R3ref )1 h i × sechη32 γ32 R3 − ζ3 (R3ref )1 h i  ref 2 tanh β3 (R3 − ζ3 (R3 )  + c33  (R3ref )2 h i sechη33 γ33 R3 − ζ3 (R3ref )2 (14)

Because the expansion for V3 is used to fit three-body interaction energies only, any constant term or two-body terms should be excluded from it. 20,21 Thus, ˜ = 3, the expansion coefficients, for M C000 , C100 , C010 , C001 , C200 , C020 , C002 , C300 , C030 , C003 are excluded as they represent diatomic energies and C000 , a constant term, denote asymptotic limit when three atoms are in-


4

Page 5 of 15

The Journal of Physical Chemistry 0.05

that anharmonicity in OH+ is roughly 3 times larger than O2+ and hence the harmonic approximation is no longer sufficient in calculating vibrational frequencies of OH stretching in the triatomic HO2+ ion. In turn, the long range interaction is dominated by the leading charge-induced dipole interaction in OH+ which shows a Morse-like behaviour, while in O2+ it is the charge-quadrupole interaction which contributes to a weak minimum in the long-range part of the potential on account of such energy stabilization.

0 -0.05 energy/Eh

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

OH+ fit

-0.10

O2 + fit OH+ data

-0.15

O2 + data

-0.20 -0.25 2

4

6

8

r(a0 )

Figure 1: CHIPR model potential energy curves for O2+ and OH+ diatomic ions. The rootmean-square-deviations relative to the calculated ab initio points are 0.4 cm-1 and 0.7 cm-1 , in the same order.

3.2

The HO2+ CHIPR PES

The raw ab inito points were calculated as a function of Jacobi coordinates at MRCI level + of theory using cc-pVTZ and cc-pVQZ basis while for for O2 it is i set and subsequently extrapolated to complete h ref basis set limit in split form with the doubley1 = c11 (s1 R1 + t1 )sechη11 γ11 R1 − ζ1 (R1 )0 22,23 adopted earlier 15 for the di h i level approach atomic molecules. The reference wave function + c12 (s2 R1 + t2 )sechη12 γ12 R1 − ζ1 (R1ref )1 h i of the CI expansion at each geometry was taken + c13 (s3 R1 + t3 )sechη13 γ13 R1 − ζ1 (R1ref )2 from the state averaged CASSCF wave func h i tion involving the five lowest electronic states + c14 (s4 R1 + t4 )sechη14 γ14 R1 − ζ1 (R1ref )3 of 3 A00 symmetry at the preceding step. The (17) magnitude of the Davidson correction ranges from −0.02 to −0.03 Eh over the entire calwhere the specific qij in Eq. (8) are in the above culation. The energies so-corrected were the equation the linear functions before the sech ones utilized in the fit. It should be noted term. The sum in the fitting expression as dethat such a correction affects only slightly the fined in Eq. (6) is fixed at m = 4 in the case agreement with experiment of the calculated of OH+ and m = 6 in the case of O2+. The remolecular properties. All computations were gionwise (repulsive, equilibrium and long range) performed using the MOLPRO 28 suite of proroot mean square deviations with arbitrary engrams. The grid details selected for each of ergy criterion falls roughly close to the overall ˚ = the Jacobi variables are given here: R(A) deviation calculated with entire data set. 0, 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, The equilibrium diatomic properties com1.6, 1.8, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4, puted with the CHIPR curve are found to ˚ = 0.8, 0.9, 0.95, 1, 1.05, 1.1, 4.5, 5, 6, 8, 10, r(A) be same as that of DMBE curve reported 1.14, 1.18, 1.2, 1.22, 1.24, 1.28, 1.3, 1.34, 1.38, 1.4, earlier. 15 For example, the equilibrium dis1.5, 1.6, 1.8, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, tance predicted from this model is 1.119 Å and 4.25, 4.5, 4.75, 5, and γ(deg) = 0, 15, 30, 45, 60, 1.028 Å for O2+ and OH+, respectively, while 75, 90. In the present work, r refers to the disthe corresponding DMBE values are 1.117 Å tance between the center-of-mass of O2 and H+, and 1.028 Å. Both set of predictions were very R refers to the diatomic distance and γ is the close to the available experimental values in the angle between R and r vectors. In valence coliterature. 24–27 Figure 1 shows the two curves ordinates, R1 is the O2 distance, hence equal to superimposed for better comparison. The difR in Jacobi coordinates, and R2 and R3 are the ference in curvature of the curves suggests two equivalent OH distances. In total, the grid ACS Paragon Plus Environment

5


Page 6 of 15

thesis: γ = 0◦ (439 cm−1 ), γ = 15◦ (300 cm−1 ), γ = 30◦ (210 cm−1 ), γ = 45◦ (201 cm−1 ), γ = 60◦ (230 cm−1 ), γ = 75◦ (240 cm−1 ), and γ = 90◦ (370 cm−1 ). The stratified rmsd of points below 10 kcal mol−1 is 78 cm-1 without weights and with weights, this reduces to 20 cm-1 . The selection of weights is based on stratified energy cuts with the minimum energy spacing of 10 kcal mol−1 and maximum of 100 kcal mol−1 near dissociation; in general, energy and weights are inversely related. In Table 2, the weights and specified energy range (in kcal mol−1 ) are given. Fig. 2 shows a simple 2D plot, where the predictive power of the final CHIPR potential function can be proven along all three Jacobi variables before examining more complicated optimized 3D surface plots. Clearly, the quality of the fit is very good, particularly near the global minima. In all three plots, the prediction from the analytical function agrees well with ab initio data within the accumulated rmsd except in the vicinity of the conical intersections which keep the general topography but smooth out the cusp. Although this could be remedied, 29–31 no attempt is here done to do so. This happens at γ = 90◦ and 0◦ , with the disagreement being relatively higher in comparison with other angular orientations. In the first plot, the curves ˚ are plotted as a function of R at fixed r = 1.24A for different γ values. The selected r value is very close to the O2 bond distance in the global minimum geometry of HO2+ predicted by the CHIPR PES. Some of the points in the perpendicular and collinear geometry do not fall in the line drawn from the model function, as it is apparent from the plot. In the middle plot, the curves are shown as a function of r at a fixed ˚ (the equilibrium distance in value of R = 1.3A CHIPR) for different γ values as indicated in the plot. The last panel is shown as a function of γ for various R values at a fixed r (1.24Å). ˚ passes through Note that the curve at R = 1.4A global minimum near 45◦ and shows a slightly ˚ near γ = 90◦ due larger deviation at R = 2A to approaching the conical intersection at Tshaped geometries. For other values of R at the same angle, this deviation is not worth notable due to the corresponding geometry lying

size amounts to 31 × 32 × 13 = 13, 888 points after replication for permutational symmetry. Of the total number of calculated points, not all were included, with an energy criterion adopted to define a cut-off limit (mostly leaving aside points in regions of high energy). Such a cut-off limit has been fixed at 0.5 Eh (314 kcal mol−1 ) relative to the minimum. This reduces the data set to 10,625 points which fall below 314 kcal mol−1 . To fit these points, we converted the Jacobi variables to interatomic ones (R1 , R2 , R3 ) using the cosine rule. On dissociation of HO2+, the diatomic species O2+ and OH+ are obtained, as discussed at length in the previous subseciton. To use Eq. (1), one then needs to compute the three-body energies of HO2+ at the above 10,625 discrete geometries. This was done by subtracting the sum of diatomic curves using the already obtained potential functions from the total interaction potential. This gives the three-body energy data set to be fitted to V3 (R1 , R2 , R3 ). Because there are three dissociation channels, two of them equivalent, we need to provide three contracted basis functions which were shown in the previous section by Eq. (12), Eq. (13), and Eq. (14) with the same number of primitives taken in each of three functions, m = 3 in the HO2+ case. The polynomial model function chosen in the contracted basis has the highest power of 15, with the total number of monomials equalling 413 parameters, in addition to the 36 equilibrium from the contracted functions. The whole set of 449 parameters obtained in the fit are provided in the supplementary information (SI). The number of monomials in each power drastically increases with increasing power, thereby enhancing the computational power required to perform the nonlinear least-squares fit. They were given up to order 6 elsewhere, 18 and listed here fully in Table 1. The accumulated root mean square deviation (rmsd) calculated for the fitted points works out to be 284 cm-1 , which is well within the accuracy of 1 kcal mol−1 . When looked carefully, Jacobi anglewise the rmsd is larger near linear and collinear geometries as well as at the minimum near the equilibrium region. It is listed inside the paren-


6

Page 7 of 15

Table 1: Number of monomials in each power of polynomial expansion power 2 3 number 2 4

4 7

5 6 10 14

7 8 18 23

9 10 11 28 34 40

12 13 47 54

14 15 62 70

Table 2: Choice of weights used in the present work range 0-10 10-20 20-40 40-60 60-80 80-100 100-150 150-250 weights 10 1 0.1 0.01 0.001 0.0001 0.0001 0.0001 far from the region of conical intersection. The dissociation occurs along two possible channels, O2+ + H and OH+ + O, and the long range interactions are dominated by electrostatic and induction forces in addition to dispersion forces. Because the data points involve distances up to 10 Å and 5 Å in (along the intermolecular and intramolecular coordinates, respectively) it is virtually impossible to account for them individually. This will be left to the CHIPR method itself, which is best suited for the purpose by taking care of all long-range forces globally through the direct fit to the ab initio data. The potential energy surface of a triatomic involves four dimensions and hence it will not be not possible to visualize it in a single plot by varying all three independent variables (shape coordinates) simultaneously. To see all features of the PES present in the entire range, one of the variables has been optimized (relaxed) leaving the other two fixed, with the procedure repeated to obtain a rectangular grid of uniformly spaced independent variables with optimized energy. In this way, all relevant topographical features can be seen in a 3D plot. In the present work, we show such plots in contour form for four distinct coordinate representations; see Fig. 3 and Fig. 4. Note that the contour lines are drawn for a set of specific values (thence not uniformly spaced with a constant increment), thus being selected such as to show even the fine details of the CHIPR PES. Clearly visible in all plots are two equivalent global-minimum geometries, which have been identified with M1 and M2 in the relaxed triangular plot in hyperspherical coordinates 32 shown in the top panel. Note that M1 and M2 are interconverted via a transitIon state geom-

Energy/Eh

-0.20

-0.25

-0.30

-0.35

γ=90 °

γ=30 °

γ=75 °

γ=15 °

γ=60 °

γ=0°

γ=45

-0.40 2

4

°

6

8

R(a0 ) 0

Energy/Eh

-0.1

-0.2

-0.3

γ=90 °

γ=30 °

γ=75 °

γ=15 °

γ=60 °

γ=0°

γ=45 °

-0.4 2

4

6 r(a0 )

-0.275

8

10

R=1.4Å

R=1.8Å

R=1.5Å

R=2Å

R=1.6Å

-0.300 Energy/Eh

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-0.325 -0.350 -0.375 -0.400 0

50

100 γ(deg)

150

Figure 2: One-dimensional cut along R, r and γ to illustrate the goodness of the fit with respect to the calculated ab initio points near equlibrium region (thence, points below the stratified cut up to 10 kcal mol−1 above the global minimum).


7


1.0

��

��

��

�

�� + +�

+�

γ*

+

0.5

��

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 15

��

0.0

��

��

-0.5

��

��

��

��

-1.0

�+�� + -1.0

0.0

-0.5

0.5

1.0

*

β

Figure 3: Countour relaxed triangular plot in hyperspherical coordinates 32 for ground state HO2+ PES. r = 2.1a◦ , R = 4.4a◦ , γ = π or 0, E = −0.2750Eh . The existence of a conical intersection along the OH stretching mode suggests that the bare proton (or H for the neutral) is sufficient enough to induce non-adiabaticity of electronic states at selected regions, rendering Born-Oppenheimer approximation invalid in those regions. Most structural features mentioned above fall below the H + O+ 2 dissociation channel, which occurs at about 0.25 Eh above the global minima. However, there are other relevant features occurring close to the other dissociation channel, OH+ + O which lies about

etry, TS, in an obvious similarity to what is found for the neutral HO2. 16–18 In close analogy to this, there are also long range minima. One is denoted M5 in the hyperspherical plot at a geometry with C2v symmetry, but cannot be seen in the other three panels. There are also three conical intersections, one at T-shaped geometries, the other two, permutationally equivalent, at collinear geometries (denoted by Ci1, Ci2, Ci3 in the same order). The estimated geometrical parameters of those conical intersections in Jacobi coordinates are r = 2.2a◦ , R = 3.5a◦ , γ = π2 , E = −0.2745Eh and


8

Page 9 of 15

8

0.35 Eh about the global minima. The most important of them in this high energy regions of the PES refer to the two equivalent H-bonded minima (denoted by M3 and M4) which occur at linear geometries and are likely to play a role in reaction dynamics at low and very low collision energies (very much like for the neutral case 17 ). They are interconnected via a low barrier saddle point whose structure can hardly be identified with confidence and hence is only marked with a X in Fig 4. Apart from this, there is a highly symmetrical linear stationary structure of D2h point group marked as M6. All minima have been confirmed through ab initio vibrational frequencies calculations and were reported previously, 15 while the frequencies obtained from the present CHIPR model potential have here been calculated 33 for the global minima and presented in Table 3 along with available experimental and previously theoretical results. The agreement for the weak OH stretching is found to be within the experimental uncertainity. The model predicts a zero point energy (ZPE) of 2866 cm-1 while prediction from ab initio calculation at the CBS limit using the uniform singlet- and triplet extrapolation 34 scheme is 2882 cm-1 . This is also another indication regarding the high quality of the ab initio model. Table 4 gathers the first 30 vibrational eigenenergies from the CHIPR PES with proper assignment of quantum states. They were computed by solving the nuclear Schrödinger equation using the discrete variable representation (3D-DVR) code of Tennyson and coworkers. 33 In turn, the geometrical parameters of the structures marked in the plots are given in Table 5 where ∆ refers to the energetic ordering of various isomers on the relative energy scale with respect to the global minimum. Naturally, the agreement of the calculated structural parameters with those obtained from the CHIPR PES supports the fact that the latter reproduces accurately the data points over the entire energy range. The previously reported 3 geometrical properties for the global minimum are r(O—O)=1.2272Å, r(O— H)=1.0102AA, α(deg)=112.743, calculated at coupled cluster with triples with enlarged aug-

7

r2 /a0

6 5

��

��

4

��

��

3 ��

2

��

2

3

��

4

5 r1 /a0

6

7

8

8

R/a0

6

4

��

��

��

2

��

��

0 0.0

��

��

0.5

1.0

1.5 2.0 γ/radian

2.5

3.0

8 6 y/a0

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4 2

��

�� 0 -8 -6 -4 -2

��

��

� ��

0 2 x/a0

��

4

6

8

Figure 4: Contour plot in various other coordinate representations to showcase topographical features. From top to bottom, the first panel is a valence coordinate plot, the secondd a plot in Jacobi coordinates, and the third a cartesian representation. In all panels, the O2 bond length is relaxed to attain minimum energy for specific values of the other two parameters.


9


Page 10 of 15

Table 3: Comparison of vibrational frequencies (cm-1 ) from CHIPR PES with experimental values of tagged as well as free HO2+. The uncertainties in experiments are given in parenthesis. Species Ne HO2+ He HO2+

HO2+

method exp:matrix exp:gas-phase exp:gas-phase exp:gas-phase exp:gas-phase exp:prediction ab initio ab initio ab initio ab initio CHIPR

OH str. 2760.8(2) 2915 2905(1) 3016.726(5) 3020(40) 3021.7 3028 3033 3030 3054

mented Dunning’s basis set. Table 5 gives almost the same values except the valence angle, with the discrepancy steming from the different ab initio methodology employed in both works. In Table 6, we compare the results from this work with past computational work for few selected equiilbrium properties. As seen, the estimates for a given property remain more or less the same in all four works, including the present one. However, the calculated can only be validated when reliable experimental estimates will be available in the future. As of now, only the fundamental vibrational frequencies are available in the high resolution infrared studies 4 which are supported by this work. Thermochemical properties such as proton affinity and enthalpy of dissociation were estimated as 103 kcal mol−1 and 261 kcal mol−1 and are in good agreement with experimental values 36,37 of 100.98 kcal mol−1 and 264 kcal mol−1 , respectively. To analyse the dissociation channels, we show in Fig. 5 contour plot as a function of r1 and r2 where r1 , r2 is defined as O—O and O—H distance. Top left is for α = 0◦ , top right for α = 90◦ and bottom for α = 180◦ . Looking at the plots it is clear that two weak minima are being developed as O—O and O—H bonds are breaking. At α = 180◦ , higher energy regions looks mirror image of OH+ + O channel and cross over the other via saddle point at 4.5 a0 . On the other side of divergence of potential energy function

OO str.

1387(7)

1379 1440 1406 1418 1399

bend 1026.5

refs.

8 7 1070(5) 4 7 4 7 1059 9 1068 3 1037 8 1040 35 1089 this work

(which are chopped off in the plot for clarity) occurs the low energy features with local minima, maxima and van der waals minima. The energy gap between these channel as predicted by α = 180◦ contour function is 1.697 eV as expected due to the requirement of enormous amount of energy in breaking the O—O bond. Therefore, O—H breaking is exothermic while O—O breaking is endothermic.

4

Conclusions

The CHIPR method has been found of great use in fitting all ab initio data points calculated for the title molecular cation with an accuracy well below 1kcal mol−1 . It has further been shown to fit well all long range features present in its ground-state adiabatic PES. Based on optimized contour plots in various coordinate systems, the topographical features of the final PES so obtained have also been fully discussed. Moreover, the spectroscopic properties obtained from the CHIPR HO+ 2 PES have been shown to agree with the available experimental data within the assigned errors. The resulting model potential energy function may thefore be commended for undertaking dynamical studies of the relevant molecular reactive scattering problems. Acknowledgement This work has the support of Fundação para a Ciência e a Tecnologia,


10

Page 11 of 15

10

r2 /a0

8

Table 4: Vibrational levels of the CHIPR PES; the quantum numbers n1 , n2 and n3 denote symmetric stretching, bending and asymetric stretching, respectively. n2 0 1 0 2 2 0 0 3 2 1 1 0 1 4 0 4 2 2 5 2 3 2 1 2 4 2 1 3 0 3

n3 0 0 1 0 1 2 0 0 1 3 0 3 1 1 1 2 4 0 1 3 2 2 4 1 1 0 2 3 0 0

2

Energy(cm−1 ) 0 1089 1399 2145 2468 2769 3054 3171 3512 3817 3955 4111 4117 4168 4451 4535 4851 5011 5128 5140 5141 5334 5426 5499 5531 5636 5816 5868 5929 6038

2

3

4

5 6 r1 /a0

7

8

9

2

3

4

5 6 r1 /a0

7

8

9

2

3

4

5 6 r1 /a0

7

8

9

10 8

r2 /a0

n1 0 0 0 0 0 0 1 1 1 0 1 0 1 0 1 0 0 2 1 0 0 0 0 1 2 0 1 0 2 2

6 4

6 4 2

10 8

r2 /a0

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6 4 2

Figure 5: Contour plot of ground state HO2+ PES at a fixed valence angle of α = 0◦ , 90◦ , 180◦ to show the two distinct dissociation channels. The x and y axes refer to O— O and O—H distances. The contour lines are equally spaced with an increment of 0.077Eh ACS Paragon Plus Environment

11


Page 12 of 15

Table 5: Structural properties of idendified isomers in ground state CHIPR PES of HO2+. PointGroup Cs C2v D2h C2v C2v

r(O—O) (Å) model ab initio 1.23 1.23 3.22 3.00 2.42 2.42 1.12 1.12 1.27 1.27

r(O—H) (Å) model ab initio 1.01 1.01 1.02 1.03 1.21 1.21 2.80 2.82 1.26 1.26

α(deg) model ab initio 114 112 0 0 0 0 79 79 59 60

∆(eV) model ab initio 0 0 4.48 4.46 5.05 5.14 2.81 2.81 1.68 1.68

Table 6: Comparison of molecular properties from CHIPR PES with other theoretical work. Properties r(O—O)(Å) r(O—H)(Å) α(deg) A (cm−1 ) B (cm−1 ) C (cm−1 ) ν1 (cm−1 ) ν2 (cm−1 ) ν3 (cm−1 )

Weaver et al. 3 1.2272 1.0102 112.743 21.9919 1.27902 1.19699 3028 1440 1068

Huang and Lee 9 1.23192 1.00973 112.461 21.5086 1.2777 1.20612 3230.4 1409.3 1096.7

Robbe et al. 12 1.237 1.007 118.8 21.5040 1.2698 1.1990 3268.25 1392.45 1113.65

this work 1.2304 1.0065 113.884 22.1963 1.2676 1.19911 3054 1399 1089

Portugal, via project UID/QUI/00313/2019. The support to A.J.C.V. from China’s Shandong Province “Double-Hundred Talent Plan” (2018) is also gratefully acknowledged.

Detection of HO2 in an Atmospheric Pressure Plasma Jet using Optical Feedback Cavity-Enhanced Absorption Spectroscopy. New J. Phys. 2016, 18, 1–9.

Supporting Information Available: Included are: 1) Analytic function used in diatomic fit of OH+; 2) Analytic function used in diatomic fit of O2; 3) Analytic function used in triatomic fit of HO2+; 4) Vibrational wave functions for the first 30 states of HO2+ in Figure S1. This material is available free of charge via the Internet at http://pubs.acs.org/.

(3) Weaver, S. L. W.; Woon, D. E.; Ruscic, B.; McCall, B. J. Is HO2+ a Detectable Interstellar Molecule? Astrophys. J 2009, 697, 601–609. (4) Kohguchi, H.; Jusko, P.; Yamada, K. M. T.; Schlemmer, S.; Asvany, O. HighResolution Infrared Spectroscopy of O2H+ in a Cryogenic Ion Trap. J. Chem. Phys. 2018, 148, 144303(1)–144303(10).

References

(5) Dyke, J. M.; Jonathan, N. B. H.; Morris, A.; Winter, M. J. Vacuum Ultraviolet Photoelectron Spectroscopy of Transient Species. Mol. Phys. 1981, 44, 1059–1066.

(1) Parise, B.; Bergman, P.; Du, F. Detection of the Hydroperoxyl Radical HO2 toward ρ Ophiuchi A Addtional Constraints on the Water Chemical Network. A&A 2012, 541, 1–4.

(6) W.C.Ho,; C.J.Pursell,; Oka, T. Infrared Spectroscopy in an H2O2+He Discharge: H3O+. J. Mol. Spectrosc. 1991, 149, 530– 541.

(2) Gianella, M.; Reuter, S.; Aguila, A. L.; Ritchie, G. A. D.; van Helden, J.-P. H.


12

Page 13 of 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(7) Nizkorodov, S. A.; Roth, D.; Olkhov, R. V.; Maier, J. P.; Dopfer, O. Infrared Predissociation Spectra of HeHO2+ and Ne-HO2+: Prediction of the ν1 Frequency of HO2+. Chem. Phys. Lett. 1997, 278, 26–30.

(16) Varandas, A. J. C. Accurate CombinedHyperbolic-Inverse-Power-Representation of Ab initio Potential Energy Surface for the Hydroperoxyl Radical and Dynamics Study of O + OH Reaction. J. Chem. Phys. 2013, 138, 134117(1)–134117(7).

(8) Jacox, M. E.; Thompson, W. E. Infrared Spectra of Products of the Reaction of H Atoms with O2 Trapped in Solid Neon: HO2, HO2+, HOHOH–, and H2O(HO). Phys. Chem. A 2013, 117, 9380–9390.

(17) Varandas, A. J. C. Combined-HyperbolicInverse-Power-Representation of Potential Energy Surfaces: A Preliminary Assessment for H3 and HO2. J. Chem. Phys. 2013, 138, 054120(1)–054120(14).

(9) Huang, X.; Lee, T. J. A Procedure for Computing Accurate Ab initio Quartic Force Fields: Application to HO2+ and H2O. J. Chem. Phys. 2008, 129, 044312– 044312.

(18) Varandas, A. J. C. Reaction Rate Constant Computations: Theories and Applications, 6th ed.; Royal Society of Chemistry: Thomas Graham House, Milton Road, United Kingdom, 2013.

(10) Vazquez, G. J.; Buenker, R. J.; Peyerimhoff, S. D. The Electronic Structure of HO2+. Mol. Phys. 1986, 59, 291–316.

(19) Varandas, A. J. C. Recalibration of a Single-Valued Double Many-Body Expansion Potential Energy Surface for Ground State HO2 and Dynamics Calculations for O2 + H Reaction. J. the O + OH Phys. Chem. 1990, 94, 8073–8080.

(11) Grimbert, D.; Lassier-govers, B.; Sidis, V. Model Potentials and Related Diabatic States for the H++O2 Collisional System. Chem. Phys. 1988, 124, 187–204.

(20) Murrell, J. N.; Carter, S.; Farantos, S. C.; Huxley, P.; Varandas, A. J. C. Molecular Potential Energy Functions, 1st ed.; Wiley: Chichester, United Kingdom, 1984.

(12) Robbe, J. M.; Monnerville, M.; Chambaud, G.; Rosmus, P.; Knowles, P. J. Theoretical Spectroscopic Data of the HO2+ ion. Chem. Phys. 2000, 252, 9–16.

(21) Varandas, A. J. C.; Murrell, J. N. A ManyBody Expansion of Polyatomic Potential Energy Surfaces: Application to Hn Systems. Faraday Discuss. Chem. Soc. 1977, 62, 92–109.

(13) F. George D. Xavier, Nonadiabatic Dynamics on the Two Coupled Electronic PESs : The H+ + O2 System. J. Phys. Chem. A 2010, 114, 10357-10366.

(22) Karton, A.; Martin, J. M. L. Comment on: "Estimating the Hartree-Fock Limit from Finite Basis Set Calculations"[Jensen F (2005) Theor Chem Acc 113:267]. Theor. Chem. Acc. 2006, 115, 330–333.

(14) F. George D. Xavier,; Kumar, S. Ab initio Adiabatic and Quasidiabatic Potential Energy Surfaces of the Lowest Four Electronic States of the H+ + O2 System. J. Chem. Phys. 2010, 133, 164304(1)– 164304(12).

(23) Varandas, A. J. C. Extrapolating to One Electron Basis Set Limit in Electronic Structure Calculations. J. Chem. Phys. 2007, 126, 244105(1)–244105(15).

(15) F. George D. Xavier,; MartinezGonzález, M.; Varandas, A. J. C. Accurate Ab initio Potential for HO2+: CBS Extrapolated Energies and DirectFit Diatomic Curves. Chem. Phys. Lett. 2018, 691, 421–430.

(24) Merer, A. J.; Malm, D. N.; Martin, R. W. The UltraViolet Emission Spectra of OH+


13


and OD+ . Rotational Structure and Perturbations in the A3 Π − X3 Σ− Transition. Can. J. Phys. 1975, 53, 251–283.

Page 14 of 15

(34) Varandas, A. J. C. Extrapolating to One Electron Basis Set Limit in Electronic Structure Calculations. J. Chem. Phys. 2007, 126, 244105(1)–244105(15).

(25) Spirko, J. A.; Mallis, J.; Hickman, A. P. Calculation of Adiabatic and Diabatic 3 Σ− State of OH+ . J. Phys. B: Mol. Opt. Phys. 2000, 33, 2395–2407.

(35) Litorja, M.; Ruscic, B. A Photoionization Study of the Hydroperoxyl Radical, HO2, and Hydrogen Peroxide, H2O2. J. Electron Spectrosc. Relat. Phenom. 1998, 97, 131– 146.

(26) Huber, K. P.; Herzberg, G. Molecular Spectra and Molecular Structure; Van Nostrand Reinhold Company: New York, 1979.

(36) Ádám Ganyecz,; Csontos, J.; Nagy, B.; Kállay, M. Theoretical and Thermochemical Network Approaches To Determine the Heats of Formation for HO2 and Its Ionic Counterparts. J. Phys. Chem. A 2015, 119, 1164–1176.

(27) Irikura, K. K. Experimental Vibrational Zero-Point Energies : Diatomic Molecules. J. Phys. Chem. Ref. Data 2007, 36, 389– 397.

(37) Ruscic, B.; Pinzon, R. E.; Morton, M. L.; Srinivasan, N. K.; Su, M.-C.; Sutherland, J. W.; Michael, J. V. Active Thermochemical Tables: Accurate Enthalpy of Formation of Hydroperoxyl Radical, HO2. J. Phys. Chem. A 2006, 110, 6592–6601.

(28) MOLPRO, Version 2012.1, a Package of Ab Initio Programs. (29) Galvão, B. R. L.; Mota, V. C.; Varandas, A. J. C. Modeling Cusps in Adiabatic Potential Energy Surfaces. J. Phys. Chem. A 2015, 119, 1415–1421. (30) Galvão, B. R. L.; Mota, V. C.; Varandas, A. J. C. Modeling Cusps in Adiabatic Potential Energy Surfaces Using a Generalized Jahn-Teller Coordinate. Chem. Phys. Lett. 2016, 660, 55–59. (31) Rocha, C. M. R.; Varandas, A. J. C. Multiple Conical Intersections in Small Linear Parameter Jahn-Teller Systems: The DMBE Potential Energy Surface of Ground-State C3 Revisited. Phys. Chem. Chem. Phys. 2018, 20, 10319–10331. (32) Varandas, A. J. C. A Useful Triangular Plot of Triatomic Potential Energy Surfaces. Chem. Phys. Lett. 1987, 138, 455– 461. (33) Tennyson, J.; Kostin, M. A.; Barletta, P.; Harris, G. J.; Polyansky, O. L.; Ramanlal, J.; Zobov, N. F. DVR3D: A Program Suite for the Calculation of RotationVibration Spectra of Triatomic Molecules. Comput. Phys. Commun. 2004, 163, 85– 116.


14

Page 15 of 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry a

c

b

d

e

Minima on the ground potential energy surface of HO2+


15

Global Potential Energy Surface for HO2+ Using the CHIPR Method

Recommend Documents