Bond Dissociation Energy Skill

Goal

Calculate the homolytic and/or heterolytic bond dissociation energy (BDE) for each single bond in a molecule using Machine Learning Interatomic Potentials (MLIPs).

Homolytic BDE (radical fragments): $$\text{BDE}_\text{homo}(A{-}B) = E(A\bullet) + E(B\bullet) - E(A{-}B)$$

Heterolytic BDE (ionic fragments, minimum over both polarity variants): $$\text{BDE}_\text{hetero}(A{-}B) = \min!\bigl(E(A^+)+E(B^-),; E(A^-)+E(B^+)\bigr) - E(A{-}B)$$

[!IMPORTANT] This skill computes BDEs by relaxing both the intact molecule and fragments with an MLIP. For purpose-trained GNN models that predict BDE directly from SMILES (MAE ~0.6 kcal/mol), consider ALFABET or BonDNet instead.

Background

BDE is a fundamental thermodynamic quantity that determines:

Drug metabolism: CYP450 enzymes abstract H from the weakest C–H bond
Electrolyte stability: Which bonds break first under electrochemical voltage
Combustion chemistry: Rate-determining bond-breaking steps in fuel oxidation
Polymer degradation: Weakest links in polymer backbone chains

A 2024 study (Zubatyuk et al., JCTC) demonstrated that MACE potentials achieve BDE RMSE of 1.37 kcal/mol for aliphatic C–H bonds in drug-like molecules, outperforming semi-empirical methods and ALFABET for BDE ranking.

1. Prerequisites

Conda Environment: mace-agent (includes RDKit, ASE, and MACE)
Input: SMILES string or structure file (.sdf, .mol2)
RDKit: Required for bond identification and molecular fragmentation

2. Choosing a Foundation Potential

Refer to the foundation-potentials skill for model selection.

[!IMPORTANT] Model requirements by cleavage mode:

| Mode | Recommended model | supports_charge_spin | Validated? | |:---|:---|:---|:---| | homolytic | MACE-OFF23-small/medium/large | Not required | ✅ | | heterolytic or both | MACE-OMOL-extra-large (env: mace-agent) | ✅ Required | ✅ | | heterolytic or both | MACE-MH-1 with omol head (env: mace-agent) | ✅ Required | ✅ | | heterolytic or both | FairChem uma-s-1p1 with --task_name omol (env: fairchem-agent) | ✅ Required | ✅ |

