Hi everyone,
I’m using pymatgen.core.surface.SlabGenerator to create surface slabs from bulk structures. My workflow is as follows:
from pymatgen.io.vasp import Poscar
from pymatgen.core.surface import SlabGenerator
from pymatgen.analysis.structure_matcher import StructureMatcher
structure = Poscar.from_file("POSCAR").structure
structure.add_oxidation_state_by_element({"Al": 3, "N": -3})
slabgen = SlabGenerator(
initial_structure=structure,
miller_index=(1, 1, 1),
min_slab_size=15,
min_vacuum_size=15,
center_slab=True,
in_unit_planes=False,
primitive=False,
lll_reduce=True
)
slabs = slabgen.get_slabs(symmetrize=False)
When I set symmetrize=True, the generated slab sometimes has a different stoichiometry than the input bulk (for example, fewer N layers in AlN (111) or fewer C layers in Mo₂C (100)).
However, when I use symmetrize=False, the slabs retain bulk stoichiometry, but since they are asymmetric, the top and bottom terminations differ, making them polar and potentially unreliable for direct surface energy comparisons (where each termination should be defined separately).
My questions are:
- Why does
symmetrize=Truesometimes alter the stoichiometry of the generated slab? - Is this behavior expected for polar surfaces where the two terminations are inherently different?
- What’s the recommended best practice if I want to keep bulk stoichiometry and study surface energies for a specific termination (without mixing terminations in the same slab)?
- Would it be better to build asymmetric, stoichiometric slabs (
symmetrize=False) and then create mirrored “double slabs” to remove the dipole before DFT calculations?
Any clarification or recommended workflow for handling this situation would be greatly appreciated.
Thanks,
Nadav