Quantum confinement of molecular deuterium clusters in carbon nanotubes: ab initio evidence for hexagonal close packing

An ab initio study of quantum confinement of deuterium clusters in carbon nanotubes is presented. First, density functional theory (DFT)-based symmetry-adapted perturbation theory is used to derive parameters for a pairwise potential model describing the adsorbate–nanotube interaction. Next, we analyze the quantum nuclear motion of N D2 molecules (N o 4) confined in carbon nanotubes using a highly accurate adsorbate-wave-function-based approach, and compare it with the motion of molecular hydrogen. We further apply an embedding approach and study zero-point energy effects on larger hexagonal and heptagonal structures of 7–8 D2 molecules. Our results show a preference for crystalline hexagonal close packing hcp of D2 molecules inside carbon nanotubes even at the cost of a reduced volumetric density within the cylindrical confinement.


Introduction
Clean sources of energy and efficient energy storage are at the core of the current research in nanoscience and nanotechnology.Hydrogen is being actively pursued as a substitute for fossil fuels since its combustion produces only heat and water, and it can be efficiently combined with oxygen in a fuel cell to produce electricity.However, realizing the promise of hydrogen as a widely used fuel still requires a substantial development of efficient storage materials, 1,2 including metal organic frameworks 3 and carbon-based nanoporous materials. 4,5he general goal in hydrogen storage methods is to pack hydrogen molecules as close as possible.There exists direct experimental evidence for solid-like packing of hydrogen at supercritical temperatures. 6Moreover, recent low-temperature (11-22 K) neutron-diffraction-based experimental measurements by Cabrillo et al. have provided evidence for the formation of a crystalline hcp structure of molecular deuterium clusters inside carbon nanotubes (CNTs). 7,8This work has raised the fundamental question of the role of nuclear quantum effects for hydrogen adsorption and confinement at low temperatures.Provided that the inter-adsorbate interactions are strong, the crystal formation of the adsorbate inside CNTs has been observed with a symmetry reflecting that of the adsorbate bulk crystal since a long time ago. 9In stark contrast, when the interadsorbate interactions are weak, both experimental observations and theoretical calculations have revealed that the structures associated to the confined (solid) adsorbate depart largely from those corresponding to the bulk crystal (for instance, considering ice 10,11 ).In this sense, the experimental observation of hcp structures for adsorbates characterized by very weak interactions such as D 2 -D 2 is relevant and needs a fundamental understanding.
Very recent studies have addressed the possible existence of either a superfluid 12 or a crystalline phase 13,14 for parahydrogen molecules inside carbon nanotubes at zero temperature 12,14 or temperatures below 4 K, 13 revealing the impact of quantum nuclear effects (see also ref. 15 for a recent review).Motivated by active experimental research, the last few years have seen marked progress in the development of methods to describe the quantum nuclear motion of H 2 and D 2 isotopes encapsulated by carbon nanotubes, including diffusion dynamics (see, e.g., ref. 16 and 17 and references therein).These theoretical studies have provided explanations of relevant confinement effects such as the quantum-induced reversed trend in H 2 and D 2 diffusion rates upon lowering the temperature, 17 confirming the experimental observations, 18 and have revealed a strong influence of the adsorbate-nanotube interaction potential on the rate values.
Modern intermolecular interaction theory allows an ab initio characterization of van der Waals-dominated interactions of molecular or atomic adsorbates with carbon surfaces such as the D 2 -CNT interaction.4][25][26][27][28][29][30][31][32][33][34][35][36][37] For the specific case of molecule-CNT interaction, an accurate treatment by Patkowski and collaborators includes the benchmarking of dispersion-including DFT functionals with composite second-order Mo ¨ller-Plesset perturbation theory and the coupled-cluster approach applied on models of carbon nanotubes. 41,42his work presents an ab initio study of the quantum nuclear motion of deuterium clusters in carbon nanotubes at zero temperature.For this purpose, we propose an additive pairwise potential model for the D 2 -nanotube interaction which relies on high-level ab initio calculations using the SAPT(DFT) approach. 19,20The latter technique assumes that the zero-order wavefunction of the total system is a product of the wavefunction of each isolated subsystem, with antisymmetry being taken into account at each order of the perturbation expansion.It is worth noting that the SAPT(DFT) approach has been previously applied to describe molecular hydrogen adsorption on planar carbon surfaces, 23,25 achieving a remarkable agreement with experimental measurements of bound state energies.The additive pairwise potential model uses the partition of van der Waals interaction into dispersion and dispersionless contributions.Specifically, the parameters of analytical dispersionless and dispersion terms are tuned up for a small representative CNT using SAPT(DFT)-based dispersionless and dispersion energies.As discussed in ref.31, the potential model relies on two basic conclusions from accurate studies of van der Waals-dominated adsorbate/surface systems (see, e.g., ref. 37 and 43) (1) the shortrange nature of the dispersionless contribution allows an appropriate tuning of model parameters using ab initio energies on small representative cluster models; (2) the dispersion is longrange, but dispersion parameters show excellent transferability properties upon increasing the size of the surface cluster models.
To describe the molecular motion of D 2 clusters, we apply a highly accurate adsorbate wave-function-based approach, an embedding scheme as well as a vibrational model to present evidence of the hcp structure observed in the experiment. 7,8he adsorbate wave-function-based approach consists in explicitly solving the nuclear Schro ¨dinger equation in the real space, using the the Discrete Variable Representation (DVR) method 44 for the discretization, and treating the molecular adsorbate as a pseudoatom.Recent studies 23,25 have shown that the behaviour of the H 2 /coronene system is not affected by the relative orientation of the molecule, suggesting that it can be treated as a pseudo-atom.Within the selected model (adiabatic formulation, reduced dimension space for D 2 , and pairwise interaction approximation), the method provides an exact solution if convergence with the DVR basis is achieved.Appealing features of the method are the absence of phenomenological terms, the inclusion of bosonic symmetry, and the possibility of calculating excited states.The main disadvantage is in the limited number of adsorbates (up to three) that can be included in well-converged solutions.It has been very recently applied to describe the quantum nuclear motion of helium and molecular nitrogen clusters in carbon nanotubes, 31 revealing very different quantum confinement effects depending on the CNT diameter, and agreeing very well with experimental measurements. 45In the next section, we briefly present the methods used in this study.Our results are reported in Section 3. Finally, the main conclusions are provided in Section 4.

