Solution to the key problem of statistical physics -- calculations of partition function of many-body systems

TL;DR

This paper introduces a novel, faster method for calculating the partition function of many-body systems, enabling accurate predictions of thermal properties in complex molecules and condensed matter, with significant improvements over existing algorithms.

Contribution

The authors present a new computational approach that significantly accelerates partition function calculations, allowing for precise thermal property predictions in large and complex systems.

Findings

Method is at least four times faster than existing algorithms.

Achieves higher precision in thermodynamic property calculations.

Successfully reproduces the isochoric equation of state for solid argon.

Abstract

The key problem of statistical physics standing over one hundred years is how to exactly calculate the partition function (or free energy) of many-body interaction systems, which severely hinders application of the theory for realistic systems. Here we present a novel approach that works at least four orders faster than state-of-the-art algorithms to the problem and can be applied to predict thermal properties of large molecules or macroscopic condensed matters via \emph{ab initio} calculations.The method was demonstrated by C$_{60}$ molecules, solid and liquid copper (up to $\sim 600$GPa), solid argon, graphene and silicene on substrate, and the derived internal energy or pressure is in a good agreement with the results of vast molecular dynamics simulations in a temperature range up to $2500$K, achieving a precision at least one order higher than previous methods. And, for the first time, the realistic isochoric equation of state for solid argon was reproduced directly from the partition function.

Equations

$$a_i^{\prime }={\int }_0^{a_i}\exp [-U(0\dots q_i\dots 0)]d{q}_i,(i=1,2\dots N)$$

$$I_{ND}\simeq \prod _{i=1}^N{a}_i^{\prime }.$$

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

$$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})].$$

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

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

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

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

$$\mathcal{Q}=e^{-\beta U_0}\prod _{I=1}^M{\left[{\mathcal{V}}_I\right]}^{N_I},$$

Full text

I Introduction

By the end of $19$th century, the born of statistical physics brought a promising prospect that products of complex chemical reactionsChipot and Pohorille (2007) and all the thermodynamic properties of macroscopic systems, such as equation of states (EOSs)D Stacey (2005); Ross and Young (1993) and phase transitionsMonson and Kofke (2000), can be thoroughly predicted without empirical data by calculating the free energy (FE) or partition function (PF). Nevertheless, it was soon realized that such an implementation is impossible for large molecules and condensed mattersUshcats et al. (2016) because PF is in principle determined by all the possible microstates over entire phase space and a $3N$-fold integral has to be solved to obtain the PF for a system consisting of $N$ particles, which goes far beyond the capability of modern supercomputers if the standard algorithm of numerical integral is appliedAllen and Tildesley (1987). For instance, a rough calculation of PF for a C60 molecule, involving a $180$-fold integral, would cost at least $10^{100}$ years by using even the fastest high-performance computing facility with $\sim 10^{16}$ FP64 operations per second.

Alternative routes are resorted to sampling approaches, that is, either following the formalism of molecular dynamics (MD) simulation to trace the trajectories by integrating Newton’s equation of motionRapaort (2004), or, in a manner of Monte Carlo (MC) way to explore the microstates with substantial contributions to PF by stochastic walk in phase spaceLandau and Kurt (2005). Although time-average-based MD and ensemble-average-based MC algorithms have been developed over half a century and proved to be impressive to gain mechanical propertiesAllen and Frenkel (1989), when it comes to thermodynamic properties, the long-standing problem in face of both schemes remains that a balance has to be reluctantly sought for between limited-length sampling and the demanding requirement of ergodicity as the size of a system increasesBallard et al. (2015). Progress has been made in evaluations of the relative difference of PF (or free energy)Jarzynski (1997); Moustafa et al. (2015), and more attentions are being paid to the density of state to calculate absolute PFWang and Landau (2001); *DOS; *nxj; *multicanonical; *transitionmatrix. Nested sampling may be state-of-the-art techniqueSkilling (2004), which aims at uniformly sampling a series of fixed fractions partitioned by potential energies (PEs) in configurational spacePártay et al. (2010); *nestsample7 and has been applied in several systems described by empirical potentialsWilson et al. (2015); Pártay et al. (2014); Do et al. (2011); Coe et al. (2009); Baldock et al. (2016); Do and Wheatley (2013, 2016); Do et al. (2012); Burkoff et al. (2012); Baldock et al. (2017); Bolhuis and Csányi (2018). Despite of its improved $2$-$3$ orders of computational efficiency, NS can hardly work with ab initio calculations because of too much computational cost, and even if pairwise interaction potentials are employed, the affordable systems are limited to a scale of hundreds particles.

