scene

A rotatable collection of spheres and edges, decoupled from chemistry.

Everything below the Scene boundary is geometry: centres, radii, colours and a camera. Nothing here knows what a crystal is, which is what makes it portable to another language later if the frame rate ever demands it.

Scene dataclass

Source code in chmpy/tui/scene.py

@dataclass
class Scene:
    centers: np.ndarray  # (N, 3)
    radii: np.ndarray  # (N,)
    colors: np.ndarray  # (N, 3)
    bonds: list = field(default_factory=list)  # [(p0, p1, rgb0, rgb1, radius)]
    edges: list = field(default_factory=list)  # [(p0, p1)] in world coords
    axes: list = field(default_factory=list)  # [(start, end, rgb)] cell vectors
    origin: np.ndarray = None  # rotation pivot

    def __post_init__(self):
        if self.origin is None:
            self.origin = (
                self.centers.mean(axis=0) if len(self.centers) else np.zeros(3)
            )

    def __len__(self):
        return len(self.centers)

    def arrow_size(self):
        """(shaft, head) radii for the cell vectors, in world units."""
        span = float(np.linalg.norm(self.centers - self.origin, axis=1).max()) if len(
            self.centers
        ) else 1.0
        shaft = max(0.09, span * 0.014)
        return shaft, shaft * 2.4

    def geometry(self):
        """Every point the frame has to contain, with its padding radius.

        A bond or an arrow is a rod of finite thickness, so its endpoints need
        the same padding an atom gets; a cell edge is a hairline and needs
        none.
        """
        points = [self.centers]
        radii = [self.radii]

        def add(endpoints, radius):
            points.append(np.asarray(endpoints, dtype=float))
            radii.append(np.full(2, radius))

        for start, end, _rgb0, _rgb1, radius in self.bonds:
            add((start, end), radius)
        for start, end in self.edges:
            add((start, end), 0.0)
        head = self.arrow_size()[1]
        for start, end, _rgb in self.axes:
            add((start, end), head)
        return np.vstack(points), np.concatenate(radii)

    def view_extent(self, rotation, margin=1.04):
        """Half-width and half-height needed about the pivot, for this view.

        Measured about the rotation pivot rather than the projected centroid,
        so the structure does not drift across the frame as it turns. The
        margin keeps the outermost atom off the frame edge, where an exact
        fit would leave it sitting in the border row.
        """
        points, radii = self.geometry()
        p = (points - self.origin) @ rotation.T
        reach = np.abs(p[:, :2]) + radii[:, None]
        ex, ey = reach.max(axis=0) * margin
        return max(float(ex), 1e-3), max(float(ey), 1e-3)

arrow_size()

(shaft, head) radii for the cell vectors, in world units.

Source code in chmpy/tui/scene.py

def arrow_size(self):
    """(shaft, head) radii for the cell vectors, in world units."""
    span = float(np.linalg.norm(self.centers - self.origin, axis=1).max()) if len(
        self.centers
    ) else 1.0
    shaft = max(0.09, span * 0.014)
    return shaft, shaft * 2.4

geometry()

Every point the frame has to contain, with its padding radius.

A bond or an arrow is a rod of finite thickness, so its endpoints need the same padding an atom gets; a cell edge is a hairline and needs none.

Source code in chmpy/tui/scene.py

def geometry(self):
    """Every point the frame has to contain, with its padding radius.

    A bond or an arrow is a rod of finite thickness, so its endpoints need
    the same padding an atom gets; a cell edge is a hairline and needs
    none.
    """
    points = [self.centers]
    radii = [self.radii]

    def add(endpoints, radius):
        points.append(np.asarray(endpoints, dtype=float))
        radii.append(np.full(2, radius))

    for start, end, _rgb0, _rgb1, radius in self.bonds:
        add((start, end), radius)
    for start, end in self.edges:
        add((start, end), 0.0)
    head = self.arrow_size()[1]
    for start, end, _rgb in self.axes:
        add((start, end), head)
    return np.vstack(points), np.concatenate(radii)

view_extent(rotation, margin=1.04)

Half-width and half-height needed about the pivot, for this view.

Measured about the rotation pivot rather than the projected centroid, so the structure does not drift across the frame as it turns. The margin keeps the outermost atom off the frame edge, where an exact fit would leave it sitting in the border row.

Source code in chmpy/tui/scene.py

def view_extent(self, rotation, margin=1.04):
    """Half-width and half-height needed about the pivot, for this view.

    Measured about the rotation pivot rather than the projected centroid,
    so the structure does not drift across the frame as it turns. The
    margin keeps the outermost atom off the frame edge, where an exact
    fit would leave it sitting in the border row.
    """
    points, radii = self.geometry()
    p = (points - self.origin) @ rotation.T
    reach = np.abs(p[:, :2]) + radii[:, None]
    ex, ey = reach.max(axis=0) * margin
    return max(float(ex), 1e-3), max(float(ey), 1e-3)

