Studying Properties of Floppy Molecules Using ... - ACS Publications

Chapter 12

Downloaded by STANFORD UNIV GREEN LIBR on May 16, 2012 | http://pubs.acs.org Publication Date (Web): February 6, 2012 | doi: 10.1021/bk-2012-1094.ch012

Studying Properties of Floppy Molecules Using Diffusion Monte Carlo Anne B. McCoy,* Charlotte E. Hinkle, and Andrew S. Petit Department of Chemistry, The Ohio State University, 100 W. 18th Ave., Columbus, Ohio 43210 *E-mail: [email protected]

Recent extensions to the Diffusion Monte Carlo approach are presented and discussed in the context of polyatomic molecules that undergo large amplitude vibrational motions, even in their ground states. The methods specifically involve the evaluation of rotationally excited state energies and wave functions, and are discussed in the context of fixed-node studies of H3O+ and CH5+.

Introduction An area in which Diffusion Quantum Monte Carlo (DMC) provides great promise and one for which the methodology has not been fully exploited is in studies of molecules or clusters that contain one or more vibrational degrees of freedom that sample large regions of the potential surface even in the vibrational ground state (1–3). The advantage of using DMC to study these systems, that is not realized in electronic structure calculations, comes from the fact that the analysis of the results of calculations of the bosonic (nodeless) ground state provides insights into properties of the molecules of interest. While most of the contributions to this volume focus on Quantum Monte Carlo approaches for electronic structure calculations, several have used DMC approaches to study rotation-vibration states. Initially this line of investigation focused on hydrogen bonded clusters (1–8). In many of these studies, the intramolecular degrees of freedom were kept fixed, and the studies focused on the intermolecular degrees of freedom (9). An attractive feature of these systems was the availability of potential surfaces, which could be readily evaluated for dimer systems and extended to larger clusters using sums of pair-wise and perhaps higher order interactions. © 2012 American Chemical Society In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.


Over the past decade, interest in using DMC in studies of vibrational properties has moved from molecular clusters to protonated systems, like CH5+ (2, 6, 7, 10) and protonated water clusters (2, 6, 11, 12). These are systems that include several large amplitude vibrations and in which many of the 3N-6 vibrational degrees of freedom are strongly coupled. Often the ground state is delocalized over multiple minima on the potential surface with significant amplitude in the region of the transition states that connect these minima. This type of system is difficult to treat by conventional methods since most of these approaches require a well-defined reference structure. On the other hand, the delocalization of the ground state wave function along the vibrational degrees of freedom is not nearly as extensive as is seen in electronic structure problems. In the most extreme case of CH5+, the ground state wave function has similar amplitude in all 120 equivalent minima and in the region of the 180 low-energy saddle points that connect these minima (10). This would lead one to expect that a particle on a sphere basis might be effective for describing the vibrational dynamics in CH5+ (13). On the other hand, while the amplitudes of the displacements sampled in the coordinates that interconnect the minima are large even in the vibrational ground state, the overall envelope of the CH stretch region of the spectrum is well-described by the average of the harmonic spectra obtained at these 300 stationary points (14). The above characteristics of CH5+, which make it difficult to treat by conventional methods also make it well-suited for studies based on Diffusion Monte Carlo approaches. Two sets of challenges are posed by these systems. First, unlike molecular clusters, the potential needs to be separately constructed for each molecule of interest. Second, DMC, in its purest form, provides a method for obtaining the zero-point energy and a Monte Carlo sampling of the ground state wave function. How can a ground state method be used to provide insights into molecules for which much of our understanding of their properties is based on spectroscopic and kinetic studies, which sample rotationally and vibrationally excited states? The first challenge has been overcome, at least in part, through sophisticated methodologies for obtaining multi-dimensional potential surfaces from numerous electronic energies, and in some cases gradients and Hessians (15–17). All of this required data can be generated by any one of a number of commonly available electronic structure programs, and can be evaluated at a variety of levels of theory, depending on the size of the system and the level of accuracy that is desired. With this information in hand, Braams, Bowman and their co-workers fit the electronic energies to high-order polynomial expansions in functions of the atom-atom distances. The fit is constrained to ensure that the potential energy is not affected by exchange of identical particles (17). This is the form of the potentials used in the studies of CH5+ (18) and H3O+ (19), described below. In an alternative approach, Collins, Jordan and co-workers GROW their potential by generating a series of quadratic expansions of the potential in terms of inverse bond-lengths and using Shephard Interpolation to obtain electronic energies at intermediate points (15, 16). This latter approach is ideally suited for DMC studies, since, while the interpolation can lead to artificial wiggles in the potential, the statistical sampling of the wave function, generated by the DMC approach, tends to be insensitive to these small oscillations. 146 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.


