Quantum mechanical (QM) + molecular mechanics (MM) models are developed to represent potential energy surfaces (PESs) for the HBr^{+} + CO_{2} → Br + HOCO^{+} reaction with HBr^{+} in the ^{2}Π_{3/2} and ^{2}Π_{1/2} spin-orbit states. The QM component is the spin-free PES and spin-orbit coupling for each state is represented by a MM-like analytic potential fit to spin-orbit electronic structure calculations. Coupled-cluster single double and perturbative triple excitation (CCSD(T)) calculations are performed to obtain “benchmark” reaction energies without spin-orbit coupling. With zero-point energies removed, the “experimental” reaction energy is 44 ± 5 meV for HBr^{+}(^{2}Π_{3/2}) + CO_{2} → Br(^{2}P_{3/2}) + HOCO^{+}, while the CCSD(T) value with spin-orbit effects included is 87 meV. Electronic structure calculations were performed to determine properties of the BrHOCO^{+} reaction intermediate and [HBr⋯OCO]^{+} van der Waals intermediate. The results of different electronic structure methods were compared with those obtained with CCSD(T), and UMP2/cc-pVTZ/PP was found to be a practical and accurate QM method to use in QM/MM direct dynamics simulations. The spin-orbit coupling calculations show that the spin-free QM PES gives a quite good representation of the shape of the PES originated by ^{2}Π_{3/2}HBr^{+}. This is also the case for the reactant region of the PES for ^{2}Π_{1/2} HBr^{+}, but spin-orbit coupling effects are important for the exit-channel region of this PES. A MM model was developed to represent these effects, which were combined with the spin-free QM PES.

## I. INTRODUCTION

There is considerable experimental interest in the chemical dynamics of ion-molecule reactions with the reactants in specific quantum states.^{1–7} Viggiano and co-workers^{1} have investigated the influence of vibrational and rotational energies on the rates of ion-molecule reactions. For exothermic reactions, they found that rotational energy has a negligible influence on the reaction efficiency, while rotational energy increases the efficiency for endothermic reactions. More detailed information regarding the role of different types of energy on the reaction dynamics is obtained by studying state-selected molecular ions.^{2} Anderson and co-workers^{3–5} have been particularly interested in the effects of reactant vibrational excitation on the dynamics of ion-molecule reactions.

In recent research, Paetow *et al.*^{6,7} used a guided ion beam apparatus to measure the rate constants for HBr^{+} and DBr^{+}, in the ^{2}Π_{3/2} and ^{2}Π_{1/2} spin-orbit (SO) states, reacting with CO_{2} to form Br + HOCO^{+}/DOCO^{+}. The mean rotational energy of the HBr^{+} and DBr^{+} ions was varied and it was found that the rate constant decreased with increase in rotational energy. Similarly, the rate constants decreased with increase in reactant collision energy. These energy and state specific effects were found for both the endothermic reaction with HBr^{+}(DBr^{+}) in the ^{2}Π_{3/2} state and the exothermic reaction with these ions in the ^{2}Π_{1/2} state. The potential energy surface (PES) for these proton and deuteron transfer reactions has been investigated at the UMP2 and coupled-cluster single double and perturbative triple excitation (CCSD(T)) levels of theory,^{6,8} without the inclusion of SO coupling,^{9} and the resulting energies are summarized in Table I. An important feature of the PES is the reaction intermediate BrHOCO^{+} (BrDOCO^{+}). The 0 K experimental reaction energetics^{10,11} are depicted in Figure 1.

Theory . | E^{QM}(BrHOCO^{+})
. | ΔE_{r}^{QM}
. |
---|---|---|

UMP2/DZP^{b} | −706(−580) | −132(25) |

UMP2/TZ2P^{c} | −840(−729) | −114(41) |

CCSD(T)/TZ2P^{b} | −135(27) | |

CCSD(T)//PMP2/TZ2P^{b} | −824(−713) | −133(22) |

CCSD(T)/cc-pVTZ/SDB | −42 | |

CCSD(T)/cc-pVQZ/SDB | −6 | |

CCSD(T)/aug-cc-pVTZ/SDB | 86 | |

CCSD(T)/aug-cc-pVQZ/SDB | 77 | |

CCSD(T)/cc-pVDZ/PP | −668 | 16 |

CCSD(T)/cc-pVTZ/PP | −698 | 67 |

CCSD(T)/cc-pVQZ/PP | 73 | |

CCSD(T)/CBS/cc-pVXZ/PP | 75 | |

CCSD(T)/aug-cc-pVDZ/PP | 152 | |

CCSD(T)/aug-cc-pVTZ/PP | 161 | |

CCSD(T)/aug-cc-pVQZ/PP | 102 | |

CCSD(T)/CBS/aug-cc-pVXZ/PP | 62 |

Theory . | E^{QM}(BrHOCO^{+})
. | ΔE_{r}^{QM}
. |
---|---|---|

UMP2/DZP^{b} | −706(−580) | −132(25) |

UMP2/TZ2P^{c} | −840(−729) | −114(41) |

CCSD(T)/TZ2P^{b} | −135(27) | |

CCSD(T)//PMP2/TZ2P^{b} | −824(−713) | −133(22) |

CCSD(T)/cc-pVTZ/SDB | −42 | |

CCSD(T)/cc-pVQZ/SDB | −6 | |

CCSD(T)/aug-cc-pVTZ/SDB | 86 | |

CCSD(T)/aug-cc-pVQZ/SDB | 77 | |

CCSD(T)/cc-pVDZ/PP | −668 | 16 |

CCSD(T)/cc-pVTZ/PP | −698 | 67 |

CCSD(T)/cc-pVQZ/PP | 73 | |

