Virtual undergraduate chemical engineering labs based on density functional theory calculations

2.1 Case 1: Calculation of the reaction rate constant for a unimolecular reaction

In this example, students are trained on how to calculate the reaction rate constant k(T) for the decomposition reaction of acetic acid (CH₃C(O)OH) in the gas phase:

The procedure requires the implementation of transition state theory (TST) and the Eyring equation (Truhlar et al., 1996):

The terms k _B (1.38 × 10⁻²³ J/K), h (1.05 × 10⁻³⁴ J/s), ∆S ^# (in J/mol K), ∆H ^# (in J/mol), and R (8.314 J/mol K) denote Boltzmann’s constant, Planks’ constant, entropy of activation, enthalpy of activation, and the universal gas constant, respectively. The equation describes how the reaction progresses at the molecular level. Thermodynamic functions are derived from partition functions of the electronic, vibrational, and rotational components, as described in statistical thermodynamics formulations (Gilson et al., 1997). The two Arrhenius parameters, the pre-exponential A factor, and the activation energy E _a can then be obtained from TST:

Experimentally, k(T) values are obtained (say, for a unimolecular irreversible elementary reaction) A → B by measuring the conversion, x, as a function of temperature:

First, students often miss the fact that the values of k(T) are only sensitive to the chemistry of the reaction itself, regardless of the type of reactor utilized. Nonetheless, a batch reactor model is often used to experimentally measure k(T) values. The concept of the transition state or the activated complex in the reaction kinetics can also be ambiguous for chemical engineering students, if not pre-covered in elementary physical chemistry. Using DFT to investigate reaction kinetics and mapping out mechanisms allows students to better grasp the physical meaning and importance of the transition state and how the underlying atomic movements determine the course of the reaction. Central to this aim is the display of atomic movements and their associated vibrational frequencies, a task that can be carried out by numerous graphical user interfaces, such as Avogadro and GaussView. Likewise, students are encouraged to appreciate the importance of the analysis of vibrational frequencies in determining the nature of a structure as a stable species (reactant or a product) or a transition state. In the latter, there should be only one imaginary frequency, which indicates the existence of a maxima energy point along the specified reaction pathway. The visualized atomic movements portray the pathway that the reactant(s) follow in their conversion into the product(s). Such important aspects could not be attained experimentally. An important point to highlight in this regard is the direct connection between macroscopic scales (i.e., moles) and microscopic scales. For example, the calculated heat capacity for a molecule from DFT applies to real-world quantities.

Figure 1 outlines the procedure for computing kinetic parameters for the dissociation reaction CH₃C(O)OH → CO₂ + CH₄. Tutorial 1.a in the Supplementary Materials (SM) provides a detailed description of the steps involved. In summary, Gaussian16 code is used to compute the energy profile of the reaction along the Reactant → Transition state → Products pathway. Vibrational frequencies and rotational constants from the Gaussian16 output files are then used to calculate thermodynamic values (namely, the enthalpies of formation ∆H _f ^o _(T) and entropies of formation ∆S _f ^o _(T)) for all involved species, including the transition state. This task will be carried out using the ChemRate code. Once the transition state is optimized, students are to view the calculated vibrational frequencies to confirm the existence of a single imaginary frequency in the obtained transition state, which amounts to 1097.7i. Visualizing atomic movements associated with this particular vibrational frequency confirms the dissociation of the CH₃C(O)OH into methane and water molecules. Figure 2 shows the potential energy diagram for the reaction. The computed rection enthalpy amounts to −27.7 kJ/mol and the activation enthalpy is 298.1 kJ/mol. As Tutorial 1 explains, the ChemRate code extracts thermodynamics values from which ΔH ^# and ΔS ^# are calculated:

where ∆H _f (TS), ∆H _f (R), ∆S _f (TS), and ∆S _f (R) denote the temperature-dependent values of the enthalpy of formation for the transition state, enthalpy of formation of reactant, entropy of formation of the transition state, and entropy of formation of the reactant, respectively. Table 1 enlists these values, whereas Figure 3 plots entropic and enthalpic values for the reaction and corresponding activation values. Finally, applying Equations (1) and (2) yield the two Arrhenius parameters as Figure 4 shows. The E _a value is computed to be 302.4 kJ/mol, a value that matches the corresponding experimental value of 305.0 kJ/mol (Mackie & Doolan, 1984).

