Dynamical continuum simulation of condensed matter from first-principles

TL;DR

This paper introduces a method to automatically generate first-principles equations of state for continuum mechanics simulations using Gaussian processes on ab initio molecular dynamics data, demonstrated on silicon shock wave simulations.

Contribution

It presents a novel approach to derive microscale material properties directly from first-principles calculations for macroscale continuum modeling.

Findings

Successfully generated EOS for silicon from DFT data

Simulated shock wave dynamics using the derived EOS

Demonstrated the method's applicability to hyperelasticity simulations

Abstract

Macroscale continuum mechanics simulations rely on material properties stemming from the microscale, which are normally described using phenomenological equations of state (EOS). A method is proposed for the automatic generation of first-principles unconstrained EOSs using a Gaussian process on a set of ab initio molecular dynamics simulations, thereby closing the continuum equations. We illustrate it on a hyperelasticity simulation of bulk silicon using density-functional theory (DFT), following the dynamics of shock waves after a cylindrical region is instantaneously heated.

Equations

$$(\rho F_{ij}{)}_t$$

$$(\rho u_i{)}_t$$

$$(\rho E{)}_t$$

$$\nabla \times {\mathbf{F}}_i=0.$$

$$\mathbf{G}={\mathbf{F}}^{\mathsf{T}}\mathbf{F}.$$

$$\mathbf{\sigma }=2\rho F_{il}{\left(\frac{\partial E}{\partial G_{lm}}\right)}_S{F}_{jm},$$

$$dE\to \frac{1}{2}V\mathrm{tr}({\mathbf{\sigma }}^{\mathsf{T}}\mathrm{d}\mathbf{G})$$

$${\hat{t}}^{\prime }={\mathbf{k}}^{\mathsf{T}}{\mathbf{C}}^{-1}\mathbf{t},$$

$$C_{ij}=C({\mathbf{x}}^{(i)},{\mathbf{x}}^{(j)})+{\nu }^2{\delta }_{ij},$$

$$C({\mathbf{x}}^{(1)},{\mathbf{x}}^{(2)})={\zeta }^2\exp \left(-\frac{1}{2}\sum _{i,j}\frac{(G_{ij}^{(1)}-G_{ij}^{(2)}{)}^2}{r_{G_{ij}}^2}-\frac{(E_0^{(1)}-E_0^{(2)}{)}^2}{r_{E_0}^2}\right).$$

$$k_i=C({\mathbf{x}}^{(i)},{\mathbf{x}}^{\ast }).$$

$$E^{\prime }(\mathbf{G},E_0)=E(\mathbf{G},E_0)-E_0$$

Full text

Continuum mechanics simulations are of great importance for the simulation of macroscopic condensed matter, from flows in the oil and gas industry Koblitz et al. (2018); Sverdrup et al. (2018) or the modelling of high strain-rate structural deformation Michael and Nikiforakis (2018), to shock waves and detonation in condensed-phase media Schoch et al. (2013a, b, c). The phenomena these simulations describe have their origin in the interactions and dynamics of electrons and nuclei at a scale much smaller than the one of the continuum simulation. Continuum mechanical theories disregard the constituent particles, however, using phenomenological equations of state (EOS) Eliezer et al. (2002), to which a substantial research effort is dedicated (see, e.g. Menikoff (2009); Heuzé (2012); Menikoff (2015); Wilkinson et al. (2017); Weck et al. (2018); Rakhel (2018)). At the atomic scale, electronic structure methods allow the accurate and predictive evaluation of the EOS of condensed matter from first-principles Swift et al. (2001); Correa et al. (2018). Indeed, this was one of its earliest applications, predicting cold pressure-volume curves for isotropic simple materials Cohen (1982). An equation of state is generally much more complex, including all deformations beyond linear, and temperature, for a total of seven dimensions in the case of a solid.

Here we describe a procedure for computing an EOS from first principles in a general way, that neither depends on the material, imposes any functional form, nor requires parameter fitting. This contrasts with the use of traditional EOS’s, such as Mie-Grüneisen Eliezer et al. (2002); Hiermaier (2008). The generality and accuracy of the former is only limited by those of the underlying first-principles theory, although the latter keeps an advantage in efficiency.

