Search for thermodynamically stable ambient-pressure superconducting hydrides in the GNoME database

TL;DR

Researchers searched for stable hydrides that could superconduct at normal pressure and found 25 with superconducting potential, up to 17 K.

Contribution

The study identifies 25 stable hydrides with superconductivity at ambient pressure using the GNoME database and ab initio methods.

Findings

25 cubic hydrides with superconducting critical temperatures above 4.2 K were identified.

The highest critical temperature found was 17 K in stable hydrides.

Stable hydrides are likely experimentally accessible and technologically relevant.

Abstract

Hydrides are considered to be one of the most promising families of compounds for achieving high temperature superconductivity. However, there are very few experimental reports of ambient-pressure hydride superconductivity, and the superconducting critical temperatures (Tc) are typically less than 10 K. At the same time several hydrides have been predicted to exhibit superconductivity around 100 K at ambient pressure but in thermodynamically unfavorable phases. In this work we aim at assessing the superconducting properties of thermodynamically stable hydride superconductors at room pressure by investigating the GNoME material database, which has been recently released and includes thousands of hydrides thermodynamically stable at 0K. To scan this large material space we have adopted a multi stage approach which combines machine learning for a fast initial evaluation and cutting edge ab initio methods to obtain a reliable estimation of Tc. Ultimately we have identified 25 cubic hydrides with Tc above 4.2 K and reach a maximum Tc of 17 K. While these critical temperatures are modest in comparison to some recent predictions, the systems where they are found, being stable, are likely to be experimentally accessible and of potential technological relevance.

The recently published GNoME material database contains a large number of new thermodynamically stable hydride compounds. In this work the authors assess their superconducting properties using ab initio methods and identify 25 cubic hydrides with a maximum critical temperature of 17 K.

Figures

The computational workflow of the screening of cubic superconducting hydrides with thermodynamic stability by combining the ab initio study with the ALIGNN model.The integers in parentheses indicate the number of hydrides.

The conventional cells of the hydrides are selected from vacancy-ordered double perovskites and fluorite-like structures.a Normal double perovskite EuDyReH_6_ with vacancy ordering. b Inverted double perovskite Ce_2_HRh_6_ with vacancy ordering. c Fluorite-like Ta_3_NbH_8_. d Fluorite-like Ta_6_NbH_16_ with vacancy. The gold spheres represent the ordered vacancies *V*_vac_.

Properties of LiZrH_6_Ru at ambient pressure.Electronic structure (a), crystal structure (b), and phonon and electron–phonon coupling (c). The gray dashed lines in the diagrams of phonon spectrum and Eliashberg spectral functions were calculated by including the ionic quantum and anharmonic effect in the stochastic self-consistent harmonic approximation (SSCHA).

Electronic, phonon, and electron–phonon coupling properties of LiZrH_6_Ru.Top panels a refer to 80 GPa and bottom panels b to 160 GPa. In either panel a or b, the left plot shows the electronic band structure and the atom projected density of states, and the right plot displays the phonon band structure, the atom projected phonon density of states, the electron-phonon spectral function *α*^2^*F*(*ω*) and its integration curve *λ*(*ω*).

Coupling properties of LiZrH_6_Ru.a Electronic bands near the Fermi level, including atomic projections on the Zr and Ru sites. b Fermi surface and the electron-phonon coupling (*λ*_k_). c Coulomb interaction matrix resolved on iso-energy surfaces. d Density of electronic states.

Full text

Introduction

The near room temperature superconductivity, experimentally observed in several superconducting hydrides such as LaH_10_^1,2^ and YH_9_^3^ hosting superconductivity above −30 and −23 °C, respectively, demonstrates that the investigation of phonon-mediated superconductors involving the lightest element H is well merited^4^. However, achieving near room temperature superconductivity in these hydrides requires subjecting them to extreme pressures over 170 billion pascals (170 GPa), making the real-world application unlikely.

