Ab initio vibrational free energies including anharmonicity for multicomponent alloys

TL;DR

This paper introduces a density-functional-theory based method combining thermodynamic integration and machine learning to accurately compute anharmonic vibrational free energies in complex multicomponent alloys.

Contribution

It presents a novel, efficient approach to include anharmonic effects in vibrational free energy calculations for chemically complex alloys.

Findings

Successfully computed anharmonic free energy for VNbMoTaW alloy

Demonstrated the method's accuracy and efficiency

Provides a new tool for thermodynamic modeling of alloys

Abstract

A density-functional-theory based approach to efficiently compute numerically exact vibrational free energies - including anharmonicity - for chemically complex multicomponent alloys is developed. It is based on a combination of thermodynamic integration and a machine-learning potential. We demonstrate the performance of the approach by computing the anharmonic free energy of the prototypical five-component VNbMoTaW refractory high entropy alloy.

Figures

Equations

$$F=F^{\mathrm{ref}}+{\int }_0^{1}d\lambda ⟨E^{\mathrm{DFT}}-E^{\mathrm{ref}}{⟩}_{\lambda },$$

$$M_{n,m}=\sum _j{f}_{n,i,j}(r_{ij})\underset{m\text{ times}}{\underbrace{{\mathbf{r}}_{ij}\otimes {\mathbf{r}}_{ij}\otimes \dots \otimes {\mathbf{r}}_{ij}}}.$$

Full text

Recent developments in the field of multicomponent alloys [high entropy alloys (HEAs), compositionally complex alloys (CCAs)] have opened new materials design perspectives Gludovatz et al. (2014); Li et al. (2016); Murty et al. (2014); Gao et al. (2016). The prediction and exploration of thermodynamic properties and phase stabilities is therefore of critical importance. To this end, parameter-free ab initio calculations, particularly using density-functional theory (DFT), are rapidly gaining popularity IGK2018 . However, the requirement Engin et al. (2008); Grabowski et al. (2011); Zhu et al. (2017) to accurately capture small free energy differences ($\approx 1$ meV/atom) poses severe challenges. Only very recently the required tools to accurately compute free energies of selected unary and ordered binary systems have been developed Duff et al. (2015a); Glensk et al. (2015); Zhu et al. (2017), while efforts to treat the immense chemical complexity of multicomponent alloys are still in their infancy. Here, we propose a highly efficient approach to compute free energies of such multicomponent alloys. We apply it to a prototypical five-component equiatomic body-centered cubic (bcc) refractory VNbMoTaW HEA in its solid solution. This alloy has attracted attention for its superior high-temperature mechanical properties Senkov et al. (2011, 2010).

The free energy is determined by various mechanisms including electronic excitations, configurational entropy and atomic vibrations (e.g., Refs. Zhang et al. (2018); IGK2018 ). The major contribution is due to atomic vibrations (e.g., for VNbMoTaW Körmann et al. (2017): $8.6\,k_{\rm B}$ at 1500 K vs. a configurational entropy of $1.6\,k_{\rm B}$), the leading term of which can be captured by the quasiharmonic approximation. However, the latter accounts only for the phonon softening due to volume expansion and misses out the temperature-dependent phonon softening and broadening. Effective harmonic Hamiltonians Souvatzis et al. (2008, 2009); Hellman et al. (2013); Errea et al. (2014); Sun et al. (2014); Carreras et al. (2017) can approximately account for the temperature induced changes. Numerically exact vibrational free energies can be obtained by thermodynamic integration Alfè et al. (2001, 2002a, 2002b); Vocadlo et al. (2003); Moustafa et al. (2017),

[TABLE]

