Article pubs.acs.org/EF
ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Lignite Depolymerization in Supercritical Methanol with LigniteRelated Model Compounds Bo Chen, Xian-Yong Wei,* Zhu-Sheng Yang, Chang Liu, Xing Fan, Yu Qing, and Zhi-Min Zong Key Laboratory of Coal Processing and Efficient Utilization (Ministry of Education), China University of Mining and Technology, Xuzhou 221008, Jiangsu, People’s Republic of China ABSTRACT: To investigate the detailed mechanisms for lignite methanolysis, we used ReaxFF reactive force field to perform a series of molecular dynamics simulations (MDSs) on a unimolecular model compound. The α-O-4 and β-O-4 types of ligniterelated model compounds were selected as representatives of linkages in lignites. The reaction products predicted by ReaxFF MDSs are consistent with those from experimental results reported. The initiation reaction observed in ReaxFF MDSs involving the ether linkage cleavage and methanol participation closely matches the results observed from previously reported experiments. The agreement of these results with available experimental observations demonstrates that ReaxFF MDSs can give an atomistic description of the initiation mechanism for methanolysis and provide useful insights into the complicated reaction processes.
1. INTRODUCTION The lignite reserve is estimated to be more than one trillion tons over the world.1 However, the low calorific value, high ash yield, and water content limit the industrial application of lignites, especially for thermoelectricity generation.2 It is necessary to develop efficient conversion processes to minimize such disadvantages. Supercritical methanolysis has been widely used for biodiesel production3−5 and polyethylene terephthalate depolymerization.6,7 Because organic matter in lignites is rich in oxygen-bridge bonds (OBBs), cleaving OBBs in lignites should be an crucial step for converting lignite to clean fuel and value-added chemicals.8 Understanding the mechanisms of lignite degradation in supercritical methanol (SCM) is important for optimization of the conversion process and efficient use of liquid products. The depolymerization of woody biomass (WBM) using model compounds in SCM9−13 is a typical example for related investigations. Lignites retain more integrate macromolecular structures of coal-forming plants and have larger amounts of oxygen-containing moieties,14,15 including OBBs, than higher rank coals. Therefore, WBM-related model compounds could be appropriate for studying lignite depolymerization in SCM. Recently, the ReaxFF reactive force field,16−20 a force field method capable of describing bond breaking, bond formation, and chemical reactivity, has been developed as an effective method for complicated chemical reactions and other complicated phenomena under extreme conditions. ReaxFF allows for molecular dynamics simulations (MDSs) for computational costs nearly as low as for simple force fields and reproduces nearly the accuracy of quantum mechanics for reaction systems, including reactants, transition states (TSs), and products, for thermal decomposition of lignite,19,20 thermal-induced and shock-induced chemistry in energetic material 1,3,5-trinitro-1,3,5-triazine,21,22 thermal initiation chemical events of triacetonetriperoxide,23 and pyrolysis and combustion of JP-10 jet fuel.24 To elucidate the mechanisms and the dominant reaction pathways (RPWs) for the © 2012 American Chemical Society
decomposition of lignites or lignin in SCM and describe the formation processes of intermediates and final products, we performed ReaxFF MDSs with N (number), V (volume), and T (temperature) constants and designated these conditions as NVT. Here, (1S,2S)-1-(4-hydroxy-3-methoxyphenyl)-2-(2methoxyphenoxy)propane-1,3-diol (1) and 2-methoxy-4-((2methoxyphenoxy)methyl)phenol (2) were selected as model compounds for studying the cleavage of β-O-4 and α-O-4 linkages, respectively, in lignites.
2. COMPUTATIONAL APPROACH Organic matter in lignites consists of various aromatic units with alkyl side chains and oxygen-containing functional groups, such as hydroxyl, methoxy, carbonyl, and carboxylic groups.14,15 These units are usually linked by OBBs. Hence, compounds 1 and 2 were used as ligniterelated model compounds. To simulate the depolymerization of model compound 1 or 2 in SCM, we placed one model compound molecule and 78 methanol molecules in a 20 × 20 × 20 Å periodic box (PB; Table 1). The use of a single model compound can obtain detailed evolution of the chemical reactions related to methanolysis.
Table 1. Simulation Conditions for the ReaxFF NVT-MDSsa system
model compound
model compound number
initial density (g cm−3)
S1 S2 S3
1 2
1 1
0.52 0.59 0.57
a
S2 includes S2i−S2iv; methanol number, 78; volume, 8000 Å3.
The temperature-programmed MDSs were processed using a velocity Verlet approach with a time step of 0.25 fs. A Berendsen thermostat with a damping constant of 100 fs was used for temperature control.25 Each reaction system (Table 1) was energyReceived: August 13, 2011 Revised: January 4, 2012 Published: January 4, 2012 984
dx.doi.org/10.1021/ef201234j | Energy Fuels 2012, 26, 984−989
Energy & Fuels
Article
Scheme 1. RPWs for Methanolysis Observed in the ReaxFF NVT-MDS
minimized and subsequently equilibrated at 800 K in a 10 ps NVTMDS. Then, the equilibrated system was performed 100 ps NVTMDS for the cook-off process, in which the system was heated from 800 to 2200 K at a rate of 100 K · ps−1 in the first 14 ps and then kept at 2200 K for another 86 ps. For all of the possible RPWs available in ReaxFF, we processed multiple simulations for each reaction system and analyzed the trajectories (every 0.25 fs), from which the system temperature, potential energy, fragment number, and molecular species, including intermediates, could be taken directly to obtain the most favorable RPWs. For S2 (Scheme 1), four RPWs, i.e., i−iv, are the most kinetically and thermodynamically favorable. The corresponding reaction systems are defined as S2i−S2iv, respectively. To observe the chemical reactions within a computationally practical simulation time, we chose an elevated temperature beyond normal experimental conditions, which could change the reaction rates rather than the preferred RPWs.26,27 The intermediates and products formed during the simulation were analyzed with a 0.3 bond order cutoff for recognition of molecules. Although the above method can be used to calculate the potential energy barrier (PEB) of hydroxyl elimination, its ensemble averaging is computationally complex because of the multiple RPWs. The sliding restraint (SR) proved to be a relatively simple and computationally inexpensive method for the ReaxFF system energy description.28,29 Thus, we applied this method to calculate the PEB of hydroxyl elimination. To derive the PEB, we augmented the ReaxFF potential with an additional SR to the α−C−O bond to force hydroxyl elimination
the product. The thermodynamic effects of the overall reaction energy and PEB can be determined by monitoring the system energy during the restrained MDS.
3. RESULTS AND DISCUSSION A NVT-MDS was employed for S1 (Table 1), and it was found that a few methanol molecules were converted to formaldehyde (Table 2) during the simulation time, which is in agreement with experimental results reported.30 However, Chenoweth et al.31 reported that no reaction occurred during a 1000 ps NVTMDS at 2000 K with just 30 methanol molecules in a 20 × 20 × 20 Å PB. This disagreement might be induced by a higher temperature and density, which resulted in methanol degradation via a methanol dimer because of a higher potential energy. With the rise in the temperature from 800 to 2200 K in the first 14 ps cook-off process, the potential energy increased continuously. However, the potential energy leveled off in the rest 86 ps of the simulations, where the system temperature was kept at 2200 K and the molecular system was stable (Figure 1). In fact, during the 100 ps NVT-MDS, the total fragment change was less than 17% and the consumed methanol molecules were no more than 12 (Figure 2). As Figure 1 shows, the potential energy was decreased by mixing either model compound 1 or 2 with methanol, suggesting that intermolecular attraction existed between the model compound and methanol. Furthermore, the decomposition of the methanol dimer was only observed in S1 (Scheme 2) because of the higher potential energy of S1 than that of S2 or S3. As Scheme 1 illustrates, the initial reaction for the depolymerization of model compounds 1 and 2 in SCM
Erestraint = k1{1 − exp[−k2(rij − R restraint)2 ]} where Rrestraint, k1, k2, and rij denote the restraint distance, the bond value (31 401 kJ mol−1), the restraint force constant (0.25 Å−2), and the distance between C and O in the α−C−O bond, respectively. In the SR formulation, Rrestraint shifted gradually along the reaction coordinate, so that rij changed from the value of the reactant to that of 985
dx.doi.org/10.1021/ef201234j | Energy Fuels 2012, 26, 984−989
Energy & Fuels
Article
Then, the nucleophile (i.e., alkoxide) attacks the positively charged atom, e.g., the carbon atom in the carbonyl group from the reactant, leading to the formation of a corresponding intermediate. For model compound 1, the cleavage of ether linkage leads to the formation of intermediates 2b and 1a. Intermediate 2b abstracts hydrogen from the hydroxyl group in methanol to produce product 7, which can further lose CH3•. The hydrogen abstraction from the methyl group in methanol leads to the formation of product 8. Phenolic compounds are the most abundant products from the methanolysis of both lignite and biomass.38−40 In comparison to intermediate 2b, intermediate 1a undergoes complex reactions via different RPWs, i.e., i−iv, to form distinct products, as shown in Scheme 1. As illustrated in Figure 4, conversion of intermediate 1a to 1b in RPW i requires ca. 35.6 kJ mol−1. Hydrogen abstraction from phenolic hydroxyl by CH3O• affords a biradical intermediate 1e and methanol via TS II (Figure 3), which has a PEB of ca. 79.5 kJ mol−1 with respect to the corresponding intermediate. The localization of unpaired electron at branched α−C makes intermediate 1e labile, so that intermediate 1e can undergo intramolecular rearrangement (IMRA) of the carbon skeleton by crossing a PEB of ca. 113.0 kJ mol−1 to produce intermediate 1g via TS 1f. Although both intermediates 1e and 1g are biradical fragments, intermediate 1g is nearly 104.7 kJ mol−1 thermodynamically stable compared to intermediate 1e by delocalizing its unpaired electrons onto adjacent conjugated π electrons in the benzene ring. The successive H migration from methanol to intermediate 1g leads to the formation of intermediate 1h via TS III (Figure 3), which is consistent with the experimental result that methanol participated in the reaction as a hydrogen donor through the radical reaction mechanism.41 TS III has a PEB of 47.3 kJ mol−1 over intermediate 1g and properly connects 1g with 1h. Subsequently, the OH• elimination from intermediate 1h involving homolytic α−C−OH cleavage requires an additional energy of ca. 143.6 kJ mol−1, leading to the formation of compound 3. In RPW ii, the elimination of γ−OH from intermediate 1a affords compound 4 via TS I (Figure 3) by crossing a PEB of ca. 116.8 kJ mol−1. This process involves the formation of the hydrogen bond between intermediate 1a and methanol and subsequent elimination of H2O. The final products are ca. 88.8 kJ mol−1 endothermic from intermediate 1a + methanol. It is noteworthy that compound 4, an isomer of compound 3, has never been identified in previous reports.12,13 Hydroxyl migration from α−C to β−C via TS 1c with a PEB of ca. 52.3 kJ mol−1 leads to the formation of intermediate 1i, as shown in RPW iii (Scheme 1). Subsequently, the IMRA of the carbon skeleton, leading to the formation of intermediate 1k, requires a PEB of ca. 158.7 kJ mol−1 with respect to intermediate 1i. The structures of intermediates 1i and 1k are similar to those of intermediates 1g and 1e, respectively, in which hydrogen is abstracted from the phenolic hydroxyl group by CH3O•. Hydrogen abstraction from γ−OH of intermediate 1k by CH3O• results in intermediate 1l + methanol via TS IV (Figure 3) by crossing a PEB of ca. 11.7 kJ mol−1. Further elongation of the C−C bond between α−C and γ−C in intermediate 1l produces compound 5 + formaldehyde as the final product. The difference in energy between the final product and intermediate 1a + CH3O• is ca. 167.5 kJ mol−1. Dehydrogenation of intermediate 1a involving homolytic γ− C−H cleavage requires ca. 251.6 kJ mol−1. Thus, in RPW iv,
Table 2. Products Observed in the Final Configuration from the ReaxFF NVT-MDS
Figure 1. Potential energy change for different systems.
without catalyst involves the cleavage of ether linkage, which is in agreement with non-catalytic experiments12,13 but different from the catalytic methanolysis involving nucleophilic attack by CH3O− generated from the catalytic dissociation of methanol.32−37 The general mechanism for the catalyzed methanolysis involves an interaction between the catalyst and methanol, leading to the formation of the corresponding alkoxide and protonated catalyst in the initial reaction step. 986
dx.doi.org/10.1021/ef201234j | Energy Fuels 2012, 26, 984−989
Energy & Fuels
Article
Figure 2. (a−c) Total fragments and (d−f) methanol consumed analyses for ReaxFF MDS of methanolysis under different systems.
carbon rearrangement has the same structure as intermediate 1h derived from hydroxyl rearrangement. Although four kinds of final products were observed in MDS, only product 3 was identified in the experiments reported.12,13 When the structure of product 3 is compared to that of intermediate 1a, it is reasonable to consider that the final product 3 can be derived from direct elimination of hydroxyl from intermediate 1a. Thus, the ReaxFF potential with an additional SR to the α−C−O bond, which was described in detail by van Duin et al.,28,29 was performed to monitor the stable intermediates as well as PEBs related to TSs. The α−OH in intermediate 1a favors ortho migration, leading to the formation of stable intermediate 1i or 1h with the lowest potential energy, as exhibited in Figure 5. For model compound 2, the cleavage of ether linkage affords intermediates 2b and 2a. Subsequently, methyl abstraction from methanol by intermediate 2a yields 4-ethyl-2-methoxyphenol (9) as the final product rather than the experimentbased reaction product 2-methoxy-4-(methoxymethyl)phenol, which may result from the methoxylation of intermediate 2a by the addition of CH3O•. In our MDS, few methoxy intermediates were observed without catalyst, whereas the formation of CH3O• from methanol was reported in Fe- and Ni-catalyzed reactions;42,43
Scheme 2. Methanol Autodecomposition Observed in the ReaxFF NVT-MDS for S1
there is no appreciable PEB for the formation of intermediate 1d and final product 6. However, benzaldehyde and phenylmethanol with similar structures to compounds 5 and 6, respectively, were identified as the products from supercritical methanolysis of dibenzyl ether.10 Except the cleavage of ether linkage, the rate-controlling steps are the elimination of α−OH from intermediate 1h in RPW i and the formation of intermediate 1j from intermediate 1i in RPW iii. In fact, the stable intermediate 1i derived from
Figure 3. Structures of the fragments observed in the ReaxFF NVT-MDS. 987
dx.doi.org/10.1021/ef201234j | Energy Fuels 2012, 26, 984−989
Energy & Fuels
Article
provide a reasonable atomistic description of the initiation mechanism for methanolysis.
■
AUTHOR INFORMATION
Corresponding Author
*Telephone: +86-516-83995916. Fax: +86-516-83884399. Email:
[email protected].
■
ACKNOWLEDGMENTS This research work was supported by National Basic Research Program of China (Grant 2011CB201302), the Fund from National Natural Science Foundation of China for Innovative Research Group (Grant 50921002), the Key Project of Natural Science Foundation of China (Grant 20936007), the Key Project of Coal Joint Fund from National Natural Science Foundation of China and Shenhua Group Corporation Limited (Grant 51134021), the National High Technology Research and Development Program of China (Grant 2007AA06Z113), National Natural Science Foundation of China (Grants 20776149 and 51074153), and the Priority Academic Program Development of Jiangsu Higher Education Institutions. The authors also thank Dr. Adri van Duin of Pennsylvania State University for the use of the ReaxFF MD code and for his continuous help on the use of that code.
■
Figure 4. ReaxFF energies for model compound 1 unimolecular dissociation pathways observed in the ReaxFF NVT-MDS.
■
NOMENCLATURE OBB = oxygen-bridge bond SCM = supercritical methanol WBM = woody biomass MDS = molecular dynamics simulation TS = transition state RPW = reaction pathway NVT = number, volume, and temperature constants PEB = potential energy barrier SR = sliding restraint IMRA = intramolecular rearrangement REFERENCES
(1) Thielemann, T.; Schmidt, S.; Peter Gerling, J. Int. J. Coal Geol. 2007, 72 (1), 1−14. (2) Nolan, P.; Shipman, A.; Rui, H. Eur. Manage. J. 2004, 22 (2), 150−164. (3) Saka, S.; Kusdiana, D.; Minami, E. J. Sci. Ind. Res. 2006, 65 (5), 420−425. (4) Demirbas, A. Energy Convers. Manage. 2007, 48 (3), 937−941. (5) Issariyakul, T.; Kulkarni, M. G.; Dalai, A. K.; Bakhshi, N. N. Fuel Process. Technol. 2007, 88 (5), 429−436. (6) Yang, Y.; Lu, Y.; Xiang, H.; Xu, Y.; Li, Y. Polym. Degrad. Stab. 2002, 75 (1), 185−191. (7) Sako, T.; Sugeta, T.; Otake, K.; Nakazawa, N.; Sato, M.; Namiki, K.; Tsugumi, M. J. Chem. Eng. Jpn. 1997, 30 (2), 342−346. (8) Dong, L.; Yuan, Q.; Yuan, H. Fuel 2006, 85 (17−18), 2402− 2407. (9) Yokoyama, C.; Nishi, K.; Nakajima, A.; Seino, K. Sekiyu Gakkaishi 1998, 41 (4), 243−250. (10) Yokoyama, C.; Nishi, K.; Otake, K.; Takahashi, S. Sekiyu Gakkaishi 1994, 37 (1), 34−44. (11) Yokoyama, C.; Nishi, K.; Takahashi, S. Sekiyu Gakkaishi 1997, 40 (6), 465−473. (12) Tsujino, J.; Kawamoto, H.; Saka, S. Wood Sci. Technol. 2003, 37 (3−4), 299−307. (13) Minami, E.; Kawamoto, H.; Saka, S. J. Wood. Sci. 2003, 49 (2), 158−165. (14) Hatcher, P. G.; Breger, I. A.; Szeverenyi, N.; Maciel, G. E. Org. Geochem. 1982, 4 (1), 9−18.
Figure 5. Energy profile of the hydroxyl elimination from intermediate 1a.
i.e., Fe and Ni catalyzed the formation of CH3O• from methanol.
4. CONCLUSION Using the ReaxFF reactive force field to perform a series of MDSs, we studied unimolecular depolymerization of the β-O-4 and α-O-4 types of model compounds in SCM. The results show that the model compound depolymerization is initiated by ether linkage cleavage, followed by complicated multi-step processes of carbon skeleton rearrangement or hydroxyl migration with several reactive intermediates along the RPWs, including biradicals and hydrogen-bonded TSs. The simulations also demonstrate that methanol donates hydrogen from the hydroxyl or methyl group to other intermediate radicals. The mechanism for initiation depolymerization of model compounds and the hydrogen-donating ability of methanol observed in ReaxFF MDSs agree with the results from the previous experiments. These results validate that ReaxFF can be used to obtain insight into complicated reaction processes and 988
dx.doi.org/10.1021/ef201234j | Energy Fuels 2012, 26, 984−989
Energy & Fuels
Article
(15) Ralph, J. P.; Catcheside, D. E. A. Fuel 1993, 72 (12), 1679− 1686. (16) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. III J. Phys. Chem. A 2001, 105 (41), 9396−9409. (17) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. III J. Phys. Chem. A 2008, 112 (5), 1040−1053. (18) Kamat, A. M.; van Duin, A. C. T.; Yakovlev, A. J. Phys. Chem. A 2010, 114 (48), 12561−12572. (19) Salmon, E.; van Duin, A. C. T.; Lorant, F.; Marquaire, P. M.; Goddard, W. A. III Org. Geochem. 2009, 40 (12), 1195−1209. (20) Salmon, E.; van Duin, A. C. T.; Lorant, F.; Marquaire, P. M.; Goddard, W. A. III Org. Geochem. 2009, 40 (3), 416−427. (21) Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A. III Phys. Rev. Lett. 2003, 91 (9), No. 098301. (22) Strachan, A.; Kober, E. M.; van Duin, A. C. T.; Oxgaard, J.; Goddard, W. A. III J. Chem. Phys. 2005, 122 (5), No. 054502. (23) van Duin, A. C. T.; Zeiri, Y.; Dubnikova, F.; Kosloff, R.; Goddard, W. A. III J. Am. Chem. Soc. 2005, 127 (31), 11053−11062. (24) Chenoweth, K.; van Duin, A. C. T.; Dasgupta, S.; Goddard, W. A. III J. Phys. Chem. A 2009, 113 (9), 1740−1746. (25) Weismiller, M. R.; van Duin, A. C. T.; Lee, J.; Yetter, R. A. J. Phys. Chem. A 2010, 114 (17), 5485−5492. (26) Mueller, J. E.; van Duin, A. C. T.; Goddard, W. A. III J. Phys. Chem. C 2010, 114 (12), 5675−5685. (27) Mueller, J. E.; van Duin, A. C. T.; Goddard, W. A. III J. Phys. Chem. C 2010, 114 (11), 4939−4949. (28) van Duin, A. C. T.; Damsté, J. S. S. Org. Geochem. 2003, 34 (4), 515−526. (29) Rahaman, O.; van Duin, A. C. T.; Goddard, W. A. III; Doren, D. J. J. Phys. Chem. B 2011, 115 (2), 249−261. (30) Bazaev, A. R.; Abdulagatov, I. M.; Bazaev, E. A.; Abdurashidova, A. A.; Ramazanova, A. E. J. Supercrit. Fluids 2007, 41 (2), 217−226. (31) Chenoweth, K.; van Duin, A. C. T.; Persson, P.; Cheng, M. J.; Oxgaard, J.; Goddard, W. A. III J. Phys. Chem. C 2008, 112 (37), 14645−14654. (32) Cerro-Alarcón, M.; Corma, A.; Iborra, S.; Gómez, J. P. Appl. Catal., A 2008, 346 (1−2), 52−57. (33) Kalita, D.; Sarma, R.; Baruah, J. B. Inorg. Chem. Commun. 2009, 12 (6), 569−571. (34) Chapuzet, J.-M.; Beauchemin, S.; Daoust, B.; Lessard, J. Tetrahedron 1996, 52 (12), 4175−4180. (35) Juárez, R.; Corma, A.; García, H. Top. Catal. 2009, 52 (12), 1688−1695. (36) Bolotin, B. M.; Kunavin, N. I.; Zeryukina, L. S.; Safina, R. U. Chem. Heterocycl. Compd. 1974, 10 (6), 658−661. (37) Isaeva, Z. G.; Bakaleinik, G. A.; Kovylyaeva, G. I.; Efremov, Y. Y.; Korshunov, R. L. Russ. Chem. Bull. 1987, 36 (12), 2606−2609. (38) Tang, S. R.; Zong, Z. M.; Zhou, L.; Zhao, W.; Li, X. B.; Peng, Y. L.; Xie, R. L.; Chen, X. F.; Gu, W. T.; Wei, X. Y. Renewable Energy 2010, 35 (5), 946−951. (39) Zhao, W.; Xu, W. J.; Lu, X. J.; Sheng, C.; Zhong, S. T.; Tang, S. R.; Zong, Z. M.; Wei, X. Y. Energy Fuels 2009, 24 (1), 136−144. (40) Chen, B.; Wei, X. Y.; Zong, Z. M.; Yang, Z. S.; Qing, Y.; Liu, C. Appl. Energy 2011, 88 (12), 4570−4576. (41) Wu, B. C.; Klein, M. T.; Sandler, S. I. AIChE J. 1990, 36 (8), 1129−1136. (42) Madix, R. J. Catal. Rev. 1984, 26 (2), 281−297. (43) Kumbhar, P. S.; Rajadhyaksha, R. A. Catal. Lett. 1991, 10 (1), 131−135.
989
dx.doi.org/10.1021/ef201234j | Energy Fuels 2012, 26, 984−989
