Bulk and Point Defect Properties in α-Zr: Uncertainty Quantification on a Semi-Empirical Potential †

Abstract

Modelling studies of irradiation defects in $\alpha$-Zr, such as point defects and their multiple clusters, often use semi-empirical potentials because of their higher computational efficiency as compared to ab initio approaches. Such potentials rely on a fixed number of parameters that need to be fitted to a reference dataset (ab initio and/or experimental), and their reliability is closely related to the uncertainty associated with their parameters, coming from both data inconsistency and model approximations. In this work, parametric uncertainties are quantified on a Second Moment Approximation (SMA) potential, focusing on bulk and point defect properties in $\alpha$-Zr. A surrogate model, based on polynomial chaos expansion, is first built for properties of interest computed from atomistics, and simultaneously allows us to analytically compute the sensitivity indices of the observed properties to the potential parameters. This additional information is then used to select a limited number of material properties for the Bayesian inference. The posterior probability distributions of the parameters are estimated through two Markov Chain Monte Carlo (MCMC) sampling algorithms. The estimated posteriors of the model parameters are finally used to estimate materials properties (not used for the inference): in any case, most of the properties are closer to the reference ab initio and experimental data than those obtained from the original potential.

1. Introduction

Zirconium-based alloys are used as a nuclear cladding material since zirconium (Zr) has a low neutron absorption cross-section combined with satisfactory thermo-mechanical properties and good corrosion resistance. Neutron irradiation triggers the creation of a large number of point defects—both of vacancy and self-interstitial type—that subsequently diffuse and cluster within the crystalline hcp (hexagonal close-packed) Zr matrix. This directly affects the structural and mechanical properties of the cladding material, up to the macroscopic scale [1]. Accurately computing the properties of the multiple possible configurations of point defects and their clusters is a first necessary step towards understanding irradiation effects in $\alpha$-Zr [2,3,4].

The atomic scale is relevant for the study of small defects, with ab initio approaches, like the Density Functional Theory (DFT), being the most popular, as they account for the electronic structure of the materials. Unfortunately, their high computational cost complicates their systematic use. Alternatively, interatomic potentials offer efficient numerical evaluations of energies and forces [5], which then would enable systematic sampling of defect configurations. These potentials are defined by a fixed number of parameters, that are typically identified by a fitting on a user-defined set of target properties. However, they carry a certain amount of uncertainties of quite diverse origins: details of the fitting procedure (cost function, weights of the cost function, etc....), choice of the reference database (experimental error, various DFT values with different exchange-correlation functionals, etc.), and model inadequacy of the interatomic potential itself. Therefore, the accuracy of the physical data predicted using these approaches directly depends on the uncertainties associated with the parameters. As pointed out by Cailliez et al. [6] very few studies quantify the uncertainties associated with the use of such potentials [7,8,9,10,11], and the used potentials are generally simply fitted on a number of properties with no assumed inherent uncertainties (see e.g., [12]).

In this work, we quantify the parametric uncertainties of a Tight-Binding Second Moment Approximation (SMA) potential [13]. We purposefully select a potential that was not adjusted on $\alpha$-Zr properties [14], and that inaccurately predicts vacancy and self-interstitial formation energies. Our methodology combines molecular static simulations, sensitivity analysis, and Bayesian inference. We first build a polynomial chaos-based surrogate model of bulk and point defect properties, which provide analytical expressions for the sensitivity indices of the observed properties to the potential parameters. Then, we identify a limited number of quantities of interest to be used for Bayesian inference, through a careful examination of indices. We calculate the posterior probability distributions of the parameters, using both ab initio and experimental data, and examine the influence of the chosen property for Bayesian inference and of the sampling algorithm. Finally, we use the obtained posteriors to propagate uncertainties from the parameters of the potential up to $\alpha$-Zr bulk and point defect properties and achieve more accurate predictions.

2. Methods & Material

2.1. Second Moment Approximation (SMA) Potential

