Quantification of Methane Hydration Energy Through Free Energy Perturbation Method

1. Introduction

Computational modeling and simulation serve as valuable complements to experimental studies of energetic materials [1,2,3,4,5,6]. Atomistic simulations can accurately determine many key material properties, helping guide or supplement laboratory work. This approach can significantly reduce the scope of investigations, lowering financial costs and minimizing physical risks. To ensure thermodynamic rigor and relevance to experimental conditions, it is crucial that an atomistic model is fully characterized within its intended application range [3,7,8,9,10,11,12].

In the context of free energy, calculations of hydration play a significant role in understanding various chemical and biochemical processes, such as solubility, solvent effects on equilibrium and reaction rates, and partition coefficients. The transfer of a solute from the gas phase to a solvent, characterized by the free energy of hydration, is a fundamental thermodynamic quantity. This calculation involves determining the difference in Gibbs free energy between the solute in the gas phase and in the solution phase at a given temperature and pressure [13,14,15,16].

Statistical perturbation theory and thermodynamic integration are the predominant computational methods used for these calculations [17,18]. Statistical perturbation theory, in particular, has shown high precision for computing the hydration free energies of small organic molecules. This method involves gradually changing the interaction parameters of a solute in a solvent and calculating the corresponding free energy changes. Despite the method's precision, it has traditionally been used for relative free energy calculations rather than absolute free energies, except in a few cases such as noble gases and small molecules like methane [19,20].

Molecular dynamics (MD) simulations provide a powerful tool for studying the microscopic details of solvation processes. Among various MD simulation packages, LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) is widely used due to its flexibility and efficiency in simulating large systems. However, LAMMPS does not have built-in functions for performing free energy calculations directly. This limitation can be addressed using the Free Energy Perturbation (FEP) method, which involves the calculating the free energy difference associated with this process [21,22].

In the context of methane hydration, FEP can be utilized to compute the free energy difference when a methane molecule is transferred from the gas phase into water. This process involves simulating the methane molecule and calculating the interactions with the surrounding water molecules. William et al. calculated the free energies of hydration of methane using FEP via Monte Carlo simulations with the BOSS program, achieving results that were in good agreement with experimental data [23]. Subsequently, Shirts et al. employed the Bennett acceptance ratio method using GROMACS, which provided an improved free energy of hydration for methane, yielding results that closely matched experimental values [24]. In this study, we aim to calculate this free energy using FEP through LAMMPS to compare its effectiveness with the previous studies.

By implementing a procedure to perform these calculations in LAMMPS, we seek to provide accurate and efficient results comparable to those obtained with other MD simulation packages. This work not only demonstrates the application of FEP in LAMMPS but also contributes to the broader understanding of solvation thermodynamics in molecular simulations. The second section offers a detailed exposition of the initial force field models, molecular dynamic setup, and the FEP method. Section 3 explores our findings, highlighting the results have been presented for the absolute free energies of hydration for methane, using the TIP4P model of water. This approach has allowed for a detailed analysis of both methodological issues and the accuracy of the intermolecular potential functions used, and finally, Section 4 summarizes our finding.

2. Materials and Methods

Here's a well-organized methods section for calculating the free energy of methane in water using the Free Energy Perturbation method in LAMMPS. This section incorporates insights from the provided articles and the LAMMPS documentation [21,25].

A key aspect of free energy calculations involves the perturbation of interaction parameters and the recalculation of the pair potential energy without changing the atomic coordinates from unperturbed system. This approach can be used to calculate free energy differences. The potential energy of the system includes three terms: the background term U_bkg (interaction sites whose parameters remain constant), the reference term U_ref (initial interactions of the atoms that will undergoes perturbation), and the term U_pert corresponding to the final interactions of these atoms. A coupling parameter λ varying between 0 and 1 connects the reference and perturbed systems is shown as below:

$$U\left(\lambda \right)=U_{bkg}+U_{ref}\left(\lambda \right)+U_{pert}\left(\lambda \right)$$

