Shock Hugoniot calculations using on-the-fly machine learned force fields with ab initio accuracy

Dynamic shock compression experiments^1–3 are commonly used for investigating material properties and behavior under extreme conditions of temperature and pressure, including those encountered in the warm dense matter (WDM) and hot dense matter (HDM) regimes. The shock Hugoniot, defined as the locus of all possible final shocked states for a given initial state, has a number of applications. These include understanding the dynamic behavior and shock-induced reactions during planetary impacts,⁴ modeling of shock initiation and detonation in polymer-bonded explosives,⁵ and investigation of implosions for inertial confinement fusion.⁶ However, given the significant cost and complexity of such experiments, Hugoniot data for different materials are limited.

Kohn–Sham density functional theory (DFT)^7,8 is among the most widely used electronic structure methods for understanding and predicting the properties of materials from the first principles of quantum mechanics, with no empirical parameters. However, Kohn–Sham DFT calculations are associated with significant computational costs. This is particularly the case for ab initio molecular dynamics (AIMD) simulations, where it is common to solve the Kohn–Sham equations many thousands of times and more. There is also a rapid increase in cost with temperature as the number of partially occupied states increases. This has limited the application of Kohn–Sham DFT for Hugoniot calculations (e.g., Refs. 9–25), with many of these restricted to a small temperature range at the lower end of the spectrum.^{10,13,14,16,19} Path integral Monte Carlo (PIMC)²⁶ is an alternative first-principle method for calculating the Hugoniot.^27–30 However, it is associated with significant computational costs that rapidly increase as the temperature decreases. This makes it impractical at the lower temperatures, particularly for systems with heavier chemical elements.

There have been recent efforts directed at developing solution strategies for Kohn–Sham DFT that are tailored to calculations at high temperature.^31–39 Though they represent significant advances, the methods still entail considerable computational cost, particularly for the AIMD simulations required in shock Hugoniot calculations. This limitation can be overcome using machine learned force field (MLFF) schemes,^40–42 where a surrogate model is developed for the potential energy surface. In particular, in their original form, MLFFs have found application for temperatures up to the lower part of the WDM regime.^43–50 These schemes have been extended to include the prediction of the internal energy,^51–56 making them amenable to the study of the WDM and HDM regimes, including the calculation of the shock Hugoniot.⁵¹ However, these methods require extensive training datasets composed of DFT data from thousands of calculations. Such a process is not only computationally and labor-intensive but also needs to be repeated for different conditions. This limitation can be overcome using on-the-fly MLFF, where training of the model is performed during the molecular dynamics (MD) simulation itself.^57–63 Such a scheme has been applied to calculate transport properties of WDM and HDM.^64,65 However, it does not provide the internal energy, as required for shock Hugoniot calculations.

In this work, we present a framework for calculating the shock Hugoniot using on-the-fly MLFF MD simulations. In particular, we employ an MLFF model based on the kernel method and Bayesian linear regression to compute the free energy, atomic forces, and pressure; along with a model based on linear regression between the internal and free energies to compute the internal energy, with all training data generated from Kohn–Sham DFT. Using this framework, we demonstrate that shock Hugoniot calculations can be accelerated by two orders of magnitude while retaining Kohn–Sham accuracy. We apply this framework to the calculation of the shock Hugoniots for 14 materials, where we find good agreement with the FPEOS database.⁶⁶ In so doing, we confirm that the inter-element interaction in compounds decreases with temperature. The developed framework represents an attractive alternative not only to DFT but also to other more approximate methods used for shock Hugoniot calculations, such as the average-atom method^67–69 and orbital-free DFT.^70,71 In particular, while the accuracy of these more approximate methods decreases with temperature and can be system dependent, the developed framework can be employed seamlessly to calculate the complete Hugoniot across a wide range of systems with ab initio accuracy.

The remainder of this paper is organized as follows: In Sec. II, we discuss the formulation and implementation of the on-the-fly MLFF scheme for shock Hugoniot calculations. In Sec. III, we verify the accuracy and performance of the framework and then apply it to compute the shock Hugoniots for various materials. Finally, we provide concluding remarks in Sec. IV.

