Isomerization and Decomposition of 2-Methylfuran with

Article Cite This: J. Chem. Inf. Model. 2019, 59, 3454−3463

pubs.acs.org/jcim

Isomerization and Decomposition of 2‑Methylfuran with External Forces Agnieszka Brzyska† and Krzysztof Wolinś ki*,‡ †

Jerzy Haber Institute of Catalysis and Surface Chemistry, Polish Academy of Sciences, Niezapominajek 8, 30-239 Krakow, Poland Department of Theoretical Chemistry, Maria Curie-Sklodowska University, pl. Maria Curie-Sklodowska 3, 20-031 Lublin, Poland

‡

Downloaded via NOTTINGHAM TRENT UNIV on September 5, 2019 at 15:07:43 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: The primary goal of this project was to evaluate the performance of the Standard and Enforced Geometry Optimization (SEGO) method which we have recently developed. The SEGO method has been designed for an automatic location of multiple minima on the molecular Potential Energy Surface (PES), and its usefulness has been demonstrated so far for three molecules only. In this project we applied the SEGO method to explore the 2-methylfuran (2MF) PES. Our choice was not accidental: this molecule recently gained a great deal of interest as a potential candidate for biofuel, and therefore its pyrolysis is extensively studied. To understand pyrolysis of 2MF a detailed knowledge about its PES is needed. In these studies we explored the 2MF PES and located a surprisingly large number of local minima corresponding to 2MF isomers and decomposition products. Some of the 2MF isomers and fragments have amazing structures which most likely were previously unknown. All structures presented in this paper were found in an automatic manner as a result of the enforced chemical reactions driven by the SEGO method. Thus, in these studies we have not only proven that the SEGO method is a very efficient tool for exploring molecular potential energy surfaces, but we also obtained very detailed knowledge about the 2MF PES.

■

methylfuran.5,9 This paper is organized as follows. In the next section a short theory of the SEGO method will be presented. Then the computational details will be given. They will be followed by the presentation of the structures found in several SEGO searches performed with external forces of a different magnitude. In the last section we will present conclusions arising from our studies.

INTRODUCTION Over the past decade there has been a growing interest in biofuels which can replace less economical and environmentally more dangerous and harmful sources of energy like coal and oil. Among others, 2-methylfuran has recently received a great deal of attention as a possible candidate for such a biofuel.1,2 For this reason the pyrolysis of 2-methylfuran has been extensively studied. The main goal of those experimental and theoretical studies3−9 was to find possible products and intermediates as well as to determine major reaction paths. This simply means exploration of the molecular Potential Energy Surface (PES) of 2-methylfuran in order to locate stationary points there. In general it is not a trivial task, but there is now an efficient tool available to do this in essentially a routine manner. Very recently we have proposed the new Standard and Enforced Geometry Optimization (SEGO) method10 which allows automatic location of multiple minima on PES. This method is based on the use of external forces which move us from one local minimum toward another and cross the separating energy barrier. Preliminary results obtained for butene, dioxane, and pyranose ring have proven the usefulness of the SEGO method. Thus, having a quite efficient tool to explore any molecular PES we applied it in these studies to the potential energy surface of 2-methylfuran. We expected to find some exotic species related to carbenes which play an important role in the decomposition of 2© 2019 American Chemical Society

■

THE SEGO METHOD The SEGO method arose from the earlier proposed Enforced Geometry Optimization (EGO) method11,12 in which geometry optimization is performed in the presence of external forces applied to selected nuclei in a molecule. The EGO method has been originally introduced to simulate Atomic Force Microscopy (AFM) experiments, and it turned out to be very successful13−15 in those applications. It was also successfully applied to drive isomerization reactions of stilbene and azobenzene.16,17 It should be noted that the EGO method is closely related to three other methods introduced in the field around the same time.18−20 In the SEGO method we begin the search for multiple local minima on PES from the first given minimum obtained in the Standard Geometry Optimization (SGO) step. The external forces are applied to all atoms in a Received: April 25, 2019 Published: June 25, 2019 3454

DOI: 10.1021/acs.jcim.9b00352 J. Chem. Inf. Model. 2019, 59, 3454−3463

Article

Journal of Chemical Information and Modeling

