Self-Adaptive Quadrature and Numerical Path Integration

In the present paper we explore the use of generalized Gaussian quadrature methods in the context of equilibrium path integral applications. Using moment techniques, we devise a compact, self-adaptive approach for use in conjunction with selected classes of interaction potentials. We demonstrate that, when applicable, the resulting approach reduces appreciably the number of potential energy evaluations required in equilibrium path integral simulations.


I. INTRODUCTION
As evidenced by both the breadth and depth of recent applications, 1,2 numerical path integral methods are effective tools for the study of finite-temperature, many-body, quantum-mechanical systems.In their most complete form these methods make it possible to obtain refinable estimates of the equilibrium properties of complex physical systems directly from specified microscopic force laws without the need to introduce untestable numerical approximations.][4][5] Numerical path integral approaches, as well as their classical Monte Carlo analogs, 6 produce estimates for desired properties via surprisingly simple stochastic procedures.These procedures unify the treatment of quantum and classical problems within a common framework.Well suited to implementation on current and anticipated computational machines, they also benefit from a ''replica'' character that renders them inherently scalable. 7 key feature of path integral formulations is their use of dimensionality as a tool.In addition to the ''physical'' variables that label the relevant ''particles'' in a molecular application, one introduces ''auxiliary'' degrees of freedom to describe the quantum-mechanical dispersion of imaginary time ''paths'' along which the physical particles move.The number and details of the physical variables vary with the application.In contrast, the number of auxiliary variables is formally infinite even for a single quantum particle in one dimension.
Path integral simulations 8 typically consist of two separate steps, one formal and one numerical.The first step consists of using path integral quantum-mechanical expressions to create integral representations of specified thermodynamic quantities.Statistical methods are then invoked to transform these formal representations into practical, numerical ''estimators'' of the corresponding properties.Differences in convergence and/or efficiency makes the design and implementation of the various steps in this overall procedure an important practical issue in actual applications.
A general strategy to improve the performance of statistically based approaches is to incorporate auxiliary information into the numerical procedure.Such information can be from ''external'' sources or it can be obtained ''internally'' in a ''boot-strap'' fashion from the simulation itself.Examples of external information would be a prior knowledge of the asymptotic behavior of the solution or perhaps of rigorous sum rules that the solution must obey.0][11] As discussed elsewhere, 6,12 such information is typically the starting point for the development of an assortment of variance reduction strategies.
In the present work, we explore the use of a special class of adaptive quadrature techniques in the context of equilibrium path integration.By building path-related information into the quadrature method, we show that it is possible to formulate the path integral problem in such a way that, for restricted classes of problems, the number of evaluations of the system's potential energy is significantly reduced relative to conventional, a priori quadrature approaches.
The organization of the remainder of the present paper is as follows: In Sec.II we present the basic formal developments and discuss the resulting computational procedures.We apply the resulting methods to two prototypical condensed phase applications in Sec.III to explore and illustrate their computational merits.Finally, we summarize our findings and offer speculation concerning their possible future developments in Sec.IV.

A. Fourier path integral methods
Feynman's ''sum over histories'' approach 13,14 provides a convenient formal starting point for deriving integral representations of equilibrium thermodynamic properties.In the present context, we implement this strategy using Fourier techniques.Such a formulation has the feature that it explicitly separates the algorithmic task of path enumeration from that of path averaging.Given that the details of the Fourier method are described elsewhere, 15,16 we limit the present discussion to those topics necessary to introduce the new developments.For simplicity, we utilize one-dimensional notation throughout, indicating multidimensional details only where necessary.
The canonical ensemble density matrix, ͗xЈ͉ ϫexp(Ϫ␤H)͉x͘, provides the information necessary to compute various equilibrium properties.This object, or a suitable average of it, is thus a representative computational objective ''typical'' of general thermodynamic simulations.In the Feynman approach 13,14 the density matrix is written formally as ͗xЈ͉exp͑Ϫ␤H͉͒x͘ϭ ͚ paths e ϪS͓path͔ , ͑2.1͒ where the summation denotes a ''sum over paths,'' x(), that begin at an ''initial'' point, x, at a ''time'' ϭ0 and end at the ''final'' point, xЈ, at a time ϭ␤ប.The statistical significance of each path is governed by the Boltzmann-type factor, exp(ϪS), in which the functional, S, is given in terms of the path by In order to implement these expressions in practice, one must devise procedures for enumerating the relevant paths and for constructing the formal sums over them.
To label the paths, it is convenient to rewrite Eq. ͑2.2͒ by switching to a scaled time variable, u, defined as u ϭ/(␤ប), and to express the paths as 15 x͑u,a ͒ϭxϩ͑ xЈϪx ͒uϩ ͚ kϭ1 ϱ a k sin͑ku ͒.

͑2.3͒
In Eq. ͑2.3͒ we are using a ''reference path,'' here taken to be a simple, linear expression, to build in the proper uϭ0 and uϭ1 boundary conditions.By construction, therefore, the ''fluctuations'' of the path about this reference must vanish at uϭ0 and uϭ1 and can be written as a ͑formally infinite͒ Fourier sine series.As emphasized by Hamming, 17 the use of such a ''reference path'' approach is a general and convenient device for accelerating the convergence rates in Fourier-based applications.Generalizations of the abovementioned results for multidimensional applications and/or for more elaborate reference paths are straightforward.
With the underlying paths labeled, the remaining task is to evaluate the path summation that appears in Eq. ͑2.1͒.Substituting Eq. ͑2.3͒ into Eq.͑2.2͒ shows that the ''action'' for a particular path, x(u,a), is given in terms of the path variables xЈ, x, and the path's Fourier coefficients, ͕a k ͖, by where and where du V͑x͑u,a͒͒.

͑2.6͒
Using Eqs.͑2.4͒-͑2.6͒, it is easy to show that the ratio of the density matrix to its free-particle counterpart, ͗xЈ͉ ϫexp(Ϫ␤H)͉x͘ FP , can be written as Equation ͑2.7͒ is an infinite-dimensional integral representation of our original computational goal, the quantummechanical density matrix.This expression, either directly or with slight embellishment, is the typical starting point for the Monte Carlo construction of general thermodynamic averages.Depending upon the application, explicit construction of the underlying density matrix elements is unnecessary.It should be noted from Eqs. ͑2.5͒ and ͑2.7͒ that Fourier methods introduce a natural internal length scale into the path integral problem.[20]

B. Quadrature convergence characteristics
In path integral applications, a principal task is the evaluation of the path average of the potential energy, Eq. ͑2.6͒.In most problems this path average is not available analytically and instead must be evaluated with numerical quadrature.It is critical to note that even when the system itself has many degrees of freedom, the path average in Eq. ͑2.6͒ remains one dimensional.Consequently, Eq. ͑2.6͒ can, in general, be evaluated using conventional quadrature.A key practical issue, and a principal concern of the present study, is the efficiency of the associated quadrature method.
To address this concern we examine the magnitude of the errors incurred by various approximate quadrature formulas.As shown in the following we find that the dependence of the error on the number of time-domain quadrature points depends on the nature of the frequency spectrum of V(x(u)).Specifically, if V(x(u)) has significant strength at all frequencies, then the convergence of the quadrature error is governed by aliasing effects 21 associated with the high frequency components of V(x(u)).This is significant because the errors introduced by these aliasing errors turn out to be largely independent of the details of the quadrature method employed.In contrast, if V(x(u)) is band limited, then the quadrature error reflects explicitly the details of the particular numerical method involved ͑i.e., is governed by standard convergence criteria͒.
We begin by evaluating the u-space integration for a specified path in Eq. ͑2.6͒ with an ordinary N-point quadrature approximation, 21 V ¯N , defined by The points and weights for this otherwise unspecified quadrature are denoted in Eq. ͑2.8͒ by ͕W n ͖ and ͕u n ͖, re- spectively.It is important to emphasize that in Eq. ͑2.8͒ the points and weights are fixed, a priori, and are independent of the Fourier variables that define the path, x(u,a).As with the quantum-mechanical path itself, the potential energy as a function of time can be written as the sum of ''reference'' and ''fluctuation'' terms.Specifically, for paths defined by Eq. ͑2.3͒, the corresponding potential energy along the path can be written as

͑2.9͒
The leading terms in Eq. ͑2.9͒ build in the proper boundary conditions at uϭ0 and uϭ1 while the remaining term describes deviations about this reference.These deviations vanish, by construction, at uϭ0 and uϭ1 and are thus naturally described by a Fourier sine series.The Fourier coefficients for the potential energy fluctuations along the path, ͕V ˆk͖, are generally nonlinear functions of the underlying path coefficients, ͕a k ͖.From Eqs. ͑2.6͒ and ͑2.9͒, the path average of the potential energy, V ¯, becomes where the time averages of the sine terms, S k , are given analytically as Combining Eqs.͑2.8͒-͑2.10͒,we see that the error, ␦V ¯N , in a general N-point quadrature approximation to the timedomain path average is given ͑assuming the quadrature errors for the reference terms vanish͒ by Here ␦S k,N is the difference between S k and the corresponding N-point approximation to the time-domain average of the kth Fourier sine term, S k,N , defined by In discussing the quadrature error in the path average of the potential it is useful to identify two limiting situations.If one dealing with a relatively ''smooth'' potential energy function and if the number of Fourier path variables that contribute to the underlying quantum-mechanical paths is finite, then the potential energy's time-domain signal, Eq. ͑2.9͒, is ''band limited.''That is, it contains no significant strength at frequencies beyond some maximum value.The asymptotic convergence of the numerical estimate for the path average in such applications is governed by the details of the underlying quadrature method.On the other hand, if one is dealing with appreciably nonlinear potentials and/or if arbitrarily high-order Fourier terms contribute significantly to the underlying quantum-mechanical paths, then the potential energy generally contains signal strength at all frequencies.For any finite, u-domain quadrature, such highfrequency components lead to inevitable ''aliasing'' effects 21 and to quadrature errors that are largely independent of the particulars of the specific method employed.To quantify these general remarks, it is useful to partition Eq. ͑2.12͒ into ''low-order'' and ''high-order'' components by writing ͑2.14͒ If we are dealing with a band-limited situation, then we can always select a value of N that is sufficiently large that for kϾN the Fourier coefficients, ͕V ˆk͖, are effectively zero.
Hence, only the first summation on the right-hand side of Eq.
͑2.14͒ is important.The quadrature errors, ␦S k,N , in the loworder terms in this summation thus reflect the details of the numerical method involved.For example, evaluating the summation in Eq. ͑2.13͒ analytically for the trapezoidal method yields From Eqs. ͑2.11͒, ͑2.14͒, and ͑2.15͒ we see that the trapezoidal quadrature error for band-limited applications ͓i.e., the low-order portion of Eq. ͑2.14͔͒ scales for large N as 1/N 2 , as expected. 21Improving this asymptotic convergence rate for band-limited applications is a simple task: One needs merely to replace the trapezoidal quadrature with a higherorder method.For non-band-limited applications, however, the situation is more complex.By definition, the time dependence of the potential energy in such situations contains frequency components of arbitrarily high order.A nontrivial portion of the signal strength in Eq. ͑2.9͒ thus lies beyond the frequency threshold imposed by the constraints of finite u-domain sampling.Unlike band-limited applications, we can not, therefore, improve the convergence rate qualitatively by simply switching to a ''better'' quadrature.The best we can generally hope to accomplish instead is to remove the low-order components of Eq. ͑2.14͒.Even if one succeeds in eliminating all such low-order terms, however, ''aliasing'' errors in the high-frequency (kϾN)␦S k,N terms of Eq. ͑2.14͒ will remain.
We can obtain a rough estimate of the magnitude of high-order aliasing errors if we assume that the corresponding quadrature errors are uncorrelated Gaussian random vari-ables of order unity.This, and the assumption that the errors in the low-order terms in Eq. ͑2.14͒ have been reduced to negligible levels, implies that the root mean square quadrature error in the path average of the potential energy is given by

