Toward Predictive Molecular Dynamics Simulations of Asphaltenes in Toluene and Heptane.

The conventional definition of asphaltenes is based on their solubility in toluene and their insolubility in heptane. We have utilized this definition to study the influence of partial charge parametrization on the aggregation behavior of asphaltenes using classical atomistic molecular dynamics simulations performed on the microsecond time scale. Under consideration here are toluene- and heptane-based systems with different partial charges parametrized using the general AMBER force field (GAFF). Systems with standard GAFF partial charges calculated by the AM1-BCC and HF/6-31G*(RESP) methods were simulated alongside systems without partial charges. The partial charges implemented differ in terms of the resulting electrical negativity of the asphaltene polyaromatic core, with the AM1-BCC method giving the greatest magnitude of the total core charge. Based on our analysis of the molecular relaxation and orientation, and on the aggregation behavior of asphaltenes in toluene and heptane, we proposed to use the partial charges obtained by the AM1-BCC method for the study of asphaltene aggregates. A good agreement with available experimental data was observed on the sizes of the aggregates, their fractal dimensions, and the solvent entrainment for the model asphaltenes in toluene and heptane. From the results obtained, we conclude that for a better predictive ability, simulation parameters must be carefully chosen, with particular attention paid to the partial charges owing to their influence on the electrical negativity of the asphaltene core and on the asphaltenes aggregation.

Introduction

Asphaltenes, the heaviest and the most polar fractions of oil, are conventionally defined as a class of compounds that are soluble in toluene and insoluble in heptane.[1] They present many industrial challenges for oil production, transportation, and refining.[2] Among these challenges are increases in density and viscosity, stabilization of water–oil emulsions, clogging of wells and pipelines, and equipment fouling.[3] The tendency of asphaltenes to self-assemble and form aggregates has become the main physical factor underlying their negative impact.[1,4,5] Therefore, investigating interactions between asphaltenes and their aggregation behavior at the molecular level is of great importance.

Significant progress has already been made toward identifying the driving forces behind asphaltene aggregation. The current assumption is that the aggregation of asphaltenes is principally governed by π–π interactions, the interplay between Coulomb interactions and dominating dispersion forces.[6−9] Nevertheless, it remains difficult to study the aggregation experimentally because of the nanometer scale of individual asphaltene molecules and even of their aggregates.[1] A limited number of reported experimental studies have focused on aggregation at the nanoscale.[5,10−22] Therefore, it has become a promising line of inquiry to utilize computer simulation methods, such as quantum-chemical calculations, and classical atomistic and coarse-grained simulations, along with experimental research, to describe asphaltene aggregation.[23]

To be predictive, any simulations being run must capture the relevant physics and chemistry of asphaltenes. As they are based on the solution of the Schrödinger equation, quantum-chemical calculations allow for the most complete description of interactions between asphaltenes. However, at this level of theory, it becomes quite difficult to study the aggregation of more than three asphaltene molecules, due to the computational expense involved. Moreover, asphaltene aggregates usually include about 3–10 molecules.[19,21,24,25] Far fewer physical details are captured in atomistic and coarse-grained simulations, since quantum effects and electronic degrees of freedom are omitted. For such simulations, the main challenges are the time scales and the description of the interactions between the asphaltenes themselves.[23,26] Currently, there is a demand to perform simulations over microseconds not only to reach an equilibrium state but also to obtain sufficient statistics.[26−28] The interactions in the atomistic and coarse-grained simulations are taken into account using sets of the potential energy functions, so-called force fields. Existing force fields provide an efficient framework to simulate tens or even hundreds of thousands of atoms over microseconds and to study the aggregation of asphaltenes.[23] Recent developments in atomistic simulations utilizing polarizable force fields may provide the means to accurately describe the interactions between asphaltenes.[29] However, the accuracy and inherent complexity of these simulations limit their ability to attain the required simulation time scales because of the increase in the computational power required.[30,31] In classical atomistic simulations, the intermolecular contributions to the system potential energy mainly include the Lennard-Jones and Coulomb terms.[32] The former is defined by the short-range repulsion and dispersion parts; this is the main term used to model π–π interactions.[33−35] The latter term describes the Coulomb interactions between the atomic partial charges when, for example, π–π interactions occur.[7,8,33,36] The parametrization of the partial charges is, however, a highly complex task that can be viewed as a mean to vary the predictive abilities of simulations.[37,38] Usually, the standard partial charges supplied with the available atomistic force fields are utilized.[23] The use of quantum-chemical calculations to obtain partial charges is a way to represent the chemical features of the asphaltenes in the atomistic simulations.[23,39,40] The replacement of groups of asphaltene atoms with “grains” is a way to increase the available time and length scales in comparison with simulations based on the atomistic models.[41−45] The coarse-graining procedure significantly simplifies the chemistry of the asphaltenes and their interactions. For example, atomic partial charges and the Coulomb interactions between them are often omitted in coarse-grained simulations, while keeping only the Lennard-Jones contribution to model π–π interactions.[41−44] However, taking partial charges into account can be important for the reproduction of the stacked geometry of asphaltene aggregates, even at the coarse-grained level.[45]

