An analytic and complete equation of state for condensed phase materials

Modeling the response of condensed phase materials at high pressure is central to fields such as aerospace, astrophysics, geophysics, and shock and detonation physics. At the continuum level, the governing partial differential equations describe the dynamics of the system where the closure conditions are provided by constitutive models. It is common in high-pressure physics to ignore the contribution of shear deformations to the free energy.¹ Once the yield condition for plastic flow is exceeded, solids begin to behave like fluids and only a thermodynamically consistent equation of state (EOS) is required.² 

Historically, the majority of high-pressure EOS are empirical in nature and rely on limited sets of experimental data to reproduce the relevant physics. Examples of such models include the Davis reactants,³ the linear shock-particle velocity, the Hayes,⁴ or the Tillotson⁵ equations of state. Many of these models are of the Mie–Grüneisen type with the full EOS surface based on a reference curve that can be a Hugoniot, an isentrope, or an isotherm. Mie–Grüneisen EOS forms are a simple linear extrapolation off the reference curve. Presuming one samples near the reference curve, an accurate EOS can be attained, but that might not be true in regions far from the reference.

One of the common deficiencies across these models is the use of a reference curve that leads to stability problems under tension.⁶ Solids and liquids can go into expanded states where the material will eventually fail under certain tensile loading conditions (spallation/cavitation). The EOS must remain qualitatively correct to allow for the dynamic evolution of failure models such as Rajendran–Dietenberger–Grove (RDG)⁷ or tensile plasticity (TEPLA).⁸ If the material enters a region where the EOS has unphysical properties—e.g., negative compressibility—the simulation will break down.

A second limitation stems from a lack of thermal data that forces the use of simplifying assumptions—such as a constant specific heat capacity, $C V$⁠, or a Grüneisen parameter, $Γ$⁠, that is only a function of the specific volume, $V$⁠. These approximations are often reasonable for atomic crystals, such as metals, above room temperature.⁹ If the flow regime is dominated by the compressive pressure, any uncertainties in the value of $C V$ and $Γ$ are usually not critical.¹⁰ Conversely, materials like polymers and energetic molecular crystals are characterized by large molecules with many internal degrees of freedom. In the presence of temperature-dependent chemical reactions, the EOS must be able to capture accurately their thermal response.

Winey et al.¹¹ recognized these limitations and developed a thermodynamically consistent equation of state for liquid nitromethane (NM) using analytic expressions for $C V$⁠, the coefficient of thermal pressure, and the isothermal bulk modulus. They report temperature measurements of NM under stepwise loading conditions using Raman spectroscopy. Despite the relatively large uncertainties (⁠ $±100$ K), they show how simplifying assumptions—such as constant $C V$ or $Γ/V$—can lead to significant inaccuracies and stress the need for additional temperature measurements.

The lack of experimental data for constraining an EOS parameterization at very high temperatures and compressions can be mitigated through the use of atomic-scale calculations. Provided an accurate description of the binding energy and interatomic forces for a material is available, the contributions to the Helmholtz free energy from the static lattice energy and the phonon densities of states can be computed using standard methodologies.^12,13 For crystalline materials, the calculation of the static lattice energy and phonon frequencies for EOS development using density functional theory (DFT)^14,15 are both routine and accurate given the availability well-validated, high-performance codes. Example equations of state include the ones developed by Cawkwell et al.¹⁶ and Sergeev et al.¹⁷ for pentaerythritol tetranitrate (PETN) and Cawkwell et al.¹⁸ for cyclotetramethylene tetranitramine (HMX).

The calculation of the static lattice energies and vibrational frequencies for systems without long-range order, such as liquids, requires instead molecular dynamics (MD) simulations that typically employ a more approximate description of the interatomic forces than that provided by DFT. Nevertheless, fast semi-empirical, empirical, and machine-learning (ML)-based interatomic potentials and force fields for metals and organic materials are widely available for use in these types of simulations.^19,20 Given the availability of high-quality calculated data, there is a need for analytic EOS models that can incorporate these advancements in a computationally efficient manner and remain qualitatively correct outside their calibration domain.

