Machine learning materials physics: Integrable deep neural networks enable scale bridging by learning free energy functions

TL;DR

This paper introduces Integrable Deep Neural Networks (IDNNs) that learn free energy functions from derivative data, enabling accurate phase field simulations and bridging scales in materials physics.

Contribution

The work presents a novel IDNN approach that can be trained on derivative data and analytically integrated to accurately represent free energy functions in materials modeling.

Findings

IDNN accurately predicts free energy and phase boundaries.

IDNN outperforms B-spline in representing complex free energy landscapes.

Enables scale bridging in materials simulations.

Abstract

The free energy of a system is central to many material models. Although free energy data is not generally found directly, its derivatives can be observed or calculated. In this work, we present an Integrable Deep Neural Network (IDNN) that can be trained to derivative data, then analytically integrated to recover an accurate representation of the free energy. The IDNN is demonstrated by training to the chemical potential values of a binary alloy with B2 ordering. The resulting DNN representation of the free energy is used in a phase field simulation and found to predict the appropriate formation of antiphase boundaries in the material. In contrast, a B-spline representation of the same data failed to represent the physics of the system with sufficient fidelity to resolve the antiphase boundaries.

Equations

$$a_i^{\ell }$$

$$z_i^{\ell }$$

$$Y_i=b_i^{n+1}+\sum _{j=1}^{m_n}{W}_{i,j}^{n+1}{a}_j^n$$

$$\displaystyle\boldsymbol{\hat{W}},\boldsymbol{\hat{b}}=\underset{\boldsymbol{W},\boldsymbol{b}}{\mathrm{arg\,min}}\,\mathrm{MSE}\left(\boldsymbol{Y}(\boldsymbol{x},\boldsymbol{W},\boldsymbol{b})\Big{|}_{\boldsymbol{\hat{x}}_{\theta}},\boldsymbol{\hat{Y}}_{\theta}\right)$$

$$\displaystyle\boldsymbol{\hat{W}},\boldsymbol{\hat{b}}=\underset{\boldsymbol{W},\boldsymbol{b}}{\mathrm{arg\,min}}\,\sum_{k=1}^{n}\mathrm{MSE}\left(\frac{\partial\boldsymbol{Y}(\boldsymbol{x},\boldsymbol{W},\boldsymbol{b})}{\partial x_{k}}\Big{|}_{\boldsymbol{\hat{x}}_{\theta}},\boldsymbol{\hat{y}}_{k_{\theta}}\right)$$

$$\frac{\partial a_i^{\ell }}{\partial x_k}$$

$$a_i^{\ell }$$

$$\frac{\partial z_i^{\ell }}{\partial x_k}$$

$${\beta }_{k_i}^{\ell }$$

$${\alpha }_i^{\ell }$$

$$y_{i,k}$$

$$Y_i$$

$${\beta }_{k_i}^1$$

$$z_i^1$$

$${\alpha }_i^1$$

$$c=\frac{x_1+x_2}{\sqrt{2}}$$

$$\eta =\frac{x_1-x_2}{\sqrt{2}}$$

$$E\left(\vec{\sigma }\right)=V_o+\sum _i{V}_i{\sigma }_i+\sum _{i,j}{V}_{i,j}{\sigma }_i{\sigma }_j+\sum _{i,j,k}{V}_{i,j,k}{\sigma }_i{\sigma }_j{\sigma }_k+...$$

$$Z({\varphi }_c,{\varphi }_{\eta },{\kappa }_c,{\kappa }_{\eta },T)=\sum _{\vec{\sigma }}\exp \left(-\frac{E(\vec{\sigma })+M{\varphi }_c(c(\vec{\sigma })-{\kappa }_c{)}^2+M{\varphi }_{\eta }(\eta (\vec{\sigma })-{\kappa }_{\eta }{)}^2}{k_{B}T}\right)$$

$$\mu_{c}=\frac{\partial g}{\partial c}\Bigr{|}_{(\langle c\rangle,\langle\eta\rangle)}=-2\phi_{c}(\langle c\rangle-\kappa_{c})$$

$$\mu_{\eta}=\frac{\partial g}{\partial\eta}\Bigr{|}_{(\langle c\rangle,\langle\eta\rangle)}=-2\phi_{\eta}(\langle\eta\rangle-\kappa_{\eta})$$

$$\Pi [c,\eta ]$$

$${\mu }_c$$

$${\mu }_{\eta }$$

$$\frac{\partial c}{\partial t}$$

$$\frac{\partial \eta }{\partial t}$$

$$M(c,\eta )=16\tilde{M}{c}^2(\sqrt{2}-c{)}^2(1/2-{\eta }^2{)}^2$$

