First Principles Calculation of the Entropy of Liquid Aluminum.

The information required to specify a liquid structure equals, in suitable units, its thermodynamic entropy. Hence, an expansion of the entropy in terms of multi-particle correlation functions can be interpreted as a hierarchy of information measures. Utilizing first principles molecular dynamics simulations, we simulate the structure of liquid aluminum to obtain its density, pair and triplet correlation functions, allowing us to approximate the experimentally measured entropy and relate the excess entropy to the information content of the correlation functions. We discuss the accuracy and convergence of the method.

1. Introduction

Let be the probability of occurrence of a state i, in thermodynamic equilibrium. The Gibbs’ and Von Neumann’s formulas for the entropy [1,2], are mathematically equivalent to the information measure defined by Shannon [3]. Entropy is thus a statistical quantity that can be calculated without reference to the underlying energetics that created the probability distribution, as recognized by Jaynes [4]. Previously we applied this concept to calculate the entropy of liquid aluminum, copper and a liquid aluminum-copper alloy binary alloy [5], using densities and correlation functions obtained from first principles molecular dynamics simulations that are nominally exact within the approximations of electronic density functional theory. In this paper we discuss the convergence and principal error sources for the case of liquid aluminum. As shown in Figure 1, we are able to reproduce the experimentally known entropy [6,7] to an accuracy of about 1 J/K/mol, suggesting that this method could provide useful predictions in cases where the experimental entropy is not known.

Figure 1

Calculated entropies compared with experimental values [6,7]. is from Equation (12), is from Equation (11), is the pair-correlation correction from Equation (14), and is from Equation (27)). We expect the best liquid state result from . In the solid state, below melting at T = 933 K, is the vibrational entropy in the quasiharmonic approximation.

In a classical fluid [8], the atomic positions and momenta ( for N atoms in volume V) take a continuum of values so that the probability becomes a function, , and the entropy becomes in the canonical ensemble. In this expression, the factor of corrects for the redundancy of configurations of identical particles, and the factors of Planck’s constant h are derived from the quantum mechanical expression. For systems whose Hamiltonians separate into additive terms for the kinetic and configurational energies, factorizes into a product of independent Maxwell-Boltzmann distributions for individual atomic momenta, times the N-body positional correlation function .

Equation (2) can be reexpressed in terms of n-body distribution functions [8,9,10,11], with , as where the n-body terms are

The subscripts N indicate that the correlation functions are defined in the canonical ensemble, with a fixed number of atoms N, and they obey the constraints

Each term can be interpreted in terms of measures of information. Briefly, is the entropy of a single particle in volume , and hence in the absence of correlations. is the difference between the information content of the pair correlation function , and the uncorrelated entropy, which must be added to . Similarly, is the difference between the information contents of the three-body correlation and the two-body correlation , which must be added to . Notice that the information content of the n-body is also contained in the -body and higher-body correlations because of the identity that expresses as a marginal distribution of .

Mutual information measures how similar a joint probability distribution is to the product of its marginal distributions [12]. In the case of a liquid structure, we may compare the two-body joint probability density [13,14] with its single-body marginal, . The mutual information per atom tells us how much information gives us concerning the positions of atoms at a distance r from another atom. Mutual information is nonnegative definite. We recognize the term in Equation (6) as the negative of the mutual information content of , with the factor of correcting for double-counting of pairs of atoms.

2. General Theory

2.1. One-Body Term

The one-body term in Equation (5) can be evaluated explicitly, yielding where is the quantum De Broglie wavelength. Both terms in Equation (11) have simple information theoretic interpretations [15]. While an infinite amount of information is required to specify the exact position of even a single particle, in practice, due to quantum mechanical uncertainty we should only specify position with a resolution of . Consider a volume . In the absence of other information, the probability that a single particle is localized within a given volume is . Summing over the ()-many such volumes yields . Similarly, the in Equation (11) is simply the entropy of the Gaussian momentum distribution, Equation (3).

Notice that resembles the absolute entropy of the ideal gas,

The difference lies in the constant term in vs. in . We shall discover that the difference is accounted for in the many-body terms . Indeed, this is clear if we place N particles in the volume . The derivation of Equation (11) generalizes to , but this must be corrected [15] by the irrelevant information, , that identifies the individual particles in each separate volume . The leading term of the Stirling approximation converts into , while the second term adds 1 to yielding .

Either or can be taken as a starting point for an expansion of the entropy in multi-particle correlations. Prior workers [11,16,17,18,19] tend to favor , while we shall find it more natural to begin with .

2.2. Two- and Three-Body Terms