Overall, to date, no consensus has been reached on the influence and role of partial charge parametrization in the aggregation of asphaltenes in computer simulations. The aim of this study was to address this issue. To this end, the aggregation behavior of the asphaltenes in toluene and heptane was analyzed using classical atomistic molecular dynamics simulations over a time scale of several microseconds. Model systems with different partial charges were considered using the general AMBER force field (GAFF).

The rest of the paper is organized as follows: in the following section, we will set out our analysis of the mutual arrangement of asphaltenes, the interactions between asphaltenes, the sizes and the fractal dimensions of the asphaltene aggregates, and the solvent entrainment inside the aggregates. We will then summarize our conclusions, and, finally, we will provide the description of the models studied and the details of our simulations.

Results and Discussion

Mutual Arrangement of Asphaltenes

Let us start by analyzing the mutual arrangement of the asphaltenes in the simulated systems. To this end, the autocorrelation function Ppair(Δt) of the centers-of-mass distances between the pairs of asphaltenes was calculated aswhere d(t) and d(t + Δt) denote the centers-of-mass distances between two asphaltenes at t and t + Δt simulation times, respectively; ⟨···⟩ is the averaging over the time intervals and over all of the asphaltenes in a system; and ⟨d⟩2 and ⟨d2⟩ are the square of the average centers-of-mass distance and the mean-square center-of-mass distance for the pair of asphaltenes, respectively.

The dependence of the autocorrelation function Ppair(Δt) for the systems under consideration is presented in Figure . Each Ppair(Δt) dependence was calculated using 5 μs long trajectories.

Figure 1

Autocorrelation function Ppair(Δt) for asphaltenes in toluene- and heptane-based systems (a) without partial charges, (b) with partial charges calculated by the HF/6-31G*(RESP) method, and (c) with partial charges calculated by the AM1-BCC method obtained using 5 μs long simulation trajectories.

As can be seen in Figure , the change in the mutual arrangement of asphaltenes was characterized by two relaxation processes in all of the systems studied: fast and slow. The characteristic times of the relaxation processes observed were evaluated using the sum of two exponentials; see the Supporting Information (Table s2). The fast process can be attributed to the fluctuations of the asphaltene positions, while the slow process may indicate the presence of asphaltene aggregation. Notably, the choice of partial charges significantly influenced the relaxation in toluene and heptane. In the case of the partial charges calculated by the HF/6-31G*(RESP) method, the relaxation of the Ppair(Δt) autocorrelation function for the asphaltenes was very similar in both solvents. In the toluene- and heptane-based systems without partial charges, the fast processes decayed at different times, whereas the slow processes appeared to be similar. The use of partial charges calculated by the AM1-BCC method revealed a drastically different trend: both relaxation processes were much slower in heptane than in toluene. This may indicate that asphaltenes aggregate more strongly in heptane than in toluene. It is worth noting that only the partial charges obtained using the AM1-BCC method led to a prominent difference in the relaxation of asphaltenes in toluene and heptane.

