A Comparative Study of Diffusion Dynamics in Methane (CH4) and Fluoroform (CHF3) using Molecular Dynamics

Abstract

This study investigated the diffusion behaviour of Methane (CH₄) and Fluoroform (CHF₃) in an aqueous environment, using a Molecular Dynamics simulation with Groningen Machine for Chemical Simulations (GROMACS). Using the OPLS-AA force field, solvated systems were simulated to understand the rate of diffusion in a polar solvent like water, both as pure water and as a saline solution. The diffusion coefficient for each moleculewas calculated using the Einstein relation. This paper offers insight into the comparative diffusivity of relatively small molecules in water and concludes that because CH₄ is lighter than CHF₃, its diffusion coefficient is higher than that of the latter. The result of the study also shows that, generally, a molecule will diffuse faster in pure water than in saline solutions. Molecular mass affects diffusivity. However, in solution, diffusion reflects solute–solvent interactions and also has contributions from effective hydrodynamic size and solvation structure. CH₄ and CHF₃ molecules are different not only in terms of mass but also in terms of polarity, which affects their interaction with water molecules and leads to the observed differences in diffusion coefficients.

Introduction

Diffusion is a process involving the net movement of a substance from a region of higher concentration to a region of lower concentration, driven by a gradient of Gibbs free energy or chemical potential¹. This movement of particles is observed everywhere in nature. For example, diffusion occurs in our lungs when we breathe and affects the rates of chemical reactions. It is understood that in order for a chemical reaction to occur, atoms/molecules that are known as the reactants must collide together with sufficient energy. The frequency of the collisions of the reactant molecules is determined by how fast those molecules diffuse, which results in the variation of the rate of reaction, as it is dependent on the collision frequency.

Molecular dynamics is a way to intuitively understand the interactions between tiny molecules, their collisions, and their net movement, known as diffusion, which will be the theme of this paper.² Molecular dynamics helps us to capture the mass and transport mechanisms at the atomic scale.

The benefit of these simulations is that they are able to predict the trajectories of molecules and use theoretical equations and principles to give ideal experimental values. This is done by applying fundamental laws of Newtonian mechanics to a large number of microscopic molecules to explain macroscopic effects such as pressure, temperature, and volume³. Molecular dynamics simulation software, such as GROMACS⁴, allows for reduced experimental error when carrying out investigations.

The two molecules that this investigation was carried out on were Methane and Fluoroform:

Methane, the simplest hydrocarbon, was chosen for the study because of its fundamental structure. Its smaller molecular size affects factors such as the diffusion coefficient⁵.
Fluoroform has a heavier molecule than Methane. This helps us compare the effect of the mass of a molecule to its diffusion rate (diffusion coefficient).⁶

Methodology

GROMACS uses the gmx msd function in order to calculate the diffusion coefficient of a molecule. MSD, a quantity commonly used in Statistical mechanics, stands for Mean square displacement. MSD is the average distance of a particle from its starting point over a short period of time. “Equation of motion is calculated using leapfrog algorithm with an integration time step”⁷.  In our case, this time step was 2 fs.“All diffusion coefficients are calculated with the mean-squared-displacement (MSD) method which is based on the Einstein relation that MSD is linear with diffusion time.”⁸.

(1)  

Einstein relation for diffusion, but , therefore:

(2)  

where is the number of dimensions. Since there are 3 dimensions in the simulation, .

The simulation times were kept short and limited sampling was carried out as this study is focused mainly on trend comparison rather than precision.⁹ The investigation comprised three different simulations:

Methane in water (referred to as CH₄), Fluoroform in water (CHF₃)and CH₄in salt-water (CH₄SW).

Throughout all the simulations, the OPLS-AA force field was used.¹⁰

The structure file of CH₄ was created using ChemSketch¹¹ which was then imported to GROMACS. The structure files can be viewed using software such as Avogadro¹² and VMD¹³, as shown in Figure 2.1 and 2.2.

Next, topology files – specifically the .top and .gro files – were generated. This was done by giving it the input of the .pdb file or SMILES code using OpenBabel¹⁴. The SMILES code for the CHF₃ molecule’s topology is: FC(F)F, as given by LigParGen¹⁵. Additionally, .mdp files were needed to configure the parameters of the simulations.

Next, a simulation box was defined using gmx editconf:

gmx editconf -f CH₄.gro -o newbox.gro -c -d 1.0 -bt cubic.

In this case, a cubic box having periodic boundary conditions was created with the dimensions 2.56nm×2.56nm×2.56nm. The two molecules were placed in relatively larger boxes to avoid possible solute–solute interactions. Additionally, using a smaller box results in greater finite size effects acting on it, causing suppressed long-range correlations and leads to an underestimation of molecular movement. This decreases the accuracy of the calculated diffusion coefficient.

