Ab Initio Vibrational Spectroscopy of cis- and trans ... - ACS Publications

Nov 16, 2016 - Ab Initio Vibrational Spectroscopy of cis- and trans-Formic Acid from a Global Potential Energy Surface. David P. Tew* and Wataru Mizuk...
1 downloads 0 Views 3MB Size
Subscriber access provided by United Arab Emirates University | Libraries Deanship

Article

Ab Initio Vibrational Spectroscopy of Cis- and TransFormic Acid From a Global Potential Energy Surface David P. Tew, and Wataru Mizukami J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b09952 • Publication Date (Web): 16 Nov 2016 Downloaded from http://pubs.acs.org on November 26, 2016

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry A 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 16

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

Ab initio Vibrational Spectroscopy of cis- and trans-Formic Acid from a Global Potential Energy Surface David P. Tew⇤ and Wataru Mizukami

School of Chemistry, University of Bristol, Bristol BS8 1TS, United Kingdom⇤ (Dated: November 16, 2016) A global analytic potential energy surface for the ro-vibrational dynamics of cis- and trans-formic acid is presented, constructed using LASSO-based regression to reproduce CCSD(T)(F12*)/ccpVTZ-F12 energies. The fit is accurate to 0.25% has an RMS deviation from the ab initio data of 9 cm 1 for the energy range 0-15000 cm 1 . Converged J = 0 vibrational eigenstates are reported, computed using vibrational configuration interaction with an internal coordinate path Hamiltonian for the torsional motion connecting the cis and trans rotamers. Methodological choices concerning the appropriate definitions of the curvilinear and diabatic bath coordinates are discussed. The zero point of the cis rotamer is 1412 cm 1 above that of the trans, which lies at 7354 cm 1 . The computed fundamentals match the bands recorded from gas-phase IR spectroscopy with an RMSD of only 3 cm 1 . A fresh assignment of the overtone spectra of both the cis and trans rotamers is presented for the energy range 0–4720 cm 1 , where 14 out of the 51 bands are reassigned on the basis of the VCI calculations. I.

INTRODUCTION

Formic acid has attracted many spectroscopic investigations since it is relevant to atmospheric chemistry, being present throughout the Earth’s troposphere,1 and is assumed to play an important role in interstellar chemistry as a possible building block for other organic compounds.2–4 Several microwave5–7 and submillimeter3,8,9 rotational spectra of trans-formic acid have been recorded, from which the structure and molecular constants have been determined as well as and effective Hamiltonian parameters and hyperfine structure. Gas-phase and matrix isolated Raman10,11 and infra red12–20 spectra have also been recorded and fundamental, overtone and combination bands listed in the mid-IR and near-IR regions, as well as high-resolution work on the low frequency ⌫7 and ⌫9 bands relevant to astrophysical detection. The trans rotamer is energetically favoured over the cis rotamer and the population at ambient conditions is approximately 1000:1. The cis rotamer was first detected by microwave spectroscopy,21 and has since been the subject of gas-phase sub-millimeter/IR3,22,23 and matrix isolated IR spectroscopy.15,24–29 From a theoretical perspective, however, far fewer investigations have been performed.15,30–34 The overtone spectra of trans-formic acid was assigned with the aid of correlation-corrected vibrational self-consistent-field calculations on a Møller–Plesset potential surface.15 Highlevel coupled cluster theory has been used to compute the equilibrium structure and semi-diagonal quartic force fields for both rotamers,9,31 which have been used in combination with vibrational perturbation theory to compute vibrationally averaged properties9 and to reassess the effective Hamiltonians used to model the ro-vibrational spectra.31 This article describes a highly accurate and comprehensive unified theoretical treatment of the midIR overtone spectrum of cis- and trans-formic acid for the purpose of a full assignment the observed bands and to provide the community with a high quality global po-

tential energy surface. The work reported builds upon a preliminary publication35 where we applied a newly developed theoretical approaches to compute the fundamental bands of trans-formic acid. For a theoretical prediction to be accurate, the quantum mechanical treatment of both the electronic and nuclear degrees of freedom must be sufficiently complete in both the level of many-body correlation included and the detail in which this is represented. For cases such as formic acid, where the ground electronic state is well separated in energy from high-lying states, accurate predictions can be obtained by first computing the Born– Oppenheimer electronic potential energy hypersurface36 at a high level of electronic structure theory, and then computing the ro-vibrational energy levels on that surface at a high-level of ro-vibrational structure theory. An accurate potential energy surface is provided by near basis set limit fc-CCSD(T) calculations37 (frozen-core coupled cluster singles, doubles with a perturbative treatment of triples), since the contributions from core-valence and post-(T) correlation, relativistic and non-adiabatic e↵ects largely cancel.38,39 With the advent of mature F12 explicitly correlated techniques,40,41 near basis set limit energies can be computed using small basis sets42–44 and the tens of thousands of energy evaluations required to map out high-dimensional surfaces are now routinely possible.45,46 Concurrent developments in variational ro-vibrational wavefunction methods have made it possible to compute the ro-vibrational energy levels of polyatomic molecules with up to 10 atoms and floppy degrees of freedom.47–56 Developments include highly efficient vibrational configuration interaction contraction schemes,57,58 scalable approximations to the n-body Hamiltonian matrix elements59–62 and general kinetic energy operators using curvilinear coordinate systems.48,63 For molecules with one or two large amplitude motions, an attractively simple representation is provided by generalisations of the semi-rigid-bender Hamiltonian approach of Hougen,

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 2 of 16 2

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

FIG. 1. Potential energy curve and ⌫9 torsional states of trans and cis formic acid.

Bunker and Johns (HBJ)64 which combine one or two curvilinear coordinates with a normal coordinate treatment of semi-rigid motions. Such approaches include the internal coordinate path,64 intrinsic reaction path,65 and reaction surface66,67 Hamiltonians and have been applied to the variational calculation of ro-vibrational energy levels of water,64 hydrogen peroxide,68,69 glyoxal,70,71 ammonia,72 methanol69,73 and malonaldehyde.67,73 In this work we present a highly accurate fulldimensional analytic potential energy surface for formic acid, computed by fitting to near basis set limit CCSD(T)(F12*) energies,44 using our recently developed fitting procedure35,74 based on the LASSO approach to linear regression with compact functional support.75 We present results of high-level VCI calculations on the J = 0 energy levels, using an incrementally expanded59,60 (multimode) internal coordinate path Hamiltonian64 with up to 5-mode coupling, and compare our predicted transition energies to those observed in the IR spectra of cisand trans-formic acid. The J=0 vibrational energy levels of formic acid have posed a significant challenge to ab initio simulation. The motion of this molecule includes a internal rotation of the OH moiety that connects the cis and trans rotamers (⌫9 ) (see Fig. 1), as well as a Fermi resonance between 2⌫9 and ⌫5 , the OH internal rotation and the C-O-H bending motion. In an earlier preliminary study we used second order vibrational active space perturbation theory, VASPT2,35 to compute the energy levels and found that the strength of the Fermi resonance and the position of the 2⌫9 and ⌫5 levels were highly sensitive to the number of basis functions in ⌫9 and could not be fully converged. This was because the calculations used the normal coordinates of the trans rotamer to represent the wavefunction, which cannot not describe the motions at the cis

rotamer without including a very high order of modecoupling away from the trans minimum. In this work, we report high-level VCI calculations using an internal coordinate path Hamiltonian, which is particularly well suited to the present case of one curvilinear coordinate. We show that the traditional choice of body-fixed axes for the torsional motion, as well as the choice of adiabatic perpendicular coordinates are in fact inappropriate and present refined choices that lead to a more rigorous theoretical treatment. Our previous calculations on formic acid were performed using a full dimensional analytic potential energy surface, which we constructed using distributed multivariate Gaussian functions, fit to reproduce CCSD(T)(F12*)/cc-pVDZ-F12 energies using a method similar to LASSO constrained optimisation.74,75 Although the potential reproduced the ab initio data to a root mean squared deviation (RMSD) of 0.25% over the energy range 0–15000 cm 1 , it was biased to reproduce the lower energy trans rotamer. In this work we construct an improved surface that is high quality for both the cis and trans rotamers and is fit to a new data set of 17076 CCSD(T)(F12*)/cc-pVTZ-F12 energies. This paper is organised as follows: The methods used to generate data and fit the electronic potential energy surface are described in section II, together with a discussion of the quality of the surface. The methods used to compute the quantum eigenstates of the ro-vibrational nuclear motion on the surface are described in section III. The results of the calculations are presented in section V, where computed fundamentals, overtones and combination bands of both trans- and cis-formic acid are listed together with a fresh assignment of the experimentally observed bands.

II.

THE POTENTIAL ENERGY SURFACE

The availability of an analytic potential energy surface greatly facilitates the quantum treatments of the ro-vibrational motion of molecules, since the number of ab initio electronic energy evaluations required to train a model potential is orders of magnitude smaller than the number of grid points required in a fully converged quantum dynamics treatment,76 and the size and nature of the quadrature grid can be increased or adapted to obtain convergence without repeating costly electronic structure calculations. Possible strategies for constructing analytic potential surfaces include interpolation schemes,77–83 linear regression on a selected function space52,84–89 and non-linear optimisation of generating functions.90–96 We have recently developed a new approach35,74 that combines a simple potential of Morse oscillators, that have qualitatively correct asymptotic behaviour, with a highly flexible, short-ranged and minimally perturbing correction surface. The use of elastic net regression, rather than plain least squares regression, is key to the success of our approach, since it both constrains the magnitude of

ACS Paragon Plus Environment

Page 3 of 16

The Journal of Physical Chemistry 3 20

200

18

150

a)

16

1.5 b)

1

c)

100

12 10 8 6

0.5

50

% Error

14 Error in cm-1

RMSD in cm-1

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 -50

0 -0.5

-100

4

-1

-150

2 0

-200 0

2500 5000 7500 10000 12500

-1.5 0

2500 5000 7500 10000 12500

2500 5000 7500 10000 12500

Energy in cm-1

Energy in cm-1

Energy in cm-1

FIG. 2. a) Cumulative RMSD error of the model PES. b) Scatter plot of errors. b) Scatter plot of percentage errors. Blue points are from the training set and red points are from the control set.

the correction surface and minimises the number of functions that contribute to the surface. We have applied our new approach to fit a model potential energy surface for the vibrational and torsional motions of formic acid. We have computed 17076 CCSD(T)(F12*) energies at random structures in the energy window 0-15000 cm 1 relative to the trans equilibrium. We used the cc-pVTZF12 basis set,97 the frozen-core approximation and a geminal exponent of 1.0 a0 1 .98 For the density fitting approximation used to accelerate the CCSD(T)(F12*) calculation, the aug-cc-VTZ fitting basis sets99 were employed for the MP2 and Fock terms. For the complementary auxiliary basis required for the F12 treatment, the OPT-RI cc-pVTZ-F12 basis set was used.100 The CCSD(T)(F12*) method has been shown to have the best cost-accuracy ratio among the available approximate CCSD(T)-F12 methods.101 The 17076 structures were generated by adding to the existing set of 7420 by random displacements from the cis and trans equilibria, as well as from the transition state for the OH torsional motion, using the same approach as applied in Ref. 35. This was further augmented with points generated from normal mode displacements away from the torsional path, using the internal coordinate path coordinate system (vide infra). The form of the analytic potential fit to the data is V (r) = V0 (r) + Vc (r) ⇣ X V0 (r) = Db2 1 e

↵b (rb eb )

b

Vc (r) =

XX i

+

dµi (e

µ i (ri

2 cµ i )

d⌫ij (e

⌫ i (ri

2 c⌫ i)

⌘2

(1) (2) sµi )

µi

XX

⌫ j (rj

2 c⌫ j)

s⌫ij )

i>j ⌫ij

+ ···

data with the primary role of ensuring that the potential behaves properly in the asymptotic regions of dissociation and atom-atom coalescence. The set b is chosen to contain all atom pairs involved in the stretching and bending internal coordinates, that is, all distances apart from rH2 H5 and rO3 H5 , where the atom numbering is according to Fig 1. The correction surface is a sum of distributed multivariate Gaussian functions of combinations of atom-atom separations ri , where 1d, 2d, 3d and 4d Gaussian functions distributed on 1d, 2d, 3d and 4d grids are included in the fit. These local functions only occupy the ranges of configuration space spanned by the ab initio data and Vc does not a↵ect the asymptotic behaviour of V . In Eq. (3) the index µ labels a Gaussian centred at cµi with width parameter iµ . The corresponding shift parameter sµi is fixed by the standardisation requirement. The linear parameters d of the correction surface Vc are determined by a method similar to LASSO regression, where the R2 quality P of the fit is optimised subject to the constraint that  |d | = . This imposes that d = 0 for all but a relatively small subset of functions and it is thus possible to use a very large basis of functions without risk of overfitting. The LASSO procedure provides a family of solutions with favourable ratios of number of functions to quality of fit, automatically selecting which functions are most needed to fit the data. The total number of functions in the fit set for the LASSO procedure was 26 1d, 706 2d, 5642 3d and 13956 4d Gaussians, which were generated by placing functions on sets of grids with spacings in the range 0.1–1.0 a0 . The width parameters of the Gaussians were selected by setting the standard deviation of the Gaussians to twice the spacing between the centres. Each ab initio data point is assigned a weight wi according to Ei , it’s energy in wavenumbers above the trans minima, with the weight function

(3) wi =

The zeroth order surface V0 is a sum of Morse functions of atom-atom separations rb that are loosely fit to the

100 Ei + 100

(4)

This ensures that the cumulative RMS error grows ap-

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 4 of 16 4

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

TABLE I. Harmonic wavenumbers of trans and cis formic acid.

ZP A0 O-H str. C-H str. C=O str. H-C-O bend C-O-H bend C-O str. O-C-O bend A00 C-H wag O-H tor. a

trans PES CCSD(T)a 7456 7454

cis PES CCSD(T)a 7386 7379

⌫1 ⌫2 ⌫3 ⌫4 ⌫5 ⌫6 ⌫7

3767 3092 1818 1412 1323 1140 632

3766 3093 1818 1410 1319 1140 633

3829 3007 1861 1428 1299 1124 664

3830 3009 1859 1422 1295 1123 663

⌫8 ⌫9

1056 673

1057 672

1038 522

1035 521

fc-CCSD(T)(F12*)/cc-pVTZ-F12 calculations

proximately linearly as a function of energy above the reference. It should be noted that the actual weight of a particular region of configuration space in the fit depends both on the assigned weight and an intrinsic weight arising from the density of data points in that region. The method for data point selection was chosen such that the density loosely approximates a Boltzmann distribution with kT ⇠ 12 !1 ⇠ 2000 cm 1 , where !1 is the highest harmonic wavenumber. 90% of the data set was used to train the potential form and 10% used as a control set to assess the quality of the fit. We found that the LASSO optimisation could be continued until the RMS error fell below 10 cm 1 and that the error distribution of the control set remained within 1 cm 1 of that of the training set in the range 0–7500 cm 1 . This resulted in 810 Gaussian functions being included in the final potential. Plots of the fitting error against the energy above the equilibrium structure are displayed in Fig. 2. Panel a) contains the cumulative RMS deviation, and scatter plots of the error distributions are displayed in panels b) and c) for total and percentage errors, respectively. In all plots, data for the training set are displayed in blue and the data for the control set in red. The 810 Gaussian functions in the correction surface are sufficient to reduce the percentage error to below 0.25% for all but 633 data points in the training set and 140 in the control set. The harmonic wavenumbers at the cis and trans minima for the model potential are compared to the corresponding data for the ab initio potential in Table I. The RMS deviation is 3 cm 1 . The structures on the model potential are compared to the ab initio values in Table II and agree to within 0.1 pm for bond lengths and 0.1 degrees for angles, which is smaller than the deviation between the ab initio values and the semi-experimental re values from Ref. 31, where ro-vibrational perturbation theory was used to remove zero-point e↵ects. The energies of the cis minima and torsional barrier are 1488 and 4422 cm 1 on the model potential compared to 1488 and 4423 cm 1 computed ab initio. Benchmark studies examining

TABLE II. Geometric parameters of trans and cis formic acid in pm and degrees. trans C=0 C O C H O H O C=O H C=O H C O C O H

Expt.a 120.3(5) 134.2(5) 109.7(5) 97.2(5) 124.8(4) 123.2(15) 112.0(15) 106.3(4)

re b 119.8 134.1 109.2 96.6 124.8 125.0 110.1 107.0

PES 119.9 134.2 109.4 96.7 124.9 125.1 110.0 106.7

CCSD(T)c 119.9 134.2 109.4 96.7 124.9 125.1 110.0 106.7

cis C=0 C O C H O H O C=O H C=O H C O C O H

Expt.d 119.5(3) 135.2(3) 110.5(4) 95.6(5) 122.1(4) 123.2(6) 114.6(6) 109.7(4)

re b 119.2 134.7 109.8 96.1 122.3 123.3 114.5 109.3

PES 119.2 134.9 110.0 96.2 122.3 124.0 113.7 109.2

CCSD(T)c 119.2 134.9 110.0 96.2 122.3 124.0 113.7 109.2

a b c d

rs structure Ref. 7 with uncertainties in parentheses Semi-experimental re values from Ref. 31 re from fc-CCSD(T)(F12*)/cc-pVTZ-F12 calculations rs structure Ref. ? with uncertainties in parentheses

ab initio contributions beyond basis set limit, frozen core CCSD(T) theory indicate that the ab initio data used for our fit has an error bar of approximately 0.2 pm for structures39 and 8 cm 1 for harmonic frequencies.38 The quality of our potential is such that the fitting errors are smaller than the intrinsic error of the data fit to. III.

THE INTERNAL COORDINATE PATH HAMILTONIAN

In this section we summarise the internal coordinate path Hamiltonian approach and describe the details of our implementation in the DYNAMOL package. In addition, we elucidate the solutions to two problems previously encountered in HBJ-based treatments, namely how decide on the appropriate definition of the body-fixed axes for the curvilinear path and how to appropriately determine the normal coordinates as a function of the path coordinate. A.

The path coordinate and path-dependent orthogonal normal coordinates

A large amplitude motion is characterised by a curvilinear path in mass-weighted Cartesian space ai↵ (s), parameterised by the arc length s, where a0i↵ a0i↵ = 1. The compound index i↵ refers to the mass-weighted Cartesian coordinate ↵ of the ith particle in the molecule. Here and throughout we use a prime to denote derivative with respect to arc length and Einstein summation convention

ACS Paragon Plus Environment

Page 5 of 16

The Journal of Physical Chemistry 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

for repeated indices. The remaining 3N 7 vibrational degrees of freedom are characterised by displacements Qk in normal coordinates with vectors li↵,k (s), defined to be orthogonal to each other and to the curvilinear path coordinate, and to be orthogonal to the vectors of instantaneous translation and rotation of the molecule in configuration ai↵ (s). The large amplitude motion path could be, for example, the intrinsic reaction path connecting transition states with minima, or an internal coordinate path obtained from constrained minimisation. These paths are usually rather similar and the latter is much more straightforward to obtain and is favoured in this work. Ab initio electronic structure calculations provide a series of points ai↵ (⌧x ). These nodes must first be oriented relative to a common body-fixed axis and an interpolation procedure must then be established to enable conversion from ⌧ to arc length s and to define the path ai↵ (s) for all values of s. Once the path is defined, the vector space of the orthonormal vibrational coordinates perpendicular to the path can then be established and a smoothly varying coordinate set li↵,k defined. In the following we discuss the formal aspects of the various choices involved in defining the internal coordinate path coordinate system and a practical implementation is presented in section III B Each geometry of the curvilinear path is placed with its centre of mass at the origin of the body-fixed axis frame so that motion in the path is orthogonal to translational motion of the molecule in the usual way. Three possibilities for orienting the nodes ai↵ (⌧x ) in the body-fixed axis present themselves: i) apply the Eckart condition ✏↵ a0i ai = 0 so that motion along the path a0i↵ (⌧x ) is orthogonal to the instantaneous rotation of the molecular geometry ai↵ (⌧x ) ii) apply the Eckart condition ✏↵ a0i a0i = 0 so that motion in the path is orthogonal to the instantaneous rotation of a special geometry a0i↵ iii) rotate each path geometry onto its principal axes. Here and throughout ✏↵ is the antisymmetric tensor. The three options di↵er in the treatment of the Coriolis coupling between motion in the path and rotation of the body-fixed axes. Option i) defines a path with zero angular momentum and this Coriolis coupling contribution is eliminated. Option i) was used in the original HBJ and subsequent RPH treatments. However, when the path circumscribes a torsional motion, the zero-angular momentum definition induces a rotation of the principal axes and violates the requirement that ai↵ (⌧ + 2⇡) 6= ai↵ (⌧ ). For torsional coordinates, either ii) or iii) must be used. In this work we select option ii). The coordinates perpendicular to the path must also be defined to complete the specification of the coordinate system. For cases where the path motion is an order of magnitude slower than the perpendicular motions, it