The slow process can take up to ∼250 ns of simulation time, depending on the partial charges; see Table s2. This estimate showed that on the scale of a few hundred nanoseconds the asphaltenes lose their correlations (memory). Therefore, to ensure that equilibration was reached in the systems, the first 300 ns were defined as the equilibration time. The production run time was deemed to be more than 4 μs, that is, 1 order of magnitude greater than the slowest relaxation process in the systems. Thus, we believe that such long simulations are necessary to carefully equilibrate the systems and to obtain sufficient statistics.

Interactions of Asphaltenes

We studied the interactions of the asphaltenes in toluene and heptane, depending on the choice of the partial charges. The pair correlation function[46]g(r) between the asphaltene centers of masses has been analyzed aswhere N0 = 50 is the number of asphaltenes in the systems, ⟨ρ⟩ is the asphaltene number density, δ is the Dirac delta function, and r is the distance between the centers of mass of the ith and jth asphaltenes (i ≠ j).

This pair correlation function g(r) describes the local density of asphaltenes and is presented in Figure for the systems with different partial charges.

Figure 2

Pair correlation function g(r) between the asphaltene centers of mass for the toluene- and heptane-based systems (a) without partial charges, (b) with partial charges calculated by the HF/6-31G* method, and (c) with partial charges calculated by the AM1-BCC method.

All g(r) dependences showed strong maxima at short distances between the asphaltene centers of mass, for all of the systems examined. This indicates the formation of asphaltene aggregates. At the same time, the number of maxima and, consequently, the microstructure of the aggregates were heavily dependent on the partial charges.

In particular, there were two maxima located at distances r ∼ 0.35 and 0.5 nm in the systems without partial charges, Figure a. These maxima may correspond to two possible distances between the adjacent asphaltene centers of masses with cofacial parallel (face-to-face) and parallel displaced (offset) stacked geometries in the aggregates, respectively; see Figure a. In the systems with partial charges calculated by the HF/6-31G*(RESP) method, the polyaromatic core carried the charge of −0.17; see the Supporting Information. This could result in a decreased overlap area between the cores.[47] Indeed, it turned out that the configuration of the asphaltenes with cofacial parallel stacked geometry was less preferable than the parallel displaced stacked geometry, since the maximum at r ∼ 0.35 nm was much lower than the maximum at r ∼ 0.5 nm; see Figure b. When the partial charges calculated by the AM1-BCC method were used, the adjacent asphaltenes adopted parallel displaced stacked geometry with r ∼ 0.6 nm; see Figures c and 3b. This result could stem from the fact that the electrical negativity of the polyaromatic core of the model asphaltenes having a charge of −0.29 was large enough to prohibit the formation of cofacial parallel stacked geometry in this case. The maxima at larger r in Figure can be attributed to the distances between the asphaltenes separated by one or several molecules in the aggregates and also to the distances between the asphaltenes in the configuration with T-shaped (edge-to-face) geometry. Thus, the presence of π–π interactions between the asphaltenes was established for all sets of the partial charges, taking into account that the center-of-mass distances between the asphaltenes were considered in the calculations.[48]

Figure 3

Aggregation of the asphaltenes with the adjacent asphaltenes adopting (a) both cofacial parallel (face-to-face) and parallel displaced (offset) stacked geometries with r ∼ 0.35 and 0.5 nm, respectively, and (b) mainly parallel displaced stacked geometry with r ∼ 0.6 nm.

We conclude that the molecular orientation of the asphaltenes in the aggregates was influenced by the partial charges. The configuration of the adjacent asphaltenes with parallel displaced stacked geometry dominated in the asphaltene aggregates, in agreement with both previous experiments and quantum-chemical calculations.[13,49] This was not the case in the systems without partial charges, where cofacial parallel stacked geometry was found to be the most likely outcome.

It should also be noted that the values of g(r) in the case of heptane solutions were always larger than the corresponding values in the case of toluene solutions, even at large values of r. It can be assumed that the strength of the asphaltene interactions and the size of the aggregates were larger in heptane than in toluene. The highest difference between the g(r) curves was observed for the partial charges calculated by the AM1-BCC method, whereas when the HF/6-31G*(RESP) method was used, the dependences almost coincide.

