Electronic Moment Tensor Potentials include both electronic and vibrational degrees of freedom

Abstract

We present the electronic moment tensor potentials (eMTPs), a class of machine-learning interatomic models and a generalization of the classical MTPs, reproducing both the electronic and vibrational degrees of freedom, up to the accuracy of ab initio calculations. Following the original polynomial interpolation idea of the MTPs, the eMTPs are defined as polynomials of vibrational and electronic degrees of freedom, corrected to have a finite interatomic cutoff. Practically, an eMTP is constructed from the classical MTPs fitted to a training set, whose energies and forces are calculated with electronic temperatures corresponding to the Chebyshev nodes on a given temperature interval. The eMTP energy is hence a Chebyshev interpolation of the classical MTPs. Using the eMTP, one can obtain the temperature-dependent vibrational free energy including anharmonicity coming from phonon interactions, the electronic free energy coming from electron interactions, and the coupling of atomic vibrations and electronic excitations. Each of the contributions can be accessed individually using the proposed formalism. The performance of eMTPs is demonstrated for two refractory systems which have a significant electronic, vibrational and coupling contribution up to the melting point—unary Nb, and a disordered TaVCrW high-entropy alloy. Highly accurate thermodynamic and kinetic quantities can now be obtained just by using eMTPs, without any further ab initio calculations. The proposed construction to include the electronic degree of freedom can also be applied to other machine-learning models.

Introduction

Ab initio modeling and prediction of thermodynamic properties of materials, especially at high temperatures, necessitates accurate prediction of free energies including all relevant excitation mechanisms¹. The vibrational free energy coming from thermal excitations of atomic motion is the major contributor². Nevertheless, other contributions such as the electronic free energy coming from the excitation of electrons with temperature can also be significant and, if neglected, can lead to falsely predicted properties. This is especially true in metallic materials ranging from simple unary systems to more complex multi-component alloys^3,4. Additionally, the electronic free energy gets affected considerably by thermal vibrations of atoms, which also eventually alters the predicted thermodynamic properties^5,6.

The most common ab initio method to calculate the free energy is density functional theory (DFT). However, it is expensive and restricted to small systems and short time scales. The vibrational free energy can be calculated using ab initio molecular dynamics (AIMD) simulations. But, without additional sampling techniques, these require several tens of thousands of steps for convergence⁷. To calculate the electronic free energy, one could use approximate methods such as the Sommerfeld approximation⁸ which depends on the effective, temperature-independent density of states (DOS) at the Fermi level. But, this approximation is far too simplified and leads to severe inaccuracies when the DOS gets affected by thermal fluctuations. The electronic free energy thus needs to be calculated from self-consistent calculations within the finite temperature DFT framework. Since the atomic vibrations have an effect on the electronic free energy, such calculations need to be performed on multiple configurations at different volumes and temperatures for an accurate and converged representation of the total free energy surface, making them computationally expensive.

Progress in machine-learning (ML)-based interatomic potential development^{9,10,11,12,13,14,15} has accelerated materials simulations and thermodynamic property predictions. The ML models represent the energy of a system as a function of the atomic positions. Hence, a well-trained model, by virtue of its construction, can accurately capture the vibrational degrees of freedom. Progress has also been made to include the effect of electronic temperature in an ML-based framework, such as within the electronic temperature dependent deep potential molecular dynamics¹⁶. More recently, data extracted from the local electronic density of states from finite-temperature DFT calculations¹⁷, and exclusive ground-state DFT calculations on configurations generated at 0 K and finite-temperature^13,15, have been incorporated into neural-network (NN)-based potentials to estimate total free energies. These ML models are parametrized using a non-linear NN regression scheme.

A complementary class of ML potentials is based on the linear expansion of the site energy with respect to non-linear tensor-based descriptors. One such ML-model is the moment tensor potential (MTP)¹⁸ whose efficiency and accuracy has been shown for systems ranging from unaries¹ to complex high-entropy alloys (HEAs)^4,19,20. MTPs are typically built with less fitting parameters and can thus be optimized with fewer input data in comparison to the non-linear layers in NN models. The performance of the MTPs has prompted their use in high-temperature thermodynamic property predictions up to melting point. Recently, Jung et al.¹ proposed an efficient methodology—direct upsampling—to accurately predict free energies of various systems. The methodology used an MTP as an intermediate model in a thermodynamic integration (TI) scheme. The MTP was able to predict accurate vibrational free energies, including anharmonic contributions, up to the melting point. However, an ‘upsampling’ step to DFT was necessary to calculate the electronic free energy, including the coupling coming from vibrations. This stage of the methodology was computationally the most expensive, since self-consistent DFT calculations were performed on several tens of configurations at each point on the volume-temperature grid.

