Torsion Coupling in H - ACS Publications - American

Article pubs.acs.org/JPCA

Rotation/Torsion Coupling in H5+, D5+, H4D+, and HD4+ Using Diffusion Monte Carlo Published as part of The Journal of Physical Chemistry A virtual special issue “Spectroscopy and Dynamics of Medium-Sized Molecules and Clusters: Theory, Experiment, and Applications”. Melanie L. Marlett, Zhou Lin,‡ and Anne B. McCoy*,§ Department of Chemistry and Biochemistry, The Ohio State University, Columbus, Ohio 43210, United States

Downloaded by UNIV OF NEBRASKA-LINCOLN on August 26, 2015 | http://pubs.acs.org Publication Date (Web): August 24, 2015 | doi: 10.1021/acs.jpca.5b05773

S Supporting Information *

ABSTRACT: Two methods for studying the rotation/torsion coupling in H5+ are described. The first involves a fixed-node treatment in which the nodal surfaces are obtained from a reduced dimensional calculation in which only the rotations of the outer H2 groups are considered. In the second, the torsion and rotation dependence of the wave function is described in state space, and the other internal coordinates are described in configuration space. Such treatments are necessary for molecules, like H5+, where there is a very low-energy barrier to internal rotation. The results of the two approaches are found to be in good agreement with previously reported energies for J = 0. The diffusion Monte Carlo treatment allows us to extend the calculations to low J, and results are reported for the three lowest energy torsion excited states with J ≤ 3. For the level of rotational and vibrational excitation investigated, only modest changes in the vibrational wave functions are found. The effects of deuteration are also investigated, focusing on D5+ and the symmetric variants of H4D+ and HD4+.

■

include H5+,9,15,16 CH5+,10,17 and the NH4+·(H2O) complex.12,13 In studying systems that contain a large amplitude internal rotor, it is common to develop an effective Hamiltonian that treats the torsion as a fourth rotational coordinate:

INTRODUCTION Rotation/torsion spectra of polyatomic molecules provide important insights into their structures. The relative intensities of transitions provide additional insights into the extent to which the electronic structure of the molecule of interest is affected by vibrational excitation. Many connections between structure and spectroscopy may be made by means of the harmonic oscillator/rigid rotor treatments of molecular vibrations and rotations. These treatments are often augmented by low-order perturbation theory.1−6 These approaches have been shown to be highly effective for molecules that undergo small amplitude displacements from their equilibrium structures upon vibrational or rotational excitation. There are many molecules and molecular ions that undergo large amplitude vibrational motions.7−14 The interpretation of their spectra is of interest due to the importance of these molecules in astrochemical and atmospheric processes or other situations where detection relies on high-resolution spectra. Challenges arise in studies of systems that display large amplitude vibrational motions due to the loss of a good zeroorder structure. This is particularly true for an ion like H5+ that undergoes large amplitude excursions from its equilibrium structure even in the ground vibrational state. Although a variety of large amplitude motions can lead to these challenges, we focus on systems with large amplitude torsional motions, in which the ground state probability amplitude is approximately independent of this coordinate due to small barriers to internal rotation. One of the first systems of this type to be investigated spectroscopically was dimethylacetylene.7,8 Other molecules and molecular ions with similarly low-energy torsion barriers © XXXX American Chemical Society

2 2 2 Ĥ = AJÂ + BJB̂ + CJĈ +

pϕ 2 2μϕ

+ V (ϕ) (1)

Here A, B, and C are the rotational constants associated with the overall rotation of the molecule, p̂ϕ provides the momentum associated with the torsion, and V(ϕ) is the potential energy associated with this coordinate. The above description maintains the separation of the rotational and torsional degrees of freedom. Often the dynamics are more easily visualized in an uncoupled representation, where the four rotations are described by the rotations of the two internal rotors (illustrated by the red and blue arrows in Figure 1) with

Figure 1. Structure of H5+ illustrating the coordinates discussed in the text. Received: June 16, 2015 Revised: July 21, 2015

A

DOI: 10.1021/acs.jpca.5b05773 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A angular momentum operators j1̂ and j2̂ along with rotation of the molecule about the two axes that are perpendicular to the internal rotor axis, JB̂ and JĈ . In this representation,


2 2 2 2 Ĥ = BJB̂ + CJĈ + b1 j1̂ + b2 j2̂ + V(ϕ)

in terms of the vibrational coordinate dependent coefficients in the expansion of the rotational state vector in a symmetric top basis. Studies using these two approaches show that for H3+ and its deuterated analogues the fixed-node treatment provides slightly greater accuracy, whereas the basis set approach for describing the rotational contribution to the wave function allows for simultaneous evaluation of multiple rotation/vibration states. This leads to a savings in time of several orders of magnitude. Additional details about both approaches are provided in the following section. In the present study, we extend the above approaches to include torsional degrees of freedom as part of the rotational Hamiltonian. We use this approach to investigate rotation/ torsion couplings in H5+ and selected deuterated analogues. This provides us further opportunities to explore these approaches. In the course of the study, we revisited the socalled recrossing correction that is used to improve the accuracy of the fixed-node DMC calculations. Results will be presented for J ≤ 3 and the first three excited states in vϕ.

(2)

where b1 and b2 are the rotational constants associated with the two rotors. The two forms of the Hamiltonian are related by the fact that j1̂ + j2̂ = JÂ , whereas the torsion angle, ϕ = (ϕ1 − ϕ2)/2. Additional details about the definition of the coordinates is provided in the Supporting Information. In the limit that the torsion potential is constant, the solution to the Hamiltonians in eqs 1 and 2 can be obtained analytically. In the present study, we focus on H5+ and its isotopologues for which b1 = b2. In this case, we can equate K = k1 + k2 and vϕ = k1 − k2, where K, k1, k2, and vϕ are the quantum numbers associated with JÂ , j1̂ , j2̂ , and p̂ϕ, respectively. Because k1 and k2 must both be integers, their sum and difference (e.g., K and vϕ) must both be even or both be odd. The implications of this for H5+ are described in a recent publication.14 In the case of asymmetric rotors, e.g., NH4+·(H2O), the decoupling of the internal rotors leads to rotational structure of the rotation/ vibration transitions that has a rotation constant relative to the A-axis that corresponds to either the NH3 or H2O end of the molecule rotating when the NH or OH stretch is excited, respectively.12,13 In a series of studies,18−22 we investigated the rotation/ vibration couplings in a series of astrochemically important molecules using diffusion Monte Carlo (DMC).23 We have described two such approaches. Both were tested for H3+ and its deuterated analogues and applied to studies of CH5+. The first approach uses vibrationally averaged rotational constants obtained from a ground state DMC calculation24 to calculate the wave functions for rotationally excited states on the basis of an effective rotational Hamiltonian.18,25 These wave functions are then used to determine the locations of the nodal surfaces for the excited states of interest. Fixed-node DMC calculations23,24,26 are performed to evaluate the energies and wave functions for the rotationally excited state of interest. In these calculations, the wave function is constrained to a single nodal region in which the sign of the wave function does not change. Though this approach was demonstrated to be quite accurate, the excited states must be evaluated one at a time, and often calculations need to be performed in several of the nodal regions. This makes the approach very expensive when information about multiple excited states is desired. Further, by employing separate DMC calculations for each state, the corresponding wave functions are obtained at different sets of Monte Carlo sampling points. This makes it difficult to evaluate quantities, such as intensities, which require integrals over products of wave functions for two different states. In a second approach,22 we partitioned the rotation/vibration Hamiltonian as 1 Ĥ = Ĥ vib + 2

■

THEORY To provide a context for the present study, we briefly review the diffusion Monte Carlo (DMC) approach. A more complete description of the theory and our implementation can be found in earlier publications.18,22,24,25,27 Briefly, the algorithm that we use is based on the studies of Anderson.23,27 The approach derives from the parallel structure of the diffusion equation, which includes a coordinate-dependent rate process dC = D∇2 C( r ⃗ ,t ) − k( r ⃗) C( r ⃗ ,t ) dt

and the time-dependent Schrödinger equation, which has been recast in terms of the imaginary time variable, τ = it/ℏ, dΨ = dτ

N

∑ i=1

ℏ2 2 ∇i Ψ( r ⃗ ,τ ) − (V ( r ⃗) − Eref )Ψ( r ⃗ ,τ ) 2mi

(5)

The long-time solution to eq 5 is the lowest energy eigenstate of the Hamiltonian for the N-atomic system of interest, whereas the zero-point energy can be obtained by shifting the value of Eref so that the amplitude of Ψ does not change with τ. In this way, the ground state solution for a Hamiltonian of interest can be obtained using the same Monte Carlo approaches that would be used to study the equilibrium solution to eq 4. On this basis, the wave function is replaced by an ensemble of δ-functions, or walkers, each of which reflects a molecular geometry of the system of interest. The simulation is carried out through a series of small, but finite, time intervals, δτ. During each time step, the values of each of the 3N coordinates of each of the walkers are displaced by a random amount that is consistent with a Gaussian distribution with a width,

∑ Jα̂ (I−1)α ,β Jβ̂ α ,β

(4)

(3)

σi =

Here, Ĥ vib represents the vibrational Hamiltonian, which depends on the 3N − 6 vibrational degrees of freedom. The second term provides the rotation/vibration interactions through the coordinate-dependent moment of inertia tensor (I). With this partitioning, the vibrational part of the wave function is described using standard DMC approaches, whereas the rotational contribution to the wave function is propagated

δτ mi

(6)

where mi is the mass associated with the ith coordinate. The potential is used to determine the relative weight of each walker within the simulation. This can be achieved in one of two ways. The first, termed branching, requires each of the walkers have a weight of 1. After the coordinates of a walker have been displaced, we evaluate B


Article

The Journal of Physical Chemistry A Pn(τ +δτ ) = exp[−(V ( rn⃗ (τ +δτ )) − Eref (τ ))δτ ]

associated with dn(τ), μn. Typically, this mass is given by 1/Gd,d, where Gd,d is the Wilson G-matrix element associated with dn(τ). As such,21,23

(7)

The integer value of Pn gives the number of walkers at that configuration that will be carried forward into the next time step. The fractional part of Pn provides the probability that one additional walker will be generated at this configuration. In the second approach, the weight of each walker can be described by a real number, wn, where wn(τ + δτ ) = Pn(τ +δτ ) wn(τ )

⎧ 2dn(τ ) dn(τ +δτ ) μ n(τ ) μ n(τ +δτ ) ⎫ ⎪ ⎪ ⎬ Pn ,recross = exp⎨− ⎪ ⎪ δτ ⎩ ⎭ (11)

(8)

The evaluation of dn(τ) in eq 11 is straightforward when the node is defined in terms of a single coordinate or a simple linear combination of multiple coordinates. When it becomes a more complicated function of multiple coordinates, or, as is the case of the rotational coordinates the space is non-Euclidian, the evaluation of dn(τ) is much less straightforward. In our previous study of excited states of asymmetric top molecules, we developed an algorithm for addressing this issue.21 Although the approach was effective, it was not straightforward to implement. In the present study, we employ a more general approach. Specifically, we recognize that both the wave function and dn are zero at the node, and the value of the wave function increases linearly with dn near the node. As we have an approximate form for the rotation/torsion wave function that comes from solving eq 1 or 2, we can use the magnitude of this trial wave function, ΨT, as a measure of the distance of the walker from the node. Because the recrossing correction will be important when a walker is near the node, we approximate our trial wave function as

In this scheme, when the weight of a walker falls below the reciprocal of the number of walkers, its coordinates are replaced by those of the walker with the largest weight, and the weight is divided equally between these two walkers. Finally, in both cases,


Eref (τ ) = V̅ (τ ) − α

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

(9)