Adsorbate-nanotube interaction
The adsorbate/CNT interaction V AÀCNT is introduced by applying an additive pairwise potential model. 31As mentioned above, our model uses different functions for the dispersionless and dispersion energy contributions.It accounts for the typical exponential growth of the dominant dispersionless contribution, the exchange-repulsion but also including a Gaussian-type 'cushion' to describe weakly attractive tails stemming from other dispersionless terms where R c is a cut-off distance, R A-C stands for the distance between the adsorbate center-of-mass and one carbon atom of the nanotube, and y C is the angle between the radial vector going from the nanotube center to one carbon atom and the vector R A-C pointing from the adsorbate center-of-mass to the same C atom.The dimensionless factor g A in the first term accounts for the anisotropy of the C-C bonds.The sum in the second term runs over all carbon atoms of the nanotube.For the dispersion part, we apply the typical C 6 /C 8 expansion with the damping functions of Tang and Toennies f n (n = 6, 8) 46 where g R is also a dimensionless anisotropy parameter.The inclusion of g A and g R anisotropy terms has been found to be important when modelling potential corrugation effects on both curved carbon and metal surfaces. 30,31,47odel parameters are derived by fitting the dispersionless and dispersion SAPT(DFT) terms, as calculated for the interaction between a single H 2 and a hydrogen-saturated nanotube CNT (5,5) tube made of 62 atoms, considering two transverse sections of the tube.Using the SAPT(DFT) MOLPRO 21,48 This journal is © the Owner Societies 2017 Phys.Chem.Chem.Phys., 2017, 19, 28621--28629 | 28623 implementation, we followed the same computational setup as reported in ref.31, using the Perdew-Burke-Ernzerhof (PBE) density functional 49 and the augmented polarized correlationconsistent triple-zeta basis 50 for all atoms but the hydrogen atoms saturating the dangling bonds, for which the correlationconsistent double-zeta basis was used (see ref. 31 and ESI † for further details).

