Molecular Dynamics Simulations of Methane Hydrate Using Polarizable Force Fields

Squeezing Oil into Water under Pressure: Inverting the Hydrophobic Effect

The molecular structure of dense homogeneous fluid water-methane mixtures has been determined for the first time using high-pressure neutron-scattering techniques at 1.7 and 2.2 GPa. A mixed state with a fully H-bonded water network is revealed. The hydration shell of the methane molecules is, however, revealed to be pressure-dependent with an increase in the water coordination between 1.7 and 2.2 GPa. In parallel, ab initio molecular dynamics simulations have been performed to provide insight into the microscopic mechanisms associated with the phenomenon of mixing. These calculations reproduce the observed phase change from phase separation to mixing with increasing pressure. The calculations also reproduce the experimentally observed structural properties. Unexpectedly, the simulations show mixing is accompanied by a subtle enhancement of the polarization of methane. Our results highlight the key role played by fine electronic effects on miscibility and the need to readjust our fundamental understanding of hydrophobicity to account for these.

Acoustic-propagation properties of methane clathrate hydrates from non-equilibrium molecular dynamics

Given methane hydrates' importance in marine sediments, as well as the widespread use of seabed acoustic-signaling methods in oil and gas exploration, the elastic characterization of these materials is particularly relevant. A greater understanding of the properties governing phonon, sound, and acoustic propagation would help to better classify methane-hydrate deposits, aiding in their discovery. Recently, we have published a new nonequilibrium molecular-dynamics (NEMD) methodology to recreate longitudinal and transverse perturbations, observing their propagation through a crystalline lattice by various metrics, to study the underlying S- and P-wave velocities (achieving excellent agreement with experiment) [Melgar et al., J. Phys. Chem. 122(5), 3006-3013 (2018); ibid.150, 084101 (2019)]. Here, we apply these NEMD methods to methane-clathrate systems to study acoustic-propagation characteristics, as well as the lattice elastic behavior. In so doing, we determine S- and P-wave velocities in excellent accord with experiment; we also ascertain the allowable magnitude range of acoustic perturbation and establish a threshold for lattice breakup and hydrate decomposition. Interestingly, upon dissociation, we observe the formation of methane nanobubbles, which agrees with previous studies on the microscopic fundamentals of hydrate dissociation by various means.

Non-equilibrium molecular-dynamics study of electromagnetic-field-induced propane-hydrate dissociation

Non-equilibrium molecular-dynamics simulations have been performed for dissolution of planar propane-hydrate/water interfaces in externally-applied electromagnetic (e/m) fields in the microwave to far infrared range (∼2.45-200 GHz) at electric-field intensities up to 2.0 V/nm and at roughly 20 K over/under temperatures vis-à-vis the zero-field propane-hydrate melting point. Upon e/m-field application, there is a field-frequency threshold above which the dissociation rate drops significantly, with a plateau therein for larger-frequencies. It was found that higher intensity and lower frequency facilitates dissociation. Except in the presence of a thermal driving-force, the 10 GHz frequency shows more substantial rate-enhancement effect vis-à-vis static electric fields or, indeed, lower-frequency e/m fields.

Molecular-dynamics study of propane-hydrate dissociation: Fluctuation-dissipation and non-equilibrium analysis

Equilibrium and non-equilibrium molecular-dynamics (MD) simulations have been performed to investigate thermal-driven break-up of planar propane-hydrate interfaces in contact with liquid water over the 260-320 K range. Two types of hydrate-surface water-lattice molecular termination were adopted, at the hydrate edge with water, for comparison: a 001-direct surface cleavage and one with completed cages. Statistically significant differences in melting temperatures and initial break-up rates were observed between both interface types. Dissociation rates were observed to be strongly dependent on temperature, with higher rates at larger over-temperatures vis-à-vis melting. A simple coupled mass and heat transfer model, developed previously, was applied to fit the observed dissociation profiles, and this helps us to identify clearly two distinct hydrate-decomposition régimes; following a highly temperature-dependent break-up phase, a second well-defined stage is essentially independent of temperature, in which the remaining nanoscale, de facto two-dimensional system's lattice framework is intrinsically unstable. Further equilibrium MD-analysis of the two-phase systems at their melting point, with consideration of the relaxation times gleaned from the auto-correlation functions of fluctuations in a number of enclathrated guest molecules, led to statistically significant differences between the two surface-termination cases; a consistent correlation emerged in both cases between the underlying, non-equilibrium, thermal-driven dissociation rates sampled directly from melting with that from an equilibrium-MD fluctuation-dissipation approach.