of the maximum force allowed, the SEGO search was executed to locate five local minima. The number of requested local minima to be found is given to the program through the input deck. Density functional theory (DFT)21,22 is the most commonly applied method23 in the routine calculations for medium- and large-sized molecular systems. In particular, the B3LYP24 exchange-correlation functional has made great success in predicting the ground-state electronic structures, molecular geometries, and reaction energetics.25 That is why in this work we used mainly the DFT/B3LYP approach with the 6-31G(d,p) basis set.26,27 Additional calculations were also performed with the 6-31+G(d,p)26−28 and 6-31++G(d,p) basis sets.27,29 However, the conventional correlation MP2 method in general can provide somewhat better results. Therefore, in order to verify stability of some suspicious structures, we also ran calculations at the MP230,31 level. In all geometry optimizations our standard default convergence criteria was used (i.e., the maximum residual force on any nucleus must be less than 0.0003 au and either the maximum atomic displacement must be less than 0.0003 au or the energy change from the previous cycle must be less than 0.000001 hartree). All calculations were performed in a parallel mode with the PQS program package,32 and structural changes were visualized and analyzed using our Graphical User Interface, PQSMol.32 Since in these studies we are interested in the ordinary potential energy surface of 2MF, therefore all energy values reported here refer to PES not to FMPES which means that the extra term in eq 2 is not included.

molecule, and the Enforced Geometry Optimization (EGO) is performed. Once this EGO step is done, external forces are switched off, and a stressed molecule is relaxed in the SGO step which hopefully locates the next local minimum on PES. For this structure new external forces are generated and applied to nuclei in the following EGO step. In order to locate, let us say, L local minima on PES, this procedure has to be repeated L times. Starting from the first minimum found in ordinary SGO, two consecutive geometry optimizations, EGO and SGO, are required to locate the next (second) local minimum. Thus, in the whole procedure a sequence of SGO, EGO, SGO, ... steps is performed, and therefore the method has been named SEGO. The crucial issue in SEGO is the choice of external forces. While in the EGO method constant forces (in magnitude) are applied to selected atoms, in the SEGO method variable forces act on all atoms. External forces in SEGO are uniquely defined for a given local minimum from which the EGO step begins. They are derived from the energy Hessian H(r0) calculated at the current minimum r0. At an arbitrary geometry r these external forces are given as f(r) = −s ∗H(r0) ∗(r − rref )

where the vectors f and r collect forces and coordinates of all N atoms in a molecule f = (f1, f2, ..., fN) and r = (r1, r2, ..., rN), respectively, with rk = (xk, yk, zk) and similarly for fk for a nucleus k. The reference geometry vector rref = (rcom, rcom, ..., rcom) refers to the center of mass with rcom = (xcom, ycom, zcom). Since in our case a molecule is always reoriented in such a way that its center of mass is placed at the origin, the reference vector is always zero, and rref will be omitted in all formulas. Thus, external forces used in SEGO can be written as f(r) = −s ∗H(r0) ∗r

■

RESULTS AND DISCUSSION In these studies we focus on two issues. The first one is to locate as many local minima on the potential energy surface of 2MF as possible because they are important for understanding the pyrolysis of 2MF. The second goal is to find possible channels and products of the decomposition of 2MF. It should be emphasized that in the present studies we do not propose any models of the 2MF pyrolysis and do not search for any transition states. We are just interested in finding local minima on the potential energy surface of 2-methylfuran, i.e., its isomers and decomposition products. There is one important issue which we also do not consider here. This is the proton transfer in the 2MF ring which leads to various cyclic carbenes, important intermediates in the 2MF pyrolysis. The proton transfer reactions in 2MF will be the subject of a separate paper which will be ready very soon. Originally we considered doing all calculations with the 6-31G(d,p) basis set only. If one is interested in finding the equilibrium geometry of an ordinary molecule, then this double-ζ polarized basis set should be good enough especially for qualitative needs such as in our project. In such a case, diffuse functions should not make much difference. However, at one point we wanted to compare our results obtained for some carbenes with those reported recently by Somers et al.9 who used the 6-31+G(d,p) basis set. So, we repeated some of our calculations with this basis set; and it turned out that indeed diffuse functions do not make much difference in optimized geometries of those carbenes, but they have significant influence on the external forces generated according to eq 1. Therefore, we performed some additional SEGO searches with the basis sets containing diffuse functions, i.e., 6-31+G(d,p) and 6-31++G(d,p). Using external forces of a different magnitude and also allowing different maximum geometry displacement during geometry optimization, we performed in this work over 40 SEGO

(1)

where the scalar parameter 1 > s > 0 limits in practice maximum external forces on atoms as given by H(r0)*r s = (maximum force allowed on any atom) /(maximum H(r0)r0 force on any atom)

where a maximum force allowed is defined by the user. External forces (1) interact with a molecule, and the total energy of a system is given as Etotal(r) = E(r) + 1/2∗s ∗r T∗H(r0) ∗r

(2)

