Diffusion Monte Carlo in Internal Coordinates - The Journal of Physical

Article pubs.acs.org/JPCA

Diffusion Monte Carlo in Internal Coordinates Andrew S. Petit and Anne B. McCoy* Department of Chemistry and Biochemistry, The Ohio State University, Columbus, Ohio 43210, United States S Supporting Information *

ABSTRACT: An internal coordinate extension of diffusion Monte Carlo (DMC) is described as a first step toward a generalized reduced-dimensional DMC approach. The method places no constraints on the choice of internal coordinates other than the requirement that they all be independent. Using H+3 and its isotopologues as model systems, the methodology is shown to be capable of successfully describing the ground state properties of molecules that undergo large amplitude, zero-point vibrational motions. Combining the approach developed here with the fixed-node approximation allows vibrationally excited states to be treated. Analysis of the ground state probability distribution is shown to provide important insights into the set of internal coordinates that are less strongly coupled and therefore more suitable for use as the nodal coordinates for the fixed-node DMC calculations. In particular, the curvilinear normal mode coordinates are found to provide reasonable nodal surfaces for the fundamentals of H2D+ and D2H+ despite both molecules being highly fluxional.

■

INTRODUCTION Over the years since its introduction to the chemical physics community by Anderson,1,2 diffusion Monte Carlo (DMC) has developed into a powerful technique for both electronic structure calculations and the study of the rotational and vibrational ground and excited state properties of molecules and clusters.3−8 The methodology is based on the propagation of the time-dependent Scrhödinger equation in imaginary time using techniques originally developed for statistical simulations.9 As a result, DMC has a much more favorable scaling with system size than variational approaches.3 This is especially true when studying the nuclear dynamics of molecules and clusters that exhibit large amplitude, zero-point vibrational motions and therefore whose vibrational wave functions sample multiple minima on the potential energy surface (PES), even in the ground state.6 Such systems lack a well-defined, zerothorder description that can be used to construct an efficient direct product basis or rapidly convergent perturbative series. This characteristic makes the application of approximate methods to fluxional systems especially challenging. Despite having several advantages over variational approaches, DMC is not without challenges. As will be discussed in detail below, DMC is fundamentally a ground state method. The extension of DMC to excited states has typically been achieved using either the correlation function approach or the fixed-node approximation. The correlation function approach involves the use of a trial function constructed with the appropriate symmetries of the system of interest to obtain the excited state energies.10,11 Although this technique has been shown to yield very accurate energies, it cannot be used to obtain the excited state wave functions. The fixed-node approximation, on the other hand, provides the lowest energy state consistent with a predetermined set of nodal surfaces obtained from either model calculations or symmetry considerations.1,2 Despite introducing a more severe approx© 2013 American Chemical Society

imation than what appears in the correlation function approach, fixed-node DMC allows the excited state wave functions to be directly determined. The analysis of these wave functions can often provide important physical insights into the effects of the excitation on the system’s structure and dynamics. We have demonstrated that fixed-node DMC provides an effective approach for the study of the vibrationally and rotationally excited states of highly fluxional molecules and clusters, and we will make further use of it in this study.6−8,12−18 A more fundamental limitation of DMC follows from constraints placed on the choice of coordinates in which to perform the calculations. Specifically, in most implementations of DMC, the kinetic energy operator is expressed in the form 3N

T̂ =

∑ k=1

pk̂ 2 2mk

(1)

where the masses, mk, are all real numbers. As will be discussed in detail below, this structure allows the use of statistical simulation techniques to evaluate the imaginary time propagation of the time-dependent Schrödinger equation. When internal coordinates are used, the kinetic energy operator contains momentum coupling terms and coordinate-dependent effective masses that prevent it from being expressed in the form of eq 1. As a result, DMC simulations are typically performed in Cartesian coordinates and therefore require the use of a full-dimensional PES expressed as a function of all 3N − 6 internal coordinates. Special Issue: Joel M. Bowman Festschrift Received: December 24, 2012 Revised: February 12, 2013 Published: February 14, 2013 7009

dx.doi.org/10.1021/jp312710u | J. Phys. Chem. A 2013, 117, 7009−7018

The Journal of Physical Chemistry A

■

Article

THEORY Diffusion Monte Carlo in Cartesian Coordinates. As our implementation of the DMC methodology for the study of highly fluxional molecules and clusters has been described elsewhere, we will only provide a general overview.6,7 The method is based on the observation that when the timedependent Schödinger equation is expressed in terms of the imaginary time variable τ = it/ℏ, it assumes the same mathematical form as a diffusion equation and therefore can be solved through a Monte Carlo simulation.1,2 In order to accomplish this, the wave function of the system of interest is expanded in an ensemble of δ-functions, or walkers, each of which is localized at a specific molecular geometry in configuration space. More specifically, the kth walker in the ensemble, with the 3N Cartesian coordinates x⃗ (k) , is represented by46

Significant advances have been made in recent years in the development of full-dimensional PESs of high accuracy.19−21 Even so, the construction of such full-dimensional surfaces remains both challenging and computationally expensive. Moreover, there are situations where either a full-dimensional treatment of the nuclear dynamics is out of reach of current computational approaches or where focusing on a smaller subsystem may simplify the analysis while still providing physical insights into the larger system. Reduced-dimensional approaches have been used by us and others to study a variety of systems that undergo large amplitude motions.18,22−27 Typically such studies involve expanding the wave function in a basis set and involve fewer than 4 degrees of freedom. As we embarked on this study, we wanted to learn whether DMC could be applied in a general way to reduced-dimensional studies. Reduced-dimensional variants of DMC are not without precedent. Rigid-body DMC, originally developed by Buch, is based on an adiabatic separation between the intermolecular and intramolecular vibrational modes in a cluster.28 Specifically, the monomers are treated as rigid bodies, and the DMC simulation is performed in terms of their center of mass locations and relative orientations.28,29 This approach has been successfully used to describe low-lying ro-vibrational excited states of clusters as well as to model the quantum effects of confinement in helium nanodroplets on the apparent rotational constants and torsional splittings of the dopant.30−37 A similar approach was used by Deskevitch et al. to perform DMC in the set of bending modes of CH4 and CH+5 , with the CH bond lengths held frozen at their equilibrium values.38 In our recent application of DMC to the study of minimized energy paths, we performed the DMC simulations in Jacobi coordinates with the distance between the centers of mass of the reacting fragments constrained.39,40 However, these approaches are currently limited to very specific forms of dimensionality reduction that satisfy the constraints imposed by requiring the kinetic energy to have the structure of eq 1. A generalization that allows reduced-dimensional DMC simulations to be performed in an arbitrary set of internal coordinates has not yet been made. As a first step toward a generalized reduced-dimensional DMC methodology, we develop in this study an internal coordinate DMC approach that allows us to use any set of independent internal coordinates. In particular, the form of the kinetic energy operator used in the DMC simulations can deviate significantly from eq 1. The approach is tested on the ground and v = 1 vibrationally excited states of H+3 and its isotopologues. These molecules are known to exhibit largeamplitude, zero-point vibrational motions and therefore represent a stringent test of the method’s underlying approximations. Additionally, several full-dimensional PESs of very high accuracy have been constructed for H3+, and converged variational calculations to spectroscopic accuracy have been reported in the literature for all four isotopologues.41−45 In addition to demonstrating the effectiveness of the internal coordinate DMC approach developed in this study, we will also consider the insights into the nodal structure of vibrationally excited states that can be obtained from reduceddimensional approaches as well as from analysis of ground state DMC probability distributions. Finally, we will discuss the use of normal modes as the nodal coordinates for fixed-node DMC calculations of fluxional molecules.

