From Single Asphaltenes and Resins to ... - ACS Publications

Oct 4, 2017 - Computing (Rice University), HPC Information Superintend- ence of ... Toulhoat, H.; Nieto-Draghi, C. Energy Fuels 2015, 29, 556−566. (...
0 downloads 0 Views 7MB Size
Article pubs.acs.org/EF

Cite This: Energy Fuels 2017, 31, 11743-11754

From Single Asphaltenes and Resins to Nanoaggregates: A Computational Study Filipe Camargo Dalmatti Alves Lima, Raphael da Silva Alvim, and Caetano Rodrigues Miranda* Departamento de Física dos Materiais e Mecânica, Instituto de Física, Universidade de São Paulo, 05508-090, São Paulo, SP, Brasil S Supporting Information *

ABSTRACT: The understanding of intermolecular interactions on asphaltene nanoaggregates has important implications over the petroleum industry. Here, we investigated the aggregation of asphaltene and resin using a hierarchical approach that combines Density Functional Theory (DFT) with dispersion corrections and molecular mechanics (MM) calculations. From molecular models already established in the literature for specific types of asphaltene (A) and resin (R) molecules, we employed MM simulations to calculate the potential energy surface to obtain the best conformations for the possible dimers A-A, A-R, and R-R, as well as the relevant combinations of trimers A-R-A, A-A-R, and A-A-A, which have been further relaxed by the DFT calculations. Indeed, the formation of the dimer A-A is energetically more favored with respect to A-R and R-R mainly due to the enhanced effect on the intermolecular interaction of the aromatic region of A. In this context, the trimers have shown to be almost 3 times more stable than the dimers. Our result suggests that the nanoaggregates have a charge density well distributed but centralized between the aromatic ring π-orbitals, whereas A molecules are added, creating a tightly packed structure. In this case, the contribution from the aliphatic chains is just steric to stabilize the aggregate and shield the aromatic center for new interactions. Although the π−π stacking can guide the nanoaggregate formation, the presence of the R molecule leads to a possible disaggregation. When R molecules are inserted, the growth of the nanoaggregate seems to be yet continued due to a dipole moment and radius of gyration increasing. However, we observed charge rearrangement from the aromatic center of π−π interaction to aliphatic chains with heteroatoms that displaces into the structure. The HOMO-1 and HOMO degeneracy is broken in this time, being more significant in trimers derived from particular dimer conformation. The formation energy in some asphaltene and resin conformations is decreased compared to the asphaltene nanoaggregates with larger aromatic islands. Therefore, R rich in heteroatoms may be seen like an inner destabilizer naturally present in oil. These findings can guide new methods to control the stability of asphaltene aggregates by external chemical agents through the degradation mechanism directly upon the aliphatic chains.



INTRODUCTION Asphaltenes are large molecules with high density and polar active in certain petroleum fractions. Although insoluble in alkanes such as n-pentane and n-heptane, they are soluble in aromatic solvents such as toluene, benzene, or pyridine. The implications of asphaltene stability within the oil matrix production affect the whole oil supply chain production. This reduces the oil recovery because of mobility changes in the reservoir, by binding to rock pores and deposits in the wells, causing obstruction of the extraction channels, corrosion, and fouling in equipment and catalysts. Moreover, destabilized asphaltenes are the main cause for arterial clogging in pipelines and wellbores, being responsible for formation and strengthening of oil and water emulsions, adherence on refining equipment, coke formation, sedimentation, and plugging during crude oil storage.1−11 Given the reduction of the light-oil reserves in the world, crude oil with higher densities may become economically viable. Therefore, the oil industry is seeking the improvement and knowledge of high-density compounds.3,5 The family of asphaltenes is quite large, which turns out to be experimentally very challenging because the techniques are often unable to distinguish the slight differences between those molecules.3 Decades ago, Dickie and Yen12 proposed a hierarchical approach consisting of two levels of aggregation: first, the © 2017 American Chemical Society

asphaltenes adhere with themselves and form aggregates; then, these aggregates form clusters. In 2010, Mullins proposed the modified Yen model, which defines the asphaltenes as a central moderately-sized polycyclic aromatic hydrocarbon ring system linked with alkane chains on the exterior, with nanoaggregates up to six molecules and clusters of eight nanoaggregates on average.13 The nanoaggregation nature of asphaltenes is a debated topic in the literature.1,3,5,13−27 The most basic information related to the binding mechanism is the π−π stacking of the aromatic groups.3 Rogel attributed the stabilization of aggregates due to interactions between molecules with higher aromaticity and higher aromatic condensation degree.28,29 Currently, it is known that the asphaltenes have molecular density between 500 and 1000 g/mL, which can vary from the source that they are extracted.3 Since asphaltenes are a collection of molecules within a given solubility, the exact molecular structure and composition of their building blocks is unknown and the system is usually characterized by bulk properties.30−35 The literature presents different approaches to Received: July 11, 2017 Revised: October 4, 2017 Published: October 4, 2017 11743

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754

Article

Energy & Fuels design models for asphaltenes,3,5,24,36−39 which most of the time are limited for specific cases. In particular, Boek et al. created asphaltene models using the Quantitative Molecular Representation (QMR),36,37 widely accepted and used in the literature.16,23,36,37 In the first step of the QMR method, Monte Carlo simulations build the threedimensional molecules from the most common heteroatoms, aromatic and aliphatic groups found in the asphaltene fraction. These candidate molecules are selected by the best match with experimental data of mass spectroscopy, nuclear magnetic resonance, and elemental analysis using nonlinear optimization techniques.36,37 The literature shows that asphaltene nanoaggregates have up to 6 nm diameters, which is remarkably small considering their damaging effects for the oil industry.14,15 Experimental studies also stablish the critical nanoaggregation concentration (CNAC) and covers the role of enthalpy and entropy effects for the micelle formation.40−42 Boek et al.,15 in a study of aggregation in capillary flows, show the planes of aromatic groups from asphaltenes preferably align and the free energy interaction reasonably fit with the inverse square of the distance. In addition, the nanoaggregation mechanism was also investigated within different solvents and physical-chemistry conditions. Molecular dynamics (MD) results for the nanoaggregation process in toluene and heptane indicate that free energies are similar for both solvents, showing evidence between the behavior of nanoaggregates and their colloidal structure.14,15 For toluene, the lifetime of the dimers/trimers of asphaltene is of the order of nanoseconds. The lifetime of up to 20 ns was obtained for dimers/trimers in heptane, as shown by Boek et al.14 Within supercritical and liquid carbon dioxide, the asphaltenes show a dense aggregate for a few nanoseconds without a preferred structure, with molecules toward the aggregate during the MD simulation.43 In 2013, Sedghi et al.16 presented a systematic MD study to evaluate the effect of heteroatoms, different sizes and numbers of alkane chains for asphaltene dimers and trimers, showing that the free energy of trimers is 50% higher compared to two monomers. In addition, they suggest that aromaticity may have larger contributions than polarity effects, in agreement with Hansen solubility analysis of asphaltenes.44 Besides, Yang et al.24 presented results for the stabilization of subfractions of asphaltenes water-in-crude oil emulsions. They indicated that asphaltenes self-associate at stacking of polyaromatic rings with heterogeneous structures that behave differently at interfaces, being influenced by heteroatoms such as sulfur and oxygen. Ruiz-Morales et al.45 show that the steric hindrance of asphaltenes in the water−oil interface with π−π stacking interaction of the aromatic cores is observed in the bulk of the oil region. Indeed, the most soluble fraction in the nanoaggregate termination (alkyl chains) has more affinity for the oil region, favoring the interaction from the influence of solvent. In 2006, Alvarez-Ramirez et al.46 established the potential curve between asphaltenes and resin using Density Functional Theory (DFT) calculations, finding interaction energies between 4.00 and 15.00 kcal/mol, respectively. Zamudio-Rivera et al.38 used experimental results and electronic structure methods to suggest mechanism steps for formation of supramolecular interactions with asphaltenes and surfactants. Buehler et al.39 employed Clar sextet theory and electronic structure calculations to design isolated asphaltene molecules with optimized electrons in the aromatic core and reduced

geometrical strain. However, none of these main works have addressed discussions at the electronic level regarding a nanoaggregation mechanism influenced by resin molecules. In this work, the role of the resin molecules in asphaltenes nanoaggregates is investigated from the electronic structure viewpoint. The energetically most favorable aggregation pathways were analyzed due to interaction between asphaltene and resin molecules from different possibilities of conformations. Using a hierarchical approach, we combined molecular mechanics (MM) and DFT to investigate nanoaggregate formation. The MM approach provides a fast and reliable way to explore the potential energy surface between the molecules, a crucial step to survey the most stable structural conformations. In addition, by DFT calculations, we analyzed the charge rearrangement toward the stability of the π−π bonds during the asphaltene + resin interaction having as reference the asphaltene aggregates. Accessing the most fundamental interactions upon these complex systems, we could infer the possible changes in the gyration radius and in the dipole moment at the molecular level. The introduction of the resin is thus seen as part of nanoaggregates of asphaltenes. However, it is expected that resins might be destabilizers on the aggregation process.



COMPUTATIONAL METHODS

Asphaltene and Resin Models. The asphaltene and resin models were proposed by Boek et al., using the QMR method,36 which were also employed in a molecular dynamics investigation for nanoaggregation models.14,15,43 In particular, the asphaltene (A) model with molecular mass 726 Da and chemical composition C53H58S was reported in the literature as one the best matched results from experimental and theoretical studies for UV fluorescence spectra.14,47,48 The model referred as resin (R) is an asphaltene with low molecular weight, 431 Da, and a small aromatic core, which is more like resin than asphaltene with a chemical composition of C26H41NS2.14 These models were chosen because they have several important organic groups, such as thioether, pyridine, polyaromatic cores, and different sizes of aliphatic chains. They form a large set of asphaltene organic groups that allow us to investigate several different interactions possibly presented in the nanoaggregates. Here, the computer simulation took two different levels: Molecular mechanics (MM) and DFT. During the MM procedure, the molecules were considered static. It was done as a first triage to obtain the most favorable energy regions for each set of possible interactions regarding the symmetry of the molecules. Several conformations were tested through a considerably cheap computational cost in each dimer and trimer, which made it possible to choose the best energetically conformations for each dimer and trimer. In a second step, the DFT procedure considered the optimal conformations obtained from MM and accurately relaxed the atomic positions under a tight convergence criterion, until forces were smaller than 0.05 eV/Å. Atomistic Modeling. The classical MM calculations were performed using the LAMMPS49 software, with the CHARMM force field50,51 that was already successfully employed for the description of the light oil model52−54 and description of asphaltenelike molecules.21 The total energies of the dimers and trimers were considered static, without temperature effects and under vacuum conditions. The analysis of the potential energy surface was carried out within the MM method, evaluating the stacking of the aromatic groups of all possible interaction conformations for each dimer combination of asphaltene (A) and resin (R). Since the molecules have a symmetry, when A (PSA) and R (SER) are rotated by 180° with respect to the vector perpendicular to the molecule, we will address the acronym A′ (ASP) and R′ (RES), respectively. Considering the interaction side, the dimer combination has the following distinct conformations, namely, A-A (PSA-PSA) = A′-A′ (ASP-ASP), A-A′ (PSA-ASP), and 11744

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754

Article

Energy & Fuels

Figure 1. (left) Top view of the molecular structures of asphaltene (A) and resin (R). (right) HOMO and LUMO orbitals for A and R. The positive and negative parts of the wave functions are represented by the red and blue spheroids near the atoms. Atomic colors: C (gray), N (blue), H (white), and S (yellow). Isosurface of −0.0005 to 0.0005 electrons/Bohr.3 A′-A (ASP-PSA) for the asphaltene dimers; A-R (PSA-SER), A-R′ (PSA-RES), A′-R (ASP-SER), and A′-R′ (ASP-RES), for the asphaltene and resin dimers; and R-R (SER-SER) = R′-R′ (RESRES), R-R′ (SER-RES), and R′-R (RES-SER) for the resin dimers. Furthermore, we adopted the following procedure for each interaction conformation shown above: two molecules are placed and kept parallel to each other at a starting distance of 5.00 Å, and then, one of the molecules is rotated each 10.00°, increasing the intermolecular distance by 0.50 Å until it reaches 14.00 Å. Finally, for each dimer combination, 540 total energy calculations were performed. In each 2D energy MM map of the potential energy surface, the intermolecular interactions were carried out evaluating the stacking of the aromatic groups of all possible conformations for A, R, combinations of both, and symmetry operations. The procedure to obtain the surface energy map for the trimer started from the dimers that presented the most stable binding energy. The trimer combination has the following conformations, namely: PSA-RES-PSA (A-R′-A), PSA-RES-ASP (A-R′-A′), ASP-PSA-RES (A′-A-R′), PSA-PSA-RES (A-A-R′), PSA-ASP-RES (A-A′-R′), PSAASP-SER (A-A′-R), PSA-ASP-PSA (A-A′-A), PSA-ASP-ASP (A-A′A′). These models were based from the more energy-stable dimers AR′ and A-A′, and they were grouped according to the final setup. The conformation of the third added molecule of asphaltene (A or A′) or resin (R or R′) is underlined. The procedure to obtain the trimer energy landscape is similar to the dimers: a third molecule was initially placed in parallel at 8.00 Å from the center of mass, rotated each 10°, and increasing the intermolecular distance by 0.50 Å until 17.00 Å. In addition, we also investigated the placement of the third molecule below the dimer, ranging from −3.00 to −12.00 Å in the 2D energy MM maps of the potential energy surface. To investigate the dimensions of the nanoaggregates around the center of mass, we have employed the radius of gyration (Rg)

R g2 =

1 M

8.00 Å between the periodic images, thus avoiding spurious interactions. For the trimer, the supercell used has the dimensions (39.54 × 43.98 × 40.60) Å3. The atomic geometry of the systems was fully relaxed with the Broyden−Fletcher−Goldfarb−Shanno (BFGS) quasiNewton algorithm until Hellmann−Feynman forces were smaller than 0.05 eV/Å. Only the Γ-point in the Brillouin zone was considered. The formation energy (Efor) is related to the enthalpy contribution of the nanoaggregation Gibbs free energy. For the aggregates, it was obtained from the difference between the total energy of the system entire (Etot) and isolated molecules (Emol), which can be Emol,1, Emol,2, and Emol,3:

Efor = Etot − Emol ,1 − Emol ,2 − Emol ,3

Since the calculations were carried out at 0 K and vacuum conditions, our results cannot infer the entropy effects that are present in the nanoaggregation aspects. From the literature, our results can be understood as the preassociation of the nanoaggregation formation, which was already discussed in the literature as an enthalpy driven effect.41 On the other hand, the DFT calculations provide unique access to the electronic properties that allow us to investigate the charge rearrangement during the aggregate formation. We evaluate the charge density at the intermolecular interface of the systems given by N

Δρ = ρtotal −

i=1

∑ ρmol , i

(3)

i=1 56

where total and mol are associated with the charge density of the entire system and in isolated components, respectively. While new aspects have been addressed for the aggregation mechanism discussion, we also calculate the electric dipole moment for the asphaltene models, which can be compared to the results from the literature,16 defined by the following equation

N

∑ mi(ri − rmean)2

(2)

N

(1)

P=

∑ qi(ri − rcm) i=1

where M is the total atomic mass of the system, ri is related to the atomic position, and the rmean is the geometric center of the molecule. The main MM molecular models were optimized under DFT calculation. DFT Calculations. For the electronic structure calculations, we used the DFT55 as implemented in the Quantum ESPRESSO suit of codes.56 The Vanderbilt’s ultrasoft pseudopotential was employed to describe the core region of the plane-wave functions.31 We choose the revPBE57 exchange correlation functional modified with vdW dispersions (vdW-DF).58−60 This is a second-order correction of the long-range term, which expresses nonlocal correlations using density− density interactions. The electronic wave functions and charge density were expanded in a plane-wave basis set with a cutoff of 47.00 and 323.00 Ry, respectively. Each dimer was placed in a supercell with dimensions (24.5 × 30.25 × 38.63) Å3, with a minimum distance of

(4)

where rcm is the center of mass of the system and the partial charges (qi) were obtained from the DFT results with the Bader charge model.61−63 The Bader method is based on the partition of the electronic charge density obtained from quantum calculations, searching for zero flux regions, which are 2D surfaces where the charge density is minimum perpendicular to its surface.



RESULTS AND DISCUSSION Asphaltene and Resin Molecules. For comparison purpose, we investigate the electronic structure of the isolated A and R molecules, referring the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). The frontier Kohn−Sham (KS) orbitals for A and R 11745

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754

