Structural Elucidation of Aqueous Methane Solutions Based on Pair Distribution Functions

Yannick Agbor1, Mahammad Khalilov2, Usama Ahmed Khand3, Jestin B. Mandumpal4*

1Department of Petroleum Engineering, University of North Dakota, USA

2Department of Computer Science, Khazar University, Baku, Republic of Azerbaijan

3Department of Petroleum and Gas Engineering, Balochistan, University of Information Technology, Engineering, and Management Sciences, Quetta, Pakistan

4College of engineering and technology, American University of the Middle East, Kuwait

*Corresponding author:Jestin B. Mandumpal, College of engineering and technology, American University of the Middle East, Kuwait, E-mail:

Citation: Mandumpal JB, Aghbor Y, Khalilov MD, Ahmed Khand U (2020) Structural Elucidation of Aqueous Methane Solutions Based on Pair Distribution Functions. J Chem Sci Chem Engg 1(2):7-18.

Received Date: May 16, 2020;Accepted Date: May 21, 2020; Published Date: May 26, 2020


We report our findings on aqueous solution of methane (at very low concentrations) based on Molecular Dynamics simulations. Temperature plays a vital role regarding solute-solvent interactions (between methane and water molecules). Increasing temperature enhances the entropy driven interactions between methane molecules (methane clustering), in particular between 270 K and 290 K. Our analysis points into the fact that the water cages expand in order to accommodate more number of water molecules. Positional disorder of hydrogen atoms increases at higher temperatures. Estimation of the number of methane – methane contacts reveals that as the concentration of methane in water increases, the solute structure shows greater stability against variation of temperatures.

Keywords: Methane; water; aqueous solutions; molecular dynamics; computer simulation


Methane (the most abundant hydrocarbon in natural gas) − water mixture is one of the simplest prototypes exhibiting hydrophobic effect, the tendency of nonpolar molecules to aggregate in aqueous solution by excluding water molecules. Methane in water at low concentrations (generally known as foams) is a precursor to an interesting class of chemical architecture known as methane hydrates in which water molecules act as host materials surrounding the guest molecule (methane) at specific temperatures and pressures. Many researchers are riveted by the properties of methane hydrates due to the impact they make on petroleum industry and energy sectors: on one hand, the growth of methane hydrates posits a threat to swift transportation of oil and gas, and on the other hand, methane hydrates offer an alternative to the already decaying traditional energy resources, fossil fuels. A molecular level understanding of how methane and water molecules interact can be helpful for optimising the oil transportation by designing highly efficient oil plants, and methane gas recovery from underground [1]

Considerable efforts have been expended over the years for procuring a detail understanding of solute−solvent interactions in aqueous solutions, and hydrophobicity of solutes in water has been at the centre of debates over decades. Franks and Evans proposed their ice−berg model for explaining the hydrophobic effect, while Kauzmann introduced the concept of a hydrophobic bond. Glew suggested the formation of clathrate like structures in order to account for the solute−solvent interactions. For a brief discussion on these propositions, please see [2]. The structural changes in solutions are accompanied by variations in thermodynamic state functions such as entropy and free energy. The latter increases when a methane molecule is added to water (solvent), which spurs on structural changes in solutions, such as the formation of molecular cages built of water molecules around the solutes [3−4]. The increase in free energy is correlated with the changes in two other thermodynamic state variables, enthalpy and entropy [5].

One of the prominent hypotheses regarding this binary mixture is the structural changes that can be attributed to ice−berg model, an increased ordering of water molecules around the “hydrophobic” methane molecule [6]. However, the ordering does not limit to this traditional view of tetragonal−coordinated water molecules [5]: in an interesting computational study, Head Gordon has identified the existence of larger structures containing up to eleven water molecules around a methane molecule, which suggests that the structural changes upon the insertion of hydrophobic molecules be very complex [7].