In this work, instead of tackling PF in the fashion of sampling, we established a direct integral approach (DIA) to solve the $3N$-fold integral on the basis of reinterpretation of original sense of integral, which works at least four orders faster than NS algorithm. Validations of the method were made in the systems including C60 clusters, condensed copper, solid argon, graphene and silicene on substrate, where the internal energies or EOSs obtained by DIA were compared to MD simulations, and for solid argon, the computed isochoric EOSs under high pressure zone were directly compared to experiments. Excellent agreements confirmed the accuracy of DIA. Ultrahigh efficiency of DIA paves a way to calculate the FE of large molecules and macroscopic condensed matters with ab initio computations.

II Direct Integral Approach to Partition Function

The original sense of one-fold ($1D$) integral $I_{1D}=\int_{0}^{a_{1}}f(x)dx$ is interpreted as the sum of infinite number of rectangles with area $A_{i}=f(x_{i})\Delta x$, and $\displaystyle I_{1D}=\lim_{\Delta x\to 0}\sum_{i}A_{i}$. Here, we interpret the integral from a different angle: The length of the $1D$ element $\Delta x$ at $x_{i}$ is modulated by $f(x_{i})$ to be a new length element $\Delta x^{\prime}_{i}=f(x_{i})\Delta x$ and $\displaystyle I_{1D}=\sum_{i}\Delta x^{\prime}_{i}$. In other words, the $1$D integral is mapped to a summation of length elements instead of area elements and equals to an effective length of $a_{1}$ (see left in Fig.1). Similarly, a two-fold integral $I_{2D}=\int_{0}^{a_{1}}\int_{0}^{a_{2}}dxdyf(x,y)$ equals to an effective area of $a_{1}\cdot a_{2}$ because the area element $ds=dxdy$ is enlarged (or shrunk) by $f(x,y)$ giving rise to an effective area element $ds^{\prime}=f(x,y)dxdy$ (see right in Fig.1). Followed by this notion, an $N$-fold integral $I_{ND}=\int_{0}^{a_{1}}\int_{0}^{a_{2}}\ldots\int_{0}^{a_{N}}dq_{1}dq_{2}\ldots dq_{N}f(q_{1},q_{2}\ldots q_{N})$ equals to an effective volume of $a_{1}\cdot a_{2}\ldots a_{N}$.

When the integrand $f(q_{1},q_{2}\ldots q_{N})$ is in a form of $\exp[-U(q_{1},q_{2}\ldots q_{N})]$ with $U(q_{1},q_{2}\ldots q_{N})$ being positive definite within the integral domain and having minimum at the origin ($U(0)=0$), the effective length of $a_{i}$ is defined as

[TABLE]

and the effective volume approximates to a product $\displaystyle\prod_{i=1}^{N}a^{\prime}_{i}$ (see proof in Supplementary Information), i.e.,

[TABLE]

Now consider the configurational integral (CI) in PF ($\mathcal{Z}$) for a continuum system consisting of $N$ particles at a given temperature $T$ (see details in Supplementary Information),

[TABLE]

where $\beta=1/k_{B}T$ with $k_{B}$ the Boltzmann constant, $q^{3N}=\{q_{1},q_{2}\ldots q_{3N}\}$ the Cartesian coordinates of particles and $U(q^{3N})$ the potential function. Although the integrand is of the same form as required by Eq.(2), it may not be positive definite or have no minimum at the origin ($q^{3N}=0$). Letting the set $Q^{3N}=\{Q_{1},Q_{2}\ldots Q_{3N}\}$ be the coordinates of particles in state of the lowest potential energy $U_{0}$, we may introduce a function