$$\begin{array}{rl}0& ={\int }_{\Omega }\left[w\frac{\partial c}{\partial t}+4\chi \left(M\nabla w\cdot \nabla g_{,c}+\chi (M{\nabla }^{2}w+\nabla M\cdot \nabla w){\nabla }^{2}c\right)\right]\mathrm{d}V\end{array}$$

$$\begin{array}{rl}0& ={\int }_{\Omega }\left[w\frac{\partial \eta }{\partial t}+\frac{1}{4}\left(Mw{g}_{,\eta }+\chi (M\nabla w+w\nabla M)\cdot \nabla \eta \right)\right]\mathrm{d}V\end{array}$$

$$Y({\xi }_1,{\xi }_2)=\sum _{i=1}^m\sum _{j=1}^n{C}_{ij}{N}_{i,p}({\xi }_1)M_{j,p}({\xi }_2)$$

$$N_{i,p}({\xi }_1)$$

$$N_{i,0}({\xi }_1)$$

$$Y({\xi }_1,{\xi }_2)=\sum _{I=0}^{mn-1}{c}_I{P}_{I,p}({\xi }_1,{\xi }_2)$$

$$Y_k$$

$$A_{kI}=P_{I,p}({\hat{\xi }}_{1_k},{\hat{\xi }}_{2_k})$$

$$B_{kI}$$

Full text

1 Introduction

An accurate description of the free energy plays a critical role in many physics-based models of materials systems. The Euler-Lagrange equations of stationary problems are obtained by requiring that the first variational derivative of the free energy functional vanish. In the example of elasticity, this leads to the equilibrium equation satisfied by the displacement field [1]. Evolution equations can also require first variations of the free energy as inputs. For example, the variational derivatives of the free energy with respect to composition and order parameters define the chemical potential fields used in phase field dynamics [2, 3, 4]. A related result is that the first derivative of the free energy density function with respect to appropriate strain measures gives the conjugate stress. Second derivatives with respect to the temperature yield the heat capacity, those with respect to strains define elasticities [5].

High-fidelity evaluations of free energy density functions are attainable by a combination of atomic scale models (including quantum mechanics and various molecular models) and statistical mechanics. The use of such free energy density functions in the partial differential equations of continuum physics at coarse-grained scales, such as in phase field models, other mass and heat transport equations, and the equations of nonlinear elasticity is a rigorous, if classical approach to scale bridging. In principle, if highly accurate free energy density functions were available, quantitative predictivity would become accessible in a specialized but significant spectrum of mechano-chemically driven phenomena, such as phase transformations, nonlinear and strain gradient elasticity.

Implementing these computations, however, can be challenging for a number of reasons. Free energy data is often computed at individual points rather than as an analytic function directly. As such, it becomes necessary to use a fitting process to create a faithful mathematical model of the free energy. The free energy density can be a function of multiple variables, including composition, temperature, strain, and order parameters, thus leading to a high-dimensional function [6, 7]. Additionally, realistic data for the free energy can contain regions with rapid fluctuations, along with other regions with very gradual slopes [8]. Finally, because of the importance of the derivatives of the free energy, it is desirable that the fitting function be smooth. All of these considerations pose challenges to the fitting technique.

Machine learning methods are readily applicable to this problem, and we specifically consider Deep Neural Networks (DNNs) [9]. DNNs have been successfully applied to problems with large numbers of inputs, such as RGB pixel values in an image [10], the log-power spectra of speech data [11], the configurational energy of multicomponent alloys [12], predicting precipitate morphologies [13] and biomarkers for predicting human age [14], among numerous other applications. Because the activation functions used in DNNs are generally global functions with one local feature, DNNs are capable of capturing local phenomena without negatively affecting longer range attributes of the data. Also, with the proper choice of activation functions, DNNs are infinitely differentiable, thus allowing all necessary derivatives to be computed.

One potential disadvantage of DNNs is that they are not, in general, analytically integrable. This is a particular challenge in the case of fitting free energy data, because the free energy is often not directly measured or computed from atomic models and statistical mechanics. Instead, the derivatives of the free energy (i.e. the chemical potentials) are first observed or computed, then integrated to find the free energy of the system [15, 16, 17, 18, 19, 20, 21, 22, 23, 8]. Such an approach is of particular importance in cases where the chemical potential representation must be integrated with respect to chemical variables to obtain the free energy density, whose derivatives with respect to strain then yield the stress for elasticity [6]. In order to preserve as much information about the derivatives as possible, it is ideal to train directly to the derivative data itself, as opposed to a numerically integrated data set. This requires an alternate form of DNN that can be trained to derivative data, then analytically integrated to represent the free energy itself. We present such an approach in this work.

