Fixed-Node, Importance-Sampling Diffusion Monte Carlo for

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

Fixed-Node, Importance-Sampling Diffusion Monte Carlo for Vibrational Structure with Accurate and Compact Trial States Ireneusz W. Bulik,*,† Michael J. Frisch,‡ and Patrick H. Vaccaro† †

Department of Chemistry, Yale University, New Haven, Connecticut 06511, United States Gaussian, Inc., 340 Quinnipiac St. Bldg. 40, Wallingford, Connecticut 06492, United States

‡

S Supporting Information *

ABSTRACT: We present a new approach to importance sampling in diffusion Monte−Carlo (DMC) simulations of vibrational excited states whereby the trial wave functions for low-energy states are incorporated into the diffusion equations so as to enforce their orthogonality. For the model systems examined here, simple variational wave functions based on the vibrational self-consistent field (VSCF) and the simplest vibrational configuration interaction (VCI) are effective in importance sampling provided that internal coordinates used in the underlying one-particle functions have been variationally optimized. The resulting model yields results comparable in accuracy to the best unguided DMC calculations without requiring an a priori choice of coordinates to specify nodal hyperplanes.

1. INTRODUCTION The prediction of complete vibrational spectra with high accuracy is an important but challenging problem. The harmonic approximation is adequate for assigning fundamentals above 1000 cm−1 in typical molecules, and improved accuracy in the energies and intensities for these features usually can be achieved using vibrational perturbation theory based on the normal-mode description.1−10 More accurate solutions within the normal-mode framework that entail Taylor-series expanded Hamiltonians have been developed using coupled-cluster theory11−13 and large-scale configuration interaction.14−22 However, when normal modes do not provide an adequate zeroth-order description, then practical models built upon such a coordinate basis become unreliable. Hindered rotations and low-frequency wagging modes of large molecules are the most common examples of these deficiencies, but systems with unusual bonding motifs, such as H5+ and H3+ (the latter of which is considered here), also can be problematic. Efforts to employ perturbation theory beyond the normalmode description have been reported,23,24 but most work using the full vibrational Hamiltonian has focused on variational treatments for small systems.25−37 These efforts have produced highly accurate results, but their computational cost has been shown to grow exponentially with system size. Less computationally demanding variational implementations of the vibrational self-consistent field (VSCF) ansatz, based upon either normal-mode or fixed-internal coordinates, have been reported.10,38−41 In previous work performed on several difficult cases, we have shown that a reasonable description can be produced by exploiting VSCF for the ground state along with minimal vibrational configuration interaction (VCI) for the excited states, as long as the model includes nonlinear © XXXX American Chemical Society

variational optimization of the internal coordinates used to produce the single-mode functions that form the product wave functions.42 The quality of predictions made by this model certainly could be improved through more extensive configuration interaction; however, since the complexity of these calculations increases exponentially with the number of modes, such an approach would limit the size of tractable molecular systems. In this paper, we examine the use of Diffusion Monte Carlo (DMC) methods to improve upon the accuracy of our previous results. DMC is an attractive alternative to large VCI because it scales much more reasonably in cost with system size. Furthermore, if VSCF and small VCI wave functions are qualitatively reasonable, then DMC using importance sampling based on them should be substantially more efficient than unguided DMC. Aside from lessening the computational load needed to achieve a given accuracy, the availability of reliable trial states should allow higher-lying overtone and combination bands, which are difficult or impossible to model otherwise, to be studied systematically. The application of DMC methods for the large manifold of states that must be considered in the vibrational problem presents new methodological issues which, to the best of our knowledge, have not been addressed previously in detail. We first present our approach to DMC using VSCF and VCI wave functions for importance sampling and describe how such trial wave functions can be exploited to apply DMC to excited states. Results obtained by bringing this formalism to bear on Received: January 5, 2018 Published: February 21, 2018 A

DOI: 10.1021/acs.jctc.8b00016 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation specific molecular systems will be presented and compared with those obtained by other methods.

, (τ ) =

2. IMPORTANCE-SAMPLING DMC FOR GROUND AND EXCITED STATES 2.1. Importance-Sampling DMC. Diffusion Monte Carlo is an established approach to the quantum many-body problem, hence the ensuing presentation is limited to defining terms and presenting those aspects which are new or specific to our implementation. In the importance sampling form of DMC, the nonnegative product f(τ,x⃗) of the DMC wave function ψ(τ,x⃗) and a fixed guiding wave function ψT(x⃗) are evolved in imaginary time according to the Hamiltonian of the system (atomic units used throughout):43−45 ∂f (τ , x ⃗) = ∂τ

N

N i

−(El(x ⃗) − ,(τ ))f (τ , x ⃗) (1)

where the Di = 1/(2mi) is the diffusion coefficient of ith atom with mass mi, F⃗i(x⃗) = 2[∇iψT(x⃗)/ψT(x⃗)] is usually referred to as the “quantum force,” El(x⃗) = [Ĥ ψT(x⃗)/ψT(x⃗)] is the local energy at x⃗, and ,(τ ) is an energy offset. For small time steps and with a second-order Trotter expansion, the appropriate Green’s function for time evolution takes the form:43,44,46 G(x ⃗ → x ′⃗ , Δτ ) ⎡ ⎛ E (x ′⃗ ) + El(x ⃗) ⎞⎤ = exp⎢ −Δτ ⎜ l − , (τ )⎟ ⎥ ⎝ ⎠⎦ 2 ⎣ × ⎛ N ⎡ (x ′ − x ⃗ − ΔτD F ⃗(x ⃗))2 ⎤⎞ i i i ⎜⎜∏ (4πDiΔτ )−3/2 exp⎢ − i⃗ ⎥⎟⎟ 4DiΔτ ⎣ ⎦⎠ ⎝ i (2)

The action of this operator on a function represented as a discrete set of Dirac delta functions (walkers) yields a set having the form of Gaussians multiplied by factors arising from the potential (branching) term. This set is sampled stochastically, yielding yet another discrete collection of delta functions, and the procedure is repeated. To be specific, at each time slice the coordinates of every atom in each walker is propagated as xi⃗′ = xi⃗ + ΔτDi Fi ⃗(x ⃗) + χi⃗

∑μ |wμ(τ )|

−α

∑μ |wμ(τ )| − W0 W0

(5)

