Revealing the Intermolecular Interactions of ... - ACS Publications

Feb 13, 2017 - ABSTRACT: Understanding the nature of non-covalent interactions (NCIs) between asphaltene molecules is not only theoretically interesti...
1 downloads 0 Views 1MB Size
Subscriber access provided by University of Colorado Boulder

Article

Revealing the Intermolecular Interactions of Asphaltene Dimers by Quantum Chemical Calculations Huanjiang Wang, Haiyan Xu, Weihong Jia, Juan Liu, and Sili Ren Energy Fuels, Just Accepted Manuscript • DOI: 10.1021/acs.energyfuels.6b02738 • Publication Date (Web): 13 Feb 2017 Downloaded from http://pubs.acs.org on February 14, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Energy & Fuels is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Revealing the Intermolecular Interactions of Asphaltene Dimers by Quantum Chemical Calculations Huanjiang Wang,†,‡ Haiyan Xu,†,‡ Weihong Jia,†,‡ Juan Liu,†,‡ Sili Ren†* †

State Key laboratory of Solid Lubrication, Lanzhou Institute of Chemical Physics, Chinese Academy of Sciences, Lanzhou, P. R. China. ‡

University of Chinese Academy of Sciences, Beijing, P. R. China.

Corresponding Author: *Email: [email protected] Address:NO.18, Tianshui Middle Road, Lanzhou City, Gansu Provience, 730000, P.R.China. Tel: +86-931-496-8051; Fax: +86-931-827-708

ABSTRACT: Understanding the nature of noncovalent interactions between asphaltene molecules is not only theoretically interesting but also important for practical application. We performed quantum chemical calculations to reveal the configuration feature and intermolecular interaction characteristics of asphaltene dimes using three representative asphaltene model compounds and their derivatives. The frontier molecular orbitals and electrostatic potential map (ESP) of the model asphaltenes were analyzed to reveal the nature of interaction between the asphaltene monomers. Calculation of binding energies

1

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

indicates that the stability of asphaltene dimers not only depend on the number of aromatic rings but also rely on the presence of heteroatoms in the aromatic core and aliphatic side chains, which could change the electrostatic charge distribution on the molecular vdW surface. In addition, noncovalent interactions (NCI) and the natural bond order (NBO) analysis method was used to identify the interactions that promote the formation of asphaltene dimers. It was found that the reduced density gradient (RDG) isosurfaces could clearly reveal the type of interactions between two asphaltene monomers in their dimers. The results indicate that various interactions possess either electrostatic or dispersive nature including hydrogen bonding, θ−θ, θ−π, and π−π interactions, among which the π−π stacking interaction is believed the major driving force for asphaltenes aggregation. 1. INTRODUCTION Asphaltenes are known as the largest, heaviest, and surface-active fraction of crude oil and bitumen. They are insoluble in n-alkanes, such as n-heptane and n-pentane, but soluble in aromatic solvents like toluene and benzene.1-2 The complexity of the molecular structure of asphaltenes brings much trouble in understanding their properties. In recent years, a variety of experimental techniques including the Small-Angle X-ray, Neutron Scattering fluorescence depolarization technique, fluorescence correlation spectroscopy and mass spectrometry are employed to analyze the structure of asphaltenes. 2-4 Now, it is generally accepted that the asphaltene molecules have an average molecular weight of 750 g/mol, with one polycyclic aromatic core consisting of 4−10 aromatic rings and aliphatic chains about 3 to 7 carbons in length. 2-3, 5 In addition, a few heteroatom of N, O and S are strewn in the molecules.

2

ACS Paragon Plus Environment

Page 2 of 23

Page 3 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

As the most complex component in crude oil, asphaltenes always brings various problems that trouble the oil industry. For examples, it could not only increase the fluid viscosity of the oil, but also contribute to produce and stabilize the oil−water emulsion. Moreover, it often gives rise to solid deposits blocking the pipeline. All these problems are closely related to the high tendency of self-association of the asphaltenes. 6-8 Over the past decades, great efforts have been made to explore the aggregation mechanism and strength of asphaltenes. Mullinset et al.

9

investigated the self-assembled process of

asphaltenes and assumed that the asphaltenes self-aggregated in a stepwise fashion: first aggregating to produce “‘nanocrystals” by π−π stacking of the polyaromatic core, following with self-association of the nanocrystals to form larger colloidal nanoaggregates. However, experimental research has inevitable shortcomings in studying the aggregation of asphaltenes at nanoscale or molecular levels. In contrast, the computer simulation techniques could provide more useful information on the interaction between asphaltene molecules. Pacheco-Sánchez et al. 10 investigated the asphaltenes’ aggregation in vacuum and found that the self-association of asphaltenes has at least three different geometric piling structures of face-to-face, π-offset and T-shaped stacking modes. Mohammad et al.

11

elucidated the asphaltenes aggregation in solvents by molecular

dynamics (MD) simulations. The results suggested that π-π interaction was the most significant force for the formation of the dimer aggregates. To study the stability of the asphaltene dimers, Alvarez-Ramirez et al.12 studied the interaction energy between asphaltene-like molecules by using molecular mechanics (MM) and quantum mechanics (QM). It was found that the interaction energy between asphaltene molecules in vacuum to be about −12 to −15 kcal/mol. They also discussed the possible conformation and their

3

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 23

energy for a total of more than 3,000 asphaltene dimers using quantum chemistry and zindo transitions.13 However, as previous reports,14 the classical MM method is difficult to correctly describe the interaction energy for the asphaltene dimers. The standard density functions used in QM are also failed to describe the long-range dispersion interactions.15 It is worthy of pointing out that accurate description of the noncovalent interactions in molecular systems is a prerequisite for well understanding the mechanism and strength of the asphaltenes aggregation. Otherwise, any attempt to model the asphaltenes aggregation would be incomplete or inaccurate. Recently, with development of the dispersion-corrected density functional (DFT-D), 16-19

it is provided a unique tool for the investigation and analysis of a wide range of

complex aromatic systems. However, only a few works deal with the intermolecular interaction of asphaltene dimes using the DFT-D method. In this study, we employ DFT-D3 method developed by Grimme et al,19-20 to study the structure-aggregation relationship of asphaltenes using three representative asphaltene molecular models and corresponding analogues. In addition, to reveal the interactions between the asphaltenes, the visualization scheme based on the noncovalent interactions (NCI) and natural bond order (NBO) analysis method

21-23

was applied to identify the noncovalent interactions

that promote the formation of asphaltene dimers. Our results will provide a clear picture for understanding the asphaltenes association and precipitation. 2. COMPUTATIONAL MTHODS All calculations were performed with the ORCA-3.0.3 program package.24 The geometries of the asphaltene molecular models and the corresponding dimers were optimized with the dispersion corrected BLYP-D3 functional

4

ACS Paragon Plus Environment

25

together with the

Page 5 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

double-ζ Def2-SVP

26

and verified by harmonic vibrational frequency computations.

Basis set superposition error (BSSE) was corrected with the geometrical counterpoise (gCP) method. 27 This combination has been tested and proved to be an effective method to provide accurate energies and structures of large systems. computational time, resolution-of-identity (RI) technique

32

28-31

In order to reduce the

was used to accelerate the

self-consistent field (SCF) process and second-order perturbative correction part. For comparison, the single point energy of the monomers and the corresponding dimers were also calculated by BLYP-NL and BLYP-D3 functionals in conjunction with Def2-TZVP basis set (listed in table s1). 33 To take into account of the solvent effects, the geometries of all the monomers and the corresponding dimers were also optimized in solvent of toluene and/or heptane using the conductor-like screening model (COSMO).34 In order to obtain a reliable signifier for the strength of interactions in asphaltene dimers, the binding energy (∆E) during dimerization was calculated by eq (1):

∆E = Edimer − 2 Emonomer ( non −bonding )

(1)

Where Emonomer(non-bonding) is the energy of isolated asphaltene monomer at its equilibrium geometry, and Edimer is the asphaltene dimer total energy. Binding energy (∆E) is indicative of the thermodynamic stability of the system. A negative value of ∆E indicates that the asphaltenes dimer could spontaneously formed and maintained stable. In addition, the NCI analysis was used to detect the interactions between the asphaltenes in their dimers. NCI is a topological tool to visualize noncovalent interactions by searching critical points in the electron density topology through the reduced density gradient (RDG). A singularity in the RDG corresponds to an electron density critical point. The eigenvalue of the electron density Hessian, along with the magnitude of the 5

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 23

density, was used to determine the strength of a particular nonbonding interaction, such as the van der Waals force, hydrogen bonds, and steric repulsion. analysis was carried out with the Multiwfn program

36

21-22, 35

The NCI

utilizing either its default

parameters or a uniform spatial grid with a step of 0.1 a.u. In order to focus on the weak noncovalent interactions, only the regions of 0.05 a.u. were analyzed. The NCI regions were visualized using Visual Molecular Dynamics software version 1.9.1. 37 NBO analysis was performed using the NBO 3.1 program as implemented in the Gaussian09 package at the BLYP/Def2-SVP level to obtain the second-order perturbative energy (E(2)), which has been used as a valuable tool to evaluate the donor-acceptor interactions. It is defined as:

E ( 2 ) = ∆Eij = qij

( F i , j )2

(2)

E j −E i

Where qi represents the occupancy of donor orbital, Ei and Ej are diagonal elements (orbital energies), and Fi,j is the off diagonal NBO Fock matrix element. 38 Table 1. Structural parameters of various model asphaltene molecules

a1 a2 b1 b2 c1 c2

Dipole moment (debye) 0.48027 0.86212 1.35755 1.54933 2.37834 1.51197

No. of aromatic carbons 19 20 26 26 25 25

No. of aliphatic carbons 4 4 21 22 32 32

Heteroatoms N N, S N N N

3. RESULTS AND DISCUSSION

6

ACS Paragon Plus Environment

No. of H/C hydrogen

Mw (g/mol)

19 20 55 57 71 63

309 308 666 648 770 762

0.83 0.83 1.17 1.19 1.25 1.11

Page 7 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