where E(r) is the energy of a molecule at the current geometry r without external forces, i.e., it defines the ordinary PES, while Etotal(r) defines the so-called19 Force-Modified PES (FMPES). The energy formula 2 is consistent with the definition of external forces (1), i.e., the negative first derivative of Etotal gives the external forces f (1). It should be emphasized that according to (1) the external force on a given atom depends on the position of ALL atoms, and, in general, it changes in every geometry optimization cycle

■

COMPUTATIONAL PROCEDURE All our SEGO calculations commenced from the 2-methylfuran (2MF) optimized ground state geometry obtained by the standard geometry optimization procedure. Maximum external forces on one atom were limited by the scalar parameter s in (1). We performed calculations with the maximum external force allowed on any atom varying from 0.100 to 1.000 au, usually in steps of 0.050 au (1 au = 82400 pN). For each value 3455


Article


Figure 1. a. Geometry optimization history for the selected external force f = 0.200 au; the example of the isomerization channel for 2MF. b. Example of the SEGO decomposition channel for 2MF. c. Example of the SEGO isomerization/decomposition channel for 2MF.

searches with the 6-31G(d,p) basis set, over 20 with the 631+G(d,p) basis set, and over 20 with the 6-31++G(d,p) basis

set. Each SEGO search started from the 2MF equilibrium geometry found in the standard geometry optimization (SGO) 3456


Article


Figure 2. a. Unique molecular systems obtained from 2MF by the enforced intramolecular rearrangement at the DFT/B3LYP/6-31G(d,p) level. b. Unique molecular systems obtained from 2MF by the enforced intramolecular rearrangement at the DFT/B3LYP/6-31+G(d,p) level. c. Unique molecular systems obtained from 2MF by the enforced intramolecular rearrangement at the DFT/B3LYP/6-31++G(d,p) level.

step. In every SEGO search we requested five local minima to be found on ordinary PES. The number of requested minima

(five) is quite arbitrary, but if one asks for more, then the output files get very long which is not very handy. The main 3457


Article


Scheme 1. a. Unique Isomerization Channels Obtained at the DFT/B3LYP Level with Different Basis Sets,a b. Unique Decomposition Channels Obtained at the DFT/B3LYP Level with Different Basis Sets,a and c. Unique Decomposition/ Isomerization Channels at the DFT/B3LYP Level with Different Basis Setsa

Carbene marked in bold, * rotamer.

a

SEGO search represents some kind of the “enforced reaction path” triggered by external forces. In our studies, in each SEGO search we want to find five local minima so each SEGO search represents the five-stage “enforced reaction”. Three types of such a path appeared in our studies. They represent channels of the following: isomerization, decomposition, and

output from the SEGO calculations is contained in the socalled geometry optimization history plot which shows the total energy changes with geometry optimization cycles. On such a plot all local minima found in the SEGO search can be seen as well as energy barriers separating two consecutive minima. Thus, the geometry optimization history plot for each 3458


Article

Journal of Chemical Information and Modeling mixed isomerization-decomposition. In the first case the reaction goes from a one 2MF isomer to another, so all local minima found in SEGO correspond to the isomers of 2MF. In the second case we deal with a sequence of decompositions of 2MF into fragments, i.e., each energy minimum on a SEGO plot represents a molecular system which consists of two or three fragments. The third case is a mixture of the previous two cases. Typical examples of these SEGO searches are shown in Figure 1a−c, respectively. We will discuss these figures later in separate sections. At the B3LYP/6-31G(d,p) level we located overall more than 150 structures on the 2MF PES. Fifty structures found in the first ten SEGO searches with external forces varying from 0.100 au to 0.650 au (with a step of 0.050 au) are shown in the Supporting Information in Figure S1. One can distinguish there local minima corresponding to the isomers of 2MF and products of the decompositions of 2MF into two or three fragments. The vast majority of these structures obtained in SEGO searches are the same. Among these we found 16 unique isomers and 22 sets of the decomposition products of 2MF. At the B3LYP/6-31+G(d,p) and also at the B3LYP/6-31+ +G(d,p) levels, we performed over 20 SEGO searches with maximum external forces varied from 0.100 to 0.900 au with a step of 0.100 au. In the same cases the structures found were the same as before with the 6-31G(d,p) basis set. However, in many cases we got new structures, some of them completely unexpected. We found 27 new unique 2MF isomers with the 631+G(d,p) basis set and 18 more with the 6-31++G(d,p) basis set. With these bigger basis sets we also found 20 more sets of decomposition of products. Thus, overall we found 61 isomers and 42 decomposition products of 2-methylfuran. Isomerization Channels. Sixteen unique isomers of 2MF obtained with the 6-31G(d,p) basis set are presented in Figure 2a. Out of 16 isomers of 2MF found, ten of them are carbenes, but only three of them are cyclic five-membered rings. Twentyseven new unique 2MF isomers found with the 6-31+G(d,p) basis set and 18 more with the 6-31++G(d,p) basis set are shown in Figure 2b and Figure 2c, respectively. Thus, overall we found 61 isomers of 2MF. In two cases transition states were located first instead of local minima. This happened for the I15 and I30 structures where transition states (TS) between two methyl group rotamers were found in the SEGO searches. Reoptimizing these TS structures with a slightly rotated −CH3 group led to the stable isomers with the same skeleton shape. Among 61 unique 2MF isomers found there were 3methylfuran (I1, 3MF) and 22 carbenes. We expected to find five-membered cyclic carbenes which were supposed to be the proton transfer products in 2MF. Indeed 3 such carbenes were found (structures: I2, I3, and I4 in Figure 2a). However, we also found some unexpected isomers with very exotic structures. In all SEGO searches we found 17 unique isomerization channels. They are presented schematically in Scheme 1a. Each of these channels corresponds to one SEGO search. Figure 1a shows an example of such a channel obtained from the SEGO B3LYP/6-31G(d,p) search where the maximum external force allowed on any atom was set up to 0.200 au. Such external forces were generated according to eq 1 on atoms in 2MF at the equilibrium geometry first. They are shown in Table 1. It can be seen that in this case the maximum external force appears on the oxygen atom. With these external forces applied to atoms the EGO step begins.