In addition, many studies have explored the orientation mechanism of water molecules around nonpolar hydrophobic solutes. Tangential Hypothesis has been proposed decades ago [8]. By comparing the distance between the central solute molecule and the hydrogen and oxygen atoms of the surrounding water molecules, Rapaport et al. proposed that the orientation of oxygen− hydrogen bonds be tangential to a hypothetical hydration sphere. The Tangential Model has also been deduced straight forward from the radial distribution functions, according to Bratos et al [9] on the basis of the “fact” that the first peaks of solute−oxygenwater and solute−hydrogenwater pair distribution functions are exactly superimposed. A team led by Buckingham has also reached the same conclusion based on their Neutron Diffraction investigations [10]. In addition, Bratos et al note that water shells around different solutes are linked via hydrogen bonding, one of the dominant intermolecular interactions occurring in condensed matter [9]. In a following article, two of the authors proposed that the tangential orientation of water molecules to methane molecule is altered significantly upon variation of temperatures: more water molecules around the solute orient towards the bulk, precisely the hydrogen atoms are found to be facing towards bulk solvent at relatively high temperatures [11]. Tani et al have proposed another interesting scenario on the orientation of water molecules around hydrophobic solutes: beyond the first hydration shell, the water molecules orient to the solute molecules with both hydrogen atoms point radially inward [12]. On the contrary, while studying about the orientations of nonpolar solutes such as neon, Stillinger et al have identified that both water protons in the first hydration shell around the solute orient outward symmetrically [2]. In the second hydrations shell, they however observed unsymmetrical orientation of water protons.

It can be noticeable from the foregoing discussions that different studies have yielded profusion of hypothesis on the solution structure of methane−water systems. Thanks to the development in computing, liquid structure can nowadays be investigated microscopically for longer timescales which cannot even be imagined several decades ago [13]. Eying on resolving the structural conundrum of methane−water foams, we conducted a series of Molecular Dynamics computer simulations, extracted Radial Distribution functions of various atomic pairs, and report herein our analysis of how methane and water molecules orient to each other. This computational procedure has been proved to be one of the most efficient tools available to computational physicists to answer important questions related to liquid structure, in particular of systems wherein water acts as the principal component [14]. The purpose of this article is twofold: (i) to investigate the effect of temperature variation, primarily in the upper part of supercooling regime of water and beyond in its normal temperature domain (from 260 K to 300 K), on the structural features of the methane− water foams; (ii) to study the preferential orientation of water molecules around the hydrophobic solutes, in particular in the first hydration shell, using Molecular Dynamics simulations. In general, we are motivated by the internal dynamics of aqueous methane solutions promoted by temperature variations partially leading to the formation of interesting structures known as methane hydrates that have profound technological significance.


A series of Molecular Dynamics simulations using 4.6.5 version of GROMACS package was performed, and its important technical details are given in this section [15]. We employed TIP4P potential (Transferable Interaction Potential Four Point) for water and OPLS AA (Optimised Potential for Liquid Simulation All Atom) forcefield for methane [16−17]. It has been reported that TIP4P water model has been found to reproduce water structure and thermodynamic quantities such as solvation entropies and therefore its choice is justified [18]. Furthermore, while comparing with Advanced Light Source (ALS) data, the water structure has been qualitatively represented by TIP4P model albeit its being slightly less accurate than TIP5P model [19]. The variants of TIP4P model, such as TIP4P/ICE [20], TIP4P/2005 [21] and TIP4P/Ewald perform better in terms of the extraction of many physical properties of water. We chose, however, the standard model because of the fact that solution properties are greatly influenced by the solute model employed in the simulation, and OPLS models have been optimised with standard Transferable Interaction Potential (TIP) protocols (for example TIP4P or TIP5P). For a brief discussion regarding the water models, please see [22] for a lengthy discussion. First, a homogenous mixture of water and methane was constructed using the PACKMOL software which ensured better mixing of methane and water molecules [23]. The initial box size was determined, taking account of water’s density at each temperature. The total number of particles was fixed to 864 with the mole fraction of methane at 10 %. Followed by the energy minimisation of 5000 steps, the system was equilibrated first in NPT ensemble by keeping pressure at 1 atm for 2 ns and a further 500 pico second (ps) in NVT ensemble at each temperature. SHAKE algorithm was utilised for constraining the bonds and angles. Periodic Boundary conditions (PBC) were applied in order to mimic the system at infinite order. A 10 Å cut off was used for van der Waals and short-range electrostatic interactions. All long-range electrostatic interactions were dealt with Particle Mesh Ewald (PME) method. The simulations were carried out at five different temperatures, from 260 K to 300 K at 10 K interval. Nose−Hover thermostat was employed for efficient temperature control[24-25]. The length of the production run was determined by the criterion of the ‘least noise possible’ in the resulting pair distribution functions. Thus, the simulations were run for several nano seconds (100 ns – 200 ns), and the trajectory was saved in every 10 fem to second (fs) to ensure greater accuracy. The first 2 nano seconds of the production trajectory was not considered for the analysis. The radial distribution functions were obtained using a built−in GROMACS tool, gird. The trajectories were divided into 250 pico second data. These were ported to MICROSOFT EXCEL platform, and the average was computed using a built−in function from which the heights of the peaks in each of the radial distribution functions were accurately found out.