where index μ runs over all walkers, El(xμ⃗ ) is the local energy at the position of the μth walker, W0 is the desired total weight (here, the number of walkers), and α is a population stabilizing parameter, which has been set equal to (Δτ)−1.46,51 The absolute value of wμ in eq 5 is necessary because we have adopted the convention of encoding the sign of the state in the sign of walkers. The energies reported in the present work were obtained by imaginary-time averaging of the local energy only and not from ,(τ ). With the Anderson correction, one can control the total weight; however, it also is necessary to avoid accumulation of the total weight in just a few walkers. Hence, whenever the weight of a walker falls under a predefined threshold (here, W−1 0 , suggested in ref 46 as a good compromise between efficiency of simulation and population control bias), the walker is annihilated and the walker with maximal weight is split in two.52 2.2. DMC for Excited States. The prescription outlined above is sufficient for the nodeless ground state of the vibrational problem, but additional constraints are needed for excited states. In particular, unconstrained propagation of the excited-state wave function according to the conventional DMC prescription will lead to its collapse into the nodeless ground-state distribution. In the case of Fermionic systems, techniques for handling the required nodal properties of the ground state are well established.44,45 Although the present work does not entail nodal surfaces arising from particleexchange symmetries, kindred methods must be developed to treat nodal-surface constraints imposed upon excited vibrational manifolds. The propagation of f(τ,x⃗) in the importance-sampling scheme implies that the nodal surface(s) of the DMC wave function are the same as those of the trial wave function. Here, each nodal surface is treated as an impenetrable barrier by eliminating any walker that attempts to change sign, thereby leading to a de facto independent ground-state simulation in each of the nodal pockets. If the nodal surface(s) were exact, then the simulation in every isolated nodal pocket converges to the exact energy of an excited state. However, with approximate nodal surface(s), the average energy in different pockets is not guaranteed to be the same,53,54 and eq 4 dictates collapse of all weight into a single nodal pocket if unconstrained propagation is allowed. The resulting single-pocket distribution describes a mixture of the desired state with lower-energy states. There have been previous approaches to handle analogous problems in electronic excited states.53−56 A comparison with previous work on vibrational excited states will be presented in section 4. Any importance-sampling DMC calculation uses the guiding wave functions to provide nodal surfaces and to accelerate convergence. The ensuing work also has exploited guiding wave functions to enforce approximate orthogonality to lower states during propagation. This prevents collapse of the wave function to a single pocket and eliminates mixing of the lower energy states (to an extent determined by the similarity of the guiding and exact wave functions of the energetically lower states). An estimate, Λk, of the overlap of the DMC excited state of interest, ψ(τ,x⃗) = f(τ,x⃗)/ψT(x⃗), with the kth variational trial state of lower energy, ψkT(x⃗), can be defined by

∑ Di∇i2 f (τ , x ⃗) − ∑ Di∇·i (f (τ , x ⃗)Fi ⃗(x ⃗)) i

∑μ |wμ(τ )|El(xμ⃗ )

(3)

where χ⃗i ∈ 3 is composed of normally distributed random numbers with zero mean and variance 2DiΔτ. The branching terms are accounted for via a continuous reweighting approach,46,47 whereby the weight wμ of the μth walker at position xμ⃗ is modified as follows: ⎡ ⎛ El(xμ′⃗ ) + El(xμ⃗ ) ⎞⎤ wμ(τ + Δτ ) ← wμ(τ ) exp⎢ −Δτ ⎜ − ,(τ )⎟⎥ ⎢⎣ 2 ⎝ ⎠⎥⎦ (4)

The energy offset ,(τ ) is taken to be the average of the local energy at the previous step plus the population-stabilizing term of Anderson:46,48−50 B


Article

Journal of Chemical Theory and Computation

Λ k (τ ) =

∑μ ∑μ

of the local energy by ∼8000 cm−1 for a reasonable Δτ value of 0.5 au. Similarly, to quench rapid diffusion due to the sampling of configurations next to nodes of the trial wave function, the forces on the ith atom are computed as

wμ(τ )

ψ k (x ⃗ ) |ψ (x )| T μ T

μ⃗

wμ(τ )

ψ k (x ⃗ ) |ψT(xμ⃗ )| T μ

∝ ⟨ψ |ψTk⟩(τ ) (6)

m ⎧ ⃗ |Fi ⃗| < i F ⎪ ⎪ i Δτ Fi ⃗ ← ⎨ ⎪ mi sgn(F ⃗) otherwise ⎪ i ⎩ Δτ

where Λk ∈ [−1,1] and the aforementioned convention that the weight of the walker encodes the sign of current guiding wave function has been adopted (necessitating the norm of ψT(x⃗) to be used in the above as the sign is carried by wμ(τ)). This metric is used to constrain the weights independently in a manner similar to the Anderson correction:57

In practice, the latter correction has been found to be less important. Unlike unguided simulations,51 a recrossing correction48−50 has not been used in our initial studies. This is in keeping with previous reports, which suggest it is less important when guiding states are implemented.46 Recrossing correction as outlined in ref 61 shall be explored in future studies. Additionally, the magnitude of ψT(x⃗) in eq 6 and eq 7 has been replaced by max(|ψT(x⃗)|, tol) where tol has been set equal to 10−5. This has been applied to the orthogonalization correction for each lower state. The regions where the wave function is small in magnitude not only contribute little to the overlap but also may be sampled poorly, hence increasing the noise in the simulation or leading to numerical instabilities. In the Supporting Information, we present tests which show that the results are insensitive to variations of this threshold by an order of magnitude. Finally, additional reweighting of each walker is performed according to

wμ(τ + Δτ ) ← ⎛ ⎛ wμ(τ ) ⎞⎞ wμ(τ ) ∏ exp⎜⎜ −λk Λk(τ ) sgn⎜⎜ ψTk(xμ⃗ )⎟⎟⎟⎟ ⎝ |ψT(xμ⃗ )| ⎠⎠ ⎝ k

(7)

