<jupyter_start><jupyter_text>Spacegroup determination: comparing spglib vs neural network based on the diffraction intensity in spherical harmonics (DISH) descriptorAuthor: <NAME> (<EMAIL>; <EMAIL>)<jupyter_code>from ai4materials.dataprocessing.preprocessing import load_dataset_from_file
from ai4materials.wrappers import load_descriptor
from ai4materials.utils.utils_config import set_configs
from ai4materials.utils.utils_config import setup_logger
from ai4materials.dataprocessing.preprocessing import load_dataset_from_file
from ai4materials.dataprocessing.preprocessing import prepare_dataset
from ai4materials.utils.utils_crystals import random_displace_atoms, create_vacancies
from ai4materials.utils.utils_crystals import get_spacegroup, get_spacegroup_analyzer
from ai4materials.utils.utils_data_retrieval import read_ase_db
from ai4materials.utils.utils_data_retrieval import read_ase_db
from ai4materials.utils.utils_crystals import create_supercell
import random 
from ase.spacegroup import get_spacegroup as ase_get_spacegroup
from asap3.analysis.localstructure import RestrictedCNA
from asap3 import FullNeighborList
from asap3.analysis import CoordinationNumbers, FullCNA, PTM
from collections import Counter
import itertools
import matplotlib.pyplot as plt
import numpy as np
import os
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix
import ase.io
%matplotlib inline  

main_folder = '/home/ziletti/Documents/calc_nomadml/rot_inv_3d/'
dataset_folder = os.path.abspath(os.path.normpath(os.path.join(main_folder, 'datasets')))
desc_folder = os.path.abspath(os.path.normpath(os.path.join(main_folder, 'desc_folder')))

ase_db_dataset_dir = '/home/ziletti/Documents/calc_nomadml/rot_inv_3d/db_ase'

configs = set_configs(main_folder=main_folder)
logger = setup_logger(configs, level='INFO', display_configs=False)
configs['io']['dataset_folder'] = dataset_folder
configs['io']['desc_folder'] = desc_folder
def plot_confusion_matrix(cm, classes,
                          normalize=False,
                          title='Confusion matrix',
                          cmap=plt.cm.Blues):
    """
    This function prints and plots the confusion matrix.
    Normalization can be applied by setting `normalize=True`.
    """
    if normalize:
        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
        print("Normalized confusion matrix")
    else:
        print('Confusion matrix, without normalization')

    print(cm)

    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    plt.title(title)
    plt.colorbar()
    tick_marks = np.arange(len(classes))
    plt.xticks(tick_marks, classes, rotation=45)
    plt.yticks(tick_marks, classes)

    fmt = '.2f' if normalize else 'd'
    thresh = cm.max() / 2.
    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
        plt.text(j, i, format(cm[i, j], fmt),
                 horizontalalignment="center",
                 color="white" if cm[i, j] > thresh else "black")

    plt.ylabel('True label')
    plt.xlabel('Predicted label')
    plt.tight_layout()<jupyter_output><empty_output><jupyter_text>1. Loading prototype database and check the classification for pristine structures<jupyter_code>db_files_prototypes_basedir = '/home/ziletti/Documents/calc_nomadml/rot_inv_3d/db_ase/'

proto_names = ['A_cP1_221_a']

prototypes_basedir = '/home/ziletti/Documents/calc_nomadml/rot_inv_3d/prototypes_aflow_new'
ase_db_files = [os.path.join(db_files_prototypes_basedir, proto_name) + '.db' for proto_name in proto_names]

db_protos = zip(proto_names, ase_db_files)

y_pred = []
for idx_db, db_proto in enumerate(db_protos):
    ase_atoms_list = read_ase_db(db_path=ase_db_files[idx_db])<jupyter_output><empty_output><jupyter_text>Here we load the pristine dataset; one can load other dataset simply uncommenting the corresponding line below.E.g. to load the dataset with 1% displacements:1. comment the line *ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_pristine' + '.db')*2. uncomment the line    * ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_displacement-1%' + '.db')*<jupyter_code>ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_pristine' + '.db')

# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_displacement-0.1%' + '.db')
# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_displacement-0.2%' + '.db')
# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_displacement-0.6%' + '.db')

# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_displacement-1%' + '.db')
# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_displacement-2%' + '.db')
# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_displacement-4%' + '.db')
# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_displacement-5%' + '.db')
# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_displacement-8%' + '.db')

# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_vacancies-1%' + '.db')
# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_vacancies-2%' + '.db')
# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_vacancies-5%' + '.db')
# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_vacancies-10%' + '.db')
# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_vacancies-20%' + '.db')
# ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_vacancies-50%' + '.db')

