Group Additivity for Estimating Thermochemical Properties of Furanic

Article pubs.acs.org/IECR

Group Additivity for Estimating Thermochemical Properties of Furanic Compounds on Pd(111) Vassili Vorotnikov, Shengguang Wang, and Dionisios G. Vlachos* Catalysis Center for Energy Innovation, Department of Chemical and Biomolecular Engineering, University of Delaware, Newark, Delaware 19716, United States S Supporting Information *

ABSTRACT: We analyze the adsorption of 101 furan chemistry-related adsorbates, including intermediates with varying levels of hydrogenation of the furan ring as well as ring-open compounds, on Pd(111). The standard heat of formation (Hf,298), standard absolute entropy (S298), and heat capacity (Cp,T) are estimated using dispersion-corrected density functional theory (DFT-D3) and statistical mechanics formulas, and all Hf,298 values are referenced to a single set of gas-phase molecules to ensure consistency. We further estimate the dispersion contribution, ΔEdisp, to the heat of adsorption. Subsequently, we develop a group additivity scheme to quickly estimate thermochemical properties of furanic molecules. For the proposed model, the vibrational and rotational contributions of gas-phase-like groups are estimated using G4-level calculations. We find that the group additivity scheme developed for open-chain molecules is inaccurate for furanics. We report 17 new groups, involving heteroatom−nearest neighbor interactions and four new corrections that account for furan ring deformation and the level of ring hydrogenation. The mean deviations from DFT-computed values are 1.8 kcal/mol in Hf,298, 0.3 kcal/mol in ΔEdisp, 0.9 cal/(mol K) in S298, and 0.7 cal/(mol K) in Cp,800. The largest deviations are observed in highly saturated adsorbates with multiple gas-phase-like surface groups. We further introduce a nine-parameter heteroatom-based model for estimating ΔEdisp, resulting in a mean deviation from DFT-computed values of 1.4 kcal/mol.

■

INTRODUCTION Biomass is considered a plausible alternative to petroleum for the production of fuels and chemicals. Its conversion involves hydrolysis to produce simple sugars and lignin derivatives or pyrolysis to produce bio-oil. Furanic compounds, such as furfural or 5-hydroxymethylfurfural (HMF), are produced as pyrolysis byproducts or through the isomerization and dehydration of sugars.1−3 Their selective catalytic upgrade can result in fuel-grade products, such as methylfuran or dimethylfuran.4 However, the mechanisms of these processes are yet to be understood, and as a result, the discovery of cheaper, more active, and more selective catalysts remains an art. Among single metal catalysts, palladium (Pd) exhibits interesting catalytic behavior. Reactor studies have shown that furfural can undergo either decarbonylation to furan or hydrogenation to furfuryl alcohol and fully saturated furanring compounds.5,6 Surface science experiments have further shown that furfural can ring-open and furfuryl alcohol can be converted to methylfuran on Pd(111).7 Understanding the reaction mechanism for the conversion of furans on Pd can therefore open doors to rational catalyst design of these processes on other metals. Computational methods are proving valuable in revealing mechanistic details in a reaction network.8 Density functional theory (DFT) provides thermodynamic and kinetic properties for molecule−metal interactions.9,10 However, the number of intermediates and possible reactions in a surface reaction network increases drastically with the number of heteroatoms in a molecule, making it impractical to perform DFT calculations for all biomass-relevant species. To illustrate this point, Figure 1 shows the number of elementary reactions and © 2014 American Chemical Society

Figure 1. Number of reactions and intermediates in furfural reaction network considering reactions on the side group only (blue), side group with furan ring hydrogenation (green), and side group with furan ring hydrogenation and opening (red). The horizontal axis corresponds to the total number of reactions (leftmost columns), the number of various reaction types, and the intermediates (rightmost columns).

Received: Revised: Accepted: Published: 11929

May 19, 2014 July 2, 2014 July 7, 2014 July 7, 2014 dx.doi.org/10.1021/ie502049a | Ind. Eng. Chem. Res. 2014, 53, 11929−11938

Industrial & Engineering Chemistry Research

Article

ab initio simulation package (VASP) version 5.2.12.28,29 The electron−electron exchange and correlation energies were computed using the Perdew, Burke, and Ernzerhof functional with the latest dispersion correction, PBE-D3.30,31 The projector augmented-wave (PAW) method was used for the electron−ion interactions.32,33 We used a plane-wave basis set with an energy cutoff of 400 eV. For bulk calculations, a tetrahedron method with Blochl corrections and 15 × 15 × 15 Monkhorst−Pack k-point mesh was used.34,35 The bulk lattice constant was obtained using the Birch−Murnaghan equation of state.36,37 The Pd fcc lattice constant was calculated to be 3.90 Å, in agreement with the experimental value of 3.89 Å.38 A set of gas-phase VASP calculations was used for referencing the adsorbates’ heats of formation. These calculations were performed with a 20 × 21 × 22 Å supercell and a gamma point. Further information on referencing for thermodynamic consistency can be found in the Supporting Information, section A. The metal slab was modeled with a 4 × 4 unit cell composed of four atomic layers. The bottom two layers were frozen. The vacuum between the slabs was set at 20 Å. The Brillouin zone was sampled with a 3 × 3 × 1 k-mesh. For accurate total energies, we used the Methfessel−Paxton method with a smearing parameter of 0.1. An RMM-DIIS quasi-Newton algorithm was used in the energy minimization.39 Surface relaxation was performed until all forces were lower than 0.05 eV/Å. In previous works, we showed that using a stricter criterion of 0.02 eV/Å does not result in significant change in the adsorbate conformation or total energy.9 The dispersion energy contribution to the heat of adsorption, ΔEAdisp * , was computed as the difference between the adsorbate− slab complex and the sum of bare slab with the adsorbate in the gas phase, similar to a typical adsorption energy calculation.

intermediates involved in the reaction of furfural over a single metal catalyst. The number of intermediates increases from 34 to 605 to 6661 (as estimated using Rule Input Network Generator (RING) software described in the Supporting Information, section E), as reactions are being added, starting from the chemistry on the side group of the furan ring to ring hydrogenation and ring-opening reactions, and the total number of reactions exceeds 104. In a recent DFT study of furfural reaction over Pd(111) excluding ring hydrogenation and ring opening,9 the computational time needed to estimate all of the thermodynamic and kinetic parameters exceeded 1 million CPU h on the TACC Ranger cluster. Even with the increased efficiency (∼5 times) of the new TACC cluster, Stampede, the requirement for DFT-based estimates of all properties in the full furfural reaction scheme over Pd(111) would exceed 39 million CPU h. For reference, each of the Dell C8220 compute nodes used on the Stampede cluster consists of two eight-core Xeon E5-2680 (Sandy Bridge) processors running at 2.7 GHz.11 Similar challenges exist in modeling other compounds, such as 5-hydroxymethylfurfural (HMF). Given this computational bottleneck, quick estimation of thermochemical properties for large molecules on metal surfaces, e.g., standard enthalpy of formation (Hf,298 ° ), standard absolute entropy (S298 ° ), and standard heat capacity at various temperatures (C°p,T) is essential. In this work, we develop a group additivity scheme for the family of furan compounds and reassess the accuracy of the group additivity scheme developed for small, open-chain molecules. We investigate the sources of deviation from DFT calculations and report new groups and corrections for the scheme applicable to furanic chemistry. We demonstrate the framework on Pd(111).

■

METHODS Group Additivity. The group additivity8,12−17 was introduced by Benson in the 1950s and applied to gas-phase thermochemistry by regressing thermochemical properties of interest (Hf,298 ° , S298 ° , and Cp,T ° ) over experimental data.12−14 The same scheme was subsequently extended to liquid-phase compounds.18 With computing advances, the use of ab initio methods allowed the development of group additivity for difficult-to-estimate species, such as radicals and methoxyfurans.19−25 Most recently, the scheme has been extended to surface adsorbates, first by Goddard and co-workers17 and then by Vlachos and co-workers.15,16 The simplest formulation of this model is shown in eq 1, where a thermochemical property of any compound, pi, is expressed as the sum of its heteroatomcentered groups’ contributions to that property, xk.

A(g) A* A + slab slab ΔEdisp = Edisp − Edisp − Edisp

Estimation of Thermochemical Properties from DFT. We use statistical mechanics formulas to obtain the values of H f,298 ° , S298 ° , and Cp° (T). Using frequencies from DFT calculations, we estimate S°298 and C°p (T) via eqs 3 and 4, respectively.

i Ngroups

pi =

∑ k=1

g i , k xk

(2)

⎛ x ⎞ ° = R∑⎜ x ST° = Svib − ln(1 − e−x)⎟ ⎝e − 1 ⎠ i

(3)

⎛ x 2 e −x ⎞ C p°(T ) = Cv°,vib = R ∑ ⎜ ⎟ ⎝ (1 − e−x)2 ⎠ i

(4)

x = Θνi /T (1)

Θνi =

Here, gi,k represents a coefficient associated with the number of times a kth group is counted in species i, and typically has an integer value. In a more complex formulation involving group corrections, gi,k can be any rational number. Nigroups represents the number of heteroatom-centered groups in compound i. For an illustration of the group additivity concept, see the Supporting Information, section F. Electronic Structure Calculations. DFT calculations pertaining to gas-phase group additivity values were carried out using the Gaussian 09 package (revision A.2)26 to obtain electronic energy and vibrational frequencies.27 Adsorbatepertaining DFT calculations were carried out using the Vienna

hνi k

where ⎛1⎞ νi [s−1] = 100⎜ ⎟c ⎝ λi ⎠

In these equations, 1/λi is the wavenumber in cm−1, c is the speed of light in m/s, and x is the reduced vibrational frequency. In the group additivity works of Salciccioli et al., the heat of formation of an adsorbate species CxHyOz was computed with reference to the most hydrogenated species CxHsatOz in the gas 11930

dx.doi.org/10.1021/ie502049a | Ind. Eng. Chem. Res. 2014, 53, 11929−11938


Article

phase.15,16,40 This method required multiple gas-phase reference molecules. To avoid varying reference molecules, we use a single set of accurate gas-phase reference species taken from the NIST database.41 Then, the heat of formation of any compound i, Hif,298, may be estimated from DFT using eqs 5−9. i i Hf,298 = H298 + n iA−1h

⎡1 0 4 ⎤−1 ⎢ ⎥ = −96416.625 + [2 1 4]⎢ 0 0 2 ⎥ ⎣0 1 2 ⎦ ⎡−17.89 + 25374.936 ⎤ ⎢ ⎥ ⎢ 0 + 730.356 ⎥ ⎢⎣−57.799 + 47909.848 ⎥⎦

(5)

i slab + i slab i H298 = E0K − E0K + E(vib,298) i i i H298 = E0K + Evib,298

CH3CHO CH3CHO Hf,298 = H298 + nCH3CHOA−1h

= −41.6 kcal/mol

(for adsorbates)

(for gas‐phase species)

compared to the experimental value of −40.8 ± 0.4 kcal/mol.43 The computed estimate is in line with the error expected from the G4 method (see Supporting Information, section A, for further details on the calculation of Hif,298).27 Much larger errors occur when using VASP for estimating gas-phase reference values. Estimation of Group Additivity Values for Furanic Thermochemistry. In order to use group additivity, we first need to parametrize each of the possible groups of the adsorbates involved in the chemistry. Similar to previous works,15,16 we use DFT parametrization in order to estimate each of the group contributions, xk, in eq 1. Given M adsorbates with a total of N groups (M ≥ N for exact or overdetermined system), it is possible to estimate each group’s contribution to a specific property by solving a system of linear equations

(6)

n i = [ NC NO NH ]

(7)

ref1 ⎤ ⎡ H ref1 − H298 ⎢ f,298,NIST ⎥ ⎢ ref2 ref2 ⎥ h = ⎢ Hf,298,NIST − H298 ⎥ ⎢ ref3 ref3 ⎥ ⎢⎣ Hf,298,NIST − H298 ⎥⎦

