A heat capacity estimator for Fourier path integral simulations

Previous heat capacity estimators useful in path integral simulations have variances that grow with the number of path variables included. In the present work a new specific heat estimator for Fourier path integral Monte Carlo simulations is derived using methods similar to those used in developing virial energy estimators. The resulting heat capacity estimator has a variance that is roughly independent of the number of Fourier coefficients (kmax) included, and the asymptotic convergence rate is shown to be proportional to 1/kmax2 when partial averaging is included. Quantum Monte Carlo simulations are presented to test the estimator using two one-dimensional models and for Lennard-Jones representations of Ne13. For finite kmax, using numerical methods, the calculated heat capacity is found to diverge at low temperatures for the potential functions studied in this work. Extrapolation methods enable useful results to be determined over a wide temperature range.


I. INTRODUCTION
Monte Carlo ͑MC͒ 1 and molecular dynamics techniques 2 have been essential in the understanding of many particle systems at nonzero temperature.These methods allow the implementation of microscopic statistical mechanics models without introducing uncontrolled approximations.4][5][6] A large number of these algorithms uses the path integral method, inspired by Feynman's formulation of statistical mechanics. 7n the present work we use the Fourier path integral ͑FPI͒ method, 3,4 in which paths are expressed as a Fourier series.In this method, the extremes of the path are two points in configuration space, and the Fourier coefficients are extra degrees of freedom associated with the quantum character of the problem.Because the number of Fourier coefficients required is infinite, for practical calculations the Fourier series is truncated by including k max terms.
Several estimators of the energy have been developed in the FPI method. 3,8Each of these estimators gives the same average energy in the limit of an infinite number of Fourier coefficients, but with differing levels of statistical noise in a simulation.One approach, often called the T-method, 3 is obtained by differentiating the canonical partition function with respect to the temperature.Another technique, called the H-method, 3 is obtained by direct application of the quantum Hamiltonian operator on the density matrix.A third method based on the virial theorem is discussed below.The T-method estimator is easy to compute, consisting primarily of terms needed during the quantum random walk.However, the variance of the energy in a T-method calculation can be shown to increase with k max .The H-method gives an energy with a variance that is, at most, weakly dependent on k max , but the terms in the estimator include derivatives of the potential and integrations with respect to the time variable that can be difficult and expensive to evaluate.
Estimators for the heat capacity of quantum systems are obtained by differentiating the expression for the average energy using either the T-method or the H-method with respect to temperature.The resulting estimators have been called the TT-method or the TH-method, 9 depending on the estimator for the energy used.The differentiation with respect to temperature results in an estimator with a variance that increases with another factor of k max .Consequently, the variance in the heat capacity of the TT-method grows as k max 2 , and the variance in the heat capacity in the TH-method grows as k max . 9s with discretized path integral methods, 6,10 it is possible to develop an energy estimator, based on the T-method result, where the terms leading to an ill-behaved variance are explicitly eliminated. 8This estimator has been called the virial energy estimator 10 ͑hereafter referred to as the V-method estimator͒, because it is based on the virial theorem.The convergence rate of the virial energy estimator is identical to the convergence rate when the T-method is used to evaluate the energy, but the variance of the virial result is only weakly dependent on k max .As with the T-method and the H-method, the heat capacity in the V-method can be evaluated by direct temperature differentiation.As a result of this additional temperature differentiation, the resulting estimator for the heat capacity has a variance that grows linearly in k max .In this paper we use methods analogous to those used in deriving the V-method estimator for the energy to derive an estimator for the heat capacity whose variance is, at most, weakly dependent on k max .We illustrate the estima-tor with applications to some one-dimensional model problems, and to quantum Ne 13 .In addition to its intrinsic interest, the Ne 13 calculation uses the first application of the extrapolation methods outlined in Ref. 8.
The remainder of this paper is organized as follows: in Sec.II we briefly review the basic expressions used in the FPI method and derive an estimator for the specific heat.We also demonstrate how to remove the terms that cause the variance of the heat capacity to increase with increasing k max .In Sec.III we present some sample results for two one-dimensional model systems and a detailed calculation of the heat capacity for quantum Ne 13 modeled by Lennard-Jones forces.We summarize our results in Sec.IV.

A. The FPI method
Although the FPI method has been explained in detail elsewhere, 3 it is useful to begin the present discussion with an outline of the key ideas and formulas.In that way the derivation of the heat capacity estimator can be made clear.In what follows, we let d be the total dimensionality of the system.For example, for an N particle system in three dimensions, dϭ3N.In the following, we assume all particles have mass m, the extension of the formulas to many particles of differing masses being straightforward.
In path integral approaches to quantum statistical mechanics, the quantum density matrix at temperature T is expressed as a path integral in imaginary time ͑x,xЈ;␤ ͒