In this first presentation of our proposed framework, we consider the problem of chemistry, postponing coupled problems for later communications. To demonstrate the method of training a DNN to derivative data of the free energy, we consider a simple binary substitutional alloy with B2 ordering on a BCC parent structure [24]. The resulting free energy, as a function of composition and an order parameter, is then used in the Cahn-Hilliard and Allen-Cahn phase field equations.

In Section 2, we present the form of an Integrable Deep Neural Network (IDNN). Section 3 describes the B2 material system and the method of calculating the chemical potential data via a combination of atomic scale and statistical mechanics computations. In this context, the IDNN training to the chemical potential data and the resulting free energy DNN are shown in Section 3.1. We describe the phase field formulation in Section 3.2 and present the computational results obtained using the DNN representation of the free energy. Section 4 describes a process for fitting the chemical potential data using B-splines and compares the resulting fit and phase field simulation with those of the DNN. Concluding discussions are presented in Section 5.

2 Integrable Deep Neural Network

For a fully connected DNN, the following equations are commonly used to define the activation value of unit $i$ in layer $\ell$, denoted here by $a^{\ell}_{i}$:

[TABLE]

where $b^{\ell}_{i}$ is the bias, $W^{\ell}_{i,j}$ is the weight, $m_{\ell-1}$ is the number of units in layer $\ell-1$, and $f(\cdot)$ is the activation function (see Figure 1). In a DNN used for regression, the output $Y_{i}$ is computed using the activation units from the final layer without an activation function:

[TABLE]

The DNN can be thought of as a function of inputs $\boldsymbol{x}$, weights $\boldsymbol{W}$, and biases $\boldsymbol{b}$, i.e. $\boldsymbol{Y}=\boldsymbol{Y}(\boldsymbol{x},\boldsymbol{W},\boldsymbol{b})$. Training the DNN consists of optimizing the weights and biases to minimize the loss function for a given dataset $\{(\hat{\boldsymbol{x}}_{\theta},\hat{\boldsymbol{Y}}_{\theta})\}_{\theta=1}^{N_{\text{data}}}$. The loss function is generally the mean square error (MSE) for regression problems. We can represent this as follows:

[TABLE]

where the sum over $\theta=1,\dots N_{\text{data}}$ is implied in the definition of the MSE.

As argued in the Introduction, it can be desirable to train a function with data describing the derivatives. While it is possible to train a standard DNN directly to the derivative data, there are at least two drawbacks. The first is that the standard DNN cannot, in general, be analytically integrated to recover the originating function. Furthermore, in the multidimensional case with multiple sets of partial derivative data, the trained standard DNNs representing the partial derivatives would not necessarily be themselves the derivatives of the same function. This inconsistency poses a problem even in situations where mathematical representations of only the partial derivatives are needed and not the integral itself.

Both of these difficulties can be overcome by differentiating the standard DNN with respect to the desired input variables, say $x_{k}$ with $k=1,\ldots,n$, then training these differentiated functions to the (partial) derivative data $\{(\hat{\boldsymbol{x}}_{\theta},\hat{\boldsymbol{y}}_{k_{\theta}})\}$, as represented by the following:

[TABLE]

The functional form that results from differentiating the standard DNN is, of course, analytically integrable. As such, it will be referred to in this work as an IDNN (Integrable Deep Neural Network). The antiderivative of the IDNN is simply a standard DNN with the weights and biases of the trained IDNN. The IDNN has the following form:

[TABLE]

where $z^{\ell}_{i}$ is given by Equation (2). In a slightly more abstracted form with two sets of activation units, $\alpha$ and $\beta_{k}$, we have:

[TABLE]

The values of the trained derivative function $y_{i,k}:=\partial Y_{i}/\partial x_{k}$ and its integral (within an integration constant) $Y_{i}$ are found as follows:

[TABLE]

Note that the following are used to compute the activation values for the first hidden layer:

[TABLE]

The activation function, $f(\cdot)$, can be chosen such that its derivative, $f^{\prime}(\cdot)$, is also a common activation function. For example, with the SoftPlus activation function, $f(x):=\ln(1+e^{x})$, the derivative, $f^{\prime}(x)$ is the commonly used logistic function, i.e. $f^{\prime}(x)=1/(1+e^{-x})$.

We emphasize that the form of the IDNN in Eq. (6)–(8) has been chosen such that its integral has the form of a standard DNN, where both the IDNN and its integral use the same weights and biases, as is clear from the optimizations (5) and (4). Thus, once the weights and biases of an IDNN have been trained using (5), the analytically integrated DNN is constructed by simply using the IDNN’s weights and biases in a standard DNN given by (1) and (2).