The challenge of evaluating excited states using a fundamentally ground state method is more fundamental. This has been a focus of our research over the past decade. Ways in which we have addressed some of these questions is the focus of the present article. In the previous volume of this series (20), we outlined our use of the fixed-node treatment, as described by Anderson (21, 22) to calculate vibrational excited states. In the present article we turn our attention to the calculation of rotationally excited states. As mentioned above, the focus of this work will be on systems that display large amplitude excursions from the equilibrium structure, even in the rotational and vibrational ground state, specifically H3O+ and CH5+.

General Background Diffusion Monte Carlo approaches and their applications to a variety of systems are described throughout this volume. In our work we use the simplest form of DMC. We base our implementation on the two early studies of Anderson (21)(22). Specifically, we take advantage of the fact that when the time-dependent Schrödinger equation is expressed in terms of an imaginary time variable, τ = it/ħ, the long-time solution converges to the lowest energy solution to the time-independent Schrödinger equation, and the amplitude of the wave function decays with a rate that is proportional to the ground state energy. If the wave function that is being propagated is expressed by an ensemble of δ-functions (or walkers), the positions of which are allowed to move in space under the influence of the imaginary-time time-dependent Schrodinger equation, the resulting equilibrium distribution of walkers represents a Monte-Carlo sampling of the ground state wave function. A more detailed description of this approach and our implementation may be found in our earlier publications (2, 23–25). As described in Refs. (20) and (22), one way to generate excited states is through a fixed-node treatment. Here we exploit the fact that if the nodal surface can be described by a simple functional form, the wave function in the vicinity of the nodal surface has similar properties to a wave function in the vicinity of an infinite potential barrier. As such, within our implementation of the fixed-node treatment, we divide configuration space into two parts using the nodal surface as the dividing surface. The two parts correspond to the regions in which the wave function has positive and negative amplitude. The “true” potential surface is used in one region, while the potential is made to be infinite in the other. Because the amplitude of the wave function is zero in regions of configuration space where the potential is infinite, any walker that moves into the region of configuration space where the potential is infinite is removed from the simulation. In practice, when the nodal surfaces divide the wave function into three or more parts, we choose to treat each region separately. Since we use finite time steps in the simulation, typically 1-10 atomic units, there is a finite probability that a walker that appears to have remained on the same side of the dividing surface could have crossed over to the forbidden region and back within a single time step. If the trajectory had been divided into smaller 147 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.


time increments, this walker would have been removed from the simulation. A statistical measure of the probability of this occurring, based on the distance of the walker from the dividing surface before and after a time step in the simulation, was developed by Anderson (22), and this recrossing correction is included in most of our simulations. We have demonstrated the effectiveness of this approach for vibrational problems through direct comparison of DMC calculations of vibrationally excited states to variational results obtained for the same surfaces for Ne2OH, Ne2SH (26) and H3O2− (12).

Direct Evaluation of Rotationally Excited States The approach described above can be readily extended to rotationally excited states since all that is required is knowledge of the functional form for the nodal surfaces, and, if the recrossing correction is to be used, an expression for the reduced mass associated with motion perpendicular to the nodal surface. The solutions to the rotational Hamiltonian for a rigid molecular system are well-known, and are expressed as functions of the three Euler angles that relate the molecule- and space-fixed Cartesian coordinate systems. Specifically, if the rotational Hamiltonian is expressed as

analytic solutions exist in three cases: • •

•

When J=0, Erot = 0 and the wave functions are isotropic. When J=1, there are three solutions with energies of A+B, A+C and B+C. The corresponding wave functions have analytical solutions that are independent of the values of A, B, and C. For symmetric (or spherical) tops with higher values of J, the wave functions can be expressed by the Tesseral harmonics, which are the set of real functions that are eigenfunctions of Eq. (1) when either A = B or B = C.

