Molecular Modelling and Simulation of Asphaltenes

Prévia do material em texto

Department of Chemical & Process Engineering
M.Eng Chemical & Process Engineering
18530
Molecular Modelling and Simulation of Asphaltenes
This project is submitted in partial fulfilment of the regulations governing the award of Degree of
MEng in Chemical Engineering at the University of Strathclyde
Author: João Luís Lemos Fontes Silva Costa
Project Supervisors: Jason Zhang & Paul Mulheran
April 27, 2015
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Summary
Products derived from crude oil, such as gasoline, kerosene and many other stocks, are currently vital
to developed countries due to their importance in the global economy, with especial focus on the
transportation sector. Along the years, petroleum industry have faced several difficulties in all its
stages, from the mapping of reservoirs to the crude oil refining procedures. The heaviest fractions of
petroleum, classified as asphaltenes, have been related to many undesirable effects in the extraction of
crude oil and also during the transportation steps, compromising the pipelines and wellbores, owing
to their tendency to form aggregates. Understanding asphaltenes is a vitally important matter for
scientists, engineers and operators in the petroleum industry, specially in the areas where enhanced oil
recovery strategies are required. Asphaltenes are complex and yet not determined molecular structures,
formed usually by a group of polyaromatic rings surrounded by varied alkane chains. Molecular
Dynamic (MD) simulations are one of the main methods used these days to link macroscopic properties
of asphaltenes with their molecular characteristics. In this project, four asphaltene (two big core and
two small core) models were built to study how MD simulations could be performed for asphaltenes in
toluene and heptane and whether the structural differences would lead to different aggregation profiles.
It was observed that the small core models did not aggregate in any case, regardless of the solvents
they were in. However the big core models indeed presented tendencies to form stacks, specially when
in heptane. Van der Waals forces were confirmed as the driving forces of the interactions between
asphaltenes and higher number of aromatic cores showed to lead to stronger aggregation. The presence
of aliphatic tails changed the shape of the aggregates favouring parallel conformation rather than
T-shaped ones.
April 27, 2015 i
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Table of contents
1 Introduction 1
2 Experimental Section 3
2.1 GROMACS simulation package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1.1 Potential Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1.2 Force Field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1.3 Periodic Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1.4 Radial Distribution Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Molecular Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Initial Simulation Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.4 Details of Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.5 Analysis tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3 Results 8
3.1 Small core models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.1.1 Small Cores in Toluene . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.1.2 Small Cores in Heptane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2 Big core models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2.1 Big Cores in Toluene . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2.2 Big Cores in Heptane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3 Aggregation events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.4 Structural analysis of models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4 Reflection & Review 19
5 Conclusions 20
6 Nomenclature 21
7 References 22
Appendix A Samples of MD files 25
A.1 Protein Database file (.pdb) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
A.2 Molecular Dynamics Parameters file (.mdp) . . . . . . . . . . . . . . . . . . . . . . . . 26
Appendix B Insertion of charges in topologies 27
Appendix C Extra simulations 28
C.1 Water-toluene box . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
C.2 Toluene-heptane mixture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
April 27, 2015 ii
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Acknowledgements
First of all, I would like to thank my supervisors Dr. Paul Mulheran and Dr. Zhenyu J. Zang for their
full commitment and support throughout the project. They showed me how to think critically about all
the aspects involved in this work and helped me to make the most of it. Secondly, I am very thankful
to the Brazilian federal government, with special reference to the National Council for Scientific and
Technological Development (CNPq) for the financial support that made my participation in this
project possible. Additionally, I would like to recognize the steadfast support of my dear friends André
Crescenzo and Vanesa Peixoto for sharing their thoughts and experiences about their own projects with
me, which helped me to overcome many challenges I have faced. PhD candidates Dorin Simionesie,
Javier Corona Amengual and Alessia Centi also deserve my sincere recognition for their advice regarding
project management skills and simulations, which enabled me to fulfil my goals, specially during my
first experiences with the simulation tools adopted for the work. Last but not least, I must thank all
my family and friends for their entirely incentive and encouragement during the project and at all times.
April 27, 2015 iii
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
1 Introduction
Technologies used in petroleum industry are currently one of the keystones of developed countries
and a mandatory characteristic of developing economies due to the undeniable importance of products
derived from crude oil, e.g. gasoline, kerosene, and Liquid Petroleum Gas (LPG), specially with
regards to the transportation sector (Mullins, Sheu, Hammani, & Marshall, 2007). Nevertheless, many
challenges arise from this demand: starting with the mapping and availability of potential sources of
oil, plus the extraction and refinement of crude oil in its different fractions at the refineries.
Petroleum industry typically split crude oil into four different fractions based on their chemical
properties: saturates, aromatics, asphaltenes and resins (SARA) (Burya, Yudin, Dechabo, Kosov, &
Anisimov, 2001). Asphaltenes are commonly defined as a solubility class of complex organic molecules
soluble in aromatic solvents, such as toluene, whilst insoluble in n-alkanes, just as n-heptane and
pentane (Wang et al., 2012). It is recognized that these compounds have a molecular weight between
500 and 750 Da, with aromatic cores of 4 to 10 aromatic rings, and aliphatic tails ranging from 3 to 7
carbons. Asphaltenes’ aromaticity, the ratio of the number of aromatic carbons to the total number
of carbons in the molecule, typically ranges between 0.35 to 0.65 (Sedghi, Goual, Welch, & J., 2013).
Asphaltenes may also present H/C atomic ratios from 1.0 to 1.2, heteroatoms such as N, O, S, and
traces of metals like Nickel and Vanadium (Spiecker, Gawrys,& Kilpatrick, 2003; Sheremata, Gray,
Dettman, & McCaffrey, 2004).
Two different types of asphaltenes’ structures may be found in crude oil. They are named as
Continental and Archipelago models (Mullins et al., 2007). The differences between these two
configurations are basically the way that the poly-aromatic cores and aliphatic tails are organized.
Continental structures (also known as ‘Island‘ structures) present only one condensed poly-aromatic
(PA) core, surrounded by aliphatic tails, while the archipelago arrangement has several smaller cores
connected by its tails (Tukhvatullina et al., 2013). These two types are exemplified in figure 1. Although
both structures can be found in different sources of crude oil, the continental structure can often
represent the major portion of asphaltenes’ properties that were studied more recently, indicating the
importance of using this structural concept as a reference for further analysis (Mullins et al., 2012).
Figure 1: Different Types of Asphaltene structure. a) Continental type; b) Archipelago type.
(Tukhvatullina et al., 2013)
Asphaltenes are specially present in the heavier fractions of petroleum and require a special attention,
in view of their tendency to form aggregates. It can often result in wellbore clogging and fouling of
pipelines (Subirana & Sheu, 2013). Also, coke formation in crude oil is being related to core aggregation
of asphaltenes in the pyrolysis process of petroleum thermal degradation (Garcia Barneto, Carmona,
& Barrón, 2015). Additionally, many studies have reported that asphaltenes play a significant role
stabilizing Water-in-crude oil emulsion.(McLean & Kilpatrick, 1997; Kokal, 2007).
Concerns related to asphaltenes are particularly important in the locations where Enhanced Oil
Recovery strategies (EOR) are required (Alvarez, Poteau, Argillier, Langevin, & Salager, 2009).
Nonetheless, petroleum industry as a whole intends to address the challenge of understanding as-
phaltenes. The main efforts in this case are either to enable the refinement of asphaltenes in lighter
April 27, 2015 1
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Figure 2: Scheme of asphaltenes’ aggregation
(Mullins et al., 2012)
fractions, which is currently not possible, or find alternatives to separate it from the petroleum liquids
before the refining process (Subirana & Sheu, 2013).
Self-aggregation of asphaltenes can be compared with the behaviour of surfactants in a micelle,
despite the fact that the driving forces of the former cannot be well established (Ibid). Van der Waals
interactions seem to be one of the key factors that explain the stacking of the polyaromatic units.
However, it is stated that the pi-electron binding competes with the steric repulsion caused by the
alkane chains, which can increase solubility in a considerable extent (Teklebrhan, Ge, Bhattacharjee,
Xu, & Sjöblom, 2012). Polar interactions may as well play an important role, depending on the
polarizability of the asphaltene’s structure. This effect can be evaluated when comparing the increased
solubility in very polar solvents such as acetone (Mullins et al., 2007). On the other hand, variables
such as pressure, temperature, and composition of environment around asphaltene molecules can not
be neglected, adding even more complexity to the subject (Diallo, Cagin, Faulon, & Goddard III,
2000).
Thus, the investigation of asphaltene‘s features, in particular its structure, is absolutely needed
to enable more accurate predictions of the potential for production in oil reservoirs. Through its
structure, it is expected to understand asphaltene‘s phase behaviour, mainly looking at the molecular
scale (Mullins et al., 2007). It is worth mentioning the advances resulting from experiments that
analysed bulk properties of asphaltene. Good estimations for molecular weight values of asphaltenes
can be obtained by means of nuclear magnetic resonance (N.M.R) (Ali, Al-Ghannam, & Al-Rawi,
1990). Moreover, small-angle neutron scattering (SANS) and small angle x-ray scattering (SAXS) are
two particularly efficient methods to study colloidal structures and provide information about the
effects of temperature, solvent and source of crude-oil in the aggregation phenomena (Spiecker et al.,
2003). To study the presence and quantity of heteroatoms, such as Nitrogen and Sulfur, x-ray methods
can also be applied. Nowadays, the X-ray Absorption Near Edge Structure (XANES) spectroscopy is
known for providing a powerful analysis of those particular atoms or functional groups (Subirana &
Sheu, 2013).
Even though there are a wide range of methods similar to the above mentioned that might help
scientists to learn more about asphaltenes, there is still no method that can precisely determine the
molecular structure of such a complex set of compounds (Mullins et al., 2007). From this perspective,
molecular simulations have a great potential in the petroleum industry, owing to the fact that the
computational approach can implement a straight route from microscopic properties of a model (atomic
composition, molecular geometry, intermolecular interactions, and so on) to macroscopic outcomes
that can be compared with real experiments (P. Allen & Tildesley, 1987). In other words, one can
simulate the influence of each variable involved in the system of interest (e.g. number of aromatic rings
within PA core) in similar conditions to those where it is usually found on site.
Molecular Dynamics (MD) is being used among other techniques by several research groups to
tentatively reproduce asphaltene’s experimental data (Headen, Boek, & Skipper, 2009; Kuznicki,
Masliyah, & Bhattacharjee, 2008; Sedghi et al., 2013). Molecular dynamics can be defined as the set of
simulations that study particle interactions through the solutions of the classical equations of motion
described by Newton. It is reported that this method was firstly used in the fifties and is continuously
being developed with the computational advances and a thoroughly study of molecular interactions
and particles movement, specifically with regards the Lennard-Jonnes potentials and electrostatic
interactions (P. Allen & Tildesley, 1987).
April 27, 2015 2
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
A Great part of the simulations in the asphaltene field intend to probe the structural models chosen,
aiming a reliable reproduction of physical and chemical properties of real asphaltenes. More than
that, it is important to ensure the reliability of the simulation package and the set of parameters
used as well. Despite the fact that there is broad range of programmes that can be employed to run
MD simulations, there is a collection of data required for those softwares, i.e. three-dimensional (3D)
coordinates, molecular geometry and topology parameters. These information are specific for each
molecule, then it might also be necessary to access MD-related applications, such as topology builder
scripts or even additional tools to perform Quantum Mechanic (QM) calculations.
Despite MD simulations might enable one to gain a whole new set of information with regards
molecular behaviour, there is still limitations due to the huge amount of data needed to compute
intermolecular interactions (M. Allen, 2004). It restrains the time and the dimensions of simulations,
currently in the house of nanoseconds (ns) and nanometers (nm) respectively. Thus, the utilization of
artifices, like implementing cut-offs and periodic boundary conditions, are commonly adopted. These
and other tactics will be further described in section 2.
In this project, different molecular geometries were adopted in the study of asphaltenes’ aggregation
behaviour in two organic solvents (toluene and heptane). Through these simulations, it is expected
to have a more clear picture of the driving forces in theaggregation phenomena. Furthermore, this
project intends to probe the path used to build the molecular structure of asphaltenes. The results of
these simulations will be compared with the literature available, concluding whether or not the routine
utilized was consistent.
A brief introduction about MD theory and equations will be made, focusing on the key concepts of
molecular dynamics. Then, some insights of how asphaltene models were built and the routine adopted
will also be given. Ultimately, the results of all the simulations will be discussed, followed by reflection
and conclusions of the project. Additional information regarding some of the parameters used within
simulations are available in the appendix A.
2 Experimental Section
2.1 GROMACS simulation package
The Groningen Machine for Chemical Simulations (GROMACS) is a free software available under
the GNU Lesser General Public License 2.1 and is designed to perform molecular dynamic simulations
and energy minimization, the main techniques in the sphere of molecular modelling and computational
chemistry (van der Spoel et al., 2013). This engine is widely used for simulation of biomolecular
systems, with special emphasis in proteins and lipids in an aqueous environment (Pronk et al., 2013).
It is reported that MD simulations can provide better ensembles for the analysis of dynamic systems
and also yield better statistic data in a certain period of time than any other existent computational
strategy. However, when preparing MD ensembles, one have to make sure that the initial configuration
of the system is not too far from equilibrium, fact that could lead to excessive forces being applied in
the molecules (Ibid).
As previously mentioned, MD simulations solve numerically the classical equations of motion, that
for a system of N interacting atoms are:
mi
∂2ri
∂t2
= Fi, i = 1...N. (1)
These forces are the negative derivatives of the potential function V (r1, r2, ..., rN ):
Fi = −∂V
∂ri
(2)
The equations are solved step-by-step, and the positions are written periodically in a output file.
The coordinates as a function of time represent the trajectory of the atoms, taking to account the
temperature and pressure of the system. After primary adjustments to reach an equilibrated state,
the software allows calculation of macroscopic properties averaging over the equilibrium trajectory
(van der Spoel et al., 2013). The main features are detailed below.
2.1.1 Potential Functions
There are three different categories of potential functions that may take part in the system’s
trajectories: Non-bonded, Bonded and Restraints.
April 27, 2015 3
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
The Non-bonded potentials comprise specially Lennard-Jones and Coulomb potentials. These
functions are computed in the system based in a neighbour list, that is, the computation of these
potentials depends on the presence of those in a specified radius. These neighbourlists automatically
exclude any bonded atoms inside that range. Bonded potentials contain predefined parameters for
covalent bonds, angles, and also proper and improper dihedrals. In the same way, the Restraint
potentials are all based on fixed lists and introduce information about restraints in position, angles,
orientations, dihedrals and distances. The values that do not depended on the positions in the
simulations are referred as fixed lists (Ibid).
2.1.2 Force Field
GROMACS package, as any other standard MD tools, makes use of the concept of Force fields to
define the set of parameters used to calculate the potential functions. Force fields are the group of
files that store information about the fixed lists that contain all the data GROMACS will require to
compute interactions between particles (van der Spoel et al., 2013). Therefore, the topology file .top for
each molecular specie should contain settings consistent with the force field adopted. The development
and application of each force field is directly related with the characteristics of the systems one might
want to study.
Force fields can be classified by the manner their parameters are set for the simulations:
� All atoms: Contain parameters for every single atom present in the system;
� United atoms: Do not contain parameters for non-polar hydrogens;
� Coarse-grained: Group several atoms into "super atoms" in order to create an abstract representation
of molecules.
The OPLS/AA (Optimized potential for liquid simulation / All Atoms) was the Force Field adopted
for the simulations due to its precision in the representation of aromatic compounds properties, such
as enthalpy of vaporization and density (Jorgensen, Maxwell, & Tirado-Rives, 1996). This same force
field has been already used in simulations of asphaltenes and is probably better suited than any other
to study the subtleties of their structure and how it affects aggregation (Headen et al., 2009; Boek,
Yakovlev, & Headen, 2009; Sedghi et al., 2013).
Non-bonded interactions in OPLS/AA are represented by (3):
Eab =
a∑
i
b∑
j
[
qiqje
2
rij
+ 4�ij
(
σ12ij
r12ij
− σ
6
ij
r6ij
)]
fij (3)
Where the first term in the sum refers to the coulombic interactions, at the same time the second term
corresponds to the Lennard-Jones potential and Eab means the interaction energy between molecules a
and b due these two factors. It is reported that the interaction between atoms of the same molecule
separated by three or more bonds also are represented by this equation (van der Spoel et al., 2013).
Bonded interactions are read from the internal information present in the fixed lists of the force field,
while the constraints on the models can be computed through different algorithms. Usually, the widely
accepted LINear Constraint Solver (LINCS) is adopted in most simulation packages and is compatible
with OPLS/AA (Hess, Bekker, Berendsen, & Fraaije, 1997).
2.1.3 Periodic Boundary Conditions
Reproduction of bulk properties in small systems, such as those commonly used in MD simulations
(10 < N < 10,000 molecules), can be achieved when periodical boundaries are used for the simulation
box. It can be described as a box that is thoroughly surrounded by copies of itself (P. Allen & Tildesley,
1987). So, when a molecule cross the edge of one box it will reappear in the opposite side of the
neighbour copy, as shown in figure 3.
2.1.4 Radial Distribution Functions
Radial Distribution functions (RDF’s), which are also known as pair correlation functions gab(r)
are widely used to characterize properties of systems, providing important information in terms of the
randomness of their structures (Hansen & McDonald, 2013; van der Spoel et al., 2013). As stated by
the latter, a RDF between particles of types A and B can be defined as follows:
gAB =
〈ρB(r)〉
〈ρB〉local =
1
〈ρB〉local
1
NA
NA∑
i�A
NB∑
j�B
δ(rij − r)
4pir2 (4)
April 27, 2015 4
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Figure 3: System with periodic boundaries in 2 dimensions
Where 〈ρB(r)〉 is the density of type B particles around A within a radial distance r, and 〈ρB〉local
the average density of all particles B in relation to the maximum possible distance (rmax) around
particles type A. Usually, rmax is given as half of the cell box length (Ibid).
Figure 4: RDF’s graphic description
Figure 5 illustrates regular radial distribution functions for liquids. As the radius increases, a
tendency for an uniform distribution is perceived, i.e. g(r) converging to 1. The top peaks will give
information about the average distance between the molecules. It can be seen that heptane has a very
distributed range of possible distances with low peaks, especially due to its non-polar characteristics.
Although toluene is also non-polar, the pi delocalized electrons favourssome specific orientations, which
affect how toluene molecules attract each other (Headen et al., 2010).
Therefore, utilization of pair correlation functions might provide very useful information in respect
to aggregation behaviour as a whole, showing clearly the many different ways that asphaltenes might
interact and form groups (Gao, Xu, Liu, & Yuan, 2014).
2.2 Molecular Models
For the purpose of this project, four asphaltene models were built , two with a small core (SC) of
aromatic rings and two models with bigger cores (BC). All aliphatic tails for these structures contain
five carbons. The SC models do not attend all the requirements to be compared with real asphaltenes
(Tang et al., 2015). Instead, these models were used to analyse whether similar geometries would
display the same behaviour even in lighter structures. The geometries of the BC structures were
based on several different publications (Mullins et al., 2007; Sedghi et al., 2013; Kuznicki et al., 2008),
differing only in terms of the presence of heteroatoms. The presence of such atoms increases the
complexity of charges insertion step, although it is included in the future plans of the project to study
their role within asphaltenes’ molecules. More details about the models are shown in table 1. The
first stage in the construction of the models consists of obtaining the chemical structure, then generate
April 27, 2015 5
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Figure 5: Example of RDF for toluene and heptane boxes
Figure 6: Structures of the four different asphaltene models
3-D coordinates for all the molecules used, and lastly, set up the parameters required by the force field
(angles, bonds, dihedrals, etc.). For each of these purposes, different softwares might be necessary. The
Accelrys Draw package was employed to design the models (Accelrys Draw 4.2 - Academic Version,
2007). This open source software is suited for academic use and can quickly draw chemical structures.
These structures were then exported to other applications where the 3D-optimized structures were
created. PRODRG online server (Schüttelkopf & van Aalten, 2004) and ArgusLab (Thompson, 2004)
could both be used in this step with the same objective of obtaining a protein database file .pdb, which
is a widely accepted input file format for MD packages like GROMACS. Once with the pdb files, the
MKTOP application was used to generate all the parametrization for the OPLS topologies, apart from
atomic charges (Ribeiro, Horta, & Alencastro, 2008).
The insertion of the charges was based in existent topologies used for simulation of aromatic
compounds using OPLS/AA as the force field (Caleman et al., 2012). In addition to that, a scheme
specific for an asphaltene case is given in the appendix B (Morimoto et al., 2014). According the latter,
there were used seven different C atom types for the asphaltene and two different types for H’s. These
April 27, 2015 6
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Table 1: Structural parameters for asphaltene models
Model no. aromatic carbons no. aliphatic tails Mw(g/mol) Aromaticity H/C
SC0 14 0 204 0.875 0.750
SC3 14 3 414 0.452 1.355
BC0 26 0 342 0.765 0.701
BC3 26 3 642 0.531 1.102
differentiation is required to give the molecule more real features, such as the planar conformation of
aromatic rings. From OPLS database, they can be described as follows:
� Opls_135, opls_136 and opls_137 represent CH3, CH2 and CH alkyl carbons respectively;
� Opls_145 and opls_147 represented the terminal and fused carbons of aromatic cores;
� Atom types opls_149 and opls_515 correspond to the next neighbours to aromatic carbon and the
sp3 carbon. These two only differ in the number of hydrogens they bond: the former has two H atoms
attached, while the latter only have one;
� Opls_140 and opls_146 were used as alkane and aromatic hydrogens respectively (Ibid).
Support for the use of empirical data in the design of topologies is given by the main paper related
with the OPLS/AA force field (Jorgensen et al., 1996). It is said that the charges of functional groups
are transferable between molecules. This eliminates the necessity of Quantum Mechanic calculations
to fit an electrostatic potential surface, which is a common procedure when charges are obtained
case-by-case.
Toluene and heptane were used as explicit solvents in the simulations. These are stated as the limits
in terms of the solubility spectrum of asphaltenes, hence they are adopted in most simulations and
experiments that study aggregation events involving these chemicals (Headen et al., 2009; Khoshandam
& Alamdari, 2010; Teklebrhan et al., 2012). Heptane structure was designed according to the routine
above mentioned, while topology and coordinates for toluene were obtained from a database of organic
liquids (Caleman et al., 2012).
2.3 Initial Simulation Configuration
The first step within each simulation was create a cubic box with suitable dimensions, and insert
six asphaltene molecules in random positions. Then the appropriate number of solvent molecules was
iteratively added to the box. The ensembles were first energy minimized to guarantee a good starting
conformation for the simulations. After Energy Minimization, the systems were equilibrated through
100 ps of NVT (temperature coupling) and 350 ps of NPT (pressure coupling) ensembles respectively.
The main goal of these equilibration steps is to ensure that temperature and pressure of the simulation
boxes are adjusted to the desired values (Sedghi et al., 2013). These ensembles are also necessary
to adjust the size of the box to the actual number of particles (Gao et al., 2014). Complementary
information about the simulation boxes are presented in table 2.
Table 2: Details of simulation boxes
Model NAsphaltene NToluene NHeptane final box sizea (nm)
Toluene systems
SC0 06 190 - 3.24 x 3.24 x 3.24
SC3 06 386 - 4.16 x 4.16 x 4.16
BC0 06 320 - 3.89 x 3.89 x 3.89
BC3 06 600 - 4.81 x 4.81 x 4.81
Heptane systems
SC0 06 - 175 3.52 x 3.52 x 3.52
SC3 06 - 355 4.46 x 4.46 x 4.46
BC0 06 - 295 4.19 x 4.19 x 4.19
BC3 06 - 553 5.18 x 5.18 x 5.18
a Initial size of SC0 simulation boxes was set at 5.00 x 5.00 x 5.00 nm3 to avoid
shrinking the box more than half of its initial size during equilibration steps.
Boxes of 7.00 x 7.00 x 7.00 nm3 were set for the remaining models.
April 27, 2015 7
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
2.4 Details of Molecular Dynamics
MD simulations were performed using the GROMACS 4.6.5 simulation package (van der Spoel
et al., 2013). The velocity rescaling algorithm was the thermostat employed (Bussi, Donadio, &
Parrinello, 2007) whilst Berendsen (Berendsen, Postma, van Gunsteren, DiNola, & Haak, 1984) and
Parrinello-Rahman (Parrinello & Rahman, 1981) barostats were used respectively. During the first 200
ps of NPT step, the Berendsen algorithm was used to relax the internal pressure in the box. Next,
Parrinello-Rahman was chosen in order to reproduce a better thermodynamic ensemble. Temperature
and pressure of all boxes were checked after the equilibration stage to ensure the systems were in the
desired temperature (298 K) and pressure (1 bar).
The Particle-Mesh Ewald (PME) algorithm was used to compute the long range electrostatic
interactions (Darden, York, & Pedersen, 1993). A cut-off of 1.0 nm was used for Van der Waals and
Coloumb interactions. This value is within the average values used in OPLS/AA and should work
fine for computing the forces between particles (Sedghi et al., 2013; Headen et al., 2009). These and
further conditions applied to the simulation boxes must be given in the molecular dynamics parameters
files (.mdp). Each one of the equilibration steps will require specific information to be contained inits
mdp file, in order to generate the correct trajectories (.trr or .xtc) and all energy terms saved in any
simulation (.edr). An example of .mdp file can be found in the appendices (see appendix A.2).
The validation of the parameters used was achieved by MD simulations with toluene at ambient
conditions. The methodology to perform these calculations was based on standard procedures (Caleman
et al., 2012). From this paper the following equations were used:
ρ = M〈V〉 (5)
and
∆Hvap = (Epot(g) + kBT )− Epot(l) (6)
Vaporization enthalpy and density were determined for this compound and the results are consistent
with experimental data, as shown in table 3.
Table 3: Enthalpy of Vaporization and density of toluene
Theoretical 1 MD simulation
enthalpy of vaporization (KJ/mol) 38.01 37.55
liquid density (Kg/m3) 864 866
1 (PubChem Compound Database; CID=1140 , 2015)
In order to observe the occurrence of aggregation for the structures designed, 20 ns simulations
were performed for all the models with the periodic boundary conditions, first in toluene, and later in
heptane. This simulation time might provide relevant information about the overall behaviour of the
asphaltene models, with the utilization of a convenient computational power.
2.5 Analysis tools
A visualization software was required to investigate the behaviour of the models during the
simulations. For this purpose, it was used the Visual Molecular Dynamics (VMD) software (Humphrey,
Dalke, & Schulten, 1996). Radial Distribution Functions (RDF) were generated to determine aggregation
phenomena by GROMACS built-in analysis tools, as well as the evaluations of temperature, pressure,
potential energies and enthalpies. It is important to state that macroscopic properties, such as pressure,
might present considerable fluctuation in a instant analysis (Pronk et al., 2013). However, in situations
like these the average value over the total simulation time is the correct parameter to be trusted.
3 Results
8 different simulations were carried out to study the importance of the structure in the aggregation
behaviour at ambient temperature and pressure conditions. Concentration of asphaltenes in each
simulation box was kept constant, around 7% wt, to simulate concentration of many sources of petroleum.
April 27, 2015 8
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
This concentration allowed simulations to occur in compact boxes, enhancing data processing steps,
which is therefore always desirable.
The criteria for aggregation used in this work was based on previous MD simulations for asphaltenes
(Sedghi et al., 2013; Gao et al., 2014). It was reported that parallel aggregation of two asphaltene
molecules usually occur when the distance between the centres of mass (COM) is around 0.5 nm for
asphaltenes with aliphatic tails, and amid 0.3 - 0.4 nm for aromatic planes. Furthermore, T shaped
aggregates could appear up to 1.3 nm. Between 0.9 and 1.3 nm, the analysis of following or previous
frames is the only way to accept or reject the aggregation. Moreover, the intensity of the peaks showed
in the radial distribution functions can also be used to determine the occurrence of aggregation. When
the heights of those peaks are above 10, indicating a local density higher than 10 times the average for
the total simulation environment, it is likely that aggregation events can be visualized (Headen et al.,
2009).
3.1 Small core models
3.1.1 Small Cores in Toluene
Figures 7 and 8 present the first and last frames of the SC0 and SC3 models respectively. The
radial distribution functions for these simulations were prepared splitting the total simulation time in
4 different sections, to identify any subtleties that could take place during the simulation. The total
average distance between the centres of mass is displayed in the RDFs as well. It can be seen that the
both SC models did not show any tendency to form nanoaggregates. The small number of PA units
did not provide enough intermolecular interactions between the models. The frames obtained for the
simulations also demonstrate that the presence of aliphatic tails in SC3 could not increase attraction
between asphaltenes to any greater extent.
Figure 7: Models SC0 in toluene: a) initial and b) final frames
Figure 8: Models SC3 in toluene: a) initial and b) final frames
The RDF outputs also show that the average distances between asphaltenes was between 0.5 and
April 27, 2015 9
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
0.9 nm during the whole simulation time. The width and height of these peaks in the RDFs are results
of the mixed forms of attraction between asphaltenes, driven by weak Van der Waals forces in SC0
and SC3 models. The latter also showing some interactions between the tails and the PA cores. The
SC0 structures remained spread around the box throughout the simulation, showing the high affinity
of this particular model with toluene.
Figure 9: RDF of SC0 models in toluene. The wide peak between 0.1 and 0.9 shows that SC0 asphaltenes are
perfectly soluble in toluene and, hence present no propensity to aggregation.
The RDF of SC3 simulation highlighted some proximity between asphaltenes amid [5-10] ns.
However the low intensity of the peak suggests that the aliphatic tails might have interacted with
the PA cores, as shown in the top left corner of figure 8.b). These interactions are not understood as
effective aggregations, but they can affect the solubility of these compounds when compared with SC0
models.
Figure 10: RDF of SC3 models in toluene. Solubility of the small aromatic cores was altered due to the addition
of alkane chains. Still, overall behaviour of those models suggested no trend in the direction of aggregation.
April 27, 2015 10
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
3.1.2 Small Cores in Heptane
Following the same procedure used for inserting asphaltene models in toluene, the SC0 and SC3
asphaltenes were solvated with heptane. First and last simulation frames are shown in figures11 and
12. In comparison with the simulations performed using toluene, an overall closeness of the molecules
in heptane can be perceived. However, it is also straightforward from the visual information that such
a reduced aromatic core would not lead to aggregation. The change in the solvent slightly increased
the probability of finding asphaltenes closer to one another, although the peaks exposed in this case
have maintained the width of simulations in toluene.
Figure 11: Model SC0 in heptane: a) initial and b) final frames
Figure 12: Model SC3 in heptane: a) initial and b) final frames
The radial distribution functions for the small core models are presented below in figures 13 and
14. One can observe that the general comportment of the RDFs were analogous during all the 20 ns of
simulation for both cases, presenting reasonable small peaks, which endorse the absence of aggregation
tendencies. SC0 models had a range between 0.6 to 0.8 nanometers, while SC3 models, exhibited a
wider range, 0.5 to 1.0 nm, owing to the interactions of alkane chains with the aromatic units. In
another words, the information obtained through the behaviour of the distribution functions match the
analysis of simulation frames and provide no evidence of stacking for the small core models in heptane.
April 27, 2015 11
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Figure 13: RDF of SC0 models in heptane. Overall arrangement of asphaltene models was the same for the
whole simulation. The low intensity of the distribution peak reached its maximum value during the last 5 ns.
Figure 14: RDF of SC3 models in heptane. Wide peak ranging from 0.5 to 0.9 nm showing different possible
orientations for SC3structures. Presence of aliphatic tails increase interactivity between asphaltene models,
whereas it does not imply the occurrence of aggregation.
3.2 Big core models
3.2.1 Big Cores in Toluene
Initial and final frames for the boxes containing the big core models in toluene are displayed in
figures 15 and 16. Radial distribution functions for different sections of simulations are presented for
BC0 (figure 17) and BC3 (figure 18) asphaltenes respectively. Big core models showed a substantial
higher tendency to stack in comparison to the small core models. During the simulations there were
observed stronger interactions between asphaltenes. However, BC0 models did not present any group
with more than pairs of molecules and only for some picoseconds, which cannot be classified as a
nanoaggregation event. For BC3 structures, groups with 3 or more molecules were formed, nonetheless
they did not remain stable for more than 2 nanoseconds.
April 27, 2015 12
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Figure 15: Model BC0 in toluene: a) initial and b) final frames
It could be perceived the stacking behaviour may change a lot depending on the number of aliphatic
tails in the molecule. The presence of three alkane chains led to different aggregation conformations,
whereas more angular stackings where observed rather parallel or T-shaped ones. One considerable
difference amid the SC and BC models is the that the same number of aliphatic tails, with the same
size, can produce much stronger interactions as the aromatic ring increases. It can be explained by the
larger distances between the tails, which can smooth the steric repulsion effects.
Figure 16: Model BC3 in toluene: a) initial and b) final frames
April 27, 2015 13
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Figure 17: RDF for BC0 models in toluene. [15-20 ns] function displays an important peak, between 0.6 and
0.8 nm, that remained almost unnoticed in the [0-20 ns] curve. This interval correspond to the distances where
T-shaped and offset parallel conformations take place.
Figure 18: RDF for BC3 models in toluene. One peak can be visualized showing some tendency to parallel
aggregation, throughout the 20 ns of simulation. A smaller peak around 0.8 highlights the presence of some
T-stackings.
3.2.2 Big Cores in Heptane
First and last frames of big core models in heptane can be seen in figures 19 and 20. A mutual
attraction behaviour was strongly pronounced in both big core models, whereas the structural differences
lead to significant disparities between the clusters. when simulated in heptane, asphaltenes got closer
to one another much more quickly, confirming the conclusions of previous works (Sedghi et al., 2013).
Through the RDFs, one can notice the presence of high and narrowed peaks near 0.5 nm, since the
earlier moments of the simulation for BC3 asphaltenes, while BC0 structures converge to a similar
behaviour in a slower pace.
The lack of aliphatic tails gives plenty space to asphaltenes to interact and aggregate in varied
ways. For BC0 models, both parallel and T-stacking aggregates were observed in our simulation. The
RDF for BC0 shows two equal high peaks near 0.6 and 0.9 nm, which means that the occurrence of
those types of aggregation was proportionally very similar. However, in the BC3 simulations, only the
parallel stackings were observed in the end, as shown in figure 20.
April 27, 2015 14
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Figure 19: Model BC0 in heptane: a) initial and b) final frames
Figure 20: Model BC3 in heptane: a) initial and b) final frames
Some interesting differences can be observed in the RDFs for each structure. For example, figure 21
shows different average behaviour during the simulation for BC0 models. The first 10 ns produce two
peaks, with the highest value for the distance between centres of mass amid 0.6 and 0.8 nm. As the
simulation progresses parallel stackings could be observed in the [10-15 ns] and [15-20 ns] curves with
even higher intensity.
April 27, 2015 15
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Figure 21: RDF for BC0 models in heptane. Parallel and T-stacks are very pronounced in the peaks around
0.4-0.5 and 0.7-0.8 nm respectively. The final curves show increase in the parallel conformations in the last 10
ns.
For BC3 asphaltenes, It can be seen that some signs of T-stacking aggregation (see figure 22) were
spotted in the [0-5] and [5-10] functions, represented by the peaks of these curves around 0.75 nm.
On the other hand, the subsequent sets show an intense presence of parallel interactions with high
peaks between 0.4 and 0.6 nm, indicating presence of slightly offsets regarding the perfect parallel
aggregation.
Figure 22: RDF for BC3 models in heptane. Intense evidences for parallel aggregation are shown in the peaks
around 0.5 nm. Smaller peaks with tops in 0.8 and 1.25 nm are related with the brief T-stacking events.
Steric effects are responsible for these offsets, since the molecules try to diminish repulsion forces
caused by the aliphatic tails. In this case, the positions where the PA cores are facing one another and
the tails are not superposed were preferred. Figure 23 shows in detail the shape of those aggregates.
April 27, 2015 16
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Figure 23: Two different aggregation forms for the same pair of asphaltenes in heptane. Figure a) shows
T-stack aggregation at 3 ns for these BC3 models, which was gradually changing for a b) parallel interaction.
The BC3 pair remained in this latter configuration until the end of the simulation, featuring the repulsion
between the tails as a more stable conformation.
3.3 Aggregation events
From all the simulations performed, there were observed effective tendencies to aggregation only in
the big core models. The shape and intensity of these events showed to be strongly dependent on the
structure of the asphaltene and the type of solvent used.
Figure 24: Aggregation compared through RDFs for big core models . Intensification in the stacking behaviour
using heptane as solvent was mainly observed for BC3 models, although BC0 structures have also followed the
same trend.
It can be noticed that the aggregation increased significantly as the properties of the structures
came nearer to those that represent real asphaltenes. Plus, the number of aromatic cores proved
to be decisive with regards interactions between the models. The presence of alkane chains is as
well established as a key factor in the shape and size of the clusters. Regarding the influence of the
solvents, it was demonstrated that aggregation behaviour is enhanced in heptane, confirming laboratory
experiments and results from other simulations (Kuznicki et al., 2008; Boek et al., 2009). The peaks
for BC0 models have doubled their height in the region related to the parallel staking and increased
almost 3 times the intensity in the distance that represents T-stacking. However, the biggest contrast
April 27, 2015 17
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
is observed between BC3 asphaltenes, which showed peaks 6 times higher in heptane, than those
exhibited in toluene, featuring a huge tendency to form parallel clusters.
3.4 Structural analysis of models
The observance of aggregation for the big core models might indicate that these molecules contain
the main properties that lead to aggregation perceived in real asphaltenes. In addition to that, the
similarities obtained in the radial distribution functions of those structures, when compared with more
complex models found in the literature, endorse the applicabilityof BC architecture for future studies.
Thus, the presence of N, O, S, Ni and Va atoms, and varied alkane chains, in a structure based on
the BC models adopted in this project might introduce worthwhile information to the investigation of
asphaltenes’ behaviour.
April 27, 2015 18
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
4 Reflection & Review
The awareness of the importance of molecular modelling for chemical engineering is probably
one of the main outcomes of this work. It was possible to see a parallel between real experiments
and simulation-based publications and observe that the present computational power can accomplish
important discovers, not only for petroleum industry, but in many other areas where chemical engineers
might be needed. For this reason it is worthwhile to seek for experience in this area and the project
was considerably handful for this purpose.
The first steps towards modelling asphaltenes were very challenging due to the lack of an introductory
literature regarding the construction of molecular models. However, one should comprehend that
very similar models may be constructed following diverse routines. It is important though, to first
understand the theory behind the softwares one want to use, to then grasp information about the
development of new models. Then it is possible to advance in the direction of more complex systems.
This project was intended to show more clearly one of those possible routines. As a result of our
approach, it is expected that this project can help students and researchers in their first contact with
Molecular Dynamics’ tools.
Computational simulations can provide an enormous range of possibilities, enabling investigation
of as many different processes as one might desire. Owing to this fact the guidance of the project
supervisors is fundamental to avoid losing the initial focus of the project. During this project many
simulation environments were constructed, the first mainly to gain experience with the MD package
and grasp what could be extracted from it (see appendix C). Later on the main simulations were
made in order to obtain the pursued results. It is worth mentioning that read through the manuals of
the simulation packages and tutorials is highly recommended and can save hours of incomprehensible
errors and alerts.
Additionally, through this research project it was possible to develop personal communication
abilities and academic skills in a massive extent. Every weekly meeting was an opportunity to deliver
ideas and points of view about the project. The preparation of short-term results for these meetings
was a great exercise to learn how to summarize key information and think critically about the project’s
targets. Furthermore, the opportunity to exchange information with other students and scientists with
regards their experiences within molecular simulations was very valuable. These moments did also
influence the manner how one may analyse his own work, learning how others have approached their
difficulties and comprehend the importance of the contribution of different individuals in the progress
of research activity.
April 27, 2015 19
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
5 Conclusions
The main goal of the project was to understand the basic concepts of Molecular Dynamics and
build a range of models that will ultimately be able to predict properties of real systems containing
asphaltenes. The routine followed to build the model structures can be seen as a simpler and reliable
alternative when compared with systems that require advanced quantum mechanic calculations. At the
same time it was a consistent approach that permitted observation of similarities between laboratory
experiments and computational models.
It was perceived that the small core structures SC0 and SC3 did not reproduce any of the tendencies
predicted for asphaltenes, which was expected given fundamental differences between general asphaltenes
and these models. These lighter molecular structures moved around the boxes behaving like liquids
especially because the deficit of characteristics that are described as the main agents of the aggregation
condition, i.e. Van der Waals forces. The absence of aggregation for these compounds express that
the routine used to build the models was not biased in conditions that would favour aggregation to
happen regardless of the models.
In another hand, the big core models showed a comportment close to what is observed in real
asphaltene samples and in a series of MD simulations. They showed the impact of the molecules’
structure to the final outcomes. Both BC0 and BC3 models exhibited higher tendencies to attract
each other in heptane than in toluene, confirming again the influence of Van der Waals forces over
the electrostatic interactions driving the aggregation tendency, and in case of simulations in toluene,
intensifying the solubility of asphaltenes.
In general terms, the results obtained in this project will serve as a start point for further analysis of
asphaltenes. Now that this ground work is established,it is possible to use MD simulations to produce
more refined information about the aggregation phenomena, such as the calculation of the Potential
of Mean Force and the dimerization energies, which are important outputs that are currently being
studied. The investigation of bigger and more complex asphaltenes is included in the future plans of the
research group to seek a thorough understanding of each variable in the overall asphaltene-asphaltene
interaction and as well of asphaltene-solvent effects deriving from structures containing heteroatoms,
varied aliphatic chains with different lengths and so forth.
Finally, the outcome of the project transcends the theoretical contributions to the petroleum
industry and to the academic field. It represents an unique opportunity for students to fully experi-
ence the atmosphere of scientific research and from this, develop communication skills in all senses,
thinking critically and sharpen the ability of taking coherent decisions. In other words, it was a very
rewarding introduction to a professional environment, producing outcomes that are highly evaluated
for corporations and research centres.
April 27, 2015 20
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
6 Nomenclature
Epot - Potential energy of simulation system;
fij - Scaling factor for inter- and intramolecular interactions;
Fi - Force on particle i;
kB - Boltzmann constant;
mi - Atomic mass;
M - Molecular mass in the bulk;
NA - Number of particles or molecules A, depending on the configuration of the radial distribution
function;
qi - Particle’s charge;
ri - Position vector of particle i;
rij - Distance between the centres of each particle;
T - Temperature;
V - Bulk volume;
ρ - Density;
σij - Van der Waals radius. Distance where intermolecular potential between particles is zero;
�ij - Dielectric constant. Measures the intensity of attraction between particles.
April 27, 2015 21
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
7 References
Accelrys draw 4.2 - academic version. (2007). San Diego: Accelrys Software Inc. Retrieved from
http://accelrys.com/products/informatics/cheminformatics/draw/no-fee.php
Ali, L. H., Al-Ghannam, K. A., & Al-Rawi, J. M. (1990). Chemical structure of asphaltenes
in heavy crude oils investigated by n.m.r. Fuel, 69 (4), 519 - 521. Retrieved from http://
www.sciencedirect.com/science/article/pii/001623619090326L
Allen, M. (2004). Introduction to molecular dynamics simulation. In Computational soft matter: from
synthetic polymers to proteins: Lecture notes. (p. 1-28). NIC.
Allen, P., & Tildesley, D. (1987). Computer simulation of liquids. Clarendon Press.
Alvarez, G.,Poteau, S., Argillier, J.-F., Langevin, D., & Salager, J.-L. (2009). Heavy oil-water
interfacial properties and emulsion stability: Influence of dilution. Energy & Fuels, 23 (1),
294-299. Retrieved from http://dx.doi.org/10.1021/ef800545k
Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F., DiNola, A., & Haak, J. R. (1984).
Molecular dynamics with coupling to an external bath. The Journal of Chemical Physics, 81 (8),
3684-3690. Retrieved from http://dx.doi.org/10.1063/1.448118
Boek, E. S., Yakovlev, D. S., & Headen, T. F. (2009). Quantitative molecular representation of
asphaltenes and molecular dynamics simulation of their aggregation†. Energy & Fuels, 23 (3),
1209-1219. Retrieved from http://dx.doi.org/10.1021/ef800876b
Burya, Y. G., Yudin, I. K., Dechabo, V. A., Kosov, V. I., & Anisimov, M. A. (2001, Aug). Light-
scattering study of petroleum asphaltene aggregation. Appl. Opt., 40 (24), 4028–4035. Retrieved
from http://ao.osa.org/abstract.cfm?URI=ao-40-24-4028
Bussi, G., Donadio, D., & Parrinello, M. (2007). Canonical sampling through velocity rescaling.
The Journal of Chemical Physics, 126 (1), 1-8. Retrieved from http://dx.doi.org/10.1063/
1.2408420
Caleman, C., van Maaren, P. J., Hong, M., Hub, J. S., Costa, L. T., & van der Spoel, D. (2012). Force
field benchmark of organic liquids: Density, enthalpy of vaporization, heat capacities, surface
tension, isothermal compressibility, volumetric expansion coefficient, and dielectric constant.
Journal of Chemical Theory and Computation, 8 (1), 61-74. Retrieved from http://dx.doi.org/
10.1021/ct200731v
Darden, T., York, D., & Pedersen, L. (1993). Particle mesh ewald: An nlog(n) method for ewald
sums in large systems. The Journal of Chemical Physics, 98 (12), 10089-10092. Retrieved from
http://dx.doi.org/10.1063/1.464397
Diallo, M., Cagin, T., Faulon, J., & Goddard III, W. (2000). Chapter 5 thermodynamic properties
of asphaltenes: A predictive approach based on computer assisted structure elucidation and
atomistic simulations. In T. F. Yen & G. V. Chilingarian (Eds.), Asphaltenes and asphalts,
2 (Vols. 40, Part B, p. 103 - 127). Elsevier. Retrieved from http://dx.doi.org/10.1016/
S0376-7361(09)70276-6
Gao, F., Xu, Z., Liu, G., & Yuan, S. (2014). Molecular dynamcis simulation: The behavior of
asphaltene in crude oil and at the oil/water interface. Energy & Fuels, 28(12), 7368–7376.
Garcia Barneto, A., Carmona, J. A., & Barrón, A. (2015). Thermogravimetric monitoring of
crude oil and its cuts in an oil refinery. Energy & Fuels, 29 (4), 2250-2260. Retrieved from
http://dx.doi.org/10.1021/ef5028795
Hansen, J.-P., & McDonald, I. R. (2013). Chapter 2 - statistical mechanics. In J.-P. Hansen &
I. R. McDonald (Eds.), Theory of simple liquids (fourth edition) (Fourth Edition ed., p. 13 - 59).
Oxford: Academic Press. Retrieved from http://dx.doi.org/10.1016/B978-0-12-387032-2
.00002-7
Headen, T. F., Boek, E. S., & Skipper, N. T. (2009). Evidence for asphaltene nanoaggregation in
toluene and heptane from molecular dynamics simulations. Energy & Fuels, 23 (3), 1220-1229.
Retrieved from http://dx.doi.org/10.1021/ef800872g
Headen, T. F., Howard, C. A., Skipper, N. T., Wilkinson, M. A., Bowron, D. T., & Soper, A. K.
(2010). Structure of pi − pi interactions in aromatic liquids. Journal of the American Chemical
Society, 132 (16), 5735-5742. Retrieved from http://dx.doi.org/10.1021/ja909084e
Hess, B., Bekker, H., Berendsen, H. J. C., & Fraaije, J. G. E. M. (1997). Lincs: A linear con-
straint solver for molecular simulations. Journal of Computational Chemistry, 18 (12), 1463-
1472. Retrieved from http://dx.doi.org/10.1002/(SICI)1096-987X(199709)18:12<1463::
AID-JCC4>3.0.CO;2-H
Humphrey, W., Dalke, A., & Schulten, K. (1996). VMD – Visual Molecular Dynamics. Journal of
Molecular Graphics, 14 , 33-38.
April 27, 2015 22
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Jorgensen, W. L., Maxwell, D. S., & Tirado-Rives, J. (1996). Development and testing of the opls
all-atom force field on conformational energetics and properties of organic liquids. Journal of
the American Chemical Society, 118 (45), 11225-11236. Retrieved from http://dx.doi.org/
10.1021/ja9621760
Khoshandam, A., & Alamdari, A. (2010). Kinetics of asphaltene precipitation in a heptane-toluene
mixture. Energy & Fuels, 24 (3), 1917-1924. Retrieved from http://dx.doi.org/10.1021/
ef9012328
Kokal, S. (2007). Chapter 12: Crude oil emulsions. In Petroleum engineering handbook. Society of
Petroleum Engineers.
Kuznicki, T., Masliyah, J. H., & Bhattacharjee, S. (2008). Molecular dynamics study of model
molecules resembling asphaltene-like structures in aqueous organic solvent systems. Energy &
Fuels, 22 , 2379-2389.
McLean, J. D., & Kilpatrick, P. K. (1997). Effects of asphaltene solvency on stability of water-in-
crude-oil emulsions. Journal of Colloid and Interface Science, 189 (2), 242 - 253. Retrieved from
http://dx.doi.org/10.1006/jcis.1997.4807
Morimoto, M., Boek, E. S., Hibi, R., Matsuoka, T., Uetani, T., Murata, S., & . . . Honda, H. (2014).
Investigation of asphaltene-asphaltene association and aggregation for compositional reservoir
simulators by quantitative molecular representations. International Petroleum Technology
Conference. Retrieved from https://www.onepetro.org/conference-paper/IPTC-18097-MS
Mullins, O. C., Sabbah, H., Eyssautier, J., Pomerantz, A. E., Barré, L., Andrews, A. B., . . . Zare,
R. N. (2012). Advances in asphaltene science and the yen–mullins model. Energy & Fuels, 26 (7),
3986-4003. Retrieved from http://dx.doi.org/10.1021/ef300185p
Mullins, O. C., Sheu, E. Y., Hammani, A., & Marshall, A. G. (2007). Asphaltenes, heavy oils, and
petroleomics. New York: Springer.
Parrinello, M., & Rahman, A. (1981). Polymorphic transitions in single crystals: A new molecular
dynamics method. Journal of Applied Physics, 52 (12), 7182-7190. Retrieved from http://
dx.doi.org/10.1063/1.328693
Pronk, S., Páll, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., . . . Lindahl, E. (2013).
Gromacs 4.5: a high-throughput and highly parallel open source molecular simulation toolkit
(Vol. 29) (No. 7). Retrieved from http://dx.doi.org/10.1093/bioinformatics/btt055
Pubchem compound database; cid=1140. (2015). National Center for Biotechnology Information.
Retrieved from http://pubchem.ncbi.nlm.nih.gov/compound/1140
Ribeiro, A. A. S. T., Horta, B. A. C., & Alencastro, R. B. d. (2008). Mktop: a program for automatic
construction of molecular topologies. Journal of the Brazilian Chemical Society, 19 (7), 1433-
1435. Retrieved from http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0103
-50532008000700031&nrm=iso
Schüttelkopf, A. W., & van Aalten, D. M. F. (2004). PRODRG: a tool for high-throughput crystallog-
raphy of protein-ligand complexes. Acta Crystallogr., D60 , 1355-1363.
Sedghi, M., Goual, L., Welch, W., & J., K. (2013). Effect of asphaltene structure on association and
aggregation using molecular dynamics. Journal of Physical Chemestry, 117 , 5765-5776.
Sheremata, J. M., Gray, M. R., Dettman, H. D., & McCaffrey, W. C. (2004). Quantitative molecular
representation and sequential optimization of athabasca asphaltenes. Energy & Fuels, 18 (5),
1377-1384. Retrieved from http://dx.doi.org/10.1021/ef049936+
Spiecker, P., Gawrys, K. L., & Kilpatrick, P. K. (2003). Aggregation and solubility behavior of
asphaltenes and their subfractions. Journal of Colloid and Interface Science, 267 (1), 178 - 193.
Retrieved from http://dx.doi.org/10.1016/S0021-9797(03)00641-6
Subirana, M., & Sheu, E. (2013). Asphaltenes: Fundamentals and applications. Springer US.
Tang, W., Hurt, M. R., Sheng, H., Riedeman, J. S., Borton, D. J., Slater, P., & Kenttämaa, H. I.
(2015). Structural comparison of asphaltenes of different origins using multi-stage tandem mass
spectrometry.Energy & Fuels, 29 (3), 1309-1314. Retrieved from http://dx.doi.org/10.1021/
ef501242k
Teklebrhan, R. B., Ge, L., Bhattacharjee, S., Xu, Z., & Sjöblom, J. (2012). Probing struc-
ture–nanoaggregation relations of polyaromatic surfactants: A molecular dynamics simulation
and dynamic light scattering study. The Journal of Physical Chemistry B, 116 (20), 5907-5918.
Retrieved from http://dx.doi.org/10.1021/jp3010184
Thompson, M. (2004). Arguslab 4.0.1. planaria software. LLC, Seattle, W.A..
Tukhvatullina, A., Barskaya, E., Kouryakov, V., Ganeeva, Y., Yusupova, T., & Romanov, G. (2013).
Supramolecular structures of oil systems as the key to regulation of oil behavior. Journal of
Petroleum & Environmental Biotechnology, 4 (4), 1-8. Retrieved from http://dx.doi.org/
April 27, 2015 23
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
10.4172/2157-7463.1000152
van der Spoel, D., Lindahl, E., Hess, B., van Buuren, A. R., Apol, E., Meulenhoff, P. J., . . . Berendsen,
H. J. C. (2013). Gromacs user manual version 4.6.5 [Computer software manual].
Wang, J., van der Tuuk Opedal, N., Lu, Q., Xu, Z., Zeng, H., & Sjöblom, J. (2012). Probing
molecular interactions of an asphaltene model compound in organic solvents using a surface
forces apparatus (sfa). Energy &/ Fuels, 26 (5), 2591-2599. Retrieved from http://dx.doi.org/
10.1021/ef201366y
April 27, 2015 24
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Appendix A Samples of MD files
A.1 Protein Database file (.pdb)
This file format is used in most molecular simulation softwares and contain the tridimensional
coordinates of each atom in the molecule. The bonds between particles are determined either by the
distances from one to the other or based on the CONECT options that might be provided, just as in
the example below.
;TITLE .PDB file for SC0 model
HETATM 1 CAO SC0 1 16.190 −6.140 −0.270 1.00 20.00 C
HETATM 2 HAO SC0 1 16.135 −6.008 −1.351 1.00 20.00 H
HETATM 3 HAP SC0 1 17.071 −5.677 0.174 1.00 20.00 H
HETATM 4 CAP SC0 1 14.930 −5.460 0.310 1.00 20.00 C
HETATM 5 HAQ SC0 1 15.071 −5.433 1.390 1.00 20.00 H
HETATM 6 HAR SC0 1 14.832 −4.469 −0.134 1.00 20.00 H
HETATM 7 CAM SC0 1 13.750 −6.190 0.050 1.00 20.00 C
HETATM 8 CAC SC0 1 13.780 −7.610 0.020 1.00 20.00 C
HETATM 9 CAB SC0 1 15.030 −8.280 0.000 1.00 20.00 C
HETATM 10 CAG SC0 1 16.230 −7.530 −0.010 1.00 20.00 C
HETATM 11 CAH SC0 1 17.470 −8.200 0.010 1.00 20.00 C
HETATM 12 HAH SC0 1 18.400 −7.630 0.050 1.00 20.00 H
HETATM 13 CAJ SC0 1 17.510 −9.610 0.030 1.00 20.00 C
HETATM 14 HAJ SC0 1 18.470 −10.120 0.050 1.00 20.00 H
HETATM 15 CAI SC0 1 16.310 −10.350 0.010 1.00 20.00 C
HETATM 16 HAI SC0 1 16.360 −11.440 0.000 1.00 20.00 H
HETATM 17 CAD SC0 1 15.060 −9.700 0.000 1.00 20.00 C
HETATM 18 CAF SC0 1 13.860 −10.440 0.000 1.00 20.00 C
HETATM 19 HAF SC0 1 13.890 −11.530 0.000 1.00 20.00 H
HETATM 20 CAE SC0 1 12.630 −9.770 0.010 1.00 20.00 C
HETATM 21 HAE SC0 1 11.700 −10.350 0.000 1.00 20.00 H
HETATM 22 CAA SC0 1 12.580 −8.360 0.020 1.00 20.00 C
HETATM 23 CAK SC0 1 11.350 −7.680 0.010 1.00 20.00 C
HETATM 24 HAK SC0 1 10.410 −8.240 0.010 1.00 20.00 H
HETATM 25 CAL SC0 1 11.310 −6.280 0.000 1.00 20.00 C
HETATM 26 HAL SC0 1 10.350 −5.760 −0.010 1.00 20.00 H
HETATM 27 CAN SC0 1 12.510 −5.530 0.020 1.00 20.00 C
HETATM 28 HAN SC0 1 12.460 −4.440 0.000 1.00 20.00 H
CONECT 1 2 3 4 10
CONECT 2 1
CONECT 3 1
CONECT 4 1 5 6 7
CONECT 5 4
CONECT 6 4
CONECT 7 4 8 27
CONECT 8 7 9 22
CONECT 9 8 10 17
CONECT 10 1 9 11
CONECT 11 10 12 13
CONECT 12 11
CONECT 13 11 14 15
CONECT 14 13
CONECT 15 13 16 17
CONECT 16 15
CONECT 17 9 15 18
CONECT 18 17 19 20
CONECT 19 18
CONECT 20 18 21 22
CONECT 21 20
CONECT 22 8 20 23
CONECT 23 22 24 25
April 27, 2015 25
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
CONECT 24 23
CONECT 25 23 26 27
CONECT 26 25
CONECT 27 7 25 28
CONECT 28 27
END
A.2 Molecular Dynamics Parameters file (.mdp)
The .mdp files contain all the simulation parameters that will be read by GROMACS. Each
equilibration ensemble and MD simulation needs information about how the interactions will be
computed. More details can be found in the GROMACS documentation (van der Spoel et al., 2013).
t i t l e = Asphaltene s imu la t i on s − MDrun step
; Run parameters
t i n i t = 0 ; i n i t i a l time
i n t e g r a t o r = md ; leap−f r o g i n t e g r a t o r
nsteps = 10000000 ; 2 ∗ 10000000 = 20000 ps , 20 ns
dt = 0.002 ; s i z e o f time s t ep s = 2 f s
; Output c on t r o l
nstxout = 1000 ; save coo rd ina t e s every 2 ps
nstvout = 1000 ; save v e l o c i t i e s every 2 ps
nstxtcout = 1000 ; xtc compressed t r a j e c t o r y output every 2 ps
nstenergy = 1000 ; save en e r g i e s every 2 ps
n s t l o g = 1000 ; update l og f i l e every 2 ps
; Bond parameters
cont inuat i on = yes ; Res ta r t ing a f t e r NPT
const ra int_a lgor i thm = l i n c s ; holonomic c on s t r a i n t s
c on s t r a i n t s = a l l−bonds ; a l l bonds ( even heavy atom−H bonds ) cons t ra ined
l i n c s_ i t e r = 1 ; accuracy o f LINCS
l inc s_orde r = 4 ; a l s o r e l a t e d to accuracy
; Ne ighborsearch ing
ns_type = gr id ; search ne ighbor ing g r id c e l l s
n s t l i s t = 10 ; 10 f s
r l i s t = 1 .0 ; short−range n e i g h b o r l i s t c u t o f f ( in nm)
rcoulomb = 1 .0 ; short−range e l e c t r o s t a t i c c u t o f f ( in nm)
rvdw = 1.0 ; short−range van der Waals c u t o f f ( in nm)
; E l e c t r o s t a t i c s
coulombtype = PME ; Pa r t i c l e Mesh Ewald f o r long−range e l e c t r o s t a t i c s
pme_order = 4 ; cub ic i n t e r p o l a t i o n
f o u r i e r s p a c i n g = 0.16 ; g r id spac ing f o r FFT
; Temperature coup l ing i s on
tcoup l = V−r e s c a l e ; modi f i ed Berendsen thermostat
tc−grps = system ; one coup l ing group − app l i e s f o r the whole system
tau_t = 0 .1 ; time constant , in ps
re f_t = 298 ; r e f e r e n c e temperature , in K
; Pressure coup l ing i s on
pcoupl = Par r i n e l l o−Rahman ; Pressure coup l ing on in NPT
pcoupltype = i s o t r o p i c ; uniform s c a l i n g o f box vec to r s
tau_p = 1 .0 ; time constant , in ps
ref_p = 1 .0 ; r e f e r e n c e pres sure , in bar
c omp r e s s i b i l i t y = 4 .5 e−5 ; Defau l t va lue f o r GROMACS, bar^−1
; Pe r i od i c boundary cond i t i on s
pbc = xyz ; 3−D PBC
; Di spe r s i on c o r r e c t i o n
DispCorr = EnerPres ; account f o r cut−o f f vdW scheme
; Ve loc i ty gene ra t i on
gen_vel = no ; Ve loc i ty gene ra t i on i s o f f
April 27, 2015 26
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Appendix B Insertion of charges in topologies
The following figure shows how charges were defined for carbon atoms in the asphaltene structure
based on the OPLS/AA parameters.
Figure 25: Scheme for charges’ insertion
(Morimoto et al., 2014)
April 27, 2015 27
Molecular Modelling and Simulation of Asphaltenes João Luís Lemos Fontes Silva Costa
Appendix C Extra simulations
C.1 Water-toluene box
The box showed below was one of the first simulations performed to study the interaction between
the solvents. From figure 26, one can notice that the solvents remained immiscible, due to their
different polarities. This simulation box is comparable with those used in the study the stabilization of
oil-water emulsion caused by asphaltenes. Hence, it could be easily used in next studies in view of
impacts of asphaltenes on that matter.
Figure 26: Simulation box containing water (red and white) and toluene (blue) molecules.
C.2 Toluene-heptane mixture
This subsequent box was designed especially to test the behaviour of the heptane structures
developed. In theory these two organic solvents should be perfectly miscible between each other,considering they are very apolar solvents. However, when in presence of asphaltenes they might take
much longer to reach an equilibrium phase than what is expected. Again, this simulation environment
could be used to study whether asphaltenes can truly affect the miscibility of these solvents. It was
observed from the simulation that the heptane models were miscible with toluene.
Figure 27: Mixed toluene (blue) and heptane (green) molecules in a simulation box.
April 27, 2015 28
	Introduction
	Experimental Section
	GROMACS simulation package
	Potential Functions
	Force Field
	Periodic Boundary Conditions
	Radial Distribution Functions
	Molecular Models
	Initial Simulation Configuration
	Details of Molecular Dynamics
	Analysis tools
	Results
	Small core models
	Small Cores in Toluene
	Small Cores in Heptane
	Big core models
	Big Cores in Toluene
	Big Cores in Heptane
	Aggregation events
	Structural analysis of models
	Reflection & Review
	Conclusions
	Nomenclature
	References
	Appendix Samples of MD files
	Protein Database file (.pdb)
	Molecular Dynamics Parameters file (.mdp)
	Appendix Insertion of charges in topologies
	Appendix Extra simulations
	Water-toluene box
	Toluene-heptane mixture