Evaluation of Matrix Elements Using Diffusion Monte Carlo Wave

Subscriber access provided by UNIV OF LOUISIANA

Article

Evaluation of Matrix Elements Using Diffusion Monte Carlo Wave Functions Victor Guang Ming Lee, Lindsey R Madison, and Anne B. McCoy J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b11213 • Publication Date (Web): 25 Apr 2019 Downloaded from http://pubs.acs.org on April 26, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Evaluation of Matrix Elements Using Di↵usion Monte Carlo Wave Functions Victor G. M. Lee,† Lindsey R. Madison,‡ and Anne B. McCoy

⇤,†

†Department of Chemistry, University of Washington, Seattle, WA 98195, USA ‡Department of Chemistry, Colby College, Waterville, ME 04901, USA E-mail: [email protected] Phone: 206-543-7464

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Approaches for using di↵usion Monte Carlo (DMC) to evaluate matrix elements involving two vibrational wave functions are explored. In the first part of this study, overlaps between a wave function obtained using DMC and one that can be calculated analytically are evaluated. In this case, the analytical wave function is used as the guiding function for an importance sampled DMC simulation. The accuracy of the calculated overlaps is found to depend strongly on the accuracy of the calculated descendant weights, which are obtained in the DMC simulation. While a single evaluation of the descendant weights is sufficient for obtaining projected probability amplitudes or expectation values of multiplicative operators, averages of multiple independent evaluations of the descendant weights are required to obtain accurate matrix elements. This approach is investigated for one-dimensional model systems as well as H2 CO, H2 D+ and D2 H+ . The approach is extended to the evaluation of matrix elements of the dipole moment operator between the ground state and states with one quantum of excitation in one of the OH stretching vibrations in H3 O2 . For these calculations, the wave functions for both the ground and excited states are evaluated using DMC. The described methodology opens the possibility of evaluating matrix elements involving two di↵erent states, both of which are obtained using DMC.

2


Page 2 of 31

Page 3 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Quantum mechanical calculations of vibrational energies and wave functions provide unique challenges for chemists. 1–12 While the potential energy for a molecular system can be described as a sum of two-body Coulombic interactions, once the Born-Oppenheimer separation of the electronic and nuclear degrees of freedom is introduced there is no longer a simple and general functional form for the potential energy for the nuclear degrees of freedom. 13 Within the Born-Oppenheimer approximation, the potential energy for the electronic degrees of freedom remains Coulombic, and the electronic wave function can be expanded as a product of functions that are expressed in terms of the three spherical coordinates of each electron. This expansion is based on the LCAO-MO model and provides important insights into the electronic structure of a broad range of chemical systems. Just as the potential energy surface for the nuclear degrees of freedom cannot be expressed in a simple generic form, it is difficult to identify a single set of coordinates and basis functions that are appropriate for all vibrational problems. If we could solve for the vibrational energies and wave functions exactly, the choice of coordinates would not a↵ect the final results. Practically, numerical evaluations of the energies and wave functions require a choice of coordinates and basis functions. The convergence properties of most vibrational calculations depend on this choice. 14 One method that does not require a carefully chosen set of coordinates and basis set is di↵usion Monte Carlo (DMC). 15–18 In this method, Monte Carlo approaches are used to evaluate the ground state wave function and zero-point energy of the system of interest. Because Monte Carlo approaches scale favorably with the system size, these calculations are typically performed in Cartesian coordinates, and the wave functions are expanded in a basis of localized functions of the 3N Cartesian coordinates that define the nuclear geometry. 19,20 The ground state energy and wave function evaluated using DMC provides a benchmark, which can be used to explore the accuracy of the energies and wave functions obtained by other approaches. Such an assessment of the quality of the wave functions requires the 3



Page 4 of 31

ability to evaluate overlap integrals between the wave functions obtained by DMC and by other approaches. In the present study we explore an approach for evaluating these integrals as well as insights that can be gained from this information. As noted, DMC can be used to evaluate the ground state energy and wave function. Because the DMC wave function,

, is expressed in a basis of localized functions (called

walkers), evaluating integrals that require

2

is not entirely straightforward. Calculating

integrals that involve products of two di↵erent wave functions is even more challenging. Several approaches have been developed and explored for evaluating expectation values of multiplicative operators. 21–24 These fall into two general categories. The first is based on finite field methods from electronic structure theory along with ideas from perturbation theory. In this approach, a small perturbation, which is proportional to the operator for which the expectation value is desired, is added to the Hamiltonian being solved. 25,26 The finite field method, while accurate, can be expensive, as separate calculations must be performed for each operator of interest. Alternatively, since the ensemble of walkers provides a Monte Carlo sampling of the ground state wave function, if we could evaluate the value of

at the

locations of each of the walkers, we would have all the information that is needed to obtain expectation values. Barnett et al. 27 showed that this information can be obtained by propagating the wave function forward in time and equating the relative number of descendants of each walker to the value of

at the locations of these walkers. This approach is often

referred to as descendant weighting. Hornik et al. 28 and Barnett et al. 29 showed that one can recast the evaluation of an o↵diagonal matrix element as an expectation value, which can be obtained using the descendant weights. That work formed the basis of several studies of properties of electronic wave functions. 30 In the present study, we explore applications of this approach to the evaluation of overlap integrals between approximate vibrational wave functions and those obtained by DMC, as well as to the evaluation of matrix elements of the dipole moment operator. The latter application allows us to evaluate intensities based on DMC simulations.

4


Page 5 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The primary challenges in evaluating matrix elements using DMC come in normalizing the wave functions. This is particularly challenging when evaluating the overlaps between DMC and approximate wave functions, where an error that is no larger than 10

4

is desired.

When we use the approaches that have been developed for evaluating expectation values p to determine h |A| i/ h | ih | i, where the Monte Carlo-sampled wave function is ,

the evaluation of the numerator and the first term in the denominator can be obtained in a

relatively straightforward manner. To obtain h | i, we equate it to h |( / )2 | i, which can be evaluated using reciprocals of the normalized descendant weights. Because the values of the descendant weights are proportional to the corresponding value of

at the coordinates

of the walker, the contributions to h | i from wakers that are in regions of configuration space where

is small are plagued by numerical instabilities. In this study, we explore how

