Role of Naphthenic Acids in Controlling Self-Aggregation of a

Article pubs.acs.org/JPCB

Role of Naphthenic Acids in Controlling Self-Aggregation of a Polyaromatic Compound in Toluene Robel B. Teklebrhan,† Cuiying Jian,‡ Phillip Choi,† Zhenghe Xu,*,†,∥ and Johan Sjöblom§ †

Department of Chemical and Materials Engineering and ‡Department of Mechanical Engineering, University of Alberta, Edmonton, Alberta, Canada T6G 1H9 § Department of Chemical Engineering, Ugelstad Laboratory, Norwegian University of Science and Technology (NTNU), Trondheim NO 7491, Norway ∥ Institute of Nuclear and New Energy Resources, Tsinghua University, Beijing 100084, China S Supporting Information *

ABSTRACT: In this work, a series of molecular dynamics simulations were performed to investigate the effect of naphthenic acids (NAs) in early stage self-assembly of polyaromatic (PA) molecules in toluene. By exploiting NA molecules of the same polar functional group but different aliphatic/cycloaliphatic nonpolar tails, it was found that irrespective of the presence of the NA molecules in the system, the dominant mode of π−π stacking is a twisted, offset parallel stacking of a slightly larger overlapping area. Unlike large NA molecules, the presence of small NA molecules enhanced the number of π−π stacked PA molecules by suppressing the hydrogen bonding interactions among the PA molecules. Smaller NA molecules were found to have a higher tendency to associate with PA molecules than larger NA molecules. Moreover, the size and distribution of π−π stacking structures were affected to different degrees by changing the size and structural features of the NA molecules in the system. It was further revealed that the association between NA and PA molecules, mainly through hydrogen bonding, creates a favorable local environment for the overlap of PA cores (i.e., π−π stacking growth) by depressing the hydrogen bonding between PA molecules, which results in the removal of some toluene molecules from the vicinity of the PA molecules. crystals,15 etc., exploring the nature of self-assembly of perylene bisimide based PA compounds for the construction of functional materials in an organic medium is of paramount importance. Perylene bisimide based PA compounds are unique compounds for their optical and electronic properties, chemical robustness, thermal stability, extended quadrupolar π system, and high fluorescence quantum yields.16−19 From a molecular perspective, it is widely believed that the noncovalent attraction between the large and fused aromatic cores of these PA molecules is the dominant driving force for their nanoaggregation.8,20−26 On the other hand, recent literature27,28 suggests that an association dominated only by aromatic−aromatic stacking interactions does not fully capture the aggregation mechanism of PA molecules in a real system, e.g., crude oil processing. This is because there is a wide range of functional groups present in the peripheral side chains of PA molecules such as carboxylic, thiophenic, sulfide, hydroxy, pyridinic, and pyrrolic groups.27 Different interactions (e.g., Brønsted acid−base interactions, dipole−dipole interactions,

1. INTRODUCTION Understanding the mechanism by which polyaromatic (PA) hydrocarbon molecules associate into nanoaggregates in solution is vital in preventing processing problems in areas such as crude oil processing, water treatment plants, oil transfer pipelines, etc.1−8 The unwanted nanoaggregates formed initially grow gradually into larger aggregates of various sizes and eventually lead to coagulation and precipitation in various sections of the oil production chains, thus blocking transport pipes and leading to system failures. Moreover, during the production and transportation processes, the nanoaggregation of these molecules greatly influences the solubility, flow behavior, viscosity, density, and other physicochemical properties of the oil under specific process conditions, thereby greatly decreasing efficiencies of oil production, transportation, and refining.1−4 It is therefore of great industrial and academic interest to understand the structural and environmental factors that influence the early stage association of such type of molecules that leads to their subsequent precipitation in organic solvent. Furthermore, due to widespread applications of PA compounds in such areas as organic thin film transistors,9 organic (opto)electronic devices and sensors,10,11 light harvesting arrays,12 pigment and laser dyes,13,14 fabrication of liquid © 2016 American Chemical Society

Received: October 15, 2015 Revised: March 6, 2016 Published: March 8, 2016 3516

DOI: 10.1021/acs.jpcb.5b10082 J. Phys. Chem. B 2016, 120, 3516−3526

Article

The Journal of Physical Chemistry B

Figure 1. Two-dimensional schematic representation of perylene bisimide based PA molecule (C5 Pe) (a) and different NA molecules (b) used in this simulation study. C5 Pe molecules incorporate a PA perylene core with two alkyl chains attached to one end and aliphatic carboxylic groups at the other end. All NA molecules possess the same number of polar carboxylic acid functional groups but different aliphatic/cycloaliphatic tails as the nonpolar part of the molecules.

cycloaliphatic tails as the nonpolar part of the molecules. The detailed methods of generating initial molecular geometries in the PDB/GRO file format, constructing and implementing the topology files, and constructing organic solvent molecules are described in our previous publications and were directly adopted in this study.20,21,25,26,29 In this work, four simulation boxes each with a dimension of 10 nm × 10 nm × 10 nm and containing a total of 18 PA molecules and 9 NA molecules were initially prepared. Each simulation box has three equally divided layers, and each layer contains nine molecules with the upper and lower layers of the simulation box containing PA molecules and the middle layer containing different types of NA molecules (see Supporting Information - S1). This kind of molecular-packing arrangements was employed to account for maximum contacts beween NA and PA molecules in the simulation systems, thus allowing us to examine the influence of NA molecules possessing different numbers of carbon (n) and rings (Z) on the aggregation of PA molecules in toluene (Figure 1b). The total number of toluene molecules added into each simulation box was within the range of 5623−5660 molecules, depending on the type of NA molecules used in the simulation box (see Supporting Information - S2). In order to discern whether toluene would have any effect on the nanoaggregation formation, we also performed an additional simulation containing no NA molecules in the system, and the data generated from this simulation are used as a baseline for the purpose of comparison. The mass concentration (wt %) of PA and NA molecules in our simulation boxes was within the range 2.2−2.4% (or ∼29.9 mmol/L) and 0.3−0.5% (or ∼14.9 mmol/ L), respectively. Furthermore, the carboxylic acid functional group of PA and NA molecules in the studied solvent (toluene) was kept protonated in all simulation systems. All MD simulations performed in this work were carried out using the GROMACS simulation package (version 4.0.5) with the GROMOS96 53a6 force field.35−38 For each system, energy minimization was first performed, followed by the isothermal− isobaric (NPT) ensemble simulation. Unless otherwise stated, the temperature and pressure of the simulation were kept constant at 298 K and 1 bar by a Nosé−Hoover thermostat39,40 and a Parrinello−Rahman barostat, respectvely.41 The isothermal compressibility of toluene (k(T) = 9.08 × 10−5 bar−1) was used for all the simulations.42 A cutoff distance of 1.4 nm

hydrogen bonding, interaction of alkyl groups that form hydrophobic pockets, etc.) could also contribute, to different extents, to the nanoaggregation.27,28 Our previous molecular dynamics (MD) simulation study revealed that, although the noncovalent attractive forces between the PA cores of the molecules play a major role in the nanoaggregate formation, various polar and nonpolar interaction sites (or moieties) on the PA molecules also exert a significant impact on initiating molecular association, the size of the nanoaggregates, 3D network formation, and the stacking mode of the individual molecules in the nanoaggregates.29,30 However, despite the significant importance of these polar and nonpolar interaction sites in nanoaggregate formation, how these active polar and nonpolar sites are influenced by various environmental factors, such as the presence of naphthenic acids (NAs), resins, aromaticity of the solvent, etc., remains poorly understood at the molecular level. Among all of these factors, of particular interest is the effect of NAs on nanoaggregate formation and growth in organic solvents. Although the amount of NAs in different oil fields varies considerably, recent studies have consistently shown increasingly higher contents of complex PA and NA molecules in several major oil field operations and water treatment plants.31−34 NAs are naturally occurring organic compounds with a generic chemical formula CnH2n+ZO2, where n indicates the number of carbons in the NA and Z specifies a homologous series (Z ranging from 0 to 12).31−34 NA molecules incorporate both polar and nonpolar functional groups in their molecular structures. Up to date, still very little is known about the possible interactions between NA and PA molecules in organic solvents. Evidently, the interaction between NA and PA molecules in the bulk oil phase can influence the selfassociation of PA molecules, thereby impacting the nanoaggregate formation and growth. The objective of this study is to investigate the influence of NAs on 3D assembly of selected PA molecules in toluene using molecular dynamics (MD) simulation.