CCSD(T)/CBS/cc-pVXZ/PP | 75 | |

CCSD(T)/aug-cc-pVDZ/PP | 152 | |

CCSD(T)/aug-cc-pVTZ/PP | 161 | |

CCSD(T)/aug-cc-pVQZ/PP | 102 | |

CCSD(T)/CBS/aug-cc-pVXZ/PP | 62 |

^{a}

Energies are in meV and are with respect to zero of energy for HBr^{+} + CO_{2}. Zero-point energies are not included in E^{QM}(BrHOCO^{+}), the energy of the BrHOCO^{+} intermediate with respect to reactants, and ΔE_{r}^{QM}, the reaction energy. Harmonic zero-point corrected energy are included in parenthesis.

Classical trajectory simulations have proven very important for interpreting experimental studies of state-selected ion-molecule reactions and determining atomistic details of their chemical dynamics. Of particular interest is the manner in which vibrational excitation of the reactants may enhance the reaction rate, which has been investigated for the reaction of H_{2}CO^{+} with D_{2} and CD_{4},^{13,14} and the NO_{2}^{+} + C_{2}H_{2}, and C_{2}H_{2}^{+} + CH_{4} reactions.^{15,16} The possibility of establishing of “Polanyi Rules” for polyatomic reactions has been investigated.^{13} Trajectories have also been used to investigate the role of state-specific vibrational excitation on collision-induced dissociation.^{17}

In the work reported here, SO coupling calculations^{9} are performed to derive PESs for the HBr^{+} + CO_{2} → Br + HOCO^{+} reaction with HBr^{+} in the ^{2}Π_{3/2} and ^{2}Π_{1/2} spin-orbit states, which in the following are denoted as ^{2}Π_{3/2} PES and ^{2}Π_{1/2} PES. An electronic structure quantum mechanical (QM) theory is used to represent the “average” ^{2}Π PES without spin-orbit coupling. Analytic molecular mechanics (MM)-like potential energy functions, which are fit to the spin-orbit coupling calculations, are then added to this QM component to give QM + MM^{18} PESs with spin-orbit coupling. These potential energy surfaces may be used in future QM + MM direct dynamics simulations^{18} to study the HBr^{+} + CO_{2} → Br + HOCO^{+} reaction with HBr^{+} in the ^{2}Π_{3/2} and ^{2}Π_{1/2} spin-orbit states. The reason for resorting to an analytic representation of the spin-orbit contribution to the PES is that different electronic structure methods are best suited for the spin-free and for the spin-orbit calculations, which would make the direct evaluation of the spin-orbit effect quite expensive. Moreover, most electronic structure packages cannot compute analytic gradients for the spin-orbit corrected PESs, as needed in direct dynamics simulations.^{19,20}

## II. COMPUTATIONAL METHODS

Without SO coupling HBr^{+} has a four-fold degenerate ^{2}Π state. With SO coupling, this state splits into a doubly degenerate ^{2}Π_{3/2} SO ground state and a doubly degenerate ^{2}Π_{1/2} SO state which is 329 meV higher in energy.^{12,21} Also considered for the SO coupling calculations described below is the HBr^{+} doubly degenerate ^{2}Σ_{1/2} state. The ^{2}P_{3/2} ground electronic state of the Br atom is four-fold degenerate, and both the ^{2}Π_{3/2} and ^{2}Π_{1/2} states of the HBr^{+} reactant correlate with this product ^{2}P_{3/2} state. The ^{2}Σ_{1/2} state of HBr^{+} correlates with the doubly degenerate ^{2}P_{1/2} excited state of the Br atom. The ^{2}P_{3/2} and ^{2}P_{1/2} states of Br are separated by 457 meV. In a BrH^{+}⋯OCO collinear arrangement and without SO coupling, the ground electronic state remains ^{2}Π as for the reactants and is four-fold degenerate. For bent planar geometries, these degenerate states split into a ^{2}A″ term and a ^{2}A′ term, still very close in energy.

### A. Calculations without spin-orbit coupling

In previous work, MP2 and CCSD(T) electronic structure calculations were performed, without SO coupling, for the HBr^{+} + CO_{2} → HOCO^{+} + Br ^{2}A″ ground state PES.^{6,8} For the work reported here, these calculations were supplemented with additional electronic structure calculations performed with the NWChem computer program.^{22}

CCSD(T)^{23} calculations, with both augmented and non-augmented correlation consistent double, triple, and quadruple zeta basis sets,^{24} were used to calculate the HBr^{+} + CO_{2} → HOCO^{+} + Br ^{2}A″ heat of reaction without SO coupling. The complete basis set (CBS) limit for these calculations was obtained using the formula of Peterson *et al.*,^{25,26} i.e.,

where *n* = 2, 3, and 4, represent X = D, T, and Q, respectively, for the cc-pVXZ basis sets.

A large number of basis functions are needed to accurately describe all the electrons for heavy atoms with many electrons, and the size of the basis set becomes important in treating both scalar and spin-orbit relativistic effects. Since the core electrons do not play an important role in chemical reactions, methods have been developed^{27,28} to replace the core electrons by effective core potentials (ECPs) or pseudopotentials (PPs).^{29–31} In the research presented here, the energy-consistent Stuttgart-Dresden-Bonn (SDB)^{32} PP and small-core PPs^{33} were used for the many electron Br atom for both accuracy and efficiency.

To compare with the benchmark CCSD(T)/CBS calculations, the reaction potential energy profile was also calculated with the three widely used density functional theory (DFT) functionals B3LYP,^{34} PBE0,^{35} and Becke98,^{36} and unrestricted MP2.^{37} A variety of basis sets were used for these calculations with the goal of ascertaining whether the reaction potential energy profile may be accurately represented by a single reference method with low computational cost.

