Energy Fuels 2010, 24, 6468–6475 Published on Web 11/09/2010
: DOI:10.1021/ef100797h
Importance of the Inclusion of Dispersion in the Modeling of Asphaltene Dimers Iain D. Mackie†,‡ and Gino A. DiLabio*,†,‡ † National Institute for Nanotechnology, National Research Council of Canada, 11421 Saskatchewan Drive, Edmonton, Alberta, Canada T6G 2M9, and ‡Department of Physics, University of Alberta, 11322 89th Avenue, Edmonton, Alberta, Canada T6G 2G7
Received June 23, 2010. Revised Manuscript Received September 14, 2010
Modeling of asphaltenes presents many challenges, not least of which is their large size and the lack of definitive experimental structure information. However, a further fundamental issue of importance is the ability of the modeling methodology to accurately predict noncovalent interactions, particularly dispersion (London) forces. The self-aggregation properties of asphaltenes are primarily driven by such interactions. Therefore, for a modeling approach to be insightful, dispersion must be accounted for. This point is illustrated in this work by examining the effect of dispersion on the conformer distribution for a series of asphaltene model dimers, including perylene bisimide-type molecules and for a model of Maya asphaltene. Inclusion of dispersion using dispersion-correcting potentials on standard density-functionals is shown to be critical to both structure and conformer population for these noncovalently bound dimers. For the Maya asphaltene, a previously postulated “open” structure is shown to be ca. 9 kcal/mol less stable than a conformer that allows for greater π-π overlap within its central archipelago-type moiety. This finding is in line with recent NMR work, which indicates that “closed” structures dominate for asphaltenes. N,N0 -dimethyl-perylene bisimide dimer has a binding energy (BE) of up to 29 kcal/mol, while the more complex N,N0 -(1-hexylheptyl)-perylene bisimide model is slightly less strongly bound (BE = 25 kcal/mol) since the large alkyl substituents restrict the ability of the extensive π-regions to overlap. Such conclusions cannot be drawn when methods that do not incorporate dispersion, e.g., the B3LYP density-functional, are used.
The self-aggregation properties of asphaltenes are primarily driven by noncovalent interactions.3-10 Despite this fact, very few electronic structure modeling studies have been published on such aggregates. The likely reason for the dearth of work in this area is that molecular orbital methods that are capable of handling large systems, viz., density-functional theory (DFT), have traditionally been unable to properly deal with noncovalent, especially dispersion, interactions. Consequently, the computational modeling of these systems has been restricted to classical molecular mechanics/dynamics, thus limiting studies of the thermochemical and electronic properties of asphaltenes and their aggregates. In order to illustrate the importance of including an accurate description of dispersion binding in molecular models suitable for application to asphaltene simulations, we apply DFT with and without dispersion corrections to demonstrate the importance of noncovalent interactions to conformers of asphaltene model dimers, and ergo, sequitur, their chemistry.
Introduction Asphaltenes are a hugely significant financial quandary for oil upgrading.1 Out of solution, they can form solid deposits in wells and pipelines, thus affecting both oil recovery and transportation. In extreme cases, asphaltene blockages can completely halt oil production. Added to this is the problem of storage since the highly polar asphaltenes can contaminate water by forming emulsions, with their oxidation often resulting in sedimentation. Their structural characterization is complicated by the fact that elemental compositions from various sources of petroleum are all different2 and are dependent upon the precipitating solvent, pressure, and temperature, among other factors. The tendency of asphaltene molecules to selfassociate severely inhibits the experimental determination of asphaltene molecular weights, with the more recently accepted range now from 500 to 1000 amu.3 *To whom correspondence should be addressed. Phone: þ1-780-6411729. E-mail:
[email protected]. (1) (a) Speight, J. G. Oil Gas Sci. Technol. 2004, 59, 467–477. (b) Speight, J. G. Oil Gas Sci. Technol. 2004, 59, 479–488. (c) Abu-Khader, M. M.; Speight, J. G. Oil Gas Sci. Technol. 2007, 62, 715–722. (2) For a synopsis of asphaltene compositions precipitated by differing flocculants from Canada and Middle East sources, see Mansoori, G. A. OPEC Rev. 1988, 12, 103–113. (3) Groenzin, H.; Mullins, O. C. Energy Fuels 2000, 14, 677–684. (4) Murgich, J.; Rodriguez, J.; Aray, Y. Energy Fuels 1996, 10, 68–76. (5) Schabron, J. F.; Speight, J. G. Pet. Sci. Technol. 1998, 16, 361–375. (6) Siddiqui, M. N. Pet. Sci. Technol. 2003, 21, 1601–1615. (7) Miura, K.; Mae, K.; Li, W.; Kusakawa, T.; Morozumi, F.; Kumano, A. Energy Fuels 2001, 15, 599–610. (8) Murgich, J. Pet. Sci. Technol. 2002, 20, 983–997. (9) Yen, T. F. Energy Sources 1974, 1, 447–451. r 2010 American Chemical Society
Models for Asphaltenes First, we calculate the binding between thiophene (1) or benzene (2) to a series of simple (small) oxygen-functionalized carbonaceous ringed molecules. Included in this series are benzoic acid (3), cyclohexanone (4), phenol (5), furan (6), 2-pyrone (7), and 4-pyrone (8); see Scheme 1. These will serve as simple, yet dramatic, examples of the importance of noncovalent interactions for binding energy and structure determination. Most density-functionals are generally very good at (10) Siffert, B.; Kuczinski, J.; Papirer, E. J. Colloid Interface Sci. 1990, 135, 107–117.
6468
pubs.acs.org/EF
Energy Fuels 2010, 24, 6468–6475
: DOI:10.1021/ef100797h
Mackie and DiLabio
Scheme 1. Binding of Thiophene (1) and Benzene (2) to 3-8 Studied in the Present Work
Scheme 2. Large Polyaromatic Dimers Studied in the Present Work
Scheme 3. Components of Maya Asphaltene Model Compound
describing hydrogen-bonding and dipole-type noncovalent interactions,11 but accurately modeling dispersion (London) forces presents more of a challenge. For this, and all subsequent work described herein, we have utilized a recently developed modification to density-functional theory (DFT),12 that allows for the efficient modeling of dispersion interactions, vide infra. By also removing the description of dispersion and relaxing the previously optimized structures, we can directly infer its importance in this context. In larger molecules, such as those more pertinent to asphaltene chemistry, dispersion will likely play an even more important role in determining complex stability. In this work, this is shown using dimers of perylene and perylene bisimide compounds (pericondensed models of asphaltenes, Scheme 2) and by a recently discussed model for Maya asphaltene (Scheme 3) that possesses an archipelago motif.13,14 Application of Dispersion-Correcting Potentials in DFT Addition of dispersion-correcting potentials (DCPs) to standard density-functionals is a convenient way to remedy the long-standing problem of describing dispersion interactions using DFT.11 The details behind the development of these potentials, which have been designed for a host of functionals and basis sets, is described elsewhere.12,15,16 For this work, the addition of two Gaussian-type functions (eq 1), which resemble atom-centered effective core potentials, reproduces a weakly attractive medium to long-range character, with the second function being tighter and repulsive to prevent overbinding. DCPs model dispersion by modifying the environment in which the valence electrons move and were developed for the 6-31þG(d,p) basis set such that basis set superposition error (BSSE) that occurs in these systems is accounted for. X 2 cli rnli e - ζli r ð1Þ Ul ðrÞ ¼ r - 2 We showed that these DCPs can significantly improve the performance of a variety of density functionals for noncovalent interactions.12,15,16 For example, B971/6-31þG(d,p) with DCPs predicts the binding energies (BEs) for a set of 22 noncovalently interacting dimers with a percent absolute
deviation (%AD) of 13.8% compared to those predicted by very high-level (and computationally intensive) wave function methods (viz., complete basis set extrapolated, coupledcluster with single and double excitations plus perturbative triples, CCSD(T)/CBS). In principle, DCPs can be formulated to modify any element but have thus far only been developed for carbon12,15,16 and silicon.17 That said, the DCPs were developed by considering a wide range of noncovalently bound dimers, including many containing heteroatoms (O, N, and F). Although DCPs are not applied to the heteroatoms under study, they are only a very minor component of these dimers. Furthermore, we have previously found that the introduction of potentials to sulfur or nitrogen atoms is not required to accurately reproduce “gold-standard”, CCSD(T)-calculated18,19 binding energies for the thiophene dimer20 and for CO2 bound to nitrogencontaining heterocycles.21 When used with complete basis set extrapolations, the CCSD(T) approach is able to produce highly accurate noncovalent binding energies. However, such an approach can only be practically applied to relatively small systems because such calculations scale to the power 7 with the size of the molecule. Therefore, the ability of a methodology to accurately reproduce high-level data, both in terms of binding
(11) Johnson, E. R.; Mackie, I. D.; DiLabio, G. A. J. Phys. Org. Chem. 2009, 22, 1127–1135. (12) DiLabio, G. A. Chem. Phys. Lett. 2008, 455, 348–353. (13) Takonahashi, T.; Sato, S.; Tanaka, R. Pet. Sci. Technol. 2004, 22, 901–914. (14) Stoyanov, S. R.; Gusarov, S.; Kovalenko, A. Mol. Simul. 2008, 34, 953–960. (15) Mackie, I. D.; DiLabio, G. A. J. Phys. Chem. A 2008, 112, 10968– 10976. (16) Mackie, I. D.; DiLabio, G. A. Phys. Chem. Chem. Phys. 2010, 12, 6092–6098.
(17) DiLabio, G. A.; Wolkow, R. A.; Johnson, E. R. J. Chem. Phys. 2005, 122, 044708. (18) Tsuzuki, S.; Honda, K.; Azumi, R. J. Am. Chem. Soc. 2002, 124, 12200–12209. (19) Vogiatzis, K. D.; Mavrandonakis, A.; Klopper, W.; Froudakis, G. E. ChemPhysChem 2009, 10, 374–383. (20) Mackie, I. D.; McClure, S. A.; DiLabio, G. A. J. Phys. Chem. A 2009, 113, 5476–5484. (21) Mackie, I. D.; DiLabio, G. A. Phys. Chem. Chem. Phys. 2010, in press.
i¼1
6469
Energy Fuels 2010, 24, 6468–6475
: DOI:10.1021/ef100797h
Mackie and DiLabio
fore, not carried out frequency analysis on structures determined without DCPs. For comparison purposes, calculations using the M06-2X29 density functional are given in the Supporting Information. We find good agreement to our DCP binding energies using this functional.
energy and structure, gives obvious benefits to both researchers without access to enormous computing resources but also to those who have an interest in larger molecules (larger than benzene dimer, for example). The DCP-method is easily applied with commonly available software and without the requirement of code modification. Moreover, the method allows for the use of relatively small basis sets such that less computing resources are required. A sample input for the Gaussian-0322 program is provided in Supporting Information. All calculations described herein were carried out using Gaussian-03, with B97123/6-31þG(d,p)-DCP or PBE24/ 6-31þG(d,p)-DCP. These combinations of functional and basis set use correcting potentials with ζ1 = 0.08, ζ2 = 0.12, c1 = -0.001438, c2 = 0.003475 (for B971) and ζ1 = 0.08, ζ2 = 0.12, c1 = -0.001550, c2 = 0.003300 (for PBE). For an illustration of how the DCP approach corrects the potential energy surface for the ring separation of thiophene dimer, see Figure 1 of ref 20. Treatment of dispersion was switched-off by removing the DCPs and reoptimizing the structures, starting from the optimized DCP geometries. It should be noted that the B971 functional does predict some long-range binding in some cases.25 B326LYP27 is the most commonly used density-functional in the literature, yet it is notoriously poor at describing van der Waals dimers.11 We therefore include calculations carried out using this method for comparison. No corrections for BSSE were made for these nondispersion-corrected approaches; therefore, binding energies reported are likely overestimated relative to those that would be obtained using a complete basis set. With the exception of dimers of 11 and 12 in Scheme 2 and the species in Scheme 3, frequency calculations for the DCPtype calculations confirmed each structure as a minimum on the potential energy surface. However, all reported binding energies do not include corrections for zero-point energies.28 The inability of conventional density-functionals to describe dispersive interactions is well established,11 we have, there-
Binding of Thiophene and Benzene to Oxygen-Functionalized Carbon A conformer search of possible dimers of thiophene and benzene bound to 1-8 was carried out using B971-DCP, with geometry optimization of up to 60 starting structures. Those that optimized to energy minimum, and without imaginary frequencies, were subsequently used for further analysis. The dimers presented, vide infra, represent a mixture of π-stacked, T-shaped, and coplanar structures. For our purposes, we label hydrogen-bonded dimers in which each monomer lies perpendicular to its neighbor as belonging to the T-shaped class, which will also include dimers where functional groups of one monomer interact with the π-region of the other such as with a CH 3 3 3 π interaction, for example. Of interest is the importance of each of these motifs to the contribution of intermolecular interaction and how properly describing noncovalent interactions influences the observed binding, both in terms of energy and structure. The data in Figure 1 show how dispersion influences the number of conformers of a particular structural motif for dimers involving monomer 1 with monomers 1-8. A histogram for dimers involving 2 shows very similar distributions to those in Figure 1 and is, therefore, relegated to the Supporting Information. Figures showing all DCP-optimized structures can also be found in the Supporting Information. The simple dimers, vide supra, ardently illustrate that the inclusion of dispersion is critical to structure and conformer
(22) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, Revision D.01; Gaussian, Inc.: Pittsburgh, PA, 2004. (23) Hamprecht, F. A.; Cohen, A. J.; Tozer, D. J.; Handy, N. C. J. Chem. Phys. 1998, 109, 6264–6271. (24) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865–3868. (25) Johnson, E. R.; DiLabio, G. A. Chem. Phys. Lett. 2006, 419, 333– 339. (26) Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652. (27) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785–789. (28) Frequency calculations are incredibly expensive computationally and, therefore, cannot be easily calculated for the largest species discussed. Inclusion of zero-point vibrational corrections reduces the predicted binding energies by ca. 15% for perylene dimer and ca. 40% for dimers of N,N0 -dimethyl-perylene bisimide. (29) M06-2X is a relatively new density-functional that has been developed with van der Waals forces in mind, see Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215–241.
Figure 1. Histogram showing the number of conformers of each structural motif predicted for dimers 1•X (X = 1-8, see Scheme 1) as predicted by dispersion-corrected and conventional densityfunctional theory methods.
6470
Energy Fuels 2010, 24, 6468–6475
: DOI:10.1021/ef100797h
Mackie and DiLabio
Table 1. Density-Functional Theory/6-31þG(d,p) Calculated Binding Energies (kcal/mol) for the Most Weakly and Most Strongly Bound Conformers of Each Dimera weakest bound dimer 1•1 1•2 2•2 1•3 2•3 1•4 2•4 1•5 2•5 1•6 2•6 1•7 2•7 1•8 2•8
B971-DCP b
0.7 2.3c 2.5d 2.6c 2.7c 3.1d 2.3b 2.3c 2.5c 2.3c 2.2c 2.2c 2.6c 2.2c 2.3b
B971 d
0.2 0.7d 1.0d 1.4d 1.4c 1.8d 1.7c 1.3c 1.4c 1.4c 1.3c 1.4c 1.6c 1.1c 1.5d
Table 2. Density-Functional Theory/6-31þG(d,p) Calculated Binding Energies (kcal/mol) and Monomer Center-of-Mass Distances (r, A˚) for Dimer Conformers of 9, 10, 11, and 12a
strongest bound
B3LYP b
0.1 0.4c 0.2d 0.6c 0.6c 0.8c 0.8c 0.4c 0.3d 0.5c 0.4c 1.4b 1.2b 0.6c 0.5d
B971-DCP
B971
c
3.0 3.2c 2.5d 4.9c 4.9c 4.4d 4.5d 5.4c 5.4c 2.8c 2.8c 4.8d 4.7d 4.5d 4.6d
c
1.9 1.9c 1.5c 4.3c 3.8c 3.0c 2.4d 4.2c 4.0c 1.7c 1.7c 3.2b 2.6c 3.3c 2.6c
with dispersion
B3LYP
without dispersion
B971
c
0.7 0.7c 0.4c 3.1c 2.6c 2.2c 1.4c 2.8c 2.6c 1.0c 0.7c 2.2b 1.4c 2.4c 1.6c
a The structural motif is also indicated. M06-2X calculated binding energies differ by an average of 0.22 kcal/mol (ca. 7.6%) from the B971DCP values, see Supporting Information. b Coplanar. c T-shaped. d π-stacked.
population for these noncovalently bound dimers. Figure 1 clearly displays that a significantly different distribution of conformer types (coplanar, T-shaped, π-stacked) are obtained when DCPs are used compared to when dispersion is not included in the calculation. When dispersion is included with DCPs, the number of π-stacked dimer conformers increases significantly compared to calculations without dispersion. For example, considering just the π-stacked conformers of 1•7, B971-DCP predicts the existence of 11 conformers and B971 without DCPs results in the prediction of only 4 conformers. B3LYP, which does very poorly for π-stacked systems, finds no conformers of this type. Of course, the energetics of binding are also greatly dependent on the inclusion of dispersion. The data in Table 1 shows the binding energies for the most weakly and most strongly bound dimers, independently determined from each method. Table 1 illustrates the stark differences in binding strength between DCP-derived structures and non-DCP dimers. Of the lowest energy dimers of each type, there exists a mean lowering of binding energy of ca. 47% (B971) and 75% (B3LYP) without dispersion, while for the maximum bound dimers the mean differences are ca. 34% and 61% for the same methods. Note too that, in several cases, the predicted dimers have different structural motifs. Such differences may have serious consequences. For example, dispersion binding is observed to play a role in orienting chemical reactants, thus affecting reaction probabilities and product distributions.30,31 If such discrepancies exist with simple (small) dimers such as these, then the problems are likely magnified with larger molecules such as those discussed in the following sections.
B971
B3LYP
dimer
BE
r
BE
r
BE
r
92-D2 92-D2d 92-C1 92-Ci 102-D2 102-C2 102-Ci (1) 102-Ci (2) 102-D2d 112-D2 112-C2 112-Ci 112-C2h 112-S4 112-C2v 122-C2 (1) 122-C2 (2) 122-Ci 122-Cs
15.4 14.6 14.5 12.6 24.5 24.1 21.2 21.1 18.2 28.8 28.1 24.6 23.4 18.0 15.7 24.8 21.8 19.4 14.7
3.4 3.5 3.8 5.0 3.4 3.5 3.8 5.0 3.4 3.4 3.4 3.8 3.7 3.4 3.7 3.4 3.5 5.0 5.1
3.0 3.2 3.3 3.4 6.3 6.1 6.0 5.5 6.3 9.1 9.0 7.3 6.3 4.5 3.0 7.2 4.4 7.3 4.7
4.1 4.1 4.8 5.5 3.9 3.9 5.1 5.5 4.0 4.0 4.0 5.0 4.6 4.0 4.2 4.1 4.6 6.5 5.8
-0.2 -0.2 0.6 0.5 0.6 0.6 0.8 1.3 0.4 2.2 2.0 1.8 1.3 0.3 -0.6 0.3 0.8 0.9 -1.7
3.5 4.9 6.2 6.2 4.4 4.3 5.5 6.6 4.4 4.3 4.4 5.5 5.4 4.3 5.0 4.9 4.8 7.7 7.2
a
The point group symmetry of the dimers is also indicated. Negative BEs represent unbound species. M06-2X calculated binding energies for the lowest energy dimers differ by an average of 0.52 kcal/mol (ca. 2.2%) from the B971-DCP values; see Supporting Information.
(11), and N,N0 -(1-hexylheptyl)-perylene bisimide (12). Dimers are henceforth labeled with subscript 2. Several symmetric conformers have been geometry optimized for each model, with the results gathered in Table 2. Given the size of these molecules, frequency calculations were only carried out on 92 and 102.28 Optimized coordinates are provided in the Supporting Information. The discrepancies between binding energies predicted with and without dispersion for the larger molecules listed in Table 2 are much more dramatic than those for the smaller species given in Table 1. B3LYP spectacularly fails to describe the structure or energetics involved in the dimers listed. Many are either unbound or have at most minimal binding, with many distances between monomers several A˚ngstroms greater than the more reliable B971-DCP structures. This is illustrated in 102, whereby the binding energies with B3LYP are up to 97% lower than those obtained using B971-DCP. B971 also performs poorly, a fact that only serves to underline the importance of dispersion to structure determination for all of the dimers studied. For example, in 112-D2, BEs of ca. 9.1 (B971) and 2.2 kcal/mol (B3LYP) pale in comparison to the DCP-value of 28.8 kcal/mol. The distance between monomer ring centers in these dimers also sees an increase by up to 0.9 A˚ when dispersion is not included. Similar differences are found for the other dimers under consideration. Perylene crystallizes32 in a herringbone fashion with a mean perpendicular distance between molecular planes of 3.46 A˚. In this gas-phase study, four different dimers have been considered. The lowest energy DCP-calculated dimer contains D2 symmetry (see Figure 2a) with a binding energy (BE) of 15.4 kcal/mol and its monomers 3.4 A˚ apart. The least bound 92 with DCPs is of Ci-symmetry (BE = 12.6 kcal/mol, r = 5.0 A˚),
Pericondensed Asphaltene Models As representative structures of pericondensed asphaltene models, we have investigated the dimerization of perylene (9), perylene bisimide (PBI, 10), N,N0 -dimethyl-perylene bisimide (30) Buelow, S.; Radhakrishnan, G.; Catanzarite, J.; Wittig, C. J. Chem. Phys. 1985, 83, 444–445. (31) DiLabio, G. A.; Piva, P. G.; Kruse, P.; Wolkow, R. A. J. Am. Chem. Soc. 2004, 126, 16048–16050.
(32) Camerman, A.; Trotter, J. Proc. R. Soc. London 1964, 279, 129– 146.
6471
Energy Fuels 2010, 24, 6468–6475
: DOI:10.1021/ef100797h
Mackie and DiLabio 34
by DFT-D. A further two energy minimum conformers contain Ci symmetry (Figure 3b,c) and have BE = 21.2 and 21.1 kcal/mol, respectively, which compare to 25.6 and 23.9 kcal/mol as calculated with DFT-D by Fink et al.33 We also found a conformer with C2 symmetry (BE = 24.1 kcal/mol) that is very similar in structure and energy to the lowest energy D2-symmetric conformer. A fifth conformer (D2d symmetry) is bound by 18.2 kcal/mol; see Figure 3d. Addition of methyl groups to the amino termini of 10 leads to N,N0 -dimethyl-perylene bisimide, 11. Six structures of 112 have been optimized (see Table 2). As is the case for 92 and 102, we find the lowest energy dimer of 11 to be of D2 symmetry (BE = 28.8 kcal/mol) with a C2-symmetric dimer ca. 0.8 kcal/mol higher in energy (BE = 28.0 kcal/mol). A S4-symmetry dimer, related to the D2d-symmetric 102 dimer, is bound by 18.0 kcal/mol, while a binding energy of 24.6 is found for a Ci-symmetry dimer. A structure close to Ci-symmetric has been the subject of recent publication by Vura-Weis et al.35 using M06-2X29/6-31þþG(d,p), but they found a BE of less than 18 kcal/mol, a difference of over 6.5 kcal/mol compared to the DCP-calculated value in this work. The M06-2X data may be explained in part by the neglect of symmetry and lack of geometry optimization in ref 35 and is reflective of the poor performance of this functional in dealing with perfectly stacked dimers.16 Dimers of 11 have the possibility of utilizing stabilizing CH 3 3 3 O interactions between monomers, as exemplified in 112-D2 and 112-C2, which have four such interactions, averaged to 2.58 and 2.55 A˚, respectively. The difference in BE compared to their 102 equivalents comes to 4.3 and 4.0 kcal/ mol, roughly 1 kcal/mol per CH 3 3 3 O interaction. For the Ci conformer of 112, only two CH 3 3 3 O interactions exist, both of which are 2.6 A˚ long, but they still serve to stabilize this dimer by 3.3 kcal/mol compared to its 102 counterpart. N,N0 -(1-Hexylheptyl)-perylene bisimide, 12, is a more complex asphaltene model with extensive alkyl-chain substituents attached to the PBI core. Before discussing the dimerization of this molecule, it is first important to note that intramolecular interactions within a monomer can have an effect on conformation stability. We therefore sampled a number of possible conformers (see the Supporting Information) and concluded that the most likely structure has the alkyl chains lying ca. perpendicular to the plane of the PBI core. This monomer was subsequently used to construct the initial structures of the dimers described hereafter, which were then subjected to minimization. Four conformers of 122 have been considered: two with C2 symmetry (Figure 4a,b), one with Ci symmetry (Figure 4c), and one with Cs symmetry (Figure 4d). Using DCPs, the difference in energy between 122-C2(1) and 122-Cs is ca. 10 kcal/mol. This is perhaps indicative of the increased importance of rotation of the monomers over displacement in dimer formation. The small difference in stability between the two C2-dimers is important since it could have implications on how these molecules would stack in a column. The lower value for 122-C2(1) suggests that trimerization, and further oligomerization, may occur sequentially from this dimer, though further calculations are required to answer this question. 122 shows similar BEs to 102 when DCPs are included. The same can be said for r. One exception to this is 122-Ci (Figure 4c), which is closest in structural motif to 102-Ci(1).
Figure 2. Perspective views of B971/6-31þG(d,p) optimized geometries for perylene dimers, 92. Structural symmetry, binding energies (kcal/mol), center of mass separation and distance between planes (A˚). (a) with DCPs, D2, 15.4, 3.4, 3.4; (b) without DCPs, Ci, 3.4, 3.9, and (c) with DCPs, C1, 14.5, 3.8, 3.5. The C atoms of the far perylene monomer are darkened, and all hydrogen atoms are omitted for clarity.
Figure 3. Perspective views of B971/6-31þG(d,p)-DCP optimized geometries for perylene bisimide dimers, 102. Structural symmetry, binding energies (kcal/mol), and center of mass separation between monomers (A˚). (a) D2, 24.5, 3.4, 3.4; (b) Ci, 21.2, 3.8, 3.5; (c) Ci, 21.1, 5.0, 3.4, and (d) D2d, 18.2, 3.4, 3.4. The C atoms of the far perylene monomer are darkened, and all hydrogen atoms are omitted for clarity. O atoms are colored red, and N atoms are colored royal blue.
yet this motif is the most strongly bound when dispersion is not included with the B971 functional (BE = 3.4 kcal/mol, r = 5.5 A˚, Figure 2b). The third motif considered possesses C1 symmetry (BE = 14.5 kcal/mol, r = 3.8 A˚, Figure 2c) and is the result of following the displacement of the largest negative frequency from dimers that have D2h (cofacial) or C2h symmetry. Similar binding (14.6 kcal/mol) is observed for a dimer with D2d symmetry (r = 3.5 A˚). Without DCPs, the C1-dimer is only bound by 3.3 kcal/mol, a drop of over 77% compared to the DCP-derived value, with the distance between monomers more than 1 A˚ longer. The influence of polar groups appended to a perylene core in 10-12 allows for further opportunities of favorable binding. Five energy minimum dimer structures of 10 have been determined, four of which are shown in Figure 3. With binding of ca. 25 kcal/mol, a D2-symmetric structure (Figure 3a) is found to be the most stable in energy. This compares to a value of 28.0 kcal/mol recently reported by Fink et al.,33 calculated (33) Fink, R. F.; Seibt, J.; Engel, V.; Renz, M.; Kaupp, M.; Lochbrunner, S.; Zhao, H.-M.; Pfister, J.; W€ urthner, F.; Engels, B. J. Am. Chem. Soc. 2008, 130, 12858–12859. (34) DFT-D refers to an alternative methodology for describing dispersion within DFT through the addition of explicit, attractive C6/ 6 R terms between all atomic pairs. For example, see Grimme, S. J. Comput. Chem. 2004, 25, 1463–1473.
(35) Vura-Weis, J.; Ratner, M. A.; Wasielewski, M. R. J. Am. Chem. Soc. 2010, 132, 1738–1739.
6472
Energy Fuels 2010, 24, 6468–6475
: DOI:10.1021/ef100797h
Mackie and DiLabio
Figure 4. Perspective side-views of B971/6-31þG(d,p)-DCP optimized geometries, symmetry, and binding energies (kcal/mol) for dimers of N,N0 -(1-hexylheptyl)-perylene bisimide, 12. (a) C2, 24.8; (b) C2, 21.8, (c) Ci, 19.4; (d) Cs, 14.6. The atom color scheme is the same as that in Figure 3.
Figure 5. Perspective views of open (top) and folded (bottom) structures of Maya asphaltene models. The folded structure is predicted to be ca. 9 kcal/mol lower in energy than the open structure. For clarity, the C atoms of component A (Scheme 3) are darkened, and all hydrogen are atoms removed. The atom coloring scheme is the same as that used in Figure 3, with S atoms colored yellow.
The large alkyl chains present in 12 reduce the ability of the π-regions of the PBI moieties to overlap, thus causing greater displacement along the short-axis of the PBI core, resulting in r = 5.0 A˚, compared to 3.8 A˚ in 102. The latter dimer is bound by 21.2 kcal/mol, but the former only has a BE of 19.4 kcal/mol. Without dispersion corrections, B971 and B3LYP favor a dimer that results from displacement of one monomer along both the short and long axes of the PBI core, rather than rotation about their centers, as found using DCPs. In the dimers where there is prominent displacement of one PBI compared to the other (e.g., 122-Ci), BEs of ca. 7.3 and 0.9 kcal/mol are calculated by B971 and B3LYP, respectively. This is favored over rotation of each monomer relative to the other, in 122-C2(1), by 0.2 (B971) and 0.6 (B3LYP) kcal/mol. In contrast, with DCPs, the C2 conformer is favored by 5.4 kcal/mol over its Ci counterpart. Larger r values are found for the non-DCP structures because the extensive alkyl chains repel each other to afford greater displacement between monomers. For example, in the Ci-case, r increases from the DCP value by ca. 1.5 and 2.7 A˚, using B971 and B3LYP, respectively. Similar trends are found for all of the remaining dimers of 12.
molecules involved is a challenge from a computing resource perspective since garnering accurate models is both expensive and time-consuming. Second, the components of asphaltenes are dependent upon the region, even oilwell, of collection. Finally, given the lack of agreement between experiment and the difficulty of data analysis and interpretation, it is of little surprise that the macro-structure of asphaltenes is much debated.3,36,37 Published descriptions of modeling attempts have primarily relied upon molecular mechanics38 or molecular dynamics,13 with starting structures based upon fitting an average to some experimental data such as spectroscopic result, elemental analysis, and molecular weight. One such study, carried out by Takanohashi et al.,13 describes efforts to model a Maya asphaltene. This model is made up of three sections, the components of which are shown in Scheme 3 and which has also been the subject of recent discussion.14 Clearly, aggregation of these Maya model components will be more complicated than the simpler perylene-type molecules discussed previously. However, the basic building blocks share common properties that will help determine appropriate (36) Strausz, O. P.; Peng, P.; Murgich, J. Energy Fuels 2002, 16, 809– 822. (37) (a) Ruiz-Morales, Y.; Mullins, O. C. Energy Fuels 2007, 21, 256– 265. (b) Alvarez-Ramirez, F.; Ramirez-Jaramillo, E.; Ruiz-Morales, Y. Energy Fuels 2006, 20, 195–204. (38) These are reviewed in Murgich, J. Mol. Simul. 2003, 29, 451–461.
Maya Asphaltene Model The difficulty in modeling asphaltenes is attributable to a number of fundamental factors. First, the extensive size of the 6473
Energy Fuels 2010, 24, 6468–6475
: DOI:10.1021/ef100797h
Mackie and DiLabio
binding motifs. Consequently, we have focused our attention to the extensive aromatic regions and how these will pack together to form aggregates. The large size of these π-regions and the noncovalent interactions between them will act as the primary driving force behind aggregation. Optimized geometries may be heavily reliant upon the starting structure used, with multiple local energy minima likely for such an aggregate. Equally, noncovalent interactions will play a key role in determining the relative energies of said local minima. Consequently, it is essential to sample enough configurational space to be confident of accurately determining the lowest energy structure. Therefore, on the basis of our findings for the other models, it occurred to us to check the folded structure of the Maya asphaltene for comparison to the layered structure proposed in ref 14. For this, we used PBE/ 6-31þG(d,p) with DCPs, thereby utilizing the computational efficiency of this density-functional.39 Our results are summarized in Figure 5, which shows calculated open and folded conformations for the Maya asphaltene model. Our calculations indicate that the folded conformation is ca. 9 kcal/mol more stable than the open structure. The folded structure benefits from stabilization gained from π-stacking in the central fragment (A in Scheme 3). Recent NMR work by Dettman40 also indicates that closed forms for asphaltenes, similar to that in Figure 5, are the dominant conformation. Without DCPs, PBE predicts the energy difference between the open and folded forms to be much smaller, with the latter being lower in energy by only 1.3 kcal/mol. On the other hand, the unphysical long-range repulsive behavior of B3LYP causes the folded Maya conformer to open significantly.
is the case with B3LYP), this would result in lack of dimer formation. Perylene bisimide (PBI) compounds have also recently been used by others42 to highlight the importance of noncovalent interactions. Nordgard et al.,42 via surface pressure-area isotherms, conclude that the presence of carboxylic acid groups is critical to the interfacial activity and film properties of PBIs. This is a result of the formation of a stable monolayer film as the acidic groups direct the molecules normal to the solvent surface. The binding of many further PBIs are reviewed by W€ urthner and co-workers, with steric interactions playing an important (de)stabilizing role.43 With this work, we have shown that neglect of dispersion seriously affects both the strength and motif of binding, even for very simple dimers. For example, considering just the π-stacked conformers of thiophene interacting with 2-pyrone, B971-DCP predicts the existence of 11 conformers but B971 without DCPs only results in the prediction of 4 conformers. B3LYP, which does very poorly for π-stacked systems, finds no conformers of this type. Similarly differential findings are found for thiophene binding to benzene, benzoic acid, cyclohexanone, phenol, furan, and 4-pyrone. Benzene bound to the same oxygen-functionalized carbon molecules shows comparable results. For larger systems, such as those more suited as asphaltene models, dispersion becomes even more critical. Indeed, dispersion interactions can lead to very stabilizing interactions. For example, the binding within a dimer of N,N0 -dimethylPBI (a pericondensed asphaltene model) can reach up to ca. 29 kcal/mol and primarily relies upon extensive π-π overlap between the perylene bisimide (PBI) cores. For the more complex molecule, N,N0 -(1-hexylheptyl)-PBI, the long alkyl chains limit the extent to which the two π-regions of the PBI monomers can overlap, and binding of a C2-symmetric dimer equates to ca. 25 kcal/mol. Without DCPs, this dimer is only bound by ca. 7 kcal/mol, while there is negligible binding when it is calculated using B3LYP. The distance between monomer ring centers also shows a marked increase with these two methods (by 0.7 and 1.5 A˚, respectively) compared to the B971-DCP value. Such discrepancies are also found for the other PBI dimers studied. Structures and energies of archipelago-type asphaltenes also suffer when dispersion is not included. This is neatly illustrated with a three-component model of Maya asphaltene. Calculated open and folded conformations for the Maya asphaltene model in this work clearly indicate that the latter is favored over the previously published13,14 motif by ca. 9 kcal/ mol. This finding is in line with recent NMR analysis of asphaltenes,40 which show that closed structures are preferentially obtained. This conclusion could not be arrived at without a method that treats dispersion. The availability of density-functional methods capable of treating dispersion11 means that large species, such as asphaltenes, can now be studied economically. Dispersion-correcting potentials (DCPs)12,15,16 are an easily implemented and convenient way of doing this since no software modification is required, and relatively small basis sets can accurately reproduce binding energies derived from vastly more expensive methodologies that cannot be applied to large species. With
Discussion and Summary The self-association of the compounds studied herein, and of similar compounds, results from a delicate balance between many different types of noncovalent interactions. The critical importance of dispersion in asphaltene model dimers has been demonstrated here. The term “dispersion-correcting potential” might lend one to think that dispersion is described with the DCP method to the detriment of other noncovalent interactions. This is not so since the underlying density-functional is itself capable of providing relatively accurate hydrogen-bonding and dipole-type BEs,11 and the introduction of DCPs has been shown to provide accurate binding for dimers containing hydrogen-bonds, dispersion interactions, dipole-type interactions, and mixtures of these.12,15,16,20 Rather, the DCP method gives accurate binding to the dispersion-dominated dimers where the bare density-functional cannot, while at the same time treating other noncovalent interactions to good effect. Such a method is, therefore, practical to systems that will contain a multitude of noncovalent interactions, including asphaltene models. Of course, aggregation will also be strongly influenced by the presence and type of solvent. Akbarzadeh et al.41 note that pyrene and dipyrenyl decane compounds possessing oxygen groups (carbonyl and hydroxyl) dimerize in o-chlorobenzene solvent. If the chosen theoretical method predicts too high a degree of repulsive forces between aggregate components (as (39) The PBE functional was used for this very large system because it is somewhat faster than B971 with only a minor tradeoff in the quality of predicted binding energies. (40) Dettman, H. personal communication. (41) Akbarzadeh, K.; Bressler, D. C.; Wang, J.; Gawrys, K. L.; Gray, M. R.; Kilpatrick, P. K.; Yarranton, H. W. Energy Fuels 2005, 19, 1268– 1271.
(42) Nordgard, E. L.; Landsem, E.; Sj€ oblom, J. Langmuir 2008, 24, 8742–8751. (43) Chen, Z.; Lohr, A.; Saha-M€ oller, C. R.; W€ urthner, F. Chem. Soc. Rev. 2009, 38, 564–584.
6474
Energy Fuels 2010, 24, 6468–6475
: DOI:10.1021/ef100797h
Mackie and DiLabio
this in mind, we aim to apply these methods to further our understanding of asphaltene structure and aggregation.
Supporting Information Available: Sample input file, optimized Cartesian coordinates, and figures for all structures. This material is available free of charge via the Internet at http:// pubs.acs.org.
Acknowledgment. We thank the Centre of Excellence for Integrated Nanotools (CEIN) and Prof. Pierre Boulanger at the University of Alberta, and WestGrid for computing resources, and the Program for Energy Research and Development (PERD) for funding. We also thank Dr. H. Dettman for helpful discussions on asphaltene structures and Drs. S. Stoyanov and A. Kovalenko for coordinates of Maya asphaltene.
Note Added after ASAP Publication. Labels in Table 2, and the title of Scheme 3 have been modified in the version of this paper published November 9, 2010. The correct version published November 12, 2010.
6475