2. METHODS 2.1. Simulation Details. The structures of the PA and NA molecules used in this study are shown in Figure 1. As can be seen, all NA molecules studied here possess the same polar carboxylic acid functional group but different aliphatic/ 3517


Article


Figure 2. Snapshots of the MD simulation cell containing either PA or PA and NA molecules in toluene. The initial molecular arrangement of PA and NA molecules in toluene is shown in simulation 1, and the final aggregated structures are depicted in simulations 2−6. PA is represented in black, while NA is represented in green. For clarity, toluene molecules are omitted.

cores are in a stacked state, indicating strong associations between PA cores. However, the association of PA molecules via side chains is also observed (e.g., red circle of simulation 2), indicating the formation of hydrogen bonds between the carboxylic acid groups of PA molecules. Moreover, the snapshots clearly revealed association between NA and PA molecules in the system. To quantitatively describe the role of NA molecules in the system, we first analyzed the nature of the stacked PA cores in the system. Demonstration of achieving dynamic equilibrium using quantities such as the velocity autocorrelation function (VACF) of PA and NA molecules is available in the Supporting Information, where all the data were generated using the last 5 ns of the simulation time (see Supporting Information - S5). 3.1. The Effect of NA Molecules on the π−π Stacked Structures of PA Molecules. To quantify the stacked structures of the PA molecules, we first calculated the radial distribution function (RDF), g(r), for the PA cores inside each simulated system. Then, the calculated RDF was further normalized with respect to the highest peak, which differs from the original RDF by a normalization factor, so that the generated RDF is confined between 0 and 1. The RDF data before normalization is available in the Supporting Information (see Supporting Information - S6) for easy comparison of different RDFs from the same system. As shown in Figure 3, the RDFs of the PA cores exhibit similar shapes that have multiple peaks located at distances of approximately 0.39, 0.77, and 1.14 nm from a reference PA core. This finding indicates that the basic structural arrangement between two PA cores during self-association is not greatly affected by the presence of NA molecules in the system. However, the corresponding peak heights are clearly different, suggesting that NA molecules do affect the prevalence of π−π stacking among the PA molecules. The results in Figure 3 show that, except for the case with NA-4 molecules, the peak heights of the RDFs of PA cores are higher

was used to compute the van der Waals interactions and shortrange electrostatic interactions. The long-range electrostatic interactions were treated using the particle mesh Ewald (PME) summation algorithm with periodic boundary conditions applied in the x, y, and z directions.43,44 The LINCS algorithm was used to constrain the bond length.45 The neighbor list was updated every 5 steps with a cutoff distance of 1.2 nm. A leapfrog Verlet algorithm implemented in GROMACS with a time step of 2 fs (fs) was used for the integration of Newton equations of motion. GROMACS built-in analytical tools and Visual Molecular Dynamics (VMD)46 were used to process, analyze, and visualize trajectories of all atoms and hence molecules in the system. The sensitivity and applicability of the united atom force field and simulation parameters described above were validated experimentaly in our previous studies on PA aggregation in toluene.20,21,25,29,47 Furthermore, our computed values on the density, self-diffusivity, and enthalpy of evaporation of toluene as well as their comparison with experimental results are available in the Supporting Information (see Supporting Information - S3).20,21,25,29,47 Therefore, given a sufficiently long period of simulation time (∼50 ns), we are confident that the force field and simulation conditions used in this study would reproduce reasonably well the experimental observations on aggregation of PA molecules in toluene at the specified temperature and pressure.