With the excpetion of the J = 0 case, the definitions of the rotational wave functions require definitions for the arguments of the Tesseral harmonics, θ and χ, which in turn require an embedding of a body-fixed axis system. At first thought the embedding of a body-fixed axis system in a highly fluxional molecule seems like an insurmountable problem. In earlier studies of CH5+ (27), NenOH (28) and H5O2+ (2), we have demonstrated that application of the Eckart conditions (29) through a series of linear algebra manipulations provides a robust internal coordinate system that minimizes the presence of rotation-vibration coupling. The procedure is relatively straight-forward and can be easily applied to the DMC wave functions. To start, we choose a reference geometry, translate it so the center of mass is at the origin and rotate the molecule to a principle axis frame by ensuring that the moment of inertia tensor is diagonal in this coordinate system. The 3N Cartesian 148 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.

coordinates that define the reference geometry are denoted by the vectors . Following Louck and Galbraith (30), for an arbitrary structure denoted by a set of N vectors

we define a set of vectors


where the α and β indices represent the Cartesian components of the vectors (x,y,z). From these we generate a set of unit vectors

where

These vectors provide the three coordinate axes in the Eckart frame so that the Cartesian coordinates of the atoms in the Eckart frame are

One of the most rigorous tests of Eckart embedding was a study of vibrationally averaged rotational constants for CH5+ and its deuterated analogues , , and for the (27). In that study, we calculated the values of ground state of these ions by averaging the values of the elements of the inverse moment of inertia tensor, evaluated using an Eckart frame defined using the procedure described above. We showed that we obtained the same values for the rotational constants independent of the reference geometry used to define the body-fixed axis system. There were examples where this broke down, but they corresponded to using stationary point structures that were not sampled by the ground state wave function as the reference geometry. With the embedding in hand, we use standard fixed-node approaches in which we assume that the nodal structure is defined by the solutions to the Hamiltonian for a rigid rotor. Specifically, the nodes are placed at the zero’s of the Tesseral Harmonics:

for chosen values of J and K. For simplicity, we require that M=0, although calculations were performed for J=1 and non-zero values of M (24). The values 149 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.


of θ and χ used to define these functions (and the nodes) are extracted from the Eckart rotation matrix (f) in Eq. (3). Specifically

Before introducing the nodal surfaces, we run a long equilibration simulation for the J = 0 wave function. This simulation needs to be propagated for sufficiently long times to ensure that the wave function is isotropic. The final ingredient in the DMC simulations is the introduction of the recrossing correction. Specifically, the probability of a walker crossing and then recrossing a nodal surface within a single time step is given by

where d(τ) and d(τ+Δτ) are the perpendicular distances of the walker from the nodal surface at the start and end of the time step (22). In our studies on H3O+, we introduced the following definitions for the effective mass for a node in θ

and in χ

where I(τ) is the instantaneous moment of inertia tensor for a given walker evaluated in the Eckart frame. Based on our studies on H3O+, we found that this correction does not affect the calculated energies for the J=1 levels, and, as a result, it was not used in the calculations involving CH5+ where the definitions of the separate rotational constants are less straight-forward (24).

Applications to H3O+ and CH5+ The approach, described above, has been recently applied to studies of H3O+ (24), CH5+ (25), and their deuterated analogues. One challenge that we encountered results from the fact that both systems explore broad regions of the potential surface. Specifically, consider the equilibrium structures of these two molecules, shown in the left column of Figure 1, as well as the low-energy saddle points, shown to the right of the corresponding minimum energy structure. In the case of CH5+, the structures on the left and right each represent one of 120 equivalent structures, while the one in the center represents one of 60 equivalent structures. The symmetry is lower for H3O+, but here there are two equivalent minima on the potential separated by a saddle point. 150 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.


