Disorder, vacancies, and attached species

A site in a real material is not always one atom. It can be an alloy site shared between two elements, a site that is only partly occupied, or an atom carrying particles whose positions nobody modelled. This notebook walks through how httk-atomistic represents all three, and — the practical part — which of its representations can carry them and which cannot.

Every code cell below is taken verbatim from the runnable script examples/disorder_and_vacancies.py. The notebook stores no outputs; run it to see them, or read the rendered example page Disorder, vacancies, and attached species: what a site can actually hold, which shows the same script.

We start with the imports and four species that between them cover the whole story: a 50/50 Fe/Ni alloy site, a titanium site that is empty one time in ten, a methyl group, and one ordinary fully ordered oxygen.

from httk.core import unwrap

from httk.atomistic import (
    Species,
    SpeciesPrimitiveView,
    Structure,
    StructurePrimitiveView,
    UnitcellStructureView,
)

CUBIC = [[4.0, 0.0, 0.0], [0.0, 4.0, 0.0], [0.0, 0.0, 4.0]]

#: A 50/50 substitutional Fe/Ni alloy site.
FE_NI = Species(name="M", chemical_symbols=("Fe", "Ni"), concentration=(0.5, 0.5))
#: A titanium site that is empty one time in ten.
TI_WITH_VACANCY = Species(name="Ti", chemical_symbols=("Ti", "vacancy"), concentration=(0.9, 0.1))
#: A carbon carrying three hydrogens whose positions are not modelled.
METHYL = Species(name="CH3", chemical_symbols=("C",), concentration=(1.0,), attached=("H",), nattached=(3,))
#: An ordinary, fully ordered species -- the only kind a primitive triple can hold.
OXYGEN = Species(name="O", chemical_symbols=("O",), concentration=(1.0,))

What a Species actually says

A Species carries a name (its label in the structure — it need not be a chemical symbol), a tuple of chemical_symbols, and a matching tuple of concentration values. Read together, those say what fraction of the time each symbol occupies the site:

  • ("Fe", "Ni") at (0.5, 0.5) is substitutional disorder — a statement about an ensemble, not about one atom.

  • ("Ti", "vacancy") at (0.9, 0.1) is a partial vacancy. "vacancy" is one of OPTIMADE’s two pseudo-symbols (the other is "X", meaning unknown), and it names the absence of an atom.

  • attached=("H",) with nattached=(3,) records attached particles: three hydrogens bound to the carbon, with no coordinates invented for them. The two fields must always be given together and share a length.

Species.is_single_element is the predicate that matters downstream: it is true only for exactly one real element symbol with nothing attached. Watch which of the four is true.

def explain(species: Species) -> None:
    """Print what a species says, and whether it survives the trip to an atomic number."""
    occupancy = ", ".join(
        f"{100 * fraction:g}% {symbol}" for symbol, fraction in zip(species.chemical_symbols, species.concentration)
    )
    attachments = ""
    if species.attached is not None and species.nattached is not None:
        attachments = " + attached " + ", ".join(
            f"{count}x{symbol}" for symbol, count in zip(species.attached, species.nattached)
        )
    print(f"  {species.name:<4} {occupancy}{attachments}")
    print(f"       is_single_element = {species.is_single_element}")


def species_are_richer_than_elements() -> None:
    print("What each species says:")
    for species in (FE_NI, TI_WITH_VACANCY, METHYL, OXYGEN):
        explain(species)

species_are_richer_than_elements()

What each species says:
  M    50% Fe, 50% Ni
       is_single_element = False
  Ti   90% Ti, 10% vacancy
       is_single_element = False
  CH3  100% C + attached 3xH
       is_single_element = False
  O    100% O
       is_single_element = True

The Unitcell representation carries all of it

A Structure — and any UnitcellStructureView over any backend — stores the Species objects themselves, so nothing is lost: the concentrations, the vacancy pseudo-symbol, and the attachments all survive intact.

Two things worth noticing in the output. First, unwrap() on the view returns the same object we put in, not a copy — views present, they do not convert. Second, SpeciesPrimitiveView renders each species as an OPTIMADE species dict with every field still present; that is unsurprising, since this model is OPTIMADE’s.