where parameters λk take the place of α in eq 5. This ansatz prevents collapse of the wave function into region(s) of a single sign. In principle, the wave function could abandon some nodal pockets in favor of others of the same sign (via a population control step), but with reasonable trial wave functions this does not occur. We recently demonstrated that compact and accurate manymode states can be constructed via VSCF and short VCI expansions using the full vibrational Hamiltonian, without evoking constraints on the functional form of the potential or approximations to the kinetic energy.42,58 While a VSCF wave function is likely to have sufficient accuracy for the ground state, excited-state wave functions based on pure VSCF (i.e., adding one quantum to one mode to produce an excited-state wave function) cannot be expected to be adequate, since the constraint of the product-state formalism from the ground state is not relaxed. As such, our focus is on using the smallest expansion which mitigates this restriction as the trial wave functionVCI(1) in the case of fundamentals, VCI(2) in the case of two-quanta combinations and overtones, etc. Of course, the use of any trial wave function imposes the approximation of fixed nodes on the DMC calculation, and both the accuracy and speed of said analyses will depend on the quality of the trial wave function. 2.3. Additional Technical Details. Importance sampling may become problematic whenever walkers enter regions where the guiding wave function takes on very small values. In these cases, the local energy and the quantum force may exhibit divergent behavior. Similar phenomena are well recognized in electronic-structure calculations and have been treated by applying the correction of DePasquale et al.,59,60 an approach that also has been employed for the present analyses. The following modification to the algorithm arises from an inexact position of the nodes in trial states, which allows walkers to sample regions that normally should be forbidden in the limit of the exact nodal surface. In particular, this prescription states ⎧ ⎪ E l (x ⃗ ) El(x ⃗) − ET < ⎪ ⎪ E l (x ⃗ ) ← ⎨ αv ⎪ ET + otherwise Δτ ⎪ ⎪ sgn(E (x ⃗) − E ) ⎩ l T

(9)

wμ(τ + Δτ ) ← wμ(τ )

|ψT(xμ′⃗ )|2 G(xμ′⃗ → xμ⃗ , Δτ ) |ψT(xμ⃗ )|2 G(xμ⃗ → xμ′⃗ , Δτ )

(10)

As discussed for example in ref 60, this prescription removes time-step errors in sampling of the correct state in the limit of a perfect guiding state. Operationally, such additional reweighing was limited to the [0.1,2.0] interval and was not used for ground-state simulations, which are less challenging for DMC when applied to the vibrational problem.

3. RESULTS Model calculations have been performed on isotopologs of water, H2O, and the trihydrogen cation, H3+, with the optimized coordinates and the atomic masses used for each species being relegated to the Supporting Information. Primitive coordinates, which form the basis of linear transformations defining their optimized counterparts, were selected to be simple bond lengths and valence angles. This implies rH(D)O and ΘH(D)OH for water and similar quantities for H3+ in which oxygen is replaced by a unique isotope of hydrogen. In what follows, pc(oc)-VCI(n) denotes vibrational configuration interaction calculations up to the nth level of excitation employing primitive (optimized) coordinates (see ref 42 for details). Additionally, pc(oc)-DMC(n) defines fixed-node importance-sampling DMC calculations with pc(oc)-VCI(n) used as the guiding state. For ground-state simulations, VCI(1) is the same as vibrational self-consistent field (VSCF). The exact solution, as computed from vibrational full configuration interaction, is labeled as VFCI. Fundamental transitions are numbered using spectroscopic labels (sorting by symmetry and energy) but appear in energy-ascending ordering as they are used in the orthogonalization procedure. The size of the systems studied allows the DMC results to be compared

αv Δτ

(8)

where ET = ⟨ψT|Ĥ |ψT⟩/⟨ψT|ψT⟩ is the energy of the trial state. The parameter αv is chosen to be 0.025, which allows variation C


Article

Journal of Chemical Theory and Computation directly with the exact solutions for a given potential energy surface. The error arises entirely from the inexact nodal surfaces and the approximated orthogonalization procedure employed in the current work. As the orthogonalization correction is needed to avoid collapse to a single nodal pocket, it is difficult to disentangle the contributions to the error from individual sources. 3.1. H2O Isotopomers. Assessment of the proposed DMC approach for modeling vibrational excited states starts with H2O and HDO as described by the potential of ref 62 but having masses consistent with those utilized in our previous work.42 All calculations in this section employed 2000 walkers (for excited states equally divided into signed subspaces) propagated for 105 time steps and an equilibration time of 3000 a.u., after which data collection started. A set of 10 statistically independent simulations was used for each time step.63 The results for zero-point energies (zpe) are presented in Table 1. Clearly, DMC with any reasonable guiding wave

Figure 1. Dependence of fundamental transition energies (with respect to VFCI) for H2O with a VCI(1) guiding wave function and primitive (left panel) or optimized (right panel) coordinates. Leastsquares regressions based upon a model linear in Δτ are presented along with 99% confidence bands. All calculations used αv = 0.025 and λ = 0.6. VFCI corresponds to VCI(10) in primitive coordinates.

Table 1. Error in Zero-Point Energies (E − EVFCI) of Water Isotopomers as a Function of Guiding Wavefunctiona H2O pc-VCI(1) 13.0 oc-VCI(1) 2.8

pc-DMC(1) 0.2 ± 1.4 oc-DMC(1) 0.1 ± 0.2

pc-VCI(1) 13.5 oc-VCI(1) 2.8

pc-DMC(1) −0.3 ± 1.0 oc-DMC(1) 0.1 ± 0.9

pc-VCI(2) 0.3 oc-VCI(2) 1.9 HDO pc-VCI(2) 0.3 oc-VCI(2) 1.9

pc-DMC(2) 0.0 ± 0.0 oc-DMC(2) −0.1 ± 0.3

VFCI 4636.8 VFCI 4636.8

pc-DMC(2) 0.0 ± 0.1 oc-DMC(2) 0.0 ± 0.2

VFCI 4021.6 VFCI 4021.6

a

Variational and DMC results based on VCI(1) and VCI(2) wavefunctions with primitive (pc) and optimized (oc) coordinates are shown. All energies are in cm−1 with error bars representing 99% confidence limits. VFCI corresponds to VCI(10) in primitive coordinates. DMC values were extrapolated to zero time step analogously to the excited states (see text for details).

function yields an exact simulation, reflecting the fact that the ground state for the problem at hand is nodeless. Nevertheless, a better reference state yields a lower statistical error for a given length of simulation. The time-step dependence of excited-states energies is shown in Figures 1 and 2 for H2O and HDO, respectively, together with a linear regression in Δτ. The results of extrapolating to zero time step are compiled in Table 2. To verify that calculations are converged with respect to simulation length and to prove that estimates are not transient, separate averages were computed for each Δτ, where the equilibration time amounted to 5 × 104 Δτ. In each case, the mean for any excited state was not affected by more than 0.3 cm−1, which was well within the reported confidence intervals. For the case of the ν3 antisymmetric-stretching fundamental in H2O, the importance-sampling DMC with either primitive or optimized coordinates delivers a virtually exact result, representing a substantial improvement over the VCI(1) energy. The result is particularly striking in the case of the oc-VCI(1) energies, where the ν3 transition has the largest error. The description of the bending mode, ν2, also benefits from combining VCI(1) with a fixed-node DMC. For optimized coordinates, the transition-energy error is decreased over 10-fold. The ν1 (symmetric-stretching) excitation energy