Given the tough conditions of these high-Tc superconducting hydrides posed by the high pressures, one may naturally wonder whether it is possible to achieve superconductivity in the hydrides at ambient pressure^5^. In earlier studies, very few hydride superconductors have been experimentally reported at ambient pressure, and their superconducting transition temperatures are typically less than 10 K^6–8^. The desire to discover high-Tc ambient-pressure superconducting hydrides has driven some ab initio theoretical predictions recently^9–16^. In particular, Cerqueira et al. combined the machine learning methods with the ab initio study to predict around 50 superconducting hydrides with Tc above 20 K from the screening of more than 1 million compounds^9^. The cubic hydride family with MH_6_ octahedra is exceptionally attractive because some of them are near the convex hull and can maintain a high Tc above the boiling point of liquid nitrogen of 77 K^10,11^. Several perovskite hydrides have also attracted some attention due to their high Tc and small enthalpy distance from the convex hull^9,17,18^. Remarkably, ab initio calculations suggest that RbPH_3_ is dynamically and kinetically stable at ambient pressure with a superconducting critical temperature of about 100 K^19^.

As discussed in refs. ^20^ and ^15^, the evidence that emerges from recent high-throughput research is that high-Tc values of the order of 100 K appear to be systematically linked to some degree of thermodynamical instability. This may be attributed to the destabilizing effect of delocalized hydrogen-derived electrons in the metallic bonding environment of multinary hydrides. Consequently, the hydrides thermodynamically stable at 0K in large materials databases such as Materials Project are generally semiconductors or insulators. Clearly, this drops a large fraction of the problem to the synthesis, as one expects that a synthesis route for metastable compounds might be very hard to find, although some recent advances in pressure-quench protocols suggest that it may be possible to retain superconductivity in the ambient-pressure metastable phases quenched from a pressure where the corresponding phases are thermodynamically stable^21^.

GNoME has used active learning and graph neural networks to identify 381 thousand crystals that are stable at 0 K according to DFT^22^. Most of these structures are identified by elemental substitutions, a commonly used method to search for 0 K stable crystals^23–25^. As mentioned in the conclusion of GNoME, these calculations are not able to consider configurational disorder. With further improvements in force fields and computational resources, simulating finite temperature stability in a high-throughput manner might be possible. Given this set of crystals that are stable at 0 K according to DFT, might there be interesting superconducting candidates? In this work, we investigate this possibility by means of first-principles simulations.

A brute force use of ab initio methods^26^ would be computationally too expensive to be directly applicable to such a large space of materials. For this reason, we adopt a recently developed machine-learning accelerated ab initio algorithm. This allows us to identify possible high-Tc hydrides in the GNoME database and use high-level ab initio methods to validate machine learning predictions. As we will report below, the search has been successful in identifying 25 cubic hydrides with Tc exceeding the boiling point of liquid helium (4.2 K). Most of these hydrides have Tc lower than NbTi alloys (10 K), nevertheless the cubic LiZrH_6_Ru hydride with vacancy-ordered double perovskite structure stands out and shows a Tc at 17 K at ambient pressure.

Results

Machine learning accelerated a high-throughput ab initio approach

To accelerate the discovery of superconducting hydrides, we have combined the machine learning methods with high-throughput ab initio calculations. The major workflow of the superconductors screening is described in Fig. 1.Fig. 1. The computational workflow of the screening of cubic superconducting hydrides with thermodynamic stability by combining the ab initio study with the ALIGNN model.The integers in parentheses indicate the number of hydrides.

We first performed data mining on thermodynamically stable hydrogen-based compounds at 0K with a cubic structure and metallic properties from the openly available dataset of GNoME project^27^, which resulted in a subset including 490 hydrides. We focused exclusively on cubic hydrides because the cubic structure is a recurring and well-established structural motif across many classes of high-Tc hydrogen-based compounds at both high pressure (e.g., LaH_10_^1,2^ and YH_9_^3^) and ambient pressure (e.g., RbPH_3_^19^ and Li_2_AgH_6_^15^). High-throughput spin-polarized density functional theory (DFT) calculations were further used to find 261 nonmagnetic metallic hydrides. Then, the Atomistic Line Graph Neural Network (ALIGNN) model^28^, which has been proven successful in the discovery of superconductors^9,10^, was used to estimate Tc of the hydrides. To narrow down the search space of superconductors and reduce the computational cost of ab initio calculations, only the superconducting hydrides with ALIGNN-predicted \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{{{{\rm{pred}}}}}$$\end{document}  ≥ 1 K were considered in the subsequent high-throughput density functional perturbation theory (DFPT) calculations.