(1)

In order to perform a free energy calculation, we need to define the end states. For solvation free energy calculations, one end state is the isolated molecule in the gas phase along with the solution phase in the background. The other end state is the solvated molecule in the solution phase. In the first state, the atoms of the isolated molecule interact only with each other and are decoupled from the surrounding environment. In the final state, the atoms of the molecule interact with all the atoms present in the system, including the solvent and other dissolved solutes, making it the coupled state. Two common paths are that “Forward” route (involves deleting methane from water) and the “Backward” route (involves introducing methane into water), in which interactions of sites that are being created or deleted are treated using Lennard-Jones and Coulomb potentials with soft cores [26].

To compute the total free energy difference, we sum the free energy differences between successive states along the chosen route. The free energy difference along each path is determined by computing the free energy change of turning on the Lennard-Jones (LJ) and Coulomb potentials interactions between the solute and the solution. The numerical estimation of the free energy difference between two states i and f, characterized, can be obtained using the Zwanzig formula as below:

$${∆G}_{i\to f }=-k_{B}T {\sum }_{j=0}^{j=n-1}\mathit{ln}<\mathit{exp}(-{\displaystyle \frac{U\left({\lambda }_{j+1}\right)-U\left({\lambda }_j\right)}{k_{B}T}})>$$

(2)

where T is the ambient temperature, k_B is the Boltzmann constant and < denotes an ensemble average by the microstates of the starting state i [17,27,28,29,30,31]. The total free energy difference is obtained by summing the individual free energy differences, which directly correlates with experimentally measurable values. To avoid singularities during free energy calculations, a soft repulsive core is used. When the parameter λ tends to 0, the pair interaction vanishes while maintaining this soft repulsive core. As λ approaches 1, the pair interaction transitions to the normal, non-softened potential [26]. Lennard-Jones and Coulomb potentials modified by a soft core as below:

$$U_{LJ}\left(r,\lambda \right)=\lambda 4\epsilon \left({\displaystyle \frac{1}{{\left[{\left(0.5{[\AA }^2\right])\left(1-\lambda \right)}^2+{\left(\frac{r}{\sigma }\right)}^6\right]}^2}}-{\displaystyle \frac{1}{{0.5\left(1-\lambda \right)}^2+{\left(\frac{r}{\sigma }\right)}^6}} \right)$$

(3)

$$U_{coul}\left(r,\lambda \right)=\lambda \kappa {\displaystyle \frac{q_i{q}_j}{{\left[\left(10{[\AA }^2\right]){(1-\lambda )}^2+r^2\right]}^{1/2}}}$$

(4)

where, r is the distance between two atoms, $\kappa =N_A/4\pi \epsilon$, and q is the charge of particle.

We perform hydration free energy calculations using the procedure laid out in the previous sections. A cubic simulation box with periodic boundary conditions was created by fftool, containing one methane molecule was solvated in 1000 TIP4P water molecules [32,33,34]. Lennard-Jones parameters and partial charges for methane and water were assigned according to standard OPLS-AA and TIP4P values. Interactions of sites that are being created or deleted are handled using Lennard-Jones and Coulomb potentials with soft cores to avoid singularities. The initial configuration was energy minimized to remove any bad contacts using the conjugate gradient algorithm. We used a time step of 1.0 fs. Bonds involving hydrogen atoms were constrained using the SHAKE algorithm. Nonbonded interactions were cut off at 1 nm, with long-range electrostatics managed by Particle-Particle-Particle-Mesh (PPPM) summation method. The NPT ensemble was simulated using the Nose–Hoover thermostat and barostat at a temperature of 298.15 K and a pressure of 1.0 atm. - Equilibration ensured the system reached a stable density and temperature, which are prerequisites for accurate free energy calculations. All systems were equilibrated for 4 ns followed by 10 ns data collection period. Furthermore, ∆λ choose 0.1, and the fix adapt/fep command is used to dynamically adjust simulation parameters (such as pair coefficients and atomic charges) over the course of the simulation according to the value of λ. This allows for a smooth transition between states (e.g., from non-interacting to fully interacting states). The simulations were run on a high-performance cluster (HPC) with LAMMPS as parallel tasks on 128 CPUs.