[TABLE]

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

[TABLE]

Clearly, $U^{\prime}(q^{\prime 3N})$ is positive definite within all the integral domain and has minimum at the origin ($U^{\prime}(0)=0$). According to Eq.(2), the integral in Eq.(5) equals to an effective $3N$-fold volume,

[TABLE]

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

[TABLE]

In this way, the $3N$-fold integral is turned into one-fold integrals.

For some 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.(6) turns into

[TABLE]

where $\mathcal{L}$ is the effective length determined by Eq.(7). Otherwise, it is needed to calculate the effective length, $\mathcal{L}_{x}$, $\mathcal{L}_{y}$, $\mathcal{L}_{z}$ by Eq.(7) of an arbitrary particle, and Eq.(6) turns into

[TABLE]

The procedure can be extended to systems composed of different particle species by calculating the effective length of each species respectively.

For inhomogeneous systems, such as defects or interfaces existed, particles may be grouped into $M$ sets numbered by $I$ with each containing $N_{I}$ equivalent particles, and CI becomes

[TABLE]

where $\mathcal{V}_{I}=\int e^{-\beta U^{\prime}(x)}dx\int e^{-\beta U^{\prime}(y)}dy\int e^{-\beta U^{\prime}(z)}dz$ denotes the effective volume of an arbitrary particle in the $I$th set with Cartesian coordinates $x$, $y$, $z$.

To implement DIA, the first step is to find the most stable structure (MSS) of the system for determining $U_{0}$, which can be accomplished in principle by several well-developed methods, such as global optimizationsOganov and Glass (2006); *glopt2; *glopt3; *glopt4 or dynamic dampingZhang and Buch (1990); Ye et al. (2009). Actually, the MSS for crystals can be immediately obtained by placing the particles right at the lattice sites. Then, we move a particle along its one degree of freedom $q^{\prime}_{i}$ to obtain $U^{\prime}(0\ldots q^{\prime}_{i}\ldots 0)$ while its other degrees of freedom $q^{\prime}_{j}$ and all the other particles are kept fixed. Clearly, this is an easy task for ab initio calculations and therefore, the PF of a $N$-particle system can be obtained even if we have no knowledge about the analytical expression of potential function $U(q^{3N})$, which is usually hard to be constructed precisely for realistic systems composed of more than one kind of particles.

For testing the DIA, we may perform ab initio calculations of $U^{\prime}(0\ldots q^{\prime}_{i}\ldots 0)$ on some realistic systems and compare the derived results with related experiments. However, results from first principle calculations may be strongly dependent on the specific algorithms, such as different types of exchange-correlation functions in density functional theory, and the experimental data are usually insufficient for extensive comparisons. In such cases, even if the results derived from the PF are in good agreement with the experimental data, it would be yet doubted of the accuracy of DIA. In order to have a stringent test, empirical potentials were used in the computations of PF, and the derived results were compared to the MD simulations using the same potentials to see if there exist some deficiency in DIA.

III Results and discussions

III.1 For isomers of C60 molecule

For a cluster consisting of $60$ carbon atoms, the most appealing impression may be the discovery of buckminsterfullerene (BF) with a football-cage structurekroto et al. (1985). Following investigations revealed that, besides the lowest potential-energy BF, there exists a bunch of C60 isomers, such as those introduced by the Stone-Wales (SW) rearrangementHeggie et al. (2016). For a long time, researchesZhao et al. (2003); Bettinger et al. (2003); Li and Ning (2004) have been attempting to answer the very questions whether the isomers are able to survive in realistic systems at finite temperatures and what the probability relative to BF is. The most reasonable argument should be the ratio of the PF of isomers ($\mathcal{Z}_{ISM}$) to that of BF ($\mathcal{Z}_{BF}$), but to our best knowledge, no such theoretical works have been tried out. Among all the isomers, the stack-$1$ SW (SW$1$) isomer (see bottom in Fig.2(a)) has the lowest PE, which is $\sim 0.6$eV higher than that of BF when using many-body Brenner potentialLi and Ning (2004). Here DIA was applied to calculate $\mathcal{P}=\mathcal{Z}_{SW1}/\mathcal{Z}_{BF}$ to determine the relative surviving probability of SW$1$ isomer in a temperature range from $100$K to $2500$K.

