Introduction

Recent developments in the field of multicomponent alloys (high entropy alloys (HEAs) and compositionally complex alloys (CCAs)) have opened new materials design perspectives.1,2,3,4 The prediction and exploration of thermodynamic properties and phase stabilities are 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 requirement6,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 and accurate approach to compute the vibrational contribution to the free energy 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 different contributions such as atomic vibrations, electronic excitations, or chemical configurations (e.g., refs. 5,13). For a fixed atomic configuration, e.g., a given chemically ordered or disordered atomic arrangements of atoms, a main contribution is due to atomic vibrations, 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 Hamiltonians14,15,16,17,18,19 can approximately account for the temperature-induced changes. Numerically exact vibrational free energies can be obtained by thermodynamic integration,20,21,22,23,24

$$F = F^{{\mathrm{ref}}} + \mathop {\int}\limits_0^1 d \lambda \left\langle {E^{{\mathrm{DFT}}} - E^{{\mathrm{ref}}}} \right\rangle _\lambda ,$$
(1)

from a reference potential Eref with free energy Fref to DFT energies EDFT, where 〈⋯〉λ denotes a thermal average on a mixed potential Eλ = λEDFT + (1 − λ)Eref. 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 λ-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 stage—containing the explicit and 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.25 Several machine-learning potentials have been proposed so far.26,27,28,29,30,31 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.32 First attempts have been put forward to describe alloys,33 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 (MTPs), a class of machine-learning potentials first proposed in ref. 34, and recently shown to perform best among different machine-learning models.35 MTP describes the atomic environment of the ith atom by the moments of inertia of the neighboring atoms,

$$M_{n,m} = \mathop {\sum}\limits_j {f_{n,i,j}} (r_{ij})\underbrace {{\boldsymbol{r}}_{ij} \otimes {\boldsymbol{r}}_{ij} \otimes \ldots \otimes {\boldsymbol{r}}_{ij}}_{m\,{\mathrm{times}}}.$$

Here the radial functions fn,i,j(rij), 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 atom. When m = 0, Mn,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 Mn,1, the tensor of second moments of inertia Mn,2, the third moments Mn,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.33,34 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 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. In particular, we apply an MTP as a reference potential within TU-TILD for the chemically complex disordered VNbMoTaW HEA and show that it is clearly superior to alternative reference potentials.

Results and discussion

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 procedure described in section “Methods”, is sufficiently stable in the relevant volume and temperature range for the application within TU-TILD (see Fig. 1a). 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 percent (see gray contour lines) becomes unstable. The results shown in Fig. 1a correspond to a “single-shot” MTP potential fitted to an initial set of DFT data. However, the MTP provides an inherent metric36 to quantify the degree of extrapolation and thus offers the possibility to actively sample configurations for fitting (as, e.g., employed in ref. 37), ensuring stability for the temperatures of interest.

Fig. 1
figure 1

Internal energy surfaces, U, as a function of lattice constant/volume and temperature, T, for a MTP and b EAM. Green dots mark the volumes (at 3000 K) used for fitting the potentials. Black dotted lines emphasize the transition to the liquid phase. The squared volume and temperature region is the region of interest (3.20–3.28 Å in terms of the lattice constant). The gray lines in a are contours at which 1 or 5% of the MTP MD runs are unstable (see Supplementary Information), however, all the runs can be made stable by using active learning36

Full size image

The performance of the MTP as a reference potential within TU-TILD can be quantified by the following: (1) the dependence of 〈EDFT−Eref〉λ on λ (Eq. (1); where Eref stands for MTP energies) should be as smooth as possible, (2) the standard deviation of the energy difference EDFT−Eref 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 previously, 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 Supplementary Information for detailed timings.

The excellent performance of the MTP is demonstrated in Figs. 2 and 3a. The MTP energies are so close to the DFT energies that the thermal average 〈EDFT−Eref〉λ is almost independent of λ and close to the targeted error of 1 meV (Fig. 2b, black curve), i.e., the resulting MTP free energy is only 1 meV/atom away from the DFT free energy (Fig. 2a). Computing this difference can be done highly efficiently because of the small standard deviation in the range of only 2 meV/atom (Fig. 2c). Consistently, the MD forces predicted by the optimized MTP show a strong correlation with the DFT forces (Fig. 3a). This good performance of the MTP is found for the whole relevant volume range (see Supplementary Information) demonstrating that an efficient study of the thermal expansion is possible as well.

Fig. 2
figure 2

Results of thermodynamic integration to DFT for VNbMoTaW using different references at 3000 K. The integral over the curves in b gives the difference in free energy between DFT and the reference shown in a. The smaller the standard deviation shown in c, the quicker the statistical convergence of the curves in b as indicated on the right axis in c. For the CPU time calculation a standard error of 1 meV/atom, a CPU time of 4 h per ionic step, and a correlation length of 15 steps were taken

