# aligfellow_xyzgraph **Repository Path**: wuping4321/aligfellow_xyzgraph ## Basic Information - **Project Name**: aligfellow_xyzgraph - **Description**: No description available - **Primary Language**: Unknown - **License**: MIT - **Default Branch**: main - **Homepage**: None - **GVP Project**: No ## Statistics - **Stars**: 0 - **Forks**: 0 - **Created**: 2026-04-29 - **Last Updated**: 2026-04-29 ## Categories & Tags **Categories**: Uncategorized **Tags**: None ## README # xyzgraph: Molecular Graph Construction from Cartesian Coordinates **xyzgraph** is a Python toolkit for building molecular graphs (bond connectivity, bond orders, formal charges, and partial charges) directly from 3D atomic coordinates in XYZ format. It provides both **cheminformatics-based** and **quantum chemistry-based** (xTB) workflows. [![PyPI Downloads](https://static.pepy.tech/badge/xyzgraph)](https://pepy.tech/projects/xyzgraph) [![License](https://img.shields.io/github/license/aligfellow/xyzgraph)](https://github.com/aligfellow/xyzgraph/blob/main/LICENSE) [![Powered by: uv](https://img.shields.io/badge/-uv-purple)](https://docs.astral.sh/uv) [![Code style: ruff](https://img.shields.io/badge/code%20style-ruff-000000.svg)](https://github.com/astral-sh/ruff) [![Typing: ty](https://img.shields.io/badge/typing-ty-EFC621.svg)](https://github.com/astral-sh/ty) [![GitHub Workflow Status](https://img.shields.io/github/actions/workflow/status/aligfellow/xyzgraph/ci.yml?branch=main&logo=github-actions)](https://github.com/aligfellow/xyzgraph/actions) [![Codecov](https://img.shields.io/codecov/c/github/aligfellow/xyzgraph)](https://codecov.io/gh/aligfellow/xyzgraph) --- ## Table of Contents 1. [Key Features](#key-features) 2. [Installation](#installation) 3. [Quick Start](#quick-start) 4. [Methodology Overview](#methodology-overview) 5. [Workflow Comparison](#workflow-comparison) 6. [CLI Reference](#cli-reference) 7. [Python API](#python-api) 8. [Visualization](#visualization) 9. [Non-Covalent Interactions](#non-covalent-interactions) 10. [Limitations & Future Work](#limitations--future-work) 11. [Examples](#examples) 12. [References](#references) 13. [Contributing & Contact](#contributing--contact) --- ## Key Features - **Distance-based initial bonding** using *consistent* van der Waals radii across *all elements* from Charry and Tkatchenko [[1]](https://doi.org/10.1021/acs.jctc.4c00784) - **Four construction methods**: - `cheminf`: Pure cheminformatics with bond order optimization - `xtb`: Semi-empirical calculation via xTB Wiberg bond orders with Mulliken charges [[2]](https://pubs.acs.org/doi/10.1021/acs.jctc.8b01176) - `rdkit`: RDKit's DetermineBonds algorithm [[3]](https://github.com/jensengroup/xyz2mol), [[4]](https://github.com/rdkit) - `orca`: Reads Mayer bond orders and Mulliken charges from ORCA outputs. - **Cheminformatics modes**: - `--quick`: Geometric detection only — no bond order assignment, formal charges, or aromatic detection - `--optimizer`: Full optimization with valence and charge minimisation - **beam**: optimization across multiple paths (slightly slower, default) - **greedy**: iterative valence adjustment - **Aromatic detection**: Hückel 4n+2 rule for 5/6-membered rings, plus a fused-perimeter fallback for non-alternant systems (azulene, acenaphthylene). Optional `--kekule` to keep Kekulé bond orders - **Charge computation**: Gasteiger (cheminf) or Mulliken (xTB/ORCA) partial charges - **Stereochemistry assignment**: R/S, E/Z, axial (biaryl, allene, metallocene), planar (metallocene, paracyclophane), and helical (helicene) chirality from 3D geometry - **RDkit/xyz2mol comparison** validation against RDKit bond perception [[3]](https://github.com/jensengroup/xyz2mol), [[4]](https://github.com/rdkit) - **Non-covalent interaction (NCI) detection**: 17 interaction types including hydrogen bonds, pi-stacking, halogen/chalcogen/pnictogen bonds, cation-pi, and more - **ASCII 2D depiction** with layout alignment for method comparison (see also [[5]](https://github.com/whitead/moltext)) --- ## Installation ### From PyPI ```bash pip install xyzgraph ``` ### From Source ```bash git clone https://github.com/aligfellow/xyzgraph.git cd xyzgraph pip install . # or simply pip install git+https://github.com/aligfellow/xyzgraph.git ``` ### Dependencies - **Core**: `numpy`, `networkx`, `rdkit` - **Optional**: [xTB binary](https://github.com/grimme-lab/xtb) (for `--method xtb`) - **Optional**: [xyz2mol_tm](https://github.com/jensengroup/xyz2mol_tm) + `scipy` (for `--compare-rdkit-tm`) To install xTB (Linux/macOS) see [here](https://github.com/grimme-lab/xtb): ```bash conda install -c conda-forge xtb # or download from GitHub releases ``` To install xyz2mol_tm (required for `--compare-rdkit-tm`): ```bash pip install "xyzgraph[rdkit-tm]" xyz2mol_tm@git+https://github.com/jensengroup/xyz2mol_tm.git ``` This installs `scipy` (via the `rdkit-tm` extra) and `xyz2mol_tm` from source in one command. This extra step is necessary because `xyz2mol_tm` is not hosted on `pypi`. --- ## Quick Start ### CLI Examples **Minimal usage** (auto-displays ASCII depiction): ```bash xyzgraph molecule.xyz # constructs graph with cheminformatics style defaults xyzgraph molecule.out # constructs graph from ORCA output ``` **Specify charge and method**: ```bash xyzgraph molecule.xyz --method xtb --charge -1 --multiplicity 2 ``` **Detailed debug output**: ```bash xyzgraph molecule.xyz --debug ``` **Compare with RDKit**: ```bash xyzgraph molecule.xyz --compare-rdkit ``` **Compare with ORCA output**: ```bash # Compare XYZ (cheminf) vs ORCA bond orders xyzgraph molecule.xyz --orca-out molecule.out # Three-way comparison: cheminf vs ORCA vs RDKit xyzgraph molecule.xyz --orca-out molecule.out --compare-rdkit ``` **Detect non-covalent interactions**: ```bash xyzgraph molecule.xyz --nci --charge 1 ``` **Multi-frame trajectory files**: ```bash # Process specific frame from trajectory (0-indexed) xyzgraph trajectory.xyz --frame 2 # Process all frames for quick topological overview xyzgraph trajectory.xyz --all-frames ``` ### Python Example **Basic usage**: ```python from xyzgraph import build_graph, build_graph_rdkit, build_graph_orca # Cheminformatics (default method) G_cheminf = build_graph("molecule.xyz", charge=0) # RDKit's DetermineBonds G_rdkit = build_graph_rdkit("molecule.xyz", charge=0) # ORCA output (Mayer bond orders) G_orca = build_graph_orca("structure.out", bond_threshold=0.5) # Print ASCII structure from xyzgraph import graph_to_ascii print(graph_to_ascii(G_cheminf, scale=3.0, include_h=False)) ``` **Multi-frame trajectory files**: ```python from xyzgraph import read_xyz_file, build_graph # Read specific frame from trajectory atoms = read_xyz_file("trajectory.xyz", frame=2) G = build_graph(atoms, charge=0) # Process all frames from xyzgraph import count_frames_and_atoms num_frames, _ = count_frames_and_atoms("trajectory.xyz") for i in range(num_frames): atoms = read_xyz_file("trajectory.xyz", frame=i) G = build_graph(atoms, charge=0) # ... analyze G ``` **Non-covalent interaction detection**: ```python from xyzgraph import build_graph from xyzgraph.nci import detect_ncis G = build_graph("molecule.xyz", charge=1) ncis = detect_ncis(G) for nci in ncis: print(nci.type, nci.site_a, nci.site_b, nci.geometry) ``` **Comparing methods**: ```python from xyzgraph import compare_with_rdkit # Build graphs G_cheminf = build_graph("molecule.xyz", charge=-1) G_rdkit = build_graph_rdkit("molecule.xyz", charge=-1) # Compare (returns formatted report) report = compare_with_rdkit(G_cheminf, G_rdkit, verbose=True, ascii=True) print(report) ``` --- ## Methodology Overview ### Design Philosophy xyzgraph offers two distinct pathways for molecular graph construction: 1. **Cheminformatics Path** (`method='cheminf'`): - Pure graph-based approach using chemical heuristics - No external quantum chemistry calls - Cached scoring, valence, edge and graph properties - Fast and suitable for both organic *and* inorganic molecules 2. **Quantum Chemistry Path** (`method='xtb'`): - Uses GFN2-xTB (extended tight-binding) calculations [[2]](https://pubs.acs.org/doi/10.1021/acs.jctc.8b01176) - Reads in Wiberg bond orders and Mulliken charges from output - Potentially more accurate for unusual bonding situations - *though, xTB may be less robust in these situations* - Requires xTB binary installation ### Cheminformatics Workflow (method='cheminf') ``` ┌─────────────────────────────────────────────────────────────────┐ │ 1. Input Processing │ │ • Parse XYZ file internally │ │ • Load reference data (VDW radii, valences, electrons) │ └────────────────────┬────────────────────────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 2. Initial Bond Graph (Two-Step Construction) │ │ │ │ Step 1: Baseline Bonds (DEFAULT thresholds) │ │ • Uses DEFAULT threshold parameters (threshold=1.0) │ │ • Builds reliable "core" connectivity │ │ • Bonds sorted by confidence: 1.0 (short) to 0.0 (at thresh) │ │ • High confidence (>0.4): added directly │ │ • Low confidence (≤0.4): geometric validation applied │ │ • Result: stable molecular scaffold │ │ • Compute rings via shortest-cycle-per-edge BFS (SSSR-like, │ │ chemically natural smallest rings; avoids 11-ring or │ │ larger artefacts from a plain cycle basis) │ │ │ │ Step 2: Extended Bonds (if using CUSTOM thresholds) │ │ • Sorted highest-confidence-first (most reliable first) │ │ • Additional bonds require geometric validation: │ │ - Acute angle check: 15° (metals) / 30° (non-metals) │ │ - Collinearity check: trans vs spurious detection │ │ - Existing ring diagonal rejection and 3-ring validation │ | - Agostic bond filtering: H-M/F-M bonds rejected if │ │ stronger H-X or F-X bond exists (2x confidence ratio) │ │ - M-L priority check: diagonal M-ligand bonds in 3-rings │ │ rejected if stronger M-donor bond exists in ring (2x) │ │ • Allows sensible elongated bonds (e.g., TS structures) │ │ │ │ • Create graph with single bonds (order = 1.0) │ └────────────────────┬────────────────────────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 3. Kekulé Initialization for Conjugated Rings │ │ • Find 5/6-membered planar rings with C/N/O/S/B/P/Se │ │ • Initialize alternating bond orders (5-ring: 2-1-2-1-1, │ │ 6-ring: 2-1-2-1-2-1) │ │ • Priority passes: Cp-like metal-bound 5-rings, heteroatom │ │ 5-rings, fused benzene blocks, isolated rings │ │ • Handle fused rings (naphthalene, anthracene) by │ │ propagating shared-edge patterns across components │ │ • Gives optimizer a good starting point; remaining locally │ │ inconsistent pairs are resolved later via Kekulé shifts │ │ • Broader atom set than aromatic detection (P, Se included) │ └────────────────────┬────────────────────────────────────────────┘ │ ┌──────────┴─────────────┐ │ │ ┌─────────▼────────────┐ ┌───────▼──────────────────────────────┐ │ 4a. Quick Mode │ │ 4b. Full Optimization │ │ • Lock metal bonds │ │ • Lock metal bonds at 1.0 │ │ • 3 iterations │ │ • Iterative BIDIRECTIONAL search: │ │ • Promote bonds │ │ - Test both +1 AND -1 changes │ │ where both atoms │ │ - Allows Kekulé structure swaps │ │ need increased │ │ • Score = f(valence_error, │ │ valence │ │ formal_charges, │ │ • Distance check │ │ electronegativity, │ │ │ │ conjugation deficit) │ │ │ │ • Optimizer choice: │ │ │ │ - Beam: parallel hypotheses │ │ │ │ - Greedy: single best change │ │ │ │ • Top-k edge candidate selection │ │ │ │ • Kekulé-shift fallback when single-│ │ │ │ edge moves stall: flips a whole │ │ │ │ alternating chain between two │ │ │ │ valence-deficient atoms at once │ │ │ │ (escapes fused-ring traps) │ └─────────┬────────────┘ └──────────┬───────────────────────────┘ └───────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 5. Aromatic Detection (Hückel 4n+2) │ │ • Per-ring: 5/6-membered rings with C/N/O/S/B │ │ - Count π electrons (sp² C → 1e, N/O/S LP → 2e, B → 0) │ │ - Lone pairs are counted in every ring an atom belongs to │ │ - Apply Hückel 4n+2; set ring bonds to 1.5 │ │ • Fused-perimeter fallback for non-alternant systems │ │ (azulene 5-7, acenaphthylene 5-6-6) where no individual │ │ ring is Hückel but the perimeter of the union is: │ │ promotes every bond inside the component to 1.5 │ │ • --kekule / kekule=True keeps Kekulé bond orders but still │ │ stores aromatic ring metadata │ │ • Other heteroatoms (e.g. P, Se) use Kekulé structures │ └────────────────────┬────────────────────────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 6. Formal Charge Assignment │ │ • For each non-metal atom: │ │ - B = 2 × Σ(bond_orders, excluding metal coord) │ │ - target = min(8, 2 × V_electrons) │ │ (octet for C/N/O/F/..., sextet for B/Al/Ga — no │ │ spurious lone pair on trivalent boron) │ │ - L = max(0, target - B) │ │ - formal = V_electrons - L - B/2 │ │ • Balance total to match system charge │ │ • Metals assigned oxidation state as formal charge │ └────────────────────┬────────────────────────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 7. Optional: Gasteiger Partial Charges │ │ • compute_gasteiger_charges(G, target_charge) │ │ • Convert bond orders to RDKit bond types │ │ • Compute Gasteiger charges │ │ • Adjust for total charge conservation │ │ • Aggregate H charges onto heavy atoms │ │ • Stored in G.nodes[i]["charges"]["gasteiger"] │ └────────────────────┬────────────────────────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 8. Optional: Non-Covalent Interaction Detection (--nci) │ │ • Classify sites: donors, acceptors, ions, halogens, etc. │ │ • Detect pi-systems: aromatic rings + conjugated domains │ │ • Enumerate candidate pairs (graph-distance filtered) │ │ • Geometry checks: distances, angles, plane alignment │ │ • 17 interaction types (H-bond, pi-stack, sigma-hole, ...) │ │ • Stored in G.graph["ncis"] as list[NCIData] │ └────────────────────┬────────────────────────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 9. Output Graph │ │ Nodes: symbol, formal_charge, valence, metal_valence, │ │ oxidation_state (metals only) │ │ Edges: bond_order, bond_type, metal_coord │ └─────────────────────────────────────────────────────────────────┘ ``` ### xTB Workflow (method='xtb') ``` ┌─────────────────────────────────────────────────────────────────┐ │ 1. Input Processing | │ • Parse XYZ file internally │ │ • Write XYZ to temporary directory │ │ • Set up xTB calculation parameters │ └────────────────────┬────────────────────────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 2. Run xTB Calculation │ │ Command: xtb .xyz --chrg --uhf │ │ • GFN2-xTB Hamiltonian │ │ • Single-point calculation │ │ • Wiberg bond order analysis │ │ • Mulliken population analysis │ └────────────────────┬────────────────────────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 3. Parse xTB Output │ │ • Read wbo file (Wiberg bond orders) │ │ • Read charges file (Mulliken atomic charges) │ │ • Threshold: bond_order > 0.5 → create edge │ └────────────────────┬────────────────────────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 4. Build Graph from xTB Data │ │ • Create nodes with Mulliken charges │ │ • Create edges with Wiberg bond orders │ │ • No further optimization needed │ └────────────────────┬────────────────────────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 5. Cleanup (optional) │ │ • Remove temporary xTB files (unless --no-clean) │ └────────────────────┬────────────────────────────────────────────┘ │ ┌────────────────────▼────────────────────────────────────────────┐ │ 6. Output Graph │ │ Nodes: symbol, charges{'mulliken': ...}, agg_charge, │ │ valence, metal_valence │ │ Edges: bond_order (Wiberg), bond_type, metal_coord │ └─────────────────────────────────────────────────────────────────┘ ``` --- ## Workflow Comparison | Feature | cheminf (quick) | cheminf (full) | xtb | |---------|----------------|----------------|-----| | **Speed** | Fastest | Fast | Moderate | | **Accuracy** | Connectivity only | Very good across various systems | Only limited by xTB performance (QM-based) | | **External deps** | None | None | Requires xTB binary | | **Bond orders** | None (all bonds set to 1) | Optimized formal charge and valency | Wiberg (fractional) | | **Charges** | None | Gasteiger | Mulliken | | **Metal complexes** | Detection only | Reasonable | Reasonable (limited by xTB metal performance) | | **Conjugated systems** | Not detected | Excellent | Excellent | | **Best for** | Topology screening, pre-processing | Most cases | Awkward bonding, validation | ### When to Use Each Method **Use `--method cheminf` (default)**: - Most use cases - No xTB installation available - Batch processing structures **Use `--method cheminf --quick`**: - Extremely large molecules or batch pre-processing - When only connectivity (which atoms are bonded) is needed **Use `--method xtb`**: - Validation of cheminf results - Unusual electronic structures - Low confidence in bonding structure ### Optimizer Algorithms (cheminf full mode only) **Beam Search Optimizer** (`--optimizer beam` default, `--beam-width 5` default): - Explores multiple optimization paths in parallel - Maintains top-k hypotheses at each iteration (of top candidates) - Bidirectional: tests both +1 and -1 bond orders for each hypothesis - **Kekulé-shift fallback**: when no single-edge move improves any beam, searches for an alternating-BO chain (1,2,1,2,…,1) between two valence-deficient atoms and flips the whole chain in one move. Escapes local-optimum traps in fused-ring systems where every single-edge change violates valence on a saturated neighbour (e.g. BN-doped PAHs, fused-ring carbanions). Only fires on stall so it can't pre-empt legitimate single-edge solutions (like promoting a double to a triple in benzyne) - More robust against local minima - Slower, but better convergence - Best for robust bonding assignment across periodic table **Greedy Optimizer** (`--optimizer greedy`): - Tests all top candidate edges, picks single best change per iteration - Bidirectional: tests both +1 and -1 bond order changes - Fast and effective for most molecules - Can get stuck in local minima (*e.g.* alpha, beta unsaturated systems) --- ## CLI Reference ### Command Syntax ```text > xyzgraph -h usage: xyzgraph [-h] [--version] [--citation] [--method {cheminf,xtb}] [--no-clean] [-c CHARGE] [-m MULTIPLICITY] [-q] [-k] [--relaxed] [-t THRESHOLD] [-d] [-a] [--json] [-as ASCII_SCALE] [--nci] [--stereo] [-H] [--show-h-idx SHOW_H_IDX] [-b] [--frame FRAME] [--all-frames] [--compare-rdkit] [--compare-rdkit-tm] [--orca-out ORCA_OUT] [--orca-threshold ORCA_THRESHOLD] [-o {greedy,beam}] [-bw BEAM_WIDTH] [--max-iter MAX_ITER] [--edge-per-iter EDGE_PER_ITER] [--bond BOND] [--unbond UNBOND] [--threshold-h-h THRESHOLD_H_H] [--threshold-h-nonmetal THRESHOLD_H_NONMETAL] [--threshold-h-metal THRESHOLD_H_METAL] [--threshold-metal-ligand THRESHOLD_METAL_LIGAND] [--threshold-sblock-ligand THRESHOLD_SBLOCK_LIGAND] [--threshold-nonmetal THRESHOLD_NONMETAL] [--allow-metal-metal-bonds] [--threshold-metal-metal-self THRESHOLD_METAL_METAL_SELF] [--period-scaling-h-bonds PERIOD_SCALING_H_BONDS] [--period-scaling-nonmetal-bonds PERIOD_SCALING_NONMETAL_BONDS] [--period-scaling-sblock-bonds PERIOD_SCALING_SBLOCK_BONDS] [input_file] Build molecular graph from XYZ or ORCA output. positional arguments: input_file Input file (XYZ or ORCA .out) options: -h, --help show this help message and exit --version Print version and exit --citation Print citation and exit Common Options: --method {cheminf,xtb} Graph construction method (default: cheminf) --no-clean Keep temporary xTB files (only for --method xtb) -c, --charge CHARGE Total molecular charge (default: 0) -m, --multiplicity MULTIPLICITY Spin multiplicity (default: auto estimation) -q, --quick Geometric detection only: skip bond order assignment, formal charges, and aromatic detection -k, --kekule Keep Kekule bond orders (do not convert aromatic rings to 1.5) --relaxed Relaxed geometric validation (for transition states) -t, --threshold THRESHOLD Global scaling for bond thresholds (default: 1.0) Output Options: -d, --debug Enable debug output -a, --ascii Show 2D ASCII depiction --json Output graph as JSON (for generating test fixtures) -as, --ascii-scale ASCII_SCALE ASCII scaling factor (default: 2.5) --nci Detect and report non-covalent interactions --stereo Assign and display stereochemistry labels -H, --show-h Include hydrogens in visualizations --show-h-idx SHOW_H_IDX Show specific H atoms (comma-separated indices) Input Options: -b, --bohr XYZ file in Bohr units (default: Angstrom) --frame FRAME Frame index for trajectory files, 0-indexed (default: 0) --all-frames Process all frames in trajectory Comparison Options: --compare-rdkit Compare with RDKit graph --compare-rdkit-tm Compare with RDKit xyz2mol_tm graph --orca-out ORCA_OUT ORCA output file for comparison --orca-threshold ORCA_THRESHOLD Min Mayer bond order for ORCA (default: 0.25) Optimizer Options: -o, --optimizer {greedy,beam} Algorithm (default: beam) -bw, --beam-width BEAM_WIDTH Beam width (default: 5) --max-iter MAX_ITER Max iterations (default: 50) --edge-per-iter EDGE_PER_ITER Edges per iteration (default: 10) Bond Constraints: --bond BOND Force bonds (e.g., --bond 0,1 2,3) --unbond UNBOND Prevent bonds (e.g., --unbond 0,1) Advanced Thresholds: --threshold-h-h THRESHOLD_H_H H-H vdW threshold (default: 0.38) --threshold-h-nonmetal THRESHOLD_H_NONMETAL H-nonmetal vdW threshold (default: 0.42) --threshold-h-metal THRESHOLD_H_METAL H-metal vdW threshold (default: 0.45) --threshold-metal-ligand THRESHOLD_METAL_LIGAND D-block metal-ligand vdW threshold (default: 0.65) --threshold-sblock-ligand THRESHOLD_SBLOCK_LIGAND S-block metal-ligand vdW threshold (default: 0.55) --threshold-nonmetal THRESHOLD_NONMETAL Nonmetal-nonmetal vdW threshold (default: 0.55) --allow-metal-metal-bonds Allow metal-metal bonds (default: True) --threshold-metal-metal-self THRESHOLD_METAL_METAL_SELF Metal-metal vdW threshold (default: 0.7) --period-scaling-h-bonds PERIOD_SCALING_H_BONDS Period scaling for H bonds (default: 0.05) --period-scaling-nonmetal-bonds PERIOD_SCALING_NONMETAL_BONDS Period scaling for nonmetal bonds (default: 0.0) --period-scaling-sblock-bonds PERIOD_SCALING_SBLOCK_BONDS Period scaling for s-block M-L bonds (default: 0.05) ``` **Method comparison**: ```bash xyzgraph molecule.xyz --debug > cheminf.txt xyzgraph molecule.xyz --method xtb --debug > xtb.txt diff cheminf.txt xtb.txt ``` **Validate against RDKit**: ```bash xyzgraph molecule.xyz --compare-rdkit ``` --- ## Python API Direct graph construction: ```python from xyzgraph import build_graph, graph_debug_report # Cheminf full optimization G_full = build_graph( atoms='molecule.xyz', charge=0, max_iter=50, # maximum iterations (normally converged <20) edge_per_iter=6, # default 10 bond=[(0,1)], # ensure a bond between 0 and 1 debug=True ) # Keep Kekule bond orders (no 1.5 aromatic conversion) G_kekule = build_graph('molecule.xyz', kekule=True) # Aromatic rings are still detected and stored in G_kekule.graph["aromatic_rings"] # Bond orders remain as optimised Kekule values (1.0/2.0) ``` Stereochemistry assignment (R/S, E/Z, axial, planar, helical): ```python from xyzgraph import annotate_stereo stereo = annotate_stereo(G_full) # assigns all stereo types at once # Also stored as G.graph["stereo"] for downstream use (e.g. rendering) # Each key maps to a list of dicts with "label" and type-specific atom references: stereo["point"] # [{"atom": 5, "label": "R"}, ...] — R/S centres stereo["ez"] # [{"bond": [4,6], "label": "E"}, ...] — geometric isomerism stereo["axial"] # [{"atoms": [9,10], "label": "Rₐ"}, ...] — biaryl, allene, metallocene stereo["planar"] # [{"ring": [...], "label": "Rₚ"}, ...] — metallocene, paracyclophane stereo["helical"] # [{"atoms": [12,0], "label": "M"}, ...] — helicenes ``` Example structures for each stereo type are provided in `examples/stereo/`: | Type | Examples | Label | |------|----------|-------| | Point (R/S) | `mnh.xyz` | R, S | | E/Z | `E_2butene.xyz`, `Z_2butene.xyz` | E, Z | | Axial (biaryl) | `R_binol.xyz`, `S_binol.xyz` | Rₐ, Sₐ | | Axial (allene) | `Ra_allene.xyz`, `Sa_allene.xyz` | Rₐ, Sₐ | | Axial (metallocene) | `Ra_ferrocene_axial.xyz`, `Sa_ferrocene_axial.xyz` | Rₐ, Sₐ | | Axial (hindered biaryl) | `Ra_hindered_biaryl.xyz`, `Sa_hindered_biaryl.xyz` | Rₐ, Sₐ | | Planar (metallocene) | `Rp_ferrocene.xyz`, `Sp_ferrocene.xyz` | Rₚ, Sₚ | | Planar (paracyclophane) | `22paracyclophane.xyz`, `22paracyclophane_F.xyz` | None, Sₚ | | Helical | `M_helicene.xyz`, `P_helicene.xyz` | M, P | **Conventions**: - **R/S, E/Z, axial**: standard CIP rules - **Planar (general)**: IUPAC pilot-atom convention (CW from pilot = Rₚ) - **Planar (metallocene)**: Schlögl convention (view from opposite the metal, CW = Rₚ). Note: this gives the opposite label from IUPAC CIP when the metal is the pilot atom - **Helical**: IUPAC P/M helix convention **Known limitations**: - Axial chirality requires sp2 junction atoms with at least one in an aromatic ring; non-conjugated restricted rotations are not detected - Ortho steric gating (≥2 non-H ortho substituents) may be overly permissive for some aryl-vinyl or aryl-amide axes where free rotation is possible in practice - Planar chirality in paracyclophanes with chemically different bridges (e.g. -CH₂CH₂- vs -CH₂O-) may not be detected on the unsubstituted deck - Helical chirality requires fused aromatic rings; large aza-helicenes with ring detection issues may not be assigned --- ## Visualization ### ASCII Depiction xyzgraph includes a built-in ASCII renderer for 2D molecular structures. This is heavily inspired by work elsewhere, *e.g.* [[5]](https://github.com/whitead/moltext) by Andrew White. ```python from xyzgraph import graph_to_ascii # Basic rendering ascii_art = graph_to_ascii(G, scale=3.0, include_h=False) print(ascii_art) ``` **Output example** (acyl isothiouronium): ```text > xyzgraph examples/isothio.xyz -a /C / /// C\ \\ \ \ \\ C\ // // O=======C =========\ C---- \ /S\ // ---C \ / \\ // \ N---- /// \\ ----C\ C \ // ---C\ \C--- \ \ \ / \\ / \\\ \ C--- // \ / C \ // ----C \\ / / C--- // \ N\------C / ----C \ /// \\\ / \ / \ ---C C-------C/ \C---- // C---- // ---C \ \ \ C ``` **Features**: - Single bonds: `-`, `|`, `/`, `\` - Double bonds: `=`, `‖` (parallel lines) - Triple bonds: `#` - Aromatic: 1.5 bond orders shown as single - Special edges: `*` (TS), `.` (NCI) if `G.edges[i,j]['TS']=True` or `G.edges[i,j]['NCI']=True` ### Layout Alignment Compare methods by aligning their ASCII depictions: ```python from xyzgraph import build_graph, graph_to_ascii # Build with both methods G_cheminf = build_graph(atoms, method='cheminf') G_xtb = build_graph(atoms, method='xtb') # Generate aligned depictions ascii_ref, layout = graph_to_ascii(G_cheminf) ascii_xtb = graph_to_ascii(G_xtb, reference_layout=layout) print("Cheminf:\n", ascii_ref) print("\nxTB:\n", ascii_xtb) ``` ### Debug Report Tabular listing of all atoms and bonds: ```python from xyzgraph import graph_debug_report report = graph_debug_report(G, include_h=False) print(report) ``` **Full example**: ```text > xyzgraph benzene_NH4-cation-pi.xyz -c 1 -a -d ================================================================================ XYZGRAPH Molecular Graph Construction from Cartesian Coordinates A. S. Goodfellow, 2025 ================================================================================ Version: xyzgraph v1.5.0 Citation: A. S. Goodfellow, xyzgraph: Molecular Graph Construction from Cartesian Coordinates, v1.5.0, 2025, https://github.com/aligfellow/xyzgraph.git. Input: benzene_NH4-cation-pi.xyz Parameters: charge=1 ================================================================================ # Building cheminf graph from examples/benzene_NH4-cation-pi.xyz... ================================================================================ BUILDING GRAPH (CHEMINF, FULL MODE) Atoms: 17, Charge: 1, Multiplicity: 1 ================================================================================ Added 17 atoms Chemical formula: C6H10N Step 1: Found 16 baseline bonds (using default thresholds) ... ... ... Step 1: 16 baseline bonds added, 0 rejected Found 1 rings from initial bonding (excluding metal cycles) Total bonds in graph: 16 Initial bonds: 16 ================================================================================ KEKULE INITIALIZATION FOR AROMATIC RINGS ================================================================================ Ring 0 (6-membered): ['C0', 'C1', 'C2', 'C3', 'C4', 'C5'] π electrons estimate: 6 -------------------------------------------------------------------------------- Valid rings for Kekulé initialization: [0] ✓ Initialized isolated 6-ring 0 -------------------------------------------------------------------------------- SUMMARY: Initialized 1 ring(s) with Kekulé pattern -------------------------------------------------------------------------------- ================================================================================ BEAM SEARCH OPTIMIZATION (width=5) ================================================================================ Initial score: 22.50 Iteration 1: No improvements found in any beam, stopping Applying best solution to graph... -------------------------------------------------------------------------------- Explored 13 states across 1 iterations Found 0 improvements Score: 22.50 → 22.50 -------------------------------------------------------------------------------- ================================================================================ FORMAL CHARGE CALCULATION ================================================================================ Initial formal charges: Sum: +1 (target: +1) Charged atoms: N12: +1 No residual charge distribution needed (sum matches target) ================================================================================ AROMATIC RING DETECTION (Hückel 4n+2) ================================================================================ Ring 1 (6-membered): ['C0', 'C1', 'C2', 'C3', 'C4', 'C5'] π electrons: 6 (C0:1, C1:1, C2:1, C3:1, C4:1, C5:1) ✓ AROMATIC (4n+2 rule: n=1) -------------------------------------------------------------------------------- SUMMARY: 1 aromatic rings, 6 bonds set to 1.5 -------------------------------------------------------------------------------- ================================================================================ GRAPH CONSTRUCTION COMPLETE ================================================================================ Constructed graph with chemical formula: C6H10N ================================================================================ # CHEMINF GRAPH DETAILS ================================================================================ # Molecular Graph: 17 atoms, 16 bonds # total_charge=1 multiplicity=1 # (C-H hydrogens hidden; heteroatom-bound hydrogens shown; valences still include all H) # [idx] Sym val=.. metal=.. formal=.. | neighbors: idx(order / aromatic flag) # (val = organic valence excluding metal bonds; metal = metal coordination bonds) [ 0] C val=4.00 metal=0.00 formal=0 | 1(1.50*) 5(1.50*) [ 1] C val=4.00 metal=0.00 formal=0 | 0(1.50*) 2(1.50*) [ 2] C val=4.00 metal=0.00 formal=0 | 1(1.50*) 3(1.50*) [ 3] C val=4.00 metal=0.00 formal=0 | 2(1.50*) 4(1.50*) [ 4] C val=4.00 metal=0.00 formal=0 | 3(1.50*) 5(1.50*) [ 5] C val=4.00 metal=0.00 formal=0 | 0(1.50*) 4(1.50*) [ 12] N val=4.00 metal=0.00 formal=+1 | 13(1.00) 14(1.00) 15(1.00) 16(1.00) [ 13] H val=1.00 metal=0.00 formal=0 | 12(1.00) [ 14] H val=1.00 metal=0.00 formal=0 | 12(1.00) [ 15] H val=1.00 metal=0.00 formal=0 | 12(1.00) [ 16] H val=1.00 metal=0.00 formal=0 | 12(1.00) # Bonds (i-j: order) (filtered) [ 0- 1]: 1.50 [ 0- 5]: 1.50 [ 1- 2]: 1.50 [ 2- 3]: 1.50 [ 3- 4]: 1.50 [ 4- 5]: 1.50 [12-13]: 1.00 [12-14]: 1.00 [12-15]: 1.00 [12-16]: 1.00 ================================================================================ # ASCII Depiction (cheminf) ================================================================================ -C------------------------C- --- --- ---- ---- --- --- C\ -C \\ // \\\ /// \\\ /// \\ // \C------------------------C/ H | | | | H------------------------N-------------------------H | | | | H ``` --- ## Non-Covalent Interactions xyzgraph includes a geometry-based NCI detection module that identifies 17 types of non-covalent interactions from the molecular graph. NCI detection runs on top of the constructed graph and requires no additional dependencies. ### CLI Usage ```bash xyzgraph molecule.xyz --nci --charge 1 ``` ### Supported Interaction Types | Type | Description | |---|---| | `hbond` | Classical hydrogen bond (D-H...A) | | `hbond_bifurcated` | Two donors sharing the same acceptor | | `halogen_bond` | Sigma-hole bond via halogen (X...A) | | `chalcogen_bond` | Sigma-hole bond via S, Se, Te | | `pnictogen_bond` | Sigma-hole bond via P, As, Sb, Bi | | `pi_pi_parallel` | Parallel-displaced pi-stacking (ring-ring) | | `pi_pi_t_shaped` | T-shaped (edge-to-face) pi-stacking | | `pi_pi_ring_domain` | Pi-stacking between ring and non-ring pi domain | | `pi_pi_domain_domain` | Pi-stacking between two non-ring pi domains | | `cation_pi` | Cation above aromatic ring | | `anion_pi` | Anion above aromatic ring | | `halogen_pi` | Halogen sigma-hole to pi-system | | `ch_pi` | C-H...pi interaction | | `hb_pi` | H-bond donor to pi-system | | `cation_lp` | Cation to lone pair donor | | `ionic` | Electrostatic cation-anion | | `salt_bridge` | H-mediated ionic (cation-H...anion) | ### Python API ```python from xyzgraph import build_graph from xyzgraph.nci import detect_ncis, NCIThresholds G = build_graph("molecule.xyz", charge=1) ncis = detect_ncis(G) for nci in ncis: print(nci.type, nci.site_a, nci.site_b, nci.geometry) # Results are also stored on the graph ncis = G.graph["ncis"] # Custom thresholds thr = NCIThresholds(hb_da_max=3.0, pii_parallel_rmax=4.0) ncis = detect_ncis(G, thresholds=thr) ``` For trajectory analysis where topology is shared across frames, use `NCIAnalyzer` to avoid repeating site classification on every frame: ```python from xyzgraph.nci import NCIAnalyzer analyzer = NCIAnalyzer(G) # topology work done once for positions in trajectory_frames: ncis = analyzer.detect(positions) # geometry checks only ``` ### Example: Acyl Isothiouronium (Chalcogen Bond) ```bash xyzgraph examples/isothio.xyz -c 1 --nci ``` ```text ================================================================================ # Non-Covalent Interactions ================================================================================ 1 interaction(s) detected: chalcogen_bond S20 ... O0 ================================================================================ # ASCII Depiction (with NCI dotted lines) ================================================================================ C\ \\ =C \\ =======\\ \C= === \\ =O === \\ =======. C= === . === . | . -C\ C---- | S---- ----- \\ / ----C | / ---C- \\ // \ N-- // | \C / \ / --- / | | C \ / --C | | \ \ / \ | | \ C---- / \\ | | \ / ---C \ ---C\ /C \ // \ N---- \\ /// C---- / \ / \\ // ----C \ // \C/ \ / C--------C / // C--- / ----C \ \ \ \ C ``` The intramolecular S...O chalcogen bond is detected via the sigma-hole on S20 directed towards the carbonyl oxygen O0, shown as a dotted line in the ASCII depiction. --- ## Current Limitations 1. **Metal Complexes** - Bond orders locked at 1.0 (no d-orbital chemistry) - Metal-metal bonds *partially* supported but not well tested (single bond allowed) - Can deal with **both** ionic *and* neutral ligands 2. **Radicals & Open-Shell Systems** - Unlikely to appropriately solve a valence structure - Not explicity dealt with currently - *May* behave, *may* be unreliable 3. **Zwitterions** - Formal charge and valence analysis does identify `-[N+](=O)(-[O-])` bonding and formal charge pattern - This is performed **without pattern matching** - *May* not always be fully robust, and does not account for delocalisation 4. **Charged Aromatics** - Hückel electron counting is simplistic - Should still solve with valence/charge optimisation 5. **Inorganic Cages** - Homogeneous clusters (≥8 atoms, same element) bypass standard ring validation - Unlikely to be fully accurately described, *e.g.* C/B cage structures --- ### Built-in Comparison xyzgraph can directly compare its output to rdkit/xyz2mol [[3]](https://github.com/jensengroup/xyz2mol), [[4]](https://github.com/rdkit) or to rdkit/xyz2mol_tm [[6]](https://github.com/jensengroup/xyz2mol_tm), [[7]](https://jcheminf.biomedcentral.com/articles/10.1186/s13321-025-01008-1): ```bash xyzgraph molecule.xyz --compare-rdkit --debug # or xyzgraph molecule.xyz --compare-rdkit-tm --debug # integrates graph building from xyz2mol_tm ``` **Output includes**: - Layout-aligned ASCII depictions - Edge differences (bonds only in one method) - Bond order differences (Δ ≥ 0.25) **Example**: ```text # Bond differences: only_in_native=1 only_in_rdkit=0 bond_order_diffs=2 # only_in_native: 4-7 # bond_order_diffs (Δ≥0.25): # 1-2 native=1.50 rdkit=1.00 Δ=+0.50 # 2-3 native=2.00 rdkit=1.50 Δ=+0.50 ``` --- ## Examples This section demonstrates xyzgraph's capabilities on real molecular systems, showcasing Kekulé initialization, aromatic detection, metal coordination analysis, and formal charge assignment. ### Example 1: Metal Complex (Ferrocene-Manganese Hydride) This example demonstrates xyzgraph's handling of organometallic complexes with multiple ligand types. **System:** [(η⁵-Cp)₂Fe][Mn(H)(CO)₂(PNN)] - Ferrocene cation with manganese hydride complex **File:** `examples/mnh.xyz` (77 atoms) **Command:** ```bash xyzgraph examples/mnh.xyz --ascii --debug ``` **Key Features:** - Detection of Cp⁻ (cyclopentadienyl) rings coordinated to Fe - Metal coordination summary (Fe²⁺, Mn¹⁺) with ligand classification - Hydride ligand (H⁻) recognition - Carbonyl (CO) ligands with triple-bonded oxygen - Aromatic Cp rings with charge contribution to π system **Output (truncated):** ```text ================================================================================ KEKULE INITIALIZATION FOR AROMATIC RINGS ================================================================================ Ring 0 (6-membered): ['N6', 'C52', 'C58', 'C57', 'C55', 'C53'] π electrons estimate: 6 Ring 1 (5-membered): ['C7', 'C8', 'C9', 'C11', 'C13'] ✓ Detected Cp-like ring (all 5 C bonded to Fe0) π electrons estimate: 6 Ring 2 (5-membered): ['C15', 'C17', 'C19', 'C21', 'C23'] ✓ Detected Cp-like ring (all 5 C bonded to Fe0) π electrons estimate: 6 Ring 3 (6-membered): ['C25', 'C26', 'C28', 'C30', 'C32', 'C34'] π electrons estimate: 6 Ring 4 (6-membered): ['C36', 'C37', 'C39', 'C41', 'C43', 'C45'] π electrons estimate: 6 -------------------------------------------------------------------------------- Valid rings for Kekulé initialization: [0, 1, 2, 3, 4] ✓ Cp-like 5-ring 1 initialized (rotation 0) ✓ Cp-like 5-ring 2 initialized (rotation 0) ✓ Initialized isolated 6-ring 0 ✓ Initialized isolated 6-ring 3 ✓ Initialized isolated 6-ring 4 -------------------------------------------------------------------------------- SUMMARY: Initialized 5 ring(s) with Kekulé pattern -------------------------------------------------------------------------------- ================================================================================ BEAM SEARCH OPTIMIZATION (width=5) ================================================================================ Locked 16 metal bonds Initial score: 392.70 Iteration 1: Generated 2 candidates, keeping top 2 ✓ New best: O3-C64 Δtotal = 81.00 score = 311.70 Iteration 2: Generated 4 candidates, keeping top 4 ✓ New best: O4-C65 Δtotal = 81.00 score = 230.70 Iteration 3: Generated 6 candidates, keeping top 5 ✓ New best: O3-C64 Δtotal = 20.00 score = 210.70 Iteration 4: Generated 5 candidates, keeping top 5 ✓ New best: O4-C65 Δtotal = 20.00 score = 190.70 Iteration 5: No improvements (single or Kekulé shift) found, stopping Applying best solution to graph... -------------------------------------------------------------------------------- Explored 198 states across 5 iterations Found 4 improvements Score: 392.70 → 190.70 -------------------------------------------------------------------------------- ================================================================================ FORMAL CHARGE CALCULATION ================================================================================ Initial formal charges: Sum: -3 (target: +0) Metal coordination summary: [ 0] Fe oxidation_state=+2 coordination=10 • 5-ring (-1) [donor: C13] • 5-ring (-1) [donor: C19] [ 1] Mn oxidation_state=+1 coordination=6 • H (-1) [donor: H67] • CO ( 0) [donor: C64] • CO ( 0) [donor: C65] • N ( 0) [donor: N6] • P ( 0) [donor: P2] • N ( 0) [donor: N5] Metal complex detected: Residual: +3 (represents metal oxidation states) Fe0: formal_charge=+2 Mn1: formal_charge=+1 ================================================================================ AROMATIC RING DETECTION (Hückel 4n+2) ================================================================================ Ring 1 (6-membered): ['N6', 'C52', 'C58', 'C57', 'C55', 'C53'] π electrons: 6 (N6:1, C52:1, C58:1, C57:1, C55:1, C53:1) ✓ AROMATIC (4n+2 rule: n=1) Ring 2 (5-membered): ['C7', 'C8', 'C9', 'C11', 'C13'] π electrons: 6 (C7:2(fc=-1), C8:1, C9:1, C11:1, C13:1) ✓ AROMATIC (4n+2 rule: n=1) Ring 3 (5-membered): ['C15', 'C17', 'C19', 'C21', 'C23'] π electrons: 6 (C15:2(fc=-1), C17:1, C19:1, C21:1, C23:1) ✓ AROMATIC (4n+2 rule: n=1) Ring 4 (6-membered): ['C25', 'C26', 'C28', 'C30', 'C32', 'C34'] π electrons: 6 (C25:1, C26:1, C28:1, C30:1, C32:1, C34:1) ✓ AROMATIC (4n+2 rule: n=1) Ring 5 (6-membered): ['C36', 'C37', 'C39', 'C41', 'C43', 'C45'] π electrons: 6 (C36:1, C37:1, C39:1, C41:1, C43:1, C45:1) ✓ AROMATIC (4n+2 rule: n=1) -------------------------------------------------------------------------------- SUMMARY: 5 aromatic rings, 28 bonds set to 1.5 -------------------------------------------------------------------------------- ================================================================================ GRAPH CONSTRUCTION COMPLETE ================================================================================ ### Selected atoms from molecular graph: [ 0] Fe val=10.00 metal=0.00 formal=+2 | 7(1.00) 8(1.00) 9(1.00) 11(1.00) 13(1.00) 15(1.00) 17(1.00) 19(1.00) 21(1.00) 23(1.00) [ 1] Mn val=6.00 metal=0.00 formal=+1 | 2(1.00) 5(1.00) 6(1.00) 64(1.00) 65(1.00) 67(1.00) [ 3] O val=3.00 metal=0.00 formal=+1 | 64(3.00) [ 4] O val=3.00 metal=0.00 formal=+1 | 65(3.00) [ 8] C val=4.00 metal=1.00 formal=-1 | 0(1.00) 7(1.50*) 9(1.50*) 47(1.00) [ 23] C val=4.00 metal=1.00 formal=-1 | 0(1.00) 15(1.50*) 21(1.50*) [ 64] C val=3.00 metal=1.00 formal=-1 | 1(1.00) 3(3.00) [ 65] C val=3.00 metal=1.00 formal=-1 | 1(1.00) 4(3.00) [ 67] H val=0.00 metal=1.00 formal=-1 | 1(1.00) ``` **ASCII Depiction:** >[!TIP] > Avert your eyes... Not good for complex molecular visualisation... ```text C---------C / \ / \ C-- / \ // ---- / \ // --C / C // | C / C | \ / | | \ / | | \ / | O | \ / | # C C--------C |# // \ #C-- // \ // ---- // \\ /C H --C C---------C C C---- \ // \ / / \ / / \ -----C --P \ / C#####/ \ / C----\--- /\ ---- \\ \ / // /####O \ / |\\\ \ --//----C-- \\ \ /// / \ / /| \ \ / ---\| \\ // ----N C---------N C----- \\\ /--- | Mn---- \ / \ | ----Fe--- | | \ / \ | ---- /|\\ ----C | \ / \ C-- /| \\---| | \ / \ \\ // |---- \\| | C---------C \ \\/---| --C\ | / C C-\\| ---- \\\ --N- // C-- \\ ---- | --- / \C-- | --- / | | -C | | | | | H | C ``` ![mnh](examples/mnh.svg) **Analysis:** - **Ferrocene fragment:** Fe(II) coordinated to two Cp⁻ ligands (η⁵ coordination) - **Cp rings:** Detected as aromatic with 6 π electrons (includes -1 charge contribution from each ring) - **Manganese center:** Mn(I) with octahedral-like coordination - Hydride (H⁻) ligand correctly identified (formal charge -1) - Two CO ligands with C≡O triple bonds (formal charges: C: -1, O: +1), net neutral ligand - Phosphine (P) and amine (N) dative bond donors - **Charge balance:** System is neutral (Fe(II) + Mn(I) - 2×Cp⁻ - H⁻ = 0) --- ### Example 2: Organic Cation (Acyl Isothiouronium) This example shows aromatic detection, formal charge assignment, and handling of heteroaromatic systems. **System:** Acyl isothiouronium cation (quaternary nitrogen) **File:** `examples/isothio.xyz` (52 atoms, +1 charge) **Command:** ```bash xyzgraph examples/isothio.xyz --charge 1 --ascii --debug ``` **Key Features:** - Benzene ring aromatic detection - 5-membered heterocycle evaluation (thiazole-like ring) - Formal charge on quaternary nitrogen (N⁺) - Beam search optimization of carbonyl bond order **Output:** ```text > xyzgraph examples/isothio.xyz -a -d -c 1 ================================================================================ BUILDING GRAPH (CHEMINF, FULL MODE) Atoms: 52, Charge: 1, Multiplicity: 1 ================================================================================ Added 52 atoms Chemical formula: C23H25N2OS Step 1: Found 55 baseline bonds (using default thresholds) Step 1: 55 baseline bonds added, 0 rejected Found 4 rings from initial bonding (excluding metal cycles) Total bonds in graph: 55 Initial bonds: 55 ================================================================================ KEKULE INITIALIZATION FOR AROMATIC RINGS ================================================================================ Ring 0 (6-membered): ['N5', 'C6', 'C13', 'C17', 'N18', 'C19'] ✗ Not planar Ring 1 (6-membered): ['C7', 'C8', 'C9', 'C10', 'C11', 'C12'] π electrons estimate: 6 Ring 2 (5-membered): ['N18', 'C19', 'S20', 'C21', 'C26'] π electrons estimate: 7 ✗ Hückel rule violated (π=7) Ring 3 (6-membered): ['C21', 'C22', 'C23', 'C24', 'C25', 'C26'] π electrons estimate: 6 -------------------------------------------------------------------------------- Valid rings for Kekulé initialization: [1, 3] ✓ Initialized isolated 6-ring 1 ✓ Initialized isolated 6-ring 3 -------------------------------------------------------------------------------- SUMMARY: Initialized 2 ring(s) with Kekulé pattern -------------------------------------------------------------------------------- ================================================================================ BEAM SEARCH OPTIMIZATION (width=5) ================================================================================ Initial score: 489.00 Iteration 1: Generated 5 candidates, keeping top 5 ✓ New best: N18-C19 Δtotal = 116.50 score = 372.50 Iteration 2: Generated 17 candidates, keeping top 5 ✓ New best: C1-C2 Δtotal = 72.00 score = 300.50 Iteration 3: Generated 3 candidates, keeping top 3 ✓ New best: O0-C1 Δtotal = 71.00 score = 229.50 Iteration 4: No improvements (single or Kekulé shift) found, stopping Applying best solution to graph... -------------------------------------------------------------------------------- Explored 152 states across 4 iterations Found 3 improvements Score: 489.00 → 229.50 -------------------------------------------------------------------------------- ================================================================================ FORMAL CHARGE CALCULATION ================================================================================ Initial formal charges: Sum: +1 (target: +1) Charged atoms: N18: +1 No residual charge distribution needed (sum matches target) ================================================================================ AROMATIC RING DETECTION (Hückel 4n+2) ================================================================================ Ring 1 (6-membered): ['N5', 'C6', 'C13', 'C17', 'N18', 'C19'] ✗ Not planar, skipping aromaticity check Ring 2 (6-membered): ['C7', 'C8', 'C9', 'C10', 'C11', 'C12'] π electrons: 6 (C7:1, C8:1, C9:1, C10:1, C11:1, C12:1) ✓ AROMATIC (4n+2 rule: n=1) Ring 3 (5-membered): ['N18', 'C19', 'S20', 'C21', 'C26'] π electrons: 6 (N18:1(fc=+1), C19:1, S20:2(LP), C21:1, C26:1) ✓ AROMATIC (4n+2 rule: n=1) Ring 4 (6-membered): ['C21', 'C22', 'C23', 'C24', 'C25', 'C26'] π electrons: 6 (C21:1, C22:1, C23:1, C24:1, C25:1, C26:1) ✓ AROMATIC (4n+2 rule: n=1) -------------------------------------------------------------------------------- SUMMARY: 3 aromatic rings, 16 bonds set to 1.5 -------------------------------------------------------------------------------- ================================================================================ GRAPH CONSTRUCTION COMPLETE ================================================================================ Constructed graph with chemical formula: C23H25N2OS ================================================================================ # CHEMINF GRAPH DETAILS ================================================================================ # Molecular Graph: 52 atoms, 55 bonds # total_charge=1 multiplicity=1 # (C-H hydrogens hidden; heteroatom-bound hydrogens shown; valences still include all H) # [idx] Sym val=.. metal=.. formal=.. | neighbors: idx(order / aromatic flag) # (val = organic valence excluding metal bonds; metal = metal coordination bonds) [ 0] O val=2.00 metal=0.00 formal=0 | 1(2.00) [ 1] C val=4.00 metal=0.00 formal=0 | 0(2.00) 2(1.00) 5(1.00) [ 2] C val=4.00 metal=0.00 formal=0 | 1(1.00) 3(2.00) [ 3] C val=4.00 metal=0.00 formal=0 | 2(2.00) 4(1.00) [ 4] C val=4.00 metal=0.00 formal=0 | 3(1.00) [ 5] N val=3.00 metal=0.00 formal=0 | 1(1.00) 6(1.00) 19(1.00) [ 6] C val=4.00 metal=0.00 formal=0 | 5(1.00) 7(1.00) 13(1.00) [ 7] C val=4.00 metal=0.00 formal=0 | 6(1.00) 8(1.50*) 12(1.50*) [ 8] C val=4.00 metal=0.00 formal=0 | 7(1.50*) 9(1.50*) [ 9] C val=4.00 metal=0.00 formal=0 | 8(1.50*) 10(1.50*) [ 10] C val=4.00 metal=0.00 formal=0 | 9(1.50*) 11(1.50*) [ 11] C val=4.00 metal=0.00 formal=0 | 10(1.50*) 12(1.50*) [ 12] C val=4.00 metal=0.00 formal=0 | 7(1.50*) 11(1.50*) [ 13] C val=4.00 metal=0.00 formal=0 | 6(1.00) 14(1.00) 17(1.00) [ 14] C val=4.00 metal=0.00 formal=0 | 13(1.00) 15(1.00) 16(1.00) [ 15] C val=4.00 metal=0.00 formal=0 | 14(1.00) [ 16] C val=4.00 metal=0.00 formal=0 | 14(1.00) [ 17] C val=4.00 metal=0.00 formal=0 | 13(1.00) 18(1.00) [ 18] N val=4.00 metal=0.00 formal=+1 | 17(1.00) 19(1.50*) 26(1.50*) [ 19] C val=4.00 metal=0.00 formal=0 | 5(1.00) 18(1.50*) 20(1.50*) [ 20] S val=3.00 metal=0.00 formal=0 | 19(1.50*) 21(1.50*) [ 21] C val=4.50 metal=0.00 formal=0 | 20(1.50*) 22(1.50*) 26(1.50*) [ 22] C val=4.00 metal=0.00 formal=0 | 21(1.50*) 23(1.50*) [ 23] C val=4.00 metal=0.00 formal=0 | 22(1.50*) 24(1.50*) [ 24] C val=4.00 metal=0.00 formal=0 | 23(1.50*) 25(1.50*) [ 25] C val=4.00 metal=0.00 formal=0 | 24(1.50*) 26(1.50*) [ 26] C val=4.50 metal=0.00 formal=0 | 18(1.50*) 21(1.50*) 25(1.50*) # Bonds (i-j: order) (filtered) [ 0- 1]: 2.00 [ 1- 2]: 1.00 [ 1- 5]: 1.00 [ 2- 3]: 2.00 [ 3- 4]: 1.00 [ 5- 6]: 1.00 [ 5-19]: 1.00 [ 6- 7]: 1.00 [ 6-13]: 1.00 [ 7- 8]: 1.50 [ 7-12]: 1.50 [ 8- 9]: 1.50 [ 9-10]: 1.50 [10-11]: 1.50 [11-12]: 1.50 [13-14]: 1.00 [13-17]: 1.00 [14-15]: 1.00 [14-16]: 1.00 [17-18]: 1.00 [18-19]: 1.50 [18-26]: 1.50 [19-20]: 1.50 [20-21]: 1.50 [21-22]: 1.50 [21-26]: 1.50 [22-23]: 1.50 [23-24]: 1.50 [24-25]: 1.50 [25-26]: 1.50 ================================================================================ # ASCII Depiction (cheminf) ================================================================================ /C / /// C\ \\ \ \ \\ C\ // // O=======C =========\ C---- \ /S\ // ---C \ / \\ // \ N---- /// \\ ----C\ C \ // ---C \C--- \ \ \ / \ / \\\ \ C--- // \ / C \ // ----C \ / / C--- // \ N-------C / ----C \ /// \\\ / \ / \ ---C C-------C/ \C---- // C---- // ---C \ \ \ C ``` ![isothiouronium](examples/isothio.svg) **Analysis:** - **Benzene rings:** Two rings correctly identified as aromatic (bond order 1.5) - **5-membered heterocycle:** N-C-S-C-C ring retains Kekulé structure with N=C double bond - **Quaternary nitrogen:** N16 assigned +1 formal charge (4 bonds, no lone pairs) - **a,b-unsaturated:** O=C and C=C double bonds correctly optimized --- ## Bond Detection Thresholds xyzgraph uses distance-based bond detection with thresholds derived from van der Waals (vdW) radii by Charry and Tkatchenko [[1]](https://doi.org/10.1021/acs.jctc.4c00784). By default, these thresholds are calibrated for different atom pair types: | Atom Pair Type | Default Threshold | Parameter Name | |---------------|-------------------|----------------| | H-H | 0.38 × (r₁ + r₂) | `threshold_h_h` | | H-nonmetal | 0.42 × (r₁ + r₂) | `threshold_h_nonmetal` | | H-metal | 0.45 × (r₁ + r₂) | `threshold_h_metal` | | S-block metal-ligand | 0.55 × (r₁ + r₂) | `threshold_sblock_ligand` | | D-block metal-ligand | 0.65 × (r₁ + r₂) | `threshold_metal_ligand` | | Nonmetal-nonmetal | 0.55 × (r₁ + r₂) | `threshold_nonmetal_nonmetal` | | Metal-Metal (same type) | 0.7 × (2r) | `threshold_metal_metal_self` | Where r₁ and r₂ are the VDW radii of the two atoms. Period-dependent scaling is applied for H-bonds (`period_scaling_h_bonds=0.05`) and s-block metal-ligand bonds (`period_scaling_sblock_bonds=0.05`) to account for heavier elements bonding at longer distances. ### Detecting Elongated Bonds (Transition States) The two-step construction allows detection of elongated bonds in transition state structures by adjusting the global threshold: ```bash # Detect elongated bonds in TS structures xyzgraph ts_structure.xyz --threshold 1.2 --debug # For more dense connectivity, one can use relaxed mode (more permissive geometric validation) xyzgraph structure.xyz --threshold 1.2 --relaxed --debug ``` **Recommended threshold ranges**: - **1.0** (default): Ground-state structures - **1.1-1.2**: Slightly elongated bonds - **1.2-1.3**: Transition states with stretched geometries - **≥1.35**: Unstable - spurious bonding likely The two-step construction with geometric validation helps reject spurious diagonals even at higher thresholds. The `--relaxed` flag can be used for more permissive angle and diagonal thresholds (**but note:** *this is likely to produce spurious structures*). **Example workflow**: See [vib_analysis](https://github.com/aligfellow/graphRC) for a complete workflow analyzing transition state vibrational modes using xyzgraph connectivity. ### Advanced Threshold Modification (Not Recommended) **Global Scaling**: - The `--threshold` (or `threshold` in Python) parameter provides a simple way to globally scale **all** thresholds. - This is safer than modifying individual thresholds. - e.g. `--threshold 1.1` - threshold_h_nonmetal × (r₁ + r₂) × **1.1** **Individual Scaling**: These parameters are exposed for users who need to: - Handle unusual bonding situations not covered by defaults - Specifically wish to obtain dense connectivity - Fine-tune bond detection for specific molecular systems - Debug or validate bond detection behavior Can be performed using the cli *e.g.* `--threshold_h_nonmetal 0.5` or directly in python within `build_graph(threshold_h_nonmetal=0.5)` > [!WARNING] > Modifying these thresholds is **not recommended** unless you have a specific reason and understand the implications > Changing values can produce *chemically invalid structures* --- ## References 1. **van der Waals Radii**: Jorge Charry and Alexandre Tkatchenko, *J. Chem. Theory Comput.*, 2024, **20**, 7469–7478. [DOI](https://doi.org/10.1021/acs.jctc.4c00784). 2. **xTB (Extended Tight Binding)**: Christoph Bannwarth, Sebastian Ehlert, and Stefan Grimme, *J. Chem. Theory Comput.* 2019, **15**, 1652–1671. [DOI](https://pubs.acs.org/doi/10.1021/acs.jctc.8b01176). [Repo](https://github.com/grimme-lab/xtb). 3. **xyz2mol**: Jan Jensen *et al.*, [xyz2mol](https://github.com/jensengroup/xyz2mol). Now integrated into RDKit as `Chem.rdDetermineBonds.DetermineBonds()`. See also Y. Kim, W. Y. Kim, *Bull. Korean Chem. Soc.*, 2015, **36**, 1769–1777. 4. **RDKit**: RDKit: Open-source cheminformatics. [https://www.rdkit.org](https://www.rdkit.org). [Repo](https://github.com/rdkit). 5. **moltext**: A. White, *moltext*. [Repo](https://github.com/whitead/moltext) 6. **xyz2mol_tm**: Jan Jensen *et al.*, [xyz2mol_tm](https://github.com/jensengroup/xyz2mol_tm). See also ref **7.**. 7. **SMILES all around: structure to SMILES conversion for transition metal complexes**: Maria H. Rasmussen, Magnus Strandgaard, Julius Seumer, Laura K. Hemmingsen, Angelo Frei, David Balcells and Jan H. Jensen, *Journal of Cheminformatics*, 2025, **17**. [DOI](https://jcheminf.biomedcentral.com/articles/10.1186/s13321-025-01008-1). --- ## Contributors - [James O'Brien (@JamesOBrien2)](https://github.com/JamesOBrien2) — stereochemistry detection ## Contributing & Contact Contributions welcome! Please open an issue or pull request and get in touch with any questions [here](https://github.com/aligfellow/xyzgraph/issues). ## Developer Setup To develop with `xyzgraph`, you can clone the repo and use [`just`](https://github.com/casey/just) and [`uv`](https://docs.astral.sh) to setup the dev environment: ```bash just setup ``` The cli can be used with: ```bash uv run xyzgraph filename.xyz ``` Run the checks using: ```bash just check ```