The EXP pair-potential system. III. Thermodynamic phase diagram

The exponential repulsive EXP pair potential is defined by

Despite its mathematical simplicity, the system of particles defined by this potential has not been investigated much on its own right. Usually, exponential functions appear in conjunction with other terms, either in pair potentials or in more sophisticated many-body potentials.^1–6 Morse in 1929 and Born and Meyer in 1932 used an exponential repulsive term in a pair potential in conjunction with a long-ranged attractive exponential or r⁻⁶ term.^1,2 In the 1960s, Kac and co-workers used a hard-sphere pair potential minus a long-ranged EXP term for rigorously deriving the van der Waals equation of state in one spatial dimension.⁷ Around that time, a number of papers considered the pure EXP pair-potential system, calculating virial coefficients^8–12 and the high-temperature equation of state.¹³ With reference to the EXP pair potential, Maimbourg and Kurchan have recently shown that for pair-potential systems with strong repulsions, the isomorph theory becomes exact in infinite dimensions.¹⁴ The EXP pair potential was also used recently by Kooij and Lerner in a study of unjamming in models with analytic pair potentials.¹⁵ 

We can think of three reasons for studying the EXP pair-potential system in detail. The first one is the above-mentioned fact that, despite the exponential function being one of the most fundamental functions of mathematics, few studies of the EXP system have been undertaken. A notable feature of the EXP pair potential is that it is finite at zero separation in contrast to popular pair potentials like the Lennard-Jones, inverse power-law, and Yukawa potentials. One might well argue that no pair potential of the real world can diverge, so in this sense, the EXP potential is more realistic than many extensively studied pair potentials. A second reason for studying the EXP pair-potential system is that it is useful for modeling certain systems. Thus, Monchick in 1959 argued that the EXP pair potential has “long been regarded as the true qualitative form of the repulsive intermolecular potential” at small distances or high temperatures,¹⁰ and Sherwood and Mason in 1965 emphasized that the EXP pair potential is more realistic than inverse-power-law pair potentials.¹¹ The EXP pair potential is a good approximation to the low-density limit of the well-known Yukawa (screened Coulomb) pair potential,¹⁶ which is used for modeling ions in solutions and ionic liquids,¹⁷ as well as dusty plasmas.¹⁸ The third and perhaps the most important reason the EXP system deserves a thorough investigation is that it may be regarded as “the mother of all pair potentials.” Thus as recently shown, the EXP pair potential provides an explanation of simple liquids’ quasiuniversality supplementing the standard hard-sphere explanation; this is because any pair-potential system is quasiuniversal if it to a good approximation can be written as a finite sum of EXP pair-potential terms with coefficients that are large compared to the temperature.^19,20

The present paper is the third in a series^21,22 investigating the physics of the EXP pair-potential system. Paper I studied structure and dynamics along the EXP system’s fluid-phase isotherms and isochores.²¹ As for any other purely repulsive pair potential, the EXP system has no liquid-gas transition. Nevertheless, in the thermodynamic phase diagram one can clearly identify a region close to the melting line of typical liquid structure and dynamics and a typical gas-phase region far from the melting line. These regions merge smoothly into one another. We refer to the fluid state points close to freezing as “liquid” states. Paper I gave an example of gas- and liquid-state quasiuniversality^19,20 by showing that the radial distribution function of the Lennard-Jones system is close to that of the EXP system at state points with the same reduced diffusion coefficient. Paper II studied the EXP system’s gas and liquid phase isomorphs,²² demonstrating the invariance of the reduced-unit structure and dynamics along the system’s isomorphs (lines of constant excess entropy) that is expected from the EXP system’s strong virial potential-energy correlations.^21–23 

This paper establishes the thermodynamic phase diagram of the EXP system. In the solid phase, we find two thermodynamically stable crystal phases, a face-centered cubic (fcc) phase at low densities and a body-centered cubic (bcc) phase at higher densities. The investigation first determines the stable crystal structures at zero temperature as a function of density and of pressure (Sec. II). Four different crystal structures are studied. As the density is varied, only the fcc and the bcc structures give minimum free energy, however. The finite-temperature fcc–bcc, fcc-liquid, and bcc-liquid phase boundaries are established in the temperature-pressure phase diagram by integration of the Clausius–Clapeyron equation, calculations that are validated by interface-pinning simulations at selected state points (Sec. III). Section IV summarizes the results of the investigation, while Sec. V gives our suggestions for how to determine, in general, a phase diagram numerically. Finally, Sec. VI gives a brief discussion.

All quantities are reported in the unit system defined by the EXP pair-potential parameters ε and σ, the “EXP unit system.”²¹ Thus, distance is reported in units of σ, particle number density in units of σ⁻³, temperature in units of ε/k_B, pressure in units of εσ⁻³, potential energy per particle, and chemical potential in units of ε, etc. These are the standard rationalized units used when reporting the results of numerical investigations of, e.g., the Lennard-Jones system. Note that EXP units differ from the “macroscopic” state-point-dependent units of isomorph theory used in this theory’s dimensionless so-called reduced quantities.^20,23,24

