Direct transition mechanism for molecular diffusion in gas hydrates

Abstract

Several mechanisms for transport of gases inside hydrates have been proposed in the literature. Previous results have pointed out the apparent impossibility of some guests passing directly through faces connecting adjacent cages without destroying the water structure. Alternative diffusion mechanisms were proposed, such the participation of an additional guest (help-gas) or a water-hopping displacement of defects in the lattice. In this work, we use dual cage explicit atomic systems, where the position of some external atoms are fixed, to demonstrate theoretically that direct transitions are feasible through hexagonal and pentagonal faces in type I hydrate, for both carbon dioxide and methane, without compromising the overall structure integrity. Our DFT calculations show that even in the case that some bonds break during the transition, all of them are recovered, because the face distortion is absorbed locally by the hydrogen bond network. This result opens the door to improve the description of transport mechanisms of guest molecules by means of direct inter-cage transitions.

---

1 Introduction

Clathrate hydrates,¹ also known as gas hydrates, are inclusion solids formed by a crystal lattice of water molecules with regular crystalline voids which may contain small guest molecules, such as CO₂ or CH₄. Topics such as methane recovery from naturally existing hydrates,^2–4 the possibility of using them to capture anthropogenic CO₂ or their use in H₂ storage deposits are examples of the wide interest in these compounds.

Understanding the transport phenomena of guest molecules inside clathrate structures is fundamental to assess the practical possibility of storage, extraction or replacement of gases in hydrates.⁵ Transport phenomena depend on multiple properties and the characteristics of both the microscopic and macroscopic scales. For example, lattice defects, inhomogeneities, or the coexistence of different crystal structures must be considered. Among them, the transition energy barrier between adjacent crystal voids has a major role in the transport of guest molecules.

The type I hydrate structure is composed of two different cages, corresponding to the crystal void types inside the water lattice: a truncated hexagonal trapezohedron consisting of 12 pentagonal and 2 hexagonal faces (called a T cell, or alternatively 5¹² 6²), and a dodecahedron consisting of 12 pentagonal faces (a D cell, or 5¹²).

Computer simulation of clathrate hydrates is performed nowadays on a regular basis, constituting a research topic in itself. With about 30 years of history, a wealth of simulation methods and strategies have been employed so far, specially those describing periodic systems. A comprehensive review by English and MacElroy⁶ detailing the most important advances is available in the recent literature. Nevertheless, the understanding of the fundamental processes remains so far incomplete, and alternative approaches, further simulations extended in system size and time, and convenient structure models are necessary.

Previously, we have studied the structure and Raman spectra of isolated T and D cells of type I hydrate, with and without CO₂ and CH₄ as guest molecules.⁷ We also performed calculations of guest molecules passing through both types of faces present, considering either a rigid or a flexible lattice. The results showed that passing through an hexagonal face was allowed, but the transition through a pentagonal face led to the breaking of the cell. Several other approaches to this question can be found in the literature. Peters et al.⁸ proposed an interesting alternative, where the guest transport is supported by a water-hopping type mechanism. According to the authors, the transport process entails breaking completely the hydrogen bonds of a water molecule in one of the faces of a cell, while the elongation of the other bonds in the ring remains constant during the transition. Their molecular dynamics calculations show that this process yields a higher transition probability through hexagonal faces.

Other studies consider additional guest species (help-gas) that produce an enhancement of main guest diffusivity by mutual interaction through the water network, in hydrates with heavier guests.⁹ Another example of the enhancement or diminution of the guest–host hydrogen bonding by tuning the interactions was studied recently by Moudrakovski el al.¹⁰ for the guest pairs isobutane–CO₂ and THF–CO₂ in sII gas hydrate. In this direction, considerable work has been carried out recently concerning especially the diffusion of H₂ in type II hydrates. Due to the discrepancies between experimental and computed values for the diffusion constants, Trinh el al.¹¹ used Born–Oppenheimer molecular dynamics simulations with the BLYP exchange-correlation functional, and also Monte Carlo simulation employing in both cases a flexible network and considering the effect of cage occupancy to obtain the barriers of diffusing H₂ between cages. A low value of 5 kJ mol⁻¹ was reported for this transition at 100 K, but further reduction of this value as temperature increases is expected, eventually reaching the experimental value of 3 kJ mol⁻¹ at 250 K.

