Searching High Temperature Superconductors with the assistance of Graph
  Neural Networks
Predicting high temperature superconductors has long been a great challenge.
A major difficulty is how to predict the transition temperature Tc of
superconductors. Recently, progress in material informatics has led to a number
of machine learning models predicting Tc, which greatly improves the efficiency
of prediction. Unfortunately, prevailing models have not shown adequate
physical rationality and generalization ability to find new high temperature
superconductors, yet. In this work, in order to give a trustable prediction on
the unexplored materials, we built a bond-sensitive graph neural network
(BSGNN), which is optimized to process the information of chemical bond and
electron interaction in the crystal lattice, to predict the Tc maximum of each
type of superconducting materials. On the basis of the domain knowledge
considered in the data preparation and algorithm design, our model revealed a
relevance between the Tc-Tc maximum and chemical bonds. The results indicate
that shorter bond length is favored by high Tc, which is in accordance with
previous human experience. Moreover, it also shows that some specific chemical
elements are favored by high Tc, which is beyond what human experts already
knew. It gives a convenient guidance for searching high temperature
superconductors in materials database, by ruling out the materials that could
never have high Tc.
Liang Gu, Yang Liu, Pin Chen, Haiyou Huang, Ning Chen, Yang Li, Yutong Lu, Yanjing Su
# Predicting Transition Temperature of Superconductors with Graph Neural Networks

###### Abstract

Predicting high temperature superconductors has long been a great challenge. The difficulty lies in how to predict the transition temperature (\(\mathbf{T_{c}}\)) of superconductors. Although recent progress in material informatics has led to a number of machine learning models predicting \(\mathbf{T_{c}}\), prevailing models have not shown adequate generalization ability and physical rationality to find new high temperature superconductors, yet. In this work, a bond sensitive graph neural network (BSGNN) was developed to predict the \(\mathbf{T_{c}}\) of various superconductors. In BSGNN, communicative message passing and graph attention methods were utilized to enhance the model's ability to process bonding and interaction information in the crystal lattice, which is crucial for the superconductivity. Consequently, our results revealed the relevance between chemical bond attributes and \(\mathbf{T_{c}}\). It indicates that shorter bond length is favored by high \(\mathbf{T_{c}}\). Meanwhile, some specific chemical elements that have relatively large van der Waals radius is favored by high \(\mathbf{T_{c}}\). It gives a convenient guidance for searching high temperature superconductors in materials database, by ruling out the materials that could never have high \(\mathbf{T_{c}}\).

**Keywords: materials informatics, machine learning, graph neural network, superconductivity, superconductors, transition temperature**
## 1 Introduction

As the practical application of superconductors is severely constrained by the low working temperature, hunting for new superconductors with higher transition temperature (\(T_{c}\)) is a long-cherished dream for generations. However, the discovery of new high temperature superconductors (HTS) is still a significant challenge. [1, 2, 3, 4, 5, 6, 7] The \(T_{c}\) of superconductors varies with a good number of interrelated factors [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23], leading to the difficulty in predicting \(T_{c}\) of unexplored materials.

In this instance, machine learning (ML) has recently been exploited in the research on superconductors with the hope of giving guidance and inspiration to find new HTS. In the latest few years, plenty of ML models were developed to predict the \(T_{c}\) of superconductors or the existence of superconductivity. [24, 25, 26, 27, 28, 29, 30, 31, 32] For example, models predicting \(T_{c}\) were trained with 10000+ data mainly derived from SuperCon database [33]. Several algorithms, including random forest (RF) [24], atom table convolutional neural network (ATCNN) [27, 30] and convolutional gradient boosting decision tree (ConvGBDT) [34], have achieved predictive scores as good as \(R^{2}\sim 0.9\). Then the ML models were used to identify HTS from the materials in open access databases such as Inorganic Crystallographic Structure Database (ICSD) [35] and Crystallography Open Database (COD) [36]. Unfortunately, high scoring models so far gave inconsistent results when proposing candidates of HTS, and most of those candidates have not been verified experimentally, yet. Therefore, it calls for more new advanced algorithms and models with better generalization ability and physical rationality.