Formation of methane clathrates in carbon nanotubes: a molecular dynamics study

In this study, we investigated the possible formation of methane clathrates in CNTs with different sizes.

A modeling study of methane hydrate decomposition in contact with the external surface of zeolites

Methane hydrate dissociates on the external surface of siliceous zeolites with methane absorbed by the solid and water forming a liquid-like phase.

Communication: Librational dynamics in water, sI and sII clathrate hydrates, and ice Ih: Molecular-dynamics insights

Equilibrium molecular-dynamics simulations have been performed for liquid water, and on metastable sI and sII polymorphs of empty hydrate lattices, in addition to ice Ih, in order to study the dynamical properties of librational motion (rotation oscillation) depicted by protons in water molecules. In particular, hydrate lattices were found to display prominent "bifurcated" features, or peaks, at circa 70 and 80-95 meV (or ∼560 and 640-760 cm(-1), respectively), also displayed by ice, in essentially quantitative agreement with experimental neutron-scattering data. However, observed differences in dispersion between these librational modes between these two structures (both hydrate polymorphs vis-à-vis ice), owing primarily to density effects, have been decomposed into contributions arising from angular-velocity dynamics about axes in the local molecular frame of water molecules, with in-plane "wagging" and "twisting" rationalising one mode at ∼70 meV, and out-of-plane motion for the higher-frequency band. This was confirmed explicitly by a type of de facto normal-mode analysis, in which only immediate layers of water molecules about the one under consideration were allowed to move. In contrast, liquid water displayed no marked preference for such local in- or out-of-plane modes characterising librational motion, owing to the marked absence of rigid, pentamers or hexamers therein.

Effects of thermodynamic inhibitors on the dissociation of methane hydrate: a molecular dynamics study

We investigate the effects of methanol and NaCl, which are known as thermodynamic hydrate inhibitors, on the dissociation kinetics of methane hydrate in aqueous solutions by using molecular dynamics simulations. It is shown that the dissociation rate is not constant but changes with time. The dissociation rate in the initial stage is increased by methanol whereas it is decreased by NaCl. This difference arises from the opposite effects of the two thermodynamic inhibitors on the hydration free energy of methane. The dissociation rate of methane hydrate is increased by the formation of methane bubbles in the aqueous phase because the bubbles absorb surrounding methane molecules. It is found that both methanol and NaCl facilitate the bubble formation. However, their mechanisms are completely different from each other. The presence of ions enhances the hydrophobic interactions between methane molecules. In addition, the ions in the solution cause a highly non-uniform distribution of dissolved methane molecules. These two effects result in the easy formation of bubbles in the NaCl solution. In contrast, methanol assists the bubble formation because of its amphiphilic character.

Energy benchmarks for methane-water systems from quantum Monte Carlo and second-order Møller-Plesset calculations

The quantum Monte Carlo (QMC) technique is used to generate accurate energy benchmarks for methane-water clusters containing a single methane monomer and up to 20 water monomers. The benchmarks for each type of cluster are computed for a set of geometries drawn from molecular dynamics simulations. The accuracy of QMC is expected to be comparable with that of coupled-cluster calculations, and this is confirmed by comparisons for the CH4-H2O dimer. The benchmarks are used to assess the accuracy of the second-order Møller-Plesset (MP2) approximation close to the complete basis-set limit. A recently developed embedded many-body technique is shown to give an efficient procedure for computing basis-set converged MP2 energies for the large clusters. It is found that MP2 values for the methane binding energies and the cohesive energies of the water clusters without methane are in close agreement with the QMC benchmarks, but the agreement is aided by partial cancelation between 2-body and beyond-2-body errors of MP2. The embedding approach allows MP2 to be applied without loss of accuracy to the methane hydrate crystal, and it is shown that the resulting methane binding energy and the cohesive energy of the water lattice agree almost exactly with recently reported QMC values.

The Use of Anisotropic Potentials in Modeling Water and Free Energies of Hydration

