Langevin Equation Path Integral Ground State - ACS Publications

Article pubs.acs.org/JPCA

Langevin Equation Path Integral Ground State Steve Constable, Matthew Schmidt, Christopher Ing, Tao Zeng, and Pierre-Nicholas Roy* Department of Chemistry, University of Waterloo, Waterloo, Ontario, Canada, N2L 3G1 ABSTRACT: We propose a Langevin equation path integral ground state (LePIGS) approach for the calculation of ground state (zero temperature) properties of molecular systems. The approach is based on a modification of the finite temperature path integral Langevin equation (PILE) method (J. Chem. Phys. 2010, 133, 124104) to the case of open Feynman paths. Such open paths are necessary for a ground state formulation. We illustrate the applicability of the method using model systems and the weakly bound water−parahydrogen dimer. We show that the method can lead to converged zero point energies and structural properties.

■

INTRODUCTION The calculation of ground state properties of many-body quantum mechanical systems is a formidable challenge. For the specific case of nuclear degrees of freedom, within the Born− Oppenheimer approximation, exact basis set methods have a computational cost that scales exponentially with the number of degrees of freedom. Some recent advances in the so-called multi configuration time dependent Hartree method have allowed to push further the size limit for certain classes of problems.1 Quantum Monte Carlo (QMC) methods are commonly used to tackle the many-body problem. The diffusion Monte Carlo (DMC)2 in which the Schrödinger equation is modeled as a diffusion process in imaginary time has been successfully applied to a variety of systems including important water containing ions3,4 and weakly bound hydrogen clusters.5 In a recent study, Boninsegni and Moroni discuss potential issues associated with the population size bias associated with DMC methods.6 One ground state method that does not suffer from the population bias problem is the path integral ground state (PIGS) approach.7,8 The approach has been successfully applied to the study of weakly bound parahydrogen clusters.9 In PIGS, the Green’s function is viewed as a density matrix that is sampled via path integral Monte Carlo (PIMC) techniques.10 This means that specialized Monte Carlo moves have to be designed for the specific types of molecular systems under study. For the case of finite temperature, the PIMC problem can be recast from a stochastic problem to a dynamical problem through the application of path integral molecular dynamics (PIMD).11,12 Use of PIMD has the advantage over PIMC that moves do not have to be designed for the system but are derived from the gradient of the potential. Furthermore, PIMD leads naturally to ring polymer molecular dynamics (RPMD),13 which can be used to obtain approximate time correlation functions. We show here that these concepts from PIMD can © 2013 American Chemical Society

be applied to PIGS by treating the path integral as the configuration component of a classical partition function describing an open polymer. It is known from previous PIMD studies that attempts to increase the accuracy of the results by increasing the Trotter number result in ergodicity problems.14 One approach to ameliorate this issue is the use of advanced thermostatting techniques, which rephrase the PIMD problem in terms of normal modes which do not suffer from the same ergodicity issues. The current paper intends to apply such a technique, namely, the path integral Langevin equation (PILE),15 to the MD treatment of PIGS. This will lead to a novel method to solve the ground state problem where the Langevin equation is used to sample open paths rather than closed paths as in the case of finite temperature studies. The approach allows one to tackle problems that are qualitatively different, since its extends the applicability of the PILE to zero temperature properties. Note that the PILE approach is based on the Fourier representation of the path integral. Fourier path integral (FPI) methods have been extensively used in the context of PIMC sampling.10,16−20 Mak has used blocking techniques to further enhance the efficiency of FPI methods.21 To summarize, the important differences between our proposed approach and the above methods are as follows: (i) Open paths are used in order to describe the ground state; (ii) PIMD sampling is used rather than PIMC. In section 2, a brief derivation of PIGS is shown, followed by incorporation of the PILE. The topology of the MD PIGS fictitious polymer requires that the system be expressed in a discrete cosine transform (DCT) representation rather than the Special Issue: Joel M. Bowman Festschrift Received: February 11, 2013 Revised: June 3, 2013 Published: June 5, 2013 7461

dx.doi.org/10.1021/jp4015178 | J. Phys. Chem. A 2013, 117, 7461−7467

The Journal of Physical Chemistry A

Article

