Multi-scale theoretical investigation of molecular hydrogen adsorption over graphene: coronene as a case study

Abstract

The physisorption of molecular hydrogen onto coronene is studied using a multi-scale theoretical approach with Density Functional Theory (DFT) calculations and Molecular Dynamics (MD) simulations. We consider two different kinds of model conformation for the approach of hydrogen towards the coronene i.e., systematic and random. For the systematic attack of hydrogen over coronene, the resulting potential energy profiles from DFT analysis are further found to resemble the Morse potential, and even the highly flexible Murrell–Sorbie (M–S) potential. The resulting M–S fitting also shows a zero-point energy correction of ∼16–17%. On the other hand, the potential energies from the random approach have been implemented into the Improved Lennard-Jones (ILJ) force field of the DL_POLY package following a prior statistical treatment. The MD simulations have been performed at different temperatures from 10 to 390 K. For the interaction of seven hydrogen molecules with coronene, the DFT method shows an average interaction energy of −3.85 kJ mol⁻¹ per H₂, which is slightly smaller than the Coupled Cluster value (CCSD(T)) of −4.71 kJ mol⁻¹ that was calculated for a single molecule in the most favorable situation. Moreover, the MD calculations reveal a mean interaction energy of −3.69 kJ mol⁻¹ per H₂ (a gross mean E_cfg of −25.98 kJ mol⁻¹ at T = 299.97 K), which is again in good agreement with the aforementioned DFT results, proving the quality of the approach used for the study of van der Waals interactions between hydrogen and graphene.

---

1. Introduction

Hydrogen energy has received significant attention as an alternative clean energy for the future post-petroleum age, since the reserves of fossil fuels are decreasing with growing consumer demand. Hydrogen is the most abundant element in the universe, being found mainly in stars. Although molecular hydrogen does not occur on Earth in significant amounts, it can be extracted from water by decomposition using chemical, thermal or electrical energy and from renewable resources such as biomass or coal by consecutive natural processes. In this context, the use of hydrogen makes its storage a crucial issue along with other branches of hydrogen technology progress^1,2 and therefore, a considerable amount of work has been dedicated to the search for efficient and technology-friendly hydrogen storage materials. Several materials have significant hydrogen storage capacities, for instance porous materials, different kinds of hydride and carbonaceous nanomaterials.³ In particular, the potential use of graphene for gas sensing has attracted the attention of the scientific community in recent years,^4–9 being considered a good candidate for gas storage.^1,2

Graphene, with its honeycomb two-dimensional structure of sp² carbon atoms, is without doubt one of the most important materials to be investigated^10–12 because of its robustness, flexibility and enormous surface area. Furthermore, the capability of graphene to adsorb different gases is also important from an environmental point of view.¹³ It could be used to store CH₄ (ref. 14 and 15) or CO₂ (ref. 13 and 16) in order to prevent global warming caused by greenhouse gases. However, so far, the storage capacity obtained at room temperature and pressure is still far from the target of the US Department of Energy (DoE), which is 9.0 wt% gravimetric density and 81 g L⁻¹ volumetric capacity by 2015.¹⁷ Nevertheless, from a technological point of view, graphene provides an ideal environment to investigate catalytic processes, surface diffusion and molecular hydrogen formation from chemisorbed atomic hydrogen.¹⁸

The interaction of atomic and molecular hydrogen with graphene has been investigated from theoretical and experimental points of view.¹⁹ The physisorption of molecular hydrogen on models of planar and curved graphene has been reported by Okamoto et al.,²⁰ Vilaplana,²¹ and Patchkovskii et al.,²² among others. The main conclusions of these works are the importance of the availability of accurate C–H₂ potentials and the fact that among the three possible adsorption sites (hollow, bridge and on-top sites), the interaction is most attractive at the hollow sites. MD studies on the adsorption and desorption behavior of hydrogen on graphite have also been reported recently.^16,23

Today, multiscale modeling and simulation (MMS) is considered as a prominent methodological approach in computational science and engineering.²⁴ The usual strategy of MMS is to apply the more accurate and expensive results from first-principle Density Functional Theory (DFT) calculations to the atomic and/or molecular scale force fields (inexpensive but less accurate) of classical simulations, such as Molecular Dynamics (MD), Monte Carlo (MC), Molecular Mechanics (MM) etc., for a particular model chemical system. For a gas–solid system, the first-principle calculations investigate the gas adsorption behavior, including specific adsorption sites, binding energies, and interaction mechanisms. As a counterpart to this, MD simulations explore the dynamic properties of the system, which subsequently produce the simulated macroscopic gas adsorption properties for the modeled material in order to complement experimental work.