For a collinear geometry, the ground state for the HBr^{+} + CO_{2} → HOCO^{+} + Br reaction is identified as ^{2}Π, which splits into ^{2}A″ and ^{2}A′ levels for bent planar geometries at which the interaction between the Br atom and CO_{2} is not negligible. To study this splitting, two distinct ROHF calculations were performed followed by MP2 calculations. These ROHF-MP2 calculations for the ^{2}A′ state had occupation …(21*a*′)^{2}(7*a*″)^{2}(22*a*′)^{1} and for the ^{2}A″ state had occupation …(21*a*′)^{2}(7*a*″)^{1}(22*a*′)^{2}.

### B. Anharmonic zero-point energy corrections without spin-orbit coupling

The experimental heats of formation for the reactants and products are CO_{2}, −393.107 ± 0.014 kJ/mol; HBr^{+}, 1097.71 ± 0.14 kJ/mol; HOCO^{+}, 600.80 ± 0.45 kJ/mol; and Br, 117.92 ± 0.06 kJ/mol.^{38} The resulting 0 K Δ*H* = Δ*E* for the HBr^{+}(^{2}Π_{3/2}) + CO_{2} → HOCO^{+} + Br(^{2}P_{3/2}) reaction is 14.12 ± 0.47 kJ/mol (146 ± 5 meV). To compare with the *ab initio* reaction energetics, an experimental value for Δ*H* = Δ*E* is required which does not include zero-point energies (ZPEs) for the reactants and products. The following procedures were used to remove ZPE from the experimental 0 K heat of reaction.

To second-order^{39} and also for the Morse potential function,^{40} the vibrational energy levels for a diatomic molecule are given by

where *v _{e}* is the harmonic frequency and

*x*is the anharmonic correction term. These parameters are known for HBr

_{e}^{+}and are

*v*= 2441.5 cm

_{e}^{−1}and

*x*= 47.4.

_{e}^{41}The HBr

^{+}ZPE is (

*hv*/2 −

_{e}*hx*/4) = 14.74 kJ/mol.

_{e}The vibrational energy levels for CO_{2}, given by second-order perturbation theory,^{39} are

which is similar to Eq. (2) with the additional parameters *d* the degeneracy of the vibrational levels, *g*_{22} the vibrational angular momentum anharmonicity constant, and *l* the vibrational angular momentum. The parameters for the CO_{2} vibrational energy levels are (all in cm^{−1}):^{39} *v*_{1,e} = 1354.0, *v*_{2,e} = 673.2, *v*_{3,e} = 2396.3, *x*_{11} = − 2.9, *x*_{22} = + 1.1, *x*_{33} = − 12.5, *x*_{12} = − 3.6, *x*_{13} = − 19.7, *x*_{23} = − 12.4, and *g*_{22} = − 0.9. For the ZPE level, the *n _{i}* and

*l*are zero, so that the CO

_{2}ZPE is 30.33 kJ/mol.

The vibrational energy levels for the product molecule HOCO^{+} have not been measured, and thus, experimental values may not be used to determine the ZPE for this molecule. Therefore, a well-tested electronic structure theory approach, with a scale factor for the anharmonic ZPE,^{42–46} was used. The specific method used is CCSD(T)/cc-pVTZ and the same optimized scale factor for each mode, determined by a least-squares minimization of the residuals between the scaled and experimental ZPEs for a chosen molecular database. The resulting ZPE for HOCO^{+} is 54.95 kJ/mol. This approach was also used to determine the ZPE for the BrHOCO^{+} reaction intermediate considered in Sec. III.

With the above values for the HBr^{+}, CO_{2}, and HOCO^{+} ZPEs, the 14.12 ± 0.47 kJ/mol (146 ± 5 meV) heat of reaction at 0 K becomes 4.24 ± 0.47 kJ/mol (44 ± 5 meV) without ZPEs. These are energies for the ground-state HBr^{+}(^{2}Π_{3/2}) + CO_{2} → Br(^{2}P_{3/2}) + HOCO^{+} pathway. Furthermore, we can assume that the splitting of the ^{2}Π multiplet of HBr^{+} and of the ^{2}P multiplet of Br is only due to the SO interaction of the states belonging to the same multiplets, because other states are well separated in energy. Then, the SO contribution to the ground state energy of HBr^{+} is half the splitting, i.e., 164.5 meV,^{12,21} and in the case of Br, it is one third of the splitting, i.e., 152.3 meV.^{47} So, the spin-orbit contribution to the reaction energy is very small, about 12 meV. By subtracting this contribution, we find that the reaction energy without ZPE and without SO is 32 ± 5 meV.

### C. Calculations with spin-orbit coupling

The spin-free states for the HBr^{+} + CO_{2} reaction are coupled by SO interactions and mix, producing a new set of states identified as ^{2}Π_{3/2} and ^{2}Π_{1/2}. SO coupling calculations were carried out with the Breit Pauli Hamiltonian, as implemented in the Molpro^{48} program using the state averaged complete active space self-consistent field (SA-CASSCF) theory,^{49} to determine PESs for these ^{2}Π_{3/2} and ^{2}Π_{1/2} states. Dunning’s cc-pVTZ basis set^{50} was used (fully uncontracted and up to p functions for H, d for C and O, and f for Br). In order to run meaningful SO coupling calculations, three degenerate 4p^{5} states of the Br atom need to be included, which correspond to the ground ^{2}Π state plus a ^{2}Σ excited state of HBr^{+}. However, if equal weights were applied to the three states to perform a SA-CASSCF calculation, the ^{2}Σ state would be higher in energy than a charge transfer state, HBr + CO_{2}^{+} (^{2}Π). To remove this latter state, we used the “dynamical weighting” ansatz, as implemented in Molpro; i.e., the weight of each state is computed as

