vibe-sketcherDocumentation
Source code for vibesketcher.document.slab
"""Slabs from a bulk cell and Miller index, and adsorbates placed on them.
The surface cell is a unimodular integer transform of the bulk cell whose
first two rows span the (hkl) plane, following W. Sun and G. Ceder, Surf.
Sci. 617, 53 (2013), doi:10.1016/j.susc.2013.05.016. Atoms are grouped into
planes by height. A "layer" is either one atomic plane, one repeat of the
oriented bulk cell, or the smallest run of planes whose composition is a
multiple of the bulk formula (a stoichiometric layer, such as one Al-O3-Al
unit of alpha-Al2O3(0001)). Polar or nonstoichiometric terminations are
reported, never neutralised or passivated silently.
"""
from __future__ import annotations
import itertools
import math
from collections import Counter
from functools import reduce
from typing import Any, Literal
import numpy as np
from vibesketcher.document import crystal
from vibesketcher.document.model import PeriodicStructure
LayerUnit = Literal["stoichiometric", "plane", "cell"]
def _reduced(miller: tuple[int, int, int]) -> tuple[int, int, int]:
if not all(isinstance(x, int) for x in miller) or not any(miller):
raise ValueError("a Miller index is three integers, not all zero")
g = reduce(math.gcd, (abs(x) for x in miller))
return (miller[0] // g, miller[1] // g, miller[2] // g)
[docs]
def surface_transform(miller: tuple[int, int, int], lattice: np.ndarray) -> np.ndarray:
"""Integer rows (v1, v2, v3), det +1: v1 and v2 span (hkl), h . v3 = 1.
Searches short integer vectors, so the in-plane cell is as compact as the
bounded search allows; a Gauss reduction then makes it reduced.
"""
h = np.asarray(_reduced(miller))
for span in (3, 6, 10):
rng = range(-span, span + 1)
vectors = [np.asarray(u) for u in itertools.product(rng, repeat=3) if any(u)]
plane = sorted(
(u for u in vectors if int(h @ u) == 0),
key=lambda u: float(np.linalg.norm(u @ lattice)),
)
pair = None
for i, v1 in enumerate(plane):
for v2 in plane[i + 1 :]:
cross = np.cross(v1, v2)
if np.array_equal(cross, h) or np.array_equal(cross, -h):
pair = (v1, v2)
break
if pair:
break
if pair is None:
continue
v1, v2 = _gauss_reduce(pair[0], pair[1], lattice)
normal = np.cross(v1 @ lattice, v2 @ lattice)
normal /= np.linalg.norm(normal)
outs = [u for u in vectors if int(h @ u) == 1]
if not outs:
continue
# The out-of-plane vector with the smallest in-plane component.
v3 = min(
outs,
key=lambda u: float(np.linalg.norm((u @ lattice) - ((u @ lattice) @ normal) * normal)),
)
matrix = np.asarray([v1, v2, v3])
det = round(float(np.linalg.det(matrix)))
if det == -1:
matrix = np.asarray([v2, v1, v3])
if round(float(np.linalg.det(matrix))) != 1:
continue
if float(np.cross(matrix[0] @ lattice, matrix[1] @ lattice) @ (matrix[2] @ lattice)) < 0:
matrix[2] = -matrix[2] # pragma: no cover - det +1 keeps handedness
return matrix.astype(int)
raise ValueError(f"no compact surface cell found for Miller index {miller}")
def _gauss_reduce(
v1: np.ndarray, v2: np.ndarray, lattice: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
"""Lagrange-Gauss reduction of a 2-D integer basis in the metric of ``lattice``."""
a, b = v1.copy(), v2.copy()
for _ in range(100):
la, lb = a @ lattice, b @ lattice
if la @ la > lb @ lb:
a, b = b, a
continue
mu = round(float((la @ lb) / (la @ la)))
if mu == 0:
break
b = b - mu * a
# Prefer an obtuse-free (<= 90 degree) setting for figures and site naming.
if (
float((a @ lattice) @ (b @ lattice)) < -1e-9
and np.linalg.norm((a + b) @ lattice) < np.linalg.norm(b @ lattice) + 1e-9
):
b = a + b
return a, b
def _bulk_formula(structure: PeriodicStructure) -> Counter[str]:
counts = Counter(s.symbol for s in structure.sites)
g = reduce(math.gcd, counts.values())
return Counter({k: v // g for k, v in counts.items()})
def _is_multiple(counts: Counter[str], formula: Counter[str]) -> bool:
if set(counts) != set(formula):
return False
ratios = {counts[k] / formula[k] for k in formula}
return len(ratios) == 1 and float(next(iter(ratios))).is_integer()
def _composition(counts: Counter[str]) -> str:
return "".join(f"{k}{v if v > 1 else ''}" for k, v in sorted(counts.items()))
[docs]
def build_slab(
bulk: PeriodicStructure,
miller: tuple[int, int, int],
*,
identifier: str,
layers: int,
layer_unit: LayerUnit = "stoichiometric",
termination: int = 0,
vacuum: float = 15.0,
repeat: tuple[int, int] | list[list[int]] = (1, 1),
fixed_layers: int = 0,
periodicity: Literal["2d", "3d"] = "2d",
plane_tolerance: float = 0.05,
surface_cell: Literal["primitive", "given"] = "primitive",
) -> PeriodicStructure:
"""Cut a slab of ``layers`` layers from a 3-D bulk structure.
``miller`` refers to the bulk structure's own cell rows, so (0, 0, 1) of a
hexagonal conventional corundum cell is the basal (0001) plane.
``termination`` selects the bottom plane among the distinct planes of one
bulk repeat. ``vacuum`` (angstrom) is the gap between the slab and its image
along the normal; ``periodicity="2d"`` gives a true 2-D model (PBC
T, T, F), ``"3d"`` a vacuum-padded 3-D cell. The bottom ``fixed_layers``
layers are listed as fixed sites for optimisers. With
``surface_cell="primitive"`` a conventional cell of a centred lattice is
first reduced to its primitive cell and the Miller index transformed with
it, so the in-plane cell is the smallest surface cell of the crystal.
"""
if sum(bulk.pbc) != 3:
raise ValueError("slabs are cut from a three-dimensional bulk structure")
if layers < 1 or layers > 200:
raise ValueError("layers must be between 1 and 200")
if not 0 <= fixed_layers < layers + (layers == 0):
raise ValueError("fixed layers must be fewer than the slab's layers")
if not math.isfinite(vacuum) or vacuum <= 0:
raise ValueError("vacuum must be a positive length")
given_miller = tuple(int(x) for x in miller)
bravais = crystal._builder(bulk).get("bravais")
if (
surface_cell == "primitive"
and bulk.cell_setting == "conventional"
and isinstance(bravais, str)
and len(crystal.centring_vectors(bravais)) > 1
):
# Miller indices transform like the cell rows: h' = P h.
scaled = np.rint(6 * crystal.primitive_transform(bravais) @ np.asarray(given_miller))
miller = _reduced(tuple(int(x) for x in scaled)) # type: ignore[arg-type]
bulk = crystal.primitive_cell(bulk, bravais, identifier=bulk.id)
lattice = crystal.lattice_rows(bulk)
transform = surface_transform(tuple(int(x) for x in miller), lattice) # type: ignore[arg-type]
oriented = transform @ lattice
a1, a2, a3 = oriented
e1 = a1 / np.linalg.norm(a1)
e3 = np.cross(a1, a2)
e3 /= np.linalg.norm(e3)
e2 = np.cross(e3, e1)
rotation = np.asarray([e1, e2, e3]).T # columns: new axes in old Cartesian
inverse = np.linalg.inv(oriented)
cell_sites = []
for site in bulk.sites:
frac = np.asarray(bulk.cartesian(site)) @ inverse
frac -= np.floor(frac + 1e-9)
cell_sites.append((site.symbol, site.id, frac))
spacing = float(a3 @ e3)
heights = sorted({round(float(f[2]) * spacing / plane_tolerance) for _, _, f in cell_sites})
# Group heights into planes (merge heights within the tolerance, wrapping).
planes: list[float] = []
for z in (h * plane_tolerance for h in heights):
if not planes or z - planes[-1] > plane_tolerance * 1.5:
planes.append(z)
if len(planes) > 1 and planes[0] + spacing - planes[-1] <= plane_tolerance * 1.5:
# The top plane is the bottom plane's periodic image: the bottom plane
# starts at its lowest atom, one repeat below (a rumpled plane split by
# the cell boundary must not be cut in two).
planes[0] = planes.pop() - spacing
if not 0 <= termination < len(planes):
raise ValueError(f"termination must be between 0 and {len(planes) - 1} for this cut")
formula = _bulk_formula(bulk)
# Enough repeats of the oriented cell to cover the requested layers.
atoms_per_cell = len(cell_sites)
repeats = layers + 2 if layer_unit != "plane" else math.ceil(layers / len(planes)) + 2
stacked = []
for n in range(-1, repeats * (2 if layer_unit == "stoichiometric" else 1) + 1):
for symbol, parent, frac in cell_sites:
f = frac + np.asarray([0, 0, n])
stacked.append((symbol, parent, f, float(f[2]) * spacing))
start = planes[termination] - plane_tolerance * 0.75
stacked = [s for s in stacked if s[3] >= start]
stacked.sort(key=lambda s: s[3])
# Planes in the stacked slab, bottom up.
groups: list[list[tuple[str, str, np.ndarray, float]]] = []
for atom in stacked:
if groups and atom[3] - groups[-1][-1][3] <= plane_tolerance * 1.5:
groups[-1].append(atom)
else:
groups.append([atom])
chosen: list[list[tuple[str, str, np.ndarray, float]]] = []
layer_of: list[int] = []
if layer_unit == "plane":
chosen = groups[:layers]
layer_of = list(range(layers))
elif layer_unit == "cell":
top = start + layers * spacing
chosen = [g for g in groups if g[0][3] < top]
layer_of = [int((g[0][3] - start) // spacing) for g in chosen]
else:
position = 0
for layer in range(layers):
counts: Counter[str] = Counter()
for index in range(position, len(groups)):
counts.update(a[0] for a in groups[index])
chosen.append(groups[index])
layer_of.append(layer)
if _is_multiple(counts, formula):
position = index + 1
break
else:
raise ValueError("no stoichiometric layer found from this termination")
if len(chosen) > 4 * atoms_per_cell * layers:
raise ValueError("stoichiometric layers did not close; use layer_unit='plane'")
if len(chosen) < (layers if layer_unit == "plane" else 1):
raise ValueError("the slab has fewer planes than requested")
atoms = [(a, layer_of[i]) for i, g in enumerate(chosen) for a in g]
xyz = np.asarray([a[2] @ oriented for a, _ in atoms]) @ rotation
xyz[:, 2] -= xyz[:, 2].min()
thickness = float(xyz[:, 2].max())
in_plane = np.asarray([a1 @ rotation, a2 @ rotation])
in_plane[:, 2] = 0.0
matrix = _repeat_matrix(repeat)
surface_rows = np.asarray(matrix, dtype=float) @ in_plane
normal_length = thickness + vacuum
cell = np.vstack([surface_rows, [0.0, 0.0, normal_length]])
xyz[:, 2] += vacuum / 2
# Replicate in plane by the integer matrix (coset representatives).
images = _in_plane_images(matrix)
inverse_cell = np.linalg.inv(cell)
sites = []
fixed: list[int] = []
layer_index: list[int] = []
parents: list[str] = []
for image in images:
shift = image[0] * in_plane[0] + image[1] * in_plane[1]
for (atom, layer), position in zip(atoms, xyz, strict=True):
frac = (position + shift) @ inverse_cell
frac[:2] -= np.floor(frac[:2] + 1e-9)
sites.append((atom[0], frac))
layer_index.append(layer)
parents.append(atom[1])
if layer < fixed_layers:
fixed.append(len(sites) - 1)
ids = crystal._site_ids([s for s, _ in sites])
bottom = Counter(a[0] for a in chosen[0])
top_plane = Counter(a[0] for a in chosen[-1])
total = Counter(s for s, _ in sites)
warnings = []
if not _is_multiple(total, formula):
warnings.append(f"slab composition {_composition(total)} is not a multiple of the bulk")
if bottom != top_plane:
warnings.append(
f"asymmetric terminations: bottom {_composition(bottom)}, "
f"top {_composition(top_plane)}; "
"the slab may carry a dipole"
)
builder = {
"kind": "slab",
"source": bulk.id,
"miller": list(given_miller),
"miller_in_cell": list(miller),
"surface_cell": surface_cell,
"surface_transform": transform,
"layers": layers,
"layer_unit": layer_unit,
"termination": termination,
"planes_per_repeat": len(planes),
"vacuum": vacuum,
"thickness": thickness,
"repeat": matrix,
"periodicity": periodicity,
"fixed_sites": [ids[i] for i in fixed],
"layer_of_site": dict(zip(ids, layer_index, strict=True)),
"top_layer": max(layer_index),
"bottom_composition": _composition(bottom),
"top_composition": _composition(top_plane),
"warnings": warnings,
"reference": "W. Sun and G. Ceder, Surf. Sci. 617, 53 (2013)",
}
pbc = (True, True, periodicity == "3d")
structure = crystal._structure(
identifier, cell, sites, pbc, "unspecified", builder, site_ids=ids
)
return structure.model_copy(
update={
"provenance": {
**structure.provenance,
"parent_sites": dict(zip(ids, parents, strict=True)),
}
}
)
def _repeat_matrix(repeat: Any) -> list[list[int]]:
if len(repeat) == 2 and all(isinstance(x, int) for x in repeat):
matrix = [[int(repeat[0]), 0], [0, int(repeat[1])]]
else:
matrix = [[int(x) for x in row] for row in repeat]
if len(matrix) != 2 or any(len(r) != 2 for r in matrix):
raise ValueError("the in-plane repeat is n1, n2 or a 2x2 integer matrix")
det = matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]
if not 1 <= det <= 400 or any(abs(x) > 50 for r in matrix for x in r):
raise ValueError("the in-plane repeat needs a positive determinant up to 400")
return matrix
def _in_plane_images(matrix: list[list[int]]) -> list[tuple[int, int]]:
m = np.asarray(matrix, dtype=float)
inverse = np.linalg.inv(m)
lo = int(min(0, m[:, 0].min(), m[:, 0].sum())) - 1
hi = int(max(0, m[:, 0].max(), m[:, 0].sum())) + 1
lo2 = int(min(0, m[:, 1].min(), m[:, 1].sum())) - 1
hi2 = int(max(0, m[:, 1].max(), m[:, 1].sum())) + 1
images = []
for i, j in itertools.product(range(lo, hi + 1), range(lo2, hi2 + 1)):
f = np.asarray([i, j]) @ inverse
if np.all(f >= -1e-9) and np.all(f < 1 - 1e-9):
images.append((i, j))
det = round(abs(float(np.linalg.det(m))))
if len(images) != det:
raise ValueError("in-plane repeat representatives could not be qualified")
return images
Orientation = Literal["flat", "upright", "as-is"]
def _orient(xyz: np.ndarray, orientation: Orientation, rotate: float) -> np.ndarray:
centred = xyz - xyz.mean(axis=0)
if orientation != "as-is" and len(xyz) >= 2:
_, _, vt = np.linalg.svd(centred)
# vt rows: largest to smallest spread; build a right-handed frame.
if orientation == "flat":
frame = np.asarray([vt[0], vt[1], np.cross(vt[0], vt[1])])
else:
frame = np.asarray([vt[1], np.cross(vt[0], vt[1]), vt[0]])
centred = centred @ frame.T
angle = math.radians(rotate)
spin = np.asarray(
[[math.cos(angle), math.sin(angle), 0], [-math.sin(angle), math.cos(angle), 0], [0, 0, 1]]
)
return np.asarray(centred @ spin)
[docs]
def adsorb(
slab: PeriodicStructure,
molecule: list[tuple[str, tuple[float, float, float]]],
*,
identifier: str,
site: Literal["top", "bridge", "hollow", "position"] = "top",
anchor: str | None = None,
position: tuple[float, float] | None = None,
height: float = 2.5,
orientation: Orientation = "flat",
rotate: float = 0.0,
attach: int | None = None,
min_distance: float = 1.2,
label: str = "adsorbate",
) -> PeriodicStructure:
"""Place a molecule above the top surface of a slab.
``site`` names the reference point: ``top`` of the ``anchor`` surface atom
(default: the top-most atom), the ``bridge`` between it and its nearest
top-layer neighbour, the ``hollow`` centroid with its two nearest
neighbours, or an explicit in-plane fractional ``position``. ``height`` is
the vertical distance (angstrom) of the molecule's ``attach`` atom (default:
its lowest atom after orientation) above that point's surface height.
``flat`` puts the molecular plane parallel to the surface (benzene lying
down), ``upright`` its long axis along the normal. Placements closer than
``min_distance`` to any slab atom or periodic image are refused.
"""
if sum(slab.pbc[:2]) != 2 or crystal._builder(slab).get("kind") not in ("slab", "adsorbate"):
raise ValueError("adsorbates go on a slab built by the slab builder")
if not molecule:
raise ValueError("the adsorbate has no atoms")
if not math.isfinite(height) or height <= 0:
raise ValueError("height must be a positive length")
lattice = crystal.lattice_rows(slab)
xyz = crystal.cartesian(slab)
surface_z = float(xyz[:, 2].max())
top_ids = [i for i in range(len(xyz)) if xyz[i, 2] > surface_z - 0.6]
ids = [s.id for s in slab.sites]
if site == "position":
if position is None or anchor is not None:
raise ValueError("site 'position' takes an in-plane fractional position and no anchor")
point = position[0] * lattice[0] + position[1] * lattice[1]
point[2] = surface_z
else:
if position is not None:
raise ValueError(f"site {site!r} is located by an anchor atom, not a position")
if anchor is None:
centre = 0.5 * (lattice[0] + lattice[1])
centre[2] = surface_z
first = min(
top_ids,
key=lambda i: float(np.linalg.norm(_nearest(xyz[i], centre, lattice) - centre)),
)
else:
if anchor not in ids:
raise ValueError(f"anchor site {anchor!r} not found")
first = ids.index(anchor)
origin = xyz[first]
neighbours = sorted(
(float(np.linalg.norm(_nearest(xyz[i], origin, lattice)[:2] - origin[:2])), i)
for i in top_ids
if i != first
)
if site == "top":
point = origin.copy()
elif site == "bridge":
if not neighbours:
raise ValueError("a bridge site needs a second top-layer atom")
point = 0.5 * (origin + _nearest(xyz[neighbours[0][1]], origin, lattice))
else:
if len(neighbours) < 2:
raise ValueError("a hollow site needs three top-layer atoms")
second = _nearest(xyz[neighbours[0][1]], origin, lattice)
third_candidates = [_nearest(xyz[i], origin, lattice) for _, i in neighbours[1:]]
third = min(
third_candidates,
key=lambda p: float(
np.linalg.norm(p[:2] - second[:2]) + np.linalg.norm(p[:2] - origin[:2])
),
)
point = (origin + second + third) / 3
point[2] = max(point[2], surface_z if site != "top" else point[2])
placed = _orient(np.asarray([p for _, p in molecule], dtype=float), orientation, rotate)
if attach is None:
attach = int(np.argmin(placed[:, 2]))
if not 0 <= attach < len(molecule):
raise ValueError("attach must index an adsorbate atom")
placed += point - placed[attach] + np.asarray([0, 0, height])
if placed[attach, 0] != point[0]:
placed[:, :2] -= placed[attach, :2] - point[:2]
# Collision checks against the slab and the adsorbate's own images.
shifts = [i * lattice[0] + j * lattice[1] for i in (-1, 0, 1) for j in (-1, 0, 1)]
closest = min(float(np.linalg.norm(p - (q + s))) for p in placed for q in xyz for s in shifts)
if closest < min_distance:
raise ValueError(f"the adsorbate comes within {closest:.2f} A of the slab")
own = [
float(np.linalg.norm(p - (q + s)))
for p in placed
for q in placed
for s in shifts
if np.linalg.norm(s) > 1e-9
]
if own and min(own) < max(min_distance, 1.6):
raise ValueError(
f"the adsorbate comes within {min(own):.2f} A of its periodic image; "
"use a larger in-plane repeat"
)
normal = lattice[2]
top = float(placed[:, 2].max())
data = slab.model_dump()
cell = lattice.copy()
if not slab.pbc[2] or top + 2.0 > normal[2]:
# Keep the molecule inside the cell; the gap above it stays at least the slab vacuum.
builder_slab = crystal._builder(slab)
recorded = builder_slab.get("vacuum", 15.0)
vacuum = float(recorded) if isinstance(recorded, (int, float)) else 15.0
cell[2] = [0.0, 0.0, max(normal[2], top + vacuum / 2)]
inverse = np.linalg.inv(cell)
counts: Counter[str] = Counter(s.symbol for s in slab.sites)
new_sites = []
for site_data in data["sites"]:
cart = np.asarray(site_data["fractional"]) @ lattice
site_data["fractional"] = tuple(crystal._clean(x) for x in cart @ inverse)
new_sites.append(site_data)
added = []
for (symbol, _), p in zip(molecule, placed, strict=True):
counts[symbol] += 1
sid = f"{symbol}{counts[symbol]}"
while sid in ids:
counts[symbol] += 1
sid = f"{symbol}{counts[symbol]}"
added.append(sid)
new_sites.append(
{
"id": sid,
"symbol": symbol,
"fractional": tuple(crystal._clean(x) for x in p @ inverse),
}
)
builder = crystal._builder(slab)
adsorbates = list(builder.get("adsorbates") or [])
adsorbates.append(
{
"label": label,
"sites": added,
"site": site,
"anchor": None if site == "position" else ids[first],
"position": list(position) if position else None,
"height": height,
"orientation": orientation,
"rotate": rotate,
"closest_slab_distance": round(closest, 4),
}
)
builder.update({"kind": "adsorbate", "adsorbates": adsorbates, "slab": slab.id})
data.update(
{
"id": identifier,
"lattice_angstrom": tuple(tuple(float(x) for x in row) for row in cell),
"sites": new_sites,
"bonds": [],
"provenance": {**slab.provenance, "builder": crystal._json(builder)},
"source_cif": None,
"source_model_sha256": None,
"charge": slab.charge,
"multiplicity": slab.multiplicity,
}
)
return PeriodicStructure.model_validate(data)
def _nearest(point: np.ndarray, reference: np.ndarray, lattice: np.ndarray) -> np.ndarray:
"""The in-plane periodic image of ``point`` nearest ``reference``."""
best = point
distance = math.inf
for i, j in itertools.product((-1, 0, 1), repeat=2):
candidate = point + i * lattice[0] + j * lattice[1]
d = float(np.linalg.norm(candidate[:2] - reference[:2]))
if d < distance:
best, distance = candidate, d
return np.asarray(best)
DefectKind = Literal["vacancy", "substitution", "interstitial"]
[docs]
def apply_defects(
structure: PeriodicStructure,
defects: list[dict[str, Any]],
*,
identifier: str,
charge: int | None = None,
multiplicity: int | None = None,
) -> PeriodicStructure:
"""Vacancies, substitutions and interstitials on named sites of a pristine structure.
Each defect is ``{"kind": "vacancy", "site": id}``,
``{"kind": "substitution", "site": id, "symbol": "Al"}`` or
``{"kind": "interstitial", "symbol": "H", "fractional": [x, y, z]}``. The
pristine structure, the defect map and the stoichiometric change are
recorded; the charge and spin are what the caller states, never inferred.
"""
from vibesketcher.elements import ATOMIC_NUMBER
if not defects:
raise ValueError("give at least one defect")
data = structure.model_dump()
sites = {s["id"]: s for s in data["sites"]}
order = [s["id"] for s in data["sites"]]
removed: list[str] = []
record = []
counts = Counter(s.symbol for s in structure.sites)
for defect in defects:
kind = defect.get("kind")
if kind in ("vacancy", "substitution"):
target = defect.get("site")
if target not in sites or target in removed:
raise ValueError(f"defect site {target!r} not found")
if kind == "vacancy":
removed.append(str(target))
record.append({"kind": kind, "site": target, "symbol": sites[target]["symbol"]})
else:
symbol = defect.get("symbol")
if symbol not in ATOMIC_NUMBER or symbol == sites[target]["symbol"]:
raise ValueError("a substitution needs a different element")
record.append(
{"kind": kind, "site": target, "from": sites[target]["symbol"], "to": symbol}
)
sites[target] = {**sites[target], "symbol": symbol}
elif kind == "interstitial":
symbol = defect.get("symbol")
frac = defect.get("fractional")
if symbol not in ATOMIC_NUMBER or not isinstance(frac, (list, tuple)) or len(frac) != 3:
raise ValueError("an interstitial needs an element and a fractional position")
counts[str(symbol)] += 1
sid = f"{symbol}i{counts[str(symbol)]}"
while sid in sites:
counts[str(symbol)] += 1
sid = f"{symbol}i{counts[str(symbol)]}"
position = np.asarray([float(x) for x in frac])
lattice = crystal.lattice_rows(structure)
cart = position @ lattice
shifts = list(itertools.product(*[(-1, 0, 1) if p else (0,) for p in structure.pbc]))
for other in structure.sites:
if other.id in removed:
continue
for shift in shifts:
d = float(
np.linalg.norm(
np.asarray(structure.cartesian(other))
+ np.asarray(shift) @ lattice
- cart
)
)
if d < 0.7:
raise ValueError(f"the interstitial sits {d:.2f} A from {other.id}")
sites[sid] = {"id": sid, "symbol": symbol, "fractional": tuple(float(x) for x in frac)}
order.append(sid)
record.append({"kind": kind, "site": sid, "symbol": symbol, "fractional": list(frac)})
else:
raise ValueError(f"unknown defect kind {kind!r}")
remaining = [sites[i] for i in order if i not in removed]
if not remaining:
raise ValueError("defects removed every site")
before = Counter(s.symbol for s in structure.sites)
after = Counter(s["symbol"] for s in remaining)
delta = {
k: after[k] - before[k] for k in sorted(set(before) | set(after)) if after[k] != before[k]
}
builder = crystal._builder(structure)
builder.update({"defects": record, "pristine": structure.id, "stoichiometry_change": delta})
data.update(
{
"id": identifier,
"sites": remaining,
"bonds": [b for b in data["bonds"] if b["a"] not in removed and b["b"] not in removed],
"charge": charge,
"multiplicity": multiplicity,
"provenance": {**structure.provenance, "builder": crystal._json(builder)},
"source_cif": None,
"source_model_sha256": None,
}
)
if structure.parent is not None:
for site in remaining:
site.setdefault("parent_site", site["id"])
if builder.get("fixed_sites"):
builder["fixed_sites"] = [s for s in builder["fixed_sites"] if s not in removed]
data["provenance"]["builder"] = crystal._json(builder)
return PeriodicStructure.model_validate(data)
#: Standard anion charges used to assign formal charges in binary ionic crystals.
_ANION_CHARGE = {"O": -2, "S": -2, "Se": -2, "Te": -2, "F": -1, "Cl": -1, "Br": -1, "I": -1,
"N": -3, "P": -3, "As": -3, "Sb": -3} # fmt: skip
[docs]
def formal_charges(structure: PeriodicStructure) -> dict[str, float] | None:
"""Formal ionic charges of a binary compound, or None where they are not defined.
The anion takes its standard charge and the cation the charge that makes
the composition neutral. Elements, ternary and covalent (no tabulated
anion) compounds give None: no polarity verdict is made for them.
"""
counts = Counter(s.symbol for s in structure.sites)
anions = [e for e in counts if e in _ANION_CHARGE]
if len(counts) != 2 or len(anions) != 1:
return None
anion = anions[0]
(cation,) = [e for e in counts if e != anion]
q = _ANION_CHARGE[anion]
return {anion: float(q), cation: -q * counts[anion] / counts[cation]}
[docs]
def polarity(made: PeriodicStructure, charges: dict[str, float] | None) -> dict[str, Any] | None:
"""Net formal charge and dipole per area along the normal (e/Å) of a slab (Tasker 1979)."""
if charges is None:
return None
xyz = crystal.cartesian(made)
rows = crystal.lattice_rows(made)
area = float(np.linalg.norm(np.cross(rows[0], rows[1])))
q = np.asarray([charges[s.symbol] for s in made.sites])
z = xyz[:, 2] - xyz[:, 2].mean()
dipole = float(q @ z) / area
charge = float(q.sum())
return {
"charge": round(charge, 6),
"dipole_e_per_angstrom": round(dipole, 6),
"polar": abs(dipole) > 1e-3 or abs(charge) > 1e-6,
"charges": charges,
"reference": "P. W. Tasker, J. Phys. C 12, 4977 (1979)",
}
[docs]
def terminations(
bulk: PeriodicStructure,
miller: tuple[int, int, int],
*,
layers: int = 1,
layer_unit: LayerUnit = "stoichiometric",
surface_cell: Literal["primitive", "given"] = "primitive",
) -> list[dict[str, Any]]:
"""Every bottom termination of a cut, with what a slab from it would be (#232).
Each entry gives the bottom and top plane compositions, whether the slab is
stoichiometric and has equal terminations, its thickness and atom count,
the builder's warnings and, for binary ionic crystals, its formal-charge
polarity (:func:`polarity`). Terminations that give the same slab (same
compositions, thickness and atom count) point to the first one with
``same_as``.
"""
if sum(bulk.pbc) != 3:
raise ValueError("slabs are cut from a three-dimensional bulk structure")
charges = formal_charges(bulk)
out: list[dict[str, Any]] = []
seen: dict[tuple[Any, ...], int] = {}
index = 0
while index < 400: # build_slab refuses an out-of-range termination first
try:
made = build_slab(bulk, miller, identifier="probe", layers=layers,
layer_unit=layer_unit, termination=index,
surface_cell=surface_cell) # fmt: skip
except ValueError as error:
if str(error).startswith("termination must be between"):
break
out.append({"termination": index, "error": str(error)})
index += 1
continue
index += 1
record = crystal._builder(made)
warnings = list(record.get("warnings", []))
entry: dict[str, Any] = {
"termination": index - 1,
"bottom": record["bottom_composition"],
"top": record["top_composition"],
"sites": len(made.sites),
"thickness": round(float(record["thickness"]), 4),
"stoichiometric": not any("not a multiple" in w for w in warnings),
"symmetric_terminations": record["bottom_composition"] == record["top_composition"],
"warnings": warnings,
"polarity": polarity(made, charges),
}
key = (entry["bottom"], entry["top"], entry["sites"], round(entry["thickness"], 2))
if key in seen:
entry["same_as"] = seen[key]
else:
seen[key] = index - 1
out.append(entry)
if len(out) == 400 and all("error" in e for e in out):
raise ValueError(out[0]["error"])
return out