Using high-throughput DFPT methods coupled with Allen–Dynes formula, we found 25 cubic hydrides from the GNoME database with a superconducting critical temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document}  ≥ 4.2 K (see Tables 1 and 2). For the 25 cubic hydrides, the DFPT-calculated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} and ALIGNN-predicted \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{{{{\rm{pred}}}}}$$\end{document} are compared in a parity plot included in Supplementary Fig. S1. The values of ALIGNN-predicted \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{{{{\rm{pred}}}}}$$\end{document} are also included in Tables 1 and 2 explicitly. The ALIGNN predicted Tc shows a strong quantitative agreement with DFPT-calculated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} , with a mean absolute error (MAE) of only 2.5 K. Given that the MAE of 2.5 K is significantly smaller than the target benchmark (the boiling point of liquid helium, 4.2 K), this predictive uncertainty is acceptable for pre-screening purposes. In addition, the conservative criterion of ALIGNN-predicted \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{{{{\rm{pred}}}}}$$\end{document}  ≥ 1 K should ensure that all materials with true superconducting critical temperature above 4.2 K are retained in the high-throughput DFPT calculations.Table 1. Vacancy-ordered double perovskitesmat_idformulaλ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{\log }$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{{{{\rm{pred}}}}}$$\end{document} Vacancy-ordered double perovskitesa1c4087396EuDyReH_6_2.0535.844.92.52523d9ecdeNdEuTcH_6_0.66166.665.13.91243cdd0e1Ce_2_HRh_6_0.67167.005.11.697ff896f6fSmEuTcH_6_0.67168.285.36.20eddd1b470EuCdFeH_6_0.47502.935.44.3cc52a3846dEuYTcH_6_0.67179.975.65.7e5427bdb32TbEuTcH_6_0.73168.496.36.884f266aecdEuDyTcH_6_0.74166.476.57.79248b5ab15EuHoTcH_6_0.77163.216.99.14bb5b78a22SmEuReH_6_0.98105.756.92.968f1be98a3LiCeMnH_6_0.65234.446.92.75e1177e128EuErTcH_6_0.81160.527.69.58d4721d779EuYReH_6_1.2884.827.73.3f68afdcbbcEuTmTcH_6_0.88155.988.49.31e099eabf5YbCeTcH_6_0.90150.528.66.11bd4b03f35YbGdTcH_6_1.00144.049.610.28bf27d9be2LiTiH_6_Ru1.12123.389.88.7e5c2a73dc0EuLuTcH_6_1.09134.3810.19.251a97ab68dEuCdH_6_Ru2.5384.4013.66.70f97a7734cLiZrH_6_Ru1.00341.5923.54.7Superconducting transition temperatures calculated with the Allen–Dynes formula ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} in K) above 4.2 K, logarithmic average of the phonon frequencies ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{\log }$$\end{document} in K), and the electron–phonon coupling constant (λ) for vacancy-ordered double perovskites. The ALIGNN-predicted \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{{{{\rm{pred}}}}}$$\end{document} (K) are included in the last column for comparison with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} .Table 2. Fluorite-like hydridesmat_idFormulaλ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{\log }$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{{{{\rm{pred}}}}}$$\end{document} Fluorite-liked8e8825efeTa_3_VH_8_0.64178.584.94.1149cfed18bTa_3_NbH_8_0.63190.375.04.434fc41965bTaNb_3_H_8_0.72185.957.07.1Vacancy-ordered fluorite-like140fb04a4bTa_6_NbH_16_0.64179.615.06.0d93c6d6bdaTa_6_MoH_16_0.75176.757.36.0Superconducting transition temperatures calculated with the Allen-Dynes formula ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} in K) above 4.2 K, logarithmic average of the phonon frequencies ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{\log }$$\end{document} in K), and the electron-phonon coupling constant (λ) for fluorite-like hydrides and vacancy-ordered fluorite-like hydrides. The ALIGNN-predicted \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{{{{\rm{pred}}}}}$$\end{document} (K) are included in the last column for comparison with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} .

By examining the crystal structures of the 25 thermodynamically stable superconducting hydrides, we find that they can be grouped into two prototype structure families, i.e., vacancy-ordered double perovskites and fluorite-like structure with or without vacancy-ordering.

