Tutorial 9: Analyse Module - 3D Molecular Structure Analysis¶

Prerequisites: Tutorial 5 — Generation Overview · You’ll learn: optimising, validating, comparing, and featurising generated structures · Next: Tutorial 10 — Generation Sweeps

At a Glance¶


Objective	Validate and characterise generated 3D molecular structures.
You need	Generated XYZ files or an ASE database, plus optional analysis dependencies for the selected command.
Main command	`MolCraftDiff analyze <subcommand> ...`
Success looks like	The requested structures, metrics, tables, or feature arrays are written without silently skipping the input set.

This tutorial covers the analyse command group, which provides tools for post-generation analysis and validation of 3D molecular structures. The CLI spelling remains MolCraftDiff analyze.

Overview¶

The analysis command group includes five subcommands:

Command	Description
`optimize`	XTB geometry optimisation
`metrics`	Validity/connectivity metrics, plus the paired `conformer` set
`xyz2mol`	XYZ to SMILES + fingerprints
`xtb-electronic`	XTB electronic properties
`featurize`	Fixed-size molecular feature vectors (SOAP / UMA / SSL3D)

Access the CLI with:

MolCraftDiff analyze --help

Part 1: XTB Geometry Optimisation¶

Optimise generated structures using xTB (GFN1, GFN2, GFN-FF) or MMFF94.

Usage¶

MolCraftDiff analyze optimize gen_xyz/ --level gfn2 --charge 0

Options¶

Option	Default	Description
`-o, --output-path`	`input_dir/optimized_xyz`	Output directory
`-l, --level`	`gfn1`	Optimisation level: gfn1, gfn2, gfn-ff, mmff94
`-c, --charge`	`0`	Molecular charge
`-t, --timeout`	`240`	Timeout per molecule (seconds)
`-s, --scale-factor`	`1.3`	Covalent radii scale factor

Output¶

Optimised XYZ files are saved to output_dir/, each with an _opt suffix (e.g. mol_0000.xyz → mol_0000_opt.xyz). Part 3’s --metrics conformer expects exactly these *_opt.xyz files.

Part 2: Validity & Connectivity Metrics¶

Compute structural validation metrics for generated molecules.

Usage¶

MolCraftDiff analyze metrics gen_xyz/ --metrics all

The input may be a directory of .xyz files or an ASE .db file. For a database, every row is scored and the results carry an ase_db_row_id column so rows can be joined back to the source database.

Metric Types¶

Type	Description
`core`	Validity, connectivity, atom stability, and set-level uniqueness/novelty/diversity
`posebuster`	Bond lengths, angles, clashes
`geom_revised`	Aromatic-aware stability metrics
`druglike`	QED, SA, LogP, fsp3, MW, HBD, HBA, Lipinski, PAINS, ring statistics (+ optional `--rdkit-rmsd`)
`similarity3d`	Shape / ESP / pharmacophore similarity to `--reference-mol`
`sbdd`	AutoDock Vina binding affinity against a protein pocket (needs `--receptor` and the `[sbdd]` extra)
`all`	All of the above except `sbdd`, which needs a receptor and is opt-in

What `core` reports¶

Two validity columns are written per molecule, because they answer different questions:

Column	Meaning
`valid`	The perceived molecule sanitises in RDKit and its SMILES round-trips. This is the definition used throughout the literature, so it is the number to quote in comparisons.
`valid_geom`	Stricter, geometry-aware check: bond counts, bond lengths, and the coordination-shape measure of every atom. Use it to catch physically distorted structures that still sanitise.
`valid_connected`	`valid` and a single connected fragment.

Note: before this release valid held what is now valid_geom. If you have a sweep config or script that depends on the old, stricter number, point it at valid_geom.

Alongside the per-molecule CSV, core writes a <output>_summary.json with the aggregate rates (validity, connectivity, atom stability, uniqueness, novelty, diversity, the molecular-size statistics, and how many files failed to process). Sweeps pick this file up automatically.