Our EOS is constructed using a machine-learning Gaussian process that probes the relevant space by a series of ab-initio molecular-dynamics (AIMD) simulations, which inform the continuum simulations. We illustrate the method with DFT for hyperelastic solids, but the procedure is applicable much more generally in systems and underlying theory. Figure 1 shows results for shock-waves emanating into bulk silicon as described by DFT using the method in this Letter, at the level of nonlinear, anisotropic hyperelasticity. The initial condition was a cylindrical inclusion of a high temperature region, described in more detail below.

For the continuum simulations, the numerical solution of a system of conservation laws is obtained in the Eulerian frame Godunov and Romenskii (1972). A general deformation is represented by the deformation gradient tensor, $\bm{F}=\partial\bm{x}/\partial\bm{X}$, where $\bm{X}$ is a material point in the undeformed configuration, and $\bm{x}$ is its displaced position. The system evolves in time according to Godunov and Romenskii (1972)

[TABLE]

$\rho$ being the mass density, ${\bm{u}}$ the velocity field, $E$ the internal energy, and $\bm{\sigma}$ the Cauchy stress, along with the initial constraint that

[TABLE]

An EOS closes the system giving the internal energy for any deformation. For the definition of a symmetric strain tensor (excluding rotations) we use the right Cauchy–Green tensor,

[TABLE]

The Cauchy stress $\bm{\sigma}$ is given in terms of $\bm{G}$ as

[TABLE]

where the derivative for stress is at constant entropy: this would suggest as most convenient a form of equation of state of $E(\bm{G},S)$. There is, however, no need for obtaining entropy values, which would be expensive in an AIMD setting. The calculation of strain derivatives of the energy for (different values of) constant entropy is what is needed. Instead, we can define the reference tag $E_{0}$ as the internal energy (kinetic energy of nuclei plus total electronic energy) of the material if it were adiabatically brought to the reference configuration from the deformed configuration. As defined, the mapping between $E_{0}$ and entropy does not depend on deformation, and, therefore, $E_{0}$ can be used to label the isentropes. The equation of state is expressed then as $E(\bm{G},E_{0})$, and one obtains $(\partial E/\partial\bm{G})_{S}$ from an entropy preserving (thermally isolated) quasistatic AIMD simulation for a given $E_{0}$. That is, we follow isentropes but do not need to know the value of $S$.

The numerical solution to eqs. 1, 2 and 3 has been extensively discussed Godunov and Romenskii (1972); Plohr and Sharp (1988); Miller and Colella (2001); Barton et al. (2009). A finite volume formulation is used here, with fluxes from the force scheme (Toro, 2013, ch.7) using the mpweno-5 reconstruction Balsara and Shu (2000). The finite-volume formulation in the Eulerian frame allow the capture of correct weak solutions (shock waves). The particular formulation of non-linear elasticity and the method of solution used here illustrate the use of a multidimensional EOS, but the same ideas can be readily ported to other situations.

The EOS used for this simulation was obtained from first-principles simulations of silicon based on DFT. AIMD simulations were performed with the Siesta method and implementation of DFT Soler et al. (2002), using the PBE Perdew et al. (1996) exchange-correlation functional. The basis functions for the valence electrons, and the pseudopotential for the Si core electrons are the same as described in Strickson and Artacho (2016). The mesh used for integrals in real space was well converged with a grid cutoff of $100\text{\,}\,\mathrm{R}\mathrm{y}$. A $2^{3}$ grid of k-points was used on the 64 atom simulations, to give an effective cutoff length of $11\text{\,}\,\mathrm{\SIUnitSymbolAngstrom}$ Moreno and Soler (1992). Thermal electronic contributions are expected to be small at the temperatures considered, and an electronic temperature of $300\text{\,}\,\mathrm{K}$ was used throughout.

Verlet integration (modified as described below) was used to follow an isentrope, with a timestep of $1\text{\,}\,\mathrm{f}\mathrm{s}$ and forces from DFT. 480 separate deformations of a 64-atom box were performed, with AIMD runs of $2\text{\,}\mathrm{p}\mathrm{s}$. For each deformation, an initial configuration was obtained by equilibrating the system using the Tersoff empirical potential Tersoff (1986) on the intended undeformed state, before switching to DFT forces and continuing the integration. The DFT dynamics was further integrated for $250\text{\,}\mathrm{f}\mathrm{s}$ before starting the deformation, and for $250\text{\,}\mathrm{f}\mathrm{s}$ after finishing it, in order to obtain averaged final quantities.