from a reference potential $E^{\rm ref}$ with free energy $F^{\rm ref}$ to DFT energies $E^{\rm DFT}$, where $\langle\cdots\rangle_{\lambda}$ denotes a thermal average on a mixed potential $E^{\lambda}=\lambda E^{\rm DFT}+(1-\lambda)E^{\rm ref}$. Using a harmonic reference would in principle give the exact anharmonic free energy including temperature-dependent phonon softening and broadening, but such a brute-force integration is computationally prohibitive in practice. The computational effort is dominated by: (1) the number of molecular-dynamics (MD) steps needed to obtain a statistically converged average, (2) the number of $\lambda$-values required to calculate the integral, and (3) the computational effort per MD step.

A state-of-the-art method, making the three steps more feasible, is the two-stage upsampled thermodynamic integration using Langevin dynamics (TU-TILD) method Duff et al. (2015a), which employs an interatomic potential as an intermediate reference in the thermodynamic integration. The thermodynamic integration is split into two stages, first from the harmonic to the reference potential, secondly from the reference potential to full DFT. The intention is to reduce the number of steps necessary to converge the thermal average in the second stage—containing the explicit, costly DFT calculations—by fitting the potential as closely as possible to the DFT data. The brunt of the statistical convergence is then relegated to the thermal average in the first stage, which does not contain explicit DFT calculations and can be thus computed highly efficiently.

The performance of this approach relies critically on the feasibility of fitting a potential that accurately interpolates the DFT data within the thermally accessible phase space. For HEAs and CCAs, where the primary goal is the exploration of the large compositional space and thus many different atomic structures, the requirement of an efficient reference for thermodynamic integration is even more critical. At the same time, for multicomponent alloys the number of fitting parameters drastically increases which results in a serious challenge to fitting reliable potentials. A priori it is not clear whether such an approach is at all feasible for HEAs and CCAs.

A possible solution to this problem could be offered by the emerging class of machine-learning techniques which have recently been developed in various scientific fields Rupp et al. (2018). Several machine-learning potentials have been proposed so far Behler and Parrinello (2007); Behler (2011); Bartók et al. (2010, 2017); Seko et al. (2015); Takahashi et al. (2017). For example, Gaussian process regression was applied to approximate the potential free energy surface of small and medium-sized molecules across the slow degrees of freedom Stecher et al. (2014). First attempts have been put forward to describe alloys Gubaev et al. (2018), focussing on the configurational degree-of-freedom of a ternary system. Whether a machine-learning approach is applicable to compute the vibrational free energy of chemically even more complex bulk materials is so far unknown.

In this work we develop a new algorithm combining the TU-TILD method with moment tensor potentials (MTP), a class of machine-learning potentials first proposed in Ref. Shapeev (2016), and recently shown to perform best among different machine-learning models comp . We demonstrate here that the TU-TILD+MTP combination is an ideal symbiosis for an efficient and accurate calculation of the full vibrational free energy of disordered multicomponent alloys. MTP describes the atomic environment of the $i$th atom by the moments of inertia of the neighboring atoms,

[TABLE]

Here the radial functions $f_{n,i,j}(r_{ij})$, $n=1,2,\ldots$, define different shells around the $i$th atom; the contribution of the $j$th atom to the $n$th shell can depend on the types of the $i$th and $j$th atoms. When $m=0$, $M_{n,0}$ is a scalar quantity interpreted as the weight of the $n$th shell. The set of these scalar descriptors is not complete. However, this set can be made complete by adding vectorial “eccentricity” of the $n$th shell, the tensor of second moments of inertia $M_{n,2}$, the third moments $M_{n,3}$, etc. Hence, MTP can approximate an arbitrary local interaction energy by forming basis functions as different ways of contracting these tensor-valued moment descriptors to a scalar and considering an arbitrary linear combination of these basis functions with parameters fitted from data Shapeev (2016); Gubaev et al. (2018). In practice, this means that we can increase the number of parameters until the fitting error stops decreasing—this would indicate that we have reached a lower error bound that a local model can achieve with a given cutoff radius.