Translational symmetry allows us to replace the double integral over positions and in Equation (6) for with the volume V times a single integral over the relative separation . A similar transformation applies to the integral for . However, the canonical ensemble constraint Equation (8) leads to long-range (large r) contributions to the remaining integrations. Nettleton and Green [16] and Raveche [17,18] recast the distribution function expansion in the grand-canonical ensemble and obtained expressions that are better convergent. We follow Baranyai and Evans [11] and apply the identity to rewrite the two-body term as

The combined integrand of falls off rapidly, so that the sum of the two integrals converges rapidly as the range of integration extends to large r. Furthermore, the combined integral is ensemble invariant, which allowed us to substitute the grand canonical ensemble radial distribution function in place of the canonical . The same trick applies to the three-body term,

The contribution of in as given by Equation (14), together with an added from the three-body Equation (17) and higher terms, reconciles the one-body entropy with the ideal gas. For consistency with previous workers [11,16,17,18,19] who omit the from and the from , and to make connection with the ideal gas, we can add the entire series to and write which is equivalent to Equation (4).

In the grand-canonical ensemble, the term in Equation (14) arise from fluctuations in the number of atoms, N, and can be evaluated in terms of the isothermal compressibility as where is the dimensionless compressibility. Note that , and hence also , are positive definite. The remaining term is the entropy reduction due to the two-body correlation. As noted above, the mutual information content of the radial distribution function reduces the entropy by

The complete two-body term is now .

The three-body fluctuation term (see Equation (16)) also relates to isothermal compressibility [18], with

The final term in Equation (17) reduces to a difference of three- and two-body entropies, and its sign is not determined. Essentially, the term adds back the two-body mutual information and then subtracts the information contained in the three-body correlation . Note that necessarily contains all the information in because of the identity Equation (9).

The pattern illustrated in Equations (21) and (24) holds for the analogous higher-body correlations as well, because integrals of the correlation function can be written in terms of integrals and density derivatives of . One limit of special interest is the incompressible limit, where the initial terms of Equations (14) and (17) vanish and only the information-derived terms survive. This limit should apply to dense fluids at low temperatures. Another limit occurs at high temperature, where the density drops and the correlation functions approach 1. In this limit all integrals involving vanish so that and all the terms sum to .

Truncation of the series of terms is accurate if higher many-body correlation functions can be approximated by products of fewer-body correlations. That is, if the higher correlation functions contain no new information. For example, the Kirkwood superposition approximation causes to vanish.

3. Results

To provide the liquid state correlation functions needed for our study we perform ab-initio molecular dynamics (AIMD) simulations based on energies and forces calculated from first principles electronic density functional theory (DFT). We apply the plane-wave code VASP [20] in the generalized gradient approximation [21]. Simulations are performed at fixed volume for each temperature. In order to determine the proper volumes (i.e., liquid densities ) we performed simulations at several volumes to identify the volume at which the pressure (including the kinetic term) vanished. Most runs were performed using Normal precision FFT grids, however the smallest system (N = 100 atoms) was found to require accurate precision.

Figure 2 shows the result of convergence studies in both volume and plane-wave cutoff energy. Briefly, we found minimal dependence on the plane wave energy cutoff, but strong and non-monotone dependence on the number of atoms. We accept atoms as a suitable target for convergence of the volume and we use the same condition for collecting our correlation functions. Our calculated density at 973 K falls below the experimentally assessed value by about 1%, similar to the discrepancy for solid Al in the limit of low temperature. From the volume-dependence of pressure we obtain estimates of the dimensionless compressibility ranging from 0.008 at T = 973 K up to 0.015 at T = 2473 K.

Figure 2

Calculated aluminum density vs. number of atoms N at various temperatures. All results hold for the default energy cutoff of 240 eV except for red squares that hold for 320 eV.

Pair correlation functions are collected as histograms in Å bins, normalized to and subsequently smeared with a Gaussian of width Å. Triplet correlation functions utilize bin widths of Å, normalized to , and are not smeared. Our run durations for data collection were 10 ps. All structures were thoroughly equilibrated prior to data collection.

Figure 3 illustrates the pair correlation function at various temperatures. Note the oscillations that extend to large r; presumably these oscillations are responsible in part for the oscillations in as a function of N. Note also the decreasing amplitude of oscillation with increasing temperature. Figure 4 illustrates the three-body correlation function for the special case of equilateral triangles with . The inset displays the ratio (see Equation (25)). Notice that is nearly a step function, with small decaying oscillations that diminish with increasing temperature.

Figure 3

Pair correlation function at various temperatures.

Figure 4

Triplet correlation function at various temperatures. Inset: Kirkwood ratio Equation (25).

3.1. One-Body Term

