Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix for correct bonds_broken in slab generation #2105

Merged
merged 2 commits into from Apr 1, 2021
Merged

Fix for correct bonds_broken in slab generation #2105

merged 2 commits into from Apr 1, 2021

Conversation

fyalcin
Copy link
Contributor

@fyalcin fyalcin commented Mar 30, 2021

Summary

SlabGenerator._get_c_ranges() takes in a bonds dictionary and returns a set of tuples of fractional z-coordinates of neighboring (nearest) atoms if they have distinct z-coordinates. This set is then used in the SlabGenerator.get_slabs() method in order to calculate the total number of broken bonds. The issue with the current version is that the tuples are stored in a set() which can only have unique elements so all but one of the bonds between say z_coord1 and z_coord2 are ignored.

I tested this by simply generating Si(111) slabs, calculating the number of bonds broken reported by pymatgen, and comparing them by what I expect from looking at the slabs in VESTA.

Here's the code I used to generate the slabs;

from pymatgen.core.surface import SlabGenerator
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen import MPRester

struct = MPRester().get_structure_by_material_id('mp-149')

bulk_conv = SpacegroupAnalyzer(struct).get_conventional_standard_structure()

SG = SlabGenerator(initial_structure = bulk_conv,
                   miller_index = [1,1,1],
                   center_slab = True,
                   primitive = True,
                   in_unit_planes = True,
                   lll_reduce = True,
                   max_normal_search= 1,
                   min_slab_size = 7,
                   min_vacuum_size = 12)

slabs = SG.get_slabs(bonds = {('Si','Si'):2.60}, ftol = 0.1, tol = 0.1, max_broken_bonds = 20,
                              symmetrize = False, repair = False)

for slab in slabs:
    print(slab.energy)

This prints 0.5 for both terminations whereas looking at the slabs in VESTA, one expects 2 and 6.

By changing c_ranges from a set() to an empty list instead and appending to it, tuples for all the bonds are preserved and the above script prints 2 and 6.

For consistency, I also changed c_ranges to an empty list if bonds is None in get_slabs().

@shyuep
Copy link
Member

shyuep commented Mar 30, 2021

Thanks for submitting this PR. Do you mind adding a unittest for the bug? You just need to add it to pymatgen/core/tests/test_surface.py.

@fyalcin
Copy link
Contributor Author

fyalcin commented Mar 31, 2021

Just added a unit test. It is the first unit test I ever wrote so I hope I was able to pull it off :)

@shyuep
Copy link
Member

shyuep commented Apr 1, 2021

Looks good to me. Thanks!

@shyuep shyuep merged commit 2f3bee0 into materialsproject:master Apr 1, 2021
@fyalcin fyalcin deleted the bonds_broken_fix branch April 1, 2021 23:21
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants