Tutorial 2: Thermochemical Property Estimation (Heat Capacity for Gases)

# Cell 1: Import Required Packages

!pip install ASE
!pip install mace-torch
!pip install rdkit

# Cell 2: Import Required Libraries

# Tools for plotting charts
import numpy as np                                # Computational library
import matplotlib.pyplot as plt                   # Plot graphs
import pandas as pd                               # Excel of python

# RDKit is for SMILES to construct a molecule
from rdkit import Chem                            # Used to build molecules (basic package)
from rdkit.Chem import AllChem, Descriptors       # Used to build molecules (advanced package)
from rdkit.Chem import rdMolDescriptors           # Used to obtain information of molecules

# MACE-OFF is the machine learning potential trained on organic molecules
from mace.calculators import mace_off             # MACE-OFF (Machine Learning Potential)

# ASE (Atomic Simulation Environment) computes information recieved from the machine learning potential
from ase import Atoms, units                      # Represents a molecule object with information and conversion for units
from ase.optimize import LBFGS                    # Optimization / energy minimization
from ase.vibrations import Vibrations             # Used to calculate vibrational modes of the Atom object
from ase.thermochemistry import IdealGasThermo    # Allows you to calculate entropy, enthalpy, and gibbs free energy

# Cell 3: Define Molecule using Simplified Molecular Input Line Entry System (SMILES) to construct it

SMILES = "O=C=O"          # Carbon Dioxide


# Other examples to try:
# SMILES = "CC(=O)O"    # Acetic acid
# SMILES = "c1ccccc1"   # Benzene
# SMILES = "CCCC"       # Butane

TEMPERATURE_K  = 298.15   # Temperature in [K] (temperature focused on in the chart below)
MODEL_SIZE     = "medium"  # 'small' (fast) | 'medium' | 'large' (A greater size generally correlates to higher accuracy with the cost of computational time)
DEVICE         = "cpu"    # 'cpu' or 'cuda' (Determines what the code will be running on, a cpu or cuda core (gpu))
FMAX           = 0.01     # eV/Å optimisation threshold use between 0.01 and 0.05 depending on how long the code runs for (higher threshold runs faster)

# Cell 4: Defines Fucntion that takes a SMILES and returns if it is valid.  Then it prints some information about the molecule

SUPPORTED_ELEMENTS = {"C", "H", "N", "O", "F", "Cl", "Br", "S", "P"} # Elements that MACE-OFF is trained on (it should get accurate results with molecules made of these elements)

def validate_smiles(smiles):
    """Validate SMILES and check element compatibility with MACE-OFF."""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        raise ValueError(f"Invalid SMILES string: '{smiles}'")

    mol = Chem.AddHs(mol)  # add explicit hydrogens
    elements = {atom.GetSymbol() for atom in mol.GetAtoms()}
    unsupported = elements - SUPPORTED_ELEMENTS

    if unsupported:
        raise ValueError(
            f"Unsupported elements for MACE-OFF: {unsupported}\n"
            f"MACE-OFF supports: {SUPPORTED_ELEMENTS}"
        )

    formula = rdMolDescriptors.CalcMolFormula(mol)
    mw = Descriptors.MolWt(mol)
    print(f"  SMILES   : {smiles}")
    print(f"  Formula  : {formula}")
    print(f"  Mol. wt  : {mw:.2f} g/mol")
    print(f"  Elements : {elements}")
    print(f"  MACE-OFF compatible: ✓")
    return mol, formula, mw

print("Validating molecule...")
rdkit_mol, formula, mol_weight = validate_smiles(SMILES)

# Cell 5: Defines Function to Create Atom layout to be optimized later by MACE-OFF

def smiles_to_ase(smiles):
    """Convert SMILES - RDKit 3D - ASE Atoms object."""
    mol = Chem.MolFromSmiles(smiles)
    mol = Chem.AddHs(mol)

    # Embed 3D coordinates
    params = AllChem.ETKDGv3()
    params.randomSeed = 42
    result = AllChem.EmbedMolecule(mol, params)
    if result != 0:
        raise RuntimeError("RDKit failed to generate 3D coordinates.")

    # Pre-optimise with MMFF94 (gives MACE-OFF a better starting point)
    AllChem.MMFFOptimizeMolecule(mol, maxIters=2000)

    # Extract positions and elements
    conf = mol.GetConformer()
    positions = conf.GetPositions()   # Angstrom
    symbols   = [atom.GetSymbol() for atom in mol.GetAtoms()]

    return Atoms(symbols=symbols, positions=positions)

