{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stereochemistry in Drug Discovery\n",
    "\n",
    "Stereochemistry is critical in drug development - enantiomers can have dramatically different:\n",
    "\n",
    "- **Pharmacological activity** (thalidomide: one enantiomer teratogenic)\n",
    "- **Potency** (often 10-1000x difference)\n",
    "- **Metabolism** (different CYP450 interactions)\n",
    "- **Toxicity** (stereoselective off-target effects)\n",
    "\n",
    "This notebook covers:\n",
    "\n",
    "1. **Stereoisomer enumeration** - generating all possible stereoisomers\n",
    "2. **Chiral center analysis** - identifying and annotating stereocenters\n",
    "3. **3D structure generation** - maintaining stereochemistry\n",
    "4. **Enantiomer comparison** - energy and geometry differences\n",
    "\n",
    "## Chemistry Background\n",
    "\n",
    "### Types of Stereoisomers\n",
    "\n",
    "| Type | Relationship | Example |\n",
    "|------|--------------|--------|\n",
    "| **Enantiomers** | Mirror images, non-superimposable | R/S alanine |\n",
    "| **Diastereomers** | Multiple stereocenters, not mirror images | threonine vs allo-threonine |\n",
    "| **Epimers** | Differ at one stereocenter | α/β-glucose |\n",
    "| **Cis/Trans** | Geometric isomers around double bonds | cis/trans-stilbene |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import tempfile\n",
    "from pathlib import Path\n",
    "\n",
    "import Auto3D\n",
    "from Auto3D import Auto3DOptions, main\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import AllChem, Descriptors, rdMolDescriptors\n",
    "from rdkit.Chem import Draw\n",
    "from rdkit.Chem.EnumerateStereoisomers import EnumerateStereoisomers, StereoEnumerationOptions\n",
    "import pandas as pd\n",
    "\n",
    "print(f\"Auto3D version: {Auto3D.__version__}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Famous Stereochemistry Examples in Drugs\n",
    "\n",
    "These examples illustrate why stereochemistry matters in drug development."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Classic examples where stereochemistry determines activity\n",
    "stereo_examples = {\n",
    "    # Thalidomide - one of the most notorious examples\n",
    "    # (R)-thalidomide: sedative\n",
    "    # (S)-thalidomide: teratogenic\n",
    "    \"thalidomide_R\": \"O=C1CC[C@H](N2C(=O)C3=CC=CC=C3C2=O)C(=O)N1\",\n",
    "    \"thalidomide_S\": \"O=C1CC[C@@H](N2C(=O)C3=CC=CC=C3C2=O)C(=O)N1\",\n",
    "    \n",
    "    # Ibuprofen - only S-enantiomer is active\n",
    "    \"ibuprofen_S\": \"CC(C)CC1=CC=C(C=C1)[C@H](C)C(=O)O\",\n",
    "    \"ibuprofen_R\": \"CC(C)CC1=CC=C(C=C1)[C@@H](C)C(=O)O\",\n",
    "    \n",
    "    # Omeprazole vs Esomeprazole (S-omeprazole)\n",
    "    \"esomeprazole_S\": \"COC1=CC2=NC([C@H](C)S(=O)C3=NC4=CC=CC=C4N3)=NC2=CC1OC\",\n",
    "    \n",
    "    # Methadone - R-enantiomer 10x more potent\n",
    "    \"methadone_R\": \"CC[C@@](C)(C(=O)CC)C(C1=CC=CC=C1)C2=CC=CC=C2\",\n",
    "    \n",
    "    # Propranolol - S-enantiomer 100x more potent at β-receptors\n",
    "    \"propranolol_S\": \"CC(C)NC[C@H](O)COC1=CC=CC2=CC=CC=C12\",\n",
    "}\n",
    "\n",
    "print(\"Stereospecific drugs:\")\n",
    "for name, smi in stereo_examples.items():\n",
    "    mol = Chem.MolFromSmiles(smi)\n",
    "    n_stereo = len(Chem.FindMolChiralCenters(mol, includeUnassigned=True))\n",
    "    print(f\"  {name}: {n_stereo} chiral center(s)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Identifying Chiral Centers\n",
    "\n",
    "Before generating stereoisomers, we need to identify all stereocenters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analyze_stereocenters(smiles):\n",
    "    \"\"\"\n",
    "    Analyze stereochemistry of a molecule.\n",
    "    \n",
    "    Returns:\n",
    "    - Number of defined stereocenters\n",
    "    - Number of undefined stereocenters\n",
    "    - Maximum possible stereoisomers (2^n)\n",
    "    \"\"\"\n",
    "    mol = Chem.MolFromSmiles(smiles)\n",
    "    \n",
    "    # Find all chiral centers\n",
    "    chiral_centers = Chem.FindMolChiralCenters(mol, includeUnassigned=True)\n",
    "    \n",
    "    defined = sum(1 for _, tag in chiral_centers if tag in ['R', 'S'])\n",
    "    undefined = sum(1 for _, tag in chiral_centers if tag == '?')\n",
    "    \n",
    "    # Check for double bond stereochemistry\n",
    "    double_bond_stereo = 0\n",
    "    for bond in mol.GetBonds():\n",
    "        if bond.GetBondType() == Chem.BondType.DOUBLE:\n",
    "            stereo = bond.GetStereo()\n",
    "            if stereo in [Chem.BondStereo.STEREOE, Chem.BondStereo.STEREOZ,\n",
    "                          Chem.BondStereo.STEREOANY]:\n",
    "                double_bond_stereo += 1\n",
    "    \n",
    "    total_undefined = undefined + double_bond_stereo\n",
    "    max_isomers = 2 ** (defined + undefined + double_bond_stereo)\n",
    "    \n",
    "    return {\n",
    "        \"defined_centers\": defined,\n",
    "        \"undefined_centers\": undefined,\n",
    "        \"double_bond_stereo\": double_bond_stereo,\n",
    "        \"max_stereoisomers\": max_isomers,\n",
    "        \"chiral_centers\": chiral_centers\n",
    "    }\n",
    "\n",
    "\n",
    "# Analyze some drug molecules\n",
    "drugs_to_analyze = {\n",
    "    \"atorvastatin\": \"CC(C)C1=C(C(=C(N1CCC(CC(CC(=O)O)O)O)C2=CC=C(C=C2)F)C3=CC=CC=C3)C(=O)NC4=CC=CC=C4\",\n",
    "    \"simvastatin\": \"CC[C@H](C)C(=C)[C@H]1CC[C@@H]2C1(CC[C@H](C2)OC(=O)C)C\",\n",
    "    \"penicillin_G\": \"CC1([C@@H](N2[C@H](S1)[C@@H](C2=O)NC(=O)CC3=CC=CC=C3)C(=O)O)C\",\n",
    "}\n",
    "\n",
    "print(\"Stereochemistry analysis:\")\n",
    "print(\"-\" * 70)\n",
    "for name, smi in drugs_to_analyze.items():\n",
    "    info = analyze_stereocenters(smi)\n",
    "    print(f\"{name}:\")\n",
    "    print(f\"  Defined centers: {info['defined_centers']}\")\n",
    "    print(f\"  Undefined centers: {info['undefined_centers']}\")\n",
    "    print(f\"  Max stereoisomers: {info['max_stereoisomers']}\")\n",
    "    print(f\"  Chiral centers: {info['chiral_centers']}\")\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Enumerating Stereoisomers\n",
    "\n",
    "For molecules with undefined stereochemistry, we need to enumerate all possible isomers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def enumerate_stereoisomers(smiles, max_isomers=16):\n",
    "    \"\"\"\n",
    "    Enumerate all stereoisomers of a molecule.\n",
    "    \n",
    "    Args:\n",
    "        smiles: Input SMILES (can have undefined stereochemistry)\n",
    "        max_isomers: Maximum number of isomers to generate\n",
    "        \n",
    "    Returns:\n",
    "        List of (SMILES, R/S assignment) tuples\n",
    "    \"\"\"\n",
    "    mol = Chem.MolFromSmiles(smiles)\n",
    "    \n",
    "    # Configure enumeration options\n",
    "    opts = StereoEnumerationOptions(\n",
    "        tryEmbedding=True,  # Verify 3D embedding is possible\n",
    "        unique=True,        # Remove duplicates\n",
    "        maxIsomers=max_isomers\n",
    "    )\n",
    "    \n",
    "    isomers = list(EnumerateStereoisomers(mol, options=opts))\n",
    "    \n",
    "    results = []\n",
    "    for iso in isomers:\n",
    "        smi = Chem.MolToSmiles(iso, isomericSmiles=True)\n",
    "        centers = Chem.FindMolChiralCenters(iso)\n",
    "        config = \"/\".join([f\"{idx}:{tag}\" for idx, tag in centers])\n",
    "        results.append((smi, config))\n",
    "    \n",
    "    return results\n",
    "\n",
    "\n",
    "# Enumerate stereoisomers for a molecule with undefined centers\n",
    "# Ephedrine - 2 chiral centers = 4 possible stereoisomers\n",
    "ephedrine_undefined = \"CC(C(C1=CC=CC=C1)O)NC\"  # No stereochemistry defined\n",
    "\n",
    "isomers = enumerate_stereoisomers(ephedrine_undefined)\n",
    "print(f\"Ephedrine stereoisomers ({len(isomers)} found):\")\n",
    "print(\"-\" * 60)\n",
    "for i, (smi, config) in enumerate(isomers):\n",
    "    print(f\"Isomer {i+1}: {config}\")\n",
    "    print(f\"  SMILES: {smi}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. 3D Structure Generation Preserving Stereochemistry\n",
    "\n",
    "Auto3D correctly generates 3D structures that preserve the input stereochemistry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate 3D for all stereoisomers\n",
    "stereoisomer_names = {\n",
    "    \"1R,2S_ephedrine\": \"C[C@@H]([C@H](C1=CC=CC=C1)O)NC\",  # (-)-ephedrine\n",
    "    \"1S,2R_ephedrine\": \"C[C@H]([C@@H](C1=CC=CC=C1)O)NC\",  # (+)-ephedrine\n",
    "    \"1R,2R_pseudoephedrine\": \"C[C@@H]([C@@H](C1=CC=CC=C1)O)NC\",  # (+)-pseudoephedrine\n",
    "    \"1S,2S_pseudoephedrine\": \"C[C@H]([C@H](C1=CC=CC=C1)O)NC\",  # (-)-pseudoephedrine\n",
    "}\n",
    "\n",
    "# Note: ephedrine and pseudoephedrine are diastereomers!\n",
    "# They have different physical properties (melting point, solubility)\n",
    "\n",
    "with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:\n",
    "    for name, smi in stereoisomer_names.items():\n",
    "        f.write(f\"{smi} {name}\\n\")\n",
    "    stereo_file = f.name\n",
    "\n",
    "print(f\"Wrote {len(stereoisomer_names)} stereoisomers to {stereo_file}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run Auto3D - it will preserve stereochemistry\n",
    "if __name__ == \"__main__\":\n",
    "    config = Auto3DOptions(\n",
    "        path=stereo_file,\n",
    "        k=1,\n",
    "        optimizing_engine=\"AIMNET\",\n",
    "        use_gpu=True,\n",
    "        # enumerate_isomer=False,  # Don't enumerate - we already have all isomers\n",
    "    )\n",
    "    \n",
    "    stereo_output = main(config)\n",
    "    print(f\"Output: {stereo_output}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Comparing Enantiomer Energies\n",
    "\n",
    "In vacuum, enantiomers have identical energies. However:\n",
    "- Diastereomers have **different** energies\n",
    "- In chiral environments (proteins, enzymes), enantiomers interact differently"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_stereoisomer_energies(sdf_path):\n",
    "    \"\"\"\n",
    "    Compare energies of stereoisomers.\n",
    "    \n",
    "    Enantiomers should have identical energies in vacuum.\n",
    "    Diastereomers will have different energies.\n",
    "    \"\"\"\n",
    "    mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))\n",
    "    \n",
    "    data = []\n",
    "    for mol in mols:\n",
    "        if mol is None:\n",
    "            continue\n",
    "            \n",
    "        name = mol.GetProp(\"_Name\")\n",
    "        \n",
    "        # Get energy\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\n",
    "        else:\n",
    "            continue\n",
    "        \n",
    "        # Get stereochemistry\n",
    "        centers = Chem.FindMolChiralCenters(mol)\n",
    "        stereo = \"/\".join([tag for _, tag in centers])\n",
    "        \n",
    "        data.append({\n",
    "            \"Name\": name,\n",
    "            \"Stereo\": stereo,\n",
    "            \"E_hartree\": e\n",
    "        })\n",
    "    \n",
    "    df = pd.DataFrame(data)\n",
    "    df[\"E_kcal\"] = (df[\"E_hartree\"] - df[\"E_hartree\"].min()) * 627.509\n",
    "    \n",
    "    return df.sort_values(\"E_kcal\")\n",
    "\n",
    "\n",
    "if 'stereo_output' in dir():\n",
    "    df_stereo = compare_stereoisomer_energies(stereo_output)\n",
    "    print(\"Stereoisomer Energy Comparison:\")\n",
    "    print(df_stereo[[\"Name\", \"Stereo\", \"E_kcal\"]].to_string(index=False))\n",
    "    print(\"\\nNote: Enantiomer pairs should have ~identical energies.\")\n",
    "    print(\"Diastereomers (ephedrine vs pseudoephedrine) have different energies.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Double Bond Stereochemistry (E/Z)\n",
    "\n",
    "Geometric isomers around double bonds are also important in drug design."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Examples of E/Z isomerism in drugs\n",
    "ez_examples = {\n",
    "    # Tamoxifen - E and Z have opposite effects!\n",
    "    \"tamoxifen_Z\": \"CC/C(=C(\\\\C1=CC=CC=C1)/C2=CC=C(OCCN(C)C)C=C2)/C3=CC=CC=C3\",\n",
    "    \"tamoxifen_E\": \"CC/C(=C(/C1=CC=CC=C1)\\\\C2=CC=C(OCCN(C)C)C=C2)/C3=CC=CC=C3\",\n",
    "    \n",
    "    # Resveratrol - trans form is biologically active\n",
    "    \"resveratrol_trans\": \"OC1=CC=C(/C=C/C2=CC(O)=CC(O)=C2)C=C1\",\n",
    "    \"resveratrol_cis\": \"OC1=CC=C(/C=C\\\\C2=CC(O)=CC(O)=C2)C=C1\",\n",
    "    \n",
    "    # Retinoids - critical for vitamin A function\n",
    "    \"retinoic_acid_all_trans\": \"CC1=C(C(CCC1)(C)C)/C=C/C(C)=C/C=C/C(C)=C/C(=O)O\",\n",
    "}\n",
    "\n",
    "with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:\n",
    "    for name, smi in ez_examples.items():\n",
    "        f.write(f\"{smi} {name}\\n\")\n",
    "    ez_file = f.name\n",
    "\n",
    "print(f\"E/Z isomers: {len(ez_examples)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate 3D structures for E/Z isomers\n",
    "if __name__ == \"__main__\":\n",
    "    config = Auto3DOptions(\n",
    "        path=ez_file,\n",
    "        k=1,\n",
    "        optimizing_engine=\"AIMNET\",\n",
    "        use_gpu=True,\n",
    "    )\n",
    "    \n",
    "    ez_output = main(config)\n",
    "    print(f\"Output: {ez_output}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare E/Z isomer energies\n",
    "if 'ez_output' in dir():\n",
    "    df_ez = compare_stereoisomer_energies(ez_output)\n",
    "    print(\"E/Z Isomer Energy Comparison:\")\n",
    "    print(df_ez[[\"Name\", \"E_kcal\"]].to_string(index=False))\n",
    "    print(\"\\nNote: E and Z isomers have different energies (unlike enantiomers).\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Auto3D Stereoisomer Enumeration\n",
    "\n",
    "Auto3D can automatically enumerate stereoisomers during conformer generation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start with undefined stereochemistry\n",
    "undefined_stereo = {\n",
    "    \"drug_candidate\": \"CC(NC(=O)C(CC1=CC=CC=C1)NC(=O)OC(C)(C)C)C(=O)O\",\n",
    "}\n",
    "\n",
    "with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:\n",
    "    for name, smi in undefined_stereo.items():\n",
    "        f.write(f\"{smi} {name}\\n\")\n",
    "    undefined_file = f.name\n",
    "\n",
    "# Count stereocenters\n",
    "info = analyze_stereocenters(undefined_stereo[\"drug_candidate\"])\n",
    "print(f\"Drug candidate has {info['defined_centers'] + info['undefined_centers']} stereocenters\")\n",
    "print(f\"Maximum possible stereoisomers: {info['max_stereoisomers']}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run with stereoisomer enumeration\n",
    "if __name__ == \"__main__\":\n",
    "    config = Auto3DOptions(\n",
    "        path=undefined_file,\n",
    "        k=1,                        # Top conformer per stereoisomer\n",
    "        enumerate_isomer=True,      # Enumerate stereoisomers\n",
    "        optimizing_engine=\"AIMNET\",\n",
    "        use_gpu=True,\n",
    "    )\n",
    "    \n",
    "    enum_output = main(config)\n",
    "    print(f\"Output: {enum_output}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Analyze enumerated stereoisomers\n",
    "if 'enum_output' in dir():\n",
    "    mols = list(Chem.SDMolSupplier(enum_output, removeHs=False))\n",
    "    print(f\"Generated {len(mols)} stereoisomers\")\n",
    "    \n",
    "    df_enum = compare_stereoisomer_energies(enum_output)\n",
    "    print(\"\\nStereoisomer Ranking by Energy:\")\n",
    "    print(df_enum[[\"Name\", \"Stereo\", \"E_kcal\"]].to_string(index=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Best Practices for Stereochemistry in Drug Discovery\n",
    "\n",
    "### Lead Optimization\n",
    "1. **Early determination**: Identify the active stereoisomer early\n",
    "2. **Synthesize pure isomers**: Test individual stereoisomers, not racemates\n",
    "3. **Consider racemization**: Some drugs racemize in vivo (ibuprofen)\n",
    "\n",
    "### Virtual Screening\n",
    "1. **Define stereochemistry**: Use explicit @/@@/E/Z notation\n",
    "2. **Enumerate if unknown**: Use Auto3D's `enumerate_isomer=True`\n",
    "3. **Dock all isomers**: Different stereoisomers may have different binding modes\n",
    "\n",
    "### QSAR/ML\n",
    "1. **Consistent representation**: Use canonical SMILES with stereochemistry\n",
    "2. **3D descriptors**: Capture stereo-dependent properties\n",
    "3. **Chiral descriptors**: Consider adding chirality-specific features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cleanup\n",
    "for f in [stereo_file, ez_file, undefined_file]:\n",
    "    if f in dir() and os.path.exists(f):\n",
    "        os.unlink(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "This tutorial covered:\n",
    "\n",
    "1. **Stereochemistry fundamentals** - enantiomers, diastereomers, E/Z isomers\n",
    "2. **Chiral center analysis** - identifying stereocenters in molecules\n",
    "3. **Stereoisomer enumeration** - generating all possible isomers\n",
    "4. **3D generation** - Auto3D preserves input stereochemistry\n",
    "5. **Energy comparison** - enantiomers identical, diastereomers different\n",
    "\n",
    "Key takeaway: **Always specify stereochemistry** in drug discovery workflows. Undefined stereochemistry leads to ambiguous results."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}