Direct mechanisms without defects or helping gases have been studied also, for example, by Alavi and Ripmeester.¹² They obtained an estimation of 0.250–0.283 eV for inter-cage transition energy barriers of H₂ in type II hydrate by using quantum mechanics methods (B3LYP and MP2 levels with 6-311++G(d,p) Pople basis) for rigid isolated cages. Román-Pérez et al. used an ab initio van der Waals density functional formalism to obtain diffusion activation energies of H₂, CO₂ and CH₄ in type H hydrate cells.¹³ They reported that CH₄ diffusion entails a substantial relaxation of the hydrate structure that can be supported only by hexagonal faces, while forcing the guest to pass through the pentagonal face destroys the hydrate structure.

Nevertheless, in spite of the results cited, the hydrogen bonds network in type I hydrate lattice allows one to think of an enhanced global structure flexibility. If one compares the estimated transition energy of CH₄ through an hexagonal face with the total energy of the hydrogen bonds of the molecules in the face, assuming as first order approximation that all bonds share symmetrically the distortion, the transition seems a priori possible. Therefore, we propose, after evaluating the results of our calculations, that type I hydrate lattice could effectively absorb the distortion imposed by the inter-cage transition of small guest molecules as CO₂ and CH₄. In the present work, that hypothesis is investigated from a theoretical perspective considering two neighbouring explicit cages using localized wave-functions, leading to an affirmative answer and also to several additional conclusions.

2 Computational methods

The systems studied were constructed by first optimizing geometric structures of isolated empty T and D cages by means of electronic Density Functional Theory (DFT),¹⁴ with the hydrogen bond connectivity corresponding to the minimum conformational energy.¹⁵ This particular hydrogen ordering is very improbable in real cages, but it is a better candidate as a reference structure, because the energy landscape depending on the H-bond connectivity is very flat, and consists of 304 3836 conformations for an isolated T cage. The difference between this structure and that immediately higher in energy¹⁶ is around 2.6 meV (0.25 kJ mol⁻¹) only.

The B3LYP/6-311+g(d,p) approximation, as implemented in Gaussian 09 was used for geometry optimizations.¹⁷ B3LYP stands for Becke¹⁸ three parameter Lee–Yang–Parr¹⁹ hybrid functional, and 6-311+g(d,p)^20,21 is the Pople-type basis set including diffuse and polarization functions. Taking into consideration the system size, and the minimum level of theory needed to obtain reliable structures and energies, the B3LYP hybrid functional was selected as the optimal procedure, as it is able to estimate accurately a number of gas hydrate properties. For example, Cappelli and Biczysko²² performed harmonic vibrational analysis on different systems using a wide variety of methods. Their study reveals that B3LYP results are the most reliable among them, combining also low computational cost. Pele et al.²³ emphasized the superiority of the B3LYP method for the calculation of Raman spectra when compared with other available alternatives. Also, our previous experience in calculating IR and Raman spectra of type I hydrates⁷ supports this choice.

Two inter-cage transitions were considered: from T to T through an hexagonal face (T6T) and from T to D through a pentagonal face (T5D). To calculate them, previously minimized T and D structures were merged to build two-cage systems: TT and TD, and re-optimized using the same level of theory (Fig. 1). Hydrogen-bonding plays a major role in the lattice structure stability, but also determines the lattice flexibility and its contribution is crucial in the host–guest interactions. Atoms in the central shared face in each double cage system, and in the following two layers, were set as fully relaxed, but the two farther atom layers were fixed during energy profile determination. This way, we guarantee that the water molecules involved in the transition at the point of maximum energy, the central face, have all their hydrogen bonds linked to neighbours, improving previous models based on isolated cages. Moreover, the first neighbours are relaxed, accounting for the local flexibility, but the second are not, reflecting the lattice large scale restraints.