In this work, extending recent developments,^62,64 the EOS is calculated using the on-the-fly MLFF MD framework outlined in Fig. 1, which is implemented in parallel within the SPARC electronic structure code,^72,73 building on a prototype version in its serial MATLAB counterpart, M-SPARC.^74,75 The MD simulation starts with a series of Kohn–Sham DFT calculations that provide the initial training data for the free energy, atomic forces, and stresses. The machine learned models are then used to predict the internal energy, forces, and stresses in subsequent MD steps, except when the Bayesian uncertainty/error in the forces so computed exceeds the threshold s_tol, at which point a DFT calculation is performed and the quantities so calculated are added to the training dataset. The threshold s_tol is dynamic, set to the maximum value of the Bayesian error in forces for the MLFF step subsequent to the DFT training step. To mitigate the cubic scaling bottleneck encountered during training, a two-step data selection method is employed: atoms with Bayesian force errors exceeding a set threshold are included in the training dataset, followed by CUR⁷⁶ for downsampling.

The Kohn–Sham training data are generated using the SPARC code, where we employ diagonalization based algorithms for simulations at temperatures less than 100 kK and employ the Spectral Quadrature (SQ) method^31,32,77,78 for higher temperatures. Diagonalization based algorithms solve the Kohn–Sham equations for the occupied orbitals, which scale cubically with the number of orbitals. Since the number of orbitals that need to be calculated increases with temperature, as determined by the Fermi-Dirac distribution, the cost of such a procedure increases rapidly with temperature. In the Gauss SQ method, the electronic density and free energy are written as bilinear forms or sums of bilinear forms, which are then approximated by Gauss quadrature rules. In so doing, the Hellmann–Feynman atomic forces and stresses also become available.^31,32,64,79 On exploiting the locality of electronic interactions,⁸⁰ which manifests itself as the exponential decay of the finite-temperature density matrix away from its diagonal,^81–84 the SQ method scales linearly with system size. The associated prefactor decreases rapidly with temperature, making it the preferred method for Kohn–Sham calculations at high temperature.^64,85–87 Note that in addition to the standard parameters in DFT calculations, there are two additional parameters in the SQ method, namely, the truncation radius and the quadrature order, both of which can be systematically converged to the accuracy desired.

The MLFF model described above only provides access to the free energy and not the internal energy, as required for the calculation of the shock Hugoniot [Eq. (1)]. In this work, we employ an internal energy model based on linear regression between the free and internal energies. We find that higher-order regression provides negligible improvement in the model, and so that linear regression suffices. As in the MLFF model, this internal energy model is trained on the steps at which a Kohn–Sham DFT calculation is performed. An alternative to this strategy is to use a model for the internal energy similar to that of the free energy [Eq. (2)]. However, each atomic configuration provides only the internal energy and not its derivatives with respect to atom positions, which limits the amount of data available to train the model. Given the limited training data available from each atomic configuration for which a DFT calculation is performed, a relatively large number of DFT steps will be required to develop an accurate model for the internal energy, which will significantly increase the cost of the MD simulation. The strategy adopted here provides the desired accuracy, as demonstrated in Sec. III, with negligible increase in the cost of the MD simulation.

We first provide the simulation details in Sec. III A. Next, we verify the accuracy and performance of the developed MLFF formalism in Sec. III B. Finally, we apply the framework to calculate the Hugoniots for various materials in Sec. III C. The data for the Hugoniots so calculated can be found in the  Appendix, while the data for all the EOS calculations used in forming each Hugoniot can be found in the supplementary material.

In calculating the Hugoniot, the pressure is computed for four densities at each of the temperatures chosen. The densities are chosen such that EOS data are available for these points from the FPEOS database, while also being close to the Hugoniot. The point on the Hugoniot is then determined by performing linear interpolation between the density and the residual of the Rankine–Hugoniot relation [Eq. (1)]. Given the proximity of the chosen densities to the corresponding one on the Hugoniot, the quality of the fit is excellent, with $R2≥0.99$ being the coefficient of determination for the linear regression in all cases.

We compute the pressure by performing isokinetic ensemble (NVK) MD simulations with a Gaussian thermostat⁹¹ for 10 000 steps, with the initial 1000 steps used for equilibration and the remaining 9000 used for calculation of the pressure. In so doing, the statistical error in the pressure is reduced to less than 0.1%. The starting atomic configurations for each MD simulation are generated through a preliminary equilibration step in which an on-the-fly MLFF MD simulation of 1000 steps is performed. For the Kohn–Sham DFT calculations performed during this step, we choose not only the grid spacings to be twice larger than those used for training but also other numerical parameters to be significantly less strict, whereby the associated computational cost is negligible compared to the cost of the production MD simulation. This procedure ensures that the machine learned models are trained only after some initial equilibration, which helps to increase the accuracy of the resulting model.

