{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Strain Energy Analysis\n", "\n", "This notebook demonstrates strain energy calculations using Auto3D's neural network potentials:\n", "\n", "1. **Ring strain** - deviation from ideal bond angles in cyclic compounds\n", "2. **Steric strain** - unfavorable van der Waals contacts\n", "3. **Torsional strain** - eclipsing interactions\n", "4. **Conformational strain** - energy penalty for non-minimum conformers\n", "\n", "## Chemistry Background\n", "\n", "### Types of Strain\n", "\n", "| Type | Cause | Energy Range |\n", "|------|-------|-------------|\n", "| **Angle strain** | Deviation from ideal bond angles | 1-25 kcal/mol |\n", "| **Torsional strain** | Eclipsed bonds | 1-6 kcal/mol |\n", "| **Steric strain** | Nonbonded repulsion | 1-10 kcal/mol |\n", "| **Ring strain** | Combined effects in cycles | 1-27 kcal/mol |\n", "\n", "### Classic Ring Strain Values (kcal/mol)\n", "\n", "| Ring | Strain | Notes |\n", "|------|--------|-------|\n", "| Cyclopropane | ~27 | Severe angle strain (60° vs 109°) |\n", "| Cyclobutane | ~26 | Angle + torsional strain |\n", "| Cyclopentane | ~6 | Modest torsional strain |\n", "| Cyclohexane | ~0 | Nearly strain-free (chair) |\n", "| Cycloheptane | ~6 | Transannular strain |" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "import tempfile\n", "from pathlib import Path\n", "\n", "import numpy as np\n", "import Auto3D\n", "from Auto3D import Auto3DOptions, main\n", "from rdkit import Chem\n", "from rdkit.Chem import AllChem, Descriptors, rdMolDescriptors\n", "import pandas as pd\n", "\n", "print(f\"Auto3D version: {Auto3D.__version__}\")\n", "\n", "HARTREE_TO_KCAL = 627.509" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Ring Strain: Cycloalkane Series\n", "\n", "Compare strain energies across ring sizes using isodesmic reactions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Cycloalkane series\n", "cycloalkanes = {\n", " \"cyclopropane\": \"C1CC1\",\n", " \"cyclobutane\": \"C1CCC1\",\n", " \"cyclopentane\": \"C1CCCC1\",\n", " \"cyclohexane\": \"C1CCCCC1\",\n", " \"cycloheptane\": \"C1CCCCCC1\",\n", " \"cyclooctane\": \"C1CCCCCCC1\",\n", " # Reference acyclic compound (for strain calculation)\n", " \"ethane\": \"CC\",\n", " \"propane\": \"CCC\",\n", " \"butane\": \"CCCC\",\n", " \"pentane\": \"CCCCC\",\n", " \"hexane\": \"CCCCCC\",\n", "}\n", "\n", "print(\"Compounds for ring strain analysis:\")\n", "for name, smi in cycloalkanes.items():\n", " mol = Chem.MolFromSmiles(smi)\n", " n_c = sum(1 for a in mol.GetAtoms() if a.GetSymbol() == 'C')\n", " print(f\" {name}: {n_c} carbons\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Write to file\n", "with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:\n", " for name, smi in cycloalkanes.items():\n", " f.write(f\"{smi} {name}\\n\")\n", " ring_file = f.name\n", "\n", "print(f\"Input file: {ring_file}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Generate optimized structures\n", "if __name__ == \"__main__\":\n", " config = Auto3DOptions(\n", " path=ring_file,\n", " k=1,\n", " optimizing_engine=\"AIMNET\",\n", " use_gpu=True,\n", " convergence_threshold=0.005,\n", " )\n", " \n", " ring_output = main(config)\n", " print(f\"Output: {ring_output}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def extract_energies(sdf_path):\n", " \"\"\"Extract energies from Auto3D output.\"\"\"\n", " mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))\n", " \n", " energies = {}\n", " for mol in mols:\n", " if mol is None:\n", " continue\n", " name = mol.GetProp(\"_Name\").split(\"_\")[0]\n", " \n", " if mol.HasProp(\"E_hartree\"):\n", " e = float(mol.GetProp(\"E_hartree\"))\n", " elif mol.HasProp(\"E_tot\"):\n", " e = float(mol.GetProp(\"E_tot\")) / 27.2114 # eV to Hartree\n", " else:\n", " continue\n", " \n", " energies[name] = e * HARTREE_TO_KCAL # Convert to kcal/mol\n", " \n", " return energies\n", "\n", "\n", "def calculate_ring_strain_homodesmotic(energies):\n", " \"\"\"\n", " Calculate ring strain using homodesmotic reactions.\n", " \n", " Cycloalkane + n*ethane → n*propane\n", " \n", " Strain = E(cycloalkane) + n*E(ethane) - n*E(propane)\n", " \"\"\"\n", " if 'ethane' not in energies or 'propane' not in energies:\n", " return {}\n", " \n", " e_ethane = energies['ethane']\n", " e_propane = energies['propane']\n", " \n", " strains = {}\n", " ring_names = ['cyclopropane', 'cyclobutane', 'cyclopentane', \n", " 'cyclohexane', 'cycloheptane', 'cyclooctane']\n", " ring_sizes = [3, 4, 5, 6, 7, 8]\n", " \n", " for name, n in zip(ring_names, ring_sizes):\n", " if name in energies:\n", " # Homodesmotic: cyclo-Cn + n*C2H6 → n*C3H8\n", " strain = energies[name] + n * e_ethane - n * e_propane\n", " strains[name] = strain\n", " \n", " return strains\n", "\n", "\n", "if 'ring_output' in dir():\n", " energies = extract_energies(ring_output)\n", " strains = calculate_ring_strain_homodesmotic(energies)\n", " \n", " print(\"\\nRing Strain Analysis (Homodesmotic):\")\n", " print(\"=\" * 50)\n", " print(f\"{'Ring':<15} {'Strain (kcal/mol)':>15} {'Per CH2':>12}\")\n", " print(\"-\" * 50)\n", " \n", " ring_sizes = {'cyclopropane': 3, 'cyclobutane': 4, 'cyclopentane': 5,\n", " 'cyclohexane': 6, 'cycloheptane': 7, 'cyclooctane': 8}\n", " \n", " for name, strain in strains.items():\n", " n = ring_sizes[name]\n", " per_ch2 = strain / n\n", " print(f\"{name:<15} {strain:>15.1f} {per_ch2:>12.1f}\")\n", " print(\"=\" * 50)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Conformational Strain: Axial vs Equatorial\n", "\n", "In substituted cyclohexanes, axial substituents experience steric strain (1,3-diaxial interactions)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Methylcyclohexane conformers\n", "# The chair flip interconverts axial ↔ equatorial\n", "\n", "# A-values (kcal/mol preference for equatorial):\n", "# -CH3: 1.74, -OH: 0.87, -F: 0.15, -Cl: 0.53, -Br: 0.48\n", "# -t-Bu: 4.9 (essentially locked equatorial)\n", "\n", "cyclohexane_derivatives = {\n", " \"methylcyclohexane\": \"CC1CCCCC1\",\n", " \"ethylcyclohexane\": \"CCC1CCCCC1\",\n", " \"isopropylcyclohexane\": \"CC(C)C1CCCCC1\",\n", " \"tert_butylcyclohexane\": \"CC(C)(C)C1CCCCC1\",\n", " \"phenylcyclohexane\": \"C1CCCCC1C2=CC=CC=C2\",\n", "}\n", "\n", "with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:\n", " for name, smi in cyclohexane_derivatives.items():\n", " f.write(f\"{smi} {name}\\n\")\n", " deriv_file = f.name\n", "\n", "print(\"Cyclohexane derivatives for conformational analysis\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Generate multiple conformers to capture axial vs equatorial\n", "if __name__ == \"__main__\":\n", " config = Auto3DOptions(\n", " path=deriv_file,\n", " k=5, # Get multiple conformers\n", " window=5.0, # Include higher energy axial forms\n", " optimizing_engine=\"AIMNET\",\n", " use_gpu=True,\n", " )\n", " \n", " deriv_output = main(config)\n", " print(f\"Output: {deriv_output}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def analyze_conformational_strain(sdf_path):\n", " \"\"\"\n", " Analyze conformational energy differences.\n", " \n", " For cyclohexane derivatives, this approximates A-values.\n", " \"\"\"\n", " mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))\n", " \n", " # Group by molecule\n", " from collections import defaultdict\n", " by_mol = defaultdict(list)\n", " \n", " for mol in mols:\n", " if mol is None:\n", " continue\n", " name_full = mol.GetProp(\"_Name\")\n", " name = name_full.split(\"_\")[0]\n", " \n", " if mol.HasProp(\"E_rel\"):\n", " e = float(mol.GetProp(\"E_rel\"))\n", " elif mol.HasProp(\"E_hartree\"):\n", " e = float(mol.GetProp(\"E_hartree\")) * HARTREE_TO_KCAL\n", " else:\n", " continue\n", " \n", " by_mol[name].append(e)\n", " \n", " # Calculate energy range for each molecule\n", " results = []\n", " for name, energies in by_mol.items():\n", " energies = np.array(energies)\n", " e_range = energies.max() - energies.min()\n", " results.append({\n", " \"Molecule\": name,\n", " \"n_conformers\": len(energies),\n", " \"E_range\": e_range,\n", " \"E_min\": energies.min(),\n", " })\n", " \n", " return pd.DataFrame(results)\n", "\n", "\n", "if 'deriv_output' in dir():\n", " df_conf = analyze_conformational_strain(deriv_output)\n", " print(\"\\nConformational Energy Analysis:\")\n", " print(df_conf.to_string(index=False))\n", " print(\"\\nE_range approximates the axial-equatorial energy difference.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Steric Strain: Gauche Interactions\n", "\n", "Butane gauche interactions are a classic example of steric strain." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def measure_torsion_angle(mol, atom_indices):\n", " \"\"\"\n", " Measure dihedral angle for 4 atoms.\n", " \"\"\"\n", " conf = mol.GetConformer()\n", " return Chem.rdMolTransforms.GetDihedralDeg(conf, *atom_indices)\n", "\n", "\n", "# For butane: anti (180°) vs gauche (60°) vs syn (0°)\n", "# gauche-butane interaction ≈ 0.9 kcal/mol\n", "\n", "print(\"Butane conformational analysis:\")\n", "print(\" Anti (180°): most stable\")\n", "print(\" Gauche (±60°): +0.9 kcal/mol\")\n", "print(\" Eclipsed (±120°): +3.4 kcal/mol (transition state)\")\n", "print(\" Syn (0°): +4.5 kcal/mol (highest energy)\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Angle Strain Analysis\n", "\n", "Compare bond angles to ideal values." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def analyze_bond_angles(mol):\n", " \"\"\"\n", " Analyze bond angles and their deviation from ideal.\n", " \n", " Ideal angles:\n", " - sp3 carbon: 109.5°\n", " - sp2 carbon: 120°\n", " - sp carbon: 180°\n", " \"\"\"\n", " from rdkit.Chem import rdMolTransforms\n", " \n", " conf = mol.GetConformer()\n", " \n", " angles = []\n", " for atom in mol.GetAtoms():\n", " idx = atom.GetIdx()\n", " neighbors = [n.GetIdx() for n in atom.GetNeighbors()]\n", " \n", " if len(neighbors) >= 2:\n", " # Calculate all angles around this atom\n", " for i in range(len(neighbors)):\n", " for j in range(i+1, len(neighbors)):\n", " angle = rdMolTransforms.GetAngleDeg(\n", " conf, neighbors[i], idx, neighbors[j]\n", " )\n", " \n", " # Determine ideal angle based on hybridization\n", " hybridization = atom.GetHybridization()\n", " if hybridization == Chem.HybridizationType.SP3:\n", " ideal = 109.5\n", " elif hybridization == Chem.HybridizationType.SP2:\n", " ideal = 120.0\n", " elif hybridization == Chem.HybridizationType.SP:\n", " ideal = 180.0\n", " else:\n", " ideal = 109.5 # Default\n", " \n", " deviation = abs(angle - ideal)\n", " angles.append({\n", " \"atom\": atom.GetSymbol(),\n", " \"idx\": idx,\n", " \"angle\": angle,\n", " \"ideal\": ideal,\n", " \"deviation\": deviation\n", " })\n", " \n", " return pd.DataFrame(angles)\n", "\n", "\n", "# Analyze cyclopropane angles\n", "if 'ring_output' in dir():\n", " mols = list(Chem.SDMolSupplier(ring_output, removeHs=False))\n", " \n", " for mol in mols:\n", " if mol is None:\n", " continue\n", " name = mol.GetProp(\"_Name\").split(\"_\")[0]\n", " if name == \"cyclopropane\":\n", " df_angles = analyze_bond_angles(mol)\n", " # Filter to carbon-centered angles\n", " df_c = df_angles[df_angles[\"atom\"] == \"C\"]\n", " \n", " print(f\"\\nAngle Analysis for {name}:\")\n", " print(f\" Average C-C-C angle: {df_c['angle'].mean():.1f}°\")\n", " print(f\" Ideal sp3 angle: 109.5°\")\n", " print(f\" Average deviation: {df_c['deviation'].mean():.1f}°\")\n", " break" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Drug Molecule Strain\n", "\n", "Analyze strain in drug-like molecules - important for binding affinity." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Drugs with strained motifs\n", "strained_drugs = {\n", " # β-lactam antibiotics (4-membered ring)\n", " \"penicillin_core\": \"O=C1CCN1\", # Azetidinone\n", " \n", " # Cyclopropane-containing drugs\n", " \"cipro_cyclopropyl\": \"C1CC1N\", # Cyclopropylamine moiety\n", " \n", " # Cubane (extreme strain, research compound)\n", " \"cubane\": \"C12C3C4C1C5C4C3C25\",\n", " \n", " # Bridged bicyclics\n", " \"norbornane\": \"C1CC2CCC1C2\", # Bicyclo[2.2.1]heptane\n", " \"adamantane\": \"C1C2CC3CC1CC(C2)C3\",\n", "}\n", "\n", "with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:\n", " for name, smi in strained_drugs.items():\n", " f.write(f\"{smi} {name}\\n\")\n", " strain_file = f.name" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if __name__ == \"__main__\":\n", " config = Auto3DOptions(\n", " path=strain_file,\n", " k=1,\n", " optimizing_engine=\"AIMNET\",\n", " use_gpu=True,\n", " )\n", " \n", " strain_output = main(config)\n", " print(f\"Output: {strain_output}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def calculate_strain_per_atom(sdf_path, reference_energy_per_ch2=-9.77):\n", " \"\"\"\n", " Calculate relative strain using energy per CH2 unit.\n", " \n", " This is a simplified measure - cyclohexane is approximately strain-free.\n", " \"\"\"\n", " mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))\n", " \n", " results = []\n", " for mol in mols:\n", " if mol is None:\n", " continue\n", " name = mol.GetProp(\"_Name\").split(\"_\")[0]\n", " \n", " if mol.HasProp(\"E_hartree\"):\n", " e = float(mol.GetProp(\"E_hartree\")) * HARTREE_TO_KCAL\n", " else:\n", " continue\n", " \n", " n_heavy = mol.GetNumHeavyAtoms()\n", " e_per_atom = e / n_heavy\n", " \n", " results.append({\n", " \"Molecule\": name,\n", " \"n_heavy\": n_heavy,\n", " \"E_total\": e,\n", " \"E_per_atom\": e_per_atom,\n", " })\n", " \n", " return pd.DataFrame(results)\n", "\n", "\n", "if 'strain_output' in dir():\n", " df_strain = calculate_strain_per_atom(strain_output)\n", " print(\"\\nStrained Molecule Analysis:\")\n", " print(df_strain.round(2).to_string(index=False))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Binding Strain\n", "\n", "When a drug binds to a receptor, it may adopt a strained conformation. The energy penalty is called **binding strain**." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def estimate_binding_strain(sdf_path, bioactive_conformer_idx=None):\n", " \"\"\"\n", " Estimate binding strain as the energy difference between\n", " the global minimum and a higher-energy conformer.\n", " \n", " In practice, the bioactive conformer would be extracted from\n", " a crystal structure or docking pose.\n", " \"\"\"\n", " mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))\n", " \n", " energies = []\n", " for mol in mols:\n", " if mol is None:\n", " continue\n", " if mol.HasProp(\"E_rel\"):\n", " energies.append(float(mol.GetProp(\"E_rel\")))\n", " \n", " if len(energies) > 1:\n", " e_min = min(energies)\n", " e_max = max(energies)\n", " return {\n", " \"n_conformers\": len(energies),\n", " \"E_min\": e_min,\n", " \"E_max\": e_max,\n", " \"max_binding_strain\": e_max - e_min\n", " }\n", " \n", " return None\n", "\n", "\n", "print(\"\\nBinding Strain Estimation:\")\n", "print(\"-\" * 50)\n", "print(\"Binding strain = E(bioactive) - E(global minimum)\")\n", "print(\"\\nTypical values:\")\n", "print(\" Low strain: 0-2 kcal/mol\")\n", "print(\" Moderate strain: 2-5 kcal/mol\") \n", "print(\" High strain: >5 kcal/mol (unfavorable)\")\n", "print(\"-\" * 50)\n", "print(\"\\nHigh binding strain reduces binding affinity:\")\n", "print(\" ΔG = ΔG_intrinsic + ΔG_strain\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7. Summary\n", "\n", "This tutorial covered:\n", "\n", "1. **Ring strain** - cycloalkanes, homodesmotic reactions\n", "2. **Conformational strain** - axial vs equatorial in cyclohexanes\n", "3. **Steric strain** - gauche interactions\n", "4. **Angle strain** - deviation from ideal geometries\n", "5. **Binding strain** - energy penalty for bioactive conformers\n", "\n", "Key insights:\n", "- Strain affects **stability**, **reactivity**, and **binding**\n", "- Neural network potentials capture strain accurately\n", "- Drug design should minimize unnecessary strain" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Cleanup\n", "for f in [ring_file, deriv_file, strain_file]:\n", " if f in dir() and os.path.exists(f):\n", " os.unlink(f)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.10.0" } }, "nbformat": 4, "nbformat_minor": 4 }