Article

Energy & Fuels

formation energy of −46.78 kcal/mol (Table 1). This particular case that shows the effect of interactions just of asphaltene molecules was already shown by Ramirez et al.46 It is also possible to find intermolecular interactions between the aliphatic chains and the aromatic groups, from the charge density analysis in the interface of the dimer shown in Figure 3. However, the frontier orbitals are degenerate (Table 2). HOMO-1 is located mainly on the aromatic group of both asphaltene molecules simultaneously as it is HOMO orbitals in Figure 3. Therefore, A-A′ displays a strong influence of π−π stacking interaction between the asphaltene molecules regardless of the interaction of the aliphatic group. This interaction keeps the electronic properties of the asphaltene molecules, since it is worth to note that the KS orbitals from the dimer A-A′ are identical compared to isolated A in Figure 1. The dimer A-R′ displays a larger number of minimum conformations at the MM energy map for the intermolecular distance in the angles 60−100° (Figure 2b, left). R became more planar with respect to the aromatic center of A after the DFT minimization, with an intermolecular distance of about 5.00 Å. Compared to the previous conformation A-A′, the Rg = 7.01 Å indicates a packed structure. Nevertheless, the relative intensity of dipole moment for A-R′ increases in comparison to the dimer A-A′ and to the isolated molecules, presenting a value of 6.50 D, with less stable formation energy (−27.80 kcal/ mol) (Table 1). The electronic changes of A and R are clearly shown after the dimerization A-R′ (Figure 3). Comparing the KS orbitals, the frontier orbitals are located on all aromatic rings for A-R′ and the LUMO states are most affected in R than in A. For this reason, in Table 2, we can see that the HOMO− LUMO energy difference is 2.30 eV after the formation of A-R′, which is very similar to the isolated A molecule (2.35 eV). R displays a HOMO−LUMO energy difference of 3.44 eV with decreased interface electronic charge from the aliphatic chains. Compared to A-A′, the dimer R′-R presents an MM surface energy map with broad minimum energy regions between 90° to 140° and 230° to 350° (Figure 2c, left). The intermolecular distances are close to the interactions between the aliphatic chains and the aromatic group by DFT optimization (Figure 2c, right). Compared to the Boek et al.14 study, there is an agreement with their MD simulations for resin dimers in toluene, which also presented a strong aggregation at the intermolecular distance of 5.00 Å. This case showed the largest amount of short-range contact regions located between the aromatic rings where are the KS orbitals (Figure 3). Rg = 6.95 Å indicates a packed structure. The formation energy of −20.54 kcal/mol is the least stable in comparison to other dimers. However, R′-R presented a very large dipole moment of 10.73 D. Asphaltene-Resin-Asphaltene, Asphaltene-Asphaltene-Resin, and Asphaltene-Asphaltene-Asphaltene Trimers. The asphaltene-resin-asphaltene trimers were built from the dimer A-R′, adding an asphaltene molecule. The surface MM energy map for A-R′-A is shown in Figure 4a, with local minimum regions displaced over almost all rotation angles. After the DFT relaxation, the molecules stayed in a planar configuration with aromatic ring distances of about 4.00 Å, Rg = 7.63 Å, and a formation energy of −64.97 kcal/mol. The HOMO state has π-character and is placed almost at R, while the LUMO is located at A. The charge density in the interfaces is mainly localized near to the heteroatoms region from the resin with dipole moment of 7.63 D.