͑2.16͒
As discussed by Lanczos, 22 the Fourier sine coefficients of a ''smooth'' function defined on a finite interval decay asymptotically as the reciprocal of the cube of the Fourier index.Equation ͑2.16͒ thus suggests ͑heuristically͒ that timedomain quadrature error for non-band-limited path averages should converge asymptotically as roughly 1/N 5/2 .Quadrature errors in the high-order sine averages, while typically of order unity, are not, as assumed in Eq. ͑2.16͒, completely uncorrelated.Nonetheless, the present argument is sufficient to suggest that in non-band-limited applications the asymptotic convergence rates for time averages of the potential do not trivially mirror the details of the quadrature involved but instead are controlled by universal aliasing errors.

C. Self-adaptive quadrature methods
A large portion of the computational effort in path integral simulations often involves the calculation of the system's potential energy.In such situations it is thus useful to minimize the number of such evaluations.a restricted, but important class of interactions, switching to an adaptive, coordinate-space quadrature as opposed to the fixed, timedomain approach used in Eq. ͑2.8͒ permits such a reduction.As an added bonus, such adaptive quadrature will permit us to break free of the convergence limits imposed by the use of time-domain methods.
If we wish to evaluate the path average of the potential energy along a path, x(u,a), arising from fixed end points, x and xЈ, and a given set of Fourier coefficients, ͕a k ͖, then it is useful to do so in the coordinate domain as opposed to a time domain.That is, rather than averaging the potential energy in time along the specific path x(u,a) in the manner indicated in Eq. ͑2.8͒, we imagine first using x(u,a) to construct P(x u ,a), the probability distribution function of positions, x u , associated with the path.In terms of this distribution of positions, we then write the path average of the potential energy as the coordinate-space average V ¯ϭ ͵ P͑x u ,a ͒V͑ x u ͒dx u .

͑2.17͒
If we actually had to construct the entire probability distribution, P(x u ,a), Eq. ͑2.17͒ would offer little practical advantage over the corresponding time-domain formulation.However, averages of the type defined by Eq. ͑2.17͒ can be evaluated efficiently using generalized Gaussian quadrature knowing only moments of the distribution rather than the distribution itself.As discussed by Gordon 23 and Press et al., 21 given 2N power moments, ͕ n ͖, defined in terms of the probability distribution function of positions by we can construct N generalized Gauss quadrature points and weights, ͕x u,n ͖ and ͕w n ͖, in terms of which the path average can be written The power moments required to implement this procedure can be expressed as and hence can be obtained via u-domain quadrature methods.
The key practical point to be made is that this formulation of the path average of the moments involves numerous evaluations of the path, x(u,a), not the potential energy, V(x(u,a)).Fast Fourier transform-sine methods can be used to evaluate the sums in Eq. ͑2.3͒ and thus provide an efficient means for obtaining the required path information.Finally, we note that modified moment approaches are available and appear to be of practical significance in stabilizing the numerics of the moment problem. 24he advantage of the coordinate-space formulation, Eq. ͑2.19͒, over the time-domain expression, Eq. ͑2.8͒, is that the coordinate dependence of the potential energy in many molecular applications is often a relatively ''simple'' function of position whereas its time variation is typically appreciably more complex.Consequently, the number of coordinatespace quadrature points required for an accurate evaluation will generally be far fewer than the corresponding number required for the time-domain evaluation.If, for example, the potential energy is a polynomial of order 2NϪ1 or less, then N coordinate-domain quadrature points are sufficient to evaluate the required path average exactly whereas with time-domain quadrature an infinite number of points would be required.
Before considering explicit numerical applications, it is useful to examine the simplest nontrivial form of the present results.If we include only a single Gauss point ͑two power moments͒ in our adaptive quadrature, then the resulting scheme approximates the path average of the potential by the value of the potential energy at the corresponding path average position.As discussed by Feynman, 13 the path average of the potential, Eq. ͑2.6͒, in this approximation becomes where, from Eq. ͑2.20͒, ͗x͘ is the path centroid.The present developments provide a systematic way to generalize Eq. ͑2.21͒ to include information about successively high-order moments into the construction of the path average.Using such expressions, we can, in fact, rewrite the original path integral in terms of integrals over the moments or cumulants of the associated paths.Evaluating the resulting expressions via molecular-dynamical methods leads to an effective dynamics reminiscent of centroid-based path integral developments [25][26][27] in which the natural variables of the problem are the path moments ͑or cumulants͒.
We close this section by noting that adaptive quadrature approaches can be utilized for the calculation of more gen-eral thermodynamic properties.To illustrate this point, we consider a representative property, the total energy, ͗E͘.Us- ing the quantum-mechanical virial estimator, 8 the total energy is expressed as ͗E͘ϭ 1 2 ͗x͑u͒VЈ͑x͑u͒͒͘ϩ͗V͑x͑u͒͒͘.

