A Lattice Distortion Theory for Promotor Containing Clathrate Hydrates

A lattice distortion theory for promotor containing clathrate hydrates is formulated using the statistical thermodynamics based model of van der Waals and Platteeuw in association with the ab initio quantum mechanics to compute the cavity potentials. Despite of high degree of lattice distortion anticipated for large and polar molecules of liquid promotors, their variable lattice energy concept is unreported. With this intention, we estimate the lattice stabilization energy from spin-component scaled second order Møller-Plesset (SCS-MP2) perturbation theory applied with the augmented correlation-consistent polarized double zeta valence (aug-cc-pVDZ) basis set. Implementing this to compute cavity potential for different promotors, the reference properties of hydrates are harvested by regressing against the phase equilibrium conditions of their binary hydrates with methane. Our study confirms the exponential relation of reference chemical potential difference with van der Waals volume of the promotors. Moreover, using the excess Gibbs free energy theory, the higher order distortions for the multiple guests are captured. The proposed lattice distortion theory is attested with phase equilibrium conditions of eight promotors containing clathrate hydrate systems, namely propylene oxide, acetone, tetrahydrofuran, pyrrolidine, iso-butanaldehyde, cyclopentane, furan and thiophene, all having methane as a co-guest.

Natural gas hydrate (NGH) has attracted the interest for its occurrence as an energy resource which is believed to have up to 15,000 billion tons of carbon compared to 5,000 billion tons of all other sources of oil and gas around the globe 1 . Stability of the hydrates is favoured by the low temperature and high-pressure conditions, which is eventually accountable for its occurrence in permafrost and marine continental slope 2,3 . This typical crystalline compound is formed by hydrogen-bonded water molecules that entrap the small gas in a cage-like structure 4 . In addition to the gaseous guests, the light organic compounds, e.g. ethers, ketones, aldehydes, refrigerants and sulphur containing compounds, are capable to form more stable hydrates that are named as clathrate hydrates in general. Moreover, if mixed in a proper proportion with water, they can form the binary hydrates at mild conditions allowing the small gaseous components to co-guest them. This apart, some organics, namely methanol and ethylene glycol, inhibit the hydrate formation and therefore increase the formation pressure. With the ability to control the formation conditions, the promotors and inhibitors are potential candidates in flow assurance, NGH energy production and hydrate based applications of desalination and gas storage 5,6 . The phase behaviour of these hydrates is quite different from each other that needs the molecular level understanding to predict their bulk level properties. Although thermodynamics of the clathrate hydrate formation is well developed in the past four decades, the stability studies of the promotor containing hydrates and their phase behaviour are one of the current attentions 7 . The hydrate cages have pentagonal dodecahedra (5 12 ) cavity as the basic building block which recombine in such a way to form different clathrate hydrate structures, amongst the sI, sII and sH are mostly found in nature 2 .
Investigating free energy for the noble gas hydrates of increasing molecular size of the guest demonstrates that their large size favours stable hydrate lattice by contributing more potential energy 8 . On the other hand, the large size of the guest molecule distorts the hydrate lattice that is evident from analysing the adsorption energies of ethane, propane and isobutane in different cages of the hydrate 9 . To model this phenomenon, the statistical thermodynamics approach of van der Waals and Platteeuw 10 is available that is originally formulated on the basis of a constant crystal lattice assumption. This leads to a unique set of reference chemical properties in the calculation of the change in chemical potential of water in empty hydrate and bulk phase 11 . A modification for spherical coordinates of the centre of mass of guest molecule with respect to the oxygen molecule of water θ φ r ( , , ) and the coordinates of all constituting atoms of the guest with its centre of mass as origin θ φ ′ ′ ′ r ( , , ). The set of radial distances ( ′ R ) of all atoms of the guest molecule is fixed for constant geometry that results in five degrees of freedoms. The constraint imposed on the angles θ and φ is [-40, 40] that accords to its location in the cavity such that it is not much close to the cage wall. The spacing for the radial distance (R) of centre of mass of guest molecule with respect to the oxygen in water is set denser near the cage wall to arrest its stiff nature in repulsive region. On revolving the guest around the water molecule, this scheme generates 6400 nodes ( Fig. 1) of PES grid. An arithmetical average of the energies computed at each radial distance ( ) E r avg results in the averaged cavity potential that can be fitted to suitable potential models.  25 , in which they used MP2/6-31 + +G(2d,2p) level of theory to compute potential energies for a total of 18000 orientations for a methane-water dimer. For improving the accuracy, they performed a high-level MP2/ cc-pVQZ at 98 selected grids from modified Plackett-Burman design to capture the effect of high-level theory to  each angular degree of freedom. The subsequent studies are reported on estimating cavity potential of the methane, argon and CO 2 hydrates used MP2 level of theory with sufficiently large basis sets 26 . However, no estimations are reported for the hydrate promotors having large molecular size. As the required computation time is considerably high for large molecular system of hydrate promotor used in the present study (say >15 atoms), this motivates us to design the ab initio methodology that is closely accurate to the coupled cluster (CC) and take computation time comparable to MPPT. In this light, we propose the spin-component scaled second order Møller-Plesset perturbation theory (SCS-MP2) by Antony and Grimme 27 . The spin-component scaled Møller-Plesset, SCS-MP2 is applied to estimate the non-covalent interaction energies of liquid organic promotors that has potential applications where the highly accurate CCSD(T) calculations cannot be performed. However, there are local correlation theories 28 developed in recent years, e.g. domain-based local pair-natural orbital (DLPNO)-CCST(T) 29,30 and local natural orbital (LNO)-CCSD(T) 31 , that can routinely compute the CCSD(T) level energies. The SCS-MP2 technique modifies the second order electronic correlation (E C ) by scaling the double excitations of electron pair in parallel (E P ) and antiparallel (E AP ) spin as Here, the scaling parameters p S and p T have default values of 6/5 and 1/3, respectively. The basis sets discretise the Schrodinger equation into the readily solvable algebraic equations. The Dunning's basis sets are designed such that the post-Hartree-Fock calculations converge systematically to the complete basis set (CBS) limit. This is utilized to correct the computed energies by analysing the basis set convergence (BSCE) and superposition (BSSE) errors. Let us take dimer of two molecules A and B. Using the counterpoise theory, the uncorrected (∆E AB raw ) and corrected (∆E AB CP ) energies are calculated in the following way E E E (4) where, the convention E A AB represents the potential energy of molecule A on the basis of AB dimer. The difference of the counterpoise corrected and uncorrected energy gives the estimation of the BSSE as follows Another error, i.e. BSCE is estimated by extrapolating the energies recorded at the increasing basis sets. Observing the expected nature of the variation of the energy values, the following equation is expected to give best fit and estimation of energy at complete basis limit 27 .
Here, B is a constant and ∆E int CBS is the desired potential energy at CBS. The symbol N represents the cardinal number that holds n=2, 3, 4, 5 and so on, values for aug-cc-pVnZ type basis sets. In this way, the BSCE is the difference of counterpoise corrected and CBS energy The can be estimated by the collective effects of these two errors as depicted in the following equation Here, = − w 1 BSCE/BSSE is the Pauling point counterpoise weight that is multiplied to the BSSE and the resultant is to be added to the raw interaction energies to get ∆E AB CBS value. This value is calculated for a total of 6400 orientations of the guest-water dimer for the cavity potential used to calculate the hydrate phase nonidealities.
Hydrate phase. The hydrate phase equilibrium occurs when the change in the chemical potential difference of water between the filled and empty hydrate ( µ ∆ β− w H ) equals to the difference between empty hydrate and liquid Here, i and j are indices for the cavity and guest. The number of cavities per water molecule (υ i ) for sII hydrate are 1/23 and 3/23 small (5 12 ) and large (5 12 6 4 ) cages, respectively. The reaction of converting empty cages into filled cages is governed by Langmuir constant (C ij ) and the fugacity ( f j ) The fugacities of the component in vapor and liquid phases are calculated using the modified PT equation of state 23 . The Langmuir constant is a measure of cavity stabilization by the effective guest-water interaction. This is estimated by volume integration of the Boltzmann weighted averaged cavity potential as Here, ω is the averaged cavity potential for which we have developed the ab initio methodology. This can be represented by a three-parameter Kihara potential model 3 . 12 6 This generalized formulation is modified for the specific hydrate cavities as follows where, ′ R and ′ z represent the radius and coordination number of the cavity, respectively. The coordination number is the count of water molecules per hydrate cavity. The Kihara potential parameters σ, ε and a are obtained by fitting Eq. (13) to the angle averaged ab initio energies ( ) E r avg estimated in Eq. (1). For hydrate equilibrium calculation, the change in chemical potential in empty and filled hydrate is equated to the change in the empty hydrate and liquid water ( µ ∆ β− w L ). For estimation of the latter quantity, the Holder's equation 32 is applied as , indicates the reference chemical potential difference, while the other three terms correct the chemical potential for operating temperature, pressure and activity, respectively. The change in specific heat (∆h w ), molar volume difference of the water in the hydrate and liquid phase (∆V w ), and activity of water (γ w ) are used for these corrections. The heat and volume corrections for the clathrate hydrates yield following expressions: 33 The values for ∆h w o and ∆V w o at standard point are −5202.2 J·mol −1 and 5.0 cm 3 ·mol −1 for sII type hydrates, respectively 33 . For the estimation of activity of water altered by presence of promotors, the modified UNIFAC 24 model is used.
Lattice distortion model formulation. Holder et al. 32 developed a method for estimating a single set of reference properties for sI hydrates irrespective of the nature of guests. The method is updated by Lee and Holder 14 for calculation of variable reference properties for both the pure and mixture hydrates 34 . Combining Eqs. (9) and (16), Rearranging Eq. (19) and adopting the ∆h w from Eq. (20), the following equation is obtained for reference properties calculation.
Let us consider the left-hand side as Y and Eventually, Eq. (21) can be written as = + Y SX I (22) When Y is plotted against X, one can obtain the slope (S) and intercept (I) that lead to the values of reference properties.
Equations (23) and (24) produce the reference properties for the specific guests in pure and mixture hydrates. The combined effect of dissimilar guests on reference properties can be explained with the excess Gibbs free energy-based approach 15 . For the binary mixture, if the value is known for one component, the same for another can be computed using the hydrate molar concentration weighted correlation as follows, , stands for its value at reference point (273.15 K and 0 MPa). The symbols Z 1 and Z 2 represent the hydrate phase compositions of component 1 and 2, respectively. The parameters A 12 and B 12 stand for the intreaction between component 1 and 2, respectively. The superscripts m and p denote the methane and promotor molecules, respectively. The first two terms in Eqs. (25) and (26) account for the primary lattice distortion, whereas the third and fourth terms account for the higher order distortions. Simulation algorithm and model identification. The quantum mechanical simulations are performed in GAMESS-US 35 (version: 2018-R1-pgi-mkl) for evaluating the guest-water interaction energies. The individual molecules are optimized using SCS-MP2 theory and aug-cc-pVDZ basis set. The initial configuration of the guest-water dimer is chosen in such a way that their dipole moments coincide with each other. A total of 6400 input files are generated using a MATLAB ® code for the different orientations of the guest with respect to the water molecule, as described in Fig. 1. For counterpoise energy calculation, we set zero point charge for each molecule separately to generate other two sets of input files. Single point energies are computed for all three sets of input files with the same level of theory with which the individual molecules are optimized. The counterpoise corrected (∆E AB CP ) and uncorrected (∆E AB raw ) energies are calculated using Eqs. (3) and (4). The energies are arithmetically averaged at each radial position and fitted to the Kihara potential function to obtain the collision diameter (σ) and energy well depth (ε).
Subsequently, the reference properties are estimated in Holder's equation for the change in chemical potential and enthalpy during phase change from water to empty hydrate at the reference point. The cavity potential parameters are employed to estimate the Langmuir constant (C ij ) using Eq. (12). The fugacity of guest in the equilibrium phases estimated using the modified equation of state model presented in Methods section. The Langmuir constant and fugacity of guest quantify the cage occupancy, θ ij using Eq. (11) and subsequently, the filled-empty hydrate chemical potential difference is estimated using Eq. (10). This is equated to the empty hydrate-water chemical potential difference in Holder's equation presented in Eq. (16). The change in the activity of water in the liquid phase due to presence of promotor is estimated using the modified UNIFAC model featuring in Methods section. In this way, the experimental hydrate formation conditions are applied to Eq. (21) and linearly fitted the cumulative effect of hydrate stabilization and destabilization factors to the operating temperature. The intercept and slope of this line are saved as reference chemical potential and enthalpy differences, respectively.