(8)

G(M × N )x (N × 1) = p(M × 1)

Here, Hi298 is the temperature-corrected DFT energy of the ith adsorbate or ith gas-phase species, niA−1h is the reference adjustment, Eslab+i 0 K is the zero-Kelvin total energy of the slab with the adsorbate, Eslab 0 K is the zero-kelvin energy of the bare slab, and Eivib,298 is the temperature contribution from the vibrational frequencies that includes the zero-point energy correction. The vector ni contains the atomic composition of compound i. The vector h contains the energetic difference between the NIST heat of formation and the one computed purely using DFT (VASP or Gaussian) for each of the reference molecules. There are only three reference molecules, because each compound of interest here is made up of only three atoms (C, O, and H). Lastly, the matrix A contains the atomic composition of each reference compound. The only criterion that needs to be satisfied is that the reference species must have linearly independent atomic compositions to ensure that any compound can be referenced. The criterion is equivalent to the rank of A being full. One simple way to accomplish this is to select the enthalpies of formation of atoms (C, O, H) as references for computing the heats of formation in assessing higher-level methods.27,42 As an example, let us compute the heat of formation of gasphase acetaldehyde (CH3CHO) using CH4, H2, and H2O as a 3CHO = −96 reference set. Gaussian (G4) calculation gives HCH 298

(1 0 4; 0 0 2; 0 1 2). The resulting calculation of CH CHO

(11)

Here G is the coefficient matrix whose elements gm,n represent the number of times that the nth group appears in the mth adsorbate. The vector p represents the thermochemical property obtained from DFT and statistical mechanics. Lastly, the unknown vector x represents each group’s contribution to the desired thermochemical property. Each element in G can be obtained from graphical and molecular symmetry considerations. In this work, we utilized the Rule Input Network Generator (RING) software to autoextract the scaling coefficients from a SMILES-formatted adsorbate connectivity information.44,45 Importantly, this software allows us to easily add corrections (see below). Material related to automated group and group correction identification is included in the Supporting Information, section E. Equation 11 is used to determine the group additivity values and corrections for surface groups involving weak or strong atom−metal interaction. As in previous works,15,16 we take the vibrational and rotational contributions of gas-phase-like groups to be those of their gas-phase group additivity values. The difference in this work is that all of the gas-phase group additivity values were obtained using Gaussian G4-level calculations rather than experimental calorimetry data used by Benson (for further details see Supporting Information, sections B and C).12,13 This allowed us to keep track of the regressed thermochemical properties, specifically the rotational and vibrational contributions to H°f,298, S°298, and C°p,T, along with estimating groups, which do not exist in the literature, in a consistent way. The details of these calculations together with a complete set of gas-phase group additivity values can be found in the Supporting Information (section C).

416.63 kcal/mol. The composition vector nCH3CHO is (2 1 4). Href 298 is −25 374.94, −730.36, and −47 909.85 kcal/mol for CH4, H2, and H2O, respectively. The corresponding Href f,298,NIST values are −17.89, 0, and −57.80 kcal/mol.43 Lastly, the r e f e r e n c e c o m p o u n d c o m p o s i t i o n m a t r i x , A , is Hf,2983

(10)

is 11931



■

Article

RESULTS Adsorption of Furanic Compounds on Pd(111). Group additivity is a multiparameter model and, as such, it requires a large number of data to obtain statistically relevant values for the various groups. We therefore started with a rather large set of calculations from previous works consisting of 20 furanic species, 15 ring-open intermediates, and 13 small molecules involved in cleaving off the side group of furfural.9,10 Following the identification of the groups involved in these surface adsorbates, we found that further DFT calculations were required to obtain group additivity values while maintaining linear independence of the groups. We therefore performed 53

intermediates in elementary reactions, such as furanic hydrogenation, dehydrogenation, dehydroxylation, decarbonylation, and ring-opening reactions. All 101 DFT-computed structures are shown in Supporting Information, section G. While adsorption conformations vary significantly, several features remain consistent with the small molecules’ studies.15,16 Where possible, the most preferred surface configurations have all C atoms sp3-hybridized. Any unsaturated C atom tends to fulfill the sp3 character by forming additional C−M bonds with the metal (M) surface. However, the sp3 hybridization is not possible for some species, especially the highly unsaturated/dehydrogenated ones, due to a mismatch between the furanic ring and the lattice. Such structures tend to result in lower stability or tilting with respect to the surface. As a result, many unsaturated species relaxed to a preferred configuration with carbon−carbon double bonds present in the adsorbate, indicating connectivity not typical of small molecules on similar metal surfaces.15,16 The presence of such sp2-hybridized C atoms was handled by assigning these heteroatoms a “Cd” notation, consistent with the original work of Benson and the gas-phase group additivity developed here (Supporting Information, section C). These groups were treated exactly the same with the addition of connectivity to a metal atom, M, where needed. As a result, we identified 17 new groups, involving heteroatom−nearest neighbor interactions, mainly arising from the presence of carbon atoms with sp2 character. Corrections Specific to Furanic Compounds. The scheme above allows standard groups, describing the heteroatom−nearest neighbor interactions only. In practice, even the original Benson group additivity scheme includes correction terms, as shown in eq 12. Here, gi,l represents a coefficient associated with the number of times a lth correction appears in species i. Nicorrections represents the number of corrections in the ith adsorbate. xl represents the contribution of a single lth correction to the property, pi.

Figure 2. Types of species used in the group additivity regression set. The top row shows typical furan-ring-containing compounds. The second row shows ring-open intermediates. The bottom row shows hydrogenated and dehydrogenated furan ring structures.

additional calculations (101 adsorbates in total) to obtain statistically significant group additivity values. The three types of species considered in this work are shown in Figure 2: furanic adsorbates with no modification to the furan ring, ringopen derivatives, and adsorbates with varying degrees of furan ring saturation. Of these, 14 species were ring-open intermediates and 39 were furanic species, all with varying levels of carbon and oxygen saturation to account for

Table 1. Corrections Introduced in the Furanic Group Additivity Scheme on Pd(111) Together with Their Corresponding Group Additivity Valuesa

