Universal machine learning interatomic potentials are ready for phonons

Dataset and its properties

To benchmark phonon properties we use the dataset developed in the MDR database³⁴. This dataset includes around 10 000 non-magnetic semiconductors, covering a wide range of elements across the periodic table. Moreover, the phonon calculations were performed with VASP, ensuring a high degree of compatibility with the training sets used in the construction of the uMLIPs. Unfortunately, this phonon dataset was originally constructed with the Perdew-Burke-Ernzerhof (PBE) for solids (PBEsol)³⁵ approximation to the exchange-correlation functional. This is certainly a very reasonable choice, as the PBEsol functional exhibits superior structural^36,37 and phonon³⁸ properties when compared to the standard PBE³⁹. However, as all uMLIPs were trained on PBE data, a direct comparison to PBEsol phonons can be ambiguous. To mitigate this problem, we recalculated the entire phonon dataset from ref. ³⁴ with the PBE functional (see Section IV). In the following, we not only present comparisons of uMLIP calculations with PBE data, but we also include the difference between PBE and PBEsol. This gives us an estimate of the variability of the results as a function of the approximation to the exchange-correlation function, that we use as an absolute scale to assess the quality of the uMLIPs.

As illustrated in Fig. 1a, the dataset contains mostly ternary and quaternary compounds. Additionally, we observe that the majority of the compounds belong to the monoclinic and orthorhombic crystal systems, followed by approximately equal proportions of trigonal and tetragonal systems. Cubic systems are less common, with hexagonal systems representing the smallest proportion. Ultimately, these characteristics are inherited from the Materials Project database¹² and the Inorganic Crystal Structure Database (ICSD)⁴⁰. Finally, triclinic systems are absent from the MDR database, likely because of the extra computational cost that arises from the reduced symmetry.

In Fig. 2 we plot the frequency of the chemical elements in the dataset. We can see that almost all the periodic table is well represented (with a few exceptions like Tc that is radioactive or Eu and Gd for which VASP has convergence problems). We also observe a significant abundance of structures containing oxygen. However, certain compounds, such as those containing Mo and W, as well as the magnetic 3d elements (from V to Ni) are underrepresented. These biases in the MDR database³⁴ are also, to some extent, inherited from the Materials Project database¹², but should not be relevant for the benchmark we present here. Although the dataset is predominated by oxides, the band gaps of the whole set still covers a large range, as illustrated in Fig. 1c.

Relative performances of uMLIPs

We start by discussing the errors in the geometry relaxations, as shown in Table 1 and Fig. 3. The “Failed” column in Table 1 indicates for how many systems a model failed to converge the forces to below 0.005 eV/Å. We can see that CHGNet and MatterSim-v1 models appear to be the most reliable, with 0.09% and 0.10% unconverged structures, respectively. The M3GNet, SevenNet-0 and MACE-MP-0 models have a similar number of unconverged structures, while the ORB and eqV2-M models exhibit a much larger failure rate. The most unreliable model for this dataset is eqV2-M, for which 0.85% structural calculations were unable to converge. In general, there are two main reasons behind the failures, either the geometry optimization path explored regions of the potential energy surface where the uMLIP yielded unphysical forces, or there were high frequency errors in the forces that prevented the relaxation algorithm to converge to the required precision. This latter reason is behind the very large failure rate for the two models where the forces are not the exact derivatives of the energy. CHGNet shows notably higher error in energy predictions, which is expected given that we did not apply the energy correction procedure typically used during CHGNet’s training.

Looking at Fig. 3 we see that, as expected, PBEsol leads to a contraction of the unit cell, correcting the underbinding that is typical of the PBE approximation. The large majority of the systems show a difference between the PBE and PBEsol volume per atom between 0 and –2 ų/atom. All uMLIPs exhibit MAE(V) that are smaller than the mean absolute difference between PBE and PBEsol. Among them, the eqV2-M model emerges as the most accurate, closely followed by ORB. Indeed, these two uMLIPs show remarkable performances for the vast majority of the compounds in the dataset, with errors that are quite small in both absolute and relative terms. MatterSim-v1 and SevenNet-0 show solid performances, although with mean errors four times larger than the two best models. Finally, M3GNet, MACE-MP-0, and CHGNet have wider error distributions, with MAE in the range of 0.4–0.5 ų/atom. These results confirm that both eqV2-M and ORB are the best models for geometry optimization, and that they can already be used to essentially replace DFT calculations for this task.

We now turn our attention to phonon related properties. We chose to look at the maximum phonon frequency (reported in Kelvin, with 1 K = 0.695 cm⁻¹), the phonon density of state (DOS), the average of the sound velocity on the 3 accoustic branches, the vibrational entropy, the Helmholtz free energy, and the heat capacity at constant volume, the last three calculated at the temperature of 300 K. The maximum phonon frequency allows us to detect systematic errors in the prediction of the concavity of the potential energy surface, especially important as it is well known that some uMLIPs have the tendency to yield too soft phonons. The phonon DOS provides information regarding the general prediction of phonon modes with respect to frequency, while the sound velocity help identify errors in the acoustic branches in the vicinity of Γ. It should be noted that for the phonon DOS we remove values below 0.1 states/THz. The vibrational entropy and the Helmholtz free energy are important properties as they are essential to determine thermodynamic stability and phase diagrams as a function of temperature. Finally, the heat capacity is an important thermal property that can be directly measured experimentally.