atoms = smiles_to_ase(SMILES)
print(f"3D structure generated: {len(atoms)} atoms")
print(f"Symbols: {atoms.get_chemical_symbols()}")

# Cell 6: Loads MACE-OFF and moves atoms (slightly) around until the forces on them fall under the input threshold below (this is done to avoid huge forces initially ruining the simulation)

print(f"Loading MACE-OFF ({MODEL_SIZE}) on {DEVICE}...")
calc = mace_off(model=MODEL_SIZE, device=DEVICE)
print("Calculator ready.\n")

# Attach calculator and optimise
atoms.calc = calc
print("Optimising geometry with MACE-OFF...")
opt = LBFGS(atoms, logfile=None)
opt.run(fmax=FMAX)

E0 = atoms.get_potential_energy()   # eV
print(f"Optimisation complete.")
print(f"  Energy        : {E0:.6f} eV = {E0 * 96.485:.6f} kJ")
print(f"  Max force     : {np.max(np.linalg.norm(atoms.get_forces(), axis=1)):.4f} eV/Å")

# Cell 7: Finds the potential energy of the molecule including vibrational energy (also filters out imaginary frequencies)

from pathlib import Path
import shutil

# Clear stale cache
shutil.rmtree("vib_cache", ignore_errors=True)

vib_dir = Path("vib_cache")
vib_dir.mkdir(exist_ok=True)

print(f"Running vibrational analysis ({3 * len(atoms)} force evaluations)...")
vib = Vibrations(atoms, name=str(vib_dir / "vib"), delta=0.01)
vib.run()

# Get all energies
vib_energies = vib.get_energies()

# Filter out imaginary and near-zero frequencies
real_vib_energies = vib_energies[vib_energies.real > 0.01]

# Check for too many imaginary modes
n_imaginary = np.sum(vib_energies.real < 0)
if n_imaginary > 2:
    print(f"⚠️  Warning: {n_imaginary} imaginary frequencies detected.")
    print("    Consider tightening FMAX for more reliable results.")
else:
    print(f"✓  {n_imaginary} imaginary mode(s) filtered — looks clean.")

print(f"\nTotal modes      : {len(vib_energies)}")
print(f"Imaginary/zero   : {len(vib_energies) - len(real_vib_energies)}")
print(f"Real modes kept  : {len(real_vib_energies)}")

# Print frequencies manually in cm⁻¹
print("\nVibrational frequencies (real modes):")
freqs_cm = real_vib_energies.real * 8065.54   # eV → cm⁻¹
for i, f in enumerate(freqs_cm):
    print(f"  Mode {i+1:>3d}: {f:>10.2f} cm⁻¹")

# Cell 8: Determines characteristics of the molecule to accurately compute chemical properties



#**** It is correct for this tutorial, but if you change molecules / atoms, make sure to change the geometry, symmetry, and spin below in this cell****



thermo = IdealGasThermo(
    vib_energies=real_vib_energies,
    potentialenergy=E0,
    atoms=atoms,
    geometry='linear',  # Change to match molecule (Linear or nonlinear)
    symmetrynumber=2,   # Change to match molecule (how many times you can turn the molecule and have symmetry)
    spin=0,             # Change to match molecule (0.5 times number of unpaired electrons)
)

# Uses the Constant Pressure Enthalpy Equation

# Cp via finite-difference dH/dT
def get_Cp_JmolK(T, dT=1.0):
    H_p = thermo.get_enthalpy(T + dT, verbose=False)
    H_m = thermo.get_enthalpy(T - dT, verbose=False)
    return (H_p - H_m) / (2 * dT) * 96485.3   # eV/K to J/(mol·K)

Cp_ref = get_Cp_JmolK(TEMPERATURE_K)