⎡ (x ⃗ − x ⃗(k))2 ⎤ ⎛ 1 ⎞3N /2 ⎥ δ(x ⃗ − x ⃗(k)) = lim ⎜ ⎟ exp⎢ − ⎥⎦ ⎢⎣ β→ 0⎝ βπ ⎠ β

(2)

The evolution of the ensemble of walkers in imaginary time is determined by ̂

̂

Ψ(x ⃗ , τ + δτ ) ≈ e−(V − Eref (τ))δτ e−Tδτ Ψ(x ⃗ , τ )

(3)

where we employ the split-operator approximation to the propagator and evaluate the energy relative to an imaginary time dependent reference energy, Eref(τ).47 Note that throughout this work, the position representation will be implicitly assumed for all operators. If Cartesian coordinates are used, then the action of the kinetic energy component of the propagator on the kth walker in the ensemble is given by ⎡ 3N δτp2̂ ⎤ i ⎥ δ(x ⃗ − x ⃗(k)) exp⎢ −∑ ⎢⎣ i = 1 2mi ⎥⎦ 3N

=

∏ i=1 3N

∝

⎡ δτp2̂ ⎤ exp⎢ − i ⎥δ(xi − xi(k)) ⎢⎣ 2mi ⎥⎦ ⎡ m (x − x (k))2 ⎤ i i i ⎥ ⎢⎣ ⎥⎦ 2δτ

∏ exp⎢− i=1

(4)

where we employ the definition of δ(x⃗ − x⃗ ) given in eq 2. It is worth emphasizing that this result follows from the fact that, in Cartesian coordinates, all of the masses, {mi}, are independent of the molecular geometry, and the kinetic energy operator contains no p̂ip̂j cross terms. Instead of allowing the walker to spread, as is implied by eq 4, we preserve the locality of the basis by randomly displacing the walker in configuration space based on the Gaussian distribution given in eq 4. After ̂ this diffusion, the integer part of e−(V−Eref(τ))δτ determines the number of walkers added to the ensemble at the walker’s location in configuration space, while the remainder gives the probability that the original walker remains in the ensemble.1,2,5 Throughout the DMC simulation, the reference energy evolves according to (k)

⎛N (τ ) − Nwalker(0) ⎞ Eref (τ ) = V̅ (τ ) − α⎜ walker ⎟ Nwalker(0) ⎠ ⎝

(5)

where V̅ (τ) is the ensemble average of the potential energy at imaginary time τ, Nwalker(τ) is the instantaneous number of walkers in the ensemble at that time, and α is a simulation 7010



Article ̂ (int)

transformation that simplifies the action of e−T δτ is coordinate dependent, and valid only in the limit of short time steps. Indeed, as δτ is increased, the changes to the values of the Gmatrix elements within a single time step will eventually become too large to be neglected, and the approximation will break down. Thus, we evaluate the validity of G(Q⃗ (τ)) ≈ G(Q⃗ (τ + δτ)) by exploring the dependence of the DMC results on the size of δτ when internal coordinates are used. In summary, we propose the following procedure for evaluating the evolution of the DMC ensemble within a given time step when internal coordinates are employed. For a walker initially at Q⃗ (k)(τ), we first construct the displacement vector δQ̃ (k) by randomly sampling a uniform Gaussian distribution with a variance of δτ. This represents the walker’s ̂ (int) random displacement due to the action of e−T δτ as expressed in the instantaneous set of mass-weighted internal coordinates defined through G−1/2(Q⃗ (k)(τ)). The location of the walker at imaginary time τ + δτ is then determined using Q⃗ (k)(τ + δτ) = → ⎯ Q⃗ (k)(τ) + G1/2(Q⃗ (k)(τ))δQ͠ (k). The potential energy component of the propagator then acts as before with the exception that V̂ is replaced with V̂ (int). The Fixed-Node Approximation. The fixed-node approximation provides one approach for extending DMC to the study of excited states.2,5−7 It is based on the observation that the local behavior of an excited state wave function in the vicinity of a nodal surface is the same as that of a ground state wave function in the vicinity of an infinite potential barrier. Specifically, the wave function goes to zero with finite slope. As a result, we modify the potential energy used in eq 3 by making it infinite everywhere outside of a specified nodal region of the excited state of interest. This has the effect of removing from the DMC ensemble any walkers that diffuse out of that nodal region and allows us to recast an excited state calculation as a series of ground state simulations, one for each nodal region. The equilibrium distribution of walkers from a fixed-node DMC calculation will provide a Monte Carlo sampling of the portion of the excited state wave function within the specified nodal region. Likewise, at equilibrium, Eref(τ) will fluctuate about the energy of the excited state, to within the accuracy of the nodal surface used to define the modified potential. An important internal check on the validity of the fixed-node DMC calculations follows from the fact that Ĥ Ψ/Ψ = E must be locally satisfied everywhere in configuration space. Specifically, the fixed-node DMC energies calculated in different nodal regions of the same state must be identical if the nodes are properly placed. Finally, because of the use of a finite time step in the fixed-node DMC calculations, it is possible for a walker to effectively cross and then recross the nodal surface within a single time step. To correct for this, Anderson showed that the probability that a walker should be removed from the DMC ensemble because it has crossed and then recrossed the nodal surface within the time step from τ to τ + δτ is given by

parameter, which is typically ∼1/δτ.1 In the long imaginary time limit, the ensemble of walkers reaches equilibrium at which point it provides a Monte Carlo sampling of the ground state wave function. Further, once equilibrium is obtained, Eref(τ) randomly fluctuates about the zero-point energy of the system of interest, to within the accuracy obtainable from the PES and the use of a finite time-step of size δτ. Probability distributions and expectation values can be calculated using the descendent weighting approach.5,6 Generalizing Diffusion Monte Carlo to Internal Coordinates. Before discussing our implementation of DMC in internal coordinates, it is useful to consider the form that the kinetic energy operator assumes when expressed in terms of internal coordinates. Specifically, the J = 0 Hamiltonian operator for an N-atom molecule can be written in terms of a set of 3N − 6 independent internal coordinates, {Qi}, as 2 3N − 6

ℏ Ĥ =− 2 = T̂

∑ i,j=1

(int)

+ V̂

∂ ∂ Gi , j(Q⃗ ) + V ′(Q⃗ ) + V (Q⃗ ) ∂Q i ∂Q j (int)

(6)

where V̂ (int) = V̂ ′ + V̂ . The Wilson G-matrix, G(Q⃗ ), arises from the use of the chain rule to convert the derivatives in the kinetic energy operator from Cartesian to internal coordinates and its matrix elements, which act as effective inverse masses in eq 6, are given by48,49 3N

Gi , j(Q⃗ ) =

∑ k=1

∂Q i 1 ∂Q j ∂xk mk ∂xk

(7)

