Dear Pymatgen Community,
I’m using pymatgen.analysis.adsorption.AdsorbateSiteFinder in Pymatgen 2019.4.11 to automatically enumerate adsorption sites on some metal chalcogenide crystal facets. My understanding is that you do this:
# Get Pymatgen structure from some place
struct = AseAtomsAdaptor.get_structure(atoms)
asf = AdsorbateSiteFinder(struct, height=2)
sites = asf.find_adsorption_sites(positions="bridge", put_inside=True, distance=1.75, no_obtuse_hollow=False)
You then get a dictionary in sites, with hollow and all keys. The hollow key is supposed to have all of the hollow sites, and the all key is supposed to just be the union of all the sites you requested (in case you specified multiple types of adsorption sites to be found, e.g., positions=["ontop", "hollow"]).
However, I’m finding that when I specify just one type of adsorption site, as done above, the all key is returning a list of adsorption sites that is longer than the hollow key. Additionally, the all subdictionary contains some sites that are not in the hollow subdictionary, and vice versa.
Here’s an example with the (100) facet of Rh3S4 (mp-29841, I believe). I’m looking for hollow sites on this surface:
# Find lists of adsorption sites
test_atoms = baresurfaces_primitive["Rh3S4"].copy()
asf = AdsorbateSiteFinder(AseAtomsAdaptor.get_structure(test_atoms), height=2)
sites_dict = asf.find_adsorption_sites(
positions=["hollow"], put_inside=True, distance=1.75, no_obtuse_hollow=False
)
sites_dict["hollow"] looks like this:
[array([ 1.64671456, 5.6360542 , 20.18552412]),
array([ 3.47302158, 0.13405058, 20.11960924]),
array([ 2.9894536 , 1.18865787, 20.02708126]),
array([ 4.71382708, 6.21485056, 20.57152668]),
array([ 6.77737094, 4.62549356, 20.76499493]),
array([ 2.93180922, 3.46316175, 20.69213354]),
array([ 3.69690319, 2.20623252, 20.47985424]),
array([ 9.89662907, 4.57747216, 20.19444695]),
array([ 1.6388581 , 3.50416913, 20.23043774]),
array([ 5.47892106, 3.89637141, 20.72833823]),
array([ 5.47892106, 2.53966717, 20.51787986]),
array([10.95784212, 2.44161235, 20.15575541]),
array([10.95784212, 6.00652148, 19.67824878]),
array([10.3723406 , 1.38303031, 20.16467823]),
array([ 0.57764506, 0.6988903 , 19.74993088])]
while sites_dict["all"] looks like this:
[array([ 2.88752006, 5.36536195, 20.69097852]),
array([ 1.06121304, 4.57747216, 20.19444695]),
array([ 1.6388581 , 3.50416913, 20.23043774]),
array([ 9.79469554, 1.75747234, 19.74100805]),
array([ 7.96838852, 1.18865787, 20.02708126]),
array([ 5.47892106, 0.50451786, 19.6123339 ]),
array([ 4.18047118, 4.62549356, 20.76499493]),
array([ 2.93180922, 3.46316175, 20.69213354]),
array([10.95784212, 5.46513697, 20.68915758]),
array([ 4.71382708, 6.21485056, 20.57152668]),
array([10.38019706, 1.3977513 , 20.13761027]),
array([ 7.26093892, 2.20623252, 20.47985424]),
array([ 5.47892106, 2.53966717, 20.51787986]),
array([ 5.47892106, 3.89637141, 20.72833823]),
array([ 9.79469554, 2.45633334, 20.12868744]),
array([ 9.31898401, 4.20303014, 20.61811713]),
array([ 3.47302158, 0.13405058, 20.11960924]),
array([ 1.64671456, 5.6360542 , 20.18552412]),
array([ 1.16314658, 5.99180049, 19.70531675]),
array([ 1.16314658, 0.33916926, 20.14653309]),
array([ 9.79469554, 3.87463641, 19.7231624 ]),
array([10.95784212, 2.44161235, 20.15575541]),
array([ 8.7334825 , 4.56275117, 20.22151491]),
array([ 1.64671456, 4.9371932 , 19.79784473]),
array([10.95784212, 6.00652148, 19.67824878]),
array([ 0.57764506, 0.6988903 , 19.74993088])]
If I sort the the sites by x, y, and then z coordinate, we see that the hollow sites and all sites dictionaries don’t cover the same set of sites:
The all dictionary of sites looks more thorough, but many of these sites are close together and would probably relax to the same position if I do DFT geometry relaxation. I believe something in the Delaunay triangulation code is adding more sites to the all list, because this doesn’t seem to happen when I request only ontop sites.
Is there a way to obtain a minimal set of unique hollow sites? And which list should I trust—hollow or all?
Thanks,
– Sam