The study presented below focuses on the relatively low-density state points defined by ρ < 1 in which ρ ≡ N/V is the particle density. Likewise, when the pressure p is the relevant variable, we focus on low-pressure states (p < 1). Classical physics is assumed throughout, ignoring the fact that the real-world thermodynamics of low-temperature crystals is of course dominated by quantum effects. The zero-temperature calculations were performed by a custom-made code using 128-bit IEEE quad precision numbers and no pair-potential cutoff; in this way, one achieves energies that are accurate to 16 significant digits. The finite-temperature molecular dynamics (MD) simulations were carried out using the Roskilde University Molecular Dynamics (RUMD) code, a graphics-processing-unit based efficient MD code.²⁵ In these simulations, a shifted-force cutoff at r = 6 was employed.

At T = 0, we investigate the full EXP pair potential, i.e., no cutoff is introduced. Because different structures may have very similar energy, a high numerical accuracy is needed to determine the “ground-state” (T = 0) crystal structures. If $rij(0)$ for a given lattice is the equilibrium distance between particles i and j, the ground-state lattice energy U₀ is given by a sum over all particle pairs,

To calculate U₀, we included images of sufficiently many unit cells that the total energy is accurate to sixteen digits. For the highest density studied (ρ = 1), it was necessary to include up to 32 images in each direction. In order to minimize floating-point roundoff errors, the list of pair energies was sorted before the energies were added.

Figure 1(a) shows the different ground-state energies as a function of the density. The barely distinguishable blue, green, and red curves are the energies of the fcc, bcc, and hexagonal close packed (hcp) lattices, respectively. These curves agree well with the “first shell + mean-field” approximation detailed below (black curve). The cubic diamond (cd) structure has significantly higher energy than the fcc, bcc, and hcp structures (light green curve).

Figure 1(b) shows the energy difference between the bcc and fcc structures relative to the fcc energy, plotted as a function of the density. The lattice energies of the fcc and bcc structures are identical at the density given by

Below this density, the fcc structure has the lowest energy; above it, the bcc structure has the lowest energy. The fact that the hcp and cd structures at no density have lower energies than this follows from the plots of their energy relative to that of the fcc structure shown, respectively, in Figs. 1(c) and 1(d).

It is instructive to compare the lattice energies to the results of a simple mean-field approximation. For a uniform and continuous distribution of particles in space with density ρ, the energy per particle u is given by $u=(1/2)∫0∞vEXP(r)4πρr2dr$ in which the factor 1/2 compensates for double counting of the pair interactions. In the EXP unit system, this becomes

which is shown as the upper dashed curve in Fig. 1(a). The mean-field prediction works best at high densities where the pair potential becomes effectively long ranged, implying that the assumption of a uniform particle distribution in space is realistic.

In Fig. 1(a), the dotted curve is the “first-shell” prediction arrived at by including only the nearest-neighbor pairs of a fcc lattice: u = 6 exp(−r_nn), where $rnn=a/2$ is the nearest-neighbor distance and a is the lattice constant that in terms of the density is given by a³ = 4/ρ. At low densities, this prediction is barely distinguishable from the exact fcc–bcc-hcp energies, but at high densities, there are significant deviations because interactions here become effectively long ranged. The short-long dashed curve gives the prediction for a fcc lattice if the nearest-neighbor pairs are treated separately, while all other interactions are taken into account by integrating the mean-field approximation from a to ∞. This leads to

This “first-shell plus mean-field” approximation works so well that it cannot be distinguished from the exact fcc, bcc, and hcp energies in Fig. 1(a). The results for this approximation, which is not used further below, show clearly that the physics is dominated by the nearest-neighbor interactions at low densities, while high densities are described well by the opposite limit of a uniform particle distribution.

The hcp structure has energy close to, but higher than, that of the fcc structure [Fig. 1(c)]. The fcc and hcp structures are both close-packed with 12 nearest neighbors at the same distance for a given density, but the next-nearest neighbors of the hcp lattice (ABABAB plane packing) are closer than those of the fcc lattice (ABCABC packing), resulting in a slightly higher energy.

At zero temperature, the thermodynamic phase diagram is one-dimensional. It may be quantified by the density ρ or by the pressure p. The above analysis shows that the ground state is fcc at low densities and bcc at high densities, with a phase transition at densities close to that given by the condition of identical potential energy [Eq. (3)]. The corresponding one-dimensional pressure phase diagram is determined as follows. The stable phase is that with the lowest chemical potential. At zero temperature, the Gibbs free energy is given by G = U + W in which W ≡ pV is the virial. For a system described by the pair potential $v$(r), if r_ij is the distance between particles i and j, for any configuration the virial is given by W = −(1/3)∑_i<jr_ij$v′$(r_ij).²⁶ If u is the potential energy and $w$ is the virial per particle, the chemical potential μ (Gibbs free energy per particle) is given by