In the Second Moment Approximation of the tight-binding scheme [13,15], the energy of an atomic site i, with $i=1\dots N_{\mathrm{at}}$ for a system of $N_{\mathrm{at}}$ atoms, can be written as the sum of a band term $E_i^{\mathrm{b}}$ and a repulsive term $E_i^{\mathrm{rep}}$, as

$$E_i^{\mathrm{tot}}=E_i^{\mathrm{b}}+E_i^{\mathrm{rep}}=-\sqrt{\xi \sum _{j\ne i}^{r_{\mathrm{cut}}}{e}^{-2q\left(\frac{r_{ij}}{r_0}-1\right)}}+A\sum _{j\ne i}^{r_{\mathrm{cut}}}{e}^{-p\left(\frac{r_{ij}}{r_0}-1\right)},$$

(1)

where $r_0$ represents a distance close to the first nearest neighbor distance, and $r_{ij}$ is the distance between atoms i and j. $\xi >0$ is an effective hop** integral, and q describes its dependence on the relative interatomic distance. $A>0$ provides the repulsive energy scale, and p accounts for the atomic distance dependence of the repulsive term. The potential is thus based on a set of four adjustable parameters $\theta =\{p,q,A,\xi \}$. The summation in Equation (1) is truncated at a cutoff distance $r_{\mathrm{cut}}$, chosen to include a suitable number of interacting atoms. Its practical implementation involves a fifth-order polynomial, smoothly linking the SMA functional to zero in an interval whose extremities must lie between two successive neighbor shells. The SMA potential described here is part of the ‘smatb’ package of the Lammps molecular dynamics software [16,17].

The set of SMA potential parameters proposed by Dufresne et al. [14] is selected. They were calibrated by fitting some bulk properties of fcc Zr quoted in Refs. [18,19,20], and their values are ${\theta }_0=\{p,q,A,\xi \}=\{7.376,2.492,0.269,2.693\}$, with $r_0=3.17$. For this potential, the $r_{cut}$ interval lies between the second and third neighbors. The nominal potential is not perfect in reproducing small defects properties, as established in Ref. [21].

2.2. Quantities of Interest (QOIs)

We choose a subset of (observed) properties that are relevant for irradiation effects in $\alpha$-Zr. More precisely, we keep elastic constants $C_{ij}$ of the hcp matrix ($Bulk$ modulus, $C_{11}$, $C_{12}$, $C_{44}$, $C_{13}$ and $C_{33}$), and formation energies $E_f$ of single point defects (vacancy V and self-interstitial SIA), that are basic blocks for point defect clusters. Four configurations of SIAs are considered: octahedral (O), basal octahedral (BO), crowdions out of the basal plane (BC’), and split dumbbells in the basal plane (BS), as described in Refs. [22,23,24]. All these are $\mathrm{T}=0K$ quantities, and simulation boxes in molecular static (MS) simulations include ∼1500 atoms, which ensures well-converged formation energy values.

2.3. Uncertainty Quantification (UQ)

A previous study [21] showed that for any given set of parameters $\theta$, the specific choice of the cutoff interval can generate large artifacts when computing the observed properties. Those $r_{\mathrm{cut}}$-induced biases can be corrected by refining—for each random draw of the parameters set—the cutoff interval until elastic properties are stable (see Ref. [21] for details). This procedure requires an additional amount of computational time. To speed up the whole process of UQ/inference, a polynomial chaos expansion-based surrogate model of the corrected QOIs is thus developed, i.e., for the elastic constants $C_{ij}$ and the formation energies of point defects $E_f$. The mathematical background of such an expansion is well described in Ref. [25], and as explained in Ref. [26], the Sobol’s sensitivity indices—that quantify the effect of the potential parameters on the QOIs—are an additional outcome of this expansion.

The uncertainty quantification is carried out using the Bayesian methodology presented in [27]. Using similar notations, we denote $Y_{\mathrm{true}}$ as the true value of a given property Y, which deviates from the experimental or ab initio observation value $Y_{\mathrm{obs}}$ through the following relation:

$$Y_{\mathrm{obs}}=Y_{\mathrm{true}}+{ϵ}_{\mathrm{obs}}.$$

(2)