2.1 Enforcing symmetries

Neural networks have the property of being uniform approximations of continuous functions over compact domains [25]. However, underlying symmetries of the domain are not guaranteed to be reproduced exactly by the DNN. For instance, consider a function, $f(x,y)$, defined on variables $x$ and $y$, such that they are symmetric under the inversion operation. Given the symmetry of the domain of $(x,y)$, the function is also symmetric under the same operation, i.e. $f(x,y)=f(\pm x,\pm y)$. Any artificial neural network that approximates $f$ must reproduce this symmetry exactly. This can be enforced on the DNN by first transforming the inputs to a set of symmetric functions [26, 27]. For example, within a neural network to approximate $f$, rather than use the values of $(x,y)$ as input, the symmetric functions $(x^{2},y^{2})$ can be used to parameterize the DNN. These functions map all symmetrically equivalent points in the $(x,y)$ space on to the same value, while also differentiating symmetrically distinguishable points. As another example consider functions $f(x,y)$ that are invariant under an inversion about the $y=0$ line. A DNN that approximates any function in this class can be guaranteed to obey the required symmetry by using $(x,y^{2})$ as the input functions.

In general, a set of invariant inputs to the neural network may be defined by generating symmetry invariant polynomial functions with algorithms described by Thomas and Van der Ven [28]. We will denote these functions as $h(\cdot)$. Since the given derivative data is likely differentiated with respect to, for example, $y$ and not $y^{2}$, we must incorporate the symmetrized inputs, $h(\cdot)$, into the loss function, e.g. $\mathrm{MSE}\left(\frac{\partial Y}{\partial h}\frac{\partial h}{\partial x_{1}},\hat{y}\right)$, where $\hat{y}=\partial\hat{Y}/\partial x_{1}$ is the derivative data.

3 Bridging atomic to continuum scales via an IDNN representation of the free energy density

The phase field model is an outgrowth of the Allen-Cahn theory that was originally developed to describe the motion of anti-phase boundaries in ordered compounds. Many binary alloys form a disordered solid solution at elevated temperature that transforms to an ordered compound at low temperature. A well-known example is the order-disorder transition in brass, where an equiatomic alloy of Cu and Zn forms a disordered mixture over the sites of a body centered cubic (bcc) lattice at high temperature, but then orders into a CsCl-type structure at low temperatures. The low temperature CuZn phase can be viewed as an ordered arrangement of Cu and Zn over the sites of the bcc lattice, as illustrated in Figure 3. This ordering on bcc is referred to as B2. Because of the translation symmetry of the bcc lattice, there are two ways to form distinct B2 orderings on the same parent bcc lattice. The Zn could, for example, occupy the corners of the conventional cubic bcc unit cell, while the Cu occupies the body centered sites. Alternatively, Cu could occupy the corners, while Zn occupies the body centered sites. The two orderings are identical, differing only by a rigid translation in space. An alloy that orders upon cooling may adopt the first variant in some regions of the solid and the second variant in other regions. When the two variants impinge, they are separated by an anti-phase boundary.

Order parameters can be used to distinguish regions where the constituents of the alloy are ordered from regions where they are disordered. The order parameters should also be able to distinguish between the different translational variants of the ordered phase. Order parameters that satisfy these conditions for the B2 ordering on bcc are well known [3, 24]. They are defined as linear combinations of sublattice concentrations $x_{1}$ and $x_{2}$ that track the average concentrations over the two sublattice sites $1$ and $2$ of the cubic unit cell shown in Figure 3. In a binary A-B alloy (e.g. A=Cu and B=Zn), each sublattice concentration $x_{i}$ ($i=1,2$) is defined as the fraction of crystal sites belonging to sublattice site $i$ that are occupied by B atoms. Convenient order parameters to measure the degree of B2 long-range ordering can be defined as [3, 24]

[TABLE]

[TABLE]

The first order parameter, $c$, is a measure of the homogeneous composition of the alloy (i.e. proportional to the fraction of all crystal sites occupied by B atoms). The second order parameter, $\eta$, measures the degree of long-range B2-like ordering, and is equal to zero in the completely disordered state as the distinction between the two sublattices disappears in the absence of long-range order and $x_{1}$=$x_{2}$. The order parameter $\eta$ is also able to distinguish between the two translational variants of B2. This becomes clear when considering an alloy containing an equal number of A and B components (i.e. when $c=1/\sqrt{2}$). In one translational variant of B2, the B atoms occupy sublattice $1$ while the A atoms occupy sublattice $2$. For this variant $x_{1}=1$ and $x_{2}=0$ such that $\eta$ becomes equal to $1/\sqrt{2}$. In the second translational variant of B2, the site occupancies are reversed and $\eta$ becomes equal to $-1/\sqrt{2}$.