In order to investigate the variation of the number of methanemethane pairs across the temperature, a simple analysis program was written in Python. The program reads the simulation trajectory and performs the distance calculations between water molecules. The methane-methane distance is averaged over total number of molecules at each trajectory


Solute structure

In figure 1, the three inter pair distribution functions of methane have been shown, and we call them type I interactions since the same trends (an increase of the intensity (heights) of the peaks) upon hike in temperature can be observed in all three pairs, a signature of hydrophobic interactions (see Figure 1, a through c). The enhanced interactions lead to the aggregation of apolar species such as methane by themselves, which is central to many natural phenomena including protein folding [26]. It has to be noted that this trend is not generic in the case of carbonmethane – hydrogenmethane pair distribution functions (shown in figure 1 b). One can see that the temperature range at which the gradual increase in the height of the peaks observed gets narrower (from 270 K through 290 K rather than the whole temperature range simulated). It can also be noted that carbonmethane− hydrogenmethane RDFs show two distinct peaks at around 4 Å and 5 Å; besides there appears a broader peak at around 7Å, which points into the existence of solvent separated pair, as reported in the case of aqueous solutions of inert gases [2]. It has been suggested that the pair of methane molecules serve as a common edge for two neighbouring hydrations shells [2] and, at various distances between two methane pairs there exist well−structured several hydration layers [4,10].

Hydrogenmethane − hydrogenmethane pair distribution functions exhibit a clear trend as well: the peak intensities show a proportional increase with respect to temperature. However, we note that the sharpness of the peaks also increases as temperature decreases, which may be due to the fact that the number of methane configurations is naturally reduced upon drop in temperature. In general, correlations between hydrogen atoms in methane molecules are found to be weakest in hydrogenmethane− hydrogenmethane pair distribution functions.

Carbonmethane − carbonmethane pair distribution functions (shown in panel a in Figure 1) also indicate strong interactions between methane molecules as temperature increases, as a further evidence to hydrophobic effect [26]. We observe that methane− methane pairing does occur at all temperatures at which the simulations were carried out; this is expected considering the fact that methane is one of the simplest hydrophobic molecules with low solubility in water, leading to phase segregation.

Solution structure

In Figure 2, peak heights versus temperature plots, derived from four cross correlation functions between water and methane molecules (carbonmethane − oxygenwater, carbonmethane − hydrogenwater, hydrogenmethane − oxygenwater, and hydrogenmethane − hydrogenwater), have been shown. These four atom−atom pair distribution functions provide vital information regarding the solvation structure around the solute, in particular hydrophobic hydration (the structure of water in the vicinity of an apolar molecule). As can be seen from the figure 2, the heights of the all peaks are found to be decreasing as temperature increases, asserting that temperature has profound impact on methane− water interactions in foams, and we call it type II interactions. It has to be noted that this trend has also been observed for neon−water systems [11].

Solvent structure

The third type of interaction is a mixed type, and we call it Type III. As one can see in the figure 3, intra−pair distribution functions of water show two contrasting behaviours. First, all first peak heights gradually diminishes as temperature increases, which in particular is noticeable at lower temperatures.

In Figures 4 and 5, we present number of methane−methane pairs over a distance of 1 nanometres (nm) at five different temperatures, from 260 K to 300 K, at two different mole fractions (10 and 20%). Although the correlation of temperature and number of methane−methane contacts is not immediately noticeable, we can observe the rate of change of the number of methane− methane contacts beyond 0.5 nm as temperature changes in dilute solution (10 % aqueous solution). This is in contrast to 20 % aqueous solution. Furthermore, the solute structure is intact within the limit of 0.5 nm regardless of the concentration, and at higher concentration (20 %) the solute structure is more ordered (regardless of the temperature) for greater distance. It would be interesting to compare this result with methane−methane pair distributions obtained for 10 % aqueous solution explained earlier. The pair distributions clearly reveal that methane−methane pairing increases with respect to temperature. The estimation of the number of methane – methane contacts also shows similar trend, except for 260 K.

Figure 1: Type I interaction. In left panels shown are peak intensity versus temperature profiles of carbonmethane − carbonmethane, carbonmethane – hydrogen methane, and hydrogenmethane − hydrogenmethane interactions of methane molecules (from top to bottom). On right, corresponding pair distribution functions have been shown. In all except carbonmethane – hydrogenmethane pair distribution function, enhanced structural effects towards higher temperatures is observed throughout the whole temperature window, while the trend is observed within a shorter temperature domain in carbonmethane−hydrogenmethane. In panel b, the peaks corresponding to 290 K and 300 K are overlapped.

Figure 2: Type II interaction. Four cross correlation pair distribution functions between methane and water molecules: (a) carbonmethane − oxygenwater (b) carbonmethane – hydrogenwater (c) hydrogenmethane – oxygenwater (d) hydrogenmethane − hydrogenwater are shown. The heights of all peaks are found to be decreasing as temperature increases, asserting that temperature has profound impact on methane−water interactions in foams.

Figure 3: Type III interaction. From top to bottom, we show the pair distribution functions of water (solvent) molecules. (a) oxygenwater –oxygenwater (b) hydrogenwater – oxygenwater (c) hydrogenwater – hydrogenwater. The heights of first peaks and second peaks demonstrate contrasting behaviour.

Figure 4: Methane-Methane distances at various temperatures at mole fraction of 10% aqueous solution.

Figure 5: Methane-Methane distances at various temperatures at mole fraction 20%.


Although rendering of a perfect understanding of any physical phenomenon is cumbersome to achieve from a crude twodimensional statistical data such as radial (pair) distribution functions, they are still useful in providing vital information regarding the structure of liquids and aqueous solutions (in terms of solute−solvent orientations) at molecular level. The differences between the pair distribution function profiles cannot be detected by merely looking at the radial distribution functions, as the changes are only by an imperceptible degree across the whole temperature domain the simulations were carried out. This warrants accurate identification and estimation of peak heights in order to understand the trend clearly. We therefore plot the heights of major peaks (mainly the first ones) at corresponding temperatures as this refers to the solution structure directly. Based on the trends, we group these interactions into three.

One can see from Figure 1 that as the temperature decreases the intensity of methane−methane pairing decreases considerably (more than 33%). However, there appears conflicting reports regarding this; for example, Stilbs et al have reported that even at high methane concentration, methane−methane pairing does not occur [27]. They compared carbonmethane−water pair distribution functions with that of pure methane system. We believe that the simulation conditions (including the concentration and production run) adopted by Stilbs et al vary drastically with the present work, which may be one of the principal reasons of this discrepancy

Beyond 1 nano metre (nm), no subsequent peaks are observed unlike in the case of solvent−solvent radial distribution pairs, which will be discussed later in this article. This is in accordance with Potential of Mean Force (PMF) calculations performed by Suzuki et al [28]. This lack of long range order can greatly be affected by the fact that in our simulations, a virtual box containing 864 molecules has been employed with the mole fraction of methane at around 10%, and according to Koh et al, around 5% of the total water molecules can be found around a methane molecule [5]. This would mean that 16−23 water molecules can be found in the first hydration shells, according to various investigation reports [10, 29]. In addition, we note that the locations of first and second peaks are indeed much shorter than those of methane hydrates, for example sI hydrates, where in the nearest methane−methane distance, and the subsequent nearest neighbour distance are found to be 6.5 Å and 10.5 Å respectively [30]. The aggregation of methane molecules in water is still open to debate. Szalewicz et al argues that the degree of aggregation between methane molecules is very much limited in terms of number of molecules associated with it [31].

Summarising the information gathered from all three pair distribution functions ( shown in figure 1), we are able to discern that methane molecules prefer to aggregate themselves, as reported elsewhere [8, 26]. This aggregation is driven by entropy upon rise in temperature. The role of entropy is explained by the fact that clustering of nonpolar solutes reduces the number of water molecules to which the solute molecules directly interacts with [8]. The first peak in carbonmethane−carbonmethane pair distribution functions, observed around 3.9 Å, show this trend more clearly than the second, observed at around 7.5 Å. The second peak corresponds to the solvent separated methane pair as argued by Skipper [26].

The peaks we obtained for carbonmethane−oxygenwater and carbonmethane−hydrogenmethane are at around 3.5 Å, consistent with earlier and recent computer simulation findings [3, 32,33]. It can be seen that the generic shapes of these pair distribution functions (corresponding to carbonmethane − oxygenwater and carbonmethane−hydrogenwater correlation motions) exhibit similar features, in particular at higher distances. The sharper carbonmethane− oxygen water peaks suggest that more number of oxygen atoms become part of the skeleton of the water cages built around methane molecules. One of the notable differences with respect to earlier experimental findings, for example with those of Scheraga’s work, is that the peak of carbonmethane− hydrogenwater radial distribution function is found to be broader which prompted them to “confirm” the tangential orientation of oxygen−hydrogen bonds in water with respect to an imaginary sphere around methane molecules [8]. The tangential arrangement offers a distinct structural feature than the traditional view of tetrahedrally coordinated water molecules, forming undistorted “perfect” hydrogen bonds. It would be interesting at this juncture to note that the first peak in methane−oxygenwater radial distribution functions is found at 4.1Å for methane hydrates, as opposed to 3.5 Å in aqueous solutions [29]. This clearly points into the fact that the water cages by which methane hydrates are built of expand at appropriate pressure and temperature in order to accommodate more number of water molecules.

From the carbonmethane – oxygenWater pair distribution functions, Vega et al have found out that the peak height is found to be around 13 % lower at 600 K than at 298 Kelvin. On the contrary, our simulations record much sharper decrease (more than 35 %) over a smaller temperature change, from 260 K to 300 K. We also note that unlike carbonmethane – oxygenwater rdfs, the second peaks in the carbonmethane− hydrogenwater radial distribution functions are shifted towards shorter distances as temperature increases. It has also to be noted that carbonmethane−oxygenwater and carbonmethane−hydrogenwater rdfs bear resemblance to earlier works [3,29]. However, the radial distribution functions obtained in our simulations show several differences with respect to some other works in the literature in addition to the already mentioned in the previous paragraph. First, the carbonmethane−oxygenwater radial distribution functions exhibit a sharper first peak than carbonmethane−hydrogenwater rdfs. This can be explained by the fact that as temperature rises, hydrogen atoms being lighter than oxygen atoms take numerous configurations. In other words, the positional disorder of hydrogen atoms increases more due to higher diffusion rate of water molecules at higher temperatures. Second, the subsequent peaks are reported to be disrupted significantly in [3], in contrast to our data wherein we observe that the second peak is found to be perturbed only marginally.

Another important feature that can be deduced from our pair distribution functions is that hydrogenmethane −oxygenwater and hydrogenmethane − hydrogenwater rdfs appear to be broader (as was observed by Scheraga et al, based on their Monte Carlo simulations) compared to carbonmethane−oxygenwater and carbonmethane − hydrogenwater pair distribution functions [6]. Scheraga et al have also concluded that a significant proportion of water hydrogens are much closer to the centre of methane molecules than water oxygen atoms, however, we haven’t been able to identify this from our respective pair distribution functions. We rather find that carbonmethane – oxygenwater peaks are sharper than carbonmethane – hydrogenwater radial distribution peaks. Interestingly the second peak is more pronounced in the latter; this has also been observed by Scheraga et al: they attributed this trend due to enhanced hydrogen bonding between first and second shells [6].

Now, it would be interesting to compare our results with experimental (Neutron Diffraction) and coarse grain simulation findings. The distribution functions derived by Buckingham et al from their Neutron Diffraction studies regarding the nature of the peaks in general tend to agree with our results [10]. However, we note that they haven’t observed any subsequent peaks beyond the first hydration shell, indicating that the solution structure is short ranged, major implication of which they attribute to the generic behaviour of water−nonpolar solutes. Buckingham et al further note that associative behaviour of methane and water is prominent in molecular simulations, reflected in well−defined second peaks. On the other hand, comparing our results with methane−water pair distribution functions derived from a coarse grain simulation, we find that our simulations in principle agree with it as well [34]. In comparison to methane hydrates, water and methane molecules are more dispersed in methane−water foams as indicated by weaker peaks [35].

Using the same water model that we employed in the current simulations, Mancera et al have performed Monte Carlo simulations on water by varying another thermodynamic variable, pressure: from 1 atm through 1600 atm. They found that by increasing pressure the strength of carbon−oxygen interactions is diminished, and they have also observed the same trend with other popular model, SPC/E [36]. However, the disruption of the first peak under pressure has not been observed in carbonmethane−oxygenwater radial distribution function [37], and the authors explain that this is due to the fact that hydration sphere around the methane molecules is undisturbed [37].

Diminishing of first peak heights as shown in figure 3 is in accordance with the findings of Skipper [26] and also a previous computer simulation finding using TIP4P water model [38]. It has also been found out that as density increases the peak intensity decreases considerably [38]. This indicates that increasing pressure disrupts the solution structure (order) to a greater extent as expected because molecules compete with each other for the available space. On the contrary, an increase in peak (second) height in hydrogenwater − hydrogenwater pair distribution functions is observed, albeit slightly, with respect to a variation in temperature. The two contrasting behaviour indicates that while the water (solvent) structure remains intact, water molecules shuttle between first and second hydration shells around each water molecule.

As expected, the first peak heights of oxygenwater – oxygenwater radial distribution function decreases as temperature increases, and this effect has been observed by several researchers, for example Vega et al [3,38]. They have observed a strong dependence of cross correlation functions of water upon temperature, by simulating at 298 K and 600 K [39].From figure 3, it can be easily discern that the locations of the first peak (at 2.7 Å) has not been changed over the interval of 40 K (from 260 K through 300 K), as in the case of corresponding pairs in the other pair distribution function as shown below (in table 1), in analogy to other systems of practical interest such as cryosolvents [40]. According to Anderson, this indicates that tetrahedral nature of water around methane has been preserved, with on average 4.26 neighbours around each water molecule [41]. Koh et al has observed that the peak in oxygenwater−oxygenwater pair distribution function has not been changed upon the mixing of methane into water [5], and they concluded that the presence of methane molecule affects water structure only marginally in terms of shell compression and that there is no evidence of structural ordering due to methane molecules [5]. Furthermore, from figure 3, it can be seen that the profiles corresponding to second and third peaks show very little changes as well across the temperatures we simulated. By contrast, the oxygenwater−oxygenwater radial distribution functions of methane hydrates are found to be more structured [42].

Skipper et al, on the other hand, have observed that lower first minimum in the oxygenwater−oxygenwater pair distribution functions has been shifted to larger distances [43]. From our analysis outlined in the Figure 3, it is clear that the temperature does have some impact upon the water structure. The first peak of oxygenwater−hydrogenwater pair distribution function can be taken as the strength of hydrogen bonding [44], and the variation of temperature has made considerable impact upon the height of the peak, and thereby reducing the hydrogen bonding interaction between water molecules. On the other hand, it has to be noted that enhancement of pressure increases miscibility of methane and water [45].


Using radial distribution functions, we were able to classify intermolecular interactions into three: type I, II and III. Clustering in methane is found to be enhanced as temperature increases, based on the peak heights derived from various pair distribution functions. The contrasting behaviour of pair distribution functions of water (alias type III interaction) indicates that the combined solution structure (first and second hydration shell around a water molecule taken together) does not diminish and water molecules shuttle between these two layers to ensure the structural stability of solvent upon the variation in temperature. A definite orientational preference in the first hydration shell of water molecules around methane molecules can be inferred from our geometrical analysis. This study, we believe that, opens up further investigations and analytical propositions on the formation of clathrates structures wherein methane molecules is trapped by water molecules that are connected by hydrogen bonds. Further theoretical work including topological analysis such as Voronoi tessellation, Potential Mean Force (PMF) calculation, and computer simulations of higher resolution (for example ab initio, and QM/MM) may reveal the perplexing nature of aqueous methane-water solutions. On the other hand, the system can be understood better experimentally using sophisticated protocols such as Vibration−Rotation− Tunneling (VRT) spectroscopy.


