vibe-sketcherDocumentation
Source code for vibesketcher.chem.frameworks
"""Metal-organic frameworks assembled from a node, a ditopic linker and a net.
Reticular construction (Yaghi et al., Nature 423, 705, 2003; O'Keeffe et al.,
Acc. Chem. Res. 41, 1782, 2008, for the RCSR net symbols): a secondary building
unit (SBU) sits on each vertex of a net, a linker on each edge, and each
carboxylate's O-C-O plane is shared by the SBU and the linker it binds. The
lattice constant follows from the construction, edge = 2 x_C + L, with x_C the
SBU centre to carboxylate-carbon distance and L the linker's carboxylate-carbon
separation. Nothing is relaxed; coordinates are an unrelaxed starting geometry.
Nodes
``zn4o``: the octahedral Zn4O(CO2)6 cluster of MOF-5 and the IRMOF series
(Li, Eddaoudi, O'Keeffe and Yaghi, Nature 402, 276, 1999; Eddaoudi et al.,
Science 295, 469, 2002) on the pcu net. Neighbouring clusters are related
by inversion so both ends of a linker see the same carboxylate plane; the
conventional cell therefore holds eight clusters (the Fm-3m cell of MOF-5),
the primitive face-centred cell two.
``paddlewheel``: the M2(CO2)4 paddlewheel (Cu: Chui et al., Science 283,
1148, 1999) as a square-grid (sql) layer, optionally joined through its
axial sites by a ditopic N-donor pillar into a pcu framework
(Zn2(bdc)2(dabco), DMOF-1: Dybtsev, Chun and Kim, Angew. Chem. Int. Ed. 43,
5033, 2004).
Linkers are preset para-oligophenylene dicarboxylates (BDC, BPDC, TPDC; rings
kept coplanar) or any dicarboxylic acid given as SMILES, embedded with RDKit's
ETKDG v3 (Wang, Witek, Landrum and Riniker, J. Chem. Inf. Model. 60, 2044, 2020)
and MMFF94 (Halgren, J. Comput. Chem. 17, 490, 1996) under a fixed seed. A
SMILES linker must be linear ditopic: both carboxylate bonds within 15 degrees
of the axis through the carboxylate carbons.
"""
from __future__ import annotations
import math
from dataclasses import dataclass, field
from typing import Any, Literal, Self
import numpy as np
from numpy.typing import NDArray
from pydantic import BaseModel, ConfigDict, Field, model_validator
from vibesketcher.chem.periodic_geometry import (
BuildRefused,
_check_collisions,
_provenance,
_structure,
)
from vibesketcher.document.model import PeriodicStructure
__all__ = ["LINKERS", "NODES", "PILLARS", "FrameworkRecipe", "build_framework", "catalog"]
_Points = NDArray[np.float64]
#: Carboxylate C-O and O-C-O (Allen et al., J. Chem. Soc. Perkin Trans. 2, S1, 1987).
_R_CO = 1.27
_OCO = math.radians(125.0)
#: Aromatic C-C, inter-ring C-C, aryl-carboxylate C-C and aromatic C-H, angstrom.
_R_AR, _R_LINK, _R_CARB, _R_CH = 1.39, 1.49, 1.50, 1.08
#: Zn-(mu4-O) and Zn-O(carboxylate) in MOF-5 (Li et al. 1999).
_ZN4O = {"M-O4": 1.94, "M-O": 1.94}
#: Paddlewheel M-M, M-O(carboxylate), M-N(axial), angstrom (Chui 1999; Dybtsev 2004).
_PADDLEWHEEL = {
"Cu": {"M-M": 2.63, "M-O": 1.96, "M-N": 2.17},
"Zn": {"M-M": 2.94, "M-O": 2.03, "M-N": 2.05},
}
_SEED = 0x5EED
NODES: dict[str, Any] = {
"zn4o": {
"description": "Zn4O(CO2)6 octahedral cluster (MOF-5/IRMOF) on the pcu net",
"metals": ["Zn"],
"nets": {"none": "pcu"},
},
"paddlewheel": {
"description": "M2(CO2)4 paddlewheel on the sql layer; pcu with a pillar",
"metals": sorted(_PADDLEWHEEL),
"nets": {"none": "sql", "pillar": "pcu"},
},
}
LINKERS: dict[str, Any] = {
"bdc": {"description": "1,4-benzenedicarboxylate (terephthalate)", "rings": 1},
"bpdc": {"description": "biphenyl-4,4'-dicarboxylate", "rings": 2},
"tpdc": {"description": "p-terphenyl-4,4''-dicarboxylate", "rings": 3},
}
PILLARS: dict[str, Any] = {
"dabco": {"description": "1,4-diazabicyclo[2.2.2]octane", "smiles": "C1CN2CCN1CC2"},
"bipyridine": {"description": "4,4'-bipyridine", "smiles": "c1cc(ccn1)-c1ccncc1"},
}
[docs]
class FrameworkRecipe(BaseModel):
"""A MOF from a node, a ditopic dicarboxylate linker and the net they imply."""
model_config = ConfigDict(extra="forbid")
kind: Literal["framework"] = "framework"
identifier: str = Field(pattern=r"^[A-Za-z][A-Za-z0-9_-]{0,63}$")
node: Literal["zn4o", "paddlewheel"]
metal: str | None = None
linker: Literal["bdc", "bpdc", "tpdc"] | None = None
linker_smiles: str | None = Field(default=None, min_length=1, max_length=500)
pillar: Literal["dabco", "bipyridine"] | None = None
cell: Literal["conventional", "primitive"] = "conventional"
vacuum: float = Field(default=15.0, ge=5, le=200, allow_inf_nan=False)
@model_validator(mode="after")
def _choices(self) -> Self:
if (self.linker is None) == (self.linker_smiles is None):
raise ValueError("give either a linker preset or linker_smiles")
if self.pillar is not None and self.node != "paddlewheel":
raise ValueError("pillars join paddlewheel layers only")
if self.cell == "primitive" and self.node != "zn4o":
raise ValueError("the primitive cell choice applies to the zn4o (pcu) framework")
return self
[docs]
def catalog() -> dict[str, Any]:
"""Framework nodes, preset linkers and pillars, for clients that present choices."""
return {
"nodes": NODES,
"linkers": {name: entry["description"] for name, entry in LINKERS.items()},
"pillars": {name: entry["description"] for name, entry in PILLARS.items()},
}
@dataclass
class _Fragment:
"""Atoms in a local frame: x along the edge, y along the shared plane direction."""
symbols: list[str]
positions: _Points
bonds: list[tuple[int, int, float]]
start: int
end: int
length: float
notes: list[str] = field(default_factory=list)
@dataclass
class _Port:
direction: _Points
plane: _Points
carbon: int
@dataclass
class _Node:
symbols: list[str]
positions: _Points
bonds: list[tuple[int, int, float]]
ports: dict[tuple[int, int], _Port]
reach: float
#: Axial metal indices (top, bottom) of a paddlewheel.
axial: tuple[int, int] | None = None
# --- linkers -------------------------------------------------------------------
def _phenylenes(rings: int) -> _Fragment:
"""Coplanar para-oligophenylene between two carboxylate carbons, axis along x."""
symbols: list[str] = []
points: list[tuple[float, float, float]] = []
bonds: list[tuple[int, int, float]] = []
previous: list[int] = []
first: list[int] = []
for ring in range(rings):
cx = ring * (2 * _R_AR + _R_LINK)
index = []
for t in range(6):
angle = math.radians(60 * t)
index.append(len(symbols))
symbols.append("C")
points.append((cx + _R_AR * math.cos(angle), _R_AR * math.sin(angle), 0.0))
for t in range(6):
bonds.append((index[t], index[(t + 1) % 6], 1.5))
for t in (1, 2, 4, 5):
angle = math.radians(60 * t)
x, y, _ = points[index[t]]
bonds.append((index[t], len(symbols), 1.0))
symbols.append("H")
points.append((x + _R_CH * math.cos(angle), y + _R_CH * math.sin(angle), 0.0))
if previous:
bonds.append((previous[0], index[3], 1.0))
else:
first = index
previous = index
start, end = first[3], previous[0]
positions = np.array(points)
positions[:, 0] -= (positions[start, 0] + positions[end, 0]) / 2
span = float(positions[end, 0] - positions[start, 0])
return _Fragment(symbols, positions, bonds, start, end, span + 2 * _R_CARB)
def _embedded(smiles: str, what: str) -> Any:
from rdkit import Chem
from rdkit.Chem import AllChem as _all
from vibesketcher.chem._rdkit import untyped
AllChem = untyped(_all) # the stubs lack the embedding functions
mol = Chem.MolFromSmiles(smiles)
if mol is None:
raise BuildRefused(f"cannot read the {what} SMILES {smiles!r}")
if mol.GetNumHeavyAtoms() > 300:
raise BuildRefused(f"the {what} has more than 300 heavy atoms")
mol = Chem.AddHs(mol)
params = AllChem.ETKDGv3()
params.randomSeed = _SEED
params.numThreads = 1
if AllChem.EmbedMolecule(mol, params) != 0:
raise BuildRefused(f"RDKit could not embed the {what} in 3D")
if AllChem.MMFFHasAllMoleculeParams(mol):
AllChem.MMFFOptimizeMolecule(mol, maxIters=2000)
else:
AllChem.UFFOptimizeMolecule(mol, maxIters=2000)
return mol
def _frame(origin: _Points, axis: _Points, toward: _Points) -> tuple[_Points, _Points, _Points]:
"""Unit x along ``axis``, y the part of ``toward`` perpendicular to it, z = x × y."""
ex = axis / np.linalg.norm(axis)
perpendicular = toward - (toward @ ex) * ex
if np.linalg.norm(perpendicular) < 1e-6:
raise BuildRefused("cannot orient the fragment: reference lies on its axis")
ey = perpendicular / np.linalg.norm(perpendicular)
return ex, ey, np.cross(ex, ey)
def _smiles_linker(smiles: str) -> _Fragment:
from rdkit import Chem
mol = _embedded(smiles, "linker")
acid = Chem.MolFromSmarts("[CX3](=[OX1])[OX2H1,OX1-]")
groups = mol.GetSubstructMatches(acid)
carbons = sorted({g[0] for g in groups})
if len(carbons) != 2:
raise BuildRefused(
f"a linker needs exactly two carboxylic acid groups; {len(carbons)} found"
)
oxygens: dict[int, list[int]] = {c: [] for c in carbons}
for c, o1, o2 in groups:
for o in (o1, o2):
if o not in oxygens[c]:
oxygens[c].append(o)
removed: set[int] = set()
attach: dict[int, int] = {}
for c in carbons:
others = [
n.GetIdx() for n in mol.GetAtomWithIdx(c).GetNeighbors() if n.GetIdx() not in oxygens[c]
]
if len(others) != 1:
raise BuildRefused("each carboxyl carbon must bond to exactly one backbone atom")
attach[c] = others[0]
for o in oxygens[c]:
removed.add(o)
removed.update(
n.GetIdx() for n in mol.GetAtomWithIdx(o).GetNeighbors() if n.GetSymbol() == "H"
)
removed.add(c)
for atom in mol.GetAtoms():
if atom.GetIdx() not in removed and atom.GetFormalCharge() != 0:
raise BuildRefused("a charged linker backbone is not supported")
xyz = mol.GetConformer().GetPositions()
c1, c2 = carbons
axis = xyz[c2] - xyz[c1]
length = float(np.linalg.norm(axis))
unit = axis / length
for c, outward in ((c1, -unit), (c2, unit)):
bond = xyz[c] - xyz[attach[c]]
angle = math.degrees(
math.acos(float(np.clip(bond @ outward / np.linalg.norm(bond), -1, 1)))
)
if angle > 15:
raise BuildRefused(
f"the linker is not linear ditopic: a carboxylate bond is {angle:.0f} degrees "
"off the axis through the carboxylate carbons"
)
def plane(c: int) -> _Points:
o1, o2 = oxygens[c]
return np.asarray(xyz[o1] - xyz[o2], dtype=float)
ex, ey, ez = _frame(xyz[c1], axis, plane(c1))
second = plane(c2) - (plane(c2) @ ex) * ex
twist = math.degrees(
math.acos(min(1.0, abs(float(second @ ey)) / float(np.linalg.norm(second))))
)
keep = [a.GetIdx() for a in mol.GetAtoms() if a.GetIdx() not in removed]
order = {old: new for new, old in enumerate(keep)}
middle = (xyz[c1] + xyz[c2]) / 2
local = np.array([[(xyz[i] - middle) @ e for e in (ex, ey, ez)] for i in keep])
bonds = []
for b in mol.GetBonds():
i, j = b.GetBeginAtomIdx(), b.GetEndAtomIdx()
if i in order and j in order:
bonds.append((order[i], order[j], float(b.GetBondTypeAsDouble())))
notes = [f"SMILES linker embedded with RDKit {_rdkit_version()} (ETKDGv3, MMFF94/UFF)"]
if twist > 5:
notes.append(
f"the embedded linker's carboxylate planes differ by {twist:.0f} degrees; "
"the nodes impose one shared plane, so relax before analysis"
)
return _Fragment(
[mol.GetAtomWithIdx(i).GetSymbol() for i in keep],
local,
bonds,
order[attach[c1]],
order[attach[c2]],
length,
notes,
)
def _pillar(name: str) -> _Fragment:
mol = _embedded(PILLARS[name]["smiles"], "pillar")
nitrogens = [a.GetIdx() for a in mol.GetAtoms() if a.GetSymbol() == "N"]
assert len(nitrogens) == 2
xyz = mol.GetConformer().GetPositions()
n1, n2 = nitrogens
heavy = next(a.GetIdx() for a in mol.GetAtomWithIdx(n1).GetNeighbors() if a.GetSymbol() != "H")
ex, ey, ez = _frame(xyz[n1], xyz[n2] - xyz[n1], xyz[heavy] - xyz[n1])
middle = (xyz[n1] + xyz[n2]) / 2
local = np.array([[(p - middle) @ e for e in (ex, ey, ez)] for p in xyz])
bonds = [
(b.GetBeginAtomIdx(), b.GetEndAtomIdx(), float(b.GetBondTypeAsDouble()))
for b in mol.GetBonds()
]
symbols = [a.GetSymbol() for a in mol.GetAtoms()]
length = float(np.linalg.norm(xyz[n2] - xyz[n1]))
return _Fragment(symbols, local, bonds, n1, n2, length)
def _rdkit_version() -> str:
from rdkit import rdBase
return str(rdBase.rdkitVersion)
# --- nodes -------------------------------------------------------------------------
def _carboxylates(
metals: list[tuple[int, _Points]],
distance: float,
directions: list[tuple[int, int]],
symbols: list[str],
points: list[_Points],
bonds: list[tuple[int, int, float]],
) -> tuple[dict[tuple[int, int], _Port], float]:
"""Bridging carboxylates on the given (axis, sign) ports; returns ports and x_C."""
h = _R_CO * math.sin(_OCO / 2)
reach = 0.0
ports: dict[tuple[int, int], _Port] = {}
for axis, sign in directions:
e = np.zeros(3)
e[axis] = sign
pair = [(index, p) for index, p in metals if p @ e > 1e-9]
if len(pair) != 2:
pair = [(index, p) for index, p in metals if abs(p @ e) <= 1e-9]
(i1, z1), (i2, z2) = sorted(pair, key=lambda item: item[0])
along = float(z1 @ e)
perpendicular = z1 - along * e
w = float(np.linalg.norm(perpendicular))
u = perpendicular / w
x_o = along + math.sqrt(distance**2 - (h - w) ** 2)
x_c = x_o + _R_CO * math.cos(_OCO / 2)
reach = x_c
carbon = len(symbols)
symbols += ["C", "O", "O"]
points += [x_c * e, x_o * e + h * u, x_o * e - h * u]
bonds += [
(carbon, carbon + 1, 1.5),
(carbon, carbon + 2, 1.5),
(i1, carbon + 1, 1.0),
(i2, carbon + 2, 1.0),
]
ports[(axis, sign)] = _Port(e, u, carbon)
return ports, reach
def _zn4o(parity: int) -> _Node:
d = _ZN4O["M-O4"] / math.sqrt(3)
corners = [(1, 1, 1), (1, -1, -1), (-1, 1, -1), (-1, -1, 1)]
symbols = ["O"] + ["Zn"] * 4
points: list[_Points] = [np.zeros(3)] + [parity * d * np.array(c, float) for c in corners]
bonds: list[tuple[int, int, float]] = [(0, k, 1.0) for k in range(1, 5)]
metals = [(k, points[k]) for k in range(1, 5)]
directions = [(axis, sign) for axis in range(3) for sign in (1, -1)]
ports, reach = _carboxylates(metals, _ZN4O["M-O"], directions, symbols, points, bonds)
return _Node(symbols, np.array(points), bonds, ports, reach)
def _paddlewheel(metal: str) -> _Node:
sizes = _PADDLEWHEEL[metal]
half = sizes["M-M"] / 2
symbols = [metal, metal]
points: list[_Points] = [np.array([0, 0, half]), np.array([0, 0, -half])]
bonds: list[tuple[int, int, float]] = []
metals = [(0, points[0]), (1, points[1])]
directions = [(0, 1), (0, -1), (1, 1), (1, -1)]
ports, reach = _carboxylates(metals, sizes["M-O"], directions, symbols, points, bonds)
return _Node(symbols, np.array(points), bonds, ports, reach, axial=(0, 1))
# --- assembly ------------------------------------------------------------------------
class _Assembly:
def __init__(self) -> None:
self.symbols: list[str] = []
self.points: list[_Points] = []
self.bonds: list[tuple[int, int, float]] = []
def node(self, node: _Node, centre: _Points) -> int:
offset = len(self.symbols)
self.symbols += node.symbols
self.points += [centre + p for p in node.positions]
self.bonds += [(i + offset, j + offset, order) for i, j, order in node.bonds]
return offset
def edge(self, fragment: _Fragment, middle: _Points, axis: _Points, plane: _Points) -> int:
ex, ey, ez = _frame(middle, axis, plane)
offset = len(self.symbols)
self.symbols += fragment.symbols
self.points += [middle + p[0] * ex + p[1] * ey + p[2] * ez for p in fragment.positions]
self.bonds += [(i + offset, j + offset, order) for i, j, order in fragment.bonds]
return offset
def _find(nodes: list[_Points], target: _Points, inverse: _Points) -> int:
"""The node at ``target`` modulo lattice translations."""
f = target @ inverse
for index, centre in enumerate(nodes):
d = f - centre @ inverse
if np.allclose(d - np.round(d), 0, atol=1e-6):
return index
raise AssertionError("net edge does not end on a node")
def _linker(p: FrameworkRecipe) -> _Fragment:
if p.linker is not None:
return _phenylenes(int(LINKERS[p.linker]["rings"]))
assert p.linker_smiles is not None
return _smiles_linker(p.linker_smiles)
def _pcu_zn4o(p: FrameworkRecipe, identifier: str) -> PeriodicStructure:
metal = p.metal or "Zn"
if metal != "Zn":
raise BuildRefused("the zn4o node is built with Zn only")
linker = _linker(p)
nodes = {1: _zn4o(1), -1: _zn4o(-1)}
edge = 2 * nodes[1].reach + linker.length
if p.cell == "conventional":
lattice = 2 * edge * np.eye(3)
grid = [(i, j, k) for i in (0, 1) for j in (0, 1) for k in (0, 1)]
else:
lattice = edge * np.array([[0.0, 1, 1], [1, 0, 1], [1, 1, 0]])
grid = [(0, 0, 0), (1, 0, 0)]
centres = [edge * np.array(g, float) for g in grid]
parities = [1 if sum(g) % 2 == 0 else -1 for g in grid]
build = _Assembly()
offsets = [build.node(nodes[s], c) for s, c in zip(parities, centres, strict=True)]
inverse = np.linalg.inv(lattice)
for index, centre in enumerate(centres):
here = nodes[parities[index]]
for axis in range(3):
port = here.ports[(axis, 1)]
other = _find(centres, centre + edge * port.direction, inverse)
there = nodes[parities[other]].ports[(axis, -1)]
if abs(float(port.plane @ there.plane)) < 0.999:
raise AssertionError("carboxylate planes of a pcu edge disagree")
start = build.edge(
linker, centre + edge / 2 * port.direction, port.direction, port.plane
)
build.bonds.append((offsets[index] + port.carbon, start + linker.start, 1.0))
build.bonds.append((offsets[other] + there.carbon, start + linker.end, 1.0))
return _finish(
p,
identifier,
build,
lattice,
(True, True, True),
net="pcu",
nodes=len(centres),
edge=edge,
linker=linker,
multiplicity=1,
cell_setting=p.cell,
)
def _paddlewheel_net(p: FrameworkRecipe, identifier: str) -> PeriodicStructure:
metal = p.metal or "Cu"
if metal not in _PADDLEWHEEL:
raise BuildRefused(
f"the paddlewheel node is built with {' or '.join(sorted(_PADDLEWHEEL))}"
)
linker = _linker(p)
node = _paddlewheel(metal)
edge = 2 * node.reach + linker.length
build = _Assembly()
offset = build.node(node, np.zeros(3))
for axis in (0, 1):
port = node.ports[(axis, 1)]
there = node.ports[(axis, -1)]
start = build.edge(linker, edge / 2 * port.direction, port.direction, port.plane)
build.bonds.append((offset + port.carbon, start + linker.start, 1.0))
build.bonds.append((offset + there.carbon, start + linker.end, 1.0))
notes: list[str] = []
if p.pillar is None:
z = np.array(build.points)[:, 2]
height = float(z.max() - z.min()) + p.vacuum
shift = np.array([0.0, 0.0, p.vacuum / 2 - float(z.min())])
build.points = [q + shift for q in build.points]
lattice = np.diag([edge, edge, height])
pbc = (True, True, False)
net = "sql"
notes.append("open axial metal sites; real layers carry axial ligands or solvent")
else:
pillar = _pillar(p.pillar)
height = _PADDLEWHEEL[metal]["M-M"] + 2 * _PADDLEWHEEL[metal]["M-N"] + pillar.length
z = np.array([0.0, 0.0, 1.0])
between = np.array([1.0, 1.0, 0.0]) / math.sqrt(2)
start = build.edge(pillar, height / 2 * z, z, between)
assert node.axial is not None
top, bottom = node.axial
build.bonds.append((offset + top, start + pillar.start, 1.0))
build.bonds.append((offset + bottom, start + pillar.end, 1.0))
lattice = np.diag([edge, edge, height])
pbc = (True, True, True)
net = "pcu"
return _finish(
p,
identifier,
build,
lattice,
pbc,
net=net,
nodes=1,
edge=edge,
linker=linker,
multiplicity=1 if metal == "Zn" else None,
cell_setting="primitive",
notes=notes,
)
def _finish(
p: FrameworkRecipe,
identifier: str,
build: _Assembly,
lattice: _Points,
pbc: tuple[bool, bool, bool],
*,
net: str,
nodes: int,
edge: float,
linker: _Fragment,
multiplicity: int | None,
cell_setting: Any,
notes: list[str] | None = None,
) -> PeriodicStructure:
cartesian = np.array(build.points)
_check_collisions(cartesian, lattice, pbc)
spin = (
"closed-shell Zn(II)"
if multiplicity == 1
else "unspecified: coupled Cu(II) spins are not predicted"
)
return _structure(
identifier,
lattice,
pbc,
build.symbols,
cartesian,
build.bonds,
cell_setting=cell_setting,
charge=0,
multiplicity=multiplicity,
provenance=_provenance(
"framework",
p,
net=net,
nodes_per_cell=nodes,
edge_angstrom=round(edge, 6),
linker_carboxylate_separation_angstrom=round(linker.length, 6),
spin=spin,
warnings=[*(notes or []), *linker.notes],
),
)
[docs]
def build_framework(p: FrameworkRecipe) -> PeriodicStructure:
"""A MOF from its node, linker and (implied) net; refuses what it cannot build."""
if p.node == "zn4o":
return _pcu_zn4o(p, p.identifier)
return _paddlewheel_net(p, p.identifier)