Reaction Thermodynamics with Auto3D

This notebook demonstrates calculating thermodynamic properties of chemical reactions using Auto3D’s neural network potentials. We cover:

1. Reaction enthalpy (ΔH) - heat released or absorbed

2. Reaction free energy (ΔG) - spontaneity prediction

3. Entropy contributions (ΔS) - disorder changes

4. Temperature dependence - ΔG at different temperatures

Chemistry Background

Gibbs Free Energy

The Gibbs free energy change determines reaction spontaneity:

\[\Delta G = \Delta H - T\Delta S\]

• ΔG < 0: Spontaneous (exergonic)

• ΔG > 0: Non-spontaneous (endergonic)

• ΔG = 0: At equilibrium

Equilibrium Constant

\[K_{eq} = e^{-\Delta G / RT}\]

At 298K:

• ΔG = -1.36 kcal/mol → K = 10

• ΔG = -2.72 kcal/mol → K = 100

• ΔG = -5.44 kcal/mol → K = 10,000

Computational Approach

\[\Delta G_{rxn} = \sum G_{products} - \sum G_{reactants}\]

Auto3D provides:

• Electronic energy (E) from neural network potentials

• Zero-point energy (ZPE) from Hessian calculations

• Thermal corrections using ideal gas thermodynamics

import os
import tempfile
from pathlib import Path

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

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

# Constants
HARTREE_TO_KCAL = 627.509
R = 1.987e-3  # kcal/(mol·K)
T_STD = 298.15  # K

1. Simple Reaction: Ester Hydrolysis

Let’s calculate the thermodynamics of ester hydrolysis:

\[\text{Ethyl acetate} + \text{Water} \rightarrow \text{Acetic acid} + \text{Ethanol}\]

This is a classic organic chemistry reaction.

# Define reactants and products
ester_hydrolysis = {
    "reactants": {
        "ethyl_acetate": "CCOC(C)=O",
        "water": "O",
    },
    "products": {
        "acetic_acid": "CC(=O)O",
        "ethanol": "CCO",
    }
}

# Combine all species
all_species = {}
for species_dict in [ester_hydrolysis["reactants"], ester_hydrolysis["products"]]:
    all_species.update(species_dict)

print("Reaction: Ethyl acetate + H2O → Acetic acid + Ethanol")
print(f"Species to calculate: {list(all_species.keys())}")

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

print(f"Wrote {len(all_species)} species to {species_file}")

# Step 1: Generate optimized 3D structures
if __name__ == "__main__":
    config = Auto3DOptions(
        path=species_file,
        k=1,                        # Lowest energy conformer
        optimizing_engine="ANI2x",  # Use ANI2x for thermodynamics
        use_gpu=True,
        opt_steps=5000,             # Ensure full convergence
        convergence_threshold=0.003, # Tight convergence for thermo
    )

    opt_output = main(config)
    print(f"Optimized structures: {opt_output}")

# Step 2: Calculate thermodynamic properties
if 'opt_output' in dir():
    thermo_output = calc_thermo(
        path=opt_output,
        model_name="ANI2x",
        gpu_idx=0,
        opt_tol=0.003,
        opt_steps=2000
    )
    print(f"Thermodynamic output: {thermo_output}")

def extract_thermo_properties(sdf_path):
    """
    Extract thermodynamic properties from Auto3D output.

    Returns dict with H, S, G in kcal/mol (S in kcal/mol/K).
    """
    mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))

    thermo_data = {}
    for mol in mols:
        if mol is None:
            continue

        name = mol.GetProp("_Name").split("_")[0]  # Remove suffix

        # Get properties (in Hartree from calc_thermo)
        if mol.HasProp("G_hartree"):
            G = float(mol.GetProp("G_hartree")) * HARTREE_TO_KCAL
            H = float(mol.GetProp("H_hartree")) * HARTREE_TO_KCAL
            S = float(mol.GetProp("S_hartree")) * HARTREE_TO_KCAL  # Already per K
            E = float(mol.GetProp("E_hartree")) * HARTREE_TO_KCAL
            T = float(mol.GetProp("T_K"))

            thermo_data[name] = {
                "E": E,
                "H": H,
                "S": S,
                "G": G,
                "T": T
            }

    return thermo_data


