Spectroscopy of a rotating hydrogen molecule in carbon nanotubes

María Pilar de Lara-Castells *a and Alexander O. Mitrushchenkov *b
aInstituto de Física Fundamental (C.S.I.C.), Serrano 123, E-28006, Madrid, Spain. E-mail: Pilar.deLara.Castells@csic.es
bLaboratoire Modélisation et Simulation Multi Echelle, MSME UMR 8208 CNRS, Université Paris-Est, 5 bd Descartes, 77454 Marne-la-Vallée, France. E-mail: Alexander.Mitrushchenkov@u-pem.fr

Received 28th June 2018 , Accepted 3rd August 2018

First published on 6th August 2018

A first-principles study of the spectroscopy of a single hydrogen molecule rotating inside and outside of carbon nanotubes is presented. Density functional theory (DFT)-based symmetry-adapted perturbation theory (SAPT) is applied to analyze the influence of the rotation in the dispersionless and dispersion energy contributions to the adsorbate–nanotube interaction. A potential model for the H2–nanotube interaction is proposed and applied to derive the molecular energy levels of the rotating hydrogen molecule. The SAPT-based analysis shows that a subtle balance between the dispersionless and dispersion contributions is key in determining the angular dependence of the H2–nanotube interaction, which is strongly influenced by the diameter of the carbon nanotubes. As a consequence, the structure of molecular energy levels is very different in wide and narrow nanotubes with the diameter above and below 1 nanometer, respectively. Strong anisotropy effects lead to a rather constrained rotation of molecular hydrogen inside narrow nanotubes.

1 Introduction

High surface area and precisely tuned pores of carbon nanotubes (CNTs) make them relevant materials for applications such as gas adsorption, selective separation of light isotopes,1 efficient hydrogen storage, and nanoreactivity, e.g., to synthesize and confine metal nanoparticles in quasi one-dimensional configurations.2 The general goal of hydrogen storage methods being to pack hydrogen as close as possible, the existing remarkable experiments provide evidence for solid-like packing of hydrogen at supercritical temperatures of industrial importance3 as well as for high-density orientationally ordered phases at 75 K.4 On a theoretical side, a very recent study showed a possibility for hexagonal-close-packing of deuterium molecules in carbon nanotubes.5

The crucial role of quantum nuclear effects in the motion of light molecular species inside carbon nanotubes has been demonstrated in experimental observations at low temperatures. Thus, a recent work by Ohba6 revealed a quenched adsorption of helium atoms in carbon nanopores with a diameter below 0.7 nm at very low temperatures of 2–5 K. The same work also showed that more molecules of nitrogen than helium atoms adsorb in these narrow nanotubes despite the larger kinetic diameter of molecular nitrogen.6 A subsequent helium density-functional theory (DFT)-based study of carbon nanotubes immersed in helium nanodroplets pointed out that the quenched adsorption can be attributed to very large zero-point energies in the motion of 4He atoms, including the formation of cavities with zero helium densities.7 A follow-up study revealed that the experimental observation can also be explained as a result of the extreme one-dimensional confinement existing along the axis of narrow nanotubes for 4He atoms and N2 molecules.8 In the particular case of molecular hydrogen inside carbon nanotubes, a quantum-induced reverse trend in H2 and D2 diffusion rates upon lowering the temperature was experimentally detected9 and theoretically confirmed.10 The impact of quantum nuclear effects has also been highlighted when characterizing the (superfluid or crystalline) phase of parahydrogen molecules inside carbon nanotubes at zero11,12 and ultra-low temperatures.13

Recent studies using high-resolution neutron spectroscopy have allowed accurate characterization of the motion of molecular hydrogen under confinement conditions provided by carbon-based nanoporous materials such as carbon nanohorns14 and aligned carbon nanotubes.15 A key ingredient when aiming at reproducing and predicting experimental results is a realistic model for the molecule–nanotube interaction. Ab initio intermolecular interaction theory allows an accurate description of van der Waals (vdW)-dominated interactions of molecules with carbon nanostructures.7,8,16–23 Recent options are mixed schemes combining second-order Möller–Plesset perturbation theory with the coupled-cluster approach,16,17 as well as non-periodic implementations of the incremental method24 applied at the coupled-cluster level25–29 with periodic calculations using dispersionless density functional theory.28,30–32 An alternative methodological protocol8 combines DFT-based symmetry-adapted perturbation theory (SAPT), which is used to derive the parameters of potential models describing gas adsorption to the carbon material, with an ad hoc-developed adsorbate wave function treatment. This method has very recently been applied in the precise characterization of the nature of bound states of atomic helium, molecular nitrogen, and molecular deuterium clusters inside carbon nanotubes.5,8

In previous works5,8 molecules of nitrogen and deuterium were treated as structureless point-like particles using the full configuration interaction Nuclear Orbital (FCI-NO) approach (see, e.g., ref. 33–35 and references therein). As a step forward, we propose here an extension of the adsorbate wave function treatment allowing explicit inclusion of the rotational mode of molecular hydrogen. This extension is also motivated by very recent experimental evidence demonstrating the appearance of a high-density orientationally ordered phase in H2 adsorbed in subnanometer pores.4 This way, our focus will be on analyzing the effect of including the H2 orientational mode upon increasing the carbon tube diameter from the subnanometer scale. The next section is dedicated to a brief description of the computational methods. The results are reported in Section 3. Finally, the main conclusions are provided in Section 4.

