Ab initio vibrational free energies including anharmonicity for multicomponent alloys

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.

Recent developments in the field of multicomponent alloys [high entropy alloys (HEAs), compositionally complex alloys (CCAs)] have opened new materials design perspectives [1][2][3][4]. 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 [5]. However, the requirement [6][7][8] to accurately capture small free energy differences (≈ 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 [8][9][10], 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 [11,12].
The free energy is determined by various mechanisms including electronic excitations, configurational entropy and atomic vibrations (e.g., Refs. [5,13]). The major contribution is due to atomic vibrations (e.g., for VNbMoTaW [14]: 8.6 k B at 1500 K vs. a configurational entropy of 1.6 k 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 temperaturedependent phonon softening and broadening. Effective harmonic Hamiltonians [15][16][17][18][19][20] can approximately account for the temperature induced changes. Numerically exact vibrational free energies can be obtained by thermodynamic integration [21][22][23][24][25], from a reference potential E ref with free energy F ref to DFT energies E DFT , where · · · λ denotes a thermal average on a mixed potential Using a harmonic reference would in principle give the exact anharmonic free energy including temperaturedependent 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 λ-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 [9], 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 stagecontaining 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. . The gray lines in (a) are contours at which 1 or 5% of the MTP MD runs are unstable [26], however, all the runs can be made stable by using active learning [27].
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 [28]. Several machine-learning potentials have been proposed so far [29][30][31][32][33][34]. 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 [35]. First attempts have been put forward to describe alloys [36], 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. [37], and recently shown to perform best among different machine-learning models [38]. 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 ith atom by the moments of inertia of the neighboring atoms, Here the radial functions f n,i,j (r ij ), n = 1, 2, . . ., define different shells around the ith atom; the contribution of the jth atom to the nth shell can depend on the types of the ith and jth atoms. When m = 0, M n,0 is a scalar quantity interpreted as the weight of the nth shell. The set of these scalar descriptors is not complete. However, this set can be made complete by adding vectorial "eccentricity" of the nth 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 [36,37]. 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 [39]. DFT calculations are performed with vasp [40,41] employing PAW [42], GGA-PBE [43], and including electron-phonon coupling [44]. See the supplement for details [26]. We consider temperatures up to 3000 K which is close to the estimated melting point [12]. In accord with our previous works [8,45] 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 [26].
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 [27] 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) (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 [26].
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 E DFT −E ref λ is almost independent of λ 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 [26] 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 [10,46,47], (2) effective harmonic [17], and (3) embedded atom method (EAM) [9,45]. We start with the 0 K harmonic potential computed for VNbMoTaW in Ref. [14]. 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, E DFT − E ref λ (where E ref stands now for the harmonic energies), on the coupling constant λ (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 [48]. Our calculations confirm this observation quantitatively: The full DFT free energy 0 K harm. 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) [10]. An even more serious issue than the strongly nonlinear dependence of E DFT − E ref λ 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 DFT − E ref , is large as shown in Fig. 2(c). The large difference in the energies is also reflected in the weak correlation between the har- monic 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 [17], 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 [26]. 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 [49] 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 0 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 E DFT − E ref λ on λ, where E 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 λ > 0.5. Thus an effective harmonic potential offers only a slightly improved reference for thermodynamic integration (≈ 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 [17] 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 [10]. 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 [50,51]. Our tests show that the number of expansion parameters has only a small influence [26]. 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 pairpotentials. 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 ≈ 8000 energies available across all volumes. This subset consists of ≈ 2000 uniformly-spaced energies, providing sufficient points per parameter to prevent overfitting. We then take the best performing potential as a starting point for a single shot conjugate gradient fit to all ≈ 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 E DFT − E ref λ on λ, where E ref stands now for the EAM energies, is more linear and the standard deviation is smaller for all λ (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 [52] 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 [10,44,53] have already indicated that longranged 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 [54]. 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 [55,56]. 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.