In a recent work we constructed a quasi-diabatic representation, Hd, of the 1, 2, 31A adiabatic states of phenol from high level multireference single and double excitation configuration interaction electronic structure data, energies, energy gradients, and derivative couplings. That Hd accurately describes surface minima, saddle points, and also regions of strong nonadiabatic interactions, reproducing the locus of conical intersection seams and the coordinate dependence of the derivative couplings. The present work determines the accuracy of Hd for describing phenol photodissociation. Additionally, we demonstrate that a modest energetic shift of two diabats yields a quantifiably more accurate Hd compared with experimental energetics. The analysis shows that in favorable circumstances it is possible to use single point energies obtained from the most reliable electronic structure methods available, including methods for which the energy gradients and derivative couplings are not available, to improve the quality of a global representation of several coupled potential energy surfaces. Our data suggest an alternative interpretation of kinetic energy release measurements near λphot ∼ 248 nm.

The photodissociation of phenol to produce H and the phenoxyl radical

C 6 H 5 OH ( X ̃ 1 A ) + h ν C 6 H 5 OH ( 2 , 3 1 A ) H + C 6 H 5 O ( X ̃ 2 B 1 , A ̃ 2 B 2 )
(1)

is understood to involve a bright S1(1π*π)←S0(1ππ) excitation, followed by the radiationless transition 1π*π1σ*π which provides a pathway for dissociation.1–4 The appreciation that this process is not restricted to phenol and is relevant to other heteronuclear cyclic molecules including azoles, nucleobases, and aromatic amino acids has motivated a large number of experimental and computational studies of this system over the past decade.1–20 

We recently reported a four, coupled, quasi-diabatic state representation, Hd, of the 1, 2, 31A states of phenol and defined the representation on a set of nuclear configurations, the domain of definition, sufficient to describe reaction (1).17 We will refer to that work as ZY. The Hd in ZY treats all 33 internal coordinates evenhandedly, while providing a quantifiably accurate representation of the nonadiabatic parameters including the locus of conical intersection seams and the geometry dependence of the derivative couplings. That work documented how accurately the Hd, reproduces the high quality ab initio multireference single and double excitation configuration interaction data (energies, energy gradients, and derivative couplings) from which it is built. Here we establish the quality of the representation, that is, how realistically the Hd describes parameters (equilibrium and saddle point structures, harmonic frequencies, excitation and dissociation energies, conical intersection seams, etc.) describing phenol and its photodissociation. We show that a simple shifting of the diabatic potential energy functions can significantly reduce deviations from experimental results and we investigate the impact of those shifts on the global representation. With an eye toward future dynamical studies, we demonstrate that full dimensional quasi-classical surface hopping trajectories can propagate on the coupled surfaces for (at least) tens of picoseconds without encountering surface anomalies. We also investigate the sensitivity of the quasi-classical nonadiabatic dynamics to the shifting process. Finally we suggest an alternative interpretation of earlier kinetic energy release spectra.1 This analysis sets the stage for detailed dynamical studies of reaction (1) based on the current Hd.

Section II describes the original representation. Particular attention is paid to the symmetry of the diabatic states and the geometry dependence of the diabatic energies. Based on this analysis Section III discusses the shifted form of the representation. Section IV presents a detailed analysis of the quality of the original and shifted representations. Section V summarizes and concludes.

The quasi diabatic representation Hd has the form

H d , s ( R ) = l = 1 N c V l [ P ˆ Ir [ u ( l ) ] , Ir [ v ( l ) ] g ( l ) ( R ) ] B u ( l ) , v ( l ) + i = 1 N state V i s B i , i H d , 0 + i = 1 N state V i s B i , i .
(2)

In Eq. (2)R denotes the nuclear coordinates, P ˆ α , β is a standard group theoretical projector for the complete nuclear permutation inversion group,21,22g(l) are the monomial functions, themselves products of basic functions, used to represent the matrix elements of Hd, and Bu,v is an Nstate × Nstate symmetric matrix with a 1 in the (u, v) and (v, u) elements and the remaining elements are 0. Here Nstate = 4. The procedure for fitting Hd,0 to ab initio data, which determines the V and g(l) in Eq. (2), is described in Ref. 10. A complete tabulation of the g(l) used to construct Hd is given in Ref. 17. The Vs are empirical shifts whose determination in the present case is described in Sec. III. Hd,0 is denoted the nascent representation and Hd,s is denoted the shifted representation. Generic representations are denoted Hd. Accompanying Hd,w (w = s or 0) is an electronic Schrödinger equation

[ H d , w ( R ) I E a , J , ( m ) ( R ) ] d J ( R ) = 0
(3)

from which the energies Ea.J,(m)(R), energy gradients ∇Ea,J(m)(R), and derivative couplings are obtained. The superscript m, for model, indicates that the energy is Hd determined.

Table I reports the basic parameters of the fit, the number of linear parameters, Nc; the number of fit data points, Npt; the number of equations, Neq; the root mean square (RMS) energy error; the mean unsigned error (MUE) for the energy; and the percent errors for the energy gradients. More details concerning Hd, its construction, and its accuracy can be found in ZY.

TABLE I.

Parameters of the fits, V 1 B 1 s = 1039 . 2 , V 1 B 2 s = 1800 . 1 cm−1.

