Simultaneous Evaluation of Multiple Rotationally ... - ACS Publications

Article pubs.acs.org/JPCA

Simultaneous Evaluation of Multiple Rotationally Excited States of H3+, H3O+, and CH5+ Using Diffusion Monte Carlo Andrew S. Petit, Jason E. Ford, and Anne B. McCoy* Department of Chemistry and Biochemistry, The Ohio State University, Columbus, Ohio 43210, United States S Supporting Information *

ABSTRACT: An extension to diffusion Monte Carlo (DMC) is proposed for simultaneous evaluation of multiple rotationally excited states of fluxional molecules. The method employs an expansion of the rotational dependence of the wave function in terms of the eigenstates of the symmetric top Hamiltonian. Within this DMC approach, each walker has a separate rotational state vector for each rotational state of interest. The values of the coefficients in the expansion of the rotational state vector associated with each walker, as well as the locations of the walkers, evolve in imaginary time under the action of a propagator based on the imaginary-time time-dependent Schrödinger equation. The approach is first applied to H3+, H2D+, and H3O+ for which the calculated energies can be compared to benchmark values. For low to moderate values of J the DMC results are found to be accurate to within the evaluated statistical uncertainty. The rotational dependence of the vibrational part of the wave function is also investigated, and significant rotation−vibration interaction is observed. Based on the successful application of this approach to H3+, H2D+, and H3O+, the method was applied to calculations of the rotational energies and wave functions for CH5+ with v = 0 and J ≤ 10. Based on these calculations, the rotational energy progression is shown to be consistent with that for a nearly spherical top molecule, and little evidence of rotation−vibration interaction is found in the vibrational wave function.

■

INTRODUCTION Gas-phase spectroscopy continues to provide a powerful experimental tool for addressing fundamental questions in chemical dynamics, atmospheric chemistry, combustion chemistry, and astrochemistry. Analysis of vibrationally resolved spectra has provided important insights into the effects of solvation on chemical reactions, the nature of hydrogen bonding, and the properties of highly reactive intermediates.1−4 Analysis of high-resolution, rotationally resolved spectra allows for the accurate determination of molecular structure as well as investigations of the strength and nature of rotation−vibration coupling and the pathways and dynamics of intramolecular vibrational energy redistribution.5,6 Moreover, these techniques have been used to study the properties and kinetics of weakly bound and/or highly reactive species that are difficult to study using other approaches.7−9 Within the context of astrochemistry, rotationally resolved spectroscopy provides the primary tool for studying the chemistry present in the interstellar medium, and the experimental determination of the spectral fingerprints of potential astrochemical species represents an ongoing effort in the field. The theoretical calculation of ro-vibrational states can aid in the interpretation and assignment of experimental spectra as well as the investigation of the unique spectral signatures and relevant spectral ranges of a given molecular species. For semirigid molecules, the harmonic oscillator/rigid-rotor model provides a reasonable zeroth-order description. Perturbation theory can then be used to develop model Hamiltonians that © XXXX American Chemical Society

contain higher-order terms, which account for the effects of anharmonicity, centrifugal distortion, and Coriolis coupling. Although the coefficients in the model Hamiltonian often cannot be expressed as simple expressions of molecular properties, systematic approaches for their evaluation have been developed.6,10−14 When large amplitude, zero-point vibrational motions are present, however, the situation becomes considerably more complicated as the harmonic oscillator/rigid-rotor model no longer provides a reasonable zeroth-order picture of the nuclear dynamics. This is particularly true when the wave function of the ground vibrational state samples multiple minima on the potential energy surface (PES). As a result of this, the perturbation theory expansion becomes slow to converge or even divergent for fluxional molecules and clusters.15,16 Variational calculations provide a powerful and systematic means of overcoming the limitations of perturbation theory.17−19 By solving the Schrödinger equation on a grid of points in the physically relevant regions of configuration space, the delocalization of the ground vibrational state wave function over multiple minima on the PES can be accurately treated. Indeed, for fluxional molecules with three or four atoms, converged, full-dimensional variational calculations that achieve Special Issue: Kenneth D. Jordan Festschrift Received: September 3, 2013 Revised: September 19, 2013

A

dx.doi.org/10.1021/jp408821a | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

spectroscopic accuracy have been reported.20−23 However, the computational cost of this approach scales very unfavorably with system size, a problem that is further exacerbated by the need for a large basis to accurately describe the degrees of freedom which undergo large amplitude vibrational motions. Extending these variational approaches to states with J > 0 is made more challenging by the significant increase in the density of states and the inclusion of additional terms in the kinetic energy operator.18 Wang and Carrington’s recent full-dimensional variational calculation of the J = 0 states of CH5+ represents a heroic example of such a calculation. Due to the computational expense involved, extending their analysis to states with J greater than one or two will be extremely challenging.24 In the present study, we use diffusion Monte Carlo (DMC) to investigate the nature of the rotation−vibration coupling present in the rotationally excited states of fluxional molecular systems. In this method, the time-dependent Schrödinger equation is propagated in imaginary time using statistical techniques. As a result DMC has a much more favorable scaling with system size than variational approaches.25−28 Because of this, DMC has been applied to a number of especially challenging problems in electronic structure theory.28−30 Moreover, a number of previous studies have demonstrated DMC to be a powerful technique for studying the rotational and vibrational ground and excited state properties of fluxional molecules and clusters as well as the rotationally excited states of molecules embedded in quantum clusters.31−44 A significant challenge comes in the application of DMC to the study of excited states as DMC is fundamentally a ground state method. The correlation function approach and the fixednode approximation represent two approaches for circumventing this limitation. In the former, trial functions of the appropriate symmetry are constructed and the imaginary-time time-dependent autocorrelation functions of these trial functions are evaluated using DMC.45,46 Such calculations allow for the accurate determination of a number of rovibrational states with the same symmetry but do not provide any of the excited state wave functions. In contrast, the fixednode approximation allows the DMC simulation to converge onto the lowest energy state consistent with a predetermined set of nodal surfaces obtained from model calculations or symmetry considerations.25,26 This predetermination of the nodal surfaces introduces a more severe approximation than is made in the correlation function approach but allows for the evaluation of the wave functions for these excited states. Analysis of these wave functions can often provide important physical insights into how the excitation affects the molecular structure and nuclear dynamics of the system. In a number of previous calculations we have demonstrated that fixed-node DMC can be effectively used to study the vibrationally and rotationally excited states of highly fluxional molecules and clusters.35−41,47−51 Despite their success, these fixed-node DMC studies were very computationally demanding, with a separate set of statistically independent DMC simulations required for each nodal region of each excited state considered in the studies. In this study, we present an approach for describing multiple rotationally excited states within a single DMC calculation. In doing so, we dress each member of the DMC ensemble with a set of rotational state vectors, each of which evolves under the action of a propagator based on the rotational Hamiltonian. Since the inverse moment of inertia tensor that is present in the

rotational Hamiltonian depends on the molecular geometry, the effects of rotation−vibration coupling are captured in the DMC simulations. By providing a separate set of weights and reference energy for each rotational state that is included in the calculation, we are able to perform correlated sampling of the rotationally excited states. Further the computer time required for these calculations scales sublinearly with the number of rotational states included in the calculation. While the present study presents the first application of such an approach to rotationally excited states of fluxional molecules, similar correlated sampling DMC approaches have been used to calculate tunneling splittings, the effects of rare gas solvation on the vibrational frequencies of the dopant, and the gradient of the electronic energy.52−55 To verify the validity of our multistate DMC method, we calculated the energies and wave functions for rotationally excited states of H3+, H2D+, and H3O+ with J as large as 10. Because of their large amplitude, zero-point vibrational motions and unusually large rotational constants, H3+ and H2D+ exhibit especially strong rotation−vibration coupling. Moreover, several full-dimensional PES’s of very high accuracy have been developed for H3+, and converged variational calculations to spectroscopic accuracy have been performed for the rotationally excited states with J ≤ 10 of both species.20,22,56−58 The third model system, H3O+, exhibits large amplitude, zeropoint vibrational motion along the umbrella inversion coordinate due to a low energy barrier in this mode of roughly 700 cm−1.59−61 Additionally, we have previously performed fixed-node DMC calculations on the rotationally excited states of these species which provide additional information with which to calibrate the multistate DMC method.48−50 As a final application, we apply the multistate DMC method to the rotationally excited states of CH5+, a system of longstanding interest due to its extreme fluxionality and speculated role as a key intermediate in astrochemical processes.62−66 Although three experimental spectra have been reported for CH5+, only the gross features in the lower resolution spectrum have been assigned. The high-resolution, rotationally resolved spectrum remains unassigned more than a decade after it was first reported.66−69 Several experimental groups are currently obtaining new, high-precision spectra of this astrochemically important molecule in an effort to simplify the assignment.70,71 The challenge in this assignment is due to the highly delocalized nature of the wave function of CH5+, with the ground state probability distribution having comparable amplitude at all 120 symmetry equivalent minima on the PES and at the 180 low-lying transition states that separate them.33−35,72,73 Bunker and co-workers have reported several studies on the rotationally excited states of CH5+ and its isotopologues, but the model Hamiltonians used in these studies appear to have underestimated the full extent of delocalization of the ground state wave function of CH5+.74,75 We have previously applied fixed-node DMC to the study of the J = 1 rotationally excited states of CH5+ but were unable to disentangle the rotational excitation from a low-lying tunneling excited state.51 In this study, we demonstrate that the multistate DMC approach is capable of describing the pure rotationally excited states of CH5+ and report the energies for all of the states with J ≤ 10.

■

THEORY Overview of Continuous Weighting DMC. Diffusion Monte Carlo (DMC) is a well-established method. In the B



Article

discussion that follows, we will provide a summary of the important aspects of DMC for the present study. More detailed accounts of the approach and its implementation can be found elsewhere.25,26,31,37 DMC provides a statistical approach for solving the timeindependent Schrödinger equation based on the isomorphism between the time-dependent Schrödinger equation in imaginary time, τ = it/ℏ, ∂ψ ℏ2 ∂ 2ψ = − Vψ ∂τ 2m ∂x 2

Wn(τ + δτ ) = exp[− (V (R⃗ n(τ + δτ )) − Eref (τ ))]Wn(τ ) (6)

This will reduce the weight of those walkers whose potential energy is higher than Eref(τ), treating these areas of configuration space as classically forbidden regions. Likewise, eq 6 increases the weight of walkers whose potential energy is lower than Eref(τ), building amplitude in classically allowed regions of configuration space. As a result, a large fraction of the total weight can quickly accumulate on the small number of walkers that occupy regions in configuration space that have low potential energies. This will lead to a poor Monte Carlo sampling of the wave function. To prevent this, the coordinates of a walker whose weight falls below a certain threshold, Wthresh, are replaced by the coordinates of the walker in the ensemble with the largest weight. The original weight is then distributed equally between the two walkers.31 After all low weight walkers are replaced, Eref(τ + δτ) is evaluated based on the location and weights of the walkers at time τ + δτ according to eq 5. In the limit of τ → ∞, an equilibrium is reached, and the distribution of walkers provides a Monte Carlo sampling of the ground state wave function while Eref(τ) fluctuates about the fully anharmonic zero-point energy of the system. The reported energies are obtained by averaging Eref(τ) over a large number of time steps. It should be noted that the accuracies of the ground state energy and wave function that are obtained by DMC are limited only by the accuracy of the potential energy surface that is used for the calculation, the size of the ensemble, the duration of the calculation, and the time step, δτ. Descendant weighting is used to obtain probability distributions, as well as expectation values.31,37 Introduction of Rotations to the DMC Simulations. In this study, the effect of rotation is introduced by dressing each of the walkers in the DMC ensemble with a rotational state vector of well-defined total rotational angular momentum, J, so that eq 3 becomes