Figure 2. Dependence of fundamental transition energies (with respect to VFCI) for HDO with a VCI(1) guiding wave function and primitive (left panel) or optimized (right panel) coordinates. Leastsquares regressions based upon a model linear in Δτ are presented along with 99% confidence bands. All calculations used αv = 0.025 and λ = 0.6. VFCI corresponds to VCI(10) in primitive coordinates.

improves by a factor of 4 when combining oc-VCI(1) with DMC.64 For the HDO isotopolog, a consistent improvement in predicted transition energies again is obained between ocVCI(1) and oc-DMC(1) calculations; however, this is not the case for the pc-VCI(1) and pc-DMC(1) analogs. Given the lower symmetry of the molecule, such behavior essentially reflects the fact that a small VCI with primitive coordinates is bound to provide less accurate nodal surface(s) than their variationally optimized counterparts. Clearly the variationally optimized coordinates are superior not only for variational calculations but also for importancesampling DMC. Of course, the oc-DMC(1) results are not perfect as the corresponding oc-VCI(1) nodal surfaces are not exact. Furthermore, inexact VCI states have been used for the orthogonalization procedure. D


Article

Journal of Chemical Theory and Computation Table 2. Error in Fundamental Excitation Energies (E − EVFCI) of Water Isotopomers for VCI(1) and DMC Guided by VCI(1) Wavefunctions for Primitive (pc) and Optimized Coordinates (oc)a H2O ν2 ν1 ν3

ν3 ν2 ν1

pc-VCI(1)

pc-DMC(1)

16.5 11.2 5.4

2.5 ± 1.6 13.1 ± 1.7 −0.2 ± 1.8

oc-VCI(1)

pc-VCI(1)

pc-DMC(1)

oc-VCI(1)

oc-DMC(1)

VFCI

16.3 12.1 9.3

3.2 ± 1.5 12.6 ± 1.5 8.6 ± 1.4

12.9 8.2 15.2

1.2 ± 1.6 3.6 ± 1.1 1.8 ± 0.9

1403.1 2723.1 3706.6

13.8 4.4 14.5 HDO

oc-DMC(1)

VFCI

1.2 ± 0.4 −1.1 ± 1.0 −1.2 ± 1.9

1594.4 3656.1 3755.0

Figure 4. Dependence of selected overtone and combination transition energies (with respect to VFCI) for H2O (left panel) and HDO (right panel) with a VCI(2) guiding wave function and optimized coordinates. Least-squares regressions based upon a model linear in Δτ are presented along with 99% confidence bands. All calculations used αv = 0.025 and λ = 0.6. VFCI corresponds to VCI(10) in primitive coordinates.

All energies are in cm−1 with error bars representing 99% confidence limits. VFCI corresponds to VCI(10) in primitive coordinates. a

lying states and to address the aforementioned issue without the expense of a full VCI(3) calculation currently are being investigated. Since our formulation of DMC does not require a priori definition of nodes, the application to combination and overtone bands is as straightforward as that of fundamental transitions. In particular, it is sufficient to populate all the nodal pockets. In the present examples, this was accomplished by randomly displacing atoms until enough walkers of a given sign were collected. This also could be done based on the walkers from ground-state diffusion, setting the signs according to the values of the trial wave function, and scaling the weights of the positive and negative regions separately so as to achieve approximate orthogonality to the ground state. The ensuing diffusion with the orthogonalization terms would then take care of the rest. 3.2. H3+ Isotopomers. The previous section showed that our algorithm offers a “black-box” improvement over VCI for a case where the latter method is quite accurate. A more stringent test is provided by the H2D+ and D2H+ ions. Previous work on these systems employing fixed-node diffusion Monte Carlo showed a strong dependence of accuracy on the choice of coordinates used to specify nodal surfaces.65 While different coordinates could be exploited to reproduce the exact (large VCI) results for a given mode, no unique choice of coordinates was satisfactory for all vibrational degrees of freedom. Hence, it is particularly interesting to examine the usefulness of trial wave functions in which the constituent internal coordinates have been determined by variational optimization. Details of VSCF and VCI analyses conducted on these systems can be found in the Supporting Information. Calculations were performed with 5000 walkers (equally divided into signed subspaces for excited states) and an equilibration time of 5000 a.u., after which data were collected for the remainder of 2 × 105 steps. Each simulation consisted of 20 independent analyses.66 The time-step dependencies for the fundamental transitions of H2D+ and D2H+ are depicted in Figures 5 and 6, respectively. The ground-state energy is not shown, but the results of a linear time-step extrapolation for the zpe is presented in Table 4, together with the corresponding (extrapolated) values for excited states.

Figure 3. Dependence of the fundamental transitions energies (with respect to VFCI) for H2O (left panel) and HDO (right panel) with a VCI(2) guiding wave function and optimized coordinates. Leastsquares regressions based upon a model linear in Δτ are presented along with 99% confidence bands. All calculations used αv = 0.025 and λ = 0.6. VFCI corresponds to VCI(10) in primitive coordinates.

Figures 3 and 4, as well as Table 3, present results obtained by using a VCI(2) description for fundamentals and selected overtone/combination bands. Of course, the oc-VCI(2) values for the fundamental excitation energies are better than those of oc-VCI(1), but DMC produces a further systematic improvement, with the largest remaining error being 1.4 cm−1. The ocVCI(2) overtone and combination transitions have significant energy errors for both water isotopologs that also are reduced by factors of 3−10 by DMC. In most cases, the accuracy of ocDMC(2) is comparable to that of oc-VCI(4) without the exponential growth in cost that accompanies increasing levels of VCI. The worst case is the ν1+ν3 feature of HDO, where ocDMC(2) still is comparable in accuracy to oc-VCI(3). It is likely that this discrepancy reflects a qualitative change in the oc-VCI(2) spectrum just below the studied combination band (e.g., a large mixing with the second overtone of the bending mode, not present in the VCI(2) basis, is observed in the exact solution). Even so, DMC reduces the error in this transition energy by a factor of 3. Systematic ways to identify such lowerE


Article


Table 3. Error in Selected Excitation Energies (E − EVFCI) of Water Isotopologs for VCI(2) and DMC Guided by VCI(2) Wavefunctions with Optimized Coordinates (oc)a H2O oc-VCI(2) oc-DMC(2) VFCI

ν2

ν1

ν3

2ν2

ν1+ν2