͑2.22͒
The ''barred'' terms in Eq. ͑2.22͒ are the time-domain averages of the various quantities along specified quantum paths while the brackets represent ensemble averages over a thermal distribution of paths and particle position.It is straightforward to express the time-domain averages in the integrands on the right-hand side of Eq. ͑2.22͒ as analogous coordinate-domain averages.The second term on the righthand side of Eq. ͑2.22͒ is, in fact, the potential energy term considered previously ͓cf.Eq. ͑2.10͔͒.The first term, the kinetic energy portion of the total energy, can be rewritten in coordinate-domain form by simply replacing the timedomain product of the coordinate and gradient of the potential by the corresponding coordinate-domain expression.
From the preceding discussion it is apparent that the present developments can be used to construct adaptive quadrature approaches for the path integral treatment of formally one-dimensional problems.The approach is also clearly viable for those portions of the interaction potential that can be written as the sum of pair-type interactions.In Sec III, for example, we demonstrate that such methods can be applied to pair-potential models of many-body systems.While not all inclusive, such models constitute an important class of practical applications.

III. NUMERICAL EXAMPLES
We now consider a number of numerical examples to illustrate the developments of Sec II.These include a documentation of the nature of aliasing errors described in Sec.II, and the application of adaptive quadrature methods to problems representative of those that arise in ''typical'' quantum fluid simulations.
As discussed in Sec.II, quadrature errors in path averages of the potential energy for band-limited applications (NϾk) reflect the details of the numerical procedure involved.Figure 1 shows quadrature errors for Fourier sine averages, ␦S k,N ͓cf.Eqs.͑2.13͒ and ͑2.11͔͒, obtained using Gauss-Legendre and trapezoidal quadrature for a fixed Fourier index, k.These results are displayed as a function of the number of quadrature points, N. We see that the errors for the trapezoidal quadrature converge to zero as the reciprocal of the square of the number of quadrature points, 1/N 2 , for NϾk.This asymptotic rate, anticipated from standard discussions, 21 would become 1/N 4 if the trapezoidal method were replaced by Simpson's rule, 1/N 6 if replaced by Bode's rule, and so on.Figure 1 documents the exceptionally rapid asymptotic convergence rates achieved with Gaussian quadrature.
Unlike the band-limited case, however, aliasing effects are inevitable with fixed, u-domain quadrature for highfrequency signal components.Figure 2 shows the quadrature errors for Gauss-Legendre Fourier sine averages, ␦S k,N , for fixed N as a function of the Fourier index, k.These errors, effectively zero for low-order terms, become an erratic function of order unity for the high-order Fourier components.As discussed in Sec.II, such a behavior of the quadrature error produces an overall convergence rate for the path average of the potential energy that is roughly 1/N 5/2 .The trapezoidal approach, as can be seen from Eqs. ͑2.11͒ and ͑2.15͒, has a poor low-order performance that, by itself, produces an overall 1/N 2 convergence rate in calculated path averages of the potential energy.High-order aliasing errors in the trapezoidal method, although more systematic than those of Gaussian quadrature, are also of order unity.While it is straightforward to improve upon the low-order quadrature performance of the basic trapezoidal approach by combining results on different mesh sizes, 21 the high-order errors in such methods remain of order unity, thus limiting the asymptotic convergence rate of the calculated path average of the potential energy regardless of the type of quadrature method used.We now consider a simple example to illustrate the effectiveness of adaptive quadrature methods for the simulation of the equilibrium behavior of a simple pair-potential model of a quantum fluid.As our example we select a Lennard-Jones dimer whose parameters chosen to be those of helium ͑⑀ϭ10.22K, ϭ2.556 Å͒.With the atoms locked in place at a specified ͑classical͒ separation distance, r, we sample the three-dimensional quantum mechanical paths for the interacting particles from an appropriate quantummechanical distribution at a given temperature, T. For each path in the resulting statistical ensemble, we evaluate the average of the potential energy, Eq. ͑2.6͒, using both conventional, time-domain trapezoidal quadrature ͓Eq.͑2.8͔͒ and coordinate-domain, adaptive Gauss quadrature ͓Eq.͑2.10͔͒.We then compare those results with the ''exact'' value, itself obtained using high-accuracy Romberg methods.
In the present work we take Tϭ51.1 K, a temperature corresponding to one of the thermodynamic states examined previously by Ceperley and Pollock. 28The required moments of the path are computed using high-order, u-domain quadrature.Adaptive Gauss points are then obtained from these path moments using Gordon's method. 23Both Gauss-Legendre and trapezoidal quadrature are used for the required power moments.Although power moment methods are sufficient for the present study, it should be noted that more stable, modified moment approaches 21,24 may prove useful in more general applications.
Table I lists the average number of potential energy evaluations required to achieve a preselected level of accuracy in the path average of the potential energy using timedomain trapezoidal and coordinate-domain, adaptive Gauss quadrature.Results are shown as a function of the classical dimer separation distance and as a function of the number of Fourier terms in the associated paths.The numbers listed in Table I are obtained using ensembles that contain up to 1000 randomly selected paths.We see that adaptive quadrature reduces significantly the number of potential energy evaluations required compared with trapezoidal approaches.At smaller separation distances, where the paths are exploring the harshly repulsive, nonlinear portions of the potential, the reduction is roughly 100-fold.At larger separation distances, where the paths visit portions of the potential energy that vary less rapidly with distance, the reduction, although measurable, is less dramatic.
A more detailed analysis of the results in Table I offers insights into the nature of the potential energy and the contrasting natures of fixed and adaptive quadrature.From the properties of Gaussian quadrature, we know that an N-point scheme is exact for the integration of polynomials of order 2NϪ1 or less.Reversing this argument, the average order of the coordinate-domain quadrature required to obtain acceptable levels of accuracy for the path average of the potential energy serves as a ''signature'' of the local polynomial order of those portions of the potential energy visited by the associated quantum-mechanical paths.From Table I, we see that the local polynomial order of the potential decreases with increasing classical separations.At small separation, between four and five adaptive Gauss points are required to obtain an accurate evaluation of path averages of the potential, an indication that the potential energy varies rapidly in these harshly repulsive regions.Near the classical equilibrium, on the other hand, the number of adaptive Gauss points required drops to less than two, implying that over the length scale of thermal fluctuations the potential is weakly anharmonic.Finally, at large separation distances, where the potential is effectively linear on the length scale of path fluctuations, a single adaptive quadrature point is sufficient.
Local polynomial order provides an explanation of the interesting k max dependence of the results of Table I.We see in Table I that increasing the number of path variables in the present study from k max ϭ1 to k max ϭ8 leads to an increase in the number of time-domain quadrature points required to obtain an accurate path average.On the other hand, we see that as we add quantum-mechanical character to the paths, the number of adaptive coordinate-domain quadrature points required decreases at small classical separation distances and increases ͑slightly͒ for larger separations.For the present problem, increasing the number of path variables from k max ϭ1 to k max ϭ8 produces paths that are spatially more extended.From a time-dependent point of view, as the path becomes more extended, the potential energy, V(x(u,a)), becomes a more complex function of time and more time-TABLE I.The average number of potential energy evaluations required to achieve a mean absolute error of less than 0.01⑀ in the path average of the potential energy for two Lennard-Jones helium atoms using u-domain trapezoidal ͑U͒ vs self-adaptive Gauss quadrature ͑S͒.Results are shown as a function of the classical separation between particles, r/.The results are calculated from averages over of up to 1000 paths chosen randomly from a quantum-mechanical distribution at Tϭ51.1 K for varying numbers of Fourier path variables (k max ).domain quadrature points are required to obtain an accurate path average.In contrast, as paths become more extended the particles involved are free to move further from their classical positions.Depending on the classical separation distance involved, this can result in either an increase or a decrease in the number of coordinate-domain quadrature points required.
For small classical separation distances, extending the path permits it to distort toward the potential minimum away from the harsher, higher-order repulsive portions of the interaction, thus reducing the number of adaptive, coordinatedomain points required.For classical separation distances near the potential minimum, on the other hand, adding quantum mechanical dispersion permits the paths greater access to repulsive, locally higher-order portions of the potential, thereby increasing ͑slightly͒ the number of adaptive quadrature points required.
Tables II and III present the average absolute errors in the calculated path averages of the Lennard-Jones dimer potential energy.Results for time-domain, trapezoidal quadrature are listed in Table II while the corresponding results for coordinate-domain adaptive quadrature are shown in Table III.Errors for both types of quadrature are presented as a function of the classical separation.As with Table I, the results summarized in Tables II and III are generated from an ensemble of up to 1000 paths chosen randomly from the appropriate quantum-mechanical distribution.As expected, in Table II we see that the average error in the trapezoidal results for a given separation distance decrease asymptotically as the square of the number of quadrature points.On the other hand, we see from Table III that the convergence of the adaptive quadrature is significantly more rapid than observed for the trapezoidal results.
As a final study of their utility, we consider the application of adaptive quadrature methods to a many-particle quantum-mechanical system, the Lennard-Jones model of the (H 2 ) 22 cluster.This system has been characterized extensively 29,30 and is thus a convenient benchmark application.
We begin by considering the quantum-mechanical density of a related problem, the simple Lennard-Jones ''cage'' potential. 31In this model a particle of mass, m, moves in the one-dimensional potential, V cage (x), given by where V(x) is a specified pair potential.In the present application, we take V(x) to be the Lennard-Jones potential appropriate for H 2 -H 2 interactions ͑⑀ϭ34.2K, ϭ2.96 Å͒ 29 As discussed elsewhere, 30,31 we have found this cage model to be useful in assessing the performance and/or convergence of various path integral algorithms.In Table IV we list the quantum-mechanical density at a particular location (xϭ0) as a function of the number of adaptive Gauss points used to evaluate the path average of the potential for varying numbers of Fourier path variables, k max .The cage parameter, a, in these simulations is taken to be 1.2, a value that produces interparticle separation distances of roughly those observed in the molecular cluster.The necessary path moments are The calculated density ratio converges, within the estimated single standard deviation statistical error of Ϯ0.6, to the exact value of 134.17, obtained using NMM methods, 33 with only five quadrature points.This rate of convergence is consistent with that observed in our dimer study ͑Table III͒ and, as discussed in Sec.II, is appreciably more rapid than would be possible with fixed, time-domain methods.Finally, to examine the utility of self-adaptive methods for more physically demanding applications, we consider their use in the study of many-body molecular clusters.As a representative simulation, we present in Table V results for the thermodynamic potential energies of the (H 2 ) 22 cluster obtained for varying numbers of path variables at Tϭ6 K.In Table V the average potential energies are obtained via a path-average estimator 32 using adaptive quadrature methods.For a given number of path variables, results obtained using various numbers of quadrature points are compared with the corresponding results obtained from direct, time-domain Fourier path integral methods. 30As with the corresponding cage simulations, the convergence of the self-adaptive results is extremely rapid.Specifically, for a given number of path variables only four self-adaptive Gauss quadrature points are required to produce convergence of the thermodynamic potential energy to its proper limit.As discussed in Sec.II, conventional Fourier methods typically require that the number of time-domain quadrature points be at least of the order of the number of path variables in the problem.