def disorder_survives_the_simple_representation() -> None:
    """The Structure, and any UnitcellStructureView, keep every field intact."""
    structure = Structure(
        cell=CUBIC,
        sites=[[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]],
        species=[TI_WITH_VACANCY, METHYL],
        species_at_sites=["Ti", "CH3"],
    )

    by_name = {species.name: species for species in structure.species}
    print("Unitcell representation (a Structure):")
    print("  Ti  chemical_symbols ", by_name["Ti"].chemical_symbols)
    print("  Ti  concentration    ", by_name["Ti"].concentration)
    print("  CH3 attached         ", by_name["CH3"].attached)
    print("  CH3 nattached        ", by_name["CH3"].nattached)

    # A view over the same structure is the same data seen again, not a lossy copy:
    # unwrap() hands back the original object.
    view = UnitcellStructureView(structure)
    print("  round-trips through UnitcellStructureView:", unwrap(view) is structure)

    # ...and so does the OPTIMADE species dict, which is where this model comes from.
    print("  as OPTIMADE dicts:")
    for species in structure.species:
        print("   ", dict(SpeciesPrimitiveView(species)))

disorder_survives_the_simple_representation()

Unitcell representation (a Structure):
  Ti  chemical_symbols  ('Ti', 'vacancy')
  Ti  concentration     (0.9, 0.1)
  CH3 attached          ('H',)
  CH3 nattached         (3,)
  round-trips through UnitcellStructureView: True
  as OPTIMADE dicts:
    {'name': 'Ti', 'chemical_symbols': ['Ti', 'vacancy'], 'concentration': [0.9, 0.1]}
    {'name': 'CH3', 'chemical_symbols': ['C'], 'concentration': [1.0], 'attached': ['H'], 'nattached': [3]}

The primitive representation rejects unsupported species

The primitive representation is the spglib-style (lattice, positions, numbers) triple that symmetry libraries speak, and numbers is a list of bare atomic numbers. There is no integer that means “90% Ti”, and none that means “carbon plus three hydrogens”.

StructurePrimitiveView does not round, drop, or approximate. It raises TypeError and names the species it could not express, preventing a modelling statement from being silently converted into a different material. A fully ordered structure goes through untouched.

def the_primitive_representation_refuses() -> None:
    """A bare atomic number cannot mean '90% Ti', so conversion raises TypeError."""
    print("Primitive representation (an spglib-style triple):")
    for label, species in (("partial vacancy", TI_WITH_VACANCY), ("attached particles", METHYL), ("alloy", FE_NI)):
        structure = Structure(
            cell=CUBIC,
            sites=[[0.0, 0.0, 0.0]],
            species=[species],
            species_at_sites=[species.name],
        )
        try:
            StructurePrimitiveView(structure)
        except TypeError as error:
            print(f"  {label:<19} -> TypeError: {error}")
        else:  # pragma: no cover - unreachable while the check above holds
            print(f"  {label:<19} -> unexpectedly accepted")

    # A fully ordered structure converts without complaint, of course:
    ordered = Structure(cell=CUBIC, sites=[[0.0, 0.0, 0.0]], species=[OXYGEN], species_at_sites=["O"])
    _lattice, _positions, numbers = StructurePrimitiveView(ordered)
    print("  fully ordered       -> atomic numbers", numbers)

the_primitive_representation_refuses()

Primitive representation (an spglib-style triple):
  partial vacancy     -> TypeError: This structure cannot be represented as a primitive structure (species 'Ti' is not a single, unattached chemical element)
  attached particles  -> TypeError: This structure cannot be represented as a primitive structure (species 'CH3' is not a single, unattached chemical element)
  alloy               -> TypeError: This structure cannot be represented as a primitive structure (species 'M' is not a single, unattached chemical element)
  fully ordered       -> atomic numbers (8,)

Checking before you convert

Since is_single_element is exactly the condition the primitive view enforces, you can ask the question yourself before handing a structure to a symmetry library, a force field, or anything else that speaks in atomic numbers — and get back the names of the offending species rather than an exception.

The rule of thumb: keep disordered structures in the Unitcell representation, and convert to primitive only at the boundary of a tool that genuinely cannot represent disorder, where the TypeError is precisely the error you want.

def guard_before_converting() -> None:
    """The idiomatic pre-flight check before handing a structure to an element-only tool."""
    mixed = Structure(
        cell=CUBIC,
        sites=[[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]],
        species=[FE_NI, OXYGEN],
        species_at_sites=["M", "O"],
    )
    offenders = [species.name for species in mixed.species if not species.is_single_element]
    print("Pre-flight check:")
    print("  species that cannot become an atomic number:", offenders)
    print("  safe to convert to a primitive triple:", not offenders)

guard_before_converting()

Pre-flight check:
  species that cannot become an atomic number: ['M']
  safe to convert to a primitive triple: False

Where to go next