#!/usr/bin/env python
# -*- coding: utf-8 -*-
from ..settings import (adsorbate_elements,
site_heights,
adsorbate_list,
adsorbate_molecule)
from ..adsorbate_coverage import (ClusterAdsorbateCoverage,
SlabAdsorbateCoverage,
enumerate_updated_sites)
from ..utilities import (custom_warning,
is_list_or_tuple,
get_close_atoms,
get_rodrigues_rotation_matrix,
get_angle_between,
get_rejection_between)
from ..labels import (get_cluster_signature_from_label,
get_slab_signature_from_label)
from ase.formula import Formula
from ase import Atoms, Atom
import numpy as np
import warnings
import re
warnings.formatwarning = custom_warning
[docs]def add_adsorbate(atoms, adsorbate, site=None, surface=None,
morphology=None, indices=None, height=None,
composition=None, orientation=None,
tilt_angle=0., subsurf_element=None,
all_sites=None, **kwargs):
"""A general function for adding one adsorbate to the surface.
Note that this function adds one adsorbate to a random site
that meets the specified condition regardless of it is already
occupied or not. 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.
adsorbate : str or ase.Atom object or ase.Atoms object
The adsorbate species to be added onto the surface.
site : str, default None
The site type that the adsorbate should be added to.
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 morphology type that the adsorbate should be added to.
Only available for surface slabs.
indices : list or tuple
The indices of the atoms that contribute to the site that
you want to add adsorbate to. This has the highest priority.
height : float, default None
The height of the added adsorbate from the surface.
Use the default settings if not specified.
composition : str, default None
The elemental of the site that should be added to.
orientation : list or numpy.array, default None
The vector that the multidentate adsorbate is aligned to.
tilt_angle: float, default 0.
Tilt the adsorbate with an angle (in degrees) relative to
the surface normal.
subsurf_element : str, default None
The subsurface element of the hcp or 4fold hollow site that
should be added to.
all_sites : list of dicts, default None
The list of all sites. Provide this to make the function
much faster. Useful when the function is called many times.
"""
composition_effect = any(v is not None for v in
[composition, subsurf_element])
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)
else:
scomp = None
if all_sites is None:
all_sites = enumerate_updated_sites(atoms, surface=surface,
morphology=morphology, occupied=False,
composition_effect=composition_effect,
label_occupied_sites=True, **kwargs)
if indices is not None:
if site is not None:
all_sites = [s for s in all_sites if s['site'] == site]
if scomp is not None:
all_sites = [s for s in all_sites if s['composition'] == scomp]
indices = indices if is_list_or_tuple(indices) else [indices]
indices = tuple(sorted(indices))
st = next((s for s in all_sites if s['indices'] == indices), None)
elif subsurf_element is None:
st = next((s for s in all_sites if s['site'] == site and
s['composition'] == scomp), None)
else:
st = next((s for s in all_sites if s['site'] == site and
s['composition'] == scomp and s['subsurf_element']
== subsurf_element), None)
if not st:
warnings.warn('No such site can be found')
else:
if height is None:
height = site_heights[st['site']]
add_adsorbate_to_site(atoms, adsorbate, st, height,
orientation, tilt_angle)
[docs]def add_adsorbate_to_site(atoms, adsorbate, site, height=None,
orientation=None, tilt_angle=0.):
"""The base function for adding one adsorbate to a site.
Site must include information of 'normal' and 'position'.
Useful for adding adsorbate to multiple sites or adding
multidentate adsorbates.
Parameters
----------
atoms : ase.Atoms object
Accept any ase.Atoms object. No need to be built-in.
adsorbate : str or ase.Atom object or ase.Atoms object
The adsorbate species to be added onto the surface.
site : dict
The site that the adsorbate should be added to.
Must contain information of the position and the
normal vector of the site.
height : float, default None
The height of the added adsorbate from the surface.
Use the default settings if not specified.
orientation : list or numpy.array, default None
The vector that the multidentate adsorbate is aligned to.
tilt_angle: float, default None
Tilt the adsorbate with an angle (in degrees) relative to
the surface normal.
"""
if height is None:
height = site_heights[site['site']]
# Make the correct position
normal = site['normal']
if np.isnan(np.sum(normal)):
warnings.warn('The normal vector is NaN, use [0., 0., 1.] instead.')
normal = np.array([0., 0., 1.])
pos = site['position'] + normal * height
# Convert the adsorbate to an Atoms object
if isinstance(adsorbate, Atoms):
ads = adsorbate
elif isinstance(adsorbate, Atom):
ads = Atoms([adsorbate])
# Or assume it is a string representing a molecule
else:
ads = adsorbate_molecule(adsorbate)
if not ads:
warnings.warn('Nothing is added.')
return
bondpos = ads[0].position
ads.translate(-bondpos)
z = -1. if adsorbate in ['CH','NH','OH','SH'] else 1.
ads.rotate(np.asarray([0., 0., z]) - bondpos, normal)
if tilt_angle > 0.:
pvec = np.cross(np.random.rand(3) - ads[0].position, normal)
ads.rotate(tilt_angle, pvec, center=ads[0].position)
if isinstance(adsorbate, str) and (adsorbate not in adsorbate_list):
# Always sort the indices the same order as the input symbol.
# This is a naive sorting which might cause H in wrong order.
# Please sort your own adsorbate atoms by reindexing as has
# been done in the adsorbate_molecule function in acat.settings.
symout = list(Formula(adsorbate))
symin = list(ads.symbols)
newids = []
for elt in symout:
idx = symin.index(elt)
newids.append(idx)
symin[idx] = None
ads = ads[newids]
if orientation is not None:
orientation = np.asarray(orientation)
oripos = next((a.position for a in ads[1:] if
a.symbol != 'H'), ads[1].position)
v1 = get_rejection_between(oripos - bondpos, normal)
v2 = get_rejection_between(orientation, normal)
theta = get_angle_between(v1, v2)
# Flip the sign of the angle if the result is not the closest
rm_p = get_rodrigues_rotation_matrix(axis=normal, angle=theta)
rm_n = get_rodrigues_rotation_matrix(axis=normal, angle=-theta)
npos_p, npos_n = rm_p @ oripos, rm_n @ oripos
nbpos_p = npos_p + pos - bondpos
nbpos_n = npos_n + pos - bondpos
d_p = np.linalg.norm(nbpos_p - pos - orientation)
d_n = np.linalg.norm(nbpos_n - pos - orientation)
if d_p <= d_n:
for a in ads:
a.position = rm_p @ a.position
else:
for a in ads:
a.position = rm_n @ a.position
ads.translate(pos - bondpos)
atoms += ads
if ads.get_chemical_formula() == 'H2':
shift = (atoms.positions[-2] - atoms.positions[-1]) / 2
atoms.positions[-2:,:] += shift
[docs]def add_adsorbate_to_label(atoms, adsorbate, label,
surface=None, height=None,
orientation=None,
tilt_angle=0.,
composition_effect=False,
all_sites=None, **kwargs):
"""Same as add_adsorbate function, except that the site type is
represented by a numerical label. 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.
adsorbate : str or ase.Atom object or ase.Atoms object
The adsorbate species to be added onto the surface.
label : int or str
The label of the site that the adsorbate should be added to.
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'.
height : float, default None
The height of the added adsorbate from the surface.
Use the default settings if not specified.
orientation : list or numpy.array, default None
The vector that the multidentate adsorbate is aligned to.
tilt_angle: float, default 0.
Tilt the adsorbate with an angle (in degrees) relative to
the surface normal.
composition_effect : bool, default False
Whether the label is defined in bimetallic labels or not.
all_sites : list of dicts, default None
The list of all sites. Provide this to make the function
much faster. Useful when the function is called many times.
"""
if composition_effect:
slab = atoms[[a.index for a in atoms if a.symbol
not in adsorbate_elements]]
metals = sorted(list(set(slab.symbols)))
else:
metals = None
if True in atoms.pbc:
signature = get_slab_signature_from_label(label, surface,
composition_effect,
metals)
else:
signature = get_cluster_signature_from_label(label,
composition_effect,
metals)
sigs = signature.split('|')
morphology, composition = None, None
if not composition_effect:
if True in atoms.pbc:
site, morphology = sigs[0], sigs[1]
else:
site, surface = sigs[0], sigs[1]
else:
if True in atoms.pbc:
site, morphology, composition = sigs[0], sigs[1], sigs[2]
else:
site, surface, composition = sigs[0], sigs[1], sigs[2]
add_adsorbate(atoms, adsorbate, site, surface,
morphology, height=height,
composition=composition,
orientation=orientation,
all_sites=all_sites, **kwargs)
[docs]def remove_adsorbate_from_site(atoms, site, remove_fragment=False):
"""The base function for removing one adsorbate from an
occupied site. The site must include information of
'adsorbate_indices' or 'fragment_indices'. Note that if
you want to remove adsorbates from multiple sites, call
this function multiple times will return the wrong result.
Please use remove_adsorbates_from_sites instead.
Parameters
----------
atoms : ase.Atoms object
Accept any ase.Atoms object. No need to be built-in.
site : dict
The site that the adsorbate should be removed from.
Must contain information of the adsorbate indices.
remove_fragment : bool, default False
Remove the fragment of a multidentate adsorbate instead
of the whole adsorbate.
"""
if site['occupied']:
if not remove_fragment:
si = list(site['adsorbate_indices'])
else:
si = list(site['fragment_indices'])
del atoms[si]
else:
warnings.warn('This site is not occupied.')
[docs]def remove_adsorbates_from_sites(atoms, sites, remove_fragments=False):
"""The base function for removing multiple adsorbates from
an occupied site. The sites must include information of
'adsorbate_indices' or 'fragment_indices'.
Parameters
----------
atoms : ase.Atoms object
Accept any ase.Atoms object. No need to be built-in.
sites : list of dicts
The site that the adsorbate should be removed from.
Must contain information of the adsorbate indices.
remove_fragments : bool, default False
Remove the fragment of a multidentate adsorbate instead
of the whole adsorbate.
"""
if not remove_fragments:
si = [i for s in sites if s['occupied'] for
i in s['adsorbate_indices']]
else:
si = [i for s in sites if s['occupied'] for
i in s['fragment_indices']]
del atoms[si]
[docs]def remove_adsorbates_too_close(atoms, adsorbate_coverage=None,
surface=None,
min_adsorbate_distance=0.5):
"""Find adsorbates that are too close, remove one set of them.
The function is intended to remove atoms that are unphysically
close. Please do not use a min_adsorbate_distace larger than 2.
The function is generalized for both periodic and non-periodic
systems (distinguished by atoms.pbc).
Parameters
----------
atoms : ase.Atoms object
The nanoparticle or surface slab onto which the adsorbates are
added. Accept any ase.Atoms object. No need to be built-in.
adsorbate_coverage : acat.adsorbate_coverage.ClusterAdsorbateCoverage object or acat.adsorbate_coverage.SlabAdsorbateCoverage object, default None
The built-in adsorbate coverage class.
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 only remove
adsorbates from 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'.
min_adsorbate_distance : float, default 0.
The minimum distance between two atoms that is not considered to
be to close. This distance has to be small.
"""
if adsorbate_coverage is not None:
sac = adsorbate_coverage
else:
if True not in atoms.pbc:
sac = ClusterAdsorbateCoverage(atoms)
else:
sac = SlabAdsorbateCoverage(atoms, surface)
dups = get_close_atoms(atoms, cutoff=min_adsorbate_distance,
mic=(True in atoms.pbc))
if dups.size == 0:
return
hsl = sac.hetero_site_list
# Make sure it's not the bond length within a fragment being too close
bond_rows, frag_id_list = [], []
for st in hsl:
if st['occupied']:
frag_ids = list(st['fragment_indices'])
frag_id_list.append(frag_ids)
w = np.where((dups == x).all() for x in frag_ids)[0]
if w:
bond_rows.append(w[0])
dups = dups[[i for i in range(len(dups)) if i not in bond_rows]]
del_ids = set(dups[:,0])
rm_ids = [i for lst in frag_id_list for i in lst if
del_ids.intersection(set(lst))]
rm_ids = list(set(rm_ids))
del atoms[rm_ids]