We note that maximum phonon frequency was calculated from the values at the q-points commensurate to the supercell matrix, whereas the DOS and thermodynamic properties were obtained on an denser q-grid by applying Fourier interpolation (see Section IV). However, as the q-grids are consistent across DFT and uMLIPs calculations, the interpolation error should be systematic and should not affect the benchmark.

We aggregated the errors for all models in Table 2 and in Fig. 4. We first notice that the deviation between the PBE and PBEsol results is small but not negligible. This observation reinforces the necessity of using a consistent functional between the training and benchmarking stages. The difference between PBE and PBEsol exhibits a rather narrow distribution in all 6 properties, especially when compared to the MAE of most of the uMLIPs. There are also systematic differences: for example, the maximum phonon frequencies in PBEsol are higher than those with PBE, which can be understood by the contraction of the cell and subsequent hardening of the force constants. PBEsol also leads to larger values of the free energy (on average of the order of 10 kJ/mol), and to smaller values of the entropy and the heat capacity.

Based on the errors we can roughly classify the seven models into three categories. The first contains ORB and eqV2-M, which have very large errors in phonon-related properties (see Fig. 4). In fact, phonon frequencies are grossly underestimated, and are often even imaginary as we will see in the following.

In the second category we have, in increasing order of accuracy, M3GNet, CHGNet, MACE-MP-0 and SevenNet-0 (see Fig. 4). The errors of these models are on average considerably larger than the difference between PBE and PBEsol. Moreover they all exhibit systematic errors, underestimating the phonon frequencies and the free energy, and overestimating the entropy and the heat capacity. From the four models, the most accurate is clearly SevenNet-0, while the older M3GNet and CHGNet show the larger errors. In spite of the difference in topologies, these four models are all trained in the same dataset, so it is not surprising that their results are somewhat similar. This again demonstrates that training data is at least as important as the representation of the crystal structure or the topology of the model to develop a uMLIP.

Finally, MatterSim-v1 stands out as the most accurate uMLIP for the calculation of phonons. Not only does it not exhibit any strong systematic error, with all distributions essentially centered at zero, but also the dispersion of the errors is extremely small, leading to values of MAE considerably smaller than the difference between PBE and PBEsol. This indicates that MatterSim-v1 can be used to calculate phonon properties of semiconductors with an accuracy comparable to DFT codes, although at a very small fraction of the computational cost. It is very interesting to note that although MatterSim-v1 is based upon the simple M3GNet, its performance exceeds much more complicated models such as SevenNet-0 or eqV2-M that are based on equivariant networks. The key in this case is the scalability of M3GNet, which allows for an increase in the number of parameters and the efficient use of larger amounts of training data.

To have a better understanding of the general behavior of the uMLIPs, we plot in Fig. 5 the distribution of the maximum frequencies predicted. Most compounds have maximum frequencies in the range of 500–2000 K, with a few containing very light elements going up to 5500 K. The softening of the phonon frequencies by M3GNet, CHGNet, MACE-MP-0 and SevenNet-0 is evident, in particular for the first two. ORB and eqV2-M, on the other hand, exhibit completely distorted distributions peaking at zero, showing that the force constants obtained with these models are unphysical.

Another important performance metric is dynamical stability, a crucial stability descriptor utilized by many high-throughput searches of inorganic materials^41,42,43,44. A compound is dynamically stable when it is in a true minimum of the potential energy surface and not in a maximum or a saddle point. In practice, it is assured by the absence of imaginary phonon frequencies in the spectrum. Unfortunately, it is well known that numerical inaccuracies often lead to small imaginary frequencies close to the Γ-point. To avoid this problem, we consider a structure to be dynamically stable if frequencies are all real across the Brillouin zone except at Γ where we allow the three acoustic modes to have small imaginary frequencies (with a threshold of −50 K). This criterion was applied to all q-points commensurate with the supercell matrix (but not to the interpolated q-points).

The elements of the confusion matrix, when compared to the PBE, are listed in Table 3. Most compounds that are stable in the PBE are also stable in PBEsol, and vice-versa, with the differences coming mostly from the difficulty associated to small imaginary frequencies as mentioned above. MatterSim-v1 and MACE-MP-0 are the most reliable with a percentage of true positives at 95%. M3GNet, SevenNet-0 and CHGNet are somewhat less accurate, especially in what concerns the percentage of true positives. Finally, the eqV2-M and ORB models perform very poorly, with more than 80% of the unstable systems being false negatives.