In the present work, we introduce a new single integrated formalism of the MTP—the electronic moment-tensor potential (eMTP)—which, in addition to capturing the vibrational free energy, is also able to predict accurate temperature-dependent electronic free energies. Analogous to the magnetic moment-tensor potential (mMTP)²¹ which includes spin degrees of freedom, here the eMTP has an additional degree of freedom to capture the electronic part. With the current formalism, the full free energy surface up to ab initio accuracy can be obtained just by using eMTPs, thereby tremendously increasing the computational efficiency of the calculations. We train eMTPs and obtain the corresponding thermodynamic properties up to melting point (T_melt) for two prototype systems—unary Nb, and a disordered four-component TaVCrW high-entropy alloy (HEA). The results are compared to the corresponding classical MTP predictions and full DFT properties that were recently calculated for Nb¹ and TaVCrW⁴. We compensate for the remaining small inaccuracies in the final thermodynamic properties coming from the eMTP training errors by upsampling to already available DFT training data, without performing any further DFT calculations. With the new formalism, in addition to accurate thermodynamic property calculations, the eMTPs can also be used to perform large-size and long-time molecular dynamics (MD) simulations to model accurate kinetic properties, for e.g., phase transformations in materials.

Results

Electronic moment tensor potential

The objective of the eMTP formalism is to accurately predict total energies of atomic configurations, including the contribution from electronic excitations. The energies and forces of the training set come from DFT calculations. Since the electronic contribution needs to be included, we utilize finite-temperature DFT calculations which are practically implemented via the Fermi broadening approach²². The DFT energies are represented as ${E}_{{{{\rm{DFT}}}}}({{{\boldsymbol{{{{\mathcal{R}}}}}}}},{T}_{{{{\rm{el}}}}})$, where ${{{\boldsymbol{{{{\mathcal{R}}}}}}}}=\{{{{{\bf{R}}}}}_{I}\}$ is the set of atomic positions R_I, and T_el corresponds to the electronic temperature. The configurations in the training set include electronic excitations and fully anharmonic vibrations and the corresponding coupling effects. Anharmonic vibrations correspond to phonon-phonon interactions, which become very strong and indispensable^6,23 at high temperatures. The eMTP formalism captures all these relevant high-temperature contributions.

The classical MTPs represent the energy E_MTP only as a polynomial of atomic positions ${{{\boldsymbol{{{{\mathcal{R}}}}}}}}$, corrected so that the local atomic energies do not depend on atoms beyond a certain cutoff. The general concept of the eMTP is to construct a polynomial of both atomic positions ${{{\boldsymbol{{{{\mathcal{R}}}}}}}}$ and electronic temperature T_el via classical MTPs. Because it is a full polynomial of both ${{{\boldsymbol{{{{\mathcal{R}}}}}}}}$ and T_el, the eMTP is capable of describing an arbitrary coupling between the two families of degrees of freedom, and accurately predicting ${E}_{{{{\rm{DFT}}}}}({{{\boldsymbol{{{{\mathcal{R}}}}}}}},{T}_{{{{\rm{el}}}}})$. The proposed eMTP represents the energy E_eMTP as an interpolation polynomial of ${{{\boldsymbol{{{{\mathcal{R}}}}}}}}$ and T_el, and is conveniently expressed as a temperature-dependent linear combination of classical MTPs for different T_el’s fitted to training sets each computed with a fixed electronic temperature. Specifically, a set of six electronic temperatures properly distributed across the relevant temperature range is chosen for each system. The distribution of the temperatures is dictated by a combination of Chebyshev nodes of different degrees between 0 K to T_melt. A more detailed description is provided in the Methods section.

Fitting results

In Table 1, we show the fitting results for the set of MTPs that make up the eMTP formalism, for both Nb and TaVCrW. For each system, each MTP is fitted to the same set of configurations but with the energies and forces corresponding to different electronic temperatures. The set of configurations is generated at different MD temperatures and volumes. For both Nb and TaVCrW, the MTPs are trained to within 2–3 meV/atom and 0.07–0.11 eV/Å root-mean-square error (RMSE) in energies and forces, respectively. The errors on the training set achieved here are in the same range as conventional training errors of classical MTPs in earlier works^1,4,23. However, in these earlier works, one classical MTP was fit to pure vibrational effects (no Fermi broadening in the DFT training/test set) at a high temperature (melting point). This single MTP showed good performance over a large temperature range as long as no electronic temperature was included. The small errors of such a single classical MTP are no longer possible across a larger temperature range while simultaneously also including the electronic degree of freedom, and one needs to resort to the eMTP formalism.

Table 1 Energy and force root-mean-square errors (RMSEs) for the set of MTPs at different electronic temperatures that make-up the eMTP formalism for Nb and TaVCrW

Figure 1 demonstrates the performance and advantage of the eMTP formalism for both Nb (upper row) and TaVCrW (bottom row). In Fig. 1a, c, we plot the RMSE in energy predicted by the corresponding eMTP model with respect to DFT at five different temperatures. These temperatures fall in the range of 0 K to T_melt, but are not exactly the temperatures used in the training (seen in Table 1). At each temperature, we have three sets of volumes. The results are compared to the predictions of a single classical MTP, the MTP from Jung et al.¹ for Nb, and the MTP from Zhou et al.⁴ for TaVCrW.

Fig. 1: Performance of the eMTP.

The RMSE in the predicted energies with respect to DFT for a validation set for (a) Nb and (c) TaVCrW at five different temperatures, for the eMTP formalism (green) and previously fitted classical MTPs^1,4 (red). The DFT energies also include electronic effects. Energy predicted by the eMTP formalism and the classical MTP with respect to DFT (including electronic effects) in (b) Nb and (d) TaVCrW, on configurations generated at 2500 K. The ‘outlier’ dots at the beginning of each volume correspond to the ideal lattice structure at each volume.