IV. SUMMARY
In the present work we have attempted to clarify the role of the choice of numerical quadrature on the convergence properties of numerical path integration algorithms.We have shown that, for restricted classes of interaction potentials, Gauss moment methods are feasible.These self-adaptive, coordinate-domain methods break free of the limits on the convergence rates of quadrature error otherwise imposed by fixed, time-domain quadrature.When applicable, these methods appear to reduce dramatically the number evaluations of the potential energy required for typical numerical path integral applications.

FIG. 1 .
FIG. 1. Quadrature errors in the approximation to ͐ 0 1 sin(ku)du.Results are shown as a function of the number of quadrature points ͑N͒ for a fixed value of k ͑taken here as 257͒ for both Gauss-Legendre ͑closed circles/dashed line͒ and trapezoidal quadrature ͑solid line͒.For reference, the horizontal line indicates an error of zero.

TABLE II .
The mean absolute errors in the path averages of the Lennard-Jones potential energy ͑in units of ⑀͒ for two helium atoms separated ͑clas-sically͒ by a distance r/.Results, shown as a function of the number of u-domain trapezoidal points, N u , are obtained by comparing approximate and numerically exact path averages for up to 1000 paths chosen at random from a quantum-mechanical distribution at 51.1 K.The number of path variables is constant for all cases (k max ϭ8).The numbers indicated in parentheses are the exponents of the corresponding entries.By construction, N u ϭ1 values are the classical results.The minus sign ͑Ϫ͒ indicates that the associated errors are less than 10 Ϫ5 ⑀.