Figure 2 shows the relative difference in chemical potential between the bcc and fcc structures plotted as a function of the pressure. The stable phase is fcc for p < 2.651 02 × 10⁻³ and bcc for p > 2.651 02 × 10⁻³. The phase transition is of first order because at the coexistence pressure, the densities of the two phases differ; they are given by ρ_fcc = 1.746 902 × 10⁻² and ρ_bcc = 1.747 118 × 10⁻², respectively. In other words, if one plots the equilibrium density as a function of the pressure, there is a discontinuous density increase at p = 2.651 02 × 10⁻³. The phase transition is characterized by a quite narrow coexistence region, since the relative density change is only 1.2 × 10⁻⁴. The phase transition is thus only weakly first order.

The full pressure-temperature and density-temperature phase diagram of the EXP system is given in Fig. 9 in Sec. IV. To arrive at it, the present section establishes the phase boundaries in the pressure-temperature phase diagram. The fcc–bcc, fcc-liquid, and bcc-liquid phase boundaries are determined by numerical integration of the Clausius–Clapeyron equation. Interface-pinning simulations²⁷ are employed for locating the fcc–bcc-liquid triple point. The triple point is the starting point for the Clausius–Clapeyron integrations of the fcc-liquid and bcc-liquid coexistence lines, and the T = 0 coexistence point found in Sec. II is the starting point for the fcc–bcc coexistence-line integration.

While we at zero temperature investigated the full EXP pair potential, at finite temperatures, the potential was truncated at r = 6 in order to reduce the computational cost. At low densities, the error from the truncation is below the numerical accuracy, but at higher densities, the truncation gives rise to a detectable error. As an example, the pressure of the bcc structure at T = 0.0022 and ρ = 0.1 is 0.1191 for the full potential and 0.1186 for the truncated potential, compare Fig. 3. The truncation error is, however, below what is visible on the figures given below.

In order to determine the fcc–bcc coexistence line in the temperature-pressure phase diagram, we integrated numerically the Clausius–Clapeyron equation for the coexistence line slope^28–30 

Here, ΔV is the volume difference and ΔS is the entropy difference between the two phases. The latter quantity is determined from the fact that the Gibbs free energies of two phases are identical at the phase boundary. Because the kinetic energy is a function only of the temperature, if the phases are denoted by 1 and 2, the condition G₁ = G₂ implies U₁ + pV₁ − TS₁ = U₂ + pV₂ − TS₂, i.e., TΔS = ΔU + pΔV. Equation (7) thus becomes

The numerical integration of Eq. (8) was performed using the fourth-order Runge-Kutta algorithm,³¹ starting from zero temperature at which the coexistence pressure is p = 2.651 × 10⁻³ (Sec. II). To evaluate ΔU and ΔV on the right-hand side of Eq. (8), we conducted constant-pressure (NpT) simulations^32,33 of each of the two phases at the state point in question.

The resulting coexistence line is shown in Fig. 4. The choice of independent variable in the integration of the Clausius–Clapeyron equation is a matter of taste, except that close to the vertical tangent in Fig. 4 pressure must be used. The integration was carried out in two steps. First, we used temperature as the integration variable, increasing temperature in each step by 10%. In the second part of the integration, pressure was the integration variable, decreasing the pressure in each step by 10%.

The fcc–bcc-liquid triple point was determined using the interface-pinning method,²⁷ which briefly works as follows. For a liquid-solid phase transition, one adds to the system’s potential-energy function a harmonic bias potential that couples to a crystalline order parameter.²⁷ The bias potential we used is given by $Ubias=κ/2(|ρk|−a)2$ in which κ is a “spring constant,” a is an “anchor point,” $ρk=∑j=1N⁡exp(−ik⋅rj)/N$⁠, and k is the wave vector of a selected Bragg peak.²⁷ The field biases the system toward two-phase configurations. As shown in Ref. 27, the chemical-potential difference between the two phases is determined by the average force the field exerts on the system.

A snapshot from an interface-pinning simulation is shown in Fig. 5. Figure 6(a) shows the coexistence pressure determined by means of the interface-pinning method plotted vs the bcc-liquid chemical potential energy difference. Figure 6(b) plots the relative difference between the bcc-liquid and the fcc-liquid coexistence pressures. The triple point is identified from the condition that the two coexistence pressures are identical.