For both Nb and TaVCrW, the overall performance of the eMTP and how it compares to the corresponding classical MTP is similar. In both cases, the eMTP formalism predicts energies to within 2–3 meV/atom accuracy at all temperatures, as seen by the green dots in Fig. 1a, c, signifying its outstanding performance. The eMTP formalism is thus able to achieve the same accuracy on the test data across the entire temperature range as the training errors of the individual MTPs that it comprises of. Contrarily, the RMSE of the single, classical MTP (red dots in Fig. 1a, c) increases considerably with temperature in both systems. The value reaches to 89 meV/atom in Nb and 93 meV/atom in TaVCrW at 2500 K. The classical MTPs were fitted at the melting point for Nb, and at 2500 K for TaVCrW, to DFT energies without the electronic contribution^1,4. Hence, although the classical MTPs predict accurate DFT energies without the electronic contribution (as shown previously^1,4), they cannot capture the temperature-dependent electronic contribution. The RMSEs of the classical MTPs in Fig. 1a, c are in fact, a few meV/atom away from the electronic free energy corresponding to that particular temperature (cf. Zhang et al.⁵ for Nb, and Zhou et al.⁴ for TaVCrW). Obviously, such high values, that significantly increase with increasing temperature, are non-negligible and have a pronounced effect on thermodynamic property predictions, as will be discussed later in the article.

Figure 1b, d further emphasize the eMTP performance by focusing only on configurations at 2500 K from three different volumes for Nb and TaVCrW. In both cases, the energies predicted by the eMTP formalism, denoted by green dots, are within a few meV/atom away from DFT. The distribution of errors has a standard deviation of 1.3 meV/atom and 2.1 meV/atom for Nb and TaVCrW, respectively. On the other hand, the energies predicted by the classical MTP, denoted by the red dots, are around 90 meV/atom higher than the DFT values for both systems, which, as mentioned above, is approximately the magnitude of the electronic free energy at 2500 K for Nb⁵ and TaVCrW⁴. The standard deviation of the classical MTP prediction is also larger, around 5.1 meV/atom and 5.9 meV/atom for Nb and TaVCrW. Although the energies are not predicted correctly by the classical MTP, the force RMSEs are similar to those of the eMTP and fall within the range of values shown in Table 1. This arises due to the fact that the electronic temperature has a negligible effect on the atomic forces in the present systems, as also observed earlier in refractory unaries⁶. Hence, a classical MTP fitted to a single temperature (melting point) without any electronic contribution is still able to reasonably predict atomic forces across the entire temperature range. This aspect will be further discussed and validated in the ‘Discussion’ section.

Thermodynamic properties

Given the predictive capability of the eMTP formalism, we use it to calculate the full free energy surface and numerically derive thermodynamic properties for Nb and TaVCrW. We focus on three thermodynamic properties—the isobaric heat capacity (C_p), the lattice expansion coefficient (α), and the adiabatic bulk modulus (B_S)—which are highly sensitive to the free-energy surface. The full free-energy surface is obtained by parametrizing total free energies calculated on a dense volume, temperature (V, T) grid. At each (V, T), we perform thermodynamic integration using Langevin dynamics (TILD). For the TILD, we use a low-temperature effective quasiharmonic (QH) model^24,25 as the reference, and perform TILD to the eMTP model. The TILD calculation parameters, the density of the grid, and the polynomials used for the free energy surface parametrization are described in the Methods section. The computational details are consistent with those used in literature for Nb¹, and for TaVCrW⁴.

Figure 2 shows the properties for Nb (top row) and TaVCrW (bottom row) computed up to the level of accuracy of the eMTPs as solid green curves. The results are compared to those calculated using classical MTPs (red dashed curves) and full high-accuracy ab initio values (blue curves). The classical MTPs and the full DFT curves are taken from previous literature^1,4. Additionally, for Nb, we also include experimental data points from literature, denoted by blue circles. No high-temperature experimental data is available yet for the TaVCrW HEA.

Fig. 2: Thermodynamic property predictions.

Isobaric heat capacity, lattice expansion coefficient, and adiabatic bulk modulus of (a) Nb up to 2750K, and (b) TaVCrW up to 2500 K, calculated using the classical MTP (red), using the eMTP (green solid) followed by an upsampling to the DFT training data (green dashed) and up to full DFT accuracy (blue solid). The classical MTPs used for the red curves and the full DFT values are taken from Jung et al.¹ and Zhou et al.⁴ for Nb and TaVCrW, respectively. For Nb, experimental values from literature^40,41,42 are shown as light blue circles for comparison.