Options¶

Option	Default	Description
`-o, --output`	None	Output CSV file
`-m, --metrics`	`all`	Metric type to compute
`--recheck-topo`	False	Recheck topology using RDKit
`--mol-converter`	`xyz2mol`	XYZ to mol converter
`-s, --split`	`1`	Deterministic splits for mean±std summary logging
`-p, --portion`	`1.0`	Fraction of XYZ files to process
`--filter` / `--filtered-output`	None	Keep only structures whose column is truthy (e.g. `--filter valid_connected`), from any metric set. Writes an XYZ directory, or an ASE `.db` when the input was a database.
`--train-smiles`	None	Reference SMILES (`.txt` one per line, or `.csv`) to score novelty against; without it novelty is reported as `n/a`
`--receptor`	None	Protein receptor for `sbdd`: a `.pdbqt`, or a `.pdb` that meeko prepares
`--ref-ligand`	None	Reference ligand `.sdf`; its own affinity becomes the bar for `high_affinity`
`--dock-mode`	`dock`	`score` (in place) / `min` (+ local optimisation) / `dock` (full redock)
`--rmsd-threshold`	`0.5`	`conformer`: RMSD cutoff (Å) for the coverage metric
`-c, --charge` / `-l, --level` / `--xtb-timeout`	`0` / `gfn2` / `120`	`conformer`: xTB settings for the strain column (ignored when `xtb` is not installed)
`--exhaustiveness`	`8`	Vina search effort for `--dock-mode dock`
`-r, --reference-mol` / `--mol-idx`	None / `0`	Reference `.pkl`/`.sdf` (and index) for `similarity3d`; `--mol-idx -1` picks a random reference per molecule
`--rdkit-rmsd` / `--rmsd-n-conf`	off / `20`	`druglike`: RMSD of the pose against UFF-optimised RDKit conformers (slow)
`--skip-atoms`	—	Atom indices to skip in validation
`-t, --timeout`	`10`	Timeout per xyz2mol conversion (s)

What `druglike` and `similarity3d` report¶

These two were one set called shepherd until they were split, because they answer different questions and need different inputs.

druglike needs nothing but the structures:

Column	Meaning
`QED`, `SA_score`, `LogP`, `fsp3`, `MW`, `HBD`, `HBA`	standard RDKit descriptors
`lipinski`	how many of the five Lipinski rules the molecule obeys, 0–5
`pains_pass`	free of PAINS-A alerts
`ring_filter_pass`	no non-aromatic ring larger than 6
`n_rings`, `n_aromatic_rings`, `n_aliphatic_rings`	ring counts
`ring_size_3` … `ring_size_9`	per molecule a boolean; the summary reports the fraction of molecules containing each size, which is the form the SBDD papers print
`rdkit_rmsd_min` / `_median` / `_max`	only with `--rdkit-rmsd`; embeds `--rmsd-n-conf` conformers per molecule, so it is far slower than the rest

similarity3d needs --reference-mol and reports shape_sim, esp_sim and pharm_sim — Gaussian-overlap scores after alignment.

Reproducibility: shape_sim and esp_sim are not deterministic. The molecular surface is sampled randomly and the overlay uses 45 random restarts, so re-running the same input moves them by up to ~0.05. Compare them within a run, and treat small differences between runs as noise. pharm_sim is deterministic.

Migration: --metrics shepherd no longer exists. Use --metrics druglike for the descriptors and --metrics similarity3d for the reference comparison. --metrics all runs druglike always, and similarity3d only when --reference-mol is given.

What `sbdd` reports¶

Binding affinity for pocket-conditioned models, the numbers the SBDD literature quotes. Requires the [sbdd] extra (pip install 'molcraftdiffusion[sbdd]').

MolCraftDiff analyze metrics gen_xyz/ --metrics sbdd \
    --receptor pocket.pdbqt --ref-ligand reference.sdf