a

Examples of groups are in circles. The number of adsorbates, in which each correction appears, is indicated in the “#” column. 11932



Article

Figure 3. Parity plots for group-additivity- vs DFT-estimated thermochemical properties of furanic derivatives adsorbed on Pd(111): (a) Hf,298, (b) S298, (c) Cp,300, and (d) ΔEdisp. In this regression, the identified group additivity corrections are included. Fitting statistics are included in green; AD stands for absolute deviation. i Ngroups

pi =

∑ k=1

tation), are highly strained, and (c) there is a large energetic difference between the simple furanic adsorption and oncehydrogenated or once-dehydrogenated ring species. The latter suggests that not only is there “surface-ring strain” in the adsorbates but also an additional effect due to the deformation of a stable furan ring. To address these differences for furanic adsorbates, we introduce four new second-order corrections: (i) weak H−M interaction, (ii) furan ring deformation, (iii) adsorbed furan ring, and (iv) C(M)2−C(M)−C(M) strain. These corrections are presented in Table 1 together with examples and the regression-estimated contributions to thermochemical properties. We ensured enough adsorbates were included to regress the newly identified groups. The number of adsorbates, “#”, in which each correction appears, is also reported in Table 1. The “weak H−M” correction is introduced to reduce the error from using gas-phase-derived group contributions and accounts for weak but additive interactions with the surface. This interaction contributes −3.5 kcal/mol to Hf,298. Antony et al. studied the adsorption of saturated n-alkanes on Pd(111), also using the DFT-D3 method, and estimated the H−Pd energetic contribution to the adsorption to be −2.4 kcal/mol,46 in reasonable agreement with the regression-based result reported here.

i Ncorrections

g i , k xk +

∑ l=1

g i , l xl

(12)