where Δ*E _{i}* is the energy difference between state

*i*and the ground state. Setting the constant

*α*to 9 a.u. enables the program to shift the charge transfer

^{2}Π state to a higher energy than the

^{2}Σ state, so that only 3 states had to be taken into account in the SA-CASSCF calculations. The active space includes 5 orbitals and 9 electrons, which are one

*σ*molecular orbital (MO) and a pair of

*π*non-bonding MO’s on HBr, and a pair of

*π*bonding MO’s on CO

_{2}for the reactants; and the three 4

*p*orbitals of Br, and one a′ MO and one a″ MO on HOCO

^{+}for the products.

## III. COMPUTATIONAL RESULTS

### A. PES without spin-orbit coupling

#### 1. CCSD(T) PES

As discussed in the Introduction, a goal of this study is to develop accurate QM+MM PESs^{18} for the HBr^{+} + CO_{2} → Br + HOCO^{+} reaction with HBr^{+} in the ^{2}Π_{1/2} and ^{2}Π_{3/2} spin-orbit states. In this model, QM represents the average ^{2}Π PES without spin-orbit coupling and MM are analytic potential energy functions fit to spin-orbit coupling calculations (see below). An accurate QM model is required for this representation of the PES. As for other QM+MM models,^{18,51–53} the terms for each PES are additive with one QM and the other MM.

In previous work,^{7} Paetow *et al.* calculated QM energies for the HBr^{+} + CO_{2} → Br + HOCO^{+} reaction with UMP2 and CCSD(T) theories and double- and triple-zeta basis sets. Values for the reaction energy, ΔE_{r}, and the energy of the reaction intermediate BrHOCO^{+}, E(BrHOCO^{+}), were reported without ZPEs included and with a harmonic ZPE correction. As shown in Table I, the calculated ΔE_{r}^{QM}, without a ZPE correction range from approximately −114 to −135 meV.

##### a. *HBr*^{+} + *CO*_{2} → *Br* + *HOCO*^{+} reaction energy.

With the objective to establish an accurate QM energy for the HBr^{+} + CO_{2} → Br + HOCO^{+} reaction, a set of CCSD(T) calculations were performed to determine the reaction energy ΔE_{r}. Two pseudo potentials were considered for Br, i.e., PP by Peterson *et al.*^{54} and SDB by Martin and Sundermann.^{55} In addition, a range of different correlation consistent cc-pVXZ and aug-cc-pVXZ basis sets^{50} were used, with X = D, T, and Q to represent double-, triple-, and quadruple zeta basis sets, respectively.

The results of these CCSD(T) calculations are listed in Table I. A similar CBS limit is found for the “cc-” and “aug-cc-” basis sets with the PP. For the former basis sets, the CBS limit for the reaction energy ΔE_{r}^{QM} is 75 meV, while 62 meV for the latter. For the calculations with the SDB pseudopotential, the aug-cc-pVTZ and aug-cc-pVQZ basis sets give ΔE_{r}^{QM} of 86 and 77 meV, respectively. A comparison of these results, with consideration of effects of the basis set on ΔE_{r} and also the CBS limit, suggests that the CCSD(T)/CBS/cc-pVXZ/PP value for ΔE_{r}^{QM} of 75 meV is the most accurate. For these calculations, the cc-pVTZ/PP and cc-pVQZ/PP basis sets give very similar ΔE_{r}^{QM} values, which are nearly identical to the CBS value. A value of 75 meV is used as the QM benchmark for the reaction energy. This value is only slightly higher than the above “experimental value” of 32 ± 5 meV, obtained by subtracting the ZPE and SO contributions.

Though complete spin-orbit calculations are given below, it is possible to obtain a meaningful approximate ΔE_{r} value for HBr^{+}(^{2}Π_{3/2}) + CO_{2} → Br (^{2}P_{3/2}) + HOCO^{+}, *in lieu* of these calculations, based on ΔE_{r}^{QM} = 75 meV. The splitting of the ^{2}P_{1/2} and ^{2}P_{3/2} state for the Br-atom is 456.9 meV. Since the ^{2}P_{3/2} state is four-fold degenerate and the ^{2}P_{1/2} state is two-fold degenerate, the QM energy for the Br-atom lies 152.3 meV above the ^{2}P_{3/2} energy.^{47} As discussed below as part of the spin-orbit calculations, to a quite good approximation, the QM energy for HBr^{+} lies midway between the energies of the ^{2}Π_{1/2} and ^{2}Π_{3/2} states, each doubly degenerate. The splitting of these two states is 329 meV, placing the QM energy for HBr^{+} ∼ 164.5 meV above the ^{2}Π_{3/2} state.^{12,21} With these QM energies for the Br atom and HBr^{+}, and the above value for ΔE_{r}^{QM}, ΔE_{r} for the HBr^{+}(^{2}Π_{3/2}) + CO_{2} → Br (^{2}P_{3/2}) + HOCO^{+} reaction without ZPE is (75 − 152.3 + 164.5) meV = 87 meV. This value is only slightly higher than the above “experimental value” of 44 ± 5 meV without ZPE, further corrected to 32 ± 5 meV by subtracting the SO contribution.

##### b. Energies of the *BrHOCO*^{+} and *HBr*⋯*OCO*^{+} intermediates.