Column / metric	Meaning
`vina_score`	Affinity of the generated pose exactly as produced, no relaxation
`vina_min`	After local minimisation in the pocket
`vina_dock`	After a full redock — the number most papers report
`high_affinity`	Fraction of molecules beating `--ref-ligand`’s own affinity
`success_rate`	QED > 0.25 ∧ SA > 0.59 ∧ `vina_dock` < −8.18

The search box is centred on each molecule’s own bounding box plus a 5 Å buffer, matching the reference implementations — generated ligands already sit in the pocket frame.

Two things worth knowing:

Poses are perceived with OpenBabel, not xyz2mol. xyz2mol returns a topology with no conformer, which would silently discard the coordinates the model generated. Structures that cannot be perceived, that lose their conformer, or that come back as several fragments are reported with an error and no affinity rather than being skipped silently — check n_scored in the summary against your input count.
sbdd is not part of --metrics all, because it needs a receptor and costs seconds per molecule (--dock-mode score is far cheaper than dock).

Part 3: Conformer Metrics (paired generated-vs-reference)¶

--metrics conformer is the home for every metric that needs a pair of structures — a generated one and the reference it should correspond to. It replaces the old analyze compare command and the removed --check-strain flag.

It is exclusive: it is not part of --metrics all, because a plain directory of samples has no reference to pair against.

Two accepted input layouts¶

1. SDF pairs — what conformer generation writes:

generated_conformers/
  conformers.csv
  mol_0000/{conformers.sdf, reference.sdf, conformer_000.xyz, ...}

MolCraftDiff generate my_conformer_config.yaml
MolCraftDiff analyze metrics generated_conformers/ --metrics conformer

Both molecules of a pair are built from the same molecular graph, so atom and bond ordering match — which is what makes stereochemistry comparable. This layout is never routed through xyz2mol/OpenBabel: re-perceiving bonds from coordinates would destroy the stereochemistry being measured.

2. xyz + optimised — the layout Part 1 produces:

MolCraftDiff analyze optimize gen_xyz/
MolCraftDiff analyze metrics gen_xyz/ --metrics conformer

Here each structure is perceived from coordinates independently, so atom ordering is not guaranteed to agree and stereo columns are not emitted (the run logs the reason).

Computed metrics¶

Column	Meaning
`rs_ok`, `ez_ok`, `n_stereocentres`	Stereochemistry preserved vs the reference (layout 1 only). All-or-nothing per molecule; the reference’s full enantiomer counts as correct, matching the convention the conformer literature reports. `None` when the reference has no stereochemistry of that kind — never 0.
`rmsd`	Read from `conformers.csv` (heavy-atom, symmetry-corrected), not recomputed
`bond_length_mean`, `bond_angle_mean`, `torsion_angle_mean`	Mean paired deviation, weighted across bond/angle/torsion types
`mmff_strain_kcal`	MMFF94 energy drop on relaxation — always available
`xtb_strain_kcal`	Same at the xTB level. Empty with a warning when the `xtb` binary is not on `PATH`; never a fabricated number

Output¶

conformer_metrics.csv (one row per conformer) plus a conformer_metrics_summary.json with rs_score, ez_score, n_rs_scored, n_skipped, rmsd_median, rmsd_best_per_mol_mean, coverage_at_threshold, and the geometry / strain means.

Part 4: XYZ to SMILES Conversion¶

Convert 3D XYZ files to 2D SMILES and extract molecular fingerprints.

Usage¶

MolCraftDiff analyze xyz2mol gen_xyz/ --bits 2048

Output Files (in `xyz_dir/2d_reprs/`)¶

File	Description
`smiles_processed.csv`	Filename → SMILES mapping
`fingerprints.npy`	Morgan fingerprints array
`scaffolds.txt`	Murcko scaffolds
`substructures.json`	Substructure counts