Figure 1:

Computational procedure underlying Case 1.

Figure 2:

Plotted energy surface for the decomposition reaction CH₃C(O)OH → CO₂ + CH₄.

Figure 3:

Computed ΔS and ΔH (a) and their corresponding ΔS ^# and ΔH ^# values (b) for the decomposition reaction CH₃C(O)OH → CO₂ + CH₄. Entropic (S) and enthalpic (H) values are given in J/mol K and kJ/mol K; respectively.

Figure 4:

Computed Arrhenius parameters for the decomposition reaction CH₃C(O)OH → CO₂ + CH₄. ^a Elwardany et al. (2015).

It is always advisable to ask students to attempt to compare computed thermo-kinetic parameters with analogous experimental findings from the literature so as to appreciate the expected accuracy margin of the DFT calculations. Figure 4 also compares computed and experimental reaction rate constants (Elwardany et al., 2015) for this reaction. As shown in Figure 4, the difference between the two sets of values remains within 0.62–1.50 between 1200 and 1500 K.

In this regard, students can also survey the effect of the deployed DFT methodology on the accuracy of the results. In this example, the M062X DFT method (Zhao & Truhlar, 2008) along with the basis set of 6-311+G(d,p) (McLean & Chandler, 1980) were deployed, for instance. Students may repeat the procedure in this case using other simple reactions, such as H + CH₄ → H₂ + CH₃ or CH₂=CHCH₂Cl → CH₂=C=CH₂ + HCl. From a chemical engineering perspective, the calculated kinetic parameters can be used in the design equation of reactors (i.e., to find the length of a plug flow reactor necessary to attain a certain conversion). The literature provides comprehensive benchmarking of DFT-calculated thermodynamic parameters against analogous experimental values. For instance, we have previously illustrated that thermo-kinetic parameters for perfluorinated compounds could be accurately predicted by DFT calculations (Altarawneh et al., 2022). As an illustrative example, Table 2 contrasts computed and experimental values (Yu-Ran, 2002) of X–H (X = C, O, N) bond dissociation enthalpies in selected simple molecules. Detailed steps are given in Tutorial 1.b. The CBS-QB3 composite chemistry model is deployed in this case. Students were able to complete this task over three sessions (1 h each). The students appreciated the practical aspects of the case and started a discussion on how to compare kinetic petameters they typically obtain the laboratory segments of the Reaction class. Some of them indicated that a major obstacle in the design project in the final year would be to locate relevant sources of thermal and kinetic values. Hence, implementation of similar DFT calculations provide an alternative source of thermo-kinetic parameters.

2.2 Case 2: Finding the equilibrium constant for an isomerization reaction

All reactions are reversible, in principle. The direction of the reaction, either in the forward or the reverse direction, depends on its equilibrium constant, K _c (T), from the Gibbs free energy ∆Gr(T) of the reaction:

Tutorial 2 (in SM) presents a procedure for computing the K _c (T) values for the isomerization reaction n-butane ↔ i-butane based on Equation (5). The procedures involve using the standard enthalpies of formation of n-butane and i-butane and computed vibrational frequencies along with their rotational constants to derive thermochemical parameters (Gibbs free energies of formation ∆G _f ^o _(T) heat capacities, standard entropies, etc.) as a function of temperatures. Figure 5 plots the computed values of ∆Gr(T) and K _c (T) between 300 and 1000 K. Lower ∆Gr(T) prefers the forward direction of the reaction as the temperature increases. Students may attempt to calculate K _c (T) for other reactions, such as the decomposition of ammonia 2NH₃ ⇄ 3H₂ + N₂ or water–gas–shift reaction CO + H₂O ⇄ CO₂ + H₂. As demonstrated in Cases 1 and 2, the applied procedure affords a facile method to obtain thermochemical parameters. For example, Figure 6 presents the heat capacities of the i-butane molecule along with a comparison with pertinent experimental measurements (Pittam & Pilcher, 1972). The two sets of values are very close to each other. Supporting information provides the values and procedures applied in this case. The reaction rate parameters and equilibrium constants computed by DFT calculations often display excellent agreement with experimentally measured values. In a series of studies (Altarawneh, 2022; Altarawneh et al., 2022; Rawadieh et al., 2021), we thoroughly compared kinetic parameters obtained by DFT calculations with analogous experimental measurements. In many cases, the two set values reside within a factor of 1.0–3.0.