In recent years, a modification of the Lennard-Jones potential energy has been proposed, the Improved Lennard-Jones (ILJ),^25,26 which is useful for describing both bonded²⁷ and nonbonded (see for instance ref. 28 and 29) non-electrostatic interaction contributions. The model contains only one additional parameter with respect to the LJ function but it is sufficient to add flexibility and to improve the description of the interaction at both long and short range simultaneously. The parameters of the interaction can be obtained by assuming the additivity of the bond polarizabilities, but in the present investigation we are interested in using the ILJ model to fit the DFT computation.

Our purpose in this investigation is to carry out a theoretical study on the adsorption of molecular hydrogen on coronene (C₂₄H₁₂), which can be considered a good prototype for graphene. To this end, the interaction of coronene with up to seven hydrogen molecules in different geometrical conformations has been analyzed at the DFT level of theory. After this, the obtained results were fitted to an analytical function and the adsorption was investigated by means of Molecular Dynamics (MD) simulations.

The paper is structured as follows: in section 2, the results of our first principles calculations are collected. In section 3, we outline the formulation of the semiempirical PES (potential energy surface), and in section 4, the results of MD simulations are presented. The most important conclusions are given in section 5.

2. First principles calculations

The geometry of coronene was optimized at the B3LYP/6-31G** level^30,31 of DFT on the Gaussian09 package.³² Then, the physisorption of molecular hydrogen onto C₂₄H₁₂ was studied by means of a supramolecular approach correcting BSSE by the counterpoise procedure.³³ The interaction energy of this closed-shell system was computed using the B97-D functional³⁴ and split-valence triple-zeta basis supplemented with two polarization functions TZV2P.³⁵ The B97-D functional recovers the attractive dispersion forces through its empirical dispersion term whereas the use of TZV-type basis sets has also been found sufficient to represent the polarization effects in similar systems.³⁶