In this study, we developed a bond sensitive graph neural network (BSGNN) algorithm with optimized graph structures and message passing method, so as to adapt for the problem of superconductivity. Being good at processing structure information and interaction information, graph neural network (GNN) models for general purpose in materials science have shown impressive performance in predicting some material properties [37, 38, 39, 40]. There are also GNN models predicting conventional superconductivity [31], but is a lack of models powerful enough to predict new HTS. Despite of the abundance of architectures of GNN algorithm [41], it is unclear yet that which architectures work well on superconductivity. Here, according to the domain knowledge about superconductivity, we focused on the information of bonding and interaction in the crystal lattice, which is known as crucial factors influencing the superconductivity. Three methods, called nearest-neighbors-only graph representation (NGR), communicative message passing (CMP) and graph attention (GAT), were integrated into our GNN algorithm. Regression models to predict the \(T_{c}\) of superconductors were trained. With the help of BSGNN, we explored the entire ICSD to find potential HTS materials. Then, the model was used to predict some compressed hydrides as well as a series of conceived binary compounds in different lattice constants. The results showed that shorter bond length of specific atom pairs in the lattice is essential for HTS, implying that our model has learned the correlation between the \(T_{c}\) and pressure. Our work suggests that GNN model can have better interpretability and generalization ability in predicting superconductors and provide substantial help to materials scientists.

## 2 Results and discussion

### Modeling

It is well known that the crystalline structure is crucial for the superconductivity. How to deal with the information of interaction between atoms or electrons in the crystal lattice, or to say, the bond information, is vital for a ML model predicting superconductors. So, three modules in our GNN code to capture the bond information. As shown in Fig. 1, for nearest-neighbors-only graph representation (NGR) module, which encodes the crystalline structures into crystal graphs, only the chemical bonds of closest neighboring atoms were considered. The nodes in the graph are each atom, while the edges in the graph are the chemical bonds of closest neighbors. There are 16 features (various atomic attributes) on each node and one feature (normalized bond length) on each edge (see more details in section \(Methods\)). In communicative message passing (CMP) module [42], the information exchange between nodes and adjacent edges are considered. And the graph attention (GAT) module [43]. All those three modules focused on the local interplay of chemical bonds.

Regression models predicting \(T_{c}\) (actually \(log_{2}T_{c}\)) were trained and tested with the split datasets. Figure 2 shows the result of the best one in five models. Those five models achieved an average predictive score of \(R^{2}=0.85\pm 0.05\).

Table 1 shows the prediction on three superconductors found in recent few years, which are not in the input data. It can be seen that our GNN model has a good performance, not only in predicting various superconducting materials, but also in predicting unexplored materials (unseen data to the model).
\begin{table}
\begin{tabular}{l l l l} \hline No. & material & \(T_{c}^{exp}\) (K) & \(T_{c}^{pred}\) (K) \\ \hline
1 & CaH\({}_{6}\) @172 GPa & 215 & 249 \\
2 & Ti @248 GPa & 26 & 16 \\
3 & CsV\({}_{3}\)Sb\({}_{5}\) & 2.3 & 11 \\ \hline \end{tabular}
\end{table}
Table 1: Three superconductors discovered recently (not in train data).

Figure 1: Three critical modules in BSGNN: nearest-neighbors-only graph representation (NGR) module, communicative message passing (CMP) module, and graph attention (GAT) module. \(v_{n}\) and \(e_{uv}\) represent the nodes and edges, respectively. \(h_{n}\) and \(h_{uv}\) represent the node features and edge feature, respectively. \(\alpha_{uv}\) is the attention coefficients. \(n\), \(u\), \(v\) are indices of nodes.
### Screening high temperature superconductors

The model shown in Fig. 2 was used to predict the materials in ICSD. The ICSD contains over \(200,000+\) entries of ordered or disordered crystal structures. Here we take all \(110,000+\) ordered structures as the data to predict. Table 2 shows the top 10 (higher \(T_{c}\)) candidates of HTS, proposed by our GNN model.

Unexpectedly, it can be found in Tab. 2, Fig. 3, and \(Additional\)\(information\) (candidates.csv) that nearly half of the predicted \(T_{c}\) is higher than 30 K, which is inconsistent with the domain knowledge. Of course, it doesn't mean our predictions are unreliable. It just means that our model can not directly predict new high temperature superconductors, because the regression model can only predict how high the \(T_{c}\) could be, but can not predict whether or not the material is a superconductor. More specifically, the model often mistakes insulators as high temperature superconductors. But, many materials in ICSD are not superconductors, even not metals, but insulators. The predicted \(T_{c}\) of insulators should be seen as the \(T_{c}\) they could have, if they could be changed into metals by chemical doping or high pressure. Unfortunately, most insulators would not change into metals, no matter how they are doped or compressed. They may have met the necessary conditions for high \(T_{c}\), but just