The free energy, $g$, of a binary A-B alloy with a low-temperature preference for B2 ordering will be a function not only of temperature, $T$, and composition, $x$ (= $c/\sqrt{2}$), but also the long-range order parameter $\eta$. The dependence of $g$ on $T$, $x$ and $\eta$ can be calculated with statistical mechanics approaches. A binary alloy can be modeled as a lattice model that tracks the configurational degrees of freedom associated with all possible ways of arranging A and B atoms over $M$ sites of a crystal. This is done by assigning an occupation variable $\sigma_{i}$ to each crystal site $i=1,...,M$, which is equal to +1 when site $i$ is occupied by a B atoms and -1 when the site is occupied by an A atom. Any arrangement of A and B atoms on the M-site crystal is then fully specified by the configuration vector $\vec{\sigma}=\left(\sigma_{1},...,\sigma_{i},...,\sigma_{M}\right)$. The energy of the crystal for any configurations $\vec{\sigma}$ can be described with a cluster expansion Hamiltonian [16, 17]

[TABLE]

where the expansion coefficients $V_{o}$, $V_{i}$, $V_{i,j}$, $V_{i,j,k}$,… are determined by the chemistry of the alloy and can be parameterized by training to a database of first-principles electronic structure calculations. Here, we use a simplified model cluster expansion Hamiltonian that contains only the constant $V_{o}$, a point term, $V_{i}$ and a term for the nearest neighbor cluster $V_{i,j}$ where ${i,j}$ represents nearest neighbor clusters. The nearest neighbor pair interaction was chosen such that the B2 ordering is a ground state [29]. The finite temperature thermodynamics associated with the order-disorder phase transition can be calculated with Monte Carlo simulations applied to the cluster expansion Hamiltonian. Figure 2 shows the temperature versus composition phase diagram of the model alloy. The B2 stability domain (green) is separated from the disordered solid-solution domain (crimson) by a second order phase transition (solid line).

Conventional grand-canonical Monte-Carlo techniques restrict the calculated free energy to only the thermodynamically stable regions. However, since the unstable parts of the free energy are critical to describing the temporal and spatial evolution of the system within a phase field model, we employed biased Monte-Carlo techniques to sample the free energy throughout the composition, order-parameter domain [24]. The biased ensemble is defined with the following partition function:

[TABLE]

where $E(\vec{\sigma})$, $c(\vec{\sigma})$ and $\eta(\vec{\sigma})$ are the energy, homogeneous composition and long-range order parameter evaluated for configuration $\vec{\sigma}$. The quantities $\phi_{c}$, $\phi_{\eta}$, $\kappa_{c}$ and $\kappa_{\eta}$ define bias potentials, with curvatures given by $\phi_{c}$ and $\phi_{\eta}$, and centers at $\kappa_{c}$ and $\kappa_{\eta}$, respectively. The biased ensemble is sampled with Metropolis Monte-Carlo, and statistical averages of the homogeneous composition $\langle c\rangle$, and order parameter $\langle\eta\rangle$ are measured for different values of $\phi$, $\kappa$, and $T$. The statistical averages can be related to the derivatives of the free energy per atom (denoted $g=G/M$, where $G$ is the total free energy) as [24]:

[TABLE]

[TABLE]

where $\mu_{c}$ and $\mu_{\eta}$ are the chemical potentials with respect to the composition and order parameter respectively. The cluster expansions, and statistical mechanics calculations were performed with the CASM code [30, 31, 32, 21]. The measured ensemble averages were then used to calculate free energy derivatives.

3.1 Training free energy DNN

An IDNN representing the chemical potential data, $\mu_{c}$ and $\mu_{\eta}$, was trained to the derivative data, i.e. the chemical potential values, as described in Section 2, such that the IDNN could be analytically integrated to recover the free energy. The IDNN was implemented as a custom Estimator using the TensorFlow library [33], and was defined by two hidden layers with 10 units per layer. The IDNN was trained for 500 epochs using the AdagradOptimizer, with learning rates of 0.1 and 0.5 applied at different stages of training. A batch size of 10 was used, with 105,061 points in the training set and 35,021 points used for cross-validation. The resulting learning curve is plotted in Figure 4, showing the decrease of both the training and cross-validation mean square errors as training progressed. Symmetry with respect to the order parameter about $\eta=0$ was enforced.