The V′(Q⃗ ) contribution to T̂ is a manifestation of the noncommutativity of position and momentum that arises when the classical expression for the vibrational kinetic energy is transformed into a quantum mechanical operator. The transformation to internal coordinates introduces both the p̂ip̂j cross terms and the coordinate dependent effective inverse masses, Gi,j(Q⃗ ), that are present in the expression for ̂ (int) T̂ (int) given in eq 6. As a result, the action of e−T δτ on a walker does not lead to the simple expressions given in eq 4 when internal coordinates are used. However, in the limit of small δτ used in the DMC simulations, G(Q⃗ (τ)) ≈ G(Q⃗ (τ + δτ)). Throughout the time step from τ to τ + δτ, we therefore approximate the G-matrix for each walker in the ensemble by its value at time τ. This effectively removes the coordinate ̂ (int) dependence from the effective inverse masses in e−T δτ, although the p̂ip̂j cross terms are still present. However, at each molecular geometry, it is possible to define a set of mass→ ⎯ weighted internal coordinates, Q͠ = G−1/2Q⃗ , in terms of which the T̂ (int) in eq 6 takes on the simple form T̂

(int)

=−

ℏ2 2

3N − 6

∑ i=1

∂2 2 ∂Q̃

⎡ 2m d(τ )d(τ + δτ ) ⎤ Precross = exp⎢ − eff ⎥ ⎣ ⎦ δτ

(8)

i

̂ (int)

Through this transformation, the action of e−T δτ on a walker again has the same structure as eq 4. Note that the essence of this approach, in which different representations are used for the evaluation of the kinetic and potential energy operators, is often used in both variational and time-dependent quantum mechanical calculations.50−52 It is worth emphasizing, though, that in the context of internal coordinate DMC, the

(9)

where d(τ) is the walker’s distance from the nodal surface at imaginary time τ, and meff is the, generally coordinatedependent, effective mass associated with the nodal coordinate.2 The central challenge in setting up a fixed-node DMC calculation is the determination of an appropriate nodal surface 7011



Article

results of five statistically independent DMC calculations were averaged for the ground state as well as for each nodal region of the excited states.16,17 The set of internal coordinates used in this study for H+3 and its isotopologues are two bond lengths, r12 and r23, as well as the bond angle between them, θ123. For the mixed isotopologues, these coordinates are defined such that atom 2 is the unique isotope. Expressions for the G-matrix elements as well as the V′ effective potential used in this study were constructed from the tabulations of Frederick and Woywod and are given in the Supporting Information.49 Finally, note that V′(Q⃗ ) diverges as cot2(θ123) when the molecule approaches linearity, and, in order to prevent numerical problems arising from this divergence, any walker with θ123 < 2.29° or θ123 > 177.62° was removed from the DMC ensemble. This constraint was consistently found to affect only a few walkers over the course of an entire DMC simulation and, as a result, does not have a significant impact on the results. The nodal coordinates, {si}, used in this study are chosen to be linear combinations of the internal coordinates and are specified in detail below. Mathematically, they are defined by the relationship Q⃗ = Ls,⃗ and the G-matrix is expressed in terms of these nodal coordinates through the transformation G(s) = L−1G(LT)−1.48 The effective mass used in eq 9 for the recrossing correction when a node is placed in sk is given by

for the excited state of interest. For vibrational problems, we often assume that the nodal surface can be expressed as g(Q⃗ ) = ηnode, where g(Q⃗ ) is the nodal coordinate appropriate for the excited state of interest, and ηnode is a constant. Symmetry can sometimes provide enough restrictions to unambiguously define the nodes as in, for example, the asymmetric stretch of H2D+. When this is not the case, adiabatic diffusion Monte Carlo (ADMC) can be used to define an optimal value for ηnode for a chosen g(Q⃗ ).29 Specifically, fixed-node DMC calculations are performed in which ηnode is increased each time step by a constant amount that is small enough for Eref(τ) to adiabatically follow the changing location of the node. In other words, the shift in the location of the node within a single time step causes a shift in Eref(τ) that is much smaller than its statistical fluctuations. The reference energy in both nodal regions is monitored throughout these calculations, and the optimal value of ηnode is identified as the value at which Eref(τ) becomes equal on both sides of the node. It is worth emphasizing that even when an ADMC optimized location for the node is used, the accuracy of a fixed node calculation is only as good as the chosen functional form for g(Q⃗ ).6,7,13−16 In cases where comparison to other types of calculations could be made, we have found that symmetry-adapted linear combinations of valence coordinates often provide good definitions of the nodal coordinates, even in systems that exhibit large amplitude, zeropoint vibrational motions.12,15,29 Finally, note that the determination of the optimal nodal surface for a fixed-node DMC calculation is independent of the set of coordinates used to perform the DMC simulation. Moreover, as demonstrated below, a given choice for g(Q⃗ ) and ηnode will produce numerically identical results for fixed-node DMC simulations performed in Cartesian or internal coordinates. An additional challenge arises when the results of separate fixed-node DMC calculations performed in different nodal regions are spliced together to yield an overall description of the excited state. The DMC probability distribution in each nodal region is independently normalized and, in the absence of symmetry, the relative weights in each of the nodal regions cannot be directly determined. In this study, we approximate these relative weights by requiring the second derivative of the probability amplitude to be continuous across the node.13 In practice, this is accomplished by projecting the DMC probability amplitudes in each of the nodal regions onto the nodal coordinate, g(Q⃗ ), and comparing local fits of these distributions in the vicinity of the node to quadratic polynomials in g(Q⃗ ).

1 = meff

Gk(s, k) (Q⃗ (τ ))Gk(s, k) (Q⃗ (τ + δτ ))

(10)

where we have taken the geometric mean of the appropriate diagonal matrix element of G(s) at imaginary times τ and τ + δτ. For excited states where the position of the node cannot be determined by symmetry, ADMC calculations performed in internal coordinates are used to define the optimal location of the node. The same value of ηnode is used for both the internal and Cartesian fixed-node DMC calculations. Details of the ADMC calculations performed in this study can be found in the Supporting Information. Two-dimensional discrete variable representation (DVR) calculations performed in the sinc-DVR basis of Colbert and Miller were used in this study to gain additional insights into the behavior of the fixed-node DMC simulations.53 In order to simplify the evaluation of the kinetic energy matrix elements, each term in the kinetic energy operator was expanded into Hermitian form according to Luckhaus,54

■

pĵ Gĵ , k pk̂ =

COMPUTATIONAL DETAILS This study was performed using the full-dimensional PES of H+3 developed by Aguado and co-workers.41 Unless otherwise specified, the DMC calculations were performed with δτ = 0.5 au. For all of the calculations, α = (2δτ)−1. The ground state calculations began with all of the walkers at the equilibrium geometry of H+3 while the excited state simulations were initialized with ensembles of walkers from well-equilibrated H+3 ground state calculations. For both the ground and excited state DMC calculations, Nwalker(0) = 20 000, and the walkers were allowed to equilibrate for 45 000 au before collecting information. For each DMC simulation, descendent weighting was used to perform 15 statistically independent calculations of |Ψ|2, each of which lasted 250 au and was separated from the next by 2500 au. The reference energy was averaged over the entire 37 000 au descendent weighting period. Finally, the

⎡ ∂ 2Gj , k ⎤ 1⎢ ⎥ pĵ pk̂ Gĵ , k + Gĵ , k pĵ pk̂ + ℏ2 ∂Q j∂Q k ⎥⎦ 2 ⎢⎣

(11)