molecules are shown in Figure 1. While the A molecule displays regions occupied by the orbitals from the center of the polyaromatic region, the R orbitals are spread between the small aromatic center in the nitrogen and the aliphatic group closest to the sulfur. Dipole moments shown in Table 1, 2.46 Table 1. Formation Energy (kcal/mol), Module of the Dipole Moment (D), and Radius of Gyration (Å) for the Isolated Models, Dimers and Trimers of Asphaltene and Resin

asphaltene (A) resin (R) A-A′ A-R′ R′-R A-R′-A A-R′-A′ A′-A-R′ A-A-R′ A-A′-R′ A-A′-R A-A′-A A-A′-A′

formation energy (kcal/mol)

dipole moment (D)

radius of gyration (Å)

−46.78 −27.80 −20.54 −64.97 −68.60 −64.00 −66.63 −76.66 −80.53 −89.25 −94.10

2.46 4.04 1.62 6.50 10.73 6.65 5.92 5.75 5.27 2.52 4.60 2.44 3.98

6.50 6.66 6.51 7.01 6.95 7.63 7.21 7.75 7.46 7.10 6.97 7.52 7.07

and 4.04 D, respectively, for A and R, are used for reference for the dimer and trimer cases. In particular, the results obtained for A are in very good agreement with results obtained from classical molecular dynamics.16 The molecules A and R presented Rg of 6.50 and 6.66 Å, respectively. Asphaltene-Asphaltene, Asphaltene-Resin, and ResinResin Dimers. Due to several interaction pathways between A and R molecules, we explore the surface energy maps of the structural conformations for each dimer and trimer studied here. This procedure allows the first energetically screening of the best candidates for dimers and trimers structures, which were posteriorly minimized by DFT calculations. Particularly, it was possible to determine the minima energy regions, intermolecular distances, and the position of the organic groups involved in the interactions. We showed in Figure 2 the MM energy maps and, from this, the more energy-stable dimerization related to each conformation A-A′, A-R′, and R′-R. The complete analysis for the other interaction energetically less stable conformations (A-A, A′-A, A-R, A′-R, A′-R′, R-R, and R′-R) is available in the Supporting Information. These results were compared with the literature in order to ensure the description of the nanoaggregates.14,43 In general, the range of the energy surface obtained for our results from classical MM is 10.00 kcal/mol on average. The results for dimers are in agreement with the results obtained by Rogel,29 which showed stabilization energies of different asphaltene and resin molecules in organic solvent, also few contributions to the stabilization energy related to structural changes. Our energy profiles are in qualitative agreement with results from the literature;14,29 although here, it is not considered solvent effects. For the dimer A-A′, the MM energy map shows distances between the aromatic islands being regions of minimum energy (Figure 2a, left). This distance is reduced to 4.00 Å on average after DFT optimization (Figure 2a, right), with Rg = 6.51 Å similar to A, but a smaller dipole moment (1.62 D) and 11746

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754