Similarly, the output of the metamodel $Y_{MM}\left(\theta \right)$ is linked to the true value via the surrogate model error ${ϵ}_{MM}$, as

$$Y_{MM}=Y_{\mathrm{true}}+{ϵ}_{MM}.$$

(3)

Our goal is to find the posterior distribution of the parameter set $\theta$ given the observed quantity $Y_{\mathrm{obs}}$, which is obtained by Bayes’ theorem through the equation:

$$P\left(\theta \right|Y_{\mathrm{obs}})\propto P\left(\theta \right)\times P\left(Y_{\mathrm{obs}}\right|\theta )$$

(4)

with $P\left(\theta \right)$ the prior distribution. To align $Y_{MM}$ closer to $Y_{\mathrm{obs}}$, the error terms ${ϵ}_{\mathrm{obs}}$ and ${ϵ}_{MM}$ are combined into a single error term $ϵ={ϵ}_{\mathrm{obs}}+{ϵ}_{MM}$. Assuming the error $ϵ$ is distributed as a Gaussian (Histograms of residuals support this assumption (not shown here).) with zero mean and variance ${\sigma }_{ϵ}^2$, the likelihood has the following form:

$$P\left(Y_{obs}\right|\theta ,{\sigma }_{ϵ}^2)\propto \exp \left(-\frac{1}{{\sigma }_{ϵ}^2}{[Y_{\mathrm{obs}}-Y_{MM}\left(\theta \right)]}^2\right).$$

(5)

The variance ${\sigma }_{ϵ}^2$ is treated here as another parameter whose posterior probability must be determined. It is assigned a prior half-normal distribution zero-centered and with ${\sigma }_{ϵ}=\frac{{ϵ}_{\mathrm{obs}}+{ϵ}_{MM}}{3}$.

To estimate the posterior distribution, a sampling using a Monte Carlo Markov Chain (MCMC) method is used, which allows us to perform random draws and thus generate drawn samples of the parameter posterior distributions (Equation (4)). Two MCMC algorithms are tested: Metropolis-Hastings (MH) [28] and Hamiltonian Monte Carlo (HMC) [29], both implemented in the Python module pymc3 [30].

The first one works by simulating a Markov Chain whose stationary distribution is the wanted posterior distribution $\pi$. Each element of the chain is randomly drawn from $\pi$ and is accepted or rejected according to a so-called acceptance probability $\alpha$, which is a function of a transition kernel Q, defining the conditional probability of moving to a new position in space given a current position. Many versions of this algorithm exist, depending on the choice of the transition kernel. For this work, we used the classical random walk type, in which Q is symmetrical. For this scheme, a normal proposal function, $\mathcal{N}(0,1)$, and a single parameter update steps are used. The acceptance rate is close to 35% for the Metropolis-Hastings scheme. The second algorithm proceeds by random draw, for each new element of the chain, of both position and momentum. We then let the system evolve according to the Hamiltonian equations up to a new position, which is theoretically always accepted, as the acceptance probability is defined by $min(1,e^{-\Delta H})$ and the dynamics conserve energy along the trajectory (i.e., $\Delta H=0$). Nonetheless, the actual rate of acceptance is close to 65%, since the hamiltonian equations are solved numerically: thus, the $\Delta H$ is not always equal to 0. The leapfrog integrator is used to solve such equations and the time step is automatically scaled through the Non-U-Turns Sampler [31]. The assigned velocity distribution for the HMC runs is $\mathcal{N}(0,1)$. For both algorithms, the convergence is tested through the Gelman-Rubin index, which needs to be lower than 1.1 [32].

3. Results

3.1. Surrogate Model and Sensitivity Indices

To develop the above-mentioned surrogate model of the corrected QOIs, we build an experimental plan for the MS simulations by varying the four potential parameters $\theta =\{p,q,A,\xi \}$ in a maximum range of $\delta =\pm 7\%$ around their nominal values. The plan is generated by randomly drawing 5000 sets of parameters, assuming a uniform probability distribution function (PDF) for all the parameters. Legendre polynomials are chosen for the polynomial expansion since they are well adapted for uniform PDFs [25]. The convergence of the surrogate model is assessed by testing the convergence of both statistical mean and variance for all QOIs. The resulting first order ($S_1$) and total sensitivity indices ($S_T$) [33], quantifying respectively the direct effect and the sum of direct and parameter interactions effect, are given in Figure 1 for the QOIs.