The matrix elements of this operator in the sinc-DVR basis, {|m,n⟩}, were then evaluated such that the momentum operators act in the opposite direction to Ĝ j,k, that is, ⟨m , n|pĵ Gĵ , k pk̂ |m′, n′⟩ =

7012

1 Gj , k (Q 1(m′) , Q 2(n′))⟨m , n|pĵ pk̂ |m′, n′⟩ 2 1 + Gj , k (Q 1(m) , Q 2(n))⟨m , n|pĵ pk̂ |m′, n′⟩ 2 2 ℏ2 ∂ Gj , k (Q (m) , Q 2(n))δm , m′δn , n′ + 2 ∂Q j∂Q k 1

(12)



Article

found to be in better agreement with the previously reported results. However, all of the differences between the DMC and variational energies range from 1 to 6 cm−1 and are within or only slightly outside of the statistical uncertainties of the DMC calculations. Besides the DMC energies, we also investigated the projections of the DMC probability distributions onto the vibrational coordinates and found no dependence on whether the DMC calculations were performed in Cartesian or internal coordinates. These results indicate that the internal coordinate DMC approach developed in this study provides an effective approach for investigating the vibrational ground state properties of molecules and clusters that undergo largeamplitude, zero-point vibrational motions. We now turn to the dependence of the DMC results on the choice of internal coordinates in which the DMC simulations are performed. Throughout this work, we have used {r12,r23,θ123}, which represents a valid set of internal coordinates as each member is defined over a range that is independent of the values of the other coordinates in the set. In order to determine the importance of using a valid set of internal coordinates, we performed a series of internal coordinate DMC calculations in the set of three bond lengths, {r12,r23,r13}. These three bond lengths are not all independent due to the geometric constraints of the triangle inequality, and, as a result, {r12,r23,r13} does not constitute a valid set of internal coordinates. For H3+, performing the internal coordinate DMC calculations in the three bond lengths instead of {r12,r23,θ123} results in a calculated zero-point energy that is over 150 cm−1 larger than the variational value. Finally, recall that the effectiveness of the internal coordinate DMC approach presented here depends on the G-matrix varying slowly enough over the range of displacements that a walker undergoes within a single time step that G(Q⃗ (τ)) provides a good approximation to the instantaneous G-matrix throughout the propagation from τ to τ + δτ. This approximation is tested by evaluating the difference between the energies from internal coordinate DMC calculations performed with δτ = 0.5 au and δτ = 5 au, that is, E(δτ=0.5) − DMC E(δτ=5) DMC . The results of this analysis are presented in Table 1, and no statistically significant time step dependence is observed in the DMC results of any of the species considered in this study. Moreover, no consistent trend is observed in even the sign of E(δτ=0.5) − E(δτ=5) DMC DMC . The v = 1 Vibrationally Excited States of H+3 and D+3 . In this study, four coordinates are used to define the nodal surfaces for the v = 1 states of H+3 and D+3 . Three of these are the linear combinations of the bond lengths which transform according to the A′ or E′ irreducible representations of the D3h symmetry group, specifically

In this way, all terms involving the explicit action of the momentum operators on Gj,k are captured by the third term in eq 12. The remaining matrix elements of p̂ip̂j in eq 12 were then evaluated using the expressions developed by Colbert and Miller when j = k and by Luckhaus and McCoy when j ≠ k.53−55 Finally, details of the grids used to perform the 2D-DVR calculations are provided in the Supporting Information.

■

RESULTS AND DISCUSSION In the analysis reported here, we aim to address three key questions regarding the effectiveness of the internal coordinate DMC approach developed in this study as well as the general application of the fixed-node approximation to the study of the vibrationally excited states of fluxional molecules. First, how accurate are the results of the internal coordinate DMC simulations, and how do they compare to the results of DMC simulations performed in Cartesian coordinates? In particular, are the approximations underlying our treatment of the kinetic energy component of the propagator valid for the range of values for δτ typically used in our DMC calculations? Second, can analysis of the ground state DMC probability distribution provide insights into the set of internal coordinates that are less strongly coupled and therefore more suitable for use as nodal coordinates for fixed-node DMC calculations? Further, when deficiencies in the fixed-node approximation are observed, can simple reduced dimensional calculations point to the properties of the excited state wave functions and nodal surfaces that are to blame? Finally, can curvilinear normal modes provide viable definitions for the nodal coordinates of fluxional molecules, especially when the symmetry adapted linear combinations of the valence coordinates are found to be strongly coupled at the harmonic level? The Ground State of H+3 and Its Isotopologues. Table 1 compares the zero-point energies of H+3 and its isotopologues obtained from DMC calculations performed using Cartesian and internal coordinates along with the results of variational calculations reported by Tennyson and co-workers and Aguado et al.41,44 For the mixed isotopologues, the DMC energies exhibit no statistically significant dependence on the choice of coordinate system. The DMC energies of H+3 and D+3 are found to be smaller when the DMC calculations were performed in internal coordinates. Furthermore, the DMC energies obtained from the calculations performed in Cartesian coordinates are Table 1. Comparison between the Zero-Point Energies, Given in cm−1, of H+3 and Its Isotopologues Obtained from DMC Calculations Performed Using Cartesian and Internal Coordinates as Well as from Variational Calculationsa species H+3 H2D+ D2H+ D+3

Evariationalb 4362.1 3979.9 3563.1 3113.7

EDMC Cartesianc,d 4361.7 3978.6 3561.2 3113.1

± ± ± ±

1.7 1.2 1.3 1.7

EDMC internalc,e 4356.3 3976.4 3560.0 3109.1

± ± ± ±

3.7 3.7 2.0 3.0

c E(δτ=0.5) − E(δτ=5) DMC DMC

s

−3.0 ± 4.4 0.6 ± 4.4 −2.2 ± 4.6 −1.9 ± 3.8

Also shown is the time step dependence, given in cm−1, of the internal coordinate DMC calculations. bThe variational benchmark values for H+3 and its isotopologues are from Aguado et al. and Miller and Tennyson, respectively.41,44 cThe reported error bars are 99% confidence intervals based on five statistically independent simulations. d The results are from DMC calculations performed with δτ = 1 au and α = 0.5 Hartree. eThe results are from DMC calculations performed with δτ = 0.5 au and α = 1 Hartree. a

=

1 (r12 + r23 + r13) 3

a1 =

1 (r12 − r23) 2

a2 =

1 (r12 + r23 − 2r13) 6

(13)

For reasons discussed below, we will also consider the excited state associated with a node placed in θ123. Note that in the limit of small amplitude vibrations, a2 and θ123 describe the same motion. Indeed, as shown in Figure 1, placing a node in the a2 coordinate results in the DMC probability distribution 7013



Article

Table 2. The Energies, Given in cm−1, of the v = 1 Vibrationally Excited States of H+3 and D+3 Obtained from Internal Coordinate Fixed-Node DMC Calculations Compared to the Results of Variational Calculations system H+3

D+3

nodal coordinate

optimized ηnodea

Evariationalb

EDMC − Evariationalc

s a1 a2 θ123 s a1 a2 θ123

1.5866 ± 0.0001 Å 0 0.0075 ± 0.0001 Å 60.044 ± 0.005° 1.5647 ± 0.0001 Å 0 0.0051 ± 0.0001 Å 60.028 ± 0.012°

3179.1 2521.1 2521.1 2521.1 2301.2 1835.1 1835.1 1835.1