In this study, three representative asphaltene molecular models and their analogues were used to study the structure−nanoaggregation relationship. As shown in Figure 1, the asphaltene models of a1, b1 and c1 were obtained according to the Indonesian bituminous coal sample, UG8 Kuwait crude oil sample, and the Yen-Mullins model, respectively.39-40 The corresponding analogues with either heteroatoms in the aromatic core and/or the aliphatic side chains or different number of aromatic rings were used to study the role of heteroatoms and aromatic rings on asphaltenes aggregation. The asphaltene models used here were obtained by optical measurements and further analyzed by molecular orbital calculations, which were also used in early theoretical calculations.41-43 The molecular structural parameters of various model molecules were calculated and tabulated in Table 1. It was shown that the dipole moment of various model asphaltenes was different, which was attributed to the fact that the heteroatoms and the number of aromatic rings altered the distributions of electron density of various asphaltenes.

Figure 1. Molecular structures of asphaltene molecules used in this study. Color coding of various elements: white – H, cyan – C, blue – N and yellow – S. 7

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 23

3.1. Electronic characteristics of various asphaltene monomers

Figure 2. Illustration of the frontier molecular orbitals of various model asphaltene molecules at the BLYP-D3/Def2-SVP level.

Figure 2 shows the orbital structure of highest occupied molecular orbitals (HOMOs) and the lowest unoccupied molecular orbitals (LUMOs) of all asphaltene molecules. It was found that the frontier molecular orbitals of all the asphaltene molecules possess π orbital features spreading over the whole aromatic cores, indicating that all asphaltenes prefer to stack by π−π interaction. The donor-acceptor HOMO-LUMO gap (∆ε) as an good index has been used to study the molecular size and structure of asphaltenes,41 or to evaluate the π−stacking interactions between two complexes which have the similar frontier orbitals in structure.44 As seen from Table 2, the values of ∆ε for all the model asphaltene molecules are small and well fall in the experimental and theoretical range (1.9065−3.0981 eV),

8

ACS Paragon Plus Environment

39, 45

Page 9 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

confirming that the calculation methodology used in this study is reliable. By comparing the calculated ∆ε of various asphaltenes with their corresponding analogues, it was found that the effect of heteroatoms on ∆ε was almost negligible, which is in good agreement with the previous report.13 Table 2. Energy levels of LUMO, HOMO and the LUMO/HOMO gap (∆ε) for various model asphaltene molecules and their analogues. LUMO(eV) HOMO(eV) ∆ε(eV)

a1 −2.06 −4.75 2.69

a2 −2.01 −4.68 2.67

b1 −1.88 −4.35 2.47

b2 −1.89 −4.35 2.46

c1 −2.27 −4.19 1.92

c2 −2.31 −4.43 2.12

The electrostatic potential map (ESP) not only conveys a variety of information including the polarity, electronegativity, and the nature of the bonds, but also is widely used to study and predict the noncovalent interactions.46 In-depth investigation of asphaltenes, ESP is helpful for understanding the nature of interaction between asphaltene monomers. Figure 3 plots the ESP-mapped vdW surface along with the surface extrema, in which the red and blue colors represent the positive and negative potential regions, respectively. It is found that the polycyclic aromatic hydrocarbon (PAH) moiety (aromatic core) and non-PAH moiety (aliphatic and/or naphthenic hydrocarbon) of asphaltenes exhibit quite different electrostatic potential characteristics. The ESP values over the PAH moiety of various model asphaltene molecules and the analogues are moderately negative (−15 to −20 kcal/mol), reflecting that the aromatic core is rich of the π-electron cloud. In contrast, the PAH hydrogens is positively charged because the elements of C and N have larger electronegativity. It is also observed that the global surface extrema is close related to the position of the heteroatoms on the molecules. For

9

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

example, the global surface minima was found to be around the nitrogen that on the aromatic rings (Figure 3 a1, c1, c2), while the surface maxima was located at hydrogen near the nitrogen that on the pyrrolic rings (Figure 3b1, b2). Such result is attributed to either the lone pair electrons or the electronegativity of the N atom.

Figure 3. ESP-mapped molecular vdW surface for various model asphaltenes and their analogues. Surface local minima and maxima of ESP are presented in cyanand orange dots, respectively.

3.2. Configurations and binding energy of the asphaltene dimers It is known that the dimers of asphaltene molecules have three geometric stacking modes of face-to-face, π-offset and T-shaped.10 In this study, the presence of heteroatoms and side chains increases the variety of configurations for the association of two asphaltene monomers. Therefore, the molecular dynamics simulations were used to get a proper initial configuration of the asphaltene dimers. Firstly, the initial molecular configurations of the asphaltene dimers were generated by using the Packmol program. 10

ACS Paragon Plus Environment

Page 10 of 23

Page 11 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Then, the coordinates of molecular structures were then exported and submitted to the LAMMPS to obtain the formation of dimers by molecular dynamics (20 ns) at 300K (MD details are given in the Supporting Information). When the systems reached their equilibrium, we found that the asphaltene molecules aggregated spontaneously and formed dimer aggregates with the aromatic ring sheets face-to-face stacked together. The obtained asphalene dimers were then used as the initial configuration for quantum chemical calculations. Table 3. Binding energy (∆E, kcal/mol) and π−π stacking interaction distances (Å) for optimized asphaltene dimers.

