Ctrl+K

Methane combustion energy

Methane combustion energy#

Consider methane combustion:

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

We calculate this reaction energy by constructing equilibrium structures of each molecule, calculating the total energy of each molecule in turn, and taking the energy difference between the reactants and products. This example is inspired by this paper, in which the authors discuss the practical aspects of calculating the methane combustion energy, taking into account the competing concerns of computational cost and accuracy. More details on the impact of electron correlation, basis set selection, and extrapolation to the complete basis set (CBS) limit can be found therein.

import numpy as np
import py3Dmol as p3d
import veloxchem as vlx

Geometry optimization#

Force field optimization#

Construct molecules using SMILES strings and perform a force field optimization.

ch4_ff = vlx.Molecule.smiles_to_xyz("C")
h2o_ff = vlx.Molecule.smiles_to_xyz("O")
co2_ff = vlx.Molecule.smiles_to_xyz("O=C=O")
o2_ff = vlx.Molecule.smiles_to_xyz("O=O")

xTB optimization#

Set up xTB driver and perform xTB optimization.

def xtb_optimization(xyz_ff):

    molecule = vlx.Molecule.read_xyz_string(xyz_ff)
    xtb_drv = vlx.XtbDriver()
    
    xtb_opt_drv = vlx.OptimizationDriver(xtb_drv)
    xtb_opt_drv.ostream.mute()
    
    xtb_opt = xtb_opt_drv.compute(molecule)
    
    return xtb_opt["final_geometry"]


ch4_xtb = xtb_optimization(ch4_ff)
h2o_xtb = xtb_optimization(h2o_ff)
co2_xtb = xtb_optimization(co2_ff)
o2_xtb = xtb_optimization(o2_ff)

Visualization#

Convert to format readable by py3Dmol and visualize structures.

viewer = p3d.view(viewergrid=(2, 2), width=400, height=300, linked=False)
viewer.addModel(ch4_xtb, "xyz", viewer=(0, 0))
viewer.addModel(o2_xtb, "xyz", viewer=(0, 1))
viewer.addModel(co2_xtb, "xyz", viewer=(1, 0))
viewer.addModel(h2o_xtb, "xyz", viewer=(1, 1))
viewer.setViewStyle({"style": "outline", "width": 0.05})
viewer.setStyle({"stick": {}, "sphere": {"scale": 0.25}})
viewer.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Energy calculation#

Calculate the total energy of each molecule using B3LYP and a 6-311G* basis set.

Note

Keep in mind that O\(_2\) has a triplet ground state, which is here set by changing the multiplicity and using unrestricted SCF.

def scf_calculation(molecule_xyz, basis_set, xc, multiplicity=None):

    molecule = vlx.Molecule.read_xyz_string(molecule_xyz)
    basis = vlx.MolecularBasis.read(molecule, basis_set, ostream=None)

    if multiplicity is not None:
        molecule.set_multiplicity(multiplicity)
        molecule.check_multiplicity()
        scf_gs = vlx.ScfUnrestrictedDriver()
        scf_gs.ostream.mute()
    else:
        scf_gs = vlx.ScfRestrictedDriver()
        scf_gs.ostream.mute()

    scf_gs.xcfun = xc
    scf_results = scf_gs.compute(molecule, basis)

    return scf_results, scf_gs


xc = "b3lyp"
basis_set = "6-311G*"

ch4_results, ch4_scf = scf_calculation(ch4_xtb, basis_set, xc)
o2_results, o2_scf = scf_calculation(o2_xtb, basis_set, xc, multiplicity=3)
co2_results, co2_scf = scf_calculation(co2_xtb, basis_set, xc)
h2o_results, h2o_scf = scf_calculation(h2o_xtb, basis_set, xc)

Reaction energy#

The reaction energy is now taken as the energy difference between the reactants and products:

e_reactants = ch4_scf.get_scf_energy() + 2 * o2_scf.get_scf_energy()
e_products = co2_scf.get_scf_energy() + 2 * h2o_scf.get_scf_energy()

print(f"The reaction energy is {627.5*(e_products - e_reactants):.1f} kcal/mol")

The reaction energy is -157.1 kcal/mol

This can be compared to the results presented here, with combustion energies ranging from approximately -157.3 kcal/mol (obtained with Hartree-Fock and a small basis set), to a high-quality estimate of -193.2 kcal/mol.

In order to investigate the combustion energy of any other alkane, we merely need to change the corresponding SMILES string and calculating energy differences according to the general reaction equation:

\[ 2\textrm{C}_n\textrm{H}_{2n+2} + (3n+1)\textrm{O}_2 \rightarrow 2n\textrm{CO}_2 + 2(n+1)\textrm{H}_2\textrm{O} \]

With this, we obtain reaction energies for the first five alkanes, considering both the total energy and the energy per carbon atom:

Hydrocarbon

kcal/mol

kcal/mol per carbon

methane

-157.1

-157.1

ethane

-286.8

-143.4

propane

-415.4

-138.8

butane

-544.9

-136.2

heptane

-673.7

-134.7
Contents