As a first step, the systematic approach of n molecules of H₂ parallel to the C₆ symmetry axis of coronene was investigated, with n ranging from 1 to 7. In all cases, one H₂ molecule attacked the center of the coronene, while the remaining ones attacked the center of an outer benzenic hexagon, and when several conformations were allowed, all of them were considered and subsequently averaged. The obtained results are depicted in Fig. 1, where the total interaction energy is represented vs. distance for the different number of hydrogen molecules (denoted as #nH₂) approaching perpendicularly to the coronene.

Fig. 1 Mapping of potential energy curves for the adsorption of a different number of molecules of H₂ approaching coronene parallel to the C₆ symmetry axis. Distances are measured with respect to the hydrogen atom of the H₂ molecule closest to the coronene surface.

In order to better interpret these profiles, we find it convenient to fit the computed values of the interaction energy per hydrogen molecule ΔE to a Murrell–Sorbie³⁷ potential of fourth degree

ΔE = D_e(1 + a₁(r − r_e) + a₂(r − r_e)² + a₃(r − r_e)³ + a₄(r − r_e)⁴)exp(−a₁(r − r_e)) (1)

which, because of its greater flexibility, generally provides better results than the well-known Morse potential

ΔE = D_e(1 − exp(−a(r − r_e)))² (2)

From the second derivative at r = r_e of the potential, the constant of the C₂₄H₁₂–H₂ van der Waals bond, k, was computed and from it the true dissociation energy, D₀ = D_e − E₀, where E₀ is the zero point energy, was obtained. The determined parameters are collected in Table 1.

Table 1 Parameters (see text for details) from the Murrell–Sorbie fitting of the average energy per hydrogen of different model adsorption systems in the systematic adsorption of H₂ over coronene. r_e refers to the distance between the center of mass of H₂ and the coronene surface

#H2	re/Å	De/kJ mol^−1	k/N m^−1	E0/kJ mol^−1	D0/kJ mol^−1
1	3.092	4.41	1.77	0.73	3.68
2	3.085	4.15	1.47	0.66	3.49
3	3.079	4.06	1.43	0.65	3.41
4	3.074	4.02	1.42	0.65	3.36
5	3.065	3.98	1.39	0.64	3.33
6	3.057	3.96	1.38	0.64	3.31
7	3.050	3.94	1.37	0.64	3.30

---

The analysis of Table 1 clearly reveals how the addition of extra hydrogens is weakening the average bond between each hydrogen and the coronene surface. This phenomenon is observed in all the computed parameters, since the different energies and the constant of strength of the bond consistently decrease while the bond length increases when the number of hydrogen molecules is increased. Remarkably, these effects are especially important after adding the second hydrogen molecule. Nevertheless, with the only exception being the adsorption of a single hydrogen molecule, the zero point energy E₀ remains almost constant. Moreover, it is worth stressing that both D_e and D₀ approach a constant value, since the slope of the curve (see Fig. 2) is obviously diminishing.

Fig. 2 D₀ (red) and D_e (blue) dissociation energies per molecule as a function of the number of H₂ molecules.

To confirm our previous results, we have also fitted our data to a Morse potential,³⁸ obtaining almost identical results to those reported using the Murrel–Sorbie potential, again with the only exception being the adsorption of a single hydrogen molecule.

At any rate, the parallel conformations considered until now – although useful for a better understanding – are evidently too idealized, as they are extremely improbable. Therefore, it is appropriate to consider more realistic conformations of the system. To this end, we have randomly generated 99 conformations for each distance between the center of mass of the hydrogen molecules and the coronene surface, and averaged them. In this way, the effects of the orientation and position of the H₂ molecules are accounted for. As one can see in Fig. 3, the interaction energy is lowered while the bond distance is increased as expected. These facts are also evident in the data presented in Table 2. Furthermore, for the case of seven molecules of hydrogen a diminution of the harmonic frequency, ν_e, can be observed, reflecting the diminution of the bond strength. However, the first anharmonic correction, x_eν_e, remains almost unaltered in all cases.

Fig. 3 Comparison of the parallel (continuous lines) and random (dashed lines) approaches of one (red) and seven (blue) hydrogen molecules to coronene.

Table 2 Spectroscopic parameters for the adsorption of 1 or 7 H₂ molecules over coronene in parallel and random approaches. r_e indicates the distance between the center of mass of H₂ and the coronene surface

Approach	re/Å	De/kJ mol^−1	νe/cm^−1	xeνe/cm^−1
Parallel_1	3.09	368.6	122.2	0.58
Parallel_7	3.05	329.2	107.4	0.57
Random_1	3.57	318.1	129.6	0.64
Random_7	3.59	271.4	99.4	0.59

---

In section 3, the DFT results are utilized to fit a different potential well suited for MD simulations, as using DFT data should improve the performance of such simulations by providing a better description of the interaction both at short and long distances.

3. The potential energy surface

The potential model assumes that the total interaction energy of coronene and molecular hydrogen, V_total, can be expressed as a sum of interactions involving the hydrogen molecules and the carbon atoms of coronene. Accordingly,

Any V_C–H2 and V_H2–H2 term is expressed by means of the ILJ function^25,26 as

withwhere β is an adjustable parameter related to the hardness of the interacting partners, which adds flexibility to the potential energy function in comparison with the Lennard-Jones function. The variable r represents the distance from the C atom to the center of mass of the hydrogen molecules for the V_C–H2 interactions or the distance between the center of mass of two hydrogen molecules for the V_H2–H2 interactions. The m parameter is taken to be equal to 6, the typical value for neutral–neutral interactions. The well depth ε and the equilibrium distance r₀ have been derived from the DFT results for different values of the β parameter. The advantages of using the ILJ function instead of the LJ one have been well-explained elsewhere.²⁵ In order to obtain the best parameters of the potential function for V_C–H2, three statistical treatments are considered:

(1). Random averaged (RA): the 99 random conformations of one hydrogen molecule and the coronene are averaged to a H₂ molecule moving along the C₆ axis of coronene. Therefore, the graphene–hydrogen interaction, , reduces because of symmetry reasons to the simple form 6V_c + 6V_m + 12V_f, where V_c, V_m and V_f represent the interactions with the three symmetric types of carbon atom (close, medium and far, respectively).

(2). Random not averaged (RNA): there is no averaging of the configurations of adsorption interactions. So, for each of the 24 carbon atoms of the coronene structure, there are a total of 2178 C–H atom–atom interactions, resulting from 99 interactions at 22 distances, fitted to the general formula of the ILJ potential according to eqn (4).

(3). Random not averaged selected (RNA-S): the fitting is done as in the case above, but with the difference that, in this case, we only include the points near the average value for each distance i.e., those in the interval x̄ ± σ.

As the final goal of these fittings is to study the interaction of H₂ with graphene, the dangling H atoms of coronene have not been included in any of the fittings. The obtained parameters of the potential defining the interaction between the center of mass of one hydrogen molecule and one C atom of coronene are given in the first part of Table 3.

Table 3 The well depth (ε), the equilibrium distances (r₀) and the β parameter for the C–H₂ and H₂–H₂ interactions obtained by fitting the DFT results as discussed in the text

Averaged data set	ε/kJ mol^−1	r0/Å	β
C–H2
RA	0.289	3.987	8.824
RNA	0.381	3.519	7.295
RNA-S	0.373	3.561	6.777
H2–H2
RNA	0.406	3.540	6.5

---

The H₂–H₂ interaction has also been calculated according to the RNA scheme, representing each hydrogen molecule only by its center of mass. The optimized ILJ parameters are given in the last line of Table 3. Finally, the fitted ILJ function, V_C24H12,H2, as obtained from the RNA data set of the DFT computation, is represented in Fig. 4.

Fig. 4 Fitted ILJ potential (blue curve) with the averaged interaction energies (red points) for the coronene–(H₂) system. Moreover, the green points, which accumulate to give a bar or a large point, evaluate the closeness to normal distribution.

4. Molecular dynamics simulations

Bearing in mind a previous study on the adsorption dynamics of molecular hydrogen on graphene, and in order to test the potential parameters, MD simulations of the coronene–(H₂)₇ prototype system have been performed, using the DL_POLY program,³⁹ and considering the microcanonical ensemble (NVE) of particles.

The coronene compound has been kept rigid in all calculations. A time step of 1 fs has been used in all the simulations presented here, integrating the MD trajectories up to a final time of 3.5 ns. The time step chosen is large enough to keep the fluctuations of E_total well below 10⁻⁵ kJ mol⁻¹. The MD simulations are preceded by an initial equilibration period of 0.5 ns (using the same time step) during which the velocities of all atoms are rescaled to match the input temperature. The equilibration period has been excluded from the statistical analysis at the end of any trajectory.

In order to decide which set of potential parameters – derived from the three different statistical basis sets considered – is the best one, preliminary MD simulations were performed using the three sets of parameters. It was found that the mean values of the potential energy attained during the 3.5 ns of simulation, represented by E_cfg (configurational energy) are somewhat different. At a mean temperature of 299.97 K it was observed that the value of E_cfg for β = 7.295 (RNA data set) is close to that of β = 6.777 (RNA-S data set). Moreover, it was found that a better convergence of the total energy is attained when using the potential parameters associated with the RNA data set. Accordingly, we have chosen the parameters given in the third row of Table 3 to calculate the C–H₂ (see Fig. 4) interaction and to study the adsorption of molecular hydrogen on coronene.

The quality of the potential parameters was tested by performing MD simulations at low temperatures of the coronene with two molecules of H₂ adsorbed. The MD results at the lowest temperature value investigated are compared with the counterpoise-corrected CCSD(T) results as shown in Table 4. As can be seen, the MD results are in good agreement with the CCSD(T) ones. The two methods give very similar values of E_cfg (ΔE_cfg = 0.97 kJ mol⁻¹). This indicates that the potential parameters are adequate to investigate the adsorption of H₂ on coronene.

Table 4 Lowest distance from the center of mass of the adsorbed molecules to the coronene surface and configurational energy (E_cfg) for the coronene–(H₂)₂ system

Type of calculation	Distances/Å	Ecfg/kJ mol^−1
CCSD(T)	3.020	−9.43
	3.020
MD (T = 9.99 K)	3.054	−8.46
	3.046

---

The adsorption of seven molecules of H₂ on coronene was investigated by increasing the temperature from 10 K to 390 K. Results of NVE simulations indicate that the configuration energy, E_cfg, remains almost constant (E_cfg ≈ −26 kJ mol⁻¹), which indicates the stability of the system with the temperature change in the mentioned range. Accordingly, the increase in the kinetic energy (due to the increase of T) causes a decrease in the total energy. In order to analyze whether any adsorption takes place, we have adopted the same criteria found from the Murrell–Sorbie and Morse potential fittings of the DFT calculation, which establish that adsorption occurs at distances lower than 3.2 Å from the hydrogen molecules to the C atoms of coronene. Previously, the radial distribution function, RDF, between the center of mass of H₂ and the C atoms has also been analyzed (see Fig. 5). From the figure it seems that the distance criterion mentioned above is adequate for MD simulations.

Fig. 5 The RDF, g(C–X), between a C atom and the center of mass of a H₂ molecule, defined by the symbol X, as a function of distance, r, for RNA (β = 7.295) at T = 299.97 K.

The results of the MD simulations show that the number of adsorbed molecules changes, depending on the initial configuration and temperature. At temperatures higher than 250 K, a total of 3–4 H₂ molecules are adsorbed. Table 5 gives the adsorption results at a mean temperature close to 300 K, which are compared with the DFT results. As can be seen, at T = 299.97 K, E_cfg = −25.98 kJ mol⁻¹, which is in good agreement with the results obtained from the DFT approach, showing an average interaction energy of −3.85 kJ mol⁻¹ (a total of −26.95 kJ mol⁻¹ for corenene interacting with seven H₂ molecules), which is again reasonably close to the Coupled Cluster value (CCSD(T)) of −4.71 kJ mol⁻¹ reported in Table 4. A snapshot of the final configuration of the coronene–(H₂)₇ system at T = 299.97 K is visualized in Fig. 6.

Table 5 Adsorption of H₂ on coronene from DFT calculations and MD simulations

Type of calculation	#H2 adsorbed	Distances/Å	Ecfg/kJ mol^−1
DFT approach	7	3.074	−26.95
MD (T = 299.97 K)	4/7	3.065	−25.98
		3.199
		3.179
		3.144

---

Fig. 6 Snapshot of the last configuration of the simulation for the coronene–(H₂)₇ system at T = 299.97 K. The random white-brown dumbbells indicate the seven H₂ molecules over the white-green coronene surface.

5. Conclusions

Physisorption of molecular hydrogen over graphene is very delicate because of the weak van der Waals interactions between the H₂ and the graphene surface. With the computational tools available to model such a system it is still challenging to achieve an accurate and realistic interaction energy. This interaction has been modeled by means of a supramolecular approach for the interaction of coronene with up to seven hydrogen molecules. The interaction energy, both at DFT and MD level, was evaluated taking all the above facts into account. The force field (FF) fitting to the simple graphene system was imposed in a reasonably acceptable manner considering all the possible types of interactions involved in this adsorption process. For the interaction of seven hydrogen molecules with coronene, the counterpoise corrected B97-D approach shows an average interaction energy of −3.85 kJ mol⁻¹ per H₂ (a total of −26.95 kJ mol⁻¹ for seven H₂ molecules), which is reasonably close to the CCSD(T) value of −4.71 kJ mol⁻¹ reported by our group. Moreover, the MD calculation reveals a mean interaction energy of −3.70 kJ mol⁻¹ per H₂ (a gross mean E_cfg of −25.98 kJ mol⁻¹ at T = 299.97 K), which once again is in good agreement with the corresponding results mentioned above.

Acknowledgements

Part of this work was supported by the Spanish Ministerio de Economía y Competividad through project CTQ2010-19738. M. Albertí acknowledges financial support from the Ministerio de Educación y Ciencia (Spain, Project CTQ2013-41307-P) and the Generalitat de Catalunya (2014 SGR 25). Thanks are also due to the Center de Supercomputació de Catalunya CESCA-C4 and Fundació Catalana per a la Recerca for the allocated supercomputing time. N. Faginas Lago acknowledges financial support from MIUR PRIN 2008 (contract 2008KJX4SN 003), Phys4entry FP72007-2013 (contract 242311) and EGI Inspire. Thanks are also due to INSTM, CINECA, IGI and the COMPCHEM virtual organization for the allocation of computing time. The research leading to part of the results presented in this paper has been made possible thanks to the Grid resources and services provided by the European Grid Infrastructure (EGI) and the Italian Grid Infrastructure (IGI). For more information, please refer to the EGI-InSPIRE paper (http://go.egi.eu/pdnon) and the IGI web server (http://www.italiangrid.it).

References

1. C. Liu, Y. Y. Fan, M. Liu, H. Cong, H. Cheng and M. Dresselhans, Science, 1999, 286, 1127–1129.
2. Q. Wang and K. Johnson, J. Chem. Phys., 1999, 110, 577–586.
3. M. Hirscher, Handbook of Hydrogen Storage: New Materials for Future Energy Storage, Wiley-VCH, 2010.
4. F. Shedin, A. Geim, S. Morozov, E. Hill, P. Blake, M. Katsnelson and S. Novoselov, Nat. Mater., 2007, 6, 652–655.
5. O. Leenaerts, B. Partoens and F. Peeters, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 125416.
6. H. Yoon, D. Jun, J. Yang, Z. Zhou and S. M.-C. C. M. Yang, Sens. Actuators, B, 2011, 157, 310–313.
7. G. Lu, L. Ocola and J. Chen, Nanotechnology, 2009, 20, 445502.
8. K. Batinac, W. Yang, S. Ringer and F. Braet, Environ. Sci. Technol., 2010, 44, 1167–1176.
9. G. Ko, H.-Y. Kim, J. Ahn, Y.-M. Park, K. Lee and J. Kim, Curr. Appl. Phys., 2010, 10, 1002–1004.
10. K. Novosolov, A. Geim, S. Morozov, Y. Jiang, D. Zhang, S. Dubonos, I. Grigorieva and A. Firsov, Science, 2004, 306, 666–669.
11. A. Geim and K. Novoselov, Nat. Mater., 2007, 6, 183–191.
12. A. Castro Neto, F. Guinea and N. Peres, Rev. Mod. Phys., 2009, 81, 109–162.
13. K. Kemp, H. Seema, M. Saleh, N. Le, K. Mahesh, V. Chandra and K. Kim, Nanoscale, 2013, 5, 3149–3171.
14. C. Thierfelder, M. Witte, S. Blankenburg, E. Rauls and W. Schmidt, Surf. Sci., 2011, 605, 746–749.
15. N. Faginas-Lago, M. Albert, A. Lagan, A. Lombardi, L. Pacifici and A. Costantini, Computational Science and Its Applications ICCSA 2014, Springer International Publishing, 2014, vol. 8579, pp. 585–600.
16. T. Trinh, J. Vlugt, M. Hägg, D. Debeaux and S. Kjelstrup, 12th Joint European Thermodynamics Conference, 2013, pp. 514–518.
17. Graphene Chemistry: Theoretical Perspectives, ed. D.-E. Jiang and Z. Chen, Springer-Verlag, John Wiley & Sons, 2013.
18. J. Petucci, C. LeBlond, M. Karimi and G. Vidali, J. Chem. Phys., 2013, 139, 044706.
19. G. Vidali, J. Low Temp. Phys., 2013, 170, 1–30.
20. Y. Okamoto and Y. Miyamoto, J. Phys. Chem. B, 2001, 105, 3470.
21. A. Ferré-Vilaplana, J. Chem. Phys., 2005, 122, 104709.
22. S. Patchkovskii, J. S. Tse, S. N. Yurchenko, L. Zhechkov, T. Heine and G. Seifert, PNAS, 2005, 102, 10439.
23. O.-E. H. J.-M. Simon and S. Kjelstrup, J. Phys. Chem. C, 2010, 114, 10212.
24. Lecture Notes in Computational Science and Engineering: Multiscale Modeling and Simulation, ed. S. Attinger and P. Koumoutsakos, Springer-Verlag, Berlin, 2004.
25. P. Pirani, S. Brizi, L. Roncaratti, P. Casavecchia, D. Cappelletti and F. Vecchiocattivi, Phys. Chem. Chem. Phys., 2008, 10, 5489.
26. N. Faginas-Lago, F. Huarte Larrañaga and M. Albertí, Eur. Phys. J. D, 2009, 55, 75.
27. M. Albertí, A. Costantini, A. Laganà and F. Pirani, J. Phys. Chem. B, 2012, 116, 4220–4227.
28. M. Albertí, A. Aguilar, J. M. Lucas and F. Pirani, Theor. Chem. Acc., 2009, 123, 21–27.
29. M. Albertí, A. Aguilar, D. Cappelletti, A. Laganà and F. Pirani, Int. J. Mass Spectrom., 2009, 280, 50–56.
30. A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
31. R. D. W. J. Hehre and J. Pople, J. Chem. Phys., 1972, 56, 2257.
32. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, V. B. G. Scalmani, 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, J. J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, 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. Jaramilloa, 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, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian09 Revision D.01, Gaussian Inc., Wallingford CT, 2009.
33. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553.
34. S. Grimme, J. Comput. Chem., 2006, 27, 1787.
35. H. A. Schaefer and R. Ahlrichs, J. Chem. Phys., 1994, 100, 5829.
36. F. Tran, J. Weber, T. Wesolowski, F. Cheikh, Y. Ellinger and F. Pauzat, J. Phys. Chem. B, 2002, 106, 8689.
37. J. Murrell and K. Sorbie, J. Chem. Soc., Faraday Trans. 2, 1974, 70, 1552.
38. P. M. Morse, Phys. Rev., 1929, 34, 57–64.
39. W. Smith, C. Yong and P. Rodger, Mol. Simul., 2002, 28, 385–471.

---

Footnote

† These authors contributed equally to this work.

---