Da1 Da2 Db1 Db2 Dc1 Dc2

Vacuum ∆E d 3.3969 −23.02 −22.68 3.4248 −43.90 3.4071 −46.65 3.4033 −40.11 3.5606 −49.08 3.5037

Heptane ∆E d 3.3974 −22.03 −21.7 3.4212 −40.83 3.4094 −43.26 3.4034 −38.26 3.5607 −47.26 3.5035

Toluene ∆E d 3.3978 −21.72 −21.39 3.4210 −39.88 3.4090 −42.19 3.4033 −37.65 3.5607 −46.68 3.5032

The structures of various monomers for the six model asphaltene molecules and their dimers were optimized in vacumm and verified by frequency calculations using the BLYP-D3/Def2-SVP method (Figure 4). The binding energy (ΔE) of dimers and vertical distance (d) between the mean planes of two parallel fused aromatic cores of asphaltenes are listed in Table 3. In order to evaluate the reliability of the computational methodology, the single point energy of the monomers and the corresponding dimers were also calculated by BLYP-NL and BLYP-D3 functionals in conjunction with a bigger basis set Def2-TZVP (listed in table s1). As shown in table 3 and table s1, the difference of the binding energies obtained in the two methods for each dimer of various model 11

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

asphaltenes is small and less than 5 kcal/mol, which validated that using a small basis set (Def2-SVP) could also obtain reliable results. To investigate the effect of heteroatoms on asphaltenes aggregation, two kinds of asphaltene molecules were considered, i.e. the a1,2 and b1,2 as shown in Figure 1. Comparing with the molecule of a1, pyridine nitrogen atom was introduced in the aromatic ring of a2, while sulfur atom was added in the lateral chains of b1 and that was not contained in b2. As seen in Table 1, the introduction of heteroatoms changes the polarity of asphaltene molecules and thus their dipole moment, which enabled us to examine the effect of heteroatoms lied either in the aromatic or in the aliphatic group on the aggregation behavior of asphaltenes.

Figure 4. The structures of asphaltene dimers optimized in vacuum at the BLYP-D3 /Def2-SVP level.

12

ACS Paragon Plus Environment

Page 12 of 23

Page 13 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

For various dimers shown in Figure 4, the Da1 and Da2 have the most stable configuration of π-offset mode, which possesses the largest overlap between the planar surfaces of asphaltenes and hence the strongest vdW attraction. Although there is the electrostatic repulsion between the surfaces due to the electron clouds, the repulsion is much smaller than the vdW attraction. Comparing with the Da2, the Da1 is easier to be formed because the heteroatom N on the aromatic core has a positive role in stabilizing the asphaltene dimers. Nitrogen has larger electronegativity than carbon, which can enhance the π-electron cloud of the ring and thus contributes to the π stacking interaction. As reported in the literature, 47 introduction of the side chains in the continental-type asphaltenes caused a steric interference for the π−π interaction, while the θ−π and/or θ−θ interactions arising from the long side chains might compensate for the decrease of the π−π force. Based on our calculations, the binding energy for Db1 and Db2 in vacuum are −43.90 and −46.65 kcal/mol, respectively. For both of the Db1 and Db2, the π-offset stacked arrangement is still preserved in their dimer’s structure. However, they have different molecular orientation. Because the global surface minima of b1 (ESP) located at the S atom in the side aliphatic chains, it can form noncovalent (Coulomb and/or vdW) interactions with hydrogens on the neighboring asphaltene. Therefore, the heteroatoms attached to the side aliphatic chains play significant role on the stacked modes of the asphaltene dimers, although the asphaltenes aggregation behavior is highly dependent on the π−π interactions between the polyaromatic cores. To further examine the effect of the number of aromatic rings on asphaltenes aggregation, the model asphaltene c1 changed to c2 by converting three cyclic aliphatic rings to aromatic rings. The calculation result shows that the binding energy of

13

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 23

dimerization for Dc1 and Dc2 are −40.11 and −49.08 kJ/mol, respectively. Such result indicates that the stability of asphaltene dimers substantially increased with the increase of the number of aromatic rings, which enhanced the π−π interactions between the asphaltenes. Besides of the chemical structure of the asphaltene molecule itself, the physical environment of medium also play important roles on the asphaltenes aggregation. Previous work showed that the asphaltenes even at low concentrations prefered to self-associate in crude oil and solvent such as toluene.9,

49

In order to interpret the

behaviors of asphaltenes aggregation in organic solvents, toluene and heptane were used. Toluene is known as a good solvent for asphaltenes, while heptane is the standard precipitant for asphaltenes. In order to interpret the behaviors of asphaltenes aggregation in different organic solvents, toluene and heptane were used. Toluene is known as a good solvent for asphaltenes, while heptane is the standard precipitant for asphaltenes. However, the results of our computations showed that the binding energy for a1, b1 and c1 and their analogues in heptane is only slightly higher than that in toluene (Table 3), reflecting that the asphaltene molecules have the tendency to form the dimers in both solvents of heptane and toluene. In order to verify the result about the dimerization of asphaltenes implicit solvent, the MD simulations were also employed (details are given in the Supporting Information). The results from the explicit-solvent simulations are given in table 1s, which are consist with the findings obtained by quantum chemical calculations. Such result indicates that the asphaltenes are in the form of small aggregates with a few molecules even in a good solvent, which is also well consistent with that of the reports.9, 50