States along an isentrope are extracted directly with molecular dynamics in a slowly deforming box. An alternative procedure is suggested in Ref. Chentsov and Levashov (2012), who use (for a liquid) a sampling in density and temperature before solving an ODE to find the internal energy as a function of density and entropy. In our direct procedure, AIMD gives an isentropically deformed state to a given target deformation, starting from an undeformed reference state at a given (randomly-sampled) $E_{0}$ value, obtained by equilibrating to a given temperature. The box is steadily deformed by slowly varying the box vectors in a linear process from the undeformed to the target deformed state. The entropy change due to varying them is made arbitrarily small by decreasing their rate of variation, since a slowly varied parameter of a Hamiltonian preserves the entropy to first-order in the rate of variation of the parameter (see e.g. Landau and Lifshitz (1980)). To demonstrate that we can follow an isentrope numerically, we show that the process is adiabatic and reversible. That is

[TABLE]

as the time derivative of the deformation vanishes, and additionally, that if the process is reversed, the energy difference between the initial and final states tends to zero. It is important to note that the quantities involved in this expression are equilibrium ensemble averages.

Figure 2 (a) shows both the difference in total energy from this process and the integrated work. For an isentropic process, both should be zero, and any difference is systematic error introduced by the process. Two cases are illustrated: a uniaxial elastic compression of 0.9 relative volume (representative of the deformations we consider), and an uniaxial compression of a liquid, to 0.8 relative volume. The latter case is more challenging since there is additional time for relaxation of the fluid to a hydrostatic stress configuration. We can therefore apply deformations to stresses of tens of $\mathrm{G}\mathrm{P}\mathrm{a}$ over \sim$$1\text{\,}\mathrm{p}\mathrm{s} on 64 atom cells and achieve relative errors in the total energy related to the strain of around 1%, and error in the strain energy computed by integrating the work of 0.2%.

Figure 3 shows slices through the equation of state for silicon used to produce fig. 1 From AIMD we have a discrete sampling of the energy surface. The points must be interpolated to evaluate the energy of a particular arbitrary deformation at a given temperature. A suitable interpolation method and a procedure for choosing the sampling points are crucial components of this scheme. We use Gaussian process regression for the interpolation MacKay (2003); Rasmussen and Williams (2006) for several reasons. First, its ability to handle multi-dimensional data. Second, the fact that (with a suitable covariance function) the interpolated function is smooth: we require the interpolant to have continuous second derivatives, since these appear in expressions for the wave speeds. We thereby avoid unphysical wave splitting. Third, it can incorporate derivate observations (e.g. from pressure) into the learning process, and predict derivatives of the interpolated function (and therefore pressures).

The Gaussian process prediction takes the form

[TABLE]

where $\bm{t}$ is the vector of observed values (total energies and their derivatives with respect to $\bm{G}$) and $\bm{C}$ is the covariance matrix, computed from the training data as

[TABLE]

where $\bm{x}=(\bm{G},E_{0})$ is the vector of inputs, and with the squared exponential covariance function between energy observations,

[TABLE]

The vector $\bm{k}$ contains the covariances of the input to predict, $\bm{x}^{*}$, with each of the observations; that is,

[TABLE]

Covariances between value and derivative observations, and between two derivative observations, are the corresponding derivatives of the $C(\bm{x}^{(1)},\bm{x}^{(2)})$ function.

The interpretation of the hyperparameters in Eq. 10 is as follows: $\zeta$ sets the scale of the inferred function, $\nu$ represents position-independent Gaussian noise in the outcomes that is independent of the inputs, and $r_{G_{ij}}$ is the length scale over which the function varies with $G_{ij}$. Larger values indicate less rapid dependence on the input. Separate noise hyperparameters are used for value and derivative observations.

