Source code for acat.adsorption_sites

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from .settings import (adsorbate_elements, 
                       CustomSurface)
from .utilities import (expand_cell, 
                        get_mic, 
                        custom_warning,
                        is_list_or_tuple,
                        hash_composition, 
                        neighbor_shell_list, 
                        get_adj_matrix,
                        sorted_by_ref_atoms)
from .labels import (get_monometallic_cluster_labels, 
                     get_monometallic_slab_labels,
                     get_bimetallic_cluster_labels, 
                     get_bimetallic_slab_labels,
                     get_multimetallic_cluster_labels, 
                     get_multimetallic_slab_labels)
from ase.constraints import ExpCellFilter
from ase.geometry import (find_mic, 
                          wrap_positions, 
                          get_layers)
from ase.optimize import BFGS, FIRE
from ase import Atoms
from scipy.spatial.distance import pdist, squareform
from collections import defaultdict, Counter
from itertools import combinations, groupby, product
from copy import deepcopy
import networkx as nx
import numpy as np
import warnings
import random
import scipy
import math
import re
warnings.formatwarning = custom_warning


[docs]class ClusterAdsorptionSites(object): """Base class for identifying adsorption sites on a nanoparticle. Support common nanoparticle shapes including: Mackay icosahedron, (truncated) octahedron and (Marks) decahedron. The information of each site is stored in a dictionary with the following keys: **'site'**: the site type, support 'ontop', 'bridge', 'longbridge', 'shortbridge', 'fcc', 'hcp', '3fold', '4fold', '5fold', '6fold'. **'surface'**: the surface of the site, support 'vertex', 'edge', 'fcc100', 'fcc111'. **'position'**: the 3D Cartesian coordinate of the site saved as a numpy array. **'normal'**: the surface normal vector of the site saved as a numpy array. **'indices'**: the indices of the atoms that constitute the site. **'composition'**: the elemental composition of the site. Always in the lowest lexicographical order. **'subsurf_index'**: the index of the subsurface atom underneath an hcp or 4fold site. **'subsurf_element'**: the element of the subsurface atom underneath an hcp or 4fold site **'label'**: the numerical label assigned to the site if label_sites is set to True. Parameters ---------- atoms : ase.Atoms object The atoms object must be a non-periodic nanoparticle. Accept any ase.Atoms object. No need to be built-in. allow_6fold : bool, default False Whether to allow the adsorption on 6-fold subsurf sites underneath fcc hollow sites. composition_effect : bool, default False Whether to consider sites with different elemental compositions as different sites. It is recommended to set composition_effect=False for monometallics. ignore_sites : str or list of strs or None, default None The site type(s) to ignore. Default is to not ignore any site. label_sites : bool, default False Whether to assign a numerical label to each site. Labels for different sites are listed in acat.labels. Use the bimetallic labels if composition_effect=True, otherwise use the monometallic labels. surrogate_metal : str, default None The code is parameterized w.r.t. pure transition metals. The generalization of the code is achieved by mapping all input atoms to a surrogate transition metal that is supported by the asap3.EMT calculator (Ni, Cu, Pd, Ag, Pt or Au). Ideally this should be the metal that defines the lattice constant of the nanoparticle. Try changing the surrogate metal when the site identification is not satisfying. If no surrogate metal is provided, but the majority of the nanoparticle is a metal supported by asap3.EMT, the surrogate metal will be set to that metal automatically. tol : float, default 0.5 The tolerance of neighbor distance (in Angstrom). Might be helpful to adjust this if the site identification is not satisfying. When the nanoparticle is small (less than 300 atoms), Cu is normally the better choice, while Au should be good for larger nanoparticles. Example ------- The following example illustrates the most important use of a `ClusterAdsorptionSites` object - getting all adsorption sites: .. code-block:: python >>> from acat.adsorption_sites import ClusterAdsorptionSites >>> from ase.cluster import Octahedron >>> atoms = Octahedron('Ni', length=7, cutoff=2) >>> for atom in atoms: ... if atom.index % 2 == 0: ... atom.symbol = 'Pt' >>> atoms.center(vacuum=5.) >>> cas = ClusterAdsorptionSites(atoms, allow_6fold=False, ... composition_effect=True, ... label_sites=True, ... surrogate_metal='Ni') >>> site = cas.get_site() # Use cas.get_sites() to get all sites >>> print(site) Output: .. code-block:: python {'site': 'bridge', 'surface': 'edge', 'position': array([12.92, 14.68, 5. ]), 'normal': array([ 0.20864239, 0.48683225, -0.84821148]), 'indices': (115, 130), 'composition': 'NiPt', 'subsurf_index': None, 'subsurf_element': None, 'label': 10} """ def __init__(self, atoms, allow_6fold=False, composition_effect=False, ignore_sites=None, label_sites=False, surrogate_metal=None, tol=.5): assert True not in atoms.pbc, 'the cell must be non-periodic' warnings.filterwarnings('ignore', category=RuntimeWarning) warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning) atoms = atoms.copy() for dim in range(3): if np.linalg.norm(atoms.cell[dim]) == 0: atoms.cell[dim][dim] = np.ptp(atoms.positions[:, dim]) + 10. self.metal_ids = [a.index for a in atoms if a.symbol not in adsorbate_elements] atoms = Atoms(atoms.symbols[self.metal_ids], atoms.positions[self.metal_ids], cell=atoms.cell, pbc=atoms.pbc) self.atoms = atoms self.positions = atoms.positions self.symbols = atoms.symbols self.numbers = atoms.numbers self.allow_6fold = allow_6fold self.composition_effect = composition_effect if (ignore_sites is None) or is_list_or_tuple(ignore_sites): self.ignore_sites = ignore_sites else: self.ignore_sites = [ignore_sites] self.surrogate_metal = surrogate_metal self.tol = tol self.ref_atoms = self.mapping(atoms) self.cell = atoms.cell self.pbc = atoms.pbc self.metals = sorted(list(set(atoms.symbols))) self.label_sites = label_sites if self.composition_effect: if len(self.metals) <= 2: if len(self.metals) == 1: self.metals *= 2 self.label_dict = get_bimetallic_cluster_labels(self.metals) else: self.label_dict = get_multimetallic_cluster_labels(self.metals) else: self.label_dict = get_monometallic_cluster_labels() self.fullCNA = {} self.make_fullCNA() self.set_first_neighbor_distance_from_rdf() self.site_dict = self.get_site_dict() self.surf_ids, self.surf_sites = self.get_surface_sites() self.site_list = [] self.populate_site_list() self.postprocessing()
[docs] def populate_site_list(self): """Find all ontop, bridge and hollow sites (3-fold and 4-fold) given an input nanoparticle based on CNA analysis of the suface atoms and collect in a site list.""" from asap3 import FullNeighborList nblist = FullNeighborList(rCut=10., atoms=self.ref_atoms) ss = self.surf_sites ssall = set(ss['all']) fcna = self.get_fullCNA() sl = self.site_list usi = set() # used_site_indices normals_for_site = dict(list(zip(ssall, [[] for _ in ssall]))) for surface, sites in ss.items(): if surface == 'all': continue for s in sites: neighbors, _, dist2 = nblist.get_neighbors(s, self.r + 0.2) for n in neighbors: si = tuple(sorted([s, n])) # site_indices if n in ssall and si not in usi: # bridge sites pos = np.average(self.positions[[n, s]], 0) site_surf = self.get_surface_designation([s, n]) site = self.new_site() site.update({'site': 'bridge', 'surface': site_surf, 'position': np.round(pos, 8), 'indices': si}) if self.composition_effect: composition = ''.join(sorted(self.symbols[list(si)])) site.update({'composition': composition}) sl.append(site) usi.add(si) # Find normal for n, m in combinations(neighbors, 2): si = tuple(sorted([s, n, m])) if n in ssall and m in ssall and si not in usi: angle = self.get_angle([s, n, m]) if self.is_eq(angle, np.pi/3.): # 3-fold (fcc or hcp) site normal = self.get_surface_normal([s, n, m]) for i in [s, n, m]: normals_for_site[i].append(normal) pos = np.average(self.positions[[n, m, s]], 0) new_pos = pos - normal * self.r * (2./3)**(.5) isubs = [a.index for a in self.ref_atoms if np.linalg.norm(a.position - new_pos) < 0.5] if not isubs: this_site = 'fcc' else: isub = isubs[0] this_site = 'hcp' site_surf = 'fcc111' site = self.new_site() site.update({'site': this_site, 'surface': site_surf, 'position': np.round(pos, 8), 'normal': np.round(normal, 8), 'indices': si}) if self.composition_effect: metals = self.metals if len(metals) == 1: composition = 3*metals[0] elif len(metals) == 2: ma, mb = metals[0], metals[1] symbols = self.symbols[list(si)] nma = np.count_nonzero(symbols==ma) if nma == 0: composition = 3*mb elif nma == 1: composition = ma + 2*mb elif nma == 2: composition = 2*ma + mb elif nma == 3: composition = 3*ma else: nodes = self.symbols[list(si)] composition = ''.join(hash_composition(nodes)) site.update({'composition': composition}) if this_site == 'hcp': site.update({'subsurf_index': isub}) if self.composition_effect: site.update({'subsurf_element': self.symbols[isub]}) sl.append(site) usi.add(si) elif self.is_eq(angle, np.pi/2.): # 4-fold hollow site site_surf = 'fcc100' l2 = self.r * math.sqrt(2) + 0.2 nebs2, _, _ = nblist.get_neighbors(s, l2) nebs2 = [k for k in nebs2 if k not in neighbors] for k in nebs2: si = tuple(sorted([s, n, m, k])) if k in ssall and si not in usi: d1 = self.ref_atoms.get_distance(n, k) if self.is_eq(d1, self.r, 0.2): d2 = self.ref_atoms.get_distance(m, k) if self.is_eq(d2, self.r, 0.2): # 4-fold hollow site found normal = self.get_surface_normal([s, n, m]) # Save the normals now and add them to the site later for i in [s, n, m, k]: normals_for_site[i].append(normal) ps = self.positions[[n, m, s, k]] pos = np.average(ps, 0) site = self.new_site() site.update({'site': '4fold', 'surface': site_surf, 'position': np.round(pos, 8), 'normal': np.round(normal, 8), 'indices': si}) if self.composition_effect: metals = self.metals if len(metals) == 1: composition = 4*metals[0] elif len(metals) == 2: ma, mb = metals[0], metals[1] symbols = self.symbols[list(si)] nma = np.count_nonzero(symbols==ma) if nma == 0: composition = 4*mb elif nma == 1: composition = ma + 3*mb elif nma == 2: opp = max(list(si[1:]), key=lambda x: np.linalg.norm(self.positions[x] - self.positions[si[0]])) if self.symbols[opp] == self.symbols[si[0]]: composition = ma + mb + ma + mb else: composition = 2*ma + 2*mb elif nma == 3: composition = 3*ma + mb elif nma == 4: composition = 4*ma else: opp = max(list(si[1:]), key=lambda x: np.linalg.norm(self.positions[x] - self.positions[si[0]])) nni = [i for i in si[1:] if i != opp] nodes = self.symbols[[si[0], nni[0], opp, nni[1]]] composition = ''.join(hash_composition(nodes)) site.update({'composition': composition}) new_pos = pos - normal * self.r * (2./3)**(.5) isub = min(range(len(self.atoms)), key=lambda x: np.linalg.norm(self.positions[x] - new_pos)) site.update({'subsurf_index': isub}) if self.composition_effect: site.update({'subsurf_element': self.symbols[isub]}) sl.append(site) usi.add(si) # ontop sites site = self.new_site() site.update({'site': 'ontop', 'surface': surface, 'position': np.round(self.positions[s], 8), 'indices': (s,)}) if self.composition_effect: site.update({'composition': self.symbols[s]}) sl.append(site) usi.add((s)) # Add 6-fold sites if allowed if self.allow_6fold: dh = 2. * self.r / 5. subsurf_ids = self.get_subsurface() for t in sl: # Add normals to ontop sites if t['site'] == 'ontop': n = np.average(normals_for_site[t['indices'][0]], 0) t['normal'] = np.round(n / np.linalg.norm(n), 8) # Add normals to bridge sites elif t['site'] == 'bridge': normals = [] for i in t['indices']: normals.extend(normals_for_site[i]) n = np.average(normals, 0) t['normal'] = np.round(n / np.linalg.norm(n), 8) # Add subsurf sites if self.allow_6fold: if t['site'] == 'fcc': site = t.copy() subpos = t['position'] - t['normal'] * dh def get_squared_distance(x): return np.sum((self.positions[x] - subpos)**2) subsi = sorted(sorted(subsurf_ids, key=get_squared_distance)[:3]) si = site['indices'] site.update({'site': '6fold', 'position': np.round(subpos, 8), 'indices': tuple(sorted(si+tuple(subsi)))}) if self.composition_effect: metals = self.metals comp = site['composition'] if len(metals) == 1: composition = ''.join([comp, 3*metals[0]]) elif len(metals) == 2: ma, mb = metals[0], metals[1] subsyms = self.symbols[subsi] nma = np.count_nonzero(subsyms==ma) if nma == 0: composition = ''.join([comp, 3*mb]) elif nma == 1: ia = subsi[subsyms.index(ma)] subpos = self.positions[ia] if self.symbols[max(si, key=get_squared_distance)] == ma: composition = ''.join([comp, 2*mb + ma]) else: composition = ''.join([comp, ma + 2*mb]) elif nma == 2: ib = subsi[subsyms.index(mb)] subpos = self.positions[ib] if self.symbols[max(si, key=get_squared_distance)] == mb: composition = ''.join([comp, 2*ma + mb]) else: composition = ''.join([comp, mb + 2*ma]) elif nma == 3: composition = ''.join([comp, 3*ma]) else: endsym = re.findall('[A-Z][^A-Z]*', comp)[-1] subpos = [self.positions[i] for i in si if self.symbols[i] == endsym][0] nodes = self.symbols[sorted(subsi, key=get_squared_distance)] composition = comp + ''.join(hash_composition(nodes)) site.update({'composition': composition}) sl.append(site) usi.add(si)
[docs] def postprocessing(self): """Postprocessing the site list.""" if self.ignore_sites is not None: self.site_list = [s for s in self.site_list if s['site'] not in self.ignore_sites] if self.label_sites: self.get_labels() for i, st in enumerate(self.site_list): st['indices'] = tuple(self.metal_ids[j] for j in st['indices']) self.site_list.sort(key=lambda x: x['indices'])
[docs] def get_site(self, indices=None): """Get information of an adsorption site. Parameters ---------- indices : list or tuple, default None The indices of the atoms that contribute to the site. Return a random site if the indices are not provided. """ if indices is None: st = random.choice(self.site_list) else: indices = indices if is_list_or_tuple(indices) else [indices] indices = tuple(sorted(indices)) st = next((s for s in self.site_list if s['indices'] == indices), None) return st
[docs] def get_sites(self, site=None, surface=None, composition=None, subsurf_element=None, return_site_indices=False): """Get information of all sites. Parameters ---------- site : str, default None Only return sites that belongs to this site type. surface : str, default None Only return sites that are on this surface. composition : str, default None Only return sites that have this composition. subsurf_element : str, default None Only return sites that have this subsurface element. return_site_indices: bool, default False Whether to return the indices of each unique site (in the site list). """ if return_site_indices: all_sites = deepcopy(self.site_list) for i, st in enumerate(all_sites): st['site_index'] = i else: all_sites = self.site_list if site is not None: all_sites = [s for s in all_sites if s['site'] == site] if surface is not None: all_sites = [s for s in all_sites if s['surface'] == surface] if composition is not None: if '-' in composition or len(list(Formula(composition))) == 6: scomp = composition else: comp = re.findall('[A-Z][^A-Z]*', composition) if len(comp) != 4: scomp = ''.join(sorted(comp)) else: if comp[0] != comp[2]: scomp = ''.join(sorted(comp)) else: if comp[0] > comp[1]: scomp = comp[1]+comp[0]+comp[3]+comp[2] else: scomp = ''.join(comp) all_sites = [s for s in all_sites if s['composition'] == scomp] if subsurf_element is not None: all_sites = [s for s in all_sites if s['subsurf_element'] == subsurf_element] if return_site_indices: all_stids = [] for st in all_sites: all_stids.append(st['site_index']) del st['site_index'] return all_stids else: return all_sites
[docs] def get_unique_sites(self, unique_composition=False, unique_subsurf=False, return_signatures=False, return_site_indices=False, about=None, site_list=None): """Get all symmetry-inequivalent adsorption sites (one site for each type). Parameters ---------- unique_composition : bool, default False Take site composition into consideration when checking uniqueness. unique_subsurf : bool, default False Take subsurface element into consideration when checking uniqueness. Could be important for surfaces like fcc100. return_signatures : bool, default False Whether to return the unique signatures of the sites instead. return_site_indices: bool, default False Whether to return the indices of each unique site (in the site list). about: numpy.array, default None If specified, returns unique sites closest to this reference position. """ if site_list is None: sl = self.site_list else: sl = site_list key_list = ['site', 'surface'] if unique_composition: if not self.composition_effect: raise ValueError('the site list does not include ' + 'information of composition') key_list.append('composition') if unique_subsurf: key_list.append('subsurf_element') else: if unique_subsurf: raise ValueError('to include the subsurface element, ' + 'unique_composition also need to be set to True') if return_signatures: sklist = sorted([[s[k] for k in key_list] for s in sl]) return sorted(list(sklist for sklist, _ in groupby(sklist))) else: seen_tuple = [] unq_sites = [] if about is not None: sl = sorted(sl, key=lambda x: np.linalg.norm(x['position'] - about)) for i, s in enumerate(sl): sig = tuple(s[k] for k in key_list) if sig not in seen_tuple: seen_tuple.append(sig) if return_site_indices: s = i unq_sites.append(s) return unq_sites
def get_labels(self): # Assign labels for st in self.site_list: if self.composition_effect: signature = [st['site'], st['surface'], st['composition']] else: signature = [st['site'], st['surface']] stlab = self.label_dict['|'.join(signature)] st['label'] = stlab def new_site(self): return {'site': None, 'surface': None, 'position': None, 'normal': None, 'indices': None, 'composition': None, 'subsurf_index': None, 'subsurf_element': None, 'label': None}
[docs] def mapping(self, atoms): """Map the nanoparticle into a surrogate nanoparticle for code versatility.""" from asap3 import EMT ref_atoms = atoms.copy() pm = self.surrogate_metal if pm is None: common_metal = Counter(self.atoms.symbols).most_common(1)[0][0] if common_metal in ['Ni', 'Cu', 'Pd', 'Ag', 'Pt', 'Au']: ref_symbol = common_metal else: ref_symbol = 'Au' if len(atoms) > 300 else 'Cu' else: ref_symbol = pm for a in ref_atoms: a.symbol = ref_symbol ref_atoms.calc = EMT() opt = FIRE(ref_atoms, logfile=None) opt.run(fmax=0.1) ref_atoms.calc = None return ref_atoms
def get_two_vectors(self, indices): p1 = self.positions[indices[1]] p2 = self.positions[indices[2]] vec1 = p1 - self.positions[indices[0]] vec2 = p2 - self.positions[indices[0]] return vec1, vec2 def is_eq(self, v1, v2, eps=0.1): if abs(v1 - v2) < eps: return True else: return False
[docs] def get_surface_normal(self, indices): """Get the surface normal vector of the plane from the indices of 3 atoms that forms to that plane. Parameters ---------- indices : list of tuple The indices of the atoms that forms the plane. """ vec1, vec2 = self.get_two_vectors(indices) n = np.cross(vec1, vec2) l = math.sqrt(n @ n.conj()) new_pos = self.positions[indices[0]] + self.r * n / l # Add support for having adsorbates on the particles already # by putting in elements to check for in the function below j = 2 * int(self.no_atom_too_close_to_pos(new_pos, (5./6)*self.r)) - 1 return j * n / l
def get_angle(self, indices): vec1, vec2 = self.get_two_vectors(indices) p = (vec1 @ vec2)/(np.linalg.norm(vec1) * np.linalg.norm(vec1)) return np.arccos(np.clip(p, -1,1))
[docs] def no_atom_too_close_to_pos(self, pos, mindist): """Returns True if no atoms are closer than mindist to pos, otherwise False. Parameters ---------- pos : numpy.array The position to be checked. mindist : float The minimum distance (in Angstrom) that is not considered as too close. """ dists = [np.linalg.norm(atom.position - pos) > mindist for atom in self.ref_atoms] return all(dists)
[docs] def get_surface_sites(self): """Returns the indices of the surface atoms and a dictionary with all the surface designations.""" surf_sites = {'all': [], 'fcc111': [], 'fcc100': [], 'edge': [], 'vertex': [],} fcna = self.get_fullCNA() site_dict = self.site_dict surf_ids = [] for i in range(len(self.atoms)): # if i in [284, 310]: # print(fcna[i]) if sum(fcna[i].values()) < 12: surf_sites['all'].append(i) if str(fcna[i]) not in site_dict: # The structure is distorted from the original, giving # a larger cutoff should overcome this problem r = self.r + self.tol fcna = self.get_fullCNA(rCut=r) if str(fcna[i]) not in site_dict: # If this does not solve the problem we probably have a # reconstruction of the surface and will leave this # atom unused this time continue surf_ids.append(i) surf_sites[site_dict[str(fcna[i])]].append(i) return surf_ids, surf_sites
[docs] def get_subsurface(self): """Returns the indices of the subsurface atoms.""" from asap3.analysis import FullCNA notsurf = [a.index for a in self.atoms if a.index not in self.surf_ids] subfcna = FullCNA(self.ref_atoms[notsurf]).get_normal_cna() return [idx for i, idx in enumerate(notsurf) if sum(subfcna[i].values()) < 12]
def make_fullCNA(self, rCut=None): if rCut not in self.fullCNA: from asap3.analysis import FullCNA self.fullCNA[rCut] = FullCNA(self.ref_atoms, rCut=rCut).get_normal_cna()
[docs] def get_fullCNA(self, rCut=None): """Get the CNA signatures of all atoms by asap3 full CNA analysis. Parameters ---------- rCut : float, default None The cutoff radius in Angstrom. If not specified, the asap3 CNA analysis will use a reasonable cutoff based on the crystalline lattice constant of the material. """ if rCut not in self.fullCNA: self.make_fullCNA(rCut=rCut) return self.fullCNA[rCut]
[docs] def get_connectivity(self): """Get the adjacency matrix.""" nbslist = neighbor_shell_list(self.ref_atoms, 0.3, neighbor_number=1) return get_adj_matrix(nbslist)
def get_site_dict(self): icosa_dict = { # Triangle sites on outermost shell -- Icosa, Cubocta, Deca, Tocta str({(3, 1, 1): 6, (4, 2, 1): 3}): 'fcc111', 'fcc111': [str({(3, 1, 1): 6, (4, 2, 1): 3})], # Edge sites on outermost shell -- Icosa str({(3, 1, 1): 4, (3, 2, 2): 2, (4, 2, 2): 2}): 'edge', 'edge': [str({(3, 1, 1): 4, (3, 2, 2): 2, (4, 2, 2): 2})], # Vertice sites on outermost shell -- Icosa, Deca str({(3, 2, 2): 5, (5, 5, 5): 1}): 'vertex', 'vertex': [str({(3, 2, 2): 5, (5, 5, 5): 1})], } cubocta_dict = { # Edge sites on outermost shell -- Cubocta, Tocta str({(2, 1, 1): 3, (3, 1, 1): 2, (4, 2, 1): 2}): 'edge', 'edge': [str({(2, 1, 1): 3, (3, 1, 1): 2, (4, 2, 1): 2})], # Square sites on outermost shell -- Cubocta, Deca, Tocta str({(2, 1, 1): 4, (4, 2, 1): 4}): 'fcc100', 'fcc100': [str({(2, 1, 1): 4, (4, 2, 1): 4})], # Vertice sites on outermost shell -- Cubocta str({(2, 1, 1): 4, (4, 2, 1): 1}): 'vertex', 'vertex': [str({(2, 1, 1): 4, (4, 2, 1): 1})], # Triangle sites on outermost shell -- Icosa, Cubocta, Deca, Tocta str({(3, 1, 1): 6, (4, 2, 1): 3}): 'fcc111', 'fcc111': [str({(3, 1, 1): 6, (4, 2, 1): 3})], } mdeca_dict = { # Edge sites (111)-(111) on outermost shell -- Deca str({(3, 1, 1): 4, (3, 2, 2): 2, (4, 2, 2): 2}): 'edge', 'edge': [str({(3, 1, 1): 4, (3, 2, 2): 2, (4, 2, 2): 2}), str({(2, 1, 1): 3, (3, 1, 1): 2, (4, 2, 1): 2}), str({(2, 0, 0): 2, (3, 1, 1): 4, (4, 2, 1): 1})], # Edge sites (111)-(100) on outermost shell -- Deca str({(2, 1, 1): 3, (3, 1, 1): 2, (4, 2, 1): 2}): 'edge', # Edge sites (111)-(111)notch on outermost shell -- Deca str({(2, 0, 0): 2, (3, 1, 1): 4, (4, 2, 1): 1}): 'edge', # Square sites on outermost shell -- Cubocta, Deca, Tocta str({(2, 1, 1): 4, (4, 2, 1): 4}): 'fcc100', 'fcc100': [str({(2, 1, 1): 4, (4, 2, 1): 4})], # Vertice sites on outermost shell -- Icosa, Deca str({(3, 2, 2): 5, (5, 5, 5): 1}): 'vertex', 'vertex': [str({(3, 2, 2): 5, (5, 5, 5): 1}), str({(2, 0, 0): 1, (2, 1, 1): 2, (3, 1, 1): 2, (4, 2, 1): 1}), str({(2, 0, 0): 2, (3, 0, 0): 1, (3, 1, 1): 2, (3, 2, 2): 1, (4, 2, 2): 1})], # Vertice sites A on outermost shell -- Mdeca str({(2, 0, 0): 1, (2, 1, 1): 2, (3, 1, 1): 2, (4, 2, 1): 1}): 'vertex', # Vertice sites B on outermost shell -- Mdeca str({(2, 0, 0): 2, (3, 0, 0): 1, (3, 1, 1): 2, (3, 2, 2): 1, (4, 2, 2): 1}): 'vertex', # Triangle (pentagon) sites on outermost shell -- Icosa, Cubocta, Deca, Tocta str({(3, 1, 1): 6, (4, 2, 1): 3}): 'fcc111', 'fcc111': [str({(3, 1, 1): 6, (4, 2, 1): 3}), str({(3, 0, 0): 2, (3, 1, 1): 4, (4, 2, 1): 2, (4, 2, 2): 2})], # Triangle (pentagon) notch sites on outermost shell -- Deca str({(3, 0, 0): 2, (3, 1, 1): 4, (4, 2, 1): 2, (4, 2, 2): 2}): 'fcc111', } tocta_dict = { # Edge sites on outermost shell -- Cubocta, Tocta str({(2, 1, 1): 3, (3, 1, 1): 2, (4, 2, 1): 2}): 'edge', 'edge': [str({(2, 1, 1): 3, (3, 1, 1): 2, (4, 2, 1): 2}), str({(2, 0, 0): 2, (3, 1, 1): 4, (4, 2, 1): 1})], # Edge sites (111)-(111) on outermost shell -- Octa str({(2, 0, 0): 2, (3, 1, 1): 4, (4, 2, 1): 1}): 'edge', # Square sites on outermost shell -- Cubocta, Deca, Tocta str({(2, 1, 1): 4, (4, 2, 1): 4}): 'fcc100', 'fcc100': [str({(2, 1, 1): 4, (4, 2, 1): 4})], # Vertice sites on outermost shell -- Tocta str({(2, 0, 0): 1, (2, 1, 1): 2, (3, 1, 1): 2, (4, 2, 1): 1}): 'vertex', 'vertex': [str({(2, 0, 0): 1, (2, 1, 1): 2, (3, 1, 1): 2, (4, 2, 1): 1})], # Triangle (pentagon) sites on outermost shell -- Icosa, Cubocta, Deca, Octa str({(3, 1, 1): 6, (4, 2, 1): 3}): 'fcc111', 'fcc111': [str({(3, 1, 1): 6, (4, 2, 1): 3})], } fcna = self.get_fullCNA() icosa_weight = cubocta_weight = mdeca_weight = tocta_weight = 0 for s in fcna: if str(s) in icosa_dict: icosa_weight += 1 if str(s) in cubocta_dict: cubocta_weight += 1 if str(s) in mdeca_dict: mdeca_weight += 1 if str(s) in tocta_dict: tocta_weight += 1 full_weights = [icosa_weight, cubocta_weight, mdeca_weight, tocta_weight] if icosa_weight == max(full_weights): return icosa_dict elif cubocta_weight == max(full_weights): return cubocta_dict elif tocta_weight == max(full_weights): return tocta_dict else: return mdeca_dict def set_first_neighbor_distance_from_rdf(self, rMax=10, nBins=200): from asap3.analysis import rdf atoms = self.ref_atoms.copy() for j, L in enumerate(list(atoms.cell.diagonal())): if L <= 10: atoms.cell[j][j] = 12 _rdf = rdf.RadialDistributionFunction(atoms, rMax, nBins).get_rdf() x = (np.arange(nBins) + 0.5) * rMax / nBins _rdf *= x**2 diff_rdf = np.gradient(_rdf) i = 0 while diff_rdf[i] >= 0: i += 1 self.r = x[i] def get_surface_designation(self, indices): fcna = self.get_fullCNA() sd = self.site_dict if len(indices) == 1: s = indices[0] return sd[str(fcna[s])] elif len(indices) == 2: if str(fcna[indices[0]]) not in sd or str(fcna[indices[1]]) not in sd: # for s in [indices[0], indices[1]]: # scna = str(fcna[s]) # if scna not in sd: # print('CNA {} is not supported.'.format(scna)) return 'unknown' s0 = sd[str(fcna[indices[0]])] s1 = sd[str(fcna[indices[1]])] ss01 = tuple(sorted([s0, s1])) if ss01 in [('edge', 'edge'), ('edge', 'vertex'), ('vertex', 'vertex')]: return 'edge' elif ss01 in [('fcc111', 'vertex'), ('edge', 'fcc111'), ('fcc111', 'fcc111')]: return 'fcc111' elif ss01 in [('fcc100', 'vertex'), ('edge', 'fcc100'), ('fcc100', 'fcc100')]: return 'fcc100' elif len(indices) == 3: return 'fcc111'
[docs] def get_graph(self, return_adj_matrix=False): """Get the graph representation of the slab. Parameters ---------- return_adj_matrix : bool, default False Whether to return adjacency matrix instead of the networkx.Graph object. """ cm = self.get_connectivity() if return_adj_matrix: return cm G = nx.Graph() symbols = self.symbols G.add_nodes_from([(i, {'symbol': symbols[i]}) for i in range(len(symbols))]) rows, cols = np.where(cm == 1) edges = zip(rows.tolist(), cols.tolist()) G.add_edges_from(edges) return G
[docs] def get_neighbor_site_list(self, dx=0.1, neighbor_number=1, radius=None, span=True, unique=False, unique_composition=False, unique_subsurf=False): """Returns the site_list index of all neighbor shell sites for each site. Parameters ---------- dx : float, default 0.1 Buffer to calculate nearest neighbor site pairs. neighbor_number : int, default 1 Neighbor shell number. radius : float, default None The radius that defines site a is a neighbor of site b. This serves as an alternative to neighbor_number. span : bool, default True Whether to include all neighbors sites spanned within the shell. unique : bool, default False Whether to discard duplicate types of neighbor sites. """ sl = self.site_list poss = np.asarray([s['position'] for s in sl]) statoms = Atoms('X{}'.format(len(sl)), positions=poss, cell=self.cell, pbc=self.pbc) if radius is not None: nbslist = neighbor_shell_list(statoms, dx, neighbor_number=1, mic=False, radius=radius, span=span) if unique: for i in range(len(nbslist)): nbsids = nbslist[i] tmpids = self.get_unique_sites(unique_composition, unique_subsurf, return_site_indices=True, site_list=[sl[ni] for ni in nbsids]) nbsids = [nbsids[j] for j in tmpids] nbslist[i] = nbsids else: D = statoms.get_all_distances(mic=False) nbslist = {} for i, row in enumerate(D): argsort = np.argsort(row) row = row[argsort] res = np.split(argsort, np.where(np.abs(np.diff(row)) > dx)[0] + 1) if span: nbsids = sorted([n for nn in range(1, neighbor_number + 1) for n in list(res[nn])]) else: nbsids = sorted(res[neighbor_number]) if unique: tmpids = self.get_unique_sites(unique_composition, unique_subsurf, return_site_indices=True, site_list=[sl[ni] for ni in nbsids]) nbsids = [nbsids[j] for j in tmpids] nbslist[i] = nbsids return nbslist
[docs] def update(self, atoms, full_update=False): """Update the position of each adsorption site given an updated atoms object. Please only use this when the indexing of the atoms object is preserved. Useful for updating adsorption sites e.g. after geometry optimization. Parameters ---------- atoms : ase.Atoms object The updated atoms object. full_update : bool, default False Whether to update the normal vectors and composition as well. Useful when the composition of the nanoalloy is not fixed. """ sl = self.site_list if full_update: ncas = ClusterAdsorptionSites(atoms, allow_6fold=self.allow_6fold, composition_effect=self.composition_effect, ignore_sites=self.ignore_sites, label_sites=self.label_sites, surrogate_metal=self.surrogate_metal, tol=self.tol) sl = ncas.site_list else: new_cluster = atoms[self.metal_ids] nsl = deepcopy(sl) for st in nsl: si = list(st['indices']) idd = {idx: i for i, idx in enumerate(self.metal_ids)} tmpsi = [idd[k] for k in si] newpos = np.average(new_cluster.positions[tmpsi], 0) st['position'] = np.round(newpos, 8) sl = nsl
[docs]def group_sites_by_facet(atoms, sites, all_sites=None): """A function that uses networkx to group one set of sites by geometrical facets of the nanoparticle. Different geometrical facets can have the same surface type. The function returns a list of lists, each contains sites on a same geometrical facet. Parameters ---------- atoms : ase.Atoms object The atoms object must be a non-periodic nanoparticle. Accept any ase.Atoms object. No need to be built-in. sites : list of dicts The adsorption sites to be grouped by geometrical facet. all_sites : list of dicts, default None The list of all sites. Provide this to make the grouping much faster. Useful when the function is called many times. Example ------- The following example shows how to group all fcc sites of an icosahedral nanoparticle by its 20 geometrical facets: .. code-block:: python >>> from acat.adsorption_sites import ClusterAdsorptionSites >>> from acat.adsorption_sites import group_sites_by_facet >>> from ase.cluster import Icosahedron >>> atoms = Icosahedron('Pt', noshells=5) >>> atoms.center(vacuum=5.) >>> cas = ClusterAdsorptionSites(atoms) >>> all_sites = cas.get_sites() >>> fcc_sites = [s for s in all_sites if s['site'] == 'fcc'] >>> groups = group_sites_by_facet(atoms, fcc_sites, all_sites) >>> print(len(groups)) Output: .. code-block:: python 20 """ # Find all indices of vertex and edge sites if not all_sites: cas = ClusterAdsorptionSites(atoms) all_sites = cas.site_list ve_indices = set(s['indices'][0] for s in all_sites if s['site'] == 'ontop' and s['surface'] in ['vertex', 'edge']) G = nx.Graph() for site in sites: indices = site['indices'] reduced_indices = tuple(i for i in indices if i not in ve_indices) if len(reduced_indices) > 0: site['reduced_indices'] = reduced_indices nx.add_path(G, reduced_indices) components = list(nx.connected_components(G)) groups = [] for component in components: group = [] for path in [s['reduced_indices'] for s in sites if 'reduced_indices' in s]: if component.issuperset(path): group.append(path) groups.append(group) grouped_sites = defaultdict(list) for site in sites: if 'reduced_indices' in site: for group in groups: if site['reduced_indices'] in group: grouped_sites[groups.index(group)] += [site] del site['reduced_indices'] return grouped_sites
[docs]class SlabAdsorptionSites(object): """Base class for identifying adsorption sites on a surface slab. Support 20 common surfaces: fcc100, fcc111, fcc110, fcc211, fcc221, fcc311, fcc322, fcc331, fcc332, bcc100, bcc111, bcc110, bcc210, bcc211, bcc310, hcp0001, hcp10m10t, hcp10m10h, hcp10m11, hcp10m12. A CustomSurface object is also supported, which can be used for any user-defined surface, including metal oxide surfaces. The information of each site is stored in a dictionary with the following keys: **'site'**: the site type, support 'ontop', 'bridge', 'longbridge', 'shortbridge', 'fcc', 'hcp', '3fold', '4fold', '5fold', '6fold'. **'surface'**: the surface type (crystal structure + Miller indices) or an acat.settings.CustomSurface object. **'morphology'**: the local surface morphology of the site. Support 'step', 'terrace', 'corner', 'sc-tc(-x)', 'tc-cc(-x)', 'sc-cc(-x)'. 'sc', 'tc' and 'cc' represents step chain, terrace chain and corner, respectively. 'x' is the Bravais lattice that connects the 2 chains, e.g. 'h' = hexagonal, 't' = 'tetragonal', 'o' = 'orthorhombic'. **'position'**: the 3D Cartesian coordinate of the site saved as a numpy array. **'normal'**: the surface normal vector of the site saved as a numpy array. **'indices'**: the indices of the atoms that constitute the site. **'composition'**: the elemental composition of the site. Always in the lowest lexicographical order. **'subsurf_index'**: the index of the subsurface atom underneath an hcp or 4fold site. **'subsurf_element'**: the element of the subsurface atom underneath an hcp or 4fold site **'label'**: the numerical label assigned to the site if label_sites is set to True. Parameters ---------- atoms : ase.Atoms object The atoms object must be a periodic surface slab with at least 3 layers (e.g. all surface atoms make up one layer). Accept any ase.Atoms object. No need to be built-in. surface : str or acat.settings.CustomSurface object The surface type (crystal structure + Miller indices). A user-defined customized surface object can also be used, but note that the identified site types will only include 'ontop', 'bridge', '3fold' and '4fold'. allow_6fold : bool, default False Whether to allow the adsorption on 6-fold subsurf sites underneath fcc hollow sites. composition_effect : bool, default False Whether to consider sites with different elemental compositions as different sites. It is recommended to set composition_effect=False for monometallics. both_sides : bool, default False Whether to consider sites on both top and bottom sides of the slab. ignore_sites : str or list of strs or None, default None The site type(s) to ignore. Default is to not ignore any site. label_sites : bool, default False Whether to assign a numerical label to each site. Labels for different sites are listed in acat.labels. Use the bimetallic labels if composition_effect=True, otherwise use the monometallic labels. surrogate_metal : str, default None The code is parameterized w.r.t. pure transition metals. The generalization of the code is achieved by mapping all input atoms to a surrogate transition metal that is supported by the asap3.EMT calculator (Ni, Cu, Pd, Ag, Pt or Au). Ideally this should be the metal that defines the lattice constant of the slab. Try changing the surrogate metal when the site identification is not satisfying. When the cell is small, Cu is normally the better choice, while the Pt and Au should be good for larger cells. If no surrogate metal is provided, but the majority of the slab is a metal supported by asap3.EMT, the surrogate metal will be set to that metal automatically. optimize_surrogate_cell : bool, default False Whether to also optimize the cell during the optimization of the surrogate slab. Useful for generalizing the code for alloy slabs of various compositions and lattice constants. Recommand to set to True if the surrogate metal is very different from the composition of the input slab. tol : float, default 0.5 The tolerance of neighbor distance (in Angstrom). Might be helpful to adjust this if the site identification is not satisfying. The default 0.5 is usually good enough. Examples -------- The following example illustrates the most important use of a `SlabAdsorptionSites` object - getting all adsorption sites. Here we use the example of an fcc211 surface (considering both top and bottom terminations): .. code-block:: python >>> from acat.adsorption_sites import SlabAdsorptionSites >>> from ase.build import fcc211 >>> atoms = fcc211('Cu', (3, 3, 4), vacuum=5.) >>> for atom in atoms: ... if atom.index % 2 == 0: ... atom.symbol = 'Au' >>> atoms.center(axis=2) >>> sas = SlabAdsorptionSites(atoms, surface='fcc211', ... allow_6fold=False, ... composition_effect=True, ... both_sides=True, ... label_sites=True, ... surrogate_metal='Cu') >>> site = sas.get_site() # Use sas.get_sites() to get all sites >>> print(site) Output: .. code-block:: python {'site': 'hcp', 'surface': 'fcc211', 'morphology': 'sc-tc-h', 'position': array([ 4.51584136, 0.63816387, 12.86014042]), 'normal': array([-0.33333333, -0. , 0.94280904]), 'indices': (0, 2, 3), 'composition': 'AuAuCu', 'subsurf_index': 9, 'subsurf_element': 'Cu', 'label': 28} The following example illustrates the another important use of `SlabAdsorptionSites` object - getting all symmetry-inequivalent adsorption sites. A customized surface is used to represent the anatase TiO2(101) surface which is not directly implemented in ACAT: .. code-block:: python >>> from acat.adsorption_sites import SlabAdsorptionSites >>> from acat.settings import CustomSurface >>> from ase.spacegroup import crystal >>> from ase.build import surface >>> a, c = 3.862, 9.551 >>> anatase = crystal(['Ti', 'O'], basis=[(0 ,0, 0), (0, 0, 0.208)], ... spacegroup=141, cellpar=[a, a, c, 90, 90, 90]) >>> anatase_101_atoms = surface(anatase, (1, 0, 1), 4, vacuum=5) >>> anatase_101 = CustomSurface(anatase_101_atoms, n_layers=4) >>> # Suppose we have a surface structure resulting from some >>> # perturbation of the reference TiO2(101) surface structure >>> atoms = anatase_101_atoms.copy() >>> atoms.rattle(stdev=0.2) >>> sas = SlabAdsorptionSites(atoms, surface=anatase_101, ... composition_effect=True, ... ignore_sites=['4fold'], ... label_sites=True) >>> unq_sites = sas.get_unique_sites(unique_composition=True) >>> print(unq_sites) Output: .. code-block:: python [{'site': 'ontop', 'surface': CustomSurface, 'morphology': 'terrace', 'position': array([ 4.39474295, 3.84711082, 15.3573688 ]), 'normal': array([0.48335322, 0. , 0.87542542]), 'indices': (36,), 'composition': 'Ti', 'subsurf_index': None, 'subsurf_element': None, 'label': 4}, {'site': 'bridge', 'surface': CustomSurface, 'morphology': 'terrace', 'position': array([ 4.39474295, 1.91611082, 15.3573688 ]), 'normal': array([0.48335322, 0. , 0.87542542]), 'indices': (36, 36), 'composition': 'TiTi', 'subsurf_index': None, 'subsurf_element': None, 'label': 12}, {'site': '3fold', 'surface': CustomSurface, 'morphology': 'sc-tc', 'position': array([ 3.62687221, 1.94117703, 15.78133729]), 'normal': array([0.48335322, 0. , 0.87542542]), 'indices': (36, 36, 38), 'composition': 'TiTiTi', 'subsurf_index': None, 'subsurf_element': None, 'label': 35}, {'site': 'bridge', 'surface': CustomSurface, 'morphology': 'sc-tc', 'position': array([ 3.24293684, 2.91921014, 15.99332154]), 'normal': array([0.3520727 , 0. , 0.92537315]), 'indices': (36, 38), 'composition': 'TiTi', 'subsurf_index': None, 'subsurf_element': None, 'label': 18}, {'site': 'ontop', 'surface': CustomSurface, 'morphology': 'step', 'position': array([ 2.09113074, 1.99130947, 16.62927427]), 'normal': array([0.48335322, 0. , 0.87542542]), 'indices': (38,), 'composition': 'Ti', 'subsurf_index': None, 'subsurf_element': None, 'label': 2}, {'site': 'bridge', 'surface': CustomSurface, 'morphology': 'step', 'position': array([ 2.09113074, 0.06030947, 16.62927427]), 'normal': array([0.48335322, 0. , 0.87542542]), 'indices': (38, 38), 'composition': 'TiTi', 'subsurf_index': None, 'subsurf_element': None, 'label': 9}] """ def __init__(self, atoms, surface, allow_6fold=False, composition_effect=False, both_sides=False, ignore_sites=None, label_sites=False, surrogate_metal=None, optimize_surrogate_cell=False, tol=.5, _allow_expand=True): assert True in atoms.pbc, 'the cell must be periodic in at least one dimension' warnings.filterwarnings('ignore', category=RuntimeWarning) warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning) atoms = atoms.copy() ptp = np.ptp(atoms.positions[:, 2]) if np.linalg.norm(atoms.cell[2]) - ptp < 10.: atoms.cell[2][2] = ptp + 10. self.metal_ids = [a.index for a in atoms if a.symbol not in adsorbate_elements] atoms = Atoms(atoms.symbols[self.metal_ids], atoms.positions[self.metal_ids], cell=atoms.cell, pbc=atoms.pbc) self.atoms = atoms self.positions = atoms.positions self.symbols = atoms.symbols self.numbers = atoms.numbers self.surface = surface self.surrogate_metal = surrogate_metal self.optimize_surrogate_cell = optimize_surrogate_cell self.ref_atoms, self.delta_positions = self.mapping(atoms) self.cell = atoms.cell self.pbc = atoms.pbc self.metals = sorted(list(set(atoms.symbols))) self.allow_6fold = allow_6fold self.composition_effect = composition_effect self.both_sides = both_sides if (ignore_sites is None) or is_list_or_tuple(ignore_sites): self.ignore_sites = ignore_sites else: self.ignore_sites = [ignore_sites] self.label_sites = label_sites if self.composition_effect: if len(self.metals) <= 2: if len(self.metals) == 1: self.metals *= 2 self.label_dict = get_bimetallic_slab_labels(self.surface, self.metals) else: self.label_dict = get_multimetallic_slab_labels(self.surface, self.metals) else: self.label_dict = get_monometallic_slab_labels(self.surface) self.tol = tol self._allow_expand = _allow_expand self.make_neighbor_list(neighbor_number=1) self.adj_matrix = self.get_connectivity() self.surf_ids, self.subsurf_ids = self.get_termination() self.site_list = [] if len(self.surf_ids) < 4: self.populate_expanded_site_list() else: self.populate_site_list() self.postprocessing()
[docs] def populate_site_list(self): """Find all ontop, bridge and hollow sites (3-fold and 4-fold) given an input slab based on Delaunay triangulation of the surface atoms in a supercell and collect in a site list. """ top_indices = self.surf_ids sl = self.site_list normals_for_site = dict(list(zip(top_indices, [[] for _ in top_indices]))) ref_cell = self.ref_atoms.cell usi = set() # used_site_indices cm = self.adj_matrix if isinstance(self.surface, CustomSurface): surface_name = 'CustomSurface' else: surface_name = self.surface # Sort top indicies by z coordinates to identify step / terrace / corner sorted_top_indices = sorted(top_indices, key=lambda x: self.positions[x,2]) ntop = len(top_indices) stepids, terraceids, cornerids = set(), set(), set() for i, s in enumerate(sorted_top_indices): if isinstance(self.surface, CustomSurface): if i < self.surface.n_corners: morphology = 'corner' elif i < self.surface.n_corners + self.surface.n_terraces: morphology = 'terrace' else: morphology = 'step' elif self.surface in ['fcc111','fcc100','bcc100','bcc110','hcp0001']: morphology = 'terrace' elif self.surface in ['fcc211','fcc322','fcc221','fcc332','bcc111','bcc210','hcp10m12']: if i < ntop / 3: morphology = 'corner' elif i >= ntop * 2 / 3: morphology = 'step' else: morphology = 'terrace' else: if i < ntop / 2: morphology = 'terrace' else: morphology = 'step' si = (s,) site = self.new_site() site.update({'site': 'ontop', 'surface': surface_name, 'morphology': morphology, 'position': np.round(self.positions[s], 8), 'indices': si}) if self.surface in ['fcc110','bcc211','hcp10m10h'] and morphology == 'terrace': site.update({'extra': np.where(cm[s]==1)[0]}) if self.composition_effect: site.update({'composition': self.symbols[s]}) if morphology == 'step': stepids.add(s) elif morphology == 'terrace': terraceids.add(s) elif morphology == 'corner': cornerids.add(s) sl.append(site) usi.add(si) if self.surface in ['fcc110','bcc211','hcp10m10h']: geo_dict = {'step': stepids, 'terrace': terraceids, 'corner': cornerids} for st in sl: if 'extra' in st: si = st['indices'] extra = st['extra'] extraids = [e for e in extra if e in self.surf_ids and e not in geo_dict[st['morphology']]] #if len(extraids) < 4: # warnings.warn('Cannot identify other 4 atoms of 5-fold site {}'.format(si)) if len(extraids) > 4: extraids = sorted(extraids, key=lambda x: get_mic( self.ref_atoms.positions[x], refpos, ref_cell, return_squared_distance=True))[:4] if self.composition_effect: metals = self.metals if len(metals) == 1: composition = 4*metals[0] elif len(metals) == 2: ma, mb = metals[0], metals[1] symbols = self.symbols[extraids] nma = np.count_nonzero(symbols==ma) if nma == 0: composition = 4*mb elif nma == 1: composition = ma + 3*mb elif nma == 2: idd = sorted(extraids[1:], key=lambda x: get_mic(self.positions[x], self.positions[extraids[0]], self.cell, return_squared_distance=True)) opp, close = idd[-1], idd[0] if self.symbols[opp] == self.symbols[extraids[0]]: composition = ma + mb + ma + mb else: if self.symbols[close] == self.symbols[extraids[0]]: composition = 2*ma + 2*mb else: composition = ma + 2*mb + ma elif nma == 3: composition = 3*ma + mb elif nma == 4: composition = 4*ma else: opp = max(list(extraids[1:]), key=lambda x: get_mic(self.positions[x], self.positions[extraids[0]], self.cell, return_squared_distance=True)) nni = [i for i in extraids[1:] if i != opp] nodes = self.symbols[[extraids[0], nni[0], opp, nni[1]]] composition = ''.join(hash_composition(nodes)) st['composition'] += '-{}'.format(composition) st['indices'] = tuple(sorted(list(si)+ extraids)) del st['extra'] meansurfz = np.average(self.positions[self.surf_ids][:,2], 0) dh = abs(meansurfz - np.average(self.positions[self.subsurf_ids][:,2], 0)) bridge_positions, fold3_positions, fold4_positions = self.delaunay_triangulation() fold4_surfaces = ['fcc100','fcc211','fcc311','fcc322','fcc331','bcc100', 'bcc210','bcc310','hcp10m10t','hcp10m11','hcp10m12'] # Complete information of each site for n, poss in enumerate([bridge_positions,fold4_positions,fold3_positions]): if not poss: continue fracs = np.stack(poss, axis=0) @ np.linalg.pinv(ref_cell) xfracs, yfracs = fracs[:,0], fracs[:,1] # Take only the positions within the periodic boundary screen = np.where((xfracs > 0 - 1e-5) & (xfracs < 1 - 1e-5) & \ (yfracs > 0 - 1e-5) & (yfracs < 1 - 1e-5))[0] reduced_poss = np.asarray(poss)[screen] # Sort the index list of surface and subsurface atoms # so that we can retrive original indices later top2_indices = np.asarray(self.surf_ids + self.subsurf_ids) top2atoms = self.ref_atoms[top2_indices] ntop1, ntop2 = len(self.surf_ids), len(top2atoms) if isinstance(self.surface, CustomSurface): pairs = np.asarray(list(product(range(len(reduced_poss)),range(len(top2atoms))))) D = find_mic(top2atoms.positions[pairs[:,1]] - reduced_poss[pairs[:,0]], ref_cell, self.pbc)[1].reshape(len(reduced_poss),len(top2atoms)) min2_rows = np.argpartition(D, 2)[:,:2] min3_rows = np.argpartition(D[:,:ntop1], 3)[:,:3] if ntop1 == 4: min4_rows = np.tile([0, 1, 2, 3], (len(reduced_poss),1)) else: min4_rows = np.argpartition(D[:,:ntop1], 4)[:,:4] else: # Make a new neighbor list including sites as dummy atoms dummies = Atoms('X{}'.format(reduced_poss.shape[0]), positions=reduced_poss, cell=ref_cell, pbc=self.pbc) testatoms = top2atoms + dummies nblist = neighbor_shell_list(testatoms, dx=self.tol, neighbor_number=1, different_species=True, mic=True) # Make bridge sites if n == 0: fold4_poss = [] for i, refpos in enumerate(reduced_poss): if isinstance(self.surface, CustomSurface): bridge_indices = min2_rows[i] bridgeids = list(top2_indices[bridge_indices]) else: bridge_indices = nblist[ntop2+i] if not bridge_indices: continue bridgeids = list(top2_indices[[j for j in bridge_indices if j < ntop1]]) if len(bridgeids) != 2: if self.surface in ['fcc100','fcc211','fcc311','fcc322','fcc331', 'bcc100','bcc210','bcc310','hcp10m11','hcp10m12']: fold4_poss.append(refpos) continue else: bridgeids = sorted(bridgeids, key=lambda x: get_mic( self.ref_atoms.positions[x], refpos, ref_cell, return_squared_distance=True))[:2] si = tuple(sorted(bridgeids)) if self.optimize_surrogate_cell or isinstance(self.surface, CustomSurface): reffrac = refpos @ np.linalg.pinv(ref_cell) pos = (reffrac + np.average(self.delta_positions[bridgeids], 0)) @ self.cell else: pos = refpos + np.average(self.delta_positions[bridgeids], 0) occurence = np.sum(cm[bridgeids], axis=0) siset = set(si) nstep = len(stepids.intersection(siset)) nterrace = len(terraceids.intersection(siset)) if isinstance(self.surface, CustomSurface): this_site = 'bridge' ncorner = len(cornerids.intersection(siset)) if nstep == 2: morphology = 'step' elif nterrace == 2: morphology = 'terrace' elif ncorner == 2: morphology = 'corner' elif nstep == 1 and nterrace == 1: morphology = 'sc-tc' elif nstep == 1 and ncorner == 1: morphology = 'sc-cc' elif nterrace == 1 and ncorner == 1: morphology = 'tc-cc' else: morphology = 'subsurf' elif self.surface in ['fcc110','hcp10m10h']: this_site = 'bridge' if nstep == 2: morphology = 'step' elif nstep == 1 and nterrace == 1: morphology = 'sc-tc-h' elif nterrace == 2: morphology = 'terrace' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif self.surface in ['fcc311','fcc331']: this_site = 'bridge' if nstep == 2: morphology = 'step' elif nstep == 1 and nterrace == 1: cto2 = np.count_nonzero(occurence==2) if cto2 == 2: morphology = 'sc-tc-t' elif cto2 in [3, 4]: morphology = 'sc-tc-h' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif nterrace == 2: morphology = 'terrace' else: morphology = 'subsurf' elif self.surface in ['fcc211','fcc322']: this_site = 'bridge' ncorner = len(cornerids.intersection(siset)) if nstep == 2: morphology = 'step' elif nstep == 1 and nterrace == 1: morphology = 'sc-tc-h' elif nstep == 1 and ncorner == 1: morphology = 'sc-cc-t' elif nterrace == 1 and ncorner == 1: morphology = 'tc-cc-h' elif ncorner == 2: morphology = 'corner' # nterrace == 2 is actually terrace bridge, # but equivalent to tc-cc-h for fcc211 elif nterrace == 2: if self.surface == 'fcc211': morphology = 'tc-cc-h' elif self.surface == 'fcc322': morphology = 'terrace' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif self.surface in ['fcc221','fcc332']: this_site = 'bridge' ncorner = len(cornerids.intersection(siset)) if nstep == 2: morphology = 'step' elif nstep == 1 and nterrace == 1: morphology = 'sc-tc-h' elif nstep == 1 and ncorner == 1: morphology = 'sc-cc-h' elif nterrace == 1 and ncorner == 1: morphology = 'tc-cc-h' elif ncorner == 2: morphology = 'corner' elif nterrace == 2: morphology = 'terrace' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif self.surface == 'bcc110': this_site, morphology = 'shortbridge', 'terrace' elif self.surface == 'bcc111': ncorner = len(cornerids.intersection(siset)) if nstep == 1 and ncorner == 1: this_site, morphology = 'longbridge', 'sc-cc-o' elif nstep == 1 and nterrace == 1: this_site, morphology = 'shortbridge', 'sc-tc-o' elif nterrace == 1 and ncorner == 1: this_site, morphology = 'shortbridge', 'tc-cc-o' elif nstep == 2 or nterrace == 2 or ncorner == 2: continue else: warnings.warn('Cannot identify site {}'.format(si)) continue elif self.surface == 'bcc210': this_site = 'bridge' ncorner = len(cornerids.intersection(siset)) if nstep == 2: morphology = 'step' elif nstep == 1 and nterrace == 1: morphology = 'sc-tc-o' elif nstep == 1 and ncorner == 1: morphology = 'sc-cc-t' elif nterrace == 1 and ncorner == 1: morphology = 'tc-cc-o' elif ncorner == 2: morphology = 'corner' # nterrace == 2 is terrace bridge and not # equivalent to tc-cc-o for bcc210 elif nterrace == 2: morphology = 'terrace' else: warnings.warn('Cannot identify site {}'.format(si)) continue if self.surface == 'bcc211': this_site = 'bridge' if nstep == 2: morphology = 'step' elif nstep == 1 and nterrace == 1: morphology = 'sc-tc-o' elif nterrace == 2: morphology = 'terrace' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif self.surface == 'bcc310': this_site = 'bridge' if nstep == 2: morphology = 'step' elif nstep == 1 and nterrace == 1: cto2 = np.count_nonzero(occurence==2) if cto2 == 2: morphology = 'sc-tc-t' elif cto2 == 3: morphology = 'sc-tc-o' else: warnings.warn('Cannot identify site {}'.format(si)) elif nterrace == 2: morphology = 'terrace' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif self.surface == 'hcp10m10t': this_site = 'bridge' if nstep == 2: morphology = 'step' elif nstep == 1 and nterrace == 1: morphology = 'sc-tc-t' elif nterrace == 2: morphology = 'terrace' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif self.surface == 'hcp10m11': this_site = 'bridge' if nstep == 2: morphology = 'step' elif nstep == 1 and nterrace == 1: cto2 = np.count_nonzero(occurence==2) if cto2 == 2: morphology = 'subsurf' isubs = [self.subsurf_ids[i] for i in np.where( occurence[self.subsurf_ids] == 2)[0]] subpos = self.positions[isubs[0]] + .5 * get_mic( self.positions[isubs[0]], self.positions[ isubs[1]], self.cell) elif cto2 == 3: morphology = 'sc-tc-h' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif nterrace == 2: morphology = 'terrace' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif self.surface == 'hcp10m12': this_site = 'bridge' ncorner = len(cornerids.intersection(siset)) if nstep == 2: morphology = 'step' elif nstep == 1 and nterrace == 1: morphology = 'sc-tc-h' elif nstep == 1 and ncorner == 1: morphology = 'sc-cc-h' elif nterrace == 1 and ncorner == 1: morphology = 'tc-cc-t' elif ncorner == 2: morphology = 'corner' elif nterrace == 2: morphology = 'terrace' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif self.surface in ['fcc111','fcc100','bcc100','hcp0001']: this_site, morphology = 'bridge', 'terrace' site = self.new_site() special = False if self.surface in ['fcc110','bcc211','hcp10m10h'] and morphology == 'terrace': special = True extraids = [xi for xi in np.where(occurence==2)[0] if xi in self.surf_ids] #if len(extraids) < 2: # warnings.warn('Cannot identify other 2 atoms of 4-fold site {}'.format(si)) if len(extraids) > 2: extraids = sorted(extraids, key=lambda x: get_mic( self.ref_atoms.positions[x], refpos, ref_cell, return_squared_distance=True))[:2] site.update({'site': this_site, 'surface': surface_name, 'morphology': morphology, 'position': np.round(pos, 8), 'indices': tuple(sorted(bridgeids + extraids))}) elif self.surface == 'hcp10m11' and morphology == 'subsurf': special = True extraids = si site.update({'site': this_site, 'surface': surface_name, 'morphology': morphology, 'position': np.round(subpos, 8), 'indices': tuple(sorted(bridgeids + isubs))}) else: site.update({'site': this_site, 'surface': surface_name, 'morphology': morphology, 'position': np.round(pos, 8), 'indices': si}) if self.composition_effect: if self.surface == 'hcp10m11' and morphology == 'subsurf': composition = ''.join(sorted(self.symbols[isubs])) else: composition = ''.join(sorted(self.symbols[list(si)])) if special: extracomp = ''.join(sorted(self.symbols[list(extraids)])) composition += '-{}'.format(extracomp) site.update({'composition': composition}) sl.append(site) usi.add(si) if self.surface in fold4_surfaces and fold4_poss: fold4atoms = Atoms('X{}'.format(len(fold4_poss)), positions=np.asarray(fold4_poss), cell=ref_cell, pbc=self.pbc) topatoms = self.ref_atoms[self.surf_ids] newatoms = topatoms + fold4atoms newnblist = neighbor_shell_list(newatoms, dx=.1, neighbor_number=2, different_species=True, mic=True) # Make 4-fold hollow sites for i, refpos in enumerate(fold4_poss): fold4_indices = newnblist[ntop1+i] fold4ids = [self.surf_ids[j] for j in fold4_indices] if len(fold4ids) < 4: continue elif len(fold4ids) > 4: fold4ids = sorted(fold4ids, key=lambda x: get_mic( self.ref_atoms.positions[x], refpos, ref_cell, return_squared_distance=True))[:4] occurence = np.sum(cm[fold4ids], axis=0) isub = np.where(occurence >= 4)[0] isub = [i for i in isub if i in self.subsurf_ids] if len(isub) == 0: continue isub = isub[0] si = tuple(sorted(fold4ids)) if self.optimize_surrogate_cell: reffrac = refpos @ np.linalg.pinv(ref_cell) pos = (reffrac + np.average(self.delta_positions[fold4ids], 0)) @ self.cell else: pos = refpos + np.average(self.delta_positions[fold4ids], 0) normal = self.get_surface_normal([si[0], si[1], si[2]]) for idx in si: normals_for_site[idx].append(normal) site = self.new_site() if self.surface in ['hcp10m10t','hcp10m11']: this_site = '5fold' site.update({'site': this_site, 'surface': surface_name, 'morphology': 'subsurf', 'position': np.round(self.positions[isub], 8), 'normal': np.round(normal, 8), 'indices': tuple(sorted(fold4ids+[isub]))}) else: if self.surface in ['fcc211','fcc322','bcc210']: morphology = 'sc-cc-t' elif self.surface in ['fcc311','fcc331','bcc310']: morphology = 'sc-tc-t' elif self.surface == 'hcp10m12': morphology = 'tc-cc-t' else: morphology = 'terrace' this_site = '4fold' site.update({'site': this_site, 'surface': surface_name, 'morphology': morphology, 'position': np.round(pos, 8), 'normal': np.round(normal, 8), 'indices': si}) if self.composition_effect: metals = self.metals if len(metals) == 1: composition = 4*metals[0] elif len(metals) == 2: ma, mb = metals[0], metals[1] symbols = self.symbols[fold4ids] nma = np.count_nonzero(symbols==ma) if nma == 0: composition = 4*mb elif nma == 1: composition = ma + 3*mb elif nma == 2: if this_site == '5fold': idd = sorted(fold4ids[1:], key=lambda x: get_mic(self.positions[x], self.positions[fold4ids[0]], self.cell, return_squared_distance=True)) opp, close = idd[-1], idd[0] if self.symbols[opp] == self.symbols[fold4ids[0]]: composition = ma + mb + ma + mb else: if self.symbols[close] == self.symbols[fold4ids[0]]: composition = 2*ma + 2*mb else: composition = ma + 2*mb + ma else: opposite = np.where(cm[si[1:],si[0]]==0)[0] opp = si[1+opposite[0]] if self.symbols[opp] == self.symbols[fold4ids[0]]: composition = ma + mb + ma + mb else: composition = 2*ma + 2*mb elif nma == 3: composition = 3*ma + mb elif nma == 4: composition = 4*ma else: if this_site == '5fold': opp = max(list(fold4ids[1:]), key=lambda x: get_mic(self.positions[x], self.positions[fold4ids[0]], self.cell, return_squared_distance=True)) nni = [i for i in fold4ids[1:] if i != opp] nodes = self.symbols[[fold4ids[0], nni[0], opp, nni[1]]] else: opposite = np.where(cm[si[1:],si[0]]==0)[0] opp = si[1+opposite[0]] nni = [i for i in si[1:] if i != opp] nodes = self.symbols[[si[0], nni[0], opp, nni[1]]] composition = ''.join(hash_composition(nodes)) site.update({'composition': composition}) if this_site == '5fold': site['composition'] = '{}-'.format( self.symbols[isub]) + site['composition'] if this_site != '5fold': site.update({'subsurf_index': isub}) if self.composition_effect: site.update({'subsurf_element': self.symbols[isub]}) sl.append(site) usi.add(si) # Make 3-fold hollow sites (differentiate fcc / hcp) if n == 2 and self.surface not in ['fcc100','bcc100','hcp10m10t']: coexist_3_4 = (self.surface in ['fcc211','fcc311','fcc322','fcc331','bcc210', 'bcc310','hcp10m12'] or isinstance(self.surface, CustomSurface)) if coexist_3_4: fold4_sites = [s for s in sl if s['site'] == '4fold'] for i, refpos in enumerate(reduced_poss): if isinstance(self.surface, CustomSurface): fold3_indices = min3_rows[i] fold3ids = [self.surf_ids[j] for j in fold3_indices] else: fold3_indices = nblist[ntop2+i] fold3ids = list(top2_indices[[j for j in fold3_indices if j < ntop1]]) if len(fold3ids) != 3: #if self.surface != 'hcp10m11': # si = tuple(sorted(fold3ids)) # warnings.warn('Cannot find the correct atoms of this 3-fold site.', # 'Find {} instead'.format(si)) continue si = tuple(sorted(fold3ids)) if self.optimize_surrogate_cell or isinstance(self.surface, CustomSurface): reffrac = refpos @ np.linalg.pinv(ref_cell) pos = (reffrac + np.average(self.delta_positions[fold3ids], 0)) @ self.cell else: pos = refpos + np.average(self.delta_positions[fold3ids], 0) # Remove redundant 3-fold sites that belongs to 4-fold sites if coexist_3_4: dup_sites = [s for s in fold4_sites if set(fold3ids).issubset(set(s['indices']))] if dup_sites: if any(get_mic(s['position'], pos, self.cell, return_squared_distance=True) < 1.5 for s in dup_sites): continue normal = self.get_surface_normal([si[0], si[1], si[2]]) for idx in si: normals_for_site[idx].append(normal) occurence = np.sum(cm[fold3ids], axis=0) if isinstance(self.surface, CustomSurface): this_site = '3fold' siset = set(si) step_overlap = stepids.intersection(siset) terrace_overlap = terraceids.intersection(siset) corner_overlap = cornerids.intersection(siset) if step_overlap and terrace_overlap and not corner_overlap: morphology = 'sc-tc' elif terrace_overlap and corner_overlap and not step_overlap: morphology = 'tc-cc' elif step_overlap and corner_overlap and not terrace_overlap: morphology = 'sc-cc' elif terrace_overlap and not step_overlap and not corner_overlap: morphology = 'terrace' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif self.surface in ['fcc211','fcc221','fcc322','fcc332', 'hcp10m11','hcp10m12']: if np.max(occurence[self.subsurf_ids]) == 3: this_site = 'hcp' else: this_site = 'fcc' siset = set(si) step_overlap = stepids.intersection(siset) corner_overlap = cornerids.intersection(siset) if step_overlap and not corner_overlap: morphology = 'sc-tc-h' elif corner_overlap and not step_overlap: morphology = 'tc-cc-h' elif step_overlap and corner_overlap: morphology = 'sc-cc-h' else: morphology = 'terrace' elif self.surface == 'bcc210': this_site = '3fold' siset = set(si) step_overlap = stepids.intersection(siset) corner_overlap = cornerids.intersection(siset) if step_overlap and not corner_overlap: morphology = 'sc-tc-o' elif corner_overlap and not step_overlap: morphology = 'tc-cc-o' else: morphology = 'sc-tc-o' elif self.surface == 'bcc110': this_site, morphology = '3fold', 'terrace' elif self.surface == 'bcc111': this_site, morphology = '3fold', 'sc-tc-cc-o' elif self.surface in ['bcc211','bcc310']: this_site, morphology = '3fold', 'sc-tc-o' elif self.surface in ['fcc111','hcp0001']: morphology = 'terrace' if np.max(occurence[self.subsurf_ids]) == 3: this_site = 'hcp' else: this_site = 'fcc' elif self.surface in ['fcc110','fcc311','fcc331','hcp10m10h']: morphology = 'sc-tc-h' if np.max(occurence[self.subsurf_ids]) == 3: this_site = 'hcp' else: this_site = 'fcc' site = self.new_site() site.update({'site': this_site, 'surface': surface_name, 'morphology': morphology, 'position': np.round(pos, 8), 'normal': np.round(normal, 8), 'indices': si}) if self.composition_effect: metals = self.metals if len(metals) == 1: composition = 3*metals[0] elif len(metals) == 2: ma, mb = metals[0], metals[1] symbols = self.symbols[list(si)] nma = np.count_nonzero(symbols==ma) if nma == 0: composition = 3*mb elif nma == 1: composition = ma + 2*mb elif nma == 2: composition = 2*ma + mb elif nma == 3: composition = 3*ma else: nodes = self.symbols[list(si)] composition = ''.join(hash_composition(nodes)) site.update({'composition': composition}) if this_site == 'hcp': isub = np.where(occurence == 3)[0] isub = [i for i in isub if i in self.subsurf_ids] if len(isub) == 0: continue isub = isub[0] spos = pos - normal * dh if get_mic(self.positions[isub], spos, self.cell, return_squared_distance=True) > 2.: site.update({'site': 'fcc'}) else: site.update({'subsurf_index': isub}) if self.composition_effect: site.update({'subsurf_element': self.symbols[isub]}) sl.append(site) usi.add(si) # Make 4-fold hollow sites if n == 1 and (self.surface in fold4_surfaces or isinstance( self.surface, CustomSurface)) and list(reduced_poss): fold4atoms = Atoms('X{}'.format(len(reduced_poss)), positions=np.asarray(reduced_poss), cell=ref_cell, pbc=self.pbc) if not isinstance(self.surface, CustomSurface): topatoms = self.ref_atoms[self.surf_ids] newatoms = topatoms + fold4atoms newnblist = neighbor_shell_list(newatoms, dx=.1, neighbor_number=2, different_species=True, mic=True) for i, refpos in enumerate(reduced_poss): if isinstance(self.surface, CustomSurface): fold4_indices = min4_rows[i] else: fold4_indices = newnblist[ntop1+i] fold4ids = [self.surf_ids[j] for j in fold4_indices] if len(fold4ids) < 4: continue elif len(fold4ids) > 4: fold4ids = sorted(fold4ids, key=lambda x: get_mic( self.ref_atoms.positions[x], refpos, ref_cell, return_squared_distance=True))[:4] occurence = np.sum(cm[fold4ids], axis=0) isub = np.where(occurence == 4)[0] isub = [i for i in isub if i in self.subsurf_ids] if len(isub) == 0: isub = None else: isub = isub[0] si = tuple(sorted(fold4ids)) if self.optimize_surrogate_cell or isinstance(self.surface, CustomSurface): reffrac = refpos @ np.linalg.pinv(ref_cell) pos = (reffrac + np.average(self.delta_positions[fold4ids], 0)) @ self.cell else: pos = refpos + np.average(self.delta_positions[fold4ids], 0) normal = self.get_surface_normal([si[0], si[1], si[2]]) for idx in si: normals_for_site[idx].append(normal) site = self.new_site() if self.surface in ['hcp10m10t','hcp10m11']: if isub is None: continue this_site = '5fold' site.update({'site': this_site, 'surface': surface_name, 'morphology': 'subsurf', 'position': np.round(self.positions[isub], 8), 'normal': np.round(normal, 8), 'indices': tuple(sorted(fold4ids+[isub]))}) else: this_site = '4fold' if isinstance(self.surface, CustomSurface): siset = set(si) step_overlap = stepids.intersection(siset) terrace_overlap = terraceids.intersection(siset) corner_overlap = cornerids.intersection(siset) if step_overlap and terrace_overlap and not corner_overlap: morphology = 'sc-tc' elif terrace_overlap and corner_overlap and not step_overlap: morphology = 'tc-cc' elif step_overlap and corner_overlap and not terrace_overlap: morphology = 'sc-cc' elif terrace_overlap and not step_overlap and not corner_overlap: morphology = 'terrace' else: warnings.warn('Cannot identify site {}'.format(si)) continue elif self.surface in ['fcc211','fcc322','bcc210']: morphology = 'sc-cc-t' elif self.surface in ['fcc311','fcc331','bcc310']: morphology = 'sc-tc-t' elif self.surface == 'hcp10m12': morphology = 'tc-cc-t' else: morphology = 'terrace' site.update({'site': this_site, 'surface': surface_name, 'morphology': morphology, 'position': np.round(pos, 8), 'normal': np.round(normal, 8), 'indices': si}) if self.composition_effect: metals = self.metals if len(metals) == 1: composition = 4*metals[0] elif len(metals) == 2: ma, mb = metals[0], metals[1] symbols = self.symbols[list(si)] nma = np.count_nonzero(symbols==ma) if nma == 0: composition = 4*mb elif nma == 1: composition = ma + 3*mb elif nma == 2: if this_site == '5fold': idd = sorted(fold4ids[1:], key=lambda x: get_mic(self.positions[x], self.positions[fold4ids[0]], self.cell, return_squared_distance=True)) opp, close = idd[-1], idd[0] if self.symbols[opp] == self.symbols[fold4ids[0]]: composition = ma + mb + ma + mb else: if self.symbols[close] == self.symbols[fold4ids[0]]: composition = 2*ma + 2*mb else: composition = ma + 2*mb + ma else: opposite = np.where(cm[si[1:],si[0]]==0)[0] opp = si[1+opposite[0]] if self.symbols[opp] == self.symbols[fold4ids[0]]: composition = ma + mb + ma + mb else: composition = 2*ma + 2*mb elif nma == 3: composition = 3*ma + mb elif nma == 4: composition = 4*ma else: if this_site == '5fold': opp = max(list(fold4ids[1:]), key=lambda x: get_mic(self.positions[x], self.positions[fold4ids[0]], self.cell, return_squared_distance=True)) nni = [i for i in fold4ids[1:] if i != opp] nodes = self.symbols[[fold4ids[0], nni[0], opp, nni[1]]] else: opposite = np.where(cm[si[1:],si[0]]==0)[0] opp = si[1+opposite[0]] nni = [i for i in si[1:] if i != opp] nodes = self.symbols[[si[0], nni[0], opp, nni[1]]] composition = ''.join(hash_composition(nodes)) site.update({'composition': composition}) if this_site == '5fold': if isub is None: continue site['composition'] = '{}-'.format( self.symbols[isub]) + site['composition'] if this_site != '5fold': site.update({'subsurf_index': isub}) if self.composition_effect and (isub is not None): site.update({'subsurf_element': self.symbols[isub]}) sl.append(site) usi.add(si) if self.surface == 'bcc110': bcc110_long_bridges = [] for t in sl: stids = t['indices'] this_site = t['site'] # Add normals to ontop sites if this_site == 'ontop': n = np.average(normals_for_site[stids[0]], 0) t['normal'] = np.round(n / np.linalg.norm(n), 8) nstids = len(stids) if nstids > 1: t['site'] = '{}fold'.format(nstids) t['indices'] = tuple(sorted(stids)) # Add normals to bridge sites elif 'bridge' in this_site: if t['morphology'] in ['terrace', 'step', 'corner']: normals = [] for i in stids[:2]: normals.extend(normals_for_site[i]) n = np.average(normals, 0) t['normal'] = np.round(n / np.linalg.norm(n), 8) nstids = len(stids) if nstids > 2: t['site'] = '{}fold'.format(nstids) t['indices'] = tuple(sorted(stids)) else: if self.surface == 'hcp10m10t': normals = [s['normal'] for s in sl if (s['site'] in ['3fold','4fold','fcc','hcp']) and (s['morphology'] == 'subsurf') and ( len(set(s['indices']).intersection( set(t['indices']))) == 2)] else: normals = [s['normal'] for s in sl if (s['site'] in ['3fold','4fold','fcc','hcp']) and ( len(set(s['indices']).intersection( set(t['indices']))) == 2)] t['normal'] = np.round(np.average(normals, 0), 8) nstids = len(stids) if nstids > 2: t['site'] = '{}fold'.format(nstids) t['indices'] = tuple(sorted(stids)) # Take care of longbridge sites on bcc110 elif self.surface == 'bcc110' and this_site == '3fold': si = t['indices'] pairs = [(si[0], si[1]), (si[0], si[2]), (si[1], si[2])] longest = max(pairs, key=lambda x: get_mic( self.ref_atoms.positions[x[0]], self.ref_atoms.positions[x[1]], ref_cell, return_squared_distance=True)) bcc110_long_bridges.append(longest) if self.surface == 'bcc110': for st in sl: if st['site'] == 'shortbridge': if st['indices'] in bcc110_long_bridges: st['site'] = 'longbridge' # Add 6-fold sites if allowed if self.allow_6fold and not isinstance(self.surface, CustomSurface): for st in sl: if self.surface in ['fcc110','hcp10m10h'] and st['site'] == 'bridge' \ and st['morphology'] == 'step': site = st.copy() subpos = st['position'] - st['normal'] * dh / 2 sid = list(st['indices']) occurence = cm[sid] surfsi = sid + list(np.where(np.sum(occurence, axis=0) == 2)[0]) if len(surfsi) > 4: surfsi = sorted(surfsi, key=lambda x: get_mic( self.positions[x], subpos, self.cell, return_squared_distance=True))[:4] subsi = [self.subsurf_ids[i] for i in np.where(np.sum( occurence[:,self.subsurf_ids], axis=0) == 1)[0]] if len(subsi) > 2: subsi = sorted(subsi, key=lambda x: get_mic( self.positions[x], subpos, self.cell, return_squared_distance=True))[:2] si = tuple(sorted(surfsi + subsi)) normal = np.asarray([0.,0.,1.]) site.update({'site': '6fold', 'position': np.round(subpos, 8), 'morphology': 'subsurf', 'normal': np.round(normal, 8), 'indices': si}) elif st['site'] == 'fcc': site = st.copy() subpos = st['position'] - st['normal'] * dh / 2 def get_squared_distance(x): return get_mic(self.positions[x], subpos, self.cell, return_squared_distance=True) subsi = sorted(sorted(self.subsurf_ids, key=get_squared_distance)[:3]) si = site['indices'] site.update({'site': '6fold', 'position': np.round(subpos, 8), 'morphology': 'subsurf', 'indices': tuple(sorted(si+tuple(subsi)))}) else: continue # Find the opposite element if self.composition_effect: metals = self.metals if self.surface in ['fcc110','hcp10m10h']: newsi = surfsi[:-1] subsi.append(surfsi[-1]) metals = self.metals if len(metals) == 1: comp = 3*metals[0] elif len(metals) == 2: ma, mb = metals[0], metals[1] symbols = self.symbols[newsi] nma = np.count_nonzero(symbols==ma) if nma == 0: comp = 3*mb elif nma == 1: comp = ma + 2*mb elif nma == 2: comp = 2*ma + mb elif nma == 3: comp = 3*ma else: nodes = self.symbols[list(si)] comp = ''.join(hash_composition(nodes)) else: comp = site['composition'] def get_squared_distance(x): return get_mic(self.positions[x], subpos, self.cell, return_squared_distance=True) if len(metals) == 1: composition = ''.join([comp, 3*metals[0]]) elif len(metals) == 2: ma, mb = metals[0], metals[1] subsyms = self.symbols[subsi] nma = np.count_nonzero(subsyms==ma) if nma == 0: composition = ''.join([comp, 3*mb]) elif nma == 1: ia = subsi[subsyms.index(ma)] subpos = self.positions[ia] if self.symbols[max(si, key=get_squared_distance)] == ma: composition = ''.join([comp, 2*mb + ma]) else: composition = ''.join([comp, ma + 2*mb]) elif nma == 2: ib = subsi[subsyms.index(mb)] subpos = self.positions[ib] if self.symbols[max(si, key=get_squared_distance)] == mb: composition = ''.join([comp, 2*ma + mb]) else: composition = ''.join([comp, mb + 2*ma]) elif nma == 3: composition = ''.join([comp, 3*ma]) else: endsym = re.findall('[A-Z][^A-Z]*', comp)[-1] subpos = [self.positions[i] for i in si if self.symbols[i] == endsym][0] nodes = self.symbols[sorted(subsi, key=get_squared_distance)] composition = comp + ''.join(hash_composition(nodes)) site.update({'composition': composition}) sl.append(site) usi.add(si)
[docs] def postprocessing(self): """Postprocessing the site list, and potentially adding more sites.""" # Check if the cell is so small that there are sites with duplicate indices small = False if not self._allow_expand: small = False elif self.surface in ['fcc211','fcc311','fcc322','fcc331','bcc210','bcc310','hcp10m12']: sorted_sets = sorted([set(s['indices']) for s in self.site_list if s['site'] in ['fcc','hcp','3fold','4fold']], key=len, reverse=True) for i in range(1, len(sorted_sets)): for s, t in zip(sorted_sets, sorted_sets[-i:]): if s.issubset(t): small = True break else: continue break elif len(set(s['indices'] for s in self.site_list)) < len(self.site_list): small = True if small: self.populate_expanded_site_list() else: if self.both_sides: self.populate_opposite_site_list() if self.ignore_sites is not None: self.site_list = [s for s in self.site_list if s['site'] not in self.ignore_sites] if self.label_sites: self.get_labels() # Wrap all sites in the box site_positions = np.asarray([s['position'] for s in self.site_list]) site_positions = wrap_positions(site_positions, self.cell, self.pbc) for i, st in enumerate(self.site_list): st['position'] = site_positions[i] st['indices'] = tuple(self.metal_ids[j] for j in st['indices']) self.site_list.sort(key=lambda x: x['indices'])
[docs] def populate_opposite_site_list(self): """Collect the sites on the opposite side of the slab.""" atoms = self.atoms.copy() atoms.positions *= [1,1,-1] if isinstance(self.surface, CustomSurface): ori_surface = deepcopy(self.surface) new_ref_atoms = self.surface.ref_atoms.copy() new_ref_atoms.positions *= [1,1,-1] new_surface = CustomSurface(new_ref_atoms, self.surface.n_layers, self.surface.tol) else: new_surface = self.surface nsas = SlabAdsorptionSites(atoms, surface=new_surface, allow_6fold=self.allow_6fold, composition_effect=self.composition_effect, both_sides=False, ignore_sites=self.ignore_sites, label_sites=self.label_sites, surrogate_metal=self.surrogate_metal, optimize_surrogate_cell=self.optimize_surrogate_cell, tol=self.tol, _allow_expand=False) # Take only the site positions within the periodic boundary bot_site_list = nsas.site_list for st in bot_site_list: if (type(st['normal']) is np.ndarray) and (st['normal'][2] > 0): st['normal'] *= [1,1,-1] st['position'] *= [1,1,-1] self.site_list += bot_site_list self.surf_ids += nsas.surf_ids self.subsurf_ids += nsas.subsurf_ids
[docs] def populate_expanded_site_list(self): """Collect the sites on the 2x2 expanded surface and cell for small unit cells and then return the sites within the original unit cell.""" atoms = self.atoms.copy() atoms.wrap() atoms *= [2,2,1] if isinstance(self.surface, CustomSurface): ori_surface = deepcopy(self.surface) new_ref_atoms = self.surface.ref_atoms * [2,2,1] new_surface = CustomSurface(new_ref_atoms, self.surface.n_layers, self.surface.tol) else: new_surface = self.surface nsas = SlabAdsorptionSites(atoms, surface=new_surface, allow_6fold=self.allow_6fold, composition_effect=self.composition_effect, both_sides=self.both_sides, ignore_sites=self.ignore_sites, label_sites=self.label_sites, surrogate_metal=self.surrogate_metal, optimize_surrogate_cell=self.optimize_surrogate_cell, tol=self.tol, _allow_expand=(len(self.surf_ids)==1)) # Take only the site positions within the periodic boundary sl = nsas.site_list poss = np.stack([s['position'] for s in sl], axis=0) poss = wrap_positions(poss, self.cell, self.pbc) D = squareform(pdist(poss)) screen, wrapped_poss = [], [] seen = set() for row, col in zip(*np.where(D < 0.15)): if row not in seen: screen.append(row) wrapped_poss.append(poss[row]) seen.update([row, col]) self.site_list = [] for i, j in enumerate(screen): s = sl[j] s['position'] = wrapped_poss[i] s['indices'] = tuple(sorted(np.mod(s['indices'], len(self.atoms)))) self.site_list.append(s)
[docs] def delaunay_triangulation(self, cutoff=5.): """Find all bridge and hollow site (3-fold and 4-fold) positions based on Delaunay triangulation of the surface atoms. Parameters ---------- cutoff : float, default 5. Radius of maximum atomic bond distance to consider. """ ext_index, ext_coords, _ = expand_cell(self.ref_atoms, cutoff) extended_top = np.where(np.in1d(ext_index, self.surf_ids))[0] ext_surf_coords = ext_coords[extended_top] # meansurfz = np.average(self.positions[self.surf_ids][:,2], 0) # surf_screen = np.where(abs(ext_surf_coords[:,2] - meansurfz) < 5.) # ext_surf_coords = ext_surf_coords[surf_screen] # Delaunay triangulation dt = scipy.spatial.Delaunay(ext_surf_coords[:,:2]) neighbors = dt.neighbors simplices = dt.simplices # Site type identification (partly borrowed from Catkit) bridge_positions, fold3_positions, fold4_positions = [], [], [] bridge_points, fold3_points, fold4_points = [], [], [] for i, corners in enumerate(simplices): cir = scipy.linalg.circulant(corners) edges = cir[:,1:] # Inner angle of each triangle corner vec = ext_surf_coords[edges.T] - ext_surf_coords[corners] uvec = vec.T / np.linalg.norm(vec, axis=2).T angles = np.sum(uvec.T[0] * uvec.T[1], axis=1) # Angle types right = np.isclose(angles, 0) obtuse = (angles < -1e-5) rh_corner = corners[right] edge_neighbors = neighbors[i] bridge = np.sum(ext_surf_coords[edges], axis=1) / 2.0 # Looping through corners allows for elimination of # redundant points, identification of 4-fold hollows, # and collection of bridge neighbors. for j, c in enumerate(corners): edge = sorted(edges[j]) if edge in bridge_points: continue # Get the bridge neighbors (for adsorption vector) neighbor_simplex = simplices[edge_neighbors[j]] oc = list(set(neighbor_simplex) - set(edge))[0] # Right angles potentially indicate 4-fold hollow potential_hollow = edge + sorted([c, oc]) if c in rh_corner: if potential_hollow in fold4_points: continue # Assumption: If not 4-fold, this suggests # no hollow OR bridge site is present. ovec = ext_surf_coords[edge] - ext_surf_coords[oc] ouvec = ovec / np.linalg.norm(ovec) oangle = np.dot(*ouvec) oright = np.isclose(oangle, 0) if oright: fold4_points.append(potential_hollow) fold4_positions.append(bridge[j]) else: bridge_points.append(edge) bridge_positions.append(bridge[j]) if not right.any() and not obtuse.any(): fold3_position = np.average(ext_surf_coords[corners], axis=0) fold3_points += corners.tolist() fold3_positions.append(fold3_position) return bridge_positions, fold3_positions, fold4_positions
[docs] def get_site(self, indices=None): """Get information of an adsorption site. Parameters ---------- indices : list or tuple, default None The indices of the atoms that contribute to the site. Return a random site if the indices are not provided. """ if indices is None: st = random.choice(self.site_list) else: indices = indices if is_list_or_tuple(indices) else [indices] indices = tuple(sorted(indices)) st = next((s for s in self.site_list if s['indices'] == indices), None) return st
[docs] def get_sites(self, site=None, morphology=None, composition=None, subsurf_element=None, return_site_indices=False): """Get information of all sites. Parameters ---------- site : str, default None Only return sites that belongs to this site type. morphology : str, default None Only return sites with this local surface morphology. composition : str, default None Only return sites that have this composition. subsurf_element : str, default None Only return sites that have this subsurface element. return_site_indices: bool, default False Whether to return the indices of each unique site (in the site list). """ if return_site_indices: all_sites = deepcopy(self.site_list) for i, st in enumerate(all_sites): st['site_index'] = i else: all_sites = self.site_list if site is not None: all_sites = [s for s in all_sites if s['site'] == site] if morphology is not None: all_sites = [s for s in all_sites if s['morphology'] == morphology] if composition is not None: if '-' in composition or len(list(Formula(composition))) == 6: scomp = composition else: comp = re.findall('[A-Z][^A-Z]*', composition) if len(comp) != 4: scomp = ''.join(sorted(comp)) else: if comp[0] != comp[2]: scomp = ''.join(sorted(comp)) else: if comp[0] > comp[1]: scomp = comp[1]+comp[0]+comp[3]+comp[2] else: scomp = ''.join(comp) all_sites = [s for s in all_sites if s['composition'] == scomp] if subsurf_element is not None: all_sites = [s for s in all_sites if s['subsurf_element'] == subsurf_element] if return_site_indices: all_stids = [] for st in all_sites: all_stids.append(st['site_index']) del st['site_index'] return all_stids else: return all_sites
[docs] def get_unique_sites(self, unique_composition=False, unique_subsurf=False, return_signatures=False, return_site_indices=False, about=None, site_list=None): """Get all symmetry-inequivalent adsorption sites (one site for each type). Parameters ---------- unique_composition : bool, default False Take site composition into consideration when checking uniqueness. unique_subsurf : bool, default False Take subsurface element into consideration when checking uniqueness. return_signatures : bool, default False Whether to return the unique signatures of the sites instead. return_site_indices: bool, default False Whether to return the indices of each unique site (in the site list). about: numpy.array, default None If specified, returns unique sites closest to this reference position. """ if site_list is None: sl = self.site_list else: sl = site_list key_list = ['site', 'morphology'] if unique_composition: if not self.composition_effect: raise ValueError('the site list does not include ' + 'information of composition') key_list.append('composition') if unique_subsurf: key_list.append('subsurf_element') else: if unique_subsurf: raise ValueError('to include the subsurface element, ' + 'unique_composition also need to be set to True') if return_signatures: sklist = sorted([[s[k] for k in key_list] for s in sl]) return sorted(list(sklist for sklist, _ in groupby(sklist))) else: seen_tuple = [] unq_sites = [] if about is not None: sl = sorted(sl, key=lambda x: get_mic(x['position'], about, self.cell, return_squared_distance=True)) for i, s in enumerate(sl): sig = tuple(s[k] for k in key_list) if sig not in seen_tuple: seen_tuple.append(sig) if return_site_indices: s = i unq_sites.append(s) return unq_sites
def get_labels(self): # Assign labels for st in self.site_list: if self.composition_effect: signature = [st['site'], st['morphology'], st['composition']] else: signature = [st['site'], st['morphology']] stlab = self.label_dict['|'.join(signature)] st['label'] = stlab def new_site(self): return {'site': None, 'surface': None, 'morphology': None, 'position': None, 'normal': None, 'indices': None, 'composition': None, 'subsurf_index': None, 'subsurf_element': None, 'label': None}
[docs] def mapping(self, atoms): """Map the slab into a surrogate reference slab for code versatility.""" if isinstance(self.surface, CustomSurface): # Sort atom indices if they deviate from the reference ref_atoms = self.surface.ref_atoms.copy() centered_atoms = ref_atoms.copy() centered_atoms.center(vacuum=5., axis=2) tmp_atoms = atoms.copy() tmp_atoms.center(vacuum=5., axis=2) centered_atoms, sorted_ids = sorted_by_ref_atoms(centered_atoms, tmp_atoms, mic=True, return_indices=True) ref_atoms = ref_atoms[sorted_ids] else: ref_atoms = atoms.copy() pm = self.surrogate_metal if pm is None: common_metal = Counter(self.atoms.symbols).most_common(1)[0][0] if common_metal in ['Ni', 'Cu', 'Pd', 'Ag', 'Pt', 'Au']: pm = common_metal area = np.linalg.norm(np.cross(atoms.cell[0], atoms.cell[1])) if self.surface in ['fcc100','fcc110','fcc311','fcc221','fcc331','fcc322', 'fcc332','bcc210','bcc211']: if area < 50.: ref_symbol = 'Cu' if pm is None else pm else: ref_symbol = 'Pt' if pm is None else pm elif self.surface in ['fcc111','fcc211','bcc111','hcp0001','hcp10m10h', 'hcp10m12']: ref_symbol = 'Cu' if pm is None else pm elif self.surface in ['bcc100','bcc110','bcc310','hcp10m10t','hcp10m11']: if area < 50.: ref_symbol = 'Cu' if pm is None else pm else: ref_symbol = 'Au' if pm is None else pm else: raise ValueError('surface {} is not implemented. '.format(self.surface) + 'Please use acat.settings.CustomSurface instead') for a in ref_atoms: a.symbol = ref_symbol try: from asap3 import EMT except: from ase.calculators.emt import EMT try: ref_atoms.calc = EMT() if self.optimize_surrogate_cell: ecf = ExpCellFilter(ref_atoms) opt = BFGS(ecf, logfile=None) else: opt = BFGS(ref_atoms, logfile=None) opt.run(fmax=0.1) except: pass centered_atoms = ref_atoms.copy() centered_atoms.center(vacuum=5., axis=2) midz = centered_atoms.cell[2][2] / 2 top_surf_indices = [] for i, a in enumerate(centered_atoms): if abs(a.position[2] - midz) < 3.: top_surf_indices = np.asarray([], dtype=int) break elif a.position[2] > midz: top_surf_indices.append(i) ref_atoms.positions[top_surf_indices] -= ref_atoms.cell[2] ref_atoms.center(vacuum=5., axis=2) ref_atoms.calc = None if self.optimize_surrogate_cell or isinstance(self.surface, CustomSurface): delta_positions = atoms.get_scaled_positions() -\ ref_atoms.get_scaled_positions() for i in range(len(delta_positions)): delta_frac = delta_positions[i] for dim in [0, 1]: if delta_frac[dim] < -0.5: delta_positions[i][dim] += 1. elif delta_frac[dim] > 0.5: delta_positions[i][dim] -= 1. else: delta_positions = atoms.positions - ref_atoms.positions return ref_atoms, delta_positions
def make_neighbor_list(self, neighbor_number=1): self.nblist = neighbor_shell_list(self.ref_atoms, self.tol, neighbor_number, mic=True)
[docs] def get_connectivity(self): """Get the adjacency matrix.""" return get_adj_matrix(self.nblist)
[docs] def get_termination(self, side='top'): """Return the indices of surface and subsurface atoms. This function relies on coordination number and the connectivity of the atoms. The top surface termination is singled out by graph adjacency using networkx. Parameters ---------- side : string, default 'top' The side of the surface termination ('top' or 'bottom'). """ assert side in ['top', 'bottom'] if isinstance(self.surface, CustomSurface): planes = get_layers(self.ref_atoms, (0,0,1), tolerance=self.surface.tol)[0] n_layers = self.surface.n_layers n_morphs = self.surface.n_morphs if side == 'top': surf = np.argwhere(planes>np.max(planes)-n_morphs).flatten() subsurf = np.argwhere((planes>np.max(planes)-2*n_morphs)& (planes<=np.max(planes)-n_morphs)).flatten() elif side == 'bottom': surf = np.argwhere(planes<np.min(planes)+n_morphs).flatten() subsurf = np.argwhere((planes<2*n_morphs)&(planes>=n_morphs)).flatten() else: cm = self.adj_matrix.copy() np.fill_diagonal(cm, 0) indices = list(range(len(self.ref_atoms))) coord = np.count_nonzero(cm, axis=1) both_surf = [] bulk = [] max_coord = np.max(coord) if self.surface == 'bcc210': max_coord -= 1 for i, c in enumerate(coord): a_s = indices[i] if c >= max_coord: bulk.append(a_s) else: both_surf.append(a_s) surfcm = cm.copy() surfcm[bulk] = 0 surfcm[:,bulk] = 0 if len(both_surf) == 2: components = [[both_surf[0]], [both_surf[1]]] else: # Use networkx to separate top layer and bottom layer rows, cols = np.where(surfcm == 1) edges = zip(rows.tolist(), cols.tolist()) G = nx.Graph() G.add_edges_from(edges) components = nx.connected_components(G) if side == 'top': surf = list(max(components, key=lambda x: np.mean( self.ref_atoms.positions[list(x),2]))) elif side == 'bottom': surf = list(min(components, key=lambda x: np.mean( self.ref_atoms.positions[list(x),2]))) subsurf = [] for a_b in bulk: for a_t in surf: if cm[a_t, a_b] > 0: subsurf.append(a_b) subsurf = list(np.unique(subsurf)) return sorted(surf), sorted(subsurf)
def get_two_vectors(self, indices): p1 = self.positions[indices[1]] p2 = self.positions[indices[2]] vec1 = get_mic(p1, self.positions[indices[0]], self.cell) vec2 = get_mic(p2, self.positions[indices[0]], self.cell) return vec1, vec2
[docs] def get_surface_normal(self, indices): """Get the surface normal vector of the plane from the indices of 3 atoms that forms to that plane. Parameters ---------- indices : list of tuple The indices of the atoms that forms the plane. """ vec1, vec2 = self.get_two_vectors(indices) n = np.cross(vec1, vec2) l = math.sqrt(n @ n.conj()) n /= l if n[2] < 0: n *= -1 return n
[docs] def get_graph(self, return_adj_matrix=False): """Get the graph representation of the nanoparticle. Parameters ---------- return_adj_matrix : bool, default False Whether to return adjacency matrix instead of the networkx.Graph object. """ cm = self.adj_matrix if return_adj_matrix: return cm G = nx.Graph() symbols = self.symbols G.add_nodes_from([(i, {'symbol': symbols[i]}) for i in range(len(symbols))]) rows, cols = np.where(cm == 1) edges = zip(rows.tolist(), cols.tolist()) G.add_edges_from(edges) return G
[docs] def get_neighbor_site_list(self, dx=0.1, neighbor_number=1, radius=None, span=True, unique=False, unique_composition=False, unique_subsurf=False): """Returns the site_list index of all neighbor shell sites for each site. Parameters ---------- dx : float, default 0.1 Buffer to calculate nearest neighbor site pairs. neighbor_number : int, default 1 Neighbor shell number. radius : float, default None The radius that defines site a is a neighbor of site b. This serves as an alternative to neighbor_number. span : bool, default True Whether to include all neighbors sites spanned within the shell. unique : bool, default False Whether to discard duplicate types of neighbor sites. """ sl = self.site_list poss = np.asarray([s['position'] for s in sl]) statoms = Atoms('X{}'.format(len(sl)), positions=poss, cell=self.cell, pbc=self.pbc) if radius is not None: nbslist = neighbor_shell_list(statoms, dx, neighbor_number=1, mic=True, radius=radius, span=span) if unique: for i in range(len(nbslist)): nbsids = nbslist[i] tmpids = self.get_unique_sites(unique_composition, unique_subsurf, return_site_indices=True, site_list=[sl[ni] for ni in nbsids]) nbsids = [nbsids[j] for j in tmpids] nbslist[i] = nbsids else: D = statoms.get_all_distances(mic=True) nbslist = {} for i, row in enumerate(D): argsort = np.argsort(row) row = row[argsort] res = np.split(argsort, np.where(np.abs(np.diff(row)) > dx)[0] + 1) if span: nbsids = sorted([n for nn in range(1, neighbor_number + 1) for n in list(res[nn])]) else: nbsids = sorted(res[neighbor_number]) if unique: tmpids = self.get_unique_sites(unique_composition, unique_subsurf, return_site_indices=True, site_list=[sl[ni] for ni in nbsids]) nbsids = [nbsids[j] for j in tmpids] nbslist[i] = nbsids return nbslist
[docs] def update(self, atoms, full_update=False): """Update the position of each adsorption site given an updated atoms object. Please only use this when the indexing of the atoms object is preserved. Useful for updating adsorption sites e.g. after geometry optimization. Parameters ---------- atoms : ase.Atoms object The updated atoms object. full_update : bool, default False Whether to update the normal vectors and composition as well. Useful when the composition of the alloy surface is not fixed. """ sl = self.site_list if full_update: nsas = SlabAdsorptionSites(atoms, surface=self.surface, allow_6fold=self.allow_6fold, composition_effect=self.composition_effect, both_sides=self.both_sides, ignore_sites=self.ignore_sites, label_sites=self.label_sites, surrogate_metal=self.surrogate_metal, optimize_surrogate_cell=self.optimize_surrogate_cell, tol=self.tol) sl = nsas.site_list else: new_slab = atoms[self.metal_ids] dvecs, _ = find_mic(new_slab.positions - self.positions, self.cell, self.pbc) nsl = deepcopy(sl) for st in nsl: si = list(st['indices']) idd = {idx: i for i, idx in enumerate(self.metal_ids)} tmpsi = [idd[k] for k in si] st['position'] += np.average(dvecs[tmpsi], 0) sl = nsl
[docs]def get_adsorption_site(atoms, indices, surface=None, return_index=False, **kwargs): """A function that returns the information of a site given the indices of the atoms that contribute to the site. The function is generalized for both periodic and non-periodic systems (distinguished by atoms.pbc). Parameters ---------- atoms : ase.Atoms object Accept any ase.Atoms object. No need to be built-in. indices : list or tuple The indices of the atoms that contribute to the site. surface : str, default None The surface type (crystal structure + Miller indices) Only required for periodic surface slabs. A user-defined customized surface object can also be used, but note that the identified site types will only include 'ontop', 'bridge', '3fold' and '4fold'. return_index : bool, default False Whether to return the site index of the site list together with the site. Example ------- This is an example of getting the site information of the (24, 29, 31) 3-fold hollow site on a fcc110 surface: .. code-block:: python >>> from acat.adsorption_sites import get_adsorption_site >>> from ase.build import fcc110 >>> atoms = fcc110('Cu', (2, 2, 8), vacuum=5.) >>> for atom in atoms: ... if atom.index % 2 == 0: ... atom.symbol = 'Au' >>> atoms.center(axis=2) >>> site = get_adsorption_site(atoms, (24, 29, 31), ... surface='fcc110', ... surrogate_metal='Cu') >>> print(site) Output: .. code-block:: python {'site': 'fcc', 'surface': 'fcc110', 'morphology': 'sc-tc-h', 'position': array([ 3.91083333, 1.91449161, 13.5088516 ]), 'normal': array([-0.57735027, 0. , 0.81649658]), 'indices': (24, 29, 31), 'composition': 'AuCuCu', 'subsurf_index': None, 'subsurf_element': None, 'label': None} """ indices = indices if is_list_or_tuple(indices) else [indices] indices = tuple(sorted(indices)) if True not in atoms.pbc: sas = ClusterAdsorptionSites(atoms, allow_6fold=True, composition_effect=True, **kwargs) else: sas = SlabAdsorptionSites(atoms, surface, allow_6fold=True, composition_effect=True, **kwargs) site_list = sas.site_list sti, site = next(((i, s) for i, s in enumerate(site_list) if s['indices'] == indices), None) if return_index: return sti, site return site
[docs]def enumerate_adsorption_sites(atoms, surface=None, morphology=None, **kwargs): """A function that enumerates all adsorption sites of the input catalyst structure. The function is generalized for both periodic and non-periodic systems (distinguished by atoms.pbc). Parameters ---------- atoms : ase.Atoms object Accept any ase.Atoms object. No need to be built-in. surface : str, default None The surface type (crystal structure + Miller indices). If the structure is a periodic surface slab, this is required. If the structure is a nanoparticle, the function enumerates only the sites on the specified surface. For periodic slabs, a user-defined customized surface object can also be used, but note that the identified site types will only include 'ontop', 'bridge', '3fold' and '4fold'. morphology : str, default None The function enumerates only the sites of the specified local surface morphology. Only available for surface slabs. Example ------- This is an example of enumerating all sites on the fcc100 surfaces of a Marks decahedral nanoparticle: .. code-block:: python >>> from acat.adsorption_sites import enumerate_adsorption_sites >>> from ase.cluster import Decahedron >>> atoms = Decahedron('Pd', p=3, q=2, r=1) >>> for atom in atoms: ... if atom.index % 2 == 0: ... atom.symbol = 'Ag' >>> atoms.center(vacuum=5.) >>> sites = enumerate_adsorption_sites(atoms, surface='fcc100', ... composition_effect=True, ... surrogate_metal='Pd') >>> print(sites[0]) Output: .. code-block:: python {'site': '4fold', 'surface': 'fcc100', 'position': array([18.86064518, 18.12221949, 11.87661345]), 'normal': array([ 0.58778525, 0.80901699, -0. ]), 'indices': (116, 117, 118, 119), 'composition': 'AgAgPdPd', 'subsurf_index': 75, 'subsurf_element': 'Pd', 'label': None} """ if True not in atoms.pbc: cas = ClusterAdsorptionSites(atoms, **kwargs) all_sites = cas.site_list if surface is not None: all_sites = [s for s in all_sites if s['surface'] == surface] else: sas = SlabAdsorptionSites(atoms, surface, **kwargs) all_sites = sas.site_list if morphology is not None: all_sites = [s for s in all_sites if s['morphology'] == morphology] return all_sites