Fig. 1 Explicit type I hydrate two-cage systems, TT and TD, considered in the present study. Cage-center to ring-center paths calculated using QTAIM theory are shown (dark green). Atoms in black were fixed during geometry optimizations.

Then, one guest molecule (CH₄ or CO₂) was added to the system (TT or TD) and the energy profile of inter-cage transition was calculated along the minimum energy path obtained using Quantum Theory of Atoms in Molecules (QTAIM).²⁴ QTAIM is based on the analysis of the topological features of the molecular electron density distribution ρ. The principal objects of molecular structures, such as atoms and bonds, can be naturally expressed using the characteristics of the system’s observable ρ. This function takes values that are particularly high around the nuclei, and decrease with the distance. Therefore, if we consider the boundaries determined by the planes of minimum ρ between nuclei in a given molecule, we can assign a definite spatial region to each nucleus, which provide us with an exact definition of the atom in the molecule. It will consist of the nucleus plus the corresponding electronic basin around, usually called Ω.²⁵ Nuclei are just one type of critical points (CP), points where the first derivative of ρ vanishes. CPs are classified according to two parameters, ω and σ, where the rank (ω) is the number of non-zero curvatures at the CP, and the signature (σ) represents the sum of the signs of the curvatures of the electronic density. In this case, over the four stable types of CPs having three non-zero eigenvalues, we are interested in the (3, −1) type. A (3, −1) point means that it has two negative curvatures and ρ presents a maximum in the plane defined by the corresponding eigenvectors, but also a minimum along the third axis, which is perpendicular to this latter plane. By analyzing the potential energy density V(r), the Lagrangian of the kinetic energy G(r), and also the Laplacian of the electron density at the bond (3, −1) CP, a classification of the bonds based on their nature^26,27 and the energy of the bond can be obtained.²⁸ This method has been used here to determine the cage-transition paths, locate the mentioned BCPs, and assess the bonding network during the processes examined. Then, an analysis of molecular structures and their geometrical variables was performed using Gabedit,²⁹ whereas topological properties were calculated with Multiwfn³⁰ over wave functions obtained using Gaussian 09.

3 Results

Series of points along the inter-cage minimum paths were considered for T6T and T5D transitions, and two different guest molecules, namely CH₄ and CO₂. For the case of CO₂, different orientations of the molecular axis with respect to the face plane were considered, and the perpendicular orientation resulted optimal. Inter-cage transition energies were directly obtained as the difference between the maximum and the minimum of the corresponding profile. In T5D, the profile is not cage-symmetrical, and therefore there are two different activation energies, T5D and D5T respectively, depending on the transition direction. This asymmetry favors in the case of CO₂ transitions towards T cages, whose occupancy ratio is important to ensure the hydrate structure mechanical stability, as shown also recently using classical molecular dynamics calculations.³¹ Fig. 2 and Table 1 summarize these results. An additional T5T transition could be considered as well, but the corresponding energy barrier is expected to be similar to that of T5D.

Fig. 2 Energy profiles of TT (hollow squares) and TD (solid circles) inter-cage transition for CO₂ (above) and CH₄ (below).

Table 1 Summary of the inter-cage transition energy barriers (E_a) applicable to the transport of CH₄ and CO₂ in type I hydrates

Guest	Transition	Ea/eV	Ea/kJ mol^−1
CH4	T6T	1.1625	112.16
	T5D	1.4119	136.23
	D5T	1.1553	111.47
CO2	T6T	0.5895	56.88
	T5D	1.1937	115.17
	D5T	0.9231	89.07

---

At first glance, as expected, CO₂ presents lower barriers than CH₄, specially in T6T transition, due to their relative size, shape, and CO₂ preferential orientation. It is noteworthy that CO₂@T6T yields a rather flat profile, corresponding to a small lattice distortion. Maximum energy values are of the order of eV, reaching 1.4 eV in CH₄@T5D transition. For CH₄, the T6T barrier is 1.16 eV, in good agreement with Román-Pérez,¹³ who reported 1.17 eV. The result obtained in that work for CO₂ is 0.42 eV, slightly lower than our value: 0.59 eV. We have also tested the system with other approximations traditionally used for hydrogen bond systems: PBE³² functional on the DEF2SV³³ fitting set and/6-311G+(2d,2pd) Pople basis set. Similar values on the interval of ±2% were obtained depending on the transition and guest molecule. The most important result of these tests is that the structure is stable in all cases during the guest transitions, in agreement with our starting hypothesis, but in contrast with previous estimations.^7,13 Assuming the network flexibility produces this difference, it makes feasible some direct transitions (e.g. T5D) that had been previously disregarded.

A direct estimation of the face deformation can be made by measuring the variation of the area of the polygon defined by the oxygen atoms. Table 2 summarizes the results obtained. As expected, deformations are considerably higher for CH₄ than for CO₂, in both T6T and T5D transitions. The minimum deformation corresponds to CO₂@T6T with less than 1% of area increment, and the maximum to CH₄@T5D, which duplicates the equilibrium area value. It is noteworthy that deformations are not completely symmetric over the face, but are related to the symmetry of the guest, so in CH₄@T6T, two groups of distances are observed between oxygen atoms in the ring, the longer corresponding to the positions of H atoms in CH₄ when passing through the face. Graphical representation of those deformations can be shown in Fig. 3 and 4.

At this point it is convenient to remember that all the calculations were done under the localized wavefunction approximation, which is the closest to the physical gas phase, but at the same time, and because of that, it has limitations for the study of lattice properties, due to the lack of a periodic explicit environment. This fact not being an issue for the geometrical stability or other similar theoretical studies, it would be of importance for the estimation of experimentally measured properties of real samples. This is particularly true for lattice energies or diffusion coefficients, which would need at least a calculation in a periodic monocrystal (see e.g. ref. 13) to compare with experiments. In this sense, the reader must be prevented from using the calculated energy values directly for any other purpose than those discussed in the preceding paragraphs.

Table 2 Ring area (A) of the face involved in the transition, and face area increment expressed in percentage

System	Guest	A/Å^2	ΔA (%)
TT	Empty	18.29
	CH4	26.87	46.95%
	CO2	18.41	0.66%
TD	Empty	12.14
	CH4	24.79	104.16%
	CO2	18.79	54.71%

---

Fig. 3 Graphical representation of the faces involved on the inter-cage migration of guest gases at the transition state. (A), (B) and (C) represent hexagonal common face of TT system in empty, CO₂ and CH₄ case, respectively. (D), (E) and (F) are the equivalent representations for the pentagonal face in TD system. The deformation caused by the guest transition is, as shown, not shared equally among all hydrogen bonds, specially in the case of CH₄ crossing the pentagonal face. Distances are expressed in Ångströms.

Fig. 4 Snapshots of the face configuration for CH₄ in TD (A) and TT (B) cases, where the face deformation process can be examined. The starting point in both cases is the minimum energy position of the guest on the T cell. The transition path followed by the guest is the minimum energy path computed with the QTAIM theory (see text for more information). The rest of water molecules that form the cell have been omitted for clarity. Distances are expressed in Ångströms.

QTAIM analysis reveals important changes in the ring bonding arrangement during the CH₄ transition (Fig. 5). Particularly, the longer elongations of the ring hydrogen bonds are accompanied by a redistribution of the corresponding BCPs: they migrate towards positions between H₂O and the guest molecule. This apparent guest–lattice stabilizing interactions are due to the structural confinement of the ring H₂O molecules, because the H₂O–guest interactions are not favourable. This result implies that the transition distortion is actually not equally shared among the hydrogen bonds in the ring, contrarily to what we supposed as a first approximation in this work. Performing the rigorous calculation reveals that the equipartition hypothesis was not exactly correct, although the structure proved to be able of assuming the distortion imposed by the guest molecule transition.

Fig. 5 Analysis of electronic density along the inter-cage ring planes in the transition states for CH₄ showing BCPs, paths connecting them, density isocontours and gradient lines. Above: T6T, below: T5D. Note that this QTAIM profiles correspond with the molecular structures (C) and (F) shown in Fig. 3.

Transitions of CO₂ produce a remarkably lower distortion in the network if compared with CH₄, as expected from the computed barrier energy values, when considering the guest in the optimal orientation (Fig. 6). In fact, transition through hexagonal face occurs without noticeable distortion of the hydrogen bond network. This fact is clearly illustrated in Fig. 7, where CO₂ is shown in the maximum energy point, surrounded by an annular symmetric region of weak non-covalent interaction with the lattice. The depicted isosurface corresponds to 0.5 of the non-dimensional reduced electronic density s expressed as

which has been reported to be a very useful tool to detect non-covalent interactions.³⁴

Fig. 6 Analysis of electronic density along the inter-cage ring planes for CO₂ transitions showing BCPs, paths connecting them, density isocontours and gradient lines. Above: T6T, below: T5D. Note that this QTAIM profiles corresponds with the molecular structures (B) and (E) shown in Fig. 3.

Fig. 7 TT system showing CO₂ passing through the hexagonal face in the point of maximum energy. The green annular region corresponds to the area of non-covalent weak interaction between guest and the water ring (oxygen atoms in black), as calculated with reduced density function (see main text for details). The blue box represents the boundaries of the region where the reduced density was computed.

During the transition through the pentagonal face, CO₂ molecule interacts significantly with H₂O lattice, inducing a particular deformation. In this case, one H₂O molecule is reoriented to favour hydrogen bonding interaction with the guest. This difference in behaviour can be attributed to the strong bond dipoles in CO₂ molecule, which are not present in CH₄. An estimation of the contribution can be obtained from the variation of binding energy of guest intra-molecular bonds (Table 3), being about 5 eV in favour of TD when compared with TT transition, where no reorientation was detected due the lack of distortion. As a guest–cage stabilizing contribution, its effect must be studied separately for a dynamic description of the system. For the purpose of the present work, though, it is not necessary, because the interaction is already present in the structures used for the energy evaluation of inter-cage barriers. Therefore, the resulting energy value is already taking into consideration this additional stabilizing contribution along the transition.

Table 3 Binding energies of the guest bonds computed using the QTAIM theory. The values of the binding energies of the guest bonds inside the minimum energy position, around the center of the T cage, are shown as a reference. T6T and T5D values correspond to the energy of the bonds at the transition states, the highest values of the profiles

Guest	System	EBond/eV	EBond/J mol^−1
CH4	T (reference)	−4.25	−410.06
	T6T	−4.26	−411.03
	T5D	−4.10	−395.59
CO2	T (reference)	−21.74	−2097.59
	T6T	−14.18	−1368.16
	T5D	−19.72	−1902.69

---

4 Conclusions

Inter-cage transition energies of CO₂ and CH₄ in type I hydrate were calculated by means of DFT at the B3LYP/6-311+g(d,p) theory level. Transition paths, bond critical points, and hydrogen bond features were determined through QTAIM.

Two-cage systems, TT and TD, were considered, fixing the position of the most external atoms, but relaxing the remaining during geometry optimizations. This way, structure flexibility is reproduced in the transition face, while global stability is attained by the fixed molecules in the more external layer.

Energy barriers found in double cages were lower than those corresponding to the rigid lattice approximation, and similar to those corresponding to relaxed single cages or isolated rings. Pentagonal faces were found to be stable under these conditions, and the lattice is then able to support the distortion due to both D5T and T5T transitions. This is in contradiction with previously reported calculations, based in flexible and isolated cells or rings, which concluded that the hydrate structure could not support the distortion, and pentagonal faces would not be stable.

As a final remark, the distance of each H-bond in the ring during inter-cage transitions depends on the relative orientation of the guest respect to the water lattice; and the total elongation is not equally shared among all ring bonds. These results indicate that the understanding of guest transport phenomena inside hydrates is not complete yet, and further studies in this direction seem to be necessary.

Acknowledgements