if 'thermo_output' in dir():
    thermo_data = extract_thermo_properties(thermo_output)

    print("Thermodynamic Properties (298 K):")
    print("-" * 60)
    print(f"{'Species':<20} {'E (kcal/mol)':>15} {'G (kcal/mol)':>15}")
    print("-" * 60)
    for name, data in thermo_data.items():
        print(f"{name:<20} {data['E']:>15.2f} {data['G']:>15.2f}")

def calculate_reaction_thermodynamics(thermo_data, reaction_def):
    """
    Calculate reaction thermodynamics from species data.

    Args:
        thermo_data: Dict of {species: {E, H, S, G}}
        reaction_def: Dict with 'reactants' and 'products' keys

    Returns:
        Dict with ΔE, ΔH, ΔS, ΔG
    """
    # Sum products
    E_prod = sum(thermo_data[name]["E"] for name in reaction_def["products"])
    H_prod = sum(thermo_data[name]["H"] for name in reaction_def["products"])
    S_prod = sum(thermo_data[name]["S"] for name in reaction_def["products"])
    G_prod = sum(thermo_data[name]["G"] for name in reaction_def["products"])

    # Sum reactants
    E_react = sum(thermo_data[name]["E"] for name in reaction_def["reactants"])
    H_react = sum(thermo_data[name]["H"] for name in reaction_def["reactants"])
    S_react = sum(thermo_data[name]["S"] for name in reaction_def["reactants"])
    G_react = sum(thermo_data[name]["G"] for name in reaction_def["reactants"])

    # Calculate deltas
    delta_E = E_prod - E_react
    delta_H = H_prod - H_react
    delta_S = S_prod - S_react
    delta_G = G_prod - G_react

    # Calculate equilibrium constant
    K_eq = np.exp(-delta_G / (R * T_STD))

    return {
        "ΔE": delta_E,
        "ΔH": delta_H,
        "ΔS": delta_S,
        "ΔG": delta_G,
        "K_eq": K_eq,
        "T": T_STD
    }


if 'thermo_data' in dir():
    rxn_thermo = calculate_reaction_thermodynamics(thermo_data, ester_hydrolysis)

    print("\nEster Hydrolysis Thermodynamics (298 K):")
    print("=" * 50)
    print(f"ΔE  = {rxn_thermo['ΔE']:>8.2f} kcal/mol (electronic)")
    print(f"ΔH  = {rxn_thermo['ΔH']:>8.2f} kcal/mol (enthalpy)")
    print(f"ΔS  = {rxn_thermo['ΔS']*1000:>8.2f} cal/(mol·K) (entropy)")
    print(f"ΔG  = {rxn_thermo['ΔG']:>8.2f} kcal/mol (free energy)")
    print(f"K_eq = {rxn_thermo['K_eq']:>8.2e}")
    print("=" * 50)

    if rxn_thermo['ΔG'] < 0:
        print("Reaction is EXERGONIC (spontaneous)")
    else:
        print("Reaction is ENDERGONIC (non-spontaneous)")

2. Isomerization Reactions

Isomerization reactions (same molecular formula, different structure) are straightforward to analyze since there’s no change in number of molecules.

# Keto-enol tautomerization of acetone
keto_enol = {
    "reactants": {"acetone_keto": "CC(=O)C"},
    "products": {"acetone_enol": "CC(O)=C"}
}

# Cyclohexane conformations (chair flip)
# Note: Both are actually the same molecule, just different conformations
# This is more about conformational thermodynamics

# cis-trans isomerization of 2-butene
butene_isomerization = {
    "reactants": {"cis_2_butene": "C/C=C\\C"},
    "products": {"trans_2_butene": "C/C=C/C"}
}

print("Isomerization reactions to analyze:")
print("1. Keto-enol tautomerization of acetone")
print("2. cis-trans isomerization of 2-butene")