Table 1. External Forces (au) Applied to Atoms in 2Methylfuran at the Equilibrium Geometrya atom

force X

force Y

force Z

norm

O1 C2 C3 C4 C5 C6 H7 H8 H9 H10 H11 H12

−0.0083668 −0.0897136 −0.0355318 0.0547744 0.1083822 −0.0243203 −0.0834378 −0.0203024 −0.0203023 −0.0682861 0.0777034 0.1094018

0.1998249 0.0640161 −0.1379468 −0.1379515 0.0865228 0.0051382 −0.0773084 0.0661709 0.0661709 −0.0997057 −0.0926286 0.0576990

0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 −0.0896040 0.0896040 0.0000000 0.0000000 0.0000000

0.2000000 0.1102116 0.1424494 0.1484279 0.1386827 0.0248572 0.1137473 0.1132239 0.1132239 0.1208480 0.1209044 0.1236848

a

The maximum force allowed on any atom is limited to 0.200 au, and the resulting scaling factor is 0.155226.

The system moves from the 2MF equilibrium geometry (initial local minimum) toward the energy barrier which is crossed at the optimization cycle number 24. External forces are switched off at the cycle 76 where the EGO step ends, and the stressed system is relaxed in the following SGO step which locates the first new local minimum at the cycle number 89. This local minimum corresponds to the cyclic carbene. The remaining four local minima were located at the optimization cycles: 196, 279, 361, and 453. All of them also correspond to carbenes but not to cyclic anymore. There was a ring opening when going from the first local minimum “89” to the second “196”. All five local minima correspond to the 2MF isomers. Decomposition Channels. In the case of the decomposition of 2MF or its isomer into fragments, the energy minimum on PES found in the SGO (no external forces) step and seen on the optimization history plot corresponds to a system of interacting fragments. These fragments are separated from each other usually by a distance of ∼3.0 Å. Such an energy minimum is not necessarily a stationary point on PES. Even if forces on all nuclei are negligible, then it may happen that such a decomposed system would have a few imaginary vibrational frequencies resulting from molecular fragment interactions. However, very often such an energy minimum corresponding to the decomposed 2MF was also a true local minimum with all real frequencies. It should be emphasized that with the same decomposition fragments we may find on PES a number of different local minima corresponding to a decomposed system with differently oriented fragments. Twenty-two sets of the unique decomposition products found in SEGO searches with the 6-31G(d,p) basis set are shown in Figure 3a. Twenty more sets were found with the 631+G(d,p) and 6-31++G(d,p) basis sets. They are presented also in Figure 3b and Figure 3c. Overall we found 42 unique main fragments of the 2MF decomposition. In 16 cases the main fragments are carbenes, both cyclic and noncyclic. The geometry of each main fragment of the 2MF decomposition was reoptimized. Some of these fragments also have very exotic structures. It should be emphasized that these 42 unique main fragments were obtained NOT directly from the decomposition of 2MF; some of them were obtained in consecutive stages of the five-stage processes. In all SEGO searches we found 13 unique decomposition channels. They are presented schematically in Scheme 1b. 3459


Article


Figure 3. a. Products of the enforced decomposition of 2MF obtained at DFT/B3LYP/6-31G(d,p). b. Products of the enforced decomposition of 2MF obtained at DFT/B3LYP/6-31+G(d,p). c. Products obtained from 2MF by the enforced decomposition with the 6-31++G(d,p) basis set.