The sampling is performed by choosing $\bm{G}$ uniformly at random over a problem-specific domain of interest, before converting it to a deformation gradient $\bm{F}$ (by a Cholesky decomposition), and thence to a target lattice $\bm{F}\bm{L}$, where $\bm{L}$ is the matrix whose columns are the lattice vectors. For larger dimensionality other samplings may be more suitable Simpson et al. (2001).

The sampling domain can be chosen generously to include the range over which the deviatoric part of the strain is expected to be less than or equal to the yield criterion, according to, for example, a continuum plasticity model, and with the isotropic part of the strain less than some bound. For the EOS given here, we sample each component independently uniformly over the range $[0.9,1.1]$ for the diagonal components and $[-0.3,0.3]$ for the off-diagonal ones. The internal energy of the undeformed state $E_{0}$ is sampled by varying the initial $T_{i}\in[0,800]$ K. Since $E_{0}$ is the dominant contribution to the internal energy, the fitting is improved by defining $E^{\prime}$

[TABLE]

as the quantity to interpolate.

The error from the reconstruction is shown in fig. 2(b), from an equation of state computed from molecular-dynamics trajectories from an empirical potential, allowing larger sampling. For the databases where gradient information is used, all six components of the gradient are included for one-sixth of the points in the database. The figure shows that this is always beneficial, but much more so for small databases, where it can reduce the error by a factor of four. In addition, if symmetry is exploited, the sampling efficiency is increased by a factor that depends on the crystal system (8 for cubic).

Validation of the multi-scale method proposed in this work is provided by the comparison shown in fig. 4 for properties of silicon shocked with a flat two-dimensional perturbation. The properties obtained here are compared with experimental results as well as with independent simulations results for the same DFT silicon obtained from an ab initio Hugoniot calculation Strickson and Artacho (2016). The agreement is highly satisfactory. In addition, a full, explicit first-principles shock wave has been simulated using AIMD with the same DFT as used here for the EOS. A 2$\times$3$\times$20 supercell with 960 Si atoms was pushed with a piston along the (001) direction with a velocity of 360 m/s. The velocity of the ensuing shock wave calculated with the method described in this Letter was 2% higher than the one obtained from the explicit AIMD simulation, offering again a satisfactory validation of the method.

It should be remembered, however, that the method described here is of a much more general applicability than flat shock waves, while the method in Ref. Strickson and Artacho (2016) is only valid for such shocks, making explicit use of the Hugoniot relations. Figures 1 and 5 illustrate a much more general case that cannot be simulated otherwise, namely, for a shock wave generated in bulk silicon from the sudden heating to $600\text{\,}\mathrm{K}$ of a cylinder of radius $R$ in zero-pressure $40\text{\,}\mathrm{K}$ bulk silicon. The figures show the behavior of the deviatoric stress (Figure 1), the radial material velocity [Figure 5 (a)] and the transverse material velocity [Figure 5 (b)] at a time $t=R/\alpha$ after the initial shock, with $\alpha=10^{4}$ m s*-1*. The initial cylindrical shock is deformed into the displayed shapes due to the anisotropy of the material. There is a scale invariance in the continuum equations that allows $R$ to be macroscopic, which is out of reach for purely atomistic simulations.

In summary, a two-scale method has been demonstrated for the generation of EOSs for macroscopic continuum simulations of condensed matter based on first-principles molecular dynamics. The AIMD simulations are performed on a sample of points selected by a machine-learning Gaussian process in the space of parameters, for the required EOS to be effectively interpolated to any other point, as requested by the continuum mechanics simulation. As a first step, it has been illustrated on complex hyperelastic shock waves in bulk silicon as obtained from DFT calculations, for which the method has been validated. Condensed matter systems of other forms or in other regimes, such as liquids, glasses, polycrystalline solids or solids under plastic deformation, would also be amenable to this method or extensions thereof, using continuum techniques (e.g. assuming yield behaviors for plastic deformation Simo and Hughes (2000)) and MD simulations at larger scales Chandra et al. (2018). The method described in this paper brings first principles to a wide range of continuum mechanics, including for materials that have not been synthesized yet.

Acknowledgments. We are grateful for the support of Dr. Alan Minchinton. OS acknowledges funding from Orica Limited through grant RG63368. Calculations were performed on the Darwin Supercomputer of the University of Cambridge High Performance Computing Service, using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England, and funding from the Science and Technology Facilities Council.