14

ACS Paragon Plus Environment

Page 15 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

3.3. NCIs analysis To further identify the characteristics of interactions between asphaltenes in their dimers, the NCI approach was applied to detect the non-covalent interactions, which is based on a 2D plot of the electron density(ρ) and its first derivatives ( s( ρ ) =

1 2 1/3

2(3π )

∇ρ

ρ 4/3

). Where ∇ is gradient operator, and ∇ρ is the norm of electron

density gradient.18 In the plots of s versus ρ, the noncovalent interactions could be reflected by the characteristic peaks appeared at low electron density, which is due to the annihilation of the density gradient at these points. However, both attractive and repulsive interactions such as the hydrogen bonds and steric repulsion may appear in the same region of the density/reduced gradient space. To distinguish these interactions, the second derivative of the density was considered. The sign of the second eigenvalue (λ2) of the electron-density Hessian matrix can be used to distinguish the bonded interactions (λ20) interactions. Analysis of the sign of λ2 thus helps to discern the different types of noncovalent interactions. Therefore, various interactions can be further differentiated from the color of gradient isosurfaces, which are plotted by the RGB (red-green-blue) scale according to the corresponding values of sign(λ2)ρ in real space. The red color of the isosurfaces indicates the repulsive interactions, while the blue corresponds to the attractive interactions and green is for very weak vdW interactions. Figure 5a shows the diagrams of s(sign(λ2)ρ) for various model asphaltene monomers and corresponding dimers. The overlapping region represents the weak interactions in the monomers, while the new generated spikes (black) refer to the noncovalent interaction in asphaltene dimers. It seems that there is no obvious difference among the diagrams for various model asphaltene molecules and their dimers. 15

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nevertheless, the interactions between asphaltenes in the dimer could be clarified with the aid of RDG isosurfaces, which provide a rich visualization of the noncovalent interactions as broad regions of real space.

Figure 5 a) diagrams of s(sign(λ2)ρ) for various model asphaltene molecules (red color) and corresponding asphaltene dimers (black color). b) RDG isosurfaces for asphaltene dimers. The isosurface at s=0.6 au was colored according to a RGB scheme over the range −0.03 < sign(λ2)ρ < 0.02 au.

In this work, the target is to investigate the aggregation behavior of the asphaltenes. Therefore, NCI visualization was carried out with emphasis on the intermolecular interactions in an asphaltene dimer and the intramolecular interaction within an asphaltene molecule was screened. As shown in Figure 5b, two attractive isosurfaces were observed in the parallel-displaced structure of Da1, in which a continuous isosurface 16

ACS Paragon Plus Environment

Page 16 of 23

Page 17 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

appears in the overlapping region between the monomers and a disk appears between the H(85) and C(22), corresponding to the π−π stacking interaction and a weak CH…C dispersion interaction, respectively. Corresponding to the s(sign(λ2)ρ) plots (Figure 5a), peaks at the ρ values over a range of −0.006 to 0.006 a.u. are associated with the π−π interactions, while peak at about ρ=−0.012 corresponds to the CH…C dispersion interaction. As to the dimer of Da2, only a continuous isosurface appears in the overlapped region between the monomers (Figure 5b), corresponding to the π−π stacking interactions. The CH…C dispersion interaction was not found in Da2 and the binding energy is therefore slightly smaller than that of Da1 (Table 3). It is further noticed that the isosurfaces for the Db1 and Db2 are significantly different (Figure 5b). Db1 has three type attractive isosurfaces, i.e. a continuous isosurface appeared in the overlapped region between the monomers corresponds to the π−π stacking interaction, a discontinuity isosurface appeared between the side aliphatic chains and the aromatic rings corresponds to the θ−π interactions, and a series of small disks appeared between the S and H in different asphaltene molecules correspond to the electrostatic interaction. In the Db2, two type attractive isosurfaces appear in the overlapping region between the monomers, corresponding to π−π stacking interaction and weaker θ−θ interaction, respectively. Although the dimers of Db1 and Db2 are stabilized by different interactions, both of them possess the π-offset stacked configuration. However, the molecular orientations in the π-offset stacked configuration are different. It is believed that the π−π stacking interaction between the poly-aromatic cores determines the aggregate configuration, while the heteroatom in the side chain has an important role on the molecular orientation. 17

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

As shown by the NCI isosurfaces of the Dc1 and Dc2 (Figure 5b), the two dimers are stabilized by π−π, θ−π, and θ−θ interactions. However, for the Dc1, no isosurface appears on the top of the aliphatic rings, suggesting that the presence of aliphatic rings might disturb the formation of π−π stacking interaction due to the steric hindrance. Comparing with that of Dc1, the isosurface area of the overlapped region for the Dc2 increased with the increase of the number of aromatic rings in the asphaltene melecules. It is believed that the π−π interaction between the aromatic cores is the major driving force for the asphaltenes aggregation. 3.4. Orbital interaction between asphaltene monomers

Figure 6. Major second-order perturbation interaction energies for asphaltene dimers obtained from the NBO analysis (isosurfaces = 0.02).