Part 5: XTB Electronic Properties¶

Compute quantum-chemical descriptors at GFN-xTB level using morfeus.

Usage¶

# Basic energy properties
MolCraftDiff analyze xtb-electronic gen_xyz/ -p energy

# All properties with JSON output
MolCraftDiff analyze xtb-electronic gen_xyz/ -p all -f json -o results.json

# ASE database for downstream analysis
MolCraftDiff analyze xtb-electronic gen_xyz/ -p all -f ase -o results.db

Property Groups¶

Molecular-level:

Group	Properties
`energy`	HOMO, LUMO, HOMO-LUMO gap
`dipole`	Dipole vector and magnitude
`reactivity`	Ionisation potential, electron affinity
`global`	Electrophilicity, nucleophilicity, fugalities
`solvation`	Solvation energy, H-bond correction (requires `--solvent`)

Atomic-level:

Group	Properties
`charges`	Mulliken atomic charges
`fukui`	Fukui indices (f⁺, f⁻, radical, dual)
`bond_orders`	Wiberg bond orders

Output Formats¶

Format	Description
`csv`	Molecular-level properties (one row per molecule)
`json`	Full data including atomic-level properties
`ase`	ASE database with properties in atoms.info/arrays
`all`	Generate all three formats

Options¶

Option	Default	Description
`-m, --method`	`2`	XTB method: 1=GFN1, 2=GFN2, ptb=PTB
`-c, --charge`	`0`	Molecular charge
`--n-unpaired`	`0`	Number of unpaired electrons
`--auto-charge`	False	For PTB neutral singlets with odd electrons, infer ±1 charge from the XYZ
`-s, --solvent`	None	Solvent for the `solvation` group (e.g. `water`, `thf`, `chcl3`)
`-p, --properties`	`energy`	Property groups to compute
`-f, --format`	`csv`	Output format
`--corrected/--no-corrected`	True	Apply empirical IP/EA correction
`--annotate-db`	False	For ASE `.db` input, annotate input rows in place with `xtb_*` results
`-j, --n-jobs`	`1`	Parallel jobs
`-t, --timeout`	`120`	Timeout per molecule (s)

Part 6: Featurise — Fixed-Size Molecular Vectors¶

Convert a directory of XYZ files into a fixed-size feature matrix for downstream machine learning (clustering, regression, dimensionality reduction, etc.).

Three backends are available:

Backend	Description	GPU needed
`soap`	SOAP descriptor via dscribe	No
`uma`	UMA backbone embeddings from pretrained fairchem model	Optional
`ssl3d`	Embeddings from a trained SSL3D checkpoint	Optional

Usage¶

# SOAP (default) — uses the built-in species list
MolCraftDiff analyze featurize gen_xyz/

# SOAP — auto-detect species from the files
MolCraftDiff analyze featurize gen_xyz/ --autodetect

# SOAP — specify species explicitly
MolCraftDiff analyze featurize gen_xyz/ --species C --species H --species N --species O

# SOAP — custom descriptor parameters
MolCraftDiff analyze featurize gen_xyz/ --n-max 12 --l-max 9 --r-cut 8.0

# UMA — CPU inference
MolCraftDiff analyze featurize gen_xyz/ --backend uma --device cpu

# UMA — GPU inference with custom checkpoint
MolCraftDiff analyze featurize gen_xyz/ --backend uma --device cuda \
    --checkpoint training_outputs/uma-s-1p2.pt

# UMA — all spherical components (higher-dimensional embedding)
MolCraftDiff analyze featurize gen_xyz/ --backend uma --all-components

# SSL3D — use a trained SSL3D checkpoint
MolCraftDiff analyze featurize gen_xyz/ --backend ssl3d \
    --ssl3d-checkpoint runs/last.ckpt --device cuda

SOAP Options¶