bonding_percolates(crystal, tolerance=0.4)

True when covalent bonding runs through the lattice.

Each component is unwrapped by walking a spanning tree and accumulating the cell translation on every edge. If a remaining edge then disagrees with the shifts that walk assigned, the component closes a loop through a lattice translation - it is a framework or a chain, not a finite molecule.

Source code in chmpy/tui/scene.py

def bonding_percolates(crystal, tolerance=0.4):
    """True when covalent bonding runs through the lattice.

    Each component is unwrapped by walking a spanning tree and accumulating the
    cell translation on every edge. If a remaining edge then disagrees with the
    shifts that walk assigned, the component closes a loop through a lattice
    translation - it is a framework or a chain, not a finite molecule.
    """
    graph, edge_cells = crystal.unit_cell_connectivity(tolerance=tolerance)
    n_comp, labels = csgraph.connected_components(csgraph=graph, directed=False)
    shifts = np.zeros((graph.shape[0], 3))
    for comp in range(n_comp):
        root = int(np.flatnonzero(labels == comp)[0])
        order, pred = csgraph.breadth_first_order(
            csgraph=graph, i_start=root, directed=False
        )
        for j in order[1:]:
            i = pred[j]
            shifts[j] = (
                shifts[i] - edge_cells[(j, i)]
                if j < i
                else shifts[i] + edge_cells[(i, j)]
            )
    return any(
        not np.allclose(shifts[j] - shifts[i], cell)
        for (i, j), cell in edge_cells.items()
    )

crystal_scene(crystal, cells=(1, 1, 1), style='ball-and-stick')

Spheres and bonds for a block of unit cells, plus the cell edges.

A molecular crystal is built from whole molecules: chopping atoms at the cell boundary would draw half a molecule here and half a cell away, which is not what the structure looks like. A framework or ionic solid has no finite molecule to complete, so it falls back to the atoms of the block with bonds found by distance.

Source code in chmpy/tui/scene.py

def crystal_scene(crystal, cells=(1, 1, 1), style="ball-and-stick"):
    """Spheres and bonds for a block of unit cells, plus the cell edges.

    A molecular crystal is built from whole molecules: chopping atoms at the
    cell boundary would draw half a molecule here and half a cell away, which
    is not what the structure looks like. A framework or ionic solid has no
    finite molecule to complete, so it falls back to the atoms of the block
    with bonds found by distance.
    """
    hi = tuple(int(c) for c in cells)
    molecular = not bonding_percolates(crystal)

    if molecular:
        pos, z, bonds = _molecular_block(crystal, hi, style)
    else:
        pos, z, bonds = _network_block(crystal, hi, style)

    if style == "ball-and-stick":
        radii = np.array([Element[int(n)].ball_stick_radius for n in z]) * 1.1
    else:
        radii = np.array([Element[int(n)].vdw_radius for n in z])

    centre = crystal.to_cartesian(np.array([[0.5, 0.5, 0.5]]) * np.array(hi))[0]
    scene = Scene(
        pos,
        radii,
        _rgb(z),
        edges=_cell_edges(crystal, hi),
        axes=_cell_axes(crystal),
        origin=centre,
    )
    if style == "ball-and-stick":
        scene.bonds = bonds
    return scene

format_direction(indices, kind)

The crystallographic spelling of a parsed direction.

Multi-digit indices are spaced out, since running them together would produce something this module could not read back.

Source code in chmpy/tui/scene.py

def format_direction(indices, kind):
    """The crystallographic spelling of a parsed direction.

    Multi-digit indices are spaced out, since running them together would
    produce something this module could not read back.
    """
    joiner = " " if max(abs(i) for i in indices) > 9 else ""
    body = joiner.join(str(i) for i in indices)
    if kind == "hkl":
        return f"({body})"
    if kind == "cartesian":
        return "xyz"[int(np.argmax(np.abs(indices)))]
    return f"[{body}]"

lattice_direction(crystal, indices, kind='uvw')

Cartesian vector for a [uvw] zone axis or an (hkl) plane normal.

[uvw] is a direction in the direct lattice, ua + vb + wc. (hkl) names a plane, and the direction of interest is its normal, which is the reciprocal-lattice vector ha + kb + lc*. For a non-orthogonal cell these are genuinely different directions, so the distinction matters.

Source code in chmpy/tui/scene.py

