vibe-sketcherDocumentation
Source code for vibesketcher.chem.rolled_tubes
"""Nanotubes rolled from any two-dimensional periodic structure.
The rolling construction generalizes the carbon chiral vector (Saito,
Dresselhaus and Dresselhaus, *Physical Properties of Carbon Nanotubes*, 1998,
ch. 3) to any layer with in-plane lattice vectors a1, a2, as CRYSTAL does for
inorganic tubes (Noel, D'Arco, Demichelis, Zicovich-Wilson and Dovesi,
J. Comput. Chem. 31, 855, 2010):
* the rolling vector is R = n a1 + m a2 and becomes the circumference;
* the tube period is the shortest lattice vector T = t1 a1 + t2 a2 with
T . R = 0. It exists only when the lattice metric makes (a1 . R) / (a2 . R)
rational; otherwise the tube has no translational period and is refused;
* the rolled cell holds abs(n t2 - m t1) layer cells per period;
* a point at in-plane fraction s along R and height h above the layer's middle
plane goes to angle 2 pi s and radius norm(R) / (2 pi) + h. Atoms outside that
middle surface are stretched along the circumference and those inside are
compressed, as for any rolled multi-plane layer, so this is an unrelaxed
starting geometry.
The cell is one-dimensional: z is the tube axis, x and y hold the tube plus a
stated vacuum gap and are not periodic.
"""
from __future__ import annotations
import math
from fractions import Fraction
from typing import Literal
import numpy as np
from pydantic import BaseModel, ConfigDict, Field
from vibesketcher.chem.periodic_geometry import BuildRefused, _check_collisions, _provenance
from vibesketcher.chem.periodic_geometry import _structure as _record
from vibesketcher.document.model import PeriodicStructure
from vibesketcher.limits import MAX_ATOMS
__all__ = ["RolledNanotubeRecipe", "roll_layer"]
[docs]
class RolledNanotubeRecipe(BaseModel):
"""Roll a 2D periodic structure of the document along n a1 + m a2."""
model_config = ConfigDict(extra="forbid")
kind: Literal["rolled_nanotube"] = "rolled_nanotube"
identifier: str = Field(pattern=r"^[A-Za-z][A-Za-z0-9_-]{0,63}$")
source: str = Field(pattern=r"^[A-Za-z][A-Za-z0-9_-]{0,63}$")
n: int = Field(ge=-100, le=100)
m: int = Field(ge=-100, le=100)
periods: int = Field(default=1, ge=1, le=1000)
vacuum: float = Field(default=10.0, ge=3, le=200, allow_inf_nan=False)
#: Roll so the layer's +normal side faces the axis instead of outwards.
inward_normal: bool = False
#: ``lattice``: n, m count the layer's own a1, a2 (as CRYSTAL's NANOTUBE
#: input does). ``nanotube``: the carbon-nanotube convention for a
#: hexagonal layer, with a1 and a2 60 degrees apart, so (n,n) is armchair
#: and (n,0) zigzag whichever hexagonal setting the layer uses.
indices: Literal["lattice", "nanotube"] = "lattice"
#: The thinnest tube built, angstrom (as for carbon nanotubes).
_MIN_DIAMETER = 4.0
def _hexagonal_indices(a1: np.ndarray, a2: np.ndarray, n: int, m: int) -> tuple[int, int]:
"""(n, m) in the 60-degree hexagonal basis as indices of the layer's own a1, a2."""
l1, l2 = float(np.linalg.norm(a1)), float(np.linalg.norm(a2))
cosine = float(a1 @ a2) / (l1 * l2)
if abs(l1 - l2) > 1e-6 * l1 or min(abs(cosine - 0.5), abs(cosine + 0.5)) > 1e-6:
raise BuildRefused(
"nanotube indices need a hexagonal layer (|a1| = |a2|, 60 or 120 degrees); "
"use lattice indices"
)
# 120 degrees: the 60-degree second vector is a1 + a2.
return (n, m) if cosine > 0 else (n + m, m)
def _perpendicular(p: float, q: float) -> tuple[int, int]:
"""Smallest integers (t1, t2) with t1 p + t2 q = 0, or a refusal."""
scale = max(abs(p), abs(q))
if abs(q) <= 1e-10 * scale:
return (0, 1)
if abs(p) <= 1e-10 * scale:
return (1, 0)
ratio = Fraction(p / q).limit_denominator(1000)
if abs(float(ratio) - p / q) > 1e-8 * max(1.0, abs(p / q)):
raise BuildRefused(
"no lattice vector of the layer is perpendicular to the rolling vector, so the "
"tube would have no translational period; choose other indices"
)
# t1 (u / v) + t2 = 0 -> (t1, t2) = (v, -u)
return (ratio.denominator, -ratio.numerator)
[docs]
def roll_layer(recipe: RolledNanotubeRecipe, layer: PeriodicStructure) -> PeriodicStructure:
"""The 1D tube cell rolled from ``layer``, or a refusal with the reason."""
if tuple(layer.pbc) != (True, True, False):
raise BuildRefused(
"only a two-dimensional structure (periodic along a1 and a2, not a3) can be rolled"
)
if recipe.n == 0 and recipe.m == 0:
raise BuildRefused("the rolling vector (0,0) is zero")
lattice = np.asarray(layer.lattice_angstrom, dtype=float)
a1, a2 = lattice[0], lattice[1]
normal = np.cross(a1, a2)
normal /= np.linalg.norm(normal)
if recipe.inward_normal:
normal = -normal
n, m = recipe.n, recipe.m
if recipe.indices == "nanotube":
n, m = _hexagonal_indices(a1, a2, n, m)
rolling = n * a1 + m * a2
t1, t2 = _perpendicular(float(a1 @ rolling), float(a2 @ rolling))
cells = n * t2 - m * t1
if cells < 0:
t1, t2, cells = -t1, -t2, -cells
if max(abs(t1), abs(t2)) > 1000:
raise BuildRefused("the tube period would span more than 1000 layer cells")
count = cells * recipe.periods * len(layer.sites)
if count > MAX_ATOMS:
raise BuildRefused(f"{count} atoms exceed the {MAX_ATOMS}-atom limit")
period = recipe.periods * (t1 * a1 + t2 * a2)
# In-plane fractions (s along R, u along the whole period) of every layer image.
cartesian = np.array([layer.cartesian(site) for site in layer.sites], dtype=float)
heights = cartesian @ normal
middle = (heights.max() + heights.min()) / 2
basis = np.array([[n, m], [recipe.periods * t1, recipe.periods * t2]], float)
inverse = np.linalg.inv(basis) # layer-cell fractions -> (s, u)
fractions = np.array([site.fractional[:2] for site in layer.sites], dtype=float)
# The four corners of the rolled cell in layer-cell units: 0, R, T, R + T.
corners = np.array([[0.0, 0.0], basis[0], basis[1], basis[0] + basis[1]])
low = np.floor(corners.min(axis=0) - 1).astype(int)
high = np.ceil(corners.max(axis=0) + 1).astype(int)
symbols: list[str] = []
points: list[list[float]] = []
parents: list[str] = []
eps = 1e-9
radius = float(np.linalg.norm(rolling)) / (2 * math.pi)
if 2 * radius < _MIN_DIAMETER:
raise BuildRefused(
f"the rolled tube would be {2 * radius / 10:.3f} nm across; tubes thinner than "
f"{_MIN_DIAMETER / 10} nm are not built"
)
for i in range(low[0], high[0] + 1):
for j in range(low[1], high[1] + 1):
for index, site in enumerate(layer.sites):
s, u = (fractions[index] + (i, j)) @ inverse
if not (-eps <= s < 1 - eps and -eps <= u < 1 - eps):
continue
r = radius + float(heights[index] - middle)
if r <= 0:
raise BuildRefused(
f"the tube radius {radius:.3f} angstrom is smaller than the layer's "
"half thickness; choose a larger rolling vector"
)
angle = 2 * math.pi * s
symbols.append(site.symbol)
points.append([r * math.cos(angle), r * math.sin(angle), u])
parents.append(site.id)
if len(symbols) != count:
raise AssertionError("rolled cell lost or duplicated layer images")
length = float(np.linalg.norm(period))
xyz = np.array(points)
outer = float(np.hypot(xyz[:, 0], xyz[:, 1]).max())
side = 2 * outer + recipe.vacuum
xyz[:, 0] += side / 2
xyz[:, 1] += side / 2
xyz[:, 2] *= length
tube = np.diag([side, side, length])
pbc = (False, False, True)
_check_collisions(xyz, tube, pbc)
structure = _record(
recipe.identifier,
tube,
pbc,
symbols,
xyz,
[],
cell_setting="unspecified",
charge=None if layer.charge is None else layer.charge * cells * recipe.periods,
provenance=_provenance(
"rolled_nanotube",
recipe,
translation=[t1, t2],
layer_cells_per_period=cells,
radius_angstrom=round(radius, 6),
axis="z",
warnings=[
"atoms off the layer's middle surface are stretched (outside) or compressed "
"(inside) along the circumference; relax before analysis"
],
),
)
return structure.model_copy(
update={
"sites": [
site.model_copy(update={"parent_site": parent})
for site, parent in zip(structure.sites, parents, strict=True)
]
}
)