TABLE III .
The mean absolute errors in the path averages of the Lennard-Jones potential energy ͑in units of ⑀͒ for two helium atoms separated ͑clas-sically͒ by a distance r/.The format is that of Table II except that selfadaptive Gauss moment quadrature methods are used to calculate the path averages.Entries designated by a minus sign ͑Ϫ͒ indicate that the associated errors are less than 10 Ϫ5 ⑀.Note the large difference in scale of the number of self-adaptive points, N s , used in the present results and the number of trapezoidal points, N u , used in TableII.

TABLE IV .
32sted is the ratio of the density at xϭ0, ͑0͒, to the freeparticle value, f p , for the Lennard-Jones cage potential of Sec.III.The interaction potential parameters are representative of the H 2 molecular interaction ͑⑀ϭ34.2K,ϭ2.96Å͒while the cage parameter, a, is taken as 1.2.As discussed in the text, the density for a given number of adaptive quadrature points, N s , is obtained by extrapolating values obtained from a sequence of calculations with varying numbers of Fourier path coefficients.On the order of 5ϫ106Monte Carlo points are used in each path integral simulation.Errors in the calculated density ratios are single standard deviation values.The exact density ratio for this system ͑computed by NMM methods͒ is 134.17.͑2.20͒ using high-accuracy Romberg methods.These results are computed from Eq. ͑2.7͒ via Fourier path integral Monte Carlo methods.From previous work32we know that the calculated density for nonpartial averaged Fourier path integral calculations converges to its exact value asymptotically as 1/k max .Results listed in Table IV are obtained from a fit of the k max dependence of the calculated density ratio to a second-order polynomial in 1/k max .Up to five million Monte Carlo configurations are used in each individual path integral simulation.

TABLE V .
Listed are the calculated average potential energies for a Lennard-Jones model of the (H 2 ) 22 cluster.The interaction potential parameters are those of Ref. 29 ͑⑀ϭ34.2K,ϭ2.96Å͒.Results for various numbers of Fourier path variables (k max ) are shown as a function of the number of adaptive quadrature points, N s used in the simulation.Convergence to the proper limit ͑FPI͒ for each k max value is rapid and is similar to that achieved for the model cage results of TableIV.One million Monte Carlo points are used in each path integral simulation.Errors in the calculated values are single standard deviation values.