In a recent work we constructed a quasi-diabatic representation, **H**^{d}, of the 1, 2, 3^{1}A adiabatic states of phenol from high level multireference single and double excitation configuration interaction electronic structure data, energies, energy gradients, and derivative couplings. That **H**^{d} 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 **H**^{d} for describing phenol photodissociation. Additionally, we demonstrate that a modest energetic shift of two diabats yields a quantifiably more accurate **H**^{d} 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.

## I. INTRODUCTION

The photodissociation of phenol to produce H and the phenoxyl radical

is understood to involve a bright *S*_{1}(^{1}*π*^{*}*π*)←*S*_{0}(^{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, **H**^{d}, of the 1, 2, 3^{1}A 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 **H**^{d} 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 **H**^{d}, 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 **H**^{d} 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 **H**^{d}.

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.

## II. QUASI-DIABATIC REPRESENTATIONS

### A. Functional form of **H**^{d}

The quasi diabatic representation **H**^{d} has the form

In Eq. (2)**R** denotes the nuclear coordinates, $ P \u02c6 \alpha , \beta $ is a standard group theoretical projector for the complete nuclear permutation inversion group,^{21,22}*g*^{(l)} are the monomial functions, themselves products of basic functions, used to represent the matrix elements of **H**^{d}, and **B**^{u,v} is an *N*^{state} × *N*^{state} symmetric matrix with a 1 in the (*u*, *v*) and (*v*, *u*) elements and the remaining elements are 0. Here *N ^{state}* = 4. The procedure for fitting

**H**

^{d,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

**H**is given in Ref. 17. The

^{d}**V**

^{s}are empirical shifts whose determination in the present case is described in Sec. III.

**H**

^{d,0}is denoted the nascent representation and

**H**

^{d,s}is denoted the shifted representation. Generic representations are denoted

**H**

^{d}. Accompanying

**H**

^{d,w}(

*w*= s or 0) is an electronic Schrödinger equation

from which the energies *E*^{a.J,(m)}(**R**), energy gradients ∇*E*^{a,J(m)}(**R**), and derivative couplings are obtained. The superscript *m*, for model, indicates that the energy is **H**^{d} determined.

Table I reports the basic parameters of the fit, the number of linear parameters, *N ^{c}*; the number of fit data points,

*N*; the number of equations,

^{pt}*N*; 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

^{eq}**H**

^{d}, its construction, and its accuracy can be found in

*ZY*.

. | H^{d,0}
. | H^{d,s}
. |
---|---|---|

N^{c} | 28 716 | 28 716 |

N^{pt} | 7 379 | 7 709 |

N^{eq} | 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 |

. | H^{d,0}
. | H^{d,s}
. |
---|---|---|

N^{c} | 28 716 | 28 716 |

N^{pt} | 7 379 | 7 709 |

N^{eq} | 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 |

### B. Domain of definition

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

### C. Permutation inversion group symmetry

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 G_{4} permutation inversion group and the C_{2v} 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 C_{s} 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) C_{2v}: for C1O7H13 collinear and along the C4C1 axis (the highest energy structures); (ii) C_{s}-planar: for C1O7H13 planar and in the benzene plane (the lowest energy structures); and (iii) C_{s}-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 C_{2v}–C_{s} correlations. Calculations in the Franck-Condon region for C_{s}-planar structures demonstrate that the first four states carry ^{1}A′, ^{1}A′, ^{1}A″, ^{1}A″ symmetry, respectively. At structures with C_{s}-bijct symmetry the first four states carry ^{1}A′, ^{1}A″, ^{1}A′, ^{1}A″ symmetry. Therefore using Table II the lowest state carries ^{1}A_{1} symmetry and the remaining states carry ^{1}B_{2}, ^{1}B_{1}, and ^{1}A_{2} symmetry in G_{4}. A similar result is obtained using an energy optimized C_{2v} 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 C_{2v} point group symmetry in the dissociation channel shows that the four low-lying states of phenoxyl carry ^{1}B_{1}, ^{1}B_{2}, ^{1}A_{2}, and ^{1}A_{1}, symmetry, respectively. Thus it is necessary to include four states to account for the change in the ^{1}A_{1} state from the ground state in the Franck-Condon region to the 4^{1}A 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).

. | E . | P(12) . | E^{∗}
. | P(12)^{∗}
. | . |
---|---|---|---|---|---|

. | E . | C_{2}
. | σ_{v}^{a}
. | σ_{v}′^{b}
. | . |

A_{1} | 1 | 1 | 1 | 1 | |

A_{2} | 1 | 1 | −1 | −1 | |

B_{1} | 1 | −1 | −1 | 1 | |

B_{2} | 1 | −1 | 1 | −1 | |

C_{s}planar(A^{′}) | 1 | 1 | (A_{1} or B_{2}) | ||

(A″) | 1 | −1 | (A_{2} or B_{1}) | ||

C_{s}bijct(A^{′}) | 1 | 1 | (A_{1} or B_{1}) | ||

(A″) | 1 | −1 | (A_{2} or B_{2}) |

. | E . | P(12) . | E^{∗}
. | P(12)^{∗}
. | . |
---|---|---|---|---|---|

. | E . | C_{2}
. | σ_{v}^{a}
. | σ_{v}′^{b}
. | . |

A_{1} | 1 | 1 | 1 | 1 | |

A_{2} | 1 | 1 | −1 | −1 | |

B_{1} | 1 | −1 | −1 | 1 | |

B_{2} | 1 | −1 | 1 | −1 | |

C_{s}planar(A^{′}) | 1 | 1 | (A_{1} or B_{2}) | ||

(A″) | 1 | −1 | (A_{2} or B_{1}) | ||

C_{s}bijct(A^{′}) | 1 | 1 | (A_{1} or B_{1}) | ||

(A″) | 1 | −1 | (A_{2} or B_{2}) |

^{a}

Molecular plane.

^{b}

Bijection plane.

This symmetry analysis is also consistent with spectroscopic observations, where the low barrier to OH torsion requires phenol to be treated with G_{4} symmetry.^{23} See also Ref. 4. The symmetries of the ground (^{1}A_{1}) and first excited state (^{1}B_{2}) have been assigned experimentally^{24} and are in accord with the symmetry designations in Fig. 2(a).

## III. **H**^{d,s} REPRESENTATION

In this section the construction of **H**^{d,s} is described. Only the ^{1}B_{1} and ^{1}B_{2} diabats will be shifted.

### A. Determining the shifts of the diabatic potential energy surfaces

The shifts introduced into **H**^{d,s} reflect: (i) the energies of the ground state minimum and first excited state minimum of phenol, denoted $ Q phenol min ( X \u0303 1 A \u2032 ) $ and $ Q phenol min ( A \u0303 1 A \u2032 ) $, respectively; (ii) the vertical excitation energy of the 3^{1}A state; and (iii) the energies of the ground and first excited state minima of phenoxyl denoted, $ Q phenoxyl min ( X \u0303 \u2009 2 B 1 ) $ and $ Q phenoxyl min ( A \u0303 \u2009 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 **H**^{d,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}

. | H^{d,s}
. | H^{d,0}
. | Lit . |
---|---|---|---|

$ Q phenol min $ ($ X \u0303 1 A \u2032 $) | 0.0 | 0.0 | |

$ Q phenol min $ ($ A \u0303 1 A \u2033 $) | 37 597.4 | 39 395.6 | |

$ Q phenoxyl min $ ($ X \u0303 \u2009 2 B 1 $) | 30 732.5 | 31 770.4 | |

$ Q phenoxyl min $ ($ A \u0303 \u2009 2 B 2 $) | 41 321.2 | 43 120.0 | |

$ Q phenol , O\u2013H ts $ (2^{1}A) | 42 749.8 | 44 031.0 | 41 459^{a} |

$ Q phenol,tor ts $ (1^{1}A) | 869.6 | 868.67 | |

$ Q phenol,tor ts $ (2^{1}A) | 39 280.3 | 41 065.9 | |

$ Q phenol mex $ (1^{1}A − 2^{1}A) | 33 500.6 | 34 417.5 | 33 877.^{b} |

$ Q phenol mex $ (2^{1}A − 3^{1}A) | 43 486.6 | 44 659.2 | 40 329.^{b}^{,}^{c} |

$ D 0 ( H - phenoxyl ( X \u0303 \u2009 2 B 1 ) ) $ | 27 888.9 | 28 920.4 | 30 005.4^{d} |

$ T 0 ( A \u0303 1 A \u2032 , Phenol ) $ | 36 405.7 | 38 176.6 | 36 350.5^{e} |

$ E 3 ( Q phenol min ( X \u0303 1 A \u2032 ) ) $ | 46 436.0 | 47 473.4 | 46 459.0^{f} |

$ T 0 ( A \u0303 \u2009 2 B 2 , Phenoxyl ) $ | 10 601.8 | 11 365.6 | 8549.9^{g} |

. | H^{d,s}
. | H^{d,0}
. | Lit . |
---|---|---|---|

$ Q phenol min $ ($ X \u0303 1 A \u2032 $) | 0.0 | 0.0 | |

$ Q phenol min $ ($ A \u0303 1 A \u2033 $) | 37 597.4 | 39 395.6 | |

$ Q phenoxyl min $ ($ X \u0303 \u2009 2 B 1 $) | 30 732.5 | 31 770.4 | |

$ Q phenoxyl min $ ($ A \u0303 \u2009 2 B 2 $) | 41 321.2 | 43 120.0 | |

$ Q phenol , O\u2013H ts $ (2^{1}A) | 42 749.8 | 44 031.0 | 41 459^{a} |

$ Q phenol,tor ts $ (1^{1}A) | 869.6 | 868.67 | |

$ Q phenol,tor ts $ (2^{1}A) | 39 280.3 | 41 065.9 | |

$ Q phenol mex $ (1^{1}A − 2^{1}A) | 33 500.6 | 34 417.5 | 33 877.^{b} |

$ Q phenol mex $ (2^{1}A − 3^{1}A) | 43 486.6 | 44 659.2 | 40 329.^{b}^{,}^{c} |

$ D 0 ( H - phenoxyl ( X \u0303 \u2009 2 B 1 ) ) $ | 27 888.9 | 28 920.4 | 30 005.4^{d} |

$ T 0 ( A \u0303 1 A \u2032 , Phenol ) $ | 36 405.7 | 38 176.6 | 36 350.5^{e} |

$ E 3 ( Q phenol min ( X \u0303 1 A \u2032 ) ) $ | 46 436.0 | 47 473.4 | 46 459.0^{f} |

$ T 0 ( A \u0303 \u2009 2 B 2 , Phenoxyl ) $ | 10 601.8 | 11 365.6 | 8549.9^{g} |

From Table III it is seen that the product channel $ T 0 ( A \u0303 \u2009 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 \u0303 1 A \u2032 ) ) $ 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 3^{1}A 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 **H**^{d,0} determined $ H I , I d , 0 ( R ) $ and adiabatic energies *E*^{a,I,(m)}(**R**), *I* = 1, …, 4 along a segmented path connecting the Franck-Condon region which includes $ Q phenol min ( X \u0303 1 A \u2032 ) $ and $ Q phenol min ( A \u0303 1 A \u2032 ) $, to the minimum energy points of conical intersection of the 2^{1}A − 3^{1}A and 1^{1}A − 2^{1}A states, $ Q phenol mex $ (2^{1}A − 3^{1}A) and $ Q phenol mex $ (1^{1}A − 2^{1}A) and finally to the product channel asymptote with phenoxyl in its equilibrium ground state geometry, $ Q phenoxyl min ( X \u0303 \u2009 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 ^{1}B_{2} diabat down by 1800.1 cm^{−1}, $ V 1 B 2 s =\u22121800.1$ cm^{−1} (to correct $ T 0 ( Q phenol min ( A \u0303 1 A \u2032 ) ) $—a 1772 cm^{−1} improvement is found in Table III) and shifting the ^{1}B_{1} diabat down by 1039.2 cm^{−1}, $ V 1 B 1 s =\u22121039.2 cm \u2212 1 $ (to correct the vertical excitation energy of the 3^{1}A state—a 1037 cm^{−1} improvement is found in Table III). $ T 0 ( A \u0303 \u2009 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 T_{0}s, since the differences in the **H**^{d,0} and **H**^{d,s} determined ZPEs are small, for phenol($ X \u0303 $), −13 cm^{−1}, phenol($ A \u0303 $), 14 cm^{−1}, phenoxyl($ X \u0303 $), 10 cm^{−1}, and phenoxyl($ A \u0303 $), 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 **H**^{d,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 **H**^{d,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.

### B. Domain of definition

**H**^{d,0} and **H**^{d,s} need not have the same domain of definition. The domain of definition for the shifted **H**^{d,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 **H**^{d,0} is redetermined. If necessary the V^{s} 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 **H**^{d,s}. However it is important to emphasize when used without the shift the revised **H**^{d} based on the extended domain of definition gives essentially the same results as those obtained with **H**^{d,0}, which is expected of a domain of definition. In this work, for simplicity the unshifted form of **H**^{d, s} is used for **H**^{d,0}.

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

## IV. COMPARISON OF **H**^{d,s} AND **H**^{d,0} REPRESENTATIONS

In this section the **H**^{d,0} and **H**^{d,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 **H**^{d,s} by the shifting procedure.

### A. Franck-Condon region

Figs. 1(a) and 1(b) report the key parameters of the **H**^{d,s} and **H**^{d,0} determined $ Q phenol min ( X \u0303 1 A \u2032 ) $ and $ Q phenol min ( A \u0303 1 A \u2032 ) $ 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 **H**^{d,s} and **H**^{d,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 material^{25} compares the **H**^{d,0} and **H**^{d,s} determined harmonic frequencies with each other and the experimental values which are available for $ Q phenol min ( X \u0303 1 A \u2032 ) $. The **H**^{d,0} and **H**^{d,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 \u0303 1 A \u2032 ) $ 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 **H**^{d,0} and **H**^{d,s} determined electronic energies for $ Q phenol min ( X \u0303 1 A \u2032 ) $, $ Q phenol min ( A \u0303 1 A \u2032 ) $ as well as for three transition states, torsions on the 1^{1}A and 2^{1}A 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 2^{1}A state, $ Q phenol , O \u2212 H ts $. The latter is discussed in detail below. The low torsional barriers on the 1^{1}A and 2^{1}A PESs necessitate the use of the G_{4} symmetry group in spectroscopic studies. A careful spectroscopic analysis of the torsional barrier on the 1^{1}A 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 **H**^{d,s} determined $ T 0 ( Q phenol min ( A \u0303 1 A \u2032 ) ) $ and the vertical excitation energy of 3^{1}A 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.

### B. Asymptotic region

Figs. 3(a) and 3(b) repeat the above analysis for the asymptotic states, analyzing the structures of $ Q phenoxyl min ( X \u0303 \u2009 2 B 1 ) $ and $ Q phenoxyl min ( A \u0303 \u2009 2 B 2 ) $. No experimental data are available for comparison. In each case the **H**^{d,s} and **H**^{d,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 **H**^{d,0} and **H**^{d,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}.

The energies of $ Q phenoxyl min ( X \u0303 \u2009 2 B 1 ) $, $ Q phenoxyl min ( A \u0303 \u2009 2 B 2 ) $ are reported in Table III, as is $ D 0 ( H - phenoxyl ( X \u0303 \u2009 2 B 1 ) ) $. The only deficiency in the asymptotic results is the energy of the $ X \u0303 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 \u0303 \u2009 2 B 1 ) ) $ and $ T 0 ( A \u0303 2 B 2 ) $.

### C. Minimum energy conical intersections

Figs. 4(a) and 4(b) extend the analysis to the minimum energy conical intersections of the 2^{1}A − 3^{1}A and 1^{1}A − 2^{1}A states, $ Q phenol mex $ (2^{1}A − 3^{1}A) and $ Q phenol mex $ (1^{1}A − 2^{1}A), respectively. For $ Q phenol mex $ (2^{1}A − 3^{1}A) the **H**^{d,0} and **H**^{d,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 **H**^{d,s} determined structure. A similar situation is obtained for $ Q phenol mex $ (1^{1}A − 2^{1}A) with the change in R(O–H) being only ∼0.01 Å.

The conical parameters ‖**g**‖, ‖**h**‖, and s_{x}, s_{y}, are reported in Table IV. Here **g** is the energy difference gradient vector;^{29} **h** is the interstate coupling^{29} vector; s_{x} and s_{y}, 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 $ (1^{1}A − 2^{1}A) and for $ Q phenol mex $ (2^{1}A − 3^{1}A) are little affected by the shifts. Also note that the $ Q phenol mex $ (2^{1}A − 3^{1}A) and $ Q phenol mex $ (1^{1}A − 2^{1}A) conical intersections are tilted along the **g** axis with s_{x} ∼ 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.

$ Q phenol mex $ (1^{1}A, 2^{1}A)
. | ||
---|---|---|

H^{d,s} : |g| = 0.125 933 9 | |h| = 0.032 427 19 | s_{x} = 0.066 061 53 |

H^{d,0} : |g| = 0.127 453 1 | |h| = 0.032 946 97 | s_{x} = 0.071 259 36 |

$ Q phenol mex $ (2^{1}A, 3^{1}A) | ||

H^{d,s} : |g| = 0.129 753 6 | |h| = 0.013 423 06 | s_{x} = 0.078 993 16 |

H^{d,0} : |g| = 0.122 603 9 | |h| = 0.013 063 28 | s_{x} = 0.083 162 09 |

$ Q phenol mex $ (1^{1}A, 2^{1}A)
. | ||
---|---|---|

H^{d,s} : |g| = 0.125 933 9 | |h| = 0.032 427 19 | s_{x} = 0.066 061 53 |

H^{d,0} : |g| = 0.127 453 1 | |h| = 0.032 946 97 | s_{x} = 0.071 259 36 |

$ Q phenol mex $ (2^{1}A, 3^{1}A) | ||

H^{d,s} : |g| = 0.129 753 6 | |h| = 0.013 423 06 | s_{x} = 0.078 993 16 |

H^{d,0} : |g| = 0.122 603 9 | |h| = 0.013 063 28 | s_{x} = 0.083 162 09 |

Table III reports the energies at $ Q phenol mex $ (2^{1}A − 3^{1}A) and $ Q phenol mex $ (1^{1}A − 2^{1}A) and compares them with experimentally inferred values. When compared to using **H**^{d,0}, the **H**^{d,s} determined the energy at $ Q phenol mex $ (2^{1}A − 3^{1}A) [$ Q phenol mex $ (1^{1}A − 2^{1}A)] 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 ^{1}B_{2} and ^{1}B_{1} shifts discussed in Sec. III. Also note that the literature value for the energy of $ Q phenol mex $ (2^{1}A − 3^{1}A) includes ZPE effects. As a first approximation one could compute the ZPE of the standard Hessian obtained by averaging the Hessians of the 2,3^{1}A states centered at $ Q phenol mex $ (2^{1}A − 3^{1}A) 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 $ (2^{1}A − 3^{1}A). This point is addressed further in Sec. IV D where the topography in the vicinity of $ Q phenol mex $ (2^{1}A − 3^{1}A) 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 $ (2^{1}A − 3^{1}A), (R(O–H), energy) = (1.12,43 486) [(1.27, 43 153)] while for $ Q phenol mex $ (1^{1}A − 2^{1}A) (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.

### D. The 2^{1}A transition state for O–H dissociation

The $ Q phenol , O \u2212 H ts ( 2 1 A ) $ transition state is pictured in Fig. 5. Its role in phenol photodissociation and its relation to the 2^{1}A − 3^{1}A 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 \u2212 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 $ (2^{1}A − 3^{1}A). In *YXZT* $ Q phenol , O \u2013 H ts ( 2 1 A ) $ was found to be geometrically close to, and energetically only somewhat below, the 2^{1}A − 3^{1}A 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 **H**^{d,s} is to lower the energy of both $ Q phenol , O \u2013 H ts ( 2 1 A ) $ and $ Q phenol mex $ (2^{1}A − 3^{1}A), the subtle relationship is preserved. The **H**^{d,s} determined separation of the saddle point and minimum energy conical intersection is 737 cm^{−1}, only slightly greater than the **H**^{d,0} determined value of 628 cm^{−1}.

Finally comparing our R(O–H) of 1.16 (1.33) Å and barrier height of $ Q phenol , O \u2212 H ts ( 2 1 A ) $ with respect to $ Q phenol min ( A \u0303 1 A \u2032 ) $, of 0.64 (0.72) eV with those of *YXZT* given parenthetically, we see that as was the case with the 2^{1}A − 3^{1}A 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 material^{25} we find the barrier height of $ Q phenol , O \u2212 H ts ( 2 1 A ) $ relative to $ Q phenol min ( X \u0303 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 (2^{1}A − 3^{1}A) conical intersection seam. The current analysis suggests that $ Q phenol , O \u2212 H ts ( 2 1 A ) $ may play a role in this observation. *YXZT* point out that tunneling through the $ Q phenol , O \u2212 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 $ (2^{1}A − 3^{1}A) with a minimum, $ Q phenol min ( A \u0303 1 A \u2032 ) $, for slightly shorter R(O–H) (approximately along the negative **g** direction) and two saddle points $ Q phenol , O \u2212 H ts ( 2 1 A ) $ symmetrically placed along **h** has the appearance of a C_{s} symmetry Jahn-Teller-like conical intersection where the brim of the Mexican hat is distorted for increasing **g**.

### E. Accuracy

The good agreement between the **H**^{d,s} determined adiabatic excitation energy of the 2^{1}A state and the vertical excitation energy of the 3^{2}A state with the available literature values strongly support the accuracy of the **H**^{d,s} representation, particularly in the Franck-Condon region. The fact that the **H**^{d,s} predicted energy of $ Q phenol , O \u2212 H ts ( 2 1 A ) $ (∼5 eV) and of the nearby conical intersection $ Q phenol mex $ (2^{1}A − 3^{1}A) 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 **H**^{d,s} in the Franck-Condon region.

### F. Nonadiabatic dynamics

In this section quasi-classical surface hopping trajectories using the ANT computer program^{31} based on Tully’s fewest switches algorithm^{32} are used to assess **H**^{d,0} and **H**^{d,s}. All trajectories originate on S_{1}. Propagations for greater than 10 ps encounter no pathologies using either **H**^{d.0} or **H**^{d,s}. Two types of trajectories were used: (i) trajectories based on a S_{1} harmonic Wigner distribution, with all **v** = 0 denoted v_{OH} = 0; and (ii) trajectories with 1 quantum in the O-H mode and all the rest of the modes in their ground states, denoted v_{OH} = 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 2^{1}A PES. In this section results are presented as **H**^{d,0} determined followed by H^{d,s} determined in parentheses. Two observations are interesting. For v_{OH} = 0, approximately 34% (33%) of the trajectories have dissociated to phenoxyl and hydrogen within 10 ps. For v_{OH} = 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 **H**^{d,s} and **H**^{d,0} trajectory results lend further support the utility of **H**^{d,s} for detailed studies of the phenol dissociation dynamics.

## V. SUMMARY AND CONCLUSIONS

Recently we reported a quasi-diabatic representation of the 1, 2, 3 ^{1}A 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 ^{1}A 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 T_{0} for the first excited state of phenol [T_{0}(S_{1})] and the vertical excitation energy to the 3^{1}A state (S_{2}) 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 S_{1}-S_{2} minimum energy conical intersection and the associated S_{1} 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 **H**^{d,s} determined coupled potential energy surfaces for nonadiabatic dynamics. Trajectories initiated on S_{1} 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 **H**^{d,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.

## Acknowledgments

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.