We now apply an MTP as a reference potential within TU-TILD for the disordered VNbMoTaW HEA. Chemical disorder is modeled by a special quasirandom structure (SQS) in a 125 atomic supercell Zunger et al. (1990). DFT calculations are performed with vasp Kresse and Furthmüller (1996a, b) employing PAW Blöchl (1994), GGA-PBE Perdew et al. (1996), and including electron-phonon coupling Zhang et al. (2017). See the supplement for details sup . We consider temperatures up to 3000 K which is close to the estimated melting point Senkov et al. (2010). In accord with our previous works Zhu et al. (2017); Gong et al. (2018) we use DFT MD simulations at several volumes at a high temperature to provide sufficient fitting data for MTP [green dots in Fig. 1(a)]. We choose a cutoff radius of $5$ Å (including the first up to the third neighbor shell). Additional tests for a smaller cutoff including two shells show a small change sup .

Since machine-learning potentials have an inherently low extrapolation capacity, stability over the relevant part of the phase space is a critical issue. Detailed tests reveal that the MTP, fitted according to the above procedure, is sufficiently stable in the relevant volume and temperature range for the application within TU-TILD [see Fig. 1(a)]. In fact the potential can be also used at extrapolated volumes and temperatures and predicts even the onset of the liquid phase. Only a small number of MD runs in the range of a few % (see gray contour lines) becomes unstable. The results shown in Fig. 1(a) correspond to a “single-shot” MTP potential fitted to an initial set of DFT data. However, the MTP provides an inherent metric Podryabinkin and Shapeev (2016) to quantify the degree of extrapolation and thus offers the possibility to actively sample configurations for fitting, ensuring stability for the temperatures of interest.

The performance of the MTP as a reference potential within TU-TILD can be quantified by: (1) The dependence of $\langle E^{\rm DFT}-E^{\rm ref}\rangle_{\lambda}$ on $\lambda$ [Eq. (1); where $E^{\rm ref}$ stands for MTP energies] should be as smooth as possible, (2) the standard deviation of the energy difference $E^{\rm DFT}-E^{\rm ref}$ should be as small as possible, and (3) the correlation in the forces should be as strong as possible. The thermodynamic integration from the harmonic potential to MTP will not be discussed since, as mentioned above, this stage of the integration can be computed highly efficiently given the fact that the MTP is more than six orders of magnitude faster than DFT. See the supplement for detailed timings sup .

The excellent performance of the MTP is demonstrated in Figs. 2 and 3(a). The MTP energies are so close to the DFT energies that the thermal average $\langle E^{\rm DFT}-E^{\rm ref}\rangle_{\lambda}$ is almost independent of $\lambda$ and close to the targeted error of 1 meV [Fig. 2(b), black curve], i.e., the resulting MTP free energy is only 1 meV/atom away from the DFT free energy [Fig. 2(a)]. Computing this difference can be done highly efficiently because of the small standard deviation in the range of only 2 meV/atom [Fig. 2(c)]. Consistently, the MD forces predicted by the optimized MTP show a strong correlation with the DFT forces [Fig. 3(a)]. This good performance of the MTP is found for the whole relevant volume range sup demonstrating that an efficient study of the thermal expansion is possible as well.

To set a baseline for the performance of the MTP we employ alternative reference potentials that have been used previously for chemically less complex unary and binary systems: (1) 0 K harmonic Grabowski et al. (2009); Glensk et al. (2015); Grabowski et al. (2015), (2) effective harmonic Hellman et al. (2013), and (3) embedded atom method (EAM) Duff et al. (2015a); Gong et al. (2018). We start with the 0 K harmonic potential computed for VNbMoTaW in Ref. Körmann et al. (2017). As mentioned above, using this reference in Eq. (1) directly provides the anharmonic contribution.