2 Method

2.1 SAPT-based decomposition of the adsorbate–nanotube interaction

The SAPT(DFT) method36,37 has been used to decompose the interaction energy for the H2/CNT system as a function of the angular orientation of the molecular adsorbate. The total interaction energy is written as a sum of first- and second-order interaction terms, namely first-order electrostatic Eelec and exchange Eexch, second-order induction Eind and dispersion Edisp terms, along with their respective exchange corrections (Eexch–ind and Eexch–disp). The dispersion and exchange–dispersion SAPT(DFT) contributions were summed to provide the total exchange–repulsion part, as done in the case of the induction and exchange–induction contributions to provide the total induction term.
2.1.1 Computational details. A carbon nanotube (CNT) is usually characterized by the so called helicity index (m,n) (see, e.g., ref. 38). Considering a (5,5) CNT as an example of a tube with subnanometer diameter, we have applied the SAPT(DFT) method as implemented within the MOLPRO electronic structure package39 including density fitting (DF).40 The dangling bonds were saturated with hydrogen atoms. We followed the same computational setup reported in previous works,5,8 using the Perdew–Burke–Ernzerhof (PBE) density functional,41 and a modified version of the augmented polarized correlation-consistent triple-zeta basis42 (aug-cc-pVTZ) for nanotube carbon atoms as well as for adsorbate hydrogen atoms. The DF of Coulomb and exchange integrals employs the auxiliary basis set developed for the aug-cc-pVTZ basis by Weigend,43 while the aug-cc-pVTZ/MP2Fit basis44 was used for fits of integrals containing virtual orbitals. The exchange–correlation PBE potential was asymptotically corrected45 using the ionization potential (IP) value reported in the NIST Chemistry Web Book for the H2 molecule46 while the IP value for the nanotube was estimated using the DFT PBE0 approach of Adamo and Barone47 in ref. 8. In total, the system consists of 64 atoms, comprising an orbital basis of 2466 basis functions and an effective basis size of 4712 functions when including density fitting (see the ESI of ref. 5 for further details).

2.2 Adsorbate wave-function approach: Hamiltonian model

The full motion of a H2 molecule inside a carbon nanotube is a six dimensional (6D) problem. In this section, we describe a series of approximations reducing the dimensionality of the problem to three degrees of freedom only.
2.2.1 Potential model of the adsorbate–nanotube interaction. As a first approximation, let us assume that the interaction of the H2 molecule with the carbon nanotube can be written as a sum of the interactions of the individual H atoms with the nanotube. Explicitly,
VH2–CNT(ρ1,ρ2) = VH–CNT(ρ1) + VH–CNT(ρ2)
where ρ1 and ρ2 stand for the radial distance of each H atom to the tube axis. Notice that an isotropic nature of the H–CNT interaction is also being assumed. Specifically, the very small corrugation of the VH–CNT potential is ignored, and the potential is considered to be almost constant along the tube axis, decaying very slowly at its borders. Hence, using cylindrical coordinates, we assume VH–CNT(ρ,z,ϕ) ≈ VH–CNT(ρ). The validity of this approximation has been assessed numerically. Using this model, the interaction VH–CNT(ρ) can be easily calculated from the known H2–CNT interaction for a parallel orientation with respect to the tube axis from ref. 5, VH2–CNT,||(ρ), as
image file: c8cp04109a-t1.tif(1)
As discussed in ref. 5, the VH2–CNT,|| values are obtained by applying an additive pairwise potential model (APPM) for the dispersionless and dispersion contributions of the adsorbate–nanotube interaction. The parameters of the APPM were extracted by fitting to SAPT(DFT)-based values of dispersionless and dispersion terms. This APPM is specially designed to characterize the interaction of adsorbates with curved carbon-based nanostructures5,8,23 and metallic solid surfaces.48 It should be stressed that the influence of the nanotube curvature-induced dipole49 in the H2–CNT interaction is implicitly accounted for in the potential model through the dependence of dispersionless and dispersion components on the angle between the radial vector going from the nanotube center to one carbon atom and the vector pointing from the H2 center-of-mass to the same C atom.5 Applying the potential model, the interaction of the whole H2 molecule with nanotubes is thus described.
2.2.2 Kinetic energy operator and variable representation. We use cylindrical coordinates R, Z, Φ for the center-of-mass (COM) H2 position and spherical coordinates r, x = cos(θ), ϕ to account for the relative motion of the two hydrogen atoms. After the corresponding Jacobian transformation of the wave-functions, the kinetic energy operator can be written as:
image file: c8cp04109a-t2.tif
where M is the total mass of the hydrogen molecule, and μ is the reduced mass. Since the H2–nanotube interaction potential does not depend on Z (i.e., on the overall translational motion), the Z degree of freedom can be simply omitted. Next, the H2 internal vibrational motion is not considered here since the characteristic vibrational frequency (about 546 meV) is much higher than the energy range studied in the present work (below 65 meV). Therefore, the coordinate r is fixed to the vibrationally averaged rigid rotor value r0 and the corresponding term is omitted from the kinetic energy operator. In principle, more accurate approaches including non-adiabatic corrections could be applied but these corrections are small and neglected in this work. For further simplifications, let us write the expression of the variables on which the adsorbate–nanotube VH2–CNT interaction depends, i.e., ρ1 and ρ2. In what follows, the ρ1 and ρ2 degrees of freedom which we referred to as ρ+ and ρ:
image file: c8cp04109a-t3.tif
This expression suggests that a suitable transformation is,
ϕ[small phi, Greek, tilde] = ϕΦ
Omitting tildes in what follows, the kinetic energy operator becomes
image file: c8cp04109a-t4.tif
image file: c8cp04109a-t5.tif(2)
Hence, the coordinate Φ can be separated out along with the corresponding integer overall rotation quantum number Λ and the associated wave function
image file: c8cp04109a-t6.tif
After this separation, the final problem becomes three dimensional (3D) with a Λ-dependent kinetic energy operator:
image file: c8cp04109a-t7.tif(3)
It can be noticed that the total 3D Hamiltonian is invariant under the symmetry transformation x → −x, and for Λ = 0, ϕ → 2π − ϕ. Moreover, since the adsorbate–nanotube interaction potential is a sum of terms with ρ1 and ρ2 as the variables, the total Hamiltonian is also invariant under the transformation ϕϕ + π (that corresponds to ρ1ρ2).
2.2.3 Calculation of energy levels. The resulting 3D Hamiltonian problem is solved using the Discrete Variable Representation (DVR) approach50 and sinc-DVR functions5,8 for radial coordinate R and analytic basis of spherical harmonics for angular coordinates. To be able to perform analytical integration for angular coordinates, the model H–CNT interaction potential VH–CNT (eqn (1)) is first expanded in a polynomial series
image file: c8cp04109a-t8.tif(4)
Using the explicit expression for ρ±2, eqn (2), this can be rewritten in the form
image file: c8cp04109a-t9.tif(5)
with coefficients PJIK being a linear combination of VL. Explicit expressions of the coefficients as well as their derivations are provided in Section S1 of the ESI.

The integration over x and ϕ can then be performed analytically using spherical harmonics as a basis, Yjm(θ,ϕ), with x = cos(θ) and jm. Using this basis, the Hamiltonian matrix can be explicitly written as

image file: c8cp04109a-t10.tif(6)
where i stands for the i-th radial DVR function, with Ri as the corresponding DVR grid point. The details for the analytical calculations of the angular matrix elements 〈jm′|(1 − x2)I[thin space (1/6-em)]cos2K(ϕ)|jm〉 are given in Section S3 of the ESI. The symmetry noted above translates into the fact that odd and even j are not mixed, as well as odd and even m. Moreover, for Λ = 0, symmetric and anti-symmetric functions of ϕ, i.e. cos() and sin(), do not mix either. Since an important part of the wave function can be close to R = 0 for nanotubes of small diameter, special attention should be paid to the explicit evaluation of image file: c8cp04109a-t11.tif.

2.2.4 Computational details in adsorbate wave-function calculations. The interaction potential VH–CNT was fitted using the polynomial expansion of eqn (4) with n values ranging from 4 to 8 upon increasing the nanotube diameter (see Section S1 of the ESI for pictures showing the accuracy of the fittings and expansion coefficients). A basis of spherical harmonics with mmax = 10 and jmax = 30 was used to account for the angular degrees of freedom. For the radial R coordinate we used from 61 to 90 Sinc-DVR functions on the intervals going from 0 to 9 Bohr (from 17 to 25 Bohr) for H2 rotating inside (outside) carbon nanotubes.

3 Results and discussion

3.1 Influence of rotational H2 motion on the H2–nanotube interaction

For a short and narrow nanotube of helicity index (5,5) and a diameter of 6.74 Å, the left-hand panel of Fig. 1 shows the total H2/sCNT(5,5) interaction potential for three orientations of the H2 molecule as a function of the radial distance R between the molecule center-of-mass and the nanotube center. The electrostatic Eelec, exchange–repulsion Eexch–rep, induction Eind, and dispersion Edisp contributions are also shown. As illustrated in the picture of the right-hand panel, θ = 0° (full lines) corresponds to a parallel orientation of the H2 to the tube axis while θ = 90° (dashed lines) is associated with a perpendicular configuration instead.
image file: c8cp04109a-f1.tif
Fig. 1 Radial scan of the total interaction energies between a single H2 molecule and the short carbon nanotube (sCNT) of helicity index (5,5) represented in the right-hand panel. Gray spheres represent carbon atoms while blue spheres stand for one H2 molecule inside. The sCNT was saturated with hydrogen atoms (not shown here). The interaction energies are represented as a function of the distance from the H2 center-of-mass to the nanotube axis, R. Full, dotted and dashed lines correspond to θ values of 0°, 45°, and 90°, respectively, with θ as the angle between the H–H internuclear axis and the nanotube axis. The different SAPT-based energy contributions to the total interaction energy are also shown. The H–H bond length was fixed to the vibrationally averaged rigid rotor value, r0 = 0.7508 Å, from ref. 18.

As can be observed in Fig. 1 (left-hand panel) and also discussed in ref. 5, the H2/CNT attractive interaction is dispersion-dominated, with the net repulsive dispersionless contribution mainly characterized by the exchange–repulsion term. The potential minimum is located at the center of the narrow nanotube because this symmetric position allows the adsorbate to benefit from the dispersion interaction with carbon atoms at both sides of the carbon cage. The exchange–repulsion contribution grows exponentially as the distance between the H2 and the carbon cage decreases although such behavior is somewhat smoothed out by the attractive electrostatic contribution. The induction term contributes very little in the relevant range of radial distances in this work.