Previous group additivity work on Pt(111) was focused on developing the simplest additivity rules to keep consistent with the Benson scheme.16 A key correction was introduced to account for the cyclic strain of a species binding with various heteroatoms to the surface (multidentate adsorption). Salciccioli et al. observed that molecular ring strain increases with decreasing number of atoms in the ring and related it to the deformation of natural bond angles. They further quantified this “surface-ring strain” correction, allowing for its estimation without complex geometric information.16 By applying this model to furan-related compounds, we found large maximum deviations in all thermochemical properties: 27.6 kcal/mol in Hf,298, 17.6 cal/(mol K) in S298, and 5.9 cal/(mol K) in Cp,800. Upon closer examination of the adsorbates with the largest deviations, we found certain connectivity patterns that distinguish these molecules from the small ones. We found that (a) furanic compounds with multiple highly saturated carbon atoms have additional weak H−metal (H−M) interactions, (b) highly dehydrogenated ring intermediates, such as C1(Pd)2CH(Pd)C(Pd)2C(Pd)(C(Pd)O)O1 or O CHC(Pd)2CH(Pd)CH(Pd)C(H)O∼(Pd) (in SMILES no11933



Article

Table 2. Surface Group Additivity Values of Furanic Chemistry on Pd(111)a group

Hf,298

ΔEdisp

S298

Cp,300

Cp,400

Cp,500

Cp,600

Cp,800

Cp,1000

Cp,1500

#

C(C)(CO)(H)(M) C(C)(CO)(M)(O) C(C)(CO)(M)2 C(C)(COwk)(H)(O) C(C)(Cd)(H)(M) C(C)(Cd)(M)2 C(C)(H)(M)2 C(C)(H)2(M) C(C)(M)3 C(C)2(H)(M) C(C)2(M)(O) C(C)2(M)2 C(CO)(H)(M)(O) C(CO)(M)2(O) C(Cd)(CO)(M)(O) C(Cd)(H)(M)(O) C(Cd)(M)2(O) C(Cd)(M)3 C(H)(M)2(O) C(H)(M)3 C(H)2(M)(O) C(H)2(M)2 C(H)3(M) C(M)3(O) C(M)4 CO(C)(M) CO(Cd)(M) CO(M)2 COwk(C)(H) COwk(C)(M) COwk(Cd)(H) COwk(H)(M) Cd(C)(M) Cd(Cd)(M) Cd(M)(O) H(M) O(C)(M) O(Cd)(M) O(H)(M) O(M)2 Owk(C)(H) Owk(C)2 Owk(H)2 group corrections

2.93 8.92 2.47 −15.33 −29.20 −21.85 6.86 −1.28 17.26 −4.04 −0.75 17.85 −10.12 9.37 −17.28 −18.76 5.34 7.49 0.12 −0.55 −8.78 1.64 −9.31 −14.49 19.56 −43.46 −9.33 −68.81 −43.88 −44.32 −24.94 −49.57 31.76 8.46 −17.65 −15.15 −30.93 −18.62 −46.43 −36.43 −52.60 −41.95 −65.85 Hf,298

−11.82 −10.84 −2.80

−1.67 −3.36 −6.58 −3.74 −2.58 −3.18 −4.96 −4.79 ΔEdisp

−0.58 −3.91 −2.68 −4.24 0.79 −2.90 0.47 3.05 −5.35 0.71 −3.75 −1.73 −0.14 0.73 9.87 9.19 −6.75 −4.93 2.12 2.54 5.42 7.08 12.36 0.60 1.66 7.89 7.33 8.94 11.84 7.64 −6.10 9.30 −18.01 1.67 −5.96 0.38 4.44 5.69 6.73 2.97 8.32 4.04 12.20 S298

7.64 7.94 5.28 6.17 4.96 4.81 6.26 6.26 7.45 4.91 3.19 4.32 2.43 2.95 7.65 5.32 3.96 −0.85 5.58 7.31 8.58 9.71 11.45 9.96 5.42 3.96 4.10 10.60 9.34 7.81 0.87 11.80 3.54 3.45 3.42 3.42 3.84 2.76 9.25 6.45 6.35 13.41 10.66

9.54 9.59 6.62 6.75 6.70 6.13 7.95 8.11 9.18 6.65 4.71 5.60 4.07 3.92 8.53 6.61 5.21 0.33 6.99 8.87 10.26 11.46 13.26 10.83 6.32 4.94 4.94 11.25 10.60 8.79 1.28 13.28 5.27 4.17 4.63 4.54 4.21 2.85 10.21 7.04 7.06 14.32 11.35

10.84 10.52 7.42 6.94 8.04 6.99 9.03 9.50 10.10 7.86 5.57 6.32 5.21 4.48 9.35 7.78 5.82 0.93 7.97 9.82 11.57 12.67 14.75 11.35 6.83 5.62 5.40 11.71 11.76 9.54 1.53 14.38 5.95 4.69 5.39 5.40 4.50 2.93 10.75 7.34 7.69 14.97 11.86

12.88 11.52 8.34 6.75 10.25 8.00 10.67 12.15 10.92 9.79 6.43 7.02 7.02 5.03 10.73 9.86 6.39 1.46 9.64 11.27 14.10 14.89 17.96 12.02 7.48 6.74 6.05 12.61 14.29 10.84 1.70 16.41 6.43 5.52 6.09 6.74 5.09 3.26 11.55 7.70 9.05 16.00 13.02 Cp,800

13.59 11.70 8.51 6.54 11.01 8.19 11.27 13.28 10.99 10.48 6.55 7.11 7.68 5.10 11.11 10.61 6.50 1.50 10.28 11.83 15.20 15.89 19.50 12.20 7.64 7.12 6.27 12.97 15.37 11.28 1.67 17.22 6.52 5.76 6.14 7.14 5.29 3.42 11.90 7.79 9.67 16.27 13.69 Cp,1000

14.54 11.74 8.59 6.07 12.05 8.28 12.15 15.02 10.84 11.44 6.59 7.08 8.66 5.14 11.54 11.71 6.57 1.49 11.25 12.69 16.89 17.51 22.04 12.37 7.81 7.60 6.56 13.42 16.89 11.84 1.59 18.39 6.48 5.99 6.06 7.57 5.49 3.67 12.57 7.87 10.66 16.44 15.07 Cp,1500

1 10 4 2 3 6 13 10 5 61 22 13 1 2 2 2 3 1 4 1 4 1 1 2 1 13 2 1 23 2 1 1 1 2 1 1 10 1 1 1 4 3 1 #

ads furan ring furan ring def surface-ring strain weak H−M C(M)2−C(M)−C(M)

−2.40 2.52 −3.14 −3.51 6.46

0.12 0.38 5.30 −0.28 0.22

2.13 −0.01 −5.52 0.90 0.59

2.28 0.03 −5.32 0.87 0.64

2.50 0.08 −4.73 0.80 0.70

7 16 66 18 27

−3.11 −2.41 −3.91 −4.10 −5.07 −4.41 −4.07 −4.13 −3.85 −2.84 −5.25 −6.55 −3.20 −2.70 −1.93 −4.52 −3.35 −4.97 −2.41 −1.68 1.51 −10.17 −4.12 −5.63 −4.81 5.12 −4.78 −10.26

1.56 −0.05 13.20 −0.90 −0.21

Cp,300

Cp,400

Cp,500

11.73 11.04 7.89 6.93 9.01 7.51 9.74 10.58 10.56 8.71 6.04 6.71 6.00 4.78 9.98 8.67 6.12 1.23 8.68 10.45 12.59 13.56 15.99 11.67 7.14 6.11 5.69 12.07 12.76 10.10 1.66 15.21 6.22 5.06 5.79 6.00 4.74 3.04 11.10 7.52 8.22 15.43 12.29 Cp,600

1.89 −0.10 −2.36 0.75 0.28

1.83 −0.11 −4.25 0.86 0.37

1.87 −0.08 −5.15 0.91 0.44

1.95 −0.06 −5.48 0.92 0.50

Common groups and second-order corrections are listed. The units are kcal/mol for Hf,298 and ΔEdisp, cal/(mol K) for S298 and Cp,T. The number of adsorbates that each group appears in, represented by “#”, is listed in the last column. a

intermediates. To account for this difference, we introduced two corrections. The first group correction, “adsorbed furan ring” applies to any adsorbate with four C−M bonds, each one being essentially sp3-hybridized, as shown in Table 1. The Cα− Cβ bond distances of gas-phase furan rings are 1.35−1.36 Å, whereas those of adsorbed furan rings are ∼1.46 Å with neartetrahedral bond angles.47 The value for surface furanics is indicative of furanic carbons’ rehybridization to sp3. Similar

In Benson’s gas-phase group additivity scheme, the gas-phase furan ring is a cyclic group correction with negative enthalpic contribution (−5.8 kcal/mol),14 effectively contributing to the stability of the related furan molecules, including furfural, 5hydroxymethylfurfural, and methylfuran. The stability of furan, arising from its aromaticity, also plays a considerable role in surface adsorption as manifested from the large energetic difference between the furans and the once (de)hydrogenated 11934



Article

Figure 4. Parity plots for heteroatom-based model vs DFT-estimated quantities. The heteroatom-based model is not adequate for (a) Hf,298 estimation, but is reasonable for (b) ΔEdisp. Fitting statistics are included in green; AD stands for absolute deviation.

trends are observed in the adsorption of linear olefins on Pd(111), with C−C bond distances of 1.45−1.47 Å for the di-σ configurations.48 Similar to its gas-phase counterpart, the “adsorbed furan ring” correction has a negative enthalpic contribution, implying that furans are highly stable as adsorbed molecules. The “adsorbed furan ring” correction may be related to the stability of other aromatic species on the surface, such as benzene; however, we cannot generalize this as other aromatics were not computed in this work. The second correction, referred to as “furan ring deformation”, accounts for the level of hydrogenation of the ring. Upon ring dehydrogenation, the dehydrogenated carbon lacks the desired sp3 hybridization due to steric effects, thus introducing some strain. The more dehydrogenated the molecule is, the greater the effect. Similarly, upon furan ring hydrogenation, the molecular plane is broken leading to intramolecular strain. However, this effect is not as pronounced at higher levels of hydrogenation as the adsorbate detaches from the surface and attains its gas-phase preferred geometry.10 In order to account for all possible combinations of the “furan ring deformation”, we (i) wrote a RING code to locate all of the fragments pertaining to deformation and (ii) assigned a group contribution to each of these special fragments. These group contribution assignments are integers used in the fitting matrix and are reported in the Supporting Information (section D). For instance, a twice-dehydrogenated furan ring would include 2 times the “furanic ring deformation” correction, whereas the one shown in Table 1 would only include one such correction. This enthalpic contribution is positive, contributing to the destabilization of a particular intermediate. The last correction accounts for “C(M)2−C(M)−C(M)” type groups. Here, the “surface-ring strain” accounts for multiple neighboring heteroatoms. Such a high-strain structural fragment has a substantial destabilizing enthalpic effect of 6.46 kcal/mol. The results of regression suggest that, in addition to the standard group-additivity-estimated reaction energy, there is an additional endothermic effect of ∼4.9 kcal/mol in either initial hydrogenation or dehydrogenation of the surface furan ring, which can be expected for particularly stable aromatic molecules.

Figure 3 shows the parity plots using these corrections. We were able to reduce the maximum deviations in all thermochemical properties to 5.5 kcal/mol in Hf,298, 1.5 kcal/ mol in ΔEdisp, 3.7 cal/(mol K) in S298, and 3.6 cal/(mol K) in Cp,800. The related group additivity values are given in Table 2. Using these values, one may estimate thermochemical properties of any furanic intermediate on Pd(111). For example, one can estimate the heat of formation of dimethylfuran. Provided it binds similarly as other furans, its adsorbed structure would correspond to a SMILES string “CC1(Pd)C(Pd)C(Pd)C(Pd)(C)O1”. With the help of the group identification code, we find the following groups: 2 C(C)(H)3, 1 O(C)2, 2 C(C)2(H)(M), 2 C(C)2(M)(O), 1 “adsorbed furan ring”, a contribution of 0.65 from the bidentate-type “surface-ring strain”, and an additional 4 “weak H−M” interactions from the methyl groups. Using the values from Table 2, we estimate the value for the heat of formation as −75.9 kcal/mol, compared to a DFTcomputed value of −75.0 kcal/mol. Heteroatom-Based Regression Model for Hf,298 and ΔEdisp. If a group additivity value is missing from the database (i.e., some specific connectivity was not a part of the regression set), one uses a similar group value in estimating a thermochemical property. Using a substitute group is commonly referred to as “remapping”. However, when the group has not been encountered in the original fitting set, group swapping results in sacrificing accuracy. To avoid making postfitting assumptions, we remapped each of the established groups to its simplest connectivity to the metal: a carbon atom can be saturated or bound to one to three metal atoms, while an oxygen atom can be saturated, weakly bound,15,16 or have one to two bonds with the surface. In doing so, we reduce the number of fitting parameters from ∼40 to 9, requiring fewer DFT calculations. We regress the same data set for both Hf,298 and ΔEdisp using this heteroatom-based model, and the resulting parity plots are shown in Figure 4. We can further compare the performance of this model to group additivity. As shown in Figure 3a,d, the performance of group additivity is satisfactory. For the full group additivity scheme, the mean deviation in Hf,298 is 1.8 kcal/mol and the maximum deviation is 5.5 kcal/mol, and the mean deviation in ΔEdisp is 0.3 kcal/mol with a maximum deviation of 1.5 kcal/mol. On the other hand, the heteroatom11935



Article

mochemical values are shown in Figure 5. Only one value for twice hydrogenated 5-methylfurfural (compound 3 in Figure 5) exceeds the maximum deviation shown in Figure 3: the Hf,298 is underestimated by 6.8 kcal/mol. This is expected since five out of seven groups in this adsorbate are gas-phase-like and, as such, they are estimated from gas-phase group additivity. All other values fall well within the deviations observed for the fitting set.

based model results in a mean deviation of 12.0 kcal/mol and a maximum deviation of over 79.5 kcal/mol in Hf,298. These results are consistent with ring strain, intraadsorbate interactions, and other local effects playing a major role in adsorption properties. Figure 4b indicates that the heteroatom-based model is adequate for the dispersion contribution to the heat of adsorption with a relatively high R2 of 0.88 and a low mean deviation of 1.4 kcal/mol. Even in the worst case, ΔEdisp is overestimated by 8.3 kcal/mol because the interaction between the saturated part of the molecule and the surface is absent (see adsorbate “F12b” in the Supporting Information, section G). The heteroatom-based model’s good performance for ΔEdisp is indicative of the dispersion contribution to the adsorption energy from the heteroatom−surface interaction. The resulting heteroatom-based model values reported in Table 3 show that the more saturated the heteroatom, the higher its contribution to the dispersion energy.

■

CONCLUSIONS In this paper, we developed a group additivity scheme for the furanics family and applied it on Pd(111). We have regressed 101 furan-related adsorbates in order to obtain group additivity values. We have found that the model used for small openchain molecules is inaccurate for furanics; the maximum deviations from DFT-computed values were 27.6 kcal/mol in Hf,298, 17.6 cal/(mol K) in S298, and 5.9 cal/(mol K) in Cp,800. We have identified 17 new groups involving heteroatom− nearest neighbor interactions and introduced four corrections specific to furanic compounds that account for ring deformation and the level of hydrogenation of the ring carbons. These corrections result in much improved accuracy of the model, with the maximum deviation from DFT-computed values being 5.5 kcal/mol in Hf,298, 3.7 cal/(mol K) in S298, and 3.6 cal/(mol K) in Cp,800. Two corrections were determined as most important in error reduction: the “weak H−M” interactions, effectively adjusting the contribution from the saturated gas-phase-derived groups, and the “furan ring deformation”, accounting for the lack of sites for the desired sp3 hybridization of the dehydrogenated ring carbons. We have further introduced two models to estimate the dispersion contribution to the heat of adsorption. The full group additivity model showed good results for ΔEdisp with the mean deviation from DFT-computed values of 0.3 kcal/mol and the maximum deviation of 1.5 kcal/mol. The nineparameter heteroatom-based model resulted in the mean deviation of 1.4 kcal/mol and a maximum deviation of 8.3 kcal/mol and is suitable for a rough estimate of ΔEdisp but poor as a group additivity scheme. Finally, assessment of the model with compounds not in the training set demonstrates the accuracy of the proposed scheme.

Table 3. Heteroatom-Based Group Additivity Values for Dispersion Contribution to the Adsorption Energya group

ΔEdisp

#

Csat C(M) C(M)2 C(M)3 Osat Owk O(M) O(M)2 H(M)

−2.81 −3.40 −1.84 −1.43 −1.90 −5.54 −2.99 −2.58 −1.67

54 80 38 9 71 33 12 1 1

a

The units are kcal/mol. The number of adsorbates each group appears in, represented by “#”, is in the last column.

Assessing the Group Additivity Model. Next we assess the accuracy of the group additivity scheme on five representative adsorbates: two furanics, two hydrogenated furanics, and a ring-open compound. All of these compounds are expected to be typical intermediates in hydrogenation or ring-opening reactions on Pd(111). The resulting differences between group-additivity-estimated and DFT-estimated ther-

Figure 5. Thermochemistry differences between group-additivity-computed and DFT-computed values for a set of adsorbates typical to furanic chemistry not used in the fitting. 11936



Article

(9) Vorotnikov, V.; Mpourmpakis, G.; Vlachos, D. G. DFT study of furfural conversion to furan, furfuryl alcohol, and 2-methylfuran on Pd(111). ACS Catal. 2012, 2496−2504. (10) Wang, S.; Vorotnikov, V.; Vlachos, D. G. A DFT study of furan hydrogenation and ring opening on Pd(111). Green Chem. 2014, 16, 736−747. (11) Stampede. https://www.tacc.utexas.edu/stampede/. (12) Benson, S. W.; Buss, J. H. Additivity rules for the estimation of molecular properties. Thermodynamic properties. J. Chem. Phys. 1958, 29, 546−572. (13) Benson, S. W.; Cruickshank, F. R.; Golden, D. M.; Haugen, G. R.; O’Neal, H. E.; Rodgers, A. S.; Shaw, R.; Walsh, R. Additivity rules for the estimation of thermochemical properties. Chem. Rev. 1969, 69, 279−324. (14) Benson, S. Thermochemical Kinetics, 2nd ed.; John Wiley & Sons, Inc.: New York, 1976. (15) Salciccioli, M.; Chen, Y.; Vlachos, D. G. Density functional theory-derived group additivity and linear scaling methods for prediction of oxygenate stability on metal catalysts: Adsorption of open-ring alcohol and polyol dehydrogenation intermediates on Ptbased metals. J. Phys. Chem. C 2010, 114, 20155−20166. (16) Salciccioli, M.; Edie, S. M.; Vlachos, D. G. Adsorption of acid, ester, and ether functional groups on Pt: Fast prediction of thermochemical properties of adsorbed oxygenates via DFT-based group additivity methods. J. Phys. Chem. C 2012, 116, 1873−1886. (17) Kua, J.; Faglioni, F.; Goddard, W. A. Thermochemistry for hydrocarbon intermediates chemisorbed on metal surfaces: CHnm(CH3)m with n = 1, 2, 3 and m ≤ n on Pt, Ir, Os, Pd, Rh, and Ru. J. Am. Chem. Soc. 2000, 122, 2309−2321. (18) Cohen, N. Revised group additivity values for enthalpies of formation (at 298 K) of carbon–hydrogen and carbon–hydrogen– oxygen compounds. J. Phys. Chem. Ref. Data 1996, 25, 1411−1481. (19) Sumathi, R.; Green, J. W. H. Oxygenate, oxyalkyl and alkoxycarbonyl thermochemistry and rates for hydrogen abstraction from oxygenates. Phys. Chem. Chem. Phys. 2003, 5, 3402−3417. (20) Sabbe, M. K.; Saeys, M.; Reyniers, M.-F.; Marin, G. B.; Van Speybroeck, V.; Waroquier, M. Group additive values for the gas phase standard enthalpy of formation of hydrocarbons and hydrocarbon radicals. J. Phys. Chem. A 2005, 109, 7466−7480. (21) Sabbe, M. K.; De Vleeschouwer, F.; Reyniers, M.-F. O.; Waroquier, M.; Marin, G. B. First principles based group additive values for the gas phase standard entropy and heat capacity of hydrocarbons and hydrocarbon radicals. J. Phys. Chem. A 2008, 112, 12235−12251. (22) Khan, S. S.; Yu, X.; Wade, J. R.; Malmgren, R. D.; Broadbelt, L. J. Thermochemistry of radicals and molecules relevant to atmospheric chemistry: Determination of group additivity values using G3//B3LYP theory. J. Phys. Chem. A 2009, 113, 5176−5194. (23) Sebbar, N.; Bozzelli, J. W.; Bockhorn, H. Thermochemical properties, rotation barriers, bond energies, and group additivity for vinyl, phenyl, ethynyl, and allyl peroxides. J. Phys. Chem. A 2004, 108, 8353−8366. (24) da Silva, G.; Kim, C.-H.; Bozzelli, J. W. Thermodynamic properties (enthalpy, bond energy, entropy, and heat capacity) and internal rotor potentials of vinyl alcohol, methyl vinyl ether, and their corresponding radicals. J. Phys. Chem. A 2006, 110, 7925−7934. (25) Hudzik, J. M.; Bozzelli, J. W. Structure and thermochemical properties of 2-methoxyfuran, 3-methoxyfuran, and their carboncentered radicals using computational chemistry. J. Phys. Chem. A 2010, 114, 7984−7995. (26) Frisch, M. J.; et al. Gaussian 09, revision A.2; Gaussian Inc.: Wallingford, CT, 2009. (27) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-4 theory. J. Chem. Phys. 2007, 126, 084108. (28) Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 1996, 6, 15.

The largest deviations occur in higher-saturated molecules, likely due to approximating gas-phase-like surface groups with the vibrational and rotational gas-phase group contributions. The proposed scheme is expected to hold for other metal catalysts with metal-specific values. The potential applications of this work include the estimation of energetic preference of various conformers, e.g., furfural adsorbed through the carbonyl functional group vs the furanic functional group, and the quick estimation of thermochemistry for furanic intermediates involved in various C−H, C−O, and C−C scission reactions on Pd(111).

■

ASSOCIATED CONTENT

* Supporting Information S

Information on referencing the heat of formation for thermodynamic consistency, the equations for obtaining gasphase thermochemistry from statistical mechanics, the ab initio based gas-phase group additivity values, the details of counting the corrections in the group additivity scheme, the specifics of using group additivity tools, and structural information for each of the adsorbates. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: (302) 831-2830. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We gratefully acknowledge the support of Air Liquide. DFT calculations were performed using TeraGrid resources provided by Texas Advanced Computing Center (TACC) of the University of Texas at Austin.

■

REFERENCES

(1) Choudhary, V.; Sandler, S. I.; Vlachos, D. G. Conversion of xylose to furfural using Lewis and Brønsted acid catalysts in aqueous media. ACS Catal. 2012, 2, 2022−2028. (2) Choudhary, V.; Mushrif, S. H.; Ho, C.; Anderko, A.; Nikolakis, V.; Marinkovic, N. S.; Frenkel, A. I.; Sandler, S. I.; Vlachos, D. G. Insights into the interplay of Lewis and Brønsted acid catalysts in glucose and fructose conversion to 5-(hydroxymethyl)furfural and levulinic acid in aqueous media. J. Am. Chem. Soc. 2013, 135, 3997− 4006. (3) Mettler, M. S.; Mushrif, S. H.; Paulsen, A. D.; Javadekar, A. D.; Vlachos, D. G.; Dauenhauer, P. J. Revealing pyrolysis chemistry for biofuels production: Conversion of cellulose to furans and small oxygenates. Energy Environ. Sci. 2012, 5, 5414−5424. (4) Lange, J.-P.; van der Heide, E.; van Buijtenen, J.; Price, R. Furfurala promising platform for lignocellulosic biofuels. ChemSusChem 2012, 5, 150−166. (5) Sitthisa, S.; Resasco, D. Hydrodeoxygenation of furfural over supported metal catalysts: A comparative study of Cu, Pd and Ni. Catal. Lett. 2011, 141, 784−791. (6) Sitthisa, S.; Pham, T.; Prasomsri, T.; Sooknoi, T.; Mallinson, R. G.; Resasco, D. E. Conversion of furfural and 2-methylpentanal on Pd/ SiO2 and Pd-Cu/SiO2 catalysts. J. Catal. 2011, 280, 17−27. (7) Pang, S. H.; Medlin, J. W. Adsorption and reaction of furfural and furfuryl alcohol on Pd(111): Unique reaction pathways for multifunctional reagents. ACS Catal. 2011, 1272−1283. (8) Salciccioli, M.; Stamatakis, M.; Caratzoulas, S.; Vlachos, D. G. A review of multiscale modeling of metal-catalyzed reactions: Mechanism development for complexity and emergent behavior. Chem. Eng. Sci. 2011, 66, 4319−4355. 11937



Article

(29) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 1996, 54, 11169. (30) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865. (31) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (32) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 1994, 50, 17953. (33) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 1999, 59, 1758. (34) Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188−5192. (35) Blochl, P. E.; Jepsen, O.; Andersen, O. K. Improved tetrahedron method for Brillouin-zone integrations. Phys. Rev. B 1994, 49, 16223− 16233. (36) Murnaghan, F. D. The compressibility of media under extreme pressures. Proc. Natl. Acad. Sci. U. S. A. 1944, 30, 244−247. (37) Birch, F. Finite elastic strain of cubic crystals. Phys. Rev. 1947, 71, 809−824. (38) Lamber, R.; Wetjen, S.; Jaeger, N. I. Size dependence of the lattice parameter of small palladium particles. Phys. Rev. B 1995, 51, 10968−10971. (39) Pulay, P. Convergence acceleration of iterative sequencesthe case of SCF iteration. Chem. Phys. Lett. 1980, 73, 393−398. (40) Salciccioli, M.; Vlachos, D. G. Kinetic modeling of Pt catalyzed and computation-driven catalyst discovery for ethylene glycol decomposition. ACS Catal. 2011, 1, 1246−1256. (41) Burgess, D. Thermochemical data. In NIST Chemistry Webbook, NIST Standard Reference Database Number 69; Linstrom, P., Mallard, W., Eds.; National Institute of Standards and Technology: Gaithersburg, MD, 2009. (42) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. Assessment of Gaussian-2 and density functional theories for the computation of enthalpies of formation. J. Chem. Phys. 1997, 106, 1063−1079. (43) http://webbook.nist.gov/chemistry/. (44) Rangarajan, S.; Bhan, A.; Daoutidis, P. Language-oriented rulebased reaction network generation and analysis: Description of ring. Comput. Chem. Eng. 2012, 45, 114−123. (45) Rangarajan, S.; Bhan, A.; Daoutidis, P. Language-oriented rulebased reaction network generation and analysis: Applications of ring. Comput. Chem. Eng. 2012, 46, 141−152. (46) Antony, A.; Hakanoglu, C.; Asthagiri, A.; Weaver, J. F. Dispersion-corrected density functional theory calculations of the molecular binding of n-alkanes on Pd(111) and PdO(101). J. Chem. Phys. 2012, 136, 054702. (47) Bradley, M. K.; Robinson, J.; Woodruff, D. P. The structure and bonding of furan on Pd(111). Surf. Sci. 2010, 604, 920−925. (48) Belelli, P. G.; Ferullo, R. M.; Castellani, N. J. Unsaturated hydrocarbons adsorbed on low coordinated Pd surface: A periodic DFT study. Surf. Sci. 2010, 604, 386−395.

11938


Group Additivity for Estimating Thermochemical Properties of Furanic

Recommend Documents