{'fraction_non_alphanumeric': 0.03312452543659833, 'fraction_numerical': 0.02050113895216401, 'mean_word_length': 4.979757850926977, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Molecular Dynamics Simulations of Volumetric Thermophysical Properties of Natural Gases

The accurate knowledge of thermophysical properties of natural gas mixtures is of great importance for practical purposes for the gas industry from exploration stages to final customer use (Jaescke et al., 2002; Wagner & Kleinrahm, 2004; Gallagher, 2006). Two main properties are required by the oil and natural gas industry: i) phase equilibria and ii) pressure–density–temperature (PρT) data. The large impact of PρT, volumetric, data on production, processing and transportation of natural gas is well–known (Hall & Holste, 1990; Husain, 1993; Wagner & Kleinrahm, 2004; Bluvshtein, 2007). Although the required accuracy for the considered properties varies depending on the purpose for which they are used (performance analysis or gas sales, Mokhatab et al. 2006), a high degree of accuracy is frequently required for most of the applications. Thermophysical properties of natural gas systems must be accurately known for national and international custody transfer considering that flowmeters measurements applied to custody transfer are used to buy and sell natural gas between pipeline companies. It should be remarked that custody transfer operations relies on accurate density data, and thus, the economic impact of accurate density measurements is very large both for the calculation of energy content of natural gas and for the flow rate obtained usually from orifice meters (which are used in about 90 % of the metering stations). Inaccuracies in density measurements may lead to very large economical losses or profits for gas producers or buyers, and thus, to conflicts between companies or even countries. Therefore, density values with accuracies below the 0.1 % level are commonly required. Experimental measurements using highly accurate densimeters, mainly magnetic suspension type, may lead to measurements with a level of uncertainty below the required 0.1 % limit (Wagner & Kleinrahm, 2004; Patil et al., 2007), and thus, this would be the best and most reliable option to obtain density data. Nevertheless, natural gas composition changes remarkably from one to other reservoir (because of age and deep), moreover, the advances in exploration and drilling technologies allow to explore and produce from nonconventional reservoirs (such as ultra-deep ones) which characteristics are very different to the traditional ones. Therefore, considering that to carry out experimental density 18

measurements is very costly, both in time and economical resources, and having in mind that it is not possible to measure all possible compositions of natural gas in the wide pressure-temperature ranges required for production, transportation and processing purposes, thus, accurate predictive models are required by the gas industry. The current industrial standard model for custody transfer purposes is AGA8-DC92 (Starling & Savidge, 1992); this is a complex multiparametric equation of state which has been used for years as an international standard. Nevertheless, in spite of the common use of AGA8-DC92 in the gas industry it presents several problems. Any predictive model must be validated with reliable experimental data obtained on a limited number of samples that have well defined compositions. None of the samples used for the AGA8-DC92 validation are in the extended region, and thus, the application of this model for mixtures with large concentrations of CO 2 , N 2 , or long alkanes (as the natural gases from non-conventional reservoirs; Babusiaux, 2004) may be problematic. Moreover, the accuracy of AGA8-DC92 model decreases remarkably on going from region 1 (265 to 335 K, 0 to 12 MPa, deviations < 0.1 %) to regions 2 ( 12 to 17 MPa, 211 to 394 K, deviation < 0.3 %) and 3 (17 to 70 MPa, 144 to 477 K, deviation < 0.5 %) of model validity, and thus the model performance decreases remarkably when pressure and temperature increases. Recent studies have reported larger deviations than the claimed ones for AGA8-DC92 (Patil et al., 2007). Therefore, two main conclusions may be extracted: i) performance of AGA8-DC92 has to be systematically analyzed using a collection of carefully selected mixtures and ii) probably a new model has to be proposed as international standard for custody transfer and transmission purposes. Properties of natural gases, including PρT behavior, are a reflection of intermolecular forces rising by the simultaneous presence of very different types of molecules in these complex multicomponent fluids, and their evolution with pressure and temperature. Therefore macroscopic properties of these fluids, such as density, are a consequence of the microscopic structure and behavior of the involved molecules in the considered phases. These relationships between microscopic and macroscopic behavior may be developed through the principles of statistical mechanics if an accurate knowledge of the forces acting between the involved molecules were available. Hence, if we want to develop reliable and accurate models to predict natural gas density, we should get a deeper insight into the microscopic structure for complex gas mixtures and their relationships with PρT behavior. For this purpose, computational chemistry methods, mainly classical molecular dynamics and Monte Carlo approaches, are very useful tools allowing i) to infer microscopic structural and energetic features and ii) to predict macroscopic relevant properties, such as density, as a function of pressure and temperature. We report in this work a computational study in which the ability of classical molecular dynamics simulation methods to predict PρT behavior of complex natural gas mixtures is analyzed. Computational predictions are compared with available highly accurate experimental density data. The possibility of using this approach to predict density values with an acceptable degree of accuracy at moderate computational costs is analyzed, the weaknesses and strengths of the method together with possible future directions are considered. www.intechopen.com