Focusing on the influence of the particular H2 orientation on the total energy, it can be observed that the parallel orientation (θ = 0°, full lines) is clearly energetically preferred over the perpendicular configuration (θ = 90°, dashed lines). It can also be seen that this preference arises from a smaller exchange–repulsion contribution when the molecule is located along the tube axis. Contrarily, the dispersion contribution becomes more pronounced as the θ value increases. Therefore, the preferential orientation arises from the counterbalance between exchange–repulsion and dispersion terms.

The subtle balance between dispersion and dispersionless contributions is very much influenced by the diameter of the nanotube. To illustrate this, we apply the potential model to longer nanotubes with helicity indexes (5,5), (10,5) and (10,10), obtaining the θ-dependence of total, dispersionless and dispersion energies presented in Table 1 (see also Fig. 2). First, it can be observed in Table 1 that the potential model provides the total energies at the potential minimum (values in parentheses) agreeing with the SAPT(DFT) values to within 6%. Second, it can also be seen that the weight of the exchange–repulsion at the potential minimum decreases as the diameter of the tube increases and, then, the dispersion becomes more relevant in determining the θ-dependence of the total energies. This is the reason why the energy difference between parallel and perpendicular orientations becomes smaller when the tube diameter increases, as can be clearly seen in Fig. 2. From Table 1, we can also notice the very weak ϕ-dependence of the total energies, reflecting the small corrugation of the carbon material.

Table 1 Total energies (Etotmin, in meV) and positions (R0, in Å) of the potential minima for the interaction between a single H2 molecule rotating in carbon nanotubes (CNTs) with helicity indexes (n,m) and a diameter (dCNT, in Å) of increasing value. The term “s(5,5)” stands for the short CNT used in the SAPT(DFT) calculations (see the right-hand panel of Fig. 1). Nanotubes with 10 cell units (a length equal or larger than 23.36 Å) were chosen when applying the potential model. Values in parentheses have been obtained using this model
CNT s(5,5)/SAPT(DFT) (5,5)/model (10,5)/model (10,10)/model
d CNT 6.7 6.7 10.4 13.6
R 0 0.0 0.0 2.0 3.6

(θ,ϕ) E disp-lessmin E dispmin E totmin E disp-lessmin E dispmin E totmin E disp-lessmin E dispmin E totmin E disp-lessmin E dispmin E totmin
(0°,0°) 90 −236 −146 (−148) 96 −279 −183 37 −127 −90 31 −103 −72
(45°,0°) 108 −247 −139 (−137) 118 −289 −171 51 −133 −82 44 −109 −65
(90°,0°) 123 −254 −131 (−124) 141 −300 −159 66 −139 −73 57 −114 −57
(45°,45°) 108 −247 −139 (−137) 118 −289 −171 46 −132 −86 32 −100 −68

image file: c8cp04109a-f2.tif
Fig. 2 Scan of relative energies EEmin as a function of θ for one H2 molecule inside carbon nanotubes (CNTs) of helicity indexes (5,5), (10,5), and (10,10).

3.2 Molecular energy levels of H2 in carbon nanotubes

As illustrative examples, we have calculated the molecular energy levels of H2 inside nanotubes with helicity indexes (5,5), (10,5), and (10,10), having diameters of 6.7, 10.4, and 13.6 Å, respectively, as well as the levels outside the (10,10) nanotube. These energy levels are assigned using three quantum numbers, n for the radial R degree of freedom, j for the molecular rotation associated with x = cos(θ), and m for the projection of the angular momentum quantum number corresponding to ϕ. For all nanotubes, the motions in R and x, ϕ are quite well separated so that the vibrational quantum number n and (jm) can be uniquely attributed to each adsorbate bound state. These states have been calculated for several values of the total angular momentum projection Λ along the tube axis. It should be noticed that, for Λ = 0 and m > 0, we can distinguish between symmetric (referred to as “ms”) and anti-symmetric (“ma”) energy levels. For larger nanotubes and Λ > 0 this symmetry is still approximately preserved, in contrast to the (5,5) nanotube where the image file: c8cp04109a-t12.tif term plays a major role so that m itself becomes a better quantum number. The energies of the most representative bound states are presented in Table 2 while the complete list of the calculated molecular energy levels is provided in Section S2 of the ESI.
Table 2 Energy levels (in cm−1) of a rotating H2 molecule inside nanotubes of helicity indexes (5,5), (10,5) and (10,10) and a diameter (dCNT, in Å) of increasing value. For the (10,10) nanotube, results for the H2 molecule placed outside the carbon cage are also given (denoted as CNT(10,10)-ext)
CNT(5,5) CNT(10,5) CNT(10,10) CNT(10,10)-ext
The energies at the potential minima (Emin, in cm−1) and the zero-point energies (Ezp, in cm−1) are also collected. The energy values are relative to that of the ground state (0 0 0). The energy of the ground state (0 0 0) is given in parentheses. For Λ > 0, ΔEΛ = (EEΛ=0)/Λ2. The energies of the lowest bound states without including the rotation are presented in brackets. Additional energy levels are provided in Section S2 of the ESI.
d CNT, Å 6.7 10.5 13.6 13.6
E min, cm−1 −1474 −725 −575 −296
E zp, cm−1 474 99 99 81

n j m E Λ=0 ΔEΛ=1 E Λ=0 ΔEΛ=1 E Λ=0 ΔEΛ=1 E Λ=0 ΔEΛ=1
0 0 0 0.0 154.088 0.0 2.905 0.0 0.727 0.0 0.081
(−997.908) (−625.627) (−475.369) (−215.515)
0 [−1147.07] [−650.433] [−502.037] [−233.840]
0 1 0 61.429 150.275 110.148 2.778 109.371 0.714 112.222 0.082
0 1 1s 340.182 231.611 140.229 4.183 138.983 0.828 130.974 0.082
0 1 1a 309.649 −144.839 112.058 1.646 109.775 0.645 112.313 0.081
0 2 0 338.445 154.303 345.202 2.072 342.789 0.680 346.558 0.081
0 2 1s 510.758 229.605 367.277 4.263 365.256 0.829 361.193 0.083
0 2 1a 488.408 −143.881 345.509 1.360 342.748 0.619 346.560 0.080
0 2 2s 825.042 284.147 382.259 11.179 377.855 −0.254 367.951 0.092
0 2 2a 835.032 −251.237 375.399 −4.388 367.170 0.129 361.452 0.071
1 0 0 522.302 169.873 117.289 5.473 123.191 0.895 93.195 0.076
[486.978] [126.178] [130.899] [100.079]
1 1 0 571.230 165.199 230.937 5.068 235.674 0.876 208.018 0.076
1 1 1s 912.150 233.919 256.851 8.366 258.877 1.070 220.864 0.076
1 1 1a 857.645 −142.973 234.331 0.745 236.068 0.748 208.069 0.075
1 2 0 877.251 151.240 467.850 3.594 470.371 0.765 443.506 0.074
1 2 1s 1058.936 239.767 486.606 8.764 487.858 1.082 453.361 0.077
1 2 1a 1019.278 −148.993 469.366 0.100 470.381 0.698 443.503 0.074
1 2 2s 1390.569 320.492 505.672 20.008 499.811 0.215 458.592 0.088
1 2 2a 1414.431 −267.805 500.797 −10.555 490.208 −0.112 453.576 0.064
2 0 0 1087.333 183.886 205.702 14.606 218.204 1.247 155.322 0.069
[1022.891] [220.287] [231.836] [166.508]
2 1 0 1123.970 177.651 321.341 13.555 332.906 1.211 271.706 0.069
2 1 1s 1521.755 221.946 353.967 18.234 350.725 1.632 279.919 0.069
2 1 1a 1435.031 −137.1229 332.991 −6.195 333.508 0.886 271.743 0.067
2 2 0 1465.285 141.410 560.816 9.167 569.158 0.498 507.901 0.067
2 2 1s 1646.533 247.978 584.240 19.849 581.976 1.679 514.236 0.071
2 2 1a 1592.477 −152.7319 568.653 −7.433 568.833 0.771 507.907 0.067
2 2 2s 1992.091 351.905 614.377 33.948 592.094 4.737 518.548 0.078
2 2 2a 2029.057 −279.984 612.508 −22.713 585.367 −1.427 514.422 0.063

From the values of zero-point energies Ezp in Table 2, it can be observed that the strong confinement of H2 upon adsorption inside the (5,5) tube leads to a considerable zero-point energy due to the steepness of the interaction potential along the radial R distance (see Fig. 1). Similar to atomic helium as the adsorbate,7 it arises from a very pronounced exchange–repulsion contribution even at the center of the nanotube. This finding is consistent with the experimental observations of Ohba,6 showing the quenched adsorption of species as light as helium atoms in nanopores with diameters below 7 Å at low (2–5 K) temperatures.

It is also interesting to compare the lowest-energy states in Table 2 with those calculated by considering only the H2 center-of-mass motion R in ref. 5, with the H2 molecule lying parallel to the tube axis. The energies of these states are shown in brackets and they would be associated with the set (j m) = (0 0). Notice that the wave function for states with j = m = 0 is quite broad in θ, which results in penalty from the less attractive interaction potential for a perpendicular configuration and a higher energy when including the H2 rotation into our treatment. The energy differences with and without including the H2 rotation are more pronounced for H2 inside the narrow nanotube where, for (j m) = (0 0), the number of vibrational bound states is reduced from three to just two. Contrarily, the less steep H2–CNT(10,10) interaction potential supports more than six vibrational states with (j m) = (0 0).

To better understand the dependence of the energy levels on the different quantum numbers, the total energies have been decomposed in the following contributions: (1) the average potential energy term 〈VH2–CNT〉; (2) the kinetic R-dependent term, 〈K(R)〉, with image file: c8cp04109a-t13.tif (i.e., the first term from the kinetic energy operator in eqn (3)); (3) the “internal” kinetic term 〈Kint〉 (i.e., the terms multiplied by the H2 rotational constant image file: c8cp04109a-t14.tif in eqn (3)); and (4) the “external” kinetic term 〈Kext〉 (i.e., the terms multiplied by the “external” rotational constant image file: c8cp04109a-t15.tif in eqn (3)). This decomposition is presented in Table 3 for the particular case of the (5,5) nanotube.

