Strain Energy Analysis

This notebook demonstrates strain energy calculations using Auto3D’s neural network potentials:

1. Ring strain - deviation from ideal bond angles in cyclic compounds

2. Steric strain - unfavorable van der Waals contacts

3. Torsional strain - eclipsing interactions

4. Conformational strain - energy penalty for non-minimum conformers

Chemistry Background

Types of Strain

Type	Cause	Energy Range
Angle strain	Deviation from ideal bond angles	1-25 kcal/mol
Torsional strain	Eclipsed bonds	1-6 kcal/mol
Steric strain	Nonbonded repulsion	1-10 kcal/mol
Ring strain	Combined effects in cycles	1-27 kcal/mol

Classic Ring Strain Values (kcal/mol)

Ring	Strain	Notes
Cyclopropane	~27	Severe angle strain (60° vs 109°)
Cyclobutane	~26	Angle + torsional strain
Cyclopentane	~6	Modest torsional strain
Cyclohexane	~0	Nearly strain-free (chair)
Cycloheptane	~6	Transannular strain

import os
import tempfile
from pathlib import Path

import numpy as np
import Auto3D
from Auto3D import Auto3DOptions, main
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors, rdMolDescriptors
import pandas as pd

print(f"Auto3D version: {Auto3D.__version__}")

HARTREE_TO_KCAL = 627.509

1. Ring Strain: Cycloalkane Series

Compare strain energies across ring sizes using isodesmic reactions.

# Cycloalkane series
cycloalkanes = {
    "cyclopropane": "C1CC1",
    "cyclobutane": "C1CCC1",
    "cyclopentane": "C1CCCC1",
    "cyclohexane": "C1CCCCC1",
    "cycloheptane": "C1CCCCCC1",
    "cyclooctane": "C1CCCCCCC1",
    # Reference acyclic compound (for strain calculation)
    "ethane": "CC",
    "propane": "CCC",
    "butane": "CCCC",
    "pentane": "CCCCC",
    "hexane": "CCCCCC",
}

print("Compounds for ring strain analysis:")
for name, smi in cycloalkanes.items():
    mol = Chem.MolFromSmiles(smi)
    n_c = sum(1 for a in mol.GetAtoms() if a.GetSymbol() == 'C')
    print(f"  {name}: {n_c} carbons")

# Write to file
with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:
    for name, smi in cycloalkanes.items():
        f.write(f"{smi} {name}\n")
    ring_file = f.name

print(f"Input file: {ring_file}")

# Generate optimized structures
if __name__ == "__main__":
    config = Auto3DOptions(
        path=ring_file,
        k=1,
        optimizing_engine="AIMNET",
        use_gpu=True,
        convergence_threshold=0.005,
    )

    ring_output = main(config)
    print(f"Output: {ring_output}")

def extract_energies(sdf_path):
    """Extract energies from Auto3D output."""
    mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))

    energies = {}
    for mol in mols:
        if mol is None:
            continue
        name = mol.GetProp("_Name").split("_")[0]

        if mol.HasProp("E_hartree"):
            e = float(mol.GetProp("E_hartree"))
        elif mol.HasProp("E_tot"):
            e = float(mol.GetProp("E_tot")) / 27.2114  # eV to Hartree
        else:
            continue

        energies[name] = e * HARTREE_TO_KCAL  # Convert to kcal/mol

    return energies


def calculate_ring_strain_homodesmotic(energies):
    """
    Calculate ring strain using homodesmotic reactions.

    Cycloalkane + n*ethane → n*propane

    Strain = E(cycloalkane) + n*E(ethane) - n*E(propane)
    """
    if 'ethane' not in energies or 'propane' not in energies:
        return {}

    e_ethane = energies['ethane']
    e_propane = energies['propane']

    strains = {}
    ring_names = ['cyclopropane', 'cyclobutane', 'cyclopentane',
                  'cyclohexane', 'cycloheptane', 'cyclooctane']
    ring_sizes = [3, 4, 5, 6, 7, 8]

    for name, n in zip(ring_names, ring_sizes):
        if name in energies:
            # Homodesmotic: cyclo-Cn + n*C2H6 → n*C3H8
            strain = energies[name] + n * e_ethane - n * e_propane
            strains[name] = strain

    return strains


if 'ring_output' in dir():
    energies = extract_energies(ring_output)
    strains = calculate_ring_strain_homodesmotic(energies)

    print("\nRing Strain Analysis (Homodesmotic):")
    print("=" * 50)
    print(f"{'Ring':<15} {'Strain (kcal/mol)':>15} {'Per CH2':>12}")
    print("-" * 50)

    ring_sizes = {'cyclopropane': 3, 'cyclobutane': 4, 'cyclopentane': 5,
                  'cyclohexane': 6, 'cycloheptane': 7, 'cyclooctane': 8}

    for name, strain in strains.items():
        n = ring_sizes[name]
        per_ch2 = strain / n
        print(f"{name:<15} {strain:>15.1f} {per_ch2:>12.1f}")
    print("=" * 50)