27.3 ± 4.0 0.1 ± 4.0 21.5 ± 4.2 3.0 ± 4.0 14.8 ± 3.8 0.2 ± 3.4 9.6 ± 4.1 1.3 ± 3.7

The optimized node locations, ηnode, are the average of the results of five ADMC calculations performed in internal coordinates, and the error bars represent one standard deviation. bThe variational benchmark values for H+3 and D+3 are from Aguado et al. and Miller et al., respectively, and are reported relative to the zero-point energy.41,44 cFive statistically independent DMC simulations were performed in each nodal region, and the reported error bars are 99% confidence intervals. a

lengths, and differences were observed between the fixed-node DMC energies of states, which should have been degenerate. As shown in Table 2, fixed-node DMC calculations performed with a node placed in θ123 yield an energy that agrees with both the literature value for the energy of the fundamental in the E′ mode and the fixed-node DMC energy of the state with a node placed in a1 to within the statistical uncertainty of the DMC calculations. Thus, even though θ123 and a2 describe the same vibrational motion in the limit of small amplitude vibrations, the fixed-node approximation is more accurate when the node is placed in θ123. In order to understand the cause of this, we constructed two-dimensional cuts of the PES along a2 and s as well as along θ123 and s, both of which are provided in the Supporting Information. We find that the potential cut in terms of a2 and s possesses a somewhat larger curvature than when the potential is expressed in terms of θ123 and s. This indicates that the potential energy coupling between a2 and s is greater than that between θ123 and s and as a result, θ123 is a better nodal coordinate for the fixed-node DMC calculations than a2. Unlike the asymmetric modes a1 and θ123, the fixed-node DMC energy of the v = 1 symmetric stretch state exhibits a statistically significant deviation from the benchmark value. Specifically, the DMC energies are too high by 27.3 ± 4.0 cm−1 for H+3 and 14.8 ± 3.8 cm−1 for D+3 , errors that amount to less than 1% of the total excitation energy. In our recent study of H+5 , we found that reduced-dimensional approaches can aid in understanding the properties of excited state wave functions and their nodal surfaces.18 Here, we perform a 2D-DVR calculation of H+3 in the symmetric and asymmetric stretch coordinates (r12 ± r23)2−1/2 at θ123 = 60°, the details of which are given in the Supporting Information. As shown in Figure 2, the wave function for the asymmetric stretch fundamental has a nodal surface that is a straight line as a result of symmetry. However, the nodal surface of the symmetric stretch v = 1 wave function displays significant curvature. The presence of curvature in a nodal surface is problematic because both the ADMC and the fixed-node DMC calculations used in this study assume that the nodal surface is a straight line defined by the nodal coordinate. Over a series of fixed-node

Figure 1. Projections of the DMC probability distributions onto a2 (a) and θ123 (b) are shown for the ground state (solid black line) as well as the asymmetric stretch excited state calculated with a node placed in a2 and the fixed-node DMC calculation performed in internal (gray solid line) and Cartesian (red dashed line) coordinates.

also exhibiting a node in θ123 and conversely, placing a node in θ123 results in a node appearing in a2. Finally, it is important to reemphasize that the determination of the nodal surface is independent of the choice of coordinates in which to perform the DMC simulation. The optimized node locations, ηnode, and the internal coordinate fixed-node DMC energies of the degenerate asymmetric states with a node in either a1 or a2 are given in Table 2. For the state with a node placed in a1, ηnode = 0, as expected for an asymmetric mode. However, the optimal location for a node placed in a2 is slightly nonzero. Further, the DMC energies of these nominally degenerate states differ by around 20 cm−1 for H+3 and 10 cm−1 for D+3 . Although the energy of the state with a node in a1 is in excellent agreement with the literature value, the energy of the state with a node in a2 is higher than the previously reported value by more than twice the statistical uncertainty of the DMC calculations. Similar behavior was seen in a recent fixed-node DMC study of the CH stretch fundamentals of CH+5 in which we used symmetry adapted linear combinations (SALCs) of the CH bond lengths as the nodal coordinates.13 Specifically, nonzero optimal node locations were found for those asymmetric stretch SALCs constructed out of more than three of the CH bond 7014



Article

performed using these definitions of the nodal coordinates are reported in Table 3 and compared to variational energies Table 3. The Energies, Given in cm−1, of the v = 1 Vibrationally Excited States of H2D+ and D2H+ Obtained from Internal Coordinate Fixed-Node DMC Calculations Compared to the Results of Variational Calculations system H2D

+

H+3

Figure 2. Wave functions are shown for the states of with one quanta in (r12 − r23)2−1/2 (a) or (r12 + r23)2−1/2 (b) obtained from a 2D-DVR calculation performed in the coordinates (r12 ± r23)2−1/2 with θ123 = 60°.

D2H+

H3O−2 ,

DMC studies of fluxional molecules, including and H+5 , we have found that the DMC energies of the v = 1 states of the XH stretches often differ from the results of other computational approaches by 10−50 cm−1.12,15,18 For the larger systems, these differences are attributed to complicated nodal patterns resulting from the fact that the molecular eigenstates have contributions from not only the zeroth-order state that has a single node in the XH stretch coordinate but also a number of zeroth-order states that contain multiple quanta of excitation in lower frequency modes. In H+3 with J = 0, the density of states is quite low in the region of the fundamentals and the situation is simpler. As H+3 is displaced along a1 or θ123, the H−H distances that result in the lowest energy configuration shift from their equilibrium value. This causes curvature in the nodal surface associated with the symmetric stretch coordinate, as illustrated in Figure 2b. This curvature is not captured in the straight line nodal surface used in the ADMC and fixed-node DMC calculations of the excited state with a node in s, leading to the deviations between the fixed-node DMC energies and the literature values. Finally, we conclude this section with a brief summary of two tests of the methodology’s treatment of the v = 1 excited states of H3+ and D+3 . First, we find no statistically significant dependence of the fixed-node DMC results on whether the calculations were performed in Cartesian or internal coordinates. Note that this is true for both the DMC energies and the projections of the DMC probability distributions onto the vibrational coordinates. An example of the latter is shown in Figure 1 while a detailed comparison between the fixed-node DMC energies obtained from calculations performed in Cartesian and internal coordinates is provided in the Supporting Information. Second, we consider the magnitude of the time step dependence of the DMC results and find −1 |E(δτ=0.5) − E(δτ=5) for H+3 DMC DMC | to be consistently less than 7 cm −1 + and 4.2 cm for D3 , values that are either within or slightly outside of the statistical uncertainties. As a result, we find no indication of a significant breakdown in the approximate treatment of either the kinetic energy propagator or the recrossing correction. Further details can be found in the Supporting Information. The v = 1 Excited States of H2D+ and D2H+. For the mixed isotopogues of H+3 , we initially used a1, θ123, and 1 s′ = (r12 + r23) (14) 2 H5O+2 ,

nodal coordinate

Evariationala

EDMC − Evariationalb

s′ θ123 a1 s′ normal mode θ123 normal mode s′ θ123 a1 s′ normal mode θ123 normal mode

2992.5 2205.9 2335.3 2992.5 2205.9 2737.3 1967.8 2079.2 2737.3 1967.8

−442.0 ± 3.9 96.4 ± 10.9 −2.7 ± 4.4 12.4 ± 4.6 9.3 ± 4.3 7.3 ± 3.8 136.1 ± 9.2 −1.9 ± 2.9 26.9 ± 2.4 22.0 ± 2.8