Overall, the analyses of the mutual arrangement and the interactions of asphaltenes complemented each other and revealed the necessity of accounting for the partial charges. Among the partial charges considered, the set calculated by the AM1-BCC method more accurately reproduced the relaxation, molecular orientation, and aggregation behavior of the asphaltenes in toluene and heptane. Therefore, this set was chosen for the further analysis of the asphaltene aggregates presented below.

Sizes of Asphaltene Aggregates

The average number of asphaltene molecules in the aggregates, in other words, the aggregation number, is one of the main characteristics used in both experimental and theoretical studies to evaluate aggregation behavior of asphaltenes in various surroundings.[14,15,25,26,50]

In the present simulations, the aggregation number is estimated by a neighbor search for asphaltenes using a geometric criterion. The simplest geometric criterion was used. Any two individual asphaltenes are neighbors if the nearest distance between their atoms does not exceed a certain value, the so-called cutoff radius rcut.[26] Here, rcut was set at 0.4 nm. Determining the groups of neighboring asphaltenes allows one to describe their aggregates and, consequently, the average number of molecules in the aggregates. Note that the individual asphaltenes, i.e., those without neighbors, are not considered as aggregates in this analysis.[26] The time dependence of the number-average aggregate sizes of asphaltenes ⟨N⟩(t) in toluene and heptane is presented in Figure .

Figure 4

Time dependence of the number-average aggregate sizes ⟨N⟩(t) of the asphaltenes in (a) toluene and (b) heptane in the systems with partial charges calculated by the AM1-BCC method. The horizontal dotted line shows the total number N0 = 50 of the asphaltenes in the systems.

It can be seen in Figure a that the time dependence of the number-average aggregate sizes of asphaltenes in toluene fluctuated around the average value of 7 ± 3 molecules. This estimate is in line with a number of experimental studies.[15,17,24] In heptane, large clusters of several aggregates were formed, which comprised almost all of the asphaltenes in the system and were stable on time scales on the order of hundreds of nanoseconds (Figure b). Large fluctuations were observed, and the clusters could easily be destroyed by thermal fluctuations.

Important information on the behavior of asphaltenes in toluene and heptane is contained in the distribution of the probability density of the aggregate sizes P(N); see Figure .

Figure 5

Distribution of the probability density of the aggregate sizes P(N) for asphaltenes in toluene and heptane.

In the case of toluene, the P(N) distribution decayed rapidly to N ∼ 35 molecules, with dimers constituting the dominant state of aggregation. This means that the asphaltenes are dissolved and uniformly distributed in toluene.[26] In turn, the distribution of the aggregate sizes was bimodal in the case of heptane, suggesting the presence of two main populations of aggregates in the system: dimers and large clusters of about 40–45 asphaltenes. These results confirm the strong tendency of asphaltenes to aggregate in heptane.[41]

Therefore, analysis of the evolution of the aggregate sizes and distribution confirms that the model asphaltenes matched their definition based on solubility: in toluene, they dissolved, forming small aggregates, while in heptane, they aggregated.

Fractal Dimension of the Aggregates

We now turn to a comparative analysis of the structure of aggregates composed of a similar number of asphaltenes in toluene and heptane. One of the characteristics widely used in experiments is the dependence of asphaltene mass M in the aggregate on its radius of gyration Rg.[14,19,22,51] The M(Rg) dependence allows one to study the influence of the solvent type on the sizes of the aggregates with similar composition, as well as to define the fractal dimension Df of the aggregates from the scaling dependence M ∼ Rg; see Figure .

Figure 6

Mass M of asphaltene aggregates versus their radius of gyration Rg in toluene and heptane. The arrows indicate the scaling dependence M ∼ Rg, where Df is the fractal dimension of the aggregates calculated as the slope of the M(Rg) curve.

Figure shows that the Rg values of the aggregates with similar masses in heptane were lower than the ones in toluene. Thus, the sizes of aggregates were influenced by the solvent type. Given a similar number of asphaltenes, the aggregates were more compact in heptane than in toluene.