2. Conformational Strain: Axial vs Equatorial

In substituted cyclohexanes, axial substituents experience steric strain (1,3-diaxial interactions).

# Methylcyclohexane conformers
# The chair flip interconverts axial ↔ equatorial

# A-values (kcal/mol preference for equatorial):
# -CH3: 1.74, -OH: 0.87, -F: 0.15, -Cl: 0.53, -Br: 0.48
# -t-Bu: 4.9 (essentially locked equatorial)

cyclohexane_derivatives = {
    "methylcyclohexane": "CC1CCCCC1",
    "ethylcyclohexane": "CCC1CCCCC1",
    "isopropylcyclohexane": "CC(C)C1CCCCC1",
    "tert_butylcyclohexane": "CC(C)(C)C1CCCCC1",
    "phenylcyclohexane": "C1CCCCC1C2=CC=CC=C2",
}

with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:
    for name, smi in cyclohexane_derivatives.items():
        f.write(f"{smi} {name}\n")
    deriv_file = f.name

print("Cyclohexane derivatives for conformational analysis")

# Generate multiple conformers to capture axial vs equatorial
if __name__ == "__main__":
    config = Auto3DOptions(
        path=deriv_file,
        k=5,                    # Get multiple conformers
        window=5.0,             # Include higher energy axial forms
        optimizing_engine="AIMNET",
        use_gpu=True,
    )

    deriv_output = main(config)
    print(f"Output: {deriv_output}")

def analyze_conformational_strain(sdf_path):
    """
    Analyze conformational energy differences.

    For cyclohexane derivatives, this approximates A-values.
    """
    mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))

    # Group by molecule
    from collections import defaultdict
    by_mol = defaultdict(list)

    for mol in mols:
        if mol is None:
            continue
        name_full = mol.GetProp("_Name")
        name = name_full.split("_")[0]

        if mol.HasProp("E_rel"):
            e = float(mol.GetProp("E_rel"))
        elif mol.HasProp("E_hartree"):
            e = float(mol.GetProp("E_hartree")) * HARTREE_TO_KCAL
        else:
            continue

        by_mol[name].append(e)

    # Calculate energy range for each molecule
    results = []
    for name, energies in by_mol.items():
        energies = np.array(energies)
        e_range = energies.max() - energies.min()
        results.append({
            "Molecule": name,
            "n_conformers": len(energies),
            "E_range": e_range,
            "E_min": energies.min(),
        })

    return pd.DataFrame(results)


if 'deriv_output' in dir():
    df_conf = analyze_conformational_strain(deriv_output)
    print("\nConformational Energy Analysis:")
    print(df_conf.to_string(index=False))
    print("\nE_range approximates the axial-equatorial energy difference.")

3. Steric Strain: Gauche Interactions

Butane gauche interactions are a classic example of steric strain.

def measure_torsion_angle(mol, atom_indices):
    """
    Measure dihedral angle for 4 atoms.
    """
    conf = mol.GetConformer()
    return Chem.rdMolTransforms.GetDihedralDeg(conf, *atom_indices)


# For butane: anti (180°) vs gauche (60°) vs syn (0°)
# gauche-butane interaction ≈ 0.9 kcal/mol

print("Butane conformational analysis:")
print("  Anti (180°): most stable")
print("  Gauche (±60°): +0.9 kcal/mol")
print("  Eclipsed (±120°): +3.4 kcal/mol (transition state)")
print("  Syn (0°): +4.5 kcal/mol (highest energy)")

4. Angle Strain Analysis

Compare bond angles to ideal values.

def analyze_bond_angles(mol):
    """
    Analyze bond angles and their deviation from ideal.

    Ideal angles:
    - sp3 carbon: 109.5°
    - sp2 carbon: 120°
    - sp carbon: 180°
    """
    from rdkit.Chem import rdMolTransforms

    conf = mol.GetConformer()

    angles = []
    for atom in mol.GetAtoms():
        idx = atom.GetIdx()
        neighbors = [n.GetIdx() for n in atom.GetNeighbors()]

        if len(neighbors) >= 2:
            # Calculate all angles around this atom
            for i in range(len(neighbors)):
                for j in range(i+1, len(neighbors)):
                    angle = rdMolTransforms.GetAngleDeg(
                        conf, neighbors[i], idx, neighbors[j]
                    )

                    # Determine ideal angle based on hybridization
                    hybridization = atom.GetHybridization()
                    if hybridization == Chem.HybridizationType.SP3:
                        ideal = 109.5
                    elif hybridization == Chem.HybridizationType.SP2:
                        ideal = 120.0
                    elif hybridization == Chem.HybridizationType.SP:
                        ideal = 180.0
                    else:
                        ideal = 109.5  # Default

                    deviation = abs(angle - ideal)
                    angles.append({
                        "atom": atom.GetSymbol(),
                        "idx": idx,
                        "angle": angle,
                        "ideal": ideal,
                        "deviation": deviation
                    })

    return pd.DataFrame(angles)