Figure 1. Equilibrium (left) and low-energy transition state (center and right) structures of CH5+ (top) and H3O+ (bottom). Without doing any calculations we can anticipate the delocalization of the ground state wave function by comparing the energies of these saddle points, relative to the potential minimum, to the harmonic zero-point energy associated with the mode that connects the minimum to the saddle point geometry. For CH5+, the energies of the two saddle points are 29.1 and 340.7 cm-1 above the global minimum, while the harmonic zero-point energy in the associated modes are 99.9 and 419.7 cm-1 as determined from the potential surface of Jin, Braams and Bowman (18). For H3O+, the energy of the saddle point is 690 cm-1 while the zero-point energy in the out-of-plane vibration is 445.5 cm-1 (19). In both cases the zero-point energy is close to or larger than the energy of the corresponding saddle point. If we analyze the fully anharmonic ground state probability amplitude obtained from DMC, we find that there is significant amplitude in all of the equivalent minima and at the saddle points that connect the minima (10, 24), as was anticipated in the harmonic analysis. This behavior is also reflected in the lowest energy vibrationally excited state in both systems being substantially non-zero (46.4 cm-1 for H3O+ (24) and 10.4 cm-1 for CH5+ (31)) but considerably smaller than the lowest calculated harmonic frequency of 891 cm-1 in H3O+ (19) and 199.9 cm-1 in CH5+ (18). Such behavior is indicative of a sizable tunneling splitting of the ground state and is consistent with the ground state wave function having amplitude at the saddle point as well as the minimum energy structures. This large delocalization of the ground state wave function, on its own, does not cause any problems for the evaluation of rotationally excited states. On the other hand, if we identify the atoms in each of these ions, as is done when we imposed a numbering scheme in Eqs. (2) to (5), the equilibrium structures are chiral, and, based on the above discussion, the ground state wave function has amplitude along the path that connects sets of chiral pairs. The reason this becomes a complicating factor in evaluating rotationally excited states is due to the fact that when a walker samples certain regions of configuration space, the algorithm for the Eckart embedding not only rotates the molecule, it also inverts the axis system. This is more easily seen for H3O+, and is illustrated in Figure 2. The top two diagrams represent the reference structure for H3O+ (left) and the structure that one would obtain with excitation of the out-of-plane (umbrella) mode. No rotation of the axis system will bring the structure on the right close to that on the left. If 151 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.


instead we first inverted the axis system then performed a 180° rotation (following the arrows through the intermediate structure), we return to the structure on the left. Numerically, the algorithm that we use for the Eckart embedding generates a rotation matrix that has a determinant of –1, and we cannot use this matrix to determine the Euler angles used to define the nodal surfaces. While the problem is easily identified, the solution is less obvious. The only one we have found to be uniformly effective is to use a reference structure that does not contain a chiral center. This can be achieved by using the transition state structure for H3O+, shown in the lower right of Figure 1. However, there are no such achiral, energetically accessible stationary point structures for CH5+.

Figure 2. An illustration of the relationship between “inversion” of H3O+ and an inversion of the axis system. Energies for CH5+ and H3O+, calculated by this method, for J=0 and J=1 are reported in Table I. For comparison, we also report the energies that are obtained for J=0 when we remove all inverted walkers. To differentiate these states, we denote this state with a superscripted (-) and the ground state with a (+). For comparison, results performed using variational approaches with the same potential surface are also reported (19, 31). Calculations for the deuterated analogues of both species have been performed and the results can be found in the original publications (24, 25). As indicated by the notation, and discussed above, since there is no achiral reference structure that could be chosen for calculations of CH5+, the rotationally excited states for CH5+ have a vibrational wave function that is antisymmetric with respect to inversion (25). We have also performed calculations for higher J levels in H3O+ (24), and, more recently on H3+ and D3+. Generally the agreement remains good at all values of J and we are investigating whether the origins of any break-downs that are observed can be attributed to problems with the model that we employ for defining the nodal surfaces or rotation-vibration interactions that are not well-described by this model. While the agreement is much less good at J=10 than for J < 6, to our knowledge these systems have more nodes in the rotation-vibration wave function than any previous fixed-node studies of nuclear dynamics. 152 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.

Table I. J = 0 and 1 energies for H3O+ (24) and CH5+ (25) calculated using DMC are compared to previously reported energies obtained using the same potential surfaces Molecule


H3O+

CH5+

Statea

DMC (cm-1)

Variationalb (cm-1)

7452±2

7450.6

44±5

46.4

23±4

22.3

20±3

17.4

19±3

17.4

10916±5

10926.5

11±2

10.4

19±1 20±1 19±1

The states are labeled as , where and correspond to the lower and upper members of the ground state tunneling doublet respectively and the subscript ± provide the parity of the rotational level. b References (19) and (31). a

Concluding Remarks In this paper, we discuss a method for describing rotationally excited states of molecules that undergo large amplitude motions that was recently developed in our group. The above approach has several advantages over correlation functionbased approaches, previously employed by Whaley and Roy and their co-workers to study rotationally excited states of molecules embedded in helium droplets (4, 5). Most important, within the context of floppy molecular systems, the present approach allows for a non-rigid molecule. Further, in contrast to approaches that are based on linear response theory, we are able to obtain the probability amplitude associated with the excited states in a straightforward manner. We are presently working to further characterize the properties of this algorithm through studies of 153 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.

H3+ and D3+ and NH3 for which rotation vibration energies for J > 5 have been evaluated using the same potential surfaces.