In the Kohn–Sham calculations performed during the on-the-fly MLFF scheme, we use the Perdew–Burke–Ernzerhof (PBE)⁹² exchange-correlation functional, unless otherwise specified. In addition, we use optimized norm-conserving Vanderbilt (ONCV) pseudopotentials,⁹³ chosen from the SPMS set⁹⁴ for temperatures lower than 100 kK, and from previous work^62,85,95 for higher temperatures. The internal energy and pressure for the initial state, i.e., u₀ and P₀, respectively, are calculated using their respective pseudopotentials to ensure that the initial and shocked states have a common reference while using Eq. (1). All the numerical parameters in DFT, including the grid spacing in the standard diagonalization calculations, and the grid spacing, truncation radius, and quadrature order in the SQ calculations, are chosen such that pressures are converged to within 1%. In particular, the truncation radius for any given temperature within the 100 kK and 2 MK range where SQ calculations are performed is relatively independent of the material system, with values of $∼2.7$ bohr and $∼0.75$ bohr at 100 kK and 2 MK, respectively. Given that a similar grid spacing of $∼0.1$ bohr is used for all materials, the quadrature order for any given temperature is also independent of the material system, with values of $∼70$ and $∼40$ at 100 kK and 2 MK, respectively. Note that the errors in the DFT calculations can be reduced as desired by choosing smaller grid spacings, and in the case of SQ, also larger truncation radii and quadrature orders. In the MLFF calculations, the hyperparameters are chosen to be the same as in previous work.^62,64 Indeed, we have found that the accuracy and performance of the MLFF scheme is insensitive to the choice of these parameters for the range of materials systems and conditions studied here.

We now compute the shock Hugoniot for carbon between temperatures of 10 kK and 10 MK, containing both the WDM and HDM regimes. In particular, we consider the following temperatures: 10 kK, 20 kK, 50 kK, 100 kK, 200 kK, 500 kK, 750 kK, 1 MK, 2 MK, 5 MK, and 10 MK. The temperature and density for the initial state are chosen to be 300 K and 3.515 g/cm³, respectively. In the MD simulations, we employ a time step of 0.63 fs for the temperature of 10 kK, adjusting it for other temperatures as the inverse square root of the temperature. We consider system sizes of 500 and 64 atoms for the lowest and highest temperatures of 10 kK and 10 MK, respectively, while linearly interpolating the system size for temperatures in between. We employ the local density approximation (LDA)⁷ for the exchange-correlation functional to facilitate comparison with previous DFT results.

In Fig. 2, we present the shock Hugoniot of carbon so computed using the on-the-fly MLFF scheme, and comparison to recent Kohn–Sham DFT results in the literature.⁸⁵ We observe that there is very good agreement between the MLFF scheme and Kohn–Sham DFT, with a maximum difference of 1.9% in the pressure and 0.4% in the density along the Hugoniot. The corresponding errors in the pressure and internal energy predicted by the MLFF scheme, as determined by averaging the error over all the DFT steps, are 1.4% and 0.005 ha/atom, respectively. In particular, the error in the internal energy is dominated by the error in the MLFF model for the free energy, given the excellent fit in the internal energy model, as shown in Fig. 3. These results indicate that the MLFF scheme can be used to reliably calculate the shock Hugoniot over a wide range of temperatures. Note that it is possible to formulate the MLFF model directly in terms of the pressure, rather than the different components of the stress tensor, as done here. We have found both strategies to provide similar accuracies, e.g., for a temperature of 500 kK, the formalisms in terms of the pressure and stress tensor result in errors of 1.1% and 0.8%, respectively, in the pressure. Note also that the internal energy model can be trained subsequent to the complete MD simulation, i.e., as a post-processing step. We have found that this strategy provides similar accuracy to that adopted here. Indeed, the current strategy of on-the-fly training of the internal energy model accounts for changes in the response of the system over the course of the MD simulation.