is appropriate to select adiabatic perpendicular coordinates. At each point along the path, the set of 3N 7 adiabatic perpendicular coordinates li↵,k (s) are the principal eigenvectors of the projected Hessian HP HP = (1

P(s))H(s)(1

P(s))

(5)

where H(s) is the mass weighted second derivative of the potential at that point on the path and P(s) projects onto the vectors for instantaneous translation and rotation of that path geometry and instantaneous motion along the path. This adiabatic definition was introduced by HBJ and adopted in subsequent reaction path treatments and the character of the perpendicular modes often change significantly as a function of the path coordinate, particularly at avoided crossings. In most cases encountered in molecular dynamics, the timescale for the path motion is commensurate with the perpendicular motions and a diabatic coordinate choice is more appropriate.102 Diabatic perpendicular coordinates can be defined by minimising l0k · lk , which minimises the Coriolis coupling constant associated with the motion along the path (vide infra). A practical implementation of this requirement is outlined in section III B.

B.

Implementation

Having obtained a set of geometries ai↵ (⌧x ) that lie on path by constrained geometry optimisation, fixing the OC-O-H dihedral angle ⌧ , these geometries are straightforwardly placed in the appropriate Eckart body-fixed axis frame by translation to the centre of mass and rotation to satisfy ✏↵ ai a0i = 0, which is an orthogonal procrustes problem and solved through SVD.103 The condition that ✏↵ a0i a0i = 0 follows automatically, as revealed by differentiation. The special reference geometry a0i↵ is chosen to be the equilibrium. Interpolation between the nodes to obtain a continuous path is performed by means of periodic cubic spline on each of the Cartesian components ai↵ independently. This interpolation approach has the property that geometries at interpolated values of ⌧ obey both the centre of mass and Eckart conditions exactly.104 60 nodes were sufficient for the interpolated geometries at the midpoints between the nodes to have a dihedral within 10 5 radians of the interpolated ⌧ . The vectors dai↵ /d⌧ and thus a0i↵ and d⌧ /ds are analytically available from the cubic spline, from which arc length as a function of ⌧ can be obtained to high precision by numerical integration. At each node, the projected hessian HP is then formed where the null-space consists of the translation and rotation of the Eckart frame and motion along the path. The eigenvectors of HP at the equilibrium structure are in fact the normal modes at that point. Diabatic vectors along the path are defined by minimising l0k (⌧ ) · lk (⌧ ), which is achieved approximately in practice. The vectors at the two nodes neighbouring the equilibrium are obtained by maximising the overlap between

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 6 of 16 6

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

li↵,k (⌧x ) and li↵,k (⌧0 ) for the set of vectors with non-zero eigenvalues. The vectors at the next nodes along are determined using their neighbours and so on along the path. Finding an orthonormal vector set with maximal overlap with a reference is another orthogonal Procrustes problem, solved by SVD. Normal coordinate vectors at values of ⌧ between the nodes are determined by periodic cubic spline interpolation in the same way as for the geometry. We find that 60 nodes are sufficient for the interpolated vectors at the midpoints to be orthogonal to better than 0.001. Evaluation of the potential energy matrix elements requires conversion between the path coordinates and the atom-atom separations used for our representation of the potential energy surface. The Cartesian coordinates for every vibrational configuration (s, Q1 , . . . , Q3N 7 ) are given by ri↵ = ai↵ (s) + li↵,k (s)Qk , from which the atomatom distances required for the potential energy function are readily obtained.

C.

Kinetic energy operator

The kinetic energy operator for Internal Coordinate Path Hamiltonian is ˆd ⇡ ˆe ˆd )µde µ 1/2 (⇧ Tˆ = 12 µ1/4 (⇧ + 1 µ1/4 Pˆk µ 1/2 Pˆk µ1/4 .

⇡ ˆe )µ1/4 (6)

2

The indices d, e 2 x, y, z, s, where x, y, z are the bodyfixed Cartesian axes and s is the path coordinate. The index k lables the 3N 7 normal modes orthogonal to the path. The angular momentum, path momentum, Coriolis and normal mode momentum operators are @ ⇧ = (Jˆx , Jˆy , Jˆz , i~ @s ) ⇡ ˆd = ⇣d,kl Qk Pˆl

(8)

Pˆk =

(9)

(7)

@ i~ @Q k

and ↵ ,k (s)

= ✏↵ " ✏↵ " (li ,k (s)ai (s) + ai (s)li ,k (s)) (15) 0 0 (16) ↵s,k (s) = ✏↵ (li ,k (s)ai (s) + ai (s)li ,k (s)) ss,k (s)

= 2 li0

1 2

k Qk )I0

1

(I0 +

1 2

l Ql )

= Y I0 1 Y.

(17)

⇣↵,kl (s) = ✏↵ li ,k (s)li ,l (s) ⇣s,kl (s) = li0 ,k (s)li ,l (s)

(18) (19)

Note that the terms in Eq. (13) are zero if the Eckart condition i) is applied, but is not zero in this work. D.

Matrix elements