Each of these channels corresponds to one SEGO search. An example of such a decomposition channel is shown in Figure 1

(no. 8) where all energy minima correspond to the systems consisting of two or three fragments: 3460


Article


results concerning stability of 61 isomers are presented in Table S1 in the Supporting Information. Among 41 main decomposition fragments stable at the B3LYP level, stability of 38 were conformed at the MP2 level with all real vibrational frequencies. The D33 fragment which is an isomer of C4H6 turned out to be a transition state with one imaginary frequency (−69.91 cm−1). A very special nature of the D27 fragment should be noted here. It is a chain isomer of C3H2. At the B3LYP level it has 3 divalent carbon atoms with two equal C−C bonds (1.294 Å) in H−C−C−C−H. However, at the MP2 level these two CC bonds were rearranged to C−C and CC with the lengths of 1.373 and 1.242 Å, respectively, and the D27 isomer became the H−C− CC−H carbene with one divalent carbon atom. The MP2 results for the main decomposition fragments of 2MF are shown in Table S2 in the Supporting Information. Finally, in order to have a uniform characteristic of the 2MF PES, we recalculated at the B3LYP/6-31++G(d,p) level energies and vibrational frequencies of the isomers obtained before with the smaller basis sets. These results are included in Table S3 in the Supporting Information. In some cases the energy and vibrational frequencies of two isomers are the same (for instance: I20, I30, etc.), but they are different enantiomers.

2MF → D4 + CH4 → D27 + CH4 + CO → D27 + CH4 + CO → D27 + CH4 + CO → D23 + CH4

Mixed Channels. In all our SEGO searches we found 18 unique channels where both isomerization and decomposition of 2MF took place. All unique mixed channels are presented schematically in Scheme 1c. Some of these cases are very interesting, because after decomposition of a given structure (2MF or its isomer) external forces put fragments together and again create an isomer. An example of such a mixed channel is shown in Figure 1c (no. 13: 2MF → I17 → D38 + H2 → I17 → I18 → I20). Verification of Stability. As mentioned before, some of the 2MF isomers and decomposition fragments found in these studies have unexpected and very exotic structures. We found cyclic structures containing one, two, or even three rings. Structures with one ring may have 3-membered (I26, I27, I44, I48, I55), 4-membered (I5, I6, I7, I28, I54), and 5-membered (I1−I4, I17, I33, I34, I56) rings. There are also structures with two rings: 3- and 4-membered (I8, I9, I18, I21 (inside carbene), I22 (inside carbene), I25, I35, I38, I45, I47, I49, I51, I52, I53), 3- and 5-membered (I16, I37, I40), and 4- and 4membered (I19, I20, I24) rings. We even found one structure with three rings: 3- and 3- and 4-membered (I46). There are noncyclic structures containing chains with 6, 5, or 4 heavy atoms. There are chains with triple CC bond moieties (I13, I29, I30, I57), with the side divalent carbon (I11, I14, I15, I31, I32) and inside divalent carbon (I12). There are also exotic species among the main decomposition fragments, for instance, the five-membered D4 ring or sixmembered D9 ring with two and three divalent carbons, respectively (see Figure 3). Some carbene-type fragments contain side divalent carbons C (D8, D22) and also −C−H (D18). Such carbenes are supposed to be even less stable than cyclic ones. Among 61 2MF isomers found in these studies 60 were predicted to be stable at the B3LYP level (local minima on PES) with all real vibrational frequencies. The I57 isomer had one imaginary frequency (−99.89 cm−1 with 6-31+ +G(d,p)), so its structure was a transition state (TS), not a local minimum on the B3LYP/6-31++G(d,p) PES. At the B3LYP level 41 main decomposition fragments were found stable with all real vibrational frequencies. One decomposed system as a whole, D8+CH2CO, was found to be a local minimum on PES with all real frequencies; however, its main fragment D8 alone was not stable and transferred into the D7 stable structure. In two other decomposed systems, namely D25+CO and D25A+CO, their main fragments D25 and D25A were different: D25 was planar, while D25A was not. However, after the geometry reoptimization, the D25 fragment went into the D25A stable structure. Thus, the SEGO method found 41 unique, stable main decomposition fragments of 2MF at the B3LYP level. Because most of the isomers and main decomposition fragments of 2MF found in these studies have very strange structures, we verified their stability at the MP2 level of theory. Molecular geometries obtained in the SEGO searches, i.e., optimized at B3LYP, were reoptimized at the MP2 level. It turned out that all 60 isomers found stable before are also stable at the MP2 level. Surprisingly, the I57 isomer which was TS at the B3LYP level went also into a stable structure with the −O−C−CC−H moiety slightly out of the skeleton plane which was planar at the B3LYP level. The MP2

