SITH: Splitting Intramolecular Tension due to StretcHing

SITH is a quantum-chemical Python framework for strain energy decomposition and bond destabilization analysis under mechanical stress. It allows computational chemists and mechanochemists to quantify how external pulling forces distribute energy across individual internal coordinates (bonds, angles, dihedrals) of molecules.

---

📚 Reference

If you use SITH in your research, please cite the following paper:

SITH: A Quantum-Chemical Framework for Predicting Bond Destabilization in Stretched Molecules
Daniel Sucerquia, Mikaela Farrugia, Andreas Dreuw, Frauke Gräter
Submitted on 28 Jan 2026
arXiv:2601.20441 [physics.chem-ph]

---

🌟 Key Features

• Anharmonic Accuracy: Unlike harmonic approximations (e.g. JEDI), SITH numerically integrates internal forces along the stretching pathway, accurately predicting energy accumulation up to bond rupture.
• Internal Coordinate Decomposition: Transforms Cartesian gradients and forces into 3N-6 internal degrees of freedom (bonds, valence angles, and dihedrals).
• COGEF & Stretching Simulations: Automated constrained geometry optimization (COGEF) pipeline for pulling molecules step-by-step.
• Multi-Backend Quantum Chemistry Interfaces: Seamlessly integrates with Psi4, Gaussian, and cclib for parsing energies, forces, and geometries.
• Ring Structure Support: Handles cyclic molecules (e.g. proline, cyclic polymers) with customizable reference coordinate sets.
• Kinetic Barrier Estimation: Integrates with the Bell-Evans model to calculate force-dependent activation barriers and relative rupture rates.
• Comprehensive Visualization: Built-in plotting tools for energy decomposition breakdown, bond-breaking rankings, and force-extension curves.

---

📖 Theoretical Background

When a molecule is stretched by external mechanical force:

1. Energy is accumulated non-uniformly across internal degrees of freedom $q_i$ (bonds $r$, angles $\theta$, dihedrals $\phi$).
2. SITH computes internal forces $F_i(q)$ at each stretching configuration step $k$ by transforming Cartesian forces using the pseudoinverse of Wilson's B-matrix: $$\mathbf{F}{\text{internal}} = (\mathbf{B}^+)^T \mathbf{F}{\text{Cartesian}}$$
3. The total strain energy stored in degree of freedom $i$ is calculated via numerical integration (trapezoidal rule): $$E_i = \int_{q_{i,0}}^{q_{i,\text{final}}} F_i(q_i) , dq_i$$

Because SITH directly integrates forces rather than assuming quadratic potential energy wells ($E = \frac{1}{2} k \Delta q^2$), it avoids overestimating bond energies at large deformations near cleavage.

---

📦 Installation

From PyPI

pip install sith-chem

With Quantum Chemistry Extras

pip install sith-chem[qm]

From Source (Editable Mode)

git clone https://github.com/username/sith.git
cd sith
pip install -e .[dev]

---

🚀 Quick Start

1. Define a Molecule

import numpy as np
from sith.core.molecule import Molecule

# Water molecule Cartesian coordinates (Angstroms)
atoms = ["O", "H", "H"]
coords = np.array([
    [0.0000, 0.0000, 0.1173],
    [0.0000, 0.7572, -0.4692],
    [0.0000, -0.7572, -0.4692]
])

mol = Molecule(atoms=atoms, coordinates=coords)
print(f"Number of atoms: {len(mol)}")
print(f"Internal coordinates (3N-6): {len(mol.degrees_of_freedom)}")

2. Perform SITH Energy Decomposition

from sith.analysis.sith_decomposition import SITHAnalysis
from sith.core.trajectory import StretchingPath

# Initialize SITH analysis on a stretching path
# (StretchingPath contains sequence of molecular configurations & forces)
analysis = SITHAnalysis(trajectory=stretching_path)
decomposition = analysis.decompose()

# Display energy accumulated per internal coordinate
bond_energies = decomposition.get_bond_energies()
for bond, energy in bond_energies.items():
    print(f"Bond {bond}: {energy:.2f} kcal/mol")

# Identify the bond most prone to rupture
ranking = decomposition.get_energy_ranking()
print(f"Highest strain stored in: {ranking[0]}")

3. Estimate Force-Modified Reaction Barriers

from sith.analysis.barriers import BellEvansBarrier

# Calculate barrier reduction under 500 pN pulling force
bell_model = BellEvansBarrier(activation_barrier_0=35.0, transition_state_distance=0.15)
reduced_barrier = bell_model.compute_barrier(force=500.0) # pN
rate_acceleration = bell_model.compute_relative_rate(force=500.0, temperature=298.15)

print(f"Reduced barrier: {reduced_barrier:.2f} kcal/mol")
print(f"Rate acceleration factor: {rate_acceleration:.2e}")

4. Plot Energy Breakdown

from sith.visualization.plots import plot_energy_decomposition

# Generate breakdown bar plot
fig, ax = plot_energy_decomposition(decomposition, unit="kcal/mol")
fig.savefig("sith_energy_decomposition.png", dpi=300)

---

🛠️ Comparison: SITH vs. JEDI

Feature	SITH (This Package)	JEDI (Harmonic)
Deformation Limit	Works up to bond rupture	Valid only near equilibrium ($\Delta q \ll 1$)
Potential Model	Anharmonic (numerical integration)	Harmonic ($E = \frac{1}{2}k\Delta q^2$)
Accuracy at Large Deformations	High (captures force saturation)	Overestimates stored energy
QM Backends	Psi4, Gaussian, cclib	Gaussian

---

🧪 Development & Running Tests

Run the test suite using pytest:

pytest

To run with coverage report:

pytest --cov=sith --cov-report=term-missing

---

📄 License

Distributed under the MIT License. See LICENSE for details.