JBM is grateful to the American University of the Middle East (AUM) for the infrastructure provided for the completion of the article. MK thanks to Prof. Hamlet Isaxanli and the dean of the School of Engineering and Applied Sciences, Dr. Hassan Niknafs for financial assistance for establishing Khazar Computational Research Center (KCRC), where a part of the work was completed.


  1. Errington JR, Boulougouris GC, Economou IG, Panagiotopoulos AZ and Theodorou DN. Molecular simulation of phase equilibria for water-methane and water-ethane mixtures. J Phys Chem B, 1988;102(44):8865-8873.
  2. Geiger A, Rahman A, Stillinger FH. Molecular Dynamics Study of the Hydration of Lennard-Jones Solutes. J Chem Phys. 1979;70:263-276.
  3. Hernandez-Cobos J, Mackie AD, Vega LF. The hyddrophobic hydration of methane as a function of temperature from histogram reweigting Monte Carlo simulations. J Chem Phys. (2001);114(17):7527-7535.
  4. Li J, Car R, Tang C and Wingreen NS. Hydrophobic interaction and hydrogen-bond network for a methane pair in liquid water, Proceedings of National Academy of Sciences. 2007;104(8):2626-2630.
  5. Buchanan P, Aldiwan N, Soper AK, Creek JL, Koh CA. Decreased structure on dissolving methane in water. Chem Phys Lett. 2005;415(1-3):89-93.
  6. Owicki JC, Scheraga HA. Monte Carlo calculations in the isothermal-isobaric ensemble. 2. Dilute aqueous solution of methane. J Amer Chem Soc. 1977;99(23):7413-7418.
  7. Head-Goron T. Is water structure around hydrophobic groups clathrate-like? Proceedings of National Academy of Sciences. 1995;92(18):8308-8312.
  8. Rapaport DC, Scheraga HA. Hydration of inert solutes. A molecular dynamic study. J Phys Chem. 1982;86(6):873-880.
  9. Guillot B, Guissani Y, Bratos S. A computer simulation study of hydrophoobic The J Chem Phys. 1991;95:3643-3648.
  10. De Jong PHK. Wilson JE, Neilson GW, Buckingham AD. Hydrophobic hydration of methane, Mol Phys. 1997;91:99-103.
  11. Guillot B, Guissani Y. A computer simulation study of the temperature dependence of the hydrophobic hydration. J Phys Chem. 1993;99:8075-8094.
  12. Alagona G, Tani A. Structure of a dilute aqueous solution of argon. A Monte Carlo simulation, The J Chem Phys. 1980;72:580-588.
  13. Walsh MR, Koh CA, Sloan ED, Sum AK and Wu DT. Microsecond simulations of spontaneous methane hydrate nucleation and growth, Science. 2009;326:1095-1098.
  14. Hirata F, Rossky PJ. A realization of ‘‘V structure’’ in liquid water, J Chem Phys. 1981:74;6867-6874.
  15. Abraham MA, Murtola T, Schultz R, Pall S, Smith JC, Hess B and Lindahl E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers Software. 2015:19-25.
  16. Jorgensen WL, Maxwell DS. Tirado-Rives, Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids J. J Amer Chem Soc. 1996;118(45):11225-11236.
  17. Tung Y, Chen L, Chen Y, Lin S. The Growth of Structure I Methane Hydrate from Molecular Dynamics Simulations, J Phys Chem. B 2010;114:10804-10813.
  18. Paschek D. Temperature dependence of the hydrophobic hydration and interaction of simple solutes: An examination of five popular water models, J Chem Phys. 2004;120:6674-6690.
  19. Sorenson JM, Hura G, Glaeser RM and Head-Goron T. What can x-ray scattering tell us about the radial distribution functions of water? J Chem Phys. 2000;113:9149-9161.
  20. Abascal LF, Sanz E, Fernández RG, Vega C. A potential model for the study of ices and amorphous water: TIP4P/Ice, J Chem Phys. 2005;122:234511.
  21. Abascal LF, Vega C. A general-purpose model for the condensed phases of water: TIP4P/2005, J Chem Phys, 2005;123:234505.
  22. Mandumpal JB. A Journey through water, a scientific exploration of the most anomalous liquid on earth. Bentham Science Dubai. 2017;288.
  23. Martínez L, Andrade R, Birgin EG, Martínez JM. PACKMOL: A package for building initial configurations for molecular dynamics simulations, J Comp Chem. 2009;30(13):2157-2164.
  24. Hoover WG. Canonical dynamics: Equilibrium phase-space distributions Phys Rev. 1985;31:1695-1697.
  25. Nose S. A united formulation of the constant temperature molecular dynamics methods, J Chem Phys. 1984;81:511-519.
  26. Skipper NT. Computer simulation of methane—water solutions. Evidence for a temperature-dependent hydrophobic attraction.  Chem Phys Lett. 1993;207:424-429.
  27. Laaksonen A, Stilbs P.  Molecular dynamics and NMR study of methane-water systems Mol Phys. 1991;74:747-764. 10.1080/00268979100102551

  29. Fukunishi Y, Tateishi T, Suzuki M. Octane/Water Interfacial Tension Calculation by Molecular Dynamics Simulation J Coll Inter Sci. 1996;180:188-192.
  30. Tse JS, Klein ML, McDonald IR. Molecular dynamics studies of ice Ic and the structure I clathrate hydrate of methane. J Phys Chem. 1983;87:4198-4203.
  31. Hawtin RW, Quigley D, Rodger PM. Gas hydrate nucleation and cage formation at a water/methane interface, Phys. Chem. Chem. Phys. 2008;10:4853-4864.
  32. Swaminathen S, Harrison SW, Beveridge DL. Monte Carlo studies on the structure of a dilute aqueous solution of methane, J. Amer. Chem. Soc. 1978;5705-5712.
  33. Akin-Ojo O, Szalewicz K. Does a pair of methane molecules aggregate in water? J Chem Phys. 2019;150:084501.
  34. Islam N, Flint M, Rick SW. Water hydrogen degrees of freedom and the hydrophobic effect, J Chem Phys. 2019;150:014502.
  35. Jacobsen LC, Molinero V. A Methane-Water model for coarse grained simulations of solutions and clathrate hydrates, J Phys Chem. B 2010;114:7302-7311.
  36. Forrisdahl O, Kvamme B, Haymet ADJ. Methane Clathrate hydrates: melting, super cooling, and phase separation from molecular dynamics computer simulations.  Mol Phys. 1996;89(3):819-834.
  37. Chau PL, Mancera RL. Computer simulation of the structural effect of pressure on the hydrophobic hydration of methane, Mol Phys. 1999;96:109-122.
  38. Koh CA, Wisbey RP, Wu X, Westacott RE, Soper AK. Water ordering around methane during hydrate formation, J Chem Phys. 2000;113:6390-6397.
  39. Sokol M, Dawid A, Dendzik Z, Gburski Z. Structure and dynamics of water—molecular dynamics study, J Mol Struct. 2004;704:341-345.
  40. Ohmori T, Kimura Y. In 14th International Conference on the Properties of Water and Steam in Kyoto, Kyoto.
  41. Swope WC, Anderson HC. A Molecular Dynamics Method for Calculating the Solubility of Gases in Liquids and the Hydrophobic Hydration of Inert-Gas Atoms in Aqueous Solution, J Phys Chem. 1984;88:6548-6556.
  42. Mandumpal JB. The molecular mechanism of solvent cryoprotection: Molecular Dynamics study of cryosolvents. LAP Lambert Academic Publishing: (2012);237.
  43. English NJ, Johson JK. Molecular-dynamics simulations of methane hydrate dissociation, The J Chem Phys. 2005;123:244503-1-12.
  44. Bridgeman CH, Buckingham AD, Skipper NT. Temperature dependence of solvent structure around a hydrophobic solute: a Monte Carlo study of methane in water, Chem Phys Lett. 1996;209-215.
  45. Okazaki S, Nakaishi K, Touhara H, Watanabe N, Adachi Y. A Monte Carlo study on the size dependence in hydrophobic hydration, J. Chem. Phys. 1981;74:5863-5871.
  46. Pruteanu CG, Marenduzzo D, Loveday JS. Pressure -induced miscibility increases of CH4 in H2O: A computational study using classical potentials, J Phys Chem. B 123, 2019;8091-8095,