■

CONCLUSIONS The first conclusion which arises from the present studies is that the 2-methylfuran potential energy surface is extremely rich. We found there were 61 stable isomers and 42 different main decomposition products. This was for us a big surprise to find so much information for such a small molecule (only 6 heavy atoms and 6 hydrogens). Moreover, we are quite convinced that our findings are not complete. We expected to find in these studies some exotic species on 2MF PES, namely carbene type isomers of 2MF. Indeed, we found such structures, but we also found a number of non-carbene structures. Some of them, both carbene and non-carbene moieties, were completely unexpected and probably previously unknown. Probably some of them would never have been found without the SEGO method. What should be emphasized is that all our findings were done automatically; we constructed “by hand” just one molecular model for 2-methylfuran, and all additional 60 isomers and 41 main decomposition fragments were obtained as the result of the enforced chemical reactions driven by the SEGO method. The enforced isomerization and decomposition of 2MF led to formation of various structures belonging to different chemical compound classes with different functional groups and different bond systems. There are bi- and tricyclic compounds with both homo and hetero 3-, 4-, 5-, and 6-membered rings as well as chain structures. One can find there a variety of functional groups (alcohol, ketone, aldehyde, ether, or epoxide) combined with the presence of double and triple bonds between carbons (alkene, alkyne, diene, and allene) and also carbenes with one, two, and even three divalent carbon atoms. Our studies show that the potential energy surface of 2-methylfuran is much more feature-rich than it was thought. With the SEGO method we found many more isomers, decomposition fragments, and many more possible isomerization and decomposition channels of 2MF than previously considered by Lifshitz et al.,5 Somers et al.,9 and Tranter et al.8 in their studies on the 2MF pyrolysis. According to those studies, the major 2MF decomposition channel begins with formation of carbenes via 3461


Article


external ones. With big external forces residing on nuclei, the geometry optimization will never converge. In our studies we had very often such cases where in the SEGO search the EGO step (with external forces on) did not finish by locating a minimum on a force-modified PES (FMPES) but was stopped before that (by limiting number of optimization cycles), after a potential energy barrier was crossed. However, the corresponding stressed structure was always good enough to begin from it a following SGO (external forces off) step and locate a next local minimum on PES. Thus, it is not necessary to go all the way down on FMPES and locate a minimum there. It is enough to cross a potential energy barrier, and it brings significant savings in SEGO calculations. It is obvious that an outcome of a SEGO search depends very much on external forces. In some sense external forces define the “enforced reaction path” on PES, and even a small change in external forces, their magnitude, and direction may result in very different paths leading to different local minima. On the one hand, this can be very useful because one may find more stationary points on PES. On the other hand, it makes the SEGO method very sensitive to essentially everything: basis set, geometry changes during optimization, choice of the coordinates in geometry optimization, accuracy of calculations, etc., and this is not very convenient. It would be better to have simpler and less sensitive external forces. Nevertheless, the SEGO method seems to be very useful and efficient in exploration of the molecular potential energy surfaces.