Option	Default	Description
`--autodetect`	False	Detect element species from files (overrides `--species`)
`--species`	See below	Element symbols; repeatable: `--species C --species H`
`--r-cut`	`6.0`	Cutoff radius in Å
`--n-max`	`8`	Radial basis functions
`--l-max`	`6`	Angular basis functions
`--sigma`	`0.1`	Gaussian smearing width
`--pooling`	`mean`	Atom pooling mode: `mean` or `sum`
`--soap-jobs`	`1`	Parallel workers

Default species list: H B C N O F Al Si P S Cl As Se Br I Hg Bi

Elements found in the files that are not in the species list are added automatically with a warning, so the run never fails silently on unseen atoms.

UMA Options¶

Option	Default	Description
`--checkpoint`	`training_outputs/uma-s-1p2.pt`	Path to UMA checkpoint `.pt` file
`--task-name`	`omol`	UMA task name
`--device`	auto	`cuda` or `cpu`
`--batch-size`	`8`	Molecules per UMA forward pass
`--charge`	`0`	Total molecular charge applied to all structures
`--spin`	`1`	Spin multiplicity applied to all structures
`--all-components`	False	Use all spherical components instead of L=0 scalars only
`--pooling`	`mean`	Atom pooling mode: `mean` or `sum`

The UMA backend requires the vendored fairchem source tree at <repo_root>/fairchem/src. If it is not found, clone it with:

git clone https://github.com/pregHosh/fairchem fairchem

Then run from the repository root, or set:

export MOLCRAFT_REPO_ROOT=/path/to/MolCraftDiffusion

SSL3D Options¶

Option	Default	Description
`--ssl3d-checkpoint`	None	Required path to a trained SSL3D `.ckpt` or `.pkl` checkpoint
`--edge-radius`	`5.0`	Radius graph cutoff in Å
`--device`	auto	`cuda` or `cpu`
`--batch-size`	`8`	Molecules per SSL3D forward pass
`--pooling`	`mean`	Atom pooling mode: `mean` or `sum`

Output Files¶

Three files are written to the output stem (default: input_dir/features):

File	Content
`features.npy`	`(N, D)` float32 feature matrix
`features.csv`	Row index → source file + frame mapping
`features_meta.json`	Backend, all parameters, feature dim, timestamp

# Custom output stem
MolCraftDiff analyze featurize gen_xyz/ -o results/soap_features
# writes: results/soap_features.npy / .csv / _meta.json

Loading the Output¶

import numpy as np
import pandas as pd

features = np.load("gen_xyz/features.npy")        # (N, D) float32
index    = pd.read_csv("gen_xyz/features.csv")    # maps row → file/frame

print(features.shape)   # e.g. (160, 81396) for SOAP, (160, 128) for UMA

Python API¶

from MolecularDiffusion.runmodes.analyze.featurize import run_featurize

# SOAP
features = run_featurize(
    input_dir="gen_xyz/",
    backend="soap",
    output_path="gen_xyz/soap_features",
    n_max=12,
    l_max=9,
)

# UMA — pass a pre-loaded list of ASE Atoms directly
from MolecularDiffusion.runmodes.analyze.uma_embeddings import get_uma_molecule_embeddings
from ase.io import read

atoms_list = [read("gen_xyz/molecule_0000.xyz"), read("gen_xyz/molecule_0001.xyz")]
results = get_uma_molecule_embeddings(
    source=atoms_list,
    checkpoint_path="training_outputs/uma-s-1p2.pt",
    device="cpu",
    charge=0,
    spin=1,
)
mol_emb = results[0]["molecule_embedding"]   # torch.Tensor, shape (128,)
node_emb = results[0]["node_embedding"]      # torch.Tensor, shape (n_atoms, 128)

Example Workflow¶

A typical post-generation analysis workflow:

# 1. Generate molecules
MolCraftDiff generate gen_config.yaml

