vibe-sketcherDocumentation
Source code for vibesketcher.render.periodic
"""Publication figures of periodic structures (#189).
A figure is drawn from the validated `PeriodicStructure`, never from a
picture: atoms (with display repeats and optional boundary completion),
bonds found from covalent radii across cell boundaries, the outline of the
cell's *periodic* edges, the lattice vectors with their lengths, a Miller
plane, the real-space Wigner-Seitz cell and, in a separate panel, the first
Brillouin zone of the reciprocal lattice with b_i . a_j = 2 pi delta_ij.
Display repeats are drawing only; they never change the model.
The camera looks along a stated direction: ``top`` (along -z), ``side``
(along +y), ``front`` (along -x), ``isometric``, a real-space direction
``[uvw]`` or the normal of a plane ``(hkl)``. Screen axes are right-handed
(right x up = towards the viewer). Atom colours follow the Jmol scheme; radii
are RDKit's covalent radii (Cordero et al., Dalton Trans. 2832 (2008)).
"""
from __future__ import annotations
import html
import itertools
import math
from typing import Literal
import numpy as np
from pydantic import BaseModel, ConfigDict, Field
from vibesketcher.document import crystal
from vibesketcher.document.model import PeriodicSite, PeriodicStructure
from vibesketcher.elements import ATOMIC_NUMBER
View = Literal["top", "side", "front", "isometric", "uvw", "hkl"]
#: Jmol element colours for common elements; others are drawn grey.
_COLOURS = {
"H": "#FFFFFF",
"He": "#D9FFFF",
"Li": "#CC80FF",
"Be": "#C2FF00",
"B": "#FFB5B5",
"C": "#909090",
"N": "#3050F8",
"O": "#FF0D0D",
"F": "#90E050",
"Na": "#AB5CF2",
"Mg": "#8AFF00",
"Al": "#BFA6A6",
"Si": "#F0C8A0",
"P": "#FF8000",
"S": "#FFFF30",
"Cl": "#1FF01F",
"K": "#8F40D4",
"Ca": "#3DFF00",
"Ti": "#BFC2C7",
"V": "#A6A6AB",
"Cr": "#8A99C7",
"Mn": "#9C7AC7",
"Fe": "#E06633",
"Co": "#F090A0",
"Ni": "#50D050",
"Cu": "#C88033",
"Zn": "#7D80B0",
"Ga": "#C28F8F",
"Ge": "#668F8F",
"As": "#BD80E3",
"Se": "#FFA100",
"Br": "#A62929",
"Sr": "#00FF00",
"Zr": "#94E0E0",
"Mo": "#54B5B5",
"Ru": "#248F8F",
"Rh": "#0A7D8C",
"Pd": "#006985",
"Ag": "#C0C0C0",
"Cd": "#FFD98F",
"Sn": "#668080",
"I": "#940094",
"Cs": "#57178F",
"Ba": "#00C900",
"La": "#70D4FF",
"Ce": "#FFFFC7",
"Pt": "#D0D0E0",
"Au": "#FFD123",
"W": "#2194D6",
"Pb": "#575961",
}
[docs]
class FigureOptions(BaseModel):
"""What a periodic figure shows; every choice is stated, none inferred."""
model_config = ConfigDict(extra="forbid")
view: View = "top"
direction: tuple[int, int, int] | None = None
repeat: tuple[int, int, int] = (1, 1, 1)
complete_boundary: bool = True
cell: bool = True
vectors: bool = True
bonds: bool = True
labels: Literal["none", "element", "site", "wyckoff"] = "none"
legend: bool = True
scale_bar: bool = True
plane: tuple[int, int, int] | None = None
wigner_seitz: bool = False
brillouin: bool = False
#: Schematic aids (#233); each is drawn from the unchanged 3D model.
#: Heights of the atoms along the view, written beside them (atoms that
#: coincide in projection share one label).
heights: bool = False
#: Atoms further from the viewer drawn smaller and paler.
depth_cue: bool = False
#: Elements whose coordination polyhedra are drawn, e.g. ["Ti"].
polyhedra: list[str] = Field(default_factory=list, max_length=6)
#: Slab layers numbered from the bottom, fixed layers marked (side views).
layers: bool = False
#: Stacking letters (A, B, C, ...) of the atomic planes of a slab (side views).
stacking: bool = False
#: The space group's symmetry elements in the projection, as in International
#: Tables diagrams (views along a, b or c of a 3D crystal).
symmetry_elements: bool = False
#: With symmetry_elements: also the two-fold axes lying in the projection plane.
inplane_axes: bool = False
#: Length tolerance (angstrom) when spglib finds the group or the Wyckoff letters.
symprec: float = Field(default=1e-3, gt=0, le=1)
atom_scale: float = Field(default=0.45, gt=0, le=2)
width: float = Field(default=520, ge=120, le=4000)
title: str | None = Field(default=None, max_length=200)
[docs]
class FigureRefused(ValueError):
"""The options ask for something this structure cannot show."""
def _radius(symbol: str) -> float:
from rdkit import Chem
return float(Chem.GetPeriodicTable().GetRcovalent(ATOMIC_NUMBER[symbol]))
def _unit(v: np.ndarray) -> np.ndarray:
n = float(np.linalg.norm(v))
if n < 1e-12:
raise FigureRefused("the view direction is zero")
return np.asarray(v / n)
[docs]
def view_direction(structure: PeriodicStructure, options: FigureOptions) -> np.ndarray:
"""The unit vector the camera looks along, in Cartesian angstrom space."""
rows = crystal.lattice_rows(structure)
if options.view in ("uvw", "hkl"):
if options.direction is None or not any(options.direction):
raise FigureRefused(f"view {options.view} needs a nonzero direction")
d = np.asarray(options.direction, dtype=float)
if options.view == "uvw":
return _unit(-(d @ rows)) # look down the direction, from its tip
reciprocal = np.linalg.inv(rows).T # rows b_i / 2 pi
return _unit(-(d @ reciprocal))
if options.direction is not None:
raise FigureRefused("direction applies to the uvw and hkl views")
return {
"top": np.asarray([0.0, 0.0, -1.0]),
"side": np.asarray([0.0, 1.0, 0.0]),
"front": np.asarray([-1.0, 0.0, 0.0]),
"isometric": _unit(np.asarray([-1.0, -0.8, -0.6])),
}[options.view]
[docs]
def camera(direction: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Screen right and up for a camera looking along ``direction``."""
d = _unit(direction)
up = np.asarray([0.0, 0.0, 1.0])
up = up - (up @ d) * d
if np.linalg.norm(up) < 1e-6:
up = np.asarray([0.0, 1.0, 0.0]) - d[1] * d
up = _unit(up)
right = np.cross(up, -d)
return _unit(right), up
def _images(
structure: PeriodicStructure, options: FigureOptions
) -> list[tuple[PeriodicSite, np.ndarray]]:
rows = crystal.lattice_rows(structure)
repeat = [n if p else 1 for n, p in zip(options.repeat, structure.pbc, strict=True)]
if any(n < 1 or n > 20 for n in repeat):
raise FigureRefused("display repeats must be between 1 and 20")
if math.prod(repeat) * len(structure.sites) > 20000:
raise FigureRefused("the figure would draw more than 20000 atoms")
out = []
for site in structure.sites:
frac = np.asarray(site.fractional, dtype=float)
shifts = []
for axis in range(3):
if not structure.pbc[axis]:
shifts.append([0])
continue
base = list(range(repeat[axis]))
# Atoms on a cell face are drawn on the opposite face too.
if options.complete_boundary and abs(frac[axis] - round(frac[axis])) < 1e-3:
base.append(repeat[axis])
shifts.append(base)
wrapped = frac.copy()
for axis in range(3):
if structure.pbc[axis] and abs(wrapped[axis] - round(wrapped[axis])) < 1e-3:
wrapped[axis] = 0.0
for shift in itertools.product(*shifts):
out.append((site, (wrapped + np.asarray(shift)) @ rows))
return out
[docs]
def plane_polygon(structure: PeriodicStructure, hkl: tuple[int, int, int]) -> np.ndarray:
"""Vertices (Cartesian) where the plane hx + ky + lz = 1 cuts the unit cell."""
h = np.asarray(hkl, dtype=float)
if not np.any(h):
raise FigureRefused("a Miller plane needs a nonzero index")
corners = np.asarray(list(itertools.product((0, 1), repeat=3)), dtype=float)
points: list[np.ndarray] = []
for i, j in itertools.combinations(range(8), 2):
a, b = corners[i], corners[j]
if np.count_nonzero(a != b) != 1:
continue
fa, fb = h @ a - 1, h @ b - 1
if abs(fa - fb) < 1e-12:
continue
t = fa / (fa - fb)
if -1e-9 <= t <= 1 + 1e-9:
p = a + t * (b - a)
if not any(np.allclose(p, q) for q in points):
points.append(p)
if len(points) < 3:
raise FigureRefused(f"the plane ({hkl[0]} {hkl[1]} {hkl[2]}) does not cut the unit cell")
pts = np.asarray(points) @ crystal.lattice_rows(structure)
centre = pts.mean(axis=0)
normal = _unit(h @ np.linalg.inv(crystal.lattice_rows(structure)).T)
e1 = _unit(pts[0] - centre)
e2 = np.cross(normal, e1)
order = np.argsort([math.atan2((p - centre) @ e2, (p - centre) @ e1) for p in pts])
return np.asarray(pts[order])
[docs]
def reciprocal_structure(structure: PeriodicStructure) -> PeriodicStructure:
"""The reciprocal lattice (rows b_i, 1/angstrom, b_i . a_j = 2 pi delta_ij) as a cell."""
rows = crystal.lattice_rows(structure)
b = 2 * math.pi * np.linalg.inv(rows).T
return PeriodicStructure(
id="reciprocal",
lattice_angstrom=tuple(tuple(float(x) for x in r) for r in b),
pbc=structure.pbc,
sites=[PeriodicSite(id="G", symbol="H", fractional=(0.0, 0.0, 0.0))],
)
class _Canvas:
"""Collects SVG elements in screen angstrom; ``u`` is the size unit for strokes and text."""
def __init__(self, right: np.ndarray, up: np.ndarray, toward: np.ndarray, u: float) -> None:
self.right, self.up, self.toward, self.u = right, up, toward, u
self.drawn: list[tuple[float, str]] = []
self.overlays: list[str] = []
self.texts: list[tuple[float, float, str, str, float]] = []
self.xs: list[float] = []
self.ys: list[float] = []
def xy(self, p: np.ndarray, *, extent: bool = True) -> tuple[float, float]:
x, y = float(p @ self.right), float(p @ self.up)
if extent:
self.xs.append(x)
self.ys.append(y)
return x, y
def depth(self, p: np.ndarray) -> float:
return float(p @ self.toward)
def line(
self,
a: np.ndarray,
b: np.ndarray,
colour: str,
width: float,
*,
dash: str | None = None,
arrow: bool = False,
depth: float | None = None,
) -> None:
(x1, y1), (x2, y2) = self.xy(a), self.xy(b)
style = f' stroke-dasharray="{dash}"' if dash else ""
svg = (
f'<line x1="{_n(x1)}" y1="{_n(y1)}" x2="{_n(x2)}" y2="{_n(y2)}" '
f'stroke="{colour}" stroke-width="{_n(width * self.u)}"{style}/>'
)
if arrow:
svg += _arrowhead(x1, y1, x2, y2, 0.5 * self.u, colour)
if depth is None:
self.overlays.append(svg)
else:
self.drawn.append((depth, svg))
def text(self, x: float, y: float, text: str, anchor: str = "start", size: float = 0.9) -> None:
self.texts.append((x, y, text, anchor, size * self.u))
def _n(value: float) -> str:
return f"{value:.3f}"
def _arrowhead(x1: float, y1: float, x2: float, y2: float, size: float, colour: str) -> str:
angle = math.atan2(y2 - y1, x2 - x1)
left = (x2 - size * math.cos(angle - 0.4), y2 - size * math.sin(angle - 0.4))
right = (x2 - size * math.cos(angle + 0.4), y2 - size * math.sin(angle + 0.4))
return (
f'<polygon points="{_n(x2)},{_n(y2)} {_n(left[0])},{_n(left[1])} '
f'{_n(right[0])},{_n(right[1])}" fill="{colour}"/>'
)
def _unit_for(points: list[tuple[float, float]]) -> float:
"""One size unit: a 26th of the drawing's larger extent (at least 4 angstrom)."""
xs = [p[0] for p in points] or [0.0]
ys = [p[1] for p in points] or [0.0]
return max(max(xs) - min(xs), max(ys) - min(ys), 4.0) / 26
def _finish(
canvas: _Canvas, header: list[tuple[str, str]], footer: list[str]
) -> tuple[str, float, float]:
"""Panel SVG in screen units (y down) with a legend row above and notes below."""
u = canvas.u
margin = 1.6 * u
minx, maxx = min(canvas.xs), max(canvas.xs)
miny, maxy = min(canvas.ys), max(canvas.ys)
top = 1.9 * u if header else 0.0
width = maxx - minx + 2 * margin
body_h = maxy - miny + 2 * margin
parts = [
f'<g transform="translate({_n(margin - minx)} {_n(top + margin + maxy)}) scale(1 -1)">'
]
parts.extend(svg for _d, svg in sorted(canvas.drawn, key=lambda item: item[0]))
parts.extend(canvas.overlays)
parts.append("</g>")
for x, y, text, anchor, size in canvas.texts:
parts.append(
f'<text x="{_n(x - minx + margin)}" y="{_n(top + maxy - y + margin)}" '
f'font-size="{_n(size)}" text-anchor="{anchor}">{html.escape(text)}</text>'
)
x = margin
for symbol, colour in header:
parts.append(
f'<circle cx="{_n(x + 0.45 * u)}" cy="{_n(0.9 * u)}" r="{_n(0.45 * u)}" '
f'fill="{colour}" stroke="#000" stroke-width="{_n(0.05 * u)}"/>'
f'<text x="{_n(x + 1.1 * u)}" y="{_n(1.25 * u)}" font-size="{_n(0.9 * u)}">'
f"{html.escape(symbol)}</text>"
)
x += (1.9 + 0.6 * len(symbol)) * u
for i, line in enumerate(footer):
parts.append(
f'<text x="{_n(margin)}" y="{_n(top + body_h + 1.2 * u * (i + 1))}" '
f'font-size="{_n(0.8 * u)}">{html.escape(line)}</text>'
)
width = max(
width,
(x + margin) if header else 0.0,
max((0.45 * u * len(t) + 2 * margin for t in footer), default=0.0),
)
height = top + body_h + (1.2 * u * len(footer) + 0.5 * u if footer else 0.0)
return "".join(parts), width, height
[docs]
def periodic_figure_svg(
structure: PeriodicStructure, options: FigureOptions | None = None
) -> bytes:
"""The figure as standalone SVG (UTF-8)."""
options = options or FigureOptions()
direction = view_direction(structure, options)
right, up = camera(direction)
panels = [_draw_structure(structure, options, right, up, -direction)]
if options.brillouin:
panels.append(_draw_zone(structure, right, up, -direction))
# Panels are scaled to a common height before they are put side by side.
target = max(h for _c, _w, h in panels)
scaled = [(c, w * target / h, h * target / h, target / h) for c, w, h in panels]
gap = 0.05 * target
total_w = sum(w for _c, w, _h, _k in scaled) + gap * (len(scaled) - 1)
scale = options.width / total_w
body = []
x = 0.0
for content, w, _h, k in scaled:
body.append(
f'<g transform="translate({_n(x * scale)} 0) scale({k * scale:.6f})">{content}</g>'
)
x += w + gap
title = html.escape(options.title or structure.id)
w_px, h_px = total_w * scale, target * scale
svg = (
'<?xml version="1.0" encoding="UTF-8"?>\n'
f'<svg xmlns="http://www.w3.org/2000/svg" width="{_n(w_px)}" height="{_n(h_px)}" '
f'viewBox="0 0 {_n(w_px)} {_n(h_px)}" font-family="Arial, Helvetica, sans-serif">'
f'<title>{title}</title><rect width="100%" height="100%" fill="#fff"/>{"".join(body)}</svg>'
)
return svg.encode("utf-8")
def _draw_structure(
structure: PeriodicStructure,
options: FigureOptions,
right: np.ndarray,
up: np.ndarray,
toward: np.ndarray,
) -> tuple[str, float, float]:
rows = crystal.lattice_rows(structure)
atoms = _images(structure, options)
positions = np.asarray([p for _, p in atoms])
periodic = [i for i in range(3) if structure.pbc[i]]
starts, ends = _edges(rows, periodic)
extent = [(float(p @ right), float(p @ up)) for p in [*positions, *starts, *ends]]
canvas = _Canvas(right, up, toward, _unit_for(extent))
u = canvas.u
radii = {s: _radius(s) for s in {site.symbol for site, _ in atoms}}
if options.bonds and len(atoms) <= 6000:
r = np.asarray([radii[site.symbol] for site, _ in atoms])
dist = np.linalg.norm(positions[:, None, :] - positions[None, :, :], axis=-1)
limit = 1.15 * (r[:, None] + r[None, :])
for i, j in zip(*np.nonzero(np.triu((dist < limit) & (dist > 0.4), k=1)), strict=True):
depth = (canvas.depth(positions[i]) + canvas.depth(positions[j])) / 2 - 1e-3
canvas.line(positions[i], positions[j], "#555", 0.12, depth=depth)
from vibesketcher.render import schematic
depths = np.asarray([canvas.depth(p) for p in positions])
span = float(depths.max() - depths.min()) if len(depths) else 0.0
record = crystal._builder(structure)
layer_of = record.get("layer_of_site") if (options.layers or options.stacking) else None
if (options.layers or options.stacking) and not isinstance(layer_of, dict):
raise FigureRefused("layer schematics need a slab made by the slab builder")
fixed = set(record.get("fixed_sites") or []) if options.layers else set()
letters = schematic.stacking_letters(structure) if options.stacking else {}
unknown = [e for e in options.polyhedra if e not in radii]
if unknown:
raise FigureRefused(f"no {', '.join(unknown)} atoms to draw polyhedra around")
for site, p in atoms:
if site.symbol not in options.polyhedra:
continue
shell = schematic.coordination(structure, p, site.symbol, {
**radii, **{s.symbol: _radius(s.symbol) for s in structure.sites}}) # fmt: skip
fill = _COLOURS.get(site.symbol, "#B0B0B0")
for face in schematic.hull_faces(shell):
corners = shell[face]
pts = " ".join(f"{_n(x)},{_n(y)}" for x, y in (canvas.xy(c) for c in corners))
canvas.drawn.append((
float(np.mean([canvas.depth(c) for c in corners])) - 1e-4,
f'<polygon class="polyhedron" points="{pts}" fill="{fill}" fill-opacity="0.35" '
f'stroke="#333" stroke-width="{_n(0.05 * u)}" stroke-linejoin="round"/>',
)) # fmt: skip
wyckoff = _wyckoff_labels(structure, options.symprec) if options.labels == "wyckoff" else {}
for index, (site, p) in enumerate(atoms):
x, y = canvas.xy(p)
rad = radii[site.symbol] * options.atom_scale
fill = _COLOURS.get(site.symbol, "#B0B0B0")
if options.depth_cue and span > 1e-6:
near = (depths[index] - depths.min()) / span # 1 nearest the viewer
rad *= 0.7 + 0.3 * near
fill = _paler(fill, 0.55 * (1 - near))
canvas.xs.extend([x - rad, x + rad])
canvas.ys.extend([y - rad, y + rad])
dash = f' stroke-dasharray="{_n(0.15 * u)} {_n(0.1 * u)}"' if site.id in fixed else ""
canvas.drawn.append(
(
depths[index],
f'<circle cx="{_n(x)}" cy="{_n(y)}" r="{_n(rad)}" '
f'data-site="{html.escape(site.id)}" fill="{fill}" stroke="#000" '
f'stroke-width="{_n(0.05 * u)}"{dash}/>',
)
)
if options.labels != "none":
canvas.text(
x,
y - rad - 0.9 * u,
{"element": site.symbol, "site": site.id}.get(options.labels)
or wyckoff.get(site.id, "?"),
"middle",
0.7,
)
height_note = None
if options.heights:
labels, height_note = schematic.heights(structure, positions, toward)
spots: dict[tuple[int, int], tuple[float, float, float, list[str]]] = {}
for (site, p), label in zip(atoms, labels, strict=True):
x, y = canvas.xy(p, extent=False)
key = (round(x / 0.05), round(y / 0.05))
rad = radii[site.symbol] * options.atom_scale
entry = spots.setdefault(key, (x, y, rad, []))
if label not in entry[3]:
entry[3].append(label)
for x, y, rad, values in spots.values():
text = ", ".join(sorted(values, key=_height_key))
canvas.text(x + rad + 0.15 * u, y - 0.25 * u, text, size=0.6)
if options.layers or options.stacking:
assert isinstance(layer_of, dict)
_draw_layers(structure, atoms, canvas, layer_of, fixed, letters, radii, options)
symmetry_note: list[str] = []
if options.symmetry_elements:
symmetry_note = _draw_symmetry(structure, canvas, rows, toward, options)
if options.cell:
for a, b in zip(starts, ends, strict=True):
canvas.line(a, b, "#000", 0.08)
if options.vectors:
placed: list[tuple[float, float]] = []
for axis, name in zip(range(3), "abc", strict=True):
if not structure.pbc[axis]:
continue
(x1, y1), (x2, y2) = canvas.xy(np.zeros(3)), canvas.xy(rows[axis])
if math.hypot(x2 - x1, y2 - y1) < 0.3 * u:
continue # seen end-on
canvas.line(np.zeros(3), rows[axis], "#000", 0.14, arrow=True)
# Keep labels of vectors that project onto each other apart.
ty = y2 + 0.3 * u
while any(abs(px - x2) < 6 * u and abs(py - ty) < 1.1 * u for px, py in placed):
ty -= 1.2 * u
placed.append((x2, ty))
canvas.text(x2 + 0.4 * u, ty, f"{name} = {float(np.linalg.norm(rows[axis])):.3f} Å")
if options.plane is not None:
if structure.pbc != (True, True, True):
raise FigureRefused("Miller planes are drawn in three-dimensional cells")
poly = plane_polygon(structure, options.plane)
pts = " ".join(f"{_n(x)},{_n(y)}" for x, y in map(canvas.xy, poly))
canvas.overlays.append(
f'<polygon points="{pts}" fill="#4dabf7" fill-opacity="0.35" '
f'stroke="#1864ab" stroke-width="{_n(0.08 * u)}"/>'
)
h, k, l_ = options.plane
cx, cy = canvas.xy(poly.mean(axis=0))
canvas.text(cx, cy, f"({h} {k} {l_})", "middle")
if options.wigner_seitz:
ws = crystal.wigner_seitz_polyhedron(structure)
verts = np.asarray(ws["vertices"], dtype=float)
for i, j in _polyhedron_edges(ws):
canvas.line(verts[i], verts[j], "#c92a2a", 0.08, dash=f"{_n(0.3 * u)} {_n(0.2 * u)}")
if options.scale_bar:
span = max(canvas.xs) - min(canvas.xs)
bar = next(b for b in (20.0, 10.0, 5.0, 2.0, 1.0) if b <= max(span / 3, 1.0))
x0, y0 = min(canvas.xs), min(canvas.ys) - 1.2 * u
canvas.ys.append(y0 - 1.0 * u)
canvas.overlays.append(
f'<line x1="{_n(x0)}" y1="{_n(y0)}" x2="{_n(x0 + bar)}" '
f'y2="{_n(y0)}" stroke="#000" stroke-width="{_n(0.15 * u)}"/>'
)
canvas.text(x0 + bar / 2, y0 - 1.0 * u, f"{bar:g} Å", "middle", 0.8)
footer = [f"{structure.id}, view {_view_text(options)}"]
if height_note:
footer.append(height_note)
footer.extend(symmetry_note)
if options.polyhedra:
footer.append("shaded: coordination polyhedra of " + ", ".join(options.polyhedra))
if options.layers and fixed:
footer.append("dashed outlines: fixed atoms")
if options.layers and isinstance(record.get("vacuum"), (int, float)):
footer.append(f"vacuum {float(record['vacuum']):g} Å between periodic images")
if options.wigner_seitz:
footer.append("dashed red: Wigner–Seitz cell of the lattice point at the origin")
if options.repeat != (1, 1, 1):
n1, n2, n3 = options.repeat
footer.append(f"display repeat {n1}×{n2}×{n3}, drawing only")
header = [(s, _COLOURS.get(s, "#B0B0B0")) for s in sorted(radii)] if options.legend else []
return _finish(canvas, header, footer)
def _wyckoff_labels(structure: PeriodicStructure, symprec: float) -> dict[str, str]:
"""Multiplicity and Wyckoff letter of each site in the standard setting (spglib)."""
from vibesketcher.document import symmetry
if not symmetry.available():
raise FigureRefused("Wyckoff labels need the symmetry extra (spglib)")
try:
report = symmetry.detect(structure, symprec=symprec)
except ValueError as error:
raise FigureRefused(str(error)) from None
lattice_points = abs(float(np.linalg.det(np.asarray(report["transformation_matrix"]))))
counts: dict[str, int] = {}
for info in report["sites"]:
counts[info["equivalent_to"]] = counts.get(info["equivalent_to"], 0) + 1
return {
site.id: f"{round(counts[info['equivalent_to']] / lattice_points)}{info['wyckoff']}"
for site, info in zip(structure.sites, report["sites"], strict=True)
}
_GLYPH_SIDES = {2: 0, 3: 3, 4: 4, 6: 6}
def _glyph(x: float, y: float, order: int, size: float, *, filled: bool) -> str:
fill = "#000" if filled else "#fff"
stroke = f'stroke="#000" stroke-width="{_n(size * 0.18)}"'
if order == 2:
return (f'<ellipse cx="{_n(x)}" cy="{_n(y)}" rx="{_n(size * 0.55)}" ry="{_n(size * 0.3)}" '
f'fill="{fill}" {stroke}/>') # fmt: skip
sides = _GLYPH_SIDES[order]
turn = math.pi / 2 if sides == 3 else math.pi / 4 if sides == 4 else 0.0
corners = " ".join(
f"{_n(x + size * 0.55 * math.cos(turn + 2 * math.pi * k / sides))},"
f"{_n(y + size * 0.55 * math.sin(turn + 2 * math.pi * k / sides))}"
for k in range(sides)
)
return f'<polygon points="{corners}" fill="{fill}" {stroke}/>'
def _draw_symmetry(
structure: PeriodicStructure,
canvas: _Canvas,
rows: np.ndarray,
toward: np.ndarray,
options: FigureOptions,
) -> list[str]:
"""Symmetry elements along the projection axis, in International Tables style."""
from vibesketcher.render import schematic, symmetry_elements
d = toward / np.linalg.norm(toward)
parallel = [
a for a in range(3) if abs(abs(float(rows[a] @ d)) / np.linalg.norm(rows[a]) - 1) < 1e-6
]
if structure.pbc != (True, True, True) or not parallel:
raise FigureRefused("symmetry elements are drawn in views along a, b or c of a 3D crystal")
axis = parallel[0]
i, j = [a for a in range(3) if a != axis]
try:
found, summary = symmetry_elements.elements(structure, axis, symprec=options.symprec)
except ValueError as error:
raise FigureRefused(str(error)) from None
u = canvas.u
def at(point: tuple[float, float]) -> tuple[float, float]:
return canvas.xy(point[0] * rows[i] + point[1] * rows[j], extent=False)
styles = {
"m": "",
"dashed": f' stroke-dasharray="{_n(0.5 * u)} {_n(0.3 * u)}"',
"dotted": f' stroke-dasharray="{_n(0.12 * u)} {_n(0.25 * u)}"',
"dashdot": (
f' stroke-dasharray="{_n(0.6 * u)} {_n(0.25 * u)} {_n(0.12 * u)} {_n(0.25 * u)}"'
), # fmt: skip
}
view_letter = "abc"[axis]
corner = at((0.0, 1.0))
stacked = 0
# Rotoinversion rings first, so an axis glyph on the same line sits inside its ring.
found = sorted(found, key=lambda e: e.kind != "rotoinversion")
for e in found:
if e.kind == "plane":
style = ("m" if e.symbol == "m" else "dotted" if e.symbol == view_letter
else "dashed" if e.symbol in ("a", "b", "c", "e") else "dashdot") # fmt: skip
(x1, y1), (x2, y2) = at(e.points[0]), at(e.points[1])
width = 0.16 if e.symbol == "m" else 0.1
canvas.overlays.append(
f'<line class="symmetry-plane" data-symbol="{e.symbol}" '
f'x1="{_n(x1)}" y1="{_n(y1)}" x2="{_n(x2)}" y2="{_n(y2)}" stroke="#000" '
f'stroke-width="{_n(width * u)}"{styles[style]}/>'
)
elif e.kind == "axis":
x, y = at(e.points[0])
order, _, screw = e.symbol.partition("_")
canvas.overlays.append(
f'<g class="symmetry-axis" data-symbol="{e.symbol}">'
f"{_glyph(x, y, int(order), 1.1 * u, filled=True)}</g>"
)
if screw:
canvas.text(x + 0.6 * u, y - 0.75 * u, screw, size=0.75)
elif e.kind == "rotoinversion":
x, y = at(e.points[0])
bar = int(e.symbol.lstrip("-"))
canvas.overlays.append(
f'<g class="symmetry-rotoinversion" data-symbol="{e.symbol}">'
f"{_glyph(x, y, bar, 1.6 * u, filled=False)}"
f'<circle cx="{_n(x)}" cy="{_n(y)}" r="{_n(0.13 * u)}" fill="#000"/></g>'
)
if e.height:
canvas.text(x + 0.95 * u, y + 0.25 * u, schematic.fraction_text(e.height), size=0.7)
elif e.kind == "inversion":
x, y = at(e.points[0])
canvas.overlays.append(
f'<circle class="symmetry-inversion" cx="{_n(x)}" cy="{_n(y)}" r="{_n(0.24 * u)}" '
f'fill="#fff" stroke="#000" stroke-width="{_n(0.05 * u)}"/>'
)
if e.height:
canvas.text(x + 0.3 * u, y + 0.3 * u, schematic.fraction_text(e.height), size=0.6)
elif e.kind == "horizontal_plane":
x, y = corner[0] - 7.5 * u, corner[1] - 1.3 * stacked * u
stacked += 1
canvas.overlays.append(
f'<polyline class="symmetry-horizontal" data-symbol="{e.symbol}" '
f'points="{_n(x)},{_n(y - 0.6 * u)} {_n(x)},{_n(y)} {_n(x + 0.6 * u)},{_n(y)}" '
f'fill="none" stroke="#000" stroke-width="{_n(0.1 * u)}"/>'
)
label = f"{e.symbol} at {schematic.fraction_text(e.height or 0)}"
canvas.text(x + 0.8 * u, y - 0.3 * u, label, size=0.75)
canvas.xs.append(x - 0.5 * u)
canvas.ys.append(y + 0.5 * u)
elif e.kind == "inplane_axis" and options.inplane_axes:
(x1, y1), (x2, y2) = at(e.points[0]), at(e.points[1])
canvas.overlays.append(
f'<line class="symmetry-inplane-axis" data-symbol="{e.symbol}" x1="{_n(x1)}" '
f'y1="{_n(y1)}" x2="{_n(x2)}" y2="{_n(y2)}" stroke="#555" '
f'stroke-width="{_n(0.05 * u)}"/>'
+ _arrowhead(x1, y1, x2, y2, 0.4 * u, "#555")
+ ("" if e.symbol == "2_1" else _arrowhead(x2, y2, x1, y1, 0.4 * u, "#555"))
)
notes = [f"symmetry elements of {summary['source']}, projected along {view_letter}",
"lines: mirror solid, glide in the plane of the drawing dashed, along the view "
"dotted, diagonal dash-dot; open circles: inversion centres"] # fmt: skip
if summary["oblique_operations_not_drawn"]:
notes.append(f"{summary['oblique_operations_not_drawn']} operations with elements "
"oblique to the view are not drawn") # fmt: skip
return notes
def _paler(colour: str, amount: float) -> str:
"""``colour`` mixed with white by ``amount`` (0 to 1)."""
rgb = [int(colour[i : i + 2], 16) for i in (1, 3, 5)]
return "#" + "".join(f"{round(c + (255 - c) * amount):02x}" for c in rgb)
def _height_key(text: str) -> float:
from fractions import Fraction
from vibesketcher.render.schematic import _GLYPHS
glyphs = {v: float(k) for k, v in _GLYPHS.items()}
try:
return glyphs[text] if text in glyphs else float(Fraction(text))
except (ValueError, ZeroDivisionError):
return 0.0
def _draw_layers(
structure: PeriodicStructure,
atoms: list[tuple[PeriodicSite, np.ndarray]],
canvas: _Canvas,
layer_of: dict[str, object],
fixed: set[object],
letters: dict[str, str],
radii: dict[str, float],
options: FigureOptions,
) -> None:
"""Layer numbers and fixed marks at the right, stacking letters at the left."""
rows = crystal.lattice_rows(structure)
normal = np.cross(rows[0], rows[1])
if abs(float(canvas.toward @ normal) / float(np.linalg.norm(normal))) > 0.1:
raise FigureRefused("layer schematics need a side view of the slab (view side or front)")
u = canvas.u
# Outside everything drawn so far: atoms, bonds and polyhedra.
left, right = min(canvas.xs) - 1.2 * u, max(canvas.xs) + 1.2 * u
by_layer: dict[int, list[float]] = {}
by_letter: dict[float, str] = {}
for site, p in atoms:
if site.id not in layer_of:
continue # adsorbate atoms
_x, y = canvas.xy(p, extent=False)
by_layer.setdefault(int(str(layer_of[site.id])), []).append(y)
if site.id in letters:
by_letter.setdefault(round(y, 1), letters[site.id])
for layer, ys in sorted(by_layer.items()) if options.layers else []:
ids = [s.id for s, _p in atoms if layer_of.get(s.id) == layer]
mark = " (fixed)" if ids and all(i in fixed for i in ids) else ""
lo, hi = min(ys), max(ys)
canvas.overlays.append(
f'<line x1="{_n(right)}" y1="{_n(lo)}" x2="{_n(right)}" y2="{_n(hi + 1e-3)}" '
f'stroke="#000" stroke-width="{_n(0.08 * u)}"/>'
)
canvas.text(right + 0.4 * u, (lo + hi) / 2 - 0.3 * u, f"layer {layer + 1}{mark}", size=0.7)
canvas.xs.append(right + 6 * u)
# Planes closer than a line of text share one label, e.g. "C/A".
merged: list[tuple[float, list[str]]] = []
for y, letter in sorted(by_letter.items()):
if merged and y - merged[-1][0] < 1.2 * u:
if letter not in merged[-1][1]:
merged[-1][1].append(letter)
else:
merged.append((y, [letter]))
for y, group in merged:
canvas.text(left - 0.3 * u, y - 0.3 * u, "/".join(group), "end", 0.75)
canvas.xs.append(left - 2.5 * u)
def _view_text(options: FigureOptions) -> str:
if options.view == "uvw" and options.direction is not None:
return "along [" + " ".join(str(x) for x in options.direction) + "]"
if options.view == "hkl" and options.direction is not None:
return "along the (" + " ".join(str(x) for x in options.direction) + ") normal"
return options.view
def _edges(rows: np.ndarray, periodic: list[int]) -> tuple[list[np.ndarray], list[np.ndarray]]:
vecs = [rows[i] for i in periodic]
n = len(vecs)
corners = [
sum((vecs[k] for k in range(n) if mask & (1 << k)), np.zeros(3)) for mask in range(1 << n)
]
starts, ends = [], []
for mask in range(1 << n):
for k in range(n):
if not mask & (1 << k):
starts.append(corners[mask])
ends.append(corners[mask | (1 << k)])
return starts, ends
def _polyhedron_edges(ws: dict[str, object]) -> set[tuple[int, int]]:
faces = ws.get("faces") or []
vertices = ws.get("vertices") or []
assert isinstance(faces, list) and isinstance(vertices, list)
edges: set[tuple[int, int]] = set()
for face in faces:
cycle = list(face)
for a, b in zip(cycle, cycle[1:] + cycle[:1], strict=True):
edges.add((min(a, b), max(a, b)))
if ws.get("dimension") == 1 and len(vertices) == 2:
edges.add((0, 1))
return edges
def _draw_zone(
structure: PeriodicStructure, right: np.ndarray, up: np.ndarray, toward: np.ndarray
) -> tuple[str, float, float]:
"""The first Brillouin zone: the Wigner-Seitz cell of the reciprocal lattice (1/angstrom)."""
reciprocal = reciprocal_structure(structure)
ws = crystal.wigner_seitz_polyhedron(reciprocal)
verts = np.asarray(ws["vertices"], dtype=float)
rows = crystal.lattice_rows(reciprocal)
tips = [rows[a] * 0.5 for a in range(3) if structure.pbc[a]]
extent = [(float(p @ right), float(p @ up)) for p in [*verts, *tips, np.zeros(3)]]
canvas = _Canvas(right, up, toward, _unit_for(extent) * 0.8)
u = canvas.u
for i, j in _polyhedron_edges(ws):
canvas.line(verts[i], verts[j], "#000", 0.08)
for axis, name in zip(range(3), ("b1", "b2", "b3"), strict=True):
if not structure.pbc[axis]:
continue
tip = rows[axis] * 0.5
(x1, y1), (x2, y2) = canvas.xy(np.zeros(3)), canvas.xy(tip)
if math.hypot(x2 - x1, y2 - y1) < 0.3 * u:
continue
canvas.line(np.zeros(3), tip, "#1864ab", 0.12, arrow=True)
canvas.text(
x2 + 0.3 * u,
y2 + 0.3 * u,
f"{name}/2, |{name}| = {float(np.linalg.norm(rows[axis])):.3f} 1/Å",
size=0.75,
)
gx, gy = canvas.xy(np.zeros(3))
canvas.overlays.append(f'<circle cx="{_n(gx)}" cy="{_n(gy)}" r="{_n(0.15 * u)}" fill="#000"/>')
canvas.text(gx - 0.3 * u, gy - 1.0 * u, "Γ", "end")
return _finish(canvas, [], ["First Brillouin zone, b_i · a_j = 2π δ_ij"])
[docs]
def periodic_figure(
structure: PeriodicStructure,
options: FigureOptions | None = None,
fmt: str = "svg",
*,
dpi: int = 300,
) -> bytes:
"""SVG, or PDF/PNG through the publication extra."""
svg = periodic_figure_svg(structure, options)
if fmt == "svg":
return svg
if fmt not in ("pdf", "png"):
raise FigureRefused(f"cannot write {fmt!r}; formats: svg, pdf, png")
from vibesketcher.render.publication import convert
return bytes(convert(svg, fmt, dpi=dpi))