# pymatgen structure manipulation gotchas: from_spacegroup, PBC wrapping, oxidation states, slabs, defects

Fix silent pymatgen structure bugs: from_spacegroup needs only the asymmetric unit, fractional coords are not auto-wrapped (1.0 vs 0.0 duplicates bonds), add_oxidation_state_by_guess returns all-zero on non-stoichiometric cells corrupting EwaldSummation, SlabGenerator min_slab_size is in Angstroms, and defect indices are cell-local.

Exact reference: {"kind":"skill_version","skill_id":"skl_yKw6ndolOZupR0cWbXSwYw","version_id":"skv_GNouELIB_PimvfyMJfT6wQ"}

Applicability: [{"constraint":"building crystal structures programmatically with pymatgen","technology":"pymatgen","version_scheme":"unknown"},{"constraint":"debugging duplicated atoms from Structure.from_spacegroup","technology":"pymatgen","version_scheme":"unknown"},{"constraint":"normalizing periodic boundary coordinates before analysis or IO","technology":"pymatgen","version_scheme":"unknown"},{"constraint":"oxidation-state decoration before Ewald summation or charge analysis","technology":"pymatgen","version_scheme":"unknown"},{"constraint":"generating slabs with SlabGenerator and interpreting empty results","technology":"pymatgen","version_scheme":"unknown"},{"constraint":"creating vacancy or defect supercells","technology":"pymatgen","version_scheme":"unknown"},{"constraint":"diagnosing wrong neighbor/bond counts from unwrapped coordinates","technology":"pymatgen","version_scheme":"unknown"}]

# pymatgen structure manipulation gotchas

Concrete pain points of building and transforming crystal structures programmatically in pymatgen: `Structure.from_spacegroup` vs manual construction, periodic-boundary wrapping with `PeriodicSite`, oxidation-state decoration before Ewald/charge analysis, and common errors generating slabs or defects.

## 1. `from_spacegroup` vs manual `Structure(lattice, species, coords)`

`Structure.from_spacegroup(sg, lattice, species, coords)` takes **only symmetrically distinct** sites and applies the spacegroup operations to generate the rest. The manual constructor takes **all** sites explicitly.

```python
from pymatgen.core import Structure, Lattice

# Only the asymmetric unit is needed here: Fm-3m expands these 2 sites
# into the full rock-salt cell.
li2o = Structure.from_spacegroup(
    "Fm-3m", Lattice.cubic(3.0),
    ["Li", "O"],
    [[0.25, 0.25, 0.25], [0, 0, 0]],
)

# Manual construction: every site is explicit.
cscl = Structure(Lattice.cubic(4.2), ["Cs", "Cl"],
                 [[0, 0, 0], [0.5, 0.5, 0.5]])
```

Gotchas:

- **Giving `from_spacegroup` the full site list duplicates atoms.** If you pass all equivalent positions, the generator creates overlapping copies and the structure is unphysical (pairs closer than 0.01 A apart). Pass only the asymmetric unit.
- **`sg` accepts an international-number int (e.g. 225 for Fm-3m) or a notation string** interpreted by `pymatgen.symmetry.groups.Spacegroup` (e.g. "R-3c", "Fm-3m"). The wrong setting/number silently yields a different cell, so always run a `SpacegroupAnalyzer(structure).get_space_group_number()` round-trip check when the origin of `sg` is ambiguous.
- **`tol` (default 1e-5) matters.** Coordinates that are distinct by less than `tol` are merged during expansion; coordinates from CIF round-off can collapse inequivalent sites. If site counts look wrong, compare with `SpacegroupAnalyzer(...).get_conventional_standard_structure()`.

## 2. Periodic boundary wrapping with `PeriodicSite`

`PeriodicSite` stores **fractional** coordinates. Distances are computed with periodic boundary conditions via `Lattice.get_distance_and_image(frac_coords1, frac_coords2)`.

Gotchas:

- **Coords are NOT wrapped into the unit cell by default.** `Structure(..., to_unit_cell=False)` leaves fractional coords exactly as given, e.g. 7.4060 stays 7.4060 (which under PBC is the same as 0.4060). Most analysis still works because distance functions use PBC, but visualizers and some neighbor routines assume coords in [0, 1). Pass `to_unit_cell=True` at construction, or rebuild with `Structure.from_sites(sites, to_unit_cell=True)`, when anything downstream indexes or bins coordinates.
- **Cartesian-instead-of-fractional is a classic silent bug.** The constructor's default `coords_are_cartesian=False` assumes fractional coords. If your numbers look like Angstroms (e.g. 7.4060 in a 7.406 A cell), set `coords_are_cartesian=True` or convert first; otherwise every distance is computed against a cell scaled by its own lattice constants.
- **Fractional coord exactly 1.0 vs 0.0 can duplicate bonds.** A site at `[1., 0.2, 0.2]` sits at the same position as `[0., 0.2, 0.2]`, but near-neighbor machinery historically floors fractional coords for the jimage, so the `1.` site is seen twice: once with jimage `(0,0,0)` and once with `(1,0,0)`. Normalize with `to_unit_cell=True` before running `CrystalNN.get_bonded_structure` or any analysis that counts neighbors.
- **Wrap after, not before, transformations.** Supercell scaling, translations, and slab shifting can push coords outside [0, 1); distance code is fine with that, but equality checks and JSON consumers may not be. When in doubt, wrap before serializing.

## 3. Oxidation-state decoration before Ewald/charge analysis

`EwaldSummation` requires every site to carry a `Species` with an oxidation state; the docs state the input structure "must have proper Species on all sites, i.e. Element with oxidation state" and point at `Structure.add_oxidation_state...`.

```python
from pymatgen.core import Structure, Lattice
from pymatgen.core.ewald import EwaldSummation

struct = Structure.from_spacegroup("Fm-3m", Lattice.cubic(4.2),
                                   ["Na", "Cl"],
                                   [[0, 0, 0], [0.5, 0.5, 0.5]])

# Decoration is in place and now returns self (post-2023), so it chains.
struct.add_oxidation_state_by_guess()
ewald = EwaldSummation(struct)
print(ewald.total_energy)
```

Gotchas:

- **Decorate on the exact cell you evaluate.** `add_oxidation_state_by_guess()` uses `Composition.oxi_state_guesses()` on the cell composition and takes the first guess. A non-stoichiometric defect or slab (e.g. a vacancy cell) can fail to find a charge-balanced combination, in which case the guess silently degrades to **0 for every element** (CuO2 guesses as Cu0+ O02-). That zeroed cell then feeds Ewald a nonsense charge distribution. Always inspect `struct[0].species_string` (e.g. "Na1+") after decoration, or explicitly pass `all_oxi_states=True` kwargs to see the ranked candidates.
- **Ewald on a neutral cell only.** `total_energy` excludes the charged-cell term, and the docs note this "is only important when the simulation cell is not charge balanced." If `struct.charge != 0`, the energy is not comparable across cells. For charged defects use `compute_sub_structure` or compare per-site energies (`get_site_energy`) instead.
- **Ewald recompute is expensive; reuse the matrix.** The docstring notes "This matrix can be used to do fast calculations of Ewald sums after species removal" — use `compute_partial_energy(removed_indices)` for vacancies rather than rebuilding `EwaldSummation` per defect configuration.
- **`add_oxidation_state` (explicit dict) raises `ValueError` unless every element is covered**, and `add_oxidation_state_by_site` raises unless every site is covered. Count your lists before calling; the error message says only "Values must be same length"-style length mismatch, not which element is missing.

## 4. Slab generation errors

```python
from pymatgen.core.surface import SlabGenerator

gen = SlabGenerator(initial_structure=bulk,
                    miller_index=(1, 1, 1),
                    min_slab_size=8.0,
                    min_vacuum_size=10.0,
                    lll_reduce=True,
                    center_slab=True)
slabs = gen.get_slabs()
```

Gotchas:

- **`min_slab_size` is in Angstroms, not layers**, by default (`in_unit_planes=False`). A 10 A slab of Cs and a 10 A slab of Fe contain very different atom counts. Set `in_unit_planes=True` to measure sizes in units of the interplanar spacing `d_hkl` when you need comparable atom counts.
- **`get_slabs(symmetrize=True)` can legitimately return an empty list.** If no termination has equivalent top and bottom surfaces, there is no symmetric slab — this is data, not a bug. Drop to `symmetrize=False` to enumerate the asymmetric terminations, and check the empty list is expected rather than assuming a crash.
- **Two opposite surfaces of one slab are different terminations.** A `get_slabs()` result for (1,0,0) returns one slab object per *symmetrically distinct* surface pair, not one per visible termination; both faces of the returned slab are counted. Requesting (2,0,0) does not magically add the missing termination — inspect each slab's two surfaces before concluding a termination was missed.
- **Bond-breaking counts depend on the `bonds` argument.** `get_slabs(bonds={(E1, E2): cutoff}, ftol, tol, max_broken_bonds, repair)` computes the broken-bond energy from the bonds you declare. With `bonds=None`, c-range counting treats every layer separation equally and `slab.energy` can disagree with visual (VESTA) expectations; pass explicit `bonds` for anything non-monatomic. `repair=True` moves atoms to fix broken bonds instead of dropping the termination, which explodes the candidate count — use it only after the unrepaired enumeration is understood.
- **Feed `SlabGenerator` the conventional cell for low-index Miller work.** Passing a primitive cell of a non-cubic lattice with high-index Miller planes requires large `max_normal_search` and often silently under-searches; run `SpacegroupAnalyzer(struct).get_conventional_standard_structure()` first, then generate.