Results
cavity potential and hydrate cage occupancy. The estimation of BSCE and BSSE is presented for propylene oxide in Fig. 2. Four different Dunning's basis sets, namely aug-cc-pVDZ, aug-cc-pVTZ, aug-cc-pVQZ and aug-cc-pV5Z having cardinal numbers n = 2, 3, 4 and 5, respectively, are utilized to compute CBS energy. By extrapolating the energies to complete basis set using Eq. 6, the value of Pauling-point correction factor is estimated as w = 0.60 for the configuration shown in inset of Fig. 2. This quantity holds value around 0.5 for different configurations and in this way, the overall effect validates the half-counterpoise method. Consequently, we choose this method for rest of the calculations.
www.nature.com/scientificreports www.nature.com/scientificreports/ To examine the accuracy of SCS-MP2 method over other schemes, we compute the energies using SCS-MP2, general AMBER force field (GAFF) and density functional theory (ωB97X-D) for methane and tetrahydrofuran system at aug-cc-pvDZ basis set (Fig. 3). The potential energies obtained from these simulations are arithmetically averaged for the angular degrees of freedom and plotted against the intermolecular distances. The nature of the energy curves suggests that the cavity potential follows the typical Lennard-Jones 12-6 theory. Fitting the ab initio energies to Eq. (13), the estimated cavity potential parameter (ε k / , σ) for CH 4 and THF using SCS-MP2, GAFF and ωB97X-D functional and resulting cage occupancies are presented in Table 2. With respect to the   experimental values of cage occupancies, the percentage absolute error using ab initio (SCS-MP2) method is observed to be insignificant as compared to GAFF and ωB97X-D. Consequently, we recommend the SCS-MP2 theory for computing cavity potential of the promotor containing clathrate hydrates (Fig. 4).
The cavity potential parameters for different promotors estimated at SCS-MP2/aug-cc-pVDZ level are reported in Table 3. The Kihara hard core radius is set as zero because it owns a negative value while fitting ab initio energies 18 . Analysing the potential curves, an important observation is made that the acetone promotes the methane hydrate formation despite of having lower energy well depth (ε k / = 124.59 K) as it is able to make more stable hydrate structure sII than the sI hydrate of pure methane. The cavity potentials are estimated with varying the distance of guest from the cage wall, thus these generalized parameters are applicable to any size of the cage and hence to every available hydrate structures.
Quantifying the extent of lattice distortion. Most of the studied sII hydrate formers are self-forming and they do not need the guest gas to form the hydrates. However, the phase equilibrium data for these promotors are mostly available as binary hydrates with methane as a co-guest. Consequently, the Holder's equation produces directly the reference properties for binary hydrate ( µ ∆ w mix o , ). The phase equilibrium data of these binary hydrates of promotors + methane are used to generate the X-Y plots (Eq. (22)), whose slope and intercept refer to the enthalpy and chemical potential difference at reference condition, respectively (Fig. 5). These values are reported in Table 4 for binary hydrates of all the eight promotors.
In case of mixture hydrates, existing lattice distortion models are formulated for the similar size of gaseous guests. They assume that the effect of distortion would be the same for all type of guests 14 . In the present case of promotor containing hydrates, the co-guest is methane which is smaller enough and in consequence, the lattice distortion effect is negligible as compared to the large molecule of promotor. In this light, we modified the theory taking into account the sole effect of lattice distortion by promotors. However, the methane-promotor interactions are considered to capture the higher order distortions. In this light, we adopted the Gibbs free energy based lattice distortion model 15,16 that leads to the estimation of reference properties for the pure promotors as well as the interaction terms listed in Table 5.
All the promotors are expected to follow the trend of decrease in hydrate formation pressures with increasing in their molecular sizes. In order to quantify the extent of lattice distortion for the promotor containing clathrate hydrates, we examined the variation of the calculated reference properties with respect to the size of the guests. For this purpose, we choose the van der Waals (vdW) volume of the promotor molecules calculated using the Bondi's group contribution method 36 ,   www.nature.com/scientificreports www.nature.com/scientificreports/   where, the volume contributions of individual atoms considered in the present study are taken from literature 36 .
The parameters N B , R A and R NR are the number of bonds present, and the number of aromatic and nonaromatic rings, respectively. The estimated vdW volumes for the promotors are shown in Table 5. The RCPD estimated using the ab initio potentials shows an exponential increment with the vdW volume of the promotor molecules (Fig. 6). This shows that the larger guest molecules destabilize the hydrate lattice more than the small guests despite of its high contribution to lattice distortion. A good exponential fit of the RCPD obtained with Kihara hard core radius for sI type of hydrate is expressed in Eq. (28). However, for the liquid phase promotors making sII hydrate, the better correlation of the computed RCPD is found with vdW volume of the promotor guest molecules as depicted in Eq. (29).
vdW where, a is the Kihara hard core diameter in pm and V vdW is in Å 3 . Equations (28) and (29) are the proposed lattice distortion models for the pure (sI) and promotor containing (sII) clathrate hydrates, respectively. This proposed lattice model for promotor containing hydrate uses vdW volume of the guest that applies to all categories of liquid sII hydrate formers. Comparing to previous lattice distortion theories in Table 6, the proposed model shows an excellent correlation between the size of the guest and RCPD as compared to Lee and Holder 16 . For sI hydrate formers, correlating the Kihara hard core radius with the estimated RCPD is observed to improve slightly from 0.90 to 0.91. The two existing models of lattice distortion consider the guest size as Kihara hard core radius 14,16 and guest diameter 21 to correlate with the estimated RCPD. The smaller gaseous guests exhibit almost spherical shape, for which the molecular diameter is reasonable to quantify their sizes 21 . However, for the highly non-spherical shape of the large molecules of promotors studied in this work, the vdW volume represents realistic property. For example, the molecular diameter of THF is 5.90 Å 37 , which gives equivalent spherical volume of 107.5 Å 3 37 , while the vdW volume is 74.19 Å 3 . This difference is attributed to the fact that the vdW volume is computed using the Bondi's group contribution method 36 that sums the volumes of individual atoms, whereas the equivalent spherical volume considers the whole molecule as a sphere. This makes the vdW volume as an effective property of non-spherical  www.nature.com/scientificreports www.nature.com/scientificreports/ guest molecules. The coefficient of correlation R 2 is improved for sII from 0.87 to 0.96. It should be noted that the previous studies only consider the nonpolar hydrate formers, while the present study investigates the polar liquid hydrocarbons. The maximum extent of lattice distortion reported in nonpolar category is 1887 J·mol −1 for isobutene, while it equals to 2034 J·mol −1 for cyclopentane in the liquid promotor category.
The hydrate formation and promoting mechanism can be analysed with the RCPDs shown in Fig. 6. The values of RCPD for sI hydrate are in the range of 1450-1700 J·mol −1 , whereas the same for sII hydrates range from 200-2000 J·mol −1 . This apparently shows that the extent of lattice distortion is more in case of sI hydrates than sII hydrates. This can be explained with the energy well depth that is observed to be greater for sII hydrate formers than the sI type. Thus, the energy contribution to stabilize the lattice is more in case of sII hydrates, which leads to lower down RCPD to 350 J·mol −1 in case of propylene oxide. phase diagrams for promotor containing clathrate hydrates. The proposed lattice distortion model is applied to predict the phase equilibrium for the binary hydrates of methane and promotors. The experimental phase equilibrium data is available for the stoichiometric proportions of the promotors to water ratio, i.e. 5.56%. Figure 7a depicts a close prediction of the hydrate formation pressures for a series of different promotor containing clathrate hydrate systems. The model performance is quantified in terms of percent absolute relative deviations (%AARD) and listed in Table 7. The AARD values hold a minimum value of 0.85% for four data points of methane + propylene oxide hydrate, while the value has a maximum value of 6.52% for 48 data points of methane + acetone hydrate. In the case of acetone, the phase predictions are shown for varying promotor dosage in Fig. 7b. The reasonable values of %AARD for the promotor containing clathrate hydrates attest the validity of the proposed lattice distortion model. Important findings of this study can be summarized as (i) formulating the ab initio methodology to compute cavity potential for promotor containing hydrates, (ii) implementing these potentials to estimate hydrate stabilization energies, (iii) regressing the reference properties against the stabilization energies and phase equilibrium conditions, (iv) obtaining a relation between reference properties with size of the promotor guests and finally, (v) using this lattice distortion model to predict the phase equilibrium conditions of promotor containing clathrate hydrates. Based on this proposed theory, it is clear that the vdW volume of the promotor molecules defines the extent of lattice distortion for hydrate promotors. As a future perspective, the study on vibrational frequencies for these promotors entrapped in the cavities can be a potential tool to envisage the relative stability of these hydrates using the loose-cage-tight-cage (LCTC) model 38 . Furthermore, this study can be extended to structure H type   www.nature.com/scientificreports www.nature.com/scientificreports/ hydrates including the shape effects of hydrate cavities. The generalized approach reported in this work can provide a basic understanding for using the promotors in hydrate based applications.