3. RESULTS AND DISCUSSION Figure 2 shows the snapshots of MD simulation boxes at the end of simulation time, which depicts the formation of nanoaggregates in the five systems (simulations 2−6). Detailed inspections of Figure 2 show that π−π stacking interaction is a dominant mode of interaction between neighboring PA cores (see Supporting Information - S4), consistent with our earlier studies.29,47,48 Particularly, in simulation 4, most of the PA 3518


Article


Table 1, the presence of NA molecules generally (except NA-4) increases the number of FLPS pairs in toluene, corresponding to a greater peak height observed in Figure 3. The second peak of the RDF is at the location of 0.77 nm away from the reference PA core, corresponding to the next shell of neighbors. At such a distance, we defined the π−π stacking as the second layer parallel stacking (SLPS) configuration, where the COM distance between the two PA cores ranges from 0.56 to 0.90 nm with cos θ larger than 0.9 (i.e., θ close to 0°). Together with the FLPS pair, we define the total stacking pair (TSP = FLPS + SLPS) as the π−π stacking configuration formed by two PA cores of COM distance less than 0.9 nm and cos θ larger than 0.9 (the first two prominent peaks in RDFs). The computed numbers of SLPS and TPS pairs are also given in Table 1. Clearly, NA molecules (except for NA-4) generally tend to enhance π−π stacking between PA molecules in toluene. Among the NA molecules studied in this work, the presence of NA-1 and NA-2 resulted in the largest number of FLPS pair and TPS pair. This trend is consistent with increasing first peak height observed in Figure 3. Compared to the first two peaks in each RDF (Figure 3), the third peak located at 1.14 nm is noticeably smaller. Therefore, such extended stacking is of less significance in our analysis and will not be discussed here. To determine whether NA molecules can affect the size and distribution of the π−π stacked structures, we further quantified the stacking type (e.g., dimer, trimer, tetramer, pentamer, etc.) of the PA molecules. Figure 4a shows the stacking types and the associated number of occurrence in each system. The varying degrees of distribution in the stacking type are clearly seen. For instance, PA molecules mixed with NA-3 molecules form π−π stacked dimer, trimer, tetramer, and pentamer, while for PA molecules mixed with NA-4 molecules only dimer and pentamer modes of stacking were seen. Further inspection of Figure 4a reveals that while the largest stacking type formed by the PA molecules without NA molecules is tetramer, a stacked PA pentamer was observed to form in all of the systems containing NA molecules, confirming again that NA molecules can enhance the stacking of PA molecules. However, compared with the other three types of NA molecules, a relatively large number of dimers and monomers formed in the system containing NA-4 molecules. This finding is consistent with the lower peak observed in the RDF (see Figure 3) and the smaller numbers of FLPS and TPS pairs shown in Table 1.

Figure 3. Normalized RDFs of PA cores. The inset shows a snapshot in the aggregate of simulation 2 (No NA), displaying the π−π stacking arrangements (t = 50 ns).

in the presence of NA molecules than in the absence of NA molecules (simulation 2), suggesting enhanced π−π stacking interactions between PA molecules in the presence of NA molecules, which is confirmed by a closer look at the aggregate structure (see Supporting Information - S4). To quantitatively compare the overall π−π stacking structures formed at each peak location (r = 0.39, 0.77, and 1.14 nm), we first define the π−π stacking formed at the first peak (r = 0.39 nm) as the first layer parallel stacking (FLPS) configuration. This configuration is quantified by two PA cores having the distance between their centers of mass (COM) less than 0.55 nm and the cosine of the angle (cos θ) between two PA ring planes larger than 0.9 (i.e., θ close to 0°). On the basis of the above distance and angle criteria, the number of FLPS pairs (first peak) in each system averaged over the last 5 ns of the 50 ns simulation is summarized in Table 1. As shown in Table 1. Number of π−π Stacking Pairs Formed in the Systems Averaged over the Last 5 ns of the 50 ns Simulation system

no NA

NA-1

NA-2

NA-3

NA-4

FLPS pairs SLPS pairs TPS pairs

11 2 13

17 5 22

15 4 19

12 2 14

10 3 13

Figure 4. (a) Type of π−π stacking and their occurrences in each system and (b) distribution (number fraction, Nf = Ni/NTotal) of stacking mode averaged over the last 5 ns of the 50 ns simulation. 3519


Article


Figure 5. Snapshots of the initial (0 ns) and final (100 ns) simulation cells. Parts a and c are the initial configurations of the simulation systems; parts b and d are the molecular configurations after 100 ns simulation time. C5 Pe is represented in black, while NA-2, in green. Solvent molecules are removed for clarity.

stacking is an offset twist parallel 3D assembly with intermediate overlapping area (see Figure 4b, inset), which is irrespective of the types of NA molecules added into the system. This finding implies that “perfect” face to face or cross stacking modes of PA molecules stacking in toluene are energetically less favorable. The results of the above structural quantifications indicate an insignificant effect of NA molecules on the overlapping area of two parallelly stacked PA molecules (i.e., FPA, CPA, or TPA). However, NA molecules can enhance the size of the π−π stacked structures, i.e., the number of PA molecules that can be involved in one stacked structure. This kind of observation is clearly visible, in particular for NA-1 and NA-2 molecules whose size is relatively small. These NA molecules are expected to affect the growth and characteristics of the stacking. To further quantify the effect of NA molecules on the size of π−π stacked structures, two additional simulations were performed for 100 ns as described in the following. Since the stacking of the PA molecules was increased mostly by NA-1 and NA-2 molecules (molecules of smaller sizes), we analyzed a single stacked structure (pentamer) shown in the inset of Figure 4a from a system of PA molecules mixed with NA-2 (simulation 4 in Figure 2). This stacked structure consisted of 5 PA (C5 Pe) molecules with 4 NA-2 molecules bounded at its active polar sites (−COOH group). The coordinates of these “NA-2 bound C5 Pe molecules” were extracted from the final trajectory in simulation 4 of Figure 2. Throughout the rest of this work, this structure is treated as a single pentamer and referred to as “pentamer C5 Pe”. Further simulation was conducted with eight identical pentamers arranged into a cubic box of 10 nm × 10 nm × 10 nm dimension (shown in Figure 5c). A reference simulation was also conducted in the same way by removing the four bounded NA-2 molecules from the active polar site of the PA pentamers (Figure 5a). Thus, the only difference between these two systems was the presence and absence of NA molecules bounded to the PA pentamers. After a 100 ns simulation run, the aggregates formed in these two systems are shown in Figure 5b and d.