Figure 5:

Calculated equilibrium constant K _eq and the free energy ΔG _r of the reaction n-butane ↔ i-butane.

Figure 6:

Obtained heat capacity for the n-butane molecule in the gas phase. ^a Pittam and Pilcher (1972).

2.3 Case 3: Estimation of the standard enthalpy of formation (∆H _f ^o ₂₉₈.₁₅) of a compound

Students often find it challenging to locate thermochemical data at specific temperatures and pressures from the literature, or that such data simply do not exist. As Cases 1 and 2 demonstrate, DFT calculations are powerful for deriving an accurate estimation of thermochemical functions. Values of ∆H _f ^o ₂₉₈.₁₅ are determined and measured from calorimetry experiments. Alternatively, isodesmic reaction schemes are readily used to estimate ∆H _f ^o ₂₉₈.₁₅ values. Isodesmic reactions are hypothetical reactions that conserve the number and type of chemical bonds on two sides of the reaction. The procedure requires obtaining the standard enthalpy of the reaction from DFT calculations (∆H _{r 298}.₁₅) and then using this value to find a missing ∆H _f ^o ₂₉₈.₁₅ that is not experimentally available. Scheme 1 presents an illustrative example in which the ∆H _f ^o ₂₉₈.₁₅ value for the 1,4-dichlorobenzene molecule is computed. As Scheme 1 shows:

Scheme 1:

Isodesmic work reactions to obtain ∆H _f ^o ₂₉₈.₁₅ value for 1,4-dichlorobenzene.

There are 10 aromatic C–H and two C–Cl bonds in the products and the reactants side. The experimental values for the chlorobenzene and benzene molecules are known to be 54.2 and 82.9 kJ/mol, respectively. These two compounds are referred to as “reference compounds”. ∆H _{r 298}.₁₅ for the reaction from the DFT calculations was calculated to be 4.2 kJ/mol. Hence, utilizing the two known ∆H _f ^o ₂₉₈.₁₅ values with the computed ∆H _{r 298}.₁₅ gives 29.7 kJ/mol as the enthalpy of the formation of the 1,4-dichlorobenzene molecule. The students can then compare this computed value with the analogous experimental value of 24.6 kJ/mol (Platonov & Simulin, 1984) where a reasonable agreement is attained. Tutorial 3 in SM provides detailed procedures involved in Case 3. Students were able to successful complete this Case within an hour lab class. The instructor may advise students to consider other DFT methodologies or composite chemistry model to assess if agreement with the experimental value could be improved.

The procedure could be applied to more complex structures of pesticides and hydrocarbons with no available ∆H _f ^o ₂₉₈.₁₅ values in the literature. The most important task is to write an Isodesmic reaction where the reference species entail known experimental ∆H _f ^o ₂₉₈.₁₅ values. It is encouraged to ask students to compare. Literature provides satisfactory comparisons between the experimental and computed values of ∆H _f ^o ₂₉₈.₁₅. For example, we have shown that the standard enthalpy of formation of o-thiophenoxy departs by less than 4.0 kJ/mol from the corresponding experimental value (Altarawneh et al., 2012).

2.4 Case 4: Calculation of redox potentials

The accurate determination of redox potentials is important in electrochemical and corrosion applications. In DFT calculations, the procedure requires implementing the Born–Haber cycle. Herein, we apply this cycle to calculate the standard one-electron redox potential E _o /V based on the Nernest equation:

The standard Gibbs free energy change for the redox half ∆G _ox,(sol) is calculated from the free energy change in the gas phase ∆G _ox,(g), and the solvation free energies of the oxidized ∆G _solv,(III) and reduced ∆G _solv,(II) forms:

This example applies Equation (6) (the Born–Harber cycle (Heinz & Suter, 2004)) to calculate the half-redox potential of nickelocene, Ni(C₅H₅)₂, as shown in Scheme 2:

Scheme 2:

Born–Harber cycle for the II/III redox of Ni(C₅H₅)₂.

Using the B3LYP/Lanl2dz level of theory, Table 3 enlists the calculated values in the thermodynamic cycle. The calculated ∆G _ox,(sol) and E _o are 4.25 and 0.18 V, respectively. Students are encouraged to see how these values compare with available experimental estimates and to re-apply the cycle for other compounds and complexes, such as nitrobenzene and Co(C₅H₅)₂. Tutorial 4 in SM provides computed energies. As a benchmark of accuracy, values pertinent to the redox potential cycle (for copper complexes) obtained by DFT calculations were found to reside within 0.15 V from their analogous experimental values (Yan et al., 2016). Utilizing nickelocene in the example also encourages students to appreciate the effect of long-range interactions and oxidation states in metal–organic systems. Students may also revisit simple systems that involve the synthesis of hydrogen, as described by Matsui et al. (2013), where an equivalent procedure applies.

2.5 Case 5: Computing material properties

The properties of materials are deployed throughout the chemical engineering curriculum (i.e., unit operations, heat transfer, corrosion, and plant design). DFT calculations have become the prime research tool for acquiring the accurate thermal, optimal, electronic, and mechanical properties of materials. Thus, it is important to train students in the use of DFT to derive the properties of materials and to investigate their structural properties. The latter are directly linked to their catalytic properties, for instance. The CASTEP program, along with the VESTA visualization code, was utilized to investigate the structural (Figure 7a), electronic (Figure 7b), and thermal (Figure 7c) properties of ceria (CeO₂). The plotted density of states curve in Figure 5b determines the nature of the material as a conductor or semiconductor. A narrow band gap, for instance (1.0–3.0 eV), indicates semi-conducting behavior. As described in the literature, DFT calculations often yield very accurate band gaps for transition metal oxides (within 0.1–0.2 eV) in reference to experimental values (Widjaja et al., 2017).

Figure 7:

Structure (a), electronic (b), and thermal properties of ceria (CeO₂).

2.6 Case 6: Mapping out a surface-catalyzed mechanism

Deriving a rate law for surface-catalyzed reactions from tabulated experimental data is an important aspect in attaining the learning outcomes of reaction courses in chemical engineering. The obtained rate law from the experimental data was then used to suggest potential reaction steps and mechanisms. This procedure produces overall kinetic parameters, that is, not for individual steps (adsorption, surface reactions, reactions between a gas phase species and an adsorbed one, and desorption). Reaction and activation energies for these initial steps in the catalyzed reaction cycle can be readily obtained by DFT calculations. This allows us to obtain a detailed reaction mechanism from which a comprehensive kinetic model can be deduced. Figure 8 provides a mechanism for the surface hydrogenation of ethene into ethane, showing all involved reactions: the dissociation of a hydrogen molecule, the adsorption of the ethene molecule, surface hydrogen transfer to adsorption species, and, finally, the desorption of the ethane molecule. The idea is to convey to students the atomic-base description of the catalyzed reactions—an aspect that has often been overlooked in the teaching of catalysis to undergraduate chemical engineering students. Other examples may include the surface oxidation of CO into CO₂ or the synthesis of syngas from water and CO₂. The literature presents insightful success stories in which DFT calculations satisfactorily account for the salient mechanistic and kinetic features of surface-mediated reactions in the interpretation of experimental observations (Li et al., 2022).

Figure 8:

Surface catalyzed hydrogenation of ethene into ethane (C₂H₄ + H₂ → C₂H₆) over Pt. All values are in kJ/mol in reference to the initial reactants. TS denotes a transition state.