4.4 0.7 ± 0.6 1594.4

4.8 1.4 ± 1.0 3656.1

3.9 −0.8 ± 0.6 3755.0

40.0 4.0 ± 1.1 3150.8

31.7 8.1 ± 1.8 5233.7

HDO oc-VCI(2) oc-DMC(2) VFCI a

ν3

ν2

ν1

2ν3

ν1+ν3

4.3 0.8 ± 0.5 1403.1

3.5 0.5 ± 0.5 2723.1

6.0 0.8 ± 0.4 3706.6

34.7 4.7 ± 1.9 2781.2

51.0 17.4 ± 3.4 5088.3

All energies are in cm−1, with error bars representing 99% confidence limits. VFCI corresponds to VCI(10) in primitive coordinates.

Table 4. Error in the Ground-State Energy (zpe) and Selected Excited-State Transition Energies (E − Eref) for H3+ Isotopologsa H2D+ oc-VCI(1) ocDMC(1) ref 68

oc-VCI(1) ocDMC(1) ref 69

zpe

ν2

ν3

ν1

49.9 −0.1 ± 1.2

148.1 −2.3 ± 2.3

151.2 −3.0 ± 1.3

61.9 16.7 ± 1.3

3978.3

2205.9 D2H+

2335.3

2992.5

zpe

ν2

ν3

ν1

61.6 0.7 ± 0.9

112.1 −0.5 ± 2.5

145.4 −3.1 ± 0.9

148.4 15.2 ± 2.2

3561.6

1967.8

2079.2

2737.3

All energies are in cm−1 with error bars representing 99% confidence intervals. Reference zpe taken from in-house VCI(11) calculations in primitive coordinates with a large single-mode basis (see Supporting Information). The fundamental transitions are taken from refs 68 and 69, which employed a different potential than the present work and ref 65, but agree to within 1 cm−1 with VCI(11) using our potential. a

Figure 5. Dependence of fundamental transition energies (with respect to VFCI) for H2D+. Least-squares analyses based upon polynomials in Δτ of indicated degree are presented with results obtained for degree 0 excluding the point at Δτ = 2 a.u. The error bars represent 99% confidence intervals. All calculations used αv = 0.025 and λ = 0.6. VFCI corresponds to VCI(11) in primitive coordinates.

systems. While a linear scaling is expected from the local energy cutoff of eq 8,59,60 a much better fit is obtained by using quadratic polynomials in Δτ .67 The energy of the ν3 fundamental has been obtained by neglecting the point at Δτ = 2 a.u. and modeling the remaining data by a constant. While the variationally optimized coordinates for the water ground state involved predominantly localized bond distances, the VSCF ground states for H2D+ and D2H+ entail symmetric and antisymmetric combinations of the analogous parameters. Other sets of optimized coordinates involving localized bonds have been identified but reside higher in energy. For the ground states of the H3+ isotopologs, Table 4 shows VSCF energieseven those involving optimized coordinates to have significantly larger deviations from exact results than those obtained for other molecules studied so far. As expected, this error is removed entirely by DMC. Indeed, the effect of using our wave function for importance sampling of the ground state is a significant reduction in the number of time steps required to reach a given statistical error, but no change in the final outcome. The oc-VCI(1) ansatz does not predict fundamental transitions quantitatively for either H2D+ or D2H+. In fact, rather slow convergence of such energies with VCI size has been observed, signaling the importance of higher-lying states for accurate simulations. Given the ∼10% errors in oc-VCI(1)

Figure 6. Dependence of fundamental transition energies (with respect to VFCI) for D2H+. Least-squares analyses based upon polynomials in Δτ of indicated degree are presented with results obtained for degree 0 excluding the point at Δτ = 2 a.u. The error bars represent 99% confidence intervals. All calculations used αv = 0.025 and λ = 0.6. VFCI corresponds to VCI(11) in primitive coordinates.

In contrast to the water isotopologs, a much stronger Δτ dependence of ν1 and ν2 transition energies is found for the H3+ F


Article


each hyperplane adjusted to make the local energies in each pocket converge to the same value51,74 (a necessary but not sufficient condition for the nodal surface to be exact). This approach sidesteps the requirement that a qualitatively reasonable guiding wave function be available, avoids the difficulties in orthogonalization in the context of an unguided method, and mitigates the exponentially growing complexity of finding arbitrary nodal surfaces. In exchange for eluding these problems, the method makes the approximation that the nodal surfaces are flat hyperplanes (when defined in terms of simple internal coordinates) and requires that a particular internal coordinate be chosen to specify each hyperplane. Of course, when the nodal surface is defined by symmetry, such methods are exact. Should this not be the case, adjustment of the hyperplane location is equivalent to a mean-field description in which the choice of coordinate/hyperplane determines the split between the degree of freedom describing the excited state and the remainder. This clearly is a more accurate mean-field approximation than VSCF because it does not constrain wave functions to be simple products, but it is not variational in nature. Table 5 compares the present results for H3+ isotopologs with the high-quality calculations of ref 65. In the latter work,

excitation energies, use of the corresponding wave functions for importance sampling really tests their most important qualitative characteristics. Indeed, Table 4 shows that importance-sampling, fixed-node DMC exploiting the ocVCI(1) wave functions dramatically enhances the accuracy of transition energies, with the ν1 mode of H2D+ improving by almost a factor of 4 and all other fundamentals improving by factors of 10 or more. This level of performance in the absence of implicit assumptions is comparable to that of ref 65, where different coordinates had to be tested for each isotopolog to produce the best results. Even for this challenging case, where oc-VCI(1) transition energies have errors of ∼100 cm−1, our wave functions are of sufficient quality to produce accurate DMC results for all states without an a priori choice of coordinates and with increased efficiency.