Adsorbate wave-function treatment
In order to calculate the adsorbate wave-functions, we have applied the same method proposed in our previous work. 31ery briefly, adsorbate molecules A are described as structureless bosonic particles.However, note that the potential felt by them is derived from the interactions of H 2 molecules in arrangements where their intramolecular axes are parallel to tube (i.e., the z) axis.The parallel orientation to the tube long axis has been chosen because it is mostly favorable energetically.Indeed any averaging over rotational motion of H 2 would lead to a less deep well due to the repulsive wall while practically leaving the potential close to the tube center unmodified.Thus, as can be seen in Fig. S2 of the ESI, † the adsorbate/nanotube potential for the parallel orientation differs only minimally from the one obtained by averaging over the J = 0 rotational state of molecular hydrogen.
Using cylindric coordinates r i (r i , f i , z i ) for the adsorbate center-of-mass position (see Fig. 1), the total Hamiltonian is written as a sum of one-particle kinetic energy terms and effective pairwise additive adsorbate/CNT (V AÀCNT ) and interadsorbate (V A-A ) interactions where R ij stands for the adsorbate-adsorbate distance, and K(r i ) is the one-particle kinetic energy term with M denoting the mass of the adsorbate.The very small corrugation along f is ignored and, considering very long nanotubes, it can be assumed that the interaction potential depends only on r (one-dimensional 1D model), V AÀCNT (r).Using this special form of the interaction potential, also the z component of the total angular momentum, L, and the total momentum is conserved and either the global rotation about z or the collective motion along z can be separated from the total adsorbate wave-function.As mentioned above, the resulting Hamiltonian problem can be discretized using the DVR approach 44 and the corresponding Schro ¨dinger equation is solved in the real space.The basis set is obtained as a direct product of functions for the different coordinates.For the adsorbate-nanotube r i degrees of freedom, we use the potentialoptimized DVR functions 51 which are built from 2D-harmonic oscillator functions, and for the relative inter-adsorbate coordinates, Sinc-DVR functions are employed instead (see ref. 31).The total Hamiltonian depends on four and seven degrees of freedom for N = 2 and 3, respectively.The dimension of the Hamiltonian for N = 3 with a proper grid basis is about half a billion so that it has been unfeasible to consider larger clusters.
To model the adsorbate-adsorbate interaction, we have used the ab initio potential energy surface reported by Hinde, 52 as restricted to our one-dimensional model.The values of the equilibrium distance (3.46 Å) and well-depth (24.72 cm À1 ) are very close to those of the well-known Silvera-Goldman potential 53 for the H 2 molecule in the J = 0 rotational state (3.41 Å and 23.84 cm À1 , respectively).