print(f"\n{'='*45}")
print(f"  GAS HEAT CAPACITY — {formula}")
print(f"  T = {TEMPERATURE_K} K")
print(f"{'='*45}")
print(f"  Cₚ = {Cp_ref:>8.2f}  J/(mol·K)")
print(f"{'='*45}")

# Cell 9: Cp vs T curve

# Rigid-Rotor Harmonic-Oscillator (Tracks atomic roational energy, not stretching at all, and treats bonds like springs)

# temps is the x-axis range of the graph and Cp_vals calls the function previously used to compute the y value or Cp for all x values, using a step up and step down to compute
temps = np.linspace(200, 1000, 80)
Cp_vals = [get_Cp_JmolK(T) for T in temps]

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(temps, Cp_vals, color="steelblue", lw=2.5)
ax.axvline(TEMPERATURE_K, color="gray", ls="--", alpha=0.7,
           label=f"T = {TEMPERATURE_K} K  →  Cₚ = {Cp_ref:.1f} J/(mol·K)")
ax.scatter([TEMPERATURE_K], [Cp_ref], color="tomato", zorder=5, s=60)
ax.set_xlabel("Temperature (K)", fontsize=12)
ax.set_ylabel("Cₚ  [J / (mol·K)]", fontsize=12)
ax.set_title(f"Gas-phase Cₚ(T) — {formula}  (MACE-OFF {MODEL_SIZE}, RRHO)", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.25)
plt.tight_layout()
plt.show()

# Cell 10: Error between MACE-OFF and NIST (CO2)

# Compute Cp at multiple temperatures
temps_limits_CO2 = [298, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200]
Cp_vals_limits_CO2 = [get_Cp_JmolK(T) for T in temps_limits_CO2]
NIST_Cp_vals_CO2 = [37.12, 37.22, 41.34, 44.61, 47.32, 49.57, 51.44, 53.00, 54.30, 55.40, 56.35]


Cp_vals_limits_CO2 = np.array(Cp_vals_limits_CO2)

data_CO2 = {
    "Temperature in [K]": temps_limits_CO2,
    "Cp [J/mol*K]": Cp_vals_limits_CO2,
    "Percent Error in [%]": np.abs((Cp_vals_limits_CO2 - NIST_Cp_vals_CO2)/NIST_Cp_vals_CO2) * 100,
}

pd.set_option("display.width", 1000)
df = pd.DataFrame(data_CO2)

display(df)

# Comparing Extremely High Temperature Values

# Compute Cp at multiple temperatures
temps_limits_ext_CO2 = [3000, 4000, 5000, 6000]
Cp_vals_limits_ext_CO2 = [get_Cp_JmolK(T) for T in temps_limits_ext_CO2]
NIST_Cp_vals_ext_CO2 = [62.23, 63.25, 64.06, 64.98]


Cp_vals_limits_ext_CO2 = np.array(Cp_vals_limits_ext_CO2)

data_1_CO2 = {
    "Temperature in [K]": temps_limits_ext_CO2,
    "Cp [J/mol*K]": Cp_vals_limits_ext_CO2,
    "Percent Error in [%]": np.abs((Cp_vals_limits_ext_CO2 - NIST_Cp_vals_ext_CO2)/NIST_Cp_vals_ext_CO2) * 100,
}

pd.set_option("display.width", 1000)
df_1 = pd.DataFrame(data_1_CO2)

print("")
display(df_1)

# Cell 11: This will compute the same information for a different gas (Carbon Monoxide) - Remember to change SMILES, symmetry number, and spin number when changing gases

# You can also change reference Temp (TEMPERATURE_K) to see Cp at different temperatures

SMILES = "[C-]#[O+]"  # Had to play with the SMILES here, C#O did not work (keep this in mind if trying another gas)

TEMPERATURE_K  = 298.15

print("Validating molecule...")
rdkit_mol, formula, mol_weight = validate_smiles(SMILES)

atoms = smiles_to_ase(SMILES)
print(f"3D structure generated: {len(atoms)} atoms")
print(f"Symbols: {atoms.get_chemical_symbols()}")