\begin{table}
\begin{tabular}{l l l l} \hline No. & material & \(T_{c}^{pred}\) (K) & category \\ \hline
1 & CF\({}_{4}\) & 271 & molecular crystal \\
2 & CaZnF\({}_{4}\) & 270 & ionic crystal \\
3 & C\({}_{2}\)F\({}_{6}\)O\({}_{3}\) & 270 & molecular crystal \\
4 & NH\({}_{4}\)F & 269 & ionic crystal \\
5 & MgF\({}_{2}\) & 265 & ionic crystal \\
6 & LiCaF\({}_{3}\) & 262 & ionic crystal \\
7 & AlF\({}_{3}\) & 254 & ionic crystal \\
8 & CaO\({}_{2}\) & 254 & ionic crystal \\
9 & CaF\({}_{2}\) & 250 & ionic crystal \\
10 & CaCO\({}_{3}\) & 248 & ionic crystal \\ \hline \end{tabular}
\end{table}
Table 2: Top 10 HTS candidates (highest \(T_{c}^{pred}\)), proposed by our GNN model.

Figure 2: Model performance. Colored points: test data (randomly 10% data).
don't have the chance to be superconductors. So, we need a materials screening to rule out the unchangeable insulators (mostly ionic crystals and molecular crystals).

Here, we used the data of energy gap (given by DFT calculation) in MATGEN [42] to identify insulators. The materials \(E_{g}>\) 2 eV were removed from the candidates list. Then, we further screened the insulators in the HTS candidate list manually. The ionic crystals and molecular crystals were removed from the candidates list. Table 3 shows the top 10 candidates after manual screening.

### Understanding the GNN model from the predictions

In order to understand what the model have learn, we looked for the patterns in the predictions. For that purpose, there is no need to box ourselves into the materials experimentally observed. So, we conceived a series of binary compounds, with the sphalerite structure (space group: F\(\bar{4}\)3m, No.216), and in different lattice parameters. We choose the structure of sphalerite because in this structure, there is only one type of chemical bond (for example, Ga-As bond is the only one in GaAs, and the Ga-Ga and As-As is too far to form a substantial bond). Then the \(T_{c}\) of these crystals of binary compound (CBC), if exists, could be considered to be solely contributed by that bond. They are purest systems to see the influence of chemical bond on the superconductivity.

Most of these CBC are not superconducting, or not metallic, or even thermodynamically unstable. But by predicting these materials don't exist, we can see the influence of bond length and elements to the \(T_{c}\). The meaning of the predictions of these CBC is how high their \(T_{c}\) maxima could be, if they were superconductors.

Figure 4a shows a result of the predictions of CBC. Each \(T_{c}\) value in Fig. 4a represents the theoretical \(T_{c}\) maxima of the CBC consisting of each two chemical elements in the heat map. It can be found that the \(T_{c}\) varies with the combination of elements. The chemical elements of the two atoms forming the bond (not each single

Figure 3: Statistical chart of the predictions on the materials in ICSD.

\begin{table}
\begin{tabular}{l l l l} \hline No. & material & \(T_{c}^{pred}\) (K) & ICSD code \\ \hline
1 & CaB\({}_{6}\)O\({}_{10}\) & 154 & 161320 \\
2 & Ca\({}_{2}\)Ga\({}_{2}\)O\({}_{5}\) & 145 & 051545 \\
3 & CaC\({}_{2}\)Bi\({}_{2}\)O\({}_{8}\) & 138 & 094741 \\
4 & Co\({}_{0.89}\)Ca\({}_{1.11}\)GeO\({}_{4}\) & 136 & 173465 \\
5 & Sn\({}_{0.976}\)Au\({}_{0.024}\)F\({}_{4}\) & 132 & 078909 \\
6 & Ca\({}_{2}\)B\({}_{5}\)O\({}_{9}\)Br & 126 & 018001 \\
7 & Ca\({}_{2}\)Nb\({}_{2}\)O\({}_{6}\)F & 121 & 063189 \\
8 & RbNi\({}_{2}\)F\({}_{6}\) & 113 & 031781 \\
9 & LiB\({}_{6}\)O\({}_{9}\)F & 111 & 420286 \\
10 & YCaGa\({}_{3}\)O\({}_{7}\) & 103 & 109448 \\ \hline \end{tabular}
\end{table}
Table 3: Top 10 HTS candidates after manual screening (may need further doping).
element, but the element combination) is vital for high \(T_{c}\). Some element combinations are able to support high \(T_{c}\), whereas some others never had a chance.

