Commit 42ba9b33 authored by Luca Massimiliano Ghiringhelli's avatar Luca Massimiliano Ghiringhelli
Browse files

Update tetradymite_PRM2020.ipynb

parent 535de8ff
......@@ -27,7 +27,7 @@
"created by: Luigi Sbailò, Thomas A. R. Purcell, Luca M. Ghiringhelli, and Matthias Scheffler \n",
" \n",
"Fritz Haber Institute of the Max Planck Society, Faradayweg 4-6, D-14195 Berlin, Germany <br>\n",
"<span class=\"nomad--last-updated\" data-version=\"v1.0.0\">[Last updated: Sep 29, 2020]</span>\n",
"<span class=\"nomad--last-updated\" data-version=\"v1.0.0\">[Last updated: Sep 15, 2020]</span>\n",
" \n",
"<div> \n",
"<img style=\"float: left;\" src=\"assets/tetradymite_PRM2020/Logo_MPG.png\" width=\"200\"> \n",
%% Cell type:markdown id: tags:
<div id="teaser" style=' background-position: right center; background-size: 00px; background-repeat: no-repeat;
padding-top: 20px;
padding-right: 10px;
padding-bottom: 170px;
padding-left: 10px;
border-bottom: 14px double #333;
border-top: 14px double #333;' >
<div style="text-align:center">
<b><font size="6.4">Discovery of new topological insulators in alloyed tetradymites
created by: Luigi Sbailò, Thomas A. R. Purcell, Luca M. Ghiringhelli, and Matthias Scheffler
Fritz Haber Institute of the Max Planck Society, Faradayweg 4-6, D-14195 Berlin, Germany <br>
<span class="nomad--last-updated" data-version="v1.0.0">[Last updated: Sep 29, 2020]</span>
<span class="nomad--last-updated" data-version="v1.0.0">[Last updated: Sep 15, 2020]</span>
<img style="float: left;" src="assets/tetradymite_PRM2020/Logo_MPG.png" width="200">
<img style="float: right;" src="assets/tetradymite_PRM2020/Logo_NOMAD.png" width="250">
%% Cell type:markdown id: tags:
### Introduction
This tutorial shows how to find descriptive parameters (short formulas) that predict whether alloyed materials are topological or trivial insulators, using the example of tetradymites. It is based on the algorithm sure independence screening and sparsifying operator (SISSO), that enables to search for optimal descriptor by scanning huge feature spaces.
<div style="padding: 1ex; margin-top: 1ex; margin-bottom: 1ex; border-style: dotted; border-width: 1pt; border-color: blue; border-radius: 3px;">R. Ouyang, S. Curtarolo, E. Ahmetcik, M. Scheffler, L. M. Ghiringhelli: <span style="font-style: italic;">SISSO: a compressed-sensing method for identifying the best low-dimensional descriptor in an immensity of offered candidates</span>, Phys. Rev. Materials 2, 083802 (2018) <a href="" target="_blank">[PDF]</a>.</div>
With the default settings, the method reproduces the same results from:
<div style="padding: 1ex; margin-top: 1ex; margin-bottom: 1ex; border-style: dotted; border-width: 1pt; border-color: blue; border-radius: 3px;">G. Cao, R. Ouyang, L. M. Ghiringhelli, M. Scheffler, H. Liu, C. Carbogno, and Z. Zhang: <span style="font-style: italic;">Artificial intelligence for high-throughput discovery of topological insulators: The example of alloyed tetradymites</span>, Phys. Rev. Materials 4, 034204 (2020) <a href="">[PDF]</a>,</div>
%% Cell type:markdown id: tags:
<div style="padding: 1ex; margin-top: 1ex; margin-bottom: 1ex; border-style: dotted; border-width: 1pt; border-color: blue; border-radius: 3px;"><b>Explanation of the method (click to expand/collapse)</b></div></summary>
We present a tool for predicting whether alloyed tetradymite are topological or trivial insulators, by using a set of descriptive parameters (a descriptor) based on free-atom data of the atomic species constituting the $AB-LMN$ materials, where $A,B \in \{ \textrm{As, Sb, Bi} \}$ and $L,M,N \in \{ \textrm{S, Se, Te} \}$. We apply a recently developed method: sure independence screening and sparsifying operator (SISSO), that allows to find an optimal descriptor in a huge feature space containing billions of features. In this tutorial an $\ell_0$-optimization is used as the sparsifying operator.
The method is described in:
<div style="padding: 1ex; margin-top: 1ex; margin-bottom: 1ex; border-style: dotted; border-width: 1pt; border-color: blue; border-radius: 3px;">
R. Ouyang, S. Curtarolo, E. Ahmetcik, M. Scheffler, L. M. Ghiringhelli: <span style="font-style: italic;">SISSO: a compressed-sensing method for identifying the best low-dimensional descriptor in an immensity of offered candidates</span>, Phys. Rev. Materials 2, 083802 (2018) <a href="" target="_blank">[PDF]</a>. <br> </div>
In this tutorial, we focus on the classification flavor of SISSO($\ell_0$).
In the space of descriptors, each category’s domain (here, topological vs trivial insulator) is approximated as
the region of space within the convex hull of the corresponding training data. SISSO finds the low-dimensional descriptor yielding the minimum overlap between these convex regions. In practice, the algorithm is iterative. At the first iteration, in the SIS step, it selects the $k$ features which yield the smallest overlap when convex regions (segments encompassing all the data in one category) over the training data are constructed. In the first iteration the feature giving the smalles overlap is already the 1D model. At each subsequent iteration $i$, in the SIS step. $k$ new features that do the same for those training points which were in the overlap regions at the previous steps (i.e., the residuals). Then, in the SO step, all $i$-tuples of features selected combining in all possible ways the features selected in the SIS steps are ranked by the size of the overlap. The $i$-tuple with the smallest overlap is the $i$D model.
In order to better identify a predictive model to classify unseen data point, at each dimension (iteration) a soft-margin support-vector machine <a href="" target="_blank">[C. Cortes & V. Vapnik, Machine learning 20, 273 (1995)]</a> is trained to define the separating hyperplanes. The resulting model is identified by the coefficents and intercept of the hyperplanes.
%% Cell type:markdown id: tags:
The idea demonstrated in this tutorial is to start from simple physical quantities ("primary features", here properties of the constituent free atoms such as Pauling electronegativity), to generate millions (or billions) of candidate formulas by applying arithmetic operations combining primary features. These candidate formulas constitute the so-called "feature space". Then, SISSO is used to select only a few of these formulas that explain the data.
By clicking directly on "Run" below, i.e., with the default selection, you can reproduce the 2D map as published in <a href="" target="_blank">PRM 2020</a>. You can also select primary features and allowed operations (by clicking the check-boxes), as well as the SISSO rung (i.e., the number of iterations in the construction of the feature space), the number of features that are selected at each iteration of the SIS step, and the max number of dimensions of the model. The materials considered here have up to 5 different atomic species in the unit cell, with the prototype formula $AB-LMN$, where the cations $A,B \in \{ \textrm{As, Sb, Bi} \}$ and the anions $L,M,N \in \{ \textrm{S, Se, Te} \}$. We have therefore grouped the features to be selected into those for cations and anions. This means that by selecting, e.g., a cation feature, such feature is added to the primary feature set for both $A$ and $B$ elements, but either is treated singularly in the feature construction and SISSO optimization. After the features' and other settings' selection, press "Run". \
After the results are shown for all models from one dimensional to the max chosen dimension, you can press "Plot interactive map" to reveal a map of tetradymites' topological vs trivial insulators, for the highest-dimensional model. If the highest-dimensional model is 2D, the support-vector-machine separation line between the two phases is shown. For higher dimensional models, the 3rd and 4th dimensions can be visualized via the size or the color of the data-point markers. Intuitive drop-down menus allow to assign axes, markers, and colors, to the descriptor components of choice.
With the selection of "PRM2020" (or default selection) as SISSO rung, a special feature space is uploaded, which contains much fewer features than in the production calculation used in <a href="" target="_blank">PRM 2020</a>. This allows to reobtain in the notebook the same result in a reasonsable time. Still, the provided feature space contains thousands of the top ranked features and SISSO finds the best nD model.
%% Cell type:code id: tags:
``` python
function code_toggle() {
if (code_show)
code_show = !code_show
$( document ).ready(code_toggle);
The Python code for this notebook is by default hidden for easier reading.
To toggle on/off the code, click <a href="javascript:code_toggle()">here</a>.
%%%% Output: display_data
%% Cell type:code id: tags:
``` python
from sissopp import Inputs, FeatureSpace, SISSOClassifier, FeatureNode, Unit
from sissopp.py_interface import read_csv
from sissopp.py_interface.import_dataframe import get_unit
from tetradymite_PRM2020.visualizer import Visualizer
import numpy as np
import pandas as pd
import os
%% Cell type:code id: tags:
``` python
# The dataset is stored in the NOMAD Archive and can be accessed with this query.
from nomad import client, config
config.client.url = ''
query = client.query_archive(query={
'dataset_id': ['BjT-NFK0QdOx81_z5TmyeQ']},
%% Cell type:code id: tags:
``` python
df_train = pd.read_pickle('./data/tetradymite_PRM2020/training_set')
%% Cell type:code id: tags:
``` python
# This piece of code is not run at initialization.
# It can create the molecular structures which are visualized.
path_structure = './data/tetradymite_PRM2020/structures/'
except OSError:
!rm ./data/tetradymite_PRM2020/structures/*
scale_factor = 10**10
alist = []
for compound in compounds:
for entry in range (1581):
labels = query[entry].section_run[0].section_system[-1].atom_labels
if (len(labels)>5):
labels_1 = str(labels[0])+'_'+str(labels[1])+'_'+str(labels[3])+'_'+str(labels[4])+'_'+str(labels[2])
labels_2 = str(labels[0])+'_'+str(labels[1])+'_'+str(labels[4])+'_'+str(labels[3])+'_'+str(labels[2])
labels_3 = str(labels[1])+'_'+str(labels[0])+'_'+str(labels[3])+'_'+str(labels[4])+'_'+str(labels[2])
labels_4 = str(labels[1])+'_'+str(labels[0])+'_'+str(labels[4])+'_'+str(labels[3])+'_'+str(labels[2])
if compound in list([labels_1, labels_2, labels_3, labels_4]):
n_atoms = len (labels)
lat_x, lat_y, lat_z = query[entry].section_run[0].section_system[-1].lattice_vectors.magnitude * scale_factor
file = open(path_structure + str(compound) +".xyz","w")
for i in [0,1,2]:
for j in [0,1,2]:
for k in [0,1,2]:
for n in range (n_atoms):
el = query[entry].section_run[0].section_system[-1].atom_labels[n]
xyz = query[entry].section_run[0].section_system[-1].atom_positions[n].magnitude * scale_factor
xyz += i*lat_x
xyz += j*lat_y
xyz += k*lat_z
file.write (el)
file.write ("\t%f\t%f\t%f\n"%(xyz[0],xyz[1],xyz[2]))
%% Cell type:code id: tags:
``` python
zeta = {'S':16, 'As':33, 'Se':34, 'Sb':51, 'Te':52, 'Bi':83}
chi = {'S':2.58, 'As':2.18, 'Se':2.55, 'Sb':2.05, 'Te':2.12, 'Bi':2.02}
lambd = {'S':0.05, 'As':0.19, 'Se':0.22, 'Sb':0.4, 'Te':0.49, 'Bi':1.25}
df_feat = pd.DataFrame(index=df_train.index, columns=[
for comp in df_train.index:
ablmn = comp.split('_')
df_feat.loc[comp] = pd.Series({
df_feat['Class'] = df_train['Class']
%% Cell type:code id: tags:
``` python
def get_feat_space_and_sisso_regressor(
selected_ops=["add", "abs_diff", "div", "sq", "exp"],
selected_features = 'all',
if default:
selected_ops = ["add", "sub", "mult", "div", "abs_diff", "sq", "cb", "sqrt", "cbrt", "inv", "abs"]
selected_features = 'all'
inputs = read_csv(
inputs = read_csv(
inputs.max_rung = max_rung
inputs.allowed_ops = selected_ops
inputs.n_sis_select = n_sis_select
inputs.n_dim = n_dim
inputs.n_residual = n_residual
inputs.n_model_store = 1
inputs.calc_type = "classification"
inputs.sample_ids_train = df_feat.index.tolist()
inputs.prop_train = df_feat["Class"].to_numpy()
inputs.prop_test = np.array([])
inputs.prop_label = "Class"
inputs.task_names = ["all_mats"]
feat_space = FeatureSpace(inputs)
sisso = SISSOClassifier(inputs, feat_space)
return feat_space, sisso
%% Cell type:code id: tags:
``` python
# In this cell interactions with buttons are defined
from ipywidgets import widgets, interactive
from IPython.display import HTML, clear_output
def handle_rung_selection(change):
if change['new'] == 'PRM2020':
default_operations = ['add', 'sub', 'abs_diff', 'mult', 'div', 'exp', 'neg_exp', 'inv', 'sq', 'cb',
'sqrt', 'cbrt', 'log', 'abs']
default_features = ['z_cations','x_cations','l_cations','z_anions','x_anions','l_anions']
for op, widget in zip(possible_operations, op_list):
widget.value = op in default_operations
widget.disabled = True
for feat, widget in zip(possible_features, feat_list):
widget.value = feat in default_features
widget.disabled = True
rung_selection.value = 'PRM2020'
feat_per_iter_selection.value = 50
dimension_selection.value = 2
for widget in op_list+feat_list:
widget.disabled = False
def plot_button_clicked(button):
with out2:
model = sisso.models[1][0]
compounds = df_train.index.to_list()
for feat in sisso.models[sisso.n_dim-1][0].feats:
classes = ['Topological insulators', 'Trivial insulators']
visualizer=Visualizer(df, sisso, classes)
def default_button_clicked(button):
rung_selection.value = 'PRM2020'
feat_per_iter_selection.value = 50
dimension_selection.value = 2
def run_button_clicked(button):
with out2:
with out1:
print('Calculating...', flush=True)
selected_features = []
allowed_operations = []
for op, widget in zip(possible_operations, op_list):
if widget.value:
for sel_feat, widget in zip(possible_features, feat_list):
if widget.value:
feat = sel_feat.split('_')[0]
typ = sel_feat.split('_')[1]
if (typ=='cations'):
selected_features.append(feat + '_'+ 'A')
selected_features.append(feat + '_'+ 'B')
if (typ=='anions'):
selected_features.append(feat + '_'+ 'L')
selected_features.append(feat + '_'+ "M")
selected_features.append(feat + '_'+ "N")
if rung_selection.value == 'PRM2020':
selected_features = "all"
tier = 0
default = True
tier = rung_selection.value
default = False
global feat_space
global sisso
feat_space, sisso = get_feat_space_and_sisso_regressor(
selected_ops = allowed_operations,
selected_features = selected_features,
max_rung = tier,
n_sis_select = feat_per_iter_selection.value,
n_dim = dimension_selection.value,
n_residual = 10,
default = default
if (dimension_selection.value>1):
print("Number of features generated: " + str(feat_space.n_feat))
for i in range(dimension_selection.value):
print(str(i+1)+'D model')
print("# misclassified: {} ".format(int(sisso.models[i][0].n_convex_overlap_train)))
string = "SVM dividing line: c0"
for nf, feat in enumerate(sisso.models[i][0].feats):
string = string + str(' + a'+str(nf)+'*'+str(feat.expr))
string = string + " = 0"
string = "c0:{:.4}".format(sisso.models[i][0].coefs[0][-1])
for j in range(i+1):
string = string + str(" | a"+str(j)+":{:.4}".format(sisso.models[i][0].coefs[0][j]))
print(string + '\n')
global df
except RuntimeError:
print("\nThe number of selected features per SIS iteration is bigger than the number of features available. Please reduce the number of selected features per SIS iteration (number of features generated / max number of dimensions) or increase the number of selected features and operations.")
print('The present selection does not lead to the creation of any derived features in the highest selected rung, please select at least one binary or power operator, or reduce the maximum rung')
%% Cell type:code id: tags:
``` python
cb_layout = widgets.Layout(width = '15px')
thin_layout = widgets.Layout(width = '100px')
mid_layout = widgets.Layout(width = '200px')
wide_layout = widgets.Layout(width = '300px')
possible_operations = ['add', 'sub', 'abs_diff', 'mult', 'div', 'exp', 'neg_exp', 'inv', 'sq', 'cb',
'sqrt', 'cbrt', 'log', 'abs']
possible_features = ['z_cations','x_cations','l_cations','z_anions','x_anions','l_anions']
tooltips = {
"z_cations" : "Atomic number",
"x_cations" : "Pauling electronegativity",
"l_cations" : "Spin orbit coupling",
"z_anions" : "Atomic number",
"x_anions" : "Pauling electronegativity",
"l_anions" : "Spin orbit coupling",
labels = {
'add' : '$x + y$', 'sub' : '$x - y$', 'abs_diff' : '$|x - y|$', 'mult' : '$x \cdot y$', 'div' : '$x / y$',
'exp' : '$\exp(x)$', 'neg_exp' : '$\exp(-x)$', 'inv' : '$1/x$', 'sq' : '$x^2$', 'cb' : '$x^3$',
'six_pow' : '$x^6$', 'sqrt' : '$\sqrt{x}$', 'cbrt' : '$\sqrt[3]{x}$', 'log' : '$\log(x)$',
'abs' : '$|x|$', 'sin' : '$\sin(x)$', 'cos' : '$\cos(x)$', 'z_cations' : '$Z_{cations}$', 'x_cations' : '$\chi_{cations}$',
'l_cations' : '$\lambda_{cations}$', 'z_anions' : '$Z_{anions}$', 'x_anions' : '$\chi_{anions}$', 'l_anions' : '$\lambda_{anions}$'
op_list = []
op_labels = []
feat_list = []
feat_labels = []
for operation in possible_operations:
op_list.append(widgets.Checkbox(description='', value=True, indent=False, layout=cb_layout))
for feature in possible_features:
feat_list.append(widgets.Checkbox(description=tooltips[feature], value=True, indent=False, layout=cb_layout))
op_box = widgets.VBox([widgets.Label()]+op_list)
op_label_box = widgets.VBox([widgets.Label(value='Operations:', layout=thin_layout)]+op_labels)
for box in op_list: box.disabled = True
feat_box = widgets.VBox([widgets.Label()]+feat_list)
feat_label_box = widgets.VBox([widgets.Label(value='Features:', layout=thin_layout)]+feat_labels)
for box in feat_list: box.disabled = True
rung_selection = widgets.Dropdown(options=['PRM2020', 1,2,3], value=2,layout=thin_layout)
rung_selection.value = 'PRM2020'
feat_per_iter_selection = widgets.BoundedIntText(value = 50, min=10, max=200, step=1, layout=thin_layout)
dimension_selection = widgets.BoundedIntText(value = 2, min=1, max=4, step=1, layout = thin_layout)
settings_box = widgets.VBox([
widgets.Label(value='Settings:', layout=wide_layout),
widgets.Label(value='SISSO rung:', layout=wide_layout),
widgets.Label(value='To unfreeze the feature selection,' , layout=wide_layout),
widgets.Label(value='please select any rung other than PRM2020.', layout=widgets.Layout(width = '300px', bottom='10px') ),
widgets.Label(value='Number of selected features per SIS iteration:', layout=wide_layout),
widgets.Label(value='Maximum number of dimensions:', layout=wide_layout),
default_button = widgets.Button(description = 'Default selection', layout=mid_layout)
run_button = widgets.Button(description = 'Run', layout=mid_layout)
plot_button = widgets.Button(description = 'Plot interactive map', disabled=True, layout=mid_layout)
button_box = widgets.VBox([default_button, run_button, plot_button])
rung_selection.observe(handle_rung_selection, names='value')
out1 = widgets.Output()
out2 = widgets.Output()
gui_box = widgets.HBox([op_box, op_label_box, feat_box, feat_label_box, settings_box, button_box])
out_box = widgets.VBox([gui_box, out1, out2])
%%%% Output: display_data
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment