Difference in stoichiometry when using symmetrize=True vs False in SlabGenerator

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:

  1. Why does symmetrize=True sometimes alter the stoichiometry of the generated slab?
  2. Is this behavior expected for polar surfaces where the two terminations are inherently different?
  3. 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)?
  4. 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

The question is what property you are interested in calculating.
I guess the symmetrize keyword actually adds atoms in order to make the surfaces symmetric (at the cost of stoichiometry).

Please note for asymmetric surfaces it is not possible to define formation energies/surface energies within periodic boundary conditions. What you can do is to calculate energy of seperation/work of adhesion.

Thank you for the clarification.
Since our goal is to calculate surface energies, we would like to ensure that the generated slabs have symmetric surfaces while also preserving the stoichiometry of the original bulk structure.

Could you please advise on the best way to achieve both, meaning to generate stoichiometric slabs with symmetric terminations suitable for surface energy calculations?

Best regards,
Nadav

Hi, would you please be able to help me figure this out?