Hd,0 Hd,s
Nc  28 716  28 716 
Npt  7 379  7 709 
Neq  883 051  907 447 
Energy error (cm−1) for electronic energies ≤ 55 000 cm−1 
RMSE  315.03  324.87 
MUE  207.8  213.26 
Energy gradient error (%) for electronic energies ≤ 55 000 cm−1 
RMSE  8.43  8.62 
MUE  6.22  6.36 
Hd,0 Hd,s
Nc  28 716  28 716 
Npt  7 379  7 709 
Neq  883 051  907 447 
Energy error (cm−1) for electronic energies ≤ 55 000 cm−1 
RMSE  315.03  324.87 
MUE  207.8  213.26 
Energy gradient error (%) for electronic energies ≤ 55 000 cm−1 
RMSE  8.43  8.62 
MUE  6.22  6.36 

A key attribute of Hd,w is its domain of definition. The domain of definition is a set of nuclear configurations R(n), 1 ≤ nNpt, sufficient to describe the nonadiabatic process in question. For Hd,0 it was determined by the requirement that more than K% of Nsh quasi classical surface hopping trajectories, in the full 33 dimensions, originating on the 31A potential energy surface (PES) stay near existing R(n). In ZY, K = 99.5 and Nsh = 1000.

A key element of this work is the symmetry of the diabatic states and the coordinate dependence of the diabatic state energies. The large amplitude motion exhibited by phenol dissociation, precludes the use point group symmetry to describe nuclear motion. Instead it is appropriate to use the permutation inversion group symmetry.21 Since only the O–H bond is broken in reaction (1) only one concerted permutation of the atoms is feasible: the simultaneous permutation of the entire benzene ring (2,6)(3,5)(8,12)(9,11), where (i, j) denotes a transposition of atoms i and j. This molecular symmetry group,22 the feasible subgroup of the complete permutation inversion group, is the direct product of inversion and the aforementioned permutation operation. It is isomorphic to the G4 permutation inversion group and the C2v point group. The characters of the diabatic states can be determined at geometries that have point group symmetry where the point group operations can be mapped to molecular symmetry group operations. For Cs geometries point group symmetry will be a proper subgroup of the molecular symmetry group.

Depending on the arrangement and orientation of the C1O7H13 moiety relative to the benzene ring (taken as planar) phenol will exhibit three point group symmetries: (i) C2v: for C1O7H13 collinear and along the C4C1 axis (the highest energy structures); (ii) Cs-planar: for C1O7H13 planar and in the benzene plane (the lowest energy structures); and (iii) Cs-bijct: for C1O7H13 planar and in a plane perpendicular to the benzene ring and containing the C4C1 axis. See Fig. 1(a) for the atom numbering used in this work. Table II provides the C2v–Cs correlations. Calculations in the Franck-Condon region for Cs-planar structures demonstrate that the first four states carry 1A′, 1A′, 1A″, 1A″ symmetry, respectively. At structures with Cs-bijct symmetry the first four states carry 1A′, 1A″, 1A′, 1A″ symmetry. Therefore using Table II the lowest state carries 1A1 symmetry and the remaining states carry 1B2, 1B1, and 1A2 symmetry in G4. A similar result is obtained using an energy optimized C2v structure in the Franck-Condon region but since these states are much higher in energy in this region, problems with root switching could be encountered. Using C2v point group symmetry in the dissociation channel shows that the four low-lying states of phenoxyl carry 1B1, 1B2, 1A2, and 1A1, symmetry, respectively. Thus it is necessary to include four states to account for the change in the 1A1 state from the ground state in the Franck-Condon region to the 41A excited state in the phenoxyl channel. To illustrate these ideas Fig. 2(a) reports the diabatic states and their energy dependence for a coplanar path explained in that figure. The corresponding adiabatic state energies are reported in Fig. 2(b).

FIG. 1.

Key geometrical parameters for phenol. Experimental values above Hd,0 determined values. Hd,s determined values differ from Hd,0 determined values by less than 0.001 Å for bond distances and 0.01° for bond angles. Plate (a) X ̃ 1 A state, experiment Ref. 38. Plate (b) A ̃ 1 A state, experiment Ref. 39.

FIG. 1.

Key geometrical parameters for phenol. Experimental values above Hd,0 determined values. Hd,s determined values differ from Hd,0 determined values by less than 0.001 Å for bond distances and 0.01° for bond angles. Plate (a) X ̃ 1 A state, experiment Ref. 38. Plate (b) A ̃ 1 A state, experiment Ref. 39.

Close modal
TABLE II.

G4 group.

E P(12) E P(12)
E C2 σva σvb
A1   
A2  −1  −1   
B1  −1  −1   
B2  −1  −1   
Csplanar(A     (A1 or B2
(A″)    −1    (A2 or B1
Csbijct(A     (A1 or B1
(A″)      −1  (A2 or B2
E P(12) E P(12)
E C2 σva σvb
A1   
A2  −1  −1   
B1  −1  −1   
B2  −1  −1   
Csplanar(A     (A1 or B2
(A″)    −1    (A2 or B1
Csbijct(A     (A1 or B1
(A″)      −1  (A2 or B2
a

Molecular plane.

b

Bijection plane.

FIG. 2.

The Hd,0 and Hd,s determined energies: plate (a) four diagonal diabatic matrix elements, and plate (b) four adiabatic state energies, along a path originating between the ground state minimum and the first excited state minimum of phenol, passing through the 21A − 31A and 11A − 21A minimum energy conical intersection points and terminating in the C6H5O + H dissociation channel at the ground state equilibrium structure of phenoxyl.

FIG. 2.

The Hd,0 and Hd,s determined energies: plate (a) four diagonal diabatic matrix elements, and plate (b) four adiabatic state energies, along a path originating between the ground state minimum and the first excited state minimum of phenol, passing through the 21A − 31A and 11A − 21A minimum energy conical intersection points and terminating in the C6H5O + H dissociation channel at the ground state equilibrium structure of phenoxyl.

Close modal

This symmetry analysis is also consistent with spectroscopic observations, where the low barrier to OH torsion requires phenol to be treated with G4 symmetry.23 See also Ref. 4. The symmetries of the ground (1A1) and first excited state (1B2) have been assigned experimentally24 and are in accord with the symmetry designations in Fig. 2(a).

In this section the construction of Hd,s is described. Only the 1B1 and 1B2 diabats will be shifted.

The shifts introduced into Hd,s reflect: (i) the energies of the ground state minimum and first excited state minimum of phenol, denoted Q phenol min ( X ̃ 1 A ) and Q phenol min ( A ̃ 1 A ) , respectively; (ii) the vertical excitation energy of the 31A state; and (iii) the energies of the ground and first excited state minima of phenoxyl denoted, Q phenoxyl min ( X ̃ 2 B 1 ) and Q phenoxyl min ( A ̃ 2 B 2 ) . These structures are reported in Figs. 1(a) and 1(b) and 3(a) and 3(b), respectively. In this work, when no confusion will result, we use the symbol Q a b (Z) to represent both the name of the structure and its geometry. Table III reports the Hd,0 determined energies of these structures and compares them with experimentally determined or inferred values. The corresponding harmonic frequencies are reported for phenol extrema in Table SI and phenoxyl extrema in Table SII in the supplementary material.25 

TABLE III.

Energetics of key structures in cm−1.

Hd,s Hd,0 Lit
Q phenol min ( X ̃ 1 A 0.0  0.0   
Q phenol min ( A ̃ 1 A 37 597.4  39 395.6   
Q phenoxyl min ( X ̃ 2 B 1 30 732.5  31 770.4   
Q phenoxyl min ( A ̃ 2 B 2 41 321.2  43 120.0   
Q phenol , O–H ts (21A)  42 749.8  44 031.0  41 459a 
Q phenol,tor ts (11A)  869.6  868.67   
Q phenol,tor ts (21A)  39 280.3  41 065.9   
Q phenol mex (11A − 21A)  33 500.6  34 417.5  33 877.b 
Q phenol mex (21A − 31A)  43 486.6  44 659.2  40 329.b,c 
D 0 ( H - phenoxyl ( X ̃ 2 B 1 ) )   27 888.9  28 920.4  30 005.4d 
T 0 ( A ̃ 1 A , Phenol )   36 405.7  38 176.6  36 350.5e 
E 3 ( Q phenol min ( X ̃ 1 A ) )   46 436.0  47 473.4  46 459.0f 
T 0 ( A ̃ 2 B 2 , Phenoxyl )   10 601.8  11 365.6  8549.9g 
Hd,s Hd,0 Lit
Q phenol min ( X ̃ 1 A 0.0  0.0   
Q phenol min ( A ̃ 1 A 37 597.4  39 395.6   
Q phenoxyl min ( X ̃ 2 B 1 30 732.5  31 770.4   
Q phenoxyl min ( A ̃ 2 B 2 41 321.2  43 120.0   
Q phenol , O–H ts (21A)  42 749.8  44 031.0  41 459a 
Q phenol,tor ts (11A)  869.6  868.67   
Q phenol,tor ts (21A)  39 280.3  41 065.9   
Q phenol mex (11A − 21A)  33 500.6  34 417.5  33 877.b 
Q phenol mex (21A − 31A)  43 486.6  44 659.2  40 329.b,c 
D 0 ( H - phenoxyl ( X ̃ 2 B 1 ) )   27 888.9  28 920.4  30 005.4d 
T 0 ( A ̃ 1 A , Phenol )   36 405.7  38 176.6  36 350.5e 
E 3 ( Q phenol min ( X ̃ 1 A ) )   46 436.0  47 473.4  46 459.0f 
T 0 ( A ̃ 2 B 2 , Phenoxyl )   10 601.8  11 365.6  8549.9g 
a

Reference 11.

b

Reference 3 includes zero point effects.

c

Reference 1.

d

Reference 4.

e

Reference 23.

f

Reference 26.

g

Reference 37.

From Table III it is seen that the product channel T 0 ( A ̃ 2 B 2 - phenoxyl ) is 11 365 (8549) cm−1 which is 2816 cm−1 greater than the experimental value given parenthetically. Similarly the T 0 ( Q phenol min ( A ̃ 1 A ) ) is 38 177 (36 350) cm−1 which is 1827 cm−1 greater than the experimental value given parenthetically. Finally the vertical excitation energy of the 31A state is 47 473.4 (46 459.0) cm−1 which is 1014 cm−1 higher than the CASPT2 estimate of Ref. 26 (their Table VIII, 5.76 eV) given in parentheses. See also Refs. 4 and 11.

Figs. 2(a) and 2(b) help put these data in perspective, plotting as dashed lines the Hd,0 determined H I , I d , 0 ( R ) and adiabatic energies Ea,I,(m)(R), I = 1, …, 4 along a segmented path connecting the Franck-Condon region which includes Q phenol min ( X ̃ 1 A ) and Q phenol min ( A ̃ 1 A ) , to the minimum energy points of conical intersection of the 21A − 31A and 11A − 21A states, Q phenol mex (21A − 31A) and Q phenol mex (11A − 21A) and finally to the product channel asymptote with phenoxyl in its equilibrium ground state geometry, Q phenoxyl min ( X ̃ 2 B 1 ) . From the energetics noted above and diabats pictured in Fig. 2(a) it is seen that the discrepancies can be reduced by shifting the 1B2 diabat down by 1800.1 cm−1, V 1 B 2 s = 1800 . 1 cm−1 (to correct T 0 ( Q phenol min ( A ̃ 1 A ) ) —a 1772 cm−1 improvement is found in Table III) and shifting the 1B1 diabat down by 1039.2 cm−1, V 1 B 1 s = 1039 . 2 cm 1 (to correct the vertical excitation energy of the 31A state—a 1037 cm−1 improvement is found in Table III). T 0 ( A ̃ 2 B 2 - phenoxyl ) will also improve by the difference in these shifts ∼760 cm−1 (764 cm−1 is reported in Table III). The changes in the electronic energies are mirrored in the zero point energy (ZPE) adjusted quantities, the T0s, since the differences in the Hd,0 and Hd,s determined ZPEs are small, for phenol( X ̃ ), −13 cm−1, phenol( A ̃ ), 14 cm−1, phenoxyl( X ̃ ), 10 cm−1, and phenoxyl( A ̃ ), 7 cm−1. See Tables S1and S2 of the supplementary material.25 Comparing Figs. 2(a) and 2(b) it is seen that the coupling of the diabats near R(O–H) ∼ 3.5. a.u. precludes shifting of the adiabatic states to correct energy errors.

In this work we have used experimental data to, as we confirm below, improve the quality of the Hd,0. However this analysis makes clear that in favorable circumstances such as these, it is possible to use single point energies obtained from the most accurate electronic structure methods available, to improve the quality of the global Hd,0. In this way accurate electronic structure methods for which derivative coupling information is not available can be used to improve coupled PESs for which derivative coupling data have been used as to advantage.

Hd,0 and Hd,s need not have the same domain of definition. The domain of definition for the shifted Hd,s is determined by running surface hopping trajectories on the shifted surfaces. If the trajectories indicate that more R(n) are needed, points are added to the current domain of definition and Hd,0 is redetermined. If necessary the Vs can be changed. The procedure is repeated until the K% threshold is reached on the shifted coupled surfaces. The revised domain of definition need only be used for Hd,s. However it is important to emphasize when used without the shift the revised Hd based on the extended domain of definition gives essentially the same results as those obtained with Hd,0, which is expected of a domain of definition. In this work, for simplicity the unshifted form of Hd, s is used for Hd,0.

Table I reports the key parameters for Hd,s. From this table it is seen that approximately 330 additional points (∼a 5% increase) were required to describe Hd.s using the same functional form, to approximately the same accuracy.

In this section the Hd,0 and Hd,s determined results are compared with each other and with experiment where available. In some cases comparisons with the results of Ref. 11, denoted YXZT below, the only other full 33 dimensional representation of phenol photodissociation, will be used. This section has two interrelated goals: (i) to demonstrate the quality of the representation and (ii) to show that no pathologies are introduced into Hd,s by the shifting procedure.

Figs. 1(a) and 1(b) report the key parameters of the Hd,s and Hd,0 determined Q phenol min ( X ̃ 1 A ) and Q phenol min ( A ̃ 1 A ) and compares them with the available experimental values. Cartesian coordinates for all structures described in this section are reported in the supplementary material.25 The Hd,s and Hd,0 determined structures in Figs. 1(a) and 1(b) are seen to be in good accord with each other, with bond distances differing by less than 0.001 Å and bond angles differing by less than 0.01°. Differences with the indicated experimental values are somewhat larger, with some bond distances differing by as much as 0.03 Å but most being a third of that, and bond angles differing by at most 1°. Table S1 of the supplementary material25 compares the Hd,0 and Hd,s determined harmonic frequencies with each other and the experimental values which are available for Q phenol min ( X ̃ 1 A ) . The Hd,0 and Hd,s determined frequencies are in excellent accord with each other. Most differences are 1-2 cm−1 with the few outliers showing differences of only 15-20 cm−1. When compared to the experimental values for Q phenol min ( X ̃ 1 A ) the computed harmonic values are too large (as expected) by <10% since anharmonic effects are neglected. Assessment of anharmonic effects will be the subject of future studies.

Table III reports the Hd,0 and Hd,s determined electronic energies for Q phenol min ( X ̃ 1 A ) , Q phenol min ( A ̃ 1 A ) as well as for three transition states, torsions on the 11A and 21A PESs, Q phenol,tor ts ( 1 1 A ) and Q phenol,tor ts ( 2 1 A ) and the transition state for O–H bond breaking on the 21A state, Q phenol , O H ts . The latter is discussed in detail below. The low torsional barriers on the 11A and 21A PESs necessitate the use of the G4 symmetry group in spectroscopic studies. A careful spectroscopic analysis of the torsional barrier on the 11A PES, which is of little consequence for the energies involved in the photodissociation considered here, is reported in Ref. 27 and suggests that our barrier is probably too low by ∼400 cm−1.

The Hd,s determined T 0 ( Q phenol min ( A ̃ 1 A ) ) and the vertical excitation energy of 31A state, quantities in the key Franck-Condon region, are within 55 cm−1 and 25 cm−1 of the literature values, respectively. As discussed in Sec. III, this good agreement reflects the shifts introduced in that section.

Figs. 3(a) and 3(b) repeat the above analysis for the asymptotic states, analyzing the structures of Q phenoxyl min ( X ̃ 2 B 1 ) and Q phenoxyl min ( A ̃ 2 B 2 ) . No experimental data are available for comparison. In each case the Hd,s and Hd,0 determined structures are indistinguishable; bond distances are within 0.001 Å and bond angles within 0.01°. The structures are largely equivalent to the CASPT2 structures of Ref. 28 with bond distances differing by at most than 0.03 Å and bond angles differing by at most 1.0°. Again the Hd,0 and Hd,s frequencies are in excellent agreement with most differences being 1-2 cm−1 and the few outliers showing differences of only 15-20 cm−1.

FIG. 3.

Key geometrical parameters for phenoxyl. Hd,0 determined. Plate (a) X ̃ 2 B 1 state. Plate (b) A ̃ 2 B 2 state. Hd,s determined results indistinguishable from Hd,0 determined values.

FIG. 3.

Key geometrical parameters for phenoxyl. Hd,0 determined. Plate (a) X ̃ 2 B 1 state. Plate (b) A ̃ 2 B 2 state. Hd,s determined results indistinguishable from Hd,0 determined values.

Close modal

The energies of Q phenoxyl min ( X ̃ 2 B 1 ) , Q phenoxyl min ( A ̃ 2 B 2 ) are reported in Table III, as is D 0 ( H - phenoxyl ( X ̃ 2 B 1 ) ) . The only deficiency in the asymptotic results is the energy of the X ̃ 2 B 1 state of phenoxyl which is predicted to be too low by ∼2000 cm−1, compared to the experimental values of D 0 ( H - phenoxyl ( X ̃ 2 B 1 ) ) and T 0 ( A ̃ 2 B 2 ) .

Figs. 4(a) and 4(b) extend the analysis to the minimum energy conical intersections of the 21A − 31A and 11A − 21A states, Q phenol mex (21A − 31A) and Q phenol mex (11A − 21A), respectively. For Q phenol mex (21A − 31A) the Hd,0 and Hd,s determined structures of the benzene ring are largely the same. The largest changes are observed for R(O–H) which is only ∼0.02 Å longer in the Hd,s determined structure. A similar situation is obtained for Q phenol mex (11A − 21A) with the change in R(O–H) being only ∼0.01 Å.

FIG. 4.

Structure of minimum energy conical intersections. Hd,0 determined parameters above Hd,s derived parameters. Hd,s shown if shown difference is >0.001 Å for distances, >0.01° for angles. Plate (a) 2,31A states intersection. Plate (b) 1,21A states intersection.

FIG. 4.

Structure of minimum energy conical intersections. Hd,0 determined parameters above Hd,s derived parameters. Hd,s shown if shown difference is >0.001 Å for distances, >0.01° for angles. Plate (a) 2,31A states intersection. Plate (b) 1,21A states intersection.

Close modal

The conical parameters ‖g‖, ‖h‖, and sx, sy, are reported in Table IV. Here g is the energy difference gradient vector;29h is the interstate coupling29 vector; sx and sy, which determine the tilt of the double cone,30 are the projections of the average energy gradient along the orthogonal g and h directions, respectively. As indicated by the comparisons in Table IV the g, h, and s vectors for Q phenol mex (11A − 21A) and for Q phenol mex (21A − 31A) are little affected by the shifts. Also note that the Q phenol mex (21A − 31A) and Q phenol mex (11A − 21A) conical intersections are tilted along the g axis with sx ∼ 0.5 ‖g‖. Correct description of the local topography, including the tilt of the double cone, as described by the conical parameters is the key to accurate nuclear dynamics in the vicinity of conical intersections.

TABLE IV.

Topographical parameters of conical intersections (sy = 0).

Q phenol mex (11A, 21A)
Hd,s : |g| = 0.125 933 9  |h| = 0.032 427 19  sx = 0.066 061 53 
Hd,0 : |g| = 0.127 453 1  |h| = 0.032 946 97  sx = 0.071 259 36 
Q phenol mex (21A, 31A) 
Hd,s : |g| = 0.129 753 6  |h| = 0.013 423 06  sx = 0.078 993 16 
Hd,0 : |g| = 0.122 603 9  |h| = 0.013 063 28  sx = 0.083 162 09 
Q phenol mex (11A, 21A)
Hd,s : |g| = 0.125 933 9  |h| = 0.032 427 19  sx = 0.066 061 53 
Hd,0 : |g| = 0.127 453 1  |h| = 0.032 946 97  sx = 0.071 259 36 
Q phenol mex (21A, 31A) 
Hd,s : |g| = 0.129 753 6  |h| = 0.013 423 06  sx = 0.078 993 16 
Hd,0 : |g| = 0.122 603 9  |h| = 0.013 063 28  sx = 0.083 162 09 

Table III reports the energies at Q phenol mex (21A − 31A) and Q phenol mex (11A − 21A) and compares them with experimentally inferred values. When compared to using Hd,0, the Hd,s determined the energy at Q phenol mex (21A − 31A) [ Q phenol mex (11A − 21A)] is lowered by 1173 [917] cm−1 improving the agreement with the estimates of Refs. 1 and 3 given in Table III. These changes are consistent with the 760 cm−1 differential in the 1B2 and 1B1 shifts discussed in Sec. III. Also note that the literature value for the energy of Q phenol mex (21A − 31A) includes ZPE effects. As a first approximation one could compute the ZPE of the standard Hessian obtained by averaging the Hessians of the 2,31A states centered at Q phenol mex (21A − 31A) in the orthogonal complement of the g and h modes. This suggests that the electronic energy reported in Table III over estimates the energy of Q phenol mex (21A − 31A). This point is addressed further in Sec. IV D where the topography in the vicinity of Q phenol mex (21A − 31A) is discussed.

These conical intersections were also determined in YXZT. For comparison it is useful to juxtapose the key parameters R(O–H) in Å and energy in cm−1 at the crossing. At Q phenol mex (21A − 31A), (R(O–H), energy) = (1.12,43 486) [(1.27, 43 153)] while for Q phenol mex (11A − 21A) (R(O–H), energy) = (1.76, 33 500) [(1.97, 33 635)] where the results of YXZT are given in square brackets. The bond lengths in YXZT are consistently longer than those reported here while the energies are quite similar.

The Q phenol , O H ts ( 2 1 A ) transition state is pictured in Fig. 5. Its role in phenol photodissociation and its relation to the 21A − 31A conical intersection has been discussed at length in YXZT. The structural relation is considered in Fig. 6, which compares g to the imaginary frequency of Q phenol , O H ts ( 2 1 A ) and reports h. Note that there are actually two transition states which are obtained largely by distortions along ±h from Q phenol mex (21A − 31A). In YXZT Q phenol , O – H ts ( 2 1 A ) was found to be geometrically close to, and energetically only somewhat below, the 21A − 31A conical intersection. It is seen from Table III, Fig. 6 and comparing Figs. 5 and 4(a) that this novel expectation is also realized in our calculation. While the effect of using Hd,s is to lower the energy of both Q phenol , O – H ts ( 2 1 A ) and Q phenol mex (21A − 31A), the subtle relationship is preserved. The Hd,s determined separation of the saddle point and minimum energy conical intersection is 737 cm−1, only slightly greater than the Hd,0 determined value of 628 cm−1.

FIG. 5.

21A saddle point for OH bond breaking. Hd,0 results above Hd,s. Hd,0 shown if difference is >0.001 Å for distances, >0.01° for angles.

FIG. 5.

21A saddle point for OH bond breaking. Hd,0 results above Hd,s. Hd,0 shown if difference is >0.001 Å for distances, >0.01° for angles.

Close modal
FIG. 6.

Comparison of Q phenol ts ( 2 1 A ) imaginary frequency and g and h vectors of Q phenol mex (21A − 31A). Note that Q phenol ts ( 2 1 A ) can be reached (in part) by motion along h.

FIG. 6.

Comparison of Q phenol ts ( 2 1 A ) imaginary frequency and g and h vectors of Q phenol mex (21A − 31A). Note that Q phenol ts ( 2 1 A ) can be reached (in part) by motion along h.

Close modal

Finally comparing our R(O–H) of 1.16 (1.33) Å and barrier height of Q phenol , O H ts ( 2 1 A ) with respect to Q phenol min ( A ̃ 1 A ) , of 0.64 (0.72) eV with those of YXZT given parenthetically, we see that as was the case with the 21A − 31A conical intersection, the energetics are similar but the R(O–H) of YXZT is appreciably longer.

Using the electronic energies in Table III and ZPEs from the supplementary material25 we find the barrier height of Q phenol , O H ts ( 2 1 A ) relative to Q phenol min ( X ̃ 1 A ) to be 5.026 eV. Reference 1 reported a distinct change in the kinetic energy release spectrum at excitation energy λphot = 248 nm (5.00 eV) and attributed it to access to the (21A − 31A) conical intersection seam. The current analysis suggests that Q phenol , O H ts ( 2 1 A ) may play a role in this observation. YXZT point out that tunneling through the Q phenol , O H ts ( 2 1 A ) barrier is expected to play an important role in the dynamics in this region. Consequently a (reduced dimensionality) fully quantum treatment is required. Such a reduced dimensionality fully quantum treatment is in progress.

Finally it is worth noting that the topography in the vicinity of Q phenol mex (21A − 31A) with a minimum, Q phenol min ( A ̃ 1 A ) , for slightly shorter R(O–H) (approximately along the negative g direction) and two saddle points Q phenol , O H ts ( 2 1 A ) symmetrically placed along h has the appearance of a Cs symmetry Jahn-Teller-like conical intersection where the brim of the Mexican hat is distorted for increasing g.

The good agreement between the Hd,s determined adiabatic excitation energy of the 21A state and the vertical excitation energy of the 32A state with the available literature values strongly support the accuracy of the Hd,s representation, particularly in the Franck-Condon region. The fact that the Hd,s predicted energy of Q phenol , O H ts ( 2 1 A ) (∼5 eV) and of the nearby conical intersection Q phenol mex (21A − 31A) are in good accord with the excitation energy, λphot = 248 nm (5.00 eV), corresponding to the increase in fast hydrogen in the experimental kinetic energy release (KER) spectrum of Ref. 1 lends further support to the accuracy of Hd,s in the Franck-Condon region.

In this section quasi-classical surface hopping trajectories using the ANT computer program31 based on Tully’s fewest switches algorithm32 are used to assess Hd,0 and Hd,s. All trajectories originate on S1. Propagations for greater than 10 ps encounter no pathologies using either Hd.0 or Hd,s. Two types of trajectories were used: (i) trajectories based on a S1 harmonic Wigner distribution, with all v = 0 denoted vOH = 0; and (ii) trajectories with 1 quantum in the O-H mode and all the rest of the modes in their ground states, denoted vOH = 1. In each case the total energy is 7000 cm−1 above the zero point energy, which is well above the 4635 (5152) cm−1 transition state barrier on the 21A PES. In this section results are presented as Hd,0 determined followed by Hd,s determined in parentheses. Two observations are interesting. For vOH = 0, approximately 34% (33%) of the trajectories have dissociated to phenoxyl and hydrogen within 10 ps. For vOH = 1, approximately 78% (76%) of the trajectories have dissociated within 10 ps. These results can be reasonably interpreted in terms of internal energy redistribution except for the fact these are classical trajectories and using 2000 cm−1 total energy 25% (20%) of trajectories have dissociated after 10 ps, which can only occur through ZPE violations.

The good agreement between the Hd,s and Hd,0 trajectory results lend further support the utility of Hd,s for detailed studies of the phenol dissociation dynamics.

Recently we reported a quasi-diabatic representation of the 1, 2, 3 1A states of phenol suitable for the description of its nonadiabatic photodissociation. The representation treats all 33 internal coordinates in an evenhanded manner and is fit solely to multireference single and double excitation configuration interaction wave function electronic structure data. That work established the accuracy of the fit, that is, its ability to reproduce the electronic structure data from which it was constructed. In this work we establish the accuracy of the representation, that is, its ability to describe the 1, 2, 3 1A states of phenol and their interactions. The representation is shown to provide an accurate description of the key Franck-Condon region where the initial nonadiabatic events transpire. In particular the error in T0 for the first excited state of phenol [T0(S1)] and the vertical excitation energy to the 31A state (S2) of phenol are each within 55 cm−1 of the literature values. The good agreement with the available data was achieved using a modest energetic shift of each of two diabats. These observations give us confidence in the predicted energies of the nearby S1-S2 minimum energy conical intersection and the associated S1 saddle point reported in Ref. 11.

Quasi-classical surface hopping trajectories based on Tully’s fewest switches algorithm demonstrated that no pathologies, holes in the representation, were introduced by the shifting procedure, although a slightly larger domain of definition was required. Quasi-classical surface hopping trajectories further document the viability of the Hd,s determined coupled potential energy surfaces for nonadiabatic dynamics. Trajectories initiated on S1 at energies both above and below the barrier (including zero point energy effects) at Q phenol , OH ts propagate for tens of picoseconds without encountering any surface anomalies, that is “holes” which can occur regions outside the domain of definition. This attests to the ability of the representation to describe the Franck-Condon regions where nonadiabatic effects are likely to be preeminent. Further long propagation times are relevant to the computational description of the transition from direct to statistical dynamics which is a problem of considerable current interest.33–36 

In future work Hd,s will be used to simulate phenol photodissociation using both full dimensional quasi classical surface hopping and reduced dimensional full quantum dynamics, tasks for which the present work shows it is quite well suited. Also in the future we will attempt to extend this fitting methodology with diabat shifting to treat even larger, more complex systems. In this regard it is useful to note the following. The use of analytic energy gradient and derivative coupling vectors data means that the amount of fitting data available at each point will increase with the size of the system and the number of electronic states. On the other hand, the system considered here, phenol, is a single dissociation channel problem. The data requirements and number of fitting parameters increase with the number of dissociative channels.

This work was supported by Department of Energy Basic Energy Sciences Grant No. DE-FG02-91ER14189 to D.R.Y. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

1.
M. G. D.
Nix
,
A. L.
Devine
,
B.
Cronin
,
R. N.
Dixon
, and
M. N. R.
Ashfold
,
J. Chem. Phys.
125
,
133318
(
2006
).
2.
O. P.
Vieuxmaire
,
Z.
Lan
,
A. L.
Sobolewski
, and
W.
Domcke
,
J. Chem. Phys.
129
,
224307
(
2008
).
3.
M. L.
Hause
,
Y. H.
Yoon
,
A. S.
Case
, and
F. F.
Crim
,
J. Chem. Phys.
128
,
104307
(
2008
).
4.
R. N.
Dixon
,
T. A. A.
Oliver
, and
M. N. R.
Ashfold
,
J. Chem. Phys.
134
,
194303
(
2011
).
5.
X.
Xu
,
K. R.
Yang
, and
D. G.
Truhlar
,
J. Chem. Theory Comput.
9
,
3612
(
2013
).
6.
G. M.
Roberts
,
A. S.
Chatterley
,
J. D.
Young
, and
V. G.
Stavros
,
J. Phys. Chem. Lett.
3
,
348
(
2012
).
7.
G. A.
Pino
,
A. N.
Oldani
,
E.
Marceca
,
M.
Fujii
,
S.
Ishiuchi
,
M.
Miyazaki
,
M.
Broquier
,
C.
Dedonder
, and
C.
Jouvet
,
J. Chem. Phys.
133
,
124313
(
2010
).
8.
M. N. R.
Ashfold
,
A. L.
Devine
,
R. N.
Dixon
,
G. A.
King
,
M. G. D.
Nix
, and
T. A. A.
Oliver
,
Proc. Natl. Acad. Sci. U. S. A.
105
,
12701
(
2008
).
9.
A.
Iqbal
,
M.
Cheung
,
M.
Nix
, and
V.
Stavros
,
J. Phys. Chem. A
113
,
8157
(
2009
).
10.
X.
Zhu
and
D. R.
Yarkony
,
J. Chem. Phys.
140
,
024112
(
2014
).
11.
K. R.
Yang
,
X.
Xu
,
J.
Zheng
, and
D. G.
Truhlar
,
Chem. Sci.
5
,
4661
(
2014
).
12.
S. G.
Ramesh
and
W.
Domcke
,
Faraday Discuss.
163
,
73
(
2013
).
13.
A.
Devine
,
M.
Nix
,
R.
Dixon
, and
M.
Ashfold
,
J. Phys. Chem. A
112
,
9563
(
2008
).
14.
G. A.
King
,
A. L.
Devine
,
M. G. D.
Nix
,
D. E.
Kelly
, and
M. N. R.
Ashfold
,
Phys. Chem. Chem. Phys.
10
,
6417
(
2008
).
15.
T. N. V.
Karsili
,
A. M.
Wenge
,
S. J.
Harris
,
D.
Murdock
,
J. N.
Harvey
,
R. N.
Dixon
, and
M. N. R.
Ashfold
,
Chem. Sci.
4
,
2434
(
2013
).
16.
M.
Nix
,
A.
Devine
,
R.
Dixon
, and
M.
Ashfold
,
Chem. Phys. Lett.
463
,
305
(
2008
).
17.
X.
Zhu
and
D. R.
Yarkony
,
J. Chem. Phys.
144
,
024105
(
2016
).
18.
H.
An
and
K.
Baeck
,
J. Phys. Chem. A
115
,
13309
(
2011
).
19.
T. N. V.
Karsili
,
A. M.
Wenge
,
B.
Marchetti
, and
M. N. R.
Ashfold
,
Phys. Chem. Chem. Phys.
16
,
588
(
2014
).
20.
Z.
Lan
,
W.
Domcke
,
V.
Vallet
,
A. L.
Sobolewski
, and
S.
Mahapatra
,
J. Chem. Phys.
122
,
224315
(
2005
).
21.
H. C.
Longuet-Higgins
,
Mol. Phys.
6
,
445
(
1963
).
22.
P. R.
Bunker
and
P.
Jensen
,
Molecular Symmetry and Spectroscopy
, 2nd ed. (
NRC Research Press
,
Ottawa
,
1998
).
23.
H. D.
Bist
,
J. C. D.
Brand
, and
D. R.
Williams
,
J. Mol. Spectrosc.
21
,
76
(
1966
).
24.
H. D.
Bist
,
J. C. D.
Brand
, and
D. R.
Williams
,
J. Mol. Spectrosc.
24
,
413
(
1967
).
25.
See supplementary material at http://dx.doi.org/10.1063/1.4944091 for tables.
26.
J.
Lorentzon
,
P.-Å.
Malmqvist
,
M.
Fülscher
, and
B. O.
Roos
,
Theor. Chem. Acc.
91
,
91
(
1995
).
27.
H. D.
Bist
,
J. C. D.
Brand
, and
D. R.
Williams
,
J. Mol. Spectrosc.
24
,
402
(
1967
).
28.
C.-W.
Cheng
,
Y.-P.
Lee
, and
H. A.
Witek
,
J. Phys. Chem. A
112
,
2648
(
2008
).
29.
D. R.
Yarkony
,
Acc. Chem. Res.
31
,
511
(
1998
).
30.
D. R.
Yarkony
,
J. Chem. Phys.
114
,
2601
(
2001
).
31.
Z. H.
Li
,
A. W.
Jasper
,
D. A.
Bonhommeau
,
R.
Valero
, and
D. G.
Truhlar
,
ANT A Computer Program
(
Univerisity of Minnesota
,
2009
).
32.
J. C.
Tully
,
J. Chem. Phys.
93
,
1061
(
1990
).
33.
H.
Tao
,
L.
Shen
,
M. H.
Kim
,
A. G.
Suits
, and
T. J.
Martinez
,
J. Chem. Phys.
134
,
054313
(
2011
).
34.
M.
Shapero
,
N. C.
Cole-Filipiak
,
C.
Haibach-Morris
, and
D. M.
Neumark
,
J. Phys. Chem. A
119
,
12349
(
2015
).
35.
N. C.
Cole-Filipiak
,
B.
Negru
,
G. M. P.
Just
,
D.
Park
, and
D. M.
Neumark
,
J. Chem. Phys.
138
,
054301
(
2013
).
36.
H.
Wang
,
S.
Yu
,
S.
Su
,
D.
Dai
,
K.
Yuan
, and
X.
Yang
,
J. Phys. Chem. A
119
,
11313
(
2015
).
37.
R. F.
Gunion
,
M. K.
Gilles
,
M. L.
Polak
, and
W. C.
Lineberger
,
Int. J. Mass Spectrom. Ion Processes
117
,
601
(
1992
).
38.
N. W.
Larsen
,
J. Mol. Struct.
51
,
175
(
1979
).
39.
D.
Spangenberg
,
P.
Imhof
, and
K.
Kleinermanns
,
Phys. Chem. Chem. Phys.
5
,
2505
(
2003
).

Supplementary Material