Since the IDNN has a more complex form than the standard DNN, it is reasonable to expect that it may require more training to achieve comparable errors. To demonstrate any differences in training, ten IDNNs and ten standard DNNs were trained to the same chemical potential data, with the same symmetry conditions imposed. All twenty neural networks consisted of two hidden layers of ten neurons, each with different initial values for the weights and biases. They were trained for 20 epochs with a learning rate of 0.1. The resulting learning curves appear in Figure 5. After 20 epochs, the average MSE for the standard DNNs was higher than the average MSE for the IDNNs. Thus, the added complexity of the IDNN does not inhibit the training.

Figure 6 shows the original chemical potential data compared with the associated IDNN that was trained to the corresponding chemical potential data. It also shows two views of the free energy surface as represented by the analytically integrated DNN. Perhaps the most significant feature of the free energy surface is the two energy wells, located at about $c=1/\sqrt{2}$, $\eta=\pm 1/\sqrt{2}$. Given that the wells exist at the same composition, the material will not separate into multiple phases, but instead form anti-phase domains. This reflects the expected physics of the system, described at the beginning of Section 3.

3.2 Phase field computation

To demonstrate the use of the DNN representation of the free energy in computations, the analytically integrated free energy DNN was used in phase field computations. The phase field model was based on the coupled Cahn-Hilliard and Allen-Cahn equations [2, 3, 34, 35, 36], solved using isogeometric analysis (IGA) [37]. The simulation was performed using the mechanoChemIGA code111Code available at github.com/mechanoChem/mechanoChemIGA, which is based on the PetIGA [38] and PETSc [39, 40, 41] libraries, and run on the XSEDE Comet HPC cluster [42].

3.2.1 Formulation

Given the homogeneous free energy density $g(c,\eta)$ as a function of concentration, $c$, and order parameter, $\eta$, we define the total free energy as the following:

[TABLE]

The corresponding chemical potentials are given by the variational derivatives of the total free energy, namely $\mu_{c}:=\delta\Pi/\delta c$ and $\mu_{\eta}:=\delta\Pi/\delta\eta$. Using standard variational methods results in the following equations for the chemical potentials:

[TABLE]

The phase field model consists of the Cahn-Hilliard and Allen-Cahn equations, given by the following, respectively:

[TABLE]

The Cahn-Hilliard equation is in conservation form, with the flux defined as $\boldsymbol{J}:=-M\nabla\mu_{c}$. It models the overall composition of the system through $c$, while conserving mass. The Allen-Cahn equation models the time evolution of the long-range ordering of the system through the non-conserved order parameter $\eta$. The two equations are coupled through the chemical potentials being derived from the same free energy. Periodic boundary conditions were applied. A degenerate mobility of the following form was applied:

[TABLE]

The weak form of the equations, as solved by the IGA formulation, takes the following form:

[TABLE]

Due to the stiffness of the equations when using the realistic free energy density, a two-stage, fourth-order Runge-Kutta time stepping scheme was used [43].

3.2.2 Phase field results

We considered a two-dimensional domain, discretized by a $200\times 200$ element mesh. The initial and boundary value problem was initialized with a uniform composition field of $c=1/\sqrt{2}$ and an order parameter field randomly perturbed about $\eta=0$ by a value up to 0.01, representing a material that had just been quenched from a higher temperature. Periodic boundary conditions were applied. The initial time step was $\Delta t=1$, with an adaptive time stepping scheme being applied to modify the time step based on the convergence of the nonlinear solver at each time step.

As shown in Figure 7, the initially disordered domain gradually forms regions of the two ordered states. The antiphase boundary has completely formed within 40 time steps. The antiphase domains take on order parameter values of $\pm 1/\sqrt{2}$, while the composition field remains nearly uniform.

As the simulation progresses, the microstructure coarsens with curvature fluctuations from a straight antiphase boundary decreasing, as shown in Figure 8.

4 Comparison with B-spline surface fit

For comparison with the IDNN’s representation properties of complex surfaces, we consider two-dimensional B-splines, motivated by their wide use in mathematics, and diverse applications in engineering [37].

4.1 Formulation

The equation for a B-spline surface $Y(\xi_{1},\xi_{2})$, with knots on a tensor grid can be written in the following form

[TABLE]

where $\xi_{1}\in[\xi^{1}_{1},\xi^{m+p+1}_{1}]$, $\xi_{2}\in[\xi^{1}_{2},\xi^{n+p+1}_{2}]$, and $N_{i,p}$ is the B-spline basis function of order $p$. The basis functions are defined by the Cox-de Boor recursion formula [44, 45]

[TABLE]

