vibe-sketcherDocumentation
Source code for vibesketcher.chem.periodic_geometry
"""Shared geometry for built periodic structures: neighbours, bond images, records.
Used by the periodic nanotube and framework builders. Cells are right-handed
row vectors in angstrom; a bond's image is the integer translation applied to
its second site, as in :class:`~vibesketcher.document.model.PeriodicStructure`.
"""
from __future__ import annotations
from collections.abc import Sequence
from typing import Any, Literal
import numpy as np
from numpy.typing import NDArray
from vibesketcher.document.model import PeriodicBond, PeriodicSite, PeriodicStructure
from vibesketcher.limits import MAX_ATOMS
__all__ = ["BUILDER_VERSION", "BuildRefused", "bond_image", "neighbour_pairs"]
#: Version of the nanotube and framework constructions, recorded in provenance.
BUILDER_VERSION = 1
_Points = NDArray[np.float64]
[docs]
class BuildRefused(ValueError):
"""Parameters a builder cannot build; the message says why."""
[docs]
def neighbour_pairs(
cartesian: _Points, lattice: _Points, pbc: tuple[bool, bool, bool], cutoff: float
) -> list[tuple[int, int, tuple[int, int, int], float]]:
"""All ``(i, j, image, distance)`` with ``|r_j + image·L - r_i| < cutoff``, i < j.
A hashed cubic grid of edge ``cutoff`` (the linked-cell method; Allen and
Tildesley, *Computer Simulation of Liquids*, 2nd ed., 2017, 5.3). Images
are the 27 neighbouring translations along periodic axes, so every cell
vector must be longer than the cutoff.
"""
for axis in range(3):
if pbc[axis] and float(np.linalg.norm(lattice[axis])) <= cutoff:
raise BuildRefused("a cell vector is shorter than the neighbour cutoff")
grid: dict[tuple[int, int, int], list[int]] = {}
keys = np.floor(cartesian / cutoff).astype(int)
for index, key in enumerate(map(tuple, keys)):
grid.setdefault(key, []).append(index)
low = cartesian.min(axis=0) - cutoff
high = cartesian.max(axis=0) + cutoff
ranges = [(-1, 0, 1) if periodic else (0,) for periodic in pbc]
images = [(i, j, k) for i in ranges[0] for j in ranges[1] for k in ranges[2]]
pairs: list[tuple[int, int, tuple[int, int, int], float]] = []
for image in images:
shift = np.array(image, dtype=float) @ lattice
# A query point q = r_i - shift lies near r_j means r_j + shift lies near r_i.
for i in range(len(cartesian)):
q = cartesian[i] - shift
if np.any(q < low) or np.any(q > high):
continue
base = np.floor(q / cutoff).astype(int)
for dx in (-1, 0, 1):
for dy in (-1, 0, 1):
for dz in (-1, 0, 1):
for j in grid.get((base[0] + dx, base[1] + dy, base[2] + dz), ()):
if j < i or (j == i and image <= (0, 0, 0)):
continue
distance = float(np.linalg.norm(cartesian[j] + shift - cartesian[i]))
if distance < cutoff:
pairs.append((i, j, image, distance))
return pairs
[docs]
def bond_image(
first: _Points, second: _Points, lattice: _Points, pbc: tuple[bool, bool, bool]
) -> tuple[tuple[int, int, int], float]:
"""The translation of ``second`` nearest ``first`` among neighbouring images."""
best: tuple[tuple[int, int, int], float] | None = None
ranges = [(-1, 0, 1) if periodic else (0,) for periodic in pbc]
for i in ranges[0]:
for j in ranges[1]:
for k in ranges[2]:
image = (i, j, k)
distance = float(np.linalg.norm(second + np.array(image) @ lattice - first))
if best is None or distance < best[1] - 1e-12:
best = (image, distance)
assert best is not None
return best
def _check_collisions(
cartesian: _Points, lattice: _Points, pbc: tuple[bool, bool, bool], minimum: float = 0.7
) -> None:
close = neighbour_pairs(cartesian, lattice, pbc, minimum)
if close:
i, j, _, distance = close[0]
raise BuildRefused(
f"{len(close)} atom pair(s) closer than {minimum} angstrom "
f"(atoms {i + 1} and {j + 1}: {distance:.3f} angstrom); "
"the requested structure overlaps itself"
)
def _structure(
identifier: str,
lattice: _Points,
pbc: tuple[bool, bool, bool],
symbols: Sequence[str],
cartesian: _Points,
bonds: Sequence[tuple[int, int, float]],
*,
cell_setting: Literal["primitive", "conventional", "unspecified"],
provenance: dict[str, Any],
charge: int | None = None,
multiplicity: int | None = None,
occupancy: Sequence[float] | None = None,
labels: Sequence[str] | None = None,
) -> PeriodicStructure:
"""Wrap periodic axes into [0, 1), number sites per element, give bonds images."""
if len(symbols) > MAX_ATOMS:
raise BuildRefused(f"{len(symbols)} sites exceed the {MAX_ATOMS}-atom limit")
fractional = cartesian @ np.linalg.inv(lattice)
for axis in range(3):
if pbc[axis]:
fractional[:, axis] -= np.floor(fractional[:, axis])
fractional[:, axis] = np.where(
fractional[:, axis] >= 1 - 1e-12, 0.0, fractional[:, axis]
)
wrapped = fractional @ lattice
counts: dict[str, int] = {}
ids: list[str] = []
for index, symbol in enumerate(symbols):
if labels is not None:
ids.append(labels[index])
continue
counts[symbol] = counts.get(symbol, 0) + 1
ids.append(f"{symbol}{counts[symbol]}")
sites = [
PeriodicSite(
id=ids[k],
symbol=symbols[k],
fractional=(float(fractional[k, 0]), float(fractional[k, 1]), float(fractional[k, 2])),
occupancy=1 if occupancy is None else occupancy[k],
)
for k in range(len(symbols))
]
periodic_bonds = []
for i, j, order in bonds:
image, distance = bond_image(wrapped[i], wrapped[j], lattice, pbc)
if distance > 3.5:
raise BuildRefused(f"bond {ids[i]}-{ids[j]} spans {distance:.2f} angstrom")
periodic_bonds.append(PeriodicBond(a=ids[i], b=ids[j], image=image, order=order))
rows = tuple(tuple(float(x) for x in row) for row in lattice)
return PeriodicStructure(
id=identifier,
lattice_angstrom=rows,
pbc=pbc,
sites=sites,
bonds=periodic_bonds,
charge=charge,
multiplicity=multiplicity,
cell_setting=cell_setting,
provenance=provenance,
)
def _provenance(kind: str, recipe: Any, **extra: Any) -> dict[str, Any]:
"""Provenance in the solid-state builders' convention: ``{"builder": {...}}``."""
return {
"builder": {
"kind": kind,
"version": BUILDER_VERSION,
"recipe": recipe.model_dump(mode="json", exclude_none=True),
"status": "unrelaxed starting geometry",
**extra,
}
}