18

ACS Paragon Plus Environment

Page 18 of 23

Page 19 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

It is believed that the π-π interaction process between two asphaltene molecules generally involves two steps. Firstly, the dispersion force plays a dominant role to bring the two monomers together. Then, with the two monomers approaching each other, Pauli repulsion will be produced and prevent them from getting closer. 48 Once the equilibrium distance reached, electrostatic interaction, Pauli repulsion, dispersion and intermolecular orbital interaction will work together. To study the intermolecular orbital interaction between the asphaltene monomers in their dimer aggregation, NBO analysis was performed, which has been employed to quantify the noncovalent interactions by analyzing the second-order orbital interaction energies E(2).49 Table 4. The Second-Order Interaction Energies (E(2), kcal/mol) for the sphaltene dimers a. dimers Da1 Da2 Db1 Db2 Dc1 Dc2

Donor (i) BD (C22-H36) π (C61-C62) π (C54-C64) π (C52-C53) π (C29-C30) LP (S174) π (C116-C117) π (C116-C117) π (C139-C140) LP (N148) π (C32-C33) LP(N140)

Acceptor (j) BD*(C67-H85) π*(C1-C18) LV(C4) π*(C5-C6) LV(C106) BD*(C85-H92) BD*(C16-H80) RY(H80) BD*(C53-H60) BD*(C69-H74) BD*(C209-H212) BD*(C60-H64)

E(2) 0.38 0.36 1.29 0.58 1.43 1.13 1.3 0.83 0.96 0.64 1.54 1.19

a

BD, BD*, LV, LP and RY stand for bonding, anti-bonding, lone vacancy, lone pair and Rydberg orbital, respectively.

The NBO analysis was carried out for all the asphaltene dimers. The results of E(2) are listed in table 4 and the corresponding NBO plots are illustrated in Figure 6. It is found that the main orbital overlap between asphaltene monomers are mainly in

19

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

π(C-C)/BD*(C-H), π(C-C)/LV(C), LP(N,S)/BD*(C-H), BD(C-H)/BD*(C-H) and π(C-C)/π*(C-C). It is noticed that the E(2) value are relative small, indicating that the intermolecular orbital interaction plays insignificant role on the aggregation of the alsphaltene dimers. In addition, the NBO analysis confirmed that hydrogen-bond was existed between the heteroatoms and hydrogen, which indicates electron transfer from a lone pair donor (N, S) to an antibonding orbital of the acceptor BD*(C-H). By comparing the NCI and NBO analyses, it was concluded that the interactions between asphaltene monomers in their dimers have both the electrostatic and dispersive nature. 4. CONCLUSION In this study, the dispersion-corrected density function BLYP-D3 was employed to study the configuration and intermolecular interaction behaviors of asphaltene dimes using a serials of representative model asphaltene molecules. The stability of the homodimers of these model compounds was evaluated in terms of the binding energy as well as the π−π interaction distances. The results suggested that the binding energies and optimized configuration of asphaltene dimers not only depend on the number of aromatic rings but also rely on the presence of heteroatoms in the aromatic core and aliphatic side chains, which can change the electrostatic charge distribution on the molecular vdW surface. Based on the NCI and NBO analyses, it is indicated that the driving forces promoting the formation of the asphaltene dimers possess either electrostatic or dispersive nature. For examples, the hydrogen bonding between hydrogen and heteroatoms in the lateral chain belongs to the electrostatic interaction, whereas the dispersion forces include θ−θ interaction between the long side chains, θ−π interaction between aliphatic side chains and aromatic rings, and π−π interaction between the PAH

20

ACS Paragon Plus Environment

Page 20 of 23

Page 21 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

moiety. It is believed that the π−π stacking interaction is the major driving force for the asphaltenes aggregation.

Supporting Information Details of molecular dynamics simulations (part1), dimerization of asphaltenes in organic solvents (part 2), π−π stacking distance for various asphaltene dimers in heptane and toluene (Table s1), binding energies (kcal/mol) of the asphaltene dimers calculated by various theoretical methods (Table s2) and coordinates of various asphaltene monomers and the corresponding dimmers in vacuum (Tables s3).

ACKNOWLEDGMENTS The authors thank the fund support from the National Science Foundation of China (Grant No. 51574217; 51374195).

REFERENCES 1. Carnahan, N. F.; Quintero, L.; Pfund, D. M.; Fulton, J. L.; Smith, R. D.; Capel, M.; Leontaritis, K., A Small-Angle X-Ray-Scattering Study of the Effect of Pressure on the Aggregation of Asphaltene Fractions in Petroleum Fluids under near-Critical Solvent Conditions. Langmuir : the ACS journal of surfaces and colloids 1993, 9 (8), 2035-2044. 2. Groenzin, H.; Mullins, O. C., Asphaltene molecular size and structure. The Journal of Physical Chemistry A 1999, 103 (50), 11237-11245. 3. Mullins, O. C.; Sabbah, H.; Eyssautier, J. l.; Pomerantz, A. E.; Barré, L.; Andrews, A. B.; Ruiz-Morales, Y.; Mostowfi, F.; McFarlane, R.; Goual, L., Advances in asphaltene science and the Yen–Mullins model. Energ Fuel 2012, 26 (7), 3986-4003. 4. Sato, S.; Takanohashi, T.; Tanaka, R., Molecular weight calibration of asphaltenes using gel permeation chromatography/mass spectrometry. Energ Fuel 2005, 19 (5), 1991-1994. 5. Sjoblom, J.; Simon, S.; Xu, Z., Model molecules mimicking asphaltenes. Advances in colloid and interface science 2015, 218, 1-16. 21

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

