A Data-Driven Statistical Model for Predicting the Critical Temperature of a Superconductor

TL;DR

This paper presents a statistical model that predicts the critical temperature of superconductors using features derived from their chemical formulas, achieving reasonable accuracy but not classifying superconductivity.

Contribution

The study introduces a data-driven approach to estimate superconducting critical temperatures based on chemical features, highlighting key factors influencing T_c.

Findings

Root-mean-squared-error of ±9.5 K in predictions

Features like thermal conductivity and atomic radius are most influential

Model predicts T_c for superconductors but not their superconducting status

Abstract

We estimate a statistical model to predict the superconducting critical temperature based on the features extracted from the superconductor's chemical formula. The statistical model gives reasonable out-of-sample predictions: $\pm 9.5$ K based on root-mean-squared-error. Features extracted based on thermal conductivity, atomic radius, valence, electron affinity, and atomic mass contribute the most to the model's predictive accuracy. It is crucial to note that our model does not predict whether a material is a superconductor or not, it only gives predictions for superconductors.

Figures

Tables

Table 1. Table 1: This table shows the properties of an element which are used for creating features to predict T c subscript 𝑇 𝑐 T_{c} .

Variable	Units	Description
Atomic Mass	atomic mass units (AMU)	total proton and neutron rest masses
First Ionization Energy	kilo-Joules per mole (kJ/mol)	energy required to remove a valence electron
Atomic Radius	picometer (pm)	calculated atomic radius
Density	kilograms per meters cubed (kg/m3)	density at standard temperature and pressure
Electron Affinity	kilo-Joules per mole (kJ/mol)	energy required to add an electron to a neutral atom
Fusion Heat	kilo-Joules per mole (kJ/mol)	energy to change from solid to liquid without temperature change
Thermal Conductivity	watts per meter-Kelvin (W/(m $\times$ K))	thermal conductivity coefficient $\kappa$
Valence	no units	typical number of chemical bonds formed by the element

Table 2. Table 2: This table summarizes the procedure for feature extraction from material’s chemical formula. The last column serves as an example; features based on thermal conductivities for Re 7 ​ Zr 1 subscript Re 7 subscript Zr 1 \mathrm{Re_{7}Zr_{1}} are derived and reported to two decimal places. Rhenium and Zirconium’s thermal conductivity coefficients are t 1 = 48 subscript 𝑡 1 48 t_{1}=48 and t 2 = 23 subscript 𝑡 2 23 t_{2}=23 W/(m × \times K) respectively. Here: p 1 = 6 7 , p 2 = 1 7 , w 1 = 48 71 , w 2 = 23 71 , A = p 1 ​ w 1 p 1 ​ w 1 + p 2 ​ w 2 ≈ 0.926 , B = p 2 ​ w 2 p 1 ​ w 1 + p 2 ​ w 2 ≈ 0.074 . formulae-sequence formulae-sequence subscript 𝑝 1 6 7 formulae-sequence subscript 𝑝 2 1 7 formulae-sequence subscript 𝑤 1 48 71 formulae-sequence subscript 𝑤 2 23 71 𝐴 subscript 𝑝 1 subscript 𝑤 1 subscript 𝑝 1 subscript 𝑤 1 subscript 𝑝 2 subscript 𝑤 2 0.926 𝐵 subscript 𝑝 2 subscript 𝑤 2 subscript 𝑝 1 subscript 𝑤 1 subscript 𝑝 2 subscript 𝑤 2 0.074 p_{1}=\frac{6}{7},p_{2}=\frac{1}{7},w_{1}=\frac{48}{71},w_{2}=\frac{23}{71},A=\frac{p_{1}w_{1}}{p_{1}w_{1}+p_{2}w_{2}}\approx 0.926,B=\frac{p_{2}w_{2}}{p_{1}w_{1}+p_{2}w_{2}}\approx 0.074.

Feature & Description	Formula	Sample Value
Mean	$=\mu =\left(t_1+t_2\right)/2$	$35.5$
Weighted mean	$=\nu =\left(p_1{t}_1\right)+\left(p_2{t}_2\right)$	$44.43$
Geometric mean	$={\left(t_1{t}_2\right)}^{1/2}$	$33.23$
Weighted geometric mean	$={\left(t_1\right)}^{p_1}{\left(t_2\right)}^{p_2}$	$43.21$
Entropy	$=-w_1\ln \left(w_1\right)-w_2\ln \left(w_2\right)$	$0.63$
Weighted entropy	$=-A\ln \left(A\right)-B\ln \left(B\right)$	$0.26$
Range	$=t_1-t_2\left(t_1>t_2\right)$	$25$
Weighted range	$=p_1{t}_1-p_2{t}_2$	$37.86$
Standard deviation	$={\left[\left(1/2\right)\left({\left(t_1-\mu \right)}^2+{\left(t_2-\mu \right)}^2\right)\right]}^{1/2}$	$12.5$
Weighted standard deviation	$=[p_1{(t_1-\nu )}^2+p_2{(t_2-\nu )}^2)]{}^{1/2}$	$8.75$

Table 3. Table 3: This table reports summary statistics on iron-based versus non-iron, and cuprate versus non-cuprate superconductors. The Size is the total number of observations of the material out of 21,263 materials. For example, 2,339 out of 21,263 materials contained iron. The rest of the columns report summary statistics for the observed critical temperatures (K): min = minimum, Q1 = first quartile, Median = median, Q3 = third quartile, Max = maximum, and SD = standard deviation.

	Size	Min	Q1	Median	Q3	Max	Mean	SD
Iron	2339	0.02	11.3	21.7	35.5	130.0	26.9	21.4
Non-Iron	18924	0.0002	4.8	19.6	68.0	185.0	35.4	35.4
Cuprate	10532	0.001	31.0	63.1	86.0	143	59.9	31.2
Non-Cuprate	10731	0.0002	2.5	5.7	12.2	185	9.5	10.7

Table 4. Table 4: This table reports the summary statistics for the critical temperatures values (K) of all 21,263 superconductors. The column headers are the min = minimum , Q1 = first quartile , median , Q3 = third quartile, Max = maximum, and SD = standard deviation of the superconducting critical temperatures (K).

Min	Q1	Median	Q3	Max	Mean	SD
0.00021	5.4	20	63	185.0	34.4	34.2

Table 5. Table 5: This figure shows the top 20 most important features based on the XGBoost gain. Here: wtd = weighted, gmean = geometric mean, std = standard deviation.

Feature	Gain
range_ThermalConductivity	0.295
wtd_std_ThermalConductivity	0.084
range_atomic_radius	0.072
wtd_gmean_ThermalConductivity	0.047
std_ThermalConductivity	0.042
wtd_entropy_Valence	0.038
wtd_std_ElectronAffinity	0.036
wtd_entropy_atomic_mass	0.025
wtd_mean_Valence	0.022
wtd_gmean_ElectronAffinity	0.021
wtd_range_ElectronAffinity	0.016
wtd_mean_ThermalConductivity	0.015
wtd_gmean_Valence	0.014
std_atomic_mass	0.013
std_Density	0.010
wtd_entropy_ThermalConductivity	0.010
wtd_range_ThermalConductivity	0.010
wtd_mean_atomic_mass	0.009
wtd_std_atomic_mass	0.009
gmean_Density	0.009

Table 6. Table 6: This table shows T c subscript 𝑇 𝑐 T_{c} predictions for a list of potential superconductors identified by Valentin et al. ( 2017 ) .

Material	Predicted $T_c$ (K)
$\mathrm{CsBe}\left({\mathrm{AsO}}_4\right)$	13.7
${\mathrm{RbAsO}}_2$	8.0
${\mathrm{KSbO}}_2$	10.2
${\mathrm{RbSbO}}_2$	11.8
${\mathrm{CsSbO}}_2$	10.1
${\mathrm{AgCrO}}_2$	53.3
$K_{0.8}\left({\mathrm{Li}}_{0.2}{\mathrm{Sn}}_{0.76}\right)O_2$	18.6
$\mathrm{Cs}\left(\mathrm{MoZn}\right)\left(O_3{F}_3\right)$	20.5
${\mathrm{Na}}_3{\mathrm{Cd}}_2\left({\mathrm{IrO}}_6\right)$	17.4
${\mathrm{Sr}}_3\mathrm{Cd}\left({\mathrm{PtO}}_6\right)$	12.8
${\mathrm{Sr}}_3\mathrm{Zn}\left({\mathrm{PtO}}_6\right)$	12.4
$\left({\mathrm{Ba}}_5{\mathrm{Br}}_2\right){\mathrm{Ru}}_2{O}_9$	17.0
${\mathrm{Ba}}_4\left({\mathrm{AgO}}_2\right)\left({\mathrm{AuO}}_4\right)$	56.7
${\mathrm{Sr}}_5{\left({\mathrm{AuO}}_4\right)}_2$	17.8
${\mathrm{RbSeO}}_{2}F$	16.7
${\mathrm{CsSeO}}_{2}F$	20.4
${\mathrm{KTeO}}_{2}F$	13.0
${\mathrm{Na}}_2{K}_4\left({\mathrm{Tl}}_2{O}_6\right)$	32.8
${\mathrm{Na}}_3{\mathrm{Ni}}_2{\mathrm{BiO}}_6$	17.1
${\mathrm{Na}}_3{\mathrm{Ca}}_2{\mathrm{BiO}}_6$	27.3
$\mathrm{CsCd}\left({\mathrm{BO}}_3\right)$	22.3
$K_2\mathrm{Cd}\left({\mathrm{SiO}}_4\right)$	17.7
${\mathrm{Rb}}_2\mathrm{Cd}\left({\mathrm{SiO}}_4\right)$	17.4
$K_2\mathrm{Zn}\left({\mathrm{SiO}}_4\right)$	19.6
$K_2\mathrm{Zn}\left({\mathrm{Si}}_2{O}_6\right)$	12.2
$K_2\mathrm{Zn}\left({\mathrm{GeO}}_4\right)$	17.6
$\left(K_{0.6}{\mathrm{Na}}_{1.4}\right)\mathrm{Zn}\left({\mathrm{GeO}}_4\right)$	25.6
$K_2\mathrm{Zn}\left({\mathrm{Ge}}_2{O}_6\right)$	10.4
${\mathrm{Na}}_6{\mathrm{Ca}}_3{\left({\mathrm{Ge}}_2{O}_6\right)}_3$	12.1
${\mathrm{Cs}}_3\left({\mathrm{AlGe}}_2{O}_7\right)$	14.8
$K_4\mathrm{Ba}\left({\mathrm{Ge}}_3{O}_9\right)$	15.1
$K_{16}{\mathrm{Sr}}_4{\left({\mathrm{Ge}}_3{O}_9\right)}_4$	13.5
$K_3\mathrm{Tb}\left[{\mathrm{Ge}}_3{O}_8{\left(\mathrm{OH}\right)}_2\right]$	11.2
$K_3\mathrm{Eu}\left[{\mathrm{Ge}}_3{O}_8{\left(\mathrm{OH}\right)}_2\right]$	11.3
${\mathrm{KBa}}_6{\mathrm{Zn}}_4\left({\mathrm{Ga}}_7{O}_{21}\right)$	30.1