Electronic structure calculations were also performed to determine the QM energy, without ZPE, for the BrHOCO^{+} reaction intermediate and the results are given in Table I. The BrHOCO^{+} energy, obtained previously by Paetow *et al.*,^{6} from UMP2 and CCSD(T) calculations, varies from −706 to −824 meV, with respect to the separated reactants HBr^{+} + CO_{2}. For the current study, an optimized BrHOCO^{+} geometry and energy were obtained with CCSD(T) calculations employing the cc-pVDZ/PP and cc-pVTZ/PP basis sets, and their respective energies are −668 and −698 meV (see Table I). Because of the extreme memory requirement for CCSD(T)/cc-pVQZ/PP calculations, it was not possible to obtain an optimized BrHOCO^{+} structure and concomitant energy at this level of theory. To approximate the CBS limiting energy for BrHOCO^{+}, the above cc-pVDZ and cc-pVTZ energies for BrHOCO^{+} were extrapolated using the well-known expression^{56,57}

in which X is the cardinal number of the basis set. The resulting CBS limit energy for BrHOCO^{+} is −710 meV. As a test of this extrapolation method, the reaction energy (ΔE_{r}^{QM}) calculated using Eq. (5), and only the double and triple zeta basis sets, was compared with that obtained in Sec. II using Eq. (1) and the double, triple, and quadruple basis sets. ΔE_{r}^{QM} from Eq. (5) is 88 meV and in a quite good agreement with the value of 75 meV obtained from Eq. (1).

In addition to the BrHOCO^{+} reaction intermediate, there is an [HBr⋯OCO]^{+} van der Waals’ intermediate in the entrance channel of the PES, whose geometry is depicted in Figure 2. The CCSD(T) energy for [HBr⋯OCO]^{+}, without ZPE, was calculated with the cc-pVDZ and cc-pVTZ basis sets and the respective values are −534 and −539 meV. The CBS extrapolated energy, using Eq. (5), is −541 meV. Thus, the potential energy for the [HBr⋯OCO]^{+} van der Waals’intermediate is ∼170 meV higher in energy than that for the BrHOCO^{+} reaction intermediate.

#### 2. An accurate QM method for direct dynamics simulations

As discussed above, CCSD(T) theory, extrapolated to the CBS limit, gives a quite accurate ΔE_{r}^{QM} = 75 meV for the HBr^{+} + CO_{2} → Br + HOCO^{+} reaction in the absence of spin-orbit coupling. However, this level of theory is impractical for the projected QM+MM direct dynamics simulations and it is important to identify a computationally practical, but sufficiently accurate, QM method for the simulations. Different QM methods and basis sets were tested and the results are summarized in Table II. An ideal QM method would not only give the accurate reaction energy ΔE_{r}^{QM} but also correctly represent the energies for the intermediates BrHOCO^{+} and [HBr⋯OCO]^{+}. In the following, QM energy values are considered for ΔE_{r}^{QM}, BrHOCO^{+}, and [HBr⋯OCO]^{+} consecutively.

. | Method . | |||
---|---|---|---|---|

Basis set . | B3LYP . | B98 . | PBE0 . | UMP2 . |

ΔE_{r}^{QM} | ||||

3-21G | −477 | |||

6-31G**/lanl2dz | −460 | |||

6-31G**/lanl2dzdp | −154 | −144 | −156 | |

6-311G | −197 | |||

6-311G** | 6 | −5 | 13 | 5 |

6-311G**/SDB | −199 | |||

cc-pVDZ/PP | −31(64) | −37 | −11 | −30 |

cc-pVTZ/PP | 83 | |||

cc-pVTZ/SDB | 10(125) | −24 | −24 | −43 |

aug-cc-pVDZ/lanl2dzdp | 102(216) | 81 | 71 | 309 |

aug-cc-pVTZ/SDB | 104(218) | 79 | 75 | 83 |

E^{QM}(BrHOCO^{+}) | ||||

3-21G | −635 | |||

6-31G**/lanl2dz | −1111 | |||

6-31G**/lanl2dzdp | −898 | −898 | −931 | |

6-311G | ||||

6-311G** | −792 | −816 | −831 | |

6-311G**/SDB | −865 | |||

cc-pVDZ/PP | −868 | −875 | −878 | −739 |

cc-pVTZ/PP | −731 | |||

cc-pVTZ/SDB | −809(−753) | −847 | −878 | −811 |

aug-cc-pVDZ/lanl2dzdp | −731(−692) | −762 | −798 | |

aug-cc-pVTZ/SDB | −741(−699) | −781 | −815 | −779 |

. | Method . | |||
---|---|---|---|---|

Basis set . | B3LYP . | B98 . | PBE0 . | UMP2 . |

ΔE_{r}^{QM} | ||||

3-21G | −477 | |||

6-31G**/lanl2dz | −460 | |||

6-31G**/lanl2dzdp | −154 | −144 | −156 | |

6-311G | −197 | |||

6-311G** | 6 | −5 | 13 | 5 |

6-311G**/SDB | −199 | |||

cc-pVDZ/PP | −31(64) | −37 | −11 | −30 |

cc-pVTZ/PP | 83 | |||

cc-pVTZ/SDB | 10(125) | −24 | −24 | −43 |

aug-cc-pVDZ/lanl2dzdp | 102(216) | 81 | 71 | 309 |

aug-cc-pVTZ/SDB | 104(218) | 79 | 75 | 83 |

E^{QM}(BrHOCO^{+}) | ||||

3-21G | −635 | |||

6-31G**/lanl2dz | −1111 | |||

6-31G**/lanl2dzdp | −898 | −898 | −931 | |

6-311G | ||||

6-311G** | −792 | −816 | −831 | |

6-311G**/SDB | −865 | |||

cc-pVDZ/PP | −868 | −875 | −878 | −739 |

cc-pVTZ/PP | −731 | |||