where V represents the ensemble average of the potential energy and Wtot(τ) provides the sum of the weights of the walkers at that time. The second term in the expression for Eref adjusts the value of Eref to ensure that the value of Wtot is roughly constant. The value of α is chosen to stabilize the simulation. We have found a value of α = 0.5/δτ to be effective for this purpose. The above description provides a method that will generate the ground state wave function and zero-point energy. Two types of approaches are taken to extend this ground state method for the evaluation of rotation/torsion excited states of H5+. They are both based on the rotation/torsion model Hamiltonians provided in eqs 1 and 2 along with the observation that the coupling between the torsion and the other vibrations in H5+ is weak and nearly isotropic. This can be seen by investigating the energy differences between the stationary points that have the two H2 units in the same and in perpendicular planes. On the basis of the potential surface used in this study,15 the torsional barrier ranges from 103 to 111 cm−1 for the six lowest energy stationary points. The nearly isotropic nature of the torsion is further evidenced by the fact that projections of the ground state probability amplitude onto the torsion angle can be fit to P(ϕ) = P0 + P4 cos 4ϕ

N

ΨT( r ⃗) ≈

∑ ∇i⃗ ΨT n ·( ri ⃗ − ri⃗(n)) i=1

(12)

where r(n) represents the coordinates of the ith atom in the nth i⃗ walker. As such, the mass associated with this distance can be evaluated by using 2 −1 ⎡N ⃗ ⎤ ∇Ψ i T ⎥ ⎢ 1 μ= = ⎢∑ ⎥ Gd , d ⎢⎣ i = 1 mi ⎥⎦

(10)

(13)

In the case of rotational motions, the Cartesian coordinates in eq 13 are replaced by the appropriate Euler angles, and the associated mass becomes the moment of inertia of the system with respect to the axis of rotation associated with the Euler angle. Expressions for these masses are provided in the Supporting Information. The results of additional tests of this generalized recrossing scheme are provided in the Supporting Information for ref 29. Because each state is described by a different nodal surface, separate simulations must be performed for each state of interest. Further, because Ĥ Ψ = EΨ must be locally satisfied in all of the nodal regions, representative calculations need to be performed in several nodal regions to test the validity of the trial wave function. A second approach improves the efficiency of the calculations by taking advantage of the approximate separation between the rotation/torsion and the vibrational degrees of freedom that result from the weak dependence of the vibrational potential on ϕ, as described above. In this model, the rotation/torsion contribution to the wave function is defined by

where P4/P0 ranges from 0.13 for H5+ to 0.25 for D5+.28 In this study, we explore two approaches for studying rotation/torsion excited states of H5+. In the first, the solutions to eq 1 or 2 are used to define nodal surfaces for a fixed-node DMC calculation that is similar to those used in studies of rotationally excited states.18,21 This approach exploits the fact that DMC is a ground state method, and within a single nodal region, the wave function behaves in the same way as the ground state solution when the potential for the system being studied is replaced by an infinite potential outside of the nodal region of interest. As the wave function in a region of infinite potential energy must be zero, any walker that moves outside of the region defined by the chosen nodal boundaries is given a weight of zero and removed from the simulation. Because finite time steps are used, a walker that remains on the same side of the node during a single time step may in fact cross the nodal boundary and cross back if a series of smaller time steps had been used. The probability of this occurring for the nth walker depends on the distance of the walker from the node, dn(τ), along with the coordinate-dependent effective mass that is C


Article

The Journal of Physical Chemistry A J, K

simulation is 5 cm−1, although the uncertainties in the relative energies are captured in the reported values. To obtain the probability amplitudes, we use the descendent weighting approach, described by Suhm and Watts.30 For these calculations descendants are gathered for 50 time steps with 3950 time steps between these evaluations of the wave function in the basis set calculations and 35 time steps spaced by 3815 time steps in the fixed-node simulations. In both cases at least 20 sets of descendants were used to evaluate the reported probability amplitudes. In the basis set calculations, states with vϕ ≤ 10 were included in the basis along with K ≤ J. This translates to 11 states for J = 0 and 31 for J = 1. In these calculations, only the diagonal elements of I are included in the Hamiltonian in eq 15. Tests were performed for various cuts through the global surface for H5+ along ϕ and it was found that for the states of interest (vϕ ≤ 6) the energies are converged to better than 0.1 cm−1. Likewise the expansion of the potential in eq 16 was based on the evaluation of the potential at 18 evenly spaced values of ϕ, and the fit was truncated at five terms. Calculations were also performed with higher order expansions, and for vϕ ≤ 4 the energies changed by less than the 5 cm−1 uncertainty of the overall simulation. We performed several calculations including the off-diagonal terms in I and found that this did not affect the calculated energies. Finally, the fixed-node calculations for J > 0 were only performed for PR = +1 as the energies do not depend on the rotational parity. For comparison, we have performed reduced dimensional calculations to obtain the energy levels. For these calculations, r1 = r2 and θ1 = θ2 = 90°, whereas the central proton is constrained to the midpoint of R⃗ in Figure 1. Operationally, these calculations correspond to solving the Hamiltonian in eq 2 in a basis with k1 ≤ 50, k 2 ≤ 50, k1 + k 2 ≤ J, and k1 − k 2 ≤ 21. The values of the parameters used for these calculations and the form of the potential are provided in the Supporting Information.

v PR ϕ PT

⎫ ⎪ [ J , K + PR (− 1)K J , − K ]⎬ ⎪ K , vϕ ⎭ ⎩ 2PR (1 + δK ,0) ⎧ ⎫ ⎪ ⎪ 1 [ vϕ + PT − vϕ ]⎬ ×⎨ ⎪ 2PT(1 + δvϕ ,0) ⎪ ⎩ ⎭ (14) =

⎧ ⎪

∑ cKJ , vϕ( r )⃗ ⎨⎪

1

where PR and PT provide the parities associated with rotation and torsion, respectively, and K and vϕ are both positive integers. The coefficients are propagated using 1 Ĥ rot,tor = 2

∑ Jα̂ (I−1)α ,β Jβ̂

2

+ b1( r ⃗)j1̂ + b2( r ⃗)j2̂

α ,β

b + b2 ̂ 2 JA + V (ϕ; r ⃗) − 1 4 Downloaded by UNIV OF NEBRASKA-LINCOLN on August 26, 2015 | http://pubs.acs.org Publication Date (Web): August 24, 2015 | doi: 10.1021/acs.jpca.5b05773