The goal is to derive a complete and thermodynamically consistent equation of state that can accurately reproduce the volumetric and thermal response of a broad class of solids and liquids in a single phase domain. The EOS is constructed from a robust cold curve⁶ and a thermal model associated with the lattice vibrations, and optionally, the contributions to the free energy from delocalized electrons in metals. Since the EOS is specifically designed to be utilized in hydrodynamic codes, it must be analytically invertible in pressure, specific internal energy, and temperature. Parametric conditions are also proposed to obtain thermodynamic stability and convexity in the phase space region of interest.

In Sec. II, we present the derivation of the equation of state using a Helmholtz free energy formalism. The model properties and relevant functional forms are proposed in Sec. III. The static DFT and MD calculations used to obtain data for the calibrations of single-crystal pentaerythritol tetranitrate (PETN), liquid nitromethane, and hexagonal close-packed (HCP) beryllium metal are described in Sec. IV. Additionally, we discuss the limitations of the Dickey–Paskin velocity autocorrelation function²¹ for obtaining the phonon DOS and propose an alternative formulation that yields phonon DOS with the appropriate normalization for EOS development. To demonstrate the qualitative features of the EOS, we show calibration examples for PETN, liquid NM, and Be in Sec. V. The three EOS calibrations have been validated against results from high pressure experiments, including single-shock and isentropic compression data, in Sec. VI. Section VII summarizes the main conclusions and the avenues for future work. The supplementary material includes a self-contained section dedicated to code developers interested in the implementation and verification of the equation of state as well as details of the Z compression experiment.

The functional forms in Eq. (25) and (26) with $n≥0$ guarantee that $Γ(V,T)≥0$ in all phase space. If $n<0$⁠, a separate check should be conducted to verify that $Γ$ takes the appropriate values in the $V$– $T$ regions of interest. In Secs. III D and III E, we use thermodynamic stability and convexity at the infinite temperature limit to derive analytic constraints for the fitting parameters.

The bounds in Eqs. (32) and (33) become useful during the EOS calibration but care must be taken to ensure thermodynamic stability in the $V$– $T$ range of interest for the general case. Section V shows a calibration example where $K T <0$ along the principal isentrope despite satisfying the inequality in Eq. (37). This is generally the case when $Γ 0 ∼ O (1)$ and indicates a deficiency of the model.

Under ambient conditions, PETN (C $5$H $8$N $4$O $12$⁠) adopts a tetragonal unit cell with space group $P 4 ¯ 2 1 c$ containing two molecules.³³ The static compression curve for hydrostatic pressures less than 10 GPa was constructed from the room temperature isothermal compression data measured by Olinger, Cady, and Halleck.³⁴ The DFT calculations used to extend the static compression curve to pressures in excess of 10 GPa and to compute the dependence of the vibrational normal modes on specific volume were described in detail in two earlier publications,^16,35 but here we provide a concise summary for completeness. All of the DFT calculations for PETN were performed using the Gaussian and Plane Waves (GPW) method in the CP2k code³⁶ with the generalized gradient approximation (GGA) exchange-correlation functional of Perdew, Burke, and Ernzerhof (PBE),³⁷ the double zeta with polarization (DZVP) basis set, Goedecker–Teter–Hutter pseudopotentials,³⁸ and the D3(BJ) dispersion correction.³⁹ The plane wave and relative cut-off energies were set to 1200 and 60 Ry, respectively. The high pressure portion was calculated in a single unit cell under periodic boundary conditions with a $2×2×3$ Monkhorst–Pack k-point mesh.⁴⁰ The phonon frequencies were calculated using a brute force (that is, without utilizing crystal symmetry) algorithm to build the force constant matrix via finite displacements of the atoms in a $2×2×2$ periodic supercell with a single k-point. The normal modes were calculated over compressions $1≥V/ V 0 ≥0.6$⁠.