using the knot vector $\Xi_{1}=\{\xi^{1}_{1},\xi^{2}_{1},\ldots,\xi^{m+p+1}_{1}\}$. $M_{j,p}(\xi_{2})$ is similarly defined using the knot vector $\Xi_{2}=\{\xi^{1}_{2},\xi^{2}_{2},\ldots,\xi^{n+p+1}_{2}\}$.

To convert the matrix $C$ to a vector (to use the standard form for least squares fitting), we use the following index conversion: $I=ni+j$, $i=0,\ldots,m-1$, $j=0,\ldots,n-1$, so that $I=0,\ldots,mn-1$. Then we can rewrite Eq. (30) as the following:

[TABLE]

where $c_{I}:=C_{ij}$ and $P_{I,p}(\xi_{1},\xi_{2}):=N_{i,p}(\xi_{1})M_{j,p}(\xi_{2})$. If evaluating multiple data points $\{(\hat{\xi}_{1_{k}},\hat{\xi}_{2_{k}})\}$, we can write the resulting vector of function evaluations using the following matrix-vector form, written in coordinate notation:

[TABLE]

where

[TABLE]

We now consider fitting to two sets of derivative data. For derivative datasets contained in the two vectors $\boldsymbol{\hat{\mu}}_{1}$ and $\boldsymbol{\hat{\mu}}_{2}$, the following matrices are defined:

[TABLE]

Then, we have the following least squares formulation, with some regularization added for numerical stability:

[TABLE]

where $\lambda$ is the regularization coefficient. Setting the gradient with respect to $\boldsymbol{c}$ equal to the zero vector leads to the following least squares solution:

[TABLE]

where $\boldsymbol{I}$ is the identity matrix. Equations (31-37) and (39) are applied, with datasets in the $(c,\eta)$ space corresponding to the $(\xi_{1},\xi_{2})$ space for B-spline surfaces. Chemical potential data for $\mu_{c}$ and $\mu_{\eta}$ are contained in the vectors $\hat{\boldsymbol{\mu}}_{1}$ and $\hat{\boldsymbol{\mu}}_{2}$, respectively.

4.2 Selection of knots

While the method presented in the previous section can be used to optimize the values of the control points for given knot vectors, the locations of the knots also can be optimized to minimize the error. A variety of approaches are available, including nonlinear least squares, bisection, and genetic algorithms [46]. Our approach was performed in two steps.

First, we divided the data into training data (75%) and validation data (25%). As before, symmetry of the function about $\eta=0$ was imposed. The optimal control points were found for a number of uniformly spaced knots, and the mean square error (MSE) using the validation data was reported, as shown in Figure 9. Higher numbers of knots were not used due to the resulting oscillatory behavior of the fit, particularly in regions of missing data. Knots placed at Chebyshev nodes also resulted in inaccurate fluctuations where data was missing when more than one interior knot was used. In comparing the resulting cross-validation error, it was found that thirteen uniformly distributed interior knots gave the lowest error with $9.25\times 10^{-4}$.

In the second step, we used Matlab’s genetic algorithm optimization routine to attempt to improve the B-spline fit using thirteen interior knots per knot vector (26 total variables). Appropriate inequality constraints were applied to maintain monotonically increasing knot vectors. The algorithm terminated after 50 generations, with a MSE of $5.903\times 10^{-3}$. Since this was not an improvement over the error with uniformly spaced knots, the B-spline fit with thirteen uniformly spaced interior knots per knot vector was taken as the best fit using B-splines.

4.3 B-spline results, and comparison with DNN

The resulting B-spline representations of the chemical potentials and the free energy density are plotted in Figure 10. Visually, these fits seem similar to those of the DNN in Figure 6, although there are notable differences.

The behavior of the dynamics related to antiphase domains is dictated primarily by the sign of $\partial g/\partial\eta$. These values are plotted in Figure 11, along with data points within $c=1/\sqrt{2}\pm 1.1\times 10^{-4}$ for comparison. Antiphase segregation occurs in a domain where the sign of $\partial g/\partial\eta$ changes from positive to negative for increasing $\eta$. In order for antiphase domains to occur in the phase field simulations, the initial conditions, which as in Section 3.2.2 consist of random perturbations about $\eta=0$ in the range (-0.01,0.01), must lie on both sides of the point where the sign of $\partial g/\partial\eta$ changes from positive to negative. While this condition holds true for the DNN representation, it does not hold for the B-spline representation. Due to oscillations in the B-spline representation, $\partial g/\partial\eta$ changes from negative to positive, rather than positive to negative, within the range of initial order parameter values. This results in a phase field solution that does not produce antiphase regions, but instead reaches an equilibrium with uniform values of $\eta=0$ and $c=1/\sqrt{2}$. It is possible, however, that other methods for knot selection or additional constraints might produce a better fit with B-splines.