In Table I, we present the performance of the on-the-fly MLFF scheme in performing the MD simulations above. We observe that the number of DFT steps constitutes less than 1.4% of the total number of MD steps, with the minimum and maximum percentages being 0.46% and 1.38%, respectively. The occurrence of these Kohn–Sham calculations within the MD simulation can be found in Fig. 4, which is representative of the other systems studied in this work. We observe that majority of the DFT calculations occur early during the MD, becoming less frequent as the simulation continues. Though similar observations have been made for ambient conditions,^57,58,62 a significantly larger number of DFT steps are performed even after a thousand MD steps here, as in previous on-the-fly MLFF MD simulations for transport properties of WDM,⁶⁴ likely due to the increased movement of atoms causing new configurations to be encountered later in the MD simulation. Even with these limited numbers of DFT calculations, they comprise on average 74% of the total CPU time over the different temperatures, with a maximum of 95% for 200 kK and a minimum of 50% for 20 kK. In terms of overall performance, the speedup provided by the on-the-fly MLFF scheme relative to Kohn–Sham DFT is up to two orders of magnitude in both CPU and wall time, with an average speedup by a factor of 62, the minimum and maximum speedups being factors of 44 and 108, respectively. Note that the speedups can be even larger while starting from the model developed for a nearby temperature and density. However, this can impact the accuracy of the resulting model, with larger errors as the difference in the temperature and density increases. Also note that the current MLFF implementation has not been developed to scale in cases where the number of processors exceeds the number of atoms. The MLFF code can be further parallelized to effectively utilize an order of magnitude more processors, by employing an additional level of parallelism over the training dataset descriptors in the covariance matrix calculation and over the radial and angular basis functions in the SOAP descriptor calculation. In addition, the SPARC electronic structure code in default operation can scale to many thousands of processors,^72,73 and the Gauss SQ method can scale to many tens of thousands of processors,^31,77,96 the speedups presented here will also be achieved in the wall time of such simulations.

We now calculate the shock Hugoniot for 14 materials from the FPEOS database, namely, 9 elements: H, He, B, C, O, N, Ne, Mg, and Si; and 5 compounds: LiF, BN, B₄C, MgO, and MgSiO₃, at temperatures between 10 kK and 2 MK, containing both the WDM and HDM regimes. To facilitate comparison, the temperatures are chosen to be the same as those in the FPEOS database.⁶⁶ The temperature for the initial state is chosen to be 300 K, and the densities for H, He, B, C, O, N, Ne, Mg, Si, LiF, BN, B₄C, MgO, and MgSiO₃ are chosen to be 0.085, 0.124, 2.465, 3.515, 0.667, 0.808, 1.507, 1.737, 2.329, 2.535, 2.258, 2.509, 3.570, and 3.208 g/cm³, again as in the FPEOS database. In the MD simulations, we employ a time step of 0.63 fs for C at 10 kK, adjusting it for other temperatures and systems as the inverse square root of the temperature and square root of the mass of the lightest element in the system, respectively. We consider system sizes of 500 and 64 atoms for the lowest and highest temperatures of 10 kK and 2 MK, respectively, while linearly interpolating the system size for temperatures in between. These system sizes are chosen to ensure finite size effects are negligible. Due to stoichiometric constraints in the case of B₄C and MgSiO₃, the smallest system size is taken as 70 atoms.

In Fig. 5, we present the shock Hugoniots for the 14 materials computed by the on-the fly MLFF scheme and compare it with those obtained from the FPEOS database. The FPEOS database has Kohn–Sham DFT data for the lower temperatures and PIMC data for the higher temperatures, generally beyond 1 MK. The script provided within the database is used to generate the FPEOS Hugoniots. We observe that there is very good agreement for both the elements as well as the compounds, with the average difference in the density and pressure across all materials being 0.8% and 2.5%, respectively. The smallest and largest differences between the Hugoniots occur for C and H, respectively, which have been plotted in Fig. 6. The larger differences in H, which occur at the lower temperatures, may be due to the fact that the FPEOS data at the lower temperatures are from PIMC simulations,^28,97 which is a many-body method with different approximations from those adopted in Kohn–Sham DFT. In particular, PIMC calculations become more challenging at the lower temperatures and, therefore, may be associated with larger error bars. We also observe that the larger differences for other materials generally occur at the higher temperatures around 1 MK, which may again be attributed to the fact that the FPEOS data at such temperatures are from PIMC simulations. To verify that the relatively large differences for H are not a consequence of the MLFF model, we have performed Kohn–Sham DFT simulations at the three lowest temperatures of 15.625, 31.25, and 62.5 kK, with corresponding densities along the Hugoniot, i.e., 0.423, 0.369, and 0.365 g/cm³, respectively. The results so obtained are also presented in Fig. 6(a). The errors in the pressure along the Hugoniot at these points are determined to be 0.62%, 0.73%, and 0.71%, respectively, which further verifies the accuracy of the on-the-fly MLFF scheme. To check the sensitivity of the results to the choice of exchange-correlation functional, we also perform on-the-fly MLFF MD simulations with LDA and find that the results agree with those from PBE to within 0.52%, 0.41%, and 0.29%, respectively. To verify the accuracy of the chosen pseudopotential, we also perform on-the-fly MLFF MD simulations with the stringent ONCV pseudopotential for H from the PseudoDOJO library,⁹⁸ and find the pressures agree with those from SPMS to within 0.48%, 0.37%, and 0.46%, respectively. To further verify the accuracy of the pseudopotential, we also perform on-the-fly MLFF MD simulations with a regularized all-electron Coulomb potential⁹⁵ and find that the pressures agree with those from SPMS to within 0.55%, 0.68%, and 0.82%, respectively. Therefore, we find that the DFT results are insensitive to both choice of pseudopotential and choice of exchange-correlation functional; hence, the relatively large differences from PIMC results for H at the lower temperatures are more likely due to the level of theory (semi-local exchange-correlation), larger error bars in the PIMC calculations, and/or other approximations such as nuclear-quantum effects.