To implement DIA, we firstly determined the MSS of the BF and the SW$1$ isomer by the Polak-Ribiere conjugate gradient algorithm. For the MSS of BF as shown in the top panel of Fig.2(a), each atom is shared by two hexagons and one pentagon, indicating that all the atoms are geometrically equivalent and thus only one atom is needed for the calculation of $\mathcal{Z}_{BF}$. For an arbitrary atom, the felt PEs along its three Cartesian coordinates are obviously not the same, so $\mathcal{L}_{x}$, $\mathcal{L}_{y}$ and $\mathcal{L}_{z}$ has to be calculated respectively and Eq.(9) was applied. The selected atom was moved step by step with an interval of $0.0001$ Å along its Cartesian $X$-axis (or $Y$, $Z$-axis, 0.5 Å in positive and negative direction) to obtain the corresponding potential-energy curves, $U^{\prime}(X)$ (or $U^{\prime}(Y)$, $U^{\prime}(Z)$), during which, the other two coordinates and all the other atoms were kept fixed (see details in Supplementary Information). Fig.2(b1) shows the internal energy $E$ of BF derived from the PF, $E=-\partial\ln\mathcal{Z}/\partial\beta$, where the atom labeled No.1 in the top panel of Fig.2(b1) was selected to calculate $\mathcal{Z}_{BF}$. As a comparison, MD simulations were performed to calculate the internal energy and the Nose-Hoover algorithm Evans and Holian (1985) for canonical ensemble was employed with a time step $0.1$ fs. The systems were allowed to relax 20 ps initially and continued to run for another 50 ps, during which averages of $E$ were recorded in every 10 fs. All the computations concerning the C60 molecules were carried out in the Large-scale Atomic/Molecular Massively Parallel Simulator software package Plimpton (1995), and the many-body Brenner potentialBrenner et al. (2002) was selected to characterize the interactions between carbon atoms. As we can see, $E$ obtained from $\mathcal{Z_{BF}}$ are in an excellent agreement with those from MD simulations and the relative difference of $E$, RDE($=|(E_{PF}-E_{MD})/E_{MD}|$), is $<0.4\%$ up to $2500$K, which verifies the accuracy of the PF obtained by DIA. It should be noted that, based on our derivations, DIA is independent of how the directions of the Cartesian coordinates are chosen. To test this, we arbitrarily selected four more atoms, labeled No.$2$-$4$ in the top panel of Fig.2(a), to calculate the $\mathcal{Z}_{BF}$ respectively. Although the potential-energy curves for different atoms are not the same, i.e., $U^{\prime}_{1}(X)$ differs a lot from $U^{\prime}_{4}(X)$, the calculated internal energies based on the five atoms are the same (see detailed data listed in Supplementary Information).

For the MSS of SW$1$ isomer shown in the bottom panel of Fig.2(a), atoms were divided into three groups, where the first group is the atoms shared by two hexagons and one pentagon (colored in brown), the second one is the atoms shared by one hexagon and two pentagons (colored in magenta), and the last one is the atoms shared by three hexagons (colored in black). Accordingly, Eq.(10) was used to calculate the $\mathcal{Z}_{SW1}$, and, an arbitrary atom in each group, surrounded by dashed lines in bottom panel of Fig.2(a), was selected to obtain the corresponding potential-energy curves respectively. The selected atoms were moved $0.5$ Å in positive and negative direction, and $10^{4}$ PEs were recorded to obtain the corresponding potential-energy curves. Fig.2(b2) shows the comparison of the internal energy derived from $\mathcal{Z}_{SW1}$ and MD simulations, where the RDE is $<0.15\%$ within the whole temperature range. To test the universality of DIA, we additionally chose other different atoms in the three groups to calculate the $\mathcal{Z}_{SW1}$ respectively and the obtained results are the same as those shown in Fig.2(b2) (see detailed data listed in Supplementary Information).

