We have performed classical molecular dynamics simulations using the fully polarizable Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) forcefield implemented within the Tinker package to determine whether a more adequate treatment of electrostatics is sufficient to correctly describe the mixing of methane with water under high pressure conditions. We found a significant difference between the ability of AMOEBA and other classical, computationally cheaper forcefields, such as TIP3P, simple point charge–extended, TIP4P, and optimized potentials for liquid simulations–all atom. While the latter models fail to detect any effect of pressure on the miscibility of methane in water, AMOEBA qualitatively captures the experimental observation of the increased solubility of methane in water with pressure. At higher temperatures, the solubility of water in methane also increases; this seems to be associated with the breakdown of the fourfold hydrogen-bonded water network structure: bonding in water is weaker, so the energy cost of solution is lowered.
I. INTRODUCTION
Water is one of the most widespread molecules in our Solar System,1 being present in both the fluid states, gaseous and liquid, and solid states, as various crystalline phases or amorphous glasses. Apart from being the main condensed matter phase on the face of the Earth (71% of the area), water is also a major component of extraterrestrial environments, being found in large quantities on Uranus, Neptune, and Titan. Water’s ubiquity is only part of the reason for ongoing interest in the subject, another part being that it is a fundamental system at the core of physics, chemistry, and biology. Water is famous for its ability to form hydrogen-bonds,2 even in the liquid state, which lead to a vast range of peculiar behaviors, such as immiscibility with certain chemical species, while being an extremely potent solvent for others.
Methane is the smallest hydrocarbon, an almost spherical molecule noted particularly for its chemical non-reactivity. It is the simplest example of a hydrophobic oil,3 which is a chemical species that does not mix with water. As such, water–methane is a model system for oil-in-water systems and their properties, on which key biological processes such as protein binding, drug-docking, and membrane transport rest. Similarly to water, methane is also widespread both on Earth and throughout extraterrestrial environments, being another major component of celestial bodies such as Neptune and Titan.4
It has long been thought that the hydrophobicity is an on/off type of property, which is either present in a system or not.5 This idea has recently been refuted,6 as experimental measurements found that under compression, methane begins to dissolve in water, saturating at a maximum of 44 mol. % at a pressure of 2 GPa and a temperature of 400 K. By contrast, there is no evidence of water dissolving in methane.
Various attempts to reproduce this unexpected behavior using widely deployed computational models, such as SPC/E (Simple Point Charge–Extended),7 TIP3P, and TIP4P8 for H2O and OPLS-AA (Optimized Potentials for Liquid Simulations–All Atom)9 for CH4, have proved fruitless, with the models failing to produce any effect of increasing pressure on the miscibility of the two chemical species.10 The reasons for the suppression of methane’s hydrophobicity were later revealed using a computationally and experimentally intensive procedure involving a combination of long-time ab initio molecular dynamics (AIMD) simulations and neutron scattering measurements.11,12 This study showed that methane acquires an induced dipole moment as pressure is increased in the presence of water, thereby becoming less hydrophobic and hence incurring a smaller energy penalty when neighbored by water. Most recently, the observed effect has been reproduced for the first time using a new generation machine-learned many-body classical forcefield.13,14 The respective authors showed that, within their forcefield, methane acquires a dipole moment, similar to the AIMD simulations, and that the mixing process takes place on a 0.5 ns timescale.
Unfortunately, this kind of approach is computationally highly expensive and is not readily generalizable to other systems, as it requires generating the necessary training data and refitting the model for every single binary mixture of interest.
One of the main unanswered questions about water and methane high pressure mixtures is what components for a model are essential to allow the two to mix. Is the mixing driven by polarization, novel chemical bonding, or increased density in the mixture? Is it sufficient to have an adequate treatment of electrostatics in the system, or is it also necessary to have the correct compressibility? More practically, is a more physically complex forcefield [such as Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA)15] better than previously tested models at reproducing the unexpected experimentally observed increase in miscibility? Therefore, can AMOEBA be used as an in silico screening tool for other potential pressure-induced qualitative changes in binary systems? In the current study, we aim to answer the two later questions.
We have performed classical molecular dynamics simulations using the AMOEBA forcefield for both methane and water. AMOEBA is a fully polarizable forcefield, computing the polarization via self-consistent induced atomic dipoles in order to produce a description of electrostatic interactions based on polarizable atomic multipoles.15 This approach renders the forcefield noticeably more computationally expensive than standard models such as TIP3P (up to 8× slower in the initial implementation, on a box of 216 water molecules15). More recent implementations of Particle-Mesh Ewald (PME) for point multipoles reduced this to only double the cost of a given calculation16 using simple point charge models. We found that the forcefield was able to correctly reproduce the mixing of methane and water when the critical temperature of water is exceeded. More surprisingly, AMOEBA predicts an increase in the miscibility of water and methane at high pressures, in contrast to other general classical forcefields previously investigated, which were completely insensitive to pressure effects. Our current simulations qualitatively reproduce the experimentally observed behavior, despite showing significant shortcomings quantitatively. The most remarkable fact is that the current simulations seem to reproduce both the pressure-enhancement of solubility and its saturation at higher pressures. Our analysis leads us to believe that the suppression of hydrophobicity in water and methane under compression is due to a fine interplay between electronic effects and thermodynamic properties such as the relative compressibility of the two species.
II. THEORY OF SOLUBILITY
ΔUe and ΔVe are the change in energy and volume on solution, respectively, and J is a local coupling between water and methane. ΔSe is the entropy change excluding the entropy of mixing.
A high pressure phase transition (separation vs mixing) at T = 0 occurs when the enthalpy of solution changes sign (PT = ΔU/ΔV). Of course, T = 0 is below the melt curve; however, the metastable phase line extends to finite T and may cross the melt line before terminating at a critical point. Whether this happens in practice depends on the relative repulsion between the two species, represented by the enthalpy J.
A. Variation of regular solution entropy
This demixing phase transition occurs when the two-phase entropy is equal to the single phase entropy. It is second order since entropy is continuous with T but heat capacity is not. This is true independent of the Boltzmann approximation to the solubility limit.
B. Upper critical solution temperature
We note that any of these critical solution points could fall in the solid region and therefore not be observable.
C. Interaction energies
Under ambient conditions, the Bragg–Williams J term can be attributed to the loss of hydrogen bonds. An isolated water molecule dissolved in methane loses all hydrogen bonds, while methane substituting for water disrupts the hydrogen bonding network. Under pressure, packing efficiency becomes a factor. It is possible for methane to be incorporated into the open water network, which can reconstruct to maintain hydrogen bonds. If the gain in density offsets the energy cost, then pressure induces mixing. More likely is that only a finite fraction of Me can be incorporated while retaining a hydrogen-bonded network. This percolation threshold places an upper bound on the methane solubility, leading to a concentration-dependent J and non-regular solution.
At higher temperatures, the thermal energy is sufficient to disrupt the hydrogen-bonded network. This corresponds to supercritical water, as in most of our simulations. Now, we return to a more typical Bragg–Williams scenario, where the energy is mainly dipole–dipole interactions in pure water, preferred to weaker dipole–induced dipole interactions for water in methane. The absence of a H-bond network means that for supercritical water, we can expect a more symmetric concentration profile.
III. METHODS
The Tinker Molecular Dynamics package20,21 was used to perform simulations using the AMOEBA forcefield,22 which contains in-built parameterizations for both water and methane. A box containing 950 H2O molecules and 470 CH4 molecules (33 mol. % methane) was constructed with the molecules initially fully mixed. This box was then relaxed to get the molecules into a realistic starting position. The desired temperatures were achieved through 250 ps NVT equilibrations using the Bussi thermostat and RESPA integrator with a 2 fs time step. To achieve the desired pressures, 250 ps NPT equilibrations were performed using the Berendsen barostat. Finally, NPT production simulations were run, each totaling 5 ns, which should be more than sufficient to allow for equilibration in solubility. In these production runs, the RESPA integrator was deployed. For all these steps, a potential energy cutoff distance of 9 Å was used along with the vdW long-range correction for distances larger than 9 Å. A cutoff distance of 7 Å was used during Ewald summation.
The radial distribution functions were calculated using the RADIAL option from the built-in analysis package of Tinker.
IV. RESULTS
A. Supercritical mixing
It is well known that as water is heated to higher temperatures, it suffers a gradual loss of hydrogen bonds.24–26 At sufficiently high temperatures, there is no H-bond network and so (supercritical) water can act as a super-solvent.27,28 To establish a conceptually easily verifiable starting point, we have first monitored this behavior at elevated pressures of 1 and 2 GPa. The pressure difference in real water tends to aid in re-forming H-bonds, leading to slightly more bonds per water molecule with increasing pressure at a given temperature. For the current study, we are computing the number of hydrogen bonds per water molecule by integrating the O–Hwater pair distribution function up to the first non-zero minimum (occurring around ∼1.8 Å), excluding the two intra-molecular covalently bound H atoms.
Similarly to our previous study,10 we have monitored the mean square displacement, the individual radial distribution functions of the mixture, and have computed the entropy of mixing (Fig. 1), presented as a ratio of Eq. (6) to the regular solution value from Eq. (3). An important note to make at this point is that we are deliberately computing the methane–water (carbon–oxygen) radial distribution function as if the system is homogeneously mixed. In demixed cases, this measure should not oscillate around the value of 1 for all distances but instead show a trend that gradually increases to a value of 1 at distances larger than the size of the demixed region.
(Top) Number of H-bonds as a function of temperature at 1 and 2 GPa and the corresponding entropy of mixing. Decrease in the number of H-bonds as a function of temperature at 1 and 2 GPa and the corresponding entropy of mixing. (Middle) Evolution of the computed entropy of mixing with increasing temperature at 1 and 2 GPa. (Bottom) Snapshot of a simulation box deemed “fully mixed,” i.e., entropy of mixing ∼1. Hydrogen atoms are removed for clarity; oxygen (from water) is red, and carbon (from methane) is turquoise.
(Top) Number of H-bonds as a function of temperature at 1 and 2 GPa and the corresponding entropy of mixing. Decrease in the number of H-bonds as a function of temperature at 1 and 2 GPa and the corresponding entropy of mixing. (Middle) Evolution of the computed entropy of mixing with increasing temperature at 1 and 2 GPa. (Bottom) Snapshot of a simulation box deemed “fully mixed,” i.e., entropy of mixing ∼1. Hydrogen atoms are removed for clarity; oxygen (from water) is red, and carbon (from methane) is turquoise.
AMOEBA predicts a uniform and gradual evolution of the system with increasing temperature. The visualization of the simulation indicates obviously non-homogeneous systems at the lowest temperatures, and increasing T steadily moves the system toward more homogeneity. This behavior is already in contrast to that observed when using the SPC/E, TIP3P, or TIP4P models for water and the OPLS-AA model for methane. In that case, there was very little change up to temperatures around 600 K (623 K being the critical temperature of real water) and a jump in the mixing state of the system above that. This view is readily reinforced once we look at the actual number of H-bonds per water molecule, which shows an almost linear decrease with increasing temperatures between 400 and 800 K, reaching almost zero at the latter pressure, suggesting a H-bond free system and hence one where homogeneous mixing with methane should be present.
As temperature is increased, the entropy gradually climbs up to this value. This is again in stark contrast to the previously deployed models, which would show very similar behaviors for 500 and 600 K followed by a noticeable jump in entropy/miscibility at 800 K.
B. Pressure-induced miscibility at 400 K
Emboldened by the more physically plausible behavior of the model when water goes supercritical, we aimed to isolate again the effect of pressure alone on miscibility. For this, we looked at the 400 K isotherm for a methane–water system along with a pure methane and pure water system along the same thermodynamic path. The equations of state followed are depicted in Fig. 2.
Equations of state for the methane–water system at 400 K and high pressures as obtained using the AMOEBA forcefield. The discontinuity in pure water between 5 and 6 GPa is due to solidification.
Equations of state for the methane–water system at 400 K and high pressures as obtained using the AMOEBA forcefield. The discontinuity in pure water between 5 and 6 GPa is due to solidification.
Pressure dependence of the average entropy of mixing for isothermal compression at 400 K.
Pressure dependence of the average entropy of mixing for isothermal compression at 400 K.
The entropy of mixing (Fig. 3) as well as the visual inspection of the trajectories in Visual Molecular Dynamics (VMD) shows that pressure has the effect of increasing the miscibility of methane in water (Fig. 4). However, our 400 K simulations at 33 mol. % CH4 remain de-mixed, while the experiments reported a single phase at 40 mol. % and comparable conditions.
Snapshot of the simulation box at 400 K for the 0.5 GPa simulation (top) and 4 GPa simulation (bottom). The phase separation and the presence of methane in water are evident.
Snapshot of the simulation box at 400 K for the 0.5 GPa simulation (top) and 4 GPa simulation (bottom). The phase separation and the presence of methane in water are evident.
The non-homogeneous mixing is readily confirmed by the inspection of the Radial Distribution Functions (RDFs) of the system (Fig. 5). However, by following a similar line of thought to that in Sec. IV A, these also reinforce the idea that pressure does increase miscibility, despite never reaching the point of turning the entire ∼33 mol. % system into a single phase. The RDFs also show a stark contrast between the lower pressures considered, where the system was completely fluid, and the highest pressures where it solidified. The distributions for water are very similar to those in the pure system, while those of methane show more visible differences. Once solidified, the pure phases have virtually identical RDFs to the frozen mixture.
Radial distribution functions for all atom pairs in the 400 K system. The distributions were calculated across the entire system, assuming that it is homogeneous; the easily visible slopes indicate that the assumption is incorrect, and the system is actually phase-separated.
Radial distribution functions for all atom pairs in the 400 K system. The distributions were calculated across the entire system, assuming that it is homogeneous; the easily visible slopes indicate that the assumption is incorrect, and the system is actually phase-separated.
Figure 6 contrasts the mixture coordination numbers (the number of nearest neighbors of a given chemical species) with increasing pressure, with pure systems. In a single component system, one expects that the coordination number in the fluid would steadily increase with pressure.29–33 In our simulation, we notice an increase in the carbon–oxygen coordination number, consistent with increased miscibility, and a corresponding decrease in the O–O and C–C coordination numbers. This is expected given the fact that even in a phase-separated system, there is always an interface present, which reduces the O–O or C–C coordination number of the molecules present on the interface; however, the accelerated decrease as a function of increasing pressure is an indication of enhanced miscibility.
First shell coordination numbers in the pure systems and the mixture as a function of pressure at 400 K.
First shell coordination numbers in the pure systems and the mixture as a function of pressure at 400 K.
A more direct means of determining the solubility, possible for in silico studies, is by analyzing the density distribution of the two species in tandem, along the z-direction of the simulation box (given that the simulation box has a slab geometry). Illustrative plots at the lowest and highest fluid pressures are presented in Fig. 7. For a phase-separated system, in which there is only one interface present within the box, so having only two phases rich in each of the chemical species present, the maximum in the z-density corresponds to the species-rich phase, while the minimum corresponds to the amount of the species dissolved in the other species. By analyzing these together, and taking the ratio of the minimum at one species to the maximum of the other (which occur at the same distance along the z-axis by definition), one can extract directly the solubility. It is worth noting that this procedure has to be performed ensuring that the minima/maxima are relatively far away from the interface in order to be accurate and meaningful. The solubilities so-obtained are depicted as a function of pressure in the bottom panel of Fig. 7, for the 400 K isotherm as well as several other higher temperature isotherms (500, 600, and 800 K).
(Top and middle) Average density of oxygen and carbon atoms, respectively, along the z-axis (long axis) of the simulation box at 400 K, showing phase separation and the increased mixing of methane in water going from lower to higher pressure. The two curves correspond to the boxes depicted in Fig. 4. (Bottom) Solubility of CH4 computed from the z-densities for all our simulations. The lines are a guide to the eye and should not be interpreted as crossing each other (due to the 400 K isotherm freezing on further pressure increase).
(Top and middle) Average density of oxygen and carbon atoms, respectively, along the z-axis (long axis) of the simulation box at 400 K, showing phase separation and the increased mixing of methane in water going from lower to higher pressure. The two curves correspond to the boxes depicted in Fig. 4. (Bottom) Solubility of CH4 computed from the z-densities for all our simulations. The lines are a guide to the eye and should not be interpreted as crossing each other (due to the 400 K isotherm freezing on further pressure increase).
It appears that within the AMOEBA forcefield, at 400 K, there is a pressure-enhanced solubility of methane in water, while at higher temperatures, this effect is nullified and even reversed, leading to a decreasing solubility with increasing pressure. This highlights the opposite effects temperature and pressure have on enhancing solubility and implies that there are possible thermodynamic paths (i.e., closely following an isochore) through p–T space where the solubility of CH4 in H2O remains constant, by carefully balancing the effects of increasing pressure against those of increasing temperature. This is in agreement with experimental studies by Strässle et al.,34 who found that increasing pressure in tandem with increasing temperature tends to preserve the H-bond network in water, constituting a key limiting factor toward the solubility of methane. Equation (5) implies that dissolving methane in water increases the average density at 400 K but decreases it at 500 and 600 K. The reasoning put forward in Sec. IV A suggests that we can associate this with the idea of the methane being located in cavities in the subcritical hydrogen-bonded water network at 400 K, but that such cavities no longer exist at higher temperatures.
C. Binary methane–water phase diagram for the AMOEBA forcefield
Analyzing the z-density curves for all the pressures and temperatures considered in the present study, we can construct the binary phase diagram for high pressure methane–water mixtures calculated using the AMOEBA forcefield (Fig. 8).
Binary phase diagrams for fluid water–methane mixtures within the AMOEBA forcefield. The circles show concentration in the two-phase region, and the crosses simply show that the fully mixed case has 33% methane.
Binary phase diagrams for fluid water–methane mixtures within the AMOEBA forcefield. The circles show concentration in the two-phase region, and the crosses simply show that the fully mixed case has 33% methane.
The top panel of Fig. 8 shows that there is a pressure-dependent UCST in the methane–water system within the AMOEBA forcefield. This appears to be strongly pressure dependent at lower pressures (0–2 GPa) and relatively pressure-insensitive above that (3–7 GPa).
In addition to the UCST, there is a Lower Critical Solution Pressure (LCSP). Pressure has the opposite effect to temperature on H-bond formation in water (the former strengthens them, while the latter breaks them). However, for lower temperatures, where pressure-induced miscibility enhancement is observed, the situation is different, with the curve resembling the extremes of an UCST curve, cut short at the top by the appearance of the solid phase, and so displaying a miscibility gap in the pressure–concentration plot.
With the advent of laser-heating methods and current high pressure equipment, it is now timely to discuss these and their implication for the phase diagrams of binary systems.
V. DISCUSSION AND CONCLUSIONS
We have shown that the AMOEBA forcefield provides a model for methane–water mixtures across a range of temperatures and pressures. The UCST, at which full mixing is observed in the fluid, is above 500 K at ambient pressure and increases with pressure. This suggests that pressure tends to reduce the tendency toward mixing. This is opposite to experiments at ambient temperature where the solubility limit of methane in water increases with pressure. We rationalize this by noting that at these temperatures, AMOEBA water is in a supercritical state. It has a coordination number around 12, typical of a simple liquid, rather than the expected number of 4 under ambient conditions. Consequently, the anomalous behavior seen in experiment is not present, and the supercritical water–methane system behaves as a typical Bragg–Williams type solution.
We can therefore associate the subcritical anomalous behavior with the hydrogen-bonded network by considering the three contributions to the Gibbs free energy. The cohesive energy comes primarily from the hydrogen bonds, which are disrupted by mixing. The entropy of mixing is always positive and so favors mixing, especially at higher pressures. The molecule-number density of the mixtures is relatively unchanged by dissolving water in methane, because the high-coordination methane fluid does not offer convenient sites for water molecules, and the water is unable to form an H-bonding network. By contrast, the open structure of water contains cages, which can accommodate some methane molecules without significantly disrupting the H-bond network. By filling these cages, dissolved methane increases both the density and entropy. This works up to a concentration of methane where the H-bond network is disrupted and this enthalpy penalty drives demixing.
The difference in the observed behavior between AMOEBA and other potentials can be ascribed to the former’s inclusion of polarizability of the methane molecule. In the presence of strongly polar water, the electrostatic dipole/induced-dipole attraction lowers the energy of methane in water, partially compensating for the disruption of the H-bond network.
Finally, we find that the AMOEBA forcefield provides a more robust and qualitatively realistic description of methane–water high pressure mixtures than standard computational forcefields, such as SPC/E, TIP3P, TIP4P, and OPLS-AA. The picture depicted by AMOEBA is closer to the available experimental results on such systems, capturing the qualitative trends (pressure acting to enhance solubility at modest temperatures), despite still falling short of reaching quantitative agreement with either experimental data, ab initio quantum mechanical calculations (DFT), or more advanced purpose-built models such as MB-pol. In spite of these shortcomings, the present study indicates that AMOEBA may be used as a computationally cheap general screening tool for unexpected behaviors of interest, such as pressure-induced solubility enhancement, where the interplay between steric and electronic effects plays a key role.
SUPPLEMENTARY MATERIAL
See the supplementary material for full simulation results: radial distribution functions for supercritical runs at 500, 600, and 800 K; entropy of mixing and radial distribution functions for all atom pairs for isothermal compression runs at 400 K, and pure systems along the same p–T path; and the influence of long-range vdW corrections on the mixing state produced by AMOEBA for representative runs.
ACKNOWLEDGMENTS
This work was supported by the ERC Advanced Grant HECATE. M.K. was supported by a Summer Studentship granted by the University of Edinburgh. We acknowledge ARCHER2 for access to computational resources.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Matthew Kerr: Formal analysis (equal); Investigation (equal); Software (equal); Visualization (equal); Writing – original draft (equal). Graeme J. Ackland: Conceptualization (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Davide Marenduzzo: Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). Giovanni B. Brandani: Methodology (equal); Software (equal). Ciprian G. Pruteanu: Conceptualization (equal); Data curation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.