Equations

$$p_1=\frac{6}{6+1}=\frac{6}{7},p_2=\frac{1}{6+1}=\frac{1}{7}.$$

$$w_1=\frac{t_1}{t_1+t_2}=\frac{48}{48+23}=\frac{48}{71},w_2=\frac{t_2}{t_1+t_2}=\frac{23}{48+23}=\frac{23}{71}.$$

$$A=\frac{p_1{w}_1}{p_1{w}_1+p_2{w}_2}\approx 0.926,B=\frac{p_2{w}_2}{p_1{w}_1+p_2{w}_2}\approx \mathrm{0.074.}$$

$$\text{out-of-sample mse}={\text{Average of (observed - predicted)}}^2$$

$${\hat{y}}_i=\sum _{k=1}^K{f}_k(x_i),$$

$$\text{Objective with respect to }f_{t}=\sum_{i=1}^{n}L\big{(}\underbrace{y_{i}}_{observed},\ \underbrace{\hat{y}^{(t-1)}_{i}+f_{t}(x_{i})}_{predicted}\big{)}+\Omega(f_{t}),$$

$$\Omega (f)=\gamma T+(1/2)\lambda \sum _{j=1}^T{w}_j^2,$$

$$\text{The Gain for the Feature}=\frac{\text{Sum of the Gains for the Feature}}{\text{Sum of the Gains for All the Features}}.$$

Full text

1 Introduction

Superconducting materials - materials that conduct current with zero resistance - have significant practical applications. Perhaps the best known application is in the Magnetic Resonance Imaging (MRI) systems widely employed by health care professionals for detailed internal body imaging. Other prominent applications include the superconducting coils used to maintain high magnetic fields in the Large Hadron Collider at CERN, where the existence of Higgs Boson was recently confirmed, and the extremely sensitive magnetic field measuring devices called SQUIDs (Superconducting Quantum Interference Devices). Furthermore, superconductors could revolutionize the energy industry as frictionless (zero resistance) superconducting wires and electrical system may transport and deliver electricity with no energy loss; see Hassenzahl (2000).

However, the wide spread applications of superconductors have been held back by two major issues: (1) A superconductor conducts current with zero resistance only at or below its superconducting critical temperature ($T_{c}$). Often impractically, a superconductor must be cooled to extremely low temperatures near or below the boiling temperature of nitrogen (77 K) before exhibiting the zero resistance property. (2) The scientific model and theory that predicts $T_{c}$ is an open problem which has been baffling the scientific community since the discovery of superconductivity in 1911 by Heike Kamerlingh Onnes, in Leiden.

In the absence of any theory-based prediction models, simple empirical rules based on experimental results have guided researchers in synthesizing superconducting materials for many years. For example, the eminent experimental physicist Matthias (1955) concluded that $T_{c}$ is related to the number of available valence electrons per atom. (A few of these rules came to be known as the Matthias’s rules.) It is now well known that many of the simple empirical rules are violated; see Conder (2016).

In this study, we take an entirely data-driven approach to create a statistical model that predicts $T_{c}$ based on its chemical formula. The superconductor data comes from the Superconducting Material Database maintained by Japan’s National Institute for Materials Science (NIMS) at http://supercon.nims.go.jp/index_en.html. After some data preprocessing, 21,263 superconductors are used.