a The variational benchmark values for H2D+ and D2H+ are from Sochi et al. and Miller et al., respectively, and are reported relative to the zero-point energy.44,45 bFive statistically independent DMC simulations were performed in each nodal region, and the reported error bars are 99% confidence intervals.

from the literature.44,45 For both species, the DMC energy for the v = 1 asymmetric stretch excited state is in excellent agreement with the benchmark value, with no statistically significant errors present. However, for H2D+, the DMC energies of the states with a single node placed in either s′ or θ123 have large errors of as much as 442 cm−1 relative to the literature values. A similar error is found in the DMC energy of the state of D2H+ with a node placed in θ123. Interestingly, the DMC energy of the symmetric stretch fundamental of D2H+ is found to be in good agreement with the benchmark value. In order to understand the cause of s′ and θ123 being poor nodal coordinates for the fixed-node DMC calculations of the v = 1 excited states of H2D+, we constructed a two-dimensional projection of the ground state DMC probability distribution onto these coordinates. This distribution is compared in Figure 3 to the two-dimensional projection of the H+3 ground state

Figure 3. Two-dimensional projections of the ground state DMC probability distributions onto s and θ123 for H+3 (a) and onto s′ and θ123 for H2D+ (b) are compared.

DMC probability distribution onto s and θ123. The principal axes of the elliptical H+3 distribution lie parallel to the s and θ123 axes, indicating that these two modes are largely uncoupled. This is consistent with s and θ123 being reasonable nodal coordinates for the v = 1 excited states of H+3 and D+3 . However, considerable coupling between s′ and θ123 is evident in the H2D+ distribution as its principle axes are significantly rotated

to define the nodal surfaces of the v = 1 excited states. Recall that in this notation, atom 2 is the unique isotope. The results of the internal coordinate fixed-node DMC calculations 7015



Article

from the s′−θ123 axis system. It is because of this pronounced coupling, present even at the harmonic level, that s′ and θ123 do not provide reasonable nodal coordinates for the v = 1 excited states of H2D+. Finally, note that the projection of the ground state DMC probability distribution of D2H+ onto s′ and θ123 is similar in appearance to that of H2D+, suggesting a comparable coupling between s′ and θ123 in D2H+. An improved set of nodal coordinates for the v = 1 states of H2D+ and D2H+ in the symmetric modes would need to take into account the coupling between s′ and θ123 shown in Figure 3. Assuming that a significant fraction of this coupling appears at the harmonic level for both species, we performed normalmode analysis in mass-weighted internal coordinates using Wilson’s FG method.48 The resulting curvilinear normal modes are tabulated in the Supporting Information. Note that the normal mode consistent with the asymmetric stretch is identical to a1 by symmetry. Specifically, because a1 has B1 symmetry while s′ and θ123 transform as A1, in the C2v symmetry group, a1 does not couple with the other coordinates at the harmonic level. The results of internal coordinate fixed-node DMC calculations of the v = 1 vibrationally excited states of H2D+ and D2H+ using the curvilinear normal modes to define the nodal surfaces are reported in Table 3, where these coordinates are identified as the s′ normal mode and the θ123 normal mode based on their zeroth-order character. Note that, except for the a1 mode where the node is fully defined by symmetry, optimized node locations were determined by ADMC, and the values of ηnode used in our calculations are given in the Supporting Information. The use of the curvilinear normal modes as the nodal coordinates instead of s′ and θ123 is found to dramatically improve the agreement of the DMC energies with the benchmark values for both excited states of H2D+ as well as the bend fundamental of D2H+. However, the DMC results of the D2H+ symmetric stretch fundamental are found to be in better agreement with the benchmark value when the nodal surface is defined in terms of s′ instead of the normal mode. In order to gain insight into why using the s′ normal mode as the nodal coordinate dramatically improves the agreement between the fixed-node DMC energy and the benchmark value for H2D+ but not for D2H+, we performed 2D-DVR calculations in the s′ and θ123 coordinates with a1 = 0 and the details are provided in the Supporting Information. The resulting wave functions for the v = 1 symmetric stretch and bend excited states of both species are shown in Figure 4. Superimposed on these wave functions are the nodal surfaces used in the fixednode DMC calculations defined in terms of the curvilinear normal modes, shown as the solid black lines, as well as in terms of either s′ or θ123, shown as the gray dashed lines. For H2D+, the nodes present in the wave functions of both excited states are found to closely match the straight line nodal surfaces defined using the normal modes as the nodal coordinates. The presence of slight curvature in the true nodal surfaces of these excited states likely accounts for the small but statistically significant differences between the DMC energies and the benchmark values. Turning now to D2H+, the wave function for the bend fundamental is found to have a similar nodal structure as is present in the analogous wave function for H2D+, with the normal mode providing a noticeably better nodal coordinate than θ123. For the v = 1 symmetric stretch wave function of D2H+, however, the two approximate nodal surfaces defined in

Figure 4. Wave functions of the bend (left column) and symmetric stretch (right column) v = 1 excited states obtained from 2D-DVR calculations are compared for H2D+ and D2H+. Also shown are the straight line nodal surfaces used in the fixed-node DMC calculations defined in terms of the curvilinear normal modes (solid black line) as well as in terms of either s′ or θ123 (dashed gray line).

terms of either s′ or the normal mode are found to be quite similar, with both providing a reasonable description of the nodal surface present in the reduced-dimensional wave function. This is consistent with the fact that when a node is placed in the normal mode associated with the symmetric stretch of D2H+, the projection of the excited state DMC wave function onto s′ also exhibits a node. Likewise, the fixed-node DMC energies calculated using the two choices of nodal surface are comparable, with a difference between them of only 19 cm−1. That the fixed-node DMC energy is in better agreement with the literature value when the node is placed in s′ instead of the corresponding normal mode likely reflects subtle differences in the cancellation of errors associated with using a straight line nodal surface defined in terms of each of these coordinates.

■

CONCLUSION In this paper, we presented an extension of DMC that allows the calculations to be performed in a set of internal, as opposed to Cartesian, coordinates. This is accomplished by, for each walker in the DMC ensemble, holding the G-matrix fixed throughout a single time step and using it to define an instantaneous representation in which the action of the kinetic energy component of the propagator is easy to evaluate. The validity of the approximations underlying this approach can be straightforwardly tested by looking for the appearance of a pronounced time step dependence in the DMC results. Finally, it is worth emphasizing that the approach places no constraints on the choice of internal coordinates other than requiring them all to be independent. We tested the internal coordinate DMC approach on H+3 and its isotopologues, molecules that exhibit large amplitude, zeropoint vibrational motions. For the ground as well as the v = 1 excited states of all four isotopologues, we found excellent agreement between the results obtained from DMC calcu7016



Article

Notes

