Direct computation of the quantum partition function by path-integral nested sampling
In the present work we introduce a computational approach to the absolute rovibrational quantum partition function using the path-integral formalism of quantum mechanics in combination with the nested sampling technique. The numerical applicability of path-integral nested sampling is demonstrated for small molecules of spectroscopic interest. The computational cost of the method is determined by the evaluation time of a point on the potential-energy surface (PES). For efficient PES implementations, the path-integral nested-sampling method can be a viable alternative to the direct Boltzmann summation technique of variationally computed rovibrational energies, especially for medium-sized molecules and at elevated temperatures.
In order to build useful databases for the modelling of radiative and chemical processes which take place under exotic conditions in the laboratory, on Earth, in the interstellar space, or in the atmosphere of distant exoplanets, we need access to accurate quantum partition functions and thermodynamic data for polyatomic molecules over a wide temperature range Mizus et al. (2017); Sousa-Silva et al. (2014a); Wenger et al. (2008); Yurchenko et al. (2014). As Ref. Sousa-Silva et al. (2014a) points out “the [absolute] partition function is necessary to establish the correct temperature dependence of spectral lines and their intensity”.
The quantum partition function of small molecular systems is usually calculated as the Boltzmann sum of variationally computed rovibrational energies including appropriate factors, which account for both the degeneracies and the spin-statistical weights. Although, the direct variational solution of the rovibrational Schrödinger equation Lauvergnat and Nauts (2002); Yurchenko et al. (2007); Mátyus et al. (2009); Fábri et al. (2011) provides all quantum dynamical information about the molecule, it is computationally feasible only for the smallest systems. And even for these, converging the value of the partition function at elevated temperatures has been considered as a challenging task Mizus et al. (2017); Sousa-Silva et al. (2014a); Wenger et al. (2008); Yurchenko et al. (2014), the variational computation of an increasing number of rovibrational energies is tedious, and near or beyond the lowest dissociation threshold it is non-trivial Szidarovszky and Császár (2015).
In order to overcome this problem, approximations to this rigorous approach can be made, e.g., the quantum partition function is very often estimated within the harmonic-oscillator-rigid-rotor model. However, the relative error of the quantum partition function from this model is several percent Furtenbacher et al. (2016) in an ideal case, moreover, this model is often qualitatively incorrect, e.g., for molecules with several torsional degrees of freedom.
Instead of looking for approximate models, we can settle for certain integral properties without explicitly computing the fully detailed quantum dynamical information. In this spirit, thermodynamic quantities of complex systems are usually calculated using the path-integral formalism of quantum mechanics Feynman et al. (2010); Ceperley (1995). For the theoretical description and understanding of many chemically interesting processes, we need only relative quantities, e.g., free-energy differences, and there are efficient approaches for computing a variety of these. In the context of the present work, we mention Refs. Garberoglio et al. (2014, 2017) in which a highly efficient path-integral Monte Carlo approach was developed for the fully quantum mechanical computation of the second virial coefficient relying on the relative properties of interacting molecules with respect to non-interacting systems. We also mention a recently published path-integral Monte Carlo approach, with a technical set-up rather different from the one proposed by us, for the direct computation of absolute rovibrational quantum partition function values and its application for the methane molecule Mielke and Truhlar (2015, 2016).
In the present work we develop a path-integral-based method which provides the absolute quantum partition function value for small molecular systems with an “arbitrary” connectivity, and thus can serve as a practical alternative to the variational-summation technique. We will focus on a temperature range (above a particular value – typically above 100 K) in which the different spatial symmetry species have the same density of states (the equality holds to a very good approximation) Sousa-Silva et al. (2014a); Yurchenko et al. (2014). In this temperature range, the effect of the spin-statistical weights (SSW) can be accounted for by a simple multiplication of the partition function computed without SSWs (more precisely, with equal SSWs), which can be straightforwardly obtained with path integrals.
The common wisdom in the field of molecular simulations is that it is not practical to directly compute the absolute partition function (neither quantum nor classical) due to the exponentially narrow regions of the configuration space which contribute to the integral value. Although there are good reasons to follow this hint, we nevertheless attempt the direct computation of the absolute partition function using the nested sampling approach Skilling (2004, 2006). Nested sampling has already been successfully used for a series of different classical systems Pártay et al. (2010, 2014); Baldock et al. (2016, 2017) to compute absolute thermodynamic quantities, and without having to have any a priori knowledge of the particular shape and location of the important basins in the configuration space.
Ii Quantum Hamiltonian and partition function
The quantum nuclear Hamiltonian of a molecule of atoms (in atomic units) is
where are the masses associated to the nuclei, labels the laboratory-frame (LF) axes, and is the potential energy surface (PES), which depends on the relative positions of the nuclei for an isolated system. The energy-spectrum of is continuous due to the overall translation of the system. In the present work, we separate off the overall translation: subtract the kinetic energy of the center of mass () and introduce Jacobi Cartesian coordinates, (). We chose Jacobi coordinates because the corresponding Hamiltonian has a simple form, similar to Eq. (1), without any derivative cross terms:
only the physical masses are replaced with the reduced masses corresponding to the selected Jacobi coordinate set (see Appendix). Eigenvalues of are the rotation-vibration energies, , and the corresponding rotation-vibration partition function is
where the last equality holds if only bound states are populated at temperature ( is the Boltzmann constant). Obviously, the partition function corresponding to the full, laboratory-fixed Hamiltonian, , is related to the rovibrational partition function, Eq. (3), as
where is the translational partition function of the center of mass McQuarrie (2000):
and is the total mass. For the present work, we shall use the Jacobi Hamiltonian, Eq. (2), which gives direct access to the rotational-vibrational partition function.
Iii Path-integral formalism
Using the path-integral formalism of quantum mechanics and the Trotter factorization of an exponential function of two non-commutative (in this case the kinetic and the potential energy) operators, we write the quantum partition function in the form:
where . Using the properties of the position eigenfunctions, the integrals can be evaluated in the usual way Ceperley (1995); Tuckerman (2002) and if the kinetic energy operator has a simple form (no cross terms and constant coefficients), we arrive to the well-known ring-polymer expression of the quantum partition function
where the cyclic “boundary condition”, , is used to shorten the notation. In Eq. (8), is the ring-polymer’s harmonic angular frequency, and is the total number of physical degrees of freedom. For a molecule in the three-dimensional space . If the molecular Hamiltonian is written in LF coordinates (see eq. (1)), is the physical mass for the nucleus, while in the case of using Jacobi coordinates (see Eq. (2)), the is replaced by the reduced mass, (see Appendix), and is the number of the Jacobi vectors.
The quantum partition function in Eq. (8) can be considered as a classical configurational integral of an extended system, which includes copies (beads) of the original system in which the neighboring beads are connected with harmonic springs of angular frequency. In short, this is an dimensional hypothetical classical system at inverse temperature with the potential energy function
For the present discussion, it is an important difference from a truly classical system that the potential energy of the ring-polymerized system, , depends on the temperature through the angular frequency of the ring-polymer springs, .
Iv Nested-sampling integration
The nested sampling technique had been originally introduced by Skilling Skilling (2004, 2006), in the field of Bayesian probability and inference to sample high-dimensional spaces and computing Bayesian evidences,
where is a given likelihood and a given prior probability density function.
The method was later adapted to be used to sample the potential energy surface of atomistic systems Pártay et al. (2010); Baldock et al. (2016); Burkoff et al. (2016); Baldock et al. (2017); Pártay (2018), thus allowing the calculation of the partition function, which is normally inaccessible except for the simplest models.
The cornerstone of the method is the reduction of the many-dimensional integral in Eq. (10) to a one-dimensional one. For this purpose, we define as the amount of prior mass with likelihood greater than some threshold ,
which transforms into the integral (see also Figure 1)
At the beginning of the sampling we take samples randomly from the prior , this is going to be the initial sample set, then a series of iterations are performed. At each iteration step, , we choose the point with the lowest likelihood value, , remove and replace it by a new point under the constraints that it has to be uniformly drawn from the phase space, where the likelihood is larger than . Now, one can approximate the evidence in Eq. (11) via a quadrature sum over the removed points, . Since the points are uniformly distributed, the value of the prior mass in the iteration can be approximated as, Skilling (2006). This procedure of the method is also summarized in Figure 2. In terms of the atomistic PES, the likelihood is the Boltzmann factor, , thus, nested sampling is a “top-down” approach, starting from the high-energy region of the PES and going towards the global minimum through a series of nested energy levels. The iteration is stopped when the evidence is converged, i.e. the contribution by the latest likelihood value is smaller than a given pre-set tolerance Skilling (2006). This stopping criterion highly depends on both the actual problem and the likelihood function, thus can be different for different purposes, however, in our case, terminating the sampling when is sufficient.
iv.1 Computational strategy to generate a new sample point
The initial samples are generated by randomly choosing points within the entire PES, i.e., generating configurations where the (Jacobi) position vectors are randomly placed within the simulation cell. Due the fast shrinkage of the phase space volume at lower energy levels, generating new configurations randomly in subsequent iterations quickly becomes impossible. To overcome this, new points are generated by cloning a randomly chosen existing sample and performing a random walk until a point is sufficiently independent from the parent configuration. The random walk is a series of Monte Carlo steps, where the coordinates of the Jacobi position vectors are changed, with each step accepted if the Boltzmann-factor of the new configuration is lower than the current limit ( in the iteration). The step size is adjusted throughout the sampling to have an acceptance ratio around 20%.
Similarly to earlier work on classical systems Baldock et al. (2016, 2017), we parallelize the nested sampling algorithm by evolving as many sample configurations as the number of processors () for an steps, instead of decorrelating the single cloned configuration alone for steps. Thus, each configuration will be evolved through steps on average before being recorded as the one with the lowest likelihood.
iv.3 Classical vs. quantum mechanical applications
For classical interatomic potential models the energy function is independent of the temperature, thus the partition function can be calculated for arbitrary values from a single nested-sampling calculation Pártay et al. (2010, 2014); Baldock et al. (2016, 2017). On the contrary, the PES of the ring-polymerized system is temperature dependent, Eq. (9), so it is necessary to carry out independent nested sampling iterations for every temperature value. However, the derivatives of the quantum partition function at a certain temperature could be computed with no additional cost, which is an aspect to be explored in a future work.
V Numerical results
v.1 Test system: A nine-dimensional integral
In order to test our implementation and check the applicability of nested sampling for integrals similar to those of the (ro)vibrational partition function of polyatomic molecules, we first computed
as an explicit nine-dimensional integral. Numerical results are shown in Table 1, calculated both with a relatively small () and with a large () sample. We shall use similar sample sizes in exploratory and production runs, respectively, for molecular systems (see next section).
Number of sample points. Each initial sample point was generated from a uniform continuous probability density over with . New sample points were generated in an MCMC iteration including 800 steps and collective moves.
and : sample mean and sample variance of the mean calculated from 20 independent runs.
v.2 Rotating-vibrating molecules
In order to demonstrate the applicability of the proposed path-integral nested sampling method for real molecular systems, we have selected three polyatomic molecules of spectroscopic interest: the parent isotopologue of magnesium hydride, water, and ammonia. For these, both accurate (and cost-efficient) potential energy surfaces as well as benchmark-quality quantum partition function values are already available in the literature.
The rovibrational partition function values, calculated at several temperatures, are collected in Table 2, and the number of sample points, , the number of Markov chain Monte Carlo (MCMC) steps, , and the number of beads, , necessary for the current accuracy are also included. The number of nested-sampling iteration steps, , was between and (generally the number of necessary iteration steps, , increases linearly with the number of sample points, ).
The overall computational cost of path-integral nested sampling is determined by the cost of a single PES call multiplied with the number of PES calls, the latter being . Both the number of sample points, , and MCMC steps, , has an effect on the accuracy, and as a rule of thumb, determines the quality of the results (for sufficiently large and values) Baldock et al. (2017).
Looking at the results it is remarkable, that path-integral nested-sampling can reproduce the Boltzmann sum of the variational energies with a less than 1 % relative error. At a given temperature, the larger the number of beads is, the more accurate the path-integral Trotter factorization becomes. However, at higher temperatures fewer beads are sufficient to achieve high accuracy, and thus, in the region where the variational technique becomes very expensive (or even unfeasible), path-integral nested sampling remains a practical alternative.
For a start, we performed exploratory computations for the three systems studied in this work using a small sample, , and only MCMC steps. The sample mean of the path-integral nested-sampling results agreed with the variational reference values within a few percent, although the sample variance was large. These exploratory calculations took a couple of minutes on a laptop for NH, for which we made use of the very fast polynomial PES of Ref. Yurchenko et al. (2011). The computations for the smaller, triatomic, HO molecule took longer than for NH due to the longer evaluation time of the water monomer PES Medders et al. (2013); Babin et al. (2014). (At the same time, we could use this NH PES only up to ca. 3000 K. Beyond this temperature value a broader configuration space give significant contribution to the integral, Eq. (8), for which unphysical regions of the polynomial PES hindered the computations).
To obtain the more accurate path-integral nested sampling results of Table 2, we used two larger parameter sets: one with sample points and MCMC steps, and a ca. four times larger one with and . The sample mean values of the two computations are in excellent agreement with each other as well as with the variational reference values. We also note that by increasing and the sample variance of the mean is significantly reduced.
In so far as it can be ascertained, the sample size—necessary to converge the integral—only slightly increases with the number of beads, and it has a very weak dependence on the number of atoms, (the “real” physical degrees of freedom). Based on nested-sampling studies of classical systems Pártay et al. (2010, 2014); Baldock et al. (2016, 2017), we think that a larger value might be necessary for systems with a large number of important basins on the PES. The efficiency of generating a new live point (now, determined by the number of MCMC steps, ) can perhaps be further improved by using total-energy Hamiltonian Monte Carlo Baldock et al. (2017) instead of all-atom Monte Carlo moves used here, which is an aspect to be explored in future work.
Classical partition function values were obtained with PI-NS using a single bead, .
Variational reference values obtained from the Boltzmann sum of variational rovibrational energy eigenvalues including the appropriate degeneracy and spin-statistical weight factors.
PI-NS: path-integral nested sampling. is the number of beads. The initial sample was generated from a uniform continuous probability density function over . The sample mean, , and the sample variance of the mean, , were calculated from 20 independent runs.
Vi Summary and conclusions
In this work we present the implementation and the first numerical applications of a direct computational approach for the absolute rovibrational quantum partition function of small, polyatomic molecules of spectroscopic interest. The approach relies on the combination of the path-integral formalism of quantum mechanics and the nested sampling technique. Nested sampling is a multi-dimensional integration technique which can be efficiently used also for integrals which have regions which are exponentially localized but give significant contributions. The computational cost of the path-integral nested-sampling method is determined by the cost of a PES call, and scales linearly with the number of the path-integral beads, the necessary number of which is approximately proportional with the inverse temperature. Thereby, we expect path-integral nested sampling to be a feasible alternative to the Boltzmann summation technique of variationally computed energy levels for small, polyatomic systems at elevated temperatures, and for larger systems for which the direct variational computation of a large number of eigenvalues is not possible.
For the versatile applicability of the method, the development of new (or a boost of existing) accurate and cost-effective potential energy surfaces is necessary. Since we anticipate the high-temperature region particularly well suited for path-integral nested sampling, the (accurate and cost-effective) PES should describe a sufficiently broad region of the configuration space to be able to account for high-temperature (and high-energy) phenomena, including dissociation.
B. Sz. and E. M. gratefully acknowledge the financial support of a PROMYS Grant (no. IZ11Z0_166525) of the Swiss National Science Foundation. B. Sz. also thanks the European Social Fund. EFOP-3.6.1-16-2016-0023. During this work we used the computer cluster ATLASZ of ELTE and the NIIF Infrastructure in Debrecen. We thank Dr. Tamás Szidarovszky for sending to us his implementation of the MgH PES based on Shayesteh et al. (2007); Szidarovszky and Császár (2015). L. B. P. acknowledges support from the Royal Society through a Dorothy Hodgkin Research Fellowship. We also thank an Instant Access Grant of the ARCHER Supercomputing Center, which allowed us to test the scalability and applicability of the parallelized implementation for large-scale production runs.
- Mizus et al. (2017) I. I. Mizus, A. Alijah, N. F. Zobov, A. A. Kyuberis, S. N. Yurchenko, J. Tennyson, and O. L. Polyansky, Mon. Not. R. Astr. Soc. 468, 1717 (2017).
- Sousa-Silva et al. (2014a) C. Sousa-Silva, N. Hesketh, S. Yurchenko, C. Hill, and J. Tennyson, J. Quant. Spectr. & Rad. Transfer 142, 66 (2014a).
- Wenger et al. (2008) C. Wenger, J. P. Champion, and V. Boudon, J. Quant. Spectr. & Rad. Transfer 109, 2697 (2008).
- Yurchenko et al. (2014) S. N. Yurchenko, J. Tennyson, J. Bailey, M. D. J. Hollis, and G. Tinetti, Proc. Nat. Acad. Sci. USA 111, 9379 (2014).
- Lauvergnat and Nauts (2002) D. Lauvergnat and A. Nauts, J. Chem. Phys. 116, 8560 (2002).
- Yurchenko et al. (2007) S. N. Yurchenko, W. Thiel, and P. Jensen, J. Mol. Spectrosc. 245, 126 (2007).
- Mátyus et al. (2009) E. Mátyus, G. Czakó, and A. G. Császár, J. Chem. Phys. 130, 134112 (2009).
- Fábri et al. (2011) C. Fábri, E. Mátyus, and A. G. Császár, J. Chem. Phys. 134, 074105 (2011).
- Szidarovszky and Császár (2015) T. Szidarovszky and A. G. Császár, J. Chem. Phys. 142, 014103 (2015).
- Furtenbacher et al. (2016) T. Furtenbacher, T. Szidarovszky, J. Hrubý, A. A. Kyuberis, N. F. Zobov, O. L. Polyansky, J. Tennyson, and A. G. Császár, J. Phys. Chem. Ref. Data 45, 043104 (2016).
- Feynman et al. (2010) R. P. Feynman, A. R. Hibbs, and D. F. Styer, Quantum Mechanics and Path Integrals (Dover Publications, Mineola, New York, 2010).
- Ceperley (1995) D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
- Garberoglio et al. (2014) G. Garberoglio, P. Jankowski, K. Szalewicz, and A. H. Harvey, J. Chem. Phys. 141, 044119 (2014).
- Garberoglio et al. (2017) G. Garberoglio, P. Jankowski, K. Szalewicz, and A. H. Harvey, J. Chem. Phys. 146, 054304 (2017).
- Mielke and Truhlar (2015) S. L. Mielke and D. G. Truhlar, J. Chem. Phys. 142, 044105 (2015).
- Mielke and Truhlar (2016) S. L. Mielke and D. G. Truhlar, J. Chem. Phys. 144, 034110 (2016).
- Skilling (2004) J. Skilling, AIP Conf. Proc. p. 395 (2004).
- Skilling (2006) J. Skilling, J. Bayesian Anal. 1, 833 (2006).
- Pártay et al. (2010) L. B. Pártay, A. P. Bartók, and G. Csányi, J. Phys. Chem. B 114, 10502 (2010).
- Pártay et al. (2014) L. B. Pártay, A. P. Bartók, and G. Csányi, Phys. Rev. E 89, 022302 (2014).
- Baldock et al. (2016) R. J. N. Baldock, L. B. Pártay, A. P. Bartók, M. C. Payne, and G. Csányi, Phys. Rev. B 93, 174108 (2016).
- Baldock et al. (2017) R. J. N. Baldock, N. Bernstein, K. M. Salerno, L. B. Pártay, and G. Csányi, Phys. Rev. E 96, 043311 (2017).
- McQuarrie (2000) D. A. McQuarrie, Statistical Mechanics (University Science Books, Sausalito, California, 2000).
- Tuckerman (2002) M. E. Tuckerman, Path Integration via Molecular Dynamics (John von Neumann Institute for Computing, Jülich, 2002), vol. Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms, NIC Series, Vol. 10, pp. 269–298.
- Burkoff et al. (2016) N. S. Burkoff, R. J. N. Baldock, C. Várnai, D. L. Wild, and G. Csányi, Comp. Phys. Comm. 201, 8 (2016).
- Pártay (2018) L. B. Pártay, Comp. Mat. Sci. 149, 153 (2018).
- Yurchenko et al. (2011) S. N. Yurchenko, R. J. Barber, J. Tennyson, W. Thiel, and P. Jensen, J. Mol. Spectrosc. 268, 123 (2011).
- Medders et al. (2013) G. R. Medders, V. Babin, and F. Paesani, J. Chem. Theory Comput. 9, 1103 (2013).
- Babin et al. (2014) V. Babin, G. R. Medders, and F. Paesani, J. Chem. Theory Comput. 10, 1599 (2014).
- Shayesteh et al. (2007) A. Shayesteh, R. D. E. Henderson, R. J. L. Roy, and P. F. Bernath, J. Phys. Chem. A 111, 12495 (2007).
- Sousa-Silva et al. (2014b) C. Sousa-Silva, N. Hesketh, S. N. Yurchenko, C. Hill, and J. Tennyson, J. Quant. Spectr. & Rad. Transfer 142, 66 (2014b).
- Suzuki and Varga (1998) Y. Suzuki and K. Varga, Stochastic Variational Approach to Quantum-Mechanical Few-Body Problems (Springer-Verlag, Berlin, 1998).
Appendix: Definition of the Jacobi vectors, Jacobi Hamiltonian, and reduced masses
The Jacobi vectors, , and the position vector of the center of mass, are constructed in a linear transformation of the laboratory-frame (LF) Cartesian coordinates, as
where the transformation matrix (see for example p. 10 of Ref. Suzuki and Varga (1998)) is
with . Upon this linear transformation of the coordinates, the Jacobi determinant is 1, and the Hamiltonian, Eq. (1), transforms to
where the reduced masses are
and is the total mass. After subtracting the kinetic energy operator of the center of mass, the Jacobi Hamiltonian is obtained as