This paper studies numerically the solid phase of a system of particles interacting by the exponentially repulsive pair potential, which is a face-centered cubic (fcc) crystal at low densities and a body-centered cubic (bcc) crystal at higher densities [U. R. Pedersen et al., J. Chem. Phys. 150, 174501 (2019)]. Structure is studied via the pair-distribution function and dynamics via the velocity autocorrelation function and the phonon density of states. These quantities are evaluated along isotherms, isochores, and three isomorphs in both crystal phases. Isomorphs are traced out by integrating the density-temperature relation characterizing configurational adiabats, starting from state points in the middle of the fcc-bcc coexistence region. Good isomorph invariance of structure and dynamics is seen in both crystal phases, which is notable in view of the large density variations studied. This is consistent with the fact that the virial potential-energy correlation coefficient is close to unity in the entire fcc phase and in most of the bcc phase (basically below the re-entrant density). Our findings confirm that the isomorph theory, developed and primarily studied for liquids, applies equally well for solids.
I. INTRODUCTION
The EXP pair potential is the purely repulsive exponentially decaying function,
The system of particles interacting via this pair potential has been studied much less than, e.g., the Lennard-Jones or inverse power-law (IPL) pair-potential systems. Nevertheless, papers reporting studies of the EXP system or, more generally, systems that involve an EXP term in the pair potential have appeared regularly ever since 1929.1–14 As argued in Ref. 15, the EXP system deserves a closer study for two reasons. First, real-world systems, e.g., the low-density limit of the Yukawa (screened Coulomb) potential system, are well described by the EXP pair potential, which also plays an important role in most potentials describing metals. Second, the EXP pair potential may be regarded as the “mother of all pair potentials” in the sense that it explains the quasiuniversality of the structure and dynamics of simple liquids that applies for a large class of pair potentials.16,17 Thus, any pair potential, which can be written as a finite sum of EXP terms with coefficients much larger than kBT, defines a system in the same quasiuniversality class as the hard-sphere and Lennard-Jones systems.16,17 The EXP system may, consequently, replace the discontinuous hard-sphere system as the generic system in liquid-state theory.
The first paper in this series studied the structure and dynamics along the EXP system’s fluid-phase isotherms and isochores.18 Paper I18 also provided an example of quasiuniversality by showing that the radial distribution function (RDF) of the Lennard-Jones system is close to that of the EXP system at state points where the two systems have the same reduced diffusion coefficient. Paper II studied the EXP system’s fluid phase isomorphs,19 demonstrating invariance of the reduced-unit structure and dynamics along isomorphs.19 Note that while there, as for any purely repulsive system, is just a single fluid phase, the EXP system has nevertheless typical gas and typical liquid state points. Moreover, in contrast to the Lennard-Jones system, the EXP system has strong virial potential-energy correlations at all low-temperature fluid-phase state points, including the system’s typical gas-phase state points. The region delimiting the liquid and gas phases defines the so-called Frenkel line, the different definitions of which20–22 are isomorph invariant.
Paper III15 established the thermodynamic phase diagram of the EXP system, which involves two crystalline phases: a face-centered cubic (fcc) phase at low densities and a body-centered (bcc) at higher densities. There is a re-entrant liquid phase at the highest densities studied, which are of order unity in the “EXP unit system” defined by σ and ε in Eq. (1). This means that above a certain temperature there is no stable solid phase, a property the EXP system has in common with other pair-potential systems such as the Gaussian-core model,23 which do not have a diverging energy at zero particle separation. Figure 1, which is reproduced from Paper III,15 shows the thermodynamic phase diagram of the EXP system in (a) the pressure-temperature representation and (b) the density-temperature representation. The phase transformations between the three phases are all of first order, i.e., associated with a density change and a latent heat. The coexistence regions in the density-temperature plot are barely visible in Fig. 1(b), however, except as a slight thickening of the phase-boundary lines. In particular, the fcc-bcc transition is only weakly first order; the relative density change is merely 0.3% at the liquid-fcc-bcc triple point and even lower (0.01%) at zero temperature.15
Thermodynamic phase diagram of the EXP pair-potential system (in the EXP unit system). [Reproduced with permission from Pedersen et al., J. Chem. Phys. 150, 174501 (2019). Copyright 2019 AIP Publishing LLC.] (a) shows the phase diagram in the pressure–temperature plane, (b) shows the phase diagram in the density-temperature plane. At low pressures/densities, the solid phase is a face-centered cubic (fcc) crystal, while at higher pressures/densities, the solid phase is body-centered cubic (bcc). The black dot marks the triple point at which the three phases meet. The star marks the re-entrant state point above which no solid phase exists. The discrete gray dashed lines in (b) are the three isomorphs studied in Sec. V.
Thermodynamic phase diagram of the EXP pair-potential system (in the EXP unit system). [Reproduced with permission from Pedersen et al., J. Chem. Phys. 150, 174501 (2019). Copyright 2019 AIP Publishing LLC.] (a) shows the phase diagram in the pressure–temperature plane, (b) shows the phase diagram in the density-temperature plane. At low pressures/densities, the solid phase is a face-centered cubic (fcc) crystal, while at higher pressures/densities, the solid phase is body-centered cubic (bcc). The black dot marks the triple point at which the three phases meet. The star marks the re-entrant state point above which no solid phase exists. The discrete gray dashed lines in (b) are the three isomorphs studied in Sec. V.
As an example of a quantity monitored in the two crystalline phases, Fig. 2 shows the excess isochoric specific heat per particle in units of kB, denoted by , in a diagram in which the colors reflect different values. The prevalent color is green in both the fcc and bcc phases, showing that the excess specific heat is close to the 1.5kB per particle expected for a perfect harmonic crystal according to the Dulong–Petit rule.24
Excess isochoric specific heat per particle in the solid phase in units of kB, denoted by . Most state points are green, showing that the specific heat is close to that of a perfect harmonic crystal (3kB per particle, implying ). The arrow marks the zero-temperature coexistence density found in Paper III: ρ = 1.747 64 · 10−2 in the EXP unit system. The line separating fcc and bcc state points found from simulations is full red; its extrapolation to the zero-temperature coexistence density15 is dashed.
Excess isochoric specific heat per particle in the solid phase in units of kB, denoted by . Most state points are green, showing that the specific heat is close to that of a perfect harmonic crystal (3kB per particle, implying ). The arrow marks the zero-temperature coexistence density found in Paper III: ρ = 1.747 64 · 10−2 in the EXP unit system. The line separating fcc and bcc state points found from simulations is full red; its extrapolation to the zero-temperature coexistence density15 is dashed.
II. METHODS
All simulations were carried out using the double-precision version of the Roskilde University Molecular Dynamics (RUMD) code, which is optimized for graphics-processing-unit (GPU) computing.25 The Newtonian molecular dynamics (MD) simulations were standard NVT simulations using the leap-frog integrator and the Nosé–Hoover thermostat. Brownian dynamics was employed to check the Nosé–Hoover NVT MD results. All findings relating to structure and dynamics reported below are based on Newtonian MD simulations at state points for which Brownian dynamics gives the same structure and specific heat. This is the case except at very low temperatures. As in Paper III,15 a cut-off at 6σ was used for all MD-simulated state points; a larger cut-off was occasionally used for the Brownian dynamics simulations (see below).
Brownian dynamics simulations were carried out using an integrator implemented in RUMD as part of this work. The integrator is given26 by (in which n is the discrete-time index, μ is the mobility, and R is the vector of all particle coordinates)
Here, Δt is the time step and N(0, 1) for each time step is a vector of random numbers drawn from a zero-mean unit-variance Gaussian distribution. Note that changing the friction coefficient affects both the magnitude of the displacements and the relative importance of the random forces compared to the forces from the gradient. The challenge is to find a friction coefficient for which the exploration of the phase space is reasonably fast.
Figure 3(a) shows the equilibrium potential-energy distributions obtained by three different simulation methods: NVT MD (black points), a standard small-step Monte Carlo method implemented for an extra consistency check (red points), and Brownian dynamics (blue points). The latter two methods result in the same Gaussian distribution, while this is not the case for the MD simulation. The reason is that MD is highly inefficient at equilibrating an almost harmonic crystal, making it virtually impossible to use MD at very low temperatures. Depending on the quantity of interest, we therefore used Brownian dynamics at low-temperature state points. Figures 3(c) and 3(d) show which state points were simulated by Brownian dynamics respectively MD, while Fig. 3(b) reports the friction coefficients used.
(a) Potential-energy distributions for three different simulation methods at a low-temperature state point. The distribution of potential energies using Brownian dynamics is given by blue symbols and that of an all-particle-move Monte Carlo algorithm by the red symbols. These two methods result in the same potential energy distribution, which as expected is Gaussian. In contrast, a bi-modal distribution is observed for Nosé–Hoover (NVT) MD simulations (black symbols). This is caused by the near-harmonic state of the system, resulting in extremely slow equilibration. NVT MD simulations are of little use at such low temperatures. (b) The friction coefficients λ ≡ 1/(mμ) used in the Brownian simulations, where m is the particle mass. The symbols mark simulated state points. (c) The Brownian dynamics cutoffs used. In order to arrive at cutoff-independent results, it was necessary to use a cutoff larger than 6σ at the highest densities and at certain low-temperature state points. (d) State points simulated by standard Newtonian Nosé–Hoover NVT MD dynamics. All results for structure and dynamics reported below were obtained by this dynamics.
(a) Potential-energy distributions for three different simulation methods at a low-temperature state point. The distribution of potential energies using Brownian dynamics is given by blue symbols and that of an all-particle-move Monte Carlo algorithm by the red symbols. These two methods result in the same potential energy distribution, which as expected is Gaussian. In contrast, a bi-modal distribution is observed for Nosé–Hoover (NVT) MD simulations (black symbols). This is caused by the near-harmonic state of the system, resulting in extremely slow equilibration. NVT MD simulations are of little use at such low temperatures. (b) The friction coefficients λ ≡ 1/(mμ) used in the Brownian simulations, where m is the particle mass. The symbols mark simulated state points. (c) The Brownian dynamics cutoffs used. In order to arrive at cutoff-independent results, it was necessary to use a cutoff larger than 6σ at the highest densities and at certain low-temperature state points. (d) State points simulated by standard Newtonian Nosé–Hoover NVT MD dynamics. All results for structure and dynamics reported below were obtained by this dynamics.
In Brownian dynamics, the configurational part of the isochoric specific heat calculated from the fluctuations should have the same value as when calculated from MD. This provides a consistency check ensuring that both methods reproduce the correct canonical ensemble probabilities. Based on this, we found it safe to use MD whenever the temperature is above 1.5 · 10−5 or the density is less five times the melting density at the temperature in question. All structure and dynamics data reported below refer to this region of phase space. At selected state points here, we compared the structure, the density-scaling exponent (Sec. III), the virial potential-energy correlation coefficient (Sec. III), and with the same quantities obtained from Brownian Dynamics and found very good agreement.
The simulations employed standard periodic boundary conditions. All fcc simulations involved 4000 particles, i.e., a cubic lattice. The bcc simulations involved 4394 particles when the cutoff used was or , i.e., a cubic lattice; at larger cutoffs a cubic lattice was used.
III. THE VIRIAL POTENTIAL-ENERGY CORRELATION COEFFICIENT AND THE DENSITY-SCALING EXPONENT
Isomorph theory applies for systems with strong correlations between the constant-volume equilibrium fluctuations of virial W and potential energy U.27–31 Recall that the average virial ⟨W⟩ gives the contribution to the pressure p deriving from particle interactions according to the general equation of state pV = NkBT + ⟨W⟩, in which N is the number of particles and T is the temperature.
Correlations between W and U are investigated numerically by standard Nosé–Hoover NVT MD simulations. For any configuration R defined as the 3N-dimensional vector of all particle positions, the potential energy U(R) is of course known from the simulation. For the EXP pair-potential system, the virial W(R) is given by the following sum over all pairs: .32 Note that the virial is always positive. This applies, in fact, throughout the phase diagram of any purely repulsive pair potential because the pressure is always larger than that of an ideal gas.
Figure 4 shows plots of virial and potential energy equilibrium NVT fluctuations as a function of time with both quantities normalized to unit variance in order to make it easier to estimate the degree of correlation.27 (a) and (b) show results for typical fcc and bcc crystalline state points, while (c) shows the fluctuations at a bcc state point close to that of the re-entrant state point. The latter case shows much weaker correlations.
Normalized plots of virial W and potential energy U deviations from average values at three different state points. The virial potential-energy correlation coefficient R [Eq. (3)] is given in each figure. (a) shows results at a low-density, low-temperature state point of the fcc phase. (b) shows results for a bcc state point. (c) shows results for a bcc state point close to the re-entrant state point [compare Fig. 1(b)]; here, the correlations are much weaker than in (a) and (b).
Normalized plots of virial W and potential energy U deviations from average values at three different state points. The virial potential-energy correlation coefficient R [Eq. (3)] is given in each figure. (a) shows results at a low-density, low-temperature state point of the fcc phase. (b) shows results for a bcc state point. (c) shows results for a bcc state point close to the re-entrant state point [compare Fig. 1(b)]; here, the correlations are much weaker than in (a) and (b).
The Pearson correlation coefficient R quantifying the correlation between W and U is defined by (in which the angular brackets denote constant-volume canonical NVT averages and Δ the deviation from the equilibrium value)
Figure 5 gives data for the virial potential-energy correlation coefficient in which (a) shows by color coding how R varies. Interestingly, the fcc-bcc phase transition does not visibly change R. In most of the phase diagram investigated the EXP system obeys the R > 0.9 criterion designating an R-simple system.17,27,28,33Figure 5(b) shows how R varies along the liquid-solid phase boundaries, plotted as a function of the pressure. Below the horizontal dashed line, R > 0.9; this applies for most of the coexistence state points, involving both fcc and bcc phases.
Virial potential-energy correlation coefficient R in the two solid phases. (a) shows that correlations are generally strong in both the fcc and the bcc crystalline phases with R depending primarily on the density, decreasing as density increases. Only at the highest densities studied, close to the re-entrant density, does the correlation coefficient go below the threshold of 0.9 that defines an R-simple system.27,28 (b) The quantity 1 − R plotted as a function of the pressure along the solid-liquid phase boundary. The green curve represents the liquid phase at coexistence, and the red and blue curves represent the fcc and bcc phases at coexistence, respectively. Note that correlations are somewhat stronger in the solid phases than in the co-existing liquid phase. At the highest pressure, which is above the re-entrant state point, R becomes negative (1 − R > 1).
Virial potential-energy correlation coefficient R in the two solid phases. (a) shows that correlations are generally strong in both the fcc and the bcc crystalline phases with R depending primarily on the density, decreasing as density increases. Only at the highest densities studied, close to the re-entrant density, does the correlation coefficient go below the threshold of 0.9 that defines an R-simple system.27,28 (b) The quantity 1 − R plotted as a function of the pressure along the solid-liquid phase boundary. The green curve represents the liquid phase at coexistence, and the red and blue curves represent the fcc and bcc phases at coexistence, respectively. Note that correlations are somewhat stronger in the solid phases than in the co-existing liquid phase. At the highest pressure, which is above the re-entrant state point, R becomes negative (1 − R > 1).
A more recent definition of an R-simple system is one for which most configurations obey the “hidden scale-invariance” condition34
This expresses the property that if one configuration Ra has a lower potential energy than another one at the same density Rb, this is also the case if the two configurations are scaled uniformly to a different density. Equation (4) only applies for all configurations in the case of perfect virial potential-energy correlations (R = 1),34 which only happens in the non-physical situation of a potential-energy function that is a constant plus an Euler-homogeneous function of the coordinates. Thus, in realistic situations, hidden scale invariance is merely an approximate symmetry.17,33
Isomorphs are defined as curves in the phase diagram of constant excess entropy Sex for an R-simple system (Sex is the entropy minus that of an ideal gas at same density and temperature,36 a quantity that is always negative). All systems have curves of constant Sex, obviously, but only for R-simple systems do these curves have approximately invariant reduced-unit structure and dynamics.30,33,34,37 Hence, the name “isomorphs” for such curves.30
The slope of an isomorph in the logarithmic density–temperature phase diagram, γ, is defined30 by
For an R-simple system, if γ were constant, there would be invariance of the reduced-unit structure and dynamics along the lines of constant ργ/T. For this reason, γ is referred to as the density-scaling exponent.30,33,38 In a computer simulation, one evaluates γ from the constant-density canonical ensemble equilibrium virial and potential-energy fluctuations via the following general statistical–mechanical identity:30
Figure 6 shows how γ varies throughout the two solid phases. We see that γ depends primarily on the density. In fact, the original 2009 version of isomorph theory30 predicted that γ is a function of the density independent of temperature.35 It was later found, however, that the generic version of isomorph theory based on Eq. (4) allows for a more general variation of γ.34 Interestingly, in the gas phase of the EXP system, γ depends more on temperature than on density.19 It is appears that γ depends primarily on the density when a system is in the condensed liquid or solid phases in which any given particle interacts strongly with several nearest neighbors.28,39 For many system, e.g., the Lennard-Jones and related systems, this corresponds to the main R-simple part of the thermodynamic phase diagram because (except at very high temperatures) the gas phase usually does not exhibit strong virial potential-energy correlations. The EXP system is an interesting exception to this that exhibits strong correlations even in the gas phase (whenever kBT ≪ ε).18
The density-scaling exponent γ defined in Eq. (5). γ depends primarily on the density and decreases with increasing density. The same behavior is observed in the condensed-liquid phase of the EXP system. In contrast, at low densities and high temperatures where the EXP system is gas-like, γ depends primarily on the temperature [compare Fig. 3(c) in Ref. 19]. The simplest version of isomorph theory predicts that γ depends only on the density.35 That treatment focused on the dense liquid phase in which each molecule interacts strongly and continuously with several nearest neighbors, which does not apply in the gas phase.
The density-scaling exponent γ defined in Eq. (5). γ depends primarily on the density and decreases with increasing density. The same behavior is observed in the condensed-liquid phase of the EXP system. In contrast, at low densities and high temperatures where the EXP system is gas-like, γ depends primarily on the temperature [compare Fig. 3(c) in Ref. 19]. The simplest version of isomorph theory predicts that γ depends only on the density.35 That treatment focused on the dense liquid phase in which each molecule interacts strongly and continuously with several nearest neighbors, which does not apply in the gas phase.
Figures 7(a) and 7(b) show γ and R as functions of density and temperature. In these figures, “predictions” refer to the theory of Ref. 28, which is summarized in Appendix A. This theory is not exact because it ignores interactions beyond nearest-neighbor particles and assumes that transverse and longitudinal fluctuations of nearest-neighbor vectors are uncorrelated. The former assumption breaks down at high densities. Figure 7(c) plots γ vs R40 for all state points studied, where blue are fcc and red are bcc state points. Paper I in this series18 derived analytically the black dashed line for the dilute gas phase and showed from simulations that all gas and liquid state points fall below this line. We now see that this also applies for the two solid phases. Interestingly, there is an almost one-to-one relation between R and γ for the solid state points, independent of the crystal phase. These data are consistent with an approximate theory ( Appendix A) according to which a single analytical expression may be derived for R and γ, which applies for both crystalline phases. Figure 7(d) shows γ vs 1 − R in a log-log plot. The arrow indicates the direction of lower temperature for data of the same density. The theoretical prediction is approached in this limit for the fcc data (red) whereas the high-density bcc data fall below the prediction. The right part of Fig. 7(d) reproduces data for the fluid phase reported in Ref. 18.
Density-scaling exponent γ and virial potential-energy correlation coefficient R of the two crystalline phases in different representations. Red symbols and lines represent fcc state points, and blue symbols and lines represent bcc state points. (a) γ (left) and R (right) plotted vs density. The dashed and full curves are the predictions of the approximate analytical theory for the T → 0 limit of perfect crystals28 ( Appendix A). The predictions work best at low densities. (b) R as a function of temperature along two fcc isochores, showing that the approximate analytical theory (horizontal lines) works best at low density and low temperature. (c) γ plotted vs R for all simulated state points. The data almost collapse to a single curve common to both crystal structures. For comparison, the black dashed line gives the exact analytical prediction for the zero-density limit of the gas phase.19 The red and blue lines are the fcc and bcc predictions of the approximate analytical theory28 ( Appendix A). The inset provides a blow up of the data with strongest correlations. (d) Log-log plot of γ vs 1 − R. The black dashed line is the gas-phase theory.18 Green lines are data for liquid-phase isochores and purple lines for gas-phase isochores. In the left part of the figure, the identical red and blue lines give the theoretical predictions, while the triangles are numerical data for different isochores with the upper ones giving data for lower densities. The arrows give the direction of decreasing temperature on a given isochore. The isochores connect to liquid-phase isochores of virtually same density (green lines, data from Ref. 18). The dashed-dotted lines represent coexistence state points. The blue (bcc) data points are below the theoretical prediction while the fcc data are above it.
Density-scaling exponent γ and virial potential-energy correlation coefficient R of the two crystalline phases in different representations. Red symbols and lines represent fcc state points, and blue symbols and lines represent bcc state points. (a) γ (left) and R (right) plotted vs density. The dashed and full curves are the predictions of the approximate analytical theory for the T → 0 limit of perfect crystals28 ( Appendix A). The predictions work best at low densities. (b) R as a function of temperature along two fcc isochores, showing that the approximate analytical theory (horizontal lines) works best at low density and low temperature. (c) γ plotted vs R for all simulated state points. The data almost collapse to a single curve common to both crystal structures. For comparison, the black dashed line gives the exact analytical prediction for the zero-density limit of the gas phase.19 The red and blue lines are the fcc and bcc predictions of the approximate analytical theory28 ( Appendix A). The inset provides a blow up of the data with strongest correlations. (d) Log-log plot of γ vs 1 − R. The black dashed line is the gas-phase theory.18 Green lines are data for liquid-phase isochores and purple lines for gas-phase isochores. In the left part of the figure, the identical red and blue lines give the theoretical predictions, while the triangles are numerical data for different isochores with the upper ones giving data for lower densities. The arrows give the direction of decreasing temperature on a given isochore. The isochores connect to liquid-phase isochores of virtually same density (green lines, data from Ref. 18). The dashed-dotted lines represent coexistence state points. The blue (bcc) data points are below the theoretical prediction while the fcc data are above it.
IV. ISOTHERMS AND ISOCHORES
This section presents results for the variation of structure and dynamics along isotherms and isochores of the fcc and bcc phases of the EXP system; the same analysis is carried out along isomorphs in Sec. V. The structure is probed by the radial distribution function (RDF), while the dynamics is probed by the velocity autocorrelation function as well as its Fourier transform giving the phonon density of states.
Figure 8 plots the RDF along four isotherms as a function of the reduced particle distance . (b) and (c) are for the same isotherm in the fcc and bcc phases, respectively; the same applies for (d) and (e). The structure changes significantly along all isotherms, with growing peak heights as density is decreased. Not surprisingly, the lowest temperature isotherm reported in (a) shows the largest and most narrow peaks, reflecting that here the thermal vibrations are smallest. Figure 9 shows data along the same four isotherms for the normalized velocity autocorrelation as a function of the reduced time defined by .30 The insets show the phonon density of states in reduced units, calculated from the velocity autocorrelation function. We conclude that the dynamics also changes significantly along isotherms.
Radial distribution functions (RDFs) along four isotherms, two of which involve both crystal phases. The RDFs are given a function of the reduced particle distance . Peak heights narrow and increase as the density is decreased; due to the use of reduced units, the peak positions do not move visibly. (a) shows results for the T = 5 · 10−7 fcc isotherm. (b) and (c) show results for the T = 3 · 10−5 fcc and bcc isotherms, respectively. (d) and (e) show results for the T = 1 · 10−4 fcc and bcc isotherms, respectively. (f) shows results for the T = 1 · 10−3 bcc isotherm.
Radial distribution functions (RDFs) along four isotherms, two of which involve both crystal phases. The RDFs are given a function of the reduced particle distance . Peak heights narrow and increase as the density is decreased; due to the use of reduced units, the peak positions do not move visibly. (a) shows results for the T = 5 · 10−7 fcc isotherm. (b) and (c) show results for the T = 3 · 10−5 fcc and bcc isotherms, respectively. (d) and (e) show results for the T = 1 · 10−4 fcc and bcc isotherms, respectively. (f) shows results for the T = 1 · 10−3 bcc isotherm.
Normalized velocity autocorrelation function and phonon density of states (inset) along the four isotherms studied in Fig. 8, plotted as a function of the reduced time. The phonon density of states was obtained from the Fourier transform of the normalized velocity autocorrelation function. There is a significant variation along all four isotherms and in both crystal phases. (a) shows results for the T = 5 · 10−7 fcc isotherm. (b) and (c) show results for the T = 3 · 10−5 fcc and bcc isotherms, respectively. (d) and (e) show results for the T = 1 · 10−4 fcc and bcc isotherms, respectively. (f) shows results for the T = 1 · 10−3 bcc isotherm.
Normalized velocity autocorrelation function and phonon density of states (inset) along the four isotherms studied in Fig. 8, plotted as a function of the reduced time. The phonon density of states was obtained from the Fourier transform of the normalized velocity autocorrelation function. There is a significant variation along all four isotherms and in both crystal phases. (a) shows results for the T = 5 · 10−7 fcc isotherm. (b) and (c) show results for the T = 3 · 10−5 fcc and bcc isotherms, respectively. (d) and (e) show results for the T = 1 · 10−4 fcc and bcc isotherms, respectively. (f) shows results for the T = 1 · 10−3 bcc isotherm.
We proceed to present a similar study of how structure and dynamics varies along the isochores, i.e., at constant volume. Figure 10 shows the reduced-unit RDFs along three isochores with (b) and (c) giving data for the fcc and bcc phases of the same isochore. At any given density, the lower the temperature is, the higher are the peaks. Figure 11 shows the normalized velocity autocorrelation function’s variation with reduced time at the same three densities. Just as along the isotherms, there is a significant variation along any given isochore of both structure and dynamics.
RDFs along three isochores as a function of the reduced particle distance. There is a significant variation of structure with higher peaks at lower temperatures. (a) shows results for the ρ = 1 · 10−3 fcc isochore. (b) and (c) show results for the ρ = 1 · 10−2 fcc and bcc isochores, respectively. (d) shows results for the ρ = 1 · 10−1 bcc isochore.
RDFs along three isochores as a function of the reduced particle distance. There is a significant variation of structure with higher peaks at lower temperatures. (a) shows results for the ρ = 1 · 10−3 fcc isochore. (b) and (c) show results for the ρ = 1 · 10−2 fcc and bcc isochores, respectively. (d) shows results for the ρ = 1 · 10−1 bcc isochore.
Normalized velocity autocorrelation function and phonon density of states (insets) along the three isochores studied in Fig. 10. The phonon density was obtained from the Fourier transform of the normalized velocity autocorrelation function. (a) shows results for the ρ = 1 · 10−3 fcc isochore. (b) and (c) show results for the ρ = 1 · 10−2 fcc and bcc isochores, respectively. (d) shows results for the ρ = 1 · 10−1 bcc isochore.
Normalized velocity autocorrelation function and phonon density of states (insets) along the three isochores studied in Fig. 10. The phonon density was obtained from the Fourier transform of the normalized velocity autocorrelation function. (a) shows results for the ρ = 1 · 10−3 fcc isochore. (b) and (c) show results for the ρ = 1 · 10−2 fcc and bcc isochores, respectively. (d) shows results for the ρ = 1 · 10−1 bcc isochore.
V. ISOMORPHS
Section IV reported the variation of structure and dynamics along the EXP system’s solid phase isotherms and isochores. Given the orders of magnitude of variation of density and temperature, not surprisingly both structure and dynamics vary a lot, even when given in reduced units in order to eliminate variations deriving from trivial scaling of time and length. As mentioned, the isomorph theory predicts the existence of lines in the phase diagram of approximately invariant reduced-unit structure and dynamics whenever virial potential-energy correlations are strong, which is the case for almost all state points of the two solid phases [compare Fig. 5(a)].
This section investigates to which degree this prediction holds. At a given temperature and pressure, the average of the fcc and bcc coexistence densities is chosen as the starting point for tracing out isomorphs in the fcc phase (going to lower densities) as well as in the bcc phase (going to higher densities). The isomorphic state points are determined by means of Eqs. (5) and (6) using a fourth-order Runge–Kutta integration in log space.41 Four simulations are needed to determine a new state point in this approach. One of these four simulations refers to the given (starting) state point, the other three are used as correctors to the prediction.41 The simulation at the state point in question was based on 2 × 106 time steps, while each of the three corrector simulations used just 50 000 time steps. The state point simulations were long because these were used also for calculating the RDF, velocity autocorrelation function, R, and γ—the Runge–Kutta integration scheme would work fine with only the short simulation time of the corrector simulations. Three isomorphs were generated by numerically integrating Eq. (5) starting from a state point on the fcc-bcc coexistence line as described above and moving, respectively, to lower and higher densities. Each direction involve one decade of density variation. The Runge–Kutta integration steps involved a density increase of ∼25% when going to higher density, i.e., into the bcc phase, and a density decrease of ∼20% when going to lower densities, i.e., into the fcc phase. The three isomorphs generated in this way and studied below are shown as gray dashed lines in Fig. 1(b), one of which is very close to the melting line. The precise state points are listed in Appendix B.
Before presenting data for the variation of structure and dynamics along selected isomorphs, we investigate the relation between melting line(s) and isomorphs. In the original paper on isomorphs from 2009,30 it was predicted that the melting line is an isomorph. This statement was qualified in 2016,42 and the melting line is now predicted to be an approximate isomorph, a fact that may nevertheless be utilized for calculating the melting-line variation of several quantities. It is generally a good approximation, however, to assume that the melting line is an isomorph. This is illustrated in Fig. 12 in which (a) shows a temperature–pressure plot of the liquid-phase isomorph through the triple point (black dashed line) and the two melting lines (fcc and bcc). The isomorph follows these lines closely, with significant deviations only at the highest densities. Here, the re-entrant point is approached, the correlation coefficient drops significantly (Fig. 5), and isomorph theory breaks down.
Melting line, configurational adiabats, and corresponding log–log slopes. The virial potential-energy correlations are strong at most points; here, the configurational adiabats are isomorphs. (a) Temperature–pressure diagram showing the configurational adiabat through the triple point as the dashed black line and the two melting lines as blue (fcc) and green (bcc) lines. (b) and (c) show as functions of pressure the log–log slopes of the melting/freezing lines and of the configurational adiabats, respectively. Note that the slopes refer to those of the temperature–density phase diagram in which the melting and freezing lines are not exactly identical. Nevertheless, the curves in (b) and (c) are very similar. (d) Density-scaling exponent γ along the two melting lines as a function of the density.
Melting line, configurational adiabats, and corresponding log–log slopes. The virial potential-energy correlations are strong at most points; here, the configurational adiabats are isomorphs. (a) Temperature–pressure diagram showing the configurational adiabat through the triple point as the dashed black line and the two melting lines as blue (fcc) and green (bcc) lines. (b) and (c) show as functions of pressure the log–log slopes of the melting/freezing lines and of the configurational adiabats, respectively. Note that the slopes refer to those of the temperature–density phase diagram in which the melting and freezing lines are not exactly identical. Nevertheless, the curves in (b) and (c) are very similar. (d) Density-scaling exponent γ along the two melting lines as a function of the density.
Figure 12(b) gives the log–log density–temperature slopes of the fcc freezing, bcc freezing, and liquid melting lines, plotted as a function of the pressure. For comparison, Fig. 12(c) gives the log–log density–temperature slopes of the isomorphs (the density-scaling exponents) in the three phases emanating from the triple point. At the fcc–bcc transition, the slope is 2.33. This is consistent with expectations from the inverse-power-law pair potential where the transition takes place around the exponent n = 7 (corresponding to the slope γ = n/3).43,44 Finally, Fig. 12(d) plots the density-scaling exponent γ in the solid phases along the melting lines as a function of the density. The overall conclusion of Fig. 12 is that both the fcc–liquid and the bcc–liquid melting lines are isomorphs to a good approximation and that, consequently, the melting-line slopes are well approximated by the density-scaling exponent. Significant deviations only occur close to the re-entrant point.
Figure 13 demonstrates invariance of the reduced-unit RDFs along six isomorphs generated by Runge–Kutta integration of Eq. (5) into the fcc and bcc phases as described above. We observe a nice collapse along all three isomorphs. A comparison with Figs. 8 and 10 shows that the invariance of structure does not arise simply because structure is invariant throughout the phase diagram.
fcc and bcc RDFs along three isomorphs. (a) and (b) show data for the fcc and bcc isomorph terminating at the ρ = 1.72 · 10−3 coexistence point. (c) and (d) show data for the fcc and bcc isomorph terminating at the ρ = 8.04 · 10−3 coexistence point. (e) and (f) show data for the fcc and bcc isomorph terminating at the ρ = 1.38 · 10−2 coexistence point.
fcc and bcc RDFs along three isomorphs. (a) and (b) show data for the fcc and bcc isomorph terminating at the ρ = 1.72 · 10−3 coexistence point. (c) and (d) show data for the fcc and bcc isomorph terminating at the ρ = 8.04 · 10−3 coexistence point. (e) and (f) show data for the fcc and bcc isomorph terminating at the ρ = 1.38 · 10−2 coexistence point.
Figure 14 shows approximate invariance of the dynamics along the same three isomorphs. The collapse is not as good as for the structure, but it should be kept in mind that a much larger density variation is considered here than in most previous studies. Thus the first data published demonstrating isomorph invariance of the dynamics were for a 3% density change of the highly viscous Kob–Andersen binary Lennard-Jones mixture,30 while each of the figures in Fig. 14 involves a density variation of a full decade. The relevant comparison is with the isotherms of Fig. 9 or the isochores of Fig. 11, revealing a considerably better collapse along the isomorphs.
Normalized velocity autocorrelation function and phonon density of states (inset) along the three isomorphs of Fig. 13. (a) and (b) show data for the fcc and bcc isomorph terminating at the ρ = 1.72 · 10−3 coexistence point. (c) and (d) show data for the fcc and bcc isomorph terminating at the ρ = 8.04 · 10−3 coexistence point. (e) and (f) show data for the fcc and bcc isomorph terminating at the ρ = 1.38 · 10−2 coexistence point.
Normalized velocity autocorrelation function and phonon density of states (inset) along the three isomorphs of Fig. 13. (a) and (b) show data for the fcc and bcc isomorph terminating at the ρ = 1.72 · 10−3 coexistence point. (c) and (d) show data for the fcc and bcc isomorph terminating at the ρ = 8.04 · 10−3 coexistence point. (e) and (f) show data for the fcc and bcc isomorph terminating at the ρ = 1.38 · 10−2 coexistence point.
Figure 15 shows how γ and R vary along the three isomorphs studied in Figs. 13 and 14. There is a substantial variation of γ, which as we have already seen is primarily a function of density. This shows that the EXP system cannot be approximated by an inverse power-law (IPL) pair-potential system. Since isomorphs are only exact for IPL systems, it is all the more remarkable that isomorph invariance of both structure and dynamics applies in the fcc and bcc phases of the EXP system. The reason is that, despite the large variation of γ, the virial potential-energy correlations are very strong (except at the largest densities studied).
γ and R along the three isomorphs studied in Figs. 13 and 14. There is a much larger variation of the density-scaling exponent γ than found for Lennard-Jones type systems; at the same time R is close to unity at most state points.
VI. OUTLOOK
This paper has investigated the two crystalline phases of the EXP system. Except at the highest densities studied, close to the density of the re-entrant transition, the virial potential-energy correlations are very strong in both the fcc and the bcc phases (Fig. 5). Consistent with this, we find excellent invariance of structure along isomorphs and, in view of the large density range simulated, also a quite good invariance of the dynamics.
Only a few studies have been performed on isomorphs in crystals,42,45–48 but the findings of this paper confirm that hidden scale invariance and the concept of isomorphs is not limited to liquids. Compared to other pair potentials, the EXP pair potential is interesting because it exhibits a large variation of the density-scaling exponent γ and has isomorphs in the gas, liquid, and solid phases. In the latter case, the large variation results in two stable crystal structures. Based on the study presented in this paper, it is now possible to include the solid phase in the proposition that the EXP system may be regarded as the mother of all simple pair-potential systems.17 Indeed, an exponential function often plays an important role in empirical as well as ab initio based pair potentials for metals.
An interesting question for future work is what happens at densities above unity, i.e., above the re-entrant point. For the Gaussian-core model, unpublished results of ours show that R may become fairly close to −1, which gives rise to “anti-isomorphs” that have negative slope but still good invariance of the reduced-unit physics. Whether a similar thing happens for the EXP system is not known. If the answer is yes, this may be utilized for interpreting observations on real systems, e.g., dense solutions of star-shaped polymers as well as of metals. In relation to the latter, interestingly it has recently has been proposed that all metals may have a reentrant transition.49
ACKNOWLEDGMENTS
This work was supported by the VILLUM Foundation’s Matter Grant (No. 16515).
APPENDIX A: APPROXIMATE ANALYTICAL THEORY FOR AND R OF CRYSTALS
Figure 7 presents predictions based on an approximate analytical theory developed some time ago.28 The theory, which assumes that all interactions beyond nearest-neighbor pairs may be ignored and that transverse and longitudinal fluctuations of nearest-neighbor vectors are uncorrelated, works as follows (compare Sec. II B.2 in Ref. 28). Let the index b represent nearest-neighbor bonds. If the nearest-neighbor bond b has length rb and the equilibrium bond length is ac, one defines
At any time, the potential energy U and the virial W can both be expressed as linear combinations of S1, S2, S3, …. At low temperatures, S1/ac and have similar variance while the fluctuations of all the other terms are insignificant.28 For the fluctuating parts of U and W, this leads to and . Here, and in which kp is defined by Taylor expanding the pair potential in question as follows, . For the EXP pair potential, if α ≡ ac/σ and β ≡ exp(−ac/σ), one finds . Thus,
If the relative displacement of the pair of particles involved in bond b, ub, is decomposed into a term parallel to and a term perpendicular to the bond vector, ub = ub,∥ + ub,⊥, it follows from Eq. (A1) that for small displacements one has S1 = ∑b|ub,⊥|2/(2ac) and S2 = ∑b|ub,∥|2. Making the Einstein-crystal-like assumption that ub,∥ and ub,⊥ are uncorrelated and assuming that the length squared of the latter is twice that of the former because there are two transverse degrees of freedom, the theory predicts28 that
and that [in Ref. 28, the below factor of 1/3 is erroneously missing from the first line of Eq. (47)]
From this, we find for the EXP pair-potential system
and
Nearest-neighbor interactions dominate in the low-density limit, i.e., when α ≫ 1; in this limit, the theory predicts R ≅ 1 and γ ≅ α/3.
Equations (A5) and (A6) give the predictions shown as full curves in Fig. 7. For the fcc crystal, the α parameter is given by , and for the bcc crystal, α is given by . According to the theory, R and γ are functions of α only. The theory thus predicts the same relation between these R and γ for both crystal phases. This is confirmed qualitatively in Figs. 7(c) and 7(d). As expected, the theory works best at low densities, though the theory is not exact even here.
APPENDIX B: STATE POINTS OF THE THREE ISOMORPHS STUDIED
The three isomorphs studied involve the state points listed in Tables I–III. The state points were determined by integrating Eqs. (5) and (6) by the fourth-order Runge–Kutta algorithm.41
State point data for the highest-temperature isomorph [compare Fig. 1(b)].
ρ − bcc . | T − bcc . | ρ − fcc . | T − fcc . |
---|---|---|---|
1.722 042 × 10−3 | 7.000 195 × 10−5 | 1.722 042 × 10−3 | 7.000 195 × 10−5 |
2.167 922 × 10−3 | 1.151 925 × 10−4 | 1.367 866 × 10−3 | 4.054 469 × 10−5 |
2.729 252 × 10−3 | 1.801 824 × 10−4 | 1.086 535 × 10−3 | 2.216 146 × 10−5 |
3.435 925 × 10−3 | 2.692 507 × 10−4 | 8.630 653 × 10−4 | 1.139 2 × 10−5 |
4.325 573 × 10−3 | 3.857 789 × 10−4 | 6.855 571 × 10−4 | 5.481 866 × 10−6 |
5.445 574 × 10−3 | 5.303 934 × 10−4 | 5.445 574 × 10−4 | 2.455 211 × 10−6 |
6.855 571 × 10−3 | 7.027 732 × 10−4 | 4.325 573 × 10−4 | 1.021 56 × 10−6 |
8.630 653 × 10−3 | 8.988 158 × 10−4 | 3.435 925 × 10−4 | 3.913 932 × 10−7 |
1.086 535 × 10−2 | 1.113 681 × 10−3 | 2.729 252 × 10−4 | 1.374 032 × 10−7 |
1.367 866 × 10−2 | 1.342 329 × 10−3 | 2.167 922 × 10−4 | 4.397 298 × 10−8 |
1.722 042 × 10−2 | 1.570 931 × 10−3 | 1.722 042 × 10−4 | 1.275 095 × 10−8 |
ρ − bcc . | T − bcc . | ρ − fcc . | T − fcc . |
---|---|---|---|
1.722 042 × 10−3 | 7.000 195 × 10−5 | 1.722 042 × 10−3 | 7.000 195 × 10−5 |
2.167 922 × 10−3 | 1.151 925 × 10−4 | 1.367 866 × 10−3 | 4.054 469 × 10−5 |
2.729 252 × 10−3 | 1.801 824 × 10−4 | 1.086 535 × 10−3 | 2.216 146 × 10−5 |
3.435 925 × 10−3 | 2.692 507 × 10−4 | 8.630 653 × 10−4 | 1.139 2 × 10−5 |
4.325 573 × 10−3 | 3.857 789 × 10−4 | 6.855 571 × 10−4 | 5.481 866 × 10−6 |
5.445 574 × 10−3 | 5.303 934 × 10−4 | 5.445 574 × 10−4 | 2.455 211 × 10−6 |
6.855 571 × 10−3 | 7.027 732 × 10−4 | 4.325 573 × 10−4 | 1.021 56 × 10−6 |
8.630 653 × 10−3 | 8.988 158 × 10−4 | 3.435 925 × 10−4 | 3.913 932 × 10−7 |
1.086 535 × 10−2 | 1.113 681 × 10−3 | 2.729 252 × 10−4 | 1.374 032 × 10−7 |
1.367 866 × 10−2 | 1.342 329 × 10−3 | 2.167 922 × 10−4 | 4.397 298 × 10−8 |
1.722 042 × 10−2 | 1.570 931 × 10−3 | 1.722 042 × 10−4 | 1.275 095 × 10−8 |
State point data for the second isomorph.
ρ − bcc . | T − bcc . | ρ − fcc . | T − fcc . |
---|---|---|---|
8.043 829 × 10−3 | 3.550 652 × 10−4 | 8.043 829 × 10−3 | 3.550 652 × 10−4 |
1.012 658 × 10−2 | 4.472 253 × 10−4 | 6.389 441 × 10−3 | 2.743 632 × 10−4 |
1.274 861 × 10−2 | 5.462 811 × 10−4 | 5.075 313 × 10−3 | 2.044 776 × 10−4 |
1.604 955 × 10−2 | 6.481 507 × 10−4 | 4.031 465 × 10−3 | 1.466 612 × 10−4 |
2.020 519 × 10−2 | 7.504 103 × 10−4 | 3.202 306 × 10−3 | 1.009 489 × 10−4 |
2.543 682 × 10−2 | 8.484 276 × 10−4 | 2.543 682 × 10−3 | 6.641 069 × 10−5 |
3.202 306 × 10−2 | 9.388 76 × 10−4 | 2.020 519 × 10−3 | 4.159 616 × 10−5 |
4.031 465 × 10−2 | 1.019 611 × 10−3 | 1.604 955 × 10−3 | 2.473 289 × 10−5 |
5.075 313 × 10−2 | 1.088 987 × 10−3 | 1.274 861 × 10−3 | 1.390 032 × 10−5 |
6.389 441 × 10−2 | 1.144 33 × 10−3 | 1.012 658 × 10−3 | 7.349 984 × 10−6 |
8.043 829 × 10−2 | 1.184 993 × 10−3 | 8.043 829 × 10−4 | 3.642 935 × 10−6 |
ρ − bcc . | T − bcc . | ρ − fcc . | T − fcc . |
---|---|---|---|
8.043 829 × 10−3 | 3.550 652 × 10−4 | 8.043 829 × 10−3 | 3.550 652 × 10−4 |
1.012 658 × 10−2 | 4.472 253 × 10−4 | 6.389 441 × 10−3 | 2.743 632 × 10−4 |
1.274 861 × 10−2 | 5.462 811 × 10−4 | 5.075 313 × 10−3 | 2.044 776 × 10−4 |
1.604 955 × 10−2 | 6.481 507 × 10−4 | 4.031 465 × 10−3 | 1.466 612 × 10−4 |
2.020 519 × 10−2 | 7.504 103 × 10−4 | 3.202 306 × 10−3 | 1.009 489 × 10−4 |
2.543 682 × 10−2 | 8.484 276 × 10−4 | 2.543 682 × 10−3 | 6.641 069 × 10−5 |
3.202 306 × 10−2 | 9.388 76 × 10−4 | 2.020 519 × 10−3 | 4.159 616 × 10−5 |
4.031 465 × 10−2 | 1.019 611 × 10−3 | 1.604 955 × 10−3 | 2.473 289 × 10−5 |
5.075 313 × 10−2 | 1.088 987 × 10−3 | 1.274 861 × 10−3 | 1.390 032 × 10−5 |
6.389 441 × 10−2 | 1.144 33 × 10−3 | 1.012 658 × 10−3 | 7.349 984 × 10−6 |
8.043 829 × 10−2 | 1.184 993 × 10−3 | 8.043 829 × 10−4 | 3.642 935 × 10−6 |
State point data for the third isomorph.
ρ − bcc . | T − bcc . | ρ − fcc . | T − fcc . |
---|---|---|---|
1.384 26 × 10−2 | 1.870 217 × 10−4 | 1.384 26 × 10−2 | 1.870 217 × 10−4 |
1.742 68 × 10−2 | 2.205 215 × 10−4 | 1.099 557 × 10−2 | 1.552 329 × 10−4 |
2.193 905 × 10−2 | 2.535 828 × 10−4 | 8.734 091 × 10−3 | 1.250 395 × 10−4 |
2.761 962 × 10−2 | 2.850 907 × 10−4 | 6.937 735 × 10−3 | 9.757 41 × 10−5 |
3.477 104 × 10−2 | 3.141 555 × 10−4 | 5.510 839 × 10−3 | 7.351 166 × 10−5 |
4.377 415 × 10−2 | 3.395 311 × 10−4 | 4.377 415 × 10−3 | 5.333 891 × 10−5 |
5.510 839 × 10−2 | 3.605 528 × 10−4 | 3.477 104 × 10−3 | 3.712 32 × 10−5 |
6.937 735 × 10−2 | 3.773 441 × 10−4 | 2.761 962 × 10−3 | 2.471 613 × 10−5 |
8.734 091 × 10−2 | 3.879 638 × 10−4 | 2.193 905 × 10−3 | 1.568 369 × 10−5 |
1.099 557 × 10−1 | 3.929 355 × 10−4 | 1.742 68 × 10−3 | 9.452 045 × 10−6 |
1.384 26 × 10−1 | 3.953 509 × 10−4 | 1.384 26 × 10−3 | 5.387 184 × 10−6 |
ρ − bcc . | T − bcc . | ρ − fcc . | T − fcc . |
---|---|---|---|
1.384 26 × 10−2 | 1.870 217 × 10−4 | 1.384 26 × 10−2 | 1.870 217 × 10−4 |
1.742 68 × 10−2 | 2.205 215 × 10−4 | 1.099 557 × 10−2 | 1.552 329 × 10−4 |
2.193 905 × 10−2 | 2.535 828 × 10−4 | 8.734 091 × 10−3 | 1.250 395 × 10−4 |
2.761 962 × 10−2 | 2.850 907 × 10−4 | 6.937 735 × 10−3 | 9.757 41 × 10−5 |
3.477 104 × 10−2 | 3.141 555 × 10−4 | 5.510 839 × 10−3 | 7.351 166 × 10−5 |
4.377 415 × 10−2 | 3.395 311 × 10−4 | 4.377 415 × 10−3 | 5.333 891 × 10−5 |
5.510 839 × 10−2 | 3.605 528 × 10−4 | 3.477 104 × 10−3 | 3.712 32 × 10−5 |
6.937 735 × 10−2 | 3.773 441 × 10−4 | 2.761 962 × 10−3 | 2.471 613 × 10−5 |
8.734 091 × 10−2 | 3.879 638 × 10−4 | 2.193 905 × 10−3 | 1.568 369 × 10−5 |
1.099 557 × 10−1 | 3.929 355 × 10−4 | 1.742 68 × 10−3 | 9.452 045 × 10−6 |
1.384 26 × 10−1 | 3.953 509 × 10−4 | 1.384 26 × 10−3 | 5.387 184 × 10−6 |