lations performed in internal and Cartesian coordinates. Furthermore, the fixed-node DMC energies of the fundamentals of H+3 and its isotopologues are in good agreement with the literature values, with the largest errors never exceeding 30 cm−1 and originating from deficiencies in the approximate, straight line nodal surfaces used to define the excited states in the fixed-node DMC calculations.41,44,45 No significant time step dependence was found in any of the DMC results. Therefore, the approximations underlying the internal coordinate DMC methodology appear to be reasonable even for a highly fluxional system like H+3 . We also demonstrated that analysis of ground state DMC probability distributions can provide important insights into what sets of internal coordinates are less strongly coupled and therefore more suitable for use as nodal coordinates for the fixed-node DMC calculations. In doing so, we found that the normal modes provide appropriate definitions for the nodal coordinates, even for a fluxional molecule like H2D+. In particular, the use of the normal modes as nodal coordinates is found to be especially important when the symmetry adapted linear combinations of the valence coordinates are strongly coupled at the harmonic level. Finally, we demonstrated the utility of performing simple reduced-dimensional calculations in order to gain insights into the properties of excited state wave functions and their nodal surfaces, especially when deficiencies in the fixed-node DMC calculations become apparent. The internal coordinate DMC methodology developed in this work provides the foundation for a generalized reduceddimensional DMC approach capable of using essentially any set of internal coordinates as the reduced-dimensional subspace. This will provide researchers with the ability to apply DMC to the specific degrees of freedom that are relevant to the physics of interest without requiring the development and use of a fulldimensional PES. Additionally, the methodology developed here can be applied to minimum energy path DMC calculations, allowing the coordinate connecting the two reacting fragments to be an internal coordinate rather than a Jacobi vector. This eliminates from the analysis complications that can arise when Jacobi vectors are used due to shifts in the centers of mass of the fragments upon deuteration.39,40 In particular, we hope to use the internal coordinate DMC approach developed here to calculate the minimum energy paths of astrochemically relevant reactions like H+3 +H2 and its isotopologues, where zero-point vibrational as well as rotational motion are expected to have an important impact on the dynamics.56−59

■

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Support through Grant CHE-1213347 from the Chemistry Division of the National Science Foundation is gratefully acknowledged. We thank Professor Octavio Roncero for providing us with the codes used to evaluate the potential used in this work. This work was supported in part by an allocation of computing resources from the Ohio Supercomputer Center.

■

ASSOCIATED CONTENT

S Supporting Information *

The analytical expressions for the G-matrix elements and V′, computational details of the reduced-dimensional calculations, two-dimensional cuts of the H+3 PES, and details of the ADMC calculations are provided. Also reported are a detailed comparison between the fixed-node DMC calculations performed in internal and Cartesian coordinates for all four isotopologues as well as a summary of the time step dependence observed in the results of the internal coordinate fixed-node DMC calculations. This material is available free of charge via the Internet at http://pubs.acs.org.

■

REFERENCES

(1) Anderson, J. B. A Random-Walk Simulation of the Schrödinger Equation: H+3 . J. Chem. Phys. 1975, 63, 1499−1503. (2) Anderson, J. B. Quantum Chemistry by Random Walk. H 2P, H+3 D3h 1A1′, H2 3Σ+u , H4 1Σ+g , Be 1S. J. Chem. Phys. 1976, 65, 4121−4127. (3) Austin, B. M.; Zubarev, D. Y.; Lester, W. A., Jr. Quantum Monte Carlo and Related Approaches. Chem. Rev. 2012, 112, 263−288. (4) 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. (5) Suhm, M. A.; Watts, R. O. Quantum Monte Carlo Studies of Vibrational States in Molecules and Clusters. Phys. Rep. 1991, 204, 293−329. (6) McCoy, A. B. Diffusion Monte Carlo Approaches for Investigating the Structure and Vibrational Spectra of Fluxional Systems. Int. Rev. Phys. Chem. 2006, 25, 77−108. (7) McCoy, A. B. Diffusion Monte Carlo Approaches for Studying Systems that Undergo Large Amplitude Vibrational Motionss. In Recent Advances in Quantum Monte Carlo Methods; Anderson, J. B., Rothstein, S. M., Eds.; ACS Symposium Series 953; American Chemical Society: Washington, DC, 2006; Vol. 3, pp 147−164. (8) McCoy, A. B.; Hinkle, C. E.; Petit, A. S. Studying Properties of Floppy Molecules Using Diffusion Monte Carlo. In Recent Advances in Quantum Monte Carlo Methods; Tanaka, S., Rothstein, S. M., William A. Lester, J., Eds.; ACS Symposium Series 1094; American Chemical Society: Washington, DC, 2012; Vol. 4, pp 145−155. (9) Metropolis, N.; Ulam, S. The Monte Carlo Method. J. Am. Stat. Assoc. 1949, 44, 335−341. (10) Ceperley, D. M.; Bernu, B. The Calculation of Excited State Properties with Quantum Monte Carlo. J. Chem. Phys. 1988, 89, 6316−6328. (11) 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. (12) McCoy, A. B.; Huang, X.; Carter, S.; Bowman, J. M. Quantum Studies of the Vibrations in H3O−2 and D3O−2 . J. Chem. Phys. 2005, 123, 064317/1−14. (13) Hinkle, C. E.; McCoy, A. B. Theoretical Investigation of Mode Mixing in Vibrationally Excited States of CH+5 . J. Phys. Chem. A 2009, 113, 4587−4597. (14) Hinkle, C. E.; McCoy, A. B. Characterizing Excited States of CH+5 with Diffusion Monte Carlo. J. Phys. Chem. A 2008, 112, 2058− 2064. (15) 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 H5O+2 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. (16) Petit, A. S.; Wellen, B. A.; McCoy, A. B. Unraveling Rotation− Vibration Mixing in Highly Fluxional Molecules Using Diffusion Monte Carlo: Application to H+3 and H3O+. J. Chem. Phys. 2012, 136, 074101/1−12. (17) Petit, A. S.; Wellen, B. A.; McCoy, A. B. Using Fixed-Node Diffusion Monte Carlo to Investigate the Effects of Rotation−

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. 7017



Article

