AdsorbateSiteFinder list all possible adsorbate sites

Hello

I am a very recent user of pymatgen. I have a 4x4 graphene supercell and would like to list all possible adsorbates (O atom) in all positions ontop, bridge, and hollow.

I first create supercell:

lattice = Lattice.hexagonal(2.46,5)
G = Structure(lattice, ["C", "C"], [[1/3,2/3,1/2],[2/3,1/3,1/2]])
G.make_supercell([4,4,1])

Then list the possible adsorbate configurations in:

adsorbate = Molecule("O", [[0, 0, 0]])
asf = AdsorbateSiteFinder(G)
asf.find_adsorption_sites()

However this only lists 4 sites:

{'ontop': [array([-2.33486982e-16,  1.42028166e+00,  4.50000000e+00])],
 'bridge': [array([-3.075     ,  7.45647873,  4.5       ]),
  array([6.15      , 6.39126748, 4.5       ]),
  array([-3.075     ,  6.03619706,  4.5       ])],
 'hollow': [],
 'all': [array([-2.33486982e-16,  1.42028166e+00,  4.50000000e+00]),
  array([-3.075     ,  7.45647873,  4.5       ]),
  array([6.15      , 6.39126748, 4.5       ]),
  array([-3.075     ,  6.03619706,  4.5       ])]}

Is there a way to list all possible sites?

Hope this makes sense
Hud

hi @hud

You should be able to set symm_reduce to 0, which will cause it not remove sites based on symmetry.

1 Like

hi @shyamd

adSites = asf.find_adsorption_sites(symm_reduce=0) worked.

len(adSites['ontop']) gives 32 atoms as expected in a 4X4 graphene supercell. However, len(adSites['bridge']) gives 158 positions while len(adSites['hollow']) gives 0 atoms.

The conditions for finding the sites are by default height=0.9 and distance=2.0. Could you advise what can be tweaked to improve the findings?

Hi @hud

adSites = asf.find_adsorption_sites(symm_reduce=0, no_obtuse_hollow=False) should return the hollow sites as well. no_obtuse_hollow is True by default, which means that hollow sites will not be assigned inside obtuse triangles, and for the graphene structure it looks like they are all obtuse.

2 Likes

Just to add a bit more, the pymatgen AdsorbateSiteFinder’s use of a triangulation algorithm to determine “hollow,” “bridge”, and “ontop” sites means that it doesn’t assign six-fold symmetric hollow sites (or four-fold, for that matter) with the appropriate designation according to human intuition, In other words, it can only consider and identify 3-atom groups for adsorbing ensembles, and doesn’t have any notion of a six-atom ensemble for the “hollow” site in graphene the way we would typically think about it. However, the six-fold site should be in the list of “bridge” sites because the delaunay triangulation of a 2-D hexagonal lattice should include one diagonal which crosses the constituent hexagons at its center. Long story short, the ASF is solid if you want to exhaustively consider symmetrically distinct adsorption sites, but it occasionally finds a few more than you would want to consider and sometimes doesn’t designate them the same way a human would.

Can you tell us more about your use case? I’m curious as to why you wouldn’t want to use the symmetry reduction.

As an aside, our designation criteria could probably be improved, I wonder if it might be better to assign site character in a post-processing step, rather than from the delaunay triangulation itself? I’ll spend some time thinking about it, but if anybody else has ideas would love to hear them.

In essence I want to optimize some property e.g. electrical conductivity based on the various configurations of adatoms using ML approach in atomate / AiiDA workflow loops. Ideally, the model should choose from minimal unique bridge, hollow or ontop adsorbate sites as how a human would but I guess manual work is still better here.

Would a post-processing step solve this problem?

Yeah, I think chances are for just a few calculations on graphene you won’t get much out of a fully-automated HT approach, but I’m glad you asked so we can spend some time thinking about how to support your use case.

For a single adsorbate, there should only be one each of “hollow”, “bridge,” and “ontop” right? Why do you need the symmetrically redundant sites?

As for postprocessing, my thinking is something along the lines of classifying site character as taking the nearest neighbor, then finding all nearest neighbors within that (bond distance) * (1 + epsilon), so basically all atoms within a certain percentage of the nearest neighbor, if there’s 1 atom, it’s “ontop”, if 2, “bridge”, if 3+, “n-fold hollow”, where is is the neighbor number.