6. Tchoukov, P.; Yang, F.; Xu, Z.; Dabros, T.; Czarnecki, J.; Sjoblom, J., Role of asphaltenes in stabilizing thin liquid emulsion films. Langmuir : the ACS journal of surfaces and colloids 2014, 30 (11), 3024-33. 7. Wang, M.; Hao, Y.; Islam, R.; Chen, C.-C., Aggregation Thermodynamics for Asphaltene Precipitation. AIChE Journal 2016, 64 (4), 1254-1263. 8. Liu, J.; Zhao, Y.; Ren, S., Molecular Dynamics Simulation of Self-Aggregation of Asphaltenes at an Oil/Water Interface: Formation and Destruction of the Asphaltene Protective Film. Energ Fuel 2015, 29, 1233−1242. 9. Mostowfi, F.; Indo, K.; Mullins, O. C.; McFarlane, R., Asphaltene Nanoaggregates Studied by CentrifugaHon†. Energ Fuel 2008, 23 (3), 1194-1200. 10. Pacheco-Sanchez, J. H.; Zaragoza, I. P.; Martinez-Magadan, J. M., Asphaltene aggregation under vacuum at different temperatures by molecular dynamics. Energ Fuel 2003, 17 (5), 1346-1355. 11. Sedghi, M.; Goual, L.; Welch, W.; Kubelka, J., Effect of Asphaltene Structure on Association and Aggregation Using Molecular Dynamics. Journal of Physical Chemistry B 2013, 117 (18), 5765-5776. 12. Alvarez-Ramirez, F.; Ramirez-Jaramillo, E.; Ruiz-Morales, Y., Calculation of the interaction potential curve between asphaltene-asphaltene, asphaltene-resin, and resin-resin systems using density functional theory. Energ Fuel 2006, 20 (1), 195-204. 13. Alvarez-Ramírez, F.; Ruiz-Morales, Y., Island versus Archipelago Architecture for Asphaltenes: Polycyclic Aromatic Hydrocarbon Dimer Theoretical Studies. Energ Fuel 2013, 27 (4), 1791-1808. 14. Mousavi, M.; Abdollahi, T.; Pahlavan, F.; Fini, E. H., The influence of asphaltene-resin molecular interactions on the colloidal stability of crude oil. Fuel 2016, 183, 262-271. 15. Rezac, J.; Hobza, P., Benchmark Calculations of Interaction Energies in Noncovalent Complexes and Their Applications. Chemical reviews 2016. 16. Grimme, S.; Hansen, A.; Brandenburg, J. G.; Bannwarth, C., Dispersion-Corrected Mean-Field Electronic Structure Methods. Chemical reviews 2016. 17. Grimme, S., Accurate description of van der Waals complexes by density functional theory including empirical corrections. Journal of computational chemistry 2004, 25 (12), 1463-73. 18. Grimme, S., Semiempirical GGA-type density functional constructed with a long-range dispersion correction. Journal of computational chemistry 2006, 27 (15), 1787-99. 19. Grimme, S.; Ehrlich, S.; Goerigk, L., Effect of the damping function in dispersion corrected density functional theory. Journal of computational chemistry 2011, 32 (7), 1456-65. 20. Grimme, S., Density functional theory with London dispersion corrections. Wiley Interdisciplinary Reviews: Computational Molecular Science 2011, 1 (2), 211-228. 21. Contreras-Garcia, J.; Johnson, E. R.; Keinan, S.; Chaudret, R.; Piquemal, J. P.; Beratan, D. N.; Yang, W., NCIPLOT: a program for plotting non-covalent interaction regions. Journal of chemical theory and computation 2011, 7 (3), 625-632. 22. Johnson, E. R.; Keinan, S.; Mori-Sanchez, P.; Contreras-Garcia, J.; Cohen, A. J.; Yang, W., Revealing noncovalent interactions. Journal of the American Chemical Society 2010, 132 (18), 6498-6506. 23. Glendening, E. D.; Landis, C. R.; Weinhold, F., Natural bond orbital methods. Wiley Interdisciplinary Reviews: Computational Molecular Science 2012, 2 (1), 1-42. 24. Neese, F., The ORCA program system. Wiley Interdisciplinary Reviews: Computational Molecular Science 2012, 2 (1), 73-78. 25. Lin, I. C.; Seitsonen, A. P.; Tavernelli, I.; Rothlisberger, U., Structure and Dynamics of Liquid Water from ab Initio Molecular Dynamics-Comparison of BLYP, PBE, and revPBE Density Functionals with and without van der Waals Corrections. Journal of chemical theory and computation 2012, 8 (10), 3902-10. 26. Weigend, F.; Ahlrichs, R., Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: design and assessment of accuracy. Physical Chemistry Chemical Physics 2005, 7 (18), 3297-3305. 27. Kruse, H.; Grimme, S., A geometrical correction for the inter- and intra-molecular basis set superposition error in Hartree-Fock and density functional theory calculations for large systems. J Chem Phys 2012, 136 (15), 154101. 28. Prampolini, G.; Livotto, P. R.; Cacelli, I., Accuracy of Quantum Mechanically Derived Force-Fields Parameterized from Dispersion-Corrected DFT Data: The Benzene Dimer as a Prototype for Aromatic Interactions. Journal of chemical theory and computation 2015, 11 (11), 5182-96. 22

