Comparison of two efficient methods for calculating partition functions

TL;DR

This paper compares a new direct integral method with nested sampling for calculating partition functions, demonstrating that the new approach is significantly faster and more precise in solid argon simulations.

Contribution

The paper introduces a direct integral approach that outperforms nested sampling in speed and accuracy for partition function calculations.

Findings

The new method is at least four orders faster than nested sampling.

The direct integral approach achieves about 10 times higher precision.

Results agree well with molecular dynamics and experimental data.

Abstract

In the long-time pursuit of the solution to calculate the partition function (or free energy) of condensed matter, Monte-Carlo-based nested sampling should be the state-of-the-art method, and very recently, we established a direct integral approach that works at least four orders faster. In present work, the above two methods were applied to solid argon at temperatures up to $300$K, and the derived internal energy and pressure were compared with the molecular dynamics simulation as well as experimental measurements, showing that the calculation precision of our approach is about 10 times higher than that of the nested sampling method.

Equations

$$\mathcal{Z}(N,V,T)=\frac{1}{N!{\Lambda }^{3N}}\int d{q}^{3N}\exp [-\beta U(q^{3N})],$$

$$\mathcal{Q}=\int d{q}^{3N}\exp [-\beta U(q^{3N})].$$

$$\mathcal{Q}=\frac{\int \exp (-\beta U)\Omega (U)dU}{\int \Omega (U)dU}$$

$$\begin{array}{rl}\mathcal{Q}& \approx \sum _m\frac{{\int }_{U_{m+1}}^{U_m}\Omega (U)dU}{\int \Omega (U)dU}\exp (-\beta ⟨U{⟩}_m)\\ & =\sum _m{\omega }_m\exp (-\beta ⟨U{⟩}_m),\end{array}$$

$$\mathcal{Q}=\sum _{m=1}^n({\alpha }^m-{\alpha }^{m+1})\exp [-\beta (U_m+U_{m+1})/2],$$

$$E=\frac{3}{2}N{k}_{B}T+\frac{\sum _{m=1}^n[(U_m+U_{m+1})/2]({\alpha }^m-{\alpha }^{m+1})\exp [-\beta (U_m+U_{m+1})/2]}{\sum _{m=1}^n({\alpha }^m-{\alpha }^{m+1})\exp [-\beta (U_m+U_{m+1})/2]}.$$

$$P\approx \frac{1}{\beta }\frac{1}{\mathcal{Q}(V)}\frac{\mathcal{Q}(V+\Delta V)-\mathcal{Q}(V)}{\Delta V}.$$

$$U^{\prime }(q^{\prime 3N})=U(q^{3N})-U_0,$$

$$\mathcal{Q}=e^{-\beta U_0}\int d{q}^{\prime 3N}\exp [-\beta U^{\prime }(q^{\prime 3N})].$$

$$Q=e^{-\beta U_0}\prod _{i=1}^{3N}{\mathcal{L}}_i,$$

$${\mathcal{L}}_i=\int e^{-\beta U^{\prime }(\mathrm{0...}{q}_i^{\prime }...0)}d{q}_i^{\prime }.$$

$$Q=e^{-\beta U_0}{\mathcal{L}}^{3N},$$

$$Q=e^{-\beta U_0}({\mathcal{L}}_x{\mathcal{L}}_y{\mathcal{L}}_z{)}^N,$$

$$E$$

$$P$$

$$\varphi (r_{ij})=4ϵ[(\frac{\sigma }{r_{ij}}{)}^{12}-(\frac{\sigma }{r_{ij}}{)}^6],$$

Full text

I Introduction