Figure 2 highlights the difficulties. Figure 2(b) shows the very nonlinear dependence of the thermodynamic average, $\langle E^{\rm DFT}-E^{\rm ref}\rangle_{\lambda}$ (where $E^{\rm ref}$ stands now for the harmonic energies), on the coupling constant $\lambda$ (green curve). Due to this nonlinear behavior the evaluation requires many sampling points. Using the MTP reference introduced above makes a highly accurate calculation of this quantity possible. We find that at 3000 K the anharmonic contribution is $-31$ meV/atom. Due to their more open structure, bcc materials are known to have a more complex temperature dependence of the anharmonic free energy than close-packed materials Fultz (2010). Our calculations confirm this observation quantitatively: The anharmonic contribution is almost twice as large as for previously investigated, close-packed face-centered cubic (fcc) structures (range of 1 to 25 meV/atom at the melting point) Glensk et al. (2015).

An even more serious issue than the strongly nonlinear dependence of $\langle E^{\rm DFT}-E^{\rm ref}\rangle_{\lambda}$ is the fact that long MD runs are required to statistically converge each of these points. The underlying reason is that the standard deviation of the energy difference, $E^{\rm DFT}-E^{\rm ref}$, is large as shown in Fig. 2(c). The large difference in the energies is also reflected in the weak correlation between the harmonic and DFT forces during an MD run as shown in Fig. 3(d).

We now investigate an effective harmonic force constant matrix constructed at finite temperature as a reference. Such a matrix is employed, e.g., in the temperature-dependent effective potential (TDEP) method Hellman et al. (2013), and provides the advantage that analytical formulas can be used to compute the vibrational free energy. Our tests show that including pair interactions up to the first- and second-nearest neighbors gives similar results as with an additional third shell sup . Results for three shells will be discussed in the following. The interactions are determined from a least-square fit of the forces from more than 1500 configurations of an MD run at 3000 K at the target lattice constant. Owing to the harmonic approximation, the fitting problem is linear Freysoldt (2017) and is solved with a standard algebraic method, namely the pseudo-inverse from singular value decomposition to avoid accidental ill-conditioning. The zero-force reference structure, where the potential attains its minimum is set to the [math] K equilibrium geometry.

The forces from such an effective harmonic potential show a slightly better correlation with the DFT MD forces at the target temperature than the 0 K harmonic forces [cf. Figs. 3(c) vs. (d)]. Correspondingly the dependence of $\langle E^{\rm DFT}-E^{\rm ref}\rangle_{\lambda}$ on $\lambda$, where $E^{\rm ref}$ now stands for the effective harmonic energies, is less nonlinear than for the 0 K harmonic energies [Fig. 2(b), red vs. green]. However, the standard deviation is smaller only at $\lambda>0.5$. Thus an effective harmonic potential offers only a slightly improved reference for thermodynamic integration ($\approx$ 1.5 times faster). The respective vibrational free energy obtained using the standard harmonic formulas including a correction of the internal energy as introduced in the TDEP method Hellman et al. (2013) gives a slightly reduced error of 19 meV/atom compared to $-31$ meV/atom for the 0 K harmonic reference [Fig. 2(a)].

The still rather large error of the effective harmonic matrix is related to strong local pairwise anharmonicity Glensk et al. (2015). The inherent asymmetry of the nearest neighbor potential when atoms move together or apart cannot be captured in general by any harmonic potential irrespective of the temperatures it is fitted to. To take the required asymmetry properly into account an asymmetric potential parametrization is required. Asymmetric potentials are offered by the MTP discussed already above or likewise by an EAM parametrization.

We thus investigate whether an EAM fit for the complex disordered VNbMoTaW HEA is possible and, if it is, how it performs for thermodynamic integration. We employ the meamfit v2 package Duff et al. (2015b); mea . Our tests show that the number of expansion parameters has only a small influence sup . Results shown in the following refer to our best EAM with 3 embedding terms per species, with 11 parameters for the electron densities, and 19 for the pair-potentials. In total there are 355 independent parameters for this chemically highly complex quinary system, rendering the extraction of an accurate potential particularly challenging. To address this challenge we fit initially to a subset of the $\approx 8000$ energies available across all volumes. This subset consists of $\approx 2000$ uniformly-spaced energies, providing sufficient points per parameter to prevent over-fitting. We then take the best performing potential as a starting point for a single shot conjugate gradient fit to all $\approx 8000$ energies.