3. Results and Discussion

In this section, we present the results of the free energy calculations for the hydration of methane using the FEP method with LAMMPS. The simulations were carried out as described in the previous section, and the following results were obtained.

Both the forward and backward routes were tested, and the results were compared. The free energy differences obtained from both pathways were consistent, demonstrating that the methodology is robust and independent of the chosen path. The free energy difference associated with the hydration of methane was calculated by gradually varying the coupling parameter λ from 0 to 1 (backward) or from 1 to 0 (forward). Figure 1 shows the Gibbs free energy (Zwanzig equation) difference with respect to the coupling parameter λ for both processes. For the x-axis, we chose the midpoint between λ_i and λ_i+1. The smooth transition between states indicates that the FEP method was successfully implemented due to the use of soft-core potentials. In Figure 1 (a), the red arrow shows the simulation direction from λ = 0 to 1, while in Figure 1(b), it shows the direction from λ = 1 to 0.

In the backward route simulation (Figure 1 (a)), the Gibbs free energy starts from a small value, approximately 0.2, when λ change from 0 to 0.1. As λ increases, $∆$G initially rises, reaching a peak at λ $\cong$0.2. This indicates a significant energy barrier in the initial stages of coupling the methane molecule to the solvent environment. Following the peak, ∆G gradually decreases, stabilizing around λ of 0.6. This suggests that once the initial interactions are established, the system stabilizes, and the energy required for further coupling decreases. From λ of 0.6 to 1.0, the ∆G varies relatively slow, indicating that the system has reached a near-equilibrium state where additional interactions contribute minimally to the overall free energy change.

In the forward route simulation (Figure 1 (b)), the process starts at λ = 1, where the methane molecule is fully coupled with the solvent. As λ decreases, ∆G further drops, reaching a minimum at λ $\cong$0.22. This suggests that the energy cost for decoupling the methane molecule is highest at this stage of the process. Beyond λ of 0.2, ∆G begins to rise, indicating that the system starts to regain some of the energy as the interactions between the methane molecule and the solvent decrease, by λ = 0, the free energy returns close to zero, indicating the decoupled state.

Both forward and backward processes display characteristic energy profiles that reflect the stability and interaction dynamics of the methane solvation process. The backward route shows an initial energy barrier, followed by stabilization, while the forward route shows an initial stabilization followed by a decrease in energy as the solvation interactions are removed.

The total free energy was compared with experimental data and previous computational studies to validate the accuracy of our simulations. Table 1 compares the total Gibbs free energy (ΔG) for the hydration of methane obtained from various methods. This study, using FEP in LAMMPS, reports ΔG_Forward of -2.218 kcal/mol with a 13% deviation and ΔG_Backward of 2.132 kcal/mol with an 8% deviation from experimental values. Previous studies using BAR in GROMACS [24] and FEP-MC [23] (Monte Carlo simulations were carried out with the BOSS program) reported ΔG_Backward of 2.27 kcal/mol with a 15% deviation and 2.90 kcal/mol with a 48% deviation, respectively, while FEP-MC also reported ΔG_Forward of -2.08 kcal/mol with a 6% deviation. However, in principle, the forward and backward free energy values should be equal in magnitude but have opposite signs. The experimental values for ΔG reported a value of $\pm$2.005 [35] and $\pm$1.912 [36] kcal/mol. For comparison with computational values, we have used the average value of $\pm$1.9585 kcal/mol. This study's results agree to the experimental values and also the previous computational studies, especially in the backward route, indicating higher accuracy and reliability of the FEP calculations performed in LAMMPS. The observed ∆G for both routes demonstrates the hysteresis effect commonly seen in such free energy perturbation calculations, indicating that the paths are not perfectly reversible, likely due to system equilibration differences. However, compared to the previous study using the FEP-MC method, which reported a significant difference of approximately 0.82 kcal/mol between the forward and backward routes, our results demonstrate a notable improvement, with a much smaller difference of only about 0.086 kcal/mol.