Literature Review
A detailed analysis of the literature shows that the available studies on the use of molecular modelling with predictive purposes for thermophysical properties of natural gas like mixtures are scarce. We will not report here the available literature on the use of Gibbss ensemble Monte Carlo methods for the prediction of phase equilibria, both for pure compounds or for the involved binary mixtures, because this is not the object of this work. Moreover, most of the studies use Monte Carlo approach instead of the classical molecular dynamics methods proposed in this work. In this section, we will analyze the relevant studies available in the open literature in which computational methods, both classical molecular dynamics and Monte Carlo approaches, are used to predict thermophysical properties of natural gas mixtures. Results of literature analysis are reported in a chronological ordering. Saager and Fischer, 1989, reported a study on NVT molecular dynamics simulations of PVT and thermal properties of pure methane up to 1000 MPa. Duan et al., 1992, reported a wide study in which PVT properties of pure methane are predicted in the 273 -2000 K and 100 -20000 bar ranges using NVE molecular dynamics simulation methods together with a united atom approach. Results show deviations within the 1.5 % range. Yoshida & Uematsu, 1996, published a study reporting the results for the prediction of PVT properties of natural gases by Monte Carlo molecular simulation. They studied light mixtures composed by methane, ethane, propane, CO2 and N 2 in the conditions of natural gas transportation in pipelines. They reported deviations lower than 1.5 % in pressure, although density prediction analysis is not carried out in a straightforward manner. Duan et al., 1996, used molecular dynamics to simulate the PVT properties of the ternary mixture methane + CO 2 + N 2 , leading to results with low deviations compared with the general equation of state proposed by the authors. Neubauer et al., 1999, reported a study in which NPT Monte Carlo method was applied for the simulation of volumetric properties of natural gas mixtures, both in the single phase and two -phase conditions. A united atom approach was used leading to density deviations up to 5 % in the high pressure region, decreasing with decreasing pressure and increasing temperature. These too large deviations, that obviously hinder the application of the proposed method with purely predictive purposes, are a consequence of the poor representation of the composition of the gas during the simulation. Errington et al., 1999, reported a NPT Monte Carlo study on the properties of the n-alkanes homologous series, developing a united atom intermolecular potential and showing its validity for alkanes up to C78. Saturated densities are predicted within a 2 % limit and liquid densities for the longer n-alkanes to 1 %. Dysthe et al., 1999, used equilibrium molecular dynamics together with the Green-Kubo formalism to predict transport coefficients of multicomponent natural gas like mixtures including alkanes up to C4 , N 2 and He, both in the gas and liquid phases. Simulations were performed in the NVT ensemble with a united atom approach for alkanes leading to viscosity deviations of 7 % and 11 % for the gaseous and liquid states, respectively. Escobedo and Chen, 2001, developed a NPT Monte Carlo study for the prediction of Joule -Thomson inversion curves for several fluids, including pure methane and a gas condensate mixture (with alkanes up to C7), leading to reliable predictions. Lagache et al., 2001, reported a study in which NPT Monte Carlo simulations were performed to compute second order derivatives of the Gibbs energy for simple alkanes up to butane and for the methane -ethane binary mixture. The authors used a united atom potential. Results show that predicted data are in fair agreement with experimental values, even for complex properties such as Joule -Thomson coefficient for which deviations below the 10 % are obtained. Ungerer, 2003, reported a wide review in which the use of Monte Carlo and molecular dynamics methods is analyzed for several relevant fields in the petroleum and gas industry. Lagache et al., 2004, reported a study in which NPT Monte Carlo was used for the prediction of density and other relevant properties (thermal expansivity, isothermal compressibility, isobaric heat capacity and Joule -Thomson coefficient) of methane, ethane and two mixtures including heavy components up to 35 carbon atoms. The authors use a united atom approach to model the involved molecules. Density, and compressibility factor, show deviations up to 3 %, whereas for the remaining studied properties deviations are lower than 10 %. The authors show the importance of the mixtures characterization and representation to obtain accurate results. Ungerer et al., 2004, reported a NPT Monte Carlo study on the prediction of relevant properties, including density, for H 2 S -rich gases, showing that although reported results provide valuable information on the understanding of the complex mixed fluids, this computational approach does not lead to the high accuracy required for process design purposes for the studied acid gases. Ungerer et al.  reported a wide and useful study on the application of Monte Carlo methods for oil and gas production and processing purposes. They showed some of the results previously reported by Lagache et al., 2004, and claimed again to the remarkable importance of an adequate characterization of studied mixtures to obtain accurate results for single phase and equilibria properties. Bessieres et al., 2006, reported a study in which NPT Monte Carlo simulations were used to predict the Joule-Thomson inversion curve of pure methane using and united atom approach. Results show very accurate predictions with deviations below the 1 % limit. Vrabec et al., 2007, carried out a study on the performance of molecular simulation methods for the prediction of Joule-Thomson inversion curves for light natural gas mixtures. Reported results show deviations usually within the 5 % range, larger for high temperatures, but being competitive with most state-of-the-art EOS in predicting Joule -Thomson inversion curves. In a review work, Ungerer et al., 2007, analyzed the weaknesses and strengths of using molecular simulation for the predication of thermophysical properties of complex fluids, including natural gas like mixtures. The authors claim that one of the main limitations of the computational approach is the availability of potentials and force field parameters tested in wide pressure-temperature ranges. The main conclusions obtained from the analysis of the available open literature are: i) Monte Carlo approach is used in an exclusive basis when thermophysical properties of natural gas mixtures are under study. ii) United atom potentials are the most common option. iii) Studies in wide pressure-temperature ranges and for multicomponent mixtures are very scarce, and thus the performance of the proposed approaches is not clear.