All the vacancy-ordered double perovskites can be generalized as defect forms of doubly cation-ordered \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$AA{\prime} BB{\prime} {X}_{6}$$\end{document} (see Fig. 2a). In the defect-free \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$AA{\prime} BB{\prime} {X}_{6}$$\end{document} , X is stabilized at 24f Wyckoff position; A and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A^{\prime}$$\end{document} are located at 4c and 4d Wyckoff positions, and are bonded to twelve equivalent X to form A**X12 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A^{\prime} {X}_{12}$$\end{document} cuboctahedra; B and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B^{\prime}$$\end{document} reside at 4a and 4b Wyckoff positions, and each B or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B^{\prime}$$\end{document} is bonded to six equivalent X ions to form B**X6 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B^{\prime} {X}_{6}$$\end{document} octahedra. In the hydrides with vacancy-ordered double perovskite structure, the number of atomic sites is reduced to nine per formula unit due to the existence of vacancies at either 4a or 4b Wyckoff positions. Figure 2a displays EuDyReH_6_ as an example of a vacancy-ordered double perovskite where the vacancy appears at the 4a position and hydrogen atoms are located at the X site. In contrast to vacancy-ordered double perovskites with X being H and space group being \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F\bar{4}3m$$\end{document} (No. 216), the hydrogen atoms are swapped with the metallic Rh atoms in Ce_2_HRh_6_, forming the inverted double perovskite with space group of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Fm\bar{3}m$$\end{document} (No. 225), where 4c and 4d positions are both occupied by Ce and the vacancy is located at the 4b position. As seen in Table 1, all these vacancy-ordered double perovskites exhibit a superconducting critical temperature ranging from 4.9 K in EuDyReH_6_ to 23.5 K in LiZrH_6_Ru, and most of them show a weak electron–phonon coupling constant λ of less than 1. In addition, these materials show a wide range of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{\log }$$\end{document} from 35.84 K in EuDyReH_6_ to 502.93 K in EuCdFeH_6_. The hydride EuCdH_6_Ru has the largest λ of 2.53, but its critical temperature remains as small as 13.6 K due to the very small \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{\log }$$\end{document} of 84.40 K.Fig. 2. The conventional cells of the hydrides are selected from vacancy-ordered double perovskites and fluorite-like structures.a Normal double perovskite EuDyReH_6_ with vacancy ordering. b Inverted double perovskite Ce_2_HRh_6_ with vacancy ordering. c Fluorite-like Ta_3_NbH_8_. d Fluorite-like Ta_6_NbH_16_ with vacancy. The gold spheres represent the ordered vacancies Vvac.

Except for vacancy-ordered double perovskites, the remaining materials are fluorite-like (CaF_2_-like) structured. For example, in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Pm\bar{3}m$$\end{document} Ta_3_NbH_8_ with 12 atomic sites in the primitive cell, each Ta at 4c position and each Nb at 1a position are bonded in a body-centered cubic geometry to eight equivalent H atoms at the 8c site. In contrast, in some other fluorite-like hydrides, such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Fm\bar{3}m$$\end{document} Ta_6_NbH_16_, the alternative appearance of vacancy and Nb atoms is observed at the 1a site, leading to 23 atomic sites in the primitive cell. We can find from Table 2 that these fluorite-like hydrides show very similar λ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{\log }$$\end{document} with averages of 0.68 and 182 K. As a result, the calculated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} are very close, all around 5.0–7.0 K. It is noted that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P\bar{4}3m$$\end{document} LuNH also features a fluorite-like crystal structure with Tc of around 16 K and λ of around 0.78 at ambient pressure, despite the fact that this phase is very thermodynamically unstable–around 0.5 eV/atom above the convex hull^29,30^.

Improved analysis of LiZrH6Ru using McMillan–Allen–Dynes formulas

The quarterly hydride LiZrH_6_Ru has the highest \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} among the studied materials. As seen in Fig. 3a, there are two bands across the Fermi level. The Li atom is fully ionized and does not contribute to the density of states (DOS) at the Fermi level. Alternatively, the Zr atom dominates the DOS at the Fermi level, with the Ru and H atoms making minor contributions.Fig. 3. Properties of LiZrH_6_Ru at ambient pressure.Electronic structure (a), crystal structure (b), and phonon and electron–phonon coupling (c). The gray dashed lines in the diagrams of phonon spectrum and Eliashberg spectral functions were calculated by including the ionic quantum and anharmonic effect in the stochastic self-consistent harmonic approximation (SSCHA).

Because the Fermi level is located on a shoulder near the peak of the density of states, we performed the improved accuracy calculations and tested the convergence against the smearing, k-grid, and q-grid. In the improved calculations with a k-grid (q-grid) of 42^3^ (6^3^), the converged λ are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{\log }$$\end{document} are 1.20 and 343.32 K, respectively. Compared to the results of the high-throughput calculations shown in Table 1 (i.e., λ ~1.0 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{\log }$$\end{document} ~341.59 K), the electron–phonon coupling constant and the logarithmic average of the phonon frequencies both increase in accurate calculations and result in a higher \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} of 30.70 K.