The NA molecules studied here are differentiated by the nonpolar aliphatic part with the size and structure being the main differences between them. For example, NA-4 molecules possess a bulky cycloaliphatic tail in comparison to NA-1, NA2, and NA-3 molecules. Since small NA molecules possess less steric hindrance compared to large NA molecules when interacting with PA molecules, the reduced stacking of PA molecules observed in the presence of NA-4 molecules is most likely due to the large steric hindrance created by the bulky cycloaliphatic tail of the NA-4 molecule. We also investigated the stability and dynamic nature of these stacking types. To do so, the stacking types and the associated number of occurrence in each system were calculated as a function of time for one of our simulation systems (PA mixed with NA-1). The result (see Supporting Information - S7) shows that the distribution of different stacking types remains essentially the same after 40−50 ns simulation time, indicating that the parallel stacking mode of PA cores is reasonably stable. In order to qualitatively determine the effect of different NA molecules on the strength of π−π binding and 3D assembly of stacked PA molecules, the overlapping area of two PA cores that form an FLPS pair (see Figure 4b) was evaluated by calculating the minimum (“MIN”) distance between the nitrogen atoms of the stacked neighboring molecules (nitrogen−nitrogen distance). For an FLPS pair of the MIN between two parallel stacked PA cores less than 0.40 nm, a face to face parallel 3D assembly (FPA) is formed to feature the largest overlapping area. For an FLPS pair of MIN distance between 0.50 and 0.60 nm, PA cores form a cross parallel 3D assembly (CPA), where the two stacked molecules have the least overlapping area. Finally, for an FLPS pair of the MIN distance between 0.40 and 0.50 nm, PA cores form a twist offset parallel 3D assembly (TOPA) with intermediate overlapping area of PA cores. The distribution (number fraction, Nf) of each type of assembly was calculated by dividing the number of stacking pairs (Ni) in each assembly mode by the total number of FLPS pairs (NTotal). As shown in Figure 4b, the dominant mode of 3520


Article


the side chains of the PA molecules were first investigated. As the PA molecule used in this study incorporates a PA perylene core with two alkyl chains attached to one end and an aliphatic carboxylic acid functional group to the other end (Figure 1a), two dominant interactions exist between NA molecules and the exteriors of the stacked PA molecules: (1) hydrogen bonding interaction between the carboxylic acid functional groups (head to head interaction) and (2) the associations of nonpolar aliphatic tails (tail to tail interaction) of interacting PA molecules. To quantify these interactions, the molecules were first divided into polar (−COOH) and nonpolar (aliphatic or cycloaliphatic tail) sites for both NA and PA molecules (see Supporting Information - S8), followed by the calculation of the MIN distance between the polar sites of PA and NA molecules as well as between their nonpolar sites. The two sites of the molecules are defined to form a “contact” if the MIN distance between them is less than 0.5 nm. For example, if the MIN distance between two carboxylic functional groups of the PA and NA molecules is less than 0.5 nm, a head to head contact is counted. Using this method, the total number of contacts in each simulation can be obtained. Figure 6 shows the total number of contacts as a function of simulation time for each system shown in Figure 2 (t = 50 ns). It can be seen from Figure 6 that irrespective of the NA molecules, the results show a sharp monotonic increase in the total number of contacts between the polar heads (−COOH) of the PA and NA molecules in comparison to the nonpolar tail to tail interactions between these molecules. It is therefore evident that the dominant interactions between NA and PA molecules are the strong head to head hydrogen bonding interactions between the carboxylic acid functional groups of the molecules. To confirm the formation of hydrogen bonding between the carboxyl functional groups of PA and NA molecules, the normalized RDF of the oxygen atom of the carboxyl functional group (O−COOH) in C5 Pe and the hydrogen atom of the carboxyl functional group (H−COOH) in NA was plotted in Figure 7a. The original RDF data is available in the Supporting Information (see Supporting Information - S9). Not surprisingly, the RDFs show very similar trends for all the systems simulated, suggesting similar binding mechanisms of different NA molecules in the simulated systems. The sharp peak located at an approximately 0.18 nm distance away from the reference oxygen atom in the normalized RDF is consistent with the

The results in Figure 5 show the formation of larger stacked structures of decamers in the system containing NA-2 molecules than in the reference system without NA-2 molecules. Adopting the same distance criteria as discussed above to compute the stacking distribution in the system, we performed the analysis every 20 ns (after 30 ns of the simulation time) to quantify the size of the stacked structure of C5 Pe molecules for these two simulated systems. The results in Table 2 reveal that after 100 ns of simulation time, almost Table 2. Number of π−π Stacked Structures with and without NA-2 in Bulk Toluene without NA-2

with NA-2

time (ns)

number of π−π pentamer

number of π−π decamer

number of π−π pentamer

number of π−π decamer

initial 35−40 55−60 75−80 95−100

8 8 8 8 8

0 0 0 0 0

8 6 6 6 4

0 1 1 1 2