The born of statistical physics laid a solid foundation to predict thermodynamic properties of macroscopic condensed matters. Phase transitions Baldock et al. (2016); Hansen and Verlet (1969), protein folding Burkoff et al. (2012) and the optimal conditions for novel material growth could be predicted theoretically as long as the partition function (PF) or free energy can be evaluated Chipot and Pohorille (2007). Nevertheless, solutions to the PF has been a lon standing problem Ushcats et al. (2016) and attempts were reluctantly turned to the help of molecular simulations Singh et al. (2012). With precedent efforts made in calculating the relative free energy, e.g., Gibbs ensemble Monte Carlo (MC) Mastny and de Pablo (2005) sampling and thermodynamic integration Mitchell and McCammon (1991), more attentions have been paid to density of states (DOSs) for computing the absolute PF Bussi et al. (2006); Hansmann (1997); Wang and Landau (2001); Li et al. (2016). The Bayesian-statistics-based nested sampling (NS) may be the state-of-the-art one Skilling (2004); Skilling et al. (2006), which aims at uniformly sampling a series of fixed fractions partitioned by potential energies in the configurational space to calculate DOS and has been applied in several systems Baldock et al. (2016); Pá?rtay et al. (2010); Do et al. (2012); Do and Wheatley (2012); Brewer et al. (2011); Nielsen (2013); Wilson et al. (2015); Do et al. (2011); Bolhuis and Csányi (2018); Baldock et al. (2017); Martiniani et al. (2014).

Very recently, we put forward a direct integral approach (DIA) to calculate the PF of condensed matters Ning et al. and the high accuracy has been proved by molecular dynamics (MD) simulations of condensed copper and argon Ning et al. , graphene and $\gamma$-graphyne materials Liu et al. (a), and silicene Liu et al. (b). Based on our reinterpreting the original sense of integral, it was shown that DIA works at least four-order faster than NS Ning et al. . On the other hand, it has not yet been confirmed whether the DIA has improved the computational precision of precedent MC methods. In this work, we carried out detailed analysis of DIA and NS in terms of the computational precision, and performed MD simulations to test the precision of internal energy and equations of state derived from the PF. It should be pointed out that the tests with MD simulations, instead of experimental data, is the most rigourous because same interatomic potentials can be used in calculations of the PF and MD simulations, which have been proved to be capable of producing very accurate results for various systems Andersen (1980); Nosé (1984); Hoover (1986). If the results derived from PF are only compared with experimental measurements, just as in most previous works, it would yet be doubted that the method for calculating the PF is accurate or not even if the agreements are excellent since it would be very likely that a deficient algorithm combined with an inappropriate empirical potential accidentally gives rise to an outcome close to the experiment.

The paper is organized as follows. In Sec. II, NS and DIA were briefly formulated, and in Sec. III, we first discussed the relationship between efficiency and accuracy of NS, and then, performed MD simulations of solid argon to test the computational precision of DIA and NS, showing that DIA has a much higher precision than NS. In addition, we found that NS works badly for the highly-condensed systems while DIA has no such a problem. A comparison with experimental data of solid argon along the melting line was presented as well, which further validates that DIA is more accurate than NS.

II Methods

PF is defined as a summation over the probabilities of all the microstates, and for a canonical ensemble consisting of N particles confined in volume $V$ at temperature $T$, it reads

[TABLE]

where $\Lambda$ is the thermal wavelength, $\beta=1/k_{B}T$ with $k_{B}$ the Boltzmann constant, $q^{3N}=\{q_{1},q_{2}...,q_{3N}\}$ the Cartesian coordinates of particles and $U(q^{3N})$ the potential energy. The $3N$-dimensional integral on the right hand of Eq.(1) is solely related to the microscopic states in configurational space, the so-called configurational integral (CI),

[TABLE]

II.1 Nested sampling

In Eq.(2), microstates in configurational space are expressed in terms of coordinates of particles. From another point of view, we may also label the microstates by their corresponding potential energy and the integral can be rewritten in terms of the DOS of potential energy

[TABLE]

where $\Omega(U)$ is the DOS of potential energy.

The strategy of NS is to partition the configurational space into a series of energy-decrease subdivisions numbered by $m$. For the $m$th subspace with upper energy limit $U_{m}$, a fixed number of configurations ($L$) with each energy $\varepsilon_{i}<U_{m}$ are generated by MC method and ordered in a sequence as $\varepsilon_{1}<\varepsilon_{2}<\ldots<\varepsilon_{L}$. The lower energy boundary $U_{m+1}$, which is the upper one for the $(m+1)$th subspace, is set to be the energy of a fixed fraction $\alpha$ of current subspace, as $U_{m+1}=\varepsilon_{I}$ with $I=\alpha L$. By the NS algorithm, Eq.(3) can be simplified as Pá?rtay et al. (2010)

[TABLE]