To our knowledge, Valentin et al. (2017) and our work are the only papers that focus on statistical models to predict $T_{c}$ for a broad class of materials. However, Owolabi et al. (2014) and Owolabi and Olatunji (2015) focus on predicting $T_{c}$ for $\mathrm{Fe}$ and $\mathrm{MgB_{2}}$ based superconductors respectively.

We derive features (or predictors) based on the superconductor’s elemental properties that could be helpful in predicting $T_{c}$. For example, consider $\mathrm{Nb_{0.8}Pd_{0.2}}$ with $T_{c}=1.98$ K. We can derive a feature based on the average thermal conductivities of the elements. Niobium and palladium’s thermal conductivity coefficients are 54 and 71 W/(m$\times$K) respectively. The mean thermal conductivity is $(54+71)/2=62.5$ W/(m$\times$K). We can treat the mean thermal conductivity variable as a feature to predict $T_{c}$. In total, we define and extract 81 features from each superconductor.

We tried various statistical models but we eventually settled on two: A multiple regression model which serves as a benchmark model, and a gradient boosted model as the main prediction model which is implemented in our software.

Our software tool to predict $T_{c}$ and the associated data are available at https://github.com/khamidieh/predict_tc and will also be available at the publisher’s complementary site. We have done our best to make the software use and access to the data as easy as possible.

Gradient boosted models create an ensemble of trees to predict a response. The trees are added in a sequential manner to improve the model by accounting for the points which are difficult to predict. Once a gradient boosted model is fitted, the weighted average of all the trees is used to give a final prediction. Gradient boosted models predict well because they are able to account for the complex interactions and correlations among the features.

The boosted models were first developed by Schapire (1990) and Freund (1995). The boosted models were generalized to gradient boosting by Friedman (2001). We use the latest improvement called XGBoost (eXtreme Gradient Boosting) by Chen and Guestrin (2016), and the associated open-source R implementation of XGBoost by Chen et al. (2018a). XGBoost is also available in other popular programming languages such as python and Julia. The full source code is at https://github.com/dmlc/xgboost.

Anthony Goldbloom, CEO of Kaggle (now a Google company), the premier data competition site, stated: “It used to be random forest that was the big winner, but over the last six months a new algorithm called XGBoost has cropped up, and it’s winning practically every competition in the structured data category.” You can see the talk at https://www.youtube.com/watch?v=GTs5ZQ6XwUM. Outside the competition realm, XGBoost has been successfully applied in disease prediction by Chen et al. (2018b), and in quantitative structure activity relationships studies by Sheridan et al. (2016).

Our XGBoost model gives reasonable predictions: an out-of-sample error of about $9.5$ K based on root-mean-squared-error (rmse), and an out-of-sample $R^{2}$ values of about $0.92$. The numbers for the multiple regression model are about $17.6$ K and $0.74$ for the out-of-sample rmse and $R^{2}$ respectively. The multiple regression serves as a benchmark model.

We are able to assess the importance of the features in prediction accuracy. Features defined based on thermal conductivity, atomic radius, valence, electron affinity, and atomic mass are the most important features in predicting $T_{c}$. On the downside, simple conclusions such as the exact nature of the relationship between the features and $T_{c}$ can’t be inferred from the XGBoost model.

Valentin et al. (2017) also create a model to predict $T_{c}$. Our approach is different than Valentin et al. (2017) in the following ways: (1) We use XGBoost versus random forests, (2) we use a larger data set, (3) we use a single large model to obtain predictions rather than a cascade of models, (4) we create a larger number features only from the elemental properties, and (5) most importantly, we quantify the out-of-sample prediction error.

2 Data Preparation

This section describes the detailed steps for the data preparation and feature extraction. Subsection (2.1) describes how the element data is obtained and processed. Subsection (2.2) describes the data preparation from NIMS Superconducting Material Database. Subsection (2.3) details how the features are extracted.

2.1 Element Data Preparation

The element data with 46 variables and 86 rows (corresponding to 86 elements) are obtained by using the ElementData function from Mathematica Version 11.1 by Wolfram and Research (2017). Appendix (A) lists the information sources for the element properties used by ElementData. The first ionization energy data came from http://www.ptable.com/ and is merged with the Mathematica data. About 12% of the entries out of the 3956 ($=46\times 86$) entries are missing.

In choosing the properties, we are guided by Conder (2016) but we also use our judgement to pick certain properties. For example, we drop the boiling point variable, and instead use the fusion heat variable which has no missing values, and is highly correlated with the boiling point variable. We had also gained some experience and insight creating some initial models for predicting $T_{c}$ of elements only. We settle on 8 properties shown in table (1).