Radial distribution functions (RDFs) provide insights into the structural arrangement of molecules around a central molecule. The Figure 2 represents the RDFs for Carbon-Oxygen of water (C-Ow), Carbon-Hydrogen of water (C-HW), Hydrogen of methane-Oxygen of water (H-OW), and Hydrogen of methane-Hydrogen of water (H-HW) interactions in a methane-water system. The g_C-OW(r) and g_C-HW(r) RDFs both exhibit a prominent first peak at around 3.5 Å, indicating strong interactions between the carbon atom of methane and the oxygen atoms of water, as well as the hydrogen atoms of water. These peaks suggest the formation of a well-defined first hydration shell around the carbon atom, with the peak heights indicating almost equal interaction strengths for these pairs. In contrast, the g_H-OW(r) and g_H-HW(r) show their first peaks at around 4.7 Å. This shift to a larger distance for the first peak suggests weaker interactions between the hydrogen atoms of methane and the surrounding water molecules, both with the oxygen and the hydrogen atoms of water. The relative heights of these peaks indicate that while these interactions are significant, they are less pronounced compared to those involving the carbon atom of methane. The differences in peak positions and heights across the RDFs highlight the varying interaction strengths and spatial arrangements within the methane-water system. The closer peaks for C-Ow and C-HW reflect stronger and more immediate interactions, forming a tighter hydration shell around the carbon atom, whereas the more distant peaks for H-OW and H-HW indicate a more diffuse and less structured interaction landscape for the hydrogen atoms of methane. This structured arrangement around methane is crucial for understanding the solvation dynamics and the energetics of methane in aqueous environments. These results in agreement with previous study [37].

4. Conclusions and Future Works

We have calculated the free energy of hydration for methane using the Free Energy Perturbation method within the LAMMPS molecular dynamics simulation framework. We were able to quantify the free energy changes associated with the solvation process, providing valuable insights into the thermodynamic behavior of methane in aqueous environments. Our simulations were meticulously designed, using the TIP4P water model and Lennard-Jones and Coulomb potentials with soft cores, to ensure accurate and stable results. The consistency observed between the forward and backward free energy calculations highlights the robustness and reliability of our methodological approach. The results demonstrated good agreement with experimental data and previous computational studies, particularly in the backward route, which highlights the accuracy of the FEP calculations performed. The total Gibbs free energy difference for methane hydration of -2.218 kcal/mol (forward) and 2.132 kcal/mol (backward), with deviations of 13% and 8% from experimental values, respectively. These results were consistent with experimental values and previous computational studies. Our study highlighted the robustness and accuracy of the FEP method, particularly in the backward route, where our results showed only an 8% deviation from experimental values. Radial distribution functions revealed distinct structural arrangements in the methane-water system. Strong interactions were observed between the carbon atom of methane and both oxygen and hydrogen atoms of water, with prominent peaks at 3.5 Å, indicating a well-defined hydration shell. In contrast, interactions involving the hydrogen atoms of methane displayed peaks at around 4.7 Å, reflecting weaker interactions and a more diffuse hydration structure. These findings highlight the varying interaction strengths and spatial arrangements, consistent with previous studies, providing key insights into the solvation dynamics of methane in water.

Our findings emphasize the importance of detailed methodological setups and parameter optimization in molecular simulations, which are critical for obtaining reliable and reproducible results. This work not only validates the use of FEP in LAMMPS for free energy calculations but also enhances our understanding of solvation thermodynamics. Future studies could extend this approach to other solute-solvent systems, further exploring the capabilities and limitations of FEP and other free energy calculation methods in molecular dynamics simulations.