Issue with calculation of broken bonds

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.

1 Like

Thanks for looking into this! I don’t know enough about how this particular method works in slab generation to say for sure, but this sounds reasonable. I’d recommend issuing a pull request to pymatgen with your proposed fix and the reasoning you’ve outlined here.

1 Like