Fig.2(b3) shows the obtained FEs of BF and SW$1$ isomer. Although the FE of BF remains smaller than that of SW$1$ isomer within the whole temperature range, the difference is quite small, which makes it difficult to conclude whether BF is the more stable one than the isomer or not. On the other hand, an unambiguous picture was presented by calculating the ratio of PF between these two clusters. Apparently, SW$1$ isomer has negligible probability to survive at low temperature zone compared to BF, while, for $T>1300$K, the surviving probability of SW$1$ prominently increases and ends up to $\sim 6\%$ at $2500$K, exhibiting that there is a considerable chance of the isomer to form at higher temperature zone. Such a finding qualitatively agrees with previous results by $40$ns-long MD simulationsLi and Ning (2004) but the exact value of $\mathcal{P}$ differs from each other, which calls for further clarifications by future experiments.

III.2 For solid and liquid Copper

Considering the huge computational cost of the MD simulations rather than DIA, the number of copper (Cu) atoms in our model is limited to $4000$, which were confined in a cubic box with periodic boundary condition (PBC) applied. The tight-binding (TB) potentialCleri and Rosato (1993) was employed to describe the interatomic interactions. For solid Cu, the MSS was found by arranging the atoms at the FCC lattice sites. In consideration of the Fm-3m symmetry of FCC lattice, the Cartesian $Z$-axis of atoms is set to [$001$] direction so that the potential $U^{\prime}$ felt by an atom moving along the $X$-axis (or $Y$-axis) is the same as the one along the $Z$-axis. To obtain $U^{\prime}(Z)$, the $Z$ coordinate of the geometry-center atom was changed step by step with its $X$ and $Y$ coordinates fixed to record the potential energies, during which all the other atoms stay fixed as well. For the liquid system, the atoms were heated up to $2.5\times 10^{4}$K in MD simulation to generate a uniform distribution and a damped trajectory methodZhang and Buch (1990) was used to determine the MSS and potential energy $U_{0}$ (see details in Supplementary Information). Different from the cases in crystal Cu, the potential $U^{\prime}$ felt by a liquid atom moving along the $X$-, $Y$- or $Z$-axis may not be the same, so $U^{\prime}(x)$, $U^{\prime}(y)$ and $U^{\prime}(z)$ were calculated respectively and Eq.(9) was applied to obtain the CI.

Common procedures for MD simulations of a canonical ensembleYe et al. (2009) was employed to produce the internal energy ($E$) and pressure ($P$) of the systems contacted with a thermal bath at given temperatures, and the Verlet algorithmVerlet (1967) was employed for integrating the equations of motion with time step $0.1$ fs and $0.01$ fs for solid and liquid respectively. We first fully relaxed the systems for $10^{5}$ steps, and ran another $10^{5}$ (or $10^{6}$) steps to record the values of $E$ and $P$ for solid (or liquid). Attentions have to be paid that the commonly used Virial equationsAllen and Tildesley (1987) are inaccurate to compute pressure in the case of many-body potential with PBC appliedLouwerse and Baerends (2006); Thompson et al. (2009), and we employed the method proposed by TsaiTsai (1979) that considers stress and momentum flux across an area for conducting the statistic of $P$.