where $\langle U\rangle_{m}$ is an averaged energy of the $m$th subspace and $\omega_{m}$ stands for the percentile of the $m$th phase space volume. It is obvious that $\omega_{m}=\alpha^{m}-\alpha^{m+1}$, and after $n$th iteration when the convergence condition is reached, CI is evaluated as

[TABLE]

where $\langle U\rangle_{m}$ is chosen to be the arithmetic average of the boundary energies of each sampled partition Do et al. (2011); Nielsen (2013). According to $E=-\frac{\partial\ln\mathcal{Z}}{\partial\beta}$, the internal energy of the N-particle system is calculated by

[TABLE]

For determining the pressure by $P=\frac{1}{\beta}\frac{\partial\ln\mathcal{Z}}{\partial V}$, another CI for the system with a volume of $V+\Delta V$ should be calculated and $P$ is obtained by

[TABLE]

II.2 Direct Integral Approach

Consider Eq.(2) and let the set $Q^{3N}=\{Q_{1},Q_{2}...Q_{3N}\}$ be the coordinates of particles in the state of the lowest potential energy $U_{0}$, we may introduce a function as

[TABLE]

where $q^{\prime}_{i}=q_{i}-Q_{i}$. By inserting Eq.(8) into Eq.(2), we obtain

[TABLE]

According to our very recent work Ning et al. , the integral can be solved as

[TABLE]

where $\mathcal{L}_{i}$ represents the effective length on the $i$th degree of freedom and is defined as

[TABLE]

For homogeneous systems with certain geometric symmetry, such as perfect one-component crystals, all the particles are equivalent and $U^{\prime}$ felt by one particle moving along $q^{\prime}_{x}$ may be the same as the one along $q^{\prime}_{y}$ (or $q^{\prime}_{z}$). In such a case, Eq.(10) turns into

[TABLE]

where $\mathcal{L}$ is determined by Eq.(11). Otherwise, it is needed to calculate the effective length, $\mathcal{L}_{x}$, $\mathcal{L}_{y}$, $\mathcal{L}_{z}$ respectively, and Eq.(10) turns into

[TABLE]

and, $E$ and $P$ are thus evaluated as

[TABLE]

III Comparisons and Discussions

The tested models are face-centered-cubic (FCC) solid Ar systems consisting of $500$ or $4000$ atoms confined in a cubic box with different sizes, and, NS and DIA were applied to calculate internal energy $E$ and pressure $P$ at different temperatures to be compared with MD simulations. The interatomic potential for solid Ar was characterized by the commonly used pairwise 12-6 Lennard-Jones (L-J) potential Allen and Tildesley (1989),

[TABLE]

where $r_{ij}$ is the distance between atoms $i$ and $j$, $\epsilon=117.05$ (K), $\sigma=3.4$ Å and the cutoff distance is $r_{cut}=12.0$ Å. The MD simulations with periodic boundary condition applied were performed by the Large-scale Atomic/Molecular Massively Parallel Simulator software package Plimpton (1995) with a time step of 0.1 fs. The Nose-Hoover constant-temperature algorithm Evans and Holian (1985) was used to produce a canonical ensemble at temperature $T$. The system was allowed to relax 20 ps at first and then continued to run for another 50 ps, during which averages of $E$ and $P$ were recorded in every 10 fs.

To implement NS, it should be at first to select appropriate values of $\alpha$ and $L$, which cooperatively balance the computational efficiency and precision of $\mathcal{Q}$ [Eq.(5)]. Apparently, the larger the values of $\alpha$ and $L$ are, the higher the calculation precision is, but the slower the computation speed is. Although the initial choice of $\alpha$ made by Pártay et al. Pá?rtay et al. (2010) is $L/(L+1)$, successive works Do et al. (2011, 2012); Do and Wheatley (2012); Wilson et al. (2015) have showed that a smaller value of $\alpha=1/2$ is sufficient enough to guarantee the calculation pricison and enables NS to be applicable to systems consisting of up to several hundred atoms, of which the computational cost is too expensive for NS with $\alpha=L/(L+1)$. Therefore, $\alpha=1/2$ was adopted in this work. Cares should be also paid to the value of $L$ because, besides the factors of efficiency and systematic errors mentioned above, fluctuations of the calculated results in NS simulations closely depend on $L$ for a fixed $\alpha$ Do et al. (2011). We performed NS with four different numbers of configuration ($L=45000,60000,75000,90000$) to calculate the pressures of the solid Ar system consisting of $500$ atoms with a density of $1.83$ g/cm3 at different temperatures, where the well-built cage model for solid systems Do and Wheatley (2012) was used. For each $L$, we ran the NS simulations $15$ times to produce the averaged value of pressure which was compared with MD simulations to see the relationship between the deviations and $L$.