We propose a novel, anisotropic rigid-body intermolecular potential model to predict the properties of water and the hydration free energies of neutral organic solutes. The electrostatic interactions of water and the solutes are modeled using atomic multipole moments up to hexadecapole; these are obtained from distributed multipole analysis of the quantum mechanically computed charge densities and include average polarization effects in solution. The repulsion-dispersion water-water interactions are modeled with a three-site, exp-6 model fitted to the experimental liquid water density and oxygen-oxygen radial distribution function at ambient conditions. The proposed water model reproduces well several water properties not used in its parametrization, including vapor-liquid coexistence densities, the maximum in liquid water density at atmospheric pressure, the structure of ordered ice polymorphs, and the liquid water heat capacity. The model is used to compute the hydration free energy of 10 neutral organic solutes using explicit-solvent free energy perturbation. The solute-solute repulsion-dispersion intermolecular potential is obtained from previous parametrizations on organic crystal structures. In order to calculate the free energies of hydration, water-solute repulsion-dispersion interactions are modeled using Lorenz-Berthelot combining rules. The root-mean-square error of the predicted hydration free energies is 1.5 kcal mol(-1), which is comparable to the error found using a continuum mean-field quantum mechanical approach parametrized using experimental free energy of hydration data. The results are also contrasted with explicit-solvent hydration free energies obtained with an atomic charge representation of the solute's charge density computed at the same level of theory used to compute the distributed multipoles. Replacing the multipole description of the solute's charge density with an atomic charge model changes the free energy of hydration by as much as 3 kcal mol(-1) and provides an estimate for the effect of the modeling quality of the intermolecular electrostatic forces in free energy of solvation calculations.

Polarizable Multipole-Based Force Field for Dimethyl and Trimethyl Phosphate

Phosphate groups are commonly observed in biomolecules such as nucleic acids and lipids. Due to their highly charged and polarizable nature, modeling these compounds with classical force fields is challenging. Using quantum mechanical studies and liquid-phase simulations, the AMOEBA force field for dimethyl phosphate (DMP) ion and trimethyl phosphate (TMP) has been developed. On the basis of ab initio calculations, it was found that ion binding and the solution environment significantly impact both the molecular geometry and the energy differences between conformations. Atomic multipole moments are derived from MP2/cc-pVQZ calculations of methyl phosphates at several conformations with their chemical environments taken into account. Many-body polarization is handled via a Thole-style induction model using distributed atomic polarizabilities. van der Waals parameters of phosphate and oxygen atoms are determined by fitting to the quantum mechanical interaction energy curves for water with DMP or TMP. Additional stretch-torsion and angle-torsion coupling terms were introduced in order to capture asymmetry in P–O bond lengths and angles due to the generalized anomeric effect. The resulting force field for DMP and TMP is able to accurately describe both the molecular structure and conformational energy surface, including bond and angle variations with conformation, as well as interaction of both species with water and metal ions. The force field was further validated for TMP in the condensed phase by computing hydration free energy, liquid density, and heat of vaporization. The polarization behavior between liquid TMP and TMP in water is drastically different.

Benchmarking the performance of density functional theory and point charge force fields in their description of sI methane hydrate against diffusion Monte Carlo

High quality reference data from diffusion Monte Carlo calculations are presented for bulk sI methane hydrate, a complex crystal exhibiting both hydrogen-bond and dispersion dominated interactions. The performance of some commonly used exchange-correlation functionals and all-atom point charge force fields is evaluated. Our results show that none of the exchange-correlation functionals tested are sufficient to describe both the energetics and the structure of methane hydrate accurately, while the point charge force fields perform badly in their description of the cohesive energy but fair well for the dissociation energetics. By comparing to ice Ih, we show that a good prediction of the volume and cohesive energies for the hydrate relies primarily on an accurate description of the hydrogen bonded water framework, but that to correctly predict stability of the hydrate with respect to dissociation to ice Ih and methane gas, accuracy in the water-methane interaction is also required. Our results highlight the difficulty that density functional theory faces in describing both the hydrogen bonded water framework and the dispersion bound methane.

Modeling organochlorine compounds and the σ-hole effect using a polarizable multipole force field

The charge distribution of halogen atoms on organochlorine compounds can be highly anisotropic and even display a so-called σ-hole, which leads to strong halogen bonds with electron donors. In this paper, we have systematically investigated a series of chloromethanes with one to four chloro substituents using a polarizable multipole-based molecular mechanics model. The atomic multipoles accurately reproduced the ab initio electrostatic potential around chloromethanes, including CCl₄, which has a prominent σ-hole on the Cl atom. The van der Waals parameters for Cl were fitted to the experimental density and heat of vaporization. The calculated hydration free energy, solvent reaction fields, and interaction energies of several homo- and heterodimer of chloromethanes are in good agreement with experimental and ab initio data. This study suggests that sophisticated electrostatic models, such as polarizable atomic multipoles, are needed for accurate description of electrostatics in organochlorine compounds and halogen bonds, although further improvement is necessary for better transferability.

A Monte Carlo simulation study of methane clathrate hydrates confined in slit-shaped pores