Article

Energy & Fuels

Figure 2. (left) 2D energy MM maps for the total energy obtained for the all the dimers (a) A-A′, (b) A-R′, (c) R′-R rotated by 180° with respect to the aromatic center. The warm to cold colors represent the high energy to lower energy interactions. The yellow color shows the region that the total energy is above the threshold specified. (right) Dimer structures optimized after the DFT calculations with aromatic ring intermolecular minimum distance. Atomic colors: C (gray), N (blue), H (white), and S (yellow).

electronic states tend to be closer to the Fermi energy, revealing the clear influence of the resin’s heteroatoms into the overall electronic structure. For both conformations, the combined orbitals are not just a superposition of the isolated ones. Particularly, the energy difference observed at HOMO-1 and HOMO for the dimer A-R′ is increased for the trimer A-R′-A (Table 2), suggesting that the frontier electronic states can be influenced by the resin molecule breaking the degeneracy. Considering the isolated molecules, the dipole moment of the trimer A-R′-A′ is smaller (5.92 D), as seen in Table 1. The asphaltene-asphaltene-resin trimers represent asphaltene and resin interaction with the dimers A-R′ and A-A′, respectively. From A-R′, Figure 4c shows the MM potential energy surface for A′-A-R′, with favored regions between 100° and 150°. The DFT relaxation led to an interplanar aromatic region distance of 3.80 and 4.00 Å for asphaltene-asphaltene and asphaltene-resin interfaces, respectively. A′-A-R′ is the least

The MM surface energy map for A-R′-A′ is shown in Figure 4b. It is seen that the local minimum of the energy is only at the intermolecular region between 200° and 230°. The major structural differences between the classical results and DFT are the following: the aromatic distances between two asphaltene and asphaltene-resin are 4.00 and 3.90 Å, respectively. The aliphatic chains kept approximately the same structure, with variations up 0.25 Å in the absence of any specific contraction or expansion. The Rg of 7.21 Å shows that the structure is very packed in comparison to the other models, increasing less than 1.00 Å in comparison to the isolated A (6.50 Å). We obtained a formation energy of −68.60 kcal/mol for this trimer A-R′-A′, almost 3 times than the energy obtained for the dimer A-R′. As seen in Figure 5, the HOMO state is located on R, indicating the pathway responsible for the formation of the nanoaggregate. From the isolated molecules to formations of the trimers A-R′-A and A-R′-A′, we observed that the unoccupied 11747

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754

Article

Energy & Fuels

Figure 3. KS HOMO and LUMO and charge density (Δρ) shown at the intermolecular interface for A-A′, A-R′, and R′-R. The positive and negative parts of the wave functions are represented by the red and blue spheroids near the atoms. HOMO and LUMO energies are represented by EH and EL, respectively. For Δρ, the regions of charge accumulation and depletion are represented by the red and blue spheroids near the atoms. The vector of dipole moment is displayed in orange on the Δρ figure. Atomic colors: C (gray), N (blue), H (white), and S (yellow). Isosurface of −0.0005 to 0.0005 electrons/Bohr.3

