{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reaction Thermodynamics with Auto3D\n",
    "\n",
    "This notebook demonstrates calculating thermodynamic properties of chemical reactions using Auto3D's neural network potentials. We cover:\n",
    "\n",
    "1. **Reaction enthalpy (ΔH)** - heat released or absorbed\n",
    "2. **Reaction free energy (ΔG)** - spontaneity prediction\n",
    "3. **Entropy contributions (ΔS)** - disorder changes\n",
    "4. **Temperature dependence** - ΔG at different temperatures\n",
    "\n",
    "## Chemistry Background\n",
    "\n",
    "### Gibbs Free Energy\n",
    "\n",
    "The Gibbs free energy change determines reaction spontaneity:\n",
    "\n",
    "$$\\Delta G = \\Delta H - T\\Delta S$$\n",
    "\n",
    "- **ΔG < 0**: Spontaneous (exergonic)\n",
    "- **ΔG > 0**: Non-spontaneous (endergonic)\n",
    "- **ΔG = 0**: At equilibrium\n",
    "\n",
    "### Equilibrium Constant\n",
    "\n",
    "$$K_{eq} = e^{-\\Delta G / RT}$$\n",
    "\n",
    "At 298K:\n",
    "- ΔG = -1.36 kcal/mol → K = 10\n",
    "- ΔG = -2.72 kcal/mol → K = 100\n",
    "- ΔG = -5.44 kcal/mol → K = 10,000\n",
    "\n",
    "### Computational Approach\n",
    "\n",
    "$$\\Delta G_{rxn} = \\sum G_{products} - \\sum G_{reactants}$$\n",
    "\n",
    "Auto3D provides:\n",
    "- Electronic energy (E) from neural network potentials\n",
    "- Zero-point energy (ZPE) from Hessian calculations\n",
    "- Thermal corrections using ideal gas thermodynamics"
   ]
  },
  {
   "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 Auto3D.ASE.thermo import calc_thermo\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import AllChem, Descriptors\n",
    "import pandas as pd\n",
    "\n",
    "print(f\"Auto3D version: {Auto3D.__version__}\")\n",
    "\n",
    "# Constants\n",
    "HARTREE_TO_KCAL = 627.509\n",
    "R = 1.987e-3  # kcal/(mol·K)\n",
    "T_STD = 298.15  # K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Simple Reaction: Ester Hydrolysis\n",
    "\n",
    "Let's calculate the thermodynamics of ester hydrolysis:\n",
    "\n",
    "$$\\text{Ethyl acetate} + \\text{Water} \\rightarrow \\text{Acetic acid} + \\text{Ethanol}$$\n",
    "\n",
    "This is a classic organic chemistry reaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define reactants and products\n",
    "ester_hydrolysis = {\n",
    "    \"reactants\": {\n",
    "        \"ethyl_acetate\": \"CCOC(C)=O\",\n",
    "        \"water\": \"O\",\n",
    "    },\n",
    "    \"products\": {\n",
    "        \"acetic_acid\": \"CC(=O)O\",\n",
    "        \"ethanol\": \"CCO\",\n",
    "    }\n",
    "}\n",
    "\n",
    "# Combine all species\n",
    "all_species = {}\n",
    "for species_dict in [ester_hydrolysis[\"reactants\"], ester_hydrolysis[\"products\"]]:\n",
    "    all_species.update(species_dict)\n",
    "\n",
    "print(\"Reaction: Ethyl acetate + H2O → Acetic acid + Ethanol\")\n",
    "print(f\"Species to calculate: {list(all_species.keys())}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write all species to file\n",
    "with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:\n",
    "    for name, smi in all_species.items():\n",
    "        f.write(f\"{smi} {name}\\n\")\n",
    "    species_file = f.name\n",
    "\n",
    "print(f\"Wrote {len(all_species)} species to {species_file}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Generate optimized 3D structures\n",
    "if __name__ == \"__main__\":\n",
    "    config = Auto3DOptions(\n",
    "        path=species_file,\n",
    "        k=1,                        # Lowest energy conformer\n",
    "        optimizing_engine=\"ANI2x\",  # Use ANI2x for thermodynamics\n",
    "        use_gpu=True,\n",
    "        opt_steps=5000,             # Ensure full convergence\n",
    "        convergence_threshold=0.003, # Tight convergence for thermo\n",
    "    )\n",
    "    \n",
    "    opt_output = main(config)\n",
    "    print(f\"Optimized structures: {opt_output}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2: Calculate thermodynamic properties\n",
    "if 'opt_output' in dir():\n",
    "    thermo_output = calc_thermo(\n",
    "        path=opt_output,\n",
    "        model_name=\"ANI2x\",\n",
    "        gpu_idx=0,\n",
    "        opt_tol=0.003,\n",
    "        opt_steps=2000\n",
    "    )\n",
    "    print(f\"Thermodynamic output: {thermo_output}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_thermo_properties(sdf_path):\n",
    "    \"\"\"\n",
    "    Extract thermodynamic properties from Auto3D output.\n",
    "    \n",
    "    Returns dict with H, S, G in kcal/mol (S in kcal/mol/K).\n",
    "    \"\"\"\n",
    "    mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))\n",
    "    \n",
    "    thermo_data = {}\n",
    "    for mol in mols:\n",
    "        if mol is None:\n",
    "            continue\n",
    "            \n",
    "        name = mol.GetProp(\"_Name\").split(\"_\")[0]  # Remove suffix\n",
    "        \n",
    "        # Get properties (in Hartree from calc_thermo)\n",
    "        if mol.HasProp(\"G_hartree\"):\n",
    "            G = float(mol.GetProp(\"G_hartree\")) * HARTREE_TO_KCAL\n",
    "            H = float(mol.GetProp(\"H_hartree\")) * HARTREE_TO_KCAL\n",
    "            S = float(mol.GetProp(\"S_hartree\")) * HARTREE_TO_KCAL  # Already per K\n",
    "            E = float(mol.GetProp(\"E_hartree\")) * HARTREE_TO_KCAL\n",
    "            T = float(mol.GetProp(\"T_K\"))\n",
    "            \n",
    "            thermo_data[name] = {\n",
    "                \"E\": E,\n",
    "                \"H\": H,\n",
    "                \"S\": S,\n",
    "                \"G\": G,\n",
    "                \"T\": T\n",
    "            }\n",
    "    \n",
    "    return thermo_data\n",
    "\n",
    "\n",
    "if 'thermo_output' in dir():\n",
    "    thermo_data = extract_thermo_properties(thermo_output)\n",
    "    \n",
    "    print(\"Thermodynamic Properties (298 K):\")\n",
    "    print(\"-\" * 60)\n",
    "    print(f\"{'Species':<20} {'E (kcal/mol)':>15} {'G (kcal/mol)':>15}\")\n",
    "    print(\"-\" * 60)\n",
    "    for name, data in thermo_data.items():\n",
    "        print(f\"{name:<20} {data['E']:>15.2f} {data['G']:>15.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_reaction_thermodynamics(thermo_data, reaction_def):\n",
    "    \"\"\"\n",
    "    Calculate reaction thermodynamics from species data.\n",
    "    \n",
    "    Args:\n",
    "        thermo_data: Dict of {species: {E, H, S, G}}\n",
    "        reaction_def: Dict with 'reactants' and 'products' keys\n",
    "        \n",
    "    Returns:\n",
    "        Dict with ΔE, ΔH, ΔS, ΔG\n",
    "    \"\"\"\n",
    "    # Sum products\n",
    "    E_prod = sum(thermo_data[name][\"E\"] for name in reaction_def[\"products\"])\n",
    "    H_prod = sum(thermo_data[name][\"H\"] for name in reaction_def[\"products\"])\n",
    "    S_prod = sum(thermo_data[name][\"S\"] for name in reaction_def[\"products\"])\n",
    "    G_prod = sum(thermo_data[name][\"G\"] for name in reaction_def[\"products\"])\n",
    "    \n",
    "    # Sum reactants\n",
    "    E_react = sum(thermo_data[name][\"E\"] for name in reaction_def[\"reactants\"])\n",
    "    H_react = sum(thermo_data[name][\"H\"] for name in reaction_def[\"reactants\"])\n",
    "    S_react = sum(thermo_data[name][\"S\"] for name in reaction_def[\"reactants\"])\n",
    "    G_react = sum(thermo_data[name][\"G\"] for name in reaction_def[\"reactants\"])\n",
    "    \n",
    "    # Calculate deltas\n",
    "    delta_E = E_prod - E_react\n",
    "    delta_H = H_prod - H_react\n",
    "    delta_S = S_prod - S_react\n",
    "    delta_G = G_prod - G_react\n",
    "    \n",
    "    # Calculate equilibrium constant\n",
    "    K_eq = np.exp(-delta_G / (R * T_STD))\n",
    "    \n",
    "    return {\n",
    "        \"ΔE\": delta_E,\n",
    "        \"ΔH\": delta_H,\n",
    "        \"ΔS\": delta_S,\n",
    "        \"ΔG\": delta_G,\n",
    "        \"K_eq\": K_eq,\n",
    "        \"T\": T_STD\n",
    "    }\n",
    "\n",
    "\n",
    "if 'thermo_data' in dir():\n",
    "    rxn_thermo = calculate_reaction_thermodynamics(thermo_data, ester_hydrolysis)\n",
    "    \n",
    "    print(\"\\nEster Hydrolysis Thermodynamics (298 K):\")\n",
    "    print(\"=\" * 50)\n",
    "    print(f\"ΔE  = {rxn_thermo['ΔE']:>8.2f} kcal/mol (electronic)\")\n",
    "    print(f\"ΔH  = {rxn_thermo['ΔH']:>8.2f} kcal/mol (enthalpy)\")\n",
    "    print(f\"ΔS  = {rxn_thermo['ΔS']*1000:>8.2f} cal/(mol·K) (entropy)\")\n",
    "    print(f\"ΔG  = {rxn_thermo['ΔG']:>8.2f} kcal/mol (free energy)\")\n",
    "    print(f\"K_eq = {rxn_thermo['K_eq']:>8.2e}\")\n",
    "    print(\"=\" * 50)\n",
    "    \n",
    "    if rxn_thermo['ΔG'] < 0:\n",
    "        print(\"Reaction is EXERGONIC (spontaneous)\")\n",
    "    else:\n",
    "        print(\"Reaction is ENDERGONIC (non-spontaneous)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Isomerization Reactions\n",
    "\n",
    "Isomerization reactions (same molecular formula, different structure) are straightforward to analyze since there's no change in number of molecules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Keto-enol tautomerization of acetone\n",
    "keto_enol = {\n",
    "    \"reactants\": {\"acetone_keto\": \"CC(=O)C\"},\n",
    "    \"products\": {\"acetone_enol\": \"CC(O)=C\"}\n",
    "}\n",
    "\n",
    "# Cyclohexane conformations (chair flip)\n",
    "# Note: Both are actually the same molecule, just different conformations\n",
    "# This is more about conformational thermodynamics\n",
    "\n",
    "# cis-trans isomerization of 2-butene\n",
    "butene_isomerization = {\n",
    "    \"reactants\": {\"cis_2_butene\": \"C/C=C\\\\C\"},\n",
    "    \"products\": {\"trans_2_butene\": \"C/C=C/C\"}\n",
    "}\n",
    "\n",
    "print(\"Isomerization reactions to analyze:\")\n",
    "print(\"1. Keto-enol tautomerization of acetone\")\n",
    "print(\"2. cis-trans isomerization of 2-butene\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine isomers\n",
    "isomers = {}\n",
    "isomers.update(keto_enol[\"reactants\"])\n",
    "isomers.update(keto_enol[\"products\"])\n",
    "isomers.update(butene_isomerization[\"reactants\"])\n",
    "isomers.update(butene_isomerization[\"products\"])\n",
    "\n",
    "with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:\n",
    "    for name, smi in isomers.items():\n",
    "        f.write(f\"{smi} {name}\\n\")\n",
    "    isomer_file = f.name\n",
    "\n",
    "print(f\"Isomers to calculate: {list(isomers.keys())}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate structures and calculate thermodynamics\n",
    "if __name__ == \"__main__\":\n",
    "    config = Auto3DOptions(\n",
    "        path=isomer_file,\n",
    "        k=1,\n",
    "        optimizing_engine=\"ANI2x\",\n",
    "        use_gpu=True,\n",
    "    )\n",
    "    \n",
    "    isomer_opt = main(config)\n",
    "    isomer_thermo_file = calc_thermo(isomer_opt, \"ANI2x\")\n",
    "    print(f\"Thermodynamic output: {isomer_thermo_file}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'isomer_thermo_file' in dir():\n",
    "    isomer_thermo = extract_thermo_properties(isomer_thermo_file)\n",
    "    \n",
    "    # Keto-enol equilibrium\n",
    "    if all(k in isomer_thermo for k in [\"acetone_keto\", \"acetone_enol\"]):\n",
    "        keto_enol_thermo = calculate_reaction_thermodynamics(isomer_thermo, keto_enol)\n",
    "        \n",
    "        print(\"Keto-Enol Tautomerization (Acetone):\")\n",
    "        print(f\"  ΔG = {keto_enol_thermo['ΔG']:.2f} kcal/mol\")\n",
    "        print(f\"  K_eq = {keto_enol_thermo['K_eq']:.2e}\")\n",
    "        print(f\"  Keto:Enol ratio = {1/keto_enol_thermo['K_eq']:.0f}:1\" \n",
    "              if keto_enol_thermo['K_eq'] < 1 \n",
    "              else f\"  Keto:Enol ratio = 1:{keto_enol_thermo['K_eq']:.0f}\")\n",
    "        print()\n",
    "    \n",
    "    # cis-trans isomerization\n",
    "    if all(k in isomer_thermo for k in [\"cis_2_butene\", \"trans_2_butene\"]):\n",
    "        butene_thermo = calculate_reaction_thermodynamics(isomer_thermo, butene_isomerization)\n",
    "        \n",
    "        print(\"cis-trans Isomerization (2-Butene):\")\n",
    "        print(f\"  ΔG = {butene_thermo['ΔG']:.2f} kcal/mol\")\n",
    "        print(f\"  trans is {'more' if butene_thermo['ΔG'] < 0 else 'less'} stable\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Combustion Reactions\n",
    "\n",
    "Combustion reactions are useful for validating computational methods against experimental heats of formation.\n",
    "\n",
    "$$\\text{CH}_4 + 2\\text{O}_2 \\rightarrow \\text{CO}_2 + 2\\text{H}_2\\text{O}$$\n",
    "\n",
    "Experimental ΔH = -212.8 kcal/mol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: For combustion, we need to handle stoichiometry\n",
    "# and triplet O2 (which NNPs may not handle well)\n",
    "\n",
    "# Let's use a simpler example - hydrogenation\n",
    "# Ethene + H2 → Ethane\n",
    "\n",
    "hydrogenation = {\n",
    "    \"reactants\": {\n",
    "        \"ethene\": \"C=C\",\n",
    "        \"hydrogen\": \"[H][H]\",\n",
    "    },\n",
    "    \"products\": {\n",
    "        \"ethane\": \"CC\",\n",
    "    }\n",
    "}\n",
    "\n",
    "# Experimental ΔH ≈ -32.7 kcal/mol\n",
    "\n",
    "print(\"Hydrogenation of Ethene:\")\n",
    "print(\"CH2=CH2 + H2 → CH3-CH3\")\n",
    "print(\"Experimental ΔH ≈ -32.7 kcal/mol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Drug Metabolism: Glucuronidation\n",
    "\n",
    "In drug metabolism, Phase II reactions like glucuronidation are important. While we can't model the full enzymatic process, we can estimate the thermodynamics of the chemical transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simplified glucuronidation model\n",
    "# Phenol + glucuronic acid → phenyl glucuronide + H2O\n",
    "\n",
    "# This is a model reaction - the real reaction involves UDP-glucuronic acid\n",
    "# and is enzyme-catalyzed\n",
    "\n",
    "phenol_metabolism = {\n",
    "    \"drug\": \"phenol\",\n",
    "    \"smiles\": \"OC1=CC=CC=C1\",\n",
    "    \"metabolite\": \"phenyl_glucuronide\",\n",
    "    \"metabolite_smiles\": \"OC1=CC=CC=C1.OC1C(O)C(O)C(O)C(C(=O)O)O1\"  # Simplified\n",
    "}\n",
    "\n",
    "print(\"Drug metabolism example:\")\n",
    "print(f\"Drug: {phenol_metabolism['drug']}\")\n",
    "print(f\"Metabolite: {phenol_metabolism['metabolite']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Temperature Dependence of ΔG\n",
    "\n",
    "The Gibbs free energy changes with temperature:\n",
    "\n",
    "$$\\Delta G(T) = \\Delta H - T\\Delta S$$\n",
    "\n",
    "This affects:\n",
    "- Equilibrium constants\n",
    "- Reaction spontaneity\n",
    "- Metabolic reactions at body temperature (310 K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_dG_at_temperature(delta_H, delta_S, T):\n",
    "    \"\"\"\n",
    "    Calculate ΔG at different temperatures.\n",
    "    \n",
    "    Assumes ΔH and ΔS are constant (valid for small T ranges).\n",
    "    \"\"\"\n",
    "    return delta_H - T * delta_S\n",
    "\n",
    "\n",
    "def plot_dG_vs_T(delta_H, delta_S, T_range=(200, 500)):\n",
    "    \"\"\"\n",
    "    Calculate ΔG over a temperature range.\n",
    "    \n",
    "    Returns DataFrame with T, ΔG, K_eq.\n",
    "    \"\"\"\n",
    "    temperatures = np.linspace(T_range[0], T_range[1], 50)\n",
    "    \n",
    "    data = []\n",
    "    for T in temperatures:\n",
    "        dG = calculate_dG_at_temperature(delta_H, delta_S, T)\n",
    "        K_eq = np.exp(-dG / (R * T))\n",
    "        data.append({\"T_K\": T, \"ΔG\": dG, \"K_eq\": K_eq})\n",
    "    \n",
    "    return pd.DataFrame(data)\n",
    "\n",
    "\n",
    "# Example: reaction where ΔH and ΔS have opposite signs\n",
    "# (entropy-driven reaction)\n",
    "if 'rxn_thermo' in dir():\n",
    "    df_temp = plot_dG_vs_T(rxn_thermo['ΔH'], rxn_thermo['ΔS'])\n",
    "    \n",
    "    print(\"Temperature Dependence of ΔG:\")\n",
    "    print(df_temp.iloc[::10].to_string(index=False))\n",
    "    \n",
    "    # Find temperature where ΔG = 0\n",
    "    if rxn_thermo['ΔS'] != 0:\n",
    "        T_eq = rxn_thermo['ΔH'] / rxn_thermo['ΔS']\n",
    "        if 0 < T_eq < 1000:\n",
    "            print(f\"\\nEquilibrium temperature (ΔG=0): {T_eq:.1f} K\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Reaction Energy Diagrams\n",
    "\n",
    "Create energy diagrams showing reactants, products, and their relative energies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_reaction_diagram_data(thermo_data, reaction_def):\n",
    "    \"\"\"\n",
    "    Create data for a reaction energy diagram.\n",
    "    \"\"\"\n",
    "    # Get reactant energies\n",
    "    reactant_names = list(reaction_def[\"reactants\"].keys())\n",
    "    reactant_G = sum(thermo_data[name][\"G\"] for name in reactant_names)\n",
    "    \n",
    "    # Get product energies  \n",
    "    product_names = list(reaction_def[\"products\"].keys())\n",
    "    product_G = sum(thermo_data[name][\"G\"] for name in product_names)\n",
    "    \n",
    "    # Set reactants as reference (0)\n",
    "    return {\n",
    "        \"states\": [\"Reactants\", \"Products\"],\n",
    "        \"G_rel\": [0, product_G - reactant_G],\n",
    "        \"labels\": [\n",
    "            \" + \".join(reactant_names),\n",
    "            \" + \".join(product_names)\n",
    "        ]\n",
    "    }\n",
    "\n",
    "\n",
    "if 'thermo_data' in dir():\n",
    "    diagram = create_reaction_diagram_data(thermo_data, ester_hydrolysis)\n",
    "    \n",
    "    print(\"Reaction Energy Diagram:\")\n",
    "    print(\"=\" * 50)\n",
    "    for state, G, label in zip(diagram[\"states\"], diagram[\"G_rel\"], diagram[\"labels\"]):\n",
    "        bar = \"█\" * max(1, int(10 + G/2))  # Simple bar representation\n",
    "        print(f\"{state:12s} {bar} {G:>8.2f} kcal/mol\")\n",
    "        print(f\"             ({label})\")\n",
    "    print(\"=\" * 50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Best Practices for Reaction Thermodynamics\n",
    "\n",
    "### Accuracy Considerations\n",
    "\n",
    "1. **Use consistent methods**: Same model and basis for all species\n",
    "2. **Tight optimization**: Use `convergence_threshold=0.003` or tighter\n",
    "3. **Verify minima**: Check for imaginary frequencies\n",
    "4. **Consider conformers**: Use lowest-energy conformer for each species\n",
    "\n",
    "### Limitations of NNPs\n",
    "\n",
    "- **Radicals and triplets**: NNPs trained on closed-shell singlets\n",
    "- **Transition metals**: Not supported by ANI/AIMNet\n",
    "- **High-energy species**: Training data bias toward stable molecules\n",
    "\n",
    "### Validation\n",
    "\n",
    "Compare with:\n",
    "- Experimental data (NIST, CRC Handbook)\n",
    "- Higher-level QM calculations (CCSD(T), DFT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cleanup\n",
    "for f in [species_file, isomer_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 demonstrated:\n",
    "\n",
    "1. **Calculating ΔH, ΔS, ΔG** for chemical reactions\n",
    "2. **Equilibrium constants** from free energy differences\n",
    "3. **Isomerization thermodynamics** - keto-enol, cis-trans\n",
    "4. **Temperature dependence** of reaction spontaneity\n",
    "5. **Reaction energy diagrams** for visualization\n",
    "\n",
    "Key workflow:\n",
    "1. Define reactants and products\n",
    "2. Optimize 3D structures with Auto3D\n",
    "3. Calculate thermodynamic properties with `calc_thermo()`\n",
    "4. Compute ΔG = ΣG(products) - ΣG(reactants)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}