The gap between first order and total indices is small for all QOIs and all parameters, meaning that interaction between parameters in the studied cases is insignificant. Parameters affect the observed properties differently: $E_f^{\mathrm{O}-\mathrm{SIA}}$ is essentially impacted by p, while $E_f^V$ is equally impacted by p, q and $\xi$, with almost no effect of A. Regarding the elastic constants, $\xi$ is the leading sensitivity index, except for $C_{44}$ that shows a stronger influence of p. We also note that for all $C_{ij}$, there is more than one significant influential parameter. As will be seen below, the sensitivity indices give valuable hints regarding which properties are relevant for performing the UQ study.

3.2. Posterior Distributions of Parameters and Uncertainty Quantification

We first evaluate the effect of the choice of properties included for Bayesian inference, and of the chosen MCMC sampling algorithm, on the estimated posterior PDFs of parameters. Then, uncertainty quantification is performed and commented on for all QOIs listed in Section 2.2.

3.2.1. Effect of Chosen QOIs and MCMC Algorithm

We first aim at choosing a subset of QOIs used to estimate the posterior PDFs of parameters. Not including all QOIs is meant to test the propagation of uncertainty for quantities that are not included in the inference, and that could, in a more general sense, be quantities for which there exist no observed values. Results of the sensitivity analysis show the selection of subsets containing three elastic constants ($Bulk$ modulus, $C_{13}$ and $C_{44}$) and $E_f^V$. Indeed, $Bulk=(2C_{11}+C_{33}+4C_{12}+2C_{13})/9$ in hcp structure, and thus encompass information regarding four QOIs at the same time. $C_{13}$ and $C_{44}$ have high p and q sensitivity indices, respectively, and the vacancy is fairly affected by all parameters but A. Input data $Y_{\mathrm{obs}}$ for the inference come from both experiments [34] and ab initio calculations ([12,23,35,36]). Posterior PDFs of parameters are calculated using each chosen property alone (UQ $Bulk$, UQ $C_{13}$, UQ $C_{44}$ and UQ $E_f^V$), and altogether (UQ Tot). For all studied cases, non-informative uniform prior distributions are assigned to the parameters $\mathcal{U}(-7\%,7\%){\theta }_0$, and the Metropolis-Hastings (MH) MCMC scheme is used.

Results are summarized in terms of the mean and standard deviation of the estimated posterior distributions of parameters in

Table 1. Means and standard deviations of posterior distributions of the parameters obtained using data of only one QOI alone (UQ $Bulk$, UQ $C_{13}$, UQ $C_{44}$, and UQ $E_f^V$) or all four QOIs together (UQ Tot.).

		UQ $\mathit{Bulk}$		UQ ${\mathit{C}}_{\mathbf{13}}$		UQ ${\mathit{C}}_{\mathbf{44}}$		UQ ${\mathit{E}}_{\mathit{f}}^{\mathit{V}}$		UQ Tot.
	Nominal	$\mathit{\mu }$	$\mathit{\sigma }$	$\mathit{\mu }$	$\mathit{\sigma }$	$\mathit{\mu }$	$\mathit{\sigma }$	$\mathit{\mu }$	$\mathit{\sigma }$	$\mathit{\mu }$	$\mathit{\sigma }$
p	7.376	7.379	0.293	7.429	0.226	7.390	0.297	7.40	0.558	7.621	0.114
q	2.492	2.488	0.986	2.485	0.1004	2.494	0.0936	2.486	0.189	2.607	0.0402
A	0.269	0.269	0.0107	0.269	0.0108	0.268	0.0106	0.269	0.0217	0.273	0.009
$\xi$	2.693	2.690	0.0866	2.690	0.107	2.700	0.0982	2.671	0.205	2.596	0.057