With the choice of the above variables, we are only missing the atomic radii of La and Ce; we replace them with their covalent radii since atomic radii and covalent radii have very high correlation ($\approx 0.95$) and approximately on the same scale and range. Some bias may be introduced into our data with this minor imputation. We add a small constant of 1.5 to the electron affinity values of all the elements to prevent issues when taking logarithm of [math].

2.2 Superconducting Material Data Preparation

Superconducting Material Database is supported by the NIMS, a public institution based in Japan. The database contains a large list of superconductors, their critical temperatures, and the source references mostly from journal articles. To our knowledge, this is the most comprehensive database of superconductors. Access to the database requires a login id and password but this is provided with a simple registration process.

We accessed the data on July 24, 2017 at http://supercon.nims.go.jp/supercon/material_menu. Once logged in, we chose “OXIDE & METALLIC” material. Figure (1) shows a screen shot of the menu. We clicked on the “search” button to get all the data. We obtained 31,611 rows of data in a comma separated file format. The key columns (variables) were “element”, the chemical formula of the material, and “Tc”, the critical temperature. Variable “num” was a unique identifier for each row. Column “refno” contained links to the referenced source. The next few steps describe the manual clean up process:

We remove columns “ma1” to “mj2”. 2. 2.

We sort the data by “Tc” from the highest to lowest. 3. 3.

The critical temperature for the following “num” variables are mistakenly shifted by one column to the right. We fix these by recording them under the “Tc” column: 31020, 31021, 31022, 31023, 31024, 31025, 153150, 153149, 42170, 42171, 30716, 30717, 30718, 30719,150001, 150002, 150003, 150004, 150005, 150006, 150007, 30712, 30713, 30714, 30715. 4. 4.

The following are removed since the critical temperatures seemed to have been misrecorded; They have critical temperatures over 203 K which as of July 2017 was the highest reliable recorded critical temperature. La0.23Th0.77Pb3 (num = 111620), Pb2C1Ag2O6 (num = 9632), Er1Ba2Cu3O7-X (num = 140) 5. 5.

All rows with “Tc” = 0 or missing are removed. 6. 6.

Columns with headings “nums”, “mo1”, “mo2”, “oz”, “str3”, “tcn”, “tcfig”, “refno” are removed. 7. 7.

We manually chang all materials with oxygen content formula such as O7-X to the best oxygen content approximation. For example, O7-X is changed to O7, O5+X is changed to O5, etc. This certainly introduces some error into our data but it is impossible to go document by document to get better estimates of the oxygen contents. At this point our data has two columns: “element” and “Tc”. 8. 8.

We use R statistical software by R Core Team (2017) and the CHNOSZ package by Dick (2008) to perform a preliminary check of the validity of the chemical formulas. The CHNOSZ package has a function makeup which reads the chemical formula in string format and breaks up the formula into the elements and their ratios. In some cases, it throws an error or a warning when the chemical formula does not make sense. For example it throws a warning message if Pb-2O is checked; Negative number of Pb does not make sense. However, the function does not check whether the material could actually exist. See figure (2) to get a sense of how this function works. With the help of the CHNOSZ package, we make the following modifications:

1. (a)

Yo975Yb0.025Ba2Cu3O, Yo975Yb0.025Ba2Cu3O, Yo975Yb0.025Ba2Cu3O are removed. There is no element with the symbol Yo. It’s likely that Y0.975 was misrecorded as Yo975 but we can’t be sure. 2. (b)

Bi1.7Pb0.3Sr2Ca1Cu2O0, La1.85Nd0Ca1.15Cu2O5.99, Bi0Mo0.33Cu2.67Sr2Y1O7.41,

Y0.5Yb0.5Ba2Sr0Cu3O7 are removed since some elements had coefficients of zero. 3. (c)

Y2C2Br0.5!1.5 is removed. The exclamation sign throws an error message. 4. (d)

Y1Ba2Cu3O6050 is removed. The coefficient of 6050 for oxygen is possibly a mistake. 5. (e)

Hg1234O10 is removed. The coefficient of 1234 for mercury is possibly a mistake. 6. (f)

Nd185Ce0.15Cu1O4 is removed. The coefficient of 185 for Neodymium is possibly a mistake. There is a Nd1.85Ce0.15Cu1O4 already in the data. 7. (g)

Bi1.6Pb0.4Sr2Cu3Ca2O1013 is changed to Bi1.6Pb0.4Sr2Cu3Ca2O10.13 since nearby rows in the data have formulas with O10.xx. 8. (h)

Y1Ba2Cu285Ni0.15O7 is changed to Y1Ba2Cu2.85Ni0.15O7 since nearby rows in the data have formulas with Cu2.xx. 9. 9.

The column headings of “Tc” and “element” are changed to “critical_temp” and “material” respectively.

6750 rows are left out because $T_{c}$ is either zero or missing. At this point we have 24,861 rows of data.