We require expectation values TLR = h L |Tˆ| R i and we currently are only interested in the J = 0 rotational state. Integrating by parts and expanding gives Z ˆ s µ1/4 L )(⇧ ˆ s µ1/4 R ) 2TLR = d⌧ µ 1/2 µss (⇧ Z ˆ s µ1/4 L )(ˆ + d⌧ µ 1/2 µse (⇧ ⇡e µ1/4 R ) Z ˆ s µ1/4 R ) + d⌧ µ 1/2 µds (ˆ ⇡d µ1/4 L )(⇧ Z + d⌧ µ 1/2 µde (ˆ ⇡e µ1/4 L )(ˆ ⇡e µ1/4 R ) Z + d⌧ µ 1/2 (Pˆk µ1/4 L )(Pˆk µ1/4 R ) (20) Watson showed that for Wilson’s semi-rigid rotor Hamiltonian, µ↵ [ˆ ⇡ , µ] = 0, but no such simplification has been demonstrated for semi-rigid-bender Hamiltonians of the kind considered here. However, the commutators ˆ s , µ] can be evaluated analytically. Intro[Pˆk , µ] and [⇧ ducing the notation i~µ1/4 q = [ i~

@ 1/4 ,µ ] @q

(21)

where q runs over both Qk and s, we have (10)

µ is the determinant of µ. The symmetric generalised reference inertia and distortion tensors are ✓ ◆ ✓ ◆ (I0 )↵ (I0 )↵s ↵ ,k ↵s,k I0 = , (11) k = (I0 )s (I0 )ss s ,k ss,k where Greek indices indicate x, y or z and (I0 )↵ (s) = ✏↵ " ✏↵ " ai (s)ai (s) (I0 )↵s (s) = ✏↵ ai (s)a0i (s) (I0 )ss (s) = a0i (s)a0i (s)

(s)

where ✏ is the unit antisymmetric tensor. Finally, the generalised Coriolis coupling constants are

µ is the inverse of a generalised instantaneous inertia tensor I, which, after application of the ⇣-sum rule factorises to I = (I0 +

0 ,k (s)ai

(12) (13) (14)

@Ide 1 I @q de 1 1 @Yde 1 @(I0 )de (I0 1 )de . 2 @q Yde + 4 @q 1 4

q = =

(22) (23)

where @Yde = 12 de,k @Qk @(I0 )de =0 @Qk 0 Yde = (I0 )0de +

ACS Paragon Plus Environment

(24) (25) 1 2

0 de,k Qk

(26)

Page 7 of 16

The Journal of Physical Chemistry 7

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

TABLE III. Zero-point and transition wavenumbers from SSVSCF calculations using the ICP Hamiltonian with m-body modal coupling and modal basis of Ms trigonometric path functions and Mk harmonic oscillator bath functions. Character ZP A0 O-H str. C-H str. C=O str. H-C-O b. O-H tor. C-O-H b. C-O str. O-C-O b. A00 C-H wag O-H tor.

m,(Ms =31,Mk =11) Ms , Mk ,(m=4) 2 3 4 27,7 29,9 31,11 Exp. 7434 7374 7388 7388 7388 7388 ⌫1 ⌫2 ⌫3 ⌫4 2⌫9 ⌫5 ⌫6 ⌫7

3600 2970 1805 1418 1259 1322 1144 659

3548 2914 1787 1380 1266 1296 1120 655

⌫8 ⌫9

1089 1036 650 654

3551 2918 1790 1383 1255 1305 1129 683

3553 2919 1790 1383 1305 1253 1129 679

3551 2918 1790 1383 1254 1305 1129 681

3551 2918 1790 1383 1255 1305 1129 683

3571 2942 1777 1380 1306 1223 1105 626

1040 1040 1040 1040 1034 648 649 649 648 641

with (I0 )0↵ (I0 )0↵s (I0 )0ss

= = =

✏↵ " ✏↵ " (a0i ✏↵ ai a00i 2a00i a0i

0 ↵ ,k

= ✏↵ " ✏↵ " (li0

0 ↵s,k 0 ss,k

= ✏↵ (ai li00 ,k + li =

2 (li00 ,k a0i

,k ai

+

+ ai li0

00 ,k ai li0 ,k a00i )

ai +

ai a0i

)

(27) (28) (29)

,k

+ li

0 ,k ai

+ a0i li

)

,k ) (30)

(31) (32)

Since both a(s) and lk (s) are cubic spline functions in our scheme, first and second order derivatives are analytic. No numerical di↵erentiation is required. The final working equations are TLR =

~2 0 C h L| Ri 2 ~2 ⇥ @ + Cq1 h @q L| Ri + h 2 ~2 2 @ @ + Cpq h @p L | @q Ri 2

@ L | @q

Ri



Cp1 2 Cpq dp

IV.

dp µde eq q

=

dp µde eq q

=

dp µde eq

=

pl ⇣d,kl Qk

+

+

Character ZP A0 O-H str. C-H str. C=O str. H-C-O b. O-H tor. C-O-H b. C-O str. O-C-O b. A00 C-H wag O-H tor.

m,(Ns , Nk =9,9) Ns , Nk ,(m=4) 2 3 4 5,5 7,7 9,9 Exp. 7427 7365 7379 7380 7379 7379 ⌫1 ⌫2 ⌫3 ⌫4 2⌫9 ⌫5 ⌫6 ⌫7

3600 2970 1805 1420 1229 1344 1143 640

3550 2914 1788 1386 1178 1118 635

⌫8 ⌫9

1091 1036 1042 1042 1043 1042 1034 640 639 633 634 635 633 641

(33)

+ k k

(34)

kp k

(35)

kp kq

(36)

ds sp

(37)

COMPUTATIONAL DETAILS

All calculations used the ICP Hamiltonian as described above and were performed using our molecular quantum

3554 2917 1790 1372 1260 1289 1127 656

3553 2917 1790 1372 1280 1290 1127 659

3553 2917 1790 1372 1266 1289 1127 657

3553 2917 1790 1372 1266 1289 1127 657

3571 2942 1777 1380 1306 1223 1105 626

dynamics package DYNAMOL. We refer readers to Ref. 35 for specifics of the implementation of VSCF, VMP2 and VCI in DYNAMOL. Trigonometric functions were used to expand the modal functions for the torsional path coordinate and Harmonic oscillator functions were used for modal expansions of the perpendicular modes, with characteristic frequencies equal to the harmonic values at the trans minimum. For the numerical integration of the Hamiltonian matrix elements, 16 point Gauss-Hermite quadrature grids were used for each coordinate with harmonic oscillator functions and a 34 point Fourier grid for the path coordinate. The parameters that define the basis used for the wavefunction expansion are Mk , the number of primitive functions used to represent the VSCF modal k and Nk , the number of VSCF modals k used for the correlation treatment. In addition to this, the Hamiltonian is approximated through a multi-mode decomposition, truncated at m-body modal interactions. V.

with C 0 = p

TABLE IV. Zero-point and transition wavenumbers from SSVMP2 calculations using the ICP Hamiltonian with m-body modal coupling and modal basis of Nt trigonometric path functions and Nh harmonic oscillator bath functions.

RESULTS

Tables III and IV display zero point, fundamental and resonance wavenumbers computed using the statespecific VSCF and VMP2 methods, respectively. These are reported to highlight the qualitative success of the mean-field model in the ICPH coordinates, but the limitation of these single-reference methods for quantitative accuracy. The VSCF and VMP2 methods are also useful for providing preliminary indications of the size of primitive and modal basis required for a converged VCI treatment. The values in Table III demonstrate the convergence of the VSCF energies with respect to the size of the primitive basis and the m-body multi-mode truncation of the Hamiltonian. Wavenumbers are listed from calculations for m=2,3,4 with (Ms , Mk )=(11,31) and from calculations for (Ms , Mk )=(7,27), (9,29) and (11,31) with m=4. We find that 12 harmonic oscillator functions

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 8 of 16 8

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

TABLE V. Zero-point and transition wavenumbers from VCI calculations using the ICP Hamiltonian with m-body modal coupling and configurations with total quanta up to Q. Character ZP A0 O-H str. C-H str. C=O str. H-C-O bend O-H tor. C-O-H bend C-O str. O-C-O bend A00 C-H wag O-H tor.

m,(Q=10) Q,(m=5) 2 3 4 5 6 8 10 Exp. Error 7428 7355 7353 7355 7355 7355 7355 ⌫1 ⌫2 ⌫3 ⌫4 2⌫9 ⌫5 ⌫6 ⌫7

3600 2970 1805 1420 1337 1236 1143 639

3559 2943 1785 1379 1314 1221 1106 629

3569 2938 1782 1376 1301 1222 1103 624

3571 2942 1777 1380 1306 1223 1105 626

4 3 6 1 1 1 3 1

⌫8 ⌫9

1091 1028 1034 1034 1036 1034 1034 1034 641 639 639 639 640 639 639 641

0 2

per mode and 31 trigonometric functions for the path coordinate is sufficient to converge the fundamentals to better than 1 cm 1 for the VSCF energies. The di↵erences between energies computed using the 3-body and 4-body Hamiltonians, however, are substantial. In an ICPH treatment the path coordinate is included in every term of the multi-mode decomposition and the 3-body Hamiltonian therefore only couples two of the perpendicular coordinates. For convergence it is necessary to extend the many-body expansion to 5-mode terms. Our implementation in DYNAMOL can treat 5-body terms in a VCI calculation, but is limited to a 4-mode expansion for VSCF and VMP2. The values in Table IV demonstrate the convergence of the VMP2 correlation treatment with respect to the size of the modal basis and the m-body multi-mode truncation of the Hamiltonian. In these calculations, the contraction from (Ms , Mk )=(11,31) to (Ns , Nk ) is performed using the VSCF ground state prior to subsequent state-specific VSCF and VMP2 calculations in the contracted modal basis. Values are presented for m=2,3,4 with (Ns , Nk )=(9,9) and for (Ns , Nk )=(5,5), (7,7) and (9,9) with m=4. For the correlation treatment it is sufficient to use 7 contracted modal functions for convergence to within 1 cm 1 . As for VSCF, the many-body expansion requires 5-mode and perhaps higher order terms for convergence. Both VSCF and VMP2 fail to describe the Fermi resonance between ⌫5 and 2⌫9 , and the energy levels change erratically as the basis changes and as the level of many-body truncation changes. In Table V we report zero point and transition energies from VCI calculations using the ICP Hamiltonian. Based on the VSCF and VMP2 results, we use a primitive basis of Ms , Mk =31,11 and a contracted modal basis of Ns , Nk =7,7. The configurations included in the VCI calculation are restricted to those where the sum of quanta of excitation are less than or equal to Q, but the number of modes simultaneously excited is otherwise unconstrained.105 We present values from calcu-

3575 2939 1783 1379 1305 1222 1108 627

3577 2940 1784 1382 1314 1234 1110 629

3575 2939 1783 1380 1306 1223 1108 627

3575 2939 1783 1379 1305 1222 1108 627

lations with Q=6,8,10 and m=5 and from calculations with m=2,3,4,5 and Q=10. With a rigorous correlation method that properly treats the resonant states, the energies are much less sensitive to the level of included mbody interaction and convergence to better than 2 cm 1 is obtained at 5-mode couplings for the zero point and better than 1 cm 1 for the fundamentals ⌫2 , ⌫3 , ⌫5 , ⌫8 and ⌫9 . For fundamentals ⌫1 , ⌫4 , ⌫6 and ⌫7 the changes from 4-body to 5-body are 6, 4, 5 and 3 cm 1 , respectively and a small contribution from 6-mode and higher couplings may reasonably be expected for these modes. The computed fundamentals are, however, in excellent agreement with the experimental values, with an RMSD of only 3 cm 1 . We proceed to compare our computed transition energies to observed frequencies of gaseous trans-formic acid and matrix isolated cis-formic acid for the fundamental, overtone and combination bands. A significant advantage of the ICPH treatment is that the nature of the vibrational energy level is easily identified from the CI coefficients, which greatly facilitates the task of assignment. Moreover, by following the stack of levels with increasing quanta of excitation in the torsional mode, we are able to identify the zero point for the cis rotamer and thus the transitions from that energy level. The calculated transitions are reported in Tables VI and VII together with available experimental data. For these final results, the number of functions Ns was increased to 9 to better converge the levels with higher quanta in the torsional mode. Of the 500 bands computed, those with up to 3 quanta of excitation and involving at most 2 modes are displayed. Some bands involving three modes or 4 quanta of excitation are displayed if they are involved in strong resonances or assignment of the experimentally observed bands. The character of the band is assigned using the CI vector. For readability only the two primary contributing states are reported and a 70% pure state is considered pure. In many cases, however, particularly at high energy, the character of the band is highly mixed in the

ACS Paragon Plus Environment

Page 9 of 16

The Journal of Physical Chemistry 9

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

CI vector. The full computed line list for the first 500 energy levels with the all the dominant states in the CI vector are provided in the SI. It was possible to identify the stack of torsional states for ⌫9 . In increasing energy (cm 1 ) relative to the trans zero point, these are 638 (⌫9 ) 1222 (2⌫9 ) 1412 (ZP cis) 1796 (3⌫9 ) 1904 (⌫9 cis) 2358 (4⌫9 ) 2400 (2⌫9 cis) 2848 (5⌫9 ). The bands belonging to the cis rotamer were thus identified and removed from Table VI on the grounds that transitions between states with character at the trans and cis are expected to be very weak. The torsional states displated in Figure 1 are those from a 2-mode VSCF calculation on the zero-point. The VCI states follow the same pattern and the first isomerising states are 6⌫9 from the trans minimum and 4⌫9 from the cis minimum. The transitions among the cis states are reported in Table VII, where energies are relative to the zero point of the cis rotamer, which is 1412 cm 1 above the zero point of the trans rotamer.

A.

trans-formic acid

A significant number of transitions of gaseous trans-formic acid have been observed in overtone spectroscopy,106 which were previously assigned by fitting the frequencies to a model e↵ective Hamiltonian. The present work allows for an independent assignment of the observed bands as well as prediction of a large number of unobserved bands (the full list is in the SI). An asterisk in the first column marks a band where the assignment from this work di↵ers from that of reference 106 or where a new assignment has been made. Since the ab initio predictions carry an uncertainty and neither rotational structure, nor absorption intensity have been computed, assignment is ambiguous at frequency regions with a high density of states. In these cases the best guess is indicated as a value in square brackets. In the following, the decisions involved in the (re)assignments are explained and justified. Band ⌫5 is in strong Fermi resonance with 2⌫9 . For the band at 1222 cm 1 , the character from the CI coefficients is 66% 2⌫9 and 26% ⌫5 , and the band at 1305 cm 1 is approximately the converse. Thus the ⌫5 fundamental is reassigned according to this work. A similar situation arises for the bands at 1855 and 1933 cm 1 , which although are in slightly poorer agreement with the observed bands at 1848 and 1931 cm 1 , are nevertheless confidently assigned from the CI coefficients. The observed bands at 2132 and 2196.3 cm 1 are in a region of low density of states and assignment is unambiguous in agreement with earlier work. The ab initio values are slightly overestimated, broadly consistent with the overestimation of 4 cm 1 per quanta in ⌫6 and 1 cm per quanta in ⌫8 expected from the accuracy of the fundamentals. The band observed at 2298.6 cm 1 was previously tentatively assigned to ⌫5 +⌫6 , which is ruled out by the

present calculations. The band ⌫6 +2⌫9 , which would replace ⌫5 +⌫6 with the reassignment of ⌫5 , is computed to be at 2312 cm 1 and is also rejected. We assign the band at 2298.6 cm 1 to ⌫7 +⌫8 +⌫9 , which is computed to be 2302 cm 1 and closely matches the expected 2 cm 1 overestimation based on the cancellation of errors in the composite fundamentals. The band observed at 2376 cm 1 is assigned to ⌫6 +⌫7 +⌫9 , computed at 2369 cm 1 , which is considered the more likely match than the neighbouring bands at 2361 (⌫6 +2⌫7 ) or 2402 (⌫3 +⌫7 ). We remark that the uncertainty in prediction is somewhat higher for bands with quanta in 3 distinct modes. The band observed at 2400.2 could be either ⌫3 +⌫7 or ⌫5 +⌫6 and is tentatively assigned to ⌫3 +⌫7 on the basis of the closest energy match. The band at 2504 cm 1 is tentatively assigned to ⌫5 +2⌫9 (2511 cm 1 ) rather than 2⌫7 +2⌫9 (2494 cm 1 ) on similar grounds. The band at 2600 cm 1 could be either ⌫4 +2⌫9 (2600 cm 1 ) or 2⌫5 (2606 cm 1 ), each of which is in strong resonance with another band through the Fermi resonance of ⌫5 and 2⌫9 . A tentative assignment of ⌫4 +2⌫9 is made based on the best energy match. The bands at 2745, 2803, 2876, 2942 and 3570.5 cm 1 are all unambiguously assigned, although we note that the ⌫2 fundamental at 2942 cm 1 is in weak resonance with 2⌫7 +⌫8 +⌫9 (10% character). We deem the accuracy of the ab initio prediction to be sufficient to assign the bands at 3057 and 3106.5 cm 1 to ⌫4 +⌫8 +⌫9 and ⌫4 +⌫6 +⌫7 rather than to ⌫3 +⌫5 (3087 cm 1 ) and 3⌫8 (3090 cm 1 ) or ⌫4 +⌫6 +⌫9 (3115 cm 1 ), respectively. The band at 3152.3 cm 1 is assigned to ⌫3 +⌫4 , which is overestimated by ab initio calculation (3161 cm 1 ) consistent with the overestimation of the ⌫3 fundamental. Similarly, the band assigned to 3⌫6 (3275 cm 1 ) is overestimated by 17 cm 1 by computation and we comment that the significant uncertainty attached to ⌫3 and ⌫6 somewhat complicates assignment of the higher energy bands. The density of states in the region of the band at 3963.6 cm 1 is high and in addition to the previously assigned band ⌫2 +2⌫8 (3952) there are other candidate bands at 3965 (2⌫4 +2⌫9 ) and 3973 cm 1 (⌫4 +2⌫5 ). We are not sufficiently confident in the predictions to change the assignment of this band. For the band at 4043.4 cm 1 , even though the computed value for ⌫2 +⌫6 is in perfect agreement, the uncertainty deduced from the ⌫2 and ⌫6 fundamentals must make this assignment tentative, with other the candidate being 2⌫4 +⌫5 . The band at 4192 cm 1 is assigned with reasonable confidence to ⌫1 +⌫7 (4199 cm 1 ) in line with previous assignment and consistent with the overestimation of the ⌫1 fundamental. The alternative plausible candidate is ⌫2 +2⌫7 at 4196 cm 1 , although this would be expected to be an underestimation based on the accuracy of the ⌫2 and ⌫7 fundamentals. The other candidate is ⌫3 +⌫4 +⌫8 at 4191 cm 1 , but this is also expected to be an underestimation. For the bands above 4200 cm 1 , although the density of states is very high, it is mostly composed of bands

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 10 of 16 10

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

TABLE VI: VCI and observed bands of trans-formic acid in cm 1 . An VCI Exp. Band State asterisk denotes newly assigned bands and square brackets a tentative 3090 3⌫8 95 assignment 3105 [3106.5] ⌫4 +⌫6 +⌫7 96 3115 ⌫4 +⌫6 +⌫9 98 VCI Exp. Band State 3159 2⌫ +⌫ (⌫ +4⌫ ) 102 5 7 7 9 7354 ZP 1 3161 [3152.3] ⌫3 +⌫4 103 627 626.16 ⌫7 2 3166 ⌫6 +2⌫8 104 638 640.72 ⌫9 3 3233 2⌫ +⌫ 109 6 8 1034 1033.47 ⌫8 4 3249 2⌫5 +⌫9 (⌫5 +3⌫9 ) 115 1108 1104.85 ⌫6 5 3288 2⌫ +2⌫ (⌫ +2⌫ ) 119 8 9 5 8 * 1222 1223 2⌫9 (⌫5 ) 6 3292 [3275] 3⌫6 121 1256 1255 2⌫7 7 3367 ⌫5 +2⌫8 (2⌫8 +2⌫9 ) 129 1268 ⌫7 +⌫9 8 3374 2⌫ +⌫ 130 4 7 * 1305 1306.2 ⌫5 (2⌫9 ) 9 3392 2⌫4 +⌫9 (⌫4 +⌫5 +⌫9 ) 131 1379 1380 ⌫4 10 3451 ⌫4 +2⌫8 (⌫3 +⌫8 +⌫9 ) 139 1661 ⌫8 +⌫7 12 3493 ⌫5 +2⌫6 (2⌫6 +2⌫9 ) 146 1672 ⌫8 +⌫9 13 3544 2⌫5 +⌫8 (⌫5 +⌫8 +2⌫9 ) 154 1733 ⌫6 +⌫7 14 3547 2⌫3 155 1739 ⌫6 +⌫9 15 3566 ⌫2 +⌫7 156 1783 1776.83 ⌫3 16 3569 ⌫4 +2⌫6 157 1796 3⌫9 (⌫5 +⌫9 ) 17 3571 ⌫ +⌫ (⌫ +3⌫ ) 159 2 9 3 9 * 1855 1847.8 ⌫7 +2⌫9 (⌫5 +⌫7 ) 18 3575 3570.5 ⌫1 161 1890 3⌫7 19 3579 ⌫ +⌫ (⌫ +3⌫ ) 162 2 9 3 9 1903 2⌫7 +⌫9 20 3588 2⌫5 +⌫6 (⌫5 +4⌫9 ) 164 * 1933 1931.1 ⌫5 +⌫7 (⌫7 +2⌫9 ) 22 3642 2⌫5 +⌫8 (⌫5 +⌫8 +2⌫9 ) 175 1947 ⌫5 +⌫9 23 3784 2⌫4 +⌫8 202 2006 ⌫4 +⌫7 24 3790 3⌫5 (⌫5 +4⌫9 ) 204 2021 ⌫4 +⌫9 25 3833 ⌫3 +2⌫8 211 2063 2⌫8 26 3836 [3826] 2⌫4 +⌫6 212 2138 2132 ⌫6 +⌫8 28 3878 ⌫4 +2⌫5 (⌫4 +⌫5 +2⌫9 ) 221 2205 2196.3 2⌫6 29 3907 2⌫5 +2⌫9 (3⌫5 ) 229 2257 ⌫8 +2⌫9 (⌫5 +⌫8 ) 30 3952 [3963.6] ⌫2 +⌫8 242 2290 2⌫7 +⌫8 31 3973 ⌫4 +2⌫5 (⌫4 +⌫5 +2⌫9 ) 247 * 2302 2298.6 ⌫7 +⌫8 +⌫9 32 3978 ⌫3 +2⌫6 250 2312 ⌫6 +2⌫9 (⌫5 +⌫6 ) 33 4040 2⌫4 +⌫5 (3⌫4 ) 263 * 2338 2338 ⌫5 +⌫8 (⌫8 +2⌫9 ) 34 4043 [4043.4] ⌫ +⌫ 265 2 6 2358 4⌫9 (⌫5 +2⌫9 ) 35 4105 3⌫4 (2⌫4 +2⌫9 ) 280 2361 ⌫6 +2⌫7 36 4160 ⌫2 +2⌫9 (⌫2 +⌫5 ) 293 * 2369 [2376] ⌫6 +⌫7 +⌫9 37 4161 2⌫3 +⌫7 294 * 2402 [2400.2] ⌫3 +⌫7 (⌫5 +⌫6 ) 39 4181 2⌫3 +⌫9 298 2406 ⌫5 +⌫6 (⌫3 +⌫7 ) 40 4191 ⌫3 +⌫4 +⌫8 (⌫6 +3⌫8 ) 302 2415 ⌫4 +⌫8 41 4196 ⌫2 +2⌫7 (⌫4 +2⌫6 +⌫7 ) 305 2420 ⌫3 +⌫9 42 4199 [4192] ⌫1 +⌫7 307 2478 ⌫4 +⌫6 45 4206 ⌫ +⌫ +⌫ 311 2 7 9 2494 2⌫7 +2⌫9 (⌫5 +2⌫7 ) 46 4212 4209 ⌫1 +⌫9 312 * 2511 [2504] ⌫5 +2⌫9 (2⌫5 ) 47 * 4237 [4242] ⌫2 +⌫5 (⌫2 +2⌫9 ) 318 2568 ⌫5 +2⌫7 (2⌫7 +2⌫9 ) 51 4256 ⌫3 +⌫4 +⌫6 322 * 2600 [2600] ⌫4 +2⌫9 (⌫4 +⌫5 ) 54 4259 2⌫6 +2⌫8 325 2608 2⌫5 (⌫5 +2⌫9 ) 55 4291 ⌫ +2⌫ (⌫ +⌫ +2⌫ ) 336 3 5 3 5 9 2636 ⌫4 +2⌫7 56 4298 4300 ⌫2 +⌫4 (⌫3 +4⌫7 ) 338 2678 ⌫4 +⌫5 (⌫4 +2⌫9 ) 59 4390 ⌫3 +2⌫5 (⌫3 +⌫5 +2⌫9 ) 372 2690 ⌫7 +2⌫8 60 * 4460 4450 ⌫3 +⌫4 +⌫5 396 2702 2⌫8 +⌫9 61 4525 4515 ⌫3 +2⌫4 412 2746 2745 2⌫4 62 4567 2⌫3 +⌫8 428 2810 2803 ⌫3 +⌫8 67 4607 4600 ⌫1 +⌫8 447 2825 2⌫6 +⌫9 (⌫6 +3⌫9 ) 68 4645 2⌫3 +⌫6 464 2829 2⌫6 +⌫7 69 4687 4670.4 ⌫ +⌫ (2⌫ +2⌫ ) 489 1 6 5 6 2886 2876 ⌫3 +⌫6 72 4713 4708 ⌫2 +⌫3 498 2938 2942 ⌫2 (2⌫7 +⌫8 +⌫9 ) 77 2940 2⌫7 +⌫8 +⌫9 (⌫2 ) 78 3003 ⌫3 +2⌫9 (⌫3 +⌫5 ) 83 involving several quanta of energy in several modes. On 3027 ⌫3 +2⌫7 86 the grounds that transitions to these bands are expected * 3058 [3057] ⌫4 +⌫8 +⌫9 91 to be very weak compared to transitions to the lower 3087 ⌫3 +⌫5 (⌫3 +2⌫9 ) 94

quanta bands, the observed and predicted bands can be

ACS Paragon Plus Environment

Page 11 of 16

The Journal of Physical Chemistry 11

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

TABLE VII: VCI and observed bands of cis-formic acid in cm 1 . An asterisk denotes newly assigned bands and square brackets a tentative assignment VCI Exp. Band State 1412 0 ZP 11 1904 492 503(6) ⌫9 21 2080 668 661(2) ⌫7 27 2400 988 980(8) 2⌫9 38 2449 1038 ⌫8 44 2514 1103 1108(4) ⌫6 48 2577 1165 ⌫7 +⌫9 52 2667 1255 1244(1) ⌫5 58 2759 1348 2⌫7 63 2806 1394 1396(2) ⌫4 66 2937 1526 ⌫8 +⌫9 76 3009 1597 ⌫6 +⌫9 85 3133 1721 ⌫7 +⌫8 99 3159 1747 ⌫5 +⌫9 101 3179 1767 ⌫6 +⌫7 106 3235 1824 1808(9) ⌫3 113 3299 1887 ⌫4 +⌫9 122 3475 2064 ⌫4 +⌫7 143 3486 2074 2⌫8 144 3536 2124 ⌫8 +2⌫9 (⌫6 +⌫8 ) 152 3570 2159 ⌫6 +⌫8 (⌫8 +2⌫9 ) 158 3595 2184 2⌫6 165 3638 2227 [2207] ⌫5 +2⌫9 174 3715 2304 ⌫5 +⌫8 189 3733 2321 ⌫3 +⌫9 193 3761 2350 2340 ⌫5 +⌫6 199 3854 2442 ⌫4 +⌫8 214 3908 2497 ⌫3 +⌫7 231 3915 2503 2⌫5 234 3918 2506 ⌫4 +⌫6 236 4061 2650 ⌫4 +⌫5 269 * 4170 2759 2760 2⌫4 (⌫2 ) 296 4274 2862 ⌫3 +⌫8 328 4292 2880 2899(11) ⌫2 (2⌫4 ) 337 4322 2911 ⌫3 +⌫5 347 4496 3084 ⌫3 +⌫6 406 4626 3214 ⌫3 +⌫4 456 4672 3261 2⌫4 +⌫9 (⌫2 +⌫9 ) 477

unambiguously matched, except for the band at 4242 cm 1 which, although assigned to ⌫2 +⌫5 (4237 cm 1 ) could also plausibly be ⌫3 +⌫4 +⌫6 . In summary, on the basis of our ab initio calculations, we have been able to assign all 40 of the observed bands below 4750 cm 1 , of which 13 are new assignments. The agreement between the ab initio values and the observed bands is remarkable, particularly since a several of the bands are resonance states.

B.

cis-formic acid

Computed transitions involving states at the cis rotamer are displayed in Table VII. Wavenumbers are listed both for the energies above the global zero point and the transition energy from the zero-point of cisformic acid, together with the character and the state

number from the ICPH-VCI calculation. Only those states with up to two quanta of excitation relative to the cis zero point are listed, plus some additional states relevant to assignment or part of Fermi doublets. This is because the states belonging to the cis rotamer, being higher in energy than those of the trans rotamer, are expected to be less well converged with respect to the size of the CI basis. The full list of states belonging to the cis rotamer can be extracted from the CI vector list in the SI. An examination of the CI coefficients revealed a high level of coupling between modes ⌫4 , ⌫5 , ⌫6 and ⌫7 , making a determination of the character of the bands more challenging than for the trans rotamer. This coupling, however appears to be an artefact of the diabatic definition of the normal coordinates as a function of the torsional internal coordinate. We conclude that although the diabatic scheme was superior in this case to the adiabatic vector definition, an optimal coordinate definition would be somewhere inbetween fully adiabatic and fully diabatic.107 We also remark that the increased coupling adds further uncertainty to the cis levels because the neglected six-body and higher coupling contributions are expected to be larger than those of the trans rotamer. As far as we know, a gas-phase IR absorption spectra has not yet been recorded for cis-formic acid due to its very low population at room temperature24 and thermal instability.108 IR spectra have been recorded for cis-formic acid trapped in an argon matrix.15,24 The observed bands are included in Table VII, taken primarily from Ref 24. The additional band at 2340 cm 1 observed in Ref 15 is also included. Where possible, an estimated experimental uncertainty is computed by taking the di↵erence between the gas-phase and argon-matrix wavenumbers for the equivalent bands of trans-formic acid. This uncertainty is given in parentheses after the experimental band energy. Since the calculations are gas phase and the experiments performed in an argon matrix, perfect agreement is not expected. This, and the anticipated lower ab initio accuracy notwithstanding, the experimental bands are readily assigned on the basis of our calculations. There are a two features of note. First, the band at 2760 cm 1 , previously unassigned, is unambiguously assigned to 2⌫4 (2759 cm 1 ), with no other plausible candidates. The band at 2⌫4 is actually part of a Fermi doublet with ⌫2 . The second feature concerns the band at 2207, which was assigned to ⌫5 +2⌫9 in Ref 24 and 2⌫6 in Ref 15. Our calculations are unfortunately inconclusive, with the observed band lying midway between the computed values for ⌫5 +2⌫9 and 2⌫6 . The ⌫5 +2⌫9 state is not likely to be tightly converged with respect to the modal basis, so is the more likely candidate from variational considerations. More information, for example from rotational structure or absorption intensities would also be helpful to resolve this and the other ambiguous assignments.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 12 of 16 12

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

VI.

CONCLUSION

This work is a successful demonstration of predictive power of modern ab initio methodology for the treatment of the ro-vibrational eigenstates of polyatomic molecules with complex dynamics including an internal rotation and strong Fermi resonances. Several methodological conclusions can be summarised from this work. The CCSD(T) method, at the basis set limit and correlating valence electrons, provides a Born–Oppenheimer electronic potential energy surface with eigenstates within 3 cm 1 of experiment for the fundamental transitions. This high level of predictive accuracy can reasonably be expected for all small non-aromatic closed-shell organic compounds. The, now standard, F12 approach makes it routinely possible to compute tens of thousands of near basis set limit CCSD(T) energies at low cost. Our calculations also demonstrate that fitting an analytic model potential energy function to the ab initio data is straightforward with our newly developed LASSO-based flexible correction surface approach. Global fits from this method are a priori free from unphysical holes and can be refined to the point where the fitting error is smaller than the intrinsic error in the CCSD(T)(F12*) ab initio data. Methodological conclusions can also be drawn from the calculation of the J = 0 vibrational eigenstates. The ICPH approach provides a natural coordinate system for molecules of this kind, with one large amplitude motion. The mean field treatment in these coordinates is qualitatively accurate, which greatly facilitates band assignment from examination of the CI coefficients. For the ICP approach with a torsional internal coordinate, the Eckart conditions for the path must exhibit the 2⇡ periodicity, which is violated if the path is defined with the zero angular momentum condition applied in previous works. The practical implementation in our program DYNAMOL, which uses a cubic spline interpolation of the coordinate system, a multi-mode expansion and numerical quadratures, is numerically robust. The hardest part of the calculation is the definition of the orthonormal coordinates perpendicular to the curvilinear path.

⇤ 1

2

3

4

[email protected] Khare, P.; Kumar, N.; Kumari, K. M.; Srivastava, S. S. Atmospheric formic and acetic acids: An overview. Rev. Geophys. 1999, 37, 227–248. Zuckerman, B.; Ball, J. A.; Gottlieb, C. A. Microwave Detection of Interstellar Formic Acid. ApJ 1971, 163, L41. Lattanzi, V.; Walters, A.; Drouin, B. J.; Pearson, J. C. Submillimeter Spectrum of Formic Acid. Astrophys. J. Suppl. S. 2008, 176, 536. Snyder, L. E. Interferometric Observations of Large Biologically Interesting Interstellar and Cometary Molecules. 2006, 103, 12243.

The standard adiabatic approach led to large many-body couplings in the kinetic energy operator due to Coriolis terms between the torsional motion and the associated change in the perpendicular coordinates. The fully diabatic scheme we adopted instead worked well, but introduced sizeable many-body coupling terms in the potential energy operator at the cis rotamer. An optimal choice of perpendicular coordinates must lie somewhere in between and should smoothly switch between fully diabatic and fully adiabatic depending on the nature of the dynamics of the system. With regards to the overtone spectra of cis- and transformic acid, a full J = 0 line list has been computed for the first 500 states, which lie from 0 to 4722 cm 1 above the zero point at 7354 cm 1 . All 40 observed transitions of the trans rotamer in that region, and all 11 of the cis rotamer have been assigned, with 14 of these being new assignments, or reassignments. A large fraction of the reassignments arise from the CI calculations predicting that the Fermi doublet bands at 1223 and 1306.2 cm 1 for the trans rotamer are 2⌫9 and ⌫5 , respectively, which is the converse assignment from earlier works, and propagates through to the reassignment of overtone and combination bands. The RMS deviation between ab initio prediction and the recorded gas phase spectrum is only 7 cm 1 . ACKNOWLEDGMENTS

D.P.T thanks the Royal Society for a University Research Fellowship. W.M. thanks the ERC for a MarieCurie fellowship (IIF-301616). This article is dedicated to the memory of Prof. N. C. Handy. SUPPORTING INFORMATION

A fortran subroutine for computing the potential energy and optionally the gradient and hessian. A file containing the full set of CI states with primary coefficients from Tables VI and VII.

5

6

7

Willemot, E.; Dangoisse, D.; Bellet, J. Microwave spectrum of formic acid and its isotopic species in D, 13C and 18O. Study of Coriolis resonances between 7 and 9 vibrational excited states. J. Mol. Spectrosc. 1978, 73, 96 – 119. Lerner, R. G.; Dailey, B. P.; Friend, J. P. Microwave Spectrum and Structure of Formic Acid. J. Chem. Phys. 1957, 26, 680–683. Davis, R.; Robiette, A.; Gerry, M.; Bjarnov, E.; Winnewisser, G. Microwave spectra and centrifugal distortion constants of formic acid containing 13C and 18O: Refinement of the harmonic force field and the molecular structure. J. Mol. Spectrosc. 1980, 81, 93 – 109.

ACS Paragon Plus Environment

Page 13 of 16

The Journal of Physical Chemistry 13

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

8

9

10

11

12

13

14

15

16

17

18

19

20

21 22

23

24

Winnewisser, M.; Winnewisser, B. P.; Stein, M.; Birk, M.; Wagner, G.; Winnewisser, G.; Yamada, K. M.; Belov, S. P.; Baskakov, O. I. Rotational Spectra of cisHCOOH, trans-HCOOH, and trans-H13COOH. J. Mol. Spectrosc. 2002, 216, 259 – 265. Cazzoli, G.,; Puzzarini, C.,; Stopkowicz, S.,; Gauss, J., Hyperfine structure in the rotational spectra of transformic acid: Lamb-dip measurements and quantumchemical calculations*. A&A 2010, 520, A64. Bertie, J. E.; Michaelian, K. H. The Raman spectra of gaseous formic acid h2 and d2. J. Chem. Phys. 1982, 76, 886–894. Olbert-Majkut, A.; Ahokas, J.; Lundell, J.; Pettersson, M. Raman spectroscopy of formic acid and its dimers isolated in low temperature argon matrices. Chem. Phys. Lett. 2009, 468, 176–183. Millikan, R. C.; Pitzer, K. S. Infrared Spectra and Vibrational Assignment of Monomeric Formic Acid. J. Chem. Phys. 1957, 27, 1305–1308. Luiz, G.; Scalabrin, A.; Pereira, D. Gas phase infrared Fourier transform spectra of H12 COOH and H13COOH. Infrared Phys. Technol. 1997, 38, 45 – 49. Tan, T.; Goh, K.; Ong, P.; Teo, H. Rovibrational Constants for the 6 and 29 Bands of HCOOD by Fourier Transform Infrared Spectroscopy. J. Mol. Spectrosc. 1999, 198, 110 – 114. Ma¸coas, ˙ E. M.; Lundell, J.; Pettersson, M.; Khriachtchev, L.; Fausto, R.; R¨ as¨ anen, M. Vibrational spectroscopy of cis- and trans-formic acid in solid argon. J. Mol. Spectrosc. 2003, 219, 70 – 80. Miyazawa, T.; Pitzer, K. S. Internal Rotation and Infrared Spectra of Formic Acid Monomer and Normal Coordinate Treatment of OutofPlane Vibrations of Monomer, Dimer, and Polymer. J. Chem. Phys. 1959, 30, 1076–1086. Redington, R. L. Vibrational spectra and normal coordinate analysis of isotopically labeled formic acid monomers. J. Mol. Spectrosc. 1977, 65, 171 – 189. Madeja, F.; Markwick, P.; Havenith, M.; Nauta, K.; Miller, R. E. Rotationally resolved infrared spectroscopy of h2- and d1-formic acid monomer in liquid He droplets. J. Chem. Phys. 2002, 116, 2870–2878. Baskakov, O. I.; Markov, I. A.; Alekseev, E. A.; Motiyenko, R. A.; Lohilahti, J.; Horneman, V.-M.; Winnewisser, B. P.; Medvedev, I. R.; De Lucia, F. C. Simultaneous analysis of rovibrational and rotational data for the 4 1, 5 1, 6 1, 7 2, 8 1, 7 1 9 1 and 9 2 states of HCOOH. J. of Mol. Struct. 2006, 795, 54–77. Baskakov, O.; Horneman, V.-M.; Alanko, S.; Lohilahti, J. FTIR spectra of the ⌫ 6 and ⌫ 8 bands of 13 C formic acid moleculeAssignment of FIR-laser lines. J. Mol. Spectrosc. 2008, 249, 60–64. Hocking, W. H. Z. Naturforsch. A 1976, 31, 1113. Baskakov, O. I.; Winnewisser, B. P.; Medvedev, I. R.; Lucia, F. C. D. The millimeter wave spectrum of cis-HCOOH in the ground state and in the v9=1 and v7=1 excited vibrational states, and cis-H13COOH in the ground state. J. of Mol. Struct. 2006, 795, 42 – 48. Baskakov, O. I.; Horneman, V.-M.; Lohilahti, J.; Alanko, S. High resolution FTIR spectra of the ⌫ 9 vibrational band of cis-rotamers HCOOH and H 13 COOH. J. of Mol. Struct. 2006, 795, 49–53. Pettersson, M.; Lundell, J.; Khriachtchev, L.; R¨ as¨ anen, M. IR Spectrum of the Other Rotamer of Formic Acid, cis-HCOOH. J. Am. Chem. Soc. 1997, 119,

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

11715–11716. Ma¸co ˆas, E.; Khriachtchev, L.; Pettersson, M.; Juselius, J.; Fausto, R.; R¨ as¨ anen, M. Reactive vibrational excitation spectroscopy of formic acid in solid argon: Quantum yield for infrared induced trans cis isomerization and solid state e↵ects on the vibrational spectrum. J. Chem. Phys. 2003, 119, 11765–11772. Pettersson, M.; Ma¸co ˆas, E.; Khriachtchev, L.; Lundell, J.; Fausto, R.; R¨ as¨ anen, M. Cis trans conversion of formic acid by dissipative tunneling in solid rare gases: Influence of environment on the tunneling rate. J. Chem. Phys. 2002, 117, 9095–9098. Pettersson, M.; Ma¸co ˆas, E. M.; Khriachtchev, L.; Fausto, R.; R¨ as¨ anen, M. Conformational isomerization of formic acid by vibrational excitation at energies below the torsional barrier. J. Am. Chem. Soc. 2003, 125, 4058–4059. Paulson, L. O.; Anderson, D. T.; Lundell, J.; Marushkevich, K.; Melavuori, M.; Khriachtchev, L. Conformation Resolved Induced Infrared Activity: trans-and cisFormic Acid Isolated in Solid Molecular Hydrogen. J. Phys. Chem. A 2011, 115, 13346–13355. Domanskaya, A.; Marushkevich, K.; Khriachtchev, L.; R¨ as¨ anen, M. Spectroscopic study of cis-to-trans tunneling reaction of HCOOD in rare gas matrices. J. Chem. Phys. 2009, 130, 154509. Bock, C. W.; Trachtman, M.; George, P. An ab initio study of the geometries, anharmonic force fields and fundamental vibration frequencies of cis-and trans-formic acid. J. Mol. Spectrosc. 1980, 80, 131–144. Demaison, J.; Herman, M.; Livin, J. Anharmonic force field of cis- and trans-formic acid from high-level ab initio calculations, and analysis of resonance polyads. J. Chem. Phys. 2007, 126 . Scribano, Y.; Benoit, D. M. Calculation of vibrational frequencies through a variational reduced-coupling approach. J. Chem. Phys. 2007, 127, 164118. Maeda, S.; Watanabe, Y.; Ohno, K. Finding important anharmonic terms in the sixth-order potential energy function by the scaled hypersphere search method: An application to vibrational analyses of molecules and clusters. J. Chem. Phys. 2008, 128, 144111. Paulson, L. O.; Kaminsky, J.; Anderson, D. T.; Bour, P.; Kubelka, J. Theoretical Study of Vibrationally Averaged Dipole Moments for the Ground and Excited C O Stretching States of trans-Formic Acid. J. Chem. Theory Comput. 2010, 6, 817–827. Mizukami, W.; Tew, D. P. A second-order multi-reference perturbation method for molecular vibrations. J. Chem. Phys. 2013, 139 . Born, M.; Oppenheimer, R. Zur Quantentheorie der Molekeln. Ann. Phys. 1927, 389, 457–484. Raghavachari, K.; Trucks, G. W.; Pople, J. A.; HeadGordon, M. A fifth-order perturbation comparison of electron correlation theories. Chem. Phys. Lett. 1989, 157, 479–483. Ruden, T. A.; Helgaker, T.; Jørgensen, P.; Olsen, J. Coupled-cluster connected quadruples and quintuples corrections to the harmonic vibrational frequencies and equilibrium bond distances of HF, N2,F2, and CO. J. Chem. Phys. 2004, 121, 5874–5884. Heckert, M.; K´ allay, M.; Gauss, J. Molecular equilibrium geometries based on coupled-cluster calculations including quadruple excitations. Mol. Phys. 2005, 103, 2109–2115.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 14 of 16 14

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

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

H¨ attig, C.; Klopper, W.; K¨ ohn, A.; Tew, D. P. Explicitly Correlated Electrons in Molecules. Chem. Rev. 2012, 112, 4–74. Kong, L.; Bischo↵, F. A.; Valeev, E. F. Explicitly Correlated R12/F12 Methods for Electronic Structure. Chem. Rev. 2012, 112, 75–107. Tew, D. P.; Klopper, W.; Neiss, C.; H¨ attig, C. Quintuple⇣ quality coupled-cluster correlation energies with triple-⇣ basis sets. Phys. Chem. Chem. Phys. 2007, 9, 1921–1930. Tew, D. P.; Klopper, W.; Neiss, C.; H¨ attig, C. Comment on Quintuple-⇣ quality coupled-cluster correlation energies with triple-⇣ basis sets by D.P. Tew, Wim Klopper, C. Neiss and C. H¨ attig, Phys. Chem. Chem. Phys., 2007, 9, 1921 [erratum]. Phys. Chem. Chem. Phys. 2008, 10, 6325–6327. H¨ attig, C.; Tew, D. P.; K¨ ohn, A. Communications: Accurate and efficient approximations to explicitly correlated coupled-cluster singles and doubles, CCSD-F12. J. Chem. Phys. 2010, 132, 231102. Ne↵, M.; Rauhut, G. Toward large scale vibrational configuration interaction calculations. J. Chem. Phys. 2009, 131, 124129. Kidwell, N. M.; Li, H.; Wang, X.; Bowman, J. M.; Lester, M. I. Unimolecular dissociation dynamics of vibrationally activated CH3CHOO Criegee intermediates to OH radical products. Nat. Chem. 2016, 8, 509 – 514. Bowman, J. M.; Carrington Jr, T.; Meyer, H.-D. Variational quantum approaches for computing vibrational energies of polyatomic molecules. Mol. Phys. 2008, 106, 2145–2182. Cs´ asz´ ar, A. G.; F´ abri, C.; Szidarovszky, T.; M´ atyus, E.; Furtenbacher, T.; Czak´ o, G. The fourth age of quantum chemistry: molecules in motion. Phys. Chem. Chem. Phys. 2011, 14, 1085–1106. Christiansen, O. Selected new developments in vibrational structure theory: potential construction and vibrational wave function calculations. Phys. Chem. Chem. Phys. 2012, 14, 6672–6687. Tennyson, J. Perspective: Accurate ro-vibrational calculations on small molecules. J. Chem. Phys. 2016, 145, 120901. Owens, A.; Yurchenko, S. N.; Yachmenev, A.; Tennyson, J.; Thiel, W. Accurate ab initio vibrational energies of methyl chloride. J. Chem. Phys. 2015, 142, 244306. Bowman, J. M.; Wang, X.; Homayoon, Z. Ab initio computational spectroscopy and vibrational dynamics of polyatomic molecules: Applications to syn and anti-CH 3 CHOO and NO 3. J. Mol. Spectrosc. 2015, 311, 2–11. Delahaye, T.; Nikitin, A. V.; Rey, M.; Szalay, P. G.; Tyuterev, V. G. Accurate 12D dipole moment surfaces of ethylene. Chem. Phys. Lett. 2015, 639, 275–282. Sarka, J.; Cs´ asz´ ar, A. G. Interpretation of the vibrational energy level structure of the astructural molecular ion H5+ and all of its deuterated isotopomers. J. Chem. Phys. 2016, 144, 154309. Wang, X.-G.; Carrington Jr, T. Calculated rotationbending energy levels of CH5+ and a comparison with experiment. J. Chem. Phys. 2016, 144, 204304. Yu, H.-G. An exact variational method to calculate rovibrational spectra of polyatomic molecules with large amplitude motion. J. Chem. Phys. 2016, 145, 084109. Cassam Chena¨ı, P.; Li´evin, J. The VMFCI method: A flexible tool for solving the molecular vibration problem.

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

J. Comput. Chem. 2006, 27, 627–640. Manthe, U. Layered discrete variable representations and their application within the multiconfigurational timedependent Hartree approach. J. Chem. Phys. 2009, 130, 054109. Jung, J. O.; Gerber, R. B. Vibrational wave functions and spectroscopy of (H2O)n, n=2,3,4,5: Vibrational selfconsistent field with correlation corrections. J. Chem. Phys. 1996, 105, 10332–10348. Carter, S.; Culik, S. J.; Bowman, J. M. Vibrational selfconsistent field method for many-mode systems: A new approach and application to the vibrations of CO adsorbed on Cu(100). J. Chem. Phys. 1997, 107, 10458– 10469. Benoit, D. M. Efficient correlation-corrected vibrational self-consistent field computation of OH-stretch frequencies using a low-scaling algorithm. J. Chem. Phys. 2006, 125, 244110. Godtliebsen, I. H.; Hansen, M. B.; Christiansen, O. Tensor decomposition techniques in the solution of vibrational coupled cluster response theory eigenvalue equations. J. Chem. Phys. 2015, 142 . Scribano, Y.; Lauvergnat, D. M.; Benoit, D. M. Fast vibrational configuration interaction using generalized curvilinear coordinates and self-consistent basis. J. Chem. Phys. 2010, 133, 094103. Hougen, J.; Bunker, P.; Johns, J. The vibration-rotation problem in triatomic molecules allowing for a largeamplitude bending vibration. J. Mol. Spectrosc. 1970, 34, 136 – 172. Miller, W. H.; Handy, N. C.; Adams, J. E. Reaction path Hamiltonian for polyatomic molecules. J. Chem. Phys. 1980, 72, 99–112. Carrington Jr., T.; Miller, W. H. Reaction surface Hamiltonian for the dynamics of reactions in polyatomic systems. J. Chem. Phys. 1984, 81, 3942–3950. Tew, D. P.; Handy, N. C.; Carter, S. A reaction surface Hamiltonian study of malonaldehyde. J. Chem. Phys. 2006, 125, 084313. Carter, S.; Handy, N. C. The vibrations of H2O2, studied by multimode, with a large amplitude motion. J. Chem. Phys. 2000, 113, 987–993. Carter, S.; Handy, N. C.; Bowman, J. M. High torsional vibrational energies of H2O2 and CH3OH studied by MULTIMODE with a large amplitude motion coupled to two e↵ective contraction schemes. Mol. Phys. 2009, 107, 727–737. Tew, D. P.; Handy, N. C.; Carter, S. Glyoxal studied with ’Multimode ’, explicit large amplitude motion and anharmonicity. Phys. Chem. Chem. Phys. 2001, 3, 1958–1964. Tew, D. P.; Handy, N. C.; Carter, S. The vibrations of glyoxal, studied by ¡80¿¡98¿Multimode¡80¿¡99¿, with a large amplitude motion, using an ab initio potential surface. Mol. Phys. 2001, 99, 393–402. L´eonard, C.; Handy, N. C.; Carter, S.; Bowman, J. The vibrational levels of ammonia. Spectrochim. acta. A 2002, 58, 825–838. Tew, D. P.; Handy, N. C.; Carter, S.; Irle, S.; Bowman, J. The internal coordinate path Hamiltonian; application to methanol and malonaldehyde. Mol. Phys. 2003, 101, 3513–3525. Mizukami, W.; Habershon, S.; Tew, D. P. A compact and accurate semi-global potential energy surface for malonaldehyde from constrained least squares regression. J.

ACS Paragon Plus Environment

Page 15 of 16

The Journal of Physical Chemistry 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

75

76

77

78

79

80

81

82

83

84

85

86 87 88

89

90

91

92

Chem. Phys. 2014, 141, –. Tibshirani, R. Regression Shrinkage and Selection via the Lasso. J. Royal Statist. Soc. B 1996, 58, 267–288. Avila, G.; Carrington, T. Using nonproduct quadrature grids to solve the vibrational Schr¨ odinger equation in 12D. J. Chem. Phys. 2011, 134, 054126–054126. Chapman, S.; Dupuis, M.; Green, S. Theoretical threedimensional potential-energy surface for the reaction of Be with HF. Chem. Phys. 1983, 78, 93–105. Bowman, J. M.; Bittman, J. S.; Harding, L. B. Abinitio calculations of electronic and vibrational energies of HCO and HOC. J. Chem. Phys. 1986, 85, 911–921. Ischtwan, J.; Collins, M. A. Molecular potential energy surfaces by interpolation. J. Chem. Phys. 1994, 100, 8080–8088. Jordan, M. J. T.; Thompson, K. C.; Collins, M. A. Convergence of molecular potential energy surfaces by interpolation: Application to the OH+H2H2O+H reaction. J. Chem. Phys. 1995, 102, 5647–5657. Collins, M. A. Molecular potential-energy surfaces for chemical reaction dynamics. Theor. Chem. Acc. 2002, 108, 313–324. Maisuradze, G. G.; Thompson, D. L.; Wagner, A. F.; Minko↵, M. Interpolating moving least-squares methods for fitting potential energy surfaces: Detailed analysis of one-dimensional applications. J. Chem. Phys. 2003, 119, 10002–10014. Guo, Y.; Kawano, A.; Thompson, D. L.; Wagner, A. F.; Minko↵, M. Interpolating moving least-squares methods for fitting potential energy surfaces: Applications to classical dynamics calculations. J. Chem. Phys. 2004, 121, 5091–5097. Braams, B. J.; Bowman, J. M. Permutationally invariant potential energy surfaces in high dimensionality. Int. Rev. Phys. Chem. 2009, 28, 577–606. Bart´ ok, A. P.; Payne, M. C.; Kondor, R.; Cs´ anyi, G. Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons. Phys. Rev. Lett. 2010, 104, 136403. Bart´ ok, A. P.; Kondor, R.; Cs´ anyi, G. On representing chemical environments. 2013, 87, 184115. Wang, X.; Carter, S.; Bowman, J. M. Fortenberry, R. C.; Yu, Q.; Mancini, J. S.; Bowman, J. M.; Lee, T. J.; Crawford, T. D.; Klemperer, W. F.; Francisco, J. S. Communication: Spectroscopic consequences of proton delocalization in OCHCO+. J. Chem. Phys. 2015, 143 . Qu, C.; Bowman, J. M. An ab initio potential energy surface for the formic acid dimer: zero-point energy, selected anharmonic fundamental energies, and ground-state tunneling splitting calculated in relaxed 1-4-mode subspaces. Phys. Chem. Chem. Phys. 2016, 18, 24835–24840. Ho, T.-S.; Rabitz, H. Reproducing kernel Hilbert space interpolation methods as a paradigm of high dimensional model representations: Application to multidimensional potential energy surface construction. J. Chem. Phys. 2003, 119, 6433–6442. Manzhos, S.; Wang, X.; Dawes, R.; Carrington, T. A Nested Molecule-Independent Neural Network Approach for High-Quality Potential Fits. J. Phys. Chem. A 2006, 110, 5295–5304. Handley, C. M.; Popelier, P. L. A. Potential Energy Surfaces Fitted by Artificial Neural Networks. J. Phys. Chem. A 2010, 114, 3371–3383.

93

94

95

96

97

98

99

100

101

102

103 104

105

106

107

108

Behler, J. Representing potential energy surfaces by highdimensional neural network potentials. J. Phys.: Condens. Matter 2014, 26, 183001. Li, J.; Carter, S.; Bowman, J. M.; Dawes, R.; Xie, D.; Guo, H. Behler, J. Neural network potential-energy surfaces in chemistry: a tool for large-scale simulations. Phys. Chem. Chem. Phys. 2011, 13, 17930–17955. Zhang, D. H.; Guo, H. Recent Advances in Quantum Dynamics of Bimolecular Reactions. Ann. Rev. Phys. Chem. 2016, 67, 135–158, PMID: 26980305. Peterson, K. A.; Adler, T. B.; Werner, H.-J. Systematically convergent basis sets for explicitly correlated wavefunctions: The atoms H, He, BNe, and AlAr. J. Chem. Phys. 2008, 128 . Tew, D. P.; Klopper, W. New correlation factors for explicitly correlated electronic wave functions. J. Chem. Phys. 2005, 123 . Weigend, F.; K¨ ohn, A.; Httig, C. Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations. J. Chem. Phys. 2002, 116 . Yousaf, K. E.; Peterson, K. A. Optimized auxiliary basis sets for explicitly correlated methods. J. Chem. Phys. 2008, 129 . K¨ ohn, A.; Tew, D. P. Explicitly correlated coupled-cluster theory using cusp conditions. I. Perturbation analysis of coupled-cluster singles and doubles (CCSD-F12). J. Chem. Phys. 2010, 133 . Luckhaus, D. 6D vibrational quantum dynamics: Generalized coordinate discrete variable representation and (a)diabatic contraction. J. Chem. Phys. 2000, 113, 1329– 1347. Strauss, H. L.; Pickett, H. M. Cubic spline interpolation retains the Eckart condition exactly at all interpollation points when the path is oriented using Eckart case ii). If instead case i), or iii) are applied, points along the path from cubic spline interpolation do not exactly obey these conditions. Carter, S.; Bowman, J. M.; Handy, N. C. Extensions and tests of “multimode”: a code to obtain accurate vibration/rotation energies of many-mode molecules. Theor. Chem. Acc. 1998, 100, 191–198. Freytes, M.; Hurtmans, D.; Kassi, S.; Li´evin, J.; Vander Auwera, J.; Campargue, A.; Herman, M. Overtone spectroscopy of formic acid. Chem. Phys. 2002, 283, 47– 61. Kraka, E. Reaction path Hamiltonian and the unified reaction valley approach. WIREs Comput. Mol. Sci. 2011, 1, 531–556. Blake, P. G.; Davies, H. H.; Jackson, G. E. Dehydration mechanisms in the thermal decomposition of gaseous formic acid. J. Chem. Soc. B 1971, 1923–1925.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 16 of 16 16

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 ACS Paragon Plus Environment