atoms.calc = calc
print("Optimising geometry with MACE-OFF...")
opt = LBFGS(atoms, logfile=None)
opt.run(fmax=FMAX)

E0 = atoms.get_potential_energy()   # eV
print(f"Optimisation complete.")
print(f"  Energy        : {E0:.6f} eV")
print(f"  Max force     : {np.max(np.linalg.norm(atoms.get_forces(), axis=1)):.4f} eV/Å")



shutil.rmtree("vib_cache", ignore_errors=True)

vib_dir = Path("vib_cache")
vib_dir.mkdir(exist_ok=True)

print(f"Running vibrational analysis ({3 * len(atoms)} force evaluations)...")
vib = Vibrations(atoms, name=str(vib_dir / "vib"), delta=0.01)
vib.run()

# Get all energies
vib_energies = vib.get_energies()

# Filter out imaginary and near-zero frequencies
real_vib_energies = vib_energies[vib_energies.real > 0.01]

# Check for too many imaginary modes
n_imaginary = np.sum(vib_energies.real < 0)
if n_imaginary > 2:
    print(f"⚠️  Warning: {n_imaginary} imaginary frequencies detected.")
    print("    Consider tightening FMAX for more reliable results.")
else:
    print(f"✓  {n_imaginary} imaginary mode(s) filtered — looks clean.")

print(f"\nTotal modes      : {len(vib_energies)}")
print(f"Imaginary/zero   : {len(vib_energies) - len(real_vib_energies)}")
print(f"Real modes kept  : {len(real_vib_energies)}")

# Print frequencies manually in cm⁻¹
print("\nVibrational frequencies (real modes):")
freqs_cm = real_vib_energies.real * 8065.54   # eV → cm⁻¹
for i, f in enumerate(freqs_cm):
    print(f"  Mode {i+1:>3d}: {f:>10.2f} cm⁻¹")


thermo = IdealGasThermo(
    vib_energies=real_vib_energies,
    potentialenergy=E0,
    atoms=atoms,
    geometry='linear',  # Change to match molecule (Linear or nonlinear)
    symmetrynumber=1,   # Change to match molecule (how many times you can turn the molecule and have symmetry)
    spin=1,             # Change to match molecule (0.5 times number of unpaired electrons)
)

Cp_ref = get_Cp_JmolK(TEMPERATURE_K)

print(f"\n{'='*45}")
print(f"  GAS HEAT CAPACITY — {formula}")
print(f"  T = {TEMPERATURE_K} K")
print(f"{'='*45}")
print(f"  Cₚ = {Cp_ref:>8.2f}  J/(mol·K)")
print(f"{'='*45}")

temps = np.linspace(200, 1000, 80)
Cp_vals = [get_Cp_JmolK(T) for T in temps]

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(temps, Cp_vals, color="steelblue", lw=2.5)
ax.axvline(TEMPERATURE_K, color="gray", ls="--", alpha=0.7,
           label=f"T = {TEMPERATURE_K} K  →  Cₚ = {Cp_ref:.1f} J/(mol·K)")
ax.scatter([TEMPERATURE_K], [Cp_ref], color="tomato", zorder=5, s=60)
ax.set_xlabel("Temperature (K)", fontsize=12)
ax.set_ylabel("Cₚ  [J / (mol·K)]", fontsize=12)
ax.set_title(f"Gas-phase Cₚ(T) — {formula}  (MACE-OFF {MODEL_SIZE}, RRHO)", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.25)
plt.tight_layout()
plt.show()

# Cell 12: Error between MACE-OFF and NIST (CO)

# Compute Cp at multiple temperatures
temps_limits_CO = [298, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300]
Cp_vals_limits_CO = [get_Cp_JmolK(T) for T in temps_limits_CO]
NIST_Cp_vals_CO = [29.15, 29.15, 29.30, 29.82, 30.47, 31.17, 31.88, 32.55, 33.18, 33.73, 34.20, 34.55]


Cp_vals_limits_CO = np.array(Cp_vals_limits_CO)

data_CO = {
    "Temperature in [K]": temps_limits_CO,
    "Cp [J/mol*K]": Cp_vals_limits_CO,
    "Percent Error in [%]": np.abs((Cp_vals_limits_CO - NIST_Cp_vals_CO)/NIST_Cp_vals_CO) * 100,
}

pd.set_option("display.width", 1000)
df = pd.DataFrame(data_CO)

display(df)

# Comparing Extremely High Temperature Values

# Compute Cp at multiple temperatures
temps_limits_ext_CO = [3000, 4000, 5000, 6000]
Cp_vals_limits_ext_CO = [get_Cp_JmolK(T) for T in temps_limits_ext_CO]
NIST_Cp_vals_ext_CO = [37.20, 37.72, 38.07, 38.37]


Cp_vals_limits_ext_CO = np.array(Cp_vals_limits_ext_CO)

data_1_CO = {
    "Temperature in [K]": temps_limits_ext_CO,
    "Cp [J/mol*K]": Cp_vals_limits_ext_CO,
    "Percent Error in [%]": np.abs((Cp_vals_limits_ext_CO - NIST_Cp_vals_ext_CO)/NIST_Cp_vals_ext_CO) * 100,
}

pd.set_option("display.width", 1000)
df_1 = pd.DataFrame(data_1_CO)

print("")
display(df_1)

# Cell 13: This will compute the same information for a different gas (Hydrogen Peroxide) - Remember to change SMILES, symmetry number, and spin number when changing gases

# You can also change reference Temp (TEMPERATURE_K) to see Cp at different temperatures

SMILES = "[H]OO[H]"

TEMPERATURE_K  = 298.15

print("Validating molecule...")
rdkit_mol, formula, mol_weight = validate_smiles(SMILES)

atoms = smiles_to_ase(SMILES)
print(f"3D structure generated: {len(atoms)} atoms")
print(f"Symbols: {atoms.get_chemical_symbols()}")

atoms.calc = calc
print("Optimising geometry with MACE-OFF...")
opt = LBFGS(atoms, logfile=None)
opt.run(fmax=FMAX)

E0 = atoms.get_potential_energy()   # eV
print(f"Optimisation complete.")
print(f"  Energy        : {E0:.6f} eV")
print(f"  Max force     : {np.max(np.linalg.norm(atoms.get_forces(), axis=1)):.4f} eV/Å")



shutil.rmtree("vib_cache", ignore_errors=True)

vib_dir = Path("vib_cache")
vib_dir.mkdir(exist_ok=True)

print(f"Running vibrational analysis ({3 * len(atoms)} force evaluations)...")
vib = Vibrations(atoms, name=str(vib_dir / "vib"), delta=0.01)
vib.run()

# Get all energies
vib_energies = vib.get_energies()

# Filter out imaginary and near-zero frequencies
real_vib_energies = vib_energies[vib_energies.real > 0.01]

# Check for too many imaginary modes
n_imaginary = np.sum(vib_energies.real < 0)
if n_imaginary > 2:
    print(f"⚠️  Warning: {n_imaginary} imaginary frequencies detected.")
    print("    Consider tightening FMAX for more reliable results.")
else:
    print(f"✓  {n_imaginary} imaginary mode(s) filtered — looks clean.")

print(f"\nTotal modes      : {len(vib_energies)}")
print(f"Imaginary/zero   : {len(vib_energies) - len(real_vib_energies)}")
print(f"Real modes kept  : {len(real_vib_energies)}")

# Print frequencies manually in cm⁻¹
print("\nVibrational frequencies (real modes):")
freqs_cm = real_vib_energies.real * 8065.54   # eV → cm⁻¹
for i, f in enumerate(freqs_cm):
    print(f"  Mode {i+1:>3d}: {f:>10.2f} cm⁻¹")


thermo = IdealGasThermo(
    vib_energies=real_vib_energies,
    potentialenergy=E0,
    atoms=atoms,
    geometry='nonlinear', # Change to match molecule (Linear or nonlinear)
    symmetrynumber=2,     # Change to match molecule (how many times you can turn the molecule and have symmetry)
    spin=0,               # Change to match molecule (0.5 times number of unpaired electrons)
)

Cp_ref = get_Cp_JmolK(TEMPERATURE_K)

print(f"\n{'='*45}")
print(f"  GAS HEAT CAPACITY — {formula}")
print(f"  T = {TEMPERATURE_K} K")
print(f"{'='*45}")
print(f"  Cₚ = {Cp_ref:>8.2f}  J/(mol·K)")
print(f"{'='*45}")

temps = np.linspace(200, 1000, 80)
Cp_vals = [get_Cp_JmolK(T) for T in temps]

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(temps, Cp_vals, color="steelblue", lw=2.5)
ax.axvline(TEMPERATURE_K, color="gray", ls="--", alpha=0.7,
           label=f"T = {TEMPERATURE_K} K  →  Cₚ = {Cp_ref:.1f} J/(mol·K)")
ax.scatter([TEMPERATURE_K], [Cp_ref], color="tomato", zorder=5, s=60)
ax.set_xlabel("Temperature (K)", fontsize=12)
ax.set_ylabel("Cₚ  [J / (mol·K)]", fontsize=12)
ax.set_title(f"Gas-phase Cₚ(T) — {formula}  (MACE-OFF {MODEL_SIZE}, RRHO)", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.25)
plt.tight_layout()
plt.show()

# Cell 14: Error between MACE-OFF and NIST (H2O2)

# Compute Cp at multiple temperatures
temps_limits_H2O2 = [298, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500]
Cp_vals_limits_H2O2 = [get_Cp_JmolK(T) for T in temps_limits_H2O2]
NIST_Cp_vals_H2O2 = [43.07, 43.20, 48.65, 52.51, 55.50, 57.92, 59.90, 61.55, 62.95, 64.17, 65.27, 66.30, 67.33, 68.42]

Cp_vals_limits_H2O2 = np.array(Cp_vals_limits_H2O2)

data_H2O2 = {
    "Temperature in [K]": temps_limits_H2O2,
    "Cp [J/mol*K]": Cp_vals_limits_H2O2,
    "Percent Error in [%]": np.abs((Cp_vals_limits_H2O2 - NIST_Cp_vals_H2O2)/NIST_Cp_vals_H2O2) * 100,
}

pd.set_option("display.width", 1000)
df = pd.DataFrame(data_H2O2)

display(df)

# Cell 15: This will compute the same information for a different gas (Methane) - Remember to change SMILES, symmetry number, and spin number when changing gases

# You can also change reference Temp (TEMPERATURE_K) to see Cp at different temperatures

SMILES = "C"

TEMPERATURE_K  = 298.15

print("Validating molecule...")
rdkit_mol, formula, mol_weight = validate_smiles(SMILES)

atoms = smiles_to_ase(SMILES)
print(f"3D structure generated: {len(atoms)} atoms")
print(f"Symbols: {atoms.get_chemical_symbols()}")

atoms.calc = calc
print("Optimising geometry with MACE-OFF...")
opt = LBFGS(atoms, logfile=None)
opt.run(fmax=FMAX)

E0 = atoms.get_potential_energy()   # eV
print(f"Optimisation complete.")
print(f"  Energy        : {E0:.6f} eV")
print(f"  Max force     : {np.max(np.linalg.norm(atoms.get_forces(), axis=1)):.4f} eV/Å")

shutil.rmtree("vib_cache", ignore_errors=True)

vib_dir = Path("vib_cache")
vib_dir.mkdir(exist_ok=True)

print(f"Running vibrational analysis ({3 * len(atoms)} force evaluations)...")
vib = Vibrations(atoms, name=str(vib_dir / "vib"), delta=0.01)
vib.run()

# Get all energies
vib_energies = vib.get_energies()

# Filter out imaginary and near-zero frequencies
real_vib_energies = vib_energies[vib_energies.real > 0.01]

# Check for too many imaginary modes
n_imaginary = np.sum(vib_energies.real < 0)
if n_imaginary > 2:
    print(f"⚠️  Warning: {n_imaginary} imaginary frequencies detected.")
    print("    Consider tightening FMAX for more reliable results.")
else:
    print(f"✓  {n_imaginary} imaginary mode(s) filtered — looks clean.")

print(f"\nTotal modes      : {len(vib_energies)}")
print(f"Imaginary/zero   : {len(vib_energies) - len(real_vib_energies)}")
print(f"Real modes kept  : {len(real_vib_energies)}")

# Print frequencies manually in cm⁻¹
print("\nVibrational frequencies (real modes):")
freqs_cm = real_vib_energies.real * 8065.54   # eV → cm⁻¹
for i, f in enumerate(freqs_cm):
    print(f"  Mode {i+1:>3d}: {f:>10.2f} cm⁻¹")


thermo = IdealGasThermo(
    vib_energies=real_vib_energies,
    potentialenergy=E0,
    atoms=atoms,
    geometry='nonlinear', # Change to match molecule (Linear or nonlinear)
    symmetrynumber=12,     # Change to match molecule (how many times you can turn the molecule and have symmetry)
    spin=0,               # Change to match molecule (0.5 times number of unpaired electrons)
)

Cp_ref = get_Cp_JmolK(TEMPERATURE_K)

print(f"\n{'='*45}")
print(f"  GAS HEAT CAPACITY — {formula}")
print(f"  T = {TEMPERATURE_K} K")
print(f"{'='*45}")
print(f"  Cₚ = {Cp_ref:>8.2f}  J/(mol·K)")
print(f"{'='*45}")

temps = np.linspace(200, 1000, 80)
Cp_vals = [get_Cp_JmolK(T) for T in temps]

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(temps, Cp_vals, color="steelblue", lw=2.5)
ax.axvline(TEMPERATURE_K, color="gray", ls="--", alpha=0.7,
           label=f"T = {TEMPERATURE_K} K  →  Cₚ = {Cp_ref:.1f} J/(mol·K)")
ax.scatter([TEMPERATURE_K], [Cp_ref], color="tomato", zorder=5, s=60)
ax.set_xlabel("Temperature (K)", fontsize=12)
ax.set_ylabel("Cₚ  [J / (mol·K)]", fontsize=12)
ax.set_title(f"Gas-phase Cₚ(T) — {formula}  (MACE-OFF {MODEL_SIZE}, RRHO)", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.25)
plt.tight_layout()
plt.show()

# Cell 16: Error between MACE-OFF and NIST (CH4)

# Compute Cp at multiple temperatures
temps_limits_CH4 = [298, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300]
Cp_vals_limits_CH4 = [get_Cp_JmolK(T) for T in temps_limits_CH4]
NIST_Cp_vals_CH4 = [35.69, 35.76, 40.63, 46.63, 52.74, 58.60, 64.08, 69.14, 73.75, 77.92, 81.68, 85.07]


Cp_vals_limits_CH4 = np.array(Cp_vals_limits_CH4)

data_CH4 = {
    "Temperature in [K]": temps_limits_CH4,
    "Cp [J/mol*K]": Cp_vals_limits_CH4,
    "Percent Error in [%]": np.abs((Cp_vals_limits_CH4 - NIST_Cp_vals_CH4)/NIST_Cp_vals_CH4) * 100,
}

pd.set_option("display.width", 1000)
df = pd.DataFrame(data_CH4)

display(df)

# Comparing Extremely High Temperature Values

# Compute Cp at multiple temperatures
temps_limits_ext_CH4 = [2000, 2500, 3000]
Cp_vals_limits_ext_CH4 = [get_Cp_JmolK(T) for T in temps_limits_ext_CH4]
NIST_Cp_vals_ext_CH4 = [101.24, 108.23, 113.55]

Cp_vals_limits_ext_CH4 = np.array(Cp_vals_limits_ext_CH4)

data_1_CH4 = {
    "Temperature in [K]": temps_limits_ext_CH4,
    "Cp [J/mol*K]": Cp_vals_limits_ext_CH4,
    "Percent Error in [%]": np.abs((Cp_vals_limits_ext_CH4 - NIST_Cp_vals_ext_CH4)/NIST_Cp_vals_ext_CH4) * 100,
}

pd.set_option("display.width", 1000)
df_1 = pd.DataFrame(data_1_CO)

print("")
display(df_1)