The authors acknowledge CESGA (http://www.cesga.es) in Santiago de Compostela, Spain, for providing access to computing facilities and Ministerio de Economía y Competitividad (FIS2012-33621 and FIS2015-68910-P, cofinanced with EU FEDER funds), Spain, for financial support.

References

1. E. D. Sloan and C. Koh, Clathrate Hydrates of Natural Gases, CRC Press, New York, 3rd edn, 2008.
2. A. K. Sum, C. A. Koh and E. D. Sloan, Ind. Eng. Chem. Res., 2009, 48, 7457–7465.
3. Natural Gas Hydrate in Oceanic and Permafrost Environments, ed. M. D. Max, Springer, 2011.
4. M. D. Max, A. H. Johnson and W. P. Dillon, Economic Geology of Natural Gas Hydrate, Springer, 2006.
5. G. Ersland, J. Husebø, A. Graue and B. Kvamme, Energy Proc., 2009, 1, 3477–3484.
6. N. J. English and J. M. D. MacElroy, Chem. Eng. Sci., 2015, 121, 133–156.
7. A. Vidal-Vidal, M. Pérez-Rodríguez, J.-P. Torré and M. M. Piñeiro, Phys. Chem. Chem. Phys., 2015, 17, 6963–6975.
8. B. Peters, N. E. R. Zimmermann, G. T. Beckham, J. W. Tester and B. L. Trout, J. Am. Chem. Soc., 2008, 130, 17342–17350.
9. S. Alavi and J. A. Ripmeester, J. Chem. Phys., 2012, 137, 054712.
10. I. L. Moudrakovski, K. A. Udachin, S. Alavi, R. I. Christopher and J. A. Ripmeester, J. Chem. Phys., 2015, 142, 074705.
11. T. T. Trinh, M. H. Waage, T. S. van Erp and S. Kjelstrup, Phys. Chem. Chem. Phys., 2015, 17, 13808–13812.
12. S. Alavi and J. Ripmeester, Angew. Chem., Int. Ed., 2007, 46, 6102–6105.
13. G. Román-Pérez, M. Moaied, J. M. Soler and F. Yndurain, Phys. Rev. Lett., 2010, 105, 145901.
14. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
15. J. D. Bernal and R. H. Fowler, J. Chem. Phys., 1933, 1, 515.
16. S. Yoo, M. V. Kirov and S. S. Xantheas, J. Am. Chem. Soc., 2009, 131, 7564–7566.
17. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09 Revision D.01, Gaussian Inc., Wallingford CT, 2009.
18. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
19. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785–789.
20. A. D. McLean and G. S. Chandler, J. Chem. Phys., 1980, 72, 5639–5648.
21. R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980, 72, 650–654.
22. C. Cappelli and M. Biczysko, in Computational Strategies for Spectroscopy, ed. V. Barone, Willey, 2012, Time-independent approach to vibrational spectroscopy, pp. 309–360.
23. I. Pele, J. Šebek, E. O. Potma and R. B. Gerber, Chem. Phys. Lett., 2011, 515, 7–12.
24. R. Bader, Chem. Rev., 1991, 91, 893.
25. C. F. Matta and R. J. Boyd, The Quantum Theory of Atoms in Molecules. From solid state to DNA and drug design, Wiley-VCH, Weinheim, 1st edn, 2007.
26. S. Emamian and S. Tayyari, J. Chem. Sci., 2013, 125, 939–948.
27. B. Shainyan, N. Chipanina, T. Aksamentova, L. Oznobikhina, G. Rosentveig and L. Rosentsveig, Tetrahedron, 2010, 66, 8551–8556.
28. E. M. E. Espinosa and C. Lecomte, J. Phys., Lett., 1998, 285, 170–173.
29. A.-R. Allouche, J. Comput. Chem., 2011, 32, 174–182.
30. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
31. J. M. Míguez, M. M. Conde, J.-P. Torré, F. J. Blas, M. M. Piñeiro and C. Vega, J. Chem. Phys., 2015, 142, 124505.
32. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
33. F. Weigend, F. Furche and R. Ahlrichs, J. Chem. Phys., 2003, 119, 12753–12762.
34. E. R. Johnson, S. Keinan, P. Mori-Sánchez, J. Contreras-García, A. J. Cohen and W. Yang, J. Am. Chem. Soc., 2010, 132, 6498–6506.

---