typical Fourier modes.22,23 This hybrid PILE−PIGS approach we name LePIGS, the Langevin equation path integral ground state method. In LePIGS, the normal modes can be propagated analytically in terms of sines and cosines, avoiding the inherent ergodicity problems of the system. Additionally, efficient thermalization for the MD simulation is realized because many of the normal modes can be assigned a Langevin friction analytically. This new LePIGS method is applied to a variety of systems in section 3. The studied systems include the simple harmonic oscillator, for which the ground state wave function with systematic error can be determined analytically, and the quartic oscillator and double-well potential, for which the ground state wave function can be derived numerically through matrix multiplication. This provides a reference for both the final result of the simulation and the convergence in parameters β and τ described below. Furthermore, results of the simulation of H2−H2O are provided, demonstrating the generality of the method. Finally, the conclusions and some discussion are included in section 4.

⟨R |exp( −βĤ )|R′⟩ = ⟨R |exp( −τĤ )|R1⟩⟨R1|exp( −τĤ )|R 2⟩... ⟨RP − 1|exp( −τĤ )|R′⟩ P−1

=

with τ = β/(P − 1) and where a high temperature density operator ρ̂ = exp(−τĤ ) is introduced. This leads to the following expression for Z0: P−1

Z0 =

N i

pi 2̂ 2mi

Z0 = lim

P →∞

⎛ ⎜

P−1

V ′(R ) =

⎧1 ⎪ if j = 1 or j = P cj = ⎨ 2 ⎪ ⎩1 otherwise

(9)

We have chosen to use the case of N = 1 particle in one spacial dimension for simplicity in the above. Generalization to many particles in higher dimensions is straightforward. The expression in eq 7 is isomorphic with that for the configurational partition function for a classical system consisting of an open polymer of P beads in some potential field. For a position dependent property A, the ground state expectation value is given by

⎟

∫ dR ∫ dR′ψT(R)⟨R|e−βĤ |R′⟩ψT(R′)

(8)

with ωP−1 = (P − 1) /βℏ and

(2)

⟨A⟩0 = lim

β , P →∞

1 Z0

∫ dR1... ∫ dRPA(R(P+1)/2)

exp{−τV ′(R )}

(10)

which is the primitive estimator of Â . Note that we have explicitly accounted for the fact that a unity trial function is considered here. Trial functions that differ from unity can be represented in exponential form and will lead to additional potential terms affecting only beads j = 1 and j = P. The above expression shows that expectation vales of operators that are diagonal in the position representation will be associated with the middle bead. If Â commutes with the density operator, exp(−βĤ ), then it is possible to shift Â all the way toward the right-hand side of eq 3 so that it is acting on |ΨT⟩. In the case of the Hamiltonian, Ĥ , we obtain the mixed estimator for the energy:

(3)

where Z0 is a pseudo partition function given by Z0 =

j=1 1/2

⎞

⟨Φ0|Â |Φ0⟩ 1 = ⟨Φ0|Â |Φ0⟩ ⟨Φ0|Φ0⟩ Z0

P

∑ [mωP − 12(R j − R j + 1)2 ] + ∑ [cjV (R j)] j=1

where ΨT(R) is some trial wave function and the parameter β/2 governs the relaxation of the trial wave function to the ground state. The choice of ΨT(R) will of course affect the convergence properties of the integral in eq 2, but in this study, it was always assumed that ΨT(R) = 1. This choice ensures that there is always some overlap between the ground state and the trial wave function, and therefore the integral should lead to the ground state. Consider the evaluation of the expectation value of some property Â in the ground state. The expression for this value is ⟨A⟩0 =

∫ dR1... ∫ dRPψT(R1)exp{−τV ′(R)}ψT(RP)

where the modified potential appearing in the expression for Z0 is defined as

(1)

∫ dR′⟨R|exp⎝− 2β Ĥ ⎠|R′⟩ΨT(R′)

mP 2πβ ℏ2

(7)

+ V̂

β→∞

(6)

The density matrix elements can be obtained using the socalled Trotter factorization24 to yield

where N is the number of atoms, p̂i and mi and the momentum and mass of atom i, and V̂ is some interaction potential. The ground state wave function Φ0(R) = ⟨R|Φ0⟩ in a position basis R is proportional to the following integral Φ0(R ) ∝ lim

∫ dR1... ∫ dRPψT(R1) ∏ [⟨R j|ρ ̂|R j+1⟩]ΨT(RP) j=1

