# 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.