2

(15)

whereas the remaining 3N − 7 vibrational coordinates are propagated using standard DMC approaches. The details of the simulation approach have been described previously.22 The most significant change is the introduction of a potential energy term in eq 15. For this study, V(ϕ) is evaluated at the position of each walker at each time step. This is accomplished by first identifying the two exterior H2 groups, labeled 1/2 and 4/5 in Figure 1. As shown, R⃗ is defined to connect the centers of mass of these two H2 groups. The H2 group that is further from hydrogen atom 3 is rotated in increments of 10° about R⃗ . The resulting points are fit to V (ϕ) = A( r ⃗) + B( r ⃗) cos 2ϕ + C( r ⃗) sin 2ϕ + D( r ⃗) cos 4ϕ + E( r ⃗) sin 4ϕ

(16)

and A replaces V in eq 7 for the diffusion of the walkers in coordinate space, whereas the remaining terms are used to propagate the coefficients in eq 14. Before moving to the results, we should note that the approaches described above rely on the fact that H5+ contains a well-defined torsion coordinate. Because it is composed of five hydrogen atoms, the potential surface for H5+ is invariant under exchange of the central and outer hydrogen atoms or between the two outer H2 units (e.g., those shown in different colors in Figure 1). If such permutations were energetically accessible for the states of interest, it would be difficult to define a torsion coordinate. On the basis of an analysis of the energetically accessible proton permutations, described in ref 29, such permutations are not expected to occur for the states of interest in the present study.

■

RESULTS AND DISCUSSION For the purpose of this study, we will divide the Results and Discussion into two parts. The first focuses on the case when J = 0, whereas the second considers the situation when J = 1, 2, and 3. Vibrations (J = 0). To test the assumptions made in developing the approaches described in the previous section, we run a standard DMC simulation for H5+, in which the potential V(r)⃗ is replaced by A(r)⃗ , the isotropic term in the expansion of the potential shown in eq 16. This results in a zero-point energy of 7238.3 ± 4.1 cm−1, which is 33 cm−1 larger than the energy obtained when the full potential is used. When higher order terms are included in the expanded potential provided in eq 16, and DMC calculations are performed using the expanded potential, the calculated zero-point energy is 7205.5 ± 2.0 cm−1. This result is in agreement with the previously reported zero-point energy obtained using the full potential, E0 = 7205.2 ± 1.129 cm−1. This gives us confidence that the form of the potential represented by eq 16 provides an accurate representation of the full potential surface. In addition to examining the energies, in Figure 2, we compare projections of the probability distributions onto several of the coordinates indicated in Figure 1. Although small differences are seen, they remain within the uncertainties of the distributions.

■

NUMERICAL DETAILS Unless stated otherwise, all of the DMC calculations were performed with ensembles of 20 000 walkers, τ = 10 au, and α = 0.5/δτ. Calculations of energies were performed for 40 000 time steps, where the system was allowed to equilibrate during the first half of the simulation, and the values of Eref from the last 20 000 time steps were used to provide the zero-point energy. Such simulations were run five times, and the reported energies and uncertainties represent the average and standard deviation of that data. As the uncertainties capture the statistical uncertainties of the simulation, they do not account for systematic errors either in the potential surface or due to the use of the finite time step. A better approximation to the uncertainties in absolute energies resulting from details in the D


Article


The Journal of Physical Chemistry A

Figure 3. Ground state probability amplitude for H5+ projected onto ϕ3, defined in Figure 1 (red curves) and fits to eq 17 (black curves) for (a) ϕ = 0, (b) ϕ = π/8, (c) ϕ = π/4, and (d) ϕ = 3π/8.

Figure 2. Ground state probability amplitudes for H5+ obtained from the basis set approach (solid green lines) and the standard DMC approach (dashed purple lines) projected onto (a) one of the H2 bond lengths, r1, (b) the length of R⃗ , (c) the projection of the central proton onto R⃗ , r3,A, and (d) the bend angle, θ1 as illustrated in Figure 1.

As is seen, there is a small but statistically significant correlation between the dependence of the probability amplitude on these two angles. Before discussing the torsion excited states, we should note that the definition of ϕ = (ϕ1 − ϕ2)/2 differs from the definition that has been used in most of the previous J = 0 studies of H5+ in which ϕ = ϕ1 − ϕ2.31−33 This difference makes our definition of vϕ twice as large as the torsion quantum number used in those studies. In Table 2, we provide the energies for states with vϕ ≤ 6. Although the zero-point energies obtained by the two DMC approaches differ by approximately 30 cm−1, the energies of the torsion excited states are very similar. Comparing these results to previously reported energies, we find that for vϕ = 2 the basis set approach is in better agreement with the MCTDH results, whereas the fixed-node treatment gives energies that are in good agreement with the FCI results. Larger deviations are seen for the states with vϕ = 4. Here the results of the basis set DMC, MCTDH, and GENIUSH calculations all show trends similar to the results of the reduced dimensional treatment. On the contrary, the fixed-node DMC calculations give lower energies and larger splittings. We have explored other choices of embedding. We have also carefully checked the probability amplitude associated with this state to check if it is consistent with a state with vϕ = 4, and not some other state with a similar nodal structure. All of our tests give us confidence that the state we calculate is the state with vϕ = 4 and its energy does not depend on the choice of embedding. We are left to conclude that the energy of this state is particularly sensitive to the way

In Table 1 we compare the zero-point energies obtained when the torsion contribution to the wave function is expanded in a basis to DMC results obtained by using δ-functions to describe all of the vibrational degrees of freedom. We also compare our results to those obtained from multiconfigurational time-dependent Hartree (MCTDH) calculations,31 vibrational configuration interaction (VCI) calculations obtained using the reaction path version of MULTIMODE,32 and full CI (FCI) calculations.33 We believe that the 33 cm−1 shift in the zero-point energy calculated using the basis set DMC approach reflects an incomplete treatment of the coupling between the torsion and the other vibrational degrees of freedom. In particular, this approach neglects couplings between the torsion and displacements of the central proton from the vector that connects the centers of mass of the outer H2 units, R⃗ in Figure 1. To assess the size of this coupling, in Figure 3 we plot projections of the probability amplitude onto ϕ3 for various values of ϕ. Here ϕ3 is defined as the torsion angle associated with the central hydrogen atom (labeled 3 in Figure 1). ϕ3 is defined to be zero when this hydrogen atom lies in the plane defined by R⃗ and the bisector of the dihedral angle between r1⃗ and r2⃗ . To aid in interpretation, the resulting distributions are fit to (0) (2) Pfit(ϕ3) = Pfit + Pfit cos(2ϕ3)

(17)

Table 1. Zero-Point Energies (cm−1) for Four Isotopologues of H5+ Calculated Using Various Methodsa DMC species H5+ D5 + H4D+ HD4+ a

standard 7205.1 5148.9 6858.1 5533.3

± ± ± ±

0.8 0.8 0.5 1.2

basis set 7236.5 5175.9 6886.5 5561.8

± ± ± ±

2.9 3.0 2.5 3.6

MCTDHb

FCIc

VCId

7202.6 5151.0 6854.1 5535.6

7214.6 5152.3 6868.8 5536.4

7244.4 5174.4

Statistical uncertainties are based on five independent DMC calculations. bIn ref 31. cIn ref 33. dIn ref 32. E


Article

The Journal of Physical Chemistry A Table 2. Energies (cm−1) for the J = 0 Torsional States of Four Isotopologues of H5+ DMC species H5

+


D5+

H4D+

HD4+

vϕ

2 2 4 4 6 6 2 2 4 4 6 6 2 2 4 4 6 6 2 2 4 4 6 6

a

PT

− +

− +

fixed node 85.6 124.0 383.0 383.2

± ± ± ±

1.8 1.4 1.9 1.2

34.6 76.8 201.7 205.8

± ± ± ±

1.0 1.1 1.2 1.1

85.0 126.0 383.1 383.5

± ± ± ±

1.3 2.3 1.9 1.8

33.6 76.4 201.7 204.8

± ± ± ±

1.6 1.4 1.4 1.9

− +

− +

− +

− +

− +

− +

− +

− +

− +

− +

basis set 91.0 130.6 439.4 440.9 986.2 986.2 35.8 78.3 220.2 223.4 493.0 493.1 91.8 130.5 441.4 442.8 988.2 988.2 35.8 79.2 221.3 224.5 493.7 493.8

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

4.2 3.6 3.7 3.7 4.3 4.3 4.4 3.8 3.8 3.8 3.7 3.7 3.5 3.5 3.8 3.9 3.2 3.2 4.9 5.3 4.9 4.9 3.7 3.7

RDb

MCTDHc

FCId

VCIe

GENIUSHf

84.0 138.0 436.8 440.1 980.2 980.3 32.2 86.0 222.4 228.9 495.4 495.6 84.0 138.0 436.8 440.1 980.2 980.3 32.2 86.0 222.4 228.9 495.4 495.6

92.6 133.3 447.1 452.3

87.3 138.7 444.0 446.8

65.6 185.4 459.6 503.0

89.9 135.8 446.3 446.9

36.9 81.5 225.0 231.6

27.6 102.7 227.5 254.9 512.7 526.6

a

See eq 14. bParameters for the reduced dimensional calculations are provided in the Supporting Information. cIn ref 31, based on the potential of Xie et al.15 dIn ref 33, based on the potential of Xie et al.15 eIn ref 32, based on the potential of Xie et al.15 fIn ref 9, based on the potential of Aguado et al.16

the calculation is performed. Interestingly, the energies of the corresponding state with J > 0 do not show this large a difference between the two DMC approaches. Given that all but the GENIUSH calculation employed the same potential, the dispersion among the results of the various approaches points to the fact that this is a molecule for which additional study would be useful. For the purpose of the present investigation, the vibrational energies obtained by both the fixed-node and basis set approaches are in very good agreement with each other and with those reported previously. This provides strong validation for the methods used in this study. As in our study of the rotational energies of H3+,22 the fixed-node treatment is more accurate than the basis set approach but is also significantly more computationally expensive. Both approaches afford the opportunity to examine how the wave functions change with torsional excitation. In Figure 4, we plot the projections of the probability distributions for the J = 0 states with vϕ = 0 and 10. As is seen, although the projections onto θ1 and r3,A are not affected by this excitation, there is a small, but significant, shift in the distribution in r1. This is expected because the angular momentum associated with r1⃗ has increased. Likewise, there is a small shift in R to smaller values. This is an indirect consequence of the rotation, as ⟨R⟩ will decrease with increased r1 or r2. We have also performed calculations of J = 0 energies for D5+, HD4+, and H4D+, where the unique atom is in the central position. An interesting and somewhat surprising observation is the similarity of the zero-point energies obtained from the reduced dimensional calculations for all four isotopologues.

Figure 4. Probability amplitude for the states of H5+ with vϕ = 0 (solid blue lines) and vϕ = 10 (dashed red lines) is projected onto (a) one of the H2 bond lengths, r1, (b) the length of R⃗ , (c) the projection of the central proton onto R⃗ , r3,A, and (d) the bend angle, θ1 as illustrated in Figure 1.

Generally, one expects to find a mass dependence as the rigid rotor energies scale as m−1 whereas the harmonic oscillator energies scale as m−1/2. The unusual mass dependence is a reflection of the fact that the potential has only a weak dependence on the torsion coordinate. If we consider ⟨H⟩ in eq 2 by assuming an isotropic ground state wave function, F


Article

The Journal of Physical Chemistry A Table 3. Energies (cm−1) J = 1 for Torsional States of Four Isotopologues of H5+ H4D+

H5+

K vϕ


0 0

a

PT +

1 1 1 1

+

0 2 0 2

+

1 3 1 3

+

0 4 0 4

+

−

−

−

−

1 5 1 5

+

0 6 0 6

+

−

fixed node 6.1 59.1 58.3 99.4 140.5 270.7 269.1 453.1 456.1 702.4 701.1

± ± ± ± ± ± ± ± ± ± ±

basis set