www.intechopen.com
Thus, in this work we report results using the molecular dynamics approach, all-atoms potential, and analysis in wide P-T ranges for selected pure and mixed fluids. This methodology was considered because of the absence of similar studies in the open literature to analyze its validity for natural gas industry production, transportation and processing purposes.

Computational Methods
Classical molecular dynamics simulations were carried out using the TINKER molecular modeling package (Ponder, 2004). All simulations were performed in the NPT ensemble; the Nosé-Hoover method (Hoover, 1985) was used to control the temperature and pressure of the simulation system. The motion equations were solved using the Verlet Leapfrog integration algorithm (Allen & Tildesley, 1989). Long-range electrostatic interactions were treated with the smooth particle mesh Ewald method (Essmann, 1995). The simulated systems consist of cubic boxes, with the number of molecules and compositions, for mixed fluids, reported in Table 1, to which periodic boundary conditions were applied in the three directions to simulate an infinite system. The composition of the mixed fluid selected to test the performance of the computational approach for multicomponent natural gas-like mixtures resembles the one reported by Patil et al., 2007, for which very accurate and reliable density data obtained through magnetic suspension densitometers are reported. The number of molecules used for the simulation of pure compounds was selected to obtain systems with 4000 -5000 total atoms leading to reasonable computing times. For the mixture we have selected a total number of molecules (1000) that allow the representation of all the involved species, even those that appear at very low mole fraction but which effect on the mixed fluid behavior is important. The simulations were performed using a cutoff radius of L/2 Å for the non bonded interactions, L being the initial box side. Initial boxes generated using the PACKMOL program (Martínez & Martínez, 2005) were minimized according to the MINIMIZE program in TINKER package to a 0.01 kcal mol -1 Å -1 rms gradient. Long simulation times are needed for computing the properties of these fluids and procedures have to be designed carefully to avoid the presence of local minima. Therefore several heating and quenching steps in the NVT ensemble up to 600 K were performed after which a 100 ps NVT equilibration molecular dynamics simulation was run at the studied temperature; finally, from the output NVT simulation configuration, a run of 1 ns (time step 1 fs) in the NPT ensemble at the studied pressure and temperature was run, from which the first 0.5 ns were used to ensure equilibration (checked through constant energy) and the remaining 0.5 ns for data collection. n-Alkanes were described according to the so called Optimized Potential for Liquid Simulations (all atom version) OPLS-AA (Jorgensen et al., 1996), eqs. 1-4. Parameters for CO 2 and N 2 were obtained from the literature (Shi & Maginn, 2008;Lagache et al., 2005), with bonds, angles and non-bonded interactions treated using eqs. 1,2 and 4. A Lennard-Jones 6-12 potential, eq. 4, was used to describe the interactions between sites which are separated more than three bonds, if they are in the same molecule (intramolecular interactions), and for interactions between sites belonging to different molecules (intermolecular interactions). Non-bonded interactions between 1-4 sites are scaled with a 0.5 factor. Lorentz-Berthelot mixing rules are applied for Lennard-Jones terms between different sites, eqs. 5-6. The used forcefield parameters are reported in Table 2.
Thermophysical properties were obtained from the statistical analysis of the results in the last 0.5 ns together with fluctuation formulae for the NPT ensemble, eqs. 7-9. We have selected five relevant properties to analyze the performance of the proposed computational approach (Lagache et al., 2004): density (ρ, because of its importance for natural gas engineering purposes), thermal expansion coefficient (α P ), isothermal compressibility (β T ), isobaric heat capacity (C P ) and Joule -Thomson coefficient (μ, calculated according to eq. 10).
In eqs. 7-9, <> stands for ensemble averages, K for kinetic energy, P for potential energy, P for pressure, T for temperature, V for volume and k B for the Boltzmann constant. The notation δX (where X stands for any of the quantities reported in eqs. 7-9) denotes X-<X>. In eq.10, v stands for the molar volume.  Table 2. Forcefield parameters used along this work obtained from the literature (Jorgensen et al., 1996;Shi & Maginn, 2008;Lagache et al., 2004)

Methane
As a first proof of the ability of the computational approach proposed in this work to reproduce accurately the properties of pure methane (the main component of natural gas mixtures) we have calculated the phase equilibria diagram of this fluid. Although the objective of this work is not devoted to the analysis of phase equilibria, these results may be considered as a previous test of the model. Phase equilibria was calculated using Monte Carlo method from the coupled-decoupled (Martin & Siepmann, 1999), dual cut-off (Vlugt et al., 1998), configurational bias (Siepmann & Frenkel, 1992) simulations, using the MF-CPN strategy (Martin & Frischknecht, 2006), which were performed using MCCCS Towhee code in the NVT-Gibbs ensemble (Errington & Panagiotopoulos, 1998). Simulated phase equilibria in comparison with literature experimental results are reported in Fig. 1. Critical properties were calculated using the Ising scaling law (Smit and Williams, 1990) with a scaling exponent of 0.32. Results reported in Fig. 1 show a fair agreement between experimental reference data (NIST webbook, Linstrom & Mallard, 2009) and simulated values. Critical temperature is just a 0.23 % lower and critical density a 5.53 % larger than the experimental values. Therefore the proposed forcefield parameterization reproduces accurately the phase equilibria of pure methane. To obtain reliable density values from molecular dynamics simulations, the evolution of calculated densities with simulation time was analyzed (together with energy values) to assure that stable configurations are obtained and to check if longer simulation times are required. Therefore in Fig. 2 we report some representative results of these analyses for pure methane as a function of pressure and temperature, we report selected results for the lowest, medium and larger densities studied.  (Linstrom & Mallard, 2009) Results reported in Fig. 2 show that density values converge rapidly, and thus no longer simulation times are required (total simulation time is 1000 ps). In some cases, as it is reported in Fig. 2 for 450 K / 20 MPa simulations (and usually appears for the lower temperature simulations), density convergence is slower, nevertheless density also converge for this temperature, and thus statistical analysis was performed with the data once convergence is reached. Once density convergence is reached, fluctuations are in the ± 1.5 % range of the average density value. Calculated densities of pure methane are reported in Table 3 in comparison with available reference data. Simulations are carried out in wide pressure-temperature ranges (250 to 450 K and 20 to 100 MPa), relevant for natural gas engineering purposes, to analyze the global performance of the approach. These pressure-temperature ranges were used for all the simulations along this work. Deviations between experimental and simulated densities are in the 1 -2 % range for the studied pressure -temperature range, and thus, slightly larger than the accuracy required for many natural gas engineering purposes.  Table 3. Comparison between reference, ρ ref (Linstrom & Mallard, 2009), and NPT molecular dynamics simulated, ρ sim , density values, for methane as a function of pressure and temperature T / K P / MPa 10 3 α P / K -1 10 3 β T / MPa -1 C P / J mol -1 K -1 μ /  Table 4. Comparison between reference, ref (Linstrom &Mallard, 2009, andLemmon et al., 2002), and NPT molecular dynamics simulated properties, sim, for methane as a function of pressure and temperature. Properties: thermal expansion coefficient (α P ), isothermal compressibility (β T ), isobaric heat capacity (C P ) and Joule-Thomson coefficient (μ) The remaining properties obtained from molecular dynamics simulations are reported in Table 4. Average absolute deviation, between reference and simulated data, are 7.70 %, 9.19 %, 2.53 % and 37.62 % for thermal expansion coefficient, isothermal compressibility, isobaric heat capacity and Joule -Thomson coefficient, respectively. Therefore, deviations obtained for these properties obtained from fluctuation analysis of simulations are too large, especially for Joule -Thomson coefficient. We report in Fig. 3 a comparison between experimental and simulated values of isobaric heat capacity and Joule -Thomson coefficient, from the reported results it may be inferred that trends of the properties with pressure and temperature are properly reproduced by simulations, with the larger deviations obtained for Joule -Thomson coefficient at the lower pressures. Deviations obtained in this work are in agreement with those reported by Lagache et al. (2004) using Monte Carlo approach with a different force field parameterization. Therefore, the main problems rising from molecular dynamics simulations are for derived properties, especially for Joule -Thomson coefficient, as we may expect. Fig. 3. Isobaric heat capacity (CP) and Joule-Thomson coefficient (μ) of pure methane as a function of pressure and temperature. Symbols, NPT molecular dynamics simulated properties, and lines, reference data (Linstrom & Mallard, 2009). Color code: black (250 K), blue (350 K) and red (450 K)

Ethane
Density predictions for ethane are reported in Table 5, deviations between experimental and simulated data are larger for ethane than for methane (Tables 3 and 5), and for methane these deviations are larger than the ones required for natural gas engineering purposes. Thus inaccuracies of simulations increase with increasing chain length of n-alkane, in spite of using an all-atoms force field approach, this could lead to increasing inaccuracies for multicomponent mixtures rich in large n-alkanes. Derived properties are reported in Table  6, large deviations are obtained for all the studied properties, especially for Joule-Thomson coefficient, as we obtained in the previous section for methane. Average absolute deviation, between reference and simulated data, are 8.86 %, 10.62 %, 3.06 % and 17.798 % for thermal expansion coefficient, isothermal compressibility, isobaric heat capacity and Joule-Thomson coefficient, respectively. www.intechopen.com  Table 5. Comparison between reference, ρ ref (Linstrom & Mallard, 2009), and NPT molecular dynamics simulated, ρ sim , density values, for ethane as a function of pressure and temperature T / K P / MPa 10 3 α P / K -1 10 3 β T / MPa -1 C P / J mol -1 K -1 μ /  Table 6. Comparison between reference, ref (Linstrom &Mallard, 2009, andLemmon et al., 2002), and NPT molecular dynamics simulated properties, sim, for ethane as a function of pressure and temperature. Properties: thermal expansion coefficient (α P ), isothermal compressibility (β T ), isobaric heat capacity (C P ) and Joule-Thomson coefficient (μ)

CO 2 and N 2
Results for carbon dioxide and nitrogen are reported in Tables 7 -10. Although these gases are minor component of common natural gas mixtures, their effect on the thermodynamic properties of natural gas mixtures is relevant, and thus, the accurate reproduction of their properties through simulations is required to obtain reliable results. Density is simulated for both gases with low deviations, average absolute deviation, between reference and simulated data, are 1.70 % and 1.71 % for carbon dioxide and nitrogen respectively. Nevertheless, these deviations are larger than 1 %, and thus, as for the results reported for nalkanes in previous sections, simulations lead to deviations too large for engineering purposes. Derived properties reported in Tables 8 and 10, show too large deviations, leading to average absolute deviations, between reference and simulated data, of 9.52 and 9.11 % (for thermal expansion coefficient, for carbon dioxide and nitrogen, respectively), 13.09 and 13.17 % (for isothermal compressibility), 3.54 and 3.32 % (for isobaric heat capacity) and 70.57 and 40.3 % (for Joule-Thomson coefficient).  Table 7. Comparison between reference, ρ ref (Linstrom & Mallard, 2009), and NPT molecular dynamics simulated, ρ sim , density values, for carbon dioxide as a function of pressure and temperature T / K P / MPa 10 3 α P / K -1 10 3 β T / MPa -1 C P / J mol -1 K -1 μ /  Table 8. Comparison between reference, ref (Linstrom &Mallard, 2009, andLemmon et al., 2002), and NPT molecular dynamics simulated properties, sim, for carbon dioxide as a function of pressure and temperature. Properties: thermal expansion coefficient (α P ), isothermal compressibility (β T ), isobaric heat capacity (C P ) and Joule-Thomson coefficient (μ)  Table 9. Comparison between reference, ρ ref (Linstrom & Mallard, 2009), and NPT molecular dynamics simulated, ρ sim , density values, for nitrogen as a function of pressure and temperature T / K P / MPa 10 3 α P / K -1 10 3 β T / MPa -1 C P / J mol -1 K -1 μ /  Table 10. Comparison between reference, ref (Linstrom &Mallard, 2009, andLemmon et al., 2002), and NPT molecular dynamics simulated properties, sim, for nitrogen as a function of pressure and temperature. Properties: thermal expansion coefficient (α P ), isothermal compressibility (β T ), isobaric heat capacity (C P ) and Joule-Thomson coefficient (μ)