͑2͒
where x and xЈ represent the particle coordinates in d dimensions, V(x) is the potential energy, and ␤ϭ1/k B T with k B the Boltzmann constant.In the FPI method the paths are expanded in a Fourier sine series about a straight line path connecting x to xЈ where u is the reduced time variable uϭ/(␤ប) that ranges between 0 and 1, and a k is a vector containing the d components of the kth Fourier coefficient associated with the d-dimensional coordinate.The result of this sine expansion is to replace the integration over all paths by an ordinary Riemann integral over all the Fourier coefficients.In practical applications it is usual to truncate the infinite summation in Eq. ͑3͒ after the first k max coefficients.The resulting approximate path is written Using Eq. ͑4͒, it is possible to show 3 the resulting expression for the density matrix, given by where S k max is the action and the line above V in Eq. ͑6͒ represents a path average The notation a k 2 implies the sum of the squares of all components of a k , and superscripts ( f p) like those appearing in Eq. ͑5͒ imply the free particle result where the potential energy is set to 0. Additionally, a denotes the truncated set of k max Fourier coefficients.
In this paper, where we focus on the energy and the heat capacity, only the diagonal elements of the density matrix are required.In such cases, it is possible to express the expectation value of any coordinate-dependent quantity A(x) by the expression To find the exact expectation value of A in Eq. ͑9͒, it is necessary to determine the limit as k max →ϱ.It has recently been shown 8 that expectation values, like those appearing in Eq. ͑9͒, converge asymptotically to the exact expression as 1/k max .To improve the convergence rate, it is possible to use partial averaging 3,8 which approximately includes contributions from the ignored terms in the Fourier series.In partial averaging, the potential energy appearing in Eq. ͑6͒ is replaced by an effective potential of the form ͑using the gradient approximation, see Ref. 8͒ Using the partial averaged action to replace the action in Eq.͑9͒, coordinate dependent expectation values converge asymptotically to the exact result as 1/k max 2 . 8In the derivations that follow, partial averaging is assumed.The corresponding primitive FPI results can be obtained by replacing V (pa) by V wherever it appears.

B. The heat capacity estimator
As mentioned in the Introduction, different energy estimators have different properties, even though their values in the limit k max →ϱ are the same.In this work we restrict attention to the V-estimator, 8 because the V-estimator expression is not costly to evaluate and has a variance that does not grow with increasing k max .In deriving the expression for the V-estimator for the energy, 8 the terms whose variance do grow with k max can be shown to cancel and are explicitly removed.In what follows, we apply the same technique to develop a heat capacity estimator that is also well-behaved.The resulting expression is valid only for systems with potential energies having well-defined derivatives through fourth order.
To begin, we define which we call the virial operator.In terms of the virial operator, we express the total energy in the V-method by 8

͑14͒
It has been proved 8 that the V-estimator of the energy converges asymptotically to the exact result as k max Ϫ2 ; i.e., ͗E͘ ϱ ϭ͗E͘ k max ϩO͑k max Ϫ2 ͒᭙T, ͑15͒ as k max →ϱ.
The expression for the specific heat C V can be obtained by temperature differentiation of the total energy V-estimator Rewriting Eq. ͑16͒, we obtain

͑17͒
where and In deriving Eq. ͑19͒, we have used the relation we find

͑25͒
The estimator for the heat capacity given in Eq. ͑24͒ has a variance that increases linearly with k max .The variance behavior is a consequence of the last two terms on the righthand side of Eq. ͑24͒.Using the same strategy as in the derivation of the V-method for the energy, we reformulate Eq. ͑24͒ in such a way that the terms that otherwise produce the linear k max growth of the variance are canceled analytically, not numerically.To make the cancellations explicit, we integrate by parts and find that

͑26͒
Substituting Eq. ͑26͒ into Eq.͑24͒, we obtain the final expression for the heat capacity

͑27͒
Equation ͑27͒ is the principal result of this paper.The convergence rate as a function of k max of the heat capacity obtained from Eq. ͑27͒ is the same as the convergence rate given in Eq. ͑24͒, but the variance of the estimator expressed in Eq. ͑27͒ is only weakly dependent on k max .To determine this convergence rate, we rewrite Eq. ͑16͒ in the form