cc-pVTZ/SDB | −809(−753) | −847 | −878 | −811 |

aug-cc-pVDZ/lanl2dzdp | −731(−692) | −762 | −798 | |

aug-cc-pVTZ/SDB | −741(−699) | −781 | −815 | −779 |

^{a}

Energies are in meV and are with respect to zero of energy for HBr^{+} + CO_{2}. Zero-point energies are not included in E^{QM}(BrHOCO^{+}), the energy of the BrHOCO^{+} intermediate with respect to reactants, and ΔE_{r}^{QM}, the reaction energy. Harmonic zero-point energy corrections are included for the energies in parenthesis.

^{b}

As benchmark results to compare with, ΔE_{r}^{QM} from CCSD(T)/CBS/cc-pVXZ (X = D, T, Q) and CCSD(T)/CBS/aug-cc-pVXZ (X = D, T, Q) are 75 and 62 meV, respectively. E^{QM}(BrHOCO^{+}) from CCSD(T)/CBS/cc-pVXZ (X = D, T) is −710 meV.

Table II shows that Pople-type basis sets give negative or small positive values of ΔE_{r}^{QM} compared to the CCSD(T)/CBS value of 75 meV. This result was illustrated previously by the calculations of Paetow *et al.*^{6} To illustrate the results in Table II, UMP2/6-311G, UMP2/6-311G*, and PBE0/6-311G* give ΔE_{r}^{QM} values of −197, 5, and 13 meV, respectively. Accurate values for ΔE_{r}^{QM} are obtained with the Dunning correlation consistent (cc) basis sets.^{49} B98 and PBE0, with the aug-cc-pVDZ/lanl2dzdp basis, give respective ΔE_{r}^{QM} values of 81 and 71 meV. UMP2 gives ΔE_{r}^{QM} = 83 meV with the aug-cc-pVTZ/SDB basis set. The B3LYP ΔE_{r}^{QM} values with the “cc” basis sets are somewhat larger than the CCSD(T)/CBS value of 75 meV, i.e., with the aug-cc-pVDZ/lanl2dzdp and aug-cc-pVDZ/SDB basis sets, the values are 102 and 104 meV, respectively.

The results in Table II show that the calculated energy for the BrHOCO^{+} reaction intermediate strongly depends on the method, basis set, and pseudopotential. DFT and MP2 using Pople-type basis sets and the lanl2dz pseudopotential substantially overestimate the HBr^{+} + CO_{2} → BrHOCO^{+} binding energy as compared to the CCSD(T)/CBS(D+T) value of −710 meV. The performance of DFT improves as the basis set gets larger, e.g., the binding energy is ∼ − 790 meV for B3LYP/6-311G**. Of the different DFT functionals, B3LYP with the aug-cc-pVTZ/SDB and aug-cc-pVDZ/land2dzdp basis sets gives the most accurate QM energies for BrHOCO^{+}, which are −741 and −731 meV, respectively. The MP2 results are best with the PP and UMP2 with cc-pVDZ/PP and cc-pVTZ/PP are ∼ − 740 and −730 meV, respectively, as compared to the CCSD(T)/CBS(D+T) value of −710 meV. Considering the performance of different electronic structure theory methods in reproducing the “accurate” ΔE_{r}^{QM} and the QM energy for BrHOCO^{+}, B3LYP/aug-cc-pVDZ/lanl2dzdp and UMP2/cc-pVTZ/PP are the best QM methods for the direct dynamics simulations.

It is important to verify that B3LYP/aug-cc-pVDZ/lanl2dzdp and UMP2/cc-pVTZ/PP are able to represent the van der Waals intermediate complex [HBr⋯OCO]^{+}. The relative energy of [HBr⋯OCO]^{+} from B3LYP/aug-cc-pVDZ/lanl2dzdp and UMP2/cc-pVTZ/PP are −705 and −532 meV, respectively, as compared to the CCDS(T)/CBS(D + T) value of −541 meV. B3LYP/aug-cc-pVDZ substantially overestimates this van der Waals’ interaction and only UMP2/cc-pVTZ/PP yields accurate energies for ΔE_{r}^{QM}, BrHOCO^{+}, and [HBr⋯OCO]^{+}.

### B. Spin-orbit coupling

The ^{2}Π_{3/2} and ^{2}Π_{1/2} spin-orbit splitting has been measured for HBr^{+} and is 329 meV.^{12,21} For the Br atom, the measured ^{2}P_{3/2} and ^{2}P_{1/2} spin-orbit splitting is 456.9 meV.^{47} The corresponding calculated splittings, as described in Sec. II C, are 295 meV for HBr^{+} and 411.7 meV for Br, both about 10% smaller than the experimental values. Hence, the calculated spin orbit contribution to the reaction energy is 10 meV, only 2 meV smaller than the experimental one (see Sec. II B).

The SO calculations were divided into two sets. One for the entrance-channel region of the PES from the reactants HBr^{+} + CO_{2} to the reaction intermediate BrHOCO^{+}, and the other for the exit-channel region from the intermediate to the products HOCO^{+} + Br. Representative potential energy scans for the entrance- and exit-channel regions of the PES are shown in Figures 3 and 4, respectively. Additional scans for the entrance-channel region are given in the supplementary material.^{58} Included in each scan is the spin-free potential energy curve and the potential energy curves for the HBr^{+} ^{2}Π_{3/2} and ^{2}Π_{1/2} spin-orbit states. The properties of these curves are discussed in the following.

#### 1. SO coupling for the ^{2}*Π*_{3/2} *HBr*^{+} PES