Discussion
The guest does distort the hydrate lattice and this phenomenon is more dominant in case of large molecule of hydrate promotors. In this view, we estimated the cavity potential for these hydrates using a novel ab initio technique featuring spin component scaled Møller-Plesset perturbation theory. This feasible and accurate methodology for computing intermolecular interaction energy is used for estimating lattice stabilization energy by encapsulation of the guest molecules. The reference chemical potential difference estimates a mixed effect of lattice distortion by differently sized molecules of methane and promotors. A negligible contribution of methane in distorting the lattice as compared to the large molecules of promotors is proposed as a modification in the existing Gibbs free energy model that is previously designed for mixtures of comparable molecular sizes. Moreover, for the liquid promotors, the van der Waals volume of the guest shows an excellent correlation coefficient of 0.96 while relating to the estimated reference chemical potential difference. This lattice distortion theory grounds the formulation of a generalized model for phase equilibrium predictions for the promotor containing clathrate hydrates.
Methods fugacity of guest. For the empty to filled hydrate phase reaction, the effective concentration of the guest is needed. We address this quantity in terms of fugacity of the guest compound in vapor and liquid phases calculated by modified Patel-Teja equation of state having the following form, This is a three parameter equation in pressure (p) and temperature (T ). Here, R is a universal gas constant, and a, b and C have the following expressions, The fugacity of the component i in the vapor phase can be estimated as Activity of water. The activity of water is influenced by the promotor introduced for hydrate formation. For this purpose, the modified UNIFAC model is adopted that expresses the molar excess Gibbs free energy (G E ) as a summation of combinatorial (G C ) and residual (G R ) parts. For multicomponent systems, the activity coefficient can be expressed as, where, γ i denotes the activity coefficient of component i. The combinatorial (γ i C ) part of the γ i has the following formulation where, j is the false index for summing over all components present in the system. The term l i is expressed as a combination of the area (q i ) and volume parameters (r i ) for species i as Here, R and Q are the area and volume parameters for the functional groups with index k. The notation v k i is the number of functional groups present of type k in component i. The parameter Z denotes the coordination number of the system having a reasonable value of 10 and it is observed to have no substantial effect on activity calculation. The parameters θ i and Ω i are the molar weighed area and volume fractions that can be calculated as, The residual part of the activity coefficient accounts for the interaction among the groups present in the system. This is represented as, The term ψ mk accounts for the interaction between the unlike-groups present in the system. The modified UNIFAC model calculates the ψ mk as follows, The values of the parameters a mk (1) , a mk (2) and a mk (3) can be found in literature.

Data availability
The simulation datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

code availability
The code which was used in the findings of this study is available from the corresponding author upon reasonable request.