Lieb Lattice Research

# notebook utils from tqdm.notebook import tqdm from IPython.display import Image # web utils import requests # data utils import pickle import itertools import numpy as np # import pymatgen import pymatgen from pymatgen.core import Lattice, Structure, Molecule from pymatgen.ext.matproj import MPRester

lieb_reference = Structure( Lattice.from_parameters(a=5, b=5, c=5, alpha=90, beta=90, gamma=90), ["Sn", "Fe", "Sn"], [[1 / 2, 0, 0], [0, 0, 0], [0, 1 / 2, 0]], ) lieb_reference.to(filename="lieb.cif")

def is_lieb_lattice(structure): distinct_species = set(structure.species) if len(structure.species) != 3 or len(distinct_species) not in (2, 3): return False # 3 total species, 2-3 distinct species if structure.lattice.angles != (90.0, 90.0, 90.0): return False # square lattice z_coords = structure.frac_coords.transpose()[-1] if not np.all(z_coords == z_coords[0]): return False # all species are at the same z-coord if len(distinct_species) == 2: # if there are 2 distinct species, make sure the edge species are the same (as opposed to the corner species and one edge species) edge_species = [structure.sites[i].species for i in np.nonzero(structure.frac_coords)[0]] if edge_species[0] != edge_species[1]: return False return True

def has_lieb_lattice(structure, max_num_lieb_species=20): if len(structure.species) < 3: return False # at least 3 total atoms (needed to form the Lieb lattice) if structure.lattice.angles.count(90.0) < 1: return False # at least one lattice angle must be 90 (the lieb lattice axis angle) for species in set(structure.species): species_idx = structure.indices_from_symbol(species.value) if len(species_idx) not in range(3, max_num_lieb_species+ 1): continue # at least three of the same element, bounded above to reduce query time if np.zeros(3) not in np.array([ structure.frac_coords[i] for i in species_idx ]): continue # at least one atom of the element is @ origin for atom_triple in itertools.combinations(species_idx, 3): # for each triple of the same element (candidate Lieb lattice atoms) nonzero_coords = np.nonzero(np.array([ structure.frac_coords[i] for i in atom_triple ])) if len(nonzero_coords[0]) != 2 or len(set(nonzero_coords[0])) != 2 or len(set(nonzero_coords[1])) != 2: return False # one @ origin, two along a lattice vector for axis, plane in enumerate(np.array([ structure.frac_coords[i] for i in atom_triple ]).transpose()): if np.all(plane[0] == plane[:]): lieb_axis, lieb_plane = axis, plane[0] if structure.lattice.angles[lieb_axis] != 90.0: return False # lieb lattice axis angle must be 90 (square lattice) if all([structure.frac_coords[i][lieb_axis] != lieb_plane for i in set(species_idx) ^ set(atom_triple)]): return True # no other species of *same element* in lieb lattice plane return False

has_lieb_lattice(lieb_reference)

def unit_vector(vector): return vector / np.linalg.norm(vector) def angle_between(v1, v2): v1_u = unit_vector(v1) v2_u = unit_vector(v2) return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0)) * 180/np.pi def get_nearest_neighbors(structure): return np.array([len(structure.get_neighbors(site, min(np.delete(structure.distance_matrix[site_idx], site_idx)) + 0.0000001 )) for site_idx, site in enumerate(structure.sites)])

def is_lieb_lattice(structure, verbose=False): if not (len(structure.species) == 3 and len(set(structure.species)) in (1, 2)): if verbose: print("# Species Constraint: # of species is %d and number of distinct species is %d." % (len(structure.species), len(set(structure.species)))) return False # exactly 3 atoms, 1 or 2 distinct atoms nearest_neighbors = get_nearest_neighbors(structure) if not (sorted(nearest_neighbors) == [2, 2, 4]): if verbose: print("Nearest Neighbors Constraint: # of nearest neighbors for each species is %s." % str(nearest_neighbors)) return False # 4 nearest neighbors for corner, 2 nearest neighbor for edge if len(set(structure.species)) == 2: edge_species = [ site.species for site in np.array(structure.sites)[np.where(nearest_neighbors == 2)[0]] ] if not (edge_species[0] == edge_species[1]): if verbose: print("Edge Species Constraint: edge species are %s and %s." % (str(edge_species[0]), str(edge_species[1]))) return False # edge species are the same element lattice_vectors = structure.cart_coords[np.where(nearest_neighbors == 2)[0]] - structure.cart_coords[np.where(nearest_neighbors == 4)[0]] if not (angle_between(lattice_vectors[0], lattice_vectors[1]) == 90.0): if verbose: print("Right Angle Constraint: angle between lattice vectors is %f." % (angle_between(lattice_vectors[0], lattice_vectors[1]) * 180/np.pi) ) return False # 90 degree angle between lattice vectors return True