the accuracy of the evaluation of the calculated descendant weights a↵ects our ability to obtain o↵-diagonal matrix elements using DMC. In addition to using this approach to assess the accuracy of wave functions obtained by other approaches, the ability to obtain o↵-diagonal matrix elements allows us to evaluate intensities of vibrational transitions using DMC approaches. In an earlier study, 31 we showed that the fixed-node treatment of excited states using DMC 20 allowed us to obtain accurate energies and wave functions for the nine states of H3 O2 with one quantum of excitation in each of the vibrational degrees of freedom. In a subsequent study, we developed an approach for obtaining both frequencies and intensities of fundamental transitions based on the ground state probability amplitude obtained from DMC, and applied this approach to a study of the spectrum of H3 O2 . 32 Since H3 O2 is an anharmonic ion that contains both large and small amplitude OH vibrations, it provides an attractive system for evaluating the accuracy of new methods for calculating vibrational spectra. Specifically, H3 O2 consists of proton that is shared by two hydroxide ions. The potential energy along the proton transfer coordinate has two minima, which are separated by a 130 cm

1

barrier. 33 As a result, the ground

state wave function has a maximum in a geometry where the proton is equidistant from the

5



two oxygen atoms, rather than in a position that corresponds to one of the minimum energy structures. This results in the shared proton stretch being very anharmonic, with a measured fundamental frequency of 697 cm

1

and a large intensity. 34 In contrast, the frequency of the

out-of-phase vibration of the hydroxide OH stretches is 3653 cm 1 , which is much closer to the frequency of the OH stretch in water. 35,36 In addition, the fundamental transition involving the hydroxide OH stretches is much less intense than the shared proton stretch. In the second part of this study, we evaluate overlap integrals when both wave functions have been obtained using DMC, and use the results of these calculations to evaluate the intensities of the fundamentals in these two OH stretch vibrations in H3 O2 . The remainder of this paper is organized as follows. In the following section, we review the DMC approach. We also describe the extensions that will be used to evaluate expectation values and matrix elements of arbitrary multiplicative operators. We explore the dependence of the accuracy of this approach on the closeness of

to the ground state solution to the

Hamiltonian through a study of a Morse oscillator with varying degrees of anharmonicity and then for a study of formaldehyde and partially deuterated H+ 3 obtained using several approaches to obtain

. 37,38 Finally, we investigate the use of DMC for evaluating intensities

of the OH stretch transitions in the H3 O2 molecular ion, where both the ground and excited state wave functions are obtained using DMC.

Evaluating Matrix Elements with Di↵usion Monte Carlo The details of the di↵usion Monte Carlo (DMC) approach and our implementation are provided in the supporting materials. In DMC, an ensemble of localized functions, which we will refer to as walkers, are propagated. At the end of the simulation, the ensemble of walkers provides a Monte Carlo sampling of the function of interest, while the average of Eref in Eq. S8 or S13 over ⌧ provides a measure of the zero-point energy. In this study, two implementations of DMC are considered. In the first, we propagate

6


Page 6 of 31

Page 7 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the ensemble of walkers based on the imaginary-time time-dependent Schr¨odinger equation, and the resulting ensemble samples the ground state wave function of the system of interest, (x). The second approach introduces a trial wave function,

T (x),

and the ensemble of

walkers provides a Monte Carlo sampling of

f (x) = (x)

T (x)

(1)

Based on this, if the DMC simulation samples f (x), we can evaluate integrals involving the product of this function and an arbitrary function, g(x) through Monte Carlo integration, as Z

f (x)g(x)dx ⇡

P

(f )

j

wj (⌧0 )g(xj (⌧0 )) W (⌧0 )

(2)

(f )

where xj (⌧0 ) and wj (⌧0 ) represent the coordinates and weight of the jth walker at a selected time in the simulation, ⌧0 . The superscript (f ) is used to indicate that the object that is being represented by the ensemble of walkers is f (x). Finally, once the simulation has equilibrated the integral in Eq. 2 does not depend on the choice of ⌧0 . Consequently, we have dropped the ⌧ -dependence of f and

in the discussion that follows.

Often other quantities are desired, for example expectation values of operators, which require the ability to evaluate | (x)|2 . When we use the Monte Carlo sampling based on f (x), the summation on the right hand side of Eq. 2 provides a way to evaluate integrals involving (x)

T (x).

Obtaining expectation values of multiplicative operators, also requires

a method for evaluating (x)/

T (x)

at the positions of the walkers. Barnett et al. and Suhm

and Watts have shown 27,39 that at the coordinates of the jth walker at a time ⌧0 , xj (⌧0 ), (x)/

T (x)

is proportional to the number of descendants this walker has after a given

number of time steps, (f )

wj (⌧ 0 ) (xj (⌧0 )) (f ) = dj (⌧0 ) = (f ) T (xj (⌧0 )) wj (⌧0 ) 7


(3)


(f )

Page 8 of 31

(f )

Here wj (⌧0 ) represents the weight associated with the jth walker at ⌧0 , and wj (⌧ 0 ) repre(f )

sents the weight of this walker at a later time. When discrete weights are used wj (⌧0 ) = 1, (f )

while wj (⌧ 0 ) represents the number of walkers at ⌧ 0 that can be traced to the jth walker at ⌧0 . For the discussion that follows, ⌧avg = ⌧ 0

⌧0 provides the amount of time that the

system is propagated to generate the descendant weights. Since at ⌧0 , f (x) is represented by an ensemble of localized functions centered at xj (⌧0 ) (f )

with an associated weight wj (⌧0 ),

R

A(x)| (x)|2 dx R | (x)|2 dx P (f ) (f ) j dj (⌧0 )wj (⌧0 )A(xj (⌧0 )) = P (f ) (f ) j dj (⌧0 )wj (⌧0 ) P (f ) 0 j wj (⌧ )A(xj (⌧0 )) = P (f ) 0 j wj (⌧ )

hAi =

(4)

for an arbitrary multiplicative operator, A.

As we consider extending this approach to evaluating matrix elements, the evaluation is complicated by the fact that we must obtain the normalization for both of the wave functions involved in the integral as R

X(x)A(x(⌧0 )) (x)dx hX|A| i = q R R ( | (x)|2 dx)( |X(x)|2 dx)

If X is an analytically determined wave function, which we require to be the same as

(5) T,

Eq. 5 becomes

R

P

(f )

wj (⌧0 )A(xj (⌧0 )) qR =s ◆ R ⌘ ✓P ⇣ ⌘2 ⇣P ( | (x)|2 dx)( | T (x)|2 dx) (f ) 0 (f ) (f ) 0 /wj (⌧ ) j wj (⌧ ) j wj (⌧0 ) T (x)A(x)