50% of the initial pentamers formed decamers in the presence of NA-2 (Figure 5d). In contrast, all of the pentamers remained in their original states for the system without NA-2 bound to its pentamers. It can therefore be concluded that the presence of NA-2 molecules can enhance the π−π stacking of the PA molecules. Considering the similar effect of NA-1 and NA-3 molecules on the stacking of PA cores, similar aggregation behavior is expected for these NA molecules (see Figures 3 and 4 as well as Table 1). It should be noted that the selection of the initial π−π stacked structure, such as dimer, trimer, etc., might influence the size of the largest stacked structures (decamers) that could form. However, the main conclusion on the effect of NA molecules on the stacking of PA molecules would remain the same. 3.2. Effect of NA Molecules on the Intermolecular Hydrogen Bonding of PA Molecules. To clarify the underlying mechanism of the observed role of NA molecules, it is crucial to examine the binding of NA molecules with PA molecules in the system. As shown in Figure 5c and d, NA molecules are mainly bonded to the exteriors of the stacked PA molecules. Thus, the interactions between NA molecules and

Figure 6. Number of tail to tail contacts between the nonpolar segments of the NA and PA molecules (a) and head to head contacts between the carboxylic acid groups of NA and PA molecules (b) as a function of simulation time. The insets show the schematic diagram depicting these interactions. 3521


Article


Figure 7. (a) The normalized RDF between the oxygen atom of C5 Pe (O− COOH of C5 Pe) and the hydrogen atom of the NA (H− COOH of NA), and (b) between hydrogen and oxygen atoms of carboxylic acid side chains of the PA molecules.

reported values in the literature of the hydrogen bond distance between two carboxylic acid functional groups in the RDF,20,21,29 confirming the formation of strong hydrogen bonding between NA−H···O−PA in the current systems. However, the sharpness of the first peak is different for different types of NA molecules, emphasizing the effect the aliphatic structure and size of NA molecules have on their binding to the PA molecules in the system (see Supporting Information S10). In fact, our normalized RDF results reveal that among all the NA molecules simulated in this study, NA-2 molecules have the highest tendency to form the hydrogen bond with the PA molecules in toluene. With the formation of hydrogen bonds between NA and PA molecules, the number of hydrogen bonds that can form between PA molecules is anticipated to reduce. To show this, the normalized RDF between the carboxylic acid groups in the side chains of the PA molecules was calculated with the results shown in Figure 7b. As anticipated, the systems containing the NA molecules exhibit a smaller first peak as compared to the case without NA molecules (i.e., reference simulation 2), indicating a suppressed degree of the PA molecular association through hydrogen bonding due to the presence of NA molecules, more so for NA molecules of small sizes. The original RDF data is available in the Supporting Information (see Supporting Information - S11). To further evaluate the formation and stability of the hydrogen bond interactions among different species in the simulated systems, the number of hydrogen bonds formed and their lifetimes were calculated. For calculating the number of hydrogen bonds, a geometric criterion of donor−acceptor distance rDA ≤ 3.5 Å and ∠acceptor−donor−hydrogen θADH ≤ 30° was used.49 The lifetime of the hydrogen bonds (τHB) was calculated from the hydrogen bond time autocorrelation function, CHB(t), at 298 K. The function CHB(t) describes the probability that a particular pair of donor and acceptor is hydrogen bonded at time t, provided that it was hydrogen bonded at reference time t = 0, and even if the bond broke at certain intermediate time.35,49 Therefore, the integration of the hydrogen bond time autocorrelation function provides an estimate of the average hydrogen bond lifetime. The numbers of hydrogen bonds and hydrogen bond lifetimes are given in Table 3. Clearly, compared with the reference simulation, the number of hydrogen bonds among PA molecules decreased upon binding of NA molecules to the PA molecules, shown as the reduced peak in the RDF of side chains in Figure 7b. Furthermore, the number of hydrogen bonds formed between NA and PA molecules depends on the type of NA molecules present. Small molecules of compact structures such as NA-2

Table 3. Average Number of Hydrogen Bonds between Various Simulated Molecules and the Hydrogen Bond Lifetime (Averaged over 50 ns Simulation Time) between PA and NA Molecules τHB (ns)

mean number of hydrogen bonds system no NA NA-1 NA-2 NA-3 NA-4

PA−PA 8.86 7.46 6.55 6.07 7.66

± ± ± ± ±

1.14 0.89 1.16 0.82 0.93

NA−PA 11.09 15.99 11.64 9.69

± ± ± ±

1.10 1.12 0.61 0.62

NA−NA

NA−PA

± ± ± ±

5.10 5.37 4.86 6.20

3.14 0.76 1.95 1.32

0.47 0.20 0.32 0.61

feature a higher tendency to form hydrogen bonds with PA molecules, leading to the formation of ∼16 hydrogen bonds. The average hydrogen bond lifetime between PA and NA molecules in toluene ranges between 4.9 and 6.2 ns (averaged over the entire 50 ns simulation time), indicating that the stability of the hydrogen bond between PA and different NA molecules is comparable despite the different numbers of hydrogen bonds formed. As discussed above, NA molecules can form hydrogen bonds with the PA molecules, thus depressing the formation of hydrogen bonds between the PA molecules. The effect of NA molecules on the hydrogen bond can be further verified by computing the number of hydrogen bonds between the PA molecules as a function of time for the two additional simulations (pentamer simulations); this is shown in Figure 8. Since molecules in these simulation systems are in an aggregated form, there exist hydrogen bonds at the beginning of the simulation time (t = 0 ns). Furthermore, as the difference between these two additional simulations is only the presence of bounded NA molecules in these PA pentamers, the numbers of hydrogen bonds among PA molecules in both simulations are initially equal. However, after starting the simulation, the calculated number of hydrogen bonds for these two simulation systems shows a dramatic change, as can be seen from Figure 8. The number of hydrogen bonds formed between the PA molecules in the system without NA molecules shows a sharp increase with time, and the number of the hydrogen bonds is doubled (i.e., increased by ∼100%) at the end of the simulation time of 100 ns. In contrast, the number of hydrogen bonds formed between the PA molecules in the system with NA molecules remains almost constant with time, confirming that formation of hydrogen bonds between PA molecules is greatly suppressed by the NA molecules in the system. 3522


Article