Formed from the dimer A-A′, the energy surface of A-A′-R′ is shown in Figure 4e, with the most favorable interaction energies placed at 200° and 11 Å from the center of the trimer. The intermolecular aromatic rings distances are 3.60 and 3.70 Å after the DFT relaxation for asphaltene-asphaltene and asphaltene-resin interfaces, respectively. The formation energy obtained is −76.66 kcal/mol, with a dipole moment of 2.52 D and a gyration radius of 7.10 Å. Compared to the previous cases, LUMO orbitals were found only at the bottom molecule A and the charge difference spread among the aromatic rings and the asphaltene aliphatic chains. The most favorable conformation found is A-A′-R, with a formation energy of −80.53 kcal/mol and the smallest radius of gyration (6.97 Å). The dipole moment, 4.60 D, is small compared to the isolated molecules A and R. The frontier orbitals also presented a strong π characteristic. They were found at R for HOMO and in the central A for LUMO. The last two trimers are composed only by asphalteneasphaltene-asphaltene from A-A′, namely, A-A′-A and A-A′-A′. The potential energy surface for both cases is shown in Figure 6, with interplanar aromatic distances of about 3.60 Å after the DFT calculations. In agreement with the literature,46 the trimers of asphaltenes showed the more stable formation energy compared to other resin trimmers studied, with −89.25 and −94.10 kcal/mol for A-A′-A and A-A′-A′, respectively. Due to the absence of the molecule R, it can be seen in Figure 7 that both systems remain the π-orbital contributions at both HOMO and LUMO. Remarkably, LUMO for A-A′-A is spread over all asphaltene molecules. While A-A′-A presents the second highest radius of gyration (7.52 Å), and the smallest

Table 2. Energy Levels in (eV) for the Orbitals HOMO-2, HOMO-1, HOMO, and LUMO of the Asphaltene and Resin Molecules, and the Dimer and Trimer Conformationsa energy level asphaltene (A) resin (R) A-A′ A-R′ R′-R A-R′-A A-R′-A′ A′-A-R′ A-A-R′ A-A′-R′ A-A′-R A-A′-A A-A′-A′ a

HOMO-2

HOMO-1

LUMO

−0.35 −0.38 −0.27 −0.31 −0.37 −0.36 −0.26 −0.14 −0.15 −0.12 −0.16 −0.10 −0.11

−0.28 −0.38 −0.03 −0.09 −0.05 −0.16 −0.08 −0.09 −0.06 −0.04 −0.10 −0.06 −0.03

2.42 3.48 2.35 2.36 3.17 2.06 2.15 2.26 2.26 2.31 2.24 2.32 2.31

HOMO was set as zero reference in all cases.

favorable trimer in terms of interaction energy (−64.00 kcal/ mol), but it displays the highest Rg = 7.75 Å. Its HOMO orbitals are located only on the R molecule, while the LUMO one stays between the A aromatic rings. The charge density spreads at the center of the trimer, as observed in Figure 5. The trimer A-A-R′ is similar to A′-A-R′ in terms of aromatic region distance, electronic orbitals, and charge density, but with a higher formation energy of −66.63 kcal/mol and smaller Rg (7.46 Å). 11748

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754

Article

Energy & Fuels

Figure 4. (left) 2D energy MM maps for the total energy obtained for the trimers (a) A-R′-A, (b) A-R′-A′, (c) A′-A-R′, (d) A-A-R′, (e) A-A′-R′, and (f) A-A′-R rotated by 180° with respect to the aromatic center. The warm to cold colors represent the high energy to lower energy interactions. The yellow color shows the region that the total energy is above the threshold specified. (right) Trimer structures optimized after the DFT calculations with aromatic ring intermolecular minimum distance. Atomic colors: C (gray), N (blue), H (white), and S (yellow).

We addressed the outcome of our findings with a different perspective on aggregation mechanism from the inclusion of R consisting of heteroatoms in the aliphatic chains in comparison to the literature of the nanoaggregation mechanism. The main

dipole moment (2.44 D), A-A′-A′ shows average values for these properties: Rg = 7.07 Å and 3.98 D. Nanoaggregation Formation. The results presented can be understood as the core information for the nanoaggregates. 11749

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754

Article

Energy & Fuels

Figure 5. KS HOMO and LUMO and charge density (Δρ) shown at the intermolecular interface for A-R′-A, A-R′-A′, A′-A-R′, A-A-R′, A-A′-R′, and A-A′-R. The positive and negative parts of the wave functions are represented by the red and blue spheroids near the atoms. HOMO and LUMO energies are represented by EH and EL, respectively. For Δρ, the regions of charge accumulation and depletion are represented by the red and blue spheroids near the atoms. Atomic colors: C (gray), N (blue), H (white), and S (yellow). Isosurface of −0.0005 to 0.0005 electrons/Bohr.3

theoretical results from the literature5,19,43,45,48 did not show the electronic contribution mainly with the presence of heteroatoms. Here, our work is more accurate by showing new possibilities in the electronic level of beginning of aggregation from both dimers asphaltene-asphaltene and asphaltene-resin. Goual et al.41 conclude that the inclusion of resins does not modify the asphaltene aggregation. Eyssautier et al.64,65 show there is not solvent penetration into the asphaltene nano-

aggregates due to the aliphatic chains. In the work of Lisitza et al.,40 the aggregation shows negative enthalpy, due to π−π stacking (including asphaltene-toluene interaction), and positive entropy, due to the entropy gain in nonpolar solutions. The rearrangement of the aliphatic chains in the aggregates is hindered with less degree of freedom, and then this could decrease the entropy. However, from the size increasing of aggregates there are two possibilities: (1) different orientations of toluene molecules around the aggregate or (2) more toluene 11750

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754

Article

Energy & Fuels

Figure 6. (left) 2D energy MM maps for the total energy obtained for the trimers (a) A-A′-A and (b) A-A′-A′ rotated by 180° with respect to the aromatic center. The warm to cold colors represent the high energy to lower energy interactions. The yellow color shows the region that the total energy is above the threshold specified. (right) Trimer structures optimized after the DFT calculations with aromatic ring intermolecular minimum distance. Atomic colors: C (gray), H (white), and S (yellow).

Figure 7. KS HOMO and LUMO and charge density (Δρ) shown at the intermolecular interface for A-A′-A and A-A′-A′. The positive and negative parts of the wave functions are represented by the red and blue spheroids near the atoms. HOMO and LUMO energies are represented by EH and EL, respectively. For Δρ, the regions of charge accumulation and depletion are represented by the red and blue spheroids near the atoms. Atomic colors: C (gray), H (white), and S (yellow). Isosurface of −0.0005 to 0.0005 electrons/Bohr.3

11751

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754

Article

Energy & Fuels