4. COMPARISON WITH PREVIOUS WORK The above results show that without coordinate optimization, VSCF and small VCI wave functions are problematic as guiding wave functions. Not only do the wave functions themselves have large errors, but the improvement from using them with DMC also is erratic. Until the present work, approximate wave functions have been unreliable as trial wave functions, especially for demanding cases such as the isotopologs of H3+. As such, it is natural that previous DMC studies on vibrational excited states have avoided using such trial states (for a notable exception, see ref 70). 4.1. Other Orthogonalization Methods. Our approach is related to other methods of explicitly enforcing orthogonality between DMC states.46 In such schemes, each state is orthogonalized to wave functions for every lower state implied by the distribution of walkers from DMC. This removes the limitations of fixed nodal surfaces(s) and of orthogonalization to approximate wave functions present in our work. However, the canonical representation of wave functions as sums of delta functions (distributions of walkers) contains far more noise than is present in the energy which accumulates as states are orthogonalized55 and leads to unacceptable statistical uncertainties for higher-lying states. Alternatively, the distribution for each state can be fit to an assumed analytic form71 that can be used subsequently in the orthogonalization of higher states.72 This is a significant approximation if a relatively simple analytical form is exploited, while the alternative of deploying a complex ansatz with many terms will tend to fit the noise as well as the underlying wave function, therefore failing to prevent the accumulation of noise in higher states. 4.2. Solution of General Nodal Surfaces. There have been efforts to determine a completely general nodal surface, at least for the case of fundamental modes, based on imposing conditions satisfied by the exact nodal surface.73 Working with an arbitrary (N − 1)-dimensional surface in an N-dimensional space has an intrinsic computational complexity that grows exponentially with the size of the problem. Consequently, whenever this method has been applied to polyatomic systems, a dimensionality reduction scheme has been implemented for evaluation of the nodal constraints. This is effectively a meanfield description in that an a priori separation of one degree of freedom is made with the rest being averaged out. 4.3. Unguided DMC with Assumed Nodal-Surface Shape. Another approach to the vibrational excited-state problem is to use unguided DMC in which the nodal surface(s) are defined by fixed hyperplane(s) in a space of internal coordinates, with the one parameter specifying the location of

Table 5. Errors in Fundamental Transitions for H2D+ and D2H+ Emerging from the Present Work (oc-DMC(1)) and from ref 65 (Qsym-DMC and Qint-DMC)a H2D+ oc-DMC(1) Qsym-DMC Qint-DMC

ν2

ν3

ν1

−2.3 ± 2.3 96.4 ± 10.9 9.3 ± 4.3

−3.0 ± 1.3 −2.7 ± 4.4

16.7 ± 1.3 −442.0 ± 3.9 12.4 ± 4.6

D2H+ oc-DMC(1) Qsym-DMC Qint-DMC

ν2

ν3

ν1

−0.5 ± 2.5 136.1 ± 9.2 22.0 ± 2.8

−3.1 ± 0.9 −1.9 ± 2.9

15.2 ± 2.2 7.3 ± 3.8 26.9 ± 2.4

All energies are in cm−1 with error bars representing 99% confidence limits. a

DMC was employed with two reasonable choices of coordinates to define the nodal surfaces, denoted here as symmetry coordinates (Qsym-DMC) and the normal modes in internal coordinates (Qint-DMC). Very accurate fundamentals were obtained with the best selection of coordinates for both H2D+ and D2H+, but different choices were necessary to obtain the best results in the two seemingly similar cases. DMC using our oc-VCI(1) wave functions, which does not require an a priori choice of coordinates, works equally well for both systems. Furthermore, the use of internal-coordinate hyperplanes as nodal surfaces becomes more problematic for overtone and combination bands, unless at least one of the nodal surfaces is constrained by symmetry. In contrast, use of oc-VCI(2) trial wave functions makes possible the application of DMC methods to any two-quantum overtone or combination feature.

5. CONCLUSIONS AND OUTLOOK A systematic approach to use compact variational wave functions for importance sampling in DMC has been presented for both ground and excited states. This makes possible the G


Article


(7) Bloino, J.; Biczysko, M.; Barone, V. General Perturbative Approach for Spectroscopy, Thermodynamics, and Kinetics: Methodological Background and Benchmark Studies. J. Chem. Theory Comput. 2012, 8, 1015−1036. (8) Jung, J. O.; Gerber, R. B. Vibrational wave functions and spectroscopy of (H2O)n, n = 2,3,4,5: Vibrational self consistent field with correlation corrections. J. Chem. Phys. 1996, 105, 10332−10348. (9) Krasnoshchekov, S. V.; Isayeva, E. V.; Stepanov, N. F. NumericalAnalytic Implementation of the Higher-Order Canonical van Vleck Perturbation Theory for the Interpretation of Medium-Sized Molecule Vibrational Spectra. J. Phys. Chem. A 2012, 116, 3691−3709. (10) Chaban, G. M.; Jung, J. O.; Gerber, R. B. Ab initio calculation of anharmonic vibrational states of polyatomic systems: Electronic structure combined with vibrational self-consistent field. J. Chem. Phys. 1999, 111, 1823−1829. (11) Christiansen, O. A second quantization formulation of multimode dynamics. J. Chem. Phys. 2004, 120, 2140−2148. (12) Christiansen, O. Vibrational coupled cluster theory. J. Chem. Phys. 2004, 120, 2149−2159. (13) Banik, S.; Pal, S.; Prasad, M. D. Calculation of vibrational energy of molecule using coupled cluster linear response theory in bosonic representation: Convergence studies. J. Chem. Phys. 2008, 129, 134111. (14) Panek, P. T.; Jacob, C. R. Efficient Calculation of Anharmonic Vibrational Spectra of Large Molecules with Localized Modes. ChemPhysChem 2014, 15, 3365−3377. (15) Avila, G.; Carrington, T., Jr. Using a pruned basis, a non-product quadrature grid, and the exact Watson normal-coordinate kinetic energy operator to solve the vibrational Schrödinger equation for C2H4. J. Chem. Phys. 2011, 135, 064101. (16) Vendrell, O.; Meyer, H.-D. Multilayer multiconfiguration timedependent Hartree method: Implementation and applications to a Henon-Heiles Hamiltonian and to pyrazine. J. Chem. Phys. 2011, 134, 044135. (17) Thomas, P. S.; Carrington, T. Using Nested Contractions and a Hierarchical Tensor Format To Compute Vibrational Spectra of Molecules with Seven Atoms. J. Phys. Chem. A 2015, 119, 13074− 13091. (18) Rakhuba, M.; Oseledets, I. Calculating vibrational spectra of molecules using tensor train decomposition. J. Chem. Phys. 2016, 145, 124101. (19) Wodraszka, R.; Carrington, T. Using a pruned, nondirect product basis in conjunction with the multi-configuration timedependent Hartree (MCTDH) method. J. Chem. Phys. 2016, 145, 044110. (20) Brown, J.; Carrington, T., Jr. Using an expanding nondirect product harmonic basis with an iterative eigensolver to compute vibrational energy levels with as many as seven atoms. J. Chem. Phys. 2016, 145, 144104. (21) Leclerc, A.; Carrington, T., Jr. Using symmetry-adapted optimized sum-of-products basis functions to calculate vibrational spectra. Chem. Phys. Lett. 2016, 644, 183−188. (22) Baiardi, A.; Stein, C. J.; Barone, V.; Reiher, M. Vibrational Density Matrix Renormalization Group. J. Chem. Theory Comput. 2017, 13, 3764−3777. (23) Sibert, E. L., III Theoretical studies of vibrationally excited polyatomic molecules using canonical Van Vleck perturbation theory. J. Chem. Phys. 1988, 88, 4378−4390. (24) McCoy, A. B.; Burleigh, D. C.; Sibert, E. L., III Rotationvibration interactions in highly excited states of SO2 and H2CO. J. Chem. Phys. 1991, 95, 7449−7465. (25) Bramley, M. J.; Carrington, T., Jr. A general discrete variable method to calculate vibrational energy levels of three- and four-atom molecules. J. Chem. Phys. 1993, 99, 8519−8541. (26) Wang, X.-G.; Carrington, T., Jr. Vibrational energy levels of CH5+. J. Chem. Phys. 2008, 129, 234102. (27) Handy, N. C.; Carter, S.; Colwell, S. M. The vibrational energy levels of ammonia. Mol. Phys. 1999, 96, 477−491.