1.4 1.3 1.7 0.9 1.6 1.3 2.2 1.6 1.5 1.5 1.8

6.3 53.2 54.9 97.2 136.8 277.5 277.5 445.6 447.2 715.4 715.4 992.5 992.5

−

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

4.1 4.3 4.5 4.3 3.6 3.4 3.9 3.7 3.7 4.9 4.9 4.3 4.3

RDb

GENIUSHc

K vϕ

7.0 56.7 56.7 91.0 145.0 280.7 280.7 443.8 447.1 712.7 712.7 987.3 987.3

6.4 57.8 57.8 96.3 142.2 283.9 283.9

0 0

+

1 1 1 1

+

0 2 0 2

+

1 3 1 3

+

0 4 0 4

+

1 5 1 5

+

0 6 0 6

+

PT

−

−

−

−

−

fixed node 6.9 58.8 58.2 98.3 141.1 267.9 269.3 451.9 455.1 701.4 701.3

± ± ± ± ± ± ± ± ± ± ±

1.0 1.1 1.8 1.0 1.5 2.1 1.1 1.4 2.4 1.1 1.7

−

D5+

K vϕ 0 0

1 1 1 1

PT +

− +

0 2 0 2

+

1 3 1 3

+

0 4 0 4

+

1 5 1 5

+

0 6 0 6

+

−

−

−

−

−

fixed node 3.0 28.4 27.4 38.9 84.7 126.2 127.3 228.4 229.0 333.1 332.9

± ± ± ± ± ± ± ± ± ± ±

1.7 1.1 1.2 1.0 1.2 1.3 2.0 1.3 2.2 1.2 0.9

basis set 6.2 53.7 53.3 98.0 136.7 278.1 278.5 447.6 449.0 716.2 716.2 994.5 994.5

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

3.5 4.1 3.5 4.7 3.5 3.5 3.6 3.8 3.8 3.2 3.2 3.2 3.2

RDb 7.0 56.7 56.7 91.0 145.0 280.7 280.7 443.8 447.1 712.7 712.7 987.3 987.3

HD4+ basis set 3.2 28.1 24.4 39.0 81.6 142.4 144.6 224.0 227.2 358.1 358.2 496.2 496.3

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

4.3 3.2 3.7 4.4 3.8 3.5 3.9 3.7 3.7 3.7 3.8 3.7 3.7

RDb

K vϕ

3.5 26.9 26.9 35.7 89.5 148.8 148.8 225.9 232.4 362.1 362.1 498.9 499.1

0 0

PT +

1 1 1 1

+

0 2 0 2

+

1 3 1 3

+

0 4 0 4

+

1 5 1 5

+

0 6 0 6

+

−

−

−

−

−

fixed node 2.6 27.4 28.4 38.8 83.5 126.8 127.6 228.5 228.2 333.2 333.3

± ± ± ± ± ± ± ± ± ± ±

2.0 1.7 1.4 1.5 1.6 1.3 1.3 1.2 1.5 1.3 1.4

−

basis set 3.2 26.8 24.7 39.0 82.4 144.6 143.7 225.4 228.6 359.6 359.6 496.9 497.0

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

5.1 4.4 4.5 4.9 5.3 4.4 5.0 5.0 5.0 4.5 4.5 3.7 3.7

RDb 3.5 26.9 26.9 35.7 89.5 148.8 148.8 225.9 232.4 362.1 362.1 498.9 499.1

a See eq 14. As noted in the text, the energies are independent of PR. As a result this quantum number is not provided. bParameters for the reduced dimensional calculations are provided in the Supporting Information. cIn ref 9, based on the potential of Aguado et al.16

ψ0 = const, we find that E0 ≈ ⟨V⟩ and ⟨T⟩ ≈ 0. This anticipates the lack of a mass dependence seen in the ground state energies obtained from the reduced dimensional calculations of the four isotopologues. Similarly, if we consider the difference between E0 obtained from the basis set and standard DMC approaches reported in Table 1, we find that for all isotopologues the difference is between 27 cm−1 for D5+ and 31 cm−1 for H5+. This is a much weaker mass dependence than one might expect, reflecting the weak mass dependence of the zero-point energy in the torsion. The results for torsion excited states for these isotopologues are reported in Table 2. Only MCTDH and VCI calculations have been performed for the D5+ isotopologue. Here both DMC approaches and the results of MCTDH calculations are in good agreement for vϕ = 2. As with H5+, larger differences among the three approaches are found when vϕ = 4. The energies for the mixed isotopologues differ from the corresponding energies for H5+ or D5+ by less than the uncertainties of the DMC energies. This is not surprising as the torsion involves primarily the outer H2 units. Rotationally Excited States. Moving to the J = 1 results, we begin by focusing on H5+. Due to torsion/rotation couplings, vϕ and K must both be even or odd. As a result, the lowest energy 4-fold degenerate K = 1 levels fall between the J = 0 energy levels with vϕ = 0 and 2. The resulting energies are reported in Table 3 for all four isotopologues, and for

vϕ ≤ 6. Energies of states with J = 2 and 3 that have been obtained using the fixed-node approach are reported in Table 4 Table 4. Energies (cm−1) for Various Levels of Torsional Excitation of H5+ with J = 2 J , K vϕ

PT

fixed-node DMC

RD

GENIUSHa

2, 0

+

0

+

20.3 ± 1.4

21.1

19.3

2, 0

+

2

−

112.8 ± 1.7

105.1

109.1

2, 0

+

2

+

153.1 ± 1.5

159.0

155.0

2, 0

+

4

−

465.6 ± 1.5

457.8

2, 0

+

4

+

469.0 ± 2.1

461.2

2, 1

+

1

−

71.3 ± 2.2

70.7

2, 1

+1

+

71.2 ± 1.2

70.7

70.4

2, 1

+

3

−

280.5 ± 1.3

294.8

296.0

2, 1

+

3

+

284.1 ± 1.2

294.8

296.0

2, 1

+

5

−

713.3 ± 1.8

726.7

2, 1

+

5

+

713.3 ± 1.8

726.7

2, 2

+

0

+

118.2 ± 1.8