def lieb_plane_is_empty(structure, lieb_species, lieb_species_idx): for stray_lieb_species in set(np.concatenate([ np.array(structure.indices_from_symbol(species.value)) for species in lieb_species ])) ^ lieb_species_idx: species_coords = np.array(structure.frac_coords)[list(lieb_species_idx.union({stray_lieb_species}))] if np.linalg.det(species_coords[1:] - species_coords[0]) == 0: return False # stray_lieb_species is coplanar w/ lieb_species return True

def has_lieb_lattice(structure, max_num_total_species=None, total_species=None, max_num_lieb_species=None, include_two_element_lattice=True): if len(structure.species) < 3 or (len(structure.species) > max_num_total_species) if max_num_total_species else False or len(structure.species) != total_species if total_species else False: return False # at least 3 atoms for species in set(structure.species): species_idx = structure.indices_from_symbol(species.value) if len(species_idx) < 2 if include_two_element_lattice else 3 or (len(species_idx) >= max_num_lieb_species) if max_num_lieb_species else False: continue # at least 2 or 3 of the same atom if include_two_element_lattice: for edge_species_pair in itertools.combinations(species_idx, 2): # every triple of (2 same atom, 1 other atom) for corner_species in set(range(len(structure.species))) ^ set(species_idx): lieb_species_idx = set(edge_species_pair).union({corner_species}) lieb_structure = Structure.from_sites(list(np.array(structure.sites)[list(lieb_species_idx)])) if is_lieb_lattice(lieb_structure) and lieb_plane_is_empty(structure, {species, structure.species[corner_species]}, lieb_species_idx): return lieb_structure if len(species_idx) >= 3: for lieb_species_idx in itertools.combinations(species_idx, 3): # every triple of same atom lieb_structure = Structure.from_sites(list(np.array(structure.sites)[list(lieb_species_idx)])) if is_lieb_lattice(lieb_structure) and lieb_plane_is_empty(structure, {species}, set(lieb_species_idx)): return lieb_structure return False

material_id = "mp-1019593"

with MPRester(API_KEY) as m: result = m.query(criteria={"material_id": material_id}, properties=["pretty_formula", "structure"])

structure = result[0]["structure"] structure

is_lieb_lattice(structure, verbose=True)

for species in set(structure.species): print(len(structure.indices_from_symbol(species.value)))

lieb_structure = has_lieb_lattice(structure) lieb_structure

is_lieb_lattice(lieb_structure)

with open('mapi_localdb.pkl', 'rb') as f: local_db = pickle.load(f) len(local_db)

lieb_lattices = { compound['material_id']: compound for compound in tqdm(local_db) if is_lieb_lattice(compound["structure"])}

len(lieb_lattices)

for id, lattice in lieb_lattices.items(): print(id, lattice['pretty_formula'])

with open('is_lieb_lattice_db.pkl', 'wb') as f: pickle.dump(lieb_lattices, f)

lieb_lattices = {} for material in tqdm(local_db): # print(material["material_id"], material["pretty_formula"]) lieb_structure = has_lieb_lattice(material["structure"], total_species=10) if lieb_structure: lieb_lattices[material["material_id"]] = { "material_properties": material, "lieb_substructure": lieb_structure }

len(lieb_lattices)

for id, lattice in lieb_lattices.items(): print(id, lattice['material_properties']['pretty_formula'], lattice['lieb_substructure'].species) # lattice['lieb_structure'].to(filename=id+"_lieb.cif")

# get bandstructure material_id = "mp-637576" response = requests.get("https://www.materialsproject.org/rest/v2/materials/%s/vasp/bandstructure" % material_id, params={'API_KEY': API_KEY})

for material_id in lieb_lattices.keys(): print("%s (%s)" % (lieb_lattices[material_id]['material_properties']['pretty_formula'], material_id)) display(Image(url='https://materialsproject.org/electronic_structure/bandstructure/plot/%s' % material_id))