Lean Natural Gas Mixtures
Results for the multicomponent natural gas mixture studied in this work are reported in Table 11. Deviations are in the 2-4 % range with an average absolute deviation from experimental results of 2.55 %, which is too large in comparison with reference methods used by the natural gas industry for density predictions. There are not available experimental data for derived properties for the studied mixture. Therefore, although thermal expansion coefficient, isothermal compressibility, isobaric heat capacity and Joule-Thomson coefficient were calculated using molecular dynamics simulations for the studied mixture, these results are not reported here because it is not possible to infer their accuracy. Nevertheless, we may expect that considering that density deviations reported in Table 11 for the studied mixture are larger than those reported in previous sections for methane, ethane, carbon dioxide and nitrogen, deviations of derived properties for the mixture should be also larger than for the studied pure fluids.  Table 11. Comparison between reference, ρ ref (Patil et al., 2007), and NPT molecular dynamics simulated, ρ sim , density values, for the natural gas like mixture (composition reported in Table 1) as a function of pressure and temperature www.intechopen.com

Conclusions
Results reported in this work show preliminary conclusions obtained in the first stages of a wide computational study that we are carrying out on the performance of molecular simulation approaches for the prediction of thermophysical properties of multicomponent natural gas like mixtures. The main conclusions that may be inferred from these initial results may be resumed in: i) Density for the main components of natural gases is predicted with deviations in the 1 -2 % range for the 250 -450 K / 20 -100 MPa studied ranges. These deviations are low considering the purely predictive character of the considered approach but are too large for natural gas engineering purposes.

ii)
Derived properties (thermal expansion coefficient, isothermal compressibility, isobaric heat capacity and Joule-Thomson coefficient), for the main components of natural gases are predicted with deviations up to 10 % (even larger for Joule-Thomson coefficient for some fluids). Therefore, although these are the deviations commonly obtained using molecular dynamics approach for many complex fluids, they are too large to use this approach for industrial purposes. Fortunately, deviations for isobaric heat capacity, which is a property with remarkable importance for natural gas engineering purposes, is predicted with deviations usually lower than 5 %.

iii)
Predictions for the studied multicomponent natural gas like mixture lead to analogous results than those mentioned in conclusions i and ii: low deviations for density which lead to probably larger deviations for the remaining properties, but for all of them too large for industrial purposes. Moreover, it should be remarked that the studied mixture is composed of n-alkanes only up to C5, and thus, heavier mixtures containing larger alkanes should lead to even larger deviations. Therefore, considering the results reported in this work, molecular dynamics approach is not able to lead to predictions that may be used for natural gas production, processing or transportation purposes in the present situation. Nevertheless, we think that results reported in this work, and by other authors in the literature (Lagache et al., 2001;Lagache et al. 2004), are very promising, and the use of molecular simulation approach to predict complex natural gas like mixtures properties should not be discarded. Therefore studies are being carried out by the authors with the next main objectives: i) Improvement of forcefield parameterizations.

ii)
Comparison of the performance of several force fields.

iii)
Improvements of all atoms force fields approach in comparison with united atoms approach, does the increase of computational time lead to a remarkable improvement in the accuracy of predictions? iv) Decrease the errors in density predictions below the 1 % limit, in wide pressure -temperature ranges, and thus, leading to results competitive with the methods commonly used by the gas industry nowadays.

v)
Analyze the predictions for heavy mixtures containing long alkanes.

vi)
Comparison of the performance of Monte Carlo and molecular dynamics approaches.

vii)
Analyze from a molecular viewpoint the structural factors that govern these complex mixtures behavior using these computational tools.