The B-spline representation is more computationally efficient than the DNN. The phase field code evaluates the free energy and all first and second derivatives at each quadrature point. Using the de Boor algorithm for B-spline evaluation [45], the FLOP count for each function evaluation is about 1700. The highest order term in the count is $\sim 165p^{2}$, with polynomial order $p$. Significantly, the FLOP count for the B-spline evaluation does not depend on the size of the knot vector. The DNN representation, on the other hand, requires about 4500 FLOPs per evaluation of the DNN and its derivatives, using a naive implementation. The highest order term of the count is $\sim 18m^{2}n$, where $m$ is the neurons per layer and $n$ is the number of layers. However, an improved evaluation algorithm could potentially reduce the FLOP count. The use of accelerators could also speed up the function evaluations, without reducing the total FLOP count.

While the FLOP count for the free energy evaluation of the DNN is more than 2.5 times that of the B-spline, this affects only the computation time for the assembly of the residual vector and tangent matrix. The wall time for the matrix-vector solution will be equivalent for the two methods, assuming both fitted functions are similar enough to give comparable condition numbers. For the DNN and B-spline fits presented here, the average computation time per time step over the first 5 time steps using B-splines was 22.8 s, while it was 24.7 s using the DNN. The total computation time using the DNN was, then, only about 1.1 times that of the computation using the B-spline to represent the free energy. For larger problems, the matrix size increases and solver time comes to dominate the average computation time. In this limit the wall times will converge.

5 Conclusions

This communication adds to our nascent, but growing body of work in machine learning and artificial intelligence targeting higher fidelity models of materials physics [12, 13]. We have explored machine learning as an approach to bridging scales, by focusing on the representation of complexity emerging from fine scale physics. Here, the fine scales come from atomic and statistical mechanics descriptions, which can parameterize thermodynamic functions, such as free energy densities with high fidelity. At the core of this work is the idea of an analytically integrable Deep Neural Network (IDNN) to represent such functions. The IDNN is of particular use in the context of mathematically representing the free energy density of a material, where only the derivative data is originally known and the trained function must be integrated to recover the free energy. It is highly relevant to multidimensional systems, where multiple sets of partial derivative data must be trained against simultaneously under the constraint that all trained functions are the partial derivatives of one common function.

Using as a prototypical case a binary alloy, an IDNN was trained to two sets of chemical potential data, which were found using first principles calculations. Since these datasets are derivatives of the same free energy density, the analytic integrability of the IDNN, which exactly preserves consistency of representation, assumes importance. Symmetry with respect to the order parameter was embedded in the IDNN. We anticipate an expansion of these approaches to embedding fundamental aspects of the physics into machine learning models as this field develops.

Phase field simulations with the analytically integrated DNN representing the free energy recovered the proper physics of the system, showing the creation and subsequent coarsening of antiphase domains in the material. This example demonstrates the ability of the IDNN to capture the relevant physics of a material system. Interestingly, the IDNN was able to represent the physics of the system more faithfully than a B-spline representation, even in the current, relatively simple case, of a two-dimensional input. We note that in earlier work, we have demonstrated that B-spline representations themselves are superior to the more traditional Redlich-Kister polynomials at resolving rapidly varying thermodynamic functions [8]. The present work is a continuation of that thread, and establishes that, even for only two-dimensional functions, the B-spline approach also may be inadequate. It is natural to expect that high-dimensional chemical potentials and free energy densities will present greater challenges.

Future work will consider systems of greater complexity, where the free energy density has such higher-dimensional dependence on variables that number $\sim\mathcal{O}(10)$. In this regime, high fidelity representations of the free energy density will be crucial to reproducing the physics through scale bridging, and the uniform approximation property of DNNs, inherited by the IDNNs, will deliver greater advantages.

6 Acknowledgements

We gratefully acknowledge the support of Toyota Research Institute, Award #849910, “Computational framework for data-driven, predictive, multi-scale and multi-physics modeling of battery materials”. This work has also been supported in part by National Science Foundation DMREF grant #1729166, “Integrated Framework for Design of Alloy-Oxide Structures”. Simulations in this work were performed using the Extreme Science and Engineering Discovery Environment (XSEDE) Comet at the San Diego Supercomputer Center through allocations TG-MSS160003 and TG-DMR180072. XSEDE is supported by National Science Foundation grant number ACI-1548562. Computing resources also were provided in part by the NSF via grant 1531752 MRI: Acquisition of Conflux, A Novel Platform for Data-Driven Computational Physics (Tech. Monitor: Ed Walker).