from pymatgen import MPRester
with MPRester() as mpr:
    struct = mpr.get_structure_by_material_id("mp-48")
# Remove bottom layer
struct = Structure.from_sites(
    [site for site in struct
     if site.frac_coords[2] == 0.75]
)
asf = AdsorbateSiteFinder(struct)
sites = asf.find_adsorption_sites()

def classify_site(coord, slab):
    site_dists = slab.get_sites_in_sphere(coord, 7)
    min_dist = min([site_dist[1] for site_dist in site_dists])
    ensemble = slab.get_sites_in_sphere(coord, min_dist * 1.15)  # probably should tweak this
    ensemble_size = len(ensemble)
    if ensemble_size == 1:
        return "ontop"
    elif ensemble_size == 2:
        return "bridge"
    else:
        return "{}-fold hollow".format(ensemble_size)

for site in sites['all']:
    character = classify_site(site, struct)
    print(character, site)
ontop [0.         0.         7.85230447]
bridge [1.23364565 1.42449128 7.85230447]
6-fold hollow [6.16822871e-09 1.42449128e+00 7.85230447e+00]
3-fold hollow [-0.61682281  1.7806141   7.85230447]

There’s still a leftover “3-fold hollow” site that I think is from the “bridge” site that’s connecting the 1,3 carbons in the delaunay triangulation, but this does reassign the “bridge” in the center of the C-hexagon to “6-fold hollow”

I naively thought not reducing symmetry would give me all the sites - now I think it’s best to follow like what you did, find adsorbate in the unit cells and manually create an array of the sites according to my supercell.

This is not elegant code, but I figured I would make a structure adsorbate so I can easily extend it to a supercell.

newsites = {}
for site in sites['all']:
    newsites[classify_site(site, struct)] = site

ontop = Structure(struct.lattice,'O',[newsites['ontop']],coords_are_cartesian=True)
hollow = Structure(struct.lattice,'O',[newsites['6-fold hollow']],coords_are_cartesian=True)
bridge = Structure(struct.lattice,'O',[newsites['bridge']],coords_are_cartesian=True)

However, it’s a bit tedious to find the rest of the bridge sites:

hollow.translate_sites([0],[struct.lattice.lengths[0]/2,0,0],frac_coords=False,to_unit_cell=False)
bridge2 = hollow

bridge.translate_sites([0],[struct.lattice.lengths[0]/2,0,0],frac_coords=False,to_unit_cell=False)
bridge1 = bridge
bridge = Structure(struct.lattice,'O',[newsites['bridge']],coords_are_cartesian=True)

bridge.insert(1,'O',list(bridge1.cart_coords[0]),coords_are_cartesian=True)
bridge.insert(2,'O',list(bridge2.cart_coords[0]),coords_are_cartesian=True)
bridge
Structure Summary
Lattice
    abc : 2.46729128 2.467291282584372 7.80307263
 angles : 90.0 90.0 119.99999996535054
 volume : 41.13742744042275
      A : 2.46729128 0.0 0.0
      B : -1.23364564 2.13673693 0.0
      C : 0.0 0.0 7.80307263
PeriodicSite: O (0.6168, 0.3561, 7.8523) [0.3333, 0.1667, 1.0063]
PeriodicSite: O (1.8505, 0.3561, 7.8523) [0.8333, 0.1667, 1.0063]
PeriodicSite: O (1.2336, 1.4245, 7.8523) [0.8333, 0.6667, 1.0063]

Is there a way to append multiple structures together?

I’m still not understanding why you need additional of the bridge sites, as all of the bridge sites in a single graphene layer are symmetrically equivalent.

Sorry for the late reply.

If I had two bridge adsorbates (in red) like these, are they different?

test

Yep, this looks like multi-adsorbate configuration of a supercell, which are more complex to generate. If you want to get all of those, you’ll have to do a bit of hacking on top of what ASF does to find combinations of adsorbate positions. You might be able to remove symmetrically redundant adsorbate sets using structure matcher, but i’ve never used it for that purpose and don’t know how well it would work. There’s also some literature on this problem from Zach Ulissi’s group that uses bayesian approaches to sampling the adsorbate configuration space.