As shown in Figures 3, 4, and S1 in the supplementary material,^{58} the spin-free potential energy curves are nearly identical to those for HBr^{+} in the ^{2}Π_{3/2} state, except the former are higher in energy. For the asymptotic HBr^{+} + CO_{2} reactants, the spin-free curve is 3.49 kcal/mol higher in energy, while for the Br + HOCO^{+} asymptotic products, this difference is 3.24 kcal/mol. Thus, using the spin-free PES to represent HBr^{+} in the ^{2}Π_{3/2} state only introduces a 0.25 kcal/mol error in the relative energies of the asymptotic reactants and products.

If the spin-free and ^{2}Π_{3/2} potential energy curves are shifted so that their potential energies are identical at the maximum distances in the plots, e.g., 20 Å in scan (a), the resulting potential energy curves are nearly identical for the scans, as shown in Figures S2 and S3 of the supplementary material.^{58} The only significant differences are for two of the entrance-channel scans in Figure 3 (the shifted curves are in Figure S2 of the supplementary material^{58}). For the shortest distance in scan (b), the spin-free potential energy is 4.11 kcal/mol higher than that for the ^{2}Π_{3/2} state; and for the repulsive region in scan (f), the spin-free curve is higher in energy. The only significant differences in the spin-free and ^{2}Π_{3/2} exit-channel scans in Figure 4 (the shifted curves are in Figure S3 of the supplementary material^{58}) are for the shortest distances in the scans (a) and (b), where the spin-free potential energies are 1.37 and 1.88 kcal/mol higher than those for the ^{2}Π_{3/2} state, respectively.

The only significant differences between the spin-free and ^{2}Π_{3/2} potential energy curves are for short-range repulsive interactions. Thus, the spin-free PES is expected to be a quite good model for both the entrance- and exit-channel regions of the PES for HBr^{+} in the ^{2}Π_{3/2} state and may be used in direct dynamics simulations for this state.

#### 2. SO coupling for the ^{2}*Π*_{1/2} *HBr*^{+} PES

Comparisons of scans for the spin-free potential energy and those for the ^{2}Π_{1/2} state are shown in Figures 3, 4, and S1 in the supplementary material.^{58} For the entrance-channel scans in Figure 3, there is a very good agreement between the spin-free and ^{2}Π_{1/2} state potential energy curves. For the asymptotic HBr^{+} + CO_{2} reactants, the ^{2}Π_{1/2} state is 3.32 kcal/mol higher than the spin-free curve. If the entrance-channel spin-free and ^{2}Π_{1/2} potential energy curves are shifted so that their potential energies are identical at the maximum distances in the plots, the resulting potential energy curves for the scans are nearly identical as shown in Figure S4 in the supplementary material.^{58} The only significant difference is for the shortest distance in Figure 3(b) (the shifted curves are in Figure S4 of the supplementary material^{58}) where the ^{2}Π_{1/2} energy is 3.97 kcal/mol higher in energy than the spin-free value. Thus, the spin-free PES is a very good representation of that for the ^{2}Π_{1/2} state in the HBr^{+} + CO_{2} entrance channel, with the latter lowered by 3.32 kcal/mol.

Figure 4 shows there are important differences between the spin-free and ^{2}Π_{1/2} potential energy curves in the exit-channel region of the PES. At the Br + HOCO^{+} product asymptotic limit, the ^{2}Π_{1/2} energy is 3.15 kcal/mol lower than the spin-free energy. Figure 5 compares the spin-free and ^{2}Π_{1/2} potential energy curves for the exit-channel scans, with the potential energy curves shifted so that they are identical at the scans’ maximum Br + HOCO^{+} separations. Important differences are clearly evident in the scans. In the following, an analytic function is developed to represent the effect of spin-orbit coupling for the ^{2}Π_{1/2} HBr^{+} + CO_{2} → Br + HOCO^{+} reaction in the product exit-channel.

#### 3. Analytic representation of the spin-orbit coupling for the ^{2}*Π*_{1/2} PES

A potential energy surface for the ^{2}Π_{1/2} state was developed by combining the spin-free potential energy surface with an analytic, i.e., MM-like, representation of the differences between the spin-free curves and those for the ^{2}Π_{1/2} state. The function for the ^{2}Π_{1/2} PES is written as

where *V _{R}* is the PES for the reactants, entrance-channel region,

*V*is the PES for the products, exit-channel region, and

_{P}*S*(

*q*) is a switching function which connects these two regions.

*V*is expressed as the $ V spin - free + V R shift $, where the former is the spin-free PES and the latter is the difference in the spin-free and the

_{R}^{2}Π

_{1/2}potential energies when the reactants are at their asymptotic separation.

*V*is expressed as

_{P} The first two terms in this equation are analogous to those for *V _{R}*, while the last is a fit to the difference between the spin-free and the

^{2}Π

_{1/2}potential energy curves, with the two sets of curves shifted so that they match at the product asymptotic separation.

The differences between the shifted ^{2}Π_{1/2} and the spin-free potential energy curves, for the product, exit-channel region of the PES, are plotted in Figure 5. An accurate fit to the differences in the curves was obtained by a sum of two-body terms between the Br-atom and the H- and O-atoms of HOCO^{+}. The two-body terms are written as

where *r* is the Br-H distance or one of the Br-O distances. Thus, *V _{fit}* was represented by a sum of the three Br-H and Br-O two-body terms. The excellent fits obtained are illustrated in Figure 6, and the fitted

*A*,

*B*,

*C*,

*D*,

*n*, and

*m*parameters for the Br-H and Br-O terms are listed in Table III.

Interaction . | A
. | B
. | C
. | D
. | n
. | m
. |
---|---|---|---|---|---|---|

Br-H | 297.464 | 2.793 98 | −9.565 44 | −13.195 6 | 8 | 2 |