The rest of the data preparation is done in R Core Team (2017). We exclude any superconductor that has an element with an atomic number greater than 86. This eliminates an additional 973 rows of data. For example, superconductors that have uranium are left out. We remove the repeating rows. It would be impossible to manually check to see whether the repeated rows are genuine independent reports from independent experiments or they are just duplicate reportings. After all the data preparation and clean up, we end up with 21,263 rows of data or about 67% of the original data we started with.

2.3 Feature Extraction

In this section, we describe the feature extraction process through a detailed example: Consider $\mathrm{Re_{7}Zr_{1}}$ with $T_{c}=6.7$ K, and focus on the features extracted based on thermal conductivity.

Rhenium and Zirconium’s thermal conductivity coefficients are $t_{1}=48$ and $t_{2}=23$ W/(m$\times$K) respectively. The ratios of the elements in the material are used to define features:

[TABLE]

The fractions of total thermal conductivities are used as well:

[TABLE]

We need a couple of intermediate values based on equations (1) and (2):

[TABLE]

Once we have obtained the values $p_{1},p_{2},w_{1},w_{2},A,$ and $B$, we can extract 10 features from Rhenium and Zirconium’s thermal conductivities as shown in table (2).

We repeat the same process above with the 8 variables listed in table (1). For example, for features based on atomic mass, just replace $t_{1}$ and $t_{2}$ with the atomic masses of Rhenium and Zirconium respectively, then carry on with the calculations of $p_{1},p_{2},w_{1},w_{2},A,B$, and finally calculate the 10 features defined in table (2). This gives us $8\times 10=80$ features. One additional features, a numeric variable counting the number of elements in the supercondutor, is also extracted. We end up with 81 features in total.

In summary: We have data with 21,263 rows and 82 columns: 81 columns corresponding to the features extracted and 1 column of the observed $T_{c}$ values.

We also considered but did not implement features that simply indicate whether an element is present in the superconductor or not. For example, we could have had a column that indicated whether say oxygen is in the material or not. However, this approach would have added a large number of indicator variables to our data, made model selection and assessment too complicated, and increased the chances of over-fitting.

3 Analysis

This section has two parts: Basic summaries of the data are given in subsection (3.1). The statistical models are described in subsection (3.2).

3.1 Descriptive Analysis

Figure (3) shows the proportions of the superconductors that had each element. For example, Oxygen is present in about 56% of the superconductors. Copper, barium, strontium, and calcium are the next most abundant elements.

Iron-based superconductors and cuprates are of particular interest in many research groups so we report some summary statistics in table (3). Iron is present in approximately 11% of the superconductors. The mean $T_{c}$ of superconductors with iron is $26.9\pm 21.4$ K. The non-iron containing superconductors’ mean is $35.4\pm 35.4$ K; the mean and standard deviations happened to be the same. A t-distribution based 95% confidence interval suggests that iron containing superconductors’ mean $T_{c}$ is lower than the non-iron’s by 7.4 to 9.5 K. Cuprates comprise approximately 49.5% of the superconductors. The cuprates’ mean $T_{c}$ is $59.9\pm 31.2$ K. The non-cuprates’ mean $T_{c}$ is $9.5\pm 10.7$ K. A t-distribution based 95% confidence interval indicates that the cuprates’ mean $T_{c}$ is higher than the non-cuprates’ mean $T_{c}$ by 49.8 to 51.0 K.

Figure (4) shows the histogram of $T_{c}$ values. The values are right skewed with a bump around 80 K. Table (4) shows the summary statistics for $T_{c}$ values.

Figure (5) shows the mean $T_{c}$ grouped by elements. Mercury containing superconductors have the highest $T_{c}$ at around 80 K on average. However, this is not the full story. Figure (6) shows the standard deviation of $T_{c}$ grouped by elements. Although mercury containing superconductors have the highest $T_{c}$ on average, these same materials show the fourth highest variability in $T_{c}$. In fact, a plot of the mean $T_{c}$ versus the standard deviation of $T_{c}$ in figure (7) shows that on average the higher the mean $T_{c}$, the higher the variability in $T_{c}$ per element.

The average absolute value of the correlation among the features is 0.35. This indicates that the features are highly correlated. Motivated by this result, we attempted to reduce the dimensionality of the data using principal component analysis (PCA). However, our PCA analysis did not show any benefits in reducing the dimensionality since a large number of principal components were needed to capture a substantial percentage of the data variation; we abandoned the PCA approach.

3.2 Model Analysis

In this section we discuss the results of the multiple regression model, and the XGBoost model. We tried a few classical models including multiple regression with interactions, principal component regression, and partial least squares but none of these make any substantial improvements to the XGBoost model. We also tried random forests but they were too slow to tune given the data size and the number of features. Scalability and speed are important advantages of using XGBoost over random forests; See Chen and Guestrin (2016).

The prediction performance of the models are compared by using out-of-sample rmse. The out-of-sample rmse is estimated by the following cross validation procedure:

Out-Of-Sample RMSE Estimation Procedure:

At random, divide the data into $2/3$ train data and $1/3$ test data. 2. 2.

Fit the model using the train data. 3. 3.

Predict $T_{c}$ of the test data. 4. 4.

Obtain an estimate of the out-of-sample mean-squared-error (mse) by using the predictions from the last step and the observed $T_{c}$ values in the test data:

[TABLE] 5. 5.

Repeat steps 1 through 4, 25 times to collect 25 out-of-sample mse’s. 6. 6.

Take the mean of the 25 collected out-of-sample mse’s and report the square root of this average as the final estimate of the out-of-sample rmse.

3.2.1 The Multiple Regression Model

The multiple regression model’s out-of-sample rmse estimated by the procedure above is about 17.6 K. The out-of-sample $R^{2}$ is about $0.74$. Figure (8) shows the predicted $T_{c}$ versus the observed $T_{c}$ when we use all the data to fit the model. The line has an intercept of zero and a slope of 1. The plot indicates that the multiple regression model under-predicts $T_{c}$ of high temperature superconductors since many predicted points are below the line for the high temperature superconductors. The model over-predicts low temperature superconductors’ $T_{c}$. The multiple regression model simply serves as a benchmark model and should not be used for prediction. There would be no use in predicting $T_{c}$ using a sophisticated model such as XGBoost, if a commonly used multiple regression model does a good job. Here, the XGBoost model vastly improves the prediction accuracy.

3.2.2 The XGBoost Model

Before we go on, we give a brief description of XGBoost set up. XGBoost is described in detail in Chen and Guestrin (2016). A readable summary is given at https://xgboost.readthedocs.io/en/latest/model.html. Hastie et al. (2009) and Izenman (2008) give general overviews on boosting as well.

The functional form of XGBoost is:

[TABLE]

where $x_{i}$ is the $i$th input feature vector, $\hat{y}_{i}$ is the predicted response, and $f_{1},\ldots,f_{K}$ is a sequence of trees. The $t$-th tree $f_{t}$ is added by minimizing the following objective function:

[TABLE]

where $L$ is the desired loss function, $n$ is the total sample size, $y_{i}$’s are the response values, $\hat{y}^{(t-1)}_{i}$ is the $i$th predicted responses at the $t-1$ step, and $\Omega$ is a penalty function. The form of $\Omega$ is:

[TABLE]

where $T$ is the number of leaves in each tree, $w_{j}$’s are the leaf weights, and $\lambda$ and $\gamma$ are regularization parameters. The goal here is to add a new tree $f_{t}$ to the overall ensemble of trees to minimizes the loss between the observed and the predicted in equation (3), while preventing over-fitting by satisfying the penalty in equation (4). The addition of this penalty function to each tree in (4) is one major XGBoost differentiator from the established method by Friedman (2001). The penalty function appears to make a big difference in practice; see Chen and Guestrin (2016). Besides the clever penalty function, Chen and Guestrin (2016) implement numerous computational tricks to make their software scalable and very fast.

In addition to the penalty function, there are a number of tuning parameters that could reduce over-fitting and enhance the model’s prediction performance; They are mainly: (1) column subsampling which means only a fraction of the features are chosen at random at each stage of adding a new tree, (2) a learning parameter $0<\eta<1$ which scales the contribution of each new tree, (3) subsample ratio which means that XGBoost only uses a small percentage of the data to grow a new tree, (4) maximum depth of a tree, and (5) minimum child weight which is the minimum number of data points needed to be in each node.

To tune XGBoost, we first split the data at random to 2/3 train and 1/3 test data. Next, we create a grid - a grid contains all the possible combination of tuning parameters - with $\eta=0.010,0.015,0.020$, column subsampling = 0.25, 0.5, 0.75, subsample ratio = 0.5, minimum node size = 1, 10, and maximum depth of a tree = 15, 16, …, 24, 25. The total gird size is 198. This means that we need 198 different XGBoost models. For each model, 750 trees are grown. The rest of the XGBoost parameters are set to the default values. (This was not our only grid; we had done some experimentations with various grids before we decided to use this grid.) Finally, we evaluate the prediction accuracy of each model based on rmse at each tree = 1, 2, …, 749, 750.

The best model (with the lowest out-of-sample rmse) turn out to be: $\eta=0.02$, maximum depth $=16$, minimum child weight $=1$, column subsampling $=0.50$, and a tree size of 374. To obtain the final out-of-sample rmse and $R^{2}$, we follow the 6 step procedure outlined at the begining of section (3.2). The procedure yield an out-of-sample rmse of $9.5$ K, and a out-of-sample $R^{2}$ of 0.92. The out-of-sample rmse of $9.5$ K has a very important interpretation: On average, the tuned XGBoost model will be off by about $9.5$ K when predicting $T_{c}$.