The static lattice energy and dependencies of the vibrational normal modes on compression of liquid NM (CH $3$NO $2$⁠) were computed using MD simulations rather than the static DFT methods because (i) the traditional force-constant matrix approach is ill-suited to liquids and (ii) the absence of long range order would make any DFT calculations prohibitively expensive. The MD simulations were performed using the LATTE quantum-MD code⁴¹ that represents the interatomic forces using semi-empirical density functional tight binding (DFTB) theory.^42,43 We used the lanl31 parameter set⁴⁴ for molecules and materials containing C, H, N, and O with the empirical dispersion correction described in Ref. 45 to improve the description of non-bonded interactions between molecules. All of the MD simulations used a time step for the integration of the equation of motion of 0.25 fs. The extended Lagrangian Born–Oppenheimer MD formalism for the propagation of the electronic degrees of freedom during DFTB-MD allowed us to use only one self-consistent field update per MD time step for the calculation of the set of Mulliken partial charges while maintaining no systematic drift in the total energy.^46,47 The DFTB density matrix was computed using $O ( N 3 )$ matrix diagonalization at zero electronic temperature.

The MD simulations used to calculate the static lattice energy and vibrational modes were performed in cubic simulation cells that contained 80 NM molecules (560 atoms) under three dimensional periodic boundary conditions. The initial atomic configuration was generated by placing NM molecules in the simulation cell at random positions, while avoiding close approaches and overlaps. The dimensions of this supercell were determined by requiring an initial density of 1.14 g cm $− 3$⁠. This configuration was then thermalized at 500 K and constant volume using a Langevin thermostat⁴⁸ for 50 000 time steps to ensure a disordered, liquid-like state was achieved. The volumetric response of liquid NM was examined by generating a series of 61 supercells with compressions $1.05≥V/ V 0 ≥0.45$ by uniformly scaling internal coordinates and lattice vectors of the thermalized, ambient density configuration. Following a short geometry optimization to eliminate any highly distorted covalent bonds, MD trajectories were computed for each system in the canonical ensemble at 298 K using a Langevin thermostat. Information from the first 50 000 time steps from each trajectory was discarded. The time histories of the pressure and atomic velocities was recorded for post-processing from the subsequent 100 000 time steps.

The time average of the instantaneous pressures from the set of MD trajectories yields a 298 K isothermal compression curve up to $V/ V 0 =0.45$⁠. A least-squares fit of Eq. (7) yields an equilibrium density of 1.18 g cm $− 3$⁠, which is slightly higher than the reported experimental values.⁴⁹ The static compression curve is generated from the 298 K isotherm by removing the contribution to the pressure from $F ion (V,T)$⁠.

We start by computing the normal mode frequencies of the nitromethane molecule in the gas phase using the dynamical matrix approach with the lanl31 DFTB model. The frequencies of the 21 normal modes are reported in Table I together with values derived from a DFT calculation at the B3LYP/cc-pVTZ level of theory⁵³ and experiment.⁵⁴ These results show that the lanl31 DFTB model yields normal mode frequencies that are in excellent agreement with experiment and more computationally expensive DFT calculations. These results also highlight that the frequencies and their groupings or degeneracies are connected to the bonding and structure of the molecule. For instance, we observe that the first six modes have imaginary frequencies in the gas phase—these correspond to trivial translations and rotations of the NM molecule. There are three modes with frequencies $> 3000 c m − 1$ that correspond to the three C–H stretches in CH $3$NO $2$⁠. The single C–N stretch corresponds to the mode at about $920 c m − 1$⁠. Indeed, all of the normal modes can be assigned to specific harmonic distortions of the NM molecule.⁵⁴ It is imperative for parameterization of the $F ion$ term in the EOS that the phonon DOS captures accurately not only the frequencies of the normal modes but also their degeneracies based on the structure and symmetry of the molecule.

The cosine transform of Eq. (45) for the condensed-phase MD trajectory computed for the system with density $1.14 g c m − 3$ is presented in Fig. 4. For clarity, the result has been normalized using Eq. (44) such that $N mode =21$⁠. The frequencies associated with each peak in the Fig. 4 correspond clearly to those reported in Table I. However, upon integrating this spectrum, which is presented in Fig. 4, we observe that the number of modes associated with each peak is not in accord with the results in presented in Table I. For instance, we observe about six rather than three C–H stretches and four rather than seven modes at very low frequencies. Hence, the cosine transform of Eq. (45) does not correspond to a phonon DOS.