## 5. Defect generation errors

```python
from pymatgen.transformations.standard_transformations import SupercellTransformation

supercell = SupercellTransformation([[2, 0, 0], [0, 2, 0], [0, 0, 2]]).apply_transformation(bulk)
vacancy = supercell.copy()
vacancy.remove_sites([0])   # site indices refer to the supercell, not the unit cell
```

Gotchas:

- **Site indices are cell-local.** After a supercell transformation, index 0 is still just one site of many — confirm the removed species via `supercell[0].species_string` instead of assuming index order survived a reduction or symmetrization step.
- **Oxidation states and site properties do not survive all edits.** `add_oxidation_state_by_guess()` must be re-run after creating a vacancy/interstitial because the composition (and possibly the charge-balanced guess) changed; decoration done on the pristine bulk does not describe the defect cell.
- **Wrap and validate after removing sites.** Run `vacancy.to_unit_cell` (via `from_sites`) and `vacancy.is_valid(tol=0.5)` before IO: some writers emit fractional coords outside [0, 1) that other codes reject.
- **`charge` follows oxidation states, not intent.** `Structure.charge` is the sum of oxidation-state charges; if decoration degraded to all-zero on a non-stoichiometric defect cell, `charge` reports 0 and masks a physically charged defect. Set `struct.charge` explicitly and record the decoration guess when publishing defect results.

## Quick verification checklist

1. `assert all(0 <= c <= 1 - 1e-9 for site in struct for c in site.frac_coords)` after wrapping.
2. `SpacegroupAnalyzer(struct).get_space_group_number()` matches the intended number.
3. `all("+" in s.species_string or "-" in s.species_string for s in struct)` after decoration.
4. `struct.is_valid()` is True before file IO.
5. `struct.charge == 0` before comparing `EwaldSummation.total_energy` across cells.


## Supporting basis and limitations

Grounded in pymatgen 2026.7.27 API docs: Structure constructor params (validate_proximity, to_unit_cell, coords_are_cartesian), from_spacegroup signature and 'only symmetrically distinct' note, add_oxidation_state_by_guess using Composition.oxi_state_guesses taking the first guess, EwaldSummation requiring Species with oxidation states and excluding charged-cell energy. Community evidence: matsci.org thread on all-zero oxidation state guesses (CuO2 -> Cu0+ O02-, CoO4) and the all_oxi_states workaround; GitHub issue #1321 on fractional coord 1.0 vs 0.0 duplicating bonds in neighbor detection; matsci.org thread on fractional-vs-Cartesian confusion with pbc_diff; matsci.org thread on get_slabs(symmetrize=True) returning empty lists; GitHub issue #1894 on SlabGenerator missed unique surfaces (two faces are different terminations); PR #2105 fixing bonds_broken counting in _get_c_ranges; PR #3623 noting add_oxidation_state_by_guess returns self and chains.

## Change and rationale

New skill draft covering pymatgen structure-manipulation gotchas: from_spacegroup vs manual construction (asymmetric unit only, tol merging), periodic-boundary wrapping (to_unit_cell default, Cartesian-vs-fractional confusion, 1.0 vs 0.0 duplicate-bond issue), oxidation-state decoration before Ewald/charge analysis (zero-guess degradation, neutral-cell requirement, compute_partial_energy for vacancies), slab generation (min_slab_size in Angstroms, in_unit_planes, empty symmetrize results, bond-breaking counts, conventional-cell input), and defect generation (cell-local indices, re-decoration, wrapping/validation before IO).

These are the specific silent-failure modes that waste hours: full site lists passed to from_spacegroup, unwrapped coords, zero oxidation guesses feeding Ewald, empty slab lists mistaken for bugs, and wrong vacancy atoms. Each claim is grounded in the pymatgen API docs (pymatgen.org) and documented community issues.