The capability of the eMTP formalism is primarily highlighted in the C_p prediction, where the eMTP-predicted curve falls almost on top of the DFT curve for both systems. The classical-MTP predicted C_p curve falls much below the DFT values, by as much as 1–1.5 k_B at the melting point, since the classical MTP does not account for the large electronic contributions that are known to significantly affect the C_p values in refractory systems^1,4,6. There are deviations in the eMTP-predicted α and B_S with respect to the full DFT values, especially near the melting point. In both systems, α is marginally underestimated by the eMTP, noticed by the green solid curves falling below the blue line. On the other hand, the bulk modulus of TaVCrW is over-estimated by the eMTP almost consistently across the entire temperature range, including at 0 K where it is higher by around 5 GPa in comparison to the DFT value. Such an over-estimation is also seen in the classical MTP prediction. Nonetheless, the overall agreement of the eMTP predictions to full DFT values is good, and better than the classical MTP predictions. It is significant to notice that the eMTP predictions also fall very close to the experimental values.

Upsampling to DFT training data

The thermodynamic properties are second derivatives of the free energy (see Methods), and extremely sensitive even to very small differences. The deviations in the estimated properties between the eMTP and the DFT (differences in the green solid and blue solid curves in Fig. 2) that were discussed above are hence a direct consequence of the RMSEs of the eMTP formalism in predicted energies and forces. To account for the difference in the free energies coming from the RMSEs, we include an ‘upsampled’ contribution to the free-energy surface. We do so by taking advantage of the existing training data that was used to build the eMTP, which includes configurations generated at a set of MD temperatures and volumes. At each of these volume-temperature points, we calculate the free-energy difference predicted by the eMTP and the DFT value (which also includes the electronic contribution) using the free-energy perturbation formula, followed by a smooth and converged parametrization. The DFT value is already available from the training data. The ‘upsampled’ surface is then added to the eMTP free-energy surface before deriving the thermodynamic properties. The procedure is described in more detail in the Methods section.

Figure 2 also shows the thermodynamic properties obtained after the eMTP upsampling, denoted as ‘eMTP^up’ by green dashed lines. The eMTP-upsampled curves fall almost on top of the full DFT curves (blue lines). This is particularly evident in α and B_S, wherein, the deviations that were observed between the eMTP and DFT become considerably smaller, all the way up to the melting point. Even the 0 K bulk modulus of TaVCrW gets rectified once we include the upsampled contribution. The marginal differences that still remain (between the green dashed and blue curves) either arise from the difference in parameters between the high-accuracy DFT and those used for the training data, or from statistical deviations in the free-energy surface parametrization. Nonetheless, we are able to reach up to full ab initio accuracy by using only the eMTP and existing training data, and without any other extensive upsampling calculations. We would like to point out that the eMTP-upsampling can only be used to correct free-energy surfaces, and the resulting thermodynamic properties, and cannot be employed while using the eMTP to study kinetic processes, e.g., phase transformations.

Different contributions using the eMTP

An added advantage of the eMTP formalism is the mutually independent calculation of the vibrational and electronic part of the total free energy. For instance, the electronic temperature T_el can be turned off (set to 0 K) in order to model purely the effect of vibrations. The effect of the electronic temperature can then be included by performing an additional set of TILD runs from the eMTP with T_el = 0 K to the eMTP with T_el = T. One can afford to do such separate calculations on large systems owing to the relative computational efficiency in comparison to performing DFT calculations. Separately calculating the vibrational (anharmonic) and the electronic contribution at each (V, T) also makes it more convenient to parametrize highly accurate individual free energy surfaces, which are later added together for thermodynamic property predictions.

Figure 3 illustrates the individual effect of atomic vibrations and electronic excitations on the predicted thermodynamic properties (C_p, α, and B_S) of Nb, calculated using the eMTP formalism. The individual curves are shown here only for Nb, since the contributions in the TaVCrW system are very similar and can also be straightforwardly obtained with the methodology. The orange curves are the thermodynamic properties obtained using the effective QH reference. The effective QH model contains the harmonic and the expansion part of the vibrational free energy, but does not include anharmonic contributions, which are significant in bcc refractory systems²³. The shift in curves due to anharmonicity is demonstrated in Fig. 3 by the pink shaded regions. The calculated set of properties including the anharmonic contribution are denoted by magenta lines.

Fig. 3: Free-energy contributions to the thermodynamic properties.

Isobaric heat capacity, lattice expansion coefficient, and adiabatic bulk modulus of Nb up to 2750K calculated using the eMTP with different contributions to the free energy. The insets show exclusively the effect of the electronic free energy including the coupling to explicit vibrations on the thermodynamic properties, as predicted by the eMTP and DFT¹.

By turning on the electronic temperature in the next stage of TILD, the eMTP formalism calculates the effect of the electronic contribution on the thermodynamic properties. The electronic free energy predicted by the eMTP at this stage includes the coupling to atomic vibrations. The resulting change in properties is denoted by black-dashed regions in Fig. 3. Evidently, the electronic free energy plays a crucial role in the final thermodynamic properties. In the insets in each of the thermodynamic property figures, we compare the change in the respective property coming purely from the electronic part, as predicted by the eMTP formalism and DFT. The curves are extremely similar, further validating the calculation of the electronic part of the free energy using the eMTP formalism. The insets also reveal the significant effect and necessity to accurately include the electronic contribution.

Speed-up