Since the ionic quantum and anharmonic effects could be crucial to the crystal structure, dynamical stability and the superconducting properties^31^, we have further fully relaxed the crystal structure with the stochastic self-consistent harmonic approximation (SSCHA) method^32–35^. Compared to the structure obtained in the standard DFT of Quantum Espresso, the lattice constants are unchanged (a = 4.604 Å, α = 60° in the primitive cell), and the atomic positions of the metal ions are preserved by the symmetry. However, the hydrogen positions are further optimized in SSCHA due to the ioinc quantum and anharmonic effect. Because only the hydrogen positions are slightly different in the two crystal structures, we only show the fully relaxed crystal structure from SSCHA calculations in Fig. 3b. As clearly seen by the zoom-in inset in Fig. 3b, there are two different types of octahedra, i.e., RuH_8_ and H_8_. The volume of the RuH8 octahedron is larger than that of the H8 octahedron because the central Ru atom in the RuH8 octahedron has a large atomic radius, whereas the H8 octahedron contains a vacant site at its center. Compared to the crystal structure in harmonic approximation, although every hydrogen octahedron remains perfect without showing any distortion in the SSCHA structure, the edges in H_8_ are slightly decreased from 2.178 Å in harmonic approximation to 2.160 Å in SSCHA. On the other hand, the edges in RuH_8_ are slightly increased from 2.426 Å in harmonic approximation to 2.444 Å in SSCHA. The optimizations of the hydrogen atoms in the H_8_ and RuH_8_ octahedra are also implied in the phonon dispersions in Fig. 3c, in which the whole high-frequency region of the pure hydrogen modes above 1250 cm^−1^ is softened by the anharmonicity renormalization. Despite the large downshift in the phonon energy, the accumulated electron-phonon coupling is not significantly changed according to the diagram of Eliashberg spectral functions in Fig. 3c. Consequently, there is only a mild improvement in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} , being 32.0 K with quantum effect.

We also studied the effect of pressure on LiZrH_6_Ru in the range from 20 to 200 GPa. LiZrH_6_Ru is found to be dynamically stable up to 180 GPa in the harmonic approximation and becomes dynamically unstable at 200 GPa due to the appearance of an imaginary mode at the L point (see Supplementary Fig. S2). Our electron-phonon coupling calculations show that the superconducting temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} can be enhanced at pressures up to 180 GPa compared to that at ambient pressure. In particular, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} is boosted to 44 K at 80 GPa. A more specific evolution of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} and λ up to 180 GPa can be found in the Supplementary Fig. S3. In Fig. 4, we show the electronic, phonon, and electron–phonon coupling properties at 80 and 160 GPa. By comparing the electronic bands at pressures of 80 and 160 GPa with those at ambient pressure (i.e., Fig. 3), the bands become more dispersive with the increase in pressure. This is because the pressure makes the atoms closer, enhancing the hopping integral between neighboring orbitals and increasing the band width. As seen from the DOS diagrams in Figs. 3 and 4, the DOS at the Fermi level is mostly contributed by Zr at all pressures. We find that only at 80GPa the Fermi level lies at the sharp peak of the projected DOS of Zr, while the Fermi level is located at the shoulder of the DOS at other pressures. In addition, the diagrams of electron-phonon spectral functions α^2^F(ω) (right plots in Figs. 3 and 4) find that the λ value is 1.22 at 80 GPa that is higher than λ of 1.19 at 0 GPa and λ of 1.11 at 160 GPa. Consequently, the combination of high DOS at the Fermi level with the large λ at 80 GPa is responsible for the highest \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{{{\rm{c}}}}}^{\,{{{\rm{AD}}}}}$$\end{document} at this pressure among the studied pressures.Fig. 4. Electronic, phonon, and electron–phonon coupling properties of LiZrH_6_Ru.Top panels a refer to 80 GPa and bottom panels b to 160 GPa. In either panel a or b, the left plot shows the electronic band structure and the atom projected density of states, and the right plot displays the phonon band structure, the atom projected phonon density of states, the electron-phonon spectral function α^2^F(ω) and its integration curve λ(ω).