͑28͒
Because each term on the right-hand side of Eq. ͑28͒ converges asymptotically as k max Ϫ2 , we can then write for the heat capacity as k max →ϱ.The asymptotic convergence characteristics expressed in Eq. ͑29͒ can also be verified directly using the methods described in Ref. 8. The derivation of the heat capacity estimator presented above needs modification for systems having a potential energy that is independent of the coordinate of the center of mass.An example of such a system is the 13-particle neon cluster examined in Sec.III.For such systems, just like the case for the virial energy estimator, 8 the free motion contribution of the center of mass is excluded by application of Eq. ͑27͒.For such systems the correction of 3k B /2 must be added explicitly.

III. APPLICATIONS A. One-dimensional systems
Simulations of single-basin one-dimensional systems can be performed in sufficient detail to provide useful results for initial investigations.In this subsection we investigate two model one-dimensional systems; the harmonic oscillator and a Lennard-Jones cage 11 often used to model the features of liquids.The thermodynamic properties of the harmonic oscillator can be evaluated analytically, and the analytic expressions provide a useful way of understanding the basic structure of the method.
For the oscillator we define the parameter ␥ϭប␤, where is the oscillator frequency.Using either Eq.͑24͒ or Eq.͑27͒ the specific heat can be shown to take the form

͑30͒
where 2␤͗V ¯kmax and the primes on the summations in Eqs.͑30͒-͑33͒ imply that only odd values of k are to be included.We have used this expression together with Metropolis MC 12 and Eq.͑27͒ to test the virial estimator.Figure 1 shows results of simulations of 10 6 passes at k max ϭ1, 2, 5, and 30.The number of quadrature points used in the u integrals has been set to max͕32,4ϫk max ͖ Gauss-Legendre points.The MC data are in excellent agreement with the analytical results ͓Eq.͑30͔͒, and the statistical errors in the simulations are smaller than the plotted points on the graph.From Fig. 1 and an analysis of Eq. ͑30͒, divergent behavior of the heat capacity at low temperatures ͑large ␥) is evident, and for a given temperature it is necessary to include sufficient Fourier coefficients to obtain meaningful results.A useful approach to the problem of determining the size of k max needed is the use of extrapolation methods, 8 which we now illustrate with the second one-dimensional model problem, the Lennard-Jones cage 11 V͑x ͒ϭ4 ͫͩ