the proton transfer in the 2MF ring, followed by elimination of carbon monoxide. As a result, two isomers of C4H6 are formed: 1-butyne and 1,3-butadiene. We found these isomers in our SEGO searches (D26 and D25), but we also found four more isomers of C4H6: D1, D2, and D30 (unsaturated 3-membered rings), and D33 (the chain carbene). We also found that carbon monoxide may be eliminated from 2MF together with another small fragment: 2MF → D35 + CO + H2 and 2MF → D27 + CO + CH4 where the main fragments D35 and D27 (carbene) are highly unsaturated chains. The 2MF decomposition models considered in the studies cited above assumed that the smaller fragment subtracted from 2MF may be one of the following: CO, H2CCO, H2C CH2, HCCH, and also some radicals. Our studies were limited to the closed-shell systems only so we could not get any radicals. We also did not find any decomposition channel with ethylene release. However, we found in SEGO searches channels with CO, H2CCO, and HCCH and also with CH4, H2CO, H2CC, H2O, and H2. In fact out of 42 unique decomposition fragments 21 were obtained in the processes like 2MF → D + H2 with the H2 release. We hope that our findings will help in developing a deeper and more detailed model of the 2MF pyrolysis. The SEGO method turned out to be very efficient in locating local minima on the 2MF PES as well as in predicting possible decomposition channels and products. All the above was achieved with SEGO in an automatic manner without much of the user’s interference. The efficiency of the SEGO method is primarily due to the nature of the external forces employed. These forces have to do the job of moving a system from its current minimum and crossing the energy barrier separating it from another minimum on PES. To do so they must have the right orientation. This is provided by the energy Hessian used to generate external forces in SEGO. These external forces at any geometry r on PES are given as a product of H(r0) and r yielding −H(r0)r. The energy Hessian H(r0) is calculated at the initial minimum r0, and therefore it contains information concerning the curvature of that minimum; and this fact ensures proper force orientations at least at the beginning of the EGO step. It should be emphasized that external forces in SEGO are very complex. Moving up from an initial minimum, toward an energy barrier, these forces have a pushing character. They squeeze a molecule until the energy barrier is crossed. Looking at the animated picture of the SEGO optimization, one can see amazing movements of nuclei especially around the potential energy maximum. We encourage the readers to take a look at the move included in the Supporting Information showing the animated picture of one of the SEGO runs. Sometimes it is even surprising that a molecule survives such severe treatment. However, very often a molecule or more precisely a molecular system under a large amount of mechanical stress does not survive, and two things may happen. A system undergoes decomposition into fragments, or the electronic structure of a system cannot be arranged and the Hartree−Fock or Kohn− Sham Self-Consistent-Field equations cannot be solved. The latter case is possible even with relatively small forces if they try to change molecular geometry in a wrong, not acceptable way. The first scenario occurs when the applied external forces are big. Too big external forces cause also another problem. If they are too big, then a molecule cannot change its geometry enough to generate internal forces which would cancel the

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.9b00352. Figure S1, Molecular Systems Obtained from 2MF (Convention: Applied External Forces-Optimization cycle number) with 6-31G(d,p) Basis Set; Table S1, Energy and Frequency Values Obtained for All SEGO Isomers at B3LYP and MP2 Levels with Different Basis Sets; Table S2, DFT/B3LYP and MP2 Energy and Frequencies for SEGO Decomposition Products Obtained at DFT/B3LYP Level; Table S3, Energy Values and Lowest/Highest Vibrational Frequency for All Isomers Calculated at DFT/B3LYP/6-31++G(d,p) (PDF) SEGO movie (MP4)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Agnieszka Brzyska: 0000-0003-3088-3425 Krzysztof Woliński: 0000-0002-0093-4271 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A.B. and K.W. acknowledge partial financial support of the project by the statutory research fund of ICSC PAS and UMCS, respectively. 3462


Article


■

(22) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140 (4A), A1133− A1138. (23) Becke, A. D. Perspective: Fifty Years of Density-Functional Theory in Chemical Physics. J. Chem. Phys. 2014, 140 (18), 18A301. (24) Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98 (7), 5648−5652. (25) Brémond, É .; Savarese, M.; Su, N. Q.; Pérez-Jiménez, Á . J.; Xu, X.; Sancho-García, J. C.; Adamo, C. Benchmarking Density Functionals on Structural Parameters of Small-/Medium-Sized Organic Molecules. J. Chem. Theory Comput. 2016, 12 (2), 459−465. (26) Hariharan, P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theor. Chim. Acta 1973, 28 (3), 213−222. (27) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. Self-consistent Molecular Orbital Methods. XXIII. A Polarization-type Basis Set for Second-row Elements. J. Chem. Phys. 1982, 77 (7), 3654−3665. (28) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Selfconsistent Molecular Orbital Methods. XX. A Basis Set for Correlated Wave Functions. J. Chem. Phys. 1980, 72 (1), 650−654. (29) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. V. R. Efficient Diffuse Function-Augmented Basis Sets for Anion Calculations. III. The 3-21+G Basis Set for First-Row Elements, LiF. J. Comput. Chem. 1983, 4 (3), 294−301. (30) Head-Gordon, M.; Pople, J. A.; Frisch, M. J. MP2 Energy Evaluation by Direct Methods. Chem. Phys. Lett. 1988, 153 (6), 503− 506. (31) Sæbø, S.; Almlöf, J. Avoiding the Integral Storage Bottleneck in LCAO Calculations of Electron Correlation. Chem. Phys. Lett. 1989, 154 (1), 83−89. (32) Baker, J.; Wolinski, K.; Janowski, T.; Saebo, S.; Pulay, P. PQS, Parallel Quantum Solutions: Box 293, Fayetteville, Arkansas 727020293, U.S.A.

REFERENCES