The predictions of the Clausius–Clapeyron integration may be checked by comparing to the results of interface-pinning simulations. This is done in Fig. 7(a) for the bcc-liquid phase boundary and in Fig. 7(b) for the fcc-liquid boundary. Note that the bcc phase is unstable above a certain temperature that defines a “re-entrant” point around the density 0.1. Above this density, the solid-liquid phase boundary has anomalous, i.e., negative slope, as found also, for instance, for sodium³⁷ and for the Gaussian-core model.³⁸ Tables I and II give interface-pinning simulation data for selected bcc-liquid and fcc-liquid coexistence state points.

To locate the re-entrant point more accurately, Fig. 8(a) shows the bcc-liquid chemical-potential difference Δμ as function of the pressure, calculated from interface-pinning simulations at a temperature close to the maximum identified in Fig. 7(a) (T = 2.3 × 10⁻³). From Δμ, the melting temperature T_m is estimated via T_m ≅ T + Δμ/Δs in which Δs is the liquid-solid entropy difference per particle. This equation is derived as follows. The expression for the entropy S = −(∂G/∂T)_p implies for the entropy per particle s = −(∂μ/∂T)_p, so for the liquid-solid difference, one has Δs = −(∂Δμ/∂T)_p. This implies T − T_m ≅ −Δμ/Δs because Δμ = 0 at coexistence. Figure 8(b) plots the melting temperature at each pressure calculated this way, allowing for a more accurate identification of the re-entrant point than merely estimating it from the maximum in Fig. 7(a).

The result of the above investigation is summarized in Fig. 9 in which (a) shows the pressure-temperature phase diagram and (b) shows the density-temperature phase diagram. The latter has regions of coexistence which are, however, so narrow that they appear merely as slightly varying line thicknesses. Figure 9(c) shows as a function of pressure the entropy of melting of the fcc-liquid and bcc-liquid transitions, as well as the phase-change entropy of the fcc–bcc transformation. Note that the latter is small compared to the melting entropy, which for a simple system liquid like the EXP system is typically around k_B per particle, i.e., unity in the EXP unit system. Finally, Fig. 9(d) shows the relative width of the three coexistence regions in the density-temperature phase diagram as a function of the pressure, confirming that the fcc–bcc transformation is only weakly first order. At high pressures, the fcc–bcc density difference becomes slightly negative; this happens at state points above the point where dp/dT changes sign, compare the Clausius–Clapeyron equation Eq. (7).

Establishing thermodynamic phase diagrams of model systems is an important objective of computational materials science and chemical physics. Phase boundaries are lines in the phase diagram where the free energies of two phases are identical. Unlike pressure or energy, it is not trivial to determine the free energy at a finite temperature state-point because entropy cannot be computed as an ensemble average. Several methods have been developed to overcome this difficulty, each with advantages and disadvantages. As an example, Widom’s insertion method³⁹ can be used to compute the free energy at a single state point, but it gives poor statistics for dense phases like the ordinary liquid and solid phases. Another possibility is to use thermodynamic integration from a state point at which the free energy is known, typically an almost ideal gas state point or a virtually harmonic crystal state point.⁴⁰ 

Integrating the Clausius–Clapeyron identity from a coexistence point³⁰ is computationally efficient since it relates to straightforward simulations of the bulk phases. A disadvantage of this method is that one coexistence point must first be determined accurately; moreover, systematic errors may accumulate during the integration. The interface-pinning method gives an accurate prediction of the coexistence point, but it is not computationally efficient because the equilibration time for interface fluctuations is often considerably larger than for the bulk phases.²⁷ 

The approach used in this paper is not specific to the EXP system and may be used for computing coexistence lines of other systems. The method we propose consists of the following three steps:

1. Determine a single coexistence point to a high precision using either a T = 0 calculation or the interface-pinning method.

2. Starting from this state point, compute the coexistence line by a fourth-order Runge-Kutta integration of the Clausius–Clapeyron equation.

3. Validate the accuracy of the computed melting line by interface-pinning simulations at selected state points.

The existence of a fcc–bcc transition is not unique to the EXP pair-potential system. For instance, the Yukawa pair potential exhibits the same transition.^41,42 This is consistent with the fact that the strong-screening limit of the Yukawa potential is dominated by the exponential “screening” term. Another interesting feature of the EXP system is the existence of a re-entrant point. This is also found in, for instance, the Gaussian-core model,³⁸ which like the EXP system has a finite value of the pair potential at zero separation.

To summarize, the thermodynamic phase diagram of the EXP pair-potential system has been determined below unit density (Fig. 9). At low densities and temperatures, the thermodynamically stable solid phase structure is fcc; at higher densities and temperatures, the solid structure is bcc. The fcc–bcc transition is weakly first order, i.e., with significantly lower volume and entropy changes than for typical melting transitions. At high densities, there is a re-entrant point above which the liquid-bcc phase boundary has a negative slope.

This work was supported by the VILLUM Foundation’s Matter (Grant No. 16515).