# Analyze cyclopropane angles
if 'ring_output' in dir():
    mols = list(Chem.SDMolSupplier(ring_output, removeHs=False))

    for mol in mols:
        if mol is None:
            continue
        name = mol.GetProp("_Name").split("_")[0]
        if name == "cyclopropane":
            df_angles = analyze_bond_angles(mol)
            # Filter to carbon-centered angles
            df_c = df_angles[df_angles["atom"] == "C"]

            print(f"\nAngle Analysis for {name}:")
            print(f"  Average C-C-C angle: {df_c['angle'].mean():.1f}°")
            print(f"  Ideal sp3 angle: 109.5°")
            print(f"  Average deviation: {df_c['deviation'].mean():.1f}°")
            break

5. Drug Molecule Strain

Analyze strain in drug-like molecules - important for binding affinity.

# Drugs with strained motifs
strained_drugs = {
    # β-lactam antibiotics (4-membered ring)
    "penicillin_core": "O=C1CCN1",  # Azetidinone

    # Cyclopropane-containing drugs
    "cipro_cyclopropyl": "C1CC1N",  # Cyclopropylamine moiety

    # Cubane (extreme strain, research compound)
    "cubane": "C12C3C4C1C5C4C3C25",

    # Bridged bicyclics
    "norbornane": "C1CC2CCC1C2",  # Bicyclo[2.2.1]heptane
    "adamantane": "C1C2CC3CC1CC(C2)C3",
}

with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:
    for name, smi in strained_drugs.items():
        f.write(f"{smi} {name}\n")
    strain_file = f.name

if __name__ == "__main__":
    config = Auto3DOptions(
        path=strain_file,
        k=1,
        optimizing_engine="AIMNET",
        use_gpu=True,
    )

    strain_output = main(config)
    print(f"Output: {strain_output}")

def calculate_strain_per_atom(sdf_path, reference_energy_per_ch2=-9.77):
    """
    Calculate relative strain using energy per CH2 unit.

    This is a simplified measure - cyclohexane is approximately strain-free.
    """
    mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))

    results = []
    for mol in mols:
        if mol is None:
            continue
        name = mol.GetProp("_Name").split("_")[0]

        if mol.HasProp("E_hartree"):
            e = float(mol.GetProp("E_hartree")) * HARTREE_TO_KCAL
        else:
            continue

        n_heavy = mol.GetNumHeavyAtoms()
        e_per_atom = e / n_heavy

        results.append({
            "Molecule": name,
            "n_heavy": n_heavy,
            "E_total": e,
            "E_per_atom": e_per_atom,
        })

    return pd.DataFrame(results)


if 'strain_output' in dir():
    df_strain = calculate_strain_per_atom(strain_output)
    print("\nStrained Molecule Analysis:")
    print(df_strain.round(2).to_string(index=False))

6. Binding Strain

When a drug binds to a receptor, it may adopt a strained conformation. The energy penalty is called binding strain.

def estimate_binding_strain(sdf_path, bioactive_conformer_idx=None):
    """
    Estimate binding strain as the energy difference between
    the global minimum and a higher-energy conformer.

    In practice, the bioactive conformer would be extracted from
    a crystal structure or docking pose.
    """
    mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))

    energies = []
    for mol in mols:
        if mol is None:
            continue
        if mol.HasProp("E_rel"):
            energies.append(float(mol.GetProp("E_rel")))

    if len(energies) > 1:
        e_min = min(energies)
        e_max = max(energies)
        return {
            "n_conformers": len(energies),
            "E_min": e_min,
            "E_max": e_max,
            "max_binding_strain": e_max - e_min
        }

    return None


print("\nBinding Strain Estimation:")
print("-" * 50)
print("Binding strain = E(bioactive) - E(global minimum)")
print("\nTypical values:")
print("  Low strain: 0-2 kcal/mol")
print("  Moderate strain: 2-5 kcal/mol")
print("  High strain: >5 kcal/mol (unfavorable)")
print("-" * 50)
print("\nHigh binding strain reduces binding affinity:")
print("  ΔG = ΔG_intrinsic + ΔG_strain")

7. Summary

This tutorial covered:

1. Ring strain - cycloalkanes, homodesmotic reactions

2. Conformational strain - axial vs equatorial in cyclohexanes

3. Steric strain - gauche interactions

4. Angle strain - deviation from ideal geometries

5. Binding strain - energy penalty for bioactive conformers

Key insights:

• Strain affects stability, reactivity, and binding

• Neural network potentials capture strain accurately

• Drug design should minimize unnecessary strain

# Cleanup
for f in [ring_file, deriv_file, strain_file]:
    if f in dir() and os.path.exists(f):
        os.unlink(f)