. For the case of posterior PDFs obtained using data of each property alone, mean values are very similar and are also fairly close to the nominal parameter values. Standard deviations are occasionally more affected (see e.g., the q parameter), but still rather similar. In contrast, the posterior distributions obtained using simultaneous data from $Bulk$ modulus, $C_{13}$, $C_{44}$, and $E_f^V$ are significantly different. Indeed, and as expected, when using data from multiple properties at the same time, more constraints have to be satisfied during the actualization of the posterior distributions. This leads to a shift of the mean and variance of the estimated posteriors against the nominal potential, especially because the QOIs used for this UQ study are different from the ones originally used to fit the nominal potential. As is well-known, the choice of observed properties is crucial for uncertainty quantification, and sensitivity analysis provides additional information about the more relevant QOIs to select. For the following, we thus chose the subset of $Bulk$ modulus, $C_{44}$, $C_{13}$ and $E_f^V$ (i.e., UQ Tot), since this combination of QOIs allows for an overall rather balanced effect of parameters.

With this subset, we now examine the influence of the MCMC algorithm on the posterior PDFs. Converged results are shown in Figure 2 for $\theta =\{p,q,A,\xi \}$. MH and HMC algorithms lead to similar results in terms of shape and standard deviations of the posteriors: p and q have slightly skewed PDFs, A has a close-to-uniform distribution, and $\xi$ is a symmetrical distribution. However, we note a small shift of the mean of the distributions for p and $\xi$ parameters when comparing, in more detail, the outcomes of the two algorithms. This will have consequences on the propagated uncertainty, as discussed below.

3.2.2. Uncertainty Quantification

Distribution of $C_{ij}$ and point defect formation energies (vacancy and four configurations of SIAs) are shown in Figure 3a,b. Results of MH and HMC algorithms are displayed using boxplot representations. Nominal parameter values are indicated with blue crosses, and as visual references, we added some experimental and ab initio values (red and dark crosses, from Refs. [23,34,36]).

In the case of elastic properties and vacancy formation energy, the choice of the algorithm has a negligible impact on the property distributions. It is quite different for the SIAs formation energies: distributions are shifted, and the spreading with HMC is larger, with many outliers. Such differences are likely due to the fact that all SIAs have essentially one strongly influential parameter, p, whose PDF is slightly shifted when comparing MH and HMC. Comparing now both UQ results to the reference values, we find that mostly, median values are closer to the reference ab initio values than the nominal parameter values. In particular, SIA $E_f$ values, that were not included in the inference, have improved predictions. Nevertheless, actual $E_f$ values and their relative ordering are not fully recovered, and this might be related to the expressivity limit of the chosen potential, having difficulties correctly describing this type of defect [21]. Since the median values of the $E_f$ related to the HMC algorithm are closer to the more accurate ab initio input data, this scheme seems to be the most appropriate for the UQ, although it is also the most computationally costly (2 to 3 times longer than the MH algorithm).

4. Conclusions

In this work, we quantified the parametric uncertainties of an SMA potential used for the study of elastic constants and point defect energetics of $\alpha$-Zr, through an approach combining molecular dynamics, polynomial chaos expansion, sensitivity analysis, and Bayesian inference. The main outcomes of this work are:

• Sensitivity analysis gives insights on the relevant properties for performing the inference of posterior PDFs parameters, here a subset of $C_{ij}$ and $E_f^V$,
• Propagation of uncertainty provides more accurate results than the nominal potential—even if not perfect—as compared to the set of reference values. This is particularly visible for SIAs $E_f$, for which predictions by the nominal potential were poor,
• MH and HMC MCMC schemes give qualitatively similar results for the PDFs and property distributions. However, mean values of $E_f$ are closer to the ab initio values for the HMC scheme.

The proposed methodology could be applied to estimate properties of point defect clusters in $\alpha$-Zr: multiple configurations having diverse energetics do exist, and are usually not perfectly captured by interatomic potentials. More generally speaking, the present UQ methodology should be applied systematically more often when using interatomic potentials (e.g., EAM, MEAM, etc.).The approach could also be used with more complex models, such as N-moment tight-binding electronic structure schemes. An application to Zr-H systems will be the topic of future work.