Next, the box was solvated by using the command gmx solvate. solvate. To solvate the CH₄ molecules with ions, Na⁺  and Cl^– ions were introduced at a concentration 0.5 M, using were added into it, using the command gmx genion as shown in Figure 2.3. A saline environment of 0.5 M was used because it is similar to that of seawater, keeping a simple composition that is often assumed in simulations.

Before the production simulations, energy minimisation (EM) was performed to relax the system to a mechanically stable configuration.

Molecular dynamics simulations were then conducted under the NVT ensemble, in which the number of particles and volume were kept constant while the temperature was controlled using a thermostat. NVT abbreviates the following: N is the number of particles; V is the Volume; T is Temperature. Pressure is not controlled, therefore the box size stays constant.

This was followed by the NPT ensemble, where both temperature and pressure were regulated, allowing the system volume to reach the target pressure. NPT abbreviates the following: N is the Number of particles; P is the Pressure; T is Temperature. The box’s volume changes to maintain the desired pressure. The thermostat coupling constant was 0.1 ps and  the barostat coupling constant was 2.0 ps.

The NVT, NPT and production molecular dynamics (Production MD) simulations were carried out at different time lengths for different molecules in different conditions. This was due to some systems stabilising before others, in terms of temperature and pressure.

Energy Minimisation

With the solvated, neutral system created, the next step was energy minimization. Energy minimisation is an important step to ensure that the system has a stable structure.¹⁶ We can initiate the simulation by using the gmx grompp

command:

gmx grompp -f minim.mdp -c solv.gro -p topol.top -o em.tpr

Figure 3.1 shows the Energy Potential Graph, obtained by the command gmx energy. As the energy potential decreases, the system becomes more stable.

Table 3.1 shows the maximum force threshold  for the energy minimisation and the no. of steps taken for this.

Parameter	Value
Integrator	steep (steepest descent)
Maximum force (emtol)	1000.0 kJmol-1nm-1
Minimization step size (emstep)	0.01 nm
Maximum steps (nsteps)	50000
Cutoff scheme	Verlet
Neighbor search type	grid
Coulomb type	PME
Coulomb cutoff	1.0 nm
van der Waals cutoff	1.0 nm
Periodic boundary conditions	xyz
Table 3.1 | EM Parameters

NVT Equilibration

The NVT Equilibration simulations utilised a thermostat to regulate and maintain the system temperature at a set point. For the NVT ensemble, a V-Rescale thermostat was used to control the velocities of each individual particle depending on the deviation of the temperature of the system from the set point (298 K). The V-Rescale thermostat adds to the Berendsen thermostat by adding accurate temperature fluctuations and it is a rather sensitive thermostat.¹⁷ The simulation can be initiated by using the gmx grompp command:

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr

Parameter	Value
Integrator	Md
Number of steps	50000
Time step (dt)	0.002 ps
Constraint algorithm	LINCS
Constraints	h-bonds
Cutoff scheme	Verlet
Coulomb cutoff	1.0 nm
van der Waals cutoff	1.0 nm
Electrostatics	PME
Thermostat	V-rescale (300 K)
Table 4.1 | NVT Equilibration Parameters

Table 4.1 shows the parameters that determine the length and accuracy of our NVT Equilibration.

Figures A1, A2, and A3 in the supplementary material show the graphs of temperature equilibration during the NVT ensemble.

NPT Equilibration

During NPT equilibration, both the  temperature and pressure were controlled using a V-Rescale thermostat and a Berendsen barostat, respectively. A barostat regulates the pressure of a system around a setpoint by changing the volume of the box gradually, while constantly checking the pressure at intervals. The Berendsen barostat is great for a quick stabilisation, but not very useful for analysis as it smooths out the fluctuations, which makes the system behave unnaturally. As a result, it does not accurately sample the NPT ensemble, and therefore it is  used only in the equilibration steps. This simulation can be initiated by using the gmx grompp command:

gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr

Parameter	Value
Integrator	md
Number of steps	50000
Time step (dt)	0.002 ps
Constraint algorithm	LINCS
Constraints	h-bonds
Cutoff scheme	Verlet
Coulomb cutoff	1.0 nm
van der Waals cutoff	1.0 nm
Electrostatics	PME
Thermostat/Barostat	V-rescale (300 K) / Parrinello-Rahman (1 bar)
Table 5.1 | NPT Equilibration Parameters

Table 5.1 shows the parameters that determine the length and accuracy of our NPT Equilibration.

