Assembling a full workflow

Assembling a full workflow#

As an example of constructing a full workflow for considering a modeling a chemical process, we 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.

import veloxchem as vlx
import veloxchem as vlx
import py3Dmol as p3d
import numpy as np
from rdkit.Chem import AllChem

Geometry optimization#

UFF optimization#

Construct molecules using SMILES strings and perform an UFF optimization.

ch4_uff = vlx.Molecule.smiles_to_xyz('C')
h2o_uff = vlx.Molecule.smiles_to_xyz('O')
co2_uff = vlx.Molecule.smiles_to_xyz('O=C=O')
o2_uff = vlx.Molecule.smiles_to_xyz('O=O')

xTB optimization#

Set up xTB driver and perform xTB optimization.

def xTX_opt(xyz_uff):
    molecule = vlx.Molecule.read_xyz_string(xyz_uff)
    xtb_drv = vlx.XtbDriver()
    method_settings = {'xtb':'gfn2'}
    xtb_drv.set_method(method_settings['xtb'].lower())
    xtb_grad_drv = vlx.XtbGradientDriver()
    xtb_opt_drv = vlx.OptimizationDriver(xtb_grad_drv)
    xtb_opt = xtb_opt_drv.compute(molecule, xtb_drv)
    return xtb_opt

ch4_xtb = xTX_opt(ch4_uff)
h2o_xtb = xTX_opt(h2o_uff)
co2_xtb = xTX_opt(co2_uff)
o2_xtb = xTX_opt(o2_uff)

Visualization#

Convert to format readable by py3Dmol and visualize structures.

ch4_xyz = ch4_xtb.get_xyz_string()
h2o_xyz = h2o_xtb.get_xyz_string()
co2_xyz = co2_xtb.get_xyz_string()
o2_xyz = o2_xtb.get_xyz_string()

viewer = p3d.view(viewergrid=(2, 2), width=400, height=300, linked=False)
viewer.addModel(ch4_xyz, 'xyz', viewer=(0, 0))
viewer.addModel(o2_xyz, 'xyz', viewer=(0, 1))
viewer.addModel(co2_xyz, 'xyz', viewer=(1, 0))
viewer.addModel(h2o_xyz, '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

Print distance matrices.

ch4_rdfit = AllChem.MolFromXYZBlock(ch4_xyz)
ch4_dm = AllChem.Get3DDistanceMatrix(ch4_rdfit)
o2_rdfit = AllChem.MolFromXYZBlock(o2_xyz)
o2_dm = AllChem.Get3DDistanceMatrix(o2_rdfit)
co2_rdfit = AllChem.MolFromXYZBlock(co2_xyz)
co2_dm = AllChem.Get3DDistanceMatrix(co2_rdfit)
h2o_rdfit = AllChem.MolFromXYZBlock(h2o_xyz)
h2o_dm = AllChem.Get3DDistanceMatrix(h2o_rdfit)
print('CH4:\n',np.around(ch4_dm,3))
print('O2:\n',np.around(o2_dm,3))
print('CO2:\n',np.around(co2_dm,3))
print('H2O:\n',np.around(h2o_dm,3))

CH4:
 [[0.    1.082 1.082 1.082 1.082]
 [1.082 0.    1.767 1.767 1.767]
 [1.082 1.767 0.    1.767 1.767]
 [1.082 1.767 1.767 0.    1.767]
 [1.082 1.767 1.767 1.767 0.   ]]
O2:
 [[0.   1.21]
 [1.21 0.  ]]
CO2:
 [[0.    1.144 2.287]
 [1.144 0.    1.144]
 [2.287 1.144 0.   ]]
H2O:
 [[0.    0.959 0.959]
 [0.959 0.    1.544]
 [0.959 1.544 0.   ]]

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(xyz, basis_set, xc, multiplicity = None):
    molecule = vlx.Molecule.read_xyz_string(xyz)
    basis = vlx.MolecularBasis.read(molecule, basis_set)
    if multiplicity:
        molecule.set_multiplicity(multiplicity)
        scf_gs = vlx.ScfUnrestrictedDriver()
    else:
        scf_gs = vlx.ScfRestrictedDriver()
    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_xyz, basis_set, xc)
o2_results, o2_scf = scf_calculation(o2_xyz, basis_set, xc, multiplicity = 3)
co2_results, co2_scf = scf_calculation(co2_xyz, basis_set, xc)
h2o_results, h2o_scf = scf_calculation(h2o_xyz, 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 in this paper, which discusses the question of computational cost versus accuracy with methane combustion as an example. They obtain combustion energies ranging from approximately -157.3 kcal/mol when using Hartree–Fock and a small basis set, which is compared to a reference estimates of -193.2 kcal/mol, obtained using high-level theory.

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