Emilija
Kohls
and
Matthias
Stein
*
Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstr. 1, 39106 Magdeburg, Germany. E-mail: matthias.stein@mpi-magdeburg.mpg.de
First published on 28th June 2017
Rhodium-catalysed hydroformylation is the major industrial process to obtain aldehyde products from olefins. The isomerization of the olefin at the catalyst is the most prominent side reaction and lowers the overall yield of the process. We here investigate whether the olefin isomer distribution obtained from batch experiments using Rh(BiPhePhos) as a catalyst and n-decene as the olefin reaches the thermodynamic equilibrium and computational quantum chemical approaches are able to accurately reproduce the experimental olefin isomer distribution. The relative energies of cis/trans configurational and double bond positional isomers of long chain n-decene were calculated using Hartree–Fock, DFT and correlated ab initio methods. Results were compared to experimental data. Electron correlation was found to be critical for the description of cis-isomer relative energies. Dispersion corrections in the DFT calculations partially compensate for deficiencies and generally improve the agreement with experiment. Adding thermodynamic corrections suffers from the neglect of certain contributions to the entropy of flexible molecules. Accounting for the entropy of mixing of multiple conformers significantly reduces the deviation from experiment. The equilibrium distribution of long chain olefins is reasonably described by correlated QM and also some density functional methods. Computational thermochemistry has thus reached a state where it provides reliable parameters for complex reaction network models and process engineering.
Apart from the desired linear aldehyde product, double bond isomerization of the olefin coordinated to the catalyst occurs as a side reaction by which undesired branched aldehydes are being produced.3–6 This lowers the overall yield of the preferred main linear product. HRh(CO)(BiPhePhos) is an advanced hydroformylation (see Fig. 1) catalyst but also known to be very active in isomerization of the olefin double bond. By choice of appropriate conditions, the formation of undesired side products can be suppressed and thus an increase of yield or selectivity can be obtained.7
Fig. 1 The hydroformylation reaction of 1-decene to the terminal undecanal as the desired product. Syngas (CO/H2) is used at 378 K, 20 bar in presence of a Rh(BiPhePhos) catalyst. |
Long chain olefins from sustainable resources are a versatile alternative to substrates from petrochemicals.8,9 Their use poses challenges to an appropriate process design in terms of the solubility of substrates and the separation of products and catalyst.10,11 Recently, Jörke et al. have conducted an isothermal batch isomerization experiment using Rh(BiPhePhos)12 as catalyst and followed the concentration profiles of olefin isomers.13,14 1-Decene was used as the starting material and was isomerizing into all internal cis- and trans-decene isomers in presence of the catalyst and absence of syngas. This was assumed to be a stepwise process (Fig. 2). All n-decene positional isomers could be detected, assigned and quantified by gas chromatography. The time-resolved concentration profiles of each isomer were followed until no further changes in the concentration profiles were observed. The experimentally measured composition of isomers (Fig. 3) was used to deduce Gibbs free energy differences of internal n-decene isomers relative to the terminal 1-decene (values given in second column in Table 1).
Fig. 2 The olefin isomerization is the major side reaction to the hydroformylation. A Rh(BiPhePhos) catalyst in absence of syngas was used to investigate the 1-decene isomerization reaction. |
Fig. 3 Experimentally measured equilibrium composition of all positional and geometrical linear decene isomers. (a) 1-Decene and internal cis-decene isomers; (b) 1-decene and internal trans-decene isomers. Reproduced with permission from the publisher from ref. 14. |
Exp14 | HF | BP86 | BP86-D3 | B3LYP | B3LYP-D3 | M06-2X | B2PLYPa | B2PLYP | MP2 | LPNO-CCSDa | |
---|---|---|---|---|---|---|---|---|---|---|---|
a Single point calculations at MP2/def2-TZVP optimized structures. | |||||||||||
1d | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.00 | 0 | 0 |
t2 | −10.07 | −9.44 | −13.94 | −11.66 | −12.46 | 10.51 | −9.80 | −10.48 | −11.11 | −10.65 | −9.22 |
t3 | −9.98 | −9.06 | −13.35 | −10.84 | −12.01 | −9.90 | −8.91 | −10.12 | −11.28 | −10.75 | −8.31 |
t4 | −10.75 | −9.40 | −13.98 | −11.41 | −12.56 | −10.37 | −9.33 | −10.46 | −11.63 | −11.24 | −9.17 |
t5 | −8.56 | −9.40 | −13.96 | −11.39 | −12.54 | −10.34 | −9.48 | −10.46 | −11.62 | −11.28 | −9.09 |
c2 | −7.28 | −2.27 | −8.64 | −7.68 | −6.57 | −5.84 | −4.93 | −5.40 | −6.48 | −5.83 | −4.75 |
c3 | −5.99 | −1.27 | −7.47 | −6.68 | −5.52 | −5.01 | −3.53 | −4.66 | −6.18 | −5.92 | −3.86 |
c4 | −6.87 | −1.62 | −8.02 | −7.59 | −5.97 | −5.72 | −3.83 | −5.02 | −7.02 | −7.00 | −4.41 |
c5 | −4.78 | −1.62 | −8.07 | −7.66 | −6.00 | −5.82 | −4.26 | −5.57 | −7.77 | −7.12 | −4.43 |
MUE | 2.73 | 2.89 | 1.33 | 1.69 | 0.91 | 1.51 | 1.07 | 1.30 | 1.07 | 1.51 |
However, it could not finally be decided whether the final distribution of n-decene isomers was fully corresponding to the equilibrium distribution.
Computational aspects of modeling the hydroformylation process have mainly focused on mechanistic investigations (for recent reviews see ref. 15–17) and kinetics18 of the main reaction but not on the side reactions.
We here aim at a quantitative description of the distribution of positional and configurational n-decene isomers using quantum chemical (QC) thermochemistry.
Knowledge of thermochemical stability of chemical species is the key to industrial applications of structure-based calculations.19 Thermochemical predictions are of crucial importance for the design and control of industrial chemical processes. There are accurate composite approaches which include the Gn methodologies20 and the Wn theory.21,22 On the other hand, quantum chemical calculations, in particular Density Functional Theory (DFT), become increasingly important in chemical industry when determining thermochemical properties of larger systems. Insight gained at the molecular level and molecular industrial thermodynamics have recently been reviewed by Deglmann et al.19 They give an overview about current state-of-the-art applications of quantum chemical calculations to problems of engineering relevance. They use the calculation of chemical reaction thermodynamics as an example to demonstrate the accuracy of mostly DFT approaches. Nevertheless, Koglin and Meijer have very early shown that DFT methods are not sufficient to obtain a quantitatively correct description of n-hexane ggg and gtg conformers which needed to be carefully benchmarked against higher level calculations.23 There is large number of publications that demonstrate that B3LYP and other popular exchange–correlation functionals tend to have significantly larger errors for larger structures (in particular hydrocarbons).24–26 The failure is assigned to the inability of most functionals to describe long-range van der Waals effects and the neglect of medium-range electron correlation.27 The efforts to overcome these shortcomings resulted in the development of highly parameterized, dispersion-corrected, as well as double hybrid functionals. M05-2X functional performed satisfactorily for main-group thermochemistry and barrier heights due to its improved description of medium-range correlation.28,29 On the other hand, the increase in error with increasing system size for hydrocarbons using B3LYP, depended on the presence of particular structural and stereoelectronic features.29
The quantum thermochemistry of long-chain n-alkanes has received attention as model system for uncoupled internal rotations,30 one-dimensional hindered rotors31 and also to assess stereoelectronic effects32 and isodesmic bond separation energies.33 There is, however, much less work dedicated to the thermochemistry of olefins.
We here use DFT and higher correlated tools such as MP2 and coupled-cluster to calculate the distribution of highly flexible unbranched decene isomers and compare the results to those obtained from experiment. We also address the problem of scaling of vibrational frequencies, solvent effects, as well as the problem of entropic contributions and the related issue of consideration of multiple conformations.
For generating a library of conformers for each of the 9 isomers the stochastic Monte Carlo multiple minimum method (MCMM)51 was used in MacroModel 10.6.52 The maximum number of steps was set to 8000 and number of steps per rotatable bond to 2000. The energy window for saving the structures was set to 10 kJ mol−1. In order to determine the molecular mechanics force field best parameterized for the decene molecules we have also performed test calculations for different force fields and checked the quality of their stretch, bend and torsion parameters (not shown). We find that OPLS_200553 is the best parameterized for this case and it was used in the subsequent conformational search. Each generated conformer was re-optimized at B3LYP-D3/TZVP level of theory using Jaguar 8.6.54
(1) |
Physicochemical data for long chain n-olefins and kinetics regarding the isomerization side reaction are hardly available in the literature but required for a qualitative and quantitative understanding and control of the process.
The Benson group contribution method is an additive empirical approach to describe relative energies of molecules from functional group energy terms fitted to experiments. There are terms for sp3–sp2 carbon–carbon single bonds (one in 1-decene and two in other n-decene isomers), sp2–sp2 carbon–carbon double bonds, sp3–sp3 carbon–carbon single bonds and so one. The Benson group contribution method, however, cannot discriminate between subtle differences of structurally complex molecules such as 3- or 4-decene isomers and conformers. The Benson group contribution method gives a −4.9 kJ mol−1 energy difference between trans olefins and their cis isomers as a result of van der Waals repulsions which agrees well with 4.5 kJ mol−1 (1 kcal mol−1) which can be found in standard organic chemistry textbooks.55 These small energy differences point out the challenges to quantum chemical methods to accurately describe the thermodynamics and relative energies of long chain olefins.
The calculated relative energies (ΔE) are summarized in Table 1 and Fig. 5.
Fig. 5 Performance of a subset of computational methods: mean unsigned errors for relative energies of all cis- and trans- n-decene isomers. Results for additional computational methods can be found in the ESI† (Table A). |
From Table 1 it becomes obvious that all computational methods (including DFT) give 1-decene to be the least stable isomer and the trans configurations to be lower in energy than their cis counterparts by 4–6.5 kJ mol−1. The experimentally determined energy differences between a cis isomer and its trans counterpart are below 4 kJ mol−1. This difference is best reproduced by the dispersion corrected functionals (BP86-D3 and B3LYP-D3). The mean unsigned error (MUE) decreases from 2.89 to 1.33 and from 1.69 to 0.91 kJ mol−1 for BP86(+D3) and B3LYP(+D3) respectively when dispersion corrections are added. This shows that dispersive corrections at least partially correct for the absence of long range dispersion in these density functionals.
From the DFs that were not corrected for vdW interactions, M06-2X with 54% of non-local, Hartree–Fock exchange performs best with a MUE of 1.51 kJ mol−1 followed by B3LYP and BP86. Since M06-2X was parameterized to account for the interpair but nonlocal electron correlations, it confirms that the medium range correlation is the missing part in the DFs for correctly describing hydrocarbon isomerization energies. The M06-2X functional, generated to overcome the correlation problem, is superior to B3LYP and other hybrid functionals but still does not provide complete correction of electron correlation. Song et al. conclude that local correlation needs to be incorporated and current exchange functionals are lacking this contribution.56 Double hybrid functionals such as B2PLYP were introduced to address the medium-range correlation issue and perform well for relative olefin double bond isomer energy differences (with a MUE of 1.30 kJ mol−1) and comparable to MP2 for single point calculations. Surprisingly, coupled cluster LPNO-CCSD calculations yield a larger deviation from experiment (1.51 kJ mol−1 on average) than MP2 (1.07 kJ mol−1) and this has to be assigned to the slower convergence with basis set size for these correlated wavefunction calculations which was also shown for the thermochemistry of the hydroformylation reaction.57
For the cis-configuration, all methods (except MP2 and B2PLYP) give cis-2-decene to be the most stable which agrees with the experimental findings. For trans-isomers, most of the methods predict either trans-2-decene or trans-4-decene (in agreement with experiment) to be the most stable one.
Hartree–Fock calculations perform well for the trans but not so for the cis isomers. This deficiency can be accounted for by inclusion of electron correlation (see MP2 values in Fig. 5). This level of treatment is necessary to correctly describe 1,4-alkyl–alkyl repulsive interactions in agreement with previous results for alkanes27 and relative stabilities of unbranched vs. branched alkene isomers.58Fig. 5 also shows that M06-2X with 54% of exact exchange exhibits larger error for the cis isomers than for the trans. This deficiency in cis stabilization arises from electron correlation and is not an exchange problem. The opposite is observed for B3LYP. It shows larger deviation for the trans isomers but this is significantly reduced by inclusion of dispersion interactions. MP2 is known to overestimate correlation effects in some chemical systems on the medium electron length scale but here we do not observe such an effect (spin-scaled MP2 yields identical results and therefore not given in Table 1).
It has to be noted that the experimental energy difference between some of the isomers is less than 1 kJ mol−1 and below the limits of computational accuracy. In several cases the difference in the calculated values is so minute that a couple of isomers can be considered to be isoenergetic and a strict energy ranking becomes impossible. According to experiment, the internal isomers (t5 and c5) are the least stable n-decene isomers. None of the computational methods, however, reproduce this finding from an electronic energy difference and thermal and entropic corrections to give Gibbs free energy differences must be accounted for (see below). Further, deviations from experimentally determined relative decene isomer stabilities may result from the neglect of solvent effects and the accessibility of multiple conformations for large, flexible long-chain molecules.33
B3LYP/def2-TZVP | ΔE (ε = 1) | ΔE (ε = 2) | ΔE (ε = 37) |
---|---|---|---|
1d | 0 | 0 | 0 |
t2 | −12.46 | −12.13 | −11.61 |
t3 | −12.01 | −11.63 | −11.03 |
t4 | −12.56 | −12.20 | −11.65 |
t5 | −12.54 | −12.19 | −11.64 |
c2 | −6.57 | −6.41 | −6.16 |
c3 | −5.52 | −5.31 | −4.99 |
c4 | −5.97 | −5.78 | −5.50 |
c5 | −6.00 | −5.83 | −5.57 |
MUE | 1.69 | 1.56 | 1.36 |
Fig. 6 Solvent effects on calculated energy differences between 1-decene and the other n-decene isomers. The deviation of B3LYP results from experiment is given. |
Inclusion of solvent effects consistently lowers the energy differences between 1-decene and the other isomers and brings them closer to experiment. The effect of solvent stabilization, however, is small and below 1 kJ mol−1 in magnitude. All isomers are almost equally stabilized by the solvent effects which can be seen in Fig. 6.
The average deviation from experiment decreases from 1.69 to 1.56 in decane and to 1.36 kJ mol−1 in DMF in comparison to gas phase B3LYP data. The effect for a decane/DMF solvent system (80:20 w/w) will be intermediate. Since COSMO is an implicit solvent model specific solvent–solute interactions are not taken into account and the solution is ideal. We have previously estimated the non-ideality of a decane/DMF solvent system and their effect on the reaction enthalpy of the hydroformylation of dodecene.57 Significant activity αi and fugacity coefficients φi were needed to rationalize the equilibrium concentrations of reactants and products. For energy differences, however, these effects are assumed to cancel out.
When vibrational frequencies are obtained with sufficient accuracy, the calculated zero point vibrational energies (ZPVE) are also reliable. On the other hand Gibbs free energies i.e. entropies cannot be calculated with the same accuracy. This has its origins in approximating vibrations as being harmonic in nature. The vibrations that are most likely to have anharmonic character are the low-frequency vibrations which contribute most to the entropy. An error of 5% in a low frequency vibration causes an error of 15–30% in entropy.59 We here
(i) use empirical scaling factors from ref. 60 to partially compensate for anharmonic effects and deficiencies of the electronic structure calculation method; and
(ii) consider the accessibility of multiple conformers.
Exp.14 | MP2a | B3LYPc | M06-2Xb | BP86d | |
---|---|---|---|---|---|
a λ = 0.9946. b λ = 0.9871. c λ = 1.0044. d λ = 1.0337. | |||||
1d | 0 | 0 | 0 | 0 | 0 |
t2 | −10.07 | −12.05 (−12.04) | −16.22 (−16.23) | −16.62 (−16.61) | −18.74 (−18.79) |
t3 | −9.98 | −11.40 (−11.40) | −15.43 (−15.43) | −14.60 (−14.58) | −17.66 (−17.69) |
t4 | −10.75 | −12.55 (−12.54) | −12.65 (−12.64) | −17.36 (−17.35) | −13.96 (−13.98) |
t5 | −8.56 | −10.46 (−10.45) | −11.31 (−11.31) | −9.05 (−9.04) | −12.95 (−12.96) |
c2 | −7.28 | −9.37 (−9.37) | −12.53 (−12.53) | −13.59 (−13.59) | −16.47 (−16.51) |
c3 | −5.99 | −9.59 (−9.59) | −8.03 (−8.02) | −8.58 (−8.57) | −10.62 (−10.64) |
c4 | −6.87 | −9.60 (−9.60) | −10.35 (−10.34) | −8.37 (−8.36) | −12.61 (−12.63) |
c5 | −4.78 | −8.16 (−8.16) | −6.75 (−6.74) | −2.62 (−2.62) | −8.81 (−8.81) |
MUE | 2.36 (2.36) | 3.62 (3.62) | 3.85 (3.84) | 5.94 (5.97) |
Enthalpies are in general less affected by the number of conformers which will be very close in energy (and thermodynamic corrections can be assumed to be identical for all conformers). Gibbs free energies, however, are much more affected by the presence of multiple conformers since they include entropic contributions. The difference in entropy between each conformer is negligible (ΔS for conformational change can be assumed to be zero; at most 2.8 J K−1 mol−161), but the entropy of mixing multiple conformers is the critical contribution. So the total molar entropy, So, is the sum of entropies of each conformer plus a correction for the mixing term (ΔSmix) (eqn (2)):
(2) |
(3) |
(4) |
ΔSmix ≈ Rln(n) | (5) |
The importance of accounting for the entropy of mixing was pointed out amongst others by Block61et al. for small alcohols radicals. Guthrie62 describes in great detail the difficulties and approaches to properly account for entropy contributions for flexible molecules. Decene molecules are highly flexible with 8 rotatable bonds and the number of possible conformers that contribute to the energy at finite temperature is large. In the entropy of mixing, only distinguishable conformations have to be considered. Rotations about bonds to symmetrical groups (i.e. CH3, CF3, C(CH3)3) do not contribute to the number of conformers, while rotations about internal bonds or bonds to unsymmetrical terminal groups do.
We have conducted an exhaustive MCMM search to identify the number of possible conformers by torsional sampling. All obtained structures were subsequently re-optimized at the B3LYP-D3/TZVP level of theory. B3LYP-D3 was the method of choice for this analysis as it showed the best performance from the DFs in terms of energy and Gibbs free energy differences. After quantum-chemical re-optimization the number of conformers in an energy window of 2RT at: 298.15 K and 378 K (that is 5 kJ mol−1 and 6.3 kJ mol−1, respectively) were counted. The number of the conformers for each of the 9 positional and configurational isomers in an energy window of 6.3 kJ mol−1 is given in Table 4. The data for 298.15 K can be found in the ESI† (Table D).
n | S (single conformer) | S corr | ΔG | Exp.14 | |
---|---|---|---|---|---|
1d | 110 | 541.0 | 580.1 | 0.00 | 0.00 |
t2 | 72 | 547.3 | 582.9 | −12.80 | −10.07 |
t3 | 102 | 548.4 | 586.9 | −13.22 | −9.98 |
t4 | 66 | 540.6 | 575.4 | −8.76 | −10.75 |
t5 | 44 | 537.0 | 568.5 | −6.25 | −8.56 |
c2 | 64 | 555.0 | 589.6 | −9.98 | −7.28 |
c3 | 37 | 555.8 | 585.9 | −7.26 | −5.99 |
c4 | 40 | 553.0 | 583.7 | −7.17 | −6.87 |
c5 | 31 | 542.0 | 570.6 | −2.08 | −4.78 |
MUE | 2.16 |
Table 4 shows that consideration of entropy of mixing of conformers improves the agreement between calculated and observed values. It reduces the average deviation from experiment from 2.79 kJ mol−1 to 2.16 kJ mol−1 in terms of Gibbs free energy differences.
The optimized conformers in an energy window of 10 kJ mol−1 were used for performing a Boltzmann-weighted energy distribution. The mole fraction of the ith conformer of n conformers in equilibrium is weighted by its probability (eqn (6)):
(6) |
The Boltzmann-weighted energy for each isomer is then obtained:
(7) |
For n-decene isomers, Boltzmann weighting did not improve the calculated relative energies (ΔE) compared to the single geometry calculations at the same level of theory (B3LYP-D3/TZVP) in Table 1.
This indicates that the very small energies of conformational change (often less than 0.1 kJ mol−1) are beyond the accuracy of current QM methods. Then, the calculated mole fractions of the conformers Boltzmann-weighted conformational energies do not reduce the deviation from experiment.
Fig. 7 Calculated binding of 1-decene to a HRh(CO)(BiPhePhos) catalyst. The steric demand by the biphenylbisoxy ligand parts poses a steric restriction to coordination of internal n-decene isomers. |
The relative stabilities of n-decene isomers at thermodynamic equilibrium can qualitatively be reproduced by QM methods to rank trans- vs. cis-configurations lower in energy by about 4 kJ mol−1. The accuracy of the computational approaches to generate reliable thermodynamic data for industrial process engineering is not fully exploited yet and requires careful inspection.19 In previous studies for other alkane and alkene isomers27,58 DFT calculations were not able to reproduce the energetic ordering of the isomers. When comparing to experiment, the MUE for all methods is below 3 kJ mol−1. Comparison of the Hartree–Fock and correlated MP2 and CCSD methods establishes the importance of electron correlation in the description of the cis isomers, but less for the trans isomers. One possible explanation is the missing correlation in the description of 1,4-alkyl–alkyl repulsive interactions. According to our study, van der Waals interactions also play a critical role and empirically correcting them significantly improves the results for the DFs. The final relative thermodynamic stabilization of trans-4- versus trans-2-decene requires a careful inclusion of electron correlation from MP2 calculations (which is also included in the B2PLYP functional).
The molecular thermochemistry of long-chain olefins still poses a challenge to current computational approaches. The RRHO approximation introduces errors in the entropic part which questions the ability of standard simplifications to the state function Q to calculate thermodynamic properties with high accuracy. An agreement with experiment can be significantly improved by correcting for the entropy of mixing for flexible molecules. This benchmark of thermochemical data of long chain olefins will help to assign parameters required in the kinetic network model for the hydroformylation process.63,64
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7nj01396e |
This journal is © The Royal Society of Chemistry and the Centre National de la Recherche Scientifique 2017 |