even where resin can be aggregate before the π−π stacking of asphaltenes (trimers A-R′-A and A-R′-A′). However, this nanoaggregate pathway seems to display an important role on the limitation in the growth of the asphaltene nanoaggregates even under the π−π stacking seen in the trimers A′-A-R′, A-A-R′, A-A′-R′, and A-A′-R. It is important observe that the tendency of HOMO-1 and HOMO to degenerate in dimers is more significantly influenced in those trimers with resin (Table 2). However, this also is related to the type of nanoaggregate conformation and not only the presence of the resin. Although the degeneracy break of the frontier orbitals above −0.09 eV observed in the asphaltene-resin dimer does not occur in all analyzed trimers (Table 2), we know that the trimers composed with R display a changing of the HOMO localization (Figure 5). The main finding of this work was to show the aggregation mechanism is related to the rearrangement of electronic charge toward the stability of the π−π stacking in the presence of resin. Besides the steric effect, the growth of the nanoaggregates was followed up by the changes toward to the degeneracy of the electronic states in Table 2 and, consequently, of the interaction energies in Table 1. The decrease or no increase linear trends of the dipole moment and radius of gyration of the aggregates with resin are compared to the strong interactions between asphaltene molecules in the asphaltene nanoaggregation, in which the increase of the radius of gyration is well-known. Even so, the asphaltene nanoaggregates are limited up to a few molecules. The tendency of reduction of the dipole moment with resin decreases the dragging force for agglomeration of new asphaltene molecules. This directly affects the most external aromatic groups that might not have electronic states available to create new strong interactions. Therefore, we can indeed infer that the upcoming of further asphaltene and resin molecules may limit the size of the nanoaggregate, but also provide a pathway to obtain nanoaggregates energetically less stable from the inclusion of destabilizing resins.

molecules are in solution instead surround the aggregate, increasing the entropy. Besides, we did not simulate toluene solvent. We anticipate that the entropic influence would be somewhat lower, but this does not affect the nanoaggregates. In the work of Goual et al.,41 the preassociation for dimers and trimers is enthalpically driven, which is exactly the situation in our study. From this idea, the solvent entropy change contributes for a positive increase in entropy upon aggregate formation. The entropic increasing would occur just for the micelle formation with CNAC. This was already shown in the literature as an important point to investigate the micelle formation.45 Our aim here is just to evaluate the enthalpy effect in the step before the formation of asphaltene clusters. By comparing of the dimers A-A′, A-R′, and R′-R (Figure 3), the A-A′ and A-R′ dimers have the most favorable formation energy, with interactions mainly between the aromatic groups. However, we see that the formation energy is inversely proportional to the dipole moment and radius of gyration (Table 1). R′-R displays a different behavior for the radius of gyration maybe due to less influence of π−π stacking in this aggregate. Our calculated radii of gyration due to the presence of resins in aggregates are comparable with ones experimentally obtained by small-angle X-ray (SAXS) and neutron (SANS) scattering. According Eyssautier et al.,64,65 the hydrodynamic radii, which could be necessary due to entropy increasing from the solvent as it is discussed in the experimental studies, may be equal to radii of gyration in larger nanoaggregates, although this represents a few percent of the asphaltene volume. As seen by calculated total energies, our work indicates a main path through interaction of asphaltenes, which tends to be sterically hindered as the aggregate increases. In particular, the work of Ruiz-Morales45 et al. by dissipative particle dynamics (DPD) theoretical simulation is based on steric hindrance of asphaltenes in the water−oil interface with π−π stacking interaction of the aromatic cores being observed in the bulk of the oil region. Indeed, the most soluble fraction in the nanoaggregate termination (alkyl chains) has more affinity for the oil region, favoring thus the interaction from the influence of toluene. Furthermore, Ruiz-Morales45 et al. do not show the electronic contribution and in the absence of heteroatoms. The fact is that, when the resin molecule is present in the beginning of the aggregation, the dimer is increasingly influenced by electronic charge, showing that the highest binding energy found on asphaltene-asphaltene dimers can be understood as the almost “pure” π−π stacking. This can be seen that this contribution comes from the charge rearrangement between the molecules that spread on along the interaction regions. Using this information as a pathway for the nanoaggregation, inclusion of a molecule A stabilizes even more the interaction energy of the formed trimer, creating thus a possible conformation also for the molecule R when it takes place in the aggregate. The formation energy of the asphaltene trimers is almost 3 times higher in comparison to the dimer of resins. Consequently, the charge rearrangement seen in the dimers also leads to the increasing of the electric dipole moment upon the trimmers, mainly ones of asphaltene-resin-asphaltene compared to asphaltene-asphaltene-asphaltene (Table 1). The formation energy follows inversely proportional to the radius of gyration in all the cases of analyzed trimmers. Resins and heteroatoms present in the interaction region of aggregates do not decrease significantly the formation energy of the aggregate,



CONCLUSIONS We presented a systematic study of the asphaltene and resin nanoaggregation mechanism using a hierarchical approach based on classical MM and DFT calculations with vdW corrections. From the atomistic results, we analyzed the potential energy surface for each asphaltene-asphaltene, asphaltene-resin, and resin-resin dimer and asphaltene-resinasphaltene, asphaltene-asphaltene-resin, and asphaltene-asphaltene-asphaltene trimer. Interestingly, the structures with more accessible regions of lower energy interactions are also the structures that displayed the lowest minimum energy. The energy difference gradient of the structures is on the order of 30.00 kcal/mol, which means they are possible upon thermal fluctuations. As solvent interactions were already discussed in the literature,14,43 here our focus is the fundamental intermolecular interactions in the nanoaggregates. Our work shows thus an important look at the electronic energy from specific interaction on the formation of the analyzed dimers and trimers independently of the possible difference of entropic energy when they are in the oil bulk. Accordingly, a temperature at 0 K was enough to determine the electronic energy that is the main goal of this work, showing further insights on the electronic structure changes upon asphaltene interaction to guide the beginning formation of nanoaggregates from the presence of the resin molecule in different association mechanisms. It is also worth noticing that 11752

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754

Energy & Fuels



