Pulling Simulations and Hydrogen Sorption Modelling on Carbon Nanotube Bundles

Abstract

Recent progress in molecular simulation technology has developed an interest in modernizing the usual computational methods and approaches. For instance, most of the theoretical work on hydrogen adsorption on carbon nanotubes was conducted a decade ago. It should be insightful to reinvestigate the field and take advantage of code improvements and features implemented in contemporary software. One example of such features is the pulling simulation modules now available in many molecular dynamics programs. We conduct pulling simulations on pairs of carbon nanotubes and measure the inter-tube distance before they dissociate in water. We use this distance to set the interval size between adjacent nanotubes as we arrange them in bundle configurations. We consider bundles with triangular, intermediate and honeycomb patterns, and armchair nanotubes with a chiral index from n = 5 to n = 10. Then, we simulate low pressure hydrogen adsorption isotherms at 77 K, using the grand canonical Monte Carlo method. The different bundle configurations adsorb great hydrogen amounts that may exceed 2% wt at ambient pressures. The computed hydrogen capacities are considered large for physisorption on carbon nanostructures and attributed to the ultra-microporous network and extraordinary high surface area of the configured models.

1. Introduction

Fossil fuels share almost the 85% of the global energy production. The rest is attributed to renewable sources and nuclear technology [1]. Presently, there is a strong policy push encouraging the transit of energy mix from fossil fuels to renewables. Important innovations and improvements are made on key technologies and relevant infrastructure. What should be surprising is that the share of non-fossil fuels to the total energy has been reduced the last decade [2]. This is due a growing aversion to nuclear power, exemplified in terms of energy safety, which negated the effect of growth in renewable sources. We have made only marginal progress on the transit of energy mix, which should be both a cause of concern and a focus for action. The evolution of fuels towards decarbonization, implies that hydrogen energy will be the last stepping stone. Hydrogen is the cleanest energy carrier and has a heating value three times higher than petroleum [3]. It can be produced from renewable sources, but it exists as low density gas under ambient conditions. Hydrogen storage is the main handicap concerning the implementation of hydrogen-based propulsion systems especially for on-board applications [4,5,6,7,8,9].

A direct method for hydrogen storage is the compression at high pressures. For instance, a recently introduced hydrogen-fuelled vehicle uses a Proton Exchange Membrane (PEM) fuel cell and a cylindrical vessel with 70 MPa $H_2$. Nevertheless, the safety risk of using on-board pressurised vessels in populated regions emerges a concern. A different method to store hydrogen is to liquefy it at 20 K. This approach requires considerable amounts of energy because of the very low storage temperature. There are also significant losses during both cooling and storing hydrogen, due to several so-called boil-off mechanisms [10,11].

A further possible way is the adsorption of hydrogen in microporous materials. Microporous materials have networks of pores with size less than 2 nm. Hydrogen molecules interact with the surface and are adsorbed inside the pores. The small size of the pores is essential for the adsorption for two reasons. First, the van der Waals forces that stem from the opposite walls may overlap at the center of the pore, likely resulting in high adsorption energies. Second, microporous solids often have great specific surface area values, and the total hydrogen capacity is proportional to this area. The relation of the capacity with the surface area holds because hydrogen is supercritical at usual adsorption conditions ($T>T_C=$ 33 K) and it is not expected to condense in wide pores (as $N_2$ does at 77 K). The adsorbed phase of a supercritical gas is attached to the surface and may grow strictly up to the thickness of a monolayer. Therefore, the adsorption capacity of supercritical gases is proportional to the number of available sites on the surface.

Different categories of microporous materials, such as carbons, zeolites and metal-Microporous carbons are lightweight materials and usually, they have large surfaces. As a consequence, they likely yield great hydrogen gravimetric capacities. Notably, carbon nanotubes were initially considered to be the solution to hydrogen storage problem but now seems clear that this is not the case [12,13]. Besides a setback on the relevant research, carbon nanotubes and their functionalized counterparts are now attracting a recursive interest as hydrogen sorbents [14,15,16,17,18,19,20,21]. This is primarily because the nanotubes are used as structured supports in advanced nanocomposites or can be embedded in membrane bilayers and mixed matrices [22,23,24,25,26,27,28]. Such systems, which have been extensively adopted in applications of biotechnology and biomedicine, are also investigated as candidate materials for energy storage.

