Molecular Dynamics simulations of Carbon Dioxide Chlathrates




Force field development

The water force field used in this study was developed in an analogous spirit Habershon et al (Habershon 2009) for H\({}_{2}\)O and Cygan et al (Cygan 2012) for CO\({}_{2}\). For both force fields the authors tried to keep the thermodynamic properties of a rigid molecule model (TIP4P/2005 in case of the water model of Habershon) while adjusting flexible intramolecular parameter so that simulated IR spectra of the liquid H\({}_{2}\)O or CO\({}_{2}\) would best match experimental IR spectra. Since the TIP4P/ICE rigid model is based on the TIP4P/2005 model, and with parameters slightly tweaked so that it reproduces the melting point of ice I\({}_{h}\) while retaining a good agreement in the temperature of maximum density of water, we chose to modify the q-TIP4P/2005F model of Habershon to incorporate characteristics of the TIP4P/ICE model rather than the TIP4P/2005. Similarly to the flexible model of Habershon the charges and Lennard-Jones parameters were to be kept strictly the same as in the TIP4P/ICE model. Next the harmonic bending potential was taken directly from q-TIP4P/2005F and the equilibrium bond angle (\(\theta_{0}\)) was modified to match the H-O-H bond angle in the TIP4P/ICE model. The location of the negatively charged virtual site was set to the same location as in TIP4P/ICE. Since this location is dynamically based on the bisector of H-O-H angle, it will fluctuate based on instanteous the O-H bond lengths and H-O-H bond angle. The O-H bonds were treated using a Morse potential with same parameters as in q-TIP4P/2005F, but with the equilibrium bond length (R\({}_{0}\)) set to that of TIP4/ICE. With this choice of R\({}_{0}\) and \(\theta_{0}\) one is guaranteed that when one takes equilibrated configurations from a simulation using the rigid TIP4P/ICE model, and changes to this flexible parameterization, the intramolecular potential terms will initially contribute zero potential energy and that the switch over from rigid to flexible models will not be abrupt. Hence this choice eliminates the necessity to an energy minimization on the configuration obtained from the rigid model simulations. However since the bond expectation value of a Morse potential is larger than the R\({}_{0}\) parameter, the flexible model will tend to have a larger dipole moment than the rigid model, hence the melting point temperature will be higher. Since empirically we have seen that the process of freezing is slower (on the order of 30-50ns) than that of melting (on the order of 1ns) this is less problematic than if the melting point of the flexible model was lowered. In practice to simulate an IR spectrum one needs on the order of 1ns of NPT equilibration with the flexible model and 1 ns of NVE simulation to obtain an IR spectrum.

If one determines the melting point using a series of phase coexistence NPT simulations on the rigid TIP4P/ICE model one obtains a melting point around 265K (since these are single MD runs error is most likely on the order of 5-10K) not far from the best estimate of the model’s melting point (271K Abscal and Vega). Allowing the H-O-H bond angle to flex doesn’t substantially change the melting point (262K with error bar on the order of 5-10K). Since the bend is described by a harmonic potential it would be expected that the average angle would be close to that of the \(\theta_{0}\) of TIP4P/ICE. Allowing flexibility in both the angle and O-H bond lengths, increase the melting point to 292K (with error bar on the order of 5-10K).