Table 3 Decomposition of total energies (Etot, in cm−1) of molecular states of H2 inside CNT(5,5) into the following contributions: average potential term 〈VH2–CNT〉, kinetic R-dependent term, 〈K(R)〉, kinetic (internal) term, 〈Kint〉, and kinetic (external) term 〈Kext〉. Energy differences ΔEtot with respect to the ground-state energies are also provided
n j m Λ E tot ΔEtot VH2–CNT K(R)〉 Kint Kext
0 0 0 0 −998 0.0 −1206 188 19 0.1
0 1 0 0 −936 61 −1246 184 125 0.0
0 1 1s 0 −658 340 −1073 168 128 120
0 1 1a 0 −688 310 −1087 160 126 113
0 2 0 0 −659 338 −1191 191 340 0.2
1 0 0 0 −476 522 −986 475 34 0.9
2 0 0 0 89 1087 −760 796 50 3.3
0 0 0 1 −844 154 −1138 159 23 112
0 1 0 1 −786 212 −1178 156 126 109
0 1 1s 1 −426 572 −971 163 129 253
0 1 1a 1 −833 165 −1153 195 125 0.3

Focusing first on the energy levels with (m Λ) = (0 0), we can notice from Tables 2 and 3 that their j-dependence deviates by just 10% from that of gas-phase H2 molecules and expressed as BH2j(j + 1). For Λ > 0, it can be observed that the Λ-dependence of the energy levels is very weak for the wide CNT(10,10) tube, being simply accounted for with the energy contribution BextΛ2. This can be explained by considering the large value of the equilibrium distance R0, of 3.6 Å (see Table 1), leading to a very small value of the “external” rotational constant image file: c8cp04109a-t16.tif (about 1 cm−1).

As opposed to the case of the wide CNT(10,10) tube, the potential minimum for the H2–CNT interaction is located at the center of narrow nanotubes (i.e., R0 = 0) so that the value of the rotational constant Bext is very large, differing by less than 13% from the value of the H2 rotational constant BH2. The quantum numbers Λ and m are thus strongly correlated and the Λ-dependence of the energies can be approximated as Bext(Λm)2. Therefore, as can be observed Table 3, the “external” kinetic term 〈Kext〉 is almost zero for levels with Λ = m.

Focusing now on the energy levels sharing the same (n j) quantum numbers and different m values, we can notice that these m-splittings decrease as the diameter of the nanotube increases. This feature is directly related to the anisotropy (i.e., the θ-dependence) of the H2–CNT interaction. For instance, the θ-dependence of the wave-functions for (0 1 0) and (0 1 1s) levels is mainly described with the functions cos(θ) and sin(θ), favoring parallel and perpendicular configurations, respectively. As can be seen in Fig. 2, the parallel configuration is energetically favored for H2 molecules inside carbon nanotubes. Consequently, the (0 1 0) level is below the (0 1 1s) levels for any nanotube size.

It is important to note that the levels (n 1 1a) and (n 1 0) are nearly degenerate with the exception of the narrow (5,5) tube. This outcome can be explained as follows: the most important term in the potential expansion splitting m = 0 and m = 1 states corresponds to that with indexes I = K = 1, see eqn (5) above (the terms with K = 0 give no splitting), which is directly related to the V2 coefficient. However, it can be demonstrated that the I = K = 1 term contributes equally to (n 1 0) and (n 1 1a) Hamiltonian matrix elements, while it is three times larger for the (n 1 1s) counterparts. A detailed analysis of the degeneracy of (n 1 1a) and (n 1 0) levels is provided in Section S4 of the ESI. For the narrow (5,5) tube, the anisotropy of the H2–CNT interaction is much more marked. In this case, the Bext(Λm)2 rotational term contributes to an increase in the energy of the (0 1 1) level by ca. 250 cm−1. The energy pattern is then determined by the whole combination of different contributions. This way, the energy differences between (0 1 0) and (0 1 1) levels are much larger such as the levels (n 1 1a) and (n 1 0) are no more degenerate. In stark contrast, the m-splitting of (n j) = (0 1) levels inside and outside of the wider (10,10) tube is below 20 cm−1 (see Table 1).

Recent experimental measurements using high-resolution neutron spectroscopy have provided precise values of m-splittings for H2 adsorbed in the outside of carbon nanohorns14 and aligned carbon nanotubes.15 Specifically, values of ca. 4 and 8 cm−1 have been reported for the energy difference between (0 1 0) and (0 1 1) levels in carbon nanohorns14 and carbon nanotubes, respectively.15 By extending our potential model to the outside of the CNT(10,10) nanotube (see Table 1), the estimated value of the m-splitting is 18 cm−1. This value is probably overestimated because our model potential is based on ab initio calculations with H2 placed inside the carbon cage. As shown in ref. 8, the repulsive core of adsorbate–CNT interaction is less marked for the adsorbate localized in the outside of the carbon cage and, then, the anisotropy of the interaction. Moreover, as discussed in ref. 49, it is well known that there is an electronic charge displacement directed toward the outside of the carbon cage in nanotubes. It is clear that this redistribution of the electronic density must affect the H2–nanotube interaction, in particular the exchange–repulsion contribution. Hence, more precise estimations would require the refitting of dispersionless contributions using ab initio values for H2 placed outside the nanotube. Using preliminary calculations of SAPT(DFT) interaction energies with the H2 molecule located outside of the carbon cage, a refitting of the pairwise potential model indicates that the m-splitting is reduced from 18 to 11 cm−1, with the global structuring of molecular energy levels being almost unmodified. Finally, a suitable modeling of assemblies of CNTs used in experimental measurements15 is probably needed.