Table 2 illustrates the speed-up achieved while using the eMTP formalism. We compare the time taken to predict thermodynamic properties of Nb, between the eMTP formalism and the recently developed direct upsampling methodology¹. The free energy calculations are performed on a 11 × 13 (V, T) grid. The table is split into two parts—one for the time taken for fitting, and the other for the free energy calculations.

Table 2 CPU-time estimate (in core hours) for different stages of the free energy calculation for Nb (PBE, 11 valence electrons) using the previous direct upsampling methodology and the current eMTP formalism

The eMTP formalism takes three times more time for the fitting in comparison to direct upsampling. In direct upsampling, a single MTP was used across the entire temperature range, whose training set involved 12 volumes at the melting point. In the eMTP formalism, we have six MTPs corresponding to six different electronic temperatures, each of which have configurations from five different volumes.

Although fitting the MTPs for the eMTP formalism takes more time, this loss in computational time is significantly smaller compared to the gain during the free energy calculations, as shown in the second part of the table. There are a couple of noteworthy points. Langevin dynamics calculations using the eMTP formalism take around five to six times more time than using a single MTP. This is noticed in the higher computational time required for calculating F^ah. Nonetheless, the most time-consuming part of the direct upsampling methodology—the upsampling—which involves DFT calculations, is now substituted instead with another round of TILD using the eMTP formalism. This leads to an increase of ×3.3 in speed in the free energy calculations. It is important to note that the second TILD using the eMTP is performed on large supercells, which is conventionally three to four times the size of the supercells used for the DFT calculations in the upsampling. By using a large supercell at this stage, the effect of long-wavelength phonons and their mutual interactions are included, albeit, only to the accuracy of the eMTP. Note also that the eMTP-upsampling to DFT that we have discussed above and demonstrated in Fig. 2 (green dashed curves) is performed on existing eMTP training data, and hence does not involve any additional computational cost.

Overall, summing up both the stages, we achieve close to a ×3 speed-up by using the eMTP formalism in highly accurate thermodynamic property predictions. Similar improvements in speed are also noticed for TaVCrW. The advantages using the newly proposed eMTP formalism are independent of the chemical complexity of the system that is studied.

Discussion

We have developed the eMTP formalism—a generalization of the classical MTPs—to calculate the total free energy of a system by taking into account both the vibrational and electronic degrees of freedom. The temperature-dependent electronic free energies are calculated by a linear combination of classical MTPs fitted at different electronic temperatures. The eMTP formalism has been applied to two refractory systems, Nb and TaVCrW, whose DFT-based thermodynamic properties are accurately reproduced in a much shorter amount of time, with no further ab initio calculations. The remaining small inaccuracies in the final properties arising from the eMTP training RMSEs are also rectified by including an upsampled free-energy contribution from the eMTP to the DFT training data. As an added benefit, the eMTP formalism can separately calculate the independent contributions (vibrational and electronic) to the total free energy without any DFT runs.

In the Results section (in Fig. 1), we have shown the superior performance of the eMTP, compared to classical MTPs, in terms of the energy RMSE with respect to DFT. We have also shown that the electronic free energy is a considerable part of the total free energy and that it significantly affects thermodynamic properties. With the eMTP formalism, we are now in a position to estimate and analyze the effect of the electronic temperature on vibrations as highlighted in Fig. 4.

Fig. 4: Nearest-neighbor distribution plots.

a Gauss-broadened first and second nearest-neighbor (1NN and 2NN) distribution function ρ (broadening parameter 0.04 Å⁻¹) for Nb at 2500 K and the corresponding volume, with and without the electronic temperature T_el, calculated using the eMTP (this work) and the classical MTP¹. The 1NN and 2NN distribution functions are plotted as a function of the neighbor distance d around the respective equilibrium distance d^eq. The corresponding 2D distribution histograms are shown in (b–e).

In contrast to energies, the electronic temperature has a negligible effect on the atomic forces, and eventually on the configurational phase-space distribution for Nb. In Fig. 4a, we plot the Gauss-broadened distribution function at 2500 K of the first nearest-neighbor (1NN) and the second nearest-neighbor (2NN), corresponding to the 2D atomic distribution histograms shown in Fig. 4b–e. We plot the 1NN and 2NN distribution with T_el ((b), (c) and green curve in (a))) and without T_el ((d), (e) and red curve in (a)). We are able to obtain a statistically relevant distribution very efficiently, with the use of the eMTP formalism and the classical MTP. Irrespective of whether we use the eMTP or the classical MTP, the distribution functions almost fall on top of each other for both the 1NN and 2NN. The similarity signifies that any difference in atomic forces that comes from the electronic temperature is smaller than the RMSE in the forces of the MTPs, making the difference practically insignificant. This also implies that the configurational space is already well-captured using a classical MTP that is trained only to the vibrational degrees of freedom. The results obtained here corroborate well with a recent study on the effect of electronic temperature on the atomic forces. For a set of bcc refractory unary systems (Mo, Nb, Ta, V and W), Forslund et al.⁶ showed that the root-mean-square deviation in atomic forces even at the melting point was less than 0.15 eV/Å, and as small as 0.06 eV/Å for Nb while including the electronic temperature, thus validating the force predictions of the eMTP formalism.