(1)

and the diffusion equation after a first-order rate term has been added, ∂C ∂ 2C = D 2 − kC ∂t ∂x

(2)

In eq 1, the diffusion coefficient, D, is replaced by ℏ /2m, while the potential acts as a coordinate-dependent rate constant. As the diffusion equation can be solved by performing a statistical simulation, it is possible to use a similar approach to solve the Schrödinger equation in imaginary time. To build a simulation out of eq 1, we first expand the wave function in an ensemble of Nwalker δ-functions, 2

Nwalker

⟨R⃗|Ψ(τ )⟩ =

∑

Wn(τ )δ 3N (R⃗ − R⃗ n(τ ))

n=1

(3)

These δ-functions are referred to as walkers, and each walker is localized at a specific molecular geometry in configuration space, R⃗ n(τ). In the continuous weighting scheme, the number of walkers is fixed, and each walker has an associated weight, Wn(τ). The weight represents the relative importance of the walker in the overall ensemble.31 The DMC ensemble evolves according to |Ψ(τ + δτ )⟩ ≈ exp[− (V̂ − Eref (τ ))δτ ] exp(−T̂δτ ) |Ψ(τ )⟩ (4)

Nwalker

where the split operator approximation to the quantum propagator is made and a reference energy, Eref(τ), has been added,76 where Eref (τ ) = V̅ (τ ) − α

Wtot(τ ) − Wtot(0) Wtot(0)

⟨R⃗|Ψ(τ )⟩ =

∑

Wn(τ ) δ 3N (R⃗ − R⃗ n(τ ))|Φn , J (τ )⟩

n=1

(7)

The rotational state vectors are each expanded in a symmetric top basis,

(5)

Equation 5 has two contributions. The first is V̅ (τ), which represents the average potential energy of the DMC ensemble at imaginary time τ. The second term adjusts the reference energy when the sum of the weights Wtot(τ) deviates from its value at τ = 0, while α is a simulation parameter that controls the size of the fluctuations in Eref(τ). As such, Eref(τ) is evaluated in a way that ensures that Wtot(τ) remains roughly constant. Once the DMC ensemble reaches equilibrium, Eref will fluctuate about the zero-point energy of the system, while the distribution of walkers will provide a Monte Carlo sampling of the ground state wave function. In the first step of implementing the simulation, the walkers diffuse according to the kinetic energy part of the propagator in eq 4. More specifically, the kth atom of each walker undergoes a displacement in each of its Cartesian coordinates by a random amount based on a Gaussian distribution of width (δτ/mk)1/2. The potential energy of each walker is then evaluated at its new geometry, and the weight of the nth walker is updated according to

J

|Φn , J (τ )⟩

∑ K =−J

Cn , J , K (τ )|J , K ⟩ (8)

with {Cn,J,K(τ)} representing the set of imaginary time dependent expansion coefficients for the nth walker. In the case of a spherical or symmetric top molecule, for example, H3+, H3O+, or CH5+, these basis functions are eigenfunctions of the rigid rotor Hamiltonian. As such, we use symmeterized linear combinations of the symmetric top eigenstates as the τ = 0 rotational state vectors for these molecular ions. The |J,K⟩ states also provide a basis in which the asymmetric top rotational wave functions can be expanded. In most of this study, K will provide the projection of J onto the body-fixed C-axis, so, unless otherwise noted K = KC. The |J,KA⟩ basis is used in some cases to investigate the sensitivity of the results to the choice of rotational basis. To account for the presence of the rotational state vectors in the propagation of the DMC ensemble, eq 4 is generalized to C



Article