Full size image
Fig. 3
figure 3

Correlation of the DFT forces versus forces from the different approximations at 3000 K. The color represents the local density. The numbers in the right lower corners represent the root mean square error of the distributions

Full size image

We have investigated the MTP-based approach also for different SQS permutations to simulate varying chemical environments for the VNbMoTaW system. We find that the total vibrational free energy varies only slightly between the four investigated supercells (standard deviation of 1.1 meV/atom). This indicates that the employed supercell size provides a converged sampling of the ideal disordered limit. Further we find a similarly good computational performance and stability for the individual MTPs for each SQS. Additional tests also show that a single MTP with a similar performance can be obtained by fitting to all SQS’ simultaneously.

To set a baseline for the performance of the MTP as compared to alternative reference potentials that have been used previously for chemically less complex unary and binary systems we study: (1) 0 K harmonic,10,38,39 (2) effective harmonic,16 and (3) embedded atom method (EAM).9,40 We start with the 0 K harmonic potential computed for VNbMoTaW in ref. 41. As mentioned in “Introduction”, using this reference in Eq. (1) directly provides the anharmonic contribution.

Figure 2 highlights the difficulties. Figure 2b shows the very nonlinear dependence of the thermodynamic average, 〈EDFT−Eref〉λ (where Eref 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.42 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–25 meV/atom at the melting point).10

An even more serious issue than the strongly nonlinear dependence of 〈EDFT−Eref〉λ for the harmonic reference 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, EDFT−Eref, is large as shown in Fig. 2c. 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. 3d.

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,16 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 (see Supplementary Information). 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 linear43 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. Fig. 3c vs. d). Correspondingly the dependence of 〈EDFT−Eref〉λ on λ, where Eref now stands for the effective harmonic energies, is less nonlinear than for the 0 K harmonic energies (Fig. 2b, 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 method16 gives a slightly reduced error of 19 meV/atom compared to −31 meV/atom for the 0 K harmonic reference (Fig. 2a).

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.44,45 Our tests show that the number of expansion parameters has only a small influence (see Supplementary Information). 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 ≈8000 energies available across all volumes. This subset consists of ≈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 ≈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. 1b) 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. 3b) than the harmonic potentials. Consistently, the dependence of 〈EDFT−Eref〉λ on λ, where Eref 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.

An attempt has been made to optimize a reference-free modified EAM (RF-MEAM) potential, however, due to the size of the potential parameter space we were unable to obtain an RF-MEAM potential, which improved on the EAM potential. It is worth noting that this is an area of active research, with recent improvements by one of us (A.I.D.) to the underlying MEAMfit algorithm as well as, e.g., consideration of preconverged binary and ternary potentials as starting points, likely to render such optimizations feasible in the near future.

Overall, the results of the present study reveal that the combination of TU-TILD with MTP represents the presently most efficient combination to compute the vibrational free energy contribution of chemically complex alloys. Ongoing investigations46 indicate that this applies not only to equiatomic compositions such as the one studied in the present study, but likewise to arbitrary nonequiatomic compositions, and further also to different crystallographic lattice types such as hcp or fcc 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 studies10,47,48 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. 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 put the here proposed method into a practical perspective, we foresee a number of potential applications such as, e.g., computing transition temperatures between ordered and disordered phases as recently performed within the quasiharmonic approximation for the V-Mo-Nb-Ta-W system by Wang et al.49 or by coupling it to the recently developed low-ranked potential method,50 which has been successfully applied to refractory multicomponent alloys.51 Other potential applications are melting temperature calculations based on the method developed in ref. 8 or accurate stacking-fault energy calculations.48 Of interest would be an extension to magnetic materials, in particular for efficiently treating vibrations in the paramagnetic state.52,53,54 The key step towards such an extension would be the development of efficient machine-trained magnetic potentials, which capture the magnetic degree of freedom. Despite recent progress employing Gaussian approximation potentials on purely ferromagnetic iron,55 the explicit inclusion of the magnetic degree of freedom has yet to be achieved in one of the current machine-learning potential frameworks.

To open the approach to a broad community, we are presently implementing it into the pyiron environment (http://pyiron.org).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.

Methods

Chemical disorder is modeled by a special quasirandom structure (SQS) in a 125 atomic supercell.57 DFT calculations are performed with VASP58,59 employing PAW,60 and GGA-PBE.61 Tests for different chemical permutations of the SQS revealed similar results. The impact of electronic excitations on the interatomic interactions has been taken into account by employing finite-temperature DFT as developed by Mermin62 (see also, e.g., ref. 47). For further details we refer to the Supplementary Information. We consider temperatures up to 3000 K, which is close to the estimated melting point.12 In accord with our previous works,8,40 we use DFT MD simulations at several volumes at a high temperature to provide sufficient fitting data for MTP (green dots in Fig. 1a). 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 (see Supplementalary Information).