The one-body term explicitly depends on density, and also depends implicitly on temperature through the De Broglie wavelength . Taking our calculated densities, and evaluating , , and , we note that and are greater than, but rather close to, the experimental liquid entropies [6,7], as shown in Figure 1. The differences drop as the temperature grows, as expected because nonideality of the liquid metal becomes less important at high temperature.

3.2. Two-Body Term

In Figure 5 we plot the terms and as defined by Equations (21) and (23), respectively, where we integrate from zero separation up to a cutoff of R. Owing to the increase of the volume differential , oscillations of are magnified at large R. The fluctuation term appears to converge towards a value close to 0, consistent with the low compressibility of the liquid metal, while the information term converges towards a negative value. Note that the oscillations are nearly opposites, so that their sum converges rapidly towards a negative value of .

Figure 5

Two-body terms and and their sum from simulated pair correlation function at T = 973 K.

Adding the entropy reduction to the single-particle entropy yields values that are close to experiment but lie slightly below, as is evident in Figure 1 (blue triangles). However, we know that liquid metals have an electronic entropy (see Section 3.3), , and when we include that term (Figure 1, orange crosses) the values lie within 1 J/K/mol of the experimental values. Had we chosen to add to instead of adding to the values would have been greater by J/K/mol, resulting in poorer agreement (Figure 1 magenta + signs). In Section 3.4 we explain why is a more suitable starting point for an expansion in multiparticle correlation functions than is.

3.3. Electronic Entropy

The electronic density of states , which comes as a byproduct of first principles calculations, determines the electronic entropy [22]. At low temperatures, all states below the Fermi energy are filled and all states above are empty. At finite temperature, single electron excitations vacate states below and occupy states above, resulting in the Fermi-Dirac occupation function

Fractional occupation probability creates an electronic contribution to the entropy,

We apply this equation to representative configurations drawn from our liquid metal simulations, with increased k-point density in order to converge the density of states.

At low temperatures, the electronic entropy approaches , which depends only on the density of states at the Fermi level. However, at the high temperatures of liquid metals the electronic entropy requires the full integral as given in Equation (27), rather than its low temperature approximation.

3.4. Three- and Higher-Body Terms

We saw in Figure 5 that the integral in Equation (21) converges slowly to the dimensionless compressibility which is a positive but very small value. Accordingly, the same must be true for the integral of the three-body fluctuation term, Equation (24), and all higher-body terms as well. Thus all fluctuation terms are essentially negligible contributions to the entropy at the temperatures considered here. This observation must break down at sufficiently high temperatures, because in the limit of very high temperature all correlation functions approach 1, so that all integrals vanish. As noted by Baranyai and Evans [11], , , and . This limit only holds at extreme high temperatures and low densities, however, the small shortfall in at T = 2473 K could reflect a need to include a small fluctuation contribution due to the many-body terms at high temperatures.

We still need to discuss the three-body information term, (Equation (17)). Previous studies have discussed this term for model Lennard-Jones and hard-sphere fluids [23,24]. This term vanishes within the Kirkwood superposition approximation, , and as seen in Figure 4 this approximation is quite accurate even at T = 973 K. Presumably the nearly free electron character of aluminum, which causes its interactions to be well described by a nearly hard-sphere pair potential [25], leads to the weak form of . The deviations of from 1 are oscillatory, both in radial dependence as seen in Figure 4, and in angle as shown in Figure 6. We lack sufficient resolution in to evaluate the complete integral, however, integrating over r at fixed angle the terms are of magnitude 0.1 or less, and they reverse sign as a function of angle, leading to further cancellation.

Figure 6

Three-body integrand for isosceles triangles of angle and sides at T = 973 K.

4. Discussion

We find that the entropy of liquid aluminum is described rather accurately using the first two terms in an expansion of the entropy in multiparticle correlations. We show in particular that it is advantageous to start the series with rather than , and in compensation to include the terms , , … within , , …, respectively, because each of these terms then becomes of the order of the small dimensionless compressibility . The remaining terms, , each have a simple information-theoretic interpretation, with being the information to specify individual particle positions with resolution , being the mutual information content of the pair correlation function, and the corresponding higher order terms reflecting the additional information contained in that is not already present in the lower order terms.

In terms of accuracy, obtaining an accurate density is important. The difference between densities predicted at different system sizes N can shift the value of by about 0.5 J/K/mol, with greater density reducing . More significant is the impact of density on , with the same difference in density increasing the mutual information by up to 6 J/K/mol. Both of these potential sources of error substantially exceed the truncation error due to neglect of multiparticle correlations, a finding that may hold generally for nearly-free-electron metals, while transition metals with angle-dependent forces may require additional terms.