# 2. Optimise geometries
MolCraftDiff analyze optimize gen_xyz/ -l gfn2 -o gen_xyz/optimized_xyz

# 3. Compute validity metrics
MolCraftDiff analyze metrics gen_xyz/optimized_xyz -o metrics.csv

# 4. Paired generated-vs-optimised metrics
MolCraftDiff analyze metrics gen_xyz/ --metrics conformer

# 5. Convert to SMILES for downstream analysis
MolCraftDiff analyze xyz2mol gen_xyz/optimized_xyz

# 6. Compute electronic properties
MolCraftDiff analyze xtb-electronic gen_xyz/optimized_xyz -p all -f ase -o electronic.db

# 7. Featurise for downstream ML (SOAP)
MolCraftDiff analyze featurize gen_xyz/optimized_xyz -o gen_xyz/soap_features

# 7b. Or use UMA embeddings (requires fairchem checkout + checkpoint)
MolCraftDiff analyze featurize gen_xyz/optimized_xyz --backend uma --device cuda \
    -o gen_xyz/uma_features

# 7c. Or use SSL3D embeddings from a trained SSL3D checkpoint
MolCraftDiff analyze featurize gen_xyz/optimized_xyz --backend ssl3d \
    --ssl3d-checkpoint runs/last.ckpt -o gen_xyz/ssl3d_features

Python API¶

The analysis functions are also available programmatically:

from MolecularDiffusion.runmodes.analyze import (
    optimize_molecule,
    get_xtb_optimized_xyz,
    compute_xtb_electronic,
    batch_xtb_electronic,
    run_xyz2mol,
)
from MolecularDiffusion.runmodes.analyze.conformer_metrics import (
    compute_conformer_metrics,
)
from MolecularDiffusion.runmodes.analyze.featurize import run_featurize
from MolecularDiffusion.runmodes.analyze.uma_embeddings import get_uma_molecule_embeddings

# Compute electronic properties for single file
result = compute_xtb_electronic(
    "molecule.xyz",
    method=2,
    properties=["energy", "charges"]
)
print(result["homo"], result["lumo"])

# Batch processing
df = batch_xtb_electronic(
    input_dir="gen_xyz/",
    output_path="results.csv",
    output_format="csv",
    properties=["energy", "reactivity"],
)

Verify the Result¶

For every analysis command, compare the reported processed count with the number of eligible input structures. Inspect missing values and per-file failures before using aggregate metrics or feature arrays downstream.

Troubleshooting¶

Use MolCraftDiff analyze <subcommand> --help to confirm optional dependencies and accepted formats.
Conversion failures often indicate unsupported elements, invalid coordinates, or disconnected structures.
xTB-backed commands require the xTB executable to be available in the active environment.

Tutorial 9: Analyse Module - 3D Molecular Structure Analysis¶

At a Glance¶

Overview¶

Part 1: XTB Geometry Optimisation¶

Usage¶

Options¶

Output¶

Part 2: Validity & Connectivity Metrics¶

Usage¶

Metric Types¶

What core reports¶

Options¶

What druglike and similarity3d report¶

What sbdd reports¶

Part 3: Conformer Metrics (paired generated-vs-reference)¶

Two accepted input layouts¶

Computed metrics¶

Output¶

Part 4: XYZ to SMILES Conversion¶

Usage¶

Output Files (in xyz_dir/2d_reprs/)¶

Part 5: XTB Electronic Properties¶

Usage¶

Property Groups¶

Output Formats¶

Options¶

Part 6: Featurise — Fixed-Size Molecular Vectors¶

Usage¶

SOAP Options¶

UMA Options¶

SSL3D Options¶

Output Files¶

Loading the Output¶

Python API¶

Example Workflow¶

Python API¶

Verify the Result¶

Troubleshooting¶

What `core` reports¶

What `druglike` and `similarity3d` report¶

What `sbdd` reports¶

Output Files (in `xyz_dir/2d_reprs/`)¶