import warnings
import numpy as np
from ase import Atom, Atoms
from ase.data import atomic_numbers, covalent_radii
from ase.neighborlist import NeighborList, natural_cutoffs
from scipy.optimize import minimize
from sklearn.cluster import MeanShift
[docs]class Adsorbate(Atoms):
"""
This is an adsorbate class which represents an adsorbate
"""
[docs] def __init__(self, symbols=None, positions=None, numbers=None, tags=None, momenta=None, masses=None, magmoms=None,
charges=None, scaled_positions=None, cell=None, pbc=None, celldisp=None, constraint=None,
calculator=None, info=None, velocities=None, host_zeotype=None, name='', description=''):
super().__init__(symbols, positions, numbers, tags, momenta, masses, magmoms, charges, scaled_positions, cell,
pbc, celldisp, constraint, calculator, info, velocities)
assert '_' not in description, 'cannot add _ to description'
if isinstance(symbols, Adsorbate):
if host_zeotype is None:
host_zeotype = symbols.host_zeotype
if description == '':
description = symbols.description
if name == '':
name = symbols.name
self.host_zeotype = host_zeotype
self.description = description
self.name = name
[docs] def min_distance(self, ads_position) -> float:
"""minimum distance from atom in a host to adsorbate :param
ads_position: np.array x,y,z of adsorbate position :return: float,
minimum distance between an atom in host MAZE-sim and adsorbate position
Args:
ads_position:
"""
assert self.host_zeotype is not None, "Cannot find min distance when host MAZE-sim is None"
dummy_atom = Atom('H', position=ads_position)
dummy_host = self.host_zeotype + dummy_atom
min_distance = min(dummy_host.get_distances(-1, [i for i in range(len(self.host_zeotype))], mic=True))
return min_distance
[docs] def avg_distance(self, ads_position):
"""average distance from position to all host atoms :param ads_position:
np.array x,y,z of adsorbate position :return: float, average distance of
host MAZE-sim atoms to adsorbate position
Args:
ads_position:
"""
assert self.host_zeotype is not None, "Cannot find average distance when host MAZE-sim is None"
dummy_atom = Atom('H', position=ads_position)
dummy_host = self.host_zeotype + dummy_atom
avg_distance = np.average(dummy_host.get_distances(-1, [i for i in range(len(self.host_zeotype))], mic=True))
return avg_distance
[docs] def find_void(self, void_position_guess):
"""finds nearest empty region in host MAZE-sim :param
void_position_guess: An initial guess for the center of the void as a
np.array with x,y,z position :return: np.array position of center of
empty void
Args:
void_position_guess:
"""
# TODO: Find a way to speed this up
assert self.host_zeotype is not None, "Cannot find void position when host MAZE-sim is None"
fun_to_min = lambda pos: -1 * self.min_distance(pos) # 1 param function for scipy.minimize
ans = minimize(fun_to_min, void_position_guess)
return ans.x
[docs] @staticmethod
def sphere_sample(radius, num_pts=500):
"""generates random positions on the surface of a sphere of certain
radius :param radius: radius of sphere surface to sample :param num_pts:
number of points to try :return: list of x,y,z positions on surface of
sphere
Args:
radius:
num_pts:
"""
position_list = []
for _ in range(int(num_pts)):
# see https://stackoverflow.com/questions/33976911/
# generate-a-random-sample-of-points-distributed-on-the-surface-of-a-unit-sphere/33977530#33977530
# for discussion on this algorithm
vec = np.random.normal(0, 1, 3) # select three random points (if normal dist no skip needed)
vec /= np.linalg.norm(vec) # normalize vector
vec *= radius # lengthen vector to desired radius
position_list.append(list(vec))
return position_list
[docs] def get_viable_positions(self, index, radius, cutoff, num_pts=None):
"""finds positions near host atom far enough from other framework atoms.
:param index: index of host atom at center :param radius: radius around
host atom to tests points :param cutoff: minimum distance from other host
atoms allowed for tests points :param num_pts: number of points to try
:return: list. positions of points which meet cutoff criteria
Args:
index:
radius:
cutoff:
num_pts:
"""
assert (radius > cutoff), "radius larger than cutoff distance"
assert self.host_zeotype is not None, "host MAZE-sim cannot be none"
guess_positions = self.sphere_sample(radius, num_pts)
host_pos = self.host_zeotype.get_positions()[index]
viable_positions = []
for pos in guess_positions:
dist = self.min_distance(pos + host_pos)
if dist > cutoff:
viable_positions.append(pos + host_pos)
return viable_positions
[docs] def get_viable_pos_cluster_centers(self, index, radius, cutoff, num_pts=None):
"""finds positions near host atom far enough from other framework atoms,
clusters these viable positions and returns the centers of these
clusters. If number of points is too small will return error :param
index: index of host atom at center :param radius: radius around host
atom to tests points :param cutoff: minimum distance from other host
atoms allowed for tests points :param num_pts: number of points to try
:return: list. center positions of clusters of points which meet
criteria
Args:
index:
radius:
cutoff:
num_pts:
"""
viable_pos = self.get_viable_positions(index, radius, cutoff, num_pts)
ms = MeanShift(bin_seeding=True)
ms.fit(np.array(viable_pos))
cluster_centers = ms.cluster_centers_
return cluster_centers
[docs] def find_best_place(self, index, radius, cutoff, num_pts=500):
"""picks the best location to place an adsorbate around the host atom
:param index: index of host atom at center :param radius: radius around
host atom to tests points :param cutoff: minimum distance from other host
atoms allowed for tests points :param num_pts: number of points to try
:return: array. x,y,z position of best location
Args:
index:
radius:
cutoff:
num_pts:
"""
best_positions = self.get_viable_pos_cluster_centers(index, radius, cutoff, num_pts)
best_avg_dist = 0
for pos in best_positions:
void = self.find_void(pos)
void_avg_dist = self.avg_distance(void) # average dist at void nearest to pos
best_pos = None
if void_avg_dist > best_avg_dist:
best_avg_dist = void_avg_dist # selects pos with largest nearby void
best_pos = pos
if best_pos is None:
assert False, 'No good positions at specified index'
return best_pos
[docs] def pick_donor(self):
"""finds atom in adsorbate most likely to bind to metal Heuristic: N > O
> P > S > X > C > H :param ads: adsorbate atoms object :return: index of
donor atom in adsorbate
"""
ads_symbols = {}
for atom in self:
if atom.symbol not in ads_symbols:
ads_symbols[atom.symbol] = atom.index # add first occurrence of atom symbol to dict
donor_symbols = ['N', 'O', 'P', 'S', 'F', 'Cl', 'Br', 'I', 'C', 'H']
for donor_symbol in donor_symbols:
if donor_symbol in ads_symbols:
ads_index = ads_symbols[donor_symbol] # picks first index of atom with specified symbol
return ads_index
else:
raise ValueError("Cannot find viable donor atom in adsorbate")
[docs] def pick_host_atom(self) -> int:
"""picks element with largest atomic number higher than 14 (Si) if one exists. If none exist,
pick Al as a donor atom.
Returns:
index of host atom
"""
sym_index_map, count = self.host_zeotype.count_elements()
max_atomic_num_sym = max(sym_index_map.keys(), key=lambda x: atomic_numbers[x])
if atomic_numbers[max_atomic_num_sym] > atomic_numbers['Si']:
return sym_index_map[max_atomic_num_sym][0]
else:
try:
return sym_index_map['Al'][0]
except KeyError:
print("No Al in host MAZE-sim")
raise
[docs] def pick_ads_position(self, donor_ind, host_ind, radius=None, cutoff=None):
"""Finds a good position to add adsorbate to host :param donor_ind:
index of donor atom on adsorbate :param host_ind: index of host binding
site :param radius: distance between host atom and donor atom :param
cutoff: minimum distance donor atom must be from host atoms :return:
vector of best position to place adsorbate
Args:
donor_ind:
host_ind:
radius:
cutoff:
"""
donor_atom_symbol, host_atom_symbol = self[donor_ind].symbol, self.host_zeotype[host_ind].symbol
donor_atom_number, host_atom_number = atomic_numbers[donor_atom_symbol], atomic_numbers[host_atom_symbol]
donor_radius, host_radius = covalent_radii[donor_atom_number], covalent_radii[host_atom_number]
if radius is None:
radius = host_radius + donor_radius
if cutoff is None:
cutoff = host_radius
pos = self.find_best_place(host_ind, radius, cutoff)
return pos
[docs] def get_donor_vec(self, donor_index):
"""finds direction of lone pair electrons on an adsorbate donor atom
:return: vector direction of donor atoms
Args:
donor_index:
"""
nl = NeighborList(natural_cutoffs(self), self_interaction=False, bothways=True)
nl.update(self)
# gets neighbors of donor atom and adds the vectors from neighbor to donor
# for most donor atoms this is roughly in the proper binding direction
donor_neighbors = nl.get_neighbors(donor_index)[0] # neighbor's index
donor_vec = np.array([0, 0, 0])
for i in donor_neighbors:
a = self.get_distance(i, donor_index, vector=True)
donor_vec = donor_vec + a
if np.linalg.norm(donor_vec) == 0:
warnings.warn("donor vector with magnitude 0 found, providing default vector")
return np.array([1, 0, 0])
donor_vec = donor_vec / np.linalg.norm(donor_vec) # normalizes donor vec
return donor_vec
[docs] def position_ads(self, donor_ind=None, host_ind=None, pos=None) -> "Adsorbate":
"""Rotates and positions adsorbate according to specified parameters if
no parameters are provided, a reasonable position is found :param pos:
vector, the position to place adsorbate's donor atom :param host_ind:
integer, index of site in host which adsorbate will be bound :param
donor_ind: integer, index of donor atom on adsorbate :return: atoms
object with adsorbate in host
Args:
donor_ind:
host_ind:
pos:
"""
if donor_ind is None:
donor_ind = self.pick_donor()
if host_ind is None:
host_ind = self.pick_host_atom()
if pos is None:
pos = self.pick_ads_position(donor_ind, host_ind)
dummy_host = self.host_zeotype + Atom('H', position=pos) # add dummy hydrogen atom to get distances to host atoms
vec = dummy_host.get_distance(-1, host_ind, mic=True, vector=True)
donor_vec = self.get_donor_vec(donor_ind) # get the direction of lone pairs on donor atom
new_self = self.__class__(self)
# rotate the ads into binding direction, move the ads to proper pos and combine
new_self.rotate(donor_vec, vec)
new_self.translate(pos - self.get_positions()[donor_ind])
return new_self