Meanwhile, Fig. 4b shows that shorter bonds usually lead to higher \(T_{c}\), which is consistent with the common sense, and indicating that our model do have learnt the dependence of \(T_{c}\) on bond length. On the basis of that, we can predict the \(T_{c}\) of superconductors at high pressure, with confidence.

Please note that what is shown in the heat-map (Fig. 4a) is not the contribution of each element to the \(T_{c}\), but of each combination of two elements. Here, it's inappropriate to discuss either one element alone. However, still, we can see the differentiation between elements by calculate the maximal or average value of each row or column in the heat-map, as shown in Fig. 4c. It suggests that some elements are more favorable than others for higher \(T_{c}\). Especially, Alkali, Alkaline earth, Chalcogen, and Halogen elements are preferred.

According to the domain knowledge of superconductivity, the factors affecting \(T_{c}\) are numerous. Bond length and elements combination are only two of them. Having learnt the influence of bond length, our model has surpassed all previous ML models, but still not goes beyond human expertise. Meanwhile, the other knowledge that our model has learnt, the influence of elements combination, is a new information that has never been noticed and interpreted.

Figure 4: Results of predicting conceived binary compounds. (a) Each \(T_{c}^{pred}\) value represents the predicted \(T_{c}\) of the conceived binary compounds consisting of two chemical elements (marked as “element 1” and “element 2”). And the bond length between the two nearest neighboring atoms is 50% of the sum of van der Waals (vdW) radii (\(d=r_{1}^{vdW}+r_{2}^{vdW}\)). (b) Dependence of \(T_{c}\) on bond length. The three data points in each line are 50%, 70%, and 90% of \(d\), respectively. (c) Black and red: the average / maximal predicted \(T_{c}\) for each element, derived from subplot (a), blue: the vdW radius of each element.
### Hunting for superior superconductors in compressed hydrides

In addition, since our model has learnt the variation of \(T_{c}\) with bond length, it is supposed to get good performance in predicting superconductors at high pressure. So we used the model shown in Fig. 2 to predict compressed hydrides. Figure 5 shows the results of several typical hydrides, in which YH\({}_{6}\), LaH\({}_{10}\), H\({}_{3}\)S and H\({}_{2}\)S are in train data, whereas CeH\({}_{9}\) and CaH\({}_{6}\) are unseen data. On this basis, our model can be used to search new compressed hydrides with higher \(T_{c}\) or lower pressure.

## 3 Conclusion

In this study, a bond sensitive GNN model (BSGNN) was developed to predict the \(T_{c}\) of various superconducting materials. Our work suggests that a well-designed GNN algorithm can achieve rational and reliable results on complex materials problems with unclear mechanism, such as superconductivity. Owing to the GAT and CMP techniques, BSGNN has effectively learnt the dependence of \(T_{c}\) on both bond length and chemical composition. The predictions of BSGNN shows that shorter bond length, as well as specific elements with relatively larger vdW radii, is in favor of high \(T_{c}\). It illustrates the necessity and importance of considering crystalline structure information when predicting superconductors. As the influence of bond length was learnt, BSGNN shall have advantages in predicting superconducting materials under different pressure. With the help of BSGNN, we predicted all the materials in ICSD, and proposed some promising candidates of HTS after further artificial screening.

## 4 Methods

### Data preparation

Our GNN models were trained with the data of 612 superconductors. The crystalline structures of those superconductors are derived from ICSD [35] (up to 2019), and the \(T_{c}\) values of them are mainly derived from SuperCon [33] (up to 2014). Some compressed hydrides and some superconductors discovered recently were added into the dataset-SuperCon. The intersection set of dataset-SuperCon and dataset-ICSD has 1526 materials. In those materials, we took the data that \(T_{c}\)\(>\) 5 K (612 entries) as the input data for our GNN model.

The 612 input data were randomly (but not fully random) divided into train/test sets in a proportion of 9:1. Before the dataset splitting, the materials family of each data was labeled. During the dataset splitting, the ratio

Figure 5: Variation of predicted \(T_{c}\) with pressure for several typical compressed hydrides. Solid filled points: experimental values, light filled points: predicted values, open filled points: DFT calculated values.
of each materials family was maintained for both the train set and test set. As a result, the data distribution is better than the case of complete random data splitting.

There are one edge feature (bond length) and 16 node features in our crystal graph (as shown Tab. 4).

### Algorithm

**NGR module**

