Prediction of Saturation Densities and Critical Properties of n-Decane, n-Pentadecane, and n-Eicosane Using Molecular Dynamics with Different Force-Fields

Prediction of thermophysical properties of heavy hydrocarbons is important because of the recent increased interest in the extraction of heavy and shale oil to meet the global energy demand. Carrying out experimental work, such as determining the critical properties of heavy hydrocarbons, is challenging due to the possibility of thermal degradation during experimentation. This study focuses on the use of molecular simulations, specifically canonical molecular dynamics, to predict the critical properties of three hydrocarbons: n-decane (n-C10), n-pentadecane (n-C15), and n-eicosane (n-C20). The method uses volume-expansion molecular dynamics (VEMD), where a single box is enlarged in one axis and the vapor–liquid equilibrium is achieved. Three different force-fields (AMBER, COMPASS, and TraPPE) were employed to compare the accuracy of the simulated results and the required computational time. The results from the simulations were compared with available experimental data, equations-of-state, and several correlations. The results indicate that TraPPE is the most accurate and efficient force-field to predict the critical properties followed by AMBER and COMPASS.


INTRODUCTION
Rapid depletion of conventional oil resources and the constant increase in world energy demand have necessitated practical and research investigations into the extraction and utilization of unconventional oil resources such as shale or heavy oil. 1,2 Unconventional oil resources are much more abundant in comparison with conventional oil resources, leading to an increased incentive in their production in the oil and gas industries. For instance, China is estimated to possess unconventional oil resources exceeding 6×10 9 tons. 3 Also, the unconventional resources are rich in Canada and Venezuela, where each country is estimated to possess 2−3 trillion barrels of unconventional oil. 4,5 Due to the abundance of unconventional oil resources and an increased interest in their production, estimation of different thermophysical properties of the individual constituents of heavy oil (i.e., heavy hydrocarbons) becomes of prime importance for the oil and gas industries. One such estimation of thermophysical properties is the determination of the critical properties of heavy hydrocarbons. Knowledge of the critical properties is essential as it allows calculation of various thermodynamic and volumetric properties that are used to design various processes for the production and refining of crude oil. 6 For example, equations-of-state (EOS) such as Peng−Robinson (PR) and Soave−Redlich−Kwong (SRK) are utilized in the petroleum industry to perform phase equilibrium calculations. However, for these EOS to be applied, accurate knowledge of the critical points of the heavy hydrocarbons is an important prerequisite. 7 Experiments involving heavy hydrocarbons are difficult to perform, especially in the case when the critical properties need to be determined. This is primarily because heavy hydrocarbons tend to decompose or thermally degrade at high temperatures and before the critical point is even attained. 6−8 Therefore, modeling techniques, such as the use of correlations or molecular simulations, are employed to predict the critical properties. In this context, molecular simulations are increasingly being used to predict the thermophysical properties of different species and mixtures at a range of temperatures and pressures. 9 Molecular simulations predict properties of a macroscopic system by following the interactions between a limited number of constituents (atoms, ions, or molecules) in a certain special region under a fixed (known a priori) external set of constraints (ensembles). The interactions between the different constituents can primarily be studied using either the molecular dynamics (MD) method (a deterministic method) or the Monte Carlo (MC) method (a stochastic approach).
In this work, molecular simulations are carried out to predict the critical properties of three n-alkanes, namely, n-decane (n-C10) as the reference compound, n-pentadecane (n-C15), and n-eicosane (n-C20). Although critical properties of n-decane and n-eicosane have been investigated using the TraPPE forcefield, 10 n-pentadecane has been studied in this work for the first time using AMBER, COMPASS, and TraPPE force-fields. Additionally, no studies exist in the open literature on the prediction of critical properties for n-decane and n-eicosane using AMBER and COMPASS force-fields. Traditionally, the Gibbs ensemble MC (GEMC) approach has been used to predict the critical properties of species. 11 However, because of their stochastic nature, MC simulations have several disadvantages compared to MD simulations. For example, it is difficult to model dense liquids using MC simulations as the acceptance ratios associated with the insertion and deletion moves for the system are very low. 12 Also, parallelization of MC simulations is not easy, which results in a very high computational time for heavy systems. 13 Hence, the simulations of the hydrocarbons in this work are mainly carried out using canonical (constant N, V, and T) MD, whereby the system has a constant number of particles, volume, and temperature throughout the simulation.

FORCE-FIELDS
An important aspect of molecular simulations that determines the accuracy of the results is the choice of the different interaction models or force-fields. Molecular simulations are derived from the formulations of statistical mechanics, and as a result, the outcome from the simulations is dependent on the total energy of the system. Calculation of kinetic energy is relatively straightforward; 14,15 however, it is the calculation of the potential energy that ultimately determines the accuracy of the results. Hence, in molecular simulations, one key factor that affects the outcome of the results is the force-field, which is used to calculate the potential energy among the different molecular species. 9 The total potential energy is usually decomposed into two main components as shown in eq 1 15 where E int is the intramolecular energy and E ext is the external or intermolecular energy. The intermolecular energy term is fundamentally the energy due to the interactions between the system's constituents. eq 1 can be further expanded according to a modern molecular potential energy model, as shown in eq 2. 9,16−18 Generally, the first four terms represent the energy from intramolecular interactions, while the last two are forms of intermolecular interactions.
In eq 2, key nonbonded interactions are the dispersionrepulsion energy between atoms of different molecules 15 or of atoms within the same molecule at a distance of three bonds or more. 12 Many potentials are available for calculating the dispersion-repulsion energy such as the Born−Mayer−Huggins potential (used when dealing with polar species), the Buckingham exp-6 model, and the most extensively used and popular model, which is the Lennard-Jones (LJ) 6−12 model. eq 3 presents the LJ 6−12 model 9 i k j j j j j j j j i k j j j j j j y { z z z z z z i k j j j j j j y { z z z z z z y { z z z z z z z z where r ij is the separation distance between force centers of particles i and j, ε ij is the energy at minimum E LJ , and σ ij is the distance at which the interparticle potential is zero. Force-fields are used in molecular simulations to carry out the potential energy calculations and contain all of the necessary information and parameters relating to atoms and groups that make up a molecule. Equation 2 represents a typical form of what is referred to as the class I force-field. Class I force-fields are typically designed based on the harmonic motion that is derived from Hooke's law. 9 This work utilizes three different force-fields to carry out the simulations and study their accuracy in predicting the critical properties of the studied hydrocarbons. These force-fields are summarized below.
2.1. AMBER. The Assisted Model Building and Energy Refinement (AMBER) force-field was developed by Cornell et al. 23 for conducting MD simulations on proteins, nucleic acids, and other organic molecules. AMBER is also an all-atom (AA) force-field where each atom of a molecule is defined. Their work showed good accuracy in reproducing the interaction energies, solvation, and conformational energies of small molecules.

COMPASS.
Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) was first developed by Sun 24 in 1998 and is also an AA force-field. The COMPASS force-field was primarily developed to study small organic and inorganic molecules and polymers.

TraPPE.
Transferable Potentials for Phase Equilibria (TraPPE) is a force-field developed for studying phase equilibria of n-alkanes, branched alkanes, ketones, aldehydes, and alcohols. 25 TraPPE is a united-atom (UA) force-field, i.e., atoms are combined to present a single "united" atom. For example, for CH 3 , instead of four atoms, a single representative atom would be used in the simulation.
Another force-field that is commonly used for alkanes is the Optimized Potential for Liquid Simulations (OPLS). However, this study does not utilizes OPLS for predicting the critical properties.

SIMULATION DETAILS
All simulations were performed on a computer system housing an Intel Core i7-6700 CPU and 8 GB of RAM. With the aid of the in-built boost system, the simulations were consistently run at a clock speed of 3.95 GHz instead of the base speed of 3.40 GHz. Additionally, the simulations were performed on a single core, thereby not taking advantage of the parallelization nature of MD simulations. Moreover, all simulations used a switching function for the calculations involving the intermolecular potential. The inner cutoff was set to 9.0 Å, while the outer cutoff was 12.0 Å. The inner cutoff was set according to the values provided by the system building software Moltemplate. 26 The outer cutoff was based on 3σ for the CH 2 group as provided by the TraPPE force-field (shown in Table 1). Although it is usually taken to be 2.5σ as recommended by Smit and Toxaevard, 27,28 the value was rounded off to provide a more accurate set of results. Furthermore, to visualize the results of the simulations, the Open Visualization Tool (OVITO) 29 was used. Table 1 depicts the Lennard-Jones parameters for the different force-fields employed.
To predict the critical properties using MD, the canonical MD was employed where the length of the simulation box was enlarged in one dimension. Furthermore, to study the effect of different potential models, three different force-fields were used to predict the critical properties of the hydrocarbons as stated before. Figure 1 below depicts the n-decane molecule for the three force-fields.
Two software packages were used to conduct the simulations. The first software used was Moltemplate. 26 This software allows the building of a molecule according to the parameters of a forcefield and, consequently, a box containing the desired number of molecules. Such a box is shown in Figure 2 for n-decane that was made using the TraPPE force-field. The second software used was LAMMPS 30 to carry out the simulations on the system that was made using Moltemplate. For all simulations, a time step of 1 fs was employed and the number of molecules was set to 1000 for cases below the boiling point of the hydrocarbon and 8000 for the cases above the boiling point. This was done to ensure the formation of a stable interface at higher temperatures. 11 As shown in Figure 2, the box not only is very large but also contains very large spacing between the molecules. To produce a more realistic system, the following steps, as provided by Moltemplate, were used as a reference. 26 1. The system was run for 5 ps using the microcanonical (constant N, V, and E) and 50 ps at canonical ensembles at a temperature of 900 K to reorient the molecules and change their direction to a more natural one. 2. To reduce the size of the box, the system was then run in the isothermal−isobaric (constant N, P, and T) ensemble with a pressure of 250 atm for 250 ps. Additionally, during this run, the temperature was decreased from 900 K to the required desired temperature. 3. After the first isothermal−isobaric run at high pressure, a second isothermal−isobaric run was performed where the pressure was reduced from 250 to 1 atm, with the desired temperature being held constant, for 250 ps. 4. Lastly, a final isothermal−isobaric run was conducted at 1 atm and the desired temperature for 500 ps. This run allowed the system to reach a constant and realistic density at the given temperature and pressure. Performing the above steps allowed the molecules to have a more natural and realistic orientation. Moreover, it allowed the system to reach an equilibrium density as shown in Figure 3. The same Figure 3 also shows the density of the system at different stages of the simulation.
After the final NPT run, the box was then suddenly enlarged in the x-dimension using a factor of 3.5 (known as volumeexpansion molecular dynamics 11 ) and the system was run at the desired temperature for 0.75 to 3 ns. The employed method allows an equilibrium to be achieved between the vapor and liquid phases. The density of the system was continuously monitored to assess if the system has reached equilibrium. After equilibrium was attained, the final liquid and vapor densities were recorded.
After completion of the simulations, the critical temperature and density were then obtained by fitting the critical scaling law (eq 4) and the law of rectilinear diameters (eq 5), respectively. 31 where ρ is the density at a particular temperature T, and the subscripts l, v, and c represent the liquid, vapor, and critical properties, respectively. λ and γ are constants that are fit with the available simulation data. Here, 0.3265 is the critical exponent and is constant. 15,32 Although molecular simulations perform well in predicting the critical densities and temperatures, it is often quite difficult to where R is the gas constant, MW is the molecular weight, and T i is the temperature at the liquid density.

VALIDITY OF RESULTS
To verify the results obtained from the molecular simulations, comparisons were made with the experimental data and various correlations available in the literature. The critical properties predicted from the different force-fields were compared with experimental data from Ambrose. 35 The experimental data from Ambrose lacks the critical density for n-eicosane; hence, eq 6 was used to calculate the critical density using the experimental critical temperatures and pressures. Moreover, experimental vapor−liquid equilibria data for n-decane 36 was compared with the results of the simulations. In addition to the experimental data, results predicted by several correlations were also compared with the results obtained from the molecular simulations. Table 2 below presents the references for the employed correlations: The above correlations only provide a single point for comparison. To study the performance of the force-fields in predicting the vapor and liquid densities, three equations-ofstate (EOS) were also used, namely, the Peng−Robinson (PR), the Soave−Redlich−Kwong (SRK), and the Wilson EOS. 43 These EOS were used to provide an insight into the accuracy of the force-fields in predicting the vapor and liquid densities.

Critical Properties.
The critical properties were predicted using the AMBER, COMPASS, and TraPPE forcefields. Several correlations and EOS were also used for crosscomparison of the molecular simulation results. The coexistence curves for n-decane and n-pentadecane are presented in Figures  4 and 5, respectively, for the AMBER, COMPASS, and TraPPE force-fields.
The results for the coexistence curves of n-eicosane are presented in Figure 6. It is important to note that experimental data for the critical density was not found in the literature. Hence, eq 6 was used along with the experimental critical temperature and pressure and density data from 36 to calculate the critical density.
The critical pressures as obtained from Vetere's equation (eq 6), as well as the experimental data, are presented in Figure 7.
Lastly, the critical temperatures, densities, and pressures obtained from the simulations are presented in Tables 3−5, respectively. In these tables, the values in the square brackets represent the percentage deviation from the experimental values.
As stated before, for each of the three hydrocarbons studied in this work (n-C10, n-C15, and n-C20), three different force-fields were used to predict the critical properties, i.e., AMBER, TraPPE, and COMPASS. The results generated by molecular simulation in this work were compared to those of several correlations, EOS, and experimental data. Referring to Tables 3−5, it can be seen that TraPPE, followed by AMBER, provides the best predictions of the critical properties. For TraPPE, all percentage deviations are less than 10% with the prediction of the critical temperatures being less than 1%. The AMBER force-field also provides accurate results for the critical properties. However, it can be seen that,     From Figure 7 and Table 5, the values of critical pressure predicted by the force-fields were compared with the experimental results. Similar to the results of critical temperature and density, the TraPPE force-field provides the most accurate results, with a maximum percentage deviation of −9.22% for npentadecane. Although the COMPASS force-field provides poor  predictions for n-decane and n-pentadecane, it gives better accuracy in the case of n-eicosane when compared with AMBER.
To compare and assess the predictive capability for the molecular simulation techniques used in this study for the equilibrium (coexisting) liquid and vapor densities, results using the different force-fields were compared with those generated by the well-known EOS, i.e., Wilson, SRK, and PR. In the case of ndecane, experimental vapor−liquid data were also used for comparison.
In the case of AMBER, there seems to be good agreement between the simulated results and the experimental data (and even with results from PR and SRK to some extent) for n-decane for the vapor density (see Figure 4). However, AMBER underpredicts the liquid density when compared to the experimental data. Rather, the liquid density more closely follows the trend set by PR, specifically at high temperatures. In comparison to n-decane, the vapor densities for n-pentadecane ( Figure 5) and n-eicosane ( Figure 6) show poor agreement for all EOS. The liquid densities for n-pentadecane and n-eicosane, on the other hand, agree well with the SRK EOS at high temperatures, similar to n-decane for PR EOS.
For the TraPPE force-field, similar trends were observed for ndecane ( Figure 4) when compared to AMBER, where the vapor density follows closely the experimental data, PR, and SRK data. The liquid density from TraPPE, on the other hand, more closely follows the liquid density predicted by PR. Additionally, for n-pentadecane ( Figure 5), the vapor density is again in good agreement with the results from PR and SRK. While for the liquid density, the predictions using TraPPE match closely with those by the SRK EOS. However, for n-eicosane (Figure 6), again similar trends are seen to that of AMBER, where the vapor density is poorly predicted by TraPPE and the liquid density at high temperatures follows the trend set by PR EOS. Lastly, the COMPASS force-field severely underpredicts both the liquid and vapor densities for all hydrocarbons and, compared to AMBER and TraPPE, provides the most inaccurate results.
Considering the results obtained from molecular simulations, it is observed that the TraPPE force-field, followed by AMBER, provides the most accurate results when compared to the experimental data. They are also in the best agreement with results generated by other correlations and EOS. This may be attributed to the fact that TraPPE is designed to have specific parameters to predict accurate phase equilibria properties of different species (such as n-alkanes and ketones). For AMBER, a clear trend is observed where the percentage deviations seem to increase with the carbon number for all critical properties. Hence, the increase in percentage deviation could be attributed to the increase in crystallization in the system, which AMBER is unable to quantify accurately.
In this work, COMPASS was found to be the least accurate of the three employed force-fields, with inaccurate results obtained for critical properties and liquid and vapor densities. This could be attributed to the fact that COMPASS does not place equal emphasis between intramolecular and intermolecular parameters. This is in agreement with the results obtained in a study by Martin,44 where the most accurate results were produced by TraPPE followed by AMBER and COMPASS. It was also observed in that study, that for smaller molecules, AMBER predicts the vapor density much more accurately with a greater set of results showing a deviation of only 1%, which is similar to the trends observed in this work where AMBER produces the vapor density for n-decane similar to the experimental data. In another study, 45 it was observed that the TraPPE force-field predicts liquid densities that are in good agreement with the Lemmon and Huber EOS for n-dodecane. These results are similar to the results obtained in this work for other hydrocarbons.
Another aspect that cannot be ignored is the systematic error that might take place during the course of this study. As    46 given the same set of simulation tasks to different users, the results appear to show significant deviations, which may not be present when systematic errors are avoided and agree within statistical uncertainties. Factors like the simulation period, system size, and even the software used for a simulation can affect the final results. 46 In the case of VLE, Stephan et al. 47 have compared the data sets of a Lennard-Jones fluid from different studies using the three different tests: compressibility factor, the Clausius− Clapeyron equation, and the outlier test. In the study, it was shown that there can be a difference of up to ±1% for saturated vapor density between a set of different studies. 47 As a result, slightly different results are to be expected between the findings presented in this and other studies. It is also important to note that the results for critical pressures were obtained from the simulated results of critical temperature and critical density by means of eq 6. Because of this, the errors for critical pressure are much higher, as the small errors and inaccuracies in critical temperatures and densities are greatly amplified.

Computational Efficiency.
The computational efficiency takes into account the number of steps completed per day by a force-field. The measurement of the steps performed provides insights into the effectiveness of each force-field, as computational time is crucial in carrying out molecular simulations that provide accurate results. Additionally, knowing the number of steps that can be computed per day allows for planning the number of simulations that can be carried out in a given period of time.
Such a comparison is provided in Figure 8 for a system of 1000 n-decane molecules in an enlarged box (in one dimension) at 520 K. Figure 8 presents data in time steps of nanoseconds per day (ns/day) for all three force-fields: AMBER, COMPASS, and TraPPE.
An important aspect of molecular simulations is the execution time required for carrying out these simulations. In an effort to reduce the computational time, different force-fields were developed (such as the united-atom TraPPE force-field) that simplify the intermolecular calculations, thereby reducing the computational time. Hence, it becomes of value to compare the computational time required for the different force-fields used in this study to gauge the computational efficiency. However, this aspect must be considered while taking into account the accuracy of the generated simulation results.
The computational time greatly depends on the model used and also on the number of particles in the system. In this work, two different models were used to design the three employed force-fields. The first model is the all-atom model (AA), which was used for the AMBER and COMPASS force-fields, while the second is the united-atom model (UA) for the TraPPE forcefield. The results for the time performance are depicted in Figure  8.
From Figure 8, two trends are observed: (1) for all hydrocarbons, the TraPPE force-field provides the highest number of time steps computed followed by AMBER and COMPASS, and (2) with the increase in carbon number, the number of time steps for all force-fields tends to decrease. These trends can be attributed to the AA and UA nature of the forcefields. In the AA model, all atoms are modeled, which in the case of n-alkanes would include hydrogen and carbon atoms. The UA model, on the other hand, only looks at a group of atoms. For example, for n-alkanes, there will be two methyl groups, while the rest will be methylene groups. This greatly reduces the number of interaction sites for a UA model when compared to the AA model, leading to better computational performance. Based on the results obtained for critical properties and those presented in Figure 8, it can be concluded that TraPPE provides the most accurate results while being efficient in terms of the required computational time.

CONCLUSIONS
The primary method of predicting the critical properties of heavy hydrocarbons, or any species in general, is through the use of the Gibbs ensemble Monte Carlo. However, a major disadvantage of using MC techniques is the lack of parallelization in a computational system, which leads to increasing computational time. Hence, the use of MD to predict the critical properties is being extensively researched. This study, therefore, uses MD simulations to predict the critical properties of n-decane, n-pentadecane, and n-eicosane, using the AMBER, COMPASS, and TraPPE force-fields. Furthermore, the computational time of the three force-fields was also compared to provide an insight into the efficiency of the simulations. From the results generated in this work, it can clearly be seen that TraPPE is the most accurate force-field for predicting the critical properties with a maximum percentage deviation of −9.22%. The second most accurate force-field is AMBER, followed by COMPASS. Moreover, when compared to other correlations, it is observed that TraPPE fits well within the results of most correlations. Moreover, compared to AMBER and COMPASS, TraPPE also provides with the most amount of time steps per day, which results in TraPPE being the most efficient in predicting the critical properties. Additionally, when compared with EOS, it is observed that when the carbon number of the hydrocarbons increases (i.e., for heavier hydrocarbons), the liquid densities predicted by TraPPE and AMBER were in good agreement with those predicted by EOS. For n-decane, both TraPPE and AMBER provided good agreement with the experimental vapor density as well. However, this is not the case for the COMPASS force-field for all species, where it generated the worst prediction for vapor and liquid densities. Based on the findings in this study, it is recommended to use TraPPE and AMBER to predict and study the critical properties of heavy hydrocarbons.