Supplimentary Info

A mode of interaction between the solute (methane) and solvent molecules (we call Linked-Torri hypothesis) may be obtained by fixing a solute and a solvent molecule viz−a−viz in the three-dimensional space. For that, we consider the location of the first peaks of the four principal cross correlation functions between methane and water molecules, namely carbonmethane− oxygenwater(3.7Å),carbonmethane−hydrogenwater(3.5Å), hydrogenmethane−oxygenwater,(4.4Å),hydrogenmethane− hydrogenwater(4.2Å), and configured the distance between the two molecules based on the location of the first peak from the corresponding RDFs. First, we place carbon and oxygen, being the centres of the two molecules (methane and water) at a distance with respect to each other, which would fix (by default) all other distances. Based on the four aforementioned interatomic distances, we obtained the resulting structure, after considering various orientational possibilities, and we found that water molecules prefer to interact with methane with both of its hydrogen atom point towards the solute, as shown in the figure 6.

Figure 6: Linked-Tori hypothesis (A) A pair of methane and water molecules aligned along carbon−oxygen axis. (B) Water molecule is rotated 90 degree, (C) schematic diagram indicating the mechanism of a possible Linked Tori interaction: hydrogen−carbon−hydrogen of methane is considered as a part of a torus (green circle), while hydrogen−oxygen−hydrogen (blue circle) axis of water can be considered as a part of another torus. The interatomic distance has been approximated to two significant figures.

Water molecules orient perpendicularly with respect to the hydrogen−carbon−hydrogen axes (there are six of such axes in methane molecule). If methane and water molecules are considered as part of two imaginary circles (one completes hydrogen−carbon− hydrogen bonds and the other hydrogen−oxygen−hydrogen bonds), the local structure can be considered as locked tori (each circles linked to the other perpendicularly).

In this scenario, both hydrogen atoms face methane molecule while the oxygen atom point outward.