ACS Paragon Plus Environment

Page 22 of 23

Page 23 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

29. Sure, R.; Grimme, S., Comprehensive Benchmark of Association (Free) Energies of Realistic Host-Guest Complexes. Journal of chemical theory and computation 2015, 11 (8), 3785-801. 30. Goerigk, L.; Kruse, H.; Grimme, S., Benchmarking density functional methods against the S66 and S66x8 datasets for non-covalent interactions. Chemphyschem 2011, 12 (17), 3421-33. 31. Wang, H.; Liu, J.; Xu, H.; Ma, Z.; Jia, W.; Ren, S., Demulsification of heavy oil-in-water emulsions by reduced graphene oxide nanosheets. RSC Adv. 2016, 6 (108), 106297-106307. 32. Weigend, F., A fully direct RI-HF algorithm: Implementation, optimised auxiliary basis sets, demonstration of accuracy and efficiency. Physical Chemistry Chemical Physics 2002, 4 (18), 4285-4291. 33. Arago, J.; Orti, E.; Sancho-Garcia, J. C., Nonlocal van der Waals Approach Merged with Double-Hybrid Density Functionals: Toward the Accurate Treatment of Noncovalent Interactions. Journal of chemical theory and computation 2013, 9 (8), 3437-43. 34. Klamt, A.; Schüürmann, G., COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. Journal of the Chemical Society, Perkin Transactions 2 1993, (5), 799-805. 35. Contreras-Garcia, J.; Yang, W.; Johnson, E. R., Analysis of hydrogen-bond interaction potentials from the electron density: integration of noncovalent interaction regions. The journal of physical chemistry. A 2011, 115 (45), 12983-90. 36. Lu, T.; Chen, F., Multiwfn: a multifunctional wavefunction analyzer. Journal of computational chemistry 2012, 33 (5), 580-92. 37. Humphrey, W.; Dalke, A.; Schulten, K., VMD: visual molecular dynamics. Journal of molecular graphics 1996, 14 (1), 33-38. 38. Reed, A. E.; Curtiss, L. A.; Weinhold, F., Intermolecular interactions from a natural bond orbital, donor-acceptor viewpoint. Chemical reviews 1988, 88 (6), 899-926. 39. Groenzin, H.; Mullins, O. C., Molecular size and structure of asphaltenes from various sources. Energ Fuel 2000, 14 (3), 677-684. 40. Mullins, O. C., The Modified Yen Model†. Energ Fuel 2010, 24 (4), 2179-2207. 41. Ruiz-Morales, Y., HOMO-LUMO gap as an index of molecular size and structure for polycyclic aromatic hydrocarbons (PAHs) and asphaltenes: A theoretical study. I. Journal of Physical Chemistry A 2002, 106 (46), 11283-11308. 42. Ruiz-Morales, Y.; Wu, X.; Mullins, O. C., Electronic absorption edge of crude oils and asphaltenes analyzed by molecular orbital calculations with optical spectroscopy. Energ Fuel 2007, 21 (2), 944-952. 43. Ruiz-Morales, Y.; Mullins, O. C., Measured and Simulated Electronic Absorption and Emission Spectra of Asphaltenes. Energ Fuel 2009, 23 (3), 1169-1177. 44. Anzellotti, A. I.; Bayse, C. A.; Farrell, N. P., Effects of nucleobase metalation on frontier molecular orbitals: potential implications for pi-stacking interactions with tryptophan. Inorganic chemistry 2008, 47 (22), 10425-31. 45. Ruizmorales, Y., Molecular Orbital Calculations and Optical Transitions of PAHs and Asphaltenes. 2007; p 95-137. 46. Murray, J. S.; Politzer, P., The electrostatic potential: an overview. Wires Comput Mol Sci 2011, 1 (2), 153-163. 47. Jian, C. Y.; Tang, T.; Bhattacharjee, S., Probing the Effect of Side-Chain Length on the Aggregation of a Model Asphaltene Using Molecular Dynamics Simulations. Energ Fuel 2013, 27 (4), 2057-2067. 48. Zhao, R.; Zhang, R. Q., A new insight into pi-pi stacking involving remarkable orbital interactions. Physical chemistry chemical physics : PCCP 2016, 18 (36), 25452-25457. 49. Zhang, S.; Wang, G.; Lu, Y.; Zhu, W.; Peng, C.; Liu, H., The Interactions between Imidazolium-Based Ionic Liquids and Stable Nitroxide Radical Species: A Theoretical Study. The journal of physical chemistry. A 2016, 120 (30), 6089-102.

23

ACS Paragon Plus Environment