Energy Fuels 2011, 25, 499–502 Published on Web 12/13/2010
: DOI:10.1021/ef1010385
)
Potential of Mean Force Calculation from Molecular Dynamics Simulation of Asphaltene Molecules on a Calcite Surface† Thomas F. Headen‡,§ and Edo S. Boek*,§,
‡ Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, United Kingdom, Schlumberger Cambridge Research, High Cross, Madingley Road, Cambridge CB3 0EL, United Kingdom, and Department of Chemical Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, United Kingdom )
§
Received August 7, 2010. Revised Manuscript Received November 6, 2010
In this paper, we present preliminary calculations of the potential of mean force (PMF) from molecular dynamics simulations of asphaltene molecules on a calcite surface. We calculate, for the first time, the PMF between an asphaltene molecule and the calcite surface directly by a constraint force method. The asphaltene molecule is obtained in a systematic way consistent with experimental information, using the quantitative molecular representation approach. The calcite surface has been defined as the {10.4} face, which is representative of the experimentally dominant crystal face. First, we calculate the time-averaged values of the constraint force as a function of the z distance between the asphaltene molecule and the calcite surface. We observe that, at a separation of 3 A˚, the value of the force is positive, which corresponds to a repulsive interaction. At 3.5 A˚, a minimum in the force is observed, corresponding to an attractive interaction. At longer distances, the average constraint force goes to zero, according to expectations. Then, we calculate the integral of the time-averaged constraint force as a function of the distance between the asphaltene molecule and the calcite surface. The difference between the minimum energy and the energy at large separations gives an estimate of the free energy of adsorption. In this case, the value of the free energy is of the order of 110 kJ/mol. This is a reasonable value for simulations in vacuo corresponding to an effective solvent.
seen by small-angle neutron scattering (SANS),4,5 and, therefore, cannot be interpreted as a single molecule. For this reason, we carry out well-controlled computer experiments to calculate the forces between a single well-defined asphaltene molecule and a well-defined mineral surface, in this case, calcite. In this paper, we report on molecular dynamics (MD) simulations of the adsorption of an asphaltene structure on a crystalline calcite surface as a model for a carbonate rock. The asphaltene molecule structure chosen for this preliminary study was calculated using the quantitative molecular representation (QMR) method by Boek et al.,6 which has shown aggregation behavior through MD simulations in toluene and heptane. These simulations may eventually help interpret AFM and surface force apparatus (SFA) experimental results, although we will not be making a direct comparison in the current paper.
Introduction It is our aim in this study to investigate the adsorption of asphaltenes on a model calcite surface. The adsorption of asphaltenes on mineral surfaces is of great importance to oil recovery applications, because it affects the wettability of the reservoir rock (e.g., see ref 1). Several atomic force microscopy (AFM) studies have appeared in recent years, measuring force-distance relationships between asphaltene aggregates and various model surfaces.2,3 However, the interpretation of the experimental data is sometimes unclear. This is due to the strong propensity for asphaltene molecules to aggregate, and therefore, forces observed may be due to aggregate-aggregate forces rather than molecule-molecule forces. For example, Long et al.2 studied the force between asphaltene aggregates and a mica surface, interpreting the resultant force-distance curve and stretching of an asphaltene aggregate. The maximum force is observed at a distance between 250 and 300 nm. This is many times the size of asphaltene nano-aggregates, as
MD Simulation Methods Classical MD uses well-defined classical intra- and intermolecular potentials to calculate interatomic forces. The system is allowed to evolve over time by stepwise integration of the equations of motion. It is important that the step size be smaller than the time period of the fastest motion in the simulation. Rigid bond lengths were used to remove the fastest moving molecular motions; a time step of 2 fs was used for all simulations. Bond
† Presented at the 11th International Conference on Petroleum Phase Behavior and Fouling. *To whom correspondence should be addressed. E-mail: e.boek@ imperial.ac.uk. (1) Buckley, J. S.; Liu, Y.; Xie, X.; Morrow, N. R. Asphaltenes and crude oil wetting;The effect of oil composition. SPE J. 1997, 2, 107–119. (2) Long, J.; Xu, Z.; Jacob, H. M. Single molecule force spectroscopy of asphaltene aggregates. Langmuir 2007, 23, 6182–6190. (3) Liu, J.; Zhang, L.; Xu, Z.; Masliyah, J. Colloidal interactions between asphaltene surfaces in aqueous solutions. Langmuir 2006, 22, 1485–1492. (4) Roux, J. N; Broseta, D.; Deme, B. SANS study of asphaltene aggregation: Concentration and solvent quality effects. Langmuir 2001, 17, 5085–5092.
r 2010 American Chemical Society
(5) Headen, T. F.; Boek, E. S.; Stellbrink, J.; Scheven., U. M. Small angle neutron scattering (SANS and V-SANS) study of asphaltene aggregates in crude oil. Langmuir 2009, 25, 422–428. (6) Boek, E. S.; Yakovlev, D. S.; Headen, T. Quantitative molecular representation of asphaltenes and molecular dynamics simulation of their aggregation. Energy Fuels 2009, 23, 1209–1219.
499
pubs.acs.org/EF
Energy Fuels 2011, 25, 499–502
: DOI:10.1021/ef1010385
Headen and Boek
angles and dihedral angles remained flexible (within appropriate force fields) to allow for molecular flexibility. For this study, the GROMACS MD code7 with the optimized potential for liquid simulations-all-atom (OPLS-AA) force-field parameters8 was used. The OPLS force field is an AA force field, which has been shown to work well for aromatic liquids in reproducing experimental data.8 AA force fields provide parameters for every atom in a system, including hydrogen, while “united-atom” (UA) force fields treat the hydrogen and carbon atoms in methyl and methylene groups as a single interaction center. A distinctive feature of the OPLS parameters is that they were optimized to fit experimental properties of liquids, such as density and heat of vaporization.8 For more details regarding the MD method and interaction potentials, we refer to our previous MD work on asphaltenes,9 GROMACS online manual,10 and textbooks, such as ref 11. One of the primary difficulties in simulating mineral-organic systems is defining a force field that is consistent with the interactions between minerals and organics. In this preliminary study, the 6-12 Leonard-Jones parameters for carbonate were approximated as the OPLS parameters for carboxylic acid; partial charges were taken from a specific calcite force field by Pavese et al.12 This approximation is satisfactory for a preliminary investigation; however, future work should include the generation of a new force field that is consistent with both mineral and organic force fields, following the example by Freeman et al., who have produced a classical force field for organics with calcite for simulation of biomineralization.13 The Freeman method uses existing potentials for the components of the system and produces cross-term potentials between these components. In this paper, we use standard combination rules for cross-term potential parameters (see ref 11). Periodic boundary conditions with the minimum image convention are used, so that a small box of ∼10 000 atoms can represent the bulk. Long-range Coulomb intermolecular forces are treated using the particle-mesh Ewald (PME) technique,14 which allows for the use of fast Fourier transforms (FFTs). One point of discussion is the use of carboxylic acid force-field parameters for carbonate. We argue that, because the OPLS force field is designed for organic molecules, there are no standard potentials for carbonate; therefore, as a starting point, carboxylic acid 6-12 potentials were used as the closest functional group. It is expected that the electrostatic component will be the strongest interaction. Therefore, the use of a similar if different 6-12 potential should not affect things too greatly.
Figure 1. Asphaltene structure used in the potential of mean force (PMF) calculations, obtained from QMR calculations in the ref 6. Oxygen is red, and sulfur is yellow. Hydrogen atoms are not shown in the structure for clarity.
possibilities. The free energy of adsorption will strongly depend upon the particular surface structure. It has been shown that the {10.4} face is the dominant crystal face on the macroscopic calcite crystal. Therefore, we have carefully made a cut through the crystal structure to obtain this representative surface.15 Subsequently, we generated a crystal cell of 4 4 4, containing 64 CaCO3 units. Then, this cell was doubled in x and y directions to obtain an 8 8 4 cell, containing 256 calcite molecules. We selected an appropriate asphaltene structure from the QMR calculations described in ref 6. This structure was found to be in excellent agreement with experimental observations, including molecular weight, nuclear magnetic resonance (NMR), and elemental analysis. For more details regarding the QMR and MD methods, we recommend refs 6 and 9, respectively. The asphaltene structure used is shown in Figure 1. This structure was placed on the calcite surface. A snapshot of the system is shown in Figure 2. PMF We calculate the PMF by a constraint force method. A series of simulations has been conducted of one asphaltene molecule in an effective medium with a dielectric constant corresponding to a hydrocarbon solvent (toluene) adjacent to the calcite surface. The z distance between the center of mass of the asphaltene and the calcite surface is constrained using the SHAKE algorithm.16 The constraint distances used were from 0.3 to 1.5 nm in steps of 0.1 nm. At each time step during the simulation, the force required to keep the molecules at the constraint distance is calculated by the SHAKE algorithm. This is averaged over the 20 ns simulation to give a mean force between the asphaltene molecule and the surface over millions of molecular conformations. Integration of the forcedistance curve yields the PMF (plus a constant of integration). This provides a method of calculating the PMF and the free energy of asphaltene adsorption on the calcite surface. It is more computationally expensive than calculation of the PMF from the radial distribution function g(r); however, it does
Calcite and Asphaltene MD Simulation Cell It is possible to generate a calcite crystal surface in many different ways, by making cuts parallel to crystal planes {hkl}. It is very important to select the correct surface out of many (7) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free J. Comput. Chem. 2005, 26, 1701–1718. (8) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118 (45), 11225–11236. (9) Headen, T; Boek, E. S.; Skipper, N. T. Evidence for asphaltene nano-aggregation in toluene and heptane from molecular dynamics simulations. Energy Fuels 2009, 23, 1220–1229. (10) www.gromacs.org. (11) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1986. (12) Pavese, A.; Catti, M.; Parker, S. C.; Wall, A. Modelling of the thermal dependence of structural and elastic properties of calcite, CaCO3. Phys. Chem. Miner. 1996, 23, 89–93. (13) Freeman, C. L.; Harding, J. H.; Cooke, D. J.; Elliott, J. A.; Lardge, J. S.; Duffy, D. M. New forcefields for modeling biomineralization processes. J. Phys. Chem. C 2007, 111 (32), 11943–11951. (14) York, D. M.; Darden, T. A.; Pedersen, L. G. The effect of longrange electrostatic interactions in simulations of macromolecular crystals: A comparison of the Ewald and truncated list methods. J. Chem. Phys. 1993, 99 (10), 8345–8348.
(15) Nygren, M. A.; Gay, D. H.; Catlow, C. R. A.; Wilson, M. P.; Rohl, A. L. Incorporation of growth-inhibiting diphosphonates into steps on the calcite cleavage plane surface. J. Chem. Soc., Faraday Trans. 1998, 94 (24), 3685–3693. (16) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23 (3), 327–341.
500
Energy Fuels 2011, 25, 499–502
: DOI:10.1021/ef1010385
Headen and Boek
Figure 5. Time-averaged value of the constraint force as a function of the z distance between the asphaltene molecule and the calcite surface.
Figure 2. Simulation snapshot of the calcite surface and asphaltene molecule.
In Figure 3, we show the three-dimensional structure of the asphaltene molecule used in the PMF MD calculations (see ref 6). In Figure 4, we show a typical transient of the constraint force as a function of time for a constraint z distance of 3 A˚. The jump in the value of the constraint force corresponds to the optimization of the free energy of the conformations near the crystal surface. Please note the step change around t = 16 000 ps. Following the asphaltene molecule as a function of time (movie not shown), we believe that this is due to a part of the aliphatic chain “finding” the surface and staying there. In Figure 5, we present the time-averaged values of the constraint force as a function of the z distance between the asphaltene molecule and the calcite surface over the course of the simulation. It can be clearly seen that, at a separation of 3 A˚, the value of the force is positive, which corresponds to a repulsive interaction. At 3.5 A˚, a minimum in the force is observed, corresponding to an attractive interaction. At longer distances, the average constraint force goes to zero, according to expectations. On the basis of the measurements shown in Figure 5, it is justified to assume that, after 20 nm, the force between the surface and the asphaltene becomes negligible, especially because there is no direct contact between the surface and the asphaltene. However, it is reasonable that more “jumps” may occur at smaller distances. Our interpretation is that features are caused by different elements of the asphaltene molecule touching the surface. One potential example is that the sulfur on the aliphatic chain may have a higher affinity for the surface than the carbon atoms because of higher electrostatic interactions. The error bars are included in the figure and are small in comparison to the value measured. In Figure 6, we present the integral of the time-averaged constraint force as a function of the distance between the asphaltene molecule and the calcite surface. The relation between Figures 5 and 6 is as follows. At each point along the reaction coordinate z, the ensemble averaged force because of the constraint is evaluated. This force is calculated as the Lagrange undetermined multiplier Æ f æ, found by solving the equation of motion of the system with the constraint. This force is shown in Figure 5. For the simple distance constraint
Figure 3. Three-dimensional structure of asphaltene molecule used in the PMF calculations (see ref 6).
Figure 4. Constraint force as a function of time for a constraint z distance of 3 A˚ between the asphaltene and the calcite surface.
have the advantage that it samples more of the ensemble in high-energy configurations, such as when the asphaltene and the calcite surface are held very close together. It can, therefore, provide a more accurate estimate for the PMF at distances that are not well-sampled during an unconstrained simulation. 501
Energy Fuels 2011, 25, 499–502
: DOI:10.1021/ef1010385
Headen and Boek
and -Æ f æ are shown in Figure 6. For more details, we refer to ref 17. The difference between the minimum energy and the energy at large separations gives an estimate of the free energy of adsorption. In this case, the value of the free energy is of the order of -110 kJ/mol. This is a reasonable value for simulations in vacuo corresponding to an effective solvent. A similar result of -130 kJ mol-1 for asphaltene-asphaltene binding energies has been observed from molecular mechanics results in vacuo.18 This value is expected to decrease in a simulation using an explicit solvent. We are planning to carry out simulations to verify this expectation and compare aggregation on the surface in a range of solvents and precipitants. Conclusions This preliminary study has described a method for the calculation of the free energy of aggregation of an asphaltene molecule on a calcite surface. The free energy of adsorption was calculated to be -110 kJ mol-1, which is reasonable for the approximate force fields used when compared to asphaltene-asphaltene aggregation in vacuo. Further work is to include the generation of more accurate force fields for organic-mineral interactions. Please note that we have neglected interactions among adsorbate molecules and this is an extension to be considered in future simulations. Also, the effect of using an explicit solvent, instead of an effective dielectric medium, has to be investigated to give more precise values for binding forces and energies between asphaltene molecules and calcite surfaces. Ultimately, this will allow for a comparison to AFM and SFA data to help interpret the nature of asphaltene aggregation on the calcite surface.
Figure 6. Integral of the time-averaged constraint force as a function of the distance between the asphaltene and the calcite surface (black curve). The green curve corresponds to the negative value of the constraint force.
that we are using in this study, the free energy change is given by dF ¼ - Æ f æz dz The free energy increments are obtained by integrating the average force along the reaction coordinate. The free energy (17) Suter, J.; Boek, E. S.; Sprik, M. Adsorption of a sodium ion on a smectite clay from constrained ab initio molecular dynamics simulations. J. Phys. Chem. C 2008, 112, 18832–18839. (18) Ortega-Rodrı´ guez, A.; Cruz, S. A.; Gil-Villegas, A.; GuevaraRodrı´ guez, F.; Lira-Galeana, C. Molecular view of the asphaltene aggregation behavior in asphaltene-resin mixtures. Energy Fuels 2003, 17, 1100–1108.
Acknowledgment. The authors thank the Natural Environment Research Council, U.K., for funding and Schlumberger Cambridge Research for use of computer resources.
502