It can also be seen that the fractal dimension of the aggregates was equal to 2 in the case of toluene solution, M ∼ Rg2. The result obtained is in an excellent agreement with the known experimental studies for asphaltene aggregates in toluene, where Df = 1.7–2.1 was reported.[14,19,22] Physically, this fractal dimension indicates that the structure of the asphaltene aggregates in toluene is spherical with high porosity or planar.[51] The M(Rg) dependence in the case of heptane had a more complex shape as compared to that for toluene. Two regions with different slopes, at small and at large Rg values, were observed. As in the case of toluene, the fractal dimension of the aggregates was equal to 2 at small Rg below 2 nm, with the dependence showing a kink at Rg ≈ 2 nm. The number N of asphaltenes forming the aggregate of such a size was about 15, while M ≈ 104 g/mol. At large Rg above 2 nm, the fractal dimension approached 3. A similar Df value was experimentally obtained for the asphaltenes in heptane.[51] According to ref (14), this result means that the model asphaltenes are insoluble in heptane. Moreover, the data obtained for the large aggregates suggested that the aggregates in heptane are rather compact and almost spherical in shape.

Therefore, the simulations reproduced not only the experimental trends of asphaltene aggregation in toluene and heptane but also the fractal dimension of the aggregates.

Solvent Entrainment inside the Aggregates

To provide a comparison with experiments, the volume fraction of the solvent inside the asphaltene aggregates was estimated. Schematically, the calculations are shown in Figure . First, an aggregate to analyze was selected, followed by building its “convex shell”, searching for the solvent molecules inside the “convex shell”, and, finally, evaluation of the volume fraction of solvent inside the aggregate. “Convex shell” can be viewed as a set of triangles that are built to surround all of the asphaltene atoms within the aggregate, the vertices of these triangles being the atoms of the aggregate surface. Qhull library[52] was used to build the “convex shells” of the asphaltene aggregates based on the atom coordinate files in the systems simulated.

Figure 7

Scheme to calculate the solvent volume fraction inside the asphaltene aggregates. (a) Instant configuration of the selected asphaltene aggregate. The example of the building of the “convex shell” (b) without and (c) with the solvent shown inside the aggregate. The aggregate is colored in blue, the “convex shell” is red, and the solvent is green. (d) Instant configuration of the selected asphaltene aggregate (blue) together with solvent molecules (red) inside the “convex shell”.

The dependence of the solvent volume fraction fv inside the aggregates on the number N of asphaltenes in their composition, both for toluene and heptane solutions, is shown in Figure .

Figure 8

Number N of asphaltenes versus the volume fraction fv of toluene and heptane inside the aggregates. The area colored in gray indicates the experimental range of toluene entrainment inside the aggregates reported by Gawrys et al.[18]

It was established that the volume fraction of the solvent inside the aggregates with the same number of asphaltenes in heptane is lower than in toluene; see Figure . The data obtained are in line with the analysis of the aggregate sizes discussed above. For the aggregates consisting of 20–35 molecules (Rg ∼ 2.6–3.3 nm) in toluene, a comparison with experimental data can be made. To the best of our knowledge, the reported fraction of toluene in large (2.5–8 nm) asphaltene aggregates lies in the range from 30 to 50 vol % in the samples at 1 wt % asphaltenes.[18] These data are in a good agreement with our results obtained at 10 wt % asphaltenes. At the same time, the fraction of heptane in the asphaltene aggregates turned out to be readily observable.

Thus, the asphaltene aggregates contained a significant amount of the solvent. This result may stem from the fact that the aggregates and their surfaces have a rather complex shape and thereby they are filled by the solvent molecules.

Conclusions

The influence of partial charge parametrization on the aggregation behavior of asphaltenes was studied using microsecond-long atomistically detailed classical molecular dynamics simulations. To this end, the asphaltenes in toluene and heptane solutions were considered according to their solubility-based definition, using the reliable and widely utilized GAFF force field together with different atomic partial charges. The standard sets of GAFF partial charges were calculated by the AM1-BCC and HF/6-31G*(RESP) methods and implemented in the simulated models. It was found that the two sets of partial charges result in different electrical negativities of the polyaromatic core of the model asphaltenes. The AM1-BCC method led to a more negative polyaromatic core than the HF/6-31G*(RESP) method. Systems without partial charges were simulated as well.