position of each of the walkers, the overall distribution of walkers, and the corresponding energy are independent of the value of the Euler angles. We now consider the procedure in which the DMC ensemble given in eq 7 evolves from τ to τ + δτ under the action of the imaginary time quantum propagator of eq 9, focusing for clarity of discussion on the nth walker. As in ground state DMC, the kinetic energy component of the propagator produces a random displacement of the walker in configuration space from R⃗ n(τ) to R⃗ n(τ + δτ). The local rotational Hamiltonian is then constructed at R⃗ n(τ + δτ) and ̂ δτ) allowed to act on the rotational state vector exp(−hrot carried by the walker according to first-order time dependent perturbation theory, so that

̂ δτ ) × |Ψ(τ + δτ )⟩ ≈ exp[− (V̂ − Eref (τ ))δτ ] exp(−hrot exp( −T̂δτ )|Ψ(τ )⟩

(9)

where the additional exp(−ĥrotδτ) to the propagator causes the evolution of the rotational state vectors. More specifically, ̂ (R ⃗ ) = hrot

1 ∑ J ̂ Iα−,1β(R⃗)Jβ̂ 2 α ,β α

(10)

is the local rotational Hamiltonian defined using the inverse moment of inertia tensor, I−1(R⃗ ), evaluated in the Eckart frame at the molecular geometry R⃗ .77−79 The operators Jα̂ and Jβ̂ , which represent the αth and βth components of the rotational angular momentum operator in the Eckart frame, are defined to only act on the rotational state vectors in eq 7. In this study, we approximate Ĥ rot−vib by the sum of Ĥ J=0 and ĥrot, where Ĥ j=0 represents the J = 0 Hamiltonian used in the ̂ is ground state DMC calculations described above, and hrot defined in eq 10. As a result, Coriolis coupling terms are not included in the Hamiltonian. The severity of this approximation is minimized by embedding the body-fixed axis system using the Eckart conditions, which are constructed to eliminate the Coriolis coupling terms in the limit of small-amplitude vibrations.77,78 Although Coriolis coupling cannot be completely eliminated in systems that undergo large-amplitude vibrational motions, in several recent DMC studies of fluxional molecules and clusters, including CH5+, we and others have shown Eckart embedding to be a reasonable approach for obtaining vibrationally averaged rotational constants. Interestingly, these rotational constants were shown to be independent of the details of how the Eckart frame was defined, and in the case of symmetric top molecules two of the vibrationally averaged rotational constants were equal to within their statistical uncertainties.32,36,39,41,79 Eckart embedding also has been shown to produce a sufficient reduction in rotation− vibration coupling to allow rotationally excited states of floppy molecules to be described using fixed-node DMC with the nodes placed in the Euler angles that relate the space-fixed and Eckart frames.38,49−51 It is worth noting that the DMC wave function and rotation−vibration Hamiltonian used in this study both include rotational degrees of freedom. Specifically, the spatial component of each walker as well as the kinetic energy operator in Ĥ J=0 include all 3N Cartesian coordinates, while the rotational state vectors and ĥrot exist in the Hilbert space of the |J,K⟩ rotational state vectors. While this appears to introduce a double-counting of the contributions from the rotational degrees of freedom, this is not the case. The spatial part of the DMC simulation provides the lowest energy state based on the sum of the potential energy and ⟨ĥrot⟩, neither of which depend on the values of the Euler angles. As a result, the wave function is separable into contributions that depend on the rotational and vibrational coordinates, and the rotational component corresponds to the J = 0 solution of the rotational Hamiltonian. This can be demonstrated by projecting the DMC wave functions onto the Euler angles, as is shown in the Supporting Information for the ground and a representative state of H3+ with J = 10. As can be seen, the DMC wave functions are isotropic to within the statistical noise of the data. While results are presented for only two states, similar distributions are obtained for all of the rotationally excited states considered in this study. As a result, while the three rotational coordinates are included in the definition of the

J

Cn , J , K (τ + δτ ) = Cn , J , K (τ ) − δτ

∑

Cn , J , K ′(τ ) ×

K ′=−J

̂ (R⃗ (τ + δτ ))|J , K ′⟩ ⟨J , K |hrot n

(11)

This approximate treatment of the rotational component of the propagator is reasonable in the limit of small δτ, and its validity in practice can be assessed by monitoring for the dependence of the DMC results on the size of δτ. Once the rotational state vector has been updated and subsequently renormalized, it can be used to determine a local rotational energy associated with the instantaneous position of the walker in configuration space and the chosen value for J according to ̂ (R⃗ (τ + δτ ))|Φ (τ + δτ )⟩ Erot(R⃗ n(τ + δτ )) = ⟨Φn , J (τ + δτ )|hrot n n,J (12)

The potential energy component of the propagator updates the weight of the walker in the DMC ensemble in the same manner as above but with Eref(τ) now compared to the sum of the potential energy and the local rotational energy, so that eq 6 becomes Wn(τ + δτ ) = exp[−(V (R⃗ n(τ + δτ )) + Erot(R⃗ n(τ + δτ )) − Eref (τ ))]Wn(τ )

(13)

The walkers with negligible weights are then identified and replaced using the same procedure as in the ground state DMC calculations described above, and the reference energy is updated according to Eref (τ ) = V̅ (τ ) + Erot ̅ (τ ) − α

Wtot(τ ) − Wtot(0) Wtot(0)

(14)

A DMC simulation performed using this approach will converge onto the lowest energy state consistent with the chosen value of J. Furthermore, projections of the DMC probability distribution onto the vibrational coordinates can provide insights into the strength and nature of the effects of rotation on the molecular structure. To see this, note that, for each walker in the DMC ensemble, ĥrot(R⃗ ) and therefore the evolution of the rotational state vector over the time step as well as the value of Erot(R⃗ ) depend on the position of the walker in configuration space. Likewise, the weight of the walker in the DMC ensemble depends on the local rotational energy through eq 13. This feedback between the spatial and the rotational state vector components of the walkers allows the DMC calculations to capture the effects of rotation−vibration coupling, to within the accuracy of the potential energy surface D



Article

the weight for each state shared evenly between the two walkers. Computational Details. Unless otherwise specified, the DMC simulations performed in this study used δτ = 5 au, α = (2δτ)−1 Hartree, and Nwalker = 40 000. The DMC calculations were initialized with all of the walkers having the same molecular geometry obtained by scaling the Cartesian coordinates of the equilibrium geometry by a factor of 1.1− 1.15. All of the weights were initially set equal to 1. The initial ̂ evaluated in rotational state vectors were the eigenstates of hrot the Eckart frame at the initial geometry. Unless specifically noted, the reference structure for the Eckart embedding was the equilibrium geometry,77,78 and the rotational state vectors were expanded in the symmetric top basis associated with projecting J ⃗ onto the C-axis of the reference structure. The calculations performed in this study included the ground state along with all of the rotationally excited states with J less than or equal to the specified Jmax. The nth walker was replaced if Wn,max fell below Wthresh = N−1 walker. In practice, this led to 1−1.5% of the walkers being replaced each time step independent of the chosen Nwalker, δτ, or Jmax. The DMC simulations were equilibrated for 45 000−50 000 au before collecting information. The reference energies were then averaged over a 37 500−38 500 au production period to obtain the energy for each of the states considered in the calculation. The average energies obtained from five statistically independent DMC calculations were then averaged together to obtain the overall average energies reported in the following sections. In doing so, the rotational excitation energies obtained in each DMC calculation were directly averaged. This led to a significant reduction in the statistical uncertainties of the overall average rotational excitation energies relative to taking the difference between the average energies of the ground and rotationally excited states, evaluated from two distinct sets of DMC calculations. Throughout the production portion of the DMC calculations, descendent weighting was used to produce 15 statistically independent calculations of |Ψ2|, each of which lasted 250 au and was separated from the next by 2250−2500 au.31 The reported average expectation values and projections of the probability distributions onto the vibrational coordinates were obtained by averaging together the 75 statistically independent calculations of |Ψ2| obtained from the five separate DMC calculations. In practice, the rotational state vectors used throughout this work were expanded in the symmetrized symmetric top basis defined via

(PES) and the molecular Hamiltonian used in the imaginary time propagator. Simultaneous Calculation of Multiple Rotationally Excited States within a Single DMC Simulation. The dressed DMC ensemble introduced in eq 7 can be generalized to simultaneously describe Nstate rotational states by dressing each walker with an array of Nstate rotational state vectors ⎛ ⎞ W n(1)(τ )|0, 0⟩ ⎜ ⎟ ⎜ ⎟ (2) (2) τ τ |Φ ⟩ W ( ) ( ) n n,J ⎜ ⎟ Nwalker ⎜ ⎟ ⟨R⃗|Ψ(τ )⟩ = ∑ (R⃗ − R⃗ n(τ )) ⎜ W n(3)(τ )|Φ(3) n , J (τ )⟩ ⎟ n=1 ⎜ ⎟ ⋮ ⎜ ⎟ ⎜⎜ (Nstate) ⎟ (Nstate) (τ )|Φn , J (τ )⟩⎟⎠ ⎝W n max (15)

These rotational state vectors are each associated with a welldefined value of J ranging from 0 to Jmax and are expanded in a symmetric top basis in the same manner as above. Note that the first state included in eq 15 represents the rotational ground state and therefore is associated with the constant rotational state vector |0,0⟩ for all of the walkers. As shown in eq 15, each walker carriers a different weight for each of the states, with W(i) n being the weight associated with the ith state for the nth walker. This allows us to obtain a Monte Carlo sampling of the wave functions for several states, as well as the corresponding energies of these states, based on a single distribution of walkers in configuration space. The numerical stability of this correlated sampling requires that the wave functions all have amplitude in the same regions of configuration space. This is necessary to ensure that a significant fraction of the DMC ensemble contributes to the description of each of the states. Likewise, to extract the energies of multiple states from a single calculation, Nstate state-dependent reference energies are evaluated at each time step in the DMC simulation. Finally, we note this is not the first application of correlated sampling to DMC simulations. Earlier studies have used similar approaches to obtain tunneling splittings in the water dimer and trimer, to study vibrational frequency shifts in HF due to rare gas solvation, and to more efficiently calculate gradients of the electronic energy.52−55 The procedure for propagating the multistate DMC ensemble given in eq 15 is very similar to that described in the previous subsection. Focusing on a single walker for convenience, after randomly displacing the walker in configuration space, the local rotational Hamiltonian is constructed and the rotational state vectors each independently evolved according to eq 11. The rotational state vectors are then orthogonalized using the modified Gram−Schmidt procedure and normalized.80 The local rotational energies are then calculated using eq 12, and the weight for each state is updated using eq 13 and the appropriate state-dependent reference energy. One complication that arises when multiple weights are associated with each walker is how to identify which walkers should be replaced because they are making a negligible contribution to the overall DMC ensemble. In this study, the maximum state-dependent weight associated with each walker is identified at the end of each time step, Wn,max. If Wn,max < Wthresh, then that walker is removed from the DMC ensemble and replaced with the walker with the largest value of Wn,max,

|J , K ⟩p =

1 2p(1 + δK ,0)

[|J , K ⟩ + p( −1)K |J , −K ⟩] (16)

where p = ±1 denotes the parity of the basis function, and the expanded rotational state vectors for the ith state in eq 15 become 81

|Φ(ni,)J (τ )⟩ =

∑ Cn(i,)J ,K ,p(τ)|J , K ⟩p K ,p

(17)

̂ are all real, With this choice of basis, the matrix elements of hrot preventing the need for complex expansion coefficients in the description of the rotational state vectors. Additionally, by using basis functions with well-defined parity, it is straightforward to calculate the average parity of the rotational states calculated in this study. Moreover, one can choose initial rotational state E



Article

Table 1. Comparison between the Multistate DMC, Rigid-Rotor, Fixed-Node DMC,49 and Benchmark Variational82,56 Energies (in cm−1) of Representative Rotationally Excited States of H3+ statea

Evar.b

|1,1⟩0.97 |1,1⟩−0.97 |1,0⟩1.00 |2,2⟩1.00 |2,2⟩−0.99 |2,1⟩0.75 |2,1⟩−0.73 |3,3⟩−0.99 |3,3⟩0.99 |3,2⟩−0.95 |3,2⟩0.97 |3,1⟩0.35 |3,1⟩−0.37 |3,0⟩0.97 |5,5⟩−0.99 |5,5⟩0.99 |5,1⟩−0.05 |5,1⟩0.01 |5,0⟩0.90 |10,10⟩0.95 |10,10⟩−0.95 |10,1⟩0.03 |10,1⟩−0.09

64.13 64.13 86.96 169.30 169.30 237.36 237.36 315.35 315.35 428.02 428.02 494.77 494.77 516.88 729.09 729.09 1250.45 1250.45 1271.41 2451.77 2451.77 4348.79 4348.79

EDMC − Evar.c −0.25 −0.12 0.14 −1.07 −1.07 −0.19 0.76 −2.53 −2.54 −0.98 −0.85 1.03 4.46 6.52 −8.44 −8.44 11.74 30.96 36.82 −36.53 −36.53 272.83 352.01

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.45 0.55 0.59 1.03 1.07 1.87 1.60 2.05 2.05 3.12 2.96 3.87 3.75 3.62 5.64 5.64 14.89 6.86 9.84 8.74 8.74 21.26 16.74

100|EDMC − Evar.|/Eexpt

ERR − Evar.d

EDMC,FN − Evar.e

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

−0.36 −0.36 −0.09 −1.09 −1.09 0.15 0.15 −2.03 −2.03 0.80 0.80 3.36 3.36 4.34 −3.55 −3.55 29.52 29.52 31.66 16.01 16.01 406.05 406.05

0.84 ± 5.23 NA 1.33 ± 2.01 NA NA NA NA 0.41 ± 4.13 NA −4.28 ± 1.88 NA −5.60 ± 2.19 NA 0.14 ± 2.02 −1.80 ± 4.06 NA NA NA NA −6.79 ± 7.55 NA 8.23 ± 6.31 NA

0.39 0.19 0.16 0.63 0.63 0.08 0.32 0.80 0.81 0.23 0.20 0.21 0.90 1.26 1.16 1.16 0.94 2.48 2.90 1.49 1.49 6.27 8.09

0.70 0.86 0.68 0.61 0.63 0.79 0.67 0.65 0.65 0.73 0.69 0.78 0.76 0.70 0.77 0.77 1.19 0.55 0.77 0.36 0.36 0.49 0.38

The states are labeled using the notation |J,K⟩⟨p⟩, where ⟨p⟩ represents the average parity of the state obtained in the multistate DMC calculation. The benchmark variational energies, Evar., are reported relative to the zero-point energy. cThe multistate DMC calculations were performed with Jmax = 3 for the states with J ≤ 3 and with Jmax = 10 for the J > 3 states. The reported error bars are the 99% confidence intervals of the DMC energies. dThe rigid-rotor energies, ERR, are the eigenvalues of eq 18 constructed using the zero-point vibrationally averaged rotational constants ⟨A⟩0 = ⟨B⟩0 = 43.44 cm−1 and ⟨C⟩0 = 20.33 cm−1.49 eThe reported fixed-node DMC energies, EDMC,FN, of the J = |K| states as well as the |10,1⟩+ state are extrapolations to the δτ = 0 limit, and the error bars are the largest combined 99% confidence intervals of the ground and rotationally excited state in the data set used to do the extrapolation. For the other states, δτ = 1 au and the error bars are the combined 99% confidence intervals.49 a b

perform calculations with Jmax = 3 or 10, corresponding to including 16 and 121 states in the calculations respectively, to a J = 0 calculation with the same number of walkers. The Jmax = 3 calculations for CH5+ took 1.085 times as long as a J = 0 calculation, while the Jmax = 10 calculation took 2.652 times as long as the corresponding J = 0 calculation. Note that this represents highly sublinear scaling with respect to the number of rotational states included in the multistate DMC calculations.

vectors that are eigenstates of both the parity operator and the local ĥrot. These initial conditions have been used in this study. We now summarize the PES’s used in this study as well as the literature benchmark data that our DMC energies are compared against. The DMC calculations of H3+ and H2D+ were performed using the PES developed by Aguado and coworkers.56 Unfortunately, variational calculations using this surface have only been reported for the states of H3+ with J = 5 and J = 10. Benchmark variational energies for the other states of H3+ considered in this study are obtained from calculations performed by Dinelli et al. on their DPT surface.82,83 The benchmark energies for all of the H2D+ states are obtained from variational calculations performed by Sochi and Tennyson using Polyansky et al.’s PES.20,22 Importantly, all three studies were performed to spectroscopic accuracy, and for those states of H3+ where the data sets overlap, the agreement between the reported energies is consistently well-within the statistical uncertainties of our DMC calculations. The H3O+ calculations were performed using the full-dimensional PES developed by Bowman and co-workers.59 Because no variational calculations have been reported for rotationally excited states of H3O+ with J > 1, the benchmark literature data used in this study is taken from a fit to high-resolution experimental data.61 Finally, the CH5+ calculations were performed using the global extension to the full-dimensional PES developed by Bowman and coworkers.34 Interestingly, within the multistate DMC calculations, the most time-consuming step is the evaluation of the potential surface. This can be seen by comparing the time required to

■

RESULTS AND DISCUSSION Using the Multistate DMC Method To Calculate the Rotational Excitation Energies of H3+, H2D+, and H3O+. Tables 1−3 provide an overview of the DMC energies of the rotationally excited states with J ≤ 3 of H3+, H2D+, and H3O+ obtained using the multistate DMC approach developed in this study with Jmax = 3 in eq 15. The agreement between the DMC energies and the literature benchmark values is excellent, with | EDMC − Ebenchmark| less than 3 cm−1 for the majority of the states and the largest such error being 6.52 cm−1 for the |3,0⟩+ state of H3+. For H3+ and H2D+, |EDMC − Ebenchmark| amounts to less than 1% of the overall rotational excitation energy for all but one state, whereas for H3O+, the absolute percent errors, 100| EDMC − Ebenchmark|/Ebenchmark, range from 1.13% to 2.51%. The 99% confidence intervals of the DMC energies are consistently less than 3.9 cm−1 for all three species and are less than 2 cm−1 for all of the H3O+ and H2D+ states. Note that, unlike H3+ and H2D+, the deviations of the H3O+ DMC energies from the F



Article

Table 2. Comparison between the Multistate DMC, Rigid-Rotor, Fixed-Node DMC,48,86 and Benchmark Experimental61 Energies (in cm−1) of Representative Rotationally Excited States of H3O+ statea

Eexptb

|1,1⟩−1.00 |1,1⟩1.00 |1,0⟩1.00 |2,2⟩1.00 |2,2⟩−1.00 |2,1⟩0.99 |2,1⟩−0.99 |2,0⟩1.00 |3,3⟩−1.00 |3,3⟩1.00 |3,2⟩−1.00 |3,2⟩1.00 |3,1⟩0.98 |3,1⟩−0.98 |3,0⟩1.00 |6,6⟩1.00 |6,6⟩−1.00 |6,1⟩0.12 |6,1⟩−0.07 |6,0⟩0.96 |10,10⟩1.00 |10,10⟩−1.00 |10,1⟩−0.04 |10,1⟩−0.02 |10,0⟩0.86

17.38 17.38 22.50 47.02 47.02 62.37 62.37 67.47 88.92 88.92 114.47 114.47 129.76 129.76 134.86 288.31 288.31 465.34 465.34 470.35 726.74 726.74 1217.28 1217.28 1222.11

EDMC − Eexptc −0.30 −0.29 −0.56 −0.63 −0.63 −1.41 −1.42 −1.66 −1.01 −1.01 −2.31 −2.31 −3.03 −3.02 −3.22 −3.12 −3.12 −10.38 −9.30 −9.23 −8.74 −8.74 −22.19 −16.17 −15.57

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

100|EDMC − Eexpt|/Evar.

ERR − Eexptd

EDMC,FN − Eexpte

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

−0.29 −0.29 −0.55 −0.61 −0.61 −1.38 −1.38 −1.63 −0.95 −0.95 −2.20 −2.20 −2.93 −2.93 −3.16 −2.28 −2.28 −9.26 −9.26 −9.42 −5.39 −5.39 −14.92 −14.92 −14.89

2.47 ± 3.34 1.45 ± 3.17 0.93 ± 3.62 1.34 ± 2.81 0.12 ± 2.41 0.52 ± 2.48 0.24 ± 2.77 2.24 ± 3.25 NA NA NA NA NA NA NA 0.96 ± 2.19 1.08 ± 2.19 NA NA 22.13 ± 2.62 −0.13 ± 2.13 0.61 ± 2.15 NA NA 84.93 ± 2.47

0.11 0.12 0.16 0.33 0.33 0.45 0.38 0.47 0.70 0.70 0.69 0.68 0.93 0.88 0.93 2.47 2.47 4.81 4.69 4.02 6.30 6.30 4.38 7.12 5.83

1.73 1.69 2.51 1.34 1.34 2.26 2.27 2.46 1.13 1.13 2.02 2.02 2.34 2.33 2.39 1.08 1.08 2.23 2.00 1.96 1.20 1.20 1.82 1.33 1.27

0.63 0.69 0.70 0.70 0.70 0.72 0.62 0.70 0.79 0.79 0.60 0.59 0.72 0.68 0.69 0.86 0.86 1.03 1.01 0.85 0.87 0.87 0.36 0.58 0.48

The states are labeled using the notation |J,K⟩⟨p⟩, where ⟨p⟩ represents the average parity of the state obtained in the multistate DMC calculation. The benchmark variational energies, Evar., are reported relative to the zero-point energy. cThe multistate DMC calculations were performed with Jmax = 3 for the states with J ≤ 3 and with Jmax = 10 for the J > 3 states. The reported error bars are the 99% confidence intervals of the DMC energies. dThe rigid-rotor energies, ERR, are the eigenvalues of eq 18 constructed using the zero-point vibrationally averaged rotational constants ⟨A⟩0 = ⟨B⟩0 = 10.975 cm−1 and ⟨C⟩0 = 6.116 cm−1. eThe reported fixed-node DMC energies, EDMC,FN, were calculated with δτ = 1 au, and the reported error bars are the combined 99% confidence intervals of the ground and rotationally excited states.48,86 a b

This is especially pronounced for the high J, low |K| states of H3+ where Coriolis coupling terms involving the doubly degenerate bending mode are known to cause significant rotation−vibration mixing.49,84,85 Indeed, many of the states of H3+ with particularly large values of |EDMC − Ebenchmark| have rotational excitation energies that are comparable to or larger than the energy of the doubly degenerate bend fundamental, 2521.42 cm−1.56 The statistical uncertainties of the DMC energies of the J > 3 states reported in Tables 1−3 are also generally larger than those of the J ≤ 3 states, with the 99% confidence intervals for the EDMC for H3+ and H2D+ being much larger than those for H3O+. However, as shown in the Supporting Information, the statistical uncertainties of the DMC energies of the J ≤ 3 states of H3+ and H2D+ are also significantly increased when Jmax is increased from 3 to 10. In order to understand this, we calculated the average fraction of walkers in the DMC ensemble that contribute to the DMC wave function of each state with weights that exceed 0.01. When Jmax = 3, this amounts to 70− 72% of the walkers contributing to a given state of H3+ or H2D+. However, when Jmax = 10, only 10−34% of the walkers have weights for a given state that are greater than 0.01. This dramatic reduction in the effective number of walkers that contribute to the DMC wave function of each state of H3+ and H2D+ results in the increased statistical uncertainties in the DMC energies. In contrast, for H3O+, 74−75% of the walkers contribute to the description of each state when Jmax = 3, while

benchmark values are consistently outside of the statistical uncertainties of the calculations. However, as the benchmark energies for H3O+ are obtained from a fit of a model Hamiltonian to experimental data, |EDMC − Ebenchmark| for this species reflect a combination of deficiencies in the DMC calculations and in the PES, and it is difficult to disentangle these two potential sources of error. It is worth noting that variational calculations of the J = 1 states of H3O+ have been reported using the same PES and these values, 17.4 cm−1 and 22.3 cm−1 for the |1,1⟩± and |1,0⟩ states, respectively, agree well with the experimental energies.59 The energies of representative states of H3+, H2D+, and H3O+ with J > 3 obtained from multistate DMC calculations performed with Jmax = 10 in eq 15 are also presented in Tables 1−3, while the energies for all states with J ≤ 6 and J = 10 can be found in the Supporting Information. For all three species, |EDMC − Ebenchmark| generally increases with rotational excitation energy. For H3O+ and to a large extent also H2D+, this increase in the errors of the J > 3 DMC energies scales with EDMC so that the absolute percent errors generally remain in the range from 1% to 2.5%. In contrast, for the high J states of H3+, the |EDMC−Ebenchmark| values are found to grow faster than EDMC so that, for example, the errors in the DMC energies of the |5,0⟩+ and |10,1⟩− states amount to 2.90% and 8.09% of the rotational excitation energy, respectively. These errors in the DMC energies are indicative of the presence of rotation−vibration coupling that is not fully accounted for in the present approach. G



Article

Table 3. Comparison between the Multistate DMC, Rigid-Rotor, Fixed-Node DMC,49,50 and Benchmark Variational22 Energies (in cm−1) of Representative Rotationally Excited States of H2D+ statea

Evar.b

10,1(1.00) 11,1(−1.00) 11,0(1.00) 20,2(1.00) 21,2(−0.99) 21,1(0.99) 22,1(−0.99) 22,0(1.00) 30,3(0.99) 31,3(−0.99) 31,2(0.99) 32,2(−0.97) 32,1(0.98) 33,1(−0.99) 33,0(0.99) 61,5(0.96) 63,3(0.92) 64,2(0.88) 65,1(0.87) 101,9(0.88) 103,7(0.77) 105,5(0.73) 107,3(0.61) 109,1(0.67) 1010,0(0.85)

45.70 60.03 72.45 131.63 138.84 175.93 218.65 223.86 251.37 254.02 326.15 354.77 376.35 458.35 459.84 991.01 1206.12 1302.33 1455.04 2301.80 2790.09 3048.45 3387.85 3905.28 4198.58

EDMC − Evar.c −0.10 −0.02 0.20 −0.51 −0.54 0.16 0.51 1.00 −1.29 −1.47 −0.18 0.04 1.59 1.69 2.42 −5.56 7.83 17.37 18.01 −19.34 17.31 70.55 112.64 108.66 84.55

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

100|EDMC − Evar.|/Evar.

ERR − Evar.d

EDMC,FN − Evar.e

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

−0.36 −0.27 −0.14 −0.99 −0.97 −0.39 0.13 0.29 −1.82 −1.84 −0.22 0.03 0.93 2.77 2.89 4.06 21.14 25.71 35.02 30.56 103.08 164.03 200.72 295.69 374.51

0.35 ± 2.02 −0.92 ± 2.04 −0.75 ± 1.95 NA NA NA NA NA 0.49 ± 5.12 NA −1.77 ± 5.79 NA −0.19 ± 5.00 NA −1.35 ± 6.17 −0.67 ± 6.73 −1.02 ± 5.58 −3.91 ± 5.30 −0.34 ± 7.10 1.83 ± 4.62 −20.88 ± 7.11 −16.08 ± 6.56 −11.54 ± 5.47 −0.47 ± 9.12 0.99 ± 6.56

0.48 0.25 0.46 1.16 0.93 1.48 1.04 1.21 1.92 1.79 2.61 2.30 2.65 2.69 2.75 9.25 10.59 7.32 12.31 9.61 7.85 11.73 7.93 12.14 10.43

0.22 0.04 0.28 0.39 0.39 0.09 0.23 0.45 0.51 0.58 0.06 0.01 0.42 0.37 0.53 0.56 0.65 1.33 1.24 0.84 0.62 2.31 3.32 2.78 2.01

1.05 0.42 0.63 0.88 0.67 0.84 0.48 0.54 0.76 0.70 0.80 0.65 0.70 0.59 0.60 0.93 0.88 0.56 0.85 0.42 0.28 0.38 0.23 0.31 0.25

The states are labeled using the notation JKA,KC(⟨p⟩), where ⟨p⟩ represents the average parity of the state obtained in the multistate DMC calculation. bThe benchmark variational energies, Evar., are reported relative to the zero-point energy. cThe multistate DMC calculations were performed with Jmax = 3 for the states with J ≤ 3 and with Jmax = 10 for the J > 3 states. The reported error bars are the 99% confidence intervals of the DMC energies. dThe rigid-rotor energies, ERR, are the eigenvalues of eq 18 constructed using the zero-point vibrationally averaged rotational constants ⟨A⟩0 = 43.361 cm−1, ⟨B⟩0 = 28.949 cm−1, and ⟨C⟩0 = 16.389 cm−1. eThe reported fixed-node DMC energies, EDMC,FN, of the J > K states are extrapolations to the δτ = 0 limit, and the error bars are the largest combined 99% confidence intervals of the ground and rotationally excited state in the data set used to do the extrapolation. For the J = 1 states, δτ = 1 au and the error bars are the combined 99% confidence intervals.49,50 a

62−71% contribute when Jmax = 10. As a result, the statistical uncertainties of the H3O+ energies are significantly less than the other species and exhibit a smaller dependence on Jmax. The greater fraction of walkers that contribute to the DMC wave functions for H3O+ reflects the smaller change in the portion of configuration space that is sampled by H3O+ upon rotational excitation compared to H3+. The reason for this as well as the observation that the multistate DMC calculations of the rotationally excited states of H3+ and H2D+ behave somewhat differently when Jmax = 3 than when Jmax = 10 can be rationalized by considering the relative importance of the local rotational energy and the potential energy in determining the evolution of the weights for each rotational state. When Jmax = 3, the rotational excitation energies are all much less than the zero-point energies of H3+ and H2D+, 4362.07 ± 3.78 cm−1 and 3977.01 ± 2.68 cm−1, respectively. The local rotational energies are therefore much less than the typical values of the potential energy sampled by the walkers and eq 13 is dominated by the potential energy, allowing a walker to have similar weights for many of the rotational states. When Jmax = 10, high J states are added to the calculation with rotational excitation energies that are comparable to or larger than the zero-point energies. For the low J states, eq 13 will still be dominated by the potential energy, whereas for the high J states, the local rotational energy and potential energy will make comparable contributions. This will lead to a greater dispersion among the state-dependent

weights carried by each walker and therefore a reduction in the number of rotational states that each walker effectively contributes to. Note that this behavior is not nearly as pronounced in H3O+ due to the fact that all of the states considered in this study have rotational excitation energies that are significantly less than the 7453.90 ± 5.84 cm−1 zero-point energy. For the symmetric top molecules H3+ and H3O+, an internal check on the reliability of the method is provided by the comparison of the DMC energies of the members of the degenerate pairs |J,K⟩±. For both species, the pairs of states with K ∼ J are consistently calculated to be degenerate to well-within the statistical uncertainties of the DMC energies. However, statistically significant differences are found between the EDMC values for the pairs of H3+ states with low values of K. The extent to which the degeneracy is spuriously broken increases with J and is most pronounced when K = 1. To identify the cause of this behavior, we performed a series of DMC calculations for H3+ in which the rotational state vectors in eq 17 were replaced by a single parity adapted symmetric top eigenstate from eq 16, and the values of the C(i) n,J,K,p were not allowed to evolve in imaginary time. Note that, by not allowing the rotational state vectors to evolve, all off-diagonal terms in the local rotational Hamiltonian that couple different symmetric top basis functions were effectively set equal to zero. The results of these DMC calculations are reported for all states with J ≤ 3 and selected states with J = 5 and 10 in the H



Article

H3O+, the rigid-rotor energies remain comparable to the multistate DMC energies even for states with J = 10, consistent with the magnitude of rotation-vibration coupling in H3O+ being much less than in H3+ and H2D+. A comparison with the results of our previously performed fixed-node DMC calculations of select rotationally excited states of H3+, H2D+, and H3O+ is also included in Tables 1−3.48−50,86 Before discussing these comparisons in detail, it is worth recalling the differences between the approximate treatments of rotational excitation in the two DMC methods. In fixed-node DMC, rotational excitation is imposed by introducing a specific set of approximate nodal surfaces consistent with the rigid-rotor wave function of the state of interest, with errors occurring when these approximate nodal surfaces differ from those of the true molecular eigenstate. In contrast, the multistate DMC approach developed in this study makes use of an approximate rotation−vibration Hamiltonian that neglects Coriolis coupling terms. Returning to the data, the fixed-node DMC energies, EDMC,FN, for the states with J ≤ 3 are found to be comparable to the multistate DMC values with neither DMC method consistently producing results in better agreement with the literature benchmark values. By directly averaging the rotational excitation energy, rather than computing the differences between the average energies for two different states that were obtained from separate DMC calculations, the multistate DMC energies tend to have somewhat smaller statistical uncertainties than are present in the fixed-node DMC values. Turning to the states with J > 3, for H3+ and H2D+, the fixed-node DMC energies are nearly consistently in better agreement with the benchmark values than the multistate DMC energies. It is worth noting, though, that this improved accuracy comes at substantial computational cost, with 20 separate fixed-node DMC calculations and an extrapolation to the δτ = 0 limit performed for each of the J > 3 states of H3+ and H2D+ reported in Tables 1 and 2. In contrast the computational time required for the multistate DMC method scales sublinearly with the number of states included in the calculation. For H3O+, whereas the |EDMC,FN − Eexpt| are smaller than the |EDMC − Eexpt| for the |6,6⟩± and |10,10⟩± states, the multistate DMC energies for the |6,0⟩+ and |10,0⟩+ states are in better agreement with the literature benchmark values than the EDMC,FN. We will return to this point below. We also investigated the dependence of the multistate DMC results on the size of the time step, the number of walkers in the DMC ensemble, and the details of the Eckart embedding. Both the split-operator approximation made in eq 9 and the use of first-order time-dependent perturbation theory to evaluate the evolution of the rotational state vectors are valid only in the limit of small δτ. To verify that these approximations are not introducing errors in our multistate DMC calculations, we reduced the time-step from δτ = 5 au to δτ = 0.5 au and found no statistically significant changes to the DMC energies, although the statistical uncertainties were somewhat reduced. Reducing the size of the DMC ensemble to Nwalker = 20 000 did not have a statistically significant affect on the calculated rotational excitation energies but did increase their statistical uncertainties. Finally, in our previously performed fixed-node DMC calculations of the rotationally excited states of H3O+, the results were found to depend strongly on the choice of static reference structure used to define the Eckart frame.48,51 Specifically, use of the C3v minimum energy structure led to the calculation of members of the manifold of rotationally excited states associated with the upper member of the

Supporting Information. They were found to be in good agreement with the multistate DMC energies presented in Table 1 and in the Supporting Information for all but the low K states of H3+. For these states, the results obtained with the simplified DMC calculations using frozen rotational state vectors were found to be in better agreement with the benchmark variational values than those calculated with evolving rotational state vectors. These results suggest that, for symmetric top molecules, the effects of the off-diagonal terms in the local rotational Hamiltonian should vibrationally average to zero. For the states of H3+ with high values of J and low values of K, the numerical vibrational averaging performed in the multistate DMC simulations is insufficient to fully accomplish this, and as a result, the degeneracy of these states is broken in the calculations. That this is particularly severe for the |J,1⟩± states reflects the fact that the degenerate |J,1⟩± basis functions are directly coupled by the term that is proportional ̂ ̂ −1 ̂ to JB̂ I−1 BCJC + JCIBCJB in eq 10. While the coefficient in front of this term is zero in the equilibrium geometry of the molecule, and must be zero on average, it deviates from zero when the molecule is displaced to a lower-symmetry structure. To test the effects of the choice of the rotational basis on the calculated energies, we repeated the calculations for H3+ with Jmax = 10, using a basis that is expanded in rotational state vectors that are eigenstates of JÂ rather than JĈ . While this may seem at first to be a poor choice for an oblate top molecule, the results of the calculations on H3+ using the A-axis to define the K quantum number are in comparable agreement with the variational results to those reported in Table 1 and do not display as large a splitting of the KC = 1 energies as was found when the C-axis defined K. The results of these calculations can be found in the Supporting Information. Before concluding this part of the discussion, we note that the absence of broken degeneracies of the DMC energies for H3O+ suggests that these problems may be relatively unique to systems, like H3+, that possess unusually large rotational constants and lack a central heavy atom to constrain the nuclear dynamics. Indeed, within the statistical uncertainties of the calculations, the DMC energies for H3O+ do not exhibit any dependence on whether the rotational state vectors were held frozen at a single symmetric top basis function or allowed to evolve during the calculation. Tables 1−3 also present the rigid rotor energies, ERR, which are the eigenvalues of the Hamiltonian 2 2 2 ĤRR = ⟨A⟩0 JÂ + ⟨B⟩0 JB̂ + ⟨C⟩0 JĈ

(18)

where ⟨A⟩0, ⟨B⟩0, and ⟨C⟩0 are the zero-point vibrationally averaged rotational constants obtained by averaging the inverse moment of inertia tensor in the Eckart frame over the DMC ground state probability distribution.36,38 For the rotationally excited states with J ≤ 3, the rigid-rotor and multistate DMC energies are found to be of comparable agreement with the literature benchmark values for all three species. Indeed, for many of the J ≤ 3 states, the ERR values are identical to the EDMC values to within the statistical uncertainties of the DMC calculations. We therefore see that, even for molecules that undergo large-amplitude, zero-point vibrational motions, eq 18 can be readily used to predict the energies of rotationally excited states with low J to within a few cm−1. For the states of H3+ and H2D+ with J > 3, the errors of the rigid-rotor energies, |ERR − Evar.|, grow rapidly with rotational energy and greatly exceed |EDMC − Evar.| for many of the states. In contrast, for I



Article

tunneling doublet. If instead, the reference structure was the D3h saddle point connecting the two wells along the umbrella inversion mode, then pure rotationally excited states were obtained from the fixed-node DMC calculations.48,51 In this study, we found that the results obtained using the multistate DMC approach did not exhibit a statistically significant dependence on whether the C3v or the D3h geometry was used as the reference structure for the Eckart embedding. Multistate DMC Probability Distributions of the Rotationally Excited States of H3+, H2D+, and H3O+. We begin this subsection by considering the extent to which the rotational state vectors obtained using the multistate DMC approach are parity eigenstates. Although eq 10 contains terms that couple basis functions with different parities, the overall effect of these terms should vibrationally average to zero as the true molecular eigenstates are also eigenstates of the parity operator. The expectation value of the parity, ⟨p⟩, is provided for each of the states reported in Tables 1−3 and in the Supporting Information. For the majority of the states with J ≤ 3, ⟨p⟩ is found to be approximately ±1. The |2,1⟩± and |3,1⟩± states of H3+ represent the only exception to this. In these cases the calculated eigenstates are admixtures involving only the |J,1⟩+ and |J,1⟩− basis states. As J is increased, the degree of mixing between the parity blocks becomes more and more pronounced and is found for states with K > 1. This is reflected in a reduction of the calculated values of |⟨p⟩|. This is especially true for the symmetric top molecules where many of the J = 10 states have |⟨p⟩| ≪ 0.5. In contrast, much less parity contamination is present in the rotational state vectors of H2D+, with all of the states with J ≤ 6 having |⟨p⟩| > 0.85 and only a single J = 10 state having |⟨p⟩| < 0.5. We now consider the effects of rotational excitation on the vibrational probability distributions of H2D+. As in our previous fixed-node DMC study of this species, rotational excitation is found to result primarily in a shift in the position of the peak in these distributions.50 To highlight the overall trends, Figure 1 focuses on the expectation values of the bond lengths and bond angles for representative J = 10 rotationally excited states along with the ground state. In the high KC limit, where the rotational excitation is primarily about the axis perpendicular to the molecular plane, centrifugal distortion shifts all three atoms away from the center of mass, leading to an increase in both bond lengths and little change to the overall molecular shape. The larger increase in the average HD bond length, ⟨rHD⟩, relative to the average HH bond length, ⟨rHH⟩, reflects the fact that the center of mass lies closer to the deuterium atom than the two hydrogen atoms. Although the calculations were performed in a |J,KC⟩p basis, H2D+ is a highly asymmetric top molecule, and in some cases it is convenient to consider the dependence of the results on the KA state labels. For the high KA states, rotational excitation leads to a much larger increase in ⟨rHH⟩ than in ⟨rHD⟩, resulting in an increase in the HDH bond angle and an associated decrease in the HHD bond angles. This is consistent with a simple picture of centrifugal distortion in which rotational motion about the A-axis; that is, the axis bisecting the HDH bond angle, will primarily lead to an increase in the H−H distance. For the states with intermediate values of KA, the manifestations of rotation−vibration coupling shown in Figure 1 do not monotonically vary between the KA = 0 and KA = J limits and are not easily rationalized from either zeroth order picture of the effects of centrifugal distortion. This can be understood by noting that the rotational state vectors

Figure 1. (a) Comparison of the expectation values ⟨rHH⟩ (red) and ⟨rHD⟩ (blue) for representative J = 10 states of H2D+ as calculated with multistate DMC (open circles) and fixed-node DMC (closed triangles). Likewise, in panel b, the red data points are for ⟨θHDH⟩, while the blue data points are for ⟨θHHD⟩. In both panels, the J = 0 expectation values obtained from a ground state DMC calculation are reported as the straight lines.50

associated with these intermediate states contain significant contributions from a larger and more diverse set of symmetric top basis functions than the rotational states at either the KA = 0 or KA = J limits. As a result, intuition gleaned from a zeroth-order model of the effects of centrifugal distortion based on having a single symmetric top basis function dominating the expansion breaks down for these intermediate states. Figure 1 also contains the expectation values obtained from our previously reported fixed-node DMC calculations.50 Although there are quantitative differences between the results of these two types of DMC calculations of rotationally excited states, particularly for the states with intermediate values of KA, the two sets of DMC calculations predict the same overall trends for the effects of rotational excitation on the molecular structure of H2D+. Therefore, even though the multistate DMC calculations yield rotational excitation energies for the J = 10 states that are in poorer agreement with the benchmark variational values than the results of the fixed node calculations, they can still provide insights into the nature and strength of rotation-vibration coupling in this fluxional molecule. We performed a similar analysis of the |10,10⟩± and |10,1⟩± states of H3+ and arrived at the same conclusion. It is worth noting that for the |10,1⟩± states of H3+, where the multistate DMC energies were spuriously nondegenerate, performing the multistate DMC calculation using frozen rotational state vectors led to a slight improvement in the agreement of the expectation values with the fixed-node DMC results, although the changes were consistently within the statistical uncertainties. J



Article

The effects of rotational excitation on the molecular structure of H3O+ are summarized in Table 4 where the expectation

symmetry axis primarily leading to a shift in the vibrational probability distribution toward more planar geometries as anticipated from a simple physical picture of centrifugal distortion. Note that the fact that the two DMC methods yield comparable expectation values for the |10,10⟩± states is consistent with the agreement between their DMC energies for these states given in Table 2. For the |10,0⟩ state, where the considerably higher rotational excitation energy and the closer spacing between adjacent rotational levels both tend to suggest stronger rotation−vibration interaction, the two DMC approaches yield qualitatively different pictures of the effects of rotation−vibration coupling on the average molecular structure. The fixed-node DMC calculations suggest that rotation−vibration coupling in the |10,0⟩ state induces similar, albeit somewhat smaller in magnitude, shifts in the vibrational probability distribution as in the |10,10⟩± states. In contrast, the multistate DMC results predict that rotational motion about an axis perpendicular to the C3v axis results in a relatively large increase in ⟨Γ⟩, indicating a shift in the vibrational probability distribution away from planar molecular geometries, along with a decrease in the average HOH bond angles and an elongation of the OH bonds. This latter picture seems more consistent with this state having very different rotational motion and nearly 500 cm−1 more rotational energy than the |10,10⟩± states. This, along with the fact that |EDMC − Eexpt| is significantly less than |EDMC,FN − Eexpt|, suggests that the multistate DMC approach is able to provide a somewhat better description of the rotation-vibration coupling present in the |10,0⟩ state of H3O+ than fixed-node DMC. Results for CH5+. Next we turn our attention to CH5+. CH5+ is a notoriously complex molecule due to the existence of 120 equivalent minima on the potential surface arising from the possible permutations of the five hydrogen atoms. Unlike most organic molecules, the barriers connecting these minima are low, and any one of the 120 minima can be connected to any other minimum through a series of saddle points that are no larger than 1 kcal mol−1 based on the potential surface used in the present calculation.34 Bunker and co-workers evaluated rotationally excited states of CH5+ using a model Hamiltonian that included the effects of torsion and the rotation, which was dressed by the 120 equivalent minima.74 In that work, the splitting due to the barriers between equivalent minima that were not explicitly included in the rotation/torsion Hamil-

Table 4. Effects of Rotational Excitation on the Molecular Structure of H3O+ statea

fixed-node DMCb,c

multistate DMCb ⟨rOH⟩ 0.995 ± 0.002 Å 0.15% 0.15% 0.39% ⟨θHOH⟩ 113.1 ± 0.3° 0.28% 0.28% −0.67% ⟨Γ⟩ 0.244 ± 0.007 Å −2.62% −2.62% 7.08%

|0,0⟩ |10,10⟩1.00 |10,10⟩−1.00 |10,0⟩0.86 |0,0⟩ |10,10⟩1.00 |10,10⟩−1.00 |10,0⟩0.86 |0,0⟩ |10,10⟩1.00 |10,10⟩−1.00 |10,0⟩0.86

0.995 ± 0.002 Å 0.09% 0.10% 0.23% 113.1 ± 0.3° 0.28% 0.28% 0.20% 0.244 ± 0.004 Å −2.73% −2.74% −1.67%

The states are labeled using the notation |J,K⟩⟨p⟩, where ⟨p⟩ represents the average parity of the state obtained in the multistate DMC calculation. bThe ground state expectation values are reported along with their standard deviations. The excited state expectation values are reported as percent differences from the ground state values. c The fixed-node DMC expectation values are from Petit and McCoy.48,86 a

values of the average OH bond length, ⟨rOH⟩, average HOH bond angle, ⟨θHOH⟩, and average planarity distance, ⟨Γ⟩, calculated in representative J = 10 states are compared to the ground state values. Note that the planarity distance, which provides a description of the umbrella inversion mode, is defined as the distance between the oxygen atom and the plane containing the three hydrogen atoms with Γ = 0 indicating a planar geometry. As in the fixed-node DMC calculations, Γ is found to be much more strongly coupled to the rotational motion than the other vibrational coordinates, a reflection of the large-amplitude zero-point motion associated with this mode. For the |10,10⟩± states, the two DMC methods yield comparable expectation values, with rotation about the C3v

Table 5. Comparison between the Multistate DMC and Rigid-Rotora Energies (in cm−1) for Representative Rotationally Excited States of CH5+ J 1 2 3 4 5 6 7 8 9 10

b Emin DMC

7.70 23.15 46.35 77.26 115.71 161.40 214.36 274.58 342.05 416.77

± ± ± ± ± ± ± ± ± ±

0.19 0.54 1.10 1.77 1.97 2.87 4.16 5.61 7.27 9.12

⟨p⟩c

⟨ρ⟩d

1.00 1.00 0.97 −0.92 0.73 0.58 −0.42 −0.31 −0.25 0.23

1.00 1.38 1.82 1.38 1.88 2.24 2.70 3.30 3.90 4.52

b Emax DMC

7.77 23.35 46.62 77.72 116.46 163.04 217.43 279.73 349.50 427.71

± ± ± ± ± ± ± ± ± ±

0.15 0.47 1.04 1.33 1.22 2.62 1.39 2.19 1.02 5.85

⟨p⟩c

⟨ρ⟩d

e Eavg DMC

Eavg RR

1.00 −0.99 −0.94 −0.79 −0.67 −0.48 0.35 0.22 0.12 0.18

1.00 1.02 1.18 1.67 2.60 3.56 4.99 5.86 6.24 6.14

7.75 23.24 46.46 77.41 116.07 162.40 216.41 278.09 347.43 424.46

7.77 23.24 46.44 77.37 116.04 162.45 216.59 278.46 348.07 425.41

The rigid-rotor energies, ERR, are the eigenvalues of eq 18 constructed using the zero-point vibrationally averaged rotational constants ⟨A⟩0 = 3.890 cm−1, ⟨B⟩0 = 3.862 cm−1, and ⟨C⟩0 = 3.849 cm−1. bThe largest and smallest energy evaluated for a given value of J using the multistate DMC calculations with Jmax = 10. The reported error bars are the 99% confidence intervals of the DMC energies. cThe average parity of the DMC state. d The effective number of symmetric top states contributing to the DMC state, see eq 20. eThe average of the calculated 2J + 1 energies for a given value of J. a

K



Article

tonian were assumed to be small, and the rotational constants used in that Hamiltonian led to symmetric top rotational progressions with A ≈ B > C. In a later work,36 we showed that the vibrationally averaged rotational constants, obtained by averaging the Eckart embedded I−1 α,β in eq 10 over the DMC ground state probability amplitude yielded three rotational constants that were effectively equal within the statistical uncertainty of the simulation. These constants provided the basis for fixed-node DMC calculations of J = 1 energies and wave functions for CH5+ and its deuterated analogues,51 using the approach that had been successfully applied to H3O+, H3+, and H2D+.48−50 In the case of CH5+ we were not able to disentangle the rotational excitation from a 10.4 cm −1 tunnelling splitting.24,38,51 Below we present our findings for rotationally excited states of CH5+, based on the multistate DMC approach, described above. The energies for the highest and lowest energy states with J between 1 and 10 are reported in Table 5. The complete list of calculated rotational energies of states with J ≤ 6 and J = 10 is provided in the Supporting Information. All reported energies were obtained from calculations in which Jmax = 10. The most striking feature of these results is the near equality of the average energies obtained from the rigid rotor and multistate DMC calculations. We will return to this point at the end of this section. Further analysis of the multistate DMC results for CH5+ shows that the energies that were obtained for a given value of J are all within the reported statistical uncertainty of the mean value of the energy for the 2J + 1 states that share a J quantum number. In other words, within the uncertainties of the DMC calculations, the v = 0, J-levels are 2J + 1-fold degenerate, as would be expected for a spherical top molecule. This contrasts the symmetric top energy pattern for states with J = 3 reported by Bunker and co-workers,74 which vary from 44.9 to 55.6 cm−1. This 11 cm−1 range is more than an order of magnitude larger than the range of energies spanned by the corresponding multistate DMC energies. While the degeneracy of the various K-levels is consistent with the nearly equal vibrationally averaged rotational constants, we should note that, because the weights for each state included in the multistate calculation are evaluated independently, there should not be any bias in the calculation toward this degeneracy. The results for H3+, H2D+, and H3O+, described above provide illustrations where such a degeneracy does not exist. The results for H3+ demonstrated more rotational state mixing of the symmetric top basis states than one would expect to find for a symmetric top molecule, and the near degeneracy of the energies for CH5+ could reflect similar state mixing. If this were the case, then the similarity of the calculated energies could result from averaging the energies of multiple nearly degenerate eigenstates. To check if this is a valid concern, we evaluated the number of symmetric top basis states that contribute to each of the calculated states in the multistate DMC calculation. Following an earlier study by Burleigh et al.16 in which they investigated the role of Coriolis coupling in rotationally excited states of formaldehyde, this is achieved by evaluating the inverse participation ratios.87 Following that study, we define a probability for the nth walker to be in the |J,K⟩p state by pn(,i)J , K , p = |Cn(i,)J , K , p|2

ρn(i) =

1 ∑K , p (pn(,i)J , K , p )2

(20)

ρ(i) n

Operationally, provides a measure of the number of rotational basis functions that contribute to a given eigenstate. The reported results reflect the average of ρ(i) over the n probability amplitude for the ith rotational state vector using descendent weighting. The fact that values reported in Table 5 and in the Supporting Information often deviate significantly from 1 is an indication that these states are highly mixed. This large deviation of many of the ⟨ρ(i)⟩ from unity supports our concern that the similarity of the rotational energies for different states with the same value of J may reflect the fact that the states that are calculated are averages of several rotational eigenstates with different energies. To test this, we repeated the evaluation of the rotation-vibration energies, but now constraining all but one of the coefficients in the expansion of the rotational state vector in eq 17 to be zero. The results of this analysis are provided in the Supporting Information. As is seen, the range of the calculated energies for a specific value of J is smaller than when the coefficients are not constrained. This provides us with confidence in the observation that in the absence of vibrational excitation the rotational energies of CH5+ are expected to follow a near-spherical top progression. Next, we turn our attention to the vibrational contribution to the wave function. We will focus on the four types of distributions we investigated in our previous studies of CH5+.35 The first two are the distributions of the CH and HH distances. These are plotted for J = 0 and one of the J = 10 states in Figure 2a and b. As is seen, the two distributions are virtually

Figure 2. Projections of the probability amplitude for CH5+ with J = 0 (red dashed lines) and J = 10 (black solid lines) onto (a) the CH distances, (b) the HH distances, (c) q, and (d) ϕ, where q and ϕ are defined in the text.

indistinguishable. The remaining two distributions require us to assign each of the walkers to one of the 120 minima in the potential. To do this, HA and HB are identified as the two hydrogen atoms that share the shortest HH bond. Of these, the hydrogen engaged in the next shortest HH bond is labeled as HB and shares that bond with HC. HD and HE are then labeled such that

(19)

rC⃗ −H2,COM·( rC⃗ −HD × rC⃗ −HE) > 0

and the inverse participation ratio is defined as L

(21)



Article

appear to be relatively insensitive to the reference structure used to evaluate the body-fixed axis system through the Eckart conditions.38,48 Based on the success of the multistate rotational DMC approach for calculations on these three well-understood molecules, we turned our attention to CH5+ where several rotationally resolved vibrational spectra have been obtained, but, as yet have not been assigned.66,69 The results of these calculations of the v = 0 rotational energies show a progression that is consistent with a spherical top molecule. Further, the deviations between the calculated rotational energies and those obtained using a rigid-rotor Hamiltonian with vibrationally averaged rotational constants are small. Given the large amplitude vibrations in this molecule and its known highly fluxional behavior, this is somewhat surprising, but not inconsistent with earlier studies.35,69 We expect that the rotational structure will become more complicated as vibrational excitation is introduced, and we are working to extend the approach to include vibrational excitation of both low- and high-frequency modes by combining this treatment with the methodologies we previously used to study vibrationally excited states,47 and to extend the rotational Hamiltonian to include the effects of Coriolis couplings, which are expected to become more important in states with both rotational and vibrational excitation. We also plan to extend these studies to partially deuterated variants of CH5+ and H5+, which will exhibit pure rotational spectra.36,41

With this notation in place, we define the coordinates that separate the minima on the potential, q and ϕ. Further descriptions of these coordinates are provided in our earlier studies of CH5+ and have not been provided here since the focus of the present study is in changes in the projected probability amplitude onto these coordinates.47 Specifically, q = rHB − HC − rHA − HD (22) The ϕ coordinate is constructed by first defining the vector that connects the carbon atom to the center of mass of the H3+ group containing HC, HD, and HE. This vector is defined to lie along the z-axis. The x-axis is defined such that rA⃗ B lies in the xz plane. The value of ϕ is given by the angle between the xz plane and the plane that contains HC, HD, or HE and the z-axis. Note that in practice the distributions associated with all three of these definitions of ϕ are averaged together. Projections of the probability amplitude onto q and ϕ are provided in Figure 2c and d. Once again the plots for the J = 0 state and one of the J = 10 states are indistinguishable. Finally, we evaluate the extent to which the various wave functions that are obtained from the multistate calculations are sampling the 120 minima on the potential surface. This is achieved by calculating the inverse participation ratio, using the approach described in eq 20, defining the probabilities as the normalized probability amplitudes in each of the 120 minima. For all of the states the value of ρ is basically constant, indicating, once again only small changes in the vibrational wave function for the v = 0 state for J ≤ 10. The largest changes are found for those states with the largest values of KC, leading us to conclude that there might be a small increase in the localization of these states, but the effect is small. Overall, for v = 0, the rotation−vibration mixing in CH5+ is modest and appears to be much smaller than what we found in H3+, H2D+, or H3O+. This is at first surprising given the lowenergy barriers and the highly fluxional behavior of this molecular ion. On the other hand, the observation is consistent with the fact that although the molecule is highly fluxional along the isomerization coordinates the distributions of the HH and CH distances, shown in Figure 2a and b, are nearly identical to those which would be obtained if a harmonic ground state probability distribution were used to generate them.35 Likewise, we demonstrated that we could reproduce the low-resolution spectrum of CH5+ by convoluting the scaled harmonic spectra obtained at the three low-lying stationary points on the potential.69

■

ASSOCIATED CONTENT

S Supporting Information *

Projections of the DMC wave function onto rotational coordinates; complete lists of calculated energies for J ≤ 6 and J = 10 for H3+, H2D+, H3O+, and CH5+ using the full multistate approach; selected calculated energies for J ≤ 10 for H3+ and CH5+ using symmetric top rotational wave functions; selected energies for H3+ calculated in a |J,KA⟩p symmetric top basis. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address

■

A.S.P.: Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104, United States.

CONCLUSION In this study, we have developed a new DMC methodology capable of simultaneously describing multiple rotationally excited states of highly fluxional molecules within a single calculation. This offers enormous computational savings over our previous fixed-node DMC calculations where each state had to be treated independently. We tested the method through calculations of the energies and wave functions for H3+, H2D+, and H3O+, where we could compare the results to reported variational or experimental energies. We find that the approach yields energies that agree with benchmark values to within the statistical uncertainties for low to moderate values of J. An unexpected observation was that better agreement with benchmark values was found for H2D+, which is the only asymmetric top molecule that was investigated in the present study. Also, in contrast to our previously reported fixed-node approach for evaluating rotationally excited states, the results

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Support through grants from the Chemistry Division of the National Science Foundation (CHE-1213347) is gratefully acknowledged. A.S.P. and J.E.F. acknowledge support from the Ohio State University through a Distinguished University Fellowship and a Mayer’s Summer Research Scholarship, respectively. This work was supported in part by an allocation of computing time from the Ohio Supercomputer Center.

■

REFERENCES

(1) Roscioli, J. R.; McCunn, L. R.; Johnson, M. A. Quantum Structure of the Intermolecular Proton Bond. Science 2007, 316, 249− 254.

M



Article

(2) Relph, R. A.; Guasco, T. L.; Elliott, B. M.; Kamrath, M. Z.; McCoy, A. B.; Steele, R. P.; Schofield, D. P.; Jordan, K. D.; Viggiano, A. A.; Ferguson, E. E.; et al. How the Shape of an H-Bonded Network Controls Proton-Coupled Water Activation in HONO Formation. Science 2010, 327, 308−312. (3) Brathwaite, A. D.; Reed, Z. D.; Duncan, M. A. Infrared Photodissociation Spectroscopy of Copper Carbonyls. J. Phys. Chem. A 2011, 115, 10461−10469. (4) Garand, E.; Kamrath, M. Z.; Jordan, P. A.; Wolk, A. B.; McCoy, A. B.; Miller, S. J.; Johnson, M. A. Direct Determination of Docking Motifs in Multiple-Contact, Non-Covalent Host-Guest Linkages Through Mass-Selective Vibrational Spectroscopy of Cold Ionic Clusters. Science 2012, 395, 694−698. (5) Keske, J.; McWhorter, D. A.; Pate, B. H. Molecular Rotation in the Presence of Intramolecular Vibrational Energy Redistribution. Int. Rev. Phys. Chem. 2000, 19, 363−407. (6) Gordy, W.; Cook, R. L. Microwave Molecular Spectra; John Wiley: New York, 1984. (7) Mollner, A. K.; Valluvadasan, S.; Feng, L.; Sprague, M. K.; Okumura, M.; Milligan, D. B.; Bloss, W. J.; Sander, S. P.; Martien, P. T.; Harley, R. A.; et al. Rate of Gas Phase Association of Hydroxyl Radical and Nitrogen Dioxide. Science 2010, 330, 646−649. (8) Murray, C.; Derro, E. L.; Sechler, T. D.; Lester, M. I. Weakly Bound Molecules in the Atmosphere: A Case Study of HOOO. Acc. Chem. Res. 2009, 42, 419−427. (9) Sharp, E. N.; Rupper, P.; Miller, T. A. The Structure and Spectra of Organic Peroxy Radicals. Phys. Chem. Chem. Phys. 2008, 10, 3955− 3981. (10) Puzzarini, C.; Stanton, J. F.; Gauss, J. Quantum-Chemical Calculation of Spectroscopic Parameters for Rotational Spectroscopy. Int. Rev. Phys. Chem. 2010, 29, 273−367. (11) Mills, I. A. In Molecular Spectroscopy: Modern Research; Rao, K. N., Mathews, C. W., Eds.; Academic Press: New York, 1972; Vol. 3; pp 115−150. (12) Barone, V. Anharmonic Vibrational Properties by a Fully Automated Second-Order Perturbative Approach. J. Chem. Phys. 2005, 122, 014108/1−10. (13) van Vleck, J. H. The Rotational Energy of Polyatomic Molecules. Phys. Rev. 1935, 47, 487−494. (14) Sibert, E. L., III. Theoretical Studies of Vibrationally Excited Polyatomic Molecules Using Canonical Van Vleck Perturbation Theory. J. Chem. Phys. 1988, 88, 4378−4390. (15) Hougen, J. T.; Kleiner, I.; Godefroid, M. Selection Rules and Intensity Calculations for a C2 Asymmetric Top Molecule Containing a Methyl Group Internal Rotor. J. Mol. Spectrosc. 1994, 163, 559−586. (16) McCoy, A. B.; Burleigh, D. C.; Sibert, E. L. Rotation-Vibration Interactions in Highly Excited States of SO2 and H2CO. J. Chem. Phys. 1991, 95, 7449−65. (17) Luckhaus, D. 6D Vibrational Quantum Dynamics: Generalized Coordinate Discrete Variable Representation and (A)diabatic Contraction. J. Chem. Phys. 2000, 113, 1329−1348. (18) Wang, X.-G.; Carrington, T., Jr. Contracted Basis Lanczos Methods for Computing Numerically Exact Rovibrational Levels of Methane. J. Chem. Phys. 2004, 121, 2937−2954. (19) Bowman, J. M.; Carrington, T.; Meyer, H.-D. Variational Quantum Approaches for Computing Vibrational Energies of Polyatomic Molecules. Mol. Phys. 2008, 106, 2145−2182. (20) Polyansky, O. L.; Tennyson, J. Ab initio Calculation of the Rotation-Vibration Energy Levels of H3+ and its Isotopomers to Spectroscopic Accuracy. J. Chem. Phys. 1999, 110, 5056−5064. (21) Ramanlal, J.; Tennyson, J. Deuterated hydrogen chemistry: partition functions, equilibrium constants, and transition intensities for the H3+ system. Mon. Not. R. Astron. Soc. 2004, 354, 161−168. (22) Sochi, T.; Tennyson, J. A computed line list for the H2D+ molecular ion. Mon. Not. R. Astron. Soc. 2010, 405, 2345−2350. (23) Huang, X.; Schwenke, D. W.; Lee, T. J. Rovibrational spectra of ammonia. I. Unprecedented accuracy of a potential energy surface used with nonadiabatic corrections. J. Chem. Phys. 2011, 134, 044320/ 1−15.

(24) Wang, X.-G.; Carrington, T., Jr. Vibrational energy levels of CH5+. J. Chem. Phys. 2008, 129, 234102/1−13. (25) Anderson, J. B. A Random-Walk Simulation of the Schrödinger Equation: H3+. J. Chem. Phys. 1975, 63, 1499−1503. (26) Anderson, J. B. Quantum Chemistry by Random Walk. H 2P, H3+ D3h1A′1, H2 3∑u+, H4 1∑g+, Be 1S. J. Chem. Phys. 1976, 65, 4121− 4127. (27) Metropolis, N.; Ulam, S. The Monte Carlo Method. J. Am. Stat. Assoc. 1949, 44, 335−341. (28) Austin, B. M.; Zubarev, D. Y.; Lester, W. A., Jr. Quantum Monte Carlo and Related Approaches. Chem. Rev. 2012, 112, 263−288. (29) Needs, R. J.; Towler, M. D.; Drummond, N. D.; López Ríos, P. Continuum Variational and Diffusion Quantum Monte Carlo Calculations. J. Phys.: Condens. Matter 2010, 22, 023201/1−15. (30) Xu, J.; Jordan, K. D. Application of the Diffusion Monte Carlo Method to the Binding Energy of Excess Electrons to Water Clusters. J. Phys. Chem. A 2009, 114, 1364−1366. (31) Suhm, M. A.; Watts, R. O. Quantum Monte Carlo Studies of Vibrational States in Molecules and Clusters. Phys. Rep. 1991, 204, 293−329. (32) Sandler, P.; Buch, V.; Sadlej, J. Ground and Excited States of the Complex of CO with Water: A Diffusion Monte Carlo study. J. Chem. Phys. 1996, 105, 10387−10397. (33) Thompson, K. C.; Crittenden, D. L.; Jordan, M. J. T. CH5+: Chemistry’s Chameleon Unmasked. J. Am. Chem. Soc. 2005, 127, 4954−4958. (34) Jin, Z.; Braams, B. J.; Bowman, J. M. An ab initio Based Global Potential Energy Surface Describing CH5+ → CH3+ + H2. J. Phys. Chem. A 2006, 110, 1569−1574. (35) Brown, A.; McCoy, A. B.; Braams, B. J.; Jin, Z.; Bowman, J. M. Quantum and Classical Studies of Vibrational Motion of CH5+ on a Global Potential Energy Surface Obtained from a Novel ab intio Direct Dynamics Approach. J. Chem. Phys. 2004, 121, 4105−4116. (36) Johnson, L. M.; McCoy, A. B. Evolution of Structure in CH5+ and its Deuterated Analogs. J. Phys. Chem. A 2006, 110, 8213−8220. (37) McCoy, A. B. Diffusion Monte Carlo for Studying Weakly Bound Complexes and Fluxional Molecules. Int. Rev. Phys. Chem. 2006, 25, 77−108. (38) McCoy, A. B.; Hinkle, C. E.; Petit, A. S. In Recent Advances in Quantum Monte Carlo Methods; Tanaka, S., Rothstein, S. M., Lester, W. A., Jr., Eds.; ACS Symposium Series 1094; American Chemical Society: Washington, DC, 2012; Vol. 4; pp 145−155. (39) Lee, H.-S.; McCoy, A. B. Quantum Monte Carlo Studies of the Structure and Spectroscopy of the NenOH (Ã 2Σ+,n = 1−4) van der Waals Complexes. J. Chem. Phys. 2001, 114, 10278−10287. (40) Hammer, N. I.; Diken, E. G.; Roscioli, J. R.; Johnson, M. A.; Myshakin, E. M.; Jordan, K. D.; McCoy, A. B.; Huang, X.; Bowman, J. M.; Carter, S. The Vibrational Predissociation Spectra of the H5O2+· RGn (RG = Ar, Ne) Clusters: Correlation of the Solvent Perturbations in the Free OH and Shared Proton Transitions of the Zundel Ion. J. Chem. Phys. 2005, 122, 244301/1−10. (41) Lin, Z.; McCoy, A. B. Investigation of the Structure and Spectroscopy of H5+ Using Diffusion Monte Carlo. J. Phys. Chem. A 2013, 10.1021/jp4014652. (42) Blume, D.; Lewerenz, M.; Whaley, K. B. Quantum Monte Carlo Methods for Rovibrational States of Molecular Systems. J. Chem. Phys. 1997, 107, 9067−78. (43) Moroni, S.; Blinov, N.; Roy, P.-N. Quantum Monte Carlo Study of Helium Clusters Doped with Nitrous Oxide: Quantum Solvation and Rotational Dynamics. J. Chem. Phys. 2004, 121, 3577−3581. (44) Suárez, A. G.; Ramilowski, J. A.; Benito, R. M.; Farrelly, D. Renormalization of the Rotational Constants of an Ammonia Molecule Seeded into a 4He Droplet. Chem. Phys. Lett. 2011, 502, 14−22. (45) Ceperley, D. M.; Bernu, B. The Calculation of Excited State Properties with Quantum Monte Carlo. J. Chem. Phys. 1988, 89, 6316−6328. (46) Bernu, B.; Ceperley, D. M.; Lester, W. A., Jr. The Calculation of Excited States with Quantum Monte Carlo. II. Vibrational Excited States. J. Chem. Phys. 1990, 93, 552−561. N



Article

Breaking Explains Infrared Spectra of CH5+ Isotopologues. Nat. Chem. 2010, 2, 298−302. (69) Huang, X.; McCoy, A. B.; Bowman, J. M.; Johnson, L. M.; Savage, C.; Dong, F.; Nesbitt, D. J. Quantum deconstruction of the infrared spectrum of CH5+. Science 2006, 311, 60−63. (70) Schlemmer, S.; Asvany, O.; Brünken, S. Low Temperature Trapping: From Reactions to Spectroscopy. 68th International Symposium on Molecular Spectroscopy, Columbus, OH, June 17−21, 2013; talk MA01. (71) Hodges, J. N.; Perry, A.; McCall, B. J. High Precision Spectroscopy of CH5+ Using NICE-OHVMS. 68th International Symposium on Molecular Spectroscopy, Columbus, OH, June 17−21, 2013; talk MG07. (72) Marx, D.; Parrinello, M. CH5+: The Cheshire Cat Smiles. Science 1999, 284, 59−61. (73) Kumar, P. P.; Marx, D. Understanding Hydrogen Scrambling and Infrared Spectrum of Bare CH5+ Based on ab initio Simulations. Phys. Chem. Chem. Phys. 2006, 8, 573−586. (74) East, A. L. L.; Kolbuszewski, M.; Bunker, P. R. Ab initio Calculation of the Rotational Spectrum of CH5+ and CD5+. J. Phys. Chem. A 1997, 101, 6746−6752. (75) Bunker, P. R.; Ostojic, B.; Yurchenko, S. A Theoretical Study of the Millimeterwave Spectrum of CH5+. J. Mol. Struct. 2004, 695, 253− 261. (76) Feit, M. D.; Fleck, J. A. Solution of the Schrödinger Equation by a Spectral Method II: Vibrational Energy Levels of Triatomic Molecules. J. Chem. Phys. 1983, 78, 301−308. (77) Eckart, C. Rotating Axes and Polyatomic Molecules. Phys. Rev. 1935, 47, 552−558. (78) Louck, J. D.; Galbraith, H. W. Eckart Vectors, Eckart Frames, and Polyatomic Molecules. Rev. Mod. Phys. 1976, 48, 69−106. (79) Ernesti, A.; Hutson, J. M. On the Rotational Constants of Floppy Molecules. Chem. Phys. Lett. 1994, 222, 257−262. (80) Rice, J. R. Experiments of Gram-Schmidt Orthogonalization. Math. Comput. 1966, 20, 325−328. (81) Huber, D. Energies of Vibrating and Rotating Molecules by Ladder Operators. Int. J. Quantum Chem. 1985, 28, 245−267. (82) Dinelli, B. M.; Neale, L.; Polyansky, O. L.; Tennyson, J. New Assignments for the Infrared Spectrum of H3+. J. Mol. Spectrosc. 1997, 181, 142−150. (83) Dinelli, B. M.; Polyansky, O. L.; Tennyson, J. Spectroscopically Determined Born-Oppenheimer and Adiabatic Durfaces for H3+, H2D+, HD2+, and D3+. J. Chem. Phys. 1995, 103, 10433−10438. (84) McNab, I. R. In Advances in Chemical Physics; Prigogine, I., Rice, S. A., Eds.; John Wiley and Sons, Inc.: New York, 1995; Vol. 89; Chapter: The spectroscopy of H3+, pp 1−87. (85) Lindsay, C. M.; McCall, B. J. Comprehensive Evaluation and Compilation of H3+ Spectroscopy. J. Mol. Spectrosc. 2001, 210, 60−83. (86) Petit, A. S.; McCoy, A. B. Correction to “Diffusion Monte Carlo Approaches for Evaluating Rotationally Excited States of Symmetric Top Molecules: Application to H3O+ and D3O+”. J. Phys. Chem. A 2011, 115, 9325−9327. (87) Bell, R. J.; Dean, P. Atomic Vibrations in Vitreous Silica. Discuss. Faraday Soc. 1970, 50, 55−61.

(47) Hinkle, C. E.; McCoy, A. B. Characterizing Excited States of CH5+ with Diffusion Monte Carlo. J. Phys. Chem. A 2008, 112, 2058− 2064. (48) Petit, A. S.; McCoy, A. B. Diffusion Monte Carlo approaches for Evaluating Rotationally Excited States of Symmetric Top Molecules: Application to H3O+ and D3O+. J. Phys. Chem. A 2009, 113, 12706− 12714. (49) Petit, A. S.; Wellen, B. A.; McCoy, A. B. Unraveling rotationvibration mixing in highly fluxional molecules using Diffusion Monte Carlo: Application to H3+ and H3O+. J. Chem. Phys. 2012, 136, 074101/1−12. (50) Petit, A. S.; Wellen, B. A.; McCoy, A. B. Using Fixed-Node Diffusion Monte Carlo to Investigate the Effects of Rotation-Vibration Coupling in Highly Fluxional Asymmetric Top Molecules: Application to H2D+. J. Chem. Phys. 2013, 138, 034105/1−11. (51) Hinkle, C. E.; Petit, A. S.; McCoy, A. B. Diffusion Monte Carlo studies of Low Energy Ro-vibrational States of CH5+ and its Deuterated Isotopologues. J. Mol. Spectrosc. 2011, 268, 189−198. (52) Gregory, J. K.; Clary, D. C. Calculations of the Tunneling Splitting in Water Dimer and Trimer Using Diffusion Monte Carlo. J. Chem. Phys. 1995, 102, 7817−7829. (53) Lewerenz, M. Quantum Monte Carlo Calculation of Argon-HF Clusters: Nonadditive Forces, Isomerization, and HF Frequency Shifts. J. Chem. Phys. 1996, 104, 1028−1039. (54) Blume, D.; Lewerenz, M.; Huisken, F.; Kaloudis, M. Vibrational Frequncy Shift of HF in Helium Clusters: Quantum Simulation and Experiment. J. Chem. Phys. 1996, 105, 8666−8683. (55) Sun, Z.; Lester, W. A., Jr.; Hammond, B. L. Correlated Sampling of Monte Carlo Derivatives With Iterative-Fixed Sampling. J. Chem. Phys. 1992, 97, 7585−7589. (56) Aguado, A.; Roncero, O.; Tablero, C.; Sanz, C.; Paniagua, M. Global potential energy surfaces for H3 + system. Analytical representation of the adiabatic ground-state 11A′ potential. J. Chem. Phys. 2000, 112, 1240−1254. (57) Cencek, W.; Rychlewski, J.; Jaquet, R.; Kutzelnigg, W. Submicrohartree Accuracy Potential Energy Surface for H3+ Including Adiabatic and Relativistic Effects. I. Calculation of the Potential Points. J. Chem. Phys. 1998, 108, 2831−2836. (58) Miller, S.; Tennyson, J. First Principles Calculation of the Molecular Constants of H3+, H2D+, D2H+ and D3+. J. Mol. Spectrosc. 1987, 126, 183−192. (59) Huang, X.; Carter, S.; Bowman, J. M. Ab initio Potential Energy Surface and Rovibrational Energies of H3O+ and its Isotopomers. J. Chem. Phys. 2003, 118, 5431−5442. (60) Rajamäki, T.; Miani, A.; Halonen, L. Six-Dimensional Ab initio Potential Energy Surfaces for H3O+ and NH3: Approaching the Subwavenumber Accuracy for the Inversion Splittings. J. Chem. Phys. 2003, 118, 10929−10938. (61) Tang, J.; Oka, T. Infrared Spectroscopy of H3O+: the ν1 Fundamental Band. J. Mol. Spectrosc. 1999, 196, 120−130. (62) Olah, A.; Rasul, G. From Kekulé’s Tetravalent Methane to Five-, Six-, and Seven-Coordinate Protonated Methanes. Acc. Chem. Res. 1997, 30, 245−250. (63) Maluendes, S. A.; McLean, A. D.; Herbst, E. Calculations of Ion−Molecule Deuterium Fractionation Reactions Involving HD. Astrophys. J. 1992, 397, 477−481. (64) Roberts, H.; Herbst, E.; Millar, T. J. Enhanced Deuterium Fractionation in Dense Interstellar Cores Resulting from Multiply Deuterated H3+. ApJ 2003, 591, L41−L44. (65) Herbst, E. Isotopic Fractionation by Ion−Molecule Reactions. Space Sci. Rev. 2003, 106, 293−304. (66) White, E. T.; Tang, J.; Oka, T. CH5+. The Infrared Spectrum Observed. Science 1999, 284, 135−137. (67) Asvany, O.; Kumar, P.; Redlich, B.; Hegeman, I.; Schlemmer, S.; Marx, D. Understanding the Infrared Spectrum of Bare CH5+. Science 2005, 309, 1219−1222. (68) Ivanov, S. D.; Asvany, O.; Witt, A.; Hugo, E.; Mathias, G.; Redlich, B.; Marx, D.; Schlemmer, S. Quantum-Induced Symmetry O


Simultaneous Evaluation of Multiple Rotationally ... - ACS Publications

Recommend Documents