In Fig. 7, we present the error in the shock Hugoniots for the compounds, i.e., LiF, BN, B₄C, MgO, and MgSiO₃, when calculated using the linear mixing approximation. In the linear mixing approximation, the internal energy of a compound at a given temperature and pressure is approximated to be the sum of the internal energies of the individual elements, while the density is approximated to be the weighted (by mass fraction) harmonic mean of the individual components. In so doing, we also compute the Hugoniots for Li and F using the on-the-fly MLFF scheme, since these have not been calculated as part of the comparison with FPEOS above. We observe that the linear mixing approximation becomes more accurate with increasing temperature, the errors being as large as 15% for the lower temperatures while becoming smaller than 1% for temperatures beyond 200 kK, which is consistent with previous results in the literature.⁹⁹ This suggests that there is a reduction in inter-element interactions with temperature for the compounds studied here, essentially becoming negligible in the HDM regime.

In this work, we have developed a framework for computing the shock Hugoniot using on-the-fly MLFF MD with ab initio accuracy. In particular, we have employed an MLFF model based on the kernel method and Bayesian linear regression to compute the free energy, atomic forces, and pressure, alongside an internal energy model based on linear regression of the internal and free energies, both trained with Kohn–Sham DFT data. We have verified the accuracy of the formalism by comparing the Hugoniot of carbon in the WDM and HDM regimes with recent Kohn–Sham DFT results in the literature. In so doing, we have demonstrated that Kohn-Sham DFT calculations of the Hugoniot can be accelerated by up to two orders of magnitude, while retaining ab initio accuracy. We have applied this framework to calculate the Hugoniots for 14 materials from the FPEOS database, comprising 9 single elements and 5 compounds, over a temperature range from 10 kK to 2 MK. We have found that the results are in good agreement with the first principles results in the database while providing tighter error bars. In addition, we have verified that the accuracy of the linear mixing approximation for compounds improves as the temperature increases.

Application of the developed framework to the calculation of the shock Hugoniot for heavier elements requires the development of accurate pseudopotentials for such elements and conditions, research we are undertaking presently. The inclusion of temperature and density as inputs to the machine learned models will make them more transferable across different conditions of temperature and density, which will help further accelerate shock Hugoniot calculations, making it a worthy subject for future research. Extending the current MLFF implementation to further increase the scalability on large-scale supercomputers, along with GPU-accelerated computation of essential kernels,¹⁰⁰ promises to decrease wall time significantly in MD simulations, making it a promising area for future research.

See the supplementary material for EOS data at various temperatures and densities for the 14 materials studied in this work, the Hugoniot initial state data, and the Hugoniot data.

S.K. and P.S. gratefully acknowledge support from Grant No. DE-NA0004128 funded by the U.S. Department of Energy (DOE), National Nuclear Security Administration (NNSA). J.E.P gratefully acknowledges support from the U.S. DOE, NNSA: Advanced Simulation and Computing (ASC) Program at Lawrence Livermore National Laboratory (LLNL). This work was performed in part under the auspices of the U.S. DOE by LLNL under Contract No. DE-AC52-07NA27344. This research was also supported by the supercomputing infrastructure provided by Partnership for an Advanced Computing Environment (PACE) through its Hive (U.S. National Science Foundation through Grant No. MRI-1828187) and Phoenix clusters at Georgia Institute of Technology, Atlanta, Georgia.

The authors have no conflicts to disclose.

Shashikant Kumar: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). John E. Pask: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Funding acquisition (equal); Investigation (supporting); Methodology (supporting); Supervision (supporting); Validation (equal); Visualization (supporting); Writing – review & editing (equal). Phanish Suryanarayana: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Software (supporting); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are available within the article and its supplementary material.

The Hugoniot data for the material systems studied in this work are summarized in Table II. The complete EOS data can be found in the supplementary material.