def lattice_direction(crystal, indices, kind="uvw"):
    """Cartesian vector for a [uvw] zone axis or an (hkl) plane normal.

    [uvw] is a direction in the direct lattice, u*a + v*b + w*c. (hkl) names a
    plane, and the direction of interest is its normal, which is the
    reciprocal-lattice vector h*a* + k*b* + l*c*. For a non-orthogonal cell
    these are genuinely different directions, so the distinction matters.
    """
    indices = np.asarray(indices, dtype=float)
    cell = crystal.unit_cell
    if kind == "hkl":
        return indices @ cell.reciprocal_lattice
    return indices @ cell.direct

look_along(direction, up_hint=None)

Rotation that puts direction out of the screen, toward the viewer.

The returned matrix has the screen basis as its rows, matching what render_scene expects. up_hint fixes the roll; the component of it perpendicular to the view direction becomes screen-up.

Source code in chmpy/tui/scene.py

def look_along(direction, up_hint=None):
    """Rotation that puts `direction` out of the screen, toward the viewer.

    The returned matrix has the screen basis as its rows, matching what
    `render_scene` expects. `up_hint` fixes the roll; the component of it
    perpendicular to the view direction becomes screen-up.
    """
    w = np.asarray(direction, dtype=float)
    norm = np.linalg.norm(w)
    if norm < 1e-12:
        raise ValueError("view direction must be non-zero")
    w = w / norm

    if up_hint is None:
        # whichever cardinal axis is least parallel to the view direction
        up_hint = np.eye(3)[int(np.argmin(np.abs(w)))]
    up_hint = np.asarray(up_hint, dtype=float)

    right = np.cross(up_hint, w)
    if np.linalg.norm(right) < 1e-9:  # hint was parallel; pick another
        up_hint = np.eye(3)[(int(np.argmin(np.abs(w))) + 1) % 3]
        right = np.cross(up_hint, w)
    right /= np.linalg.norm(right)
    return np.stack([right, np.cross(w, right), w])

molecule_scene(molecule, style='ball-and-stick')

Spheres, and bonds when the style calls for them.

Source code in chmpy/tui/scene.py

def molecule_scene(molecule, style="ball-and-stick"):
    """Spheres, and bonds when the style calls for them."""
    # the longest dimension of the molecule ends up across the screen
    pos = molecule.positions_in_molecular_axis_frame()
    z = molecule.atomic_numbers

    if style == "space-filling":
        radii = np.array([Element[int(n)].vdw_radius for n in z])
        return Scene(pos, radii, _rgb(z), origin=np.zeros(3))

    if molecule.bonds is None:
        molecule.guess_bonds()
    radii = np.array([Element[int(n)].ball_stick_radius for n in z]) * 0.9
    return Scene(
        pos,
        radii,
        _rgb(z),
        bonds=_bond_segments(pos, z, molecule.unique_bonds),
        origin=np.zeros(3),
    )

parse_direction(text)

Parse a view direction into (indices, kind).

Accepts crystallographic notation - '[uvw]' for a zone axis, '(hkl)' for a plane normal, bare digits defaulting to a zone axis - as well as the named axes a/b/c, their reciprocals a/b/c*, and the cartesian x/y/z.

Indices may be run together in the usual way, with a minus sign binding to the digit that follows ('1-10' is [1,-1,0]), or separated by spaces or commas when they need more than one digit ('10 -2 1').

Source code in chmpy/tui/scene.py

def parse_direction(text):
    """Parse a view direction into (indices, kind).

    Accepts crystallographic notation - '[uvw]' for a zone axis, '(hkl)' for a
    plane normal, bare digits defaulting to a zone axis - as well as the named
    axes a/b/c, their reciprocals a*/b*/c*, and the cartesian x/y/z.

    Indices may be run together in the usual way, with a minus sign binding to
    the digit that follows ('1-10' is [1,-1,0]), or separated by spaces or
    commas when they need more than one digit ('10 -2 1').
    """
    raw = text.strip().lower()
    if not raw:
        raise ValueError("no direction given")
    if raw in NAMED:
        return NAMED[raw]

    kind = "uvw"
    if raw.startswith("(") and raw.endswith(")"):
        kind, raw = "hkl", raw[1:-1]
    elif raw.startswith("[") and raw.endswith("]"):
        kind, raw = "uvw", raw[1:-1]
    elif raw.startswith("{") and raw.endswith("}"):
        kind, raw = "hkl", raw[1:-1]
    raw = raw.strip()

    if any(sep in raw for sep in " ,;"):
        parts = [p for p in re.split(r"[\s,;]+", raw) if p]
    else:
        # run-together form: a minus sign belongs to the digit after it
        parts = re.findall(r"-?\d", raw)
        if "".join(parts) != raw.replace(" ", ""):
            raise ValueError(f"cannot read direction {text!r}")

    if len(parts) != 3:
        raise ValueError(f"need three indices, got {len(parts)} from {text!r}")
    try:
        indices = tuple(int(p) for p in parts)
    except ValueError as exc:
        raise ValueError(f"cannot read direction {text!r}") from exc
    if not any(indices):
        raise ValueError("direction must be non-zero")
    return indices, kind