Br-O_{1}^{b} | 722.662 | 1.672 22 | −552.752 | −145.305 | 5 | 3 |

Br-O_{2}^{c} | 3298.71 | 18.095 0 | −1515.71 | −87 367.9 | 8 | 13 |

Interaction . | A
. | B
. | C
. | D
. | n
. | m
. |
---|---|---|---|---|---|---|

Br-H | 297.464 | 2.793 98 | −9.565 44 | −13.195 6 | 8 | 2 |

Br-O_{1}^{b} | 722.662 | 1.672 22 | −552.752 | −145.305 | 5 | 3 |

Br-O_{2}^{c} | 3298.71 | 18.095 0 | −1515.71 | −87 367.9 | 8 | 13 |

^{a}

Equation (7) defines the potential energy function for the ^{2}Π_{1/2} PES. The MM fit is to Eq. (8) where units of the parameters are *A* in kcal/mol, *B* in Å^{−1}, *C* in kcal Å^{n}/mol, and *D* in kcal Å^{m}/mol.

^{b}

O_{1} is the oxygen atom of HOCO^{+} which is attached to the H and C atoms.

^{c}

O_{2} is the Oxygen atom of HOCO^{+} which is attached only to the C atom.

The remaining component needed for the potential energy surface, *V*(^{2}Π_{1/2}), is the switching function *S*(*q*) connecting the reactant, entrance-channel and product, exit-channel regions of the potential. *S*(*q*), *q* = *r*_{H−Br} − *r*_{H−O}, is given by

where *q _{o}*,

*a*, and

*n*are parameters, determined by fitting Eq. (6) to all the points in the potential energy scans given in Figures 3 and 4 and the supplementary material

^{58}for the

^{2}Π

_{1/2}HBr

^{+}+ CO

_{2}reaction. The resulting values for the parameters are

*q*= 0.6215 Å,

_{o}*a*= 1263.0 Å

^{−2}, and

*n*= 2. Shown in Figure 7 are illustrations of the excellent fit by this switching function to similar potential energy points for the reactant

*V*and product

_{R}*V*regions of the PES.

_{P}## IV. SUMMARY

As described in the following, extensive calculations and analyses were made to develop accurate QM+MM PES models for the HBr^{+} + CO_{2} → Br + HOCO^{+} reaction with HBr^{+} in the ^{2}Π_{3/2} and ^{2}Π_{1/2} spin-orbit states:

Accurate anharmonic ZPE corrections were made to obtain an accurate “experimental” energy without ZPE for the ground state HBr

^{+}(^{2}Π_{3/2}) + CO_{2}→ HOCO^{+}+ Br(^{2}P_{3/2}) reaction, to compare with the results of electronic structure calculations.CCSD(T) electronic structure calculations were performed to determine “benchmark” spin-free QM energies for the HBr

^{+}+ CO_{2}→ HOCO^{+}+ Br reaction.With zero-point energies removed, the “experimental” reaction energy is 44 ± 5 meV for HBr

^{+}(^{2}Π_{3/2}) + CO_{2}→ Br(^{2}P_{3/2}) + HOCO^{+}, while the CCSD(T) value with spin-orbit effects included is 87 meV.Electronic structure calculations were performed to determine structures, vibrational frequencies, and energies for the intermediates BrHOCO

^{+}and [HBr⋯OCO]^{+}.To determine a spin-free QM model for direct dynamics simulations, calculations were performed with a broad range of electronic structure methods and their results were compared with those obtained with CCSD(T). Only UMP2/cc-pVTZ/PP was found to be a practical and accurate QM method to use in QM+MM direct dynamics simulations.

The spin-free states are coupled by their SO interactions and mix, producing a new set of states

^{2}Π_{3/2}and^{2}Π_{1/2}. The SO coupling calculations were performed to determine PESs for these states.The SO coupling calculations show that the spin-free QM PES gives a quite good representation of the PES for

^{2}Π_{3/2}HBr^{+}.The spin-free QM PES accurately describes the reactant, entrance-channel region of the PES for

^{2}Π_{1/2}HBr^{+}reaction. However, spin-orbit coupling effects are important for the product, exit-channel region of this PES. A MM model was developed to represent these effects, which were combined with the spin-free QM PES to form a QM+MM model of the PES for the^{2}Π_{1/2}HBr^{+}reaction.In principle, the PESs for these

^{2}Π_{3/2}and^{2}Π_{1/2}HBr^{+}states are non-adiabatically coupled. However, there are no crossings of these PESs. The PESs are almost parallel in their reactant and intermediate regions and gradually approach in forming the products. Therefore, non-adiabatic transitions between the^{2}Π_{3/2}and^{2}Π_{1/2}states are expected to be unimportant for the HBr^{+}+ CO_{2}→ Br + HOCO^{+}reaction dynamics.

In future work, the PESs determined and developed here will be used in direct dynamics simulations of the HBr^{+} + CO_{2} → Br + HOCO^{+} reaction with HBr^{+} in the ^{2}Π_{3/2} and ^{2}Π_{1/2} spin-orbit states.

## Acknowledgments

The calculations reported here are also based upon work supported by the National Science Foundation under the Partnership in International Research and Education (PIRE) Grant No. OISE-0730114, and the Robert A. Welch Foundation under Grant No. D-0005. Support was also provided by the High-Performance Computing Center (HPCC) at Texas Tech University, under the direction of Philip W. Smith. The authors thank Karl-Michael Weitzel and Lisa Paetow for very helpful suggestions.

## REFERENCES

*NIST Atomic Spectra Database*, version 5.2, National Institute of Standards and Technology, Gaithersburg, MD, 2012; available online at http://physics.nist.gov/asd.