During the optimization, energies are computed relative to the 0 K relaxed structure for the corresponding volume. A cutoff radius of 5 Å (as for MTP) is imposed for the pairwise terms, and negative ‘electron’ densities are allowed—although positive background densities are required overall—to provide maximum variational flexibility. The resulting potential is stable across a wide range of volumes and temperatures [see Fig. 1(b)] and predicts the onset of the liquid phase.

The forces obtained with the fitted EAM potential show a better correlation with DFT MD forces [Fig. 3(b)] than the harmonic potentials. Consistently, the dependence of $\langle E^{\rm DFT}-E^{\rm ref}\rangle_{\lambda}$ on $\lambda$, where $E^{\rm ref}$ stands now for the EAM energies, is more linear and the standard deviation is smaller for all $\lambda$ (blue curves in Fig. 2). Using the EAM as a reference is more than three times faster than an effective harmonic potential. However, the EAM cannot compete with the MTP as a reference potential which further increases the efficiency by about an order of magnitude as compared to the EAM.

These results reveal that the combination of TU-TILD with MTP represents the presently most efficient combination to compute the vibrational free energy contribution. A preliminary study prelim shows that this does not only apply to equiatomic compositions such as the one studied in the present manuscript but likewise to arbitrary non-equiatomic compositions, and further also to different crystallographic lattice types such as hcp or fcc, magnetic systems, and even to the liquid phase.

The underlying physical reason for the excellent performance of the TU-TILD+MTP combination is the fact that the vibrational free energy is determined by a rather well-defined, sufficiently smooth, and local—although strictly anharmonic—part of the phase space. Several other studies Glensk et al. (2015); Zhang et al. (2017); SFE have already indicated that long-ranged interactions that maybe present at $T=0$ K due to quantum-mechanical interference effects, vanish when explicit vibrations are introduced at finite temperatures due to the breaking of the crystal symmetries. Effective interactions at finite temperatures are thus localized and can be well fitted by a local approach. These interactions are strongly anharmonic requiring an anharmonic description as provided by the MTP or EAM, with a greater flexibility offered by the MTP foot . It should be stressed that this greater flexibility comes with a lower extrapolation capability (see Fig. 1 and corresponding discussion). A main achievement of the present work is having shown that, for free energy calculations within the TU-TILD+MTP approach, the stability of the MTP is sufficient and that providing a set of well distributed fitting data renders the sampling of the thermally accessible phase space an interpolation task; a task optimally suited for a machine learning approach.

To open the approach to a broad community we are presently implementing it into the pyiron environment Janssen et al. (2018); py . This, together with the performance of the new approach, paves the route to compute vibrational free energies not only highly accurately but with a computational performance adequate for high-throughput screening of multicomponent alloys.

Acknowledgments We thank Jan Janssen and Konstantin Gubaev for fruitful discussions. Funding by the Deutsche Forschungsgemeinschaft (SPP 2006) and the European Research Council (ERC) under the EU’s Horizon 2020 Research and Innovation Programme (Grant No. 639211) is gratefully acknowledged. F.K. was supported by the NWO/STW (VIDI grant 15707). A.S. was supported by the Russian Science Foundation (Grant No 18-13-00479). A.I.D. acknowledges support from the STFC Hartree Center program, Innovation: Return on Research, funded by the UK Department for Business, Energy & Industrial Strategy.

Competing interests The Authors declare no Competing Financial or Non-Financial Interests.

Data availability The authors declare that all data supporting the findings of this study are available within the paper and its supplementary information files.

Contributions B.G. performed the TU-TILD approach, C.F. the implementation of the effective harmonic force constant fitting procedure, A.I.D. the MEAM-fits, A.S. developed and fitted the MTPs. All authors designed the project, discussed the results, and wrote the manuscript.