The phonon DOS, $g(f(V))$⁠, for liquid NM obtained using Eq. (46) from our series of DFTB-MD trajectories over the compression ratios $1.05≥V/ V 0 ≥0.45$ are presented in Fig. 5. The set of DOS show the shift and broadening of the low frequency peak associated with the harmonic motion of the NM molecules as well as the mode-specific positive shifts in the frequencies of the intramolecular vibrational modes as compression is increased.

The static lattice energy and phonon DOS of HCP Be were calculated using the plane wave DFT code VASP.^55–58 All calculations used the PBE exchange-correlation functional, a plane wave cut-off energy of 400 eV with the PREC=ACCURATE setting, and the projector augmented wave pseudopotential for Be.^59,60 The occupancies in the vicinity of the Fermi energy were smeared using the Methfessel–Paxton method with a 0.2 eV smearing width.⁶¹ The static compression curve was obtained by performing a series of lattice parameter optimizations under hydrostatic stresses between $−$20 and 240 GPa with a $30×30×20$ $Γ$-centered Monkhorst–Pack k-point mesh. The lattice parameters predicted by the static DFT calculations at zero pressure, $a=2.265$ Å, $c=3.564$ Å, and $c/a=1.573$ are slightly smaller than experimental measurements at room temperature, $a=2.287$ Å, $c=3.583$ Å, and $c/a=1.567$⁠.⁶² These results were encouraging because they imply that in the case of Be, the calculated static lattice energy is close to the true one and that a correction by the thermal pressure is not required.

The phonon DOS, $g(f(V))$⁠, of HCP Be were calculated via the dynamical matrix method with a $6×6×6$ supercell over compression ratios $0.51≤V/ V 0 ≤1.0$ with a $5×5×3$ k-point mesh and are presented in Fig. 6. These calculations were managed using the phonopy package, which uses the space group symmetry of the crystal to construct the full dynamical matrix from a minimal number of atomic displacements.⁶³ This reduces the computational burden of building the dynamical matrix from $O ( N 4 )$ for a brute force algorithm to $O ( N 3 )$⁠. As compression increases the phonon DOS shift systematically to higher frequencies. The Be pressure-temperature phase diagram shows that HCP Be is thermodynamically stable at low temperature under hydrostatic compression to at least 250 GPa.⁶⁴ The computed phonon DOS are fully consistent with this result because they exhibit no signs of any soft or imaginary phonon modes, which would indicate a structural instability in the HCP phase, at pressures up to 240 GPa.

In this section, we describe the global optimization steps for the thermal and cold components of the model, as well as the different features of each calibration.The devised EOS is intended to capture all three terms of the total free energy in Eq. (41) in the single phase $V$– $T$ range of interest. The calibration is conducted in two sequential steps corresponding to the thermal and cold components of the model. To overcome well-known shortcoming of DFT and MD with regard to its ability to accurately predict the ambient-pressure densities of materials, we constrain the ambient state (⁠ $P a =0$ GPa, $T a =298$ K) during the optimization to match the experimental values of density $ρ a$⁠, thermal expansion coefficient $β a$⁠, and Grüneisen parameter $Γ a$⁠. The values and references for each thermodynamic property are included in Table II. Note that we use the subscript “ $a$” to refer to the ambient state—as opposed to the more common “ $0$”—to avoid confusion with the reference state at $T 0 =0$ K.

All three calibrations are conducted using the basin-hopping global optimization algorithm of Wales and Doye⁷¹ with the implementation included in the LMFIT⁷² Python package. The fits to the normalized specific heat for (a) single-crystal PETN, (b) liquid nitromethane, and (c) HCP beryllium are shown in the top row of Fig. 7 along with their root-mean-square (RMS) deviation. For beryllium, a similar plot for $Γ$ is given in Fig. 8. The markers for the thermal part indicate the DFT data plotted every 5 points in $V$ and $T$ and the solid lines represent the EOS model. The colormap represents different isochores ranging from low (blue) to high (red) compression ratios. Note that the results are only used to demonstrate the properties of the equation of state model and do not account for any uncertainties in the DFT or the experimental results.