The eMTP formalism thus primarily rectifies the missing free energies that the classical MTP cannot capture. Hence, in cases where the electronic contribution can be neglected, or while using a more-expensive direct upsampling approach, one can resort to classical MTPs to obtain full thermodynamic properties. One should also note that, while upsampling to DFT, the higher order terms of the free energy perturbation also account for the almost negligible effect of the electronic temperature on the vibrations. In combination with the other generalization of MTPs—the mMTPs²¹, one can also include magnetic effects to the free energy.

The here developed eMTP formalism brings us closer towards entirely modeling the free energy and thermodynamics of a system using an interatomic potential. Since the formalism calculates the electronic contribution as an interpolation of energies at different electronic temperatures, a relatively straight-forward extension of the idea can be applied to other machine-learning potentials in literature. In addition to the vibrational and electronic free energies, defects such as vacancies that can alter the free energies and thermodynamic properties⁶ can also be included. Since the training is on configurations from a reasonably dense volume-temperature grid, the eMTP formalism can also conveniently handle thermodynamic calculations in an NpT ensemble. The benefit of using a machine-learning based model is the ability to improve the robustness of the model by actively improving the training set. For instance, in a recent work, a single MTP was trained to both ordered and disordered phases of the TaVCrW alloy²⁶. The eMTP formalism can be extended to such a training set, thereby simultaneously including also the configurational degree of freedom, in addition to the electronic and vibrational degrees of freedom.

The eMTP formalism can also be used to study kinetically-driven effects with near-DFT accuracy including the effect of electronic excitations, by performing MD simulations for a much longer time on much larger systems. For instance, a single eMTP can predict highly accurate phase transformation behavior, where even the electronic free energy of the different phases affects transformation properties. One can also extend the formalism to predict more accurate melting points for various systems by using the eMTP, for example, in the TOR-TILD scheme²⁷.

Methods

Formalism of the eMTP

The classical MTP potentials represent the energy E_MTP as a polynomial of atomic positions ${{{\boldsymbol{{{{\mathcal{R}}}}}}}}$, corrected so that the local atomic energies do not depend on atoms beyond a certain cutoff. Following this idea, we define ${E}_{{{{\rm{eMTP}}}}}({{{\boldsymbol{{{{\mathcal{R}}}}}}}},{T}_{{{{\rm{el}}}}})$ as a polynomial of both ${{{\boldsymbol{{{{\mathcal{R}}}}}}}}$ and T_el in the following way:

$$\begin{array}{r}{E}_{{{{\rm{eMTP}}}}}({{{\boldsymbol{{{{\mathcal{R}}}}}}}},{T}_{{{{\rm{el}}}}})=\mathop{\sum }\limits_{i=0}^{n}{E}_{{{{{\rm{MTP}}}}}_{i}}({{{\boldsymbol{{{{\mathcal{R}}}}}}}})\,{P}_{i}({T}_{{{{\rm{el}}}}}),\end{array}$$

(1)

where P_i(T_el) are Chebyshev polynomials of the electronic temperature on the interval $[0,{T}_{{{{\rm{el}}}}}^{\max }]$ and ${E}_{{{{{\rm{MTP}}}}}_{i}}({{{\boldsymbol{{{{\mathcal{R}}}}}}}})$ are the energies predicted via a classical MTP for the atomic configuration ${{{\boldsymbol{{{{\mathcal{R}}}}}}}}$. In this form, the construction of an eMTP requires classical MTPs that are trained via training sets with fixed electronic temperatures. The chosen set of electronic temperatures needs to be well-distributed in the relevant temperature range in order to obtain an accurate interpolation. A pragmatic way is to choose temperatures corresponding to Chebyshev nodes on the interval $[0,{T}_{{{{\rm{el}}}}}^{\max }]$, where each ${E}_{{{{{\rm{MTP}}}}}_{i}}$ (i = 0…n) corresponds to one of the n + 1 Chebyshev nodes. Nonetheless, the final results are independent of the chosen temperature set as long as enough well-distributed temperature values are included. Thus, the process of calculating the potential energy and interatomic forces of a given configuration ${{{\boldsymbol{{{{\mathcal{R}}}}}}}}$ involves calculating such quantities via each ${E}_{{{{{\rm{MTP}}}}}_{i}}({{{\boldsymbol{{{{\mathcal{R}}}}}}}})$ and then interpolating these quantities for a given electronic temperature.

Training set generation

The training data for the MTPs that make-up the eMTP formalism come from DFT calculations. The DFT calculations are performed using VASP^28,29, with the projector augmented wave (PAW) potentials³⁰ and the generalized gradient approximation (GGA)^31,32,33. For Nb, we use 54-atom and 125-atom supercells, and for TaVCrW, we use 128-atom supercells. We use an energy cut-off of 350 eV and sample the reciprocal space with a 2 × 2 × 2 k-point mesh and a Fermi broadening scheme²², with the smearing width corresponding to the electronic temperature of the training data.