(x)dx

8

j


(6)

Page 9 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The other case we will consider involves evaluating matrix elements between two states, both of which are determined using di↵usion Monte Carlo. The wave functions for the two states of interest are

(X)

(x) and X(x), with associated weights wj

( )

and wj . As such, in

this case

R

X(x)A(x) (x)dx qR =s R ⇣P 2 2 ( | (x)| dx)( |X(x)| dx)

P

(X)

wj (⌧ 0 )A(xj (⌧0 )) ◆ ⌘ ✓P ⇣ ⌘2 ( ) 0 (X) 0 ( ) 0 /wj (⌧ ) j wj (⌧ ) j wj (⌧ ) j

(7)

To obtain the descendant weights for two states based on the same set of walkers, we first equilibrate the ensemble based on a ground state simulation. Starting from a distribution of walkers at ⌧0 , the ground state descendant weights are evaluated by continuing this propagation until ⌧ 0 , and evaluating the number of descendants each walker generates over ⌧avg = ⌧ 0

⌧0 . The descendant weights for the excited state are obtained by propagating the

same ensemble of walkers forward in time using the fixed-node approach to obtain the state of interest and evaluating the number of descendants. The details of the fixed-node procedure, and its application to H3 O2 can be found elsewhere. 20,31 In a fixed-node calculation, the potential is replaced by an infinite potentialon one side of the nodal surface that divides the wave function into regions of positive and negative amplitudes. This ensures that any walker that crosses the node is removed from this simulation. Due to this additional mechanism for removing walkers in the excited fixed-node DMC simulations, more walkers will have small ( )

weights for the excited state than have small weights for the ground state. Since 1/wj (⌧ 0 ) (X)

is needed for the calculations of the intensity and not 1/wj (⌧ 0 ), the calculations will be more stable if

(x) in Eq. 7 represents the ground state, while X(x) represents an excited

state. As we consider Eqs. 6 and 7 we note that both contain a term that is proportional to ⌘ 1 ⇣ ( ) ( ) the wj (⌧ 0 ) , where represents f or . When wj (⌧ 0 ) is small, calculations of the overlaps can su↵er from numerical instabilities. To improve the accuracy of the calculations, 9



Page 10 of 31

rather than performing the calculation of the weights a single time, as is often done when expectation values are evaluated, weights are obtained several times and these values are ( )

averaged. In this way, wj (⌧ 0 ) in Eqs. 6 and 7 is replaced by

( ) wj (⌧ 0 )

=

Pnavg

( )

wm,j (⌧ 0 ) navg

m=1

(8) ( )

where navg reflects the number of evaluations of the descendant weights, wm,j (⌧ 0 ), that are at ⌧0 . Each of the evaluations of wm,j (⌧ 0 ) uses a

obtained using a single snapshot of f or

di↵erent series of random numbers in the DMC propagation, and this procedure will result ( )

in navg distinct sets of values of the wj (⌧ 0 ). In the discussion that follows, we will consider this averaging process and its influence on the accuracy of the final results. Finally, as noted ( )

above, the value of wj (⌧ 0 ) should not depend on the choice of ⌧0 and only the amount of time over which the descendants are collected, ⌧avg . For this reason, in the discussion that ( )

follows, we replace the ⌧ 0 argument of wj

with ⌧avg .

Results and Discussion Overlap Integrals for One-Dimensional Systems To start, we focus on one-dimensional systems to better understand the convergence properties of the overlap integrals with respect to various parameters. These include the anharmonicity of the system, the amount of time the descendant weights are collected, ⌧avg , and the number of independent evaluations of the descendant weights that are used in the evaluation of the average weights, navg . We focus on the one-dimensional Morse oscillator, described above. For these calculations, we use the harmonic ground state wave function for the trial wave function,

T,

and evaluate the normalized overlap integral between

and the ground state wave function obtained from DMC, represented by

10


T (x)

(x). This overlap, which will be

Page 11 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


O = h

T|

= q R

i R

= r⇣ P

(x)) dx R 2 2 (x)dx T (x)dx

(

T (x)

W (⌧0 ) ⌘ ⇣P ⌘ 1 w (⌧ ) w (⌧ ) avg j avg j j j

(9)

is simply the normalized overlap between the unnormalized harmonic and anharmonic ground state solutions for this model system. For comparison, we also obtained the ground state wave function for the Morse oscillator described above variationally. Based on this evaluation O(Exact) = 0.9963. As seen in Eqs. 5 and 6, calculating O requires the evaluation of three integrals. The R numerator, T (x) (x)dx, is trivially equated to the number of walkers that are used to represent

(x) at ⌧0 since discrete weighting is used for this part of the propagation. In

Figure 1, we explore the convergence properties of the terms in the denominator of O, scaled by the number of walkers at ⌧0 , with green squares/dashed line and purple diamonds/dashed line, as well as O itself, which is shown with the filled red squares/solid line. These three quantities are shown as functions of navg . To illustrate the convergence behavior, each of these integrals is subtracted from its value when navg = 3500, and the results are plotted on a log10 scale. Based on the results shown in Figure 1, the value of O shows fairly slow convergence behavior. Breaking O down into its three contributions, the slow convergence R 2 can be traced to the convergence properties of the evaluation of T (x)dx. The relaR 2 (x)dx to the value of navg (for this integral tive insensitivity to the calculated value of |f (3500)

f (navg )| < 10

6

independent of the value of navg ) is consistent with earlier obser(f )

vations that a single evaluation of the descendant weights, wj (⌧avg ), is generally sufficient for calculations of expectation values. The reason for the higher sensitivity that is found ⇣ ⌘ 1 R 2 (f ) (x)dx comes from its dependence on w (⌧ ) rather than for the evaluation of avg T j 11



(f )

Page 12 of 31

(f )

wj (⌧avg ). In geometries where wj (⌧avg ) is small, a small absolute error can manifest itself ⌘ 1 ⇣ (f ) . as a large error in wj (⌧avg )

The sensitivity of the accuracy of the calculated overlap integrals to the size of navg is

further explored in Figure 2. In panel (a), the value of O is subtracted from the value obtained variationally. In these plots, the red circles/solid line provide the same results as were plotted in Figure 1, subtracted from O(Exact) rather than O(navg = 3500). The blue diamonds/dashed line show the results obtained when we do not employ importance sampling for the evaluation of O. As is seen by comparing the magnitude of the errors represented by these two curves, the quality of the results is greatly improved when importance sampling is used. On the right side of Figure 2, the values of the average weights wj (⌧avg ) for 1000 individual walkers are plotted as functions of the corresponding xj (⌧0 ) values for the three simulation conditions identified by arrows in Figure 2(a). The purple curves in panels (ii) and (iii) provide the expected values of

(x)/

T (x).

When importance sampling is not used

T (x)

is a constant, and the purple curve in panel (i) provides the ground state wave function for the Morse potential. When importance sampling is used, the quality of the results improves with navg , and the average weights lie close to the purple curve when navg = 1120. On the other hand, using the same value of navg without importance sampling leads to a much greater spread in the values of the wj (⌧avg )’s when compared to the value of the ground state wave function at the same location. It is for this reason that we will focus on the results obtained using importance sampling in the remainder of the discussion. Before continuing, it should be noted that while the individual evaluations of the average weights do not generally lie along the purple curve, as is seen in Figure S1, their average values reproduce 2 (x) for this R 2 dx in Figure 1 and the use of descendant Morse oscillator. This is why the evaluation of

weighting to obtain expectation values of multiplicative operators are not as sensitive to the accuracy of the individual weights as the overlap integrals are found to be. Finally, the red curve in Figure 2(a) has a minimum when navg ⇡ 1600. This reflects the 12


Page 13 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


fact that when the weights are evaluated by propagating for ⌧avg = 8000 a.u., the sign of O(navg )

O(Exact) changes at this value of navg . This is seen in Figure 2(b). Looking at

the dependence of O on the amount of time the DMC wave function is propagated to obtain the descendant weights, we find it to decrease, although the rate of this decrease becomes smaller for large values of navg . This reflects a buildup of the weight onto a small number of walkers when the descendant weighting is propagated to long times. This problem is even more dramatic when importance sampling is not employed. A propagation time of 8000 a.u. is chosen for the remainder of this part of the study as it is in the range of values of ⌧avg over which the calculated value of O in Eq. 9 is numerically stable. 18,40 We have also explored the time step dependence on the accuracy of the evaluation of O. The results are plotted in Figure S2. Overall, convergence behavior of the overlaps with navg shown in Figure 1 does not depend on the size of the time step. While in the limit of large navg there is some dependence of the asymptotic value of O on the value of the time step, for time steps ranging from one to ten a.u., the error in the asymptotic value is on the order of 10

4

to 10 5 .

The final quantity that is explored is the e↵ect of the quality of the trial wave function as an approximation to the ground state wave function on the accuracy of the evaluation of both the overlap and the zero-point energy. In particular, errors can be introduced if the trial wave function does not sample the range of configurations that are sampled by the true ground state. The results of this analysis are provided in Table 1. For this part of the study, we use the Morse oscillator described above. Adjusting the reduced mass has the e↵ect of tuning the anharmonicity of the system, with smaller values of µ corresponding to the ground state wave function sampling more anharmonic regions of the potential. As is seen in the results reported in Table 1, when the mass is decreased, the overlap decreases, and the error in the calculated zero-point energy increases. For overlaps exceeding roughly 0.995, the energies obtained using DMC with and without importance sampling di↵er by less than 1 cm 1 . When the overlap becomes smaller, the error in both the overlap and the

13



Page 14 of 31

energy obtained when importance sampling is employed increase, with the calculated overlap consistently providing an underestimate of the expected value. The reason for these trends can be found in panels (ii) and (iii) on the right side of Figure 2. When importance sampling is used, the walkers are localized in the region of configuration space where the trial wave function has amplitude. Regions where

T

is small

will be poorly sampled. 41–43 It is this loss of walker density in the tails of the trial wave function that leads to the underestimation of O and the errors in the energy obtained from a DMC simulation with importance sampling. As a result, the approach described in the previous section will provide a more accurate measure of the overlaps as they approach unity, while the calculated overlaps will become more qualitative when these values are smaller. The quality of the calculated value of O can be correlated to the agreement between the zeropoint energy obtained with and without importance sampling. In fact for the vibrational problems with nodeless ground states it is encouraging how well the simple unguided DMC calculation is able to reproduce the exact zero-point energy.

Overlap Integrals for Formaldehyde and H+ 3 While it is interesting to explore the accuracy of DMC for evaluating overlap integrals for model systems, the focus of this work is to use this approach as a tool for assessing the quality of approximate wave functions for polyatomic systems. Here we focus on H2 CO, a modestly anharmonic molecule, as well as H2 D+ and D2 H+ , which sample much more anharmonic regions of the potential in their ground states than H2 CO. For all three molecules, we calculate overlaps between the DMC ground states and wave functions that were evaluated using recently developed VSCF and VSCF/CI approaches, 37,38 based on the potential surface of Yacmenev et al for H2 CO 44 and the potential of Aguado et al for H2 D+ and D2 H+ . 45 We start by considering H2 CO, which is the least anharmonic of the three molecules being considered. A DMC calculation using the VSCF wave function for

T

provides an energy

that is nearly identical to those obtained using DMC without importance sampling and the 14


Page 15 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


reported ground state energy based on the potential surface used for the calculation. 37,44 In this case h

T|

i= 0.9996(2). Based on the analysis provided above, this result is expected

to be quantitatively accurate. Also, as is seen in the results reported in Table 2, h

T |H|

Ti

deviates from the literature value by 4 cm 1 , although when it is used as a guiding function for DMC, the zero-point energy obtained by the two approaches di↵er by less than 1 cm 1 , which is the statistical uncertainty of these results. While the ground state wave function for H2 CO is expected to be reasonably welldescribed by harmonic treatments, H2 D+ and D2 H+ are considerably more anharmonic. This is reflected in the larger deviation between the values of h

T |H|

Ti

obtained from har-

monic and the VSCF or VCI(1) wave functions provided in Table 2. Despite the increased anharmonicity the zero-point energies for H2 D+ and D2 H+ , calculated using all three trial wave functions, di↵er by less than 1.5 cm 1 . Additionally, energies evaluated using the normal mode wave functions based on linear combinations of internal displacement coordinates di↵er from the exact results for this potential by less than 0.1 cm 1 . Examining the calculated overlaps, we find that the overlaps between the VCI(1) wave trial wave function and the exact wave function obtained using DMC exceed 0.995. The overlaps are somewhat smaller when the harmonic ground state wave functions are used. When we compare the overlaps obtained using harmonic wave functions based on normal modes, which have been constructed as linear combinations of displacements of Cartesian and internal coordinates, we find that the overlap is larger when the harmonic trial wave function is based on Cartesian displacements. This was at first surprising as it is expected that normal modes based on internal displacement coordinates should provide a better zero-order description of bending vibrations in such anharmonic systems. 47,48 To better understand these trends, we also evaluated hHi for these trial wave functions numerically, and the results are reported in Table 2. Consistent with the above results, we find that hHi is smaller when the normal modes are expressed as linear combinations of Cartesian displacements. Based on this, we conclude that the vibrations in H2 D+ (and D2 H+ ) are unusual. In this ion, the natural vibrations

15



are somewhere between those of a triatomic molecule, where the vibrations are described as two stretching motions and a bend, and a complex of H2 with D+ where the bend is better described as a H-H stretch, the symmetric stretch becomes the D-H2 stretch and the asymmetric stretch is now the hindered rotation of H2 in the presence of D+ . As noted above, the introduction of a guiding function can introduce errors if the exact wave function has amplitude in regions of configuration space where the guiding function is vanishingly small. 42 This is why the forms of the guiding functions often focus on the long-range behavior (e.g. the use of Jastrow factors), 49–51 or good guesses to the exact wave function based on prior variational Monte Carlo calculations. 15,41,52 This deterioration in the accuracy of the DMC results as |O| in Eq. 9 decreases will also lead to challenges as we use DMC to assess the accuracy of the approximate ground state wave functions that are also being used as the guiding functions. As noted above, we can use the errors in the zero-point energy as an additional measure of the robustness of the calculation of the overlaps. Based on the good agreement between the energies for these molecular systems, obtained with and without importance sampling, we expect that the reported overlaps are accurate within the uncertainties of the DMC results.

Calculating O↵ Diagonal Matrix Elements We turn now to the application of this methodology to the calculation of matrix elements involving two vibrational states in order to determine matrix elements of the dipole moment operator, which are used to evaluate intensities of vibrational transitions. We first explore the accuracy of this approach using a one dimensional harmonic oscillator and then consider the molecular system, H3 O2 . To determine the accuracy and convergence behavior of the evaluation of matrix elements, the expectation values and o↵ diagonal matrix elements of the position operators, q, and q 2 , for a harmonic oscillator are used. The harmonic oscillator that is used for this part of the study has a mass 3.5 amu and a harmonic frequency of 913.65 cm 1 . These values are the 16


Page 16 of 31

Page 17 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


same as those used to define the Morse oscillator in the above analysis. To validate the use of sampling of geometries based on the ground state probability amplitude, hq 2 i for the first vibrational excited state is compared with

h

1 |q

where the superscript (

2

|

0)

1i

(

0)

=

*

=

P

0

✓

q2

1 (q) 0 (q)

◆2

0

+

⌘2 ⇣ ( 1) ( ) 2 q w (⌧ ) /wj 0 (⌧avg ) avg j j j ⌘2 P ⇣ ( 1) ( ) (⌧avg ) /wj 0 (⌧avg ) j wj

(10)

on the matrix element in the left side of the equation indicates

that the sampled walkers are based on a ground state calculation. In contrast, the descendent weights are obtained using the potential that is appropriate for a fixed-node calculation for the first excited state, indicated by the (

1)

superscript on wj (⌧avg ). The results for the

harmonic treatments are provided in Tables S2 and S3. The results obtained for navg = 200 and 2000 deviate from the expected value by 3% when navg = 200 and are within 0.5% of the exact value when navg is increased to 2000. To explore the accuracy of the approach applied to an anharmonic molecular system, vibrations of H3 O2 are considered using the potential surface developed by Huang et al. 53 In these calculations, q represents the asymmetric stretch of the flanking OH bonds or the shared proton stretch. The evaluations of h

0 |q|

1i

(

0)

, h

1 |q

2

|

1i

(

0)

are plotted as

functions of 1/navg for the asymmetric stretch and shared proton displacement in Figures S3 and S4, respectively. The value of h

1 |q

2

|

1i

(

1)

included as red X’s in both of these plots

for comparison. The numerical data that corresponds to these plots is provided in Table S4. As is seen, the results obtained when navg = 100 and 200 are very close to the corresponding exact values. With the approach validated, we turn our attention to matrix elements of the dipole moment, and the intensities obtained from this evaluation. The results of these calculations

17



Page 18 of 31

are provided in Table S6 and plotted in Figure S5. The results are also summarized in Table 3. In this work, the fixed node results refer to those obtained using the approach described above, where the intensities are either obtained using navg = 200 in the first row or the extrapolated value based on the fits shown in Figure S5. The corresponding energies are obtained from the converged fixed-node calculation used to obtain the h

1 |q

2

|

1i

(

1)

values

shown with red X’s in Figures S3 and S4. While the values of the intensities di↵er for the two approaches, their ratio remains nearly constant. When navg = 2000, the results are in good agreement with the extrapolated values. The intensity for the shared proton displacement is in reasonable agreement with our previous calculations in which we evaluated spectra based on ground state probability amplitudes. 32 Larger di↵erences are found for the intensity of the asymmetric stretch of the flanking OH bonds, which appear to be borrowing intensity from the shared proton stretch to a greater extent than is anticipated by the localized excitation model implicit in the results reported in Ref. 32. This larger intensity in the flanking OH stretch is consistent with the analogous protonated water cluster, H5 O+ 2 , where harmonic calculations place the intensity of the OH stretches of the outer water molecules an order of magnitude lower than the intensity of the shared proton stretch.

Conclusions In this study, we have explored approaches for obtaining matrix elements involving two di↵erent wave functions using di↵usion Monte Carlo. This enables us to evaluate the accuracy of wave functions obtained using other approximate approaches. We find that when the approximate wave function is used as a trial wave function, the combination of the accuracy of the calculated energy and the value of the overlap integral, h |

T i,

allow us to assess

the accuracy of the wave function. In cases where the wave function is of high quality, such evaluations of the overlap integral may be made with accuracy of 10 4 . For less accurate wave

18


Page 19 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


functions, the accuracy of both the zero-point energy and the wave function deteriorate. We have also shown how this approach may be applied to the evaluation of matrix elements of the dipole moment operator, which are used in the calculation of the intensities. High accuracy is achieved when the approach is tested for one-dimensional systems, and the approach is also demonstrated for a highly anharmonic molecular system.

Acknowledgement Support from the Chemistry Division of the National Science Foundation (CHE-1619660) is gratefully acknowledged. Parts of this work were performed using the Ilahie cluster, which was purchased using funds from a MRI grant from the National Science Foundation (CHE-1624430), and using Colby College’s NSCC cluster and supported by the NSF CC* grant (OAC-1659142). Additionally, this work was facilitated through the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system and funded by the STF at the University of Washington. We also thank I. W. Bulik and M. J. Frisch for providing us with the wave functions used to evaluate the VSCF wave function for H2 CO and the VCI(1) wave function for partially deuterated forms of H+ 3.

Supporting Information Available A description of Di↵usion Monte Carlo and our implementation; Comparisons of descendant weights and the ground state wave function for a Morse oscillator; Convergence properties of O for the Morse oscillator; Convergence properties of matrix elements of the displacements of the shared proton stretch and antisymmetric OH stretch of H3 O2 ; Convergence properties of the intensities of fundamental transitions in the shared proton and antisymmetric OH stretches in H3 O2 ; A table summarizing parameters used in the DMC calculation; Tables providing results for calculations of matrix elements obtained for the ground and first excited state of a harmonic oscillator; Tables providing the data used to generate figures for H3 O2 . 19



References (1) Wilson, E. B.; Decius, J. C.; Cross, P. C. Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra; Courier Corporation, 1955. (2) Scott, A. P.; Radom, L. Harmonic Vibrational Frequencies: An Evaluation of HartreeFock, Møller-Plesset, Quadratic Configuration Interaction, Density Functional Theory, and Semiempirical Scale Factors. J. Phys. Chem. 1996, 100, 16502–16513. (3) Whitehead, R.; Handy, N. Variational Calculation of Vibration-Rotation Energy Levels for Triatomic Molecules. J. Mol. Spectrosc. 1975, 55, 356–373. (4) Bowman, J. M. Self-Consistent Field Energies and Wavefunctions for Coupled Oscillators. J. Chem. Phys. 1978, 68, 608–610. (5) Willetts, A.; Handy, N. C.; Green, W. H.; Jayatilaka, D. Anharmonic Corrections to Vibrational Transition Intensities. J. Phys. Chem. 1990, 94, 5608–5616. (6) Barone, V. Anharmonic Vibrational Properties by a Fully Automated Second-Order Perturbative Approach. J. Chem. Phys. 2005, 122, 014108. (7) Bowman, J. M.; Carter, S.; Huang, X. MULTIMODE: A Code to Calculate Rovibrational Energies of Polyatomic Molecules. Int. Rev. Phys. Chem. 2003, 22, 533–549. (8) Christiansen, O. Coupled Cluster Theory with Emphasis on Selected New Developments. Theor. Chem. Acc 2006, 116, 106–123. (9) Christiansen, O. Vibrational Structure Theory: New Vibrational Wave Function Methods for Calculation of Anharmonic Vibrational Energies and Vibrational Contributions to Molecular Properties. Phys. Chem. Chem. Phys. 2007, 9, 2942–2953. (10) Bowman, J. M.; Carrington, T.; Meyer, H.-D. Variational Quantum Approaches for Computing Vibrational Energies of Polyatomic Molecules. Mol. Phys. 2008, 106, 2145– 2182. 20


Page 20 of 31

Page 21 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(11) Roy, T. K.; Gerber, R. B. Vibrational Self-Consistent Field Calculations for Spectroscopy of Biological Molecules: New Algorithmic Developments and Applications. Phys. Chem. Chem. Phys. 2013, 15, 9468–9492. (12) Bloino, J.; Baiardi, A.; Biczysko, M. Aiming at an Accurate Prediction of Vibrational and Electronic Spectra for Medium-to-Large Molecules: An Overview. Int. J. Quantum Chem. 2016, 116, 1543–1574. (13) Born, M.; Oppenheimer, R. Zur Quantentheorie der Molekeln. Ann. Phys. 1927, 389, 457–484. (14) Podolsky, B. Quantum-Mechanically Correct Form of Hamiltonian Function for Conservative Systems. Phys. Rev. 1928, 32, 812. (15) Austin, B. M.; Zubarev, D. Y.; Lester Jr, W. A. Quantum Monte Carlo and Related Approaches. Chem. Rev. 2011, 112, 263–288. (16) McCoy, A. B. Di↵usion Monte Carlo Approaches for Investigating the Structure and Vibrational Spectra of Fluxional Systems. Int. Rev. Phys. Chem. 2006, 25, 77–107. (17) Umrigar, C.; Nightingale, M.; Runge, K. A Di↵usion Monte Carlo Algorithm with Very Small Time-Step Errors. J. Chem. Phys. 1993, 99, 2865–2890. (18) Hammond, B. L.; Lester, W. A.; Reynolds, P. J. Monte Carlo Methods in Ab Initio Quantum Chemistry; World Scientific, 1994; Vol. 1. (19) Anderson, J. B. A Random-Walk Simulation of the Schr¨odinger Equation: H+ 3 . J. Chem. Phys. 1975, 63, 1499–1503. 1 0 3 (20) Anderson, J. B. Quantum Chemistry by Random Walk. H 2 P, H+ 3 D3h A1 , H2 ⌃u , 1 H4 1 ⌃ + g , Be S. J. Chem. Phys. 1976, 65, 4121–4127.

(21) Liu, K.; Kalos, M.; Chester, G. Quantum Hard Spheres in a Channel. Phys. Rev. A 1974, 10, 303. 21



(22) Reynolds, P.; Barnett, R.; Hammond, B.; Grimes, R.; Lester Jr, W. Quantum Chemistry by Quantum Monte Carlo: Beyond Ground-State Energy Calculations. Int. J. Quantum Chem. 1986, 29, 589–596. (23) Langfelder, P.; Rothstein, S. M.; Vrbik, J. Di↵usion Quantum Monte Carlo Calculation of Nondi↵erential Properties for Atomic Ground States. J. Chem. Phys. 1997, 107, 8525–8535. (24) Warren, G. L.; Hinde, R. J. Population Size Bias in Descendant-Weighted Di↵usion Quantum Monte Carlo Simulations. Phys. Rev. E 2006, 73, 056706. (25) Sandler, P.; Buch, V.; Clary, D. Calculation of Expectation Values of Molecular Systems Using Di↵usion Monte Carlo in Conjunction with the Finite Field Method. J. Chem. Phys. 1994, 101, 6353–6355. (26) Lee, H.-S.; Herbert, J. M.; McCoy, A. B. Adiabatic Di↵usion Monte Carlo Approaches for Studies of Ground and Excited State Properties of Van der Waals Complexes. J. Chem. Phys 1999, 110, 5481–5484. (27) Barnett, R.; Reynolds, P.; Lester, W. Monte Carlo Algorithms for Expectation Values of Coordinate Operators. J. Comput. Phys. 1991, 96, 258–276. (28) Hornik, M.; Snajdr, M.; Rothstein, S. M. Estimating the Overlap of an Approximate with the Exact Wave Function by Quantum Monte Carlo Methods. J. Chem. Phys. 2000, 113, 3496–3498. (29) Barnett, R.; Reynolds, P.; Lester Jr, W. Computation of Transition Dipole Moments by Monte Carlo. J. Chem. Phys. 1992, 96, 2141–2154. (30) Barnett, R.; Johnson, E.; Lester Jr, W. Quantum Monte Carlo Determination of the Lithium 2 2 S ! 2 P Oscillator Strength: Higher Precision. Phys. Rev. A 1995, 51, 2049. 22


Page 22 of 31

Page 23 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(31) McCoy, A. B.; Huang, X.; Carter, S.; Bowman, J. M. Quantum Studies of the Vibrations in H3 O2 and D3 O2 . J. Chem. Phys. 2005, 123, 064317. (32) McCoy, A. B.; Diken, E. G.; Johnson, M. A. Generating Spectra from Ground-State Wave Functions: Unraveling Anharmonic E↵ects in the OH ·H2 O Vibrational Predissociation Spectrum. J. Phys. Chem. A 2009, 113, 7346–7352. (33) Diken, E. G.; Headrick, J. M.; Roscioli, J. R.; Bopp, J. C.; Johnson, M. A.; McCoy, A. B. Fundamental Excitations of the Shared Proton in the H3 O2 and H5 O+ 2 Complexes. J. Phys. Chem. A 2005, 109, 1487–1490. (34) Roscioli, J. R.; Diken, E. G.; Johnson, M. A.; Horvath, S.; McCoy, A. B. Prying Apart a Water Molecule with Anionic H-Bonding: A Comparative Spectroscopic Study of the X ·H2 O (X = OH, O, F, Cl, and Br) Binary Complexes in the 600-3800 cm

1

Region.

J. Phys. Chem. A 2006, 110, 4943–4952. (35) Price, E. A.; Robertson, W. H.; Diken, E. G.; Weddle, G. H.; Johnson, M. A. Argon Predissociation Infrared Spectroscopy of the Hydroxide-Water Complex (OH ·H2 O). Chem. Phys. Lett. 2002, 366, 412–416. (36) Robertson, W. H.; Diken, E. G.; Price, E. A.; Shin, J.-W.; Johnson, M. A. Spectroscopic Determination of the OH Solvation Shell in the OH ·(H2 O)n Clusters. Science 2003, 299, 1367–1372. (37) Bulik, I. W.; Frisch, M. J.; Vaccaro, P. H. Vibrational Self-Consistent Field Theory Using Optimized Curvilinear Coordinates. J. Chem. Phys. 2017, 147, 044110. (38) Bulik, I. W.; Frisch, M. J.; Vaccaro, P. H. Fixed-Node, Importance-Sampling Di↵usion Monte Carlo for Vibrational Structure with Accurate and Compact Trial States. J. Chem. Theory Comput. 2018, 14, 1554–1563.

23



(39) Suhm, M. A.; Watts, R. O. Quantum Monte Carlo Studies of Vibrational States in Molecules and Clusters. Phys. Rep. 1991, 204, 293–329. (40) Blume, D. Fermionization of a Bosonic Gas Under Highly Elongated Confinement: A Di↵usion Quantum Monte Carlo Study. Phys. Rev. A 2002, 66, 053613. (41) Umrigar, C.; Wilson, K.; Wilkins, J. Optimized Trial Wave Functions for Quantum Monte Carlo Calculations. Phys. Rev. Lett. 1988, 60, 1719. (42) Bernu, B.; Ceperley, D.; Lester Jr, W. The Calculation of Excited States with Quantum Monte Carlo. II. Vibrational Excited States. J. Chem. Phys. 1990, 93, 552–561. (43) Bernu, B.; Ceperley, D.; Lester Jr, W. Erratum: The Calculation of Excited States with Quantum Monte Carlo. II. Vibrational Excited States [J. Chem. Phys. 93, 552 (1990)]. J. Chem. Phys. 1991, 95, 7782–7782. (44) Yachmenev, A.; Yurchenko, S. N.; Jensen, P.; Thiel, W. A New Spectroscopic Potential Energy Surface for Formaldehyde in its Ground Electronic State. J. Chem. Phys. 2011, 134, 244307. (45) 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 11 A0 Potential. J. Chem. Phys. 2000, 112, 1240–1254. (46) Miller, S.; Tennyson, J. First Principles Calculation of the Molecular Constants of H+ 3, H2 D+ , D2 H+ , and D+ 3 . J. Mol. Spectrosc. 1987, 126, 183–192. (47) Sibert III, E. L.; Hynes, J. T.; Reinhardt, W. P. Fermi Resonance from a Curvilinear Perspective. J. Phys. Chem. 1983, 87, 2032–2037. (48) McCoy, A. B. The Role of Electrical Anharmonicity in the Association Band in the Water Spectrum. J. Phys. Chem. B 2014, 118, 8286–8294.

24


Page 24 of 31

Page 25 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(49) Jastrow, R. Many-Body Problem with Strong Forces. Phys. Rev. 1955, 98, 1479. (50) Drummond, N.; Towler, M.; Needs, R. Jastrow Correlation Factor for Atoms, Molecules, and Solids. Phys. Rev. B 2004, 70, 235119. (51) Bouab¸ca, T.; Bra¨ıda, B.; Ca↵arel, M. Multi-Jastrow Trial Wavefunctions for Electronic Structure Calculations with Quantum Monte Carlo. J. Chem. Phys. 2010, 133, 044111. (52) Toulouse, J.; Umrigar, C. J. Optimization of Quantum Monte Carlo Wave Functions by Energy Minimization. J. Chem. Phys. 2007, 126, 084102. (53) Huang, X.; Braams, B. J.; Carter, S.; Bowman, J. M. Quantum Calculations of Vibrational Energies of H3 O2 on an Ab Initio Potential. J. Am. Chem. Soc. 2004, 126, 5042–5043.

25



FIGURES

Figure 1: The convergence of theR integrals that contribute to Eq. 9 scaled by the R number 2 2 fo walkers at time ⌧0 , f , for f = (x)dx (purple diamonds/dashed line), f = T (x)dx (green squares/dashed line) and f = O = h | T i (red filled circles/solid line) and plotted as functions of the number of independent evaluations of the weights (navg ), which is used R 2 (f ) (x)dx = to determine R 2 the wj (⌧avg ). The results are subtracted from f (3500), where 0.9999, T (x)dx = 1.0073 and h | T i = 0.9964. These calculations are based on the onedimensional Morse oscillator potential described in the text, T is the correspond harmonic ground state wave function, and weights are collected for ⌧avg =8000 a.u.

26


Page 26 of 31

Page 27 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(f )

Figure 2: (a) The convergence of O = h | T i (see Eq. 9), evaluated using wj (⌧avg ) values that are obtained from DMC simulations, which are propagated with importance sampling and ⌧avg = 8000 a.u. (red filled circles/solid line), and without (blue diamonds/dashed line) importance sampling and ⌧avg = 600 a.u. based on the Morse oscillator model described in the text. Panels (i)-(iii) provide the distribution of the average weights shown for 1000 walkers. These are obtained (i) without importance sampling and navg = 1120 and with importance sampling and (ii) navg = 4 and (iii) navg = 1120, as indicated by corresponding colored arrows in panel (a). The solid purple curves in these plots show the expected value of (i) (x) and (ii) and (iii) (x)/ T (x). (b) The dependence of the error in O on ⌧avg is plotted as a function of ⌧avg for specified values of navg .

27



Table 1: Ground State Properties for Morse Oscillators Using DMC With and Without Importance Sampling.a µb h Morse | HO i h Morse | HO ic E0 /cm 1 E0 /cm 1 E0 /cm 1 /amu Exact DMC,ISd DMC,IS DMCe Exact 28 0.9987 0.9987 (0.2) 161.17 (0.04) 161.17 (0.13) 161.16 14 0.9981 0.9982 (0.3) 227.77 (0.03) 227.66 (0.18) 227.72 7 0.9974 0.9974 (0.5) 321.70 (0.06) 321.87 (0.16) 321.63 3.5 0.9963 0.9963 (0.5) 454.37 (0.04) 454.13 (0.26) 454.04 1.725 0.9947 0.9947 (0.9) 641.28 (0.1) 640.31 (0.23) 640.48 0.8625 0.9925 0.9910 (1.6) 904.63 (0.04) 902.69 (0.32) 902.51 0.43125 0.9893 0.9802 (4.1) 1275.28 (0.11) 1269.81 (0.27) 1269.82 a Numbers in parentheses provide one standard deviation based on five calculations using statistically independent wave functions. b The reduced mass. c The standard deviations have been multiplied by 104 . d Evaluated using importance sampling with navg = 2100 and ⌧avg = 8000 a.u. e Evaluated using DMC without importance sampling.

28


Page 28 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecule

T

h

T |H|

ˆ T ib h | T ic E0 /cm 1 E0 /cm 1 /cm 1 DMC,ISd DMC,IS DMCe f CH2 O VCI(1) 5774.1 0.9996 (0.1) 5769.8 (0.1) 5770.0 (1.1) g N.M. Cart. 5854.9 0.9946 (1.5) 5770.9 (0.2) + h H2 D VCI(1) 4028.2 0.9967 (1.3) 3978.9 (0.1) 3978.4 (0.7) N.M. Cart.g 4217.6 0.9918 (3.3) 3978.4 (0.1) N.M. Inti 4259.5 0.9909 (3.2) 3979.9 (0.3) + h D2 H VCI(1) 3623.2 0.9957 (1.8) 3562.7 (0.2) 3561.7 (1.5) N.M. Cart.g 3751.3 0.9932 (2.9) 3562.0 (0.2) N.M. Inti 3816.8 0.9873 (5.1) 3563.2 (0.2) a Numbers in parentheses provide one standard deviation based on five calculations using statistically independent wave functions. b Average energy based on the T . The literature values for these zeropoint energies are 5769.78 cm 1 (H2 CO), 44 3979.905 cm 1 (H2 D+ ) and 3563.086 cm 1 (D2 H+ ). 46 c The standard deviation has been multiplied by 104 . d Evaluated with DMC using importance sampling with navg =560 and ⌧avg =150 a.u. e Evaluated with DMC without importance sampling. f Ref. 37. g Normal modes defined as a linear combination of Cartesian coordinates. h Ref. 38. i Normal modes defined as a linear combination of internal coordinates.

Table 2: Comparison of Overlaps and Zero-Point Energies For H2 CO, H2 D+ and D2 H+ .a

Page 29 of 31 The Journal of Physical Chemistry


29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Asymmetric Stretch Shared Proton Displacement Method E (cm 1 ) Intensity (km mol 1 ) E (cm 1 ) Intensity (km mol 1 ) Fixed Nodea 3607 (8) 50 (6) 649 (9) 580 (20) b c Fixed Node 95 (7) 940 (40)c d Ground State 3650 (20) 8 (5) 810 (40) 1900 (100) a Fixed node approximation described in the text with ⌧avg = 750 a.u.. ⌧ = 10 a.u., and navg = 200 averaged over calculates based on 10 statistically independent wave functions. b The extrapolated values as the limit of large navg , obtained from the fits to the data plotted in Figures S3 and S4 in the supporting information. c Extrapolated values are in good agreement with simulations that have longer equilibration times (⌧ =60 000 a.u.) and navg = 2000. Intensities of the asymmetric stretch and shared proton displacement are 80 (10) km mol 1 and 970 (40) km mol 1 , respectively. d Evaluated using the approach described in Ref. 32.

Table 3: Intensities for Fundamental Transitions in H3 O2 Evaluated Using Two Approaches



30

Page 30 of 31

Page 31 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC graphic:

31


Evaluation of Matrix Elements Using Diffusion Monte Carlo Wave

Recommend Documents