Embedding approach
To confirm that the hcp structure is preferred over a facecentered cubic fcc structure, we have applied an embedding approach as follows: first, we choose a nanotube with a length of 23.36 Å (10 cell units).Next, we define an embedding region comprising 14 D 2 molecules in two C 6 hexagonal structures with circumradius of value a as well as two D 2 molecules located at their centers (see TOC graphics).The two C 6 hexagonal structures are located in a plane perpendicular to the tube axis, separated by a distance of c/2 from z = 0.Then, the embedding potential acting on a D 2 molecule is simply with r and z denoting the cylindric coordinates of a reference molecule (i.e., center-of-mass position) and V D 2 -D 2 as the interaction potential between two isolated D 2 molecules.Note that the cylindrical symmetry of the problem is preserved.Finally, we obtain the adsorbate wave-functions of the embedded molecule as described in the previous section, but after adding the embedding potential to the adsorbate-nanotube interaction potential.In this case, however, the z-dependence of the adsorbate-nanotube potential is also included.Explicitly, the Schro ¨dinger equation of a single embedded D 2 molecule, This journal is © the Owner Societies 2017 needs to be solved for the Hamiltonian where m stands for z component of the total angular momentum along the tube long axis fixed to zero in this work.

Vibrational model
To further test whether the hcp structure is preferred for larger clusters, we have applied a vibrational model to single layers of D 2 molecules with either C 6 or C 7 symmetries inside the CNT(10,10) nanotube.We reduce our analysis to radial motions of the D 2 molecules in a plane perpendicular to the tube axis.
Using our additive pairwise potential model, an analytical expression of the Hessian matrix is obtained by calculating the second derivatives.This matrix is diagonalized to obtain 2N eigenvalues.Three of them are zero, corresponding to the two planar translations and one in-plane rotation.The remaining eigenvalues are converted into vibrational frequencies o i and the zero-point energy (ZPE) is determined via a summation, Note that internal vibration of the molecules and their rotational degrees of freedom are not taken into consideration.However, since they are barely affected by the actual arrangement of the H 2 molecules within the nanotube their contribution to the zeropoint energy can be considered constant and therefore irrelevant in the direct comparison.
3 Results and discussion

Adsorbate/nanotube interaction
The left-hand panel of Fig. 2 shows the total D 2 /CNT interaction potential as a function of the radial distance r between the molecule center-of-mass and the nanotube center, along with the electrostatic E elec , exchange-repulsion E exch-rep , induction E ind , and dispersion contributions E disp .Upon applying our pairwise potential model, we obtain the radial scans of interaction energies presented in the right-hand panel of Fig. 2 for nanotubes of increasing diameter (see also Table 1).As can be observed in Fig. 2 (left-hand panel), the D 2 /CNT attractive interaction is dispersion-dominated, with the net repulsive dispersionless contribution mainly determined by the exchange-repulsion term.As typically found, the exchangerepulsion grows exponentially as the molecule-surface distance decreases although such behaviour is somewhat smoothed out by the attractive electrostatic contribution.The induction term contributes very little at the potential minimum but considerably at the carbon cage.The potential minimum is located at the center of the narrow nanotube (diameter of 6.74 Å), where the adsorbates benefit from the dispersion interaction with carbon atoms at both sides of the carbon cage.However, upon increasing the nanotube diameter (see right-hand panel of Fig. 2), the dispersion becomes very small at the nanotube center and the potential minima shift towards the carbon cage (see also Table 2).For a nanotube of diameter of 27.21 Å, the value of the well-depth (À56 meV) is already close to that obtained for a graphene sheet and measured for the H 2 molecule adsorbed onto graphite (À51.7 meV from ref. 54).For a (5,5) nanotube, the energy at the potential minimum (À183) is about 15% larger than the well-depth of the potential proposed by Gordillo et al. 55 as a sum of Lennard-Jones-type C-H 2 interactions.From Fig. 2 it can also be observed that the potential model closely reproduces the ab initio results.
Fig. 2 Left-hand panel: Radial scan of the interaction energies between a single D 2 molecule and a short carbon nanotube (sCNT) of helicity index (5,5).Our pairwise potential model is compared with reference ab initio calculations using the SAPT(DFT) approach.Right-hand panel: Radial scan of the total interaction energies (full lines) and dispersion contributions (dashed lines) between a single D 2 molecule and carbon nanotubes with helicity indexes of increasing value.

Adsorbate wave-functions
As a model system in previous studies 12,13,56 and to simulate realistic experimental conditions, 7,8 we have performed adsorbate wave-function calculations considering a long carbon nanotube of helicity index (10,10) and diameter of 13.56 Å.Moreover, this diameter fits very well with the expected distance in the hcp structure of the bulk crystal.We focus here on the lowest-energy adsorbate (bosonic) wave-function with L = 0.The adsorption energies of (D 2 ) N and (H 2 ) N clusters to the CNT(10,10) nanotube are shown in Table 2.The one-particle density distributions D(r) are plotted in Fig. 3 as functions of the radial r distance from the adsorbate to the tube center.
As previously found for helium and molecular parahydrogen clusters confined by either an attractive dopant [57][58][59] or carbon nanotubes, 31 the total energy (see Table 2) can be expressed as N times the energy of the ground-state wave-function e 0 plus a small term hV A-A i accounting for the weak inter-adsorbate attractive interaction In fact, the adsorbate molecules tend to occupy the same nuclear orbital so that the radial densities differ only by the normalization factor N (see Fig. 3).These density distributions are extended by about 1 Å only.Along the radial direction, the molecular adsorbates are rather confined at the adsorbatenanotube potential minimum which, for the CNT(10,10) tube, is located at the lateral region close to the carbon wall with the peak at about 3.6 Å for D 2 .Notice that this value matches perfectly with the nearest-neighbour distance in the crystal hcp structure of deuterium at low temperatures (3.6 Å from ref. 60  and 61).The position of the peak for the H 2 isotope (about 3.55 Å) is also consistent with those determined by Rossi and Ancilotto 12 and Del Maestro and Boninsegni 13 for parahydrogen molecules adsorbed inside the (10,10) nanotube (about 3.70 Å from ref. 13).
A strong confinement at the adsorbate-nanotube potential minimum is also found in classical descriptions of molecular adsorption.This is consistent with the small zero-point energy for the adsorbate-nanotube vibration (about 7 meV) as compared with the well-depth (about À72.5 meV).In contrast, the D 2 -D 2 vibrational motion is characterized by an extremely large zero-point energy (2.20 meV for the pure D 2 dimer) as compared with the corresponding potential minimum (À3.065 meV), with an even larger value for the lighter H 2 dimer (2.70 meV).The associated pair density distribution is very delocalized and extends over more than 8 Å.This is clearly seen in Fig. 4 which presents plots of the pair densities D(R 12 ) as a function of the adsorbate-adsorbate distance for (D 2 ) 2 and (H 2 ) 2 dimers considering both isolated as well as confined clusters.We observe that the pair densities become more localized when going from isolated to confined clusters, and upon augmenting the mass of the hydrogen isotope.However, the carbon cage is barely affecting the adsorbate-adsorbate interaction of the confined dimers (by about À0.1 meV only, see Table 2).
Let us now consider the adsorbate wave-function dependence on relative inter-adsorbate degrees of freedom.Fig. 5 plots the density distributions of the (D 2 ) 2 and (H 2 ) 2 dimers as a function of the relative angle f 12 = f 2 À f 1 (left-hand panel) as well as the relative inter-adsorbate distance along the tube axis z 12 = z 2 À z 1 (right-hand panel).The peak at f 12 = 601 in the angular distribution of the (D 2 ) 2 dimer (left-hand panel) corresponds to a C 6 -type arrangement of the two molecules at the plane with a zero value of z 12 (right-hand panel), which is obviously consistent with the hcp structure in the extended system (see TOC graphics).Notice that the pair distribution for the lighter (H 2 ) 2 dimer is more extended and delocalized due to the larger zero-point energy.
From the right-hand panel of Fig. 5, it can be observed that the pair density distribution D(z 12 ) for the (D 2 ) 2 dimer is very broad and develops a shoulder at about the position of the potential minimum for the adsorbate-adsorbate interaction.The bimodal structure of the D(z 12 ) distribution also reflects the multi-configurational character of the adsorbate wave-function and a strong coupling between inter-adsorbate radial and angular vibrational modes due to the two-dimensional confinement induced by the carbon cage.Interestingly, the doublepeaked distribution in the (D 2 ) 2 dimer becomes disentangled in two pair distributions peaked at zero and 3.6 Å.Clearly, the density profiles along the tube axis are being shaped by the inter-adsorbate interaction.Similarly, notice that the angular distribution peaked at about 601 become more localized for the (D 2 ) 3 trimer (see left-hand panel of Fig. 5), with average of  This journal is © the Owner Societies 2017 the relative angles equal to 551 and 271.These values are close to those in the extended hcp structure (601 and 301).The adsorbate wave-function dependence on relative angles is also determined by the inter-adsorbate interaction.
Taken together, our adsorbate wave-function results show that the structuring of molecular deuterium clusters inside CNTs is determined by a delicate balance between attractive adsorbate-nanotube and inter-adsorbate interactions, with the extremely large zero-point energy corrections of the latter in a key role.This structuring is such that the attractive adsorbateadsorbate interaction is optimized, avoiding an increase of the zero-point energy cost for the adsorbate-adsorbate vibrational mode, but also taking advantage of the attractive adsorbatenanotube potential minimum.On one hand, the magnitude of the adsorbate-nanotube potential minimum (about À80 meV) is much larger than the inter-adsorbate counterpart (about À3 meV).On the other hand, as in doped helium and molecular hydrogen clusters, 58,59,62 the average adsorbate-adsorbate contribution to the total energy (about À1 meV, see Table 2) scales approximately with the number of adsorbate pairs (i.e., as N(N À 1)/2), while the one-particle term (about À65 meV) scales linearly with N. Therefore, a significant influence of the two-body term would have been expected for larger clusters than those considered here.Moreover, accurate adsorbate wave-function calculations of small doped clusters have indicated a cyclic arrangement of N (N o 6) hydrogen molecules around the dopant center-of-mass 58 which depart largely from the structuring under confinement by the carbon walls here revealed.
We notice that the findings in this work are consistent with our previous results on two-dimensional confinement of helium clusters, 31 where the main peak of the pair correlation function was located at a position very close to that measured for superfluid liquid helium 63 when considering two or three 4 He atoms. 31Thus, carbon nanotubes seem to enhance nuclear quantum effects in cold few-atom clusters at the nanoscale such as, for example, signaling the transition to a two-dimensional quantum fluid with just a few 4 He atoms.Very intriguing quantum phenomena have also been revealed when immersing carbon nanotubes in helium droplets such as an incomplete flooding of the quantum fluid. 29

Embedding approach
As illustrated in the TOC graphics, the hcp structure contains alternating C 6 and C 3 symmetry arrangements.To provide an evidence of the C 3 structure, let us now consider how the radial density distribution of a D 2 molecule is modified by two C 6 symmetry arrangements of D 2 molecules.For this purpose, as described in the second section, we calculate the wavefunction of a single D 2 in an effective potential incorporating the effect of the two hexagonal D 2 shells.As shown in Fig. 6, the addition of these shells as an embedding potential causes the shift of radial density peak from the lateral region at r = 3.6 Å to about r = 1.8 Å.The reasons of this behaviour is revealed in Fig. 6: at the tube center (z = 0), the repulsive barrier of the embedding potential is located at r = 3.6 Å as the peaks of the radial densities for the D 2 molecules from the hexagonal shells, while the largest attractive interaction between the embedded and the D 2 shells occurs at about r = 1.8 Å.This value differs from that of the exact hcp structure (for a = 3.6 Å) by 0.2 Å only for the C 3 symmetry arrangement.This deviation is small when considering that the D 2 molecules from the hexagonal shells are not fixed particles but display density distributions which, along r, extend by about 1 Å.
The effect of the embedding potential is also manifested in the binding energy of an embedded D 2 molecule (À53.21meV).Remarkably, this value differs by just 18% from the ''unembedded'' counterpart (À64.89,see Table 2).Hence, the reduced attractive D 2 -nanotube interaction upon getting further from the carbon cage in the C 3 symmetry structure is almost compensated by the attractive interaction with D 2 molecules in the adjacent C 6 shells.Once again, this is in contrast with the expected classical behaviour, i.e., the adhesion of D 2 molecules to the carbon walls whatever its number be, and highlights a delicate balance between D 2 -D 2 and D 2 -nanotube interactions.
To close our adsorbate wave-function analysis, we note that the packing of D 2 molecules on the triangular-type arrangement is also consistent with the central column of D 2 molecules distributed along the z axis is in the hcp structure (see TOC graphics).The existence of the central column has been reported by Rossi and Ancilotto 12 as well as Del Maestro and Boninsegni. 13As detailed in the ESI, † an effective model explains that the state of the dimer formed along the tube axis is extremelly weakly bound (by about À0.03 meV) for a D 2 molecule at a C 6 -type arrangement with r = 3.6 Å (i.e., the average position for a single D 2 inside the (10,10) nanotube).However, the formation of a cluster of three D 2 molecules in a C 3 arrangement inside the same nanotube makes the effective inter-adsorbate interaction much deeper, supporting much more strongly bound vibrational states.

Vibrational model
This section is dedicated to the vibrational analysis of single layers of D 2 molecules with either C 6 or C 7 symmetries inside a CNT(10,10) nanotube.Following the outline discussed in Section 2.4, we first build a total PES model by combining our pairwise potential model for the interaction between D 2 molecules and the carbon structure with an analytical expression for the D 2 -D 2 potential reported by Hinde. 52A convenient one-dimensional cut through this still complicated potential energy surface, a function of 2N À 3 vibrational degrees of freedom with N = 7, is obtained by scans over the energy with respect to the ''breathing mode'' of the deuterium cluster.Results are presented in Fig. 7.
We use these energy curves to derive an approximate analytical potential surface for planar motions of the D 2 molecules inside the tube and determine the zero-point energy in the corresponding modes from the Hessian matrix of this potential.Zero energy is set to the geometry where all D 2 molecules are at infinite distance from the tube and from each other.As can be seen, the confinement leads to approximately harmonic shapes of both potentials.The curves become steeper at larger distances, indicating that size extension is more costly than a compression of the D 2 arrangement due to the stronger repulsion effects from the wall of the carbon nanotubes than from inter-adsorbate interactions.The different packing of the D 2 molecules in either C 6 or C 7 symmetric structures has a significant impact on the curvature of the potential: adding an extra D 2 molecule to the hcp structure leads to a much steeper curvature.However, note that the minimum of energy is barely affected by the addition of an extra D 2 molecule, which shows that the total gain and costs of a particle addition are comparable to each other.A better overview is obtained if we compare electronic binding energies per D 2 molecule as presented in Table 3, which provides almost the same electronic energies per particle in both geometries.However, zero point effects need to be taken into consideration in the low energy regime (Fig. 7).
Assuming a simple model of interconnected springs (see ESI †), the potential energy of the confined D 2 molecules is fully defined by only two parameters, a spring constant f and a minimum distance r min between the interconnected D 2 molecules.Within this approach of coupled harmonic oscillators, justified by the small deviations from the optimum intermolecular distance, we obtain a force constant f = 798 cm À1 Å À2 and an equilibrium distance of 3.578 Å in the C 6 symmetry for CNT (10,10).In C 7 symmetry we obtain a force constant f = 1311 cm À1 Å À2 and an equilibrium spring distance of 3.627 Å. Diagonalization of the analytical Hessian with these parameters yields a set of vibrational frequencies o i , from which the total zero point energy of each system can be derived.For convenience, these zero point energies are then divided by the number of D 2 molecules and added to the corresponding electronic energies per particle listed in Table 3 for each symmetry.Such a comparison reveals that the C 6 symmetric arrangement, which emulates the hcp structure in the extended   tube, is now clearly preferred over a C 7 symmetric arrangement due to its lower zero point energy.

Conclusions
In summary, a new potential model for the hydrogen-nanotube interaction has been proposed which is based on accurate ab initio determinations through DFT-based symmetry adapted perturbation theory.This potential model is used in adsorbate wave-function calculations of molecular deuterium clusters confined by the CNT(10,10) nanotube, considering up to three molecules.Our results show a preference of hcp over other structural arrangements.They highlight that the intermolecular structuring is ultimately determined by a subtle balance between attractive adsorbate-nanotube and adsorbate-adsorbate interactions, with a prominent role played by the extremely large zeropoint energy associated to the latter.An embedding approach shows that the C 3 symmetry arrangement in the hcp structure can also be explained.Finally, a vibrational model demonstrates that the hcp structure is preferred over an heptagonal geometry due to quantum nuclear effects stemming from zero-point energy corrections.
The pairwise additive potential model applied in this work allows for an efficient treatment of molecular hydrogen physisorption on carbon-based materials, which is of special interest to researchers in the field of hydrogen storage and other applications concerning gas storage or gas separation.Moreover, as has been demonstrated, the accuracy of this potential model is high enough to even allow the discussion of quantum motion at very low temperatures within the frameworks of adsorbate wavefunction theory, which might become relevant in future studies of neutral molecules in weak confinement potentials for new technologies such as quantum information processing and molecular data storage.

Fig. 1
Fig. 1 Figure illustrating the atomic structural models for the H 2 /CNT system.Gray spheres represent carbon atoms while blue spheres stand for one H 2 molecule inside a carbon nanotube with diameter of helicity index(10,10).Cylindrical coordinates (r, f, z) of the H 2 center-of-mass are also indicated.The H-H bond is oriented along the z axis, with the bond length fixed to the vibrationally averaged rigid rotor value, r eq = 0.7508 Å from ref.23.

Fig. 3
Fig. 3 Radial density distributions D(r) of (D 2 ) N clusters and the (H 2 ) molecule inside the CNT(10,10) nanotube, normalized as Ð drDðrÞ ¼ N. The adsorbate-nanotube interaction potential is also shown as well as the position of the potential minimum (red-colored arrow).

Fig. 4 Fig. 5
Fig. 4 Pair density distributions D(R 12 ) as a function of the adsorbateadsorbate distance R 12 .

Fig. 6
Fig. 6 Radial density distributions D(r) of embedded and unembedded D 2 molecules inside the CNT(10,10) nanotube.The embedded molecule is confined by two symmetry C 6 arrangements of D 2 molecules located at z = AEc/2 (see TOC graphics).The corresponding distribution D(z) is peaked at z = 0 and extend up to z = AE0.5 Å.The embedding and adsorbatenanotube interaction potentials at z = 0 are also shown.

Table 1
Energies (E min , in meV) and positions (r min , in Å) of the potential minima for the interaction between a single D 2 molecule and carbon nanotubes CNT with helicity indexes (n,n) and diameter (d CNT , in Å) of increasing value This journal is © the Owner Societies 2017 Phys.Chem.Chem.Phys., 2017, 19, 28621--28629 | 28625

Table 2
Total energies and inter-adsorbate contribution hV A-A i of the lowest-energy bound states of molecular deuterium clusters in the inside of the CNT(10,10)nanotube.The bound-state energies of the isolated dimers are indicated in parenthesis

Table 3
Energies (in cm À1 ) per D 2 molecule in CNT(10,10), without and with zero-point energy (ZPE) corrections This journal is © the Owner Societies 2017