115.5

117.7

2, 2

+2

−

211.0 ± 1.0

199.5

205.2

2, 2

+2

+

252.3 ± 1.1

253.5

2, 2

+

4

−

563.4 ± 1.5

552.3

552.6b

2, 2

+

4

+

567.7 ± 1.7

555.6

552.6b

70.4

a In ref 9, based on the potential of Aguado et al.16 Comparisons are made on the basis of the energy and the value of K. bIt is unclear which of the 2, 2 + 4 ± states this state correlates to.

G


Article


were used for the two calculations, the overall agreement is generally very good. As we found for J = 1, for states with energies below 150 cm−1, the differences are no more than twice the uncertainty in the DMC calculations. For the higher energy states, the differences become larger.

and Table S1 of the Supporting Information, respectively. To date, the only reported energies for J > 0 were obtained using a different potential surface,9,16 making comparisons with previous calculations a bit more challenging. That said, when we compare the energies for the same states of H5+, for states with energies below 150 cm−1, the energies calculated by DMC and GENIUSH agree to within the reported uncertainties in Table 3 and twice the uncertainties in Table 4. Larger differences are seen for the higher energy state. This likely reflects differences in the potential surfaces used in the two calculations. Comparing the results of the fixed-node and the basis set approaches, we find the agreement to be very good, in general. It should be noted that the J > 0 energies are degenerate. For even values of vϕ and K the K ± vϕ levels are doubly

■

CONCLUSIONS In this study we proposed two approaches for studying the coupled torsion/rotation energy levels of H5+ and its deuterated analogues. The first involves a fixed-node treatment, where the nodes are based on the wave function from a reduceddimensional calculation in which the molecule is treated as two independent rotors coupled by a one-dimensional torsion potential. The second propagates the torsion and rotations in state space and the remaining degrees of freedom in coordinate space. The results of the two approaches are compared to each other and to previously reported energies for these ions. Overall, the agreement is very good, providing further evidence for the relatively weak coupling between the torsion motion and the remaining vibrational degrees of freedom. The results of this work provide a path forward for studying the rotation/ torsion spectrum of this molecular ion that exhibits very large amplitude vibrational motions.


PT