application of DMC to fundamental modes having neither symmetry constraints nor an a priori choice of coordinates, and to overtone and combination bands even when nodal surfaces are unconstrained by symmetry. The combination of oc-VCI and DMC is a “black-box” method which has both high accuracy and low scaling with system size. The scaling of the computational cost of DMC is lower than that for the solution of the trial wave functions. The latter involve numerical quadrature of matrix elements as well as optimization of coordinates for each state, so while cost scales as a polynomial, the order of the polynomial is higher for the trial state optimizations. The application of these methods to larger systems, as well as broader tests of their accuracy, will be the subject of future papers.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b00016. Details of VSCF calculations, atomic masses, and coordinates used and further details regarding presented algorithm (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Ireneusz W. Bulik: 0000-0002-0401-4376 Patrick H. Vaccaro: 0000-0001-7178-7638 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Dr. Anne McCoy for numerous helpful discussions about the details of vibrational DMC and Victor Lee for his comments on the manuscript. I.W.B. also thanks Dr. Giovanni Scalmani and Dr. Fernando Clemente for helpful discussions. Financial support provided by Gaussian Inc. and computational time provided by Yale Center for Research Computing are gratefully acknowledged. P.H.V. gratefully acknowledges support of the National Science Foundation under grant CHE-1464957.

■

REFERENCES

(1) 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. (2) Nielsen, H. H. The Vibration-Rotation Energies of Molecules. Rev. Mod. Phys. 1951, 23, 90−136. (3) Willetts, A.; Handy, N. C.; Green, W. H.; Jayatilaka, D. Anharmonic corrections to vibrational transition intensities. J. Phys. Chem. 1990, 94, 5608−5616. (4) Barone, V. Anharmonic vibrational properties by a fully automated second-order perturbative approach. J. Chem. Phys. 2005, 122, 014108. (5) Bloino, J.; Barone, V. A second-order perturbation theory route to vibrational averages and transition properties of molecules: General formulation and application to infrared and vibrational circular dichroism spectroscopies. J. Chem. Phys. 2012, 136, 124108. (6) Bloino, J. A VPT2 Route to Near-Infrared Spectroscopy: The Role of Mechanical and Electrical Anharmonicity. J. Phys. Chem. A 2015, 119, 5269−5287. H


Article


(50) Anderson, J. B. Quantum chemistry by random walk: Higher accuracy. J. Chem. Phys. 1980, 73, 3897−3899. (51) McCoy, A. B. Diffusion Monte Carlo approaches for investigating the structure and vibrational spectra of fluxional systems. Int. Rev. Phys. Chem. 2006, 25, 77−107. (52) In the case of signed walkers, the population control is performed within a given sign space. (53) Ceperley, D. M. Fermion nodes. J. Stat. Phys. 1991, 63, 1237− 1267. (54) Foulkes, W. M. C.; Hood, R. Q.; Needs, R. J. Symmetry constraints and variational principles in diffusion quantum Monte Carlo calculations of excited-state energies. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 60, 4558−4570. (55) Hammond, B. L.; Lester, W. A., Jr.; Reynolds, P. J. Monte Carlo Methods in Ab Initio Quantum Chemistry; World Scientific, 1994; Chapter 6, pp 181−197. (56) Williamson, A. J.; Hood, R. Q.; Needs, R. J.; Rajagopal, G. Diffusion quantum Monte Carlo calculations of the excited states of silicon. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 57, 12140− 12144. (57) While other ansätze are possible, the one presented is simple and similar to the Anderson stabilizing term, which works well. (58) For the vibrational problem, both the wave functions and the potentials are smooth; hence the VSCF and VCI can be used directly as guiding states without numerical issues and additional steps such as optimizing Jastrow factors etc. (59) DePasquale, M. F.; Rothstein, S. M.; Vrbik, J. Reliable diffusion quantum Monte Carlo. J. Chem. Phys. 1988, 89, 3629−3637. (60) Umrigar, C. J.; Nightingale, M. P.; Runge, K. J. A diffusion Monte Carlo algorithm with very small timestep errors. J. Chem. Phys. 1993, 99, 2865−2890. (61) Marlett, M. L.; Lin, Z.; McCoy, A. B. Rotation/Torsion Coupling in H5+, D5+, H4D+, and HD4+ Using Diffusion Monte Carlo. J. Phys. Chem. A 2015, 119, 9405−9413. (62) Partridge, H.; Schwenke, D. W. The determination of an accurate isotope dependent potential energy surface for water from extensive ab initio calculations and experimental data. J. Chem. Phys. 1997, 106, 4618−4639. (63) Each simulation consisted of a separately seeded calculation with its own equilibration period. No information was shared between different runs. The average energies from each simulation were then used to produce an average energy and error estimates. (64) The situation is slightly different with primitive coordinates and the symmetric stretch, ν1, because the pc-VCI(1) excitation energy happens to benefit from a cancellation of errors for ground- and excited-state energies. Since the ground-state energy is essentially exact, this cancellation is removed in the pc-DMC(1) excitation energy. (65) Petit, A. S.; McCoy, A. B. Diffusion Monte Carlo in Internal Coordinates. J. Phys. Chem. A 2013, 117, 7009−7018. (66) For the ground state, data were collected in 10 independent runs having 105 steps each, and for the excited states, with Δτ parameters of 1.0 a.u. and 2.0 a.u., 10 independent runs were performed. These cases converged much more rapidly than the excited states using smaller time steps. (67) The Δτ time-step dependence was used in ref 43, justified by the higher-order correction to the propagator, as well as in ref 60, where such scaling was argued to result from the killing of walkers that cross the node. In the latter, the energy affected by such behavior was computed in a different way than in the present work. Nonetheless, the current study of time-step dependence suggests that fractional powers in Δτ are needed for our incarnation of the DMC algorithm. The Δτ dependence is observed for transitions in the H3+ systems where the trial state of VCI(1) quality is not expected to be very accurate and the nodal surface is not determined by symmetry alone. (68) Sochi, T.; Tennyson, J. A computed line list for the H2D+ molecular ion. Mon. Not. R. Astron. Soc. 2010, 405, 2345−2350.

(28) Carter, S.; Pinnavaia, N.; Handy, N. C. The vibrations of formaldehyde. Chem. Phys. Lett. 1995, 240, 400−408. (29) Vendrell, O.; Gatti, F.; Lauvergnat, D.; Meyer, H.-D. Fulldimensional (15-dimensional) quantum-dynamical simulation of the protonated water dimer. I. Hamiltonian setup and analysis of the ground vibrational state. J. Chem. Phys. 2007, 127, 184302. (30) Vendrell, O.; Gatti, F.; Meyer, H.-D. Full dimensional (15dimensional) quantum-dynamical simulation of the protonated water dimer. II. Infrared spectrum and vibrational dynamics. J. Chem. Phys. 2007, 127, 184303. (31) Lauvergnat, D.; Nauts, A. Quantum dynamics with sparse grids: A combination of Smolyak scheme and cubature. Application to methanol in full dimensionality. Spectrochim. Acta, Part A 2014, 119, 18−25. (32) Avila, G.; Carrington, T., Jr. Nonproduct quadrature grids for solving the vibrational Schrödinger equation. J. Chem. Phys. 2009, 131, 174103. (33) Yurchenko, S. N.; Thiel, W.; Jensen, P. Theoretical ROVibrational Energies (TROVE): A robust numerical approach to the calculation of rovibrational energies for polyatomic molecules. J. Mol. Spectrosc. 2007, 245, 126−140. (34) Mátyus, E.; Czakó, G.; Császár, A. G. Toward black-box-type full- and reduced-dimensional variational (ro)vibrational computations. J. Chem. Phys. 2009, 130, 134112. (35) Fábri, C.; Sarka, J.; Császár, A. G. Communication: Rigidity of the molecular ion H5+. J. Chem. Phys. 2014, 140, 051101. (36) Owens, A.; Yurchenko, S. N.; Yachmenev, A.; Tennyson, J.; Thiel, W. Accurate ab initio vibrational energies of methyl chloride. J. Chem. Phys. 2015, 142, 244306. (37) Sarka, J.; Császár, A. G. Interpretation of the vibrational energy level structure of the astructural molecular ion H5+ and all of its deuterated isotopomers. J. Chem. Phys. 2016, 144, 154309. (38) Bowman, J. M. Selfconsistent field energies and wavefunctions for coupled oscillators. J. Chem. Phys. 1978, 68, 608−610. (39) Wright, N. J.; Gerber, R. B. Direct calculation of anharmonic vibrational states of polyatomic molecules using potential energy surfaces calculated from density functional theory. J. Chem. Phys. 2000, 112, 2598−2604. (40) 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. (41) Christiansen, O. Selected new developments in vibrational structure theory: potential construction and vibrational wave function calculations. Phys. Chem. Chem. Phys. 2012, 14, 6672−6687. (42) Bulik, I. W.; Frisch, M. J.; Vaccaro, P. H. Vibrational selfconsistent field theory using optimized curvilinear coordinates. J. Chem. Phys. 2017, 147, 044110. (43) Reynolds, P. J.; Ceperley, D. M.; Alder, B. J.; Lester, W. A., Jr. Fixednode quantum Monte Carlo for molecules. J. Chem. Phys. 1982, 77, 5593−5603. (44) Mitas, L. In Quantum Monte Carlo Methods in Physics and Chemistry; Nightingale, M., Umrigar, C., Eds.; Nato Science Series C, Springer: Dordrecht, The Netherlands, 1998; Chapter 9. (45) Pederiva, F.; Roggero, A.; Schmidt, K. E. In An Advanced Course in Computational Nuclear Physics: Bridging the Scales from Quarks to Neutron Stars; Hjorth-Jensen, M., Lombardo, M. P., van Kolck, U., Eds.; Springer International Publishing: Cham, Switzerland, 2017; pp 401−476. (46) Suhm, M. A.; Watts, R. O. Quantum Monte Carlo studies of vibrational states in molecules and clusters. Phys. Rep. 1991, 204, 293− 329. (47) The continuous reweighting approach was chosen due to ease of implementation, especially for the orthogonalization correction implemented in eq 7, as discussed in the following text. (48) Anderson, J. B. A randomwalk simulation of the Schrödinger equation: H3+. J. Chem. Phys. 1975, 63, 1499−1503. (49) 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. I


Article

Journal of Chemical Theory and Computation (69) Miller, S.; Tennyson, J. First principles calculation of the molecular constants of H3+, H2D+, D2H+, and D3+. J. Mol. Spectrosc. 1987, 126, 183−192. (70) Broude, S.; Jung, J. O.; Gerber, R. Combined diffusion quantum Monte Carlo-vibrational self-consistent field (DQMCVSCF) method for excited vibrational states of large polyatomic systems. Chem. Phys. Lett. 1999, 299, 437−442. (71) Coker, D. F.; Watts, R. O. Structure and vibrational spectroscopy of the water dimer using quantum simulation. J. Phys. Chem. 1987, 91, 2513−2518. (72) Sun, H.; Watts, R. O. Diffusion Monte Carlo simulations of hydrogen fluoride dimers. J. Chem. Phys. 1990, 92, 603−616. (73) Sandler, P.; Buch, V.; Sadlej, J. Ground and excited states of the complex of CO with water: A diffusion Monte Carlo study. J. Chem. Phys. 1996, 105, 10387−10397. (74) Cheng, T. C.; Bandyopadyay, B.; Wang, Y.; Carter, S.; Braams, B. J.; Bowman, J. M.; Duncan, M. A. Shared-Proton Mode Lights up the Infrared Spectrum of Fluxional Cations H5+ and D5+. J. Phys. Chem. Lett. 2010, 1, 758−762.

J


Fixed-Node, Importance-Sampling Diffusion Monte Carlo for

Recommend Documents