the behavior of the trimer formation and dipole moment reported here are in agreement with literature results.16 On one hand, our results cannot address further discussions into the CNAC,40,41 since this is related to larger aggregates between 6 and 8 molecules as suggested in the literature. On the other hand, the growth and termination of the nanoaggregates were followed up by the changes toward to the degeneracy of the electronic states, dipole moment, and radius of gyration. Accordingly, our results at the electronic level cover the fundamental aspects of how the nanoaggregation mechanism is affected by resin, taking into account mainly the dimers and trimers of asphaltenes and resins. Whereas the aromatic region dominates almost all interactions, the DFT results showed that π-orbitals are the majority even in the presence of resin and particularly degenerated in the frontier of KS orbitals on the dimers. The trimer of asphaltenes indeed presented the most stable formation energy even in the absence of degeneracy (−94.10 kcal/mol). This is almost 2 times the interaction energy obtained for the dimer of asphaltenes from the isolated asphaltene molecules, carrying a possible association for such composition. However, this energetically stable growth of aggregates by asphaltene molecules limits the dipole moment and the radius of gyration due to centralization of charge at the π−π interaction. Although the nanoaggregation mechanism is related to the aromatic groups in the asphaltene molecules, the results of electronic charge density show charge rearrangement upon the interaction regions in the trimers with resin, equally spread toward the aliphatic chains by the heteroatoms. For this reason, the asphaltene-resin-asphaltene trimers with higher dipole moment and radius of gyration display smaller interaction energy, as observed from conformation A-R′-A. In summary, the nanoaggregates show a strong competition of the π-orbitals, leading to a degeneracy of the electronic states. The role of the heteroatoms is important to influence the charge rearrangement. It can enhance the interaction among the aliphatic chains but reducing the interaction on aromatic groups. However, this effect does not reduce the effective dipole moment. The driving force to attract further asphaltene molecules into the nanoaggregate can be decreased by aggregation outside of the π−π stacking with the breakdown of degeneracy, limiting the size into a few nanometers. The destabilization mechanism of the nanoaggregates seems to be thus related to agents that pull the charges to the border again. This can be achieved for instance with surfactants that have selectivity for the aliphatic chains of the resin molecules. These results address very important findings regarding the electronic competition on the dimer/trimer formation, providing new insights to the nanoaggregates fundamental aspects. The findings reported here might be extrapolated for different types of asphaltene and resin molecules able to aggregate: the key factor is the fundamental interaction of the heteroatoms in the aliphatic chains, system size, and the organic functional groups rather than the design of the molecules. From an industrial point of view, the knowledge of the aggregation mechanism may allow the design of new surfactants or compounds with specific functional groups to control the formation and destabilization of the asphaltene nanoaggregates.

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.energyfuels.7b02002. Molecular mechanics results obtained for the interacting dimers. Table S1: minimum total energy and torsion angle obtained for the dimers from molecular mechanics (MM) calculations. Figure S1: 2D energy MM maps for the dimers A-A, A′-A, A-R, A′-R, A′-R′, R-R, R′-R (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Filipe Camargo Dalmatti Alves Lima: 0000-0001-7062-5450 Caetano Rodrigues Miranda: 0000-0002-8008-4907 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This project is financially supported by the PETROBRAS and Brazilian agencies CNPq and FAPESP. The authors acknowledge the computational time provided by the Blue Gene/Q supercomputer, supported by the Center of Research Computing (Rice University), HPC Information Superintendence of Universidade de São Paulo, UFABC, CENAPAD-SP and IFUSP SAMPA group.



REFERENCES

(1) Mullins, O. C.; Sabbah, H.; Eyssautier, J.; Pomerantz, A. E.; Barré, L.; Andrews, a. B.; Ruiz-Morales, Y.; Mostowfi, F.; McFarlane, R.; Goual, L.; Lepkowicz, R.; Cooper, T.; Orbulescu, J.; Leblanc, R. M.; Edwards, J.; Zare, R. N. Energy Fuels 2012, 26 (7), 3986−4003. (2) Mullins, O. C.; Pomerantz, a. E.; Dong, C.; Petro, D.; Seifert, D. J.; Zuo, J. Y.; et al. Oilfield Rev. 2013, 24, 14−25. (3) Adams, J. J. Energy Fuels 2014, 28 (5), 2831−2856. (4) Zhang, D.; Creek, J.; Jamaluddin, a J.; Marshall, A. G.; Rodgers, R. P.; Mullins, O. C.; et al. Oilfield Rev. 2007, 19, 22−43. (5) Jover, J. F.; Müller, E. a.; Haslam, A. J.; Galindo, A.; Jackson, G.; Toulhoat, H.; Nieto-Draghi, C. Energy Fuels 2015, 29, 556−566. (6) Higaki, Y.; Kobayashi, M.; Murakami, D.; Takahara, A. Polym. J. 2016, 48 (4), 325−331. (7) Meirer, F.; Kalirai, S.; Morris, D.; Soparawalla, S.; Liu, Y.; Mesu, G.; Andrews, J. C.; Weckhuysen, B. M. Sci. Adv. 2015, 1, e1400199. (8) Stachowiak, C.; Viguié, J. R.; Grolier, J. P. E.; Rogalski, M. Langmuir 2005, 21 (11), 4824−4829. (9) He, L.; Lin, F.; Li, X.; Sui, H.; Xu, Z. Chem. Soc. Rev. 2015, 44 (15), 5446−5494. (10) Higaki, Y.; Hatae, K.; Ishikawa, T.; Takanohashi, T.; Hayashi, J.; Takahara, A. ACS Appl. Mater. Interfaces 2014, 6 (22), 20385−20389. (11) Elsner, M.; Hoelzer, K. Environ. Sci. Technol. 2016, 50 (7), 3290−3314. (12) Dickie, J. P.; Yen, T. F. Anal. Chem. 1967, 39 (14), 1847−1852. (13) Mullins, O. C. Energy Fuels 2010, 24 (4), 2179−2207. (14) Headen, T. F.; Boek, E. S.; Skipper, N. T. Energy Fuels 2009, 23, 1220−1229. (15) Boek, E. S.; Wilson, A. D.; Padding, J. T.; Headen, T. F.; Crawshaw, J. P. Energy Fuels 2010, 24 (4), 2361−2368. (16) Sedghi, M.; Goual, L.; Welch, W.; Kubelka, J. J. Phys. Chem. B 2013, 117 (18), 5765−5776. (17) Hoepfner, M. P.; Fogler, H. S. Langmuir 2013, 29 (49), 15423− 15432. 11753

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754

Article

Energy & Fuels (18) Zabala, R.; Mora, E.; Cespedes, C.; Guarin, L.; Acuna, H.; Botero, O.; Patino, J. E.; Cortes, F. B. In: Paper OTC-24310, presented at the Offshore Technology Conference, Rio de Janeiro, Brazil, 29−31 October, 2013; DOI: 10.4043/24310-MS. (19) Mikami, Y.; Liang, Y.; Matsuoka, T.; Boek, E. S. Energy Fuels 2013, 27 (4), 1838−1845. (20) Wang, S.; Xu, J.; Wen, H. Comput. Phys. Commun. 2014, 185 (12), 3069−3078. (21) Rane, K. S.; Kumar, V.; Wierzchowski, S.; Shaik, M.; Errington, J. R. Ind. Eng. Chem. Res. 2014, 53 (45), 17833−17842. (22) Jian, C.; Tang, T.; Bhattacharjee, S. Energy Fuels 2014, 28 (6), 3604−3613. (23) Liu, J.; Zhao, Y.; Ren, S. Energy Fuels 2015, 29 (2), 1233−1242. (24) Yang, F.; Tchoukov, P.; Dettman, H.; Teklebrhan, R. B.; Liu, L.; Dabros, T.; Czarnecki, J.; Masliyah, J.; Xu, Z. Energy Fuels 2015, 29 (8), 4783−4794. (25) Forte, E.; Taylor, S. E. Adv. Colloid Interface Sci. 2015, 217, 1− 12. (26) Mehranfar, M.; Gaikwad, R.; Das, S.; Mitra, S. K.; Thundat, T. Langmuir 2014, 30 (3), 800−804. (27) Natarajan, A.; Xie, J.; Wang, S.; Masliyah, J.; Zeng, H.; Xu, Z. J. Phys. Chem. C 2011, 115 (32), 16043−16051. (28) Rogel, E. Colloids Surf., A 1995, 104 (1), 85−93. (29) Rogel, E. Energy Fuels 2000, 14 (3), 566−574. (30) Schuler, B.; Meyer, G.; Peña, D.; Mullins, O. C.; Gross, L. J. Am. Chem. Soc. 2015, 137 (31), 9870−9876. (31) Pomerantz, A. E.; Hammond, M. R.; Morrow, A. L.; Mullins, O. C.; Zare, R. N. J. Am. Chem. Soc. 2008, 130 (23), 7216−7217. (32) Sanabria-Ortega, G.; Pécheyran, C.; Bérail, S.; Donard, O. F. X. Anal. Chem. 2012, 84 (18), 7874−7880. (33) Gaspar, A.; Zellermann, E.; Lababidi, S.; Reece, J.; Schrader, W. Anal. Chem. 2012, 84 (12), 5257−5267. (34) Ramachandran, V.; van Tol, J.; McKenna, A. M.; Rodgers, R. P.; Marshall, A. G.; Dalal, N. S. Anal. Chem. 2015, 87 (4), 2306−2313. (35) Mochida, I.; Okuma, O.; Yoon, S.-H. Chem. Rev. 2014, 114 (3), 1637−1672. (36) Boek, E. S.; Yakovlev, D. S.; Headen, T. F. Energy Fuels 2009, 23 (3), 1209−1219. (37) Al Halwachi, H. K.; Yakovlev, D. S.; Boek, E. S. Energy Fuels 2012, 26 (10), 6177−6185. (38) Pons-Jiménez, M.; Cartas-Rosado, R.; Martínez-Magadán, J. M.; Oviedo-Roa, R.; Cisneros-Dévora, R.; Beltrán, H. I.; Zamudio-Rivera, L. S. Colloids Surf., A 2014, 455, 76−91. (39) Martín-Martínez, F. J.; Fini, E. H.; Buehler, M. J. RSC Adv. 2015, 5 (1), 753−759. (40) Lisitza, N. V.; Freed, D. E.; Sen, P. N.; Song, Y.-Q. Energy Fuels 2009, 23 (3), 1189−1193. (41) Goual, L.; Sedghi, M.; Zeng, H.; Mostowfi, F.; McFarlane, R.; Mullins, O. C. Fuel 2011, 90 (7), 2480−2490. (42) Wu, Q.; Pomerantz, A. E.; Mullins, O. C.; Zare, R. N. Energy Fuels 2014, 28 (1), 475−482. (43) Headen, T. F.; Boek, E. S. Energy Fuels 2011, 25 (2), 503−508. (44) Hansen, C. M. Hansen Solubility Parameters: A User’s Handbook; CRC Press: Boca Raton, FL, 1999. (45) Ruiz-Morales, Y.; Mullins, O. C. Energy Fuels 2015, 29 (3), 1597−1609. (46) Alvarez-Ramirez, F.; Ramirez-jaramillo, E.; Ruiz-morales, Y. Energy Fuels 2006, 20 (1), 195−204. (47) Ruiz-Morales, Y.; Mullins, O. C. Energy Fuels 2007, 21 (1), 256−265. (48) Ruiz-Morales, Y.; Wu, X.; Mullins, O. C. Energy Fuels 2007, 21 (2), 944−952. (49) Plimpton, S. J. Comput. Phys. 1995, 117 (1), 1−19. (50) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; Mackerell, A. D. J. Comput. Chem. 2009, NA−NA. (51) Yu, W.; He, X.; Vanommeslaeghe, K.; MacKerell, A. D. J. Comput. Chem. 2012, 33 (31), 2451−2468.

(52) Kunieda, M.; Nakaoka, K.; Liang, Y.; Miranda, C. R.; Ueda, A.; Takahashi, S.; Okabe, H.; Matsuoka, T. J. Am. Chem. Soc. 2010, 132 (51), 18281−18286. (53) de Lara, L. S.; Michelon, M. F.; Miranda, C. R. J. Phys. Chem. B 2012, 116 (50), 14667−14676. (54) de Lara, L. S.; Rigo, V. A.; Miranda, C. R. J. Phys. Chem. C 2016, 120 (12), 6787−6795. (55) Kohn, W.; Sham, L. J. J. Phys. Rev. 1965, 140 (4A), A1133− A1138. (56) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; De Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M. J. Phys.: Condens. Matter 2009, 21 (39), 395502. (57) Zhang, Y.; Yang, W. Phys. Rev. Lett. 1998, 80 (4), 890−890. (58) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Phys. Rev. Lett. 2004, 92 (24), 246401. (59) Thonhauser, T.; Cooper, V. R.; Li, S.; Puzder, A.; Hyldgaard, P.; Langreth, D. C. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76 (12), 125112. (60) Román-Pérez, G.; Soler, J. M. Phys. Rev. Lett. 2009, 103 (9), 096102. (61) Sanville, E.; Kenny, S. D.; Smith, R.; Henkelman, G. J. Comput. Chem. 2007, 28 (5), 899−908. (62) Tang, W.; Sanville, E.; Henkelman, G. J. Phys.: Condens. Matter 2009, 21 (8), 084204. (63) Yu, M.; Trinkle, D. R. J. Chem. Phys. 2011, 134 (6), 064111. (64) Eyssautier, J.; Frot, D.; Barré, L. Langmuir 2012, 28 (33), 11997−12004. (65) Eyssautier, J.; Levitz, P.; Espinat, D.; Jestin, J.; Gummel, J.; Grillo, I.; Barré, L. J. Phys. Chem. B 2011, 115 (21), 6827−6837.

11754

DOI: 10.1021/acs.energyfuels.7b02002 Energy Fuels 2017, 31, 11743−11754