The strategy of the graph representation in our BSGNN algorithm is called nearest-neighbors-only. Each node in the graph is an atom in the crystal lattice. Each edge in the graph is a chemical bond between nearest neighbors. Here, "nearest neighbors" means a set of neighboring atoms for a centering atom. There is an edge between two neighboring atoms when \(d_{neighbor}<1.2\) * \(d_{nearest}\), where \(d_{neighbor}\) is the distance between those two neighboring atoms, \(d_{nearest}\) is the distance to the nearest neighbor for either of those two atoms.

**CMP module**

Inspired by the work on molecular graph with enhanced interatomic interactions [45], we apply the CMPNN framework to perform message passing on the aforementioned crystal graphs. The message generation and message update are iterated through the following equation.

\[\left\{\begin{aligned} m_{v}^{t+1}&=\sum_{k\in N(v)}h _{kv}^{t}\\ h_{v}^{t+1}&=h_{v}^{t}+m_{v}^{t+1}\\ m_{vw}^{t+1}&=h_{v}^{t+1}-h_{wv}^{t}\\ h_{vw}^{t+1}&=ReLU(h_{vw}^{0}+W_{m}m_{vw}^{t+1}) \end{aligned}\right.\]

Where \(h_{kv}^{t}\) represents the directed edge vector pointing to node \(v\) at time \(t\), which generates the message \(m_{v}^{t+1}\) for node \(v\) at time \(t+1\). Then, the message is used to update the hidden layer representation of node \(v\) at time \(t+1\). Next, the updated node information can be used to update all the directed edge vectors in reverse. Specifically,

\begin{table}
\begin{tabular}{l l} No. & feature name \\ \hline
1 & atomic mass 1  \\
2 & vdW radius 1  \\
3 & atomic radius 1  \\
4 & electronegativity (Pauling) 2  \\
5 & orbital energy of highest valence electrons 3  \\
6 & orbital energy of lowest valence electrons 4  \\
7 & number of unfilled valence orbitals \\
8 & number of filled valence orbitals \\
9 & number of unfilled valence \(s\) orbitals \\
10 & number of filled valence \(s\) orbitals \\
11 & number of unfilled valence \(p\) orbitals \\
12 & number of filled valence \(p\) orbitals \\
13 & number of unfilled valence \(d\) orbitals \\
14 & number of filled valence \(d\) orbitals \\
15 & number of unfilled valence \(f\) orbitals \\
16 & number of filled valence \(f\) orbitals \\ \end{tabular}
\end{table}
Table 4: Node features.
the message \(m_{wv}^{t+1}\) generated by each directed edge is obtained by subtracting the vector \(h_{wv}^{t}\) of its reverse edge from the starting node vector \(h_{v}^{t+1}\). The model prevents gradient vanishing by using residual connections between the generated edge messages and the initial edge vectors. Such completes the update of all edge vectors, where \(W_{m}\) represents the learnable weight matrix.

#### Attention module

The attention layer consists of local attention layer and global attention layer. The local attention layer helps the central atom identify the neighbor atoms with greater effect as follows [46]:

\[\boldsymbol{h}^{{}^{\prime}}(u)=\boldsymbol{W}_{1}\boldsymbol{h}(u)+\sum_{v \in N(u)}\alpha_{u,v}\boldsymbol{W}_{2}\boldsymbol{h}(v) \tag{1}\]

where the attention coefficients \(\alpha_{u,v}\) are computed as:

\[\alpha_{u,v}=softmax(\frac{(\boldsymbol{W}_{3}\boldsymbol{h}(u))^{T}( \boldsymbol{W}_{4}\boldsymbol{h}(v))}{\sqrt{d}}) \tag{2}\]

while the global attention layer is used to identify which atom of the material contributes more to \(T_{c}\) as follows [47]:

\[\boldsymbol{z}=\sum_{u\in U}Set2Set(\boldsymbol{h}^{{}^{\prime}}(u)) \tag{3}\]

#### Additional information

Data and code can be found at [https://github.com/GLinustb/BSGNN](https://github.com/GLinustb/BSGNN).

#### Author contributions

The study was planned and designed by YL, PC and HH, machine learning programs were performed by LG and PC, and checked by YL, the manuscript prepared by LG, YL, HH and PC. All authors discussed the results and commented on the manuscript.

#### Competing interests

The authors declare no competing interests.

#### Acknowledgements

The authors gratefully acknowledge the financial support of Guangdong Province Key Area R&D Program (2019B010940001).