In-depth characterization of LiZrH6Ru beyond McMillan–Allen–Dynes formulas

So far, we have estimated the superconducting properties using the McMillan–Allen–Dynes formulas in the μ* approximation and assuming a conventional repulsion strength μ^*^ = 0.1. This approach is usually reliable when dealing with conventional superconductors^26,36^, however, one should consider that it has been developed using data from simple superconducting systems, and it is known to fail in describing systems with complex electronic properties, e.g., Magnesium diboride^37,38^.

LiZrH_6_Ru has several features that hint at the possibility that the McMillan estimation of Tc might not be very reliable. The first thing is that the Fermi level is located close to a peak in the density of states, but close to a band gap. This might lead to a poor Coulomb renormalization^39–41^ and an overestimation of Tc. In simple terms this point can be explained by using the Morel-Anderson formula^42^, which links μ^*^ to the Coulomb repulsion at the Fermi level (μ):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mu }^{* }=\frac{\mu }{1+\mu \,{{{\rm{ln}}}}\left(\frac{{E}_{F}}{{\omega }_{c}}\right)},$$\end{document}

where ωc is a cutoff frequency (arbitrary in Eliashberg theory, close to the maximum frequency of the phonon spectrum in McMillan–Allen–Dynes theory) and EF is the Fermi energy setting the bandwidth of the valence region. In LiZrH_6_Ru, the valence band starts about 10 eV below the Fermi level; however, the presence of a band gap between −1 and −3 eV complicates the matter. In practice, it is not obvious if one should assume EF = 1 or EF = 10. In fact, inspecting the character of the electronic bands, one can see that bands close to the Fermi level are almost exclusively transition metal states, while low-lying states have a dominant H contribution, which is small at the Fermi level.

Second, we see that the Fermi level states have a strong Ru and a Zr component. It is not uncommon that superconductivity in transition metals is affected by incipient magnetism. For example, the critical temperature of Nb (right next to Zr in the periodic table) is significantly reduced by the effect of spin fluctuations^43–45^. Last but not least, as shown in Fig. 5a, there are two electronic bands crossing the Fermi level with different orbital character. Multi-band effect in superconductivity might lead to an increase in Tc as compared to a single band approximation^38,46–48^. The k-resolved electron–phonon coupling (λk) on the Fermi surface is shown^49^ in Fig. 5, showing a substantial coupling anisotropy, although neither as strong or simple as in MgB_2_.Fig. 5. Coupling properties of LiZrH_6_Ru.a Electronic bands near the Fermi level, including atomic projections on the Zr and Ru sites. b Fermi surface and the electron-phonon coupling (λk). c Coulomb interaction matrix resolved on iso-energy surfaces. d Density of electronic states.

The first point can be addressed by computing the critical temperature by solving the Eliashberg equations, including the Coulomb interaction from first principles^50,51^. This approach includes the full electronic bandwidth of the problem and the dielectric screening in the random phase approximation. The Coulomb matrix is shown in Fig. 5. This 2D matrix ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W(\xi ,\xi ^{\prime} )$$\end{document} ) is defined by averaging the Coulomb matrix elements between iso-energetic electronic states^40,45^. By definition, μ = W(0, 0)NF, where NF is the DOS at the Fermi level. In the random phase approximation (RPA), we have μ = 0.66, which is quite large, but not uncommon in transition metal compounds^45^. As one can see, this matrix has approximately a block diagonal structure. The main valence states (panel c4) are strongly repulsive in the −10 to −4 eV range, but they interact weakly with the Fermi level states (panels c3 and c1). This fact, the large value of μ and the presence of a large band gap, lead to an extremely disruptive Coulomb contribution to the superconducting state. The calculated Tc solving the isotropic Eliashberg in RPA is 14.2 K. This estimation is about half of that obtained by the standard McMillan–Allen–Dynes approach, an uncommon deviation which reminds the case of Chevrel phases discussed in ref. ^41^.

A second aspect we would like to further investigate is the possible role of magnetic effects in superconductivity. LiZrH_6_Ru is non-magnetic in the sense that we could not stabilize any ferromagnetic or antiferromagnetic ground state order. However, transition metals might have an incipient tendency toward magnetism, which results in a large spin susceptibility, which enhances the Coulomb repulsion. In extreme cases, this might lead to unconventional superconductivity^52,53^ and in more conventional cases to a reduced value of Tc^43^. This can be included in our ab initio Eliashberg formalism using the Kukkonen-Overhauser approach discussed in refs. ^45,54^. This leads to a Tc estimation of 11.3 K. Considering that, the ab initio inclusion of the spin susceptibility tends to give a stronger Coulomb repulsion even in the absence of significant magnetic instabilities, which is consistent with a marginal role of magnetic effects.