͑34͒
This potential has been used to model simple quantum liquids in previous work. 11The potential curve has one single minimum between two strongly repulsive walls located at Ϯ␣.The parameters , , and ␣ have been chosen to model liquid argon with (ϭ119.4K, ϭ3.405 Å, ␣ϭ1.2 , and mϭ39.948amu͒.Metropolis MC simulations have been performed using k max ϭ1, 5, 7, 10, 15, 20, 30, and 50, with data accumulation over 10 6 passes after a warmup period of 10 5 passes.The number of quadrature points in the u integrations has been chosen to be max͕32,4ϫk max ͖ using Gauss-Legendre quadrature.As found in our study of the harmonic oscillator, the heat capacity is found to be divergent numerically for fixed k max in the limit of low temperatures.In the case of the Lennard-Jones cage potential, as well as the Lennard-Jones cluster discussed in the next section, the limiting form of this divergence takes the form of a power law, ␤ n .The exponent n of the power law is not universal and depends on the form of the potential energy function.For the three potentials examined in this work, numerical studies suggest that n is an integer that depends on the potential.
The heat capacity data as a function of k max Ϫ2 for several temperature values are shown in Fig. 2. As predicted by Eq. ͑29͒, the estimator is linear in k max Ϫ2 , and the linear behavior can be used to extrapolate the heat capacity at a particular temperature.We use this extrapolation method to determine the heat capacity of Ne 13 as discussed in the next section.

B. Ne 13
4][15][16] For example, much can be learned about nucleation rates from the study of the free energy of formation of the clusters. 17The heat capacity, the principal quantity of interest in this paper, can provide valuable information about cluster melting phenomena. 15In contrast to bulk systems, when a cluster ''melts'' there is a range of temperatures over which the cluster can be considered simultaneously to occupy liquid-like and solid-like structures. 14Associated with such cluster melting regions are maxima in the heat capacity as a function of temperature. 15Previous heat capacity estimators 9 have been so costly to evaluate, especially at low temperatures where many Fourier coefficients are needed, that prior quantum simulations have been limited to Lennard-Jones clusters of size 7 or less.We now illustrate the utility of the heat capacity estimator derived in this work with a simulation of the quantum properties of Ne 13 , and contrast the calculated heat capacity with that obtained from a classical simulation.Classical 13-particle Lennard-Jones clusters have well-defined heat capacity maxima, 15 and neon clusters are known to have significant quantum contributions. 9,18,19he potential energy used in this work is the pairwise additive Lennard-Jones potential modified with a constraining potential where ϭ35.6 K and ϭ2.749 Å , and r i j is the distance between particles i and j in the cluster.The constraining potential V c (x) is introduced because at finite temperatures the cluster has a finite vapor pressure, and the constraining potential insures the cluster contains 13 atoms over the course of the entire simulation.The form of the constraining potential used in this work is identical to that employed elsewhere, 20 where x i is the coordinate of particle i and x c.m. is the position of the center of mass.In some previous simulations of Lennard-Jones clusters of this size, 15 the constraining radius R c has been set to 4 .In the current studies, we have encountered difficulties insuring an ergodic simulation with the constraining radius set to 4 owing to transitions to gas-like structures at higher temperatures.These ergodicity problems have been solved by setting the constraining radius to R c ϭ2 .For the classical system, the heat capacity melting peak location and height have been found to be indistinguishable in simulations using these two values for the constraining radius.Differences in the heat capacity have been seen only at temperatures higher than the melting peak maximum owing to differences in the analog of boiling behavior.The simulations have been performed using the parallel version 21 of the J-walking Monte Carlo algorithm. 22The details of the parallel J-walking algorithm and the application of J-walking to quantum systems are presented elsewhere, 9,21 and we make no effort to review the approach here.In the parallel J-walking application to Ne 13 , the initial high temperature set of generating processors has been fixed at 13.1 K.The simulations have consisted of 10 6 warmup passes followed by 4ϫ10 6 passes with data accumulation.The maximum number of generating processors has been set to eight and the size of the configuration array has been set to 5000 ͑see Ref. 21͒.Generating processor temperatures have been included whenever the jump acceptance probabilities have dropped below 10%.The generating processors pass the configurations of both the coordinates and the Fourier coefficients to the running processors.For the Fourier representation of the paths, k max has been set to 1, 2, 4, 8, and 16.
A graph of C V as a function of T is found in Fig. 3 for each value of k max studied.From the highest calculated temperature to temperatures just below the melting peak, the behavior of the heat capacity is only weakly dependent on k max .An exception is the height of the heat capacity maximum, which shows a small dependence on the number of Fourier coefficients included.At temperatures below 8 K the dependence of the calculated heat capacity on k max is strong.To determine the heat capacity at these lower temperatures, it is convenient to extrapolate to infinite k max by using the asymptotic convergence rate of k max Ϫ2 .As an example, the heat capacity calculated as a function of k max Ϫ2 is presented in Fig. 4 for Tϭ6 K. Using such extrapolated data, the full quantum heat capacity curve is presented in Fig. 5 along with the classical heat capacity for the same system, included for comparison.The classical heat capacity data have been generated using the standard serial J-walking algorithm. 22The quantum heat capacity approaches 0 with decreasing temperatures as it must, while the classical heat capacity approaches the equipartition limit at low temperatures.As expected, quantum effects shift the melting peak to a lower temperature than the corresponding classical maximum.The calculation of the quantum heat capacity at temperatures lower than those reported in Fig. 5 requires increased computational effort, because the relative error for the diminishing heat capacity becomes large and increasing numbers of Fourier coefficients are required as T→0.

IV. SUMMARY
We have developed a new heat capacity estimator based on the virial energy estimator used in Fourier path integral Monte Carlo simulations.The estimator of the heat capacity is constructed so that the variance of C V is, at most, weakly dependent on the number of Fourier coefficients included.In calculations using partial averaging, the heat capacity calculated with the estimator converges asymptotically to the exact result as k max Ϫ2 .For finite k max the heat capacity calculated with the estimator is found numerically to be divergent as T→0, but extrapolation techniques based on the asymptotic convergence rate enable the determination of useful results over a wide temperature range.The method has been applied successfully to some model systems as well as to Lennard-Jones representations of Ne 13 .

FIG. 1 .
FIG.1.Specific heat virial estimator vs ␥ϭប␤ for several values of k max in the one-dimensional harmonic oscillator.The lines are the analytical expressions for the V-estimator at ͑from left to right͒ k max ϭ1, 2, and 5.The thick line at the right is the exact expression for C V at the infinite k max limit.The data points are taken from Monte Carlo simulations.The error bars are smaller than the resolution of the graph.

FIG. 2 .
FIG. 2. Specific heat virial estimator as a function of k maxϪ2 at fixed values of T for the one-dimensional Lennard-Jones cage potential.The circles are the results from simulations and the lines are the best linear fit for the data points.

FIG. 3 .
FIG. 3. Specific heat virial estimator as a function of T for several values of k max for the Lennard-Jones representation of Ne 13 .The errors bars represent one standard deviation.