vicinity of the PA molecules (see Figure 9). The RDF results also indicate that more toluene molecules are displaced from the vicinity of the PA molecules in the presence of NA-2 molecules, which suggests the prevailing effect of small NA molecules over the π−π stacking of PA molecules in the system. This trend is consistent with the cumulative number of toluene molecules around the PA core and hydrogen bonding sites (see Supporting Information S12), suggesting that the hydrogen bonds formed between NA and PA molecules reduce the number of solvent (toluene) molecules around the solute (PA) molecules due to the increase in π−π stacked PA structures in the system. From the above discussion, it is of great importance to point out that, while NA molecules enhanced the population of π−π stacking configurations and displaced more toluene molecules from the vicinity of the PA molecules, such an effect does not necessarily translate to a smaller overall size of the aggregates than the case in the absence of small NA molecules. To have a complete picture, other aggregation mechanisms such as hydrogen bonding need to be considered. In two additional simulations, for instance, the largest aggregates determined by a single linkage algorithm involved 35 and 40 PA molecules, respectively, for systems with and without NA-2 molecules. Clearly, the aggregates formed in the system without NA-2 molecules were even larger. In the presence of NA-2 molecules, however, two decamers were involved in the largest aggregate. Without NA-2 molecules, the largest aggregate has only stacked pentamers linked together through other types of interactions, i.e., hydrogen bonding in this particular case (see Figure 10). Therefore, the major role played by NA-2 molecules during the aggregation of PA molecules is its influence on aggregation mechanisms to lead to distinct structures, rather than just simply changing the size of the aggregates. 3.4. Thermodynamic Implications of PA and NA Interactions in Toluene. Our previous studies20,21,29,47 showed that in addition to aromatic π−π stacking, the intermolecular hydrogen bonding between PA molecules plays a crucial role in the formation and growth of nanoaggregates. While PA molecules can still aggregate via π−π stacking in the absence of NA molecules, hydrogen bonding interactions among PA molecules are clearly noticeable. Indeed, the snapshots in Figures 2 and 5 actually show the hydrogen bonding linkage between PA molecules, which is consistent with the results from our previous studies.20,21,29,47 However, as shown in the previous section, the presence of NA molecules clearly blocks the hydrogen

Figure 8. Number of hydrogen bonds among PA molecules as a function of simulation time with and without NA bounded to the PA pentamers of the system. The insets show snapshots of the initial structures for the two simulated systems.

3.3. Solvent Distribution around PA Nanoaggregates. The intermolecular binding of NA to PA molecules not only depresses the hydrogen bonding between PA molecules but also affects the solvent distribution around the PA aggregates. To illustrate this effect, the local distribution of toluene solvent molecules around the PA core and carboxylic acid functional groups of the PA molecules was calculated and compared with each other in detail. First of all, it is worth mentioning that toluene molecules can interact with PA molecules via noncovalent PA core−toluene π−π interactions and to some extent the weak hydrogen bonding between the carboxylic acid functional group and toluene. In fact, these two types of intermolecular interactions make toluene a “good” solvent for PA molecules.8,20−24 Hence, two types of RDFsthe RDFs of the toluene molecules around the PA core and around the hydrogen bonding sites of the PA moleculeswere calculated and the results are shown in parts a and b of Figure 9, respectively. It can be seen that the presence of NA molecules in the system creates a favorable local environment for the overlap of PA cores (i.e., π−π stacking growth) by depressing the hydrogen bonding interactions between PA molecules, which forces some toluene molecules to migrate away from the

Figure 9. RDFs of toluene molecules near the oxygen atom of the carboxylic acid functional group of the PA molecules (a) and near a PA core (b), averaged over the last 5 ns of the 50 ns simulation. 3523


Article


Figure 10. Schematic diagram showing the possible association mechanism of stacked PA molecules in the absence and presence of small NA molecules such as NA-2 in toluene. C5 Pe is represented in black, while NA-2, in green.

4. CONCLUSIONS The effect of NA molecules on the configurations of the nanoaggregates formed by PA molecules in toluene was studied using MD simulation. The results from this study show that while NA molecules do not significantly affect the preferred 3D assembly mode of two neighboring parallel stacked PA molecules, they can form hydrogen bonding with PA molecules, reducing the hydrogen bonds between the PA molecules. Therefore, the size of π−π stacked structures (i.e., the number of PA molecules involved in one stacked structure) was enhanced to compensate for the loss of hydrogen bond interactions among the PA molecules. The presence of NA molecules can reduce the number of solvent (toluene) molecules around the PA molecules, in particular around the hydrogen bonding sites of the PA molecules. Subtle differences were observed among different NA molecules. Depending on the size and type of the NA molecules present in the system, the stacked structures feature various distributions in size. Larger NA molecules such as NA-4 can introduce significant steric hindrance when they are bounded via hydrogen bonds to the exterior of the PA molecules. As a result, they tend to favor the formation of dimers when compared with smaller NA-1 and NA-2 molecules which promotes the formation of more trimers, tetramers, and so on in the system. It is clear from the current study that to accurately predict the physicochemical behavior of asphaltene molecules in crude oil, the critical impact of other oil components such as NAs needs to be carefully considered.

bonding among PA molecules, which reduces the possibility of forming a linkage between the carboxylic functional groups of the PA molecules and therefore increases the π−π stacking configurations. Actually, the configuration of aggregated structures is driven by the direction of free energy minimization. While π−π binding strength may decrease free energy to a larger extent in comparison with other types of interactions, such as hydrogen bonding, ordering of the PA cores decreases the entropy and does the opposite. As a result, the favored configuration formed may not necessarily be the one of higher binding energy, i.e., π−π stacking, when the effect of entropy is taken into consideration. Therefore, in the absence of NA molecules, PA molecules will aggregate in such a way to balance all types of interactions, such as π−π stacking and hydrogen bonding, to minimize the free energy of the system. For instance, Sedghi et al.50 reported that the free energy of dimerization of PA molecules in aliphatic solvent via π−π stacking is −14 kJ/mol while the corresponding value to include hydrogen bonding is −17.7 kJ/mol. These values are based on the formation of dimers. For the nanoaggregation of PA molecules, especially, for the further growth of nanoaggregates, hydrogen bonding may play a more significant role. Since the hydrogen bonding between PA molecules is depressed in the presence of NA molecules, the parallel stacking becomes more predominant in the aggregation of PA molecules. Therefore, as shown in Figure 10, two PA pentamers can further stack into a one-dimensional rod-like decamer in the presence of NA-2 molecules, while in the reference simulation, side-by-side linkage was observed during the further association of pentamers. Finally, it is of great importance to point out here that the above effect of NA molecules can only exist when the concentration of the PA molecules is high enough such that aggregation will certainly happen. In the case that PA molecules do not aggregate, such an effect may not be as significant.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b10082. Figures and tables showing the initial simulation setup, RDF before normalization, validation of the system 3524


Article


■

simulation force field, and snapshots of the simulation box that tend to show the aggregated molecules and the partitioning of polar and nonpolar active sites of the molecules (PDF)

(16) Fischer, M. K. R.; Kaiser, T. E.; Würthner, F.; Bäuerle, P. Dendritic oligothiophene-perylene bisimide hybrids: synthesis, optical and electrochemical properties. J. Mater. Chem. 2009, 19, 1129−1141. (17) Reghu, R. R.; Bisoyi, H. K.; Grazulevicius, J. V.; Anjukandi, P.; Gaidelis, V.; Jankauskas, V. Air stable electron-transporting and ambipolar bay substituted perylene bisimides. J. Mater. Chem. 2011, 21, 7811−7819. (18) Ilhan, F.; Tyson, D. S.; Stasko, D. J.; Kirschbaum, K.; Meador, M. A. Twisted, Z-Shaped Perylene Bisimide. J. Am. Chem. Soc. 2006, 128, 702−703. (19) Rybtchinski, B. Adaptive Supramolecular Nanomaterials Based on Strong Noncovalent Interactions. ACS Nano 2011, 5, 6791−6818. (20) Kuznicki, T.; Masliyah, J. H.; Bhattacharjee, S. Molecular Dynamics Study of Model Molecules Resembling Asphaltene-like Structures in Aqueous Organic Solvent Systems. Energy Fuels 2008, 22, 2379−2389. (21) Kuznicki, T.; Masliyah, J. H.; Bhattacharjee, S. Aggregation and Partitioning of Model Asphaltenes at Toluene-Water Interfaces: Molecular Dynamics Simulations. Energy Fuels 2009, 23, 5027−5035. (22) Boek, E. S.; Yakovlev, D. S.; Headen, T. F. Quantitative Molecular Representation of Asphaltenes and Molecular Dynamics Simulation of Their Aggregations. Energy Fuels 2009, 23, 1209−1219. (23) Headen, T. F.; Boek, E. S.; Skipper, N. T. Evidence for Asphaltene Nanoaggregation in Toluene and Heptane from Molecular Dynamics Simulations. Energy Fuels 2009, 23, 1220−1229. (24) McLean, L. D.; Kilpatrick, P. K. Effects of Asphaltene Aggregation and Resin-Asphaltene Interactions on Stability of Emulsions Generated from Water-Heptane-Toluene-Resin-Asphaltene Mixtures. J. Colloid Interface Sci. 1997, 196, 23−34. (25) Jian, C.; Tang, T. One-dimensional Self-assembly of Polyaromatic Compounds Revealed by Molecular Dynamics Simulations. J. Phys. Chem. B 2014, 118, 12772−12780. (26) Jian, C.; Tang, T. Molecular Dynamics Simulations Reveal Inhomogeneity-Enhanced Stacking of Violanthrone-78-Based Polyaromatic Compounds in n-Heptane−Toluene Mixtures. J. Phys. Chem. B 2015, 119, 8660−8668. (27) Sjöblom, J.; Soderlund, H.; Lindblad, S.; Johansen, E. J.; Skjarvo, I. M. Water−in− Crude Oil Emulsions from the Norwegian Continental Shelf. Part II. Chemical Destabilization and Interfacial Tensions. Colloid Polym. Sci. 1990, 268, 389−398. (28) Gray, M. R.; Tykwinski, R. R.; Stryker, J. M.; Tan, X. Supramolecular Assembly Model for Aggregation of Petroleum Asphaltenes. Energy Fuels 2011, 25, 3125−3134. (29) Teklebrhan, R. B.; Ge, L.; Bhattacharjee, S.; Xu, Z.; Sjöblom, J. Probing Structure-Nanoaggregation Relations of Polyaromatic Surfactants: a Molecular Dynamics Simulation and Dynamic Light Scattering Study. J. Phys. Chem. B 2012, 116, 5907−5918. (30) Jian, C.; Tang, T.; Bhattacharjee, S. Probing the Effect of Side Chain Length on the Aggregation of a Model Asphaltene Using Molecular Dynamics Simulations. Energy Fuels 2013, 27, 2057−2067. (31) Auflem, I. H.; Havre, T. E.; Sjöblom, J. Near-IR Study on the Dispersive Effects of Amphiphiles and Naphthenic Acids on Asphaltenes in Model Heptane-Toluene Mixtures. Colloid Polym. Sci. 2002, 280, 695−700. (32) Havre, T. E.; Sjöblom, J.; Vindstad, J. E. Oil-Water Partitioning and Interfacial Behaviour of Naphthenic Acids. J. Dispersion Sci. Technol. 2003, 24, 789−801. (33) Havre, T. E.; Ese, M. H.; Sjöblom, J.; Blokhus, A. M. Langmuir Films of Naphthenic Acids at Different pH and Electrolyte Concentrations. Colloid Polym. Sci. 2002, 280, 647−652. (34) Knudsen, D. K.; Simon, S.; Qassym, L.; Gao, B.; Sjöblom, J. Mixed Micelles of Tetrameric Acids and Naphthenic Acids in Water. Energy Fuels 2014, 28, 4469−4479. (35) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701−1718. (36) Berendsen, H.; van der Spoel, D.; Drunen, V. R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge the financial support from the Canada Research Chairs (CRC) program and the NSERCIndustrial Research Chair in Oil Sands Engineering. We are also grateful for the Oil Sands and Coal Interfacial Engineering Facility (OSCIEF) for assistance and computational resources. Financial support through the Petromaks program (Norwegian Research Council) and JIP 1 (Ugelstad Laboratory) and National Natural Science Foundation of China (Grant 21333005) is gratefully acknowledged. We would also like to express our special thanks to Dr. Subir Bhattacharjee for his active involvement in the project and insightful discussions on the scientific matter of the subject.