Setting charge/spin on MACE models: use atoms.info["charge"] and atoms.info["spin"] (the calculator's default info_keys maps "charge" → total_charge / "spin" → total_spin). Both MACE-OMOL and MACE-MH use joint_embedding to condition the network on these scalars.

If you request --cleavage both with a model that does not support charge/spin, the skill will log a warning and silently fall back to homolytic-only. Using --cleavage heterolytic with an unsupported model raises an error.

Note on single-atom fragments: When a bond produces a bare H (or other single atom), heterolytic BDE is automatically skipped — neither MACE nor FairChem UMA has signed single-atom energies (only neutral H, C, N, O… are in the reference tables).

3. Calculation Workflow

Step 1: Provide a molecule

# SMILES input (most common)
--smiles "CCO"

# Or from a structure file
--structure molecule.sdf

Step 2: Run BDE calculation

Homolytic only (default, no charge/spin needed):

# Env: mace-agent
python .agents/skills/chem-bond-dissociation/scripts/calculate_bde.py \
    --smiles CCO \
    --all_bonds \
    --cleavage homolytic \
    --model_type mace \
    --model_name MACE-OFF23-small \
    --output_dir research/my_folder/bde_results

Both homolytic and heterolytic (MACE-OMOL):

# Env: mace-agent
python .agents/skills/chem-bond-dissociation/scripts/calculate_bde.py \
    --smiles CCO \
    --all_bonds \
    --cleavage both \
    --model_type mace \
    --model_name MACE-OMOL-extra-large \
    --output_dir research/my_folder/bde_results_both

Both homolytic and heterolytic (FairChem UMA omol):

# Env: fairchem-agent
python .agents/skills/chem-bond-dissociation/scripts/calculate_bde.py \
    --smiles CCO \
    --all_bonds \
    --cleavage both \
    --model_type fairchem \
    --model_name uma-s-1p1 \
    --task_name omol \
    --output_dir research/my_folder/bde_results_both

Heterolytic only with FairChem UMA:

# Env: fairchem-agent
python .agents/skills/chem-bond-dissociation/scripts/calculate_bde.py \
    --smiles CCO \
    --all_bonds \
    --cleavage heterolytic \
    --model_type fairchem \
    --model_name uma-s-1p1 \
    --task_name omol \
    --output_dir research/my_folder/bde_hetero

Key Parameters

| Argument | Default | Description | |:---|:---|:---| | --smiles | — | SMILES string of the molecule | | --structure | — | Path to structure file (.sdf, .mol2) | | --bond | — | Specific bond as atom indices "i-j" (0-indexed) | | --all_bonds | True | Compute BDE for all single bonds | | --include_h_bonds | False | Include X–H bonds | | --cleavage | homolytic | homolytic, heterolytic, or both | | --model_type | mace | MLIP backend (mace, fairchem) | | --model_name | auto | Model checkpoint (default: MACE-OFF23-small for homolytic; uma-s-1p1 for hetero/both) | | --task_name | — | Task head for multi-task models (e.g. omol for FairChem UMA) | | --fmax | 0.01 | Force convergence for relaxation (eV/Å) | | --output_dir | required | Output directory |

4. Output Files

bde_results.json — Full results including:
- metadata: model name, cleavage mode, supports_charge_spin, SMILES, etc.
- intact_energy_eV: Energy of the relaxed intact molecule
- bonds: List of per-bond results:
  - bde_eV, bde_kJ_mol, bde_kcal_mol: Homolytic BDE (if computed)
  - heterolytic_bde_eV, heterolytic_bde_kJ_mol, heterolytic_bde_kcal_mol: Best heterolytic BDE (if computed)
  - heterolytic_best_variant: Which polarity won ("frag1+ / frag2-" or "frag1- / frag2+")
  - heterolytic_variants: Raw results for both polarity variants
- weakest_bond_homolytic, weakest_bond_heterolytic: Summary of weakest bonds
- bonds_ranked_by_homolytic_bde, bonds_ranked_by_heterolytic_bde: Sorted tables
intact_relaxed.xyz: Relaxed intact molecule
frag_bond{N}_homo_{1,2}.xyz: Homolytic radical fragments
frag_bond{N}_hetero_pos_neg_{1,2}.xyz: Heterolytic cation/anion fragments (variant A)
frag_bond{N}_hetero_neg_pos_{1,2}.xyz: Heterolytic anion/cation fragments (variant B)

5. Examples

| Example | Model | Cleavage | Notes | |:---|:---|:---|:---| | examples/ethanol_mace_off23_small/ | MACE-OFF23-small | homolytic | Standard homolytic BDE for ethanol; includes H bonds | | examples/methanol_mace_omol_both/ | MACE-OMOL-extra-large | both | Homo + heterolytic for methanol; C–O hetero = 90 vs homo = 143 kcal/mol | | examples/methanol_uma_omol_both/ | FairChem UMA omol | both | Homo + heterolytic for methanol; C–O hetero = 157 vs homo = 126 kcal/mol |

Ethanol — Homolytic BDE (MACE-OFF23)

# Env: mace-agent
python .agents/skills/chem-bond-dissociation/scripts/calculate_bde.py \
    --smiles CCO \
    --all_bonds \
    --include_h_bonds \
    --cleavage homolytic \
    --model_type mace \
    --model_name MACE-OFF23-small \
    --output_dir .agents/skills/chem-bond-dissociation/examples/ethanol_mace_off23_small

Methanol — Both Homo and Heterolytic BDE (FairChem UMA omol)

# Env: fairchem-agent
python .agents/skills/chem-bond-dissociation/scripts/calculate_bde.py \
    --smiles CO \
    --all_bonds \
    --include_h_bonds \
    --cleavage both \
    --model_type fairchem \
    --model_name uma-s-1p1 \
    --task_name omol \
    --output_dir .agents/skills/chem-bond-dissociation/examples/methanol_uma_omol_both

Experimental BDEs for ethanol (Blanksby & Ellison, 2003): | Bond | Experimental BDE (kcal/mol) | |:---|:---| | O–H | ~104 | | C–H (methyl) | ~101 | | C–H (methylene) | ~95 | | C–C | ~85 | | C–O | ~92 |

6. Constraints

Radical spin states: For homolytic BDE, MLIPs are generally "electron-agnostic" and treat fragments as neutral regardless of spin state. BDE ranking is typically more reliable than absolute values.
Ionic states: Heterolytic BDE requires a charge/spin-aware model (supports_charge_spin=True). Validated: MACE-OMOL, MACE-MH (omol head), and FairChem UMA omol. These models use atoms.info["charge"] and atoms.info["spin"] to condition on the ionic state. Models without this flag raise an error for --cleavage heterolytic.
Ring bonds: Breaking bonds in rings produces a single open-chain diradical. The script will warn and skip ring bonds.
Accuracy: Expect ~2–5 kcal/mol error for homolytic BDEs with MACE-OFF23. Heterolytic accuracy is less benchmarked with current MLIPs.
Environments:
- mace-agent for MACE models
- fairchem-agent for FairChem/UMA models

References

Blanksby & Ellison, "Bond Dissociation Energies of Organic Molecules", Acc. Chem. Res. 2003, 36, 255.
St. John et al., "Prediction of organic homolytic bond dissociation enthalpies at near chemical accuracy with sub-second computational cost", Nat. Commun. 2020, 11, 2328. (ALFABET)
Zubatyuk et al., "A Transferable MACE Potential for Open- and Closed-Shell Drug-Like Molecules", J. Chem. Theory Comput. 2024.

Author: Bowen Deng Contact: GitHub @learningmatter-mit