The final aspect we wish to address is the possibility of multiband / anisotropic superconductivity. With this motive, we have computed the two-band resolved α^2^F function and solved the multiband Eliashberg equation. This approach gives an enhanced Tc of about 20%. A nearly identical enhancement was computed using the fully anisotropic superconducting density functional theory (SCDFT) approach^55^, proving that the two-band effect is the leading contribution to the anisotropic enhancement of Tc. As shown in Fig. 5a, the two portions of the Fermi surface are chemically different and spatially disjoint. The inner Fermi surface has a mix of Ru–Zr character, while the outer one has a pure Zr projection. The components of the electron-phonon coupling from the inner and outer bands are λ ~ 0.84 and λ ~ 0.3, respectively. While these values alone could suggest strong anisotropy effects, similar to the case of MgB_2_, unlike MgB_2_, the inter-band coupling is also very large, scattering from the outer band to the inner one has λ = 0.86. In this situation, the band anisotropy is washed out by the electron-phonon coupling, leading to a minor anisotropy contribution to Tc. In the language of the Suhl, Matthias, Walker paper^46^, anisotropy has a small effect on Tc because the maximum eigenvalue of the λ matrix (1.32) is close to the isotropic coupling (1.26). Our best estimation of Tc, including multiband effects and the ab-initio RPA Coulomb interaction, is 17 K.

It is natural to compare this relatively low Tc predicted for LiZrH_6_Ru and those predicted for a similar family of compounds, Mg_2_MH_6_ (M = Rh, Ir, Pd, or Pt) in ref. ^10^ and refs. ^11,14^. The key difference is that in the Mg_2_MH_6_ family, hydrogenic states are present at the Fermi level and provide a significant contribution to the superconducting coupling. On the other hand, in the LiZrH6Ru system, the hydrogen contribution is negligible on the Fermi surface. A generally identical situation has been discussed in detail in ref. ^56^ to explain the different superconducting properties of Mg_2_IrH_6_ and Ca_2_IrH_6_. The fact that the Mg_2_MH_6_ family has high Tc but is thermodynamically unstable while LiZrH_6_Ru is thermodynamically stable is consistent with the recent analysis reported in ref. ^15^ about the limits of Tc in conventional superconductors. On a positive note, the synthesis of LiZrH_6_Ru might be easier owing to the favorable thermodynamics, while that of Mg_2_IrH_6_ is expected to be challenging^57^.

A note about ab initio predictions and their reliability

The predictive power of computational methods in superconductivity research has evolved substantially over the past decades; however, several challenging limitations remain that may affect the confidence of ab initio predictions. Therefore, it is important to discuss what we consider to be the main sources of uncertainty that may lead to discrepancies between computational predictions and experimental realizations.

The first major source of uncertainty lies in the identification of thermodynamically stable phases. The typical error in energy estimations, while strongly material dependent, is of the order of 50 meV atom^−1^^58^. Therefore, if the energy cost for a thermodynamic decomposition is less than this, theoretical predictions should be considered extremely uncertain.

A second source of uncertainty derives from the fact that, even if a phase is truly thermodynamically stable, its synthesis might not achieve a sufficiently pure form so that predictions for a perfect bulk phase still apply. This is mostly a matter of experimental work, when superconductors are of great technological interest (e.g., high temperature superconductors), synthesis methods are investigated with utmost care to achieve samples of the highest quality. In other cases, there might not be a sufficient motivation to invest resources in optimizing the synthesis process. This consideration also applies to metastable phases, which may be accessible only through complex synthesis pathways and are therefore worth investigating primarily when a compound is expected to exhibit exceptional properties. This may be the case for Mg_2_IrH_6_, which is predicted to have an extremely high Tc^9,10^; however, its slight metastability led, upon attempted synthesis, to the formation of Mg_2_IrH_5_, which is unfortunately a semiconductor^57^.

A third source of uncertainty stems from the theory of superconductivity itself and from the state-of-the-art approximations we employ. A precise error estimate is difficult to obtain; however, unless Tc is very low, modern superconductivity methods typically carry an expected uncertainty on the order of 30%^26^. It should also be noted that some studies focus primarily on structural predictions and do not apply sufficient rigor to the prediction of Tc.

