From b56ded5fa2033b4a7cc65028bd631efde01c5964 Mon Sep 17 00:00:00 2001 From: Carlo Antonio Pignedoli Date: Fri, 13 Dec 2024 08:53:13 +0000 Subject: [PATCH 1/3] using dim finder --- surfaces_tools/widgets/analyze_structure.py | 614 ++++--- surfaces_tools/widgets/lowdimfinder.py | 1620 +++++++++++++++++++ surfaces_tools/widgets/nn.py | 468 ++++++ 3 files changed, 2375 insertions(+), 327 deletions(-) create mode 100644 surfaces_tools/widgets/lowdimfinder.py create mode 100644 surfaces_tools/widgets/nn.py diff --git a/surfaces_tools/widgets/analyze_structure.py b/surfaces_tools/widgets/analyze_structure.py index c060fdd..63f4e33 100644 --- a/surfaces_tools/widgets/analyze_structure.py +++ b/surfaces_tools/widgets/analyze_structure.py @@ -1,58 +1,216 @@ import copy import itertools - import aiidalab_widgets_base as awb -import ase + import numpy as np +from ase import Atoms, neighborlist +from scipy import sparse +from scipy.signal import find_peaks +from scipy.spatial import ConvexHull +from .lowdimfinder import LowDimFinder import traitlets as tr -from ase import neighborlist -from scipy import signal, sparse, spatial - - -def to_ranges(iterable): - iterable = sorted(set(iterable)) - for _key, group in itertools.groupby(enumerate(iterable), lambda t: t[1] - t[0]): - lgroup = list(group) - yield lgroup[0][1], lgroup[-1][1] - - -def mol_ids_range(ismol): - # Shifts the list by +1. - range_string = "" - shifted_list = [i + 1 for i in ismol] - ranges = list(to_ranges(shifted_list)) - for i in range(len(ranges)): - if ranges[i][1] > ranges[i][0]: - range_string += f"{ranges[i][0]}..{ranges[i][1]}" - else: - range_string += f"{ranges[i][0]} " - return range_string - - -def conne_matrix(atoms): - cutoff = neighborlist.natural_cutoffs(atoms) - neighbor_list = neighborlist.NeighborList( - cutoff, self_interaction=False, bothways=False +def gaussian(x, sig): + return ( + 1.0 + / (sig * np.sqrt(2.0 * np.pi)) + * np.exp(-np.power(x, 2.0) / (2 * np.power(sig, 2.0))) ) - neighbor_list.update(atoms) - - return neighbor_list.get_connectivity_matrix() +def boxfilter(x, thr): + return np.asarray([1 if i < thr else 0 for i in x]) + # Piero Gasparotto +def get_types( frame): # classify the atmos in: + # 0=molecule + # 1=slab atoms + # 2=adatoms + # 3=hydrogens on the surf + # 5=unknown + # 6=metalating atoms + # frame=ase frame + # thr=threashold in the histogram for being considered a surface layer + nat = frame.get_global_number_of_atoms() + + # all atom types set to 5 (unknown) + atype = np.zeros(nat, dtype=np.int16) + 5 + area = frame.cell[0][0] * frame.cell[1][1] + minz = np.min(frame.positions[:, 2]) + maxz = np.max(frame.positions[:, 2]) + + if maxz - minz < 1.0: + maxz += (1.0 - (maxz - minz)) / 2 + minz -= (1.0 - (maxz - minz)) / 2 + + # Which values should we use below? + sigma = 0.2 # thr + peak_rel_height = 0.5 + layer_tol = 1.0 * sigma + + # Quick estimate number atoms in a layer: + nbins = int(np.ceil((maxz - minz) / 0.15)) + hist, _ = np.histogram(frame.positions[:, 2], density=False, bins=nbins) + max_atoms_in_a_layer = max(hist) + + lbls = frame.get_chemical_symbols() + n_intervals = int(np.ceil((maxz - minz + 3 * sigma) / (0.1 * sigma))) + z_values = np.linspace(minz - 3 * sigma, maxz + 3 * sigma, n_intervals) # 1000 + atoms_z_pos = frame.positions[:, 2] + + # OPTION 1: generate 2d array to apply the gaussian on + z_v_exp, at_z_exp = np.meshgrid(z_values, atoms_z_pos) + arr_2d = z_v_exp - at_z_exp + atomic_density = np.sum(gaussian(arr_2d, sigma), axis=0) + + # OPTION 2: loop through atoms + # atomic_density = np.zeros(z_values.shape) + # for ia in range(len(atoms)): + # atomic_density += gaussian(z_values - atoms.positions[ia,2], sigma) + + peaks = find_peaks( + atomic_density, + height=None, + threshold=None, + distance=None, + prominence=None, + width=None, + wlen=None, + rel_height=peak_rel_height, + ) + layersg = z_values[peaks[0].tolist()] -def clusters(matrix): - nclusters, idlist = sparse.csgraph.connected_components(matrix) - return nclusters, [np.where(idlist == i)[0].tolist() for i in range(nclusters)] + # Check top and bottom layers should be documented better + found_top_surf = False + while not found_top_surf: + iz = layersg[-1] + two_d_atoms = [ + frame.positions[i, 0:2] + for i in range(nat) + if np.abs(frame.positions[i, 2] - iz) < layer_tol + ] + coverage = 0 + if len(two_d_atoms) > max_atoms_in_a_layer / 4: + hull = ConvexHull(two_d_atoms) + coverage = hull.volume / area + if coverage > 0.3: + found_top_surf = True + else: + layersg = layersg[0:-1] + + found_bottom_surf = False + while not found_bottom_surf: + iz = layersg[0] + two_d_atoms = [ + frame.positions[i, 0:2] + for i in range(nat) + if np.abs(frame.positions[i, 2] - iz) < layer_tol + ] + coverage = 0 + if len(two_d_atoms) > max_atoms_in_a_layer / 4: + hull = ConvexHull(two_d_atoms) + coverage = hull.volume / area + if coverage > 0.3 and len(two_d_atoms) > max_atoms_in_a_layer / 4: + found_bottom_surf = True + else: + layersg = layersg[1:] -def molecules(ismol, atoms): - if ismol: - nmols, ids = clusters(conne_matrix(atoms[ismol])) - return [[ismol[i] for i in ids[j]] for j in range(nmols)] - return [] + bottom_z = layersg[0] + top_z = layersg[-1] + # check if there is a bottom layer of H + found_layer_of_h = True + for i in range(nat): + iz = frame.positions[i, 2] + if (layer_tol + iz) > bottom_z and iz < (bottom_z + layer_tol): + if lbls[i] == "H": + atype[i] = 3 + else: + found_layer_of_h = False + break + if found_layer_of_h: + layersg = layersg[1:] + # bottom_z=layersg[0] + + layers_dist = [] + iprev = layersg[0] + for inext in layersg[1:]: + layers_dist.append(abs(iprev - inext)) + iprev = inext + + for i in range(nat): + iz = frame.positions[i, 2] + if (layer_tol + iz) > bottom_z and iz < (top_z + layer_tol): + if not (atype[i] == 3 and found_layer_of_h): + atype[i] = 1 + else: + if np.min([np.abs(iz - top_z), np.abs(iz - bottom_z)]) < np.max( + layers_dist + ): + if not (atype[i] == 3 and found_layer_of_h): + atype[i] = 2 + + # assign the other types + #metalatingtypes = ("Au", "Ag", "Cu", "Ni", "Co", "Zn", "Mg") + #moltypes = ("H", "N", "B", "O", "C", "F", "S", "Br", "I", "Cl") + #possible_mol_atoms = [ + # i for i in range(nat) if atype[i] == 2 and lbls[i] in moltypes + #] + #possible_mol_atoms += [i for i in range(nat) if atype[i] == 5] + # identify separate molecules + # all_molecules=self.molecules(mol_atoms,atoms) + #all_molecules = [] + + return atype, layersg + + +def transform_vector(pos1, pos2, v1): + """ + Transform vector v1 using the transformation that maps pos2 to pos1. + + Args: + pos1 (np.ndarray): Target positions, shape (n, 3). + pos2 (np.ndarray): Source positions, shape (n, 3). + v1 (np.ndarray): Vector to transform, shape (3,). + + Returns: + np.ndarray: Transformed vector, shape (3,). + """ + # Ensure inputs are numpy arrays + pos1 = np.array(pos1) + pos2 = np.array(pos2) + v1 = np.array(v1) + + # Compute centroids of both point sets + centroid1 = np.mean(pos1, axis=0) + centroid2 = np.mean(pos2, axis=0) + + # Center the points + pos1_centered = pos1 - centroid1 + pos2_centered = pos2 - centroid2 + + # Compute the covariance matrix + H = np.dot(pos2_centered.T, pos1_centered) + + # Singular Value Decomposition (SVD) for rotation + U, S, Vt = np.linalg.svd(H) + R = np.dot(Vt.T, U.T) + + # Ensure a proper rotation matrix (det(R) = 1) + if np.linalg.det(R) < 0: + Vt[-1, :] *= -1 + R = np.dot(Vt.T, U.T) + + # Translation vector + t = centroid1 - np.dot(R, centroid2) + + # Transform the vector + transformed_v1 = -np.dot(R, v1) #+ t + transformed_v1 = transformed_v1 / np.linalg.norm(transformed_v1) + + + return transformed_v1 class StructureAnalyzer(tr.HasTraits): - structure = tr.Instance(ase.Atoms, allow_none=True) + structure = tr.Instance(Atoms, allow_none=True) details = tr.Dict() def __init__(self, only_sys_type=False, who_called_me="Boh"): @@ -60,217 +218,21 @@ def __init__(self, only_sys_type=False, who_called_me="Boh"): self.who_called_me = who_called_me super().__init__() - def gaussian(self, x, sig): - return ( - 1.0 - / (sig * np.sqrt(2.0 * np.pi)) - * np.exp(-np.power(x, 2.0) / (2 * np.power(sig, 2.0))) - ) - - def boxfilter(self, x, thr): - return np.asarray([1 if i < thr else 0 for i in x]) - - def get_types(self, frame, thr): # Piero Gasparotto - # classify the atmos in: - # 0=molecule - # 1=slab atoms - # 2=adatoms - # 3=hydrogens on the surf - # 5=unknown - # 6=metalating atoms - # frame=ase frame - # thr=threashold in the histogram for being considered a surface layer - nat = frame.get_number_of_atoms() - - # all atom types set to 5 (unknown) - atype = np.zeros(nat, dtype=np.int16) + 5 - area = frame.cell[0][0] * frame.cell[1][1] - minz = np.min(frame.positions[:, 2]) - maxz = np.max(frame.positions[:, 2]) - - if maxz - minz < 1.0: - maxz += (1.0 - (maxz - minz)) / 2 - minz -= (1.0 - (maxz - minz)) / 2 - - # Which values should we use? - sigma = 0.2 - peak_rel_height = 0.5 - layer_tol = 1.0 * sigma - - # Quick estimate number atoms in a layer. - nbins = int(np.ceil((maxz - minz) / 0.15)) - hist, bin_edges = np.histogram(frame.positions[:, 2], density=False, bins=nbins) - max_atoms_in_a_layer = max(hist) - - lbls = frame.get_chemical_symbols() - n_intervals = int(np.ceil((maxz - minz + 3 * sigma) / (0.1 * sigma))) - z_values = np.linspace(minz - 3 * sigma, maxz + 3 * sigma, n_intervals) # 1000 - atoms_z_pos = frame.positions[:, 2] - - # Generate 2d array to apply the gaussian on. - z_v_exp, at_z_exp = np.meshgrid(z_values, atoms_z_pos) - arr_2d = z_v_exp - at_z_exp - atomic_density = np.sum(self.gaussian(arr_2d, sigma), axis=0) - - peaks = signal.find_peaks( - atomic_density, - height=None, - threshold=None, - distance=None, - prominence=None, - width=None, - wlen=None, - rel_height=peak_rel_height, - ) - layersg = z_values[peaks[0].tolist()] - len(layersg) - layersg[-1] - - # Check top and bottom layers should be documented better. - found_top_surf = False - while not found_top_surf: - iz = layersg[-1] - two_d_atoms = [ - frame.positions[i, 0:2] - for i in range(nat) - if np.abs(frame.positions[i, 2] - iz) < layer_tol - ] - coverage = 0 - if len(two_d_atoms) > max_atoms_in_a_layer / 4: - hull = spatial.ConvexHull(two_d_atoms) - coverage = hull.volume / area - if coverage > 0.3: - found_top_surf = True - else: - layersg = layersg[0:-1] - - found_bottom_surf = False - while not found_bottom_surf: - iz = layersg[0] - two_d_atoms = [ - frame.positions[i, 0:2] - for i in range(nat) - if np.abs(frame.positions[i, 2] - iz) < layer_tol - ] - coverage = 0 - if len(two_d_atoms) > max_atoms_in_a_layer / 4: - hull = spatial.ConvexHull(two_d_atoms) - coverage = hull.volume / area - if coverage > 0.3 and len(two_d_atoms) > max_atoms_in_a_layer / 4: - found_bottom_surf = True - else: - layersg = layersg[1:] - - bottom_z = layersg[0] - top_z = layersg[-1] - - # Check if there is a bottom layer of H. - found_layer_of_h = True - for i in range(nat): - iz = frame.positions[i, 2] - if iz > bottom_z - layer_tol and iz < bottom_z + layer_tol: - if lbls[i] == "H": - atype[i] = 3 - else: - found_layer_of_h = False - break - if found_layer_of_h: - layersg = layersg[1:] - # bottom_z=layersg[0] - - layers_dist = [] - iprev = layersg[0] - for inext in layersg[1:]: - layers_dist.append(abs(iprev - inext)) - iprev = inext - - for i in range(nat): - iz = frame.positions[i, 2] - if iz > bottom_z - layer_tol and iz < top_z + layer_tol: - if not (atype[i] == 3 and found_layer_of_h): - atype[i] = 1 - else: - if np.min([np.abs(iz - top_z), np.abs(iz - bottom_z)]) < np.max( - layers_dist - ): - if not (atype[i] == 3 and found_layer_of_h): - atype[i] = 2 - - # assign the other types - metalatingtypes = ("Au", "Ag", "Cu", "Ni", "Co", "Zn", "Mg") - moltypes = ("H", "N", "B", "O", "C", "F", "S", "Br", "I", "Cl") - possible_mol_atoms = [ - i for i in range(nat) if atype[i] == 2 and lbls[i] in moltypes - ] - possible_mol_atoms += [i for i in range(nat) if atype[i] == 5] - - # Identify separate molecules. - all_molecules = [] - if len(possible_mol_atoms) > 0: - fragments = molecules(possible_mol_atoms, frame) - all_molecules = copy.deepcopy(fragments) - - # Remove isolated atoms. - for frag in fragments: - if len(frag) == 1: - all_molecules.remove(frag) - else: - for atom in frag: - if lbls[atom] in metalatingtypes: - atype[atom] = 6 - else: - atype[atom] = 0 - - return atype, layersg, all_molecules - - def all_connected_to(self, id_atom, atoms, exclude): - cov_radii = [ase.data.covalent_radii[a.number] for a in atoms] - - atoms.set_pbc([False, False, False]) - nl_no_pbc = neighborlist.NeighborList( - cov_radii, bothways=True, self_interaction=False - ) - nl_no_pbc.update(atoms) - atoms.set_pbc([True, True, True]) - - tofollow = [] - followed = [] - isconnected = [] - tofollow.append(id_atom) - isconnected.append(id_atom) - while len(tofollow) > 0: - indices, offsets = nl_no_pbc.get_neighbors(tofollow[0]) - indices = list(indices) - followed.append(tofollow[0]) - for i in indices: - if (i not in isconnected) and (atoms[i].symbol not in exclude): - tofollow.append(i) - isconnected.append(i) - for i in followed: - if i in tofollow: # do not remove this check - tofollow.remove(i) - - return isconnected - - def string_range_to_list(self, a): - singles = [int(s) - 1 for s in a.split() if s.isdigit()] - ranges = [r for r in a.split() if ".." in r] - for r in ranges: - t = r.split("..") - to_add = [i - 1 for i in range(int(t[0]), int(t[1]) + 1)] - singles += to_add - return sorted(singles) - @tr.observe("structure") def _observe_structure(self, _=None): with self.hold_trait_notifications(): self.details = self.analyze() - def analyze(self): + def analyze( + self, + ): if self.structure is None: return {} atoms = self.structure + + + sys_size = np.ptp(atoms.positions, axis=0) no_cell = ( atoms.cell[0][0] < 0.1 or atoms.cell[1][1] < 0.1 or atoms.cell[2][2] < 0.1 @@ -280,6 +242,35 @@ def analyze(self): atoms.cell = sys_size + 10 atoms.set_pbc([True, True, True]) + + # LowDimFinder section + low_dim_finder = LowDimFinder( + # This is called aiida_structure but it's wrong, it's an ASE structure + #aiida_structure=conventional_asecell, + aiida_structure=atoms, + vacuum_space=40.0, + radii_offset=-0.75, # [-0.75 -0.7, -0.65, -0.6, -0.55] + bond_margin=0.0, + max_supercell=3, + min_supercell=3, + rotation=True, + full_periodicity=False, + radii_source="alvarez", + orthogonal_axis_2D=True, + ) + res = low_dim_finder.get_group_data() + # End LowDimFinder section + is_a_bulk = np.amax(res['dimensionality']) == 3 + is_a_molecule = np.amax(res['dimensionality']) == 0 + is_a_wire = np.amax(res['dimensionality']) == 1 + is_a_slab = np.amax(res['dimensionality']) == 2 + max_dim_portion = res['dimensionality'].index(max(res['dimensionality'])) + # orientation + if not is_a_bulk and not is_a_molecule: + direction = transform_vector(atoms.positions[res['unit_cell_ids'][max_dim_portion]],np.array(res['positions'][max_dim_portion]),np.array([0,0,1])) + dir_short = [f"{x:.2f}" for x in direction] + # + total_charge = np.sum(atoms.get_atomic_numbers()) bottom_h = [] @@ -290,116 +281,85 @@ def analyze(self): unclassified = [] slabatoms = [] slab_layers = [] - all_molecules = None - is_a_bulk = False - is_a_molecule = False - is_a_wire = False + all_molecules = [res['unit_cell_ids'][i] for i,j in enumerate(res['dimensionality']) if j == 0] spins_up = [the_a.index for the_a in atoms if the_a.tag == 1] spins_down = [the_a.index for the_a in atoms if the_a.tag == 2] other_tags = [the_a.index for the_a in atoms if the_a.tag > 2] - # Check if there is vacuum otherwise classify as bulk and skip. - vacuum_x = sys_size[0] + 4 < atoms.cell[0][0] - vacuum_y = sys_size[1] + 4 < atoms.cell[1][1] - vacuum_z = sys_size[2] + 4 < atoms.cell[2][2] + # Do not use a set in the following line list(set(atoms.get_chemical_symbols())) - # Need ALL atoms and elements for spin guess and for cost calculation + # need ALL atoms and elements for spin guess and for cost calculation all_elements = atoms.get_chemical_symbols() - summary = "" cases = [] if len(spins_up) > 0: - summary += "spins_up: " + awb.utils.list_to_string_range(spins_up) + "\n" + summary += "spins_up: " + awb.utils.llist_to_string_range(spins_up) + "\n" if len(spins_down) > 0: - summary += ( - "spins_down: " + awb.utils.list_to_string_range(spins_down) + "\n" - ) + summary += "spins_down: " + awb.utils.llist_to_string_range(spins_down) + "\n" if len(other_tags) > 0: - summary += ( - "other_tags: " + awb.utils.list_to_string_range(other_tags) + "\n" - ) - if (not vacuum_z) and (not vacuum_x) and (not vacuum_y): - is_a_bulk = True + summary += "other_tags: " + awb.utils.llist_to_string_range(other_tags) + "\n" + + if is_a_bulk : + sys_type = "Bulk" cases = ["b"] summary += "Bulk contains: \n" slabatoms = list(range(len(atoms))) bulkatoms = slabatoms - if vacuum_x and vacuum_y and vacuum_z: - is_a_molecule = True + if is_a_molecule : + sys_type = "Molecule" - if not self.only_sys_type: - summary += "Molecule: \n" - all_molecules = molecules(list(range(len(atoms))), atoms) - com = np.average(atoms.positions, axis=0) - summary += ( - "COM: " - + str(com) - + ", min z: " - + str(np.min(atoms.positions[:, 2])) - + "\n" - ) - if vacuum_x and vacuum_y and (not vacuum_z): - is_a_wire = True - sys_type = "Wire" - cases = ["w"] - if not self.only_sys_type: - summary += "Wire along z contains: \n" - wireatoms = list(range(len(atoms))) - if vacuum_y and vacuum_z and (not vacuum_x): - is_a_wire = True - sys_type = "Wire" - cases = ["w"] - if not self.only_sys_type: - summary += "Wire along x contains: \n" - wireatoms = list(range(len(atoms))) - if vacuum_x and vacuum_z and (not vacuum_y): - is_a_wire = True + summary += "Molecule: \n" + #all_molecules = molecules(list(range(len(atoms))), atoms) + com = np.average(atoms.positions, axis=0) + summary += ( + "COM: " + + str(com) + + ", min z: " + + str(np.min(atoms.positions[:, 2])) + + "\n" + ) + if is_a_wire : + print("Wire") sys_type = "Wire" cases = ["w"] - if not self.only_sys_type: - summary += "Wire along y contains: \n" - wireatoms = list(range(len(atoms))) - # wireatoms = slabatoms - is_a_slab = not (is_a_bulk or is_a_molecule or is_a_wire) - if self.only_sys_type: - if is_a_slab: - return {"system_type": "Slab"} - else: - return {"system_type": sys_type} + summary += f"Wire along {dir_short} \n" + slabatoms = list(range(len(atoms))) + wireatoms = slabatoms + # END check - elif is_a_slab: + if is_a_slab: cases = ["s"] - tipii, layersg, all_molecules = self.get_types(atoms, 0.1) - if vacuum_x: + tipii, layersg = get_types(atoms) + if np.abs(np.dot(direction, [1, 0, 0])) > 0.9: slabtype = "YZ" - elif vacuum_y: + elif np.abs(np.dot(direction, [0, 1, 0])) > 0.9: slabtype = "XZ" else: slabtype = "XY" sys_type = "Slab" + slabtype - mol_atoms = np.where(tipii == 0)[0].tolist() - metalatings = np.where(tipii == 6)[0].tolist() - mol_atoms += metalatings + mol_atoms = [idatom for mol in all_molecules for idatom in mol] #np.where(tipii == 0)[0].tolist() + # mol_atoms=extract_mol_indexes_from_slab(atoms) + #metalatings = np.where(tipii == 6)[0].tolist() + #mol_atoms += metalatings bottom_h = np.where(tipii == 3)[0].tolist() - - # Unclassified - unclassified = np.where(tipii == 5)[0].tolist() - slabatoms = np.where(tipii == 1)[0].tolist() adatoms = np.where(tipii == 2)[0].tolist() - - # Slab layers. + unclassified = [idatom for idatom in list(range(atoms.get_global_number_of_atoms())) if idatom not in set(bottom_h+slabatoms+mol_atoms)] # np.where(tipii == 5)[0].tolist() + + + #MOVE this up to get_types slab_layers = [[] for i in range(len(layersg))] for ia in slabatoms: idx = (np.abs(layersg - atoms.positions[ia, 2])).argmin() slab_layers[idx].append(ia) + # End slab layers summary += "Slab " + slabtype + " contains: \n" summary += ( "Cell: " + " ".join([str(i) for i in atoms.cell.diagonal().tolist()]) + "\n" @@ -410,43 +370,43 @@ def analyze(self): slab_elements = set(atoms[slabatoms].get_chemical_symbols()) if len(bottom_h) > 0: - summary += "bottom H: " + mol_ids_range(bottom_h) + "\n" + summary += "bottom H: " + awb.utils.list_to_string_range(bottom_h) + "\n" if len(slabatoms) > 0: - summary += "slab atoms: " + mol_ids_range(slabatoms) + "\n" - for nlayer in range(len(slab_layers)): + summary += "slab atoms: " + awb.utils.list_to_string_range(slabatoms) + "\n" + for nlayer, the_layer in enumerate(slab_layers): summary += ( "slab layer " + str(nlayer + 1) + ": " - + mol_ids_range(slab_layers[nlayer]) + + awb.utils.list_to_string_range(the_layer) + "\n" ) if len(adatoms) > 0: cases.append("a") - summary += "adatoms: " + mol_ids_range(adatoms) + "\n" + summary += "adatoms: " + awb.utils.list_to_string_range(adatoms) + "\n" if all_molecules: cases.append("m") summary += "#" + str(len(all_molecules)) + " molecules: " - for nmols in range(len(all_molecules)): - summary += str(nmols + 1) + ") " + mol_ids_range(all_molecules[nmols]) + for nmols, the_mol in enumerate(all_molecules): + summary += str(nmols + 1) + ") " + awb.utils.list_to_string_range(the_mol) summary += " \n" if len(metalatings) > 0: + metalating_str = awb.utils.list_to_string_range(metalatings) summary += ( - "metal atoms inside molecules (already counted): " - + mol_ids_range(metalatings) - + "\n" + f"metal atoms inside molecules (already counted): {metalating_str}\n" ) if len(unclassified) > 0: cases.append("u") - summary += "unclassified: " + mol_ids_range(unclassified) + summary += "unclassified: " + awb.utils.list_to_string_range(unclassified) # Indexes from 0 if mol_ids_range is not called. + cell_str = " ".join([str(i) for i in itertools.chain(*atoms.cell.tolist())]) return { "total_charge": total_charge, "system_type": sys_type, - "cell": " ".join([str(i) for i in itertools.chain(*atoms.cell.tolist())]), + "cell": cell_str, "slab_layers": slab_layers, "bottom_H": sorted(bottom_h), "bulkatoms": sorted(bulkatoms), diff --git a/surfaces_tools/widgets/lowdimfinder.py b/surfaces_tools/widgets/lowdimfinder.py new file mode 100644 index 0000000..dd235da --- /dev/null +++ b/surfaces_tools/widgets/lowdimfinder.py @@ -0,0 +1,1620 @@ +# -*- coding: utf-8 -*- +""" +Low dimensionality atom group finder, which analyses a bulk crystal and returns +groups of atoms which are held together by weak (van der Waals) forces as +separate structures. +""" + +__copyright__ = ( + "Copyright (c), 2014-2022, École Polytechnique Fédérale de Lausanne (EPFL), Switzerland, " + "Laboratory of Theory and Simulation of Materials (THEOS). All rights reserved." +) +__license__ = ( + "Non-Commercial, End-User Software License Agreement, see LICENSE.txt file" +) +__version__ = "0.3.0" + +import numpy as np + +from ase import Atoms +from numpy import isscalar +from numbers import Number + + +## Source: http://chemwiki.ucdavis.edu/Reference/Reference_Tables/Atomic_and_Molecular_Properties/A3%3A_Covalent_Radii +# http://dx.doi.org/10.1039/b801115j, checked in paper +_map_atomic_number_radii_cordero = { + 1: 0.31, + 2: 0.28, + 3: 1.28, + 4: 0.96, + 5: 0.84, + 6: 0.76, + 7: 0.71, + 8: 0.66, + 9: 0.57, + 10: 0.58, + 11: 1.66, + 12: 1.41, + 13: 1.21, + 14: 1.11, + 15: 1.07, + 16: 1.05, + 17: 1.02, + 18: 1.06, + 19: 2.03, + 20: 1.76, + 21: 1.7, + 22: 1.6, + 23: 1.53, + 24: 1.39, + 25: 1.39, + 26: 1.32, + 27: 1.26, + 28: 1.24, + 29: 1.32, + 30: 1.22, + 31: 1.22, + 32: 1.2, + 33: 1.19, + 34: 1.2, + 35: 1.2, + 36: 1.16, + 37: 2.2, + 38: 1.95, + 39: 1.9, + 40: 1.75, + 41: 1.64, + 42: 1.54, + 43: 1.47, + 44: 1.46, + 45: 1.42, + 46: 1.39, + 47: 1.45, + 48: 1.44, + 49: 1.42, + 50: 1.39, + 51: 1.39, + 52: 1.38, + 53: 1.39, + 54: 1.4, + 55: 2.44, + 56: 2.15, + 57: 2.07, + 58: 2.04, + 59: 2.03, + 60: 2.01, + 61: 1.99, + 62: 1.98, + 63: 1.98, + 64: 1.96, + 65: 1.94, + 66: 1.92, + 67: 1.92, + 68: 1.89, + 69: 1.9, + 70: 1.87, + 71: 1.87, + 72: 1.75, + 73: 1.7, + 74: 1.62, + 75: 1.51, + 76: 1.44, + 77: 1.41, + 78: 1.36, + 79: 1.36, + 80: 1.32, + 81: 1.45, + 82: 1.46, + 83: 1.48, + 84: 1.4, + 85: 1.5, + 86: 1.5, + 87: 2.6, + 88: 2.21, + 89: 2.15, + 90: 2.06, + 91: 2, + 92: 1.96, + 93: 1.9, + 94: 1.87, + 95: 1.8, + 96: 1.69, +} + +# http://dx.doi.org/10.1002/chem.200901472 +_map_atomic_number_radii_pyykko = { + 1: 0.32, + 2: 0.46, + 3: 1.33, + 4: 1.02, + 5: 0.85, + 6: 0.75, + 7: 0.71, + 8: 0.63, + 9: 0.64, + 10: 0.67, + 11: 1.55, + 12: 1.39, + 13: 1.26, + 14: 1.16, + 15: 1.11, + 16: 1.03, + 17: 0.99, + 18: 0.96, + 19: 1.96, + 20: 1.71, + 21: 1.48, + 22: 1.36, + 23: 1.34, + 24: 1.22, + 25: 1.19, + 26: 1.16, + 27: 1.11, + 28: 1.1, + 29: 1.12, + 30: 1.18, + 31: 1.24, + 32: 1.21, + 33: 1.21, + 34: 1.16, + 35: 1.14, + 36: 1.17, + 37: 2.1, + 38: 1.85, + 39: 1.63, + 40: 1.54, + 41: 1.47, + 42: 1.38, + 43: 1.28, + 44: 1.25, + 45: 1.25, + 46: 1.2, + 47: 1.28, + 48: 1.36, + 49: 1.42, + 50: 1.4, + 51: 1.4, + 52: 1.36, + 53: 1.33, + 54: 1.31, + 55: 2.32, + 56: 1.96, + 57: 1.8, + 58: 1.63, + 59: 1.76, + 60: 1.74, + 61: 1.73, + 62: 1.72, + 63: 1.68, + 64: 1.69, + 65: 1.68, + 66: 1.67, + 67: 1.66, + 68: 1.65, + 69: 1.64, + 70: 1.7, + 71: 1.62, + 72: 1.52, + 73: 1.46, + 74: 1.37, + 75: 1.31, + 76: 1.29, + 77: 1.22, + 78: 1.23, + 79: 1.24, + 80: 1.33, + 81: 1.44, + 82: 1.44, + 83: 1.51, + 84: 1.45, + 85: 1.47, + 86: 1.42, + 87: 2.23, + 88: 2.01, + 89: 1.86, + 90: 1.75, + 91: 1.69, + 92: 1.7, + 93: 1.71, + 94: 1.72, + 95: 1.66, + 96: 1.66, + 97: 1.68, + 98: 1.68, + 99: 1.65, + 100: 1.67, + 101: 1.73, + 102: 1.76, + 103: 1.61, + 104: 1.57, + 105: 1.49, + 106: 1.43, + 107: 1.41, + 108: 1.34, + 109: 1.29, + 110: 1.28, + 111: 1.21, + 112: 1.22, + 113: 1.36, + 114: 1.43, + 115: 1.62, + 116: 1.75, + 117: 1.65, + 118: 1.57, +} + +# source: Batsanov, S. S. "Van der Waals radii of elements." Inorganic materials 37.9 (2001): 871-885. +# For the elements missing there, we used the radii of Pyykko et al (see above) +# to which we added 0.8 Angstrom (rule of thumb mentioned in Batsanov paper, +# working quite well on the elements present in both references) +_map_atomic_number_radii_van_der_Waals = { + 1: 1.1, + 2: 1.26, + 3: 2.2, + 4: 1.9, + 5: 1.8, + 6: 1.7, + 7: 1.6, + 8: 1.55, + 9: 1.5, + 10: 1.4700000000000002, + 11: 2.4, + 12: 2.2, + 13: 2.1, + 14: 2.1, + 15: 1.95, + 16: 1.8, + 17: 1.8, + 18: 1.76, + 19: 2.8, + 20: 2.4, + 21: 2.3, + 22: 2.15, + 23: 2.05, + 24: 2.05, + 25: 2.05, + 26: 2.05, + 27: 2.0, + 28: 2.0, + 29: 2.0, + 30: 2.1, + 31: 2.1, + 32: 2.1, + 33: 2.05, + 34: 1.9, + 35: 1.9, + 36: 1.97, + 37: 2.9, + 38: 2.55, + 39: 2.4, + 40: 2.3, + 41: 2.15, + 42: 2.1, + 43: 2.05, + 44: 2.05, + 45: 2.0, + 46: 2.05, + 47: 2.1, + 48: 2.2, + 49: 2.2, + 50: 2.25, + 51: 2.2, + 52: 2.1, + 53: 2.1, + 54: 2.1100000000000003, + 55: 3.0, + 56: 2.7, + 57: 2.5, + 58: 2.4299999999999997, + 59: 2.56, + 60: 2.54, + 61: 2.5300000000000002, + 62: 2.52, + 63: 2.48, + 64: 2.49, + 65: 2.48, + 66: 2.4699999999999998, + 67: 2.46, + 68: 2.45, + 69: 2.44, + 70: 2.5, + 71: 2.42, + 72: 2.25, + 73: 2.2, + 74: 2.1, + 75: 2.05, + 76: 2.0, + 77: 2.0, + 78: 2.05, + 79: 2.1, + 80: 2.05, + 81: 2.2, + 82: 2.3, + 83: 2.3, + 84: 2.25, + 85: 2.27, + 86: 2.2199999999999998, + 87: 3.0300000000000002, + 88: 2.8099999999999996, + 89: 2.66, + 90: 2.4, + 91: 2.49, + 92: 2.3, + 93: 2.51, + 94: 2.52, + 95: 2.46, + 96: 2.46, + 97: 2.48, + 98: 2.48, + 99: 2.45, + 100: 2.4699999999999998, + 101: 2.5300000000000002, + 102: 2.56, + 103: 2.41, + 104: 2.37, + 105: 2.29, + 106: 2.23, + 107: 2.21, + 108: 2.14, + 109: 2.09, + 110: 2.08, + 111: 2.01, + 112: 2.02, + 113: 2.16, + 114: 2.23, + 115: 2.42, + 116: 2.55, + 117: 2.45, + 118: 2.37, +} + +# source: Alvarez, S., "A cartography of the van der Waals territories", +# Dalton Trans., 2013, 42, 8617 +_map_atomic_number_radii_van_der_Waals_alvarez = { + 1: 1.2, + 2: 1.43, + 3: 2.12, + 4: 1.98, + 5: 1.91, + 6: 1.77, + 7: 1.66, + 8: 1.5, + 9: 1.46, + 10: 1.58, + 11: 2.5, + 12: 2.51, + 13: 2.25, + 14: 2.19, + 15: 1.9, + 16: 1.89, + 17: 1.82, + 18: 1.83, + 19: 2.73, + 20: 2.62, + 21: 2.58, + 22: 2.46, + 23: 2.42, + 24: 2.45, + 25: 2.45, + 26: 2.44, + 27: 2.4, + 28: 2.4, + 29: 2.38, + 30: 2.39, + 31: 2.32, + 32: 2.29, + 33: 1.88, + 34: 1.82, + 35: 1.86, + 36: 2.25, + 37: 3.21, + 38: 2.84, + 39: 2.75, + 40: 2.52, + 41: 2.56, + 42: 2.45, + 43: 2.44, + 44: 2.46, + 45: 2.44, + 46: 2.15, + 47: 2.53, + 48: 2.49, + 49: 2.43, + 50: 2.42, + 51: 2.47, + 52: 1.99, + 53: 2.04, + 54: 2.06, + 55: 3.48, + 56: 3.03, + 57: 2.98, + 58: 2.88, + 59: 2.92, + 60: 2.95, + 62: 2.9, + 63: 2.87, + 64: 2.83, + 65: 2.79, + 66: 2.87, + 67: 2.81, + 68: 2.83, + 69: 2.79, + 70: 2.8, + 71: 2.74, + 72: 2.63, + 73: 2.53, + 74: 2.57, + 75: 2.49, + 76: 2.48, + 77: 2.41, + 78: 2.29, + 79: 2.32, + 80: 2.45, + 81: 2.47, + 82: 2.6, + 83: 2.54, + 89: 2.8, + 90: 2.93, + 91: 2.88, + 92: 2.71, + 93: 2.82, + 94: 2.81, + 95: 2.83, + 96: 3.05, + 97: 3.4, + 98: 3.05, + 99: 2.7, +} + + +class LowDimFinderExc(Exception): + pass + + +class WeirdStructureExc(LowDimFinderExc): + """ + Raised when a weird structure, which the LowDimfinder cannot handle, is found. + """ + + +class GroupOverTwoCellsNoConnectionsExc(LowDimFinderExc): + """ + Raised when the group expands over at least two cells, but no connection is found. + """ + + +class NoRadiiListExc(LowDimFinderExc): + """ + Raised when no valid radii list is defined. + """ + + +class LowDimFinder: + """ + Take an aiida_structure, analyse the structure and return all bonded groups of atoms. + + :param bond_margin: percentage which is added to the + bond length (default: 0). + :param radii_offset: distance (in Angstrom) that is added to + each radius (before using the margin above) (default: 0). + In the end, the criterion which defines if atoms are bonded or not is + that two atoms X1 and X2 are bonded if their distance is smaller than + :math:`[R(X1) + R(X2) + 2\\text{radii_offset}] \\cdot (1 + \\text{bond_margin})`, + where :math:`R(X)` denotes the radius of a species X. + :param radii_source: Can be either ``cordero`` (default), ``pyykko``, + ``vdw`` (for Batsanov van der Walls radii), or ``alvarez`` + (for Alvarez van der Waals radii), the lowdimfinder uses + the radii list from the corresponding papers. + Another possibility is ``custom``, in that case the list of radii has to be + defined with the parameter custom_radii. + :param custom_radii: If the radii_list is set to ``custom``. The dictionary passed with this + parameters is used for the radii. + :param vacuum_space: The space added around the atoms of the structures with lower + dimensionality. (default: 20 angstrom) + :param min_supercell: size of the supercell at which the search for lower dimensionality groups + is started. (default: 3) + :param max_supercell: size at which the search for lower dimensionality groups is stopped. The search + always starts at a supercell size of 2. If the dimensionality of the largest group is 0, + the supercell size is increased by 1, to increase the chance to catch weird structured groups. + (default: 3) + :param full_periodicity: If True, it sets the periodic boundary conditions or the reduced + structures to ``[True, True, True]``, otherwise if False (default) the pbc are set according + the dimensionality. + 0D: ``[False,False,False]``, 1D: ``[False,False,True]``, + 2D: ``[True, True, False]``, 3D: ``[True, True, True]``. + :param rotation: If True, it rotates the reduced 1D structures into the z-axis and 2D structures in + the x-y plane (default: False). + :param orthogonal_axis_2D: If True (default), define the 3rd vector orthogonal to the layer, + otherwise set it in direction of the closest unconnected periodic site. + """ + + def __init__(self, aiida_structure, **kwargs): + + self.params = { + "bond_margin": 0.0, + "radii_offset": 0.0, + "vacuum_space": 20, + "max_supercell": 3, + "min_supercell": 3, + "full_periodicity": False, + "rotation": False, + "radii_source": "cordero", + "orthogonal_axis_2D": True, + "custom_radii": None, + } + self.setup_builder(**kwargs) + + # starting with a min_supercell size and increase it until + # reaching max_supercell + self.supercell_size = self.params["min_supercell"] + + self.aiida_structure = aiida_structure + self.structure = self.aiida_structure + self.n_unit = len(self.structure) # number of atoms in input structure + self.supercell = self.structure.repeat(self.supercell_size) + + self._group_number = None # index of group that is analysed + self._low_dim_index = ( + 0 # index of reduced structure that is currently calculated + ) + self._found_unit_cell_atoms = set( + [] + ) # set of all the input structure atoms that have been attributed to a group + + # independent variables, add those you want to have as output to get_group_data() + self._groups = [] + self._unit_cell_groups = [] + self._atoms_in_one_cell = [] + self._connection_counter = [] + self._connected_positions = [] + self._unconnected_positions = [] + self._vectors = [] + self._shortest_unconnected_vector = [] + + # to be returned by get_group_data() + self._dimensionality = [] + self._chemical_formula = [] + self._positions = [] + self._chemical_symbols = [] + self._cell = [] + self._tags = [] + + self._reduced_ase_structures = [] + self._rotated_structures = [] + self.reduced_aiida_structures = [] + self._3D_structures_with_layer_lattice = [] + + if self.params["radii_source"] == "custom": + self._map_atomic_number_radii_custom = self.params["custom_radii"] + else: + self._map_atomic_number_radii_custom = None + + def setup_builder(self, **kwargs): + """ + Change the builder parameters. + Pass as kwargs the internal parameters to change. + + .. note:: no checks are performed on the validity of the passed parameters. + """ + for k, v in kwargs.items(): + if k in self.params.keys(): + self.params[k] = v + + def _define_groups(self): # pylint: disable=too-many-locals + """ + Return the different groups found in the cell. + """ + # import time + # start_time = time.time() + n_sc = len(self.supercell) + # get distances between all the atoms + coords = self.supercell.get_positions() + + _map_atomic_number_radii = { + "cordero": _map_atomic_number_radii_cordero, + "pyykko": _map_atomic_number_radii_pyykko, + "vdw": _map_atomic_number_radii_van_der_Waals, + "alvarez": _map_atomic_number_radii_van_der_Waals_alvarez, + "custom": self._map_atomic_number_radii_custom, + } + radii_source = self.params["radii_source"] + try: + table = _map_atomic_number_radii[radii_source] + except KeyError: + raise NoRadiiListExc() + atomic_numbers = self.supercell.get_atomic_numbers() + radii = [table[atomic_number] for atomic_number in atomic_numbers] + radii_array = np.array(radii) + self.params["radii_offset"] + pairs = get_bonded_pairs( + coords, radii_array, factor=1.0 + self.params["bond_margin"], use_scipy=True + ) + + set_list = [] + for _ in range(n_sc): + set_list.append(set()) + for i, j in pairs: + set_list[i].add(j) + set_list[j].add(i) + + # At the beginning, set_list contains only first neighbours + # At the end, set_list will contain all neighbours at any + # hopping distance + + # Atoms we already visited + visited = set() + groups = [] + + # I run over all atoms + for atom in range(len(set_list)): # pylint: disable=consider-using-enumerate + # atom not in visited --> new atom group + if atom not in visited: + # fresh set of connected atoms. + # Initially we add the analysed atom. + connected_atoms = set([atom]) + new_neighs = set_list[atom] - connected_atoms + + all_neighbours = connected_atoms.copy() + # Loop stops when no new neighbours are found. + while new_neighs: + # conneted atoms are a set of atoms that we know are already + # connected to the atom (or atom group) + # that we are considering. + connected_atoms.update(new_neighs) + # I update the set of all neighbours of the connected atoms + all_neighbours.update(set.union(*(set_list[a] for a in new_neighs))) + # new_neigbours are the all the neighbours of the connected atoms + # minus the connected atoms. + new_neighs = all_neighbours - connected_atoms + + # Add the connected atoms to the visited list + visited.update(connected_atoms) + groups.append(sorted(connected_atoms)) + + # Returns list of groups sorted by number of atoms + self._groups = sorted(groups, key=len, reverse=True) + + def _define_unit_cell_group(self): + """ + Analyse a group and define which atoms form the unit cell of the reduced group. + unit_cell_sites are the corresponding periodic sites in the unit cell + .. note:: the unit_cell_sites do not have to be equivalent, to the unit cell atoms of + the reduced group. The group can be dispersed over two or more cells + """ + group = self.groups[self._group_number] + + unit_cell_sites = set([]) + reduced_unit_cell_group = set([]) + for i in group: + if i % self.n_unit not in unit_cell_sites: + unit_cell_sites.add(i % self.n_unit) + reduced_unit_cell_group.add(i) + + self._unit_cell_groups.append(sorted(list(reduced_unit_cell_group))) + + self._found_unit_cell_atoms.update(unit_cell_sites) + + if unit_cell_sites == reduced_unit_cell_group: + self._atoms_in_one_cell.append(True) + else: + self._atoms_in_one_cell.append(False) + + def _define_number_of_connections(self): + """ + The connections of a periodic site with its images in the other :math:`N^3-1` + cells are checked. + """ + + # current analysed group + # search the site with the maximum periodic connection in group + group = self.groups[self._group_number] + + # translate all idx back to the unit cell + group_in_unit_cell = [idx % self.n_unit for idx in group] + + # find idx with maximum occurence in list + counter_dict = {} + + for idx in group_in_unit_cell: + counter_dict[idx] = counter_dict.get(idx, 0) + 1 + # make a list of tuples with (counter, idx) + tuple_list = [(counter, idx) for idx, counter in counter_dict.items()] + + # periodic site in unit cell with highest occurrency + counter, periodic_site_in_unit_cell = max(tuple_list) + + self._connection_counter.append(counter) + + # list of connected periodic sites + pos = [] + # list of unconnected periodic sites + unpos = [] + + # get the positions of all connected and unconnected sites + for i in range(self.supercell_size**3): + if i * self.n_unit + periodic_site_in_unit_cell in group: + pos.append( + self.supercell.get_positions()[ + i * self.n_unit + periodic_site_in_unit_cell + ] + ) + else: + unpos.append( + self.supercell.get_positions()[ + i * self.n_unit + periodic_site_in_unit_cell + ] + ) + + self._connected_positions.append(pos) + self._unconnected_positions.append(unpos) + + if ( + self._connection_counter == 1 + and self.supercell_size < self.params["max_supercell"] + and self._low_dim_index == 0 + ): + + self.supercell_size = self.supercell_size + 1 + # reset all the group data + self._groups = [] + self._unit_cell_groups = [] + self._atoms_in_one_cell = [] + self._connection_counter = [] + self._connected_positions = [] + self._vectors = [] + self._dimensionality = [] + + self.supercell = self.structure.repeat(self.supercell_size) + + self._define_number_of_connections() + + def _define_dimensionality(self): # pylint: disable=too-many-locals + """ + Vectors are defined between the connected periodic positions of the group. + The rank + of the matrix made by all those vectors gives the dimensionality of the group. + 0D: an array of three orthonormal vectors is created. + 1D: The shortest vector between periodic sites is taken plus two orthonormal + vectors are created. + 2D: Two shortest linear independent vectors taken plus an orthonormal + vector created + 3D: Takes the vectors of input structure + """ + + def l_are_equals(a, b): + # function to compare lengths + return abs(a - b) <= 1e-5 + + def a_are_equals(a, b): + # function to compare angles (actually, cosines) + return abs(a - b) <= 1e-5 + + # vector to connected positions + vectors = [] + for i in self._get_connected_positions()[self._low_dim_index]: + if i is not self._get_connected_positions()[self._low_dim_index][0]: + vectors.append(i - self._connected_positions[self._low_dim_index][0]) + self._dimensionality.append(np.linalg.matrix_rank(vectors)) + + # vectors of first connected site to unconnected periodic sites + unconnected_vectors = [] + for i in self._get_unconnected_positions()[self._low_dim_index]: + unconnected_vectors.append( + i - self._connected_positions[self._low_dim_index][0] + ) + + if self._dimensionality[self._low_dim_index] == 0: + self._vectors.append( + [np.array([1, 0, 0]), np.array([0, 1, 0]), np.array([0, 0, 1])] + ) + + elif self._dimensionality[self._low_dim_index] == 1: + idx = find_shortest_vector(vectors) + normal_vector1 = np.cross( + vectors[idx], [1, 0, 0] # pylint: disable=invalid-sequence-index + ) + if np.linalg.norm(normal_vector1) < 0.000001: + normal_vector1 = np.cross( + vectors[idx], [0, 1, 0] # pylint: disable=invalid-sequence-index + ) + normal_vector2 = np.cross( + vectors[idx], normal_vector1 # pylint: disable=invalid-sequence-index + ) + normal_vector1 = normal_vector1 / np.linalg.norm(normal_vector1) + normal_vector2 = normal_vector2 / np.linalg.norm(normal_vector2) + self._vectors.append([normal_vector1, normal_vector2, vectors[0]]) + + elif self._dimensionality[self._low_dim_index] == 2: + idx = find_shortest_vector(vectors) + vector1 = vectors.pop(idx) + idx = find_shortest_vector(vectors) + vector2 = vectors.pop(idx) + if self.params["orthogonal_axis_2D"]: + vector3 = np.cross(vector1, vector2) + + while np.linalg.norm(vector3) < 0.0000001: + idx = find_shortest_vector(vectors) + vector2 = vectors.pop(idx) + vector3 = np.cross(vector1, vector2) + + else: + # to prevent case where in a supercell the two shortest vectors + # are linearly dependent + while np.linalg.norm(np.cross(vector1, vector2)) < 0.0000001: + idx = find_shortest_vector(vectors) + vector2 = vectors.pop(idx) + + # do same for 3rd vector !!! + idx = find_shortest_vector(unconnected_vectors) + vector3 = unconnected_vectors.pop(idx) + + while ( + np.linalg.norm(np.dot(np.cross(vector1, vector2), vector3)) + < 0.0000001 + ): + idx = find_shortest_vector(unconnected_vectors) + vector3 = unconnected_vectors.pop(idx) + + self._shortest_unconnected_vector.append(vector3) + + vector3 = vector3 / np.linalg.norm(vector3) + + # check the hexagonal case + norm_vector1 = np.linalg.norm(vector1) + norm_vector2 = np.linalg.norm(vector2) + cosphi = np.dot(vector1, vector2) / norm_vector1 / norm_vector2 + if l_are_equals(norm_vector1, norm_vector2) and a_are_equals(cosphi, 0.5): + # change the lattice vectors to get a "standard" hexagonal cell + # (compliant with the KpointsData class) + vector2 = vector2 - vector1 + + self._vectors.append([vector1, vector2, vector3]) + + elif self._dimensionality[self._low_dim_index] == 3: + self._vectors.append(self.structure.cell) + + def _define_reduced_ase_structure( + self, + ): # pylint: disable=too-many-branches,too-many-statements,too-many-locals + """ + Append the reduced group of atoms corresponding to the currently analysed group + to the list of reduced ase structures. + """ + min_z = None + max_z = None + min_x = None + min_y = None + max_x = None + max_y = None + + symbols = [] + positions = [] + tags = [] + + # 3 Dimensional --> return input structure + if self._get_dimensionality()[self._low_dim_index] == 3: + + self._reduced_ase_structures.append(self.structure) + self._chemical_formula.append(self.structure.get_chemical_formula()) + self._positions.append(self.structure.get_positions().tolist()) + self._chemical_symbols.append(self.structure.get_chemical_symbols()) + self._cell.append(self.structure.cell.tolist()) + self._tags.append(self.structure.get_tags().tolist()) + + self._vectors[self._low_dim_index][0] = self.structure.cell.tolist()[0] + self._vectors[self._low_dim_index][1] = self.structure.cell.tolist()[1] + self._vectors[self._low_dim_index][2] = self.structure.cell.tolist()[2] + # has to be called in case a structure is made of a 3D compound + # and a lower D compound. + self._get_unit_cell_groups() + + return + + if self._get_dimensionality()[self._low_dim_index] == 2: + # get positions and chemical symbols + for i in self._get_unit_cell_groups()[self._low_dim_index]: + positions.append( + [ + self.supercell.positions[i][0], + self.supercell.positions[i][1], + self.supercell.positions[i][2], + ] + ) + symbols.append(self.supercell.get_chemical_symbols()[i]) + tags.append(self.supercell.get_tags()[i]) + # min and max z + normal_vec_z = np.cross( + self._vectors[self._low_dim_index][0], + self._vectors[self._low_dim_index][1], + ) + current_z = np.dot( + positions[-1], normal_vec_z / np.linalg.norm(normal_vec_z) + ) + + if min_z is None or min_z > current_z: + min_z = current_z + if max_z is None or max_z < current_z: + max_z = current_z + + positions = np.array(positions) + + if self.params["orthogonal_axis_2D"]: + + positions = ( + positions + - (min_z - self.params["vacuum_space"] / 2) + * self._vectors[self._low_dim_index][2] + ) + self._vectors[self._low_dim_index][2] = self._vectors[ + self._low_dim_index + ][2] * (self.params["vacuum_space"] + max_z - min_z) + + else: + # in order to keep the distance between the layers equal to the vacuum space + # get vector orthogonal to layer + normal_vec = np.cross( + self._vectors[self._low_dim_index][0], + self._vectors[self._low_dim_index][1], + ) + + # cosinus between third vector and normal vector + cos = abs( + np.dot(normal_vec, self._vectors[self._low_dim_index][2]) + / np.linalg.norm(normal_vec) + / np.linalg.norm(self._vectors[self._low_dim_index][2]) + ) + + # translate the lengths into direction of the third vector + vacuum_space_in_third_direction = ( + float(self.params["vacuum_space"]) / cos + ) + max_z = max_z / cos + min_z = min_z / cos + + positions = ( + positions + - (min_z - vacuum_space_in_third_direction / 2) + * self._vectors[self._low_dim_index][2] + ) + self._vectors[self._low_dim_index][2] = self._vectors[ + self._low_dim_index + ][2] * (vacuum_space_in_third_direction + max_z - min_z) + + pbc = [True, True, False] + + elif self._get_dimensionality()[self._low_dim_index] == 1: + for i in self._get_unit_cell_groups()[self._low_dim_index]: + positions.append( + [ + self.supercell.positions[i][0], + self.supercell.positions[i][1], + self.supercell.positions[i][2], + ] + ) + symbols.append(self.supercell.get_chemical_symbols()[i]) + tags.append(self.supercell.get_tags()[i]) + + if min_x is None or min_x > np.dot( + positions[-1], self._vectors[self._low_dim_index][0] + ): + min_x = np.dot(positions[-1], self._vectors[self._low_dim_index][0]) + if max_x is None or max_x < np.dot( + positions[-1], self._vectors[self._low_dim_index][0] + ): + max_x = np.dot(positions[-1], self._vectors[self._low_dim_index][0]) + if min_y is None or min_y > np.dot( + positions[-1], self._vectors[self._low_dim_index][1] + ): + min_y = np.dot(positions[-1], self._vectors[self._low_dim_index][1]) + if max_y is None or max_y < np.dot( + positions[-1], self._vectors[self._low_dim_index][1] + ): + max_y = np.dot(positions[-1], self._vectors[self._low_dim_index][1]) + + positions = np.array(positions) + positions = ( + positions + - (min_x - self.params["vacuum_space"] / 2) + * self._vectors[self._low_dim_index][0] + ) + positions = ( + positions + - (min_y - self.params["vacuum_space"] / 2) + * self._vectors[self._low_dim_index][1] + ) + + self._vectors[self._low_dim_index][0] = self._vectors[self._low_dim_index][ + 0 + ] * (self.params["vacuum_space"] + max_x - min_x) + self._vectors[self._low_dim_index][1] = self._vectors[self._low_dim_index][ + 1 + ] * (self.params["vacuum_space"] + max_y - min_y) + + pbc = [False, False, True] + + elif self._get_dimensionality()[self._low_dim_index] == 0: + + for i in self._get_unit_cell_groups()[self._low_dim_index]: + positions.append( + [ + self.supercell.positions[i][0], + self.supercell.positions[i][1], + self.supercell.positions[i][2], + ] + ) + symbols.append(self.supercell.get_chemical_symbols()[i]) + tags.append(self.supercell.get_tags()[i]) + + if min_x is None or min_x > np.dot( + positions[-1], self._vectors[self._low_dim_index][0] + ): + min_x = np.dot(positions[-1], self._vectors[self._low_dim_index][0]) + if max_x is None or max_x < np.dot( + positions[-1], self._vectors[self._low_dim_index][0] + ): + max_x = np.dot(positions[-1], self._vectors[self._low_dim_index][0]) + if min_y is None or min_y > np.dot( + positions[-1], self._vectors[self._low_dim_index][1] + ): + min_y = np.dot(positions[-1], self._vectors[self._low_dim_index][1]) + if max_y is None or max_y < np.dot( + positions[-1], self._vectors[self._low_dim_index][1] + ): + max_y = np.dot(positions[-1], self._vectors[self._low_dim_index][1]) + if min_z is None or min_z > np.dot( + positions[-1], self._vectors[self._low_dim_index][2] + ): + min_z = np.dot(positions[-1], self._vectors[self._low_dim_index][2]) + if max_z is None or max_z < np.dot( + positions[-1], self._vectors[self._low_dim_index][2] + ): + max_z = np.dot(positions[-1], self._vectors[self._low_dim_index][2]) + + pbc = [False, False, False] + + positions = np.array(positions) + positions = ( + positions + - (min_x - self.params["vacuum_space"] / 2) + * self._vectors[self._low_dim_index][0] + ) + positions = ( + positions + - (min_y - self.params["vacuum_space"] / 2) + * self._vectors[self._low_dim_index][1] + ) + positions = ( + positions + - (min_z - self.params["vacuum_space"] / 2) + * self._vectors[self._low_dim_index][2] + ) + + self._vectors[self._low_dim_index][0] = self._vectors[self._low_dim_index][ + 0 + ] * (self.params["vacuum_space"] + max_x - min_x) + self._vectors[self._low_dim_index][1] = self._vectors[self._low_dim_index][ + 1 + ] * (self.params["vacuum_space"] + max_y - min_y) + self._vectors[self._low_dim_index][2] = self._vectors[self._low_dim_index][ + 2 + ] * (self.params["vacuum_space"] + max_z - min_z) + + else: + raise WeirdStructureExc("No dimensionality") + + if self.params["full_periodicity"]: + pbc = [True, True, True] + + reduced_ase_structure = Atoms( + cell=self._vectors[self._low_dim_index], + pbc=pbc, + symbols=symbols, + positions=positions, + tags=tags, + ) + reduced_ase_structure.set_positions( + reduced_ase_structure.get_positions(wrap=True) + ) + + if self.params["rotation"]: + + original_struc = self.structure.copy() + rotated_original_structure, reduced_ase_structure = self._rotate_structures( + original_struc, reduced_ase_structure + ) + + self._rotated_structures.append(rotated_original_structure) + + self._reduced_ase_structures.append(reduced_ase_structure) + self._chemical_formula.append(reduced_ase_structure.get_chemical_formula()) + self._positions.append(reduced_ase_structure.get_positions().tolist()) + self._chemical_symbols.append(reduced_ase_structure.get_chemical_symbols()) + self._cell.append(reduced_ase_structure.cell.tolist()) + self._tags.append(reduced_ase_structure.get_tags().tolist()) + + def _define_reduced_ase_structures(self): + """ + Define the ASE structures of groups of atoms, starting with the largest group + until all atoms of the original structure are assigned to an ASE structure. + .. note:: Groups containing atoms corresponding to sites in the original structure + that have already been assigned to an ASE structure, are skipped. Therefore + every atom in the original structure will be only present in one ASE structure. + """ + + # all the atoms of the input structure, which have to be allocated to groups + test = set(range(self.n_unit)) + + for idx, group in enumerate(self.groups): + # break as soon as all atoms of the input structure are in a group + if self._found_unit_cell_atoms == test: + break + self._group_number = idx + + # set of the periodic sites of the input structure + unit_cell_group = {i % self.n_unit for i in group} + + # if the atoms corresponding to the group are only a periodic + # replica of an existing low dimensionality + # group, the group is skipped and the next one is analysed.. + if unit_cell_group.issubset(self._found_unit_cell_atoms): + continue + self._define_reduced_ase_structure() + self._low_dim_index = self._low_dim_index + 1 + self._low_dim_index = self._low_dim_index - 1 + + def _define_reduced_aiida_structures(self): + """ + Converts all ASE structures into AiiDA structures + """ + # from aiida.plugins import DataFactory + if not self._reduced_ase_structures: + self._define_reduced_ase_structures() + # S = DataFactory("structure") + + self.reduced_aiida_structures = [] + for i in self._reduced_ase_structures: + self.reduced_aiida_structures.append(i) + + def _rotate_structures(self, originalstruc, asestruc): + """ + Rotates the third vector in z-axis and the first vector in x-axis. + Layer (2D) into xy-plane, wire (1D) into z-direction + + :param originalstruc: unrotated original ASE structure + :param asestruc: unrotated reduced ASE structure + :return: rotated ASE original and reduced structures + """ + if self.params["orthogonal_axis_2D"]: + originalstruc.rotate( + v=asestruc.cell[2], a=[0, 0, 1], center=(0, 0, 0), rotate_cell=True + ) + asestruc.rotate( + v=asestruc.cell[2], a=[0, 0, 1], center=(0, 0, 0), rotate_cell=True + ) + + else: + # turn the normal vector into z direction + normal_vec = np.cross(asestruc.cell[0], asestruc.cell[1]) + originalstruc.rotate( + v=normal_vec, a=[0, 0, 1], center=(0, 0, 0), rotate_cell=True + ) + asestruc.rotate( + v=normal_vec, a=[0, 0, 1], center=(0, 0, 0), rotate_cell=True + ) + # turn it into positive z direction (if needed) + if asestruc.cell[2][2] < 0.0: + originalstruc.rotate( + v=[0, 0, np.sign(asestruc.cell[2][2])], + a=[0, 0, 1], + center=(0, 0, 0), + rotate_cell=True, + ) + asestruc.rotate( + v=[0, 0, np.sign(asestruc.cell[2][2])], + a=[0, 0, 1], + center=(0, 0, 0), + rotate_cell=True, + ) + + # finally rotate the first cell vector of asestruc into x (this is a + # rotation inside the x-y plane so it will not impact its normal) + originalstruc.rotate( + v=asestruc.cell[0], a=[1, 0, 0], center=(0, 0, 0), rotate_cell=True + ) + asestruc.rotate( + v=asestruc.cell[0], a=[1, 0, 0], center=(0, 0, 0), rotate_cell=True + ) + + return originalstruc, asestruc + + def _rotate_ase_structure(self, asestruc): + """ + Rotates the third vector in z-axis and the first vector in x-axis. + Layer (2D) into xy-plane, wire (1D) into z-direction + + :param asestruc: unrotated reduced ASE structure + :return: rotated ASE structure + """ + if self.params["orthogonal_axis_2D"]: + asestruc.rotate( + v=asestruc.cell[2], a=[0, 0, 1], center=(0, 0, 0), rotate_cell=True + ) + + else: + # turn the normal vector into z direction + normal_vec = np.cross(asestruc.cell[0], asestruc.cell[1]) + asestruc.rotate( + v=normal_vec, a=[0, 0, 1], center=(0, 0, 0), rotate_cell=True + ) + # turn it into positive z direction (if needed) + if asestruc.cell[2][2] < 0.0: + asestruc.rotate( + v=[0, 0, np.sign(asestruc.cell[2][2])], + a=[0, 0, 1], + center=(0, 0, 0), + rotate_cell=True, + ) + + # finally rotate the first cell vector of asestruc into x (this is a + # rotation inside the x-y plane so it will not impact its normal) + asestruc.rotate( + v=asestruc.cell[0], a=[1, 0, 0], center=(0, 0, 0), rotate_cell=True + ) + + return asestruc + + def _get_dimensionality(self): + """ + Add dimensionality and normalized lattice vectors of the currently + analyzed group to the corresponding lists + + Return the dimensionality list. + """ + if len(self._dimensionality) <= self._low_dim_index: + self._define_dimensionality() + return self._dimensionality + + @property + def groups(self): + """ + Return all atoms in the supercell which are only held together by vdW + forces as separate groups + + :return: groups of bonded atoms + """ + if not self._groups: + self._define_groups() + return self._groups + + def _get_unit_cell_groups(self): + """ + + Return atom ids used to build the reduced ase structures. + """ + + if len(self._unit_cell_groups) <= self._low_dim_index: + self._define_unit_cell_group() + + return self._unit_cell_groups + + def _get_number_of_connections(self): + """ + Return the number of connected periodic positions in the supercell. + """ + if len(self._connection_counter) <= self._low_dim_index: + self._define_number_of_connections() + return self._connection_counter + + def _get_connected_positions(self): + """ + Return the connected periodic positions in the supercell. + """ + if len(self._connection_counter) <= self._low_dim_index: + self._define_number_of_connections() + return self._connected_positions + + def _get_unconnected_positions(self): + """ + Return the unconnected periodic positions in the supercell. + """ + if len(self._connection_counter) <= self._low_dim_index: + self._define_number_of_connections() + return self._unconnected_positions + + def _get_reduced_ase_structures(self): + """ + Return a list with all the lower dimensionality structures, as ASE atoms. + """ + if not self._reduced_ase_structures: + self._define_reduced_ase_structures() + return self._reduced_ase_structures + + def get_reduced_aiida_structures(self): + """ + Return a list with all the lower dimensionality structures found, + as AiiDA structures. + """ + if not self.reduced_aiida_structures: + self._define_reduced_aiida_structures() + return self.reduced_aiida_structures + + def get_group_data(self): + """ + Return a dictionary with list of the dimensionality, chemical_formula, + cell parameters, positions and chemical symbols of the atoms of the + extracted structures. + """ + _ = self.get_reduced_aiida_structures() + return { + "dimensionality": self._dimensionality, + "chemical_formula": self._chemical_formula, + "positions": self._positions, + "chemical_symbols": self._chemical_symbols, + "unit_cell_ids":self._unit_cell_groups, + "cell": self._cell, + "tags": self._tags, + } + + def _get_rotated_structures(self): + """ + Return rotated original structures if rotation is true. + """ + + _ = self.get_reduced_aiida_structures() + + # from aiida.plugins import DataFactory + # S = DataFactory("structure") + + rotated_aiida_structures = [] + + for asestruc in self._rotated_structures: + rotated_aiida_structures.append(asestruc) + + return rotated_aiida_structures + + def _get_3D_structures_with_layer_lattice(self): + """ + Return list of 3D structures with the same lattice as the layer, + for all 2D structures found. + """ + if not self._3D_structures_with_layer_lattice: + self._define_3D_structures_with_layer_lattice() + + return self._3D_structures_with_layer_lattice + + def _define_3D_structures_with_layer_lattice(self): + """ + Define 3D structures with the same lattice as the layers for all + 2D reduced structures found. + """ + from aiida.plugins import ( # pylint:disable=import-error,import-outside-toplevel + DataFactory, # pylint:disable=import-error,import-outside-toplevel + ) + + S = DataFactory("structure") + _ = self.get_reduced_aiida_structures() + + # 2D structure counter + counter = 0 + + for idx, dim in enumerate(self._dimensionality): + # build 3D structure with layer lattice + if dim == 2 and not self.params["orthogonal_axis_2D"]: + + positions = self.structure.get_positions() + positions = positions - self.structure.cell[2] + + # get the two first vectors from the reduced ase structure and the shortest unconnected vector. + asestruc = Atoms( + cell=[ + self._vectors[idx][0], + self._vectors[idx][1], + self._shortest_unconnected_vector[counter], + ], + positions=self.structure.get_positions(), + symbols=self.structure.get_chemical_symbols(), + tags=self.structure.get_tags(), + pbc=[True, True, True], + ) + + asestruc.set_positions(asestruc.get_positions(wrap=True)) + + if int(round(asestruc.get_volume() / self.structure.get_volume())) != 1: + # we need to do something else if the volume changes (it's a supercell) + # -> we use the supercell, wrap everything in, and remove + # the overlapping stuff + the_asestruc = Atoms( + cell=[ + self._vectors[idx][0], + self._vectors[idx][1], + self._shortest_unconnected_vector[counter], + ], + positions=self.supercell.get_positions(), + symbols=self.supercell.get_chemical_symbols(), + pbc=[True, True, True], + ) + + pos_and_symbols = zip( + the_asestruc.get_scaled_positions().tolist(), + the_asestruc.get_chemical_symbols(), + ) + the_pos, the_symbols = zip(*objects_set(pos_and_symbols)) + asestruc = Atoms( + cell=[ + self._vectors[idx][0], + self._vectors[idx][1], + self._shortest_unconnected_vector[counter], + ], + scaled_positions=the_pos, + symbols=the_symbols, + pbc=[True] * 3, + ) + + if self.params["rotation"] is True: + asestruc = self._rotate_ase_structure(asestruc) + + counter = counter + 1 + + self._3D_structures_with_layer_lattice.append(S(ase=asestruc)) + + +def find_shortest_vector(array): + """ + Takes an array of vectors and finds the shortest one. + + :param array: array of vectors + :return idx: the index of the shortest vector in the array + """ + idx = np.array([np.linalg.norm(vector) for vector in array]).argmin() + return idx + + +def get_bonded_pairs( + coords, radii_array, factor=1.0, use_scipy=True +): # pylint: disable=too-many-locals + """ + Get the list of pairs of connected atoms in the cell. + (i, j) in the list means atom i and atom j are connected. + + :param coords: n*3 array with the positions to compare (n=number of + atoms in the cell) + :param radii_array: array of length n with the radii for + each of the n atoms) + :param factor: constant factor multiplied to all the bonds + :param use_scipy: True to use the optimized scipy.spatial.cKDTree + + .. note:: if use_scipy is True, it requires scipy version >= 0.16.0 + + :return: pairs of connected atoms + """ + if use_scipy: + try: + import scipy.spatial # pylint: disable=import-outside-toplevel + + tree = scipy.spatial.cKDTree(coords) + except (ImportError, AttributeError): + raise ValueError( + "To use scipy and speed-up the lowdimfinder, you " + "need scipy with a version >= 0.16.0" + ) + max_radius = np.max(radii_array) * factor + + # prefiltering the pairs by twice the maximum radius + pairs = list(tree.query_pairs(2.0 * max_radius)) + + # among those find the pairs that are really connected + if len(pairs) > 0: + atom1_idx, atom2_idx = zip(*pairs) # list of atoms within those pairs + coords1 = coords[np.array(atom1_idx)] + coords2 = coords[np.array(atom2_idx)] + # radii of these atoms + rad1 = radii_array[np.array(atom1_idx)] * factor + rad2 = radii_array[np.array(atom2_idx)] * factor + pairs_idx = np.nonzero( + ((coords1 - coords2) ** 2).sum(axis=1) - (rad1 + rad2) ** 2 < 0.0 + ) + return np.array(pairs)[pairs_idx[0]] + return [] + + # Get the squared_distances between to lists of cartesian coordinates. + # matrix(i, j) contains the squared distance between atom i and atom j. + matrix = (coords[:, None, :] - coords[None, :, :]) ** 2 + dist_squared = np.sum(matrix, axis=2) + radii_matrix = (radii_array + radii_array[:, None]) * factor + pairs = zip(*np.nonzero(dist_squared - radii_matrix**2 < 0)) + return [pair for pair in pairs if pair[0] < pair[1]] + + +def numbers_are_equal(a, b, epsilon=1e-14): + """ + Compare two numbers a and b within epsilon. + :param a: float, int or any scalar + :param b: float, int or any scalar + :param epsilon: threshold above which a is considered different from b + :return: boolean + """ + return abs(a - b) < epsilon + + +def scalars_are_equal(a, b, **kwargs): + """ + Compare two objects of any type, except list, arrays or dictionaries. + Numbers are compared thanks to the ``numbers_are_equal`` function. + :param a: any scalar object + :param b: any scalar object + :param kwargs: parameters passed to the function numbers_are_equal + :return: boolean + :raise: NonscalarError if either a or b is a list, an array or a dictionary + """ + if isscalar(a) and isscalar(b): + if isinstance(a, Number): + return isinstance(b, Number) and numbers_are_equal(a, b, **kwargs) + return a == b + if (a is None) or (b is None): + return a == b + raise TypeError("a and b must be scalars") + + +def objects_set(objects, **kwargs): + """ + Return a set made of objects compared between them using the + 'objects_are_equal' function + :param objects: an iterable containing any kind of objects that can + be compared using the 'objects_are_equal' function (list, dict, etc.) + :param kwargs: additional keyword arguments to be passed to the + 'objects_are_equal' function (e.g. precision for floating point number) + :return: a set of non-equal objects + """ + the_set = [] + for obj in objects: + if len([o for o in the_set if objects_are_equal(obj, o, **kwargs)]) == 0: + the_set.append(obj) + + return the_set + + +def objects_are_equal(obj1, obj2, **kwargs): + """ + Recursive function. + Return True if obj1 is the same as obj2. Scalars are + compared using the function ``scalars_are_equal``. + Handles strings, floats, ints, booleans, as well as lists, arrays and + dictionaries containing such objects (possibly nested). + :param obj1: any object + :param obj2: any object + :param kwargs: parameters passed to the function scalars_are_equal + :return: boolean + """ + if isinstance(obj1, dict): + if not isinstance(obj2, dict): + return False + obj1_keys = sorted(obj1.keys()) + obj2_keys = sorted(obj2.keys()) + if not objects_are_equal(obj1_keys, obj2_keys, **kwargs): + return False + if not objects_are_equal( + [obj1[k] for k in obj1_keys], [obj2[k] for k in obj2_keys], **kwargs + ): + return False + return True + + if isinstance(obj1, (list, np.ndarray, tuple)): + if not isinstance(obj2, (list, np.ndarray, tuple)): + return False + if len(obj1) != len(obj2): + return False + for e1, e2 in zip(obj1, obj2): + if np.isscalar(e1): + if not np.isscalar(e2): + return False + if not scalars_are_equal(e1, e2, **kwargs): + return False + else: + if not objects_are_equal(e1, e2, **kwargs): + return False + return True + + try: + return scalars_are_equal(obj1, obj2, **kwargs) + except TypeError: + raise TypeError("Type of obj1 and obj2 not recognized") diff --git a/surfaces_tools/widgets/nn.py b/surfaces_tools/widgets/nn.py new file mode 100644 index 0000000..c060fdd --- /dev/null +++ b/surfaces_tools/widgets/nn.py @@ -0,0 +1,468 @@ +import copy +import itertools + +import aiidalab_widgets_base as awb +import ase +import numpy as np +import traitlets as tr +from ase import neighborlist +from scipy import signal, sparse, spatial + + +def to_ranges(iterable): + iterable = sorted(set(iterable)) + for _key, group in itertools.groupby(enumerate(iterable), lambda t: t[1] - t[0]): + lgroup = list(group) + yield lgroup[0][1], lgroup[-1][1] + + +def mol_ids_range(ismol): + # Shifts the list by +1. + range_string = "" + shifted_list = [i + 1 for i in ismol] + ranges = list(to_ranges(shifted_list)) + for i in range(len(ranges)): + if ranges[i][1] > ranges[i][0]: + range_string += f"{ranges[i][0]}..{ranges[i][1]}" + else: + range_string += f"{ranges[i][0]} " + return range_string + + +def conne_matrix(atoms): + cutoff = neighborlist.natural_cutoffs(atoms) + neighbor_list = neighborlist.NeighborList( + cutoff, self_interaction=False, bothways=False + ) + neighbor_list.update(atoms) + + return neighbor_list.get_connectivity_matrix() + + +def clusters(matrix): + nclusters, idlist = sparse.csgraph.connected_components(matrix) + return nclusters, [np.where(idlist == i)[0].tolist() for i in range(nclusters)] + + +def molecules(ismol, atoms): + if ismol: + nmols, ids = clusters(conne_matrix(atoms[ismol])) + return [[ismol[i] for i in ids[j]] for j in range(nmols)] + return [] + + +class StructureAnalyzer(tr.HasTraits): + structure = tr.Instance(ase.Atoms, allow_none=True) + details = tr.Dict() + + def __init__(self, only_sys_type=False, who_called_me="Boh"): + self.only_sys_type = only_sys_type + self.who_called_me = who_called_me + super().__init__() + + def gaussian(self, x, sig): + return ( + 1.0 + / (sig * np.sqrt(2.0 * np.pi)) + * np.exp(-np.power(x, 2.0) / (2 * np.power(sig, 2.0))) + ) + + def boxfilter(self, x, thr): + return np.asarray([1 if i < thr else 0 for i in x]) + + def get_types(self, frame, thr): # Piero Gasparotto + # classify the atmos in: + # 0=molecule + # 1=slab atoms + # 2=adatoms + # 3=hydrogens on the surf + # 5=unknown + # 6=metalating atoms + # frame=ase frame + # thr=threashold in the histogram for being considered a surface layer + nat = frame.get_number_of_atoms() + + # all atom types set to 5 (unknown) + atype = np.zeros(nat, dtype=np.int16) + 5 + area = frame.cell[0][0] * frame.cell[1][1] + minz = np.min(frame.positions[:, 2]) + maxz = np.max(frame.positions[:, 2]) + + if maxz - minz < 1.0: + maxz += (1.0 - (maxz - minz)) / 2 + minz -= (1.0 - (maxz - minz)) / 2 + + # Which values should we use? + sigma = 0.2 + peak_rel_height = 0.5 + layer_tol = 1.0 * sigma + + # Quick estimate number atoms in a layer. + nbins = int(np.ceil((maxz - minz) / 0.15)) + hist, bin_edges = np.histogram(frame.positions[:, 2], density=False, bins=nbins) + max_atoms_in_a_layer = max(hist) + + lbls = frame.get_chemical_symbols() + n_intervals = int(np.ceil((maxz - minz + 3 * sigma) / (0.1 * sigma))) + z_values = np.linspace(minz - 3 * sigma, maxz + 3 * sigma, n_intervals) # 1000 + atoms_z_pos = frame.positions[:, 2] + + # Generate 2d array to apply the gaussian on. + z_v_exp, at_z_exp = np.meshgrid(z_values, atoms_z_pos) + arr_2d = z_v_exp - at_z_exp + atomic_density = np.sum(self.gaussian(arr_2d, sigma), axis=0) + + peaks = signal.find_peaks( + atomic_density, + height=None, + threshold=None, + distance=None, + prominence=None, + width=None, + wlen=None, + rel_height=peak_rel_height, + ) + layersg = z_values[peaks[0].tolist()] + len(layersg) + layersg[-1] + + # Check top and bottom layers should be documented better. + found_top_surf = False + while not found_top_surf: + iz = layersg[-1] + two_d_atoms = [ + frame.positions[i, 0:2] + for i in range(nat) + if np.abs(frame.positions[i, 2] - iz) < layer_tol + ] + coverage = 0 + if len(two_d_atoms) > max_atoms_in_a_layer / 4: + hull = spatial.ConvexHull(two_d_atoms) + coverage = hull.volume / area + if coverage > 0.3: + found_top_surf = True + else: + layersg = layersg[0:-1] + + found_bottom_surf = False + while not found_bottom_surf: + iz = layersg[0] + two_d_atoms = [ + frame.positions[i, 0:2] + for i in range(nat) + if np.abs(frame.positions[i, 2] - iz) < layer_tol + ] + coverage = 0 + if len(two_d_atoms) > max_atoms_in_a_layer / 4: + hull = spatial.ConvexHull(two_d_atoms) + coverage = hull.volume / area + if coverage > 0.3 and len(two_d_atoms) > max_atoms_in_a_layer / 4: + found_bottom_surf = True + else: + layersg = layersg[1:] + + bottom_z = layersg[0] + top_z = layersg[-1] + + # Check if there is a bottom layer of H. + found_layer_of_h = True + for i in range(nat): + iz = frame.positions[i, 2] + if iz > bottom_z - layer_tol and iz < bottom_z + layer_tol: + if lbls[i] == "H": + atype[i] = 3 + else: + found_layer_of_h = False + break + if found_layer_of_h: + layersg = layersg[1:] + # bottom_z=layersg[0] + + layers_dist = [] + iprev = layersg[0] + for inext in layersg[1:]: + layers_dist.append(abs(iprev - inext)) + iprev = inext + + for i in range(nat): + iz = frame.positions[i, 2] + if iz > bottom_z - layer_tol and iz < top_z + layer_tol: + if not (atype[i] == 3 and found_layer_of_h): + atype[i] = 1 + else: + if np.min([np.abs(iz - top_z), np.abs(iz - bottom_z)]) < np.max( + layers_dist + ): + if not (atype[i] == 3 and found_layer_of_h): + atype[i] = 2 + + # assign the other types + metalatingtypes = ("Au", "Ag", "Cu", "Ni", "Co", "Zn", "Mg") + moltypes = ("H", "N", "B", "O", "C", "F", "S", "Br", "I", "Cl") + possible_mol_atoms = [ + i for i in range(nat) if atype[i] == 2 and lbls[i] in moltypes + ] + possible_mol_atoms += [i for i in range(nat) if atype[i] == 5] + + # Identify separate molecules. + all_molecules = [] + if len(possible_mol_atoms) > 0: + fragments = molecules(possible_mol_atoms, frame) + all_molecules = copy.deepcopy(fragments) + + # Remove isolated atoms. + for frag in fragments: + if len(frag) == 1: + all_molecules.remove(frag) + else: + for atom in frag: + if lbls[atom] in metalatingtypes: + atype[atom] = 6 + else: + atype[atom] = 0 + + return atype, layersg, all_molecules + + def all_connected_to(self, id_atom, atoms, exclude): + cov_radii = [ase.data.covalent_radii[a.number] for a in atoms] + + atoms.set_pbc([False, False, False]) + nl_no_pbc = neighborlist.NeighborList( + cov_radii, bothways=True, self_interaction=False + ) + nl_no_pbc.update(atoms) + atoms.set_pbc([True, True, True]) + + tofollow = [] + followed = [] + isconnected = [] + tofollow.append(id_atom) + isconnected.append(id_atom) + while len(tofollow) > 0: + indices, offsets = nl_no_pbc.get_neighbors(tofollow[0]) + indices = list(indices) + followed.append(tofollow[0]) + for i in indices: + if (i not in isconnected) and (atoms[i].symbol not in exclude): + tofollow.append(i) + isconnected.append(i) + for i in followed: + if i in tofollow: # do not remove this check + tofollow.remove(i) + + return isconnected + + def string_range_to_list(self, a): + singles = [int(s) - 1 for s in a.split() if s.isdigit()] + ranges = [r for r in a.split() if ".." in r] + for r in ranges: + t = r.split("..") + to_add = [i - 1 for i in range(int(t[0]), int(t[1]) + 1)] + singles += to_add + return sorted(singles) + + @tr.observe("structure") + def _observe_structure(self, _=None): + with self.hold_trait_notifications(): + self.details = self.analyze() + + def analyze(self): + if self.structure is None: + return {} + + atoms = self.structure + sys_size = np.ptp(atoms.positions, axis=0) + no_cell = ( + atoms.cell[0][0] < 0.1 or atoms.cell[1][1] < 0.1 or atoms.cell[2][2] < 0.1 + ) + if no_cell: + # set bounding box as cell + atoms.cell = sys_size + 10 + + atoms.set_pbc([True, True, True]) + + total_charge = np.sum(atoms.get_atomic_numbers()) + bottom_h = [] + adatoms = [] + bulkatoms = [] + wireatoms = [] + metalatings = [] + unclassified = [] + slabatoms = [] + slab_layers = [] + all_molecules = None + is_a_bulk = False + is_a_molecule = False + is_a_wire = False + + spins_up = [the_a.index for the_a in atoms if the_a.tag == 1] + spins_down = [the_a.index for the_a in atoms if the_a.tag == 2] + other_tags = [the_a.index for the_a in atoms if the_a.tag > 2] + + # Check if there is vacuum otherwise classify as bulk and skip. + vacuum_x = sys_size[0] + 4 < atoms.cell[0][0] + vacuum_y = sys_size[1] + 4 < atoms.cell[1][1] + vacuum_z = sys_size[2] + 4 < atoms.cell[2][2] + + # Do not use a set in the following line list(set(atoms.get_chemical_symbols())) + # Need ALL atoms and elements for spin guess and for cost calculation + all_elements = atoms.get_chemical_symbols() + + summary = "" + cases = [] + if len(spins_up) > 0: + summary += "spins_up: " + awb.utils.list_to_string_range(spins_up) + "\n" + if len(spins_down) > 0: + summary += ( + "spins_down: " + awb.utils.list_to_string_range(spins_down) + "\n" + ) + if len(other_tags) > 0: + summary += ( + "other_tags: " + awb.utils.list_to_string_range(other_tags) + "\n" + ) + if (not vacuum_z) and (not vacuum_x) and (not vacuum_y): + is_a_bulk = True + sys_type = "Bulk" + cases = ["b"] + summary += "Bulk contains: \n" + slabatoms = list(range(len(atoms))) + bulkatoms = slabatoms + + if vacuum_x and vacuum_y and vacuum_z: + is_a_molecule = True + sys_type = "Molecule" + if not self.only_sys_type: + summary += "Molecule: \n" + all_molecules = molecules(list(range(len(atoms))), atoms) + com = np.average(atoms.positions, axis=0) + summary += ( + "COM: " + + str(com) + + ", min z: " + + str(np.min(atoms.positions[:, 2])) + + "\n" + ) + if vacuum_x and vacuum_y and (not vacuum_z): + is_a_wire = True + sys_type = "Wire" + cases = ["w"] + if not self.only_sys_type: + summary += "Wire along z contains: \n" + wireatoms = list(range(len(atoms))) + if vacuum_y and vacuum_z and (not vacuum_x): + is_a_wire = True + sys_type = "Wire" + cases = ["w"] + if not self.only_sys_type: + summary += "Wire along x contains: \n" + wireatoms = list(range(len(atoms))) + if vacuum_x and vacuum_z and (not vacuum_y): + is_a_wire = True + sys_type = "Wire" + cases = ["w"] + if not self.only_sys_type: + summary += "Wire along y contains: \n" + wireatoms = list(range(len(atoms))) + # wireatoms = slabatoms + is_a_slab = not (is_a_bulk or is_a_molecule or is_a_wire) + if self.only_sys_type: + if is_a_slab: + return {"system_type": "Slab"} + else: + return {"system_type": sys_type} + + elif is_a_slab: + cases = ["s"] + tipii, layersg, all_molecules = self.get_types(atoms, 0.1) + if vacuum_x: + slabtype = "YZ" + elif vacuum_y: + slabtype = "XZ" + else: + slabtype = "XY" + + sys_type = "Slab" + slabtype + mol_atoms = np.where(tipii == 0)[0].tolist() + metalatings = np.where(tipii == 6)[0].tolist() + mol_atoms += metalatings + + bottom_h = np.where(tipii == 3)[0].tolist() + + # Unclassified + unclassified = np.where(tipii == 5)[0].tolist() + + slabatoms = np.where(tipii == 1)[0].tolist() + adatoms = np.where(tipii == 2)[0].tolist() + + # Slab layers. + slab_layers = [[] for i in range(len(layersg))] + for ia in slabatoms: + idx = (np.abs(layersg - atoms.positions[ia, 2])).argmin() + slab_layers[idx].append(ia) + + summary += "Slab " + slabtype + " contains: \n" + summary += ( + "Cell: " + " ".join([str(i) for i in atoms.cell.diagonal().tolist()]) + "\n" + ) + if len(slabatoms) == 0: + slab_elements = set() + else: + slab_elements = set(atoms[slabatoms].get_chemical_symbols()) + + if len(bottom_h) > 0: + summary += "bottom H: " + mol_ids_range(bottom_h) + "\n" + if len(slabatoms) > 0: + summary += "slab atoms: " + mol_ids_range(slabatoms) + "\n" + for nlayer in range(len(slab_layers)): + summary += ( + "slab layer " + + str(nlayer + 1) + + ": " + + mol_ids_range(slab_layers[nlayer]) + + "\n" + ) + if len(adatoms) > 0: + cases.append("a") + summary += "adatoms: " + mol_ids_range(adatoms) + "\n" + if all_molecules: + cases.append("m") + summary += "#" + str(len(all_molecules)) + " molecules: " + for nmols in range(len(all_molecules)): + summary += str(nmols + 1) + ") " + mol_ids_range(all_molecules[nmols]) + + summary += " \n" + if len(metalatings) > 0: + summary += ( + "metal atoms inside molecules (already counted): " + + mol_ids_range(metalatings) + + "\n" + ) + if len(unclassified) > 0: + cases.append("u") + summary += "unclassified: " + mol_ids_range(unclassified) + + # Indexes from 0 if mol_ids_range is not called. + + return { + "total_charge": total_charge, + "system_type": sys_type, + "cell": " ".join([str(i) for i in itertools.chain(*atoms.cell.tolist())]), + "slab_layers": slab_layers, + "bottom_H": sorted(bottom_h), + "bulkatoms": sorted(bulkatoms), + "wireatoms": sorted(wireatoms), + "slabatoms": sorted(slabatoms), + "adatoms": sorted(adatoms), + "all_molecules": all_molecules, + "metalatings": sorted(metalatings), + "unclassified": sorted(unclassified), + "numatoms": len(atoms), + "all_elements": all_elements, + "slab_elements": slab_elements, + "spins_up": spins_up, + "spins_down": spins_down, + "other_tags": other_tags, + "sys_size": sys_size, + "cases": cases, + "summary": summary, + } From 03b7a6b9142a46420fef845bc4e6aade7458ec1a Mon Sep 17 00:00:00 2001 From: Carlo Antonio Pignedoli Date: Mon, 6 Jan 2025 10:54:15 +0000 Subject: [PATCH 2/3] bug fix --- surfaces_tools/widgets/constraints.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/surfaces_tools/widgets/constraints.py b/surfaces_tools/widgets/constraints.py index 1169c34..11bfb0d 100644 --- a/surfaces_tools/widgets/constraints.py +++ b/surfaces_tools/widgets/constraints.py @@ -1,9 +1,9 @@ import ipywidgets as ipw import traitlets as tr from aiida_nanotech_empa.workflows.cp2k import cp2k_utils +import aiidalab_widgets_base as awb from ase import Atoms -from .analyze_structure import mol_ids_range class OneColvar(ipw.HBox): @@ -141,7 +141,7 @@ def add_constraint(self, b=None): + self.details["slab_layers"][1] ) self.constraints.children[0].constraint_widget.value = ( - "fixed xyz " + mol_ids_range(to_fix) + "fixed xyz " + awb.utils.list_to_string_range(to_fix) ) def _observe_help(self, change): From e30150f4e0525d529c978034b15abf1a26280b97 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 6 Jan 2025 11:13:19 +0000 Subject: [PATCH 3/3] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- surfaces_tools/widgets/analyze_structure.py | 141 +++++++++++--------- surfaces_tools/widgets/constraints.py | 3 +- surfaces_tools/widgets/lowdimfinder.py | 22 ++- 3 files changed, 91 insertions(+), 75 deletions(-) diff --git a/surfaces_tools/widgets/analyze_structure.py b/surfaces_tools/widgets/analyze_structure.py index 63f4e33..782d5fd 100644 --- a/surfaces_tools/widgets/analyze_structure.py +++ b/surfaces_tools/widgets/analyze_structure.py @@ -1,14 +1,16 @@ -import copy import itertools -import aiidalab_widgets_base as awb +import aiidalab_widgets_base as awb import numpy as np +import traitlets as tr from ase import Atoms, neighborlist from scipy import sparse from scipy.signal import find_peaks from scipy.spatial import ConvexHull + from .lowdimfinder import LowDimFinder -import traitlets as tr + + def gaussian(x, sig): return ( 1.0 @@ -16,10 +18,13 @@ def gaussian(x, sig): * np.exp(-np.power(x, 2.0) / (2 * np.power(sig, 2.0))) ) + def boxfilter(x, thr): return np.asarray([1 if i < thr else 0 for i in x]) # Piero Gasparotto -def get_types( frame): # classify the atmos in: + + +def get_types(frame): # classify the atmos in: # 0=molecule # 1=slab atoms # 2=adatoms @@ -149,28 +154,28 @@ def get_types( frame): # classify the atmos in: atype[i] = 2 # assign the other types - #metalatingtypes = ("Au", "Ag", "Cu", "Ni", "Co", "Zn", "Mg") - #moltypes = ("H", "N", "B", "O", "C", "F", "S", "Br", "I", "Cl") - #possible_mol_atoms = [ + # metalatingtypes = ("Au", "Ag", "Cu", "Ni", "Co", "Zn", "Mg") + # moltypes = ("H", "N", "B", "O", "C", "F", "S", "Br", "I", "Cl") + # possible_mol_atoms = [ # i for i in range(nat) if atype[i] == 2 and lbls[i] in moltypes - #] - #possible_mol_atoms += [i for i in range(nat) if atype[i] == 5] + # ] + # possible_mol_atoms += [i for i in range(nat) if atype[i] == 5] # identify separate molecules # all_molecules=self.molecules(mol_atoms,atoms) - #all_molecules = [] + # all_molecules = [] return atype, layersg - + def transform_vector(pos1, pos2, v1): """ Transform vector v1 using the transformation that maps pos2 to pos1. - + Args: pos1 (np.ndarray): Target positions, shape (n, 3). pos2 (np.ndarray): Source positions, shape (n, 3). v1 (np.ndarray): Vector to transform, shape (3,). - + Returns: np.ndarray: Transformed vector, shape (3,). """ @@ -200,15 +205,15 @@ def transform_vector(pos1, pos2, v1): R = np.dot(Vt.T, U.T) # Translation vector - t = centroid1 - np.dot(R, centroid2) + centroid1 - np.dot(R, centroid2) # Transform the vector - transformed_v1 = -np.dot(R, v1) #+ t + transformed_v1 = -np.dot(R, v1) # + t transformed_v1 = transformed_v1 / np.linalg.norm(transformed_v1) - return transformed_v1 + class StructureAnalyzer(tr.HasTraits): structure = tr.Instance(Atoms, allow_none=True) details = tr.Dict() @@ -230,9 +235,7 @@ def analyze( return {} atoms = self.structure - - - + sys_size = np.ptp(atoms.positions, axis=0) no_cell = ( atoms.cell[0][0] < 0.1 or atoms.cell[1][1] < 0.1 or atoms.cell[2][2] < 0.1 @@ -242,35 +245,38 @@ def analyze( atoms.cell = sys_size + 10 atoms.set_pbc([True, True, True]) - - # LowDimFinder section + + # LowDimFinder section low_dim_finder = LowDimFinder( - # This is called aiida_structure but it's wrong, it's an ASE structure - #aiida_structure=conventional_asecell, - aiida_structure=atoms, - vacuum_space=40.0, - radii_offset=-0.75, # [-0.75 -0.7, -0.65, -0.6, -0.55] - bond_margin=0.0, - max_supercell=3, - min_supercell=3, - rotation=True, - full_periodicity=False, - radii_source="alvarez", - orthogonal_axis_2D=True, - ) + # This is called aiida_structure but it's wrong, it's an ASE structure + # aiida_structure=conventional_asecell, + aiida_structure=atoms, + vacuum_space=40.0, + radii_offset=-0.75, # [-0.75 -0.7, -0.65, -0.6, -0.55] + bond_margin=0.0, + max_supercell=3, + min_supercell=3, + rotation=True, + full_periodicity=False, + radii_source="alvarez", + orthogonal_axis_2D=True, + ) res = low_dim_finder.get_group_data() # End LowDimFinder section - is_a_bulk = np.amax(res['dimensionality']) == 3 - is_a_molecule = np.amax(res['dimensionality']) == 0 - is_a_wire = np.amax(res['dimensionality']) == 1 - is_a_slab = np.amax(res['dimensionality']) == 2 - max_dim_portion = res['dimensionality'].index(max(res['dimensionality'])) + is_a_bulk = np.amax(res["dimensionality"]) == 3 + is_a_molecule = np.amax(res["dimensionality"]) == 0 + is_a_wire = np.amax(res["dimensionality"]) == 1 + is_a_slab = np.amax(res["dimensionality"]) == 2 + max_dim_portion = res["dimensionality"].index(max(res["dimensionality"])) # orientation if not is_a_bulk and not is_a_molecule: - direction = transform_vector(atoms.positions[res['unit_cell_ids'][max_dim_portion]],np.array(res['positions'][max_dim_portion]),np.array([0,0,1])) + direction = transform_vector( + atoms.positions[res["unit_cell_ids"][max_dim_portion]], + np.array(res["positions"][max_dim_portion]), + np.array([0, 0, 1]), + ) dir_short = [f"{x:.2f}" for x in direction] # - total_charge = np.sum(atoms.get_atomic_numbers()) bottom_h = [] @@ -281,14 +287,16 @@ def analyze( unclassified = [] slabatoms = [] slab_layers = [] - all_molecules = [res['unit_cell_ids'][i] for i,j in enumerate(res['dimensionality']) if j == 0] + all_molecules = [ + res["unit_cell_ids"][i] + for i, j in enumerate(res["dimensionality"]) + if j == 0 + ] spins_up = [the_a.index for the_a in atoms if the_a.tag == 1] spins_down = [the_a.index for the_a in atoms if the_a.tag == 2] other_tags = [the_a.index for the_a in atoms if the_a.tag > 2] - - # Do not use a set in the following line list(set(atoms.get_chemical_symbols())) # need ALL atoms and elements for spin guess and for cost calculation all_elements = atoms.get_chemical_symbols() @@ -297,23 +305,27 @@ def analyze( if len(spins_up) > 0: summary += "spins_up: " + awb.utils.llist_to_string_range(spins_up) + "\n" if len(spins_down) > 0: - summary += "spins_down: " + awb.utils.llist_to_string_range(spins_down) + "\n" + summary += ( + "spins_down: " + awb.utils.llist_to_string_range(spins_down) + "\n" + ) if len(other_tags) > 0: - summary += "other_tags: " + awb.utils.llist_to_string_range(other_tags) + "\n" - - if is_a_bulk : - + summary += ( + "other_tags: " + awb.utils.llist_to_string_range(other_tags) + "\n" + ) + + if is_a_bulk: + sys_type = "Bulk" cases = ["b"] summary += "Bulk contains: \n" slabatoms = list(range(len(atoms))) bulkatoms = slabatoms - if is_a_molecule : - + if is_a_molecule: + sys_type = "Molecule" summary += "Molecule: \n" - #all_molecules = molecules(list(range(len(atoms))), atoms) + # all_molecules = molecules(list(range(len(atoms))), atoms) com = np.average(atoms.positions, axis=0) summary += ( "COM: " @@ -322,7 +334,7 @@ def analyze( + str(np.min(atoms.positions[:, 2])) + "\n" ) - if is_a_wire : + if is_a_wire: print("Wire") sys_type = "Wire" cases = ["w"] @@ -342,18 +354,23 @@ def analyze( slabtype = "XY" sys_type = "Slab" + slabtype - mol_atoms = [idatom for mol in all_molecules for idatom in mol] #np.where(tipii == 0)[0].tolist() + mol_atoms = [ + idatom for mol in all_molecules for idatom in mol + ] # np.where(tipii == 0)[0].tolist() # mol_atoms=extract_mol_indexes_from_slab(atoms) - #metalatings = np.where(tipii == 6)[0].tolist() - #mol_atoms += metalatings + # metalatings = np.where(tipii == 6)[0].tolist() + # mol_atoms += metalatings bottom_h = np.where(tipii == 3)[0].tolist() slabatoms = np.where(tipii == 1)[0].tolist() adatoms = np.where(tipii == 2)[0].tolist() - unclassified = [idatom for idatom in list(range(atoms.get_global_number_of_atoms())) if idatom not in set(bottom_h+slabatoms+mol_atoms)] # np.where(tipii == 5)[0].tolist() - - - #MOVE this up to get_types + unclassified = [ + idatom + for idatom in list(range(atoms.get_global_number_of_atoms())) + if idatom not in set(bottom_h + slabatoms + mol_atoms) + ] # np.where(tipii == 5)[0].tolist() + + # MOVE this up to get_types slab_layers = [[] for i in range(len(layersg))] for ia in slabatoms: idx = (np.abs(layersg - atoms.positions[ia, 2])).argmin() @@ -388,7 +405,9 @@ def analyze( cases.append("m") summary += "#" + str(len(all_molecules)) + " molecules: " for nmols, the_mol in enumerate(all_molecules): - summary += str(nmols + 1) + ") " + awb.utils.list_to_string_range(the_mol) + summary += ( + str(nmols + 1) + ") " + awb.utils.list_to_string_range(the_mol) + ) summary += " \n" if len(metalatings) > 0: diff --git a/surfaces_tools/widgets/constraints.py b/surfaces_tools/widgets/constraints.py index 11bfb0d..4d4a0fc 100644 --- a/surfaces_tools/widgets/constraints.py +++ b/surfaces_tools/widgets/constraints.py @@ -1,11 +1,10 @@ +import aiidalab_widgets_base as awb import ipywidgets as ipw import traitlets as tr from aiida_nanotech_empa.workflows.cp2k import cp2k_utils -import aiidalab_widgets_base as awb from ase import Atoms - class OneColvar(ipw.HBox): def __init__(self, cvtype="distance"): style = {"description_width": "initial"} diff --git a/surfaces_tools/widgets/lowdimfinder.py b/surfaces_tools/widgets/lowdimfinder.py index dd235da..409c6fc 100644 --- a/surfaces_tools/widgets/lowdimfinder.py +++ b/surfaces_tools/widgets/lowdimfinder.py @@ -1,4 +1,3 @@ -# -*- coding: utf-8 -*- """ Low dimensionality atom group finder, which analyses a bulk crystal and returns groups of atoms which are held together by weak (van der Waals) forces as @@ -14,12 +13,11 @@ ) __version__ = "0.3.0" -import numpy as np +from numbers import Number +import numpy as np from ase import Atoms from numpy import isscalar -from numbers import Number - ## Source: http://chemwiki.ucdavis.edu/Reference/Reference_Tables/Atomic_and_Molecular_Properties/A3%3A_Covalent_Radii # http://dx.doi.org/10.1039/b801115j, checked in paper @@ -558,8 +556,8 @@ def __init__(self, aiida_structure, **kwargs): self._low_dim_index = ( 0 # index of reduced structure that is currently calculated ) - self._found_unit_cell_atoms = set( - [] + self._found_unit_cell_atoms = ( + set() ) # set of all the input structure atoms that have been attributed to a group # independent variables, add those you want to have as output to get_group_data() @@ -651,7 +649,7 @@ def _define_groups(self): # pylint: disable=too-many-locals if atom not in visited: # fresh set of connected atoms. # Initially we add the analysed atom. - connected_atoms = set([atom]) + connected_atoms = {atom} new_neighs = set_list[atom] - connected_atoms all_neighbours = connected_atoms.copy() @@ -683,8 +681,8 @@ def _define_unit_cell_group(self): """ group = self.groups[self._group_number] - unit_cell_sites = set([]) - reduced_unit_cell_group = set([]) + unit_cell_sites = set() + reduced_unit_cell_group = set() for i in group: if i % self.n_unit not in unit_cell_sites: unit_cell_sites.add(i % self.n_unit) @@ -1350,7 +1348,7 @@ def get_group_data(self): "chemical_formula": self._chemical_formula, "positions": self._positions, "chemical_symbols": self._chemical_symbols, - "unit_cell_ids":self._unit_cell_groups, + "unit_cell_ids": self._unit_cell_groups, "cell": self._cell, "tags": self._tags, } @@ -1387,8 +1385,8 @@ def _define_3D_structures_with_layer_lattice(self): Define 3D structures with the same lattice as the layers for all 2D reduced structures found. """ - from aiida.plugins import ( # pylint:disable=import-error,import-outside-toplevel - DataFactory, # pylint:disable=import-error,import-outside-toplevel + from aiida.plugins import ( + DataFactory, # pylint:disable=import-error,import-outside-toplevel; pylint:disable=import-error,import-outside-toplevel ) S = DataFactory("structure")