The mutual arrangement of the asphaltenes and the interactions between them were analyzed, and a significant influence of the partial charge parametrization on the relaxation, the molecular orientation, and the aggregation behavior of the asphaltenes was established. In the case of the partial charges calculated by the AM1-BCC method, it was shown that the relaxation of asphaltene aggregates is much slower in heptane than in toluene. In turn, the adjacent asphaltenes in the aggregates mainly adopted a parallel displaced stacked geometry. The strength of the asphaltene interactions and the size of the aggregates were greater in heptane than in toluene. By contrast, the other set of partial charges failed to reproduce these characteristics in toluene and heptane, possibly due to the insufficient electrical negativity of the polyaromatic core of the model asphaltenes. Thus, we found that the set of partial charges obtained by the AM1-BCC method is more promising, and this set was selected for all subsequent simulations.

The simulated systems were checked against experimental data on asphaltene aggregation behavior in toluene and heptane, the fractal dimension, and the solvent entrainment of the asphaltene aggregates. In particular, the asphaltenes in toluene mainly formed small aggregates composed of only a few molecules, while a strong tendency toward aggregation was observed for the asphaltenes in heptane. Thus, the validity of the solubility-based definition of the asphaltenes was confirmed. The sizes of the aggregates with the same number of asphaltenes were lower in heptane than in toluene. Moreover, the fractal dimension of the large aggregates Df = 2 suggested the formation of planar or spherical porous structures in toluene. The large aggregates having the radius of gyration greater than 2 nm formed rather compact spherical structures in heptane, since their fractal dimension was close to Df = 3. It was also established that the volume fraction of the solvent inside the aggregates consisting of a similar number of asphaltene molecules in heptane is lower than in toluene. Overall, a good agreement was found between the simulation results and the experimental data.

To summarize, we considered a number of experimentally valuable characteristics of asphaltenes in toluene and heptane. We conclude that the choice of partial charges may have a significant impact on the electrical negativity of the asphaltene core and the aggregation behavior of the asphaltenes in toluene and heptane. The GAFF force field together with the partial charges calculated by the AM1-BCC method could be successfully used to study the nanoscale structure of asphaltenes in toluene and heptane.

Model and Computational Methods

Asphaltenes naturally occur in the environment as a mixture of molecules, and thousands of different chemical structures of asphaltenes can potentially exist.[3] This is the reason why solubility in toluene and heptane, rather than a specific chemical structure, is the basis for their definition. An important step in the simulations of asphaltenes is the choice of their model mixture and chemical structure. The choice of the mixture of asphaltenes can never be universal, given that only a limited number of different molecules can be examined in the simulations.[26] Moreover, taking different molecules into account can significantly increase the computational resources required for the simulations.[53] Nevertheless, the definition of asphaltenes assumes that a monodispersed sample can represent the solubility of the asphaltenes in toluene and heptane.[26] Such an assumption is common in simulational research.[23] Therefore, a monodispersed sample of asphaltenes was considered in this study.

In terms of their structure, asphaltenes are complex polycyclic aromatic hydrocarbons. However, the most common characteristics of the behavior of asphaltenes in toluene and heptane can be investigated using a model molecule with a “typical” chemical structure.[26] Therefore, following the studies of Mullins[1] and Li and Greenfield,[54] we considered an asphaltene molecule with the “typical” chemical structure shown in Figure , which has been successfully used in the past to describe the structure of asphalt[55−57] and the water–oil interface.[58] This asphaltene model has a molecular weight of 707.1 g/mol. It contains seven aromatic rings, several long and short linear aliphatic side groups, and one sulfur heteroatom in the polyaromatic core of the molecule. It should be pointed out that none of the previous studies checked whether the model molecule can be considered as an asphaltene based on its solubility as the main definition; this check is performed in the simulations reported in this paper.

Figure 9

Chemical structure of the model asphaltene molecule chosen for this study; see refs (1, 54−58).