Conclusion

We have investigated whether, within the recently released GNoME dataset of predicted stable (at 0 K) compounds, there is any exceptional hydride superconductor with a high critical temperature. Our high-throughput calculations were able to identify 25 superconductors with Tc ranging from the boiling point of liquid helium (4.2 K) to 23.5 K. For the highest-Tc compound, advanced ab initio methods refined Tc down to 17 K. While not extremely high, as compared to some of the other predictions which have appeared in literature, these stand out because of the extremely accurate estimation of thermodynamic stability at 0K (due to a more densely sampled 0K convex hull). At this stage, experimental feedback that can validate both the crystal structures and superconductivity of these DFT-based candidates would greatly benefit the future search for superconductors.

Methods

All the ab initio calculations were computed by the Vienna ab initio simulation package (VASP)^59,60^ and Quantum Espresso (QE)^61,62^. In both VASP and Quantum Espresso calculations, the PBEsol functional, a revised version of the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) for solids^63^, was employed for exchange-correlation effects in all the ab initio calculations. In the QE calculations, the PBEsol pseudo potentials from the strict set of PseudoDojo were employed^64^. For elements with f electrons, it is noted that they are frozen as the core states in the PseudoDojo potentials.

To prescreen the nonmagnetic and metallic hydrides, spin-polarized density functional theory (DFT) methods implemented in the VASP were used, and those calculations were managed and run by the Quacc code^65^. The kinetic energy cutoff of 550 eV was used for the plane-wave basis. In the structure optimizations, a grid density of 1000 per \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{A}}}}{{{\rm{n}}}}{{{\rm{g}}}}{{{\rm{s}}}}{{{\rm{t}}}}{{{\rm{r}}}}{{{\rm{o}}}}{{{\rm{m}}}}}^{-3}$$\end{document} ) of reciprocal cell was used to sample the k-grid. Electronic self-consistency was achieved with an energy convergence threshold of 10^−6^ eV. The relaxations were carried out until the forces on all the atoms were smaller than 0.001 eV/Å.

We trained the ALIGNN^28^ model using as targets, simultaneously, λ, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{\log }$$\end{document} , and Tc, with the error for each property weighted equally, as these are the choices yielding the best results. We used the default hyperparameters. Data for the training dataset were obtained from refs. ^9,10,66^, and can be downloaded from https://alexandria.icams.rub.de. The used training dataset ensures that all elements are included except the lanthanides and actinides, with La and Lu as the only exceptions. It is noted that the Tc expands a wide range up to 80 K, providing some coverage of high-Tc regimes. In addition, all crystal systems except monoclinic and triclinic are represented in the training set, including cubic hydride systems, which are relevant for our target materials in this study.

The other ab initio calculations were all carried out by using Quantum ESPRESSO^61,62^. In the high-throughput Quantum Espresso calculations, geometry optimizations employed uniform Γ-centered k-point grids with a density of 3000 k-points per reciprocal atom. When this resulted in an odd number of k-points in any direction, the next even number was used. Energy, force, and stress convergence thresholds were set to 1 × 10^−8^ a.u., 1 × 10^−6^ a.u., and 5 × 10^−2^ kbar, respectively. For electron–phonon coupling calculations, we utilized a double-grid technique where the k-grid from lattice optimization was doubled in each direction for the coarse grid and quadrupled for the fine grid. Phonon q-sampling used half of the aforementioned k-point grid. The Eliashberg function was obtained through double δ-integration using a Methfessel-Paxton smearing of 0.05 Ry.

Ab initio Eliashberg simulations including Coulomb interactions are performed using the approach of ref. ^51^, where Coulomb matrix elements and average Coulomb interaction functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W(\xi ,\xi ^{\prime} )$$\end{document} are computed using the elk code^67^. SCDFT anisotropic simulations are done with a Monte-Carlo algorithm using 35 K k-points accumulated with higher probability near the Fermi surface^26,68^. The same algorithm is used to compute the electron-phonon spectral functions for multiband simulations. Electronic and coupling parameters are linearly interpolated on the random mesh, which is used for the solution of the SCDFT Kohn-Sham gap equation. We adopt the function from ref. ^55^. Fermi surfaces are plotted with the Fermisurfer code^49^.

Supplementary information

Supplementary Information