Monte Carlo simulations are used to study the structure, stability, and dissociation mechanisms of methane hydrate crystals inside carbon-like slit-shaped pores. The simulation conditions used mimic experimental studies of the dissociation of methane and propane hydrates in mesoporous silica gels (Handa, Y. P.; Stupin, D. J. Phys. Chem. 1992, 96, 8599). Simulations are performed under conditions of fixed methane pressure and fixed water loading, with the temperature increased in steps, with long equilibrations at each temperature. The initial structures of the confined hydrates are taken to be bulk-like, and pore widths chosen to accommodate integer or half-integer numbers of hydrate unit cells. Density profiles and orientational order parameter profiles are obtained and used to understand the structural changes associated with hydrate dissociation. Three different common water models, SPC/E, TIP4P, and TIP4P/2005, are used and the results compared. For water modeled using either the TIP4P or TIP4P/2005 potentials, dissociation temperatures are depressed proportionally to the inverse pore width, as predicted by the macroscopic Gibbs-Thomson equation. This behavior is observed for pores small enough that only half-cages of the clathrate structure are present. Experimental work has verified Gibbs-Thomson behavior for pores as small as 2 nm (Seshadri, K.; Wilder, J. W.; Smith, D. H. J. Phys. Chem. B 2001, 105, 2627); micropores of the size studied here have not yet been studied by experiment. Interestingly, the dissociation of hydrates modeled using the SPC/E water potential does not display the predicted pore-size dependence, and the dissociation mechanisms in this model seem to be quite different than those in the TIP4P-type models. In the SPC/E hydrates, with increasing temperature, cage dissocation occurs before methane desorption. In TIP4P-type hydrates, these processes occur either at the same temperature (to within the resolution of this study) or with dissociation occurring at higher temperatures than desorption. These simulations show that a variety of interesting clathrate structures and phase behaviors may be accessed in suitably designed microporous materials, with potentially useful applications in gas storage or separations.

The mobility of water molecules through gas hydrates

Crystal growth simulations of gas hydrates have suggested that hydrate cages may occasionally be occupied by H(2)O rather than guest molecules, leaving interstitial defects within the hydrate crystal. Further inspection of the behavior of these interstitial H(2)O molecules has revealed that they are relatively highly mobile entities within a gas hydrate. In this paper, we report these observations and examine the molecular mechanisms responsible for the transport of these interstitial molecules through hydrate crystals. Four distinct pathways for the H(2)O molecule transport between cages are found, each facilitated by the presence of empty cages. The relative richness of the observed behavior of interstitial defects suggests that interstitial diffusion could be an important mechanism for the mass transport of H(2)O molecules through gas hydrates.

Multipole electrostatics in hydration free energy calculations

Hydration free energy (HFE) is generally used for evaluating molecular solubility, which is an important property for pharmaceutical and chemical engineering processes. Accurately predicting HFE is also recognized as one fundamental capability of molecular mechanics force field. Here, we present a systematic investigation on HFE calculations with AMOEBA polarizable force field at various parameterization and simulation conditions. The HFEs of seven small organic molecules have been obtained alchemically using the Bennett Acceptance Ratio method. We have compared two approaches to derive the atomic multipoles from quantum mechanical calculations: one directly from the new distributed multipole analysis and the other involving fitting to the electrostatic potential around the molecules. Wave functions solved at the MP2 level with four basis sets (6-311G*, 6-311++G(2d,2p), cc-pVTZ, and aug-cc-pVTZ) are used to derive the atomic multipoles. HFEs from all four basis sets show a reasonable agreement with experimental data (root mean square error 0.63 kcal/mol for aug-cc-pVTZ). We conclude that aug-cc-pVTZ gives the best performance when used with AMOEBA, and 6-311++G(2d,2p) is comparable but more efficient for larger systems. The results suggest that the inclusion of diffuse basis functions is important for capturing intermolecular interactions. The effect of long-range correction to van der Waals interaction on the hydration free energies is about 0.1 kcal/mol when the cutoff is 12Å, and increases linearly with the number of atoms in the solute/ligand. In addition, we also discussed the results from a hybrid approach that combines polarizable solute with fixed-charge water in the HFE calculation.

Gaussian Multipole Model (GMM)

An electrostatic model based on charge density is proposed as a model for future force fields. The model is composed of a nucleus and a single Slater-type contracted Gaussian multipole charge density on each atom. The Gaussian multipoles are fit to the electrostatic potential (ESP) calculated at the B3LYP/6-31G* and HF/aug-cc-pVTZ levels of theory and tested by comparing electrostatic dimer energies, inter-molecular density overlap integrals, and permanent molecular multipole moments with their respective ab initio values. For the case of water, the atomic Gaussian multipole moments Q(lm) are shown to be a smooth function of internal geometry (bond length and bond angle), which can be approximated by a truncated linear Taylor series. In addition, results are given when the Gaussian multipole charge density is applied to a model for exchange-repulsion energy based on the inter-molecular density overlap.