Figures B1, B2, and B3 in the supplementary material show the graphs of pressure equilibration during the NPT equilibration.

Production MD Run

Once the equilibration phase was complete, the simulation moved to the Production Molecular Dynamics stage. This was the final stage of the simulation that produced the final result for data collection¹⁸. The simulation was initiated by using the gmx grompp command:

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr

The simulation was subsequently initiated with a small runtime and limited sampling was carried out. This is because the study is focused mainly on comparing trends rather than obtaining precise data.

Figure 6.1 shows the molecule CH₄ in Avogadro at the end of the production MD simulation. It is the visual representation of the final simulation.

Table 6.1 contains the different simulation parameters for the MD Run, obtained from a simulation file. These parameters are the same to all 3 simulations except for the simulation lengths.

Parameter	Value
Integrator	md
Number of steps	500000
Timestep	0.002 ps
Constraints	h-bonds
Cutoff scheme	Verlet
Coulomb cutoff	1.0 nm
van der Waals cutoff	1.0 nm
Electrostatics	PME
Thermostat type	V-rescale
Barostat type	Parrinello-Rahman
Table 6.1 | MD Run Parameters

Figures C1, C2, and C3 in the supplementary material show the graphs of the temperature equilibration, while figures D1, D2 and D3 show the graphs of pressure equilibration during the MD production run.

Results and Discussion

The diffusion coefficient for each of the simulations was then calculated using the function gmx msd. From the Einstein relation for diffusion, it can be deduce that:

(3)  

The mean square displacement (MSD) graphs were plotted and the diffusion coefficients were determined using the gmx msd function:

 gmx msd -f md.xtc -s md.tpr -n index.ndx -o msd.xvg

A log-log graph of MSD against time was used to determine the linear portions of the graph. This is shown by the graphs in Figure 7.1 with non-linear scales.

After plotting these graphs, the flags -beginfit and -endfit were used to add the limits of the linear portion and redo the simulation. This made the calculated diffusion coefficient more accurate. The diffusion coefficients are obtained from GROMACS are shown in Table 7.1.

Molecule Name	Diffusion Coefficient	Literature Diffusion Coefficient
CH₄
CHF₃		(estimated)
CH₄ SW		(estimated)
Table 7.1 | Diffusion Coefficients

Table 7.1 compares the diffusion coefficients achieved by the simulations being run compared to literature. However, due to the limited literature values on CHF₃ in water and CH₄ in salt water, they are given according to widely accepted or estimated values. The literature value for CH₄ in water given by Chen et al. 2018.¹⁹

Running simulations using various initial velocity seeds and calculating block average for the mean square displacement will definitely increase the statistical reliability and reduce uncertainty. In this paper, only single trajectories for calculation were used due to computational limitations. Even with single trajectories, the simulations show expected qualitative trends for molecular diffusion. These include slower diffusion for heavier solutes and less mobile behaviour in saline solutions. The lack of block averaging and uncertainty analysis is recognised as a limitation of this study. The manuscript clarifies that high uncertainty in the values of the diffusion coefficients arises from trajectory length and sampling and contains plan for next stage of work. This will also involve block averaging among other error analysis methods.

Molecular mass affects diffusivity, however, in solution, diffusion reflects solute–solvent interactions and also has contributions from effective hydrodynamic size and solvation structure. CH₄ and CHF₃ molecules are different not only in terms of mass but also in terms of polarity. This could affect their interaction with water molecules and lead to different results being obtained. The current findings are therefore to be considered as a result of the effect of both molecular mass and solvation-related interactions.

NVT and NPT equilibration, and production MD:

During NVT equilibration, NPT equilibration and production MD, the temperature is stabilised to around 298K and the pressure stabilised to about 1 bar. This is demonstrated by the graphs in the supplementary material, which have a 10ps running average that fluctuates about the set-point at an acceptable amount to be considered as stabilised.

Mean square displacement graphs:

When analysing the simulations to measure diffusion, the mean square displacement (MSD) of the molecules in the simulations must be considered. In order to obtain accurate values of the diffusion coefficients, a graph of MSD against time with non-linear scales, must be plotted. This is done by plotting the logarithm of MSD against the logarithm of time. Note that the logarithms are to the base 10. In Figure 7.1, the linear portions of the log MSD – log time graphs can be observed in order to choose the limits for the accurate calculation of the diffusion coefficient. The MSD graphs  were also plotted with a linear scale.²⁰

Diffusion coefficients:

The diffusion coefficients were calculated by GROMACS, as listed in Table 7.1. They indicate that CH₄ in water has the highest rate of diffusion out of the three simulations, followed by CHF₃in water, and lastly, CH₄ in salt water. This is due to the larger mass of a CHF₃molecule as compared to CH₄ and the slower diffusion of molecules in salt water as compared to pure water. CH₄ is a molecule that is much lighter than CHF₃. CH₄ has a mass of approximately 16.0425 gmol^-1, while CHF₃ has a mass of approximately 70.0138 gmol^-1. The diffusion coefficient of CH₄, as given in Table 7.1, is higher than the diffusion coefficient of CHF₃_. This means that the less the (molar) mass of a molecule, the faster it can diffuse in water. If NaCl salt is added to water of concentration 0.5M, we observe another decrease in the diffusion coefficient of CH₄. This means that, generally, a molecule will diffuse faster in pure water than in water with salt.²¹

Conclusion

This study investigated the diffusion behaviour of Methane (CH₄) and Fluoroform (CHF₃) in a water environment, using a molecular dynamics simulation. There were 3 separate sets of molecular dynamics simulations that were carried out. The first two simulations were a comparison of the mass of a molecule to its rate of diffusion, while the last and first simulations were a comparison of the rate of diffusion of a molecule to the presence of ions. The diffusion coefficient for each of CH₄ and CHF₃ was calculated using the Einstein relation for diffusion.²² The simulations carried out confirm the claim that due to CH₄ being much lighter than CHF₃,the diffusion coefficient of CH₄ is higher than the diffusion coefficient of CHF₃_. The simulations also indicate that, generally, a molecule will diffuse faster in pure water than in salt water. The findings are to be considered as the result of both molecular mass and solvation-related interactions. Future work will involve block-averaging with several simulations that run for longer times in order to increase the statistical reliability.

Supplementary Information

References

