When generating slabs via SlabGenerator.get_slabs()
, you can specify a bonds
dictionary and it will generate all the slabs with # of bonds broken < max_broken_bonds. The issue I found with this is that the number of broken bonds seem inaccurate as far as I can tell and the reason for this seems to be the _get_c_ranges()
method.
This method iterates over all the sites and its nearest neighbors in the oriented unit cell and if a site and its NN belong to different layers (z-coordinates in this case), the fractional z-coordinates of the site and its neighbor get added to a set()
as a sorted tuple.
The problem as I see it stems from the fact that since sets only contain unique elements, this approach ignores all the other ‘bonds’ between two layers as long as 1 bond (tuple) is already added to the set.
I tested this with asymmetric Si(111) which has two distinct terminations and the number of bonds broken should be 2 and 6 looking at the structures in VESTA. Due to this issue however, the number of broken bonds for these two terminations are assigned as 0.5 for each.
If I change the set()
in SlabGenerator._get_c_ranges()
to an empty list instead and append to it, I get the expected 2 and 6 bonds broken as a result.
I was wondering if there’s something I am overlooking here or this is indeed an issue with the calculation of the number of bonds broken.