The eMTP formalism has MTPs trained at 6 different electronic temperatures in the range from 0K to the melting point (${T}_{{{{\rm{el}}}}}^{\max }={T}_{{{{\rm{melt}}}}}$). For Nb, we choose the experimental melting point of 2750K³⁴. For the TaVCrW HEA, the experimental melting point is unavailable. Therefore, we use the eMTP formalism to predict an over-heated melting point by performing a solid-liquid heating simulation in MD using LAMMPS³⁵. We then proceed to use 2500K as the approximate melting point—a value that is ≈150K lower than the over-heated value. The 6 electronic temperatures chosen for Nb and TaVCrW are mentioned in Table 1; they fall in the [0, T_melt] range. For Nb, we choose the roots of the Chebyshev polynomials of degrees 2 and 3. Additionally, we also include 2000K to the temperature set in order to obtain a better distribution. For TaVCrW, we use the roots of the Chebyshev polynomial of degree 6 as the set of electronic temperatures.

The configurations for the training set are generated from MD simulations using previously fitted classical MTPs^1,4. MD simulations are run using LAMMPS at different temperatures—1500K, 2000K, and 2750K for Nb, and the set of six temperatures in Table 1 for TaVCrW. At each temperature, we choose five to eight different volumes to generate the configurations: (i) the equilibrium volume at the highest temperature (2750K for Nb and 2500K for TaVCrW) from previously calculated DFT lattice expansion coefficients^1,4, (ii) three to six smaller volumes, corresponding to lattice constants that are 0.02 to 0.12 Å smaller than the lattice constant at the highest temperature, and (iii) one larger volume corresponding to a lattice constant that is 0.02 Å larger than highest-temperature lattice constant. Each MD simulation with the classical MTP is run for 10,000 steps. Following this, configurations are selectively added to the initial training set that the classical MTP was trained to (refer to Jung et al.¹ and Zhou et al.⁴ for the training set of the classical MTPs). At this stage, the actual accuracy of the model used to generate the configurations is not crucial as long as sufficient configurations have been chosen that span the necessary configurational space. In total, we end up with 1150 configurations for Nb and 1620 configurations for TaVCrW. The energies and forces of all configurations are calculated at each electronic temperature using the above-mentioned DFT parameters, and they form the training set for the MTPs that make-up the eMTP.

Electronic MTP fitting

At each electronic temperature, we fit an MTP to the training data using the MLIP code^18,36. We choose a level-24 MTP for both systems. The maximum cut-off is 5Å, above which interactions are not considered. The fitting weights for energy and force are 1 and 0.01 Å ², respectively. The RMSEs in energy and forces of the individual MTPs at each electronic temperature are shown in Table 1. Once the MTPs are obtained for the eMTP formalism, we use them to perform TILD calculations to calculate the total free energy.

eMTP free-energy surface generation

The reference for the TILD calculations is an effective harmonic potential^24,25 fitted to configurations from a low-temperature MD, whose energies and forces are recalculated with stricter DFT parameters. For Nb, we use a 20 K effective harmonic and for TaVCrW, we use a 500 K effective harmonic potential. The effective harmonic potentials are fit using the HIPHIVE³⁷ and SPHINX³⁸ code for Nb and TaVCrW, respectively. The dynamical matrix is obtained by minimizing the error on forces coming from DFT calculations at the corresponding temperature (20 K and 500 K for Nb and TaVCrW) on 512-atom and 128-atom supercells. During the fitting, a cut-off radius that includes the first three neighbor shells is used to restrict the interactions. The effective harmonic matrices are then used as reference for the TILD calculations. Further details about the effective harmonic potential fitting can be found in previous articles^1,4.

We perform TILD calculations on a highly dense (V, T) grid (11 × 13 for Nb and 11 × 10 for TaVCrW) in the relevant volume-temperature range. For Nb, the temperature points in kelvin are 2 750/7 × {0, 0.33, 0.67, 1, 1.5, 2, 3, 4, 5, 5.5, 6, 6.5, 7, 7.5}, and the volumes are equidistant between 18.18 and 20.37 ų/atom. For TaVCrW, we choose every 250 K up to 2750K and volumes equidistant between 14.75 and 17.48 ų/atom. To make comparisons easier, we choose the same grid as in previous calculations^1,4. Such a high-density grid is a prerequisite for parametrizing the free energy surface in order to obtain well-converged thermodynamic properties. At each (V, T) point in the grid, we perform TILD using the eMTP formalism in two steps. First, we perform TI from the effective harmonic model to the eMTP, but with the electronic temperature turned off (T_el = 0 K). The free energy difference at this stage captures the anharmonic vibrational free energy F^ah(V, T) that is not captured by the effective harmonic reference. Next, we perform TI from the eMTP without T_el to the eMTP with T_el set to the ionic temperature T_ion. The free energy difference at this stage accounts for the electronic free energy F^el(V, T), which includes the coupling to vibrations. The TILD runs in both steps are performed in 432-atom supercells for Nb, and 128-atom supercells for TaVCrW, at 22 λ-points from 0 to 1, for 30,000 steps with a time-step of 5 fs. The free energy difference obtained in both steps is summed with the reference effective harmonic energy F^qh, and the 0 K energy-volume curve E_0K(V), in order to obtain the full free energy, as given by

$$F(V,T)={E}_{0{{{\rm{K}}}}}(V)+{F}^{{{{\rm{qh}}}}}(V,T)+{F}^{{{{\rm{ah}}}}}(V,T)+{F}^{{{{\rm{el}}}}}(V,T),$$