1. P. Naeiji, T. K. Woo, S. Alavi and R. Ohmura, 2020, “Molecular dynamics simulations of interfacial properties of the CO2–water and CO2–CH4–water systems”, The Journal of Chemical Physics, Vol. 153, pp. 044701-1 to 044701-13, DOI: 10.1063/5.0008114, (see p. 044701-5). [↩]
2. G. Battimelli and G. Ciccotti, 2018, “Berni Alder and the pioneering times of molecular simulation”, The European Physical Journal H, Vol. 43, pp. 303–335, DOI: 10.1140/epjh/e2018-90027-5, (see p. 304). [↩]
3. A. Astuti and A. Mutiara, 2009, “Performance Analysis on Molecular Dynamics Simulation of Protein Using GROMACS”, December 2009 (conference or preprint), DOI: 10.13140/RG.2.1.3245.2885, (see p. 2). [↩]
4. M. J. Abraham et al., GROMACS – Groningen Machine for Chemical Simulations, Zenodo (2019). [https://doi.org/10.5281/zenodo.14846130] [↩]
5. M. Wahlen, 1993, “The Global Methane Cycle”, Annual Review of Earth and Planetary Sciences, Vol. 21, pp. 407–426, DOI: 10.1146/annurev.ea.21.050193.002203, (see p. 407). [↩]
6. P. A. Witherspoon and D. N. Saraf, 1965, “Diffusion of Methane, Ethane, Propane, and n-Butane in Water from 25 to 43°”, The Journal of Physical Chemistry, Issue 11, Vol. 69, pp. 3752–3755, DOI: 10.1021/j100895a017, (see p. 3752). [↩]
7. Y.-A. Chen, C.-K. Chu, Y.-P. Chen, L.-S. Chu, S.-T. Lin and L.-J. Chen, 2018, “Measurements of diffusion coefficient of methane in water/brine under high pressure”, Terrestrial, Atmospheric and Oceanic Sciences, Vol. 29, pp. 577–587, DOI: 10.3319/TAO.2018.02.23.02, (see p. 4 of the accepted manuscript). [↩]
8. X. Zhao and H. Jin, 2020, “Correlation for diffusion coefficients of H2, CH4, CO, O2 and CO2 in supercritical water from molecular dynamics simulation”, Applied Thermal Engineering, Vol. 171, pp. 114941, DOI: 10.1016/j.applthermaleng.2020.114941, (see p. 6 of the accepted manuscript). [↩]
9. V. Marry, P. Turq, T. Cartailler and D. Levesque, 2002, “Microscopic simulation of structure and dynamics of water and counterions in a monohydrated montmorillonite”, Journal of Chemical Physics, Issue 7, Vol. 117, pp. 3454–3463, DOI: 10.1063/1.1493186, (see p. 3461). [↩]
10. W. L. Jorgensen and J. Tirado-Rives, 2005, “Potential energy functions for atomic-level simulations of water and organic and biomolecular systems”, Proceedings of the National Academy of Sciences of the United States of America, Issue 19, Vol. 102, pp. 6665–6670, DOI: 10.1073/pnas.0408037102, (see p. 6665). [↩]
11. ChemSketch, Version 2022.1.2, Advanced Chemistry Development, Inc., Toronto, Canada, 2022. Available at: [www.acdlabs.com]. [↩]
12. M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek and G. R. Hutchison, 2012, “Avogadro: An advanced semantic chemical editor, visualization, and analysis platform”, Journal of Cheminformatics, Vol. 4, pp. 17, DOI: 10.1186/1758-2946-4-17, (see p. 1 of the article).;Avogadro, Version 1.2.0, [http://avogadro.cc/] (accessed June 2025). [↩]
13. W. Humphrey, A. Dalke and K. Schulten, 1996, “VMD: Visual Molecular Dynamics”, Journal of Molecular Graphics, Vol. 14, pp. 33–38, DOI: 10.1016/0263-7855(96)00018-5, (see p. 33). [↩]
14. N. M. O’Boyle, M. Banck, C. A. James, C. Morley, T. Vandermeersch and G. R. Hutchison, 2011, “Open Babel: An open chemical toolbox”, Journal of Cheminformatics, Vol. 3, pp. 33, DOI: 10.1186/1758-2946-3-33, (see p. 2); Open Babel, Version 3.1.1. [http://openbabel.org] [↩]
15. L. S. Dodda, I. Cabeza de Vaca, J. Tirado-Rives and W. L. Jorgensen, 2017, “LigParGen web server: an automatic OPLS-AA parameter generator for organic ligands”, Nucleic Acids Research, Vol. 45, pp. W331–W336, DOI: 10.1093/nar/gkx312, (see p. 4); L. S. Dodda, J. Z. Vilseck, J. Tirado-Rives and W. L. Jorgensen, 2017, “1.14*CM1A-LBCC: Localized Bond-Charge Corrected CM1A Charges for Condensed-Phase Simulations”, The Journal of Physical Chemistry B, Vol. 121, pp. 3864–3870, DOI: 10.1021/acs.jpcb.7b00272, (see p. 3864). [↩]
16. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, 1987, “The missing term in effective pair potentials”, The Journal of Physical Chemistry, Issue 24, Vol. 91, pp. 6269–6271, DOI: 10.1021/j100308a038, (see p. 6271). [↩]
17. I. Shvab and R. J. Sadus, 2014, “Thermodynamic properties and diffusion of water + methane binary mixtures”, The Journal of Chemical Physics, Vol. 140, pp. 104505-1 to 104505-10, DOI: 10.1063/1.4867282, (see p. 104505-3). [↩]
18. I. Shvab and R. J. Sadus, 2014, “Thermodynamic properties and diffusion of water + methane binary mixtures”, The Journal of Chemical Physics, Vol. 140, pp. 104505-1 to 104505-10, DOI: 10.1063/1.4867282, (see p. 104505-3 of the published article). [↩]
19. Y.-A. Chen, C.-K. Chu, Y.-P. Chen, L.-S. Chu, S.-T. Lin and L.-J. Chen, 2018, “Measurements of diffusion coefficient of methane in water/brine under high pressure”, Terrestrial, Atmospheric and Oceanic Sciences, Vol. 29, pp. 577–587, DOI: 10.3319/TAO.2018.02.23.02, (see p. 6 of the accepted manuscript). [↩]
20. H. Moradi, H. Azizpour, H. Bahmanyar, M. Mohammadi and M. Akbari, 2020, “Prediction of methane diffusion coefficient in water using molecular dynamics simulation”, Heliyon, Issue 11, Vol. 6, pp. e05385, DOI: 10.1016/j.heliyon.2020.e05385, (see p. 3). [↩]
21. T. Funazukuri, 2018, “Concerning the determination and predictive correlation of diffusion coefficients in supercritical fluids and their mixtures”, The Journal of Supercritical Fluids, Vol. 134, pp. 28–32, DOI: 10.1016/j.supflu.2017.11.035, (see p. 7 of the accepted manuscript). [↩]
22. M. J. W. Frank, J. A. M. Kuipers and W. P. M. van Swaaij, 1996, “Diffusion Coefficients and Viscosities of CO2 + H2O, CO2 + CH3OH, NH3 + H2O, and NH3 + CH3OH Liquid Mixtures”, Journal of Chemical and Engineering Data, Issue 2, Vol. 41, pp. 297–302, DOI: 10.1021/je950157k, (see p. 298). [↩]