4 Conclusions

We presented a new potential model and adsorbate wave-function method based on accurate vdW ab initio calculations and an explicit treatment of the relevant nuclear degrees of freedom. The method has been applied to provide accurate energy levels of molecular hydrogen rotating inside and outside carbon nanotubes of various sizes, starting from the subnanometer regime. An analysis based on SAPT-DFT theory has been performed to estimate the influence of the different physical contributions to the H2–nanotube interaction energy on the rotating molecule. This analysis highlights the dramatic modification of the interaction energy landscapes and molecular energy patterns in going from carbon tubes of sub-nanometer diameter to nanometer sizes. In particular, strong anisotropy effects constrain the rotation of molecular hydrogen in narrow nanotubes. The new computational tool presented in this work will allow interpretation of future experimental measurements using high-resolution neutron spectroscopy of H2 molecules inside carbon nanotubes. Work is in progress to refine our potential model and to provide more accurate estimations of molecular energy levels for H2 molecules adsorbed outside carbon nanotubes and nanotube assemblies. The application of our method to other diatomic molecules in carbon nanotubes is also envisaged.

Conflicts of interest

There are no conflicts to declare.


This work has been partly supported by the Spanish Agencia Estatal de Investigación (AEI) and the Fondo Europeo de Desarrollo Regional (FEDER, UE) under Grant No. MAT2016-75354-P and COST Action CM1405 “Molecules in Motion” (MOLIM). MPdLC is greatly thankful to the CTI (CSIC) and CESGA super-computer facilities (Spain) for the provided computational resources.


  1. P. Kowalczyk, A. P. Terzyk, P. A. Gauden, S. Furmaniak, K. Kaneko and T. F. Miller III, J. Phys. Chem. Lett., 2015, 6, 3367–3372 CrossRef CAS PubMed.
  2. T. T. Nguyen and P. Serp, ChemCatChem, 2013, 5, 3595–3603 CrossRef.
  3. V. P. Ting, A. J. Ramirez-Cuesta, N. Bimbo, J. E. Sharpe, A. Noguera-Diaz, V. Presser, S. Rudic and T. J. Mays, ACS Nano, 2015, 9, 8249–8254 CrossRef CAS.
  4. R. J. Olsen, A. K. Gillespie, C. I. Contescu, J. W. Taylor, P. Pfeifer and J. R. Morris, ACS Nano, 2017, 11, 11617–11631 CrossRef CAS.
  5. M. P. de Lara-Castells, A. W. Hauser, A. O. Mitrushchenkov and R. Fernández-Perea, Phys. Chem. Chem. Phys., 2017, 19, 28621–28629 RSC.
  6. T. Ohba, Sci. Rep., 2016, 6, 28992 CrossRef CAS.
  7. A. W. Hauser and M. P. de Lara-Castells, J. Phys. Chem. Lett., 2016, 7, 4929–4935 CrossRef CAS.
  8. A. W. Hauser, A. O. Mitrushchenkov and M. P. de Lara-Castells, J. Phys. Chem. C, 2017, 121, 3807–3821 CrossRef CAS.
  9. T. X. Nguyen, H. Jobic and S. K. Bhatia, Phys. Rev. Lett., 2010, 105, 085901 CrossRef CAS.
  10. M. Mondelo-Martell and F. Huarte-Larrañaga, J. Phys. Chem. A, 2016, 120, 6501–6512 CrossRef CAS.
  11. M. Rossi and F. Ancilotto, Phys. Rev. B, 2016, 94, 100502 CrossRef.
  12. G. Ferré, M. C. Gordillo and J. Boronat, Phys. Rev. B, 2017, 95, 064502 CrossRef.
  13. A. Del Maestro and M. Boninsegni, Phys. Rev. B, 2017, 95, 054517 CrossRef.
  14. F. Fernández-Alonso, F. J. Bermejo, C. Cabrillo, R. O. Loutfy, V. Leon and M. L. Saboungi, Phys. Rev. Lett., 2007, 98, 215503 CrossRef.
  15. H. G. Schimmel, G. J. Kearley and F. M. Mulder, ChemPhysChem, 2004, 5, 1053–1055 CrossRef CAS.
  16. D. G. A. Smith and K. Patkowski, J. Phys. Chem. C, 2014, 118, 544–550 CrossRef CAS.
  17. D. G. A. Smith and K. Patkowski, J. Phys. Chem. C, 2015, 119, 4934–4948 CrossRef CAS.
  18. M. P. de Lara-Castells and A. O. Mitrushchenkov, J. Phys. Chem. A, 2015, 119, 11022–11032 CrossRef CAS.
  19. M. d. R. Cuevas-Flores, M. A. García Revilla and M. Bartolomei, J. Comput. Chem., 2018, 39, 71–80 CrossRef CAS.
  20. M. Bartolomei, R. Perez de Tudela, K. Arteaga, T. González Lezana, M. I. Hernandez, J. Campos-Martínez, P. Villarreal, J. Hernández Rojas, J. Breton and F. Pirani, Phys. Chem. Chem. Phys., 2017, 19, 26358–26368 RSC.
  21. M. Bartolomei, E. Carmona-Novillo, M. I. Hernández, J. Campos-Martínez and F. Pirani, J. Phys. Chem. C, 2013, 117, 10512–10522 CrossRef CAS.
  22. A. Hesselmann and T. Korona, Phys. Chem. Chem. Phys., 2011, 13, 732–743 RSC.
  23. A. W. Hauser and M. P. de Lara-Castells, Phys. Chem. Chem. Phys., 2017, 19, 1342–1351 RSC.
  24. H. Stoll, J. Chem. Phys., 1992, 97, 8449–8454 CrossRef CAS.
  25. E. Voloshina, D. Usvyat, M. Schütz, Y. Dedkov and B. Paulus, Phys. Chem. Chem. Phys., 2011, 13, 12041–12047 RSC.
  26. M. P. de Lara-Castells, H. Stoll, B. Civalleri, M. Causà, E. Voloshina, A. O. Mitrushchenkov and M. Pi, J. Chem. Phys., 2014, 141, 151102 CrossRef.
  27. M. P. de Lara-Castells, A. O. Mitrushchenkov and H. Stoll, J. Chem. Phys., 2015, 143, 102804 CrossRef.
  28. M. P. de Lara-Castells, N. F. Aguirre, H. Stoll, A. O. Mitrushchenkov, D. Mateo and M. Pi, J. Chem. Phys., 2015, 142, 131101 CrossRef PubMed.
  29. M. P. de Lara-Castells, M. Bartolomei, A. O. Mitrushchenkov and H. Stoll, J. Chem. Phys., 2015, 143, 194701 CrossRef PubMed.
  30. K. Pernal, R. Podeszwa, K. Patkowski and K. Szalewicz, Phys. Rev. Lett., 2009, 103, 263201 CrossRef PubMed.
  31. R. Podeszwa and K. Szalewicz, J. Chem. Phys., 2012, 136, 161102 CrossRef.
  32. R. Podeszwa, K. Pernal, K. Patkowski and K. Szalewicz, J. Phys. Chem. Lett., 2010, 1, 550–555 CrossRef CAS.
  33. M. P. de Lara-Castells, G. Delgado-Barrio, P. Villarreal and A. O. Mitrushchenkov, J. Chem. Phys., 2009, 131, 194101 CrossRef CAS.
  34. M. P. de Lara-Castells and A. O. Mitrushchenkov, J. Phys. Chem. Lett., 2011, 2, 2145–2151 CrossRef CAS.
  35. N. F. Aguirre, P. Villarreal, G. Delgado-Barrio, A. O. Mitrushchenkov and M. P. de Lara-Castells, Phys. Chem. Chem. Phys., 2013, 15, 10126–10140 RSC.
  36. A. J. Misquitta, B. Jeziorski and K. Szalewicz, Phys. Rev. Lett., 2003, 91, 033201 CrossRef.
  37. A. Heßelmann and G. Jansen, Chem. Phys. Lett., 2003, 367, 778–784 CrossRef.
  38. Z. Qiuchen and Z. Jin, Small, 2014, 10, 4586–4605 CrossRef.
  39. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, et al., MOLPRO, version2012.1, a package of ab initio programs, see http://www.molpro.net Search PubMed.
  40. A. Heßelmann, G. Jansen and M. Schütz, J. Chem. Phys., 2005, 122, 014103 CrossRef.
  41. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868 CrossRef CAS.
  42. D. E. Woon and T. H. Dunning, Jr., J. Chem. Phys., 1994, 100, 2975–2988 CrossRef CAS.
  43. F. Weigend, Phys. Chem. Chem. Phys., 2002, 4, 4285–4291 RSC.
  44. F. Weigend, A. Köhn and C. Hättig, J. Chem. Phys., 2002, 116, 3175–3183 CrossRef CAS.
  45. M. Grüning, O. V. Gritsenko, S. V. A. van Gisbergen and E. J. Baerends, J. Chem. Phys., 2001, 114, 652–660 CrossRef.
  46. S. Lias, Ionization Energy Evaluation, in NIST Chemistry WebBook, NIST Standard Reference Database Number 69, ed. P. Linstrom and W. Mallard, National Institute of Standards and Technology, Gaithersburg, 2005 Search PubMed.
  47. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170 CrossRef CAS.
  48. M. P. de Lara-Castells, R. Fernández-Perea, F. Madzharova and E. Voloshina, J. Chem. Phys., 2016, 144, 244707 CrossRef.
  49. T. Dumitrică, C. M. Landis and B. I. Yakobson, Chem. Phys. Lett., 2002, 360, 182–188 CrossRef.
  50. Z. Bačic and J. C. Light, Annu. Rev. Phys. Chem., 1989, 40, 469–498 CrossRef.


Electronic supplementary information (ESI) available: Details of the fittings to potential models, four tables with complete lists of calculated molecular energy levels, complementary expressions for the calculation of Hamiltonian matrix elements, and analysis of degenerate energy levels appearing in wide nanotubes. See DOI: 10.1039/c8cp04109a

This journal is © the Owner Societies 2019