Clathrate hydrates with hydrogen-bonding guests

Clathrate hydrates (CHs) are inclusion compounds in which "tetrahedrally" bonded H(2)O forms a crystalline host lattice composed of a periodic array of cages. The structure is stabilized by guest particles which occupy the cages and interact with cage walls via van der Waals interactions. A host of atoms or small molecules can act as guests; here the focus is on guests that are capable of strong to intermediate H-bonding to water (small ethers, H(2)S, etc.) but nevertheless "choose" this hydrate crystal form in which H-bonding is absent from the equilibrium crystal structure. These CHs can form by exposure of ice to guest molecules at temperatures as low as 100-150 K, at the (low) guest saturation pressure. This is in contrast to the "normal" CHs whose formation typically requires temperatures well above 200 K and at least moderate pressures. The experimental part of this study addresses formation kinetics of CHs with H-bonding guests, as well as transformation kinetics between different CH forms, studied by CH infrared spectroscopy. The accompanying computational study suggests that the unique properties of this family of CHs are due to exceptional richness of the host lattice in point defects, caused by defect stabilization by H-bonding of water to the guests.

A multipole-based water potential with implicit polarization for biomolecular simulations

A new water potential, DMIP (distributed multipoles, implicit polarization), is constructed using distributed multipoles to describe the electrostatic interactions, while accounting for polarization implicitly. In this procedure, small clusters are randomly sampled from atomistic simulations of bulk water using the AMOEBA (Ren and Ponder, J Comput Chem 2002, 23, 1497) potential. The multipole moments of the central water in each cluster are obtained from ab initio densities for each cluster, and the moments are then averaged over all clusters. Properties of bulk water calculated using DMIP compare favorably with existing data from AMOEBA simulations and experiment, with a conservative estimate of reduction in compute time of roughly 40%. The implicit force-field is also shown to work compatibly with existing polarizable multipole-based force-fields for biomolecules.

Calculation of protein-ligand binding free energy by using a polarizable potential

The binding of charged ligands benzamidine and diazamidine to trypsin was investigated by using a polarizable potential energy function and explicit-water molecular dynamics simulations. The binding free energies were computed from the difference between the free energies of decoupling the ligand from water and protein environments. Both the absolute and the relative free energies from the perturbation simulations agree with experimental measurements to within 0.5 kcal.mol(-1). Comparison of free-energy components sampled from different thermodynamic paths indicates that electrostatics is the main driving force behind benzamidine recognition of trypsin. The contribution of electronic polarization to binding appears to be crucial. By computing the free-energy contribution caused by the polarization between the ligand and its surroundings, we found that polarization has the opposite effect in dissimilar environments. Although polarization favors ligand solvation in water, it weakens the protein-ligand attraction by screening the electrostatic interaction between trypsin and benzamidine. We also examined the relative binding free energies of a benzamidine analog diazamidine to trypsin. The changes in free energy on benzamidine-diazamidine substitution were tens of kilocalories in both water and trypsin environments; however, the change in the total binding free energy is <2 kcal.mol(-1) because of cancellation, consistent with the experimental results. Overall, our results suggest that the use of a polarizable force field, given adequate sampling, is capable of achieving chemical accuracy in molecular simulations of protein-ligand recognition.

Heterogeneous crystal growth of methane hydrate on its sII [001] crystallographic face

This paper presents a systematic molecular simulation study of the heterogeneous crystal growth of methane hydrate sII from supersaturated aqueous methane solutions. The growth of sII hydrate on the [001] crystallographic face is achieved through utilization of a recently proposed methodology, and rates of crystal growth of 1 A/ns were sustained for the molecular models and specific conditions employed in this work. Characteristics of the crystals grown as well as properties and structure of the interface are examined. Water cages with a 5(12)6(3) arrangement, which are improper to both sI and sII structures, are identified during the heterogeneous growth of sII methane hydrate. We show that the growth of a [001] face of sII hydrate can produce an sI crystalline structure, confirming that cross-nucleation of methane hydrate structures is possible. Defects consisting of two methane molecules trapped in large 5(12)6(4) cages and water molecules trapped in small and large cages are observed, where in one instance we have found a large 5(12)6(4) cage containing three water molecules.