■

REFERENCES

(1) Speight, J. G. Petroleum Asphaltenes - Part 1: Asphaltenes, Resins and the Structure of Petroleum. Oil Gas Sci. Technol. 2004, 59, 467− 477. (2) Rogel, E. Studies on Asphaltene Aggregation via Computational Chemistry. Colloids Surf., A 1995, 104, 85−93. (3) Pierre, C.; Barre, L.; Pina, A.; Moan, M. Composition and Heavy Oil Rheology. Oil Gas Sci. Technol. 2004, 59, 489−501. (4) Speight, J. G.; Long, R. B. The Concept of Asphaltenes Revisited. Fuel Sci. Technol. Int. 1996, 14, 1−12. (5) Spiecker, P. M.; Gawrys, K. L.; Trail, C. B.; Kilpatrick, P. K. Effects of Petroleum Resins on Asphaltene Aggregation and Water-inOil Emulsion Formation. Colloids Surf., A 2003, 220, 9−27. (6) Marshall, G. A.; Rodgers, P. R. Petroleomics: The Next Grand Challenge for Chemical Analysis. Acc. Chem. Res. 2004, 37, 53−59. (7) Marshall, A. G.; Blakney, G. T.; Beu, S. C.; Hendrickson, C. L.; McKenna, A. M.; Purcell, J. M.; Rodgers, R. P.; Xian, F. Petroleomics: A Test Bed for Ultra-High-Resolution Fourier Transform Ion Cyclotron Resonance Mass Spectrometry. Eur. Mass Spectrom. 2010, 16, 367−71. (8) Mullins, O. C. The Modified Yen Model. Energy Fuels 2010, 24, 2179−2207. (9) Facchetti, A. π-Conjugated Polymers for Organic Electronics and Photovoltaic Cell Applications. Chem. Mater. 2011, 23, 733−758. (10) Huang, Y.; Quan, B.; Wei, Z.; Liu, G.; Sun, L. Self-Assembled Organic Functional Nanotubes and Nanorods and Their Sensory Properties. J. Phys. Chem. C 2009, 113, 3929−3933. (11) Li, C.; Liu, M.; Pschirer, N. G.; Baumgarten, M.; Müllen, K. Polyphenylene-Based Materials for Organic Photovoltaics. Chem. Rev. 2010, 110 (94), 6817−6855. (12) Hippius, C.; Schlosser, F.; Vysotsky, O. M.; Böhmer, V.; Würthner, F. Energy Transfer in Calixarene-Based Cofacial-Positioned Perylene Bisimide Arrays. J. Am. Chem. Soc. 2006, 128, 3870−3871. (13) Law, K.-Y. Organic photoconductive materials: recent trends and developments. Chem. Rev. 1993, 93, 449−486. (14) Díaz-García, M. A.; Calzado, E. M.; Villalvilla, J. M.; Boj, P. G.; Quintana, J. A.; Céspedes-Guirao, F. J.; Fernández-Lázaro, F.; SastreSantos, A. Synth. Met. 2009, 159, 2293−2295. (15) Mizoshita, N.; Tani, T.; Inagaki, S. Highly Conductive Organosilica Hybrid Films Prepared from a Liquid-Crystal Perylene Bisimide Precursor. Adv. Funct. Mater. 2011, 21, 3291−3296. 3525


Article

The Journal of Physical Chemistry B (37) Lindahl, E.; Hess, B.; van der Spoel, D. GROMACS 3.0: A Package for Molecular Simulation and Trajectory Analysis. J. Mol. Model. 2001, 7, 306−317. (38) Oostenbrink, C.; Villa, A.; Mark, A.; van Gunsteren, W. A Biomolecular Force Field Based on the Free Enthalpy of Hydration and Solvation: The GROMOS Force-Field Parameter Sets 53A5 and 53A6. J. Comput. Chem. 2004, 25, 1656−1676. (39) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (40) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (41) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (42) Knovel Critical Tables, 2nd ed.; Knovel: New York, 2008; online: http://www.knovel.com/web/portal/browse/display?_EXT_ KNOVEL_DISPLAY_bookid=761&VerticalID=0 (accessed on 7 June 2014). (43) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8592. (44) van der Spoel, D.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A.; Feenstra, K. The GROMACS Manual (version 3.3); Groningen University: Groningen, The Netherlands, 2004. (45) Hess, B.; Bekker, H.; Berendsen, C. J. H.; Fraaije, M. E. G. J. A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (46) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (47) Teklebrhan, R. B.; Ge, L.; Bhattacharjee, S.; Xu, Z.; Sjöblom, J. Initial Partition and Aggregation of Uncharged Polyaromatic Molecules at the Oil-Water Interface: A Molecular Dynamics Simulation Study. J. Phys. Chem. B 2014, 118, 1040−1051. (48) Jian, C.; Tang, T.; Bhattacharjee, S. Molecular Dynamics Investigation on the Aggregation of Violanthrone78-based Model Asphaltenes in Toluene. Energy Fuels 2014, 28, 3604−3613. (49) Jana, B.; Pal, S.; Bagchi, B. Hydration Dynamics of Protein Molecules in Aqueous Solution: Unity Among Diversity. J. Chem. Sci. 2012, 124, 317−325. (50) Sedghi, M.; Goual, L.; Welch, W.; Kubelka, J. Effect of Asphaltene Structure on Association and Aggregation Using Molecular Dynamics. J. Phys. Chem. B 2013, 117, 5765−5776.

3526


Role of Naphthenic Acids in Controlling Self-Aggregation of a

Recommend Documents