# Combine isomers
isomers = {}
isomers.update(keto_enol["reactants"])
isomers.update(keto_enol["products"])
isomers.update(butene_isomerization["reactants"])
isomers.update(butene_isomerization["products"])

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

print(f"Isomers to calculate: {list(isomers.keys())}")

# Generate structures and calculate thermodynamics
if __name__ == "__main__":
    config = Auto3DOptions(
        path=isomer_file,
        k=1,
        optimizing_engine="ANI2x",
        use_gpu=True,
    )

    isomer_opt = main(config)
    isomer_thermo_file = calc_thermo(isomer_opt, "ANI2x")
    print(f"Thermodynamic output: {isomer_thermo_file}")

if 'isomer_thermo_file' in dir():
    isomer_thermo = extract_thermo_properties(isomer_thermo_file)

    # Keto-enol equilibrium
    if all(k in isomer_thermo for k in ["acetone_keto", "acetone_enol"]):
        keto_enol_thermo = calculate_reaction_thermodynamics(isomer_thermo, keto_enol)

        print("Keto-Enol Tautomerization (Acetone):")
        print(f"  ΔG = {keto_enol_thermo['ΔG']:.2f} kcal/mol")
        print(f"  K_eq = {keto_enol_thermo['K_eq']:.2e}")
        print(f"  Keto:Enol ratio = {1/keto_enol_thermo['K_eq']:.0f}:1"
              if keto_enol_thermo['K_eq'] < 1
              else f"  Keto:Enol ratio = 1:{keto_enol_thermo['K_eq']:.0f}")
        print()

    # cis-trans isomerization
    if all(k in isomer_thermo for k in ["cis_2_butene", "trans_2_butene"]):
        butene_thermo = calculate_reaction_thermodynamics(isomer_thermo, butene_isomerization)

        print("cis-trans Isomerization (2-Butene):")
        print(f"  ΔG = {butene_thermo['ΔG']:.2f} kcal/mol")
        print(f"  trans is {'more' if butene_thermo['ΔG'] < 0 else 'less'} stable")

3. Combustion Reactions

Combustion reactions are useful for validating computational methods against experimental heats of formation.

\[\text{CH}_4 + 2\text{O}_2 \rightarrow \text{CO}_2 + 2\text{H}_2\text{O}\]

Experimental ΔH = -212.8 kcal/mol

# Note: For combustion, we need to handle stoichiometry
# and triplet O2 (which NNPs may not handle well)

# Let's use a simpler example - hydrogenation
# Ethene + H2 → Ethane

hydrogenation = {
    "reactants": {
        "ethene": "C=C",
        "hydrogen": "[H][H]",
    },
    "products": {
        "ethane": "CC",
    }
}

# Experimental ΔH ≈ -32.7 kcal/mol

print("Hydrogenation of Ethene:")
print("CH2=CH2 + H2 → CH3-CH3")
print("Experimental ΔH ≈ -32.7 kcal/mol")

4. Drug Metabolism: Glucuronidation

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.

# Simplified glucuronidation model
# Phenol + glucuronic acid → phenyl glucuronide + H2O

# This is a model reaction - the real reaction involves UDP-glucuronic acid
# and is enzyme-catalyzed

phenol_metabolism = {
    "drug": "phenol",
    "smiles": "OC1=CC=CC=C1",
    "metabolite": "phenyl_glucuronide",
    "metabolite_smiles": "OC1=CC=CC=C1.OC1C(O)C(O)C(O)C(C(=O)O)O1"  # Simplified
}

print("Drug metabolism example:")
print(f"Drug: {phenol_metabolism['drug']}")
print(f"Metabolite: {phenol_metabolism['metabolite']}")

5. Temperature Dependence of ΔG

The Gibbs free energy changes with temperature:

\[\Delta G(T) = \Delta H - T\Delta S\]

This affects:

• Equilibrium constants

• Reaction spontaneity

• Metabolic reactions at body temperature (310 K)