THEORY AND METHODOLOGY Here we demonstrate how the ground state can be found from a trial function, and review the PIGS method via the discrete path integral formulation. A normal mode representation is adopted in imaginary time, leading to incorporation of the PILE. The advantages of this approach are described in detail. Path Integral Ground State. The PIGS method provides an estimate of the nuclear ground state properties of a chemical system described by some Hamiltonian

∑

(5)

j=1

■

Ĥ =

∏ ⟨R j|ρ|R j + 1⟩

E0 =

(4)

In order to compute the quantity ⟨R|exp(−βĤ )|R′⟩, known as the imaginary time propagator, a factorization scheme can be adopted:

Ĥ ΨT(RP) ΨT(RP)

0

(11)

where E0 is the ground state energy and R(P) is the position of the rightmost bead. Equivalently, the leftmost bead may be used by shifting Ĥ all the way to the left and acting on ⟨ΨT|, 7462



Article

since Ĥ is Hermitian. Throughout this paper, the trial wave function is assumed to be one, so the expression for the mixed estimator of the energy may be reduced to

E0 = ⟨V (RP)⟩0

with ωk = 2ωP sin(kπ/2P). With this simplified Hamiltonian, the kinetic spring terms have been replaced with a sum of independent harmonic oscillators, whose equations of motion are known analytically in terms of sines and cosines. This implies a two-step integration scheme, where individual beads are propagated via MD in a scaled classical potential and then propagated analytically in normal mode space to capture the effects of the kinetic spring terms. In order to efficiently thermostat the integration of HP0 (p̃, q̃), the PILE is employed. The PILE dictates that, before and after the propagation of the kinetic springs in normal mode space, the momenta are scaled thusly:

(12)

The Langevin Equation Path Integral Ground State Method. It is our objective to combine the PIGS formulation within a PIMD framework with the efficient sampling features of the PILE thermostat.15,25 This would form a hybrid PILE/ PIGS method, which we term the Langevin equation path integral ground state (LePIGS). However, the PILE is derived from the PIMD Hamiltonian which describes closed paths, while the Hamiltonian for PIGS is isomorphic with open polymers (paths). It is not immediately obvious how to proceed. Equation 3 can be evaluated in a canonical ensemble using a modified N particle Hamiltonian of the form HP(p , q) = N

VP(q) =

HP0(p , 1 P

q) + VP(q)

P

N P−1

HP0(p , q) =

(14)

j=1

⎛ [p(j) ]2

∑ ∑ ⎜⎜ i=1 j=1

i

⎝ 2mi

+

■

RESULTS The LePIGS method was implemented in the freely available Molecular Modeling Tool Kit (MMTK),26 based on the implementation of PIMD-PILE by Ing et al.27 The DCT code was provided by the FFTW28 package. Model Systems. Arguably, the simplest system with which to benchmark the performance of the LePIGS method is that of an electron trapped in a 1D harmonic potential. The Hamiltonian for this system is HHO =

∑ Ckjqi(j)

where Ck is a vector of coefficients for normal mode k whose elements are defined as 1 cos(πk(j − 1/2)/P) k = 0 P

Under these transformations, N P−1

q )̃ =

∑∑ i=1 k=0

H0P

[pi(̃ k) ]2 2mi

(18)

takes on the desired form:

+

1 miωk 2[qĩ (k)]2 2

p2 mω 2 2 + x 2m 2

(24)

and when expressed in atomic units, m = ω = 1. A time step dt = 0.001 fs was for the integration of the equations of motion along with a friction coefficient of γ0 = 104 ps−1. The total simulation time was 4 ps, and data was collected every 0.56 fs to avoid correlations. The ground state density for this system, ρP, is known analytically including corrections for the systematic error of using finite P.7 Both the analytical ground state density and the histogram of the position of the middle bead calculated by LePIGS (with β = 3 and P = 15) are shown in Figure 1. Ground state energy, E0, convergence diagrams in β and τ are provided in Figures 2 and 3, respectively. Figure 2 shows that the energy convergence with respect to β is exponential as expected. For a fixed β value, Figure 3 confirms the expected parabolic convergence associated with the primitive Trotter error.

(17)

2 cos(πk(j − 1/2)/P) k > 0 P

(23)

The final piece of information needed to determine the frictions is the value τ0, which is a free parameter to the simulation and is used to optimize convergence time.

(16)

j=1

HP0(p ̃ ,

(22)

1 − [c1(k)]2

⎪

P−1

⎧ ⎪ ⎪ Ckj = ⎨ ⎪ ⎪ ⎩

c 2(k) =

⎪

P−1

qĩ (k) =

(21)

⎧1/τ0 k = 0 γ (k) = ⎨ ⎩ 2ωk k > 0

where (p, q) are conjugate pairs of momentum and position vectors. There are now P − 1 kinetic spring terms, and no circularity condition is enforced. Additionally, there would be included in VP(q) two terms corresponding to the interactions between the trial wave function, |ΨT⟩, and the first and last beads, but in our case the trial wave function is always unity, so no such term is included in eq 14. It can be noted that the simplicity of the PILE thermostat stems from the form of the Hamiltonian when expressed in terms of normal mode coordinates. An open path can be expressed in terms of modified normal mode coordinates through the use of the discrete cosine transform type II (DCTII) matrix rather than the typical Fourier matrix. Now, we define new normal mode coordinates for PIGS as

j=1

c1(k) = exp( −(Δt /2)γ (k))

where γ is the Langevin friction on the kth normal mode and ξ(k) i is an independent Gaussian number with a mean of 0 and a standard deviation of 1. Statistically optimal values for the Langevin frictions can be analytically determined from the autocorrelation function of HP0 (p̃, q̃) in the case of a free polymer, and are given by the following equation:

⎞ 1 miωP 2[qi(j + 1) − qi(j)]2 ⎟⎟ 2 ⎠

∑ Ckjpi(j)

(20)

(k)

(15)

pi(̃ k) =

mi (k) (k) c 2 ξi β

with

(13)

∑ ∑ V (qi(j)) i=1

pi(̃ k) ← c1(k)pi(̃ k) +

(19) 7463



Article

̀ Vdw(x) = 0.1x 4 − 0.5x 2

(26)

The ground states for these systems at finite P is not known analytically but can be generated numerically by explicitly forming the density matrix ⟨xi|e−βĤ |xi+1⟩ and performing matrix multiplication P times onto an initial ψT vector. This is equivalent to using the numerical matrix multiplication methods.29 In this manner, a reference for the expected convergence in β and τ may be prepared. The LePIGS convergence in β and τ is shown for the quartic oscillator in Figures 4 and 5 and for the double-well in Figures 6 and 7,

Figure 1. Harmonic oscillator ground state position density obtained from LePIGS simulation with β = 3 and P = 15 (filled circles). The exact result (solid line) is also presented for reference.

Figure 4. Exponential convergence of the LePIGS E0 with respect to β for the quartic oscillator (τ = 0.05 hartree−1). The exact value obtained from explicit diagonalization is shown as a dashed line.

Figure 2. Exponential convergence of the LePIGS E0 with respect to β for the harmonic oscillator (τ = 0.05 hartree−1). The exact value of 0.5 hartree is shown as a dashed line.

Figure 5. Convergence of the LePIGS E0 of the quartic oscillator with respect to τ for β = 3. A parabolic fit to the data is shown as a solid line. The exact value obtained from explicit diagonalization is shown as a dashed line.

respectively. The energy converges to the correct value as expected, in both cases. This implies that the code is functioning as intended, and that the methodology is sound. A Weakly Bound Complex: the Water−Hydrogen Dimer. The first application to a “real” system using this technique is the water−hydrogen dimer. The water molecule is treated as a fixed molecule and provides an external potential to the hydrogen molecule, treated as a point-like particle. The water−hydrogen potential we use is from ref 30. This potential has recently been used by our group to study a parawater molecule embedded in superfluid hydrogen clusters.31 The Hamiltonian for a system consisting of a parahydrogen particle of mass m is

Figure 3. Convergence of the LePIGS E0 of the harmonic oscillator with respect to τ for β = 3. A parabolic fit to the data is shown as a solid line. The exact value of 0.5 hartree is shown as a dashed line.

Two additional models are considered next, the quartic oscillator with potential Vq and the double-well potential with potential Vdw. The potentials have the following form (in atomic units): 1 Vq(x) = x 4 (25) 4 7464



Article

chosen on the basis of experimentation and previous low temperature path integral molecular dynamics simulations.27 Several steps are needed to determine the centroid friction value. First, we perform an energy minimization to get our initial configuration, followed by a short LePIGS simulation without friction, but skipping enough steps to ensure uncorrelated configurations. Using each of those uncorrelated configurations as initial configurations, a microcanonical ensemble simulation is run (holding constant number of atoms, volume, energy). The friction value chosen is the averaged time correlation of structural properties over all the microcanonical ensemble simulations. Once the centroid friction is determined, we run another short LePIGS simulation using that friction and not skipping any steps. The correlation times are then determined for each of the structural properties. The longest correlation time is the one used for the simulation. The “exact” energy, radial distribution, and angular distributions were obtained using a basis set calculation for a parahydrogen interacting with a fixed water molecule using the approach of ref 30. As in the previous section, the energy convergence diagram of β is shown in Figure 9. However, the

Figure 6. Exponential convergence of the LePIGS E0 with respect to β for the double well (τ = 0.05 hartree−1). The exact value obtained from explicit diagonalization is shown as a dashed line.

Figure 7. Convergence of the LePIGS E0 of the double well with respect to τ for β = 3. A parabolic fit to the data is shown as a solid line. The exact value obtained from explicit diagonalization is shown as a dashed line.

HH 2 =

p2 + VH2 − H2O(R) 2m

Figure 9. Exponential convergence of the LePIGS E0 with respect to β for the water−parahydrogen dimer (τ = 0.007 K−1).

(27)

where R is the position of the point-like hydrogen molecule in the water molecule frame. The coordinate system is depicted in Figure 8 where the vector R is decomposed into its magnitude, R, and polar and azimuthal angles, θ and χ, respectively. In order to perform the simulations, three parameters need to be determined, the time step dt, the centroid friction γ(0), and the correlation time of the system which is needed to get time-independent data points. A time step dt = 2.0 fs was

energy will not converge to the “exact” value, since our τ value is not ideal. This is not a concern, since the energy will just converge to a lower energy value than the exact ground state energy. A β value of 1.00 K−1 is chosen, since it is safely in the region of the converged energy value. Using this β value, the energy convergence with τ is shown in Figure 10. The black points are the LePIGS simulation energies, and the red dashed line is the “exact” energy. The blue curve is the quadratic fit with the blue point as the energy extrapolated to the τ = 0 limit with its associated error. The LePIGS simulation obtains the exact energy within error at τ = 0.001 and in the limit as τ = 0. The convergence of the radial and angular distributions with β and τ are also shown below and compared to the “exact” distributions. Figures 11, 12, and 13 show the structural distribution dependence on τ. When varying β, there was only a significant difference in the θ distribution shown in Figure 14. In all cases, convergence is systematically achieved when β (τ) is sufficiently large (small). Our results show that the LePIGS method is practical and accurate.

■

CONCLUDING REMARKS We have developed a new approach, LePIGS, that can sample the ground state properties of an arbitrary many-particle

Figure 8. Coordinate system definition for the water−parahydrogen dimer. 7465



Article

Figure 10. Convergence of the LePIGS E0 of the water−parahydrogen dimer with respect to τ for β = 1.00 K−1. A parabolic fit to the data is shown as a solid line. The exact value obtained from explicit diagonalization is shown as a dashed line.

Figure 13. Angular distribution of χ for β = 1.00 K−1 at two different τ values. The circles represent τ = 0.050 K−1, the triangles represent τ = 0.001 K−1, and the solid line represents the exact distribution.

Figure 11. Radial distribution for β = 1.00 K−1 at two different τ values. The circles represent τ = 0.050 K−1, the triangles represent τ = 0.001 K−1, and the solid line represents the exact distribution.

Figure 14. Angular distribution of θ for τ = 0.007 K−1 at two different β values. The circles represent β = 0.400 K−1, the triangles represent β = 1.0 K−1, and the solid line represents the exact distribution.

in the spirit of the PILE approach used for PIMD simulations of ring polymers. The PILE allows one to separate the sampling of the free particle action from the potential. It therefore allows one to avoid issues associated with the coupling between the normal modes of the free particle action, and interactions. These benefits have been observed in the original finite temperature PILE work15 and have now been confirmed for the case of open paths in the present LePIGS approach. These conclusions are substantiated by simulation results for model systems with single and double well interactions, as well as for realistic and accurate van der Waals interactions between water and molecular hydrogen. Of particular interest for LePIGS is the effect of γ(0), the centroid friction, which is a free parameter for the simulation that affects the rate of convergence to exact values by impacting the exploration of phase space. In principle, if γ(0) is too small, the canonical distribution will not be sampled, and if it is too large, then the exploration of phase space will be inefficient due to constant interactions with the random force. In practice, it is found that γ(0) = 0 still reproduces a canonical distribution of velocities for the centroid mode in the quartic oscillator, whereas for the free particle and harmonic oscillator it does not. This can be explained by the presence of the potential V(R), which acts to couple the energies of individual beads. This

Figure 12. Angular distribution of θ for β = 1.00 K−1 at two different τ values. The circles represent τ = 0.050 K−1, the triangles represent τ = 0.001 K−1, and the solid line represents the exact distribution.

system. The dynamics of the open polymer path variables is completed in normal mode space, taking advantage of the properties of the DCT to provide partial analytical propagation as well as analytical friction values for use in Langevin dynamics 7466



Article

(15) Ceriotti, M.; Parrinello, M.; Markland, T. E.; Manolopoulos, D. E. Efficient Stochastic Thermostatting of Path Integral Molecular Dynamics. J. Chem. Phys. 2010, 133, 124104. (16) Coalson, R. D.; Freeman, D. L.; Doll, J. Partial Averaging Approach to Fourier Coefficient Path Integration. J. Chem. Phys. 1986, 85, 4567. (17) Chakravarty, C. Particle Exchange in the Fourier Path-Integral Monte Carlo Technique. J. Chem. Phys. 1993, 99, 8038. (18) Chakravarty, C.; Gordillo, M.; Ceperley, D. A Comparison of the Efficiency of Fourier- and Discrete Time-Path Integral Monte Carlo. J. Chem. Phys. 1998, 109, 2123. (19) Doll, J.; Freeman, D. Comment on “A Comparison of the Efficiency of Fourier- and Discrete Time-Path Integral Monte Carlo” [J. Chem. Phys. 109, 2123 (1998)]. J. Chem. Phys. 1999, 111, 7685. (20) Chakravarty, C.; Gordillo, M.; Ceperley, D. Response to “Comment on ‘A Comparison of the Efficiency of Fourier- and Discrete Time-Path Integral Monte Carlo’” [J. Chem. Phys. 111, 7685 (1999)]. J. Chem. Phys. 1999, 111, 7687. (21) Mak, C. A Multigrid Algorithm for Sampling Imaginary-Time Paths in Quantum Monte Carlo Simulations. J. Phys. Chem. B 2004, 108, 6760. (22) Garcia, P.; Peinado, A. Diagonalizing Properties of the Discrete Cosine Transforms. Signal Process. 1995, 43, 2631. (23) Burnham, C. J.; Reiter, G. F.; Mayers, J.; Abdul-Redah, T.; Reichert, H.; Dosch, H. On the Origin of the Redshift of the OH Stretch in Ice Ih: Evidence from the Momentum Distribution of the Protons and the Infrared Spectral Density. Phys. Chem. Chem. Phys. 2006, 8, 3966. (24) Trotter, H. F. On the Product of Semi-Groups of Operators. Proc. Am. Math. Soc. 1959, 10, 545−551. (25) Bussi, G.; Parrinello, M. Accurate Sampling Using Langevin Dynamics. Phys. Rev. E 2007, 75, 5. (26) Hinsen, K. The Molecular Modeling Toolkit: a New Approach to Molecular Simulations. J. Comput. Chem. 2000, 21, 79−85. (27) Ing, C.; Hinsen, K.; Yang, J.; Zeng, T.; Li, H.; Roy, P.-N. A PathIntegral Langevin Equation Treatment of Low-Temperature Doped Helium Clusters. J. Chem. Phys. 2012, 136, 224309. (28) Frigo, M.; Johnson, S. G. The Design and Implementation of FFTW3. Proc. IEEE 2005, 93, 216−231 (special issue on “Program Generation, Optimization, and Platform Adaptation”). (29) Thirumalai, D.; Buskin, E. J.; Berne, B. J. An Iterative Scheme for the Evaluation of Discretized Path Integrals. J. Chem. Phys. 1983, 79, 5063. (30) Zeng, T.; Li, H.; Le Roy, R. J.; Roy, P.-N. Adiabatic-HinderedRotor Treatment of the Parahydrogen-Water Complex. J. Chem. Phys. 2011, 135, 094304. (31) Zeng, T.; Li, H.; Roy, P.-N. Simulating Asymmetric Top Impurities in Superfluid Clusters: a para-Water Dopant in paraHydrogen. J. Phys. Chem. Lett. 2013, 4, 18. (32) Cuervo, J.; Roy, P.-N. On the Solid and Liquid-Like Nature of Quantum Clusters in Their Ground State. J. Chem. Phys. 2008, 128, 224509. (33) Cuervo, J. E.; Roy, P.-N. Weakly Bound Complexes Trapped in Quantum Matrices: Structure, Energetics, and Isomer Coexistence in (para-H2)N(ortho-D2)3 Clusters. J. Chem. Phys. 2009, 131, 114302.

implies that a centroid friction of zero may be used to produce acceptable values for some systems. Work is currently underway to take into account a variety of trial wave functions. A properly optimized trial wave function should allow one to substantially reduce the projection time β/ 2, and hence the number of path integral beads. Simulations are also currently being performed to test the approach for parahydrogen clusters such as those studied in refs 32 and 33. The LePIGS method is general and now implemented in the MMTK. It can therefore be applied to a variety of systems for which ground state properties such as zero-point energies and structural distributions are sought.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), Ministry of Research and Innovation (MRI), Ontario, and the Canada Foundation for Innovation (CFI).