Figure (9) shows the predicted $T_{c}$ versus the observed $T_{c}$. Except for lower observed $T_{c}$ vlues, no severe bias is discernable. There are are a number outliers visible.

3.2.3 Feature Importance

Feature importance in XGBoost is measured by gain. The gain for a feature is defined as follows: Whenever a tree is split on a feature, the improvement in the objective function is recorded. The gain for the feature is then:

[TABLE]

Features with higher gain are more important.

Table (5) shows the top 20 most important features. Features extracted based on thermal conductivity, atomic radius, valence, electron affinity, and atomic mass appear to be the most important features. Also observe that features defined based on thermal conductivity, valence, electron affinity, and atomic mass appear most often on the list. This may suggest that these properties could be more important than other properties in predicting $T_{c}$.

4 Prediction Software

We have put the code for prediction at https://github.com/khamidieh/predict_tc. The software is created using R Statistical programming language, R Core Team (2017). The data could also be directly downloaded from our github site.

We demonstrate some examples using the software. Figure (10) shows the predictions for three materials: $\mathrm{Ba_{0.2}La_{1.8}CuO_{4}}$, $\mathrm{MgB_{2}}$, and $\mathrm{Hg}$. The “verbose” option uses the cosine similarity measure to pull data with similar chemical formulas. The multiple entries for $\mathrm{Ba_{0.2}La_{1.8}CuO_{4}}$ are obtained. The default value for verbose is false so no superconductors similar to $\mathrm{MgB_{2}}$ and $\mathrm{Hg}$ are shown.

We had obtained the data on July 24, 2017. We like to see what sort of predictions we could obtain for some new superconductors reported since. Nishiyama et al. (2017) report a $T_{c}$ of around 3 K for $\mathrm{Ca_{0.5}Sr_{0.5}C_{6}}$. Goto et al. (2017) report a $T_{c}$ of 1.3 K for $\mathrm{NaSn_{2}As_{2}}$. Figure (11) shows the prediction results. The XGBoost model over-predicts but it is within the $\pm 9.5$ K out-of-sample rmse. The message “Not able to find match(es)” indicates that nothing in the training data is similar to these two new superconductors. We should not expect good predictions for completely new superconductors.

Figure (12) shows what can go wrong when the XGBoost model predicts badly or when the inputs do not make sense. The prediction for $\mathrm{H_{2}S}$, which has a $T_{c}$ of 203 K under extremely high pressures, is way off. (Note that $\mathrm{H_{2}S}$ with $T_{c}$ of 203 is not in the train data.) This is perhaps expected since there is no feature that captures the dependence of $T_{c}$ on pressure. The model gives a prediction for $\mathrm{FCl}$ but this is a non-sense; The prediction model can’t check for the existence of solids. The model gives an error message for $\mathrm{mgB{2}}$ since it does not recognize $\mathrm{mg}$ with the lower case m as an element.

Next, we predict $T_{c}$ for materials identified by Valentin et al. (2017) as potential superconductors. The results are shown in table (6). None of the superconductors in table (6) are found to be (cosine) similar to the superconductors in our train data.

5 Conclusion

We have shown that a statistical model using only the superconductors’ chemical formula can predict $T_{c}$ reasonably well. We have also made the software and the data easily available. There are practical uses for our model: (1) Researchers interested in finding high temperature superconductors may use the model to narrow their search, and (2) researchers could use the cleaned data along with new data (such as pressure or crystal structure) to make better models.

6 Acknowledgements

We like to thank Dr. Allan Macdonald, professor of physics at the University of Texas at Austin, for many useful suggestions.

7 Bibliography

Appendix A Mathematica ElementData

Below is the list of sources Mathematica has used to obtained the element property data. It is directly copied from:

http://reference.wolfram.com/language/note/ElementDataSourceInformation.html.

• •

Atomic Mass Data Center. “NUBASE.” 2003. http://amdc.in2p3.fr/web/nubase_en.html

• •

Cardarelli, F. Materials Handbook: A Concise Desktop Reference. Springer, 2000.

• •

Lide, D. R. (Ed.). CRC Handbook of Chemistry and Physics. 87th ed. CRC Press, 2006.

• •

Speight, J. Lange’s Handbook of Chemistry. McGraw-Hill, 2004.

• •

United Kingdom National Physical Laboratory. “Kaye and Laby Tables of Physical and Chemical Constants.” http://www.kayelaby.npl.co.uk/

• •

United States National Institute of Standards and Technology. “Atomic Weights and Isotopic Compositions Elements.”

https://www.nist.gov/pml/atomic-weights-and-isotopic-compositions-relative-atomic-masses

• •

United States National Institute of Standards and Technology. “NIST Chemistry Webbook.” http://webbook.nist.gov/chemistry/

• •

Winter, M. ”WebElements.” 2007. https://www.webelements.com/