def calculate_dG_at_temperature(delta_H, delta_S, T):
    """
    Calculate ΔG at different temperatures.

    Assumes ΔH and ΔS are constant (valid for small T ranges).
    """
    return delta_H - T * delta_S


def plot_dG_vs_T(delta_H, delta_S, T_range=(200, 500)):
    """
    Calculate ΔG over a temperature range.

    Returns DataFrame with T, ΔG, K_eq.
    """
    temperatures = np.linspace(T_range[0], T_range[1], 50)

    data = []
    for T in temperatures:
        dG = calculate_dG_at_temperature(delta_H, delta_S, T)
        K_eq = np.exp(-dG / (R * T))
        data.append({"T_K": T, "ΔG": dG, "K_eq": K_eq})

    return pd.DataFrame(data)


# Example: reaction where ΔH and ΔS have opposite signs
# (entropy-driven reaction)
if 'rxn_thermo' in dir():
    df_temp = plot_dG_vs_T(rxn_thermo['ΔH'], rxn_thermo['ΔS'])

    print("Temperature Dependence of ΔG:")
    print(df_temp.iloc[::10].to_string(index=False))

    # Find temperature where ΔG = 0
    if rxn_thermo['ΔS'] != 0:
        T_eq = rxn_thermo['ΔH'] / rxn_thermo['ΔS']
        if 0 < T_eq < 1000:
            print(f"\nEquilibrium temperature (ΔG=0): {T_eq:.1f} K")

6. Reaction Energy Diagrams

Create energy diagrams showing reactants, products, and their relative energies.

def create_reaction_diagram_data(thermo_data, reaction_def):
    """
    Create data for a reaction energy diagram.
    """
    # Get reactant energies
    reactant_names = list(reaction_def["reactants"].keys())
    reactant_G = sum(thermo_data[name]["G"] for name in reactant_names)

    # Get product energies
    product_names = list(reaction_def["products"].keys())
    product_G = sum(thermo_data[name]["G"] for name in product_names)

    # Set reactants as reference (0)
    return {
        "states": ["Reactants", "Products"],
        "G_rel": [0, product_G - reactant_G],
        "labels": [
            " + ".join(reactant_names),
            " + ".join(product_names)
        ]
    }


if 'thermo_data' in dir():
    diagram = create_reaction_diagram_data(thermo_data, ester_hydrolysis)

    print("Reaction Energy Diagram:")
    print("=" * 50)
    for state, G, label in zip(diagram["states"], diagram["G_rel"], diagram["labels"]):
        bar = "█" * max(1, int(10 + G/2))  # Simple bar representation
        print(f"{state:12s} {bar} {G:>8.2f} kcal/mol")
        print(f"             ({label})")
    print("=" * 50)

7. Best Practices for Reaction Thermodynamics

Accuracy Considerations

1. Use consistent methods: Same model and basis for all species

2. Tight optimization: Use convergence_threshold=0.003 or tighter

3. Verify minima: Check for imaginary frequencies

4. Consider conformers: Use lowest-energy conformer for each species

Limitations of NNPs

• Radicals and triplets: NNPs trained on closed-shell singlets

• Transition metals: Not supported by ANI/AIMNet

• High-energy species: Training data bias toward stable molecules

Validation

Compare with:

• Experimental data (NIST, CRC Handbook)

• Higher-level QM calculations (CCSD(T), DFT)

# Cleanup
for f in [species_file, isomer_file]:
    if f in dir() and os.path.exists(f):
        os.unlink(f)

Summary

This tutorial demonstrated:

1. Calculating ΔH, ΔS, ΔG for chemical reactions

2. Equilibrium constants from free energy differences

3. Isomerization thermodynamics - keto-enol, cis-trans

4. Temperature dependence of reaction spontaneity

5. Reaction energy diagrams for visualization

Key workflow:

1. Define reactants and products

2. Optimize 3D structures with Auto3D

3. Calculate thermodynamic properties with calc_thermo()

4. Compute ΔG = ΣG(products) - ΣG(reactants)