■

REFERENCES

(1) Meyer, H.-D. Studying Molecular Quantum Dynamics with the Multiconfiguration Time-Dependent Hartree Method. Comput. Mol. Sci. 2012, 2, 351. (2) Anderson, J. B. A Random-Walk Simulation of the Schrodinger Equation: H+3. J. Chem. Phys. 1975, 63, 1499. (3) Hammer, N. I.; Diken, E. G.; Roscioli, J. R.; Johnson, M. A.; Myshakin, E. M.; Jordan, K. D.; McCoy, A. B.; Huang, X.; Bowman, J. M.; Carter, S. The Vibrational Predissociation Spectra of the H5O2+RGn(RG = Ar,Ne) Clusters: Correlation of the Solvent Perturbations in the Free OH and Shared Proton Transitions of the Zundel Ion. J. Chem. Phys. 2005, 122, 244301. (4) McCoy, A. B.; Huang, X.; Carter, S.; Bowman, J. M. Quantum Studies of the Vibrations in H3O2- and D3O2-. J. Chem. Phys. 2005, 123, 064317. (5) Guardiola, R.; Navarro, J. A Diffusion Monte Carlo Study of Small Para-Hydrogen Clusters. Cent. Eur. J. Phys. 2008, 6, 33−37. (6) Boninsegni, M.; Moroni, S. Population Size Bias in Diffusion Monte Carlo. Phys. Rev. E 2012, 86, 056712. (7) Hetenyi, B.; Rabani, E.; Berne, B. J. Path-Integral Diffusion Monte Carlo: Calculation of Observables of Many-Body Systems in the Ground State. J. Chem. Phys. 1999, 110, 6143. (8) Sarsa, A.; Schmidt, K. E.; Magro, W. R. A Path Integral Ground State Method. J. Chem. Phys. 2000, 113, 1366. (9) Cuervo, J. E.; Roy, P.-N. Path Integral Ground State Study of Finite-Size Systems: Application to Small (Parahydrogen)N (N = 2− 20) Clusters. J. Chem. Phys. 2006, 125, 124314−124316. (10) Ceperley, D. M. Path Integrals in the Theory of Condensed Helium. Rev. Mod. Phys. 1995, 67, 279−355. (11) Parrinello, M.; Rahman, A. Study of an F Center in Molten KCl. J. Chem. Phys. 1984, 80, 860. (12) Tuckerman, M. E. In Quantum Simulations of Complex ManyBody Systems: From Theory to Algorithms; Grotendorst, J., Marx, D., Muramatsu, A. Eds.; NIC Series, John von Neumann Institute for Computing: Jülich, 2002; Vol. 10, pp 269−298. (13) Craig, I. R.; Manolopoulos, D. E. Quantum Statistics and Classical Mechanics: Real Time Correlation Functions from Ring Polymer Molecular Dynamics. J. Chem. Phys. 2004, 121, 3368. (14) Hall, R. W.; Berne, B. J. Nonergodicity in Path Integral Molecular Dynamics. J. Chem. Phys. 1984, 81, 3641. 7467


Langevin Equation Path Integral Ground State - ACS Publications

Recommend Documents