Acknowledgments


Support for this work by the National Science Foundation (CHE-0848242) is gratefully acknowledged, as is allocation of computational resources from the Ohio Supercomputing Center with which many of the calculations described above were performed.

References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.

17. 18. 19. 20.

21.

Suhm, M. A.; Watts, R. O. Phys. Rep. 1991, 204, 293. McCoy, A. B. Int. Rev. Phys. Chem. 2006, 25, 77. Sandler, P.; Buch, V.; Sadlej, J. J. Chem. Phys. 1996, 105, 10387. Moroni, S.; Blinov, N.; Roy, P.-N. J. Chem. Phys. 2004, 121, 3577. Blume, D.; Lewerenz, M.; Whaley, K. B. J. Chem. Phys. 1997, 107, 9067. Cho, H. M.; Singer, S. J. J. Phys. Chem. A 2004, 108, 8691. Thompson, K. S.; Crittenden, D. L.; Jordan, M. J. T. J. Am. Chem. Soc. 2005, 127, 4954. Severson, M. W.; Buch, V. J. Chem. Phys. 1999, 111, 10866. Buch, V. J. Chem. Phys. 1992, 97, 726. Brown, A.; McCoy, A. B.; Braams, B. J.; Jin, Z.; Bowman, J. M. J. Chem. Phys. 2004, 121, 4105. McCoy, A. B.; Huang, X.; Carter, S.; Landeweer, M. Y.; Bowman, J. M. J. Chem. Phys. 2005, 122, 061101. McCoy, A. B.; Huang, X.; Carter, S.; Bowman, J. M. J. Chem. Phys. 2005, 123, 064317. Deskevich, M. P.; McCoy, A. B.; Hutson, J. M.; Nesbitt, D. J. J. Chem. Phys. 2008, 128, 094306. Huang, X.; McCoy, A. B.; Bowman, J. M.; Johnson, L. M.; Savage, C.; Dong, F.; Nesbitt, D. J. Science 2006, 311, 60. Ischtwan, J.; Collins, M. A. J. Chem. Phys. 1994, 100, 8080. Jordan, M. J. T.; Crittenden, D. L.; Thomson, K. C. Quantum Effects in Loosely Bound Complexes. In Advances in Quantum Monte Carlo; Anderson, J. B., Rothstein, S. M., Eds.; ACS Symposium Series 953; American Chemical Society: Washington, DC, 2006; pp 101−142. Braams, B. J.; Bowman, J. M. Int. Rev. Phys. Chem. 2009, 28, 577. Jin, Z.; Braams, B. J.; Bowman, J. M. J. Phys. Chem. A 2006, 110, 1569. Huang, X.; Carter, S.; Bowman, J. J. Chem. Phys. 2003, 118, 5431. McCoy, A. B. Diffusion Monte Carlo Approaches for Studying Systems That Undergo Large Amplitude Vibrational Motions. In Advances in Quantum Monte Carlo; Anderson, J. B., Rothstein, S. M., Eds.; ACS Symposium Series 953; American Chemical Society: Washington, DC, 2007; pp 147−164. Anderson, J. B. J. Chem. Phys. 1975, 63, 1499. 154 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.


22. Anderson, J. B. J. Chem. Phys. 1976, 65, 4121. 23. McCoy, A. B. Chem. Phys. Lett. 2000, 321, 71. 24. (a) Petit, A. S.; McCoy, A. B. J. Phys. Chem. A 2009, 113, 12706. (b) Petit, A. S.; McCoy, A. B. J. Phys. Chem. A 2011, 115, 9325. 25. Hinkle, C. E.; Petit, A. S.; McCoy, A. B. J. Mol. Spectrosc. 2011, 268, 189–198. 26. Lee, H.-S.; McCoy, A. B. J. Chem. Phys. 2002, 116, 9677. 27. Johnson, L. M.; McCoy, A. B. J. Phys. Chem. A 2006, 110, 8213. 28. Lee, H.-S.; McCoy, A. B. J. Chem. Phys. 2001, 114, 10278. 29. Eckart, C. Phys. Rev. 1935, 47, 552. 30. Louck, J. D.; Galbraith, H. W. Rev. Mod. Phys. 1976, 48, 69. 31. Wang, X.-G.; Carrington, T., Jr. J. Chem. Phys. 2008, 129, 234102.

155 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.

Studying Properties of Floppy Molecules Using ... - ACS Publications

Recommend Documents