"""Carbon nanotubes, graphene and fullerenes generated from their parameters (#137).
A nanostructure is made, not drawn: its atoms come from the honeycomb
lattice in Cartesian coordinates (angstrom), so a tube of thousands of atoms
never goes through 2D layout and the ``MAX_LAYOUT_ATOMS`` bound does not
apply; :data:`vibesketch.limits.MAX_ATOMS` does, checked before any atom is
made.
Nanotubes follow the chiral-vector construction (Saito, Dresselhaus and
Dresselhaus, *Physical Properties of Carbon Nanotubes*, Imperial College
Press 1998, ch. 3; Hamada, Sawada and Oshiyama, Phys. Rev. Lett. 68, 1579,
1992). With graphene's lattice vectors a1, a2 (``|a1| = |a2| = a = sqrt(3) a_CC``)
the chiral vector is Ch = n a1 + m a2, and
* the diameter is d = a sqrt(n^2 + nm + m^2) / pi;
* the chiral angle theta, between Ch and a1, has
cos(theta) = (2n + m) / (2 sqrt(n^2 + nm + m^2)), 0 <= theta <= 30 degrees;
* the translation vector is T = t1 a1 + t2 a2 with t1 = (2m + n)/d_R,
t2 = -(2n + m)/d_R, d_R = gcd(2m + n, 2n + m), and ``|T| = sqrt(3) |Ch| / d_R``;
* one translation period holds N = 2(n^2 + nm + m^2)/d_R hexagons, 2N carbons.
The lattice points inside the rectangle Ch x (periods T) are rolled onto a
cylinder of radius ``|Ch| / (2 pi)``. Walls of a multi-walled tube are concentric
and should be 0.34 nm apart, the graphite interlayer spacing.
Fullerenes are the truncated icosahedron C60 (Kroto et al., Nature 318, 162,
1985), from the even permutations of (0, +-1, +-3 phi), (+-1, +-(2 + phi),
+-2 phi) and (+-phi, +-2, +-phi^3). Seen along a five-fold axis C60 is two
30-atom hemispheres ending in rings of ten, which are the rings of a (5,5)
tube; C70 inserts one ring at the equator (D5h), and a capped (5,5) tube
inserts two rings per period. Stone-Wales defects turn one C-C bond by 90
degrees (Stone and Wales, Chem. Phys. Lett. 128, 501, 1986).
Coordinates are a starting geometry: bond lengths are the lattice's
a_CC = 1.42 angstrom, groups take standard bond lengths, and sp3 carbons stay
where the lattice put them. Optimise before reading geometry from them.
"""
from __future__ import annotations
import hashlib
import itertools
import json
import math
from collections import Counter
from collections.abc import Callable, Iterable, Sequence
from dataclasses import dataclass, field
from typing import Any
import numpy as np
from numpy.typing import NDArray
from vibesketch.document.model import (
Chirality,
Defect,
Fullerene,
Graphene,
NanoGroup,
Nanoribbon,
Nanostructure,
NanostructureParameters,
Nanotube,
PeriodicCell,
SidewallGroup,
)
__all__ = [
"A_CC",
"INTERLAYER_NM",
"AtomModel",
"NanostructureRefused",
"build",
"chiral_angle_deg",
"concentric_walls",
"diameter_nm",
"hexagons_per_period",
"hill_formula",
"identity_of",
"is_metallic",
"translation_angstrom",
"unrolled",
"wall_properties",
]
#: The C-C distance in graphene, angstrom.
A_CC = 1.42
#: The lattice constant a = sqrt(3) a_CC, angstrom.
A_LATTICE = math.sqrt(3) * A_CC
#: Graphite's interlayer spacing, the wall spacing of a multi-walled tube, nm.
INTERLAYER_NM = 0.34
#: How far from 0.34 nm two walls may be (Saito et al. report 0.34 to 0.39 nm).
_WALL_GAP_NM = (0.30, 0.40)
#: The thinnest tube generated; smaller ones are not stable sp2 cylinders.
_MIN_DIAMETER_NM = 0.4
#: Generator version, part of the identity: a change of construction is a new identity.
_VERSION = 1
#: Bonds are found between carbons closer than this, angstrom.
_CC_CUTOFF = 1.7
_Points = NDArray[np.float64]
[docs]
class NanostructureRefused(ValueError):
"""Parameters this generator cannot build; the message says why."""
[docs]
@dataclass
class AtomModel:
"""A generated structure: symbols, positions in angstrom, bonds ``(i, j, order)``."""
symbols: list[str]
positions: list[list[float]]
bonds: list[tuple[int, int, float]]
description: str
periodic: PeriodicCell | None = None
#: Diameter and chiral angle of each wall, for a nanotube.
walls: list[dict[str, Any]] = field(default_factory=list)
[docs]
def record(self, parameters: NanostructureParameters) -> Nanostructure:
identity, key = identity_of(parameters)
return Nanostructure(
parameters=parameters,
identity=identity,
key=key,
description=self.description,
formula=self.formula,
bonds=self.bonds,
periodic=self.periodic,
)
# --- the chiral-vector formulas (Saito, Dresselhaus and Dresselhaus 1998, ch. 3) ---
def _norm2(c: Chirality) -> int:
return c.n * c.n + c.n * c.m + c.m * c.m
[docs]
def diameter_nm(c: Chirality) -> float:
"""d = a sqrt(n^2 + nm + m^2) / pi, in nm."""
return A_LATTICE * math.sqrt(_norm2(c)) / math.pi / 10
[docs]
def chiral_angle_deg(c: Chirality) -> float:
"""theta from cos(theta) = (2n + m) / (2 sqrt(n^2 + nm + m^2)); 0 zigzag, 30 armchair."""
cosine = (2 * c.n + c.m) / (2 * math.sqrt(_norm2(c)))
return math.degrees(math.acos(min(1.0, cosine)))
[docs]
def wall_properties(c: Chirality) -> dict[str, Any]:
"""Diameter, chiral angle, translation period and electronic character of one wall."""
return {
"n": c.n,
"m": c.m,
"diameter_nm": round(diameter_nm(c), 4),
"chiral_angle_deg": round(chiral_angle_deg(c), 3),
"translation_nm": round(translation_angstrom(c) / 10, 4),
"electronic": "metallic" if is_metallic(c) else "semiconducting",
}
def _d_r(c: Chirality) -> int:
return math.gcd(2 * c.m + c.n, 2 * c.n + c.m)
[docs]
def translation_angstrom(c: Chirality) -> float:
"""``|T| = sqrt(3) |Ch| / d_R``, angstrom."""
return math.sqrt(3) * A_LATTICE * math.sqrt(_norm2(c)) / _d_r(c)
[docs]
def hexagons_per_period(c: Chirality) -> int:
"""N = 2 (n^2 + nm + m^2) / d_R; a period holds 2N carbons."""
return 2 * _norm2(c) // _d_r(c)
[docs]
def concentric_walls(inner: Chirality, count: int) -> list[Chirality]:
"""``count`` walls about ``inner``, each about 0.34 nm outside the last, of its type.
Armchair walls step (n + 5, n + 5) and zigzag walls (n + 9, 0), the
nearest whole steps to 0.34 nm; a chiral wall is scaled to the radius and
rounded, so its chiral angle stays close to the inner wall's.
"""
if count < 1:
raise NanostructureRefused("a tube has at least one wall")
walls = [inner]
for _ in range(count - 1):
last = walls[-1]
if last.n == last.m:
walls.append(Chirality(n=last.n + 5, m=last.m + 5))
elif last.m == 0:
walls.append(Chirality(n=last.n + 9, m=0))
else:
scale = (diameter_nm(last) + 2 * INTERLAYER_NM) / diameter_nm(last)
n, m = round(last.n * scale), round(last.m * scale)
walls.append(Chirality(n=min(n, 100), m=min(m, n, 100)))
return walls
# --- identity ---
def _canonical(parameters: NanostructureParameters) -> dict[str, Any]:
"""The parameters with a tube's length resolved to whole periods per wall."""
data = parameters.model_dump(mode="json")
if isinstance(parameters, Nanotube):
data.pop("length_nm", None)
data.pop("periods", None)
data["periods"] = _wall_periods(parameters)
data["ends"] = list(parameters.ends)
return data
[docs]
def identity_of(parameters: NanostructureParameters) -> tuple[str, str]:
"""The canonical identity of a nanostructure and its key, ``NANO-`` and 20 hex digits.
It depends on the parameters alone, a tube's length as whole periods of
each wall, so the same structure asked for in nm or in periods has one
identity. It replaces the InChI, which has no notation for these.
"""
identity = json.dumps(
{"generator": _VERSION, **_canonical(parameters)}, sort_keys=True, separators=(",", ":")
)
digest = hashlib.sha256(identity.encode()).hexdigest()[:20].upper()
return identity, f"NANO-{digest}"
# --- lattice and graph helpers ---
def _bonds_within(points: _Points, cutoff: float = _CC_CUTOFF) -> list[tuple[int, int]]:
"""Pairs closer than ``cutoff``, by a cell list: linear in the number of points."""
cells: dict[tuple[int, int, int], list[int]] = {}
keys = np.floor(points / cutoff).astype(int)
for index, key in enumerate(map(tuple, keys)):
cells.setdefault(key, []).append(index)
pairs: list[tuple[int, int]] = []
limit = cutoff * cutoff
for (x, y, z), members in cells.items():
for dx, dy, dz in itertools.product((-1, 0, 1), repeat=3):
others = cells.get((x + dx, y + dy, z + dz))
if not others:
continue
for i in members:
for j in others:
if j > i:
d = points[i] - points[j]
if float(d @ d) < limit:
pairs.append((i, j))
return sorted(pairs)
def _neighbours(count: int, pairs: Iterable[tuple[int, int]]) -> list[set[int]]:
around: list[set[int]] = [set() for _ in range(count)]
for i, j in pairs:
around[i].add(j)
around[j].add(i)
return around
def _prune(points: _Points, pairs: list[tuple[int, int]]) -> tuple[_Points, list[tuple[int, int]]]:
"""Drop carbons with fewer than two carbon neighbours, repeatedly: no dangling chains."""
keep = np.ones(len(points), dtype=bool)
around = _neighbours(len(points), pairs)
changed = True
while changed:
changed = False
for i in range(len(points)):
if keep[i] and len(around[i]) < 2:
keep[i] = False
changed = True
for j in around[i]:
around[j].discard(i)
around[i] = set()
index = np.cumsum(keep) - 1
kept = [(int(index[i]), int(index[j])) for i, j in pairs if keep[i] and keep[j]]
return points[keep], kept
@dataclass
class _Skeleton:
"""A carbon framework being decorated: positions, bonds and each carbon's surface normal."""
points: _Points
pairs: list[tuple[int, int]]
normal: Callable[[_Points], _Points]
#: Map from a site to its carbon, for groups and defects.
locate: Callable[[Any], int] | None = None
sp3: set[int] = field(default_factory=set)
removed: set[int] = field(default_factory=set)
used: set[int] = field(default_factory=set)
groups: list[tuple[int, NanoGroup, bool]] = field(default_factory=list)
def around(self) -> list[set[int]]:
return _neighbours(len(self.points), self.pairs)
def _rotation(axis: _Points, angle: float) -> _Points:
axis = axis / np.linalg.norm(axis)
x, y, z = axis
k = np.array([[0.0, -z, y], [z, 0.0, -x], [-y, x, 0.0]])
turn: _Points = np.eye(3) + math.sin(angle) * k + (1 - math.cos(angle)) * (k @ k)
return turn
# --- groups ---
#: Bond lengths, angstrom: to an sp2 (edge) carbon, then to an sp3 (sidewall) one.
_C_H = (1.09, 1.10)
_C_O = (1.36, 1.43)
_C_C = (1.49, 1.52)
def _group_atoms(
group: NanoGroup, at: _Points, out: _Points, side: _Points, sp3: bool
) -> tuple[list[str], list[_Points], list[tuple[int, int, float]]]:
"""A group's atoms and bonds, ``-1`` standing for the carbon it is bonded to.
The group lies in the plane of ``out`` (away from the carbon) and ``side``
(perpendicular to it): the sheet's own plane at an edge, so a carboxyl is
conjugated with it.
"""
k = 1 if sp3 else 0
if group == "H":
return ["H"], [at + _C_H[k] * out], [(-1, 0, 1.0)]
if group == "OH":
o = at + _C_O[k] * out
tilt = math.radians(71)
h = o + 0.97 * (math.cos(tilt) * out + math.sin(tilt) * side)
return ["O", "H"], [o, h], [(-1, 0, 1.0), (0, 1, 1.0)]
c = at + _C_C[k] * out
up = math.cos(math.radians(60)) * out
across = math.sin(math.radians(60)) * side
o_double = c + 1.21 * (up + across)
o_single = c + 1.35 * (up - across)
h = o_single + 0.97 * out
return (
["C", "O", "O", "H"],
[c, o_double, o_single, h],
[(-1, 0, 1.0), (0, 1, 2.0), (0, 2, 1.0), (2, 3, 1.0)],
)
def _unit(vector: _Points) -> _Points:
length = float(np.linalg.norm(vector))
return vector / length if length > 1e-9 else vector
def _perpendicular(out: _Points, prefer: _Points) -> _Points:
"""The unit part of ``prefer`` perpendicular to ``out``; another axis when they are parallel."""
for axis in (prefer, np.array([0.0, 0.0, 1.0]), np.array([1.0, 0.0, 0.0])):
side = axis - float(axis @ out) * out
if float(np.linalg.norm(side)) > 0.3:
return _unit(side)
return _unit(np.cross(out, np.array([0.0, 1.0, 0.0])))
def _terminate(skeleton: _Skeleton, edge_group: Callable[[int], NanoGroup]) -> None:
"""Ask for a group on every two-coordinated carbon not removed or decorated already.
Carbons next to a vacancy are left as they are; :func:`_assemble` decides
whether a hydroxyl or carboxyl fits.
"""
around = skeleton.around()
taken = {i for i, _, _ in skeleton.groups}
for i in range(len(skeleton.points)):
if len(around[i]) == 2 and i not in skeleton.used and i not in taken:
skeleton.groups.append((i, edge_group(i), False))
#: The closest two atoms of different groups, or a group and the lattice, may come, angstrom.
_CONTACT = {True: 1.7, False: 2.3}
#: No atom of a group reaches further than this from its carbon, with contact, angstrom.
_REACH = 8.0
def _assemble(skeleton: _Skeleton, description: str) -> AtomModel:
"""Carbons, then each group's atoms, with removed carbons dropped and bond orders set.
Groups go where they touch nothing: sidewall groups first (one that
cannot is refused), then the edge hydroxyls and carboxyls, each in the
sheet's plane or turned about their bond in 30-degree steps to the twist
with the most room, checked against the lattice, the groups placed so far
and a hydrogen on every other edge carbon; one that fits no way becomes
hydrogen. Then the hydrogens.
"""
points = skeleton.points
keep = [i for i in range(len(points)) if i not in skeleton.removed]
new_index = {old: new for new, old in enumerate(keep)}
symbols = ["C"] * len(keep)
positions: list[list[float]] = [[round(float(v), 6) for v in points[i]] for i in keep]
bonds: list[tuple[int, int, float]] = []
for i, j in skeleton.pairs:
if i in new_index and j in new_index:
order = 1.0 if i in skeleton.sp3 or j in skeleton.sp3 else 1.5
bonds.append((new_index[i], new_index[j], order))
around = skeleton.around()
normals = skeleton.normal(points)
axis = np.array([0.0, 0.0, 1.0])
def frame(carbon: int, sp3: bool) -> tuple[_Points, _Points]:
"""The direction away from ``carbon`` and one perpendicular to it in the sheet's plane."""
normal = _unit(normals[carbon])
if sp3:
return normal, _perpendicular(normal, axis)
pull = sum((_unit(points[j] - points[carbon]) for j in around[carbon]), np.zeros(3))
out = _unit(-pull)
return out, _perpendicular(out, np.cross(normal, out))
edge_h = {
carbon: points[carbon] + _C_H[0] * frame(carbon, False)[0]
for carbon, _group, sp3 in skeleton.groups
if not sp3
}
placed: list[tuple[int, NanoGroup, list[str], list[_Points], list[tuple[int, int, float]]]] = []
occupied: list[tuple[_Points, bool]] = []
kept_points = points[keep]
def surroundings(carbon: int) -> tuple[_Points, _Points, list[tuple[_Points, bool]]]:
"""What a group on ``carbon`` could touch: lattice, edge hydrogens, placed atoms nearby."""
here = points[carbon]
near = np.linalg.norm(kept_points - here, axis=1) < _REACH
for i in (carbon, *around[carbon]):
if i in new_index:
near[new_index[i]] = False
hydrogens = np.array(
[
h
for c, h in edge_h.items()
if c != carbon and float(np.linalg.norm(h - here)) < _REACH
]
).reshape(-1, 3)
placed_near = [
(p, light) for p, light in occupied if float(np.linalg.norm(p - here)) < _REACH
]
return kept_points[near], hydrogens, placed_near
def clearance(
around_it: tuple[_Points, _Points, list[tuple[_Points, bool]]],
atoms: list[str],
where: list[_Points],
) -> float:
"""How far the closest contact of the group is beyond the allowed one; < 0 overlaps."""
lattice, hydrogens, placed_near = around_it
worst = math.inf
for symbol, at in zip(atoms, where, strict=True):
light = symbol == "H"
if len(lattice):
worst = min(
worst, float(np.linalg.norm(lattice - at, axis=1).min()) - _CONTACT[light]
)
if len(hydrogens):
worst = min(
worst, float(np.linalg.norm(hydrogens - at, axis=1).min()) - _CONTACT[True]
)
for p, other_light in placed_near:
worst = min(worst, float(np.linalg.norm(at - p)) - _CONTACT[light or other_light])
return worst
# Sidewall groups first, then the bulky edge groups, then the hydrogens.
for carbon, group, sp3 in sorted(skeleton.groups, key=lambda g: (not g[2], g[1] == "H")):
out, in_plane = frame(carbon, sp3)
chosen: tuple[list[str], list[_Points], list[tuple[int, int, float]]] | None = None
if group == "H" and not sp3:
chosen = _group_atoms("H", points[carbon], out, in_plane, False)
else:
# In the sheet's plane first, then twisted about the bond in 30-degree steps.
twists = [_rotation(out, math.radians(30 * k)) @ in_plane for k in range(12)]
nearby = surroundings(carbon)
scored = [
(clearance(nearby, *_group_atoms(group, points[carbon], out, side, sp3)[:2]), k)
for k, side in enumerate(twists)
]
best, k = max(scored, key=lambda item: (round(item[0], 6), -item[1]))
if best >= 0:
chosen = _group_atoms(group, points[carbon], out, twists[k], sp3)
if chosen is None and sp3:
raise NanostructureRefused(
f"a sidewall {group} does not fit at carbon {carbon}: it would touch "
"another group or the edge; choose a site further away"
)
if chosen is None:
group = "H"
chosen = _group_atoms("H", points[carbon], out, in_plane, False)
placed.append((carbon, group, *chosen))
if group != "H" or sp3:
occupied.extend((p, s == "H") for s, p in zip(chosen[0], chosen[1], strict=True))
edge_h.pop(carbon, None)
for carbon, _group, group_symbols, group_points, group_bonds in sorted(
placed, key=lambda g: g[0]
):
start = len(symbols)
symbols.extend(group_symbols)
positions.extend([[round(float(v), 6) for v in p] for p in group_points])
for a, b, bond_order in group_bonds:
first = new_index[carbon] if a == -1 else start + a
bonds.append((min(first, start + b), max(first, start + b), bond_order))
bonds.sort()
return AtomModel(symbols, positions, bonds, description)
# --- groups and defects at sites ---
def _modify(
skeleton: _Skeleton, sidewall: Sequence[SidewallGroup], defects: Sequence[Defect]
) -> None:
"""Place sidewall groups and defects at the carbons nearest their sites.
Each site is resolved on the undecorated lattice, so one site does not
move another; two that touch the same carbons are refused.
"""
if not sidewall and not defects:
return
if skeleton.locate is None:
raise NanostructureRefused("this structure takes no sidewall groups or defects")
around = skeleton.around()
def claim(atoms: set[int], what: str) -> None:
clash = atoms & skeleton.used
if clash:
raise NanostructureRefused(
f"{what} is too close to another group or defect (carbon {min(clash)}): "
"choose sites further apart"
)
skeleton.used |= atoms
for item in sidewall:
carbon = skeleton.locate(item.site)
claim({carbon}, f"the sidewall {item.group} at {item.site.model_dump(exclude_none=True)}")
skeleton.sp3.add(carbon)
skeleton.groups.append((carbon, item.group, True))
for defect in defects:
carbon = skeleton.locate(defect.site)
what = f"the {defect.kind} at {defect.site.model_dump(exclude_none=True)}"
if defect.kind == "vacancy":
claim({carbon, *around[carbon]}, what)
skeleton.removed.add(carbon)
else:
_stone_wales(skeleton, carbon, around, claim, what)
def _stone_wales(
skeleton: _Skeleton,
carbon: int,
around: list[set[int]],
claim: Callable[[set[int], str], None],
what: str,
) -> None:
"""Turn the bond from ``carbon`` most nearly along the length by 90 degrees about its middle.
The two carbons swap one neighbour each, which makes the pentagon-
heptagon pairs of the defect; positions stay on the surface.
"""
points = skeleton.points
along = np.array([0.0, 0.0, 1.0])
partner = max(
sorted(around[carbon]), key=lambda j: abs(float(_unit(points[j] - points[carbon]) @ along))
)
first, second = sorted((carbon, partner))
old = (around[first] | around[second]) - {first, second}
claim({first, second, *old}, what)
middle = (points[first] + points[second]) / 2
normal = _unit(skeleton.normal(middle[None, :])[0])
turn = _rotation(normal, math.pi / 2)
moved = {i: middle + turn @ (points[i] - middle) for i in (first, second)}
ends = sorted(old)
# Each turned carbon keeps the two old neighbours now nearest it.
by_first = sorted(ends, key=lambda j: float(np.linalg.norm(points[j] - moved[first])))
to_first = set(by_first[:2])
to_second = set(ends) - to_first
if len(to_second) != 2 or any(
float(np.linalg.norm(points[j] - moved[second])) > 2.0 for j in to_second
):
raise NanostructureRefused(f"{what} cannot be made at the edge; choose an inner site")
for i in (first, second):
points[i] = moved[i]
pairs = [p for p in skeleton.pairs if not (set(p) & {first, second})]
pairs.append((first, second))
pairs.extend((min(first, j), max(first, j)) for j in to_first)
pairs.extend((min(second, j), max(second, j)) for j in to_second)
skeleton.pairs = sorted(pairs)
# --- nanotubes ---
def _wall_periods(tube: Nanotube) -> list[int]:
"""Whole translation periods per wall: those of the length, or of the innermost wall's."""
if tube.periods is not None:
length = tube.periods * translation_angstrom(tube.walls[0])
else:
assert tube.length_nm is not None
length = tube.length_nm * 10
return [max(1, round(length / translation_angstrom(w))) for w in tube.walls]
def _wall(c: Chirality, periods: int) -> _Points:
"""The rolled lattice of one wall, ``periods`` translation periods long, along +z."""
a1 = A_LATTICE * np.array([math.sqrt(3) / 2, 0.5])
a2 = A_LATTICE * np.array([math.sqrt(3) / 2, -0.5])
dr = _d_r(c)
ch = c.n * a1 + c.m * a2
t = ((2 * c.m + c.n) // dr) * a1 - ((2 * c.n + c.m) // dr) * a2
span = periods * t
basis = np.column_stack([a1, a2])
corners = np.array([[0.0, 0.0], ch, span, ch + span])
lattice = np.linalg.solve(basis, corners.T).T
low = np.floor(lattice.min(axis=0)).astype(int) - 1
high = np.ceil(lattice.max(axis=0)).astype(int) + 1
i, j = np.meshgrid(np.arange(low[0], high[0] + 1), np.arange(low[1], high[1] + 1))
sites = (i.reshape(-1, 1) * a1 + j.reshape(-1, 1) * a2).astype(float)
sites = np.vstack([sites, sites + (a1 + a2) / 3])
eps = 1e-7
s = sites @ ch / float(ch @ ch)
u = sites @ span / float(span @ span)
inside = (s >= -eps) & (s < 1 - eps) & (u >= -eps) & (u < 1 - eps)
s, u = s[inside], u[inside]
radius = float(np.linalg.norm(ch)) / (2 * math.pi)
angle = 2 * math.pi * s
z = u * float(np.linalg.norm(span))
points = np.column_stack([radius * np.cos(angle), radius * np.sin(angle), z])
return points[np.lexsort((angle, np.round(z, 6)))]
[docs]
def unrolled(
c: Chirality, periods: int = 1, margin: float = 0.12
) -> tuple[_Points, list[tuple[int, int]]]:
"""The graphene sheet a wall is rolled from, Ch along +x and T along +y, in angstrom.
The lattice covers the rectangle Ch x (periods T) and ``margin`` of it
around; returns the points (z = 0) and the C-C bonds between them.
"""
a1 = A_LATTICE * np.array([math.sqrt(3) / 2, 0.5])
a2 = A_LATTICE * np.array([math.sqrt(3) / 2, -0.5])
dr = _d_r(c)
ch = c.n * a1 + c.m * a2
t = (((2 * c.m + c.n) // dr) * a1 - ((2 * c.n + c.m) // dr) * a2) * periods
basis = np.column_stack([a1, a2])
corners = np.array([[0.0, 0.0], ch, t, ch + t]) * (1 + 2 * margin) - margin * (ch + t)
lattice = np.linalg.solve(basis, corners.T).T
low = np.floor(lattice.min(axis=0)).astype(int) - 1
high = np.ceil(lattice.max(axis=0)).astype(int) + 1
i, j = np.meshgrid(np.arange(low[0], high[0] + 1), np.arange(low[1], high[1] + 1))
sites = (i.reshape(-1, 1) * a1 + j.reshape(-1, 1) * a2).astype(float)
sites = np.vstack([sites, sites + (a1 + a2) / 3])
s = sites @ ch / float(ch @ ch)
u = sites @ t / float(t @ t)
inside = (s >= -margin) & (s <= 1 + margin) & (u >= -margin) & (u <= 1 + margin)
flat = np.column_stack(
[s[inside] * float(np.linalg.norm(ch)), u[inside] * float(np.linalg.norm(t))]
)
points = np.column_stack([flat, np.zeros(len(flat))])
return points, _bonds_within(points)
def _check_tube(tube: Nanotube) -> list[int]:
from vibesketch.limits import check_atoms
for wall in tube.walls:
if diameter_nm(wall) < _MIN_DIAMETER_NM:
raise NanostructureRefused(
f"({wall.n},{wall.m}) is {diameter_nm(wall):.2f} nm across; tubes thinner "
f"than {_MIN_DIAMETER_NM} nm are not generated"
)
for inner, outer in itertools.pairwise(tube.walls):
gap = (diameter_nm(outer) - diameter_nm(inner)) / 2
if not _WALL_GAP_NM[0] <= gap <= _WALL_GAP_NM[1]:
raise NanostructureRefused(
f"walls ({inner.n},{inner.m}) and ({outer.n},{outer.m}) are {gap:.3f} nm "
f"apart; concentric walls must be {_WALL_GAP_NM[0]} to {_WALL_GAP_NM[1]} nm "
f"apart (graphite: {INTERLAYER_NM} nm), innermost first"
)
if "capped" in tube.ends and [(w.n, w.m) for w in tube.walls] != [(5, 5)]:
raise NanostructureRefused(
"a capped end is generated for a single-walled (5,5) tube only, which a C60 "
"hemisphere closes; other chiralities need a cap search vibe-sketch does not do"
)
periods = _wall_periods(tube)
carbons = sum(2 * hexagons_per_period(w) * p for w, p in zip(tube.walls, periods, strict=True))
check_atoms(carbons, "the nanotube")
return periods
def _tube(tube: Nanotube) -> AtomModel:
periods = _check_tube(tube)
if "capped" in tube.ends:
points, pairs = _capped_55(periods[0], tube.ends)
else:
walls: list[_Points] = []
middle = max(p * translation_angstrom(w) for w, p in zip(tube.walls, periods, strict=True))
for wall, count in zip(tube.walls, periods, strict=True):
points = _wall(wall, count)
points[:, 2] += (middle - count * translation_angstrom(wall)) / 2
walls.append(points)
points = np.vstack(walls)
pairs = _bonds_within(points)
points, pairs = _prune(points, pairs)
outer = tube.walls[-1]
outer_radius = diameter_nm(outer) * 5
z_low, z_high = float(points[:, 2].min()), float(points[:, 2].max())
def normal(p: _Points) -> _Points:
centre = np.zeros_like(p)
centre[:, 2] = np.clip(p[:, 2], z_low + outer_radius, z_high - outer_radius)
if "capped" not in tube.ends:
centre[:, 2] = p[:, 2]
return np.array([_unit(v) for v in p - centre])
skeleton = _Skeleton(points, pairs, normal)
radii = np.hypot(points[:, 0], points[:, 1])
on_outer = radii > outer_radius - 0.5 if len(tube.walls) > 1 else np.ones(len(points), bool)
around = skeleton.around()
candidates = [i for i in range(len(points)) if on_outer[i] and len(around[i]) == 3]
def locate(site: Any) -> int:
if site.around is None:
raise NanostructureRefused("a site on a nanotube gives 'around' in degrees")
z = z_low + site.along * (z_high - z_low)
angle = math.radians(site.around)
target = np.array([outer_radius * math.cos(angle), outer_radius * math.sin(angle), z])
return min(candidates, key=lambda i: float(np.linalg.norm(points[i] - target)))
skeleton.locate = locate
_modify(skeleton, tube.sidewall, tube.defects)
lower = (z_low + z_high) / 2
_terminate(
skeleton,
lambda i: _end_group(tube.ends[0] if points[i, 2] < lower else tube.ends[1]),
)
model = _assemble(skeleton, _tube_description(tube, periods))
model.walls = [
{**wall_properties(w), "periods": p} for w, p in zip(tube.walls, periods, strict=True)
]
if len(tube.walls) == 1 and "capped" not in tube.ends:
wall = tube.walls[0]
model.periodic = PeriodicCell(
translation_angstrom=round(translation_angstrom(wall), 6),
carbons_per_period=2 * hexagons_per_period(wall),
)
return model
def _end_group(end: str) -> NanoGroup:
if end == "OH":
return "OH"
if end == "COOH":
return "COOH"
return "H"
def _count(number: int, noun: str) -> str:
return f"{number} {noun}" + ("" if number == 1 else "s")
def _tube_description(tube: Nanotube, periods: list[int]) -> str:
walls = "@".join(f"({w.n},{w.m})" for w in tube.walls)
count = {1: "single", 2: "double", 3: "triple"}.get(len(tube.walls))
kind = f"{count}-walled" if count else f"{len(tube.walls)}-walled"
length = periods[0] * translation_angstrom(tube.walls[0]) / 10
ends = tube.ends[0] if tube.ends[0] == tube.ends[1] else "/".join(tube.ends)
finish = "capped" if ends == "capped" else f"{ends}-terminated"
extra = []
if tube.sidewall:
extra.append(_count(len(tube.sidewall), "sidewall group"))
if tube.defects:
extra.append(_count(len(tube.defects), "defect"))
tail = f", {', '.join(extra)}" if extra else ""
return f"{walls} {kind} carbon nanotube, {length:.2f} nm, {finish}{tail}"
# --- fullerenes and capped (5,5) tubes ---
def _c60_along_axis() -> _Points:
"""C60 with a five-fold axis along z, bonds a_CC long."""
phi = (1 + math.sqrt(5)) / 2
seeds = [(0.0, 1.0, 3 * phi), (1.0, 2 + phi, 2 * phi), (phi, 2.0, phi**3)]
found: set[tuple[float, ...]] = set()
for x, y, w in seeds:
for p in ((x, y, w), (y, w, x), (w, x, y)):
for signs in itertools.product((1, -1), repeat=3):
found.add(tuple(round(a * s, 9) for a, s in zip(p, signs, strict=True)))
points = np.array(sorted(found))
points *= A_CC / 2.0 # the truncated icosahedron of these seeds has edges of 2
axis = _unit(np.array([0.0, 1.0, phi]))
z = np.array([0.0, 0.0, 1.0])
turn = _rotation(np.cross(axis, z), math.acos(float(axis @ z)))
return points @ turn.T
def _rings_inserted(rings: int, open_top: bool) -> _Points:
"""C60 split at its equator with ``rings`` rings of ten inserted (Ih C60, D5h C70, ...).
The two equatorial rings of C60 are rings of a (5,5) tube, a/2 apart and
turned 36 degrees about the axis from each other; inserted rings
continue that pattern. ``open_top`` leaves out the upper hemisphere.
"""
c60 = _c60_along_axis()
z = np.round(c60[:, 2], 6)
levels = sorted(set(z))
lower = c60[z < 0]
upper = c60[z > 0]
ring = c60[z == levels[3]] # the lower equatorial ring
step = float(levels[4] - levels[3])
turn = _rotation(np.array([0.0, 0.0, 1.0]), math.radians(36))
parts = [lower]
for k in range(1, rings + 1):
copy = ring @ np.linalg.matrix_power(turn, k).T
copy[:, 2] += k * step
parts.append(copy)
if not open_top:
top = upper @ (turn.T if rings % 2 else np.eye(3))
top[:, 2] += rings * step
parts.append(top)
else:
parts.append(upper[np.round(upper[:, 2], 6) == levels[4]] @ np.eye(3))
parts[-1] = parts[-1] @ (turn.T if rings % 2 else np.eye(3))
parts[-1][:, 2] += rings * step
return np.vstack(parts)
def _capped_55(periods: int, ends: tuple[str, str]) -> tuple[_Points, list[tuple[int, int]]]:
"""A (5,5) tube closed by C60 hemispheres: C60 with two rings per period inserted."""
points = _rings_inserted(2 * periods, open_top=ends != ("capped", "capped"))
if ends[0] != "capped":
points[:, 2] *= -1 # the open end is the lower one
points = points[np.lexsort((np.arctan2(points[:, 1], points[:, 0]), points[:, 2]))]
pairs = _bonds_within(points)
return points, pairs
def _fullerene(f: Fullerene) -> AtomModel:
points = _rings_inserted(0 if f.size == 60 else 1, open_top=False)
pairs = _bonds_within(points)
around = _neighbours(len(points), pairs)
if len(points) != f.size or any(len(a) != 3 for a in around):
raise AssertionError(f"C{f.size} was not built as a closed three-connected cage")
skeleton = _Skeleton(points, pairs, lambda p: np.array([_unit(v) for v in p]))
name = "buckminsterfullerene C60 (Ih)" if f.size == 60 else "fullerene C70 (D5h)"
return _assemble(skeleton, name)
# --- graphene and ribbons ---
def _sheet_lattice(width: float, length: float, x0: float = 0.0, y0: float = 0.0) -> _Points:
"""Honeycomb points in [x0, x0 + width] x [y0, y0 + length]; zigzag rows along x."""
a1 = np.array([A_LATTICE, 0.0])
a2 = np.array([A_LATTICE / 2, A_LATTICE * math.sqrt(3) / 2])
rows = math.ceil((y0 + length) / a2[1]) + 2
columns = math.ceil((x0 + width) / A_LATTICE) + rows + 2
i, j = np.meshgrid(np.arange(-rows - 2, columns), np.arange(-2, rows))
sites = i.reshape(-1, 1) * a1 + j.reshape(-1, 1) * a2
sites = np.vstack([sites, sites + (a1 + a2) / 3])
eps = 1e-6
inside = (
(sites[:, 0] >= x0 - eps)
& (sites[:, 0] < x0 + width - eps)
& (sites[:, 1] >= y0 - eps)
& (sites[:, 1] <= y0 + length + eps)
)
flat = sites[inside]
flat = flat[np.lexsort((flat[:, 0], np.round(flat[:, 1], 6)))]
return np.column_stack([flat, np.zeros(len(flat))])
def _flat(
points: _Points, edges: NanoGroup, item: Graphene | Nanoribbon, description: str
) -> AtomModel:
from vibesketch.limits import check_atoms
check_atoms(len(points), "the graphene structure")
pairs = _bonds_within(points)
points, pairs = _prune(points, pairs)
if len(points) < 6:
raise NanostructureRefused("these dimensions hold no complete hexagon")
skeleton = _Skeleton(points, pairs, lambda p: np.tile([0.0, 0.0, 1.0], (len(p), 1)))
around = skeleton.around()
inner = [i for i in range(len(points)) if len(around[i]) == 3]
lows, highs = points.min(axis=0), points.max(axis=0)
def locate(site: Any) -> int:
if site.across is None:
raise NanostructureRefused("a site on a sheet or ribbon gives 'across', a fraction")
target = np.array(
[
lows[0] + site.across * (highs[0] - lows[0]),
lows[1] + site.along * (highs[1] - lows[1]),
0.0,
]
)
if not inner:
raise NanostructureRefused("this structure has no inner carbon for a site")
return min(inner, key=lambda i: float(np.linalg.norm(points[i] - target)))
skeleton.locate = locate
_modify(skeleton, item.sidewall, item.defects)
_terminate(skeleton, lambda _i: edges)
return _assemble(skeleton, description)
def _graphene(g: Graphene) -> AtomModel:
width, length = g.width_nm * 10, g.length_nm * 10
if width * length / 2.62 > 50_000: # 2.62 A^2 per carbon
from vibesketch.limits import check_atoms
check_atoms(round(width * length / 2.62), "the graphene sheet")
points = _sheet_lattice(width, length)
# Long side vertical: 'along' runs over the length.
description = f"graphene flake, {g.width_nm:g} x {g.length_nm:g} nm, {g.edges}-terminated"
return _flat(points, g.edges, g, description)
def _ribbon(r: Nanoribbon) -> AtomModel:
from vibesketch.limits import check_atoms
check_atoms(2 * r.width * r.periods, "the nanoribbon")
if r.edge == "zigzag":
# Zigzag chains run along x, a_CC * 3/2 apart; a period is a along x.
across = (r.width - 1) * A_LATTICE * math.sqrt(3) / 2 + A_LATTICE * math.sqrt(3) / 6
flat = _sheet_lattice(r.periods * A_LATTICE, across)
points = np.column_stack([flat[:, 1], flat[:, 0], flat[:, 2]]) # length along y
period = A_LATTICE
else:
# Dimer lines run along y, a/2 apart; a period is sqrt(3) a along y.
across = (r.width - 1) * A_LATTICE / 2
points = _sheet_lattice(across + 1e-3, r.periods * math.sqrt(3) * A_LATTICE - 1e-3)
period = math.sqrt(3) * A_LATTICE
name = "AGNR" if r.edge == "armchair" else "ZGNR"
description = (
f"{r.width}-{name} graphene nanoribbon, {r.periods} periods "
f"({r.periods * period / 10:.2f} nm), {r.edges}-terminated"
)
model = _flat(points, r.edges, r, description)
if not r.sidewall and not r.defects:
model.periodic = PeriodicCell(
translation_angstrom=round(period, 6), carbons_per_period=2 * r.width
)
return model
[docs]
def build(
parameters: NanostructureParameters | dict[str, Any],
) -> tuple[NanostructureParameters, AtomModel]:
"""The atom model of ``parameters``; raises :class:`NanostructureRefused` with the reason.
A dict is validated first (a ``pydantic.ValidationError`` names the field).
"""
if isinstance(parameters, dict):
from pydantic import TypeAdapter
parameters = TypeAdapter(NanostructureParameters).validate_python(parameters)
from vibesketch.limits import LimitExceeded
try:
if isinstance(parameters, Nanotube):
return parameters, _tube(parameters)
if isinstance(parameters, Graphene):
return parameters, _graphene(parameters)
if isinstance(parameters, Nanoribbon):
return parameters, _ribbon(parameters)
return parameters, _fullerene(parameters)
except LimitExceeded as error:
raise NanostructureRefused(str(error)) from None