degenerate, whereas the energy should also be independent of PT when K and vϕ are odd. As is seen in the results reported in Table 3, there are small splittings among these energies, but the splittings are smaller than the reported uncertainties. Although these parities define the nodal surfaces used for the fixed-node calculations, the parities are not enforced for the basis set calculations. In the case of states with even values of vϕ, ⟨PR⟩ and ⟨PT⟩ are close to ±1, whereas for states with odd values of vϕ, PT ≈ 0, whereas PR = ±1. The loss of well-defined parity reflects difficulties due to mixing of degenerate states in the basis set calculations. Before discussing the energies, a few comments on the fixednode results need to be made. In calculating these states we used the nodal structure from several reduced dimensional calculations. Specifically, we used the nodal structures obtained from using Eckart embeddings based on reference structures with D2h and D2d symmetries as well as a flexible embedding in which the xz-plane is defined to bisect the dihedral angle between r1⃗ and r2⃗ .9,14 When using these embeddings, we ran into two difficulties. First the calculated energies depended on the choice of embedding, and second, the nodal structure of the calculated states does not always match the nodal structure of the reference state. For example, the calculated states based on the nodal structures for the states with J = 2, K = 2, and vϕ = 0 had energies and nodal structures that corresponded to vϕ = 1 or 2. This resulted from the fact that the definitions of ϕ and χ are not independent in a body-fixed axis system. To circumvent this problem, the nodal surfaces used for the calculations reported here are based on the solutions to the reduced dimensional Hamiltonian in an uncoupled representation, e.g., eq 2, using the D2h saddle point structure as the reference geometry. The resulting energies for J = 1 are provided in Table 3 for vϕ ≤ 6. The agreement between the two DMC approaches is very good for vϕ ≤ 6. Interestingly, the basis set energies are lower when vϕ ≤ 2 or vϕ = 4. At other levels of torsional excitation, the fixed-node treatment provides lower energies. In contrast to the J = 0 results, the energies of the states with vϕ = 4 and 5 that are calculated with the fixed-node approach are similar to those obtained from the basis set and reduced dimensional approaches. Energies for J = 2 for the fixed-node calculations are provided in Table 4, whereas the results for J = 3 are reported in the Supporting Information. We compare the fixed-node energies for states with J = 2 with those reported by Fabri et al.9 Although differences in the quantum number assignments make such comparisons challenging, and different potential surfaces

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.5b05773. Definitions of coordinates in eqs 1 and 2. Reduced dimensional Hamiltonian and masses for the rotational coordinates used for the recrossing correction in eq 13. Selected J = 3 energies for H5+ obtained using the fixed node treatment (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*A. B. McCoy. E-mail: [email protected]. Present Addresses ‡

Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA 02139, USA. § Department of Chemistry, University of Washington, Seattle, WA 98195, USA. 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. We thank Professor Joel M. Bowman for providing us with the codes used to evaluate the potential surface for H5+. This work was supported in part by an allocation of computing time from the Ohio Supercomputer Center. Z.L. thanks the Graduate School at The Ohio State University for fellowship support.

■

REFERENCES

(1) Hoy, A.; Mills, I.; Strey, G. Anharmonic Force Constant Calculations. Mol. Phys. 1972, 24, 1265−1290. (2) Townes, C. H.; Schawlow, A. L. Microwave Spectroscopy; McGraw-Hill Book Company, Inc.: New York, 1955. (3) Nielsen, H. H. The Vibration-Rotation Energies of Molecules. Rev. Mod. Phys. 1951, 23, 90−136.

H


Article



Society: Washington, DC, 2012; Vol. 4, Chapter Studying Properties of Floppy Molecules Using Diffusion Monte Carlo, pp 145−155. (26) McCoy, A. B. 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; Chapter Diffusion Monte Carlo Approaches for Studying Systems that Undergo Large Amplitude Vibrational Motions, pp 147−164. (27) Anderson, J. B. A Random-Walk Simulation of the Schrödinger Equation: H+3 . J. Chem. Phys. 1975, 63, 1499−1503. (28) Lin, Z.; McCoy, A. B. Investigation of the Structure and Spectroscopy of H+5 Using Diffusion Monte Carlo. J. Phys. Chem. A 2013, 117, 11725−11736. (29) Lin, Z.; McCoy, A. B. Probing the Relationship between LargeAmplitude Motions of H5+ and Proton Exchange Between H3+ and H2 J. Phys. Chem. A 2015, DOI: 10.1021/acs.jpca.5b05774. (30) Suhm, M. A.; Watts, R. O. Quantum Monte Carlo Studies of Vibrational States in Molecules and Clusters. Phys. Rep. 1991, 204, 293−329. (31) Valdès, A.; Prosmiti, R.; Delgado-Barrio, G. Vibrational Dynamics of the H5+ and its Isotopologues from Multiconfiguration Time-Dependent Hartree Calculations. J. Chem. Phys. 2012, 137, 214308. (32) Cheng, T. C.; Jiang, L.; Asmis, K. R.; Wang, Y.; Bowman, J. M.; Ricks, A. M.; Duncan, M. A. Mid- and Far-IR Spectra of H+5 and D+5 Compared to the Predictions of Anharmonic Theory. J. Phys. Chem. Lett. 2012, 3, 3160−3166. (33) Song, H.; Lee, S.-Y.; Yang, M.; Lu, Y. Full-Dimensional Quantum Calculations of the Vibrational States of H5+. J. Chem. Phys. 2013, 138, 124309.

(4) McCoy, A. B.; Sibert, E. L. Perturbative Approaches to Highly Excited Molecular Vibrations in H2O, D2O, and HDO. J. Chem. Phys. 1990, 92, 1893−1901. (5) Barone, V. Anharmonic Vibrational Properties by a Fully Automated Second-Order Perturbative Approach. J. Chem. Phys. 2005, 122, 014108. (6) Vazquez, J.; Stanton, J. F. Simple(r) Algebraic Equation for Transition Moments of Fundamental Transitions in Vibrational Second-Order Perturbation Theory. Mol. Phys. 2006, 104, 377−388. (7) Mills, I. M.; Thompson, H. W. Internal Rotation in Dimethyl Acetylene. Proc. R. Soc. London, Ser. A 1954, 226, 306−314. (8) Bunker, P. R.; Longuest-Higgins, H. C. The Infra-Red Spectrum of Dimethylacetylene and the Torsional Barrier. Proc. R. Soc. London, Ser. A 1964, 280, 340−352. (9) Fábri, C.; Sarka, J.; Császár, A. G. Communication: Rigidity of the Molecular Ion H5+. J. Chem. Phys. 2014, 140, 051101. (10) Kolbuszewski, M.; Bunker, P. R. Potential Barriers, Tunneling Splittings, and the Predicted J = 1← 0 Spectrum of CH+5 . J. Chem. Phys. 1996, 105, 3649−3653. (11) East, A. L. L.; Bunker, P. R. A General Rotation-Contortion Hamiltonian with Structure Relaxation: Application to the Processing Internal Rotor Model. J. Mol. Spectrosc. 1997, 183, 157−162. (12) Pankewitz, T.; Lagutschenkov, A.; Niedner-Schatteburg, G.; Xantheas, S. S.; Lee, Y.-T. Infrared Spectrum of NH4+(H2O): Evidence for Mode Specific Fragmentation. J. Chem. Phys. 2007, 126, 074307. (13) Kelleher, P. J.; Johnson, C. J.; Fournier, J. A.; Johnson, M. A.; McCoy, A. B. Persistence of Dual Free Internal Rotation in NH+4 (H2O)·Hen=0−3 Ion−Molecule Complexes: Expanding the Case for Quantum Delocalization in He Tagging. J. Phys. Chem. A 2015, 119, 4170−4176. (14) Sarka, J.; Fabri, C.; Szidarovszky, T.; Csaszar, A. G.; Lin, Z.; McCoy, A. B. Modeling Rotations, Vibrations and Rotational Couplings in Astructural Molecules: A Case Study Based on the H+5 Molecular Ion. Mol. Phys. 2015, 113, 1873−1883. (15) Xie, Z.; Braams, B. J.; Bowman, J. M. Ab initio Global PotentialEnergy Surface for H5+ →H3+ + H2. J. Chem. Phys. 2005, 122, 224307. (16) Aguado, A.; Barragan, P.; Prosmiti, R.; Delgado-Barrio, G.; Villarreal, P.; Roncero, O. A new Accurate and Full Dimensional Potential Energy Surface of H5+ based on a Triatomics-in-Molecules Analytic Functional Form. J. Chem. Phys. 2010, 133, 024306. (17) Schreiner, P. R.; Kim, S. J.; Schaefer, I.; H, F.; Schleyer, P. v. R. M. CH5+: The Never-Ending Story or the Final Word? J. Chem. Phys. 1993, 99, 3716−3720. (18) Petit, A.; 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. (19) Petit, A. S.; Wellen, B. A.; McCoy, A. B. Unraveling RotationVibration Mixing in Highly Fluxional Molecules Using Diffusion Monte Carlo: Applications to H3+ and H3O+. J. Chem. Phys. 2012, 136, 074101. (20) 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. (21) 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. (22) Petit, A. S.; Ford, J. E.; McCoy, A. B. Simultaneous Evaluation of Multiple Rotationally Excited States of H+3 , H3O+, and CH+5 Using Diffusion Monte Carlo. J. Phys. Chem. A 2014, 118, 7206−7220. (23) 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. (24) McCoy, A. B. Diffusion Monte Carlo for Studying Weakly bound Complexes and Fluxional Molecules. Int. Rev. Phys. Chem. 2006, 25, 77−108. (25) McCoy, A. B.; Hinkle, C. E.; Petit, A. S. In Recent Advances in Quantum Monte Carlo. Methods; Tanaka, S., Rothstein, S. M., William A. Lester, J., Eds.; ACS Symposium Series 1094; American Chemical I


Torsion Coupling in H - ACS Publications - American

Recommend Documents