render_scene(scene, rotation=None, width=88, zoom=1.0, edge_color=None, depth_cue=True, height=None, theme=None, shading='lit', extent=None)

Rasterise a Scene into a Pixels image under the given rotation.

Source code in chmpy/tui/scene.py

def render_scene(
    scene,
    rotation=None,
    width=88,
    zoom=1.0,
    edge_color=None,
    depth_cue=True,
    height=None,
    theme=None,
    shading="lit",
    extent=None,
):
    """Rasterise a Scene into a Pixels image under the given rotation."""
    from .canvas import theme_for

    theme = theme or theme_for(None)
    edge_color = edge_color or theme.edge
    r = np.eye(3) if rotation is None else rotation
    pts = (scene.centers - scene.origin) @ r.T

    height = height or width
    # default to a tight fit for this view; callers that rotate interactively
    # pass their own extent so the framing does not breathe frame to frame
    ex, ey = extent if extent is not None else scene.view_extent(r)
    z = max(zoom, 1e-3)
    ex, ey = ex / z, ey / z
    # both dimensions have to fit, or a frame wider than it is tall - which
    # every terminal is - crops the structure top and bottom
    scale = min((width - 1) / (2 * ex), (height - 1) / (2 * ey))
    origin = (-0.5 * width / scale, 0.5 * height / scale)

    px = Pixels(height, width)
    cue = None
    if depth_cue and len(pts):
        cue = (pts[:, 2].max(), pts[:, 2].min() - 1.0)

    shared = {"depth_cue": cue, "fade_to": theme.background,
              "shading": shading}
    zbuffer = draw_spheres(
        px, pts, scene.radii, scene.colors, scale, origin, **shared
    )
    if scene.bonds:
        starts, ends, c0, c1, radii = _bond_arrays(scene, r)
        draw_cylinders(
            px, starts, ends, radii, c0, c1, scale, origin,
            zbuffer=zbuffer, **shared,
        )

    for p0, p1 in scene.edges:
        a = (np.asarray(p0) - scene.origin) @ r.T
        b = (np.asarray(p1) - scene.origin) @ r.T
        _line3(px, a, b, scale, origin, edge_color, zbuffer=zbuffer)

    if scene.axes:
        _draw_axes(px, scene, r, scale, origin, zbuffer, shared)
    return px

rotation_matrix(yaw, pitch, roll=0.0)

Rotation from yaw (about screen y), pitch (about screen x), roll (z).

Source code in chmpy/tui/scene.py

def rotation_matrix(yaw, pitch, roll=0.0):
    """Rotation from yaw (about screen y), pitch (about screen x), roll (z)."""
    cy, sy = np.cos(yaw), np.sin(yaw)
    cp, sp = np.cos(pitch), np.sin(pitch)
    cr, sr = np.cos(roll), np.sin(roll)
    ry = np.array([[cy, 0, sy], [0, 1, 0], [-sy, 0, cy]])
    rx = np.array([[1, 0, 0], [0, cp, -sp], [0, sp, cp]])
    rz = np.array([[cr, -sr, 0], [sr, cr, 0], [0, 0, 1]])
    return rz @ rx @ ry

view_rotation(indices, kind='uvw', crystal=None)

Rotation looking along a direction, resolved against a lattice if given.

For a crystal the roll is fixed by whichever cell vector is closest to perpendicular to the view, so an axis-aligned view comes out the way it would be drawn on paper rather than at an arbitrary twist.

Source code in chmpy/tui/scene.py

def view_rotation(indices, kind="uvw", crystal=None):
    """Rotation looking along a direction, resolved against a lattice if given.

    For a crystal the roll is fixed by whichever cell vector is closest to
    perpendicular to the view, so an axis-aligned view comes out the way it
    would be drawn on paper rather than at an arbitrary twist.
    """
    if crystal is None or kind == "cartesian":
        return look_along(indices)
    w = lattice_direction(crystal, indices, kind)
    direct = crystal.unit_cell.direct
    cosine = np.abs(direct @ w) / (
        np.linalg.norm(direct, axis=1) * np.linalg.norm(w)
    )
    return look_along(w, direct[int(np.argmin(cosine))])