The majority of theoretical work on nanotubes and their interactions with hydrogen has aged about 10 years [13,29,30,31,32,33]. Since then, molecular simulation technology has been massively advanced. Modern molecular dynamics programs support a rich set of calculation types, preparation and analysis tools and a flexibility to interplay with different software. Carbon nanotube modelling can be redesigned using contemporary computational workflows and take advantage of the new technology. For example, we can employ pulling simulations on pairs of nanotubes to configure high surface area bundle models of such nanoparticles and simulate hydrogen adsorption. Pulling simulations describe the process that places a molecular substance or a structure at increasing center of mass distance from a reference point with its position maintained by a biasing potential. Pull codes can now be routinely called by most molecular dynamics programs. The method is increasingly being applied to better understand a range of molecular interactions and allow one to describe the free energy profile along a particular reaction coordinate [34,35,36,37,38,39,40,41,42,43,44,45,46,47,48].

We conduct pulling simulations on model duplicates of armchair carbon nanotubes with a chiral index from n = 5 to n = 10. We embed the models in salt water and pull the one nanotube along a reaction coordinate using a force constant. We measure the inter tube distance, the moment before the pulled nanotube is being detached. We assign this distance to the side length of two closest neighbour nanotubes and build bundle configurations of such nanoparticles. We consider bundles with three different arrangements. Namely we build triangular, intermediate and honeycomb bundle configurations. We simulate low pressure hydrogen adsorption isotherms at 77 K on the models using the grand canonical Monte Carlo method. We attribute the discrepancies of the resulted isotherm curves to the porous characteristics of the structures such as the surface area and pore size distribution.

2. Materials and Methods

We generated six armchair nanotubes CNT(n,n) using the BuildCstruct (http://chembytes.wikidot.com/buildcstruct) script and by varying the chiral index from n = 5 to n = 10 [49]. We set the length of the nanotubes to 15 nm. We defined the bonded and the non bonded parameters of the atoms on the basis of the OPLS all-atom (OPLSAA) force field [50]. Specifically, we modelled the $s{p}^2$ carbons using the references of naphthalene and aliphatic carbons. We modelled the terminated hydrogens of nanotubes as benzene hydrogens. We performed the simulations using the GROMACS package, version 2018 [51].

We conducted pulling simulations using two identical armchair nanotubes. One of the nanotubes served as a reference while the other is placed at increasing center of mass (COM) distance from the reference. Before pulling, the two nanotubes were placed at a fixed distance with parallel orientations. We solvated the simulation box with simple point charge (SPC) water and we added 100 mM NaCl. We performed steepest descents minimization to remove the overlaps. We applied position restraints on the reference nanotube. The system was then equilibrated in two consecutive steps. The first phase involved a 50 ps simulation under a constant volume (NVT) and the second another 50 ps under constant pressure (NPT). In both equilibration steps the two nanotubes were coupled to separate coupling paths. We used weak coupling to maintain the temperature at 310 K in NVT and the pressure isotropically to 1.0 bar in NPT. Production molecular dynamics (MD) simulations were conducted for 100 ps in the absence of any restraints. For these runs we used the Nosé-Hoover thermostat to maintain temperature and the Parinello–Rahman barostat to isotropically regulate the pressure.

The final configuration of the last NPT run was used as starting configuration for the actual pulling simulation. We edited the simulation box and increased appropriately the dimensions to satisfy the minimum image convention and provide space for pulling simulations along the coordinate axis which, in this case, was the x-axis. As before SPC water was added to represent the solvent and 100 mM NaCl was present in the simulation cell. Equilibration was performed for 100 ps under the NPT ensemble, using the same control parameters described above. We set the first nanotube immobile by restraining the position of its carbon atoms. The other nanotube was pulled away along the x-axis over 500 ps, using a spring constant of 1000 kJ mol${}^{-1}$nm${}^{-2}$ and pull rate 0.01 nm ps${}^{-1}$. A final center of mass distance between the nanotubes of approximately 5.5 nm was achieved. The outcome of interest in the particular simulations was the distance, s, between the nanotubes the instant before the separation. The length, s, was then assigned to the length of the unit vector of triangular lattice template, which we used to build the bundle models.

We built three bundle configurations for each armchair nanotube. We used templates for triangular, intermediate and honeycomb arrangements. Each bundle was placed in a cubic simulation box. We edited the box edges at a distance at least 0.5 nm from the outermost nanotube of the bundle. To ensure that we have processed correctly the topologies, we performed a steepest descents minimization and then vacuum molecular dynamics simulation for 50 ps in NPT at 300 K. The vacuum molecular dynamics also served to make the box ready for adsorption simulations. In adsorption simulations, a fixed initial configuration is used, which is the simulation box containing the adsorbent without the presence of adsorbate molecules (i.e., in vacuum). This is equivalent to the adsorption experiments where the sample is degassed in high vacuum before the measurement. Some critical size metrics for the different bundle configurations are given in Table 1.

We conducted detailed pore size analysis for the different bundle configurations after the vacuum simulations. We computed the total pore volume (TPV) of the structures using Monte Carlo integration. That is we iteratively inserted a single helium molecule (${ϵ}_{He}=$ 10.22 K and ${\sigma }_{He}=$ 0.258 nm) at random coordinates inside the simulation cell and counted the times we get an attractive potential [52,53]. The integration was performed in ${10}^7$ iterations. We computed the specific surface area (SSA) of the bundles based on the rolling over method of Düren [54]. In this step we used nitrogen as a probe, modelled like a sphere with ${\sigma }_{N_2}$ = 0.33 nm and cross sectional area $A_{N_2}$ = 0.169 nm${}^2$. Finally, we computed the pore size distributions (PSD) of the bundles, using the poreblazer code [55]. The code resulted a histogram for the cumulative pore volume, $V_p\left(r\right)$, of the structure for a series of equally sized bins, $dr$. The PSD was given by the numerical differentiation $dVp\left(r\right)/dr$. Different codes related to the computational pore size analysis of microporous solids have been reported in the literature [56,57].

We simulated hydrogen adsorption in the nanotube bundles using the grand canonical Monte Carlo (GCMC) method. We applied full periodic boundary conditions on the simulation box. We sampled configurations by attempting random moves of hydrogen molecules under constant chemical potential, volume and temperature ($\mu$,V,T). The moves were molecular displacements, insertions and deletions. The probability to select each type of move was $1/3$. The acceptance or the rejection of a new configuration depended on the energy change that the move created in the system and the criteria of Metropolis [58]. We attempted ${10}^7$–${10}^8$ Monte Carlo steps and collected statistics at the last 40% of the total configurations. We performed the simulations at 77 K. Each simulated isotherm consisted of 110 fugacity points. The fugacity points were distributed exponentially in the range [0,100] kPa, being more concentrated at low (or more disperced at high) pressures. The full isotherm was submitted in a single file in the form of an array of simulation jobs, listing all the 110 fugacity steps. The isotherms were expressed as excess adsorbate mass. The excess particles $n_{ex}$, were given by $n_{ex}=n_{abs}-{\rho }_{bulk}{V}_p$, where $n_{abs}$ was the absolute number of hydrogen molecules inside the simulation box, ${\rho }_{bulk}$ was the bulk hydrogen density given from the Peng Robinson equation of state and $V_p$ was the total pore volume (TPV) of the bundle computed with (helium) Monte Carlo integration [54]. The excess particles were then converted to weight per cent capacities, $wt$ % = $M_{H_2}{n}_{ex}/(100$ $m_b)$. Where, $m_b$ was the molecular weight, in g/mol, of the bundle in the simulation cell and $M_{H_2}=2$ g/mol.

We used the Feynmann–Hibbs potential to compute the intermolecular interactions of hydrogen. This is a Taylor expansion of the classical (6,12) Lennard Jones potential, $U_{LJ}$, in respect to the distance, r, of the interacting particles [31,59].

$$U_{FH}\left(r\right)=U_{LJ}\left(r\right)+\frac{{\hslash }^2}{24\mu k_{B}T}{\nabla }^2\left(U_{LJ}\left(r\right)\right)$$

where ${\nabla }^2\left(U\left(r\right)\right)=U^{″}\left(r\right)+2U^{\prime }\left(r\right)/r$, T = 77 K, is the temperature of the GCMC simulations, $k_B$ is the Boltzmann constant, ℏ is the reduced Planck constant and $\mu =(m_i^{-1}+m_j^{-1})$ is the sum of the inverse masses of the interacting particles i and j. Special feature of the Feynmann–Hibbs potential is that it incorporates the molecular masses in the calculations. This is highly critical when simulating light mass adsobates as it provides a coarse description for quantum behaviour which can be computed inline. This attribute makes the expression of Feynamnn–Hibbs (FH) potential able to calculate explicitly the interactions of hydrogen isotopes. In this regard, FH has been extensively used in quantum sieving simulations [60,61,62,63]. We used ${ϵ}_{H_2}$= 36.7 K, ${\sigma }_{H_2}$= 0.296 nm, $m_{H_2}$= 2 g/mol for molecular hydrogen and ${ϵ}_C$ = 28 K, ${\sigma }_C$ = 0.34 nm, $m_C$ = 12 g/mol for the carbon atoms.

3. Results

We conduct pulling simulations on duplicate models of armchair nanotubes with a chiral index from $n=5$ to $n=10$. The two nanotubes are placed in close contact and with parallel orientation. We restrain the position of the first nanotube and pull the second along the x axis. The applied force and the distance between the nanotubes over time are shown in Figure 1. By pulling on the nanotube, force builds up until a breaking point is reached, at which time the van der Waals interactions of the two nanoparticles are disrupted allowing the second nanotube to move. After this point, the force drops abruptly and levels off to zero. The time required to achieve the maximum force, as well as the value of the maximum force depends on the chiral diameter of the nanotube. The narrowest nanotubes request weak pulling forces, because of the high convex curvature of the surfaces. The nanotubes with chiral index $n=7,9$ and 10, require comparable pulling forces and times to begin the displacement. Beyond the breaking point and for a short time range, the velocity of the pulled nanoparticles is fast. After the applied force drops to zero, the nanotubes move at a constant speed that is equal to the pulling rate, 0.01 $nm/ps$. The distance of the two nanotubes before the breaking point, is assigned to the side length, s, of the adjacent nanotubes in a bundle configuration. That is the distance between the closest neighbour nanotubes in the arrangement.

In Figure 2, we present a set of configurations from the pulling simulations. The left nanotube (black color) is constrained at a fixed position. The right nanotube (red) is being pulled away by applying a force at the center of mass. When the force becomes sufficient to overcome the restoring forces, the red nanotube moves towards the specified coordinate. The position and the orientation of the pulled nanotube is drawn in four consecutive time steps in a range of 0.05 ns after the moment of separation. Although the reaction coordinate and the driving force constant is the same, the nanotubes are being pulled in different pathways. The time range of the illustrated configurations is highlighted in the bottom panel of Figure 1.

In Figure 3, we present three bundle configurations of the armchair carbon nanotube CNT(9,9). The triangular bundle corresponds to the most densely packed tiling of the nanoparticles. The intermediate bundle is configured by deleting two vertically aligned nanotubes from the triangular arrangement. The honeycomb bundle is structured by deleting two more nanotubes from the central horizontal axis of the intermediate bundle. The latter bundle consists of two hexagonal tiles of nanotubes. We use the arrows in Figure 3 to describe three critical size metrics of the bundle configurations. These are the side, s, of the adjacent nanotubes, the distance of the second closest nanotubes (width), $w=\sqrt{3}s$ and of the third closest nanotubes (height), $h=2s$, in the hexagonal tile. As noted, the side, s, is computed using pulling simulations on model duplicates of the specified structures. The width and the height are measured from the center of mass (COM) of the nanoparticles. The term, $\sqrt{3}$, comes from $sin\left({60}^{°}\right)$. In porosity measurements, it is useful to express the interparticle lengths and wall-to-wall distances. Therefore, if d is the chiral diameter of the nanotube, we use the expression $w^{\prime }=w-d$ for the width and $h^{\prime }=h-d$ for the height of the hexagonal tiles. Accordingly, the effective diameter of the nanotubes is given by $d^{\prime }=d-{\sigma }_C$, where, ${\sigma }_C$, is the Lennard Jones diameter of carbon. We list the pore size parameters of the different bundles in Table 1.

In Figure 4, we show the computed pore size distributions for the different bundles. We consider structures with triangular, intermediate and honeycomb arrangements of nanotubes. The peaks on the distributions indicate the diameter of the pores. The height of the peaks indicates the contribution of the pores to the total pore volume (TPV). We detect three distinct peaks on the distributions which correlate with some principle size metrics of the bundles. We use three coloured barriers inside the panels to spot these sizes. The sizes are the effective diameter, $d^{\prime }$, of the nanotubes, the width, $w^{\prime }$ and the height $h^{\prime }$ of the hexagonal tiles. The values are listed in Table 1. Apart from the main peaks on the distributions, other value fluctuations are due to the convex surface of the neighbouring nanotubes, which may intervene the space of $w^{\prime }$ or $h^{\prime }$ pores. The first peak of the distributions is related to the effective diameter, $d^{\prime }$, of the nanotube. The height of the first peak increases with the chiral index of the armchair nanotube since the internal volume also increases. The distributions of CNT(5,5) do not show relevant first peak, because the effective diameter of CNT(5,5) is smaller than the size (0.4 nm) of the finite cubes we used in the grid. The second and the third peak of the distributions locate close to the sizes, $w^{\prime }$ and $h^{\prime }$ respectively. Most of the cavities in the triangular bundles have a diameter comparable to $w^{\prime }$. Deleting the nanotubes from a triangular bundle, produces configurations with more open space than the initial structure. In this regard, in intermediate and honeycomb bundles the pores with size $w^{\prime }$ are reduced in number and new cavities with size $h^{\prime }$ are configured.

We note that the distributions peak slightly more on the left than the coloured barriers. This is expected due to the vacuum simulations, which we performed after configuring the bundles. The lengths, $w^{\prime }$ and $h^{\prime }$, are computed with the nanotubes immersed in water. These variables are proportional to the side, s, so that they are affected by a possible formation of water layer between the tubes during the pulling simulations. With vacuum simulations the adjacent nanotubes may approach each other closer than s, which is depicted to the resulted distributions. The vacuum simulations change the geometric properties of the as made bundles, similar to what is experimentally anticipated when evacuating the solid or extracting the solvent from the micropores.

In Figure 5, we present the isotherms of hydrogen adsorption simulations on the different CNT models and bundle configurations at 77 K. The isotherms refer to the excess $H_2$ uptake using the formula for $n_{ex}$ and the total pore volume, $V_p$, of the structures, listed in

Table 2. Total pore volume (helium) of the CNT bundles arranged in triangular, intermediate and honeycomb patterns and the change (%) in respect to the pore volume of the triangular bundle.

	Total (Helium) Pore Volume (cm${}^3$/g)			Change (% from Triangular)
	Triangular	Intermediate	Honeycomb	Intermediate	Honeycomb
CNT(5,5)	1.37	1.67	2.08	22.3	52.2
CNT(6,6)	1.52	1.85	2.27	21.7	49.7
CNT(7,7)	1.50	1.80	2.21	20.1	47.5
CNY(8,8)	1.68	2.03	2.50	21.4	49.2
CNT(9,9)	1.73	2.10	2.58	21.1	48.6
CNT(10,10)	1.80	2.17	2.66	20.5	47.8

. Then we convert the output to the weight per cent, ($wt$ %) values, using the molecular weights of the bundles, $m_b$, listed in

Table 3. Molecular weight of the CNT bundles arranged in triangular, intermediate and honeycomb patterns.

	Weight (g/mol) in the Cell
	Triangular	Intermediate	Honeycomb
CNT(5,5)	23,800	20,400	17,000
CNT(6,6)	28,560	24,480	20,400
CNT(7,7)	33,320	28,560	23,800
CNT(8,8)	38,080	32,640	27,200
CNT(9,9)	42,840	36,720	30,600
CNT(10,10)	47,600	40,800	34,000

. In each panel, the three adsorption curves almost coincide at low pressures. This denotes that the different bundle arrangements of a single nanotube have similar micropore content. Hydrogen is either adsorbed inside the cavity of the nanotube or at the open space between the adjacent nanotubes. Both of these confinements are less than a few molecular diameters wide. The pore network of the CNT(5,5) and CNT(6,6) bundles consists exclusively of sub-nanometer pores. Such ultra microporous structures favour the adsorption of hydrogen at extremely low pressures. For instance, the uptake of $H_2$ in the bundles of CNT(5,5) and CNT(6,6) is ≃0.5% wt at 0.01 kPa. By increasing the diameter of the nanotube, the pore network of the bundles shifts to wider sizes. This results in the substantial decrease of $H_2$ uptake at low pressures. Near ambient pressures, hydrogen uptake is proportional to the specific surface area of the samples. The surface area values of the different bundle configurations are listed in

Table 4. Surface area of the CNT bundles arranged in triangular, intermediate and honeycomb patterns and the change (%) in respect to the surface area of the triangular bundle.

	Surface Area (m${}^2$/g)			Change (% from Triangular)
	Triangular	Intermediate	Honeycomb	Intermediate	Honeycomb
CNT(5,5)	1601.25	2081.72	2364.15	30.0	47.6
CNT(6,6)	1750.66	2221.26	2522.03	26.9	44.1
CNT(7,7)	1874.01	2291.26	2566.58	22.3	37.0
CNT(8,8)	2092.28	2609.34	2878.45	24.7	37.6
CNT(9,9)	2201.74	2679.08	2937.47	21.7	33.4
CNT(10,10)	2262.70	2760.00	3005.34	22.0	32.8

. The honeycomb and intermediate bundles have larger surface areas and adsorb greater densities than the triangular bundles. The surface area of the bundles increases with the chiral diameter of the nanotube. Accordingly, the corresponding bundle structures of the wider nanotubes adsorb greater densities at 100 kPa. In some cases, the computed hydrogen uptake at ambient pressures and 77 K may exceed 2 % wt.

Table 1. Size variables of bundle models for different armchair carbon nanotubes. The side, s, is computed with pulling simulations. The width $w^{\prime }$ and the height $h^{\prime }$ of the bundles are the wall to wall distances of the variables, w and h, shown in Figure 3.

	Diameter	Diameter (eff.)	Side	Width	Height
	$\mathit{d}$  $\left(\mathbf{nm}\right)$	${\mathit{d}}^{\prime }=\mathit{d}-{\mathit{\sigma }}_{\mathit{C}}$	$\mathit{s}$  $\left(\mathbf{nm}\right)$	${\mathit{w}}^{\prime }=\sqrt{\mathbf{3}}\mathit{s}-\mathit{d}$	${\mathit{h}}^{\prime }=\mathbf{2}\mathit{s}-\mathit{d}$
CNT(5,5)	0.678	0.338	0.963	0.991	1.249
CNT(6,6)	0.814	0.474	1.106	1.103	1.399
CNT(7,7)	0.949	0.609	1.180	1.095	1.411
CNT(8,8)	1.085	0.745	1.372	1.291	1.658
CNT(9,9)	1.221	0.881	1.485	1.352	1.750
CNT(10,10)	1.356	1.016	1.595	1.407	1.835

Table 1. Size variables of bundle models for different armchair carbon nanotubes. The side, s, is computed with pulling simulations. The width $w^{\prime }$ and the height $h^{\prime }$ of the bundles are the wall to wall distances of the variables, w and h, shown in Figure 3.

	Diameter	Diameter (eff.)	Side	Width	Height
	$\mathit{d}$  $\left(\mathbf{nm}\right)$	${\mathit{d}}^{\prime }=\mathit{d}-{\mathit{\sigma }}_{\mathit{C}}$	$\mathit{s}$  $\left(\mathbf{nm}\right)$	${\mathit{w}}^{\prime }=\sqrt{\mathbf{3}}\mathit{s}-\mathit{d}$	${\mathit{h}}^{\prime }=\mathbf{2}\mathit{s}-\mathit{d}$
CNT(5,5)	0.678	0.338	0.963	0.991	1.249
CNT(6,6)	0.814	0.474	1.106	1.103	1.399
CNT(7,7)	0.949	0.609	1.180	1.095	1.411
CNT(8,8)	1.085	0.745	1.372	1.291	1.658
CNT(9,9)	1.221	0.881	1.485	1.352	1.750
CNT(10,10)	1.356	1.016	1.595	1.407	1.835

4. Conclusions

We discuss how to produce carbon nanotube bundles in saline solution and study hydrogen adsorption properties after solvent removal and evacuation. We interplay a relatively new simulation technique, namely the pulling simulation, with an old one, the grand canonical Monte Carlo. These techniques have modules with different dependencies that we encode and execute in a single workflow. We conduct pulling simulations on pairs of nanotubes and configure high surface area bundles of such nanostructures to simulate hydrogen adsorption. We use the pull codes to compute appropriate interstitial distances for six armchair carbon nanotubes. We detect comparable puling forces and times for the nanotubes to detach. As they are being pulled, the nanotubes follow separate dissociation pathways. We measure the interstitial distances the instant before the nanotubes detach and use them as unit sizes in bundle templates. We consider three templates, each with different porous characteristics. We perform vacuum molecular dynamics on the resulted bundles which slightly change the initially set interstitial distances. All the resulted bundle structures are ultra microporous with high specific surface area. The surface area and the micropore content are crucial parameters for the materials performance for hydrogen physisorption. We simulate hydrogen adsorption on the configured bundles at 77 K. We compute and analyse the isotherms on the basis of the porosity and the chiral index of the nanotubes.