(40) Hinkle, C. E.; McCoy, A. B. Isotopic Effects on the Dynamics of the CH+3 + H2 → CH+5 → CH+3 + H2 Reaction. J. Phys. Chem. A 2012, 116, 4687−4694. (41) Aguado, A.; Roncero, O.; Tablero, C.; Sanz, C.; Paniagua, M. Global Potential Energy Surfaces for the H+3 System. Analytical Representation of the Adiabatic Ground-State 11A′ Potential. J. Chem. Phys. 2000, 112, 1240−1254. (42) Cencek, W.; Rychlewski, J.; Jaquet, R.; Kutzelnigg, W. Submicrohartree Accuracy Potential Energy Surface for H+3 Including Adiabatic and Relativistic Effects. I. Calculation of the Potential Points. J. Chem. Phys. 1998, 108, 2831−2836. (43) Polyansky, O. L.; Tennyson, J. Ab Initio Calculation of the Rotation−Vibration Energy Levels of H+3 and Its Isotopomers to Spectroscopic Accuracy. J. Chem. Phys. 1999, 110, 5056−5064. (44) Miller, S.; Tennyson, J. First Principles Calculation of the Molecular Constants of H+3 , H2D+, D2H+ and D+3 . J. Mol. Spectrosc. 1987, 126, 183−192. (45) Sochi, T.; Tennyson, J. A Computed Line List for the H2D+ Molecular Ion. Mon. Not. R. Astron. Soc. 2010, 405, 2345−2350. (46) Dennery, P.; Krzywicki, A. Mathematics for Physicists; Dover: New York, 1996. (47) 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. (48) Wilson, E. B.; Decius, J. C.; Cross, P. C. Molecular Vibrations; Dover: New York, 1955. (49) Frederick, J. H.; Woywod, C. General Formulation of the Vibrational Kinetic Energy Operator in Internal Bond-Angle Coordinates. J. Chem. Phys. 1999, 111, 7255−7271. (50) Kosloff, D.; Kosloff, R. A Fourier Method Solution for the Time Dependent Schrödinger Equation as a Tool in Molecular Dynamics. J. Comput. Phys. 1983, 52, 35−53. (51) Charron, G.; Carrington, T., Jr. A Fourier-Lanczos Method for Calculating Energy Levels Without Storing or Calculating Matrices. Mol. Phys. 1993, 79, 13−23. (52) Tannor, D. J. Introduction to Quantum Mechanics: A Time Dependent Perspective; University Science Press: Sausalito CA, 2007. (53) Colbert, D. T.; Miller, W. H. A Novel Discrete Variable Representation for Quantum Mechanical Reactive Scattering via the SMatrix Kohn Method. J. Chem. Phys. 1992, 96, 1982−1991. (54) Luckhaus, D. 6D Vibrational Quantum Dynamics: Generalized Coordinate Discrete Variable Representation and (A)diabatic Contraction. J. Chem. Phys. 2000, 113, 1329−1348. (55) Gardenier, G. H.; Johnson, M. A.; McCoy, A. B. Spectroscopic Study of the Ion-Radical H-Bond in H4O+2 . J. Phys. Chem. A 2009, 113, 4772−4779. (56) Crabtree, K. N.; Indriolo, N.; Kreckel, H.; Tom, B. A.; McCall, B. J. On the Ortho:Para Ratio of H+3 in Diffuse Molecular Clouds. Astrophys. J. 2011, 729, 15/1−12. (57) Oka, T.; Epp, E. The Nonthermal Rotational Distribution of H+3 . Astrophys. J. 2004, 613, 349−354. (58) Roberts, H.; Herbst, E.; Millar, T. J. Enhanced Deuterium Fractionation in Dense Interstellar Cores Resulting from Multiply Deuterated H+3 . Astrophys. J. Lett. 2003, 591, L41−L44. (59) Roberts, H.; Herbst, E.; Millar, T. J. The Chemistry of Multiply Deuterated Species in Cold, Dense Interstellar Cores. Astron. Astrophys. 2004, 424, 905−917.

Vibration Coupling in Highly Fluxional Asymmetric Top Molecules: Application to H2D+. J. Chem. Phys. 2013, 138, 034105/1−11. (18) Lin, Z.; McCoy, A. B. Signatures of Large-Amplitude Vibrations in the Spectra of H+5 and D+5 . J. Phys. Chem. Lett. 2012, 3, 3690−3696. (19) Bowman, J. M.; Czakó, G.; Fu, B. High-Dimensional Ab Initio Potential Energy Surfaces for Reaction Dynamics Calculations. Phys. Chem. Chem. Phys. 2011, 13, 8094−8111. (20) Braams, B. J.; Bowman, J. M. Permutationally Invariant Potential Energy Surfaces in High Dimensionality. Int. Rev. in Phys. Chem. 2009, 28, 577−606. (21) Collins, M. A. Molecular Potential-Energy Surfaces for Chemical Reaction Dynamics. Theor. Chim. Acta 2002, 108, 313−324. (22) Barnes, G. L.; Sibert, E. L., III. The Effects of Asymmetric Motions on the Tunneling Splittings in Formic Acid Dimer. J. Chem. Phys. 2008, 129, 164317/1−9. (23) Nagesh, J.; Sibert, E. L., III. Vibrational Dynamics Around the Conical Intersection: A Study of Methoxy Vibrations on the X̃ 2E Surface. Phys. Chem. Chem. Phys. 2010, 12, 8250−8259. (24) Kjaergaard, H. G.; Garden, A. L.; Chaban, G. M.; Gerber, R. B.; Matthews, D. A.; Stanton, J. F. Calculation of Vibrational Transition Frequencies and Intensities in Water Dimer: Comparison of Different Vibrational Approaches. J. Phys. Chem. A 2008, 112, 4324−4335. (25) Howard, D. L.; Jørgensen, P.; Kjaergaard, H. G. Weak Intramolecular Interactions in Ethylene Glycol Identified by Vapor Phase OH-Stretching Overtone Spectroscopy. J. Am. Chem. Soc. 2005, 127, 17096−17103. (26) Leforestier, C.; Szalewicz, K.; van der Avoird, A. Spectra of Water Dimer from a New Ab Initio Potential with Flexible Monomers. J. Chem. Phys. 2012, 137, 014305/1−17. (27) Partanen, L.; Hänninen, V.; Halonen, L. Ab Initio Structural and Vibrational Investigation of Sulfuric Acid Monohydrate. J. Phys. Chem. A 2012, 116, 2867−2879. (28) Buch, V. Treatment of Rigid Bodies by Diffusion Monte Carlo: Application to para-H2···H2O and ortho-H2···H2O Clusters. J. Chem. Phys. 1992, 97, 726−729. (29) Lee, H.-S.; Herbert, J. M.; McCoy, A. B. Structure and Spectroscopy of NenSH (2Σ+) Complexes Using Adiabatic Diffusion Monte Carlo (ADMC). J. Chem. Phys. 1999, 111, 9203−9212. (30) 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. (31) Gregory, J. K.; Clary, D. C. Calculations of the Tunneling Splittings in Water Dimer and Trimer Using Diffusion Monte Carlo. J. Chem. Phys. 1995, 102, 7817−7829. (32) Lee, E.; Farrelly, D.; Whaley, K. B. Rotational Level Structure of SF6-Doped 4HeN Clusters. Phys. Rev. Lett. 1999, 83, 3812−3815. (33) Severson, M. W.; Buch, V. Quantum Monte Carlo Simulation of Intermolecular Excited Vibrational States in the Cage Water Hexamer. J. Chem. Phys. 1999, 111, 10866−10875. (34) Paesani, F.; Gianturco, F. A.; Whaley, K. B. Microsolvation and Vibrational Shifts of OCS in Helium Clusters. J. Chem. Phys. 2001, 115, 10225−10238. (35) Viel, A.; Whaley, K. B.; Wheatley, R. J. Blueshift and Intramolecular Tunneling of NH3 Umbrella Mode in 4Hen Clusters. J. Chem. Phys. 2007, 127, 194303/1−15. (36) 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. (37) Ramilowski, J. A.; Mikosz, A. A.; Farrelly, D. Rotational Structure of Small 4He Clusters Seeded with HF, HCl, and HBr Molecules. J. Phys. Chem. A 2007, 111, 12275−12288. (38) Deskevich, M. P.; McCoy, A. B.; Hutson, J. M.; Nesbitt, D. J. Large-Amplitude Quantum Mechanics in Polyatomic Hydrides. II. A Particle-on-a-Sphere Model for XHn (n = 4,5). J. Chem. Phys. 2008, 128, 094306/1−13. (39) Hinkle, C. E.; McCoy, A. B. Minimum Energy Path Diffusion Monte Carlo Approach for Investigating Anharmonic Quantum Effects: Applications to the CH+3 + H2 Reaction. J. Phys. Chem. Lett. 2010, 1, 562−567. 7018


Diffusion Monte Carlo in Internal Coordinates - The Journal of Physical

Recommend Documents