The interactions in the simulated systems studied were described using the GAFF force field.[59] This force field has been found to be reliable in simulations of both asphaltenes[39,40,60] and polyaromatic compounds.[61] The simulation input files, containing information on the atom parameters, as well as the bonded and Lennard-Jones interactions between the atoms, were generated using the ACPYPE tool, as recommended by the GAFF developers.[59,62] It is also recommended to use atomic partial charges calculated by the AM1-BCC or HF/6-31G*(RESP) methods to account for the Coulomb interactions in GAFF.[32,59] Therefore, in the present simulations, both standard sets of atomic partial charges were utilized to analyze the influence of their parametrization on the asphaltene aggregation behavior. The ACPYPE tool was used to perform calculations by the AM1-BCC method. Another set of partial charges was calculated by the HF/6-31G*(RESP) method with the aid of the Gaussian 09 software.[63] The calculated sets of partial charges can be found in Table S1 reported in the Supporting Information. It should be mentioned that the two sets of partial charges are quite different: when calculating the total charge of the polyaromatic core of the model asphaltene molecule, we found that the AM1-BCC and HF/6-31G*(RESP) methods gave charges of −0.29 and −0.17, respectively. Thus, the former method resulted in an electrically more negative polyaromatic core in comparison with the latter. Therefore, it is reasonable to analyze which set of partial charges would allow us to reproduce the aggregation behavior of asphaltenes. Additionally, systems without partial charges were also considered for comparison, i.e., zero partial charges were used.

The main goal of the present simulations is twofold: (i) to choose the most promising atomic partial charges on the basis of the solubility of the model asphaltenes in toluene and heptane by studying their aggregation behavior; (ii) to address the experimentally valuable characteristics of the asphaltenes in toluene and heptane, such as their sizes, fractal dimensions, and solvent entrainment in the asphaltene aggregates. Therefore, two types of solutions were considered, and a total of six systems were studied, which differed in terms of partial charges and solvent type. Analysis of the simulation approaches of asphaltenes in toluene and heptane using atomistic molecular dynamics allowed us to adapt the existing protocols; see the review in ref (23). At the first stage, a simulation cell containing the number N0 = 50 of asphaltenes was generated. For this purpose, the asphaltenes and the solvent molecules were randomly added into the periodic cubic cell. An example of the configuration obtained for the heptane-based system filled with asphaltenes is presented in Figure . The asphaltenes were simulated at the concentration of 10 wt % corresponding to their usual concentration in oil.[3,64] Accordingly, the system contained 3453 toluene or 3175 heptane molecules. To reach experimental densities, the obtained cells were compressed for 5 ns at a temperature of 300 K and at a pressure of 150 bar.

Figure 10

Snapshot of the instant configuration of the system composed of asphaltenes after the addition of the heptane molecules. For visual clarity, lines have been used to draw heptane molecules, while the atoms inside the asphaltene molecules are depicted by spheres. The blue lines indicate the simulation cell.

After the compression stage, systems with various sets of atomic partial charges were simulated at a temperature of 300 K and a pressure of 1 bar. Time scales on the order of 100 ns are the state of the art for atomistic simulations.[23] However, the displacement of asphaltenes to a distance comparable with their size can take hundreds of nanoseconds,[65] suggesting that a lower estimate for the equilibration times is required. It is therefore important to describe the aggregation behavior of asphaltenes on the microsecond time scale. Thus, the simulation was run for 5 μs for each system. The total time of the simulations performed in this study was 30 μs.

The Gromacs 5.1.4 package was used to perform the molecular dynamics simulations.[66] The equations of motion were integrated with the time step of 2 fs. Temperature and pressure were coupled using the Nosé–Hoover and Parrinello–Rahman algorithms, with time constants of 0.4 and 1 ps, respectively.[67−69] The nonbonded interactions were truncated at 1 nm, similar to that in ref (70). The particle mesh Ewald algorithm was used to deal with the electrostatic interactions in the systems with partial charges.[71] The simulations were analyzed using Python scripts powered by the MDAnalysis library,[72,73] while the VMD software[74] was utilized for the visualization of the trajectories.