# read the ASE atom database containing the structures to classify
ase_atoms = read_ase_db(db_path=ase_db)<jupyter_output><empty_output><jupyter_text>2. Polyhedral template matching benchmarkingThis modifier classifies the local structural environment of particles using the Polyhedral Template Matching (PTM) method [P.M. Larsen et al., Modelling Simul. Mater. Sci. Eng. 24, 055007 (2016)] and allows identifying the local crystalline structure of simple condensed phases (FCC, BCC, etc.) formed by atoms.The PTM approach promises greater reliability than e.g. Common Neighbor Analysis (CNA) in the presence of strong thermal fluctuations and strain. Furthermore, it provides a direct calculation of the local (per-atom) orientation, elastic deformation, strain and alloy type.The modifier assigns the Structure Type property to each particle with the assigned structural type encoded as an integer value:- 0 = Other, unknown coordination structure- 1 = FCC, face-centered cubic- 2 = HCP, hexagonal close-packed- 3 = BCC, body-centered cubic- 4 = ICO, icosahedral coordination- 5 = SC, simple cubicWe use the ASAP implementation of PTM: (https://wiki.fysik.dtu.dk/asap/Running%20Polyhedral%20Template%20Matching)<jupyter_code>#ptm.keys()
#classes = ['None', 'FCC', 'HCP', 'BCC', 'Ic', 'SC']
classes = ['None', '225', '194', '229', 'Ic', '221']
y_pred_ptm = []
y_true = []

for idx, atoms in enumerate(ase_atoms):
    if idx%1000 == 0: 
        print("{}/{}".format(idx, len(ase_atoms)))
    
    if str(atoms.info['target']) == '227':
        # skip the diamond structures because PTM cannot classify them
        pass
    else:
        # make a supercell
        atoms=atoms*(2, 2, 2)
        ptm = PTM(atoms=atoms)
        atom_class_ptm = [classes[item] for item in ptm['structure']]
        atom_class_true = [str(atoms.info['target'])]*len(atoms)
        y_true.extend(atom_class_true)
        y_pred_ptm.extend(atom_class_ptm)

print("Total number of atoms in the structures considered: {}".format(len(y_true)))
print(Counter(y_true), Counter(y_pred_ptm))
print('Accuracy: {}'.format(accuracy_score(y_true, y_pred_ptm)))
cnf_matrix = confusion_matrix(y_true, y_pred_ptm)
np.set_printoptions(precision=2)
#print(cnf_matrix)
plot_confusion_matrix(cnf_matrix, classes=['194', '221', '225', '229', 'Ic', 'None'],
                      normalize=False, title='Confusion matrix, without normalization')
# del y_true
# del y_pred_ptm<jupyter_output><empty_output><jupyter_text>3. Other benchmarkingsOther methods were also benchmarks: 3.1 Conventional CNA (with fixed cutoff)Here, a threshold distance criterion is used to determine whether a pair of atoms is bonded or not. The cutoff distance must be chosen according to the crystal structure at hand. For face-centered cubic (FCC) and hexagonal close-packed (HCP) structures the cutoff radius must lie midway between the first and the second shell of neighbors. For body-centered cubic (BCC) structures the cutoff radius should be positioned between the second and the third neighbor shell. OVITO provides a list of optimal cutoff distances for FCC and BCC crystal structures formed by common pure elements. These optimal radii can be found in the Presets drop-down list. 3.2 Adaptive CNA (with variable cutoff)Sometimes it may be difficult to choose the right cutoff radius for the conventional CNA, in particular in the case of multiphase systems. This is why an adaptive version of the CNA has been developed that works without a fixed cutoff. The Adaptive Common Neighbor Analysis (a-CNA) method [Stukowski, Modell. Simul. Mater. Sci. Eng. 20, 045021] determines the optimal cutoff radius automatically for each individual particle. 3.3. Bond-angle analysis (a.k.a. Ackland-Jones analysis)This modifier implements a method for identifying common crystalline and other structures based on an analysis of the distribution of angles formed by the pairs of neighbors of a central atom. The method is known as Ackland-Jones bond-angle method [Ackland and Jones, Phys. Rev. B 73, 054104]. The algorithm assigns a structural type to each particle having a local environment that matches one of the known structures (FCC, BCC, HCP, icosahedral).The text is taken from the Ovito documentation:     https://www.ovito.org/manual/particles.modifiers.common_neighbor_analysis.html, https://ovito.org/manual_testing/particles.modifiers.bond_angle_analysis.html    We use the Ovito implementation of these methods. 3.4 Symmetry-based spacegroup determination (with Spglib)Spglib is a library for finding and handling crystal symmetries written in C. There are the python, fortran, and ruby interfaces. Version history is found at ChangeLog.The detailed algorithm of spglib is describted the following text:"𝚂𝚙𝚐𝚕𝚒𝚋: a software library for crystal symmetry search", <NAME> and <NAME>, https://arxiv.org/abs/1808.01590 (written at version 1.10.4)Website: https://atztogo.github.io/spglib/  4. Classification results spglib<jupyter_code># read the ASE atom database containing the structures to classify
ase_db = os.path.join(ase_db_dataset_dir, 'hcp-sc-fcc-diam-bcc_pristine' + '.db')
ase_atoms = read_ase_db(db_path=ase_db)

# uncomment the line below to select a subset of structures
# ase_atoms = ase_atoms[:100]
# y_pred_spglib = []
# y_true = []

# for idx, atoms in enumerate(ase_atoms):
#     if idx%1000 == 0: 
#         print("{}/{}".format(idx, len(ase_atoms)))
        
#     atom_class_true = [str(atoms.info['target'])] 
#     y_true.extend(atom_class_true)
#     y_pred_spglib.append(str(ase_get_spacegroup(atoms, symprec=1e-1).no))
    
# print('Accuracy: {}'.format(accuracy_score(y_true, y_pred_spglib)))
# print(Counter(y_true), Counter(y_pred_spglib))

# # print confusion matrix
# cnf_matrix = confusion_matrix(y_true, y_pred_spglib)
# np.set_printoptions(precision=2)
# print(cnf_matrix)<jupyter_output>0/100<jupyter_text>4.1 Spglib benchmarking on defective structures<jupyter_code># true labels
y_true = []

# classification from spglib with different thresholds
y_pred_symprec1 = []
y_pred_symprec2 = []
y_pred_symprec3 = []

ase_atoms_supercell = []

for idx, atoms in enumerate(ase_atoms):
    if idx%1000 == 0: 
        print("{}/{}".format(idx, len(ase_atoms)))
    
    # retrieve the correct class
    atom_class_true = [str(atoms.info['target'])] 
    y_true.extend(atom_class_true)

    # from the pristine supercell, create displaced structures
#     supercell =  random_displace_atoms(atoms, displacement_scaled=0.02, 
#                                       create_replicas_by='user-defined', 
#                                       noise_distribution='uniform_scaled', 
#                                       target_replicas=[1,1,1],
#                                       optimal_supercell=False)
    
    supercell =  create_vacancies(atoms, target_vacancy_ratio=0.02, 
                                      create_replicas_by='user-defined',
                                      target_replicas=[1,1,1],
                                      optimal_supercell=False)
    
    alpha = random.random() * 360.0
    atoms.rotate(alpha, 'x', rotate_cell=True, center='COU')

    beta = random.random() * 360.0
    atoms.rotate(beta, 'y', rotate_cell=True, center='COU')

    gamma = random.random() * 360.0
    atoms.rotate(gamma, 'z', rotate_cell=True, center='COU')
    
    ase_atoms_supercell.append(supercell)
    
    get_spacegroup(supercell, symprec=[1e-1], angle_tolerance=5.0)
    y_pred_symprec1.append(str(supercell.info['spacegroup_nb_0.1']))
#     y_pred_symprec3.append(str(supercell.info['spacegroup_nb_0.001']))


    # uncomment the lines below to calculate the Spglib classification for different thresholds
    # get_spacegroup(supercell, symprec=[1e-1, 1e-2,1e-3], angle_tolerance=5.0)
    # y_pred_symprec1.append(str(supercell.info['spacegroup_nb_0.1']))
    # y_pred_symprec2.append(str(supercell.info['spacegroup_nb_0.01']))
    # y_pred_symprec3.append(str(supercell.info['spacegroup_nb_0.001']))

print("Accuracy (spglib symprec=1e-1, angle_tolerance=5.0): {}".format(accuracy_score(y_true, y_pred_symprec1)))
# print("Accuracy (spglib symprec=1e-2, angle_tolerance=5.0): {}".format(accuracy_score(y_true, y_pred_symprec2)))
# print("Accuracy (spglib symprec=1e-3, angle_tolerance=5.0): {}".format(accuracy_score(y_true, y_pred_symprec3)))
Counter(y_true), Counter(y_pred_symprec1), Counter(y_pred_symprec2), Counter(y_pred_symprec3)<jupyter_output><empty_output><jupyter_text>5. Visualization of crystal structures<jupyter_code>from ase.visualize import view
idx=1
# change idx if you want to visualize another structure 
# idx=0 visualizes the first structure in the list
# Note: 0<=idx<len(ase_atoms_supercell)
view(ase_atoms_supercell[idx]*(1, 1, 1), viewer='x3d')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11149947570779448, 'fraction_numerical': 0.018245368752184552, 'mean_word_length': 4.612397018438603, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 34, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '615327', 'license': 'Apache-2.0', 'n_tokens_mistral': 5081, 'n_tokens_neox': 4792, 'n_words': 1141, 'path': 'examples/single_crystal_classification.ipynb', 'repo_name': 'hpleva/ai4materials', 'size': '14359'}