As shown in Fig.3, for the solid systems with atomic volume ranging from $11.85$ to $6.07$ Å3/atom, the internal energy $E_{PF}$ and pressure $P_{PF}$ derived from the PF are in excellent agreement with those ($E_{MD}$ and $P_{MD}$) obtained by MD simulations at temperatures from $10$K to $2500$K. For $T\leq 1000$K, the relative difference of internal energy RDE$=|(E_{PF}-E_{MD})/E_{MD}|$) is less than $0.09\%$ ((Fig.3a, see data in Supplementary Information). As the temperature rises up to $2500$K, the difference gets a bit larger (the maximum is $0.69\%$), which may be attributed to the statistical fluctuations of MD simulations because the relative standard deviations of internal energy (RSDE), $\sim 0.2\%$, increase by about three times of those ($\sim 0.08\%$) below $1000$K. As shown in Fig.3b, the relative difference of pressure (RDP$=|(P_{PF}-P_{MD})/P_{MD}|$) is on the level of $\sim 2\%$ or less, which is larger than the RDE. The reason should be that the MD method for statistic of pressure needs more simulation time since the relative standard deviations of pressure (RSDP), $\sim 10\%$, are about two orders larger than the RSDE. In order to confirm this conjecture, we did similar computations except that the TB potential was replaced with a Lennard-Jones (L-J) potentialAgrawal et al. (2002), for which the Virial theorem can be applied in the MD simulations to calculate the pressure exactly. For the volumes of the system ranging from $11.85$ to $7.54$ Å3/atom at $300$K, the RSDP reduces down to $\sim 0.03\%$, and correspondingly, the RDP gets smaller to a level of $\sim 0.1\%$ (see data in Supplementary Information).

Fig.4 shows that the internal energy and pressure of liquid Cu obtained by DIA are in good agreement with the MD simulations. The RDE for all the systems is less than $2\%$, which is about ten times larger than that for solid Cu. This difference may stem from the fact that the $U^{\prime}(0\ldots q^{\prime}_{i}\ldots 0)$ felt by a liquid atom differs a little from that felt by other atoms, which is not the same as the situation in the solid systems. So, we can perform DIA for more liquid atoms to obtain more accurate results. As to the pressure, the RDP is less than $2\%$ for most systems except for the one with a density of $8.9$ g/cc, which displays larger RDP and can be understood because the RSDP is apparently larger than others (see Fig.4(b3)).

It should be noted that the precision of DIA is so high that it has reached limit of the MD simulations. For instance, the internal energy obtained by DIA of solid Cu with the atomic volume of $11.85$Å3/atom is $-13894.76$eV for $300$K and $-13383.00$eV for $800$K, which are almost the same as $-13894.96$eV and $-13383.11$eV obtained by the MD simulations with RSDE of $\sim 0.4$‰. The gradual increases of RDE and RDP should be attributed to the fluctuation rise of the MD simulations (see Fig.3(a3) and (b3)), for which smaller time step should be applied to integrate the equations of motion at high temperatures to ensure the precision.

As shown in Fig.5, the isothermal EOS of solid Cu derived from the PF using TB potentialCleri and Rosato (1993) exhibits the same trend as the experimental results though the pressures are about $10\%$ larger than the measured valuesKraus et al. (2016); Dewaele et al. (2004). This discrepancy should be attributed to the inaccuracy of the empirical potential because the pressures obtained by DIA and MD simulations coincide with each other quite well. By contrast, we also used the L-J potentialAgrawal et al. (2002) to calculate the EOS by DIA. Although the outcome coincides well with the corresponding MD simulations, it deviates much more from the experiments. Accordingly, it calls for a more accurate interatomic potential, which may resort to ab initio calculations in the future.

III.3 For Solid Argon

Considering that solid argon (Ar) has been extensively studied experimentally and is believed to be well characterized by L-J pair potentialAllen and Tildesley (1987), we calculated the PF of solid Ar by DIA to produce the EOS and compared the results with experiments. To our best knowledge, no accurate EOS has been put forward for solid Ar by directly solving the PF, though various EOSs towards Ar systems have been put forward either based on a fitting-parameter procedure (see detailed reviews in Refs.Tegeler et al. (1999); Young et al. (2016)), or in a MC sampling way for disordered liquid and gas statesDesgranges and Delhommelle (2012); Baldock et al. (2016).

In our work, the system consists of $32000$ argon atoms placed at FCC sites and the computational procedure of DIA was the same as described in the system of solid Cu according to Eq.(8) except for the interatomic potential being replaced with L-J potentialAllen and Tildesley (1987). As shown in Fig.6(a) (see data in Supplementary Information), the pressure ($P_{PF}$) derived from the PF coincides quite well with the experimental measurements ($P_{EXP}$)Lewis et al. (1974). For the atomic volume ranging from $35.54$ to $38.75$ Å3/atom, the relative difference (RDP$=|(P_{PF}-P_{EXP})/P_{EXP}|$) between the theoretical pressures and the experimental ones is within $5\%$ in most conditions as demonstrated in Fig.6(b). The maximum RDP of $6.63\%$ occurs in the system with an atomic volume of $V=35.54$ Å3/atom at $80$K, and may result from the fact, according to the original referenceLewis et al. (1974), that the uncertainty of the measured pressure ranges from $\pm 50$ bars above $80$K and becomes larger for lower temperatures. Considering the possible influence from the empirical L-J potentials, we may reach a conclusion that the isochoric EOS obtained by DIA is in a good agreement with the experiment.

III.4 Comparisons with Nested Sampling

It is very necessary and interesting to compare DIA with NS, the state-of-the-art approach to PF, in terms of computational efficiency and precision. On the efficiency, the computational cost of the two methods is determined by the number of times, $N_{c}$, to calculate the total potential energy. In previous work employing NS in solid Al (or NiTi) systems of $64$ atomsBaldock et al. (2016), $N_{c}$ is about $10^{9}$ (or $10^{10})$, while $N_{c}$ for DIA applied in the solid Cu of $4000$ atoms is $1.5\times 10^{4}$, showing that DIA is at least four orders faster. The real computer time of DIA to calculate the PF of $4000$ Cu atoms characterized by the many-body TB potentialCleri and Rosato (1993) is about $5$ minutes with full use of a desktop 8-core AMD Ryzen $1800$X CPU ($3.6$GHz per core), and, for the $32000$ Ar atoms characterized by the pairwise L-J potentialAllen and Tildesley (1987), is about $2$ minutes using one physical core of the CPU. For a solid argon system of $500$ atoms described by L-J potential at temperatures ranging from $80$K to $300$K, we ran a NS algorithm with $N_{c}\sim 10^{9}$ and DIA with $N_{c}\sim 10^{4}$ (real computer time is about $4$ hours for NS and about $5$ seconds for DIA respectively by using one physical core of the CPU) to calculate internal energy and pressure, which were compared to the MD simulations using the same potential function, demonstrating that the precision of DIA is about $10$ times higherGong et al. (2019). The accuracy of DIA has also been proved by calculating the internal energy of graphene or $\gamma$-graphyne materials on Cu substrate using Brener potential functionLiu et al. (2019a), and silicene on Ag substrate using Tersoff potential functionLiu et al. (2019b). Due to the ultrahigh efficiency, DIA has been successfully applied to predict the optimal conditions for silicene growth on Ag substrate with ab initio calculationsLiu et al. (2019b), which should be the first time to calculate the absolute free energy at finite temperatures up to thousands Kelvins with first-principle laws.

IV Conclusion

In summary, by our reinterpretation of integral, DIA to PF of large molecules and condensed systems was established. The accuracy of DIA was strictly validated by vast MD simulations for C60 clusters, condensed Cu, solid argon, graphene and silicene on substrate, and by experiments for solid Ar respectively. Compared to state-of-the-art method for PF, DIA works at least four orders faster with about one order more precise. The new approach will find its vast applications in investigating the products of complex chemical reactions, thermodynamic properties of large molecules and macroscopic systems, which highly relates to designing novel material, predicting various phase transitions and parameter-free EOS under extreme conditions.

V acknowledgement

TCW acknowledges the support from National Natural Science Foundation of China under Grant No.21727801.