As shown in Fig.1 sm , the pressures obtained by the NS are gradually approaching to those of the MD simulations as $L$ increases and the corresponding fluctuations of NS is relatively larger with the smallest $L$. On the other hand, it should be noted that the fluctuations does not monotonically decrease with the increase of $L$, which was also observed in previous works Do et al. (2012). The fluctuations for $L=60000$ and $90000$ are almost the same, which are about $30\%$ smaller than those with $L=75000$, though the pressures with $L=90000$ are slightly closer to the MD simulations. Considering that the computational time with $L=90000$ is twice as much as that with $L=60000$, we chose $L=60000$ in the following work and conducted the NS simulations at each $(N,V,T)$ conditions for 15 times to calculate the averaged values of internal energy by Eq.(6) and pressure by Eq.(7), where the volume difference $\Delta V$ was made by changing the length of the box by $1\%$ because our calculations showed that smaller volume difference would produce very unphysical results.

Relatively, systematic parameters are much fewer for implementation of DIA. For the solid Ar system, the atoms were placed right at the FCC sites to produce $U_{0}$, and $U^{\prime}(0...q^{\prime}_{i}...0)$ in Eq.(11) was obtained by moving the center atoms along its $Z$-axis ($[100]$ direction) by 2Å while the coordinates of its $X$-axis, $Y$-axis, and of all the other particles were kept fixed. $2\times 10^{4}$ potential energies were recorded to calculate the $\mathcal{L}$ by Eq.(11), and, the internal energy and pressure were subsequently calculated by Eqs.(14) and (15), where the volume difference was made by changing the length of the box by $10^{-3}\%$.

For the argon system of $500$ atoms with different densities ($1.83$, $2.13$, $2.43$ and $2.98$ g/cm3) at temperatures from $25$K to $300$K, $E$ and $P$ obtained by DIA and the NS are shown in Fig.2, where the corresponding quantities of $E_{MD}$ and $P_{MD}$ of MD simulations are also presented as comparisons sm . For the systems with a density of $1.83$ g/cm3 and $2.13$ g/cm3, the averaged relative difference of internal energy, RDE ($=|\frac{E-E_{MD}}{E_{MD}}|$), of DIA is less than $4.1\%$, which is about four times smaller than that, $16.6\%$, of NS NOT . As the density increases up to $2.43$ g/cm3 and $2.98$ g/cm3, the averaged $RDE$ of DIA decrease to $5.51\%$ and $0.48\%$ respectively, while the averaged $RDE$ of the NS climbs up to $36.44\%$ for the density of $2.43$ g/cm3 and the NS fails to work for the system with density of $2.98$ g/cm3. As to precision of the pressure, the averaged relative difference, RDP ($=|\frac{P-P_{MD}}{P_{MD}}|$), of DIA is $2.48\%$, $1.69\%$, $0.17\%$ and $0.63\%$ for the densities of $1.83$, $2.13$, $2.43$ and $2.98$ g/cm3 respectively, while the corresponding RDP of the NS is $10.22\%$, $9.18\%$, $4.54\%$ and $\infty$.

The above comparisons show that the calculation precision of DIA is much higher than that of the NS. Furthermore, DIA works better with increase of the density while the NS can hardly work when the density is higher than $2.98$ g/cm3. The difficulty should be attributed to numerical calculations of Eq.(5), where the factor ($\alpha^{m}-\alpha^{m+1}$) approaches to zero ($\alpha=1/2$) as $m$ approaches to larger number, meanwhile, the factor $e^{-\beta(U_{m}+U_{m+1})/2}$ increases quickly when $U_{m}<0$, which is the common case for the Ar systems with lower density and the product ($(\alpha^{m}-\alpha^{m+1})\cdot e^{-\beta(U_{m}+U_{m+1})/2}$) is not too large (or small) for the 16 bit number of computer to describe. However, when the density is large enough that the $U_{m}>0$, both ($\alpha^{m}-\alpha^{m+1}$) and ($e^{-\beta(U_{m}+U_{m+1})/2}$) approach to zero as $m$ getting larger, and the product ($(\alpha^{m}-\alpha^{m+1})\cdot e^{-\beta(U_{m}+U_{m+1})/2}$) gets to be so small (but not exactly ”0”) that the output of computer is exact ”0”, which makes the denominator in Eq.(6) be zero easily. For this reason, we failed to apply the NS to calculate $E$ and $P$ of the Ar system with a density of $2.98$ g/cm3. A larger value of $\alpha$ might be helpful while the computational efficiency would be slowed down. By contrast, DIA has no such a problem because the largest part of the potential energy, $U_{0}$ of the MSS, has been extracted in Eqs.(8) and (9), and the left part $U^{\prime}$ is small enough to guarantee the precision of the integral for high density systems.

The lower precision for the NS calculating the pressure can be understood as follows. The pressure is determined by Eq.(7), where the volume difference $\Delta V$ should be set as small as possible to achieve high precision. However, the integral $\mathcal{Q}$ of Eq.(5) is not very sensitive to the small changes of the volume $V$ because of the random characteristic of MC simulations, leading to large fluctuations of $\mathcal{Q}(V+\Delta V)-\mathcal{Q}(V)$ for each running of the MC simulation. Our calculations showed that the large fluctuations would produce unphysical pressures when the $\Delta V$ is smaller than $10^{-4}\%V$, which corresponds to the length of the cubic box changed by $1\%$ adopted in our calculations. In DIA for calculating the pressure [Eq.(15)], the involved quantities $U_{0}$ and $\mathcal{L}$ determined by Eq.(11) are all sensitive to volume of the system, so the volume difference in Eq.(15) can be set much smaller. We tried several values of the box length difference in the range of $10^{-1}\%-10^{-6}\%$ and confirmed that the obtained pressures converges at the volume difference of $10^{-13}\%$ ($10^{-3}\%$ box length difference).

The computational efficiency of the NS and DIA depends on the number of the total potential calculation. For the NS running, the MC algorithm has to work $6\times 10^{3}-8\times 10^{3}$ times each producing $60000$ configurations to reach the convergence, so $3.6\times 10^{8}-4.8\times 10^{8}$ times of potential energy calculations must be performed to produce the $U_{m}$ in Eq.(5). Because of the fluctuations, the NS was run $15$ times for a given system to produce the averaged results, thus the number of the total potential calculations is larger than $5.4\times 10^{9}$, which is about five orders of magnitude larger than the one, $2\times 10^{4}$, for running DIA in the same system.

Because of the ultra-high efficiency, DIA was applied to calculate the internal energy and pressure of solid argon composed of $4000$ atoms, on which the NS costs too much computer hours and we have to give up the calculations, and we performed MD simulations to give comparisons. As shown in Fig.3, both $E$ and $P$ obtained by DIA coincide well with MD simulations where both $RDE$ and $RDP$ of DIA are almost the same as those calculated in the $500$-atom systemsm .

Finally, a comparison was made of DIA and the NS with experimental data of solid Ar along melting line Crawford et al. (1976). Considering the lower efficiency of the NS, the simulated system for both DIA and NS consists of $500$ atoms and the computational procedures are the same as described above. As shown in Fig.4, the internal energy and pressure obtained by DIA are significantly better than those of the NS. The averaged relative deviation of internal energy and pressure to the experimental data is $5.34\%$ and $5.5\%$ for DIA, which are about 6 times smaller than the ones, $39.12\%$ and $28.72\%$, for the NS.

IV Conclusion

In summary, by comparisons with MD simulations as well as experimental data, we confirmed that the accuracy of DIA outperforms the NS. The precision of DIA is about four times higher than that of NS for low-density systems and about one order higher in high-density situations. We also analyzed the intrinsic deficiency of NS in calculations of systems under highly condensed situations. Since the efficiency of DIA is at least five orders faster than that of the NS at the same time, DIA paves a better way to investigate thermodynamic properties of condensed matters, especially the ones with high density under extreme conditions.

V Acknowledgement

TCW acknowledges the support by Nation Natural Science Foundation under Grant No.21727801.