(2)

with

$${F}^{{{{\rm{ah}}}}}(V,T)=\int\nolimits_{0}^{1}d\lambda {\left\langle {E}^{{{\rm{eMTP}}}}({T}_{{{{\rm{el}}}}} = 0{{{\rm{K}}}})-{E}^{{{\rm{qh}}}}\right\rangle }_{\lambda },$$

(3)

and

$${F}^{{{{\rm{el}}}}}(V,T)=\int\nolimits_{0}^{1}d\lambda {\left\langle {E}^{{{\rm{eMTP}}}}({T}_{{{{\rm{el}}}}} = {T}_{{{{\rm{ion}}}}})-{E}^{{{\rm{eMTP}}}}({T}_{{{{\rm{el}}}}} = 0{{{\rm{K}}}})\right\rangle }_{\lambda }.$$

(4)

Here, E^eMTP and E^qh are the corresponding internal energies using the eMTP and the effective harmonic potential. The integral runs over the coupling parameter λ. Once we obtain the total free energy contributions at the grid points, we parametrize individual contributions separately. The reference energy F^qh(V, T) is parametrized using a third-order polynomial in volume. The anharmonic free energy F^ah(V, T) is written as a function of renormalized frequencies that are parametrized using a fourth-order polynomial in volume and temperature. The electronic free energy F^el(V, T) is directly parametrized using a fourth-order polynomial in volume and temperature. The contributions are summed up to obtain the total free energy surface F(V, T) in Eq. 2.

Upsampling from eMTP to DFT training data

In addition to calculating thermodynamic properties predicted up to the level of accuracy of the eMTP (green solid lines in Fig. 2), we also account for the free-energy difference coming from the RMSE in energies and forces of the eMTP (green dashed lines in Fig. 2). We achieve this without performing any additional DFT calculations, but instead using existing data that was used for training the eMTP. As discussed above in the subsection ‘Training set generation’, the configurations in the training set are generated at a set of (V, T) points (in particular, three T’s and five V’s for Nb, and six T’s and eight V’s for TaVCrW). For a given T_el, at each (V, T), we calculate an ‘upsampled’ contribution ΔF^eMTP,up as given by,

$$\begin{array}{r}\Delta {F}^{{{{\rm{eMTP}}}},{{{\rm{up}}}}}(V,T)=-{k}_{B}T\ln \left\langle \exp \left(-\frac{{E}^{{{{\rm{DFT}}}}}-{E}^{{{{\rm{eMTP}}}}}}{{k}_{B}T}\right) \right\rangle .\end{array}$$

(5)

Here, E^DFT and E^eMTP are DFT energies (including electronic temperature), and eMTP energies, respectively, and the averaging is performed over the configurations in the training set. The DFT energies are already available from the training set. The upsampled contribution ΔF^eMTP,up is parametrized using a second-order polynomial in volume and temperature, followed by a parametrization in T_el. A parametrization in T_el can be avoided if we have generated configurations for the training data at the set of chosen electronic temperatures (as we do for TaVCrW). Nonetheless, we obtain an eMTP-upsampled free-energy surface that is added to the previously calculated free-energy surface F(V, T) in Eq. (2).

Thermodynamic properties

A Legendre transformation of F(V, T) provides the total Gibbs energy at ambient pressure p as given by G(p, T) = F(V, T) + pV. Necessary numerical derivatives are taken to obtain the thermodynamic properties using the following formulas:

• Isobaric heat capacity C_p

  $${C}_{p}(p,T)=T{\left(\frac{\partial S(p,T)}{\partial T}\right)}_{p},$$

  (6)
  $$S(p,T)=-{\left(\frac{\partial G(p,T)}{\partial T}\right)}_{p}.$$

  (7)

• Lattice expansion coefficient α

  $$\alpha (p,T)=\frac{1}{3{V}^{{{{\rm{eq}}}}}(p,T)}\left(\frac{\partial {V}^{{{{\rm{eq}}}}}(p,T)}{\partial T}\right),$$

  (8)
  $${V}^{{{{\rm{eq}}}}}(p,T)={\left(\frac{\partial G(p,T)}{\partial p}\right)}_{T}.$$

  (9)

• Adiabatic bulk modulus B_S

  $${B}_{S}(p,T)=\left\{\begin{array}{ll}{B}_{T}(p,T)\dfrac{{C}_{p}}{{C}_{V}}&{C}_{V}(p,T)\, >\, 0.01\,{k}_{B},\\ {B}_{T}(p,T)&{C}_{V}(p,T)\, <\, 0.01\,{k}_{B}.\end{array}\right.$$

  (10)
  $${B}_{T}(p,T)=\frac{1}{{\kappa }_{T}},\quad{\kappa }_{T}=-\frac{1}{{V}^{{{{\rm{eq}}}}}(p,T)}{\left(\frac{\partial {V}^{{{{\rm{eq}}}}}(p,T)}{\partial p}\right)}_{T}.$$

  (11)

Here, S is the entropy, V^eq is the equilibrium volume, B_T is the isothermal bulk modulus, C_V is the isochoric heat capacity, and κ_T is the compressibility.