(1) Xiao, H.; Yang, X.; Hou, B.; Wang, R.; Xue, Q.; Ju, H. Combustion Performance and Pollutant Emissions Analysis of a Diesel Engine Fueled with Biodiesel and Its Blend with 2Methylfuran. Fuel 2019, 237, 1050−1056. (2) Hoppe, F.; Burke, U.; Thewes, M.; Heufer, A.; Kremer, F.; Pischinger, S. Tailor-Made Fuels from Biomass: Potentials of 2Butanone and 2-Methylfuran in Direct Injection Spark Ignition Engines. Fuel 2016, 167, 106−117. (3) Davis, A. C.; Sarathy, S. M. Computational Study of the Combustion and Atmospheric Decomposition of 2-Methylfuran. J. Phys. Chem. A 2013, 117 (33), 7670−7685. (4) Cheng, Z.; He, S.; Xing, L.; Wei, L.; Li, W.; Li, T.; Yan, B.; Ma, W.; Chen, G. Experimental and Kinetic Modeling Study of 2Methylfuran Pyrolysis at Low and Atmospheric Pressures. Energy Fuels 2017, 31 (1), 896−903. (5) Lifshitz, A.; Tamburu, C.; Shashua, R. Decomposition of 2Methylfuran. Experimental and Modeling Study. J. Phys. Chem. A 1997, 101 (6), 1018−1029. (6) Urness, K. N.; Guan, Q.; Golan, A.; Daily, J. W.; Nimlos, M. R.; Stanton, J. F.; Ahmed, M.; Ellison, G. B. Pyrolysis of Furan in a Microreactor. J. Chem. Phys. 2013, 139 (12), 124305. (7) Grela, M. A.; Amorebieta, V. T.; Colussi, A. J. Very Low Pressure Pyrolysis of Furan, 2-Methylfuran and 2,5-Dimethylfuran. The Stability of the Furan Ring. J. Phys. Chem. 1985, 89 (1), 38−41. (8) Tranter, R. S.; Lynch, P. T.; Randazzo, J. B.; Lockhart, J. P. A.; Chen, X.; Goldsmith, C. F. High Temperature Pyrolysis of 2-Methyl Furan. Phys. Chem. Chem. Phys. 2018, 20 (16), 10826−10837. (9) Somers, K. P.; Simmie, J. M.; Metcalfe, W. K.; Curran, H. J. The Pyrolysis of 2-Methylfuran: A Quantum Chemical, Statistical Rate Theory and Kinetic Modelling Study. Phys. Chem. Chem. Phys. 2014, 16 (11), 5349. (10) Wolinski, K. Exploring Potential Energy Surface with External Forces. J. Chem. Theory Comput. 2018, 14 (12), 6306−6316. (11) Wolinski, K.; Baker, J. Theoretical Predictions of Enforced Structural Changes in Molecules. Mol. Phys. 2009, 107 (22), 2403− 2417. (12) Wolinski, K.; Baker, J. Geometry Optimization in the Presence of External Forces: A Theoretical Model for Enforced Structural Changes in Molecules. Mol. Phys. 2010, 108 (14), 1845−1856. (13) Cybulska, J.; Brzyska, A.; Zdunek, A.; Woliński, K. Simulation of Force Spectroscopy Experiments on Galacturonic Acid Oligomers. PLoS One 2014, 9 (9), e107896. (14) Brzyska, A.; Płaziński, W.; Woliński, K. Force-Induced Structural Changes in Non-Sulfated Carrageenan Based Oligosaccharides − a Theoretical Study. Soft Matter 2018, 14 (30), 6264− 6277. (15) Brzyska, A.; Woliński, K. Enforced Conformational Changes in the Structural Units of Glycosaminoglycan (Non-Sulfated HeparinBased Oligosaccharides). RSC Adv. 2014, 4 (69), 36640−36648. (16) Baker, J.; Wolinski, K. Isomerization of Stilbene Using Enforced Geometry Optimization. J. Comput. Chem. 2011, 32 (1), 43−53. (17) Baker, J.; Wolinski, K. Kinetically Stable High-Energy Isomers of C14H12 and C12H10N2 Derived from Cis-Stilbene and CisAzobenzene. J. Mol. Model. 2011, 17 (6), 1335−1342. (18) Beyer, M. K. The Mechanical Strength of a Covalent Bond Calculated by Density Functional Theory. J. Chem. Phys. 2000, 112 (17), 7307−7312. (19) Ong, M. T.; Leiding, J.; Tao, H.; Virshup, A. M.; Martínez, T. J. First Principles Dynamics and Minimum Energy Pathways for Mechanochemical Ring Opening of Cyclobutene. J. Am. Chem. Soc. 2009, 131 (18), 6377−6379. (20) Ribas-Arino, J.; Shiga, M.; Marx, D. Understanding Covalent Mechanochemistry. Angew. Chem., Int. Ed. 2009, 48 (23), 4190− 4193. (21) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136 (3B), B864−B871. 3463


Isomerization and Decomposition of 2-Methylfuran with

Recommend Documents