Importance of Van der Waals Interactions in Hydrogen Adsorption on a Silicon-carbide Nanotube Revisited with vdW-DFT and Quantum Monte Carlo

Density functional theory (DFT) is a valuable tool for calculating adsorption energies toward designing materials for hydrogen storage. However, dispersion forces being absent from the local/semi-local theory, it remains unclear as to how the consideration of van der Waals (vdW) interactions affects such calculations. For the first time, we applied diffusion Monte Carlo (DMC) to evaluate the adsorption characteristics of a hydrogen molecule on a (5,5) armchair silicon-carbide nanotube (H2-SiCNT). Within the DFT framework, we benchmarked various exchange-correlation functionals, including those recently developed for treating dispersion or vdW interactions. We found that the vdW-corrected DFT methods agree well with DMC, whereas the local (semilocal) functional significantly over (under)-binds. Furthermore, we fully optimized the H2-SiCNT geometry within the DFT framework and investigated the correlation between the structure and charge density. The vdW contribution to the adsorption was found to be non-negligible at ∼1 kcal/mol per hydrogen molecule, which amounts to 9–29% of the ideal adsorption energy required for hydrogen storage applications.


INTRODUCTION
Hydrogen energy is a promising energy resource for reducing greenhouse gas emissions. 1−3 To realize the industrial use of hydrogen energy, particularly in the transportation sector, one of the most important developmental challenges is addressing the related storage issuessafety and capacity. 2 Several material-based strategies to store hydrogen have been proposed, which involve choices such as the form of the stored hydrogen (physical vs. chemical storage) and the structure of the storage material (e.g., nanostructures and metal hydrides). To achieve adsorption-based room-temperature storage, the ideal interaction energy between the stored hydrogen and storage material has been estimated as ∼3.5− 11.5 kcal/mol. 4,5 At the moment, this technology is limited to low-temperature storage because the interaction is too weak to resist, being overpowered by thermal energy at the desired higher temperatures. Computational material design would be immensely helpful in further exploring appropriate storage materials from the massive material space, provided hydrogen adsorption energies on candidate materials can be accurately predicted at reasonable computational costs.
Silicon-carbide nanotubes (SiCNTs) are a typical nanostructure studied for the above purpose owing to their enhanced molecular interactions compared with those of the structurally related (and more common) carbon nanotubes (CNTs). This feature has been linked to the SiCNT polarized surface originating from the Si−C bonds. 6−9 Despite this, density functional studies based on local (LDA), semilocal (GGA), and hybrid (B3LYP) exchange-correlation (XC) functionals on pristine nanotubes estimated the adsorption energy of hydrogen between 0.7 and 1.98 kcal/mol, 10−13 which is still too small for storing hydrogen at the desired ambient temperature. This has motivated further investigations into doping schemes with dopants such as alkali and transition metals 12,14 or with vacancies. 11,15 These studies, however, did not incorporate van der Waals (vdW) corrections into the conventional DFT scheme, and thus they may have underestimated the true adsorption potential. From another viewpoint, the basis sets adopted for DFT calculations also matter, since Gaussian and numerical basis sets can give rise to overestimated interaction energies owing to the basis set superposition error. 16 Thus, the true adsorption potential of SiCNTs remains unclear, even when the results seem plausible and are consistent with the experimentally observed higher hydrogen uptake than that of CNTs. 17,18 In this study, we demonstrate the importance of incorporating vdW interactions for a quantitative description of H 2 adsorption on SiCNTs and related systems. We benchmarked leading vdW-corrected functionals for DFT against diffusion Monte Carlo (DMC) and conventional XCs in terms of binding curves and charge densities. DMC 19 is a stochastic method in which an "exact" many-body Hamiltonian (without any approximation) is used to project a set of electron positions into their ground-state configuration, allowing dispersion forces and dynamical correlations to be rigorously captured without the use of additional corrections. In this approach, DFT is only used to construct the initial Slater determinant from its Kohn−Sham orbitals and otherwise has only minimal effect on the final results. 20 DMC has been proved to be highly accurate for various noncovalent systems 20−31 and is comparable to the best correlated methods, while offering a better computational cost and minimal systematic error. 20 Finally, we investigated the fully optimized H 2 -SiCNT geometries for our benchmark set of XC functionals and discussed their relationship to the charge densities at a fixed geometry, along with the implications of their difference in terms of the charge densities.

COMPUTATIONAL DETAILS
2.1. Structural Modeling. We selected the type-I (5,5) armchair SiCNT with 20 atoms per primitive cell. Consisting only of Si−C bonds, type-I was found to be more stable than those containing Si−Si and C−C bonds. 7,32,33 The basic structure was built using a generic nanotube generation tool 34 and later optimized in DFT. A 14 × 14 × a Å simulation cell was used to minimize spurious interactions with periodic images, and the periodic unit length (a) was optimized along with the other geometries. For benchmark, we considered only the vertical mode of a single hydrogen molecule over the hollow site of the SiCNT (Figure 1). We discuss the effect of performing full optimization of hydrogen orientation in the later part of the results. The distance was varied to obtain a binding curve for which the energies were least-square-fitted to a Morse potential. The binding energy corresponding to the energy change from the free to the bound configuration is defined as Each component on the right-hand side was computed independently using the same simulation cell. Since the diameter of the SiCNT is not singular due to lengthwise buckling, we define r as the distance from the hydrogen molecule to the center of the nanotube subtracted by the largest diameter within the primitive cell. 2.2. Density Functional Theory (DFT). The PWSCF binary in the QUANTUM ESPRESSO 35 package was used for all DFT calculations based on plane wave basis sets with pseudopotentials. The kinetic energy cutoff and k-point grid size were converged at 150 hartree and a 1 × 1 × 6 Monkhorst−Pack grid, 36 respectively, to achieve chemical accuracy. To reduce the timestep error in the later DMC stage, the nonsingular energy-consistent pseudopotential reported by Burkatzki et al. 37 was utilized for all calculations. All geometry optimizations were performed within PBE with total force and energy thresholds of 10 −10 and 10 −4 a.u., respectively.
We compared the vdW approaches DFT-D2, DFT-D3, Tkatchenko−Scheffler (TS), Exchange-Dipole Model (XDM), vdW-DF2, and rVV10; 38−45 the first four are based on an explicit R −6 pairwise potential operating on atomic coordinates, whereas the latter two have an exchange-correlation potential modified by the addition of a nonlocal correlation term. The additional term in the exchange-correlation potential allows for changes in the charge density, providing more information for analyzing the binding formation. To make a comparison with previous studies, 10−13 we contrast their performance relative to the LDA from Perdew−Zunger (PZ) 46 and GGA from Perdew−Burke−Ernzerhof (PBE). 47 2.3. Diffusion Monte Carlo (DMC). QMCPACK 48 was used for all quantum Monte Carlo (QMC) calculations, the main DMC calculations, and for the variational Monte Carlo (VMC) calculations used in preparation of the trial wavefunction. The trial wavefunction was a Slater−Jastrow type, comprising a single determinant with PBE-DFT orbitals and a Jastrow factor with optimized parameters at the VMC level. Parameters consisting of one-and two-body interaction terms were used, each in the form of B-splines with 10 optimizable parameters for each atom type. For efficient statistical accumulation, Jastrow parameters were optimized by minimizing VMC energy using a hybrid method mixing the linear method algorithm 49 and an accelerated descent method as described by Otis and Neuscamann. 50 The orbitals were preprojected into B-splines 51 to increase the computational efficiency. The pseudopotential parts were evaluated with the T-move scheme. 52 The timestep and finite size errors were eliminated through linear extrapolation, first by the timestep and then by the real space supercell size. The DMC timesteps were 0.0025, 0.01, and 0.04 a.u, and the supercells were 2, 6, and 8 times duplicates of the primitive unit cell in the cylinder length, each with twist-averaged boundary conditions 53 on regular grids of 1 × 1 × 8, 1 × 1 × 3, and 1 × 1 × 2, respectively. The mixed boundary conditions were applied with open boundaries on the 14 Å sides and the target walker population was set to 4096. In the Supporting Information, these potential biases in DMC are discussed and proved to be small enough for the purpose of this study.

RESULTS AND DISCUSSION
3.1. Structural Properties. We found that the geometry of the (5,5) SiCNT is relatively unaffected by the choice of the XC functional (Table 1). This result is unsurprising due to its construction of only single covalent bonds involving s and p electrons. The obtained bond lengths were also consistent with those reported in prior hybrid works, 32,54 albeit slightly shorter than those obtained utilizing cluster models. 12,33 It is not unusual to have a slight distortion near cluster terminations that does not exist in periodic models such as the one used in this work. Thus, we conclude that the conventional XCs are sufficient to describe the SiCNT model used in this work.

H 2 Adsorption on SiCNT.
The main interest is how well the adsorption curves from the vdW corrections agree with the DMC results and among themselves. As shown in Figure 2, it is clear that PZ (PBE) severely over (under)-binds the DMC target value by 0.6 (1.2) kcal/mol. This result is consistent with the established behavior of LDA and GGA in vdW-dominated systems. The nonempirical vdW-DF2 and rVV10 functionals return virtually identical binding energies but their H 2 separations differ by ∼0.2 Å. Interestingly, the D2 energetics agree more with the nonempirical methods relative to D3. This was accompanied by a slight underestimation of the H 2 separation distance, which was corrected by D3. Considering they both start from PBE energies, this is satisfactory. In general, all vdW corrections underbind relative to DMC but with more reasonable adsorption minima. This suggests that while the geometries derived from these functionals can generally be trusted, a more careful consideration of their energetics is necessary. Overall, PBE + TS produces the closest binding energy to the DMC, while the H 2 distance is best reproduced by PBE + D3. LDA is known to produce spurious covalent bonding between noncovalent molecules due to the self-interaction error. 22,26 This gives rise to ca. 0.6 kcal/mol overbinding. GGA improves the self- The subscripts "p" and "d" denote perpendicular and diagonal bonds relative to the cylinder axis, respectively. interaction error, but at the cost of weak intermolecular interactions as the vdW interaction is not inherently accounted for. To illustrate this point, we plotted the charge density difference between the whole and isolated components, Δρ(r) = ρ SiCNT+H 2 (r) − ρ SiCNT (r) − ρ H 2 (r), in Figure 3. The charge accumulation between the H 2 molecule and SiCNT, i.e., spurious covalent bonding, is prominent for LDA-PZ as shown in Figure 3(a). This feature is not present in the GGA-PBE data, as shown in Figure 3(b), and is instead replaced by a much weaker redistribution of charge toward the SiCNT surface. Note that while the pairwise corrections (D2, D3, TS, and XDM) properly address the binding curve, they do not affect the charge density at a given geometry. In contrast, vdW-DF2 and rVV10 achieve their corrections by deforming the charge density (as a side effect of changes at the wavefunction level) using nonlocal perturbations in the correlation integral. As shown in Figure 3c,d for vdW-DF2 and rVV10, respectively, there is a slight dip in charge density in the shape of a bow; this is associated with the existence of noncovalent-type interactions, 55 although it is much weaker in the rVV10 case. Interestingly, the rVV10 distribution is much closer to the PZ distribution as opposed to the distinct distribution of vdW-DF2. The lack of proper vdW interactions in PZ and PBE results in misestimation by more than 0.5 kcal/mol, with further potential for errors arising from incorrect geometry predictions. Since the target adsorption energy is 3.5−11.5 kcal/mol, this misestimation is not negligible. We thus claim that the inclusion of vdW corrections in the XC functional is essential for the quantitative evaluation of adsorption energies, particularly in the vdW-heavy H 2 -SiCNT system, and more generally for further computational exploration of storage materials ( Table 2).

Fully Optimized Geometries.
In Section 3.2, we varied only the distance between H 2 and SiCNT for our benchmark purpose of evaluating the vdW interactions. Here, we consider the fully optimized geometries to evaluate a more realistic adsorption energy. Table 3 gives the structural and energetics information obtained from all of the DFT methods adopted in the present study. Similar to the findings discussed in Section 3.2, the vdW-corrected functionals (vdW-DF2 and rVV10) clearly give highly similar trends, while LDA (GGA) over (under)-binds, which is well-known for noncovalent systems. 26 Within the framework of the vdW-corrected functionals, they agree well with each other in terms of the structural parameters and adsorption energy. Compared to the vdW-corrected functionals, PBE gives poor description of all of the structural and energetics properties. In contrast, LDA-PZ gives a surface angle θ s closer to the nonlocal vdW corrections (vdW-DF2 and rVV10), but looking at the other properties, they are significantly different from each other; therefore, we may regard this coincidence in PZ as being accidental. These structural properties are associated with the charge densities as shown in Figure 3. Overall, LDA and the two vdW-DFT functionals (vdW-DF2 and rVV10) give different shapes than PBE. As the charge density at the bonding region between H 2 and SiCNT increases (LDA > rVV10 > vdW-DF2), R H 2 and θ s both decrease (LDA < rVV10 < vdW-DF2). As mentioned in Section 3.2, the LDA bonding is not noncovalent but spurious. It is evident from the correlation between the charge density and structure that the LDA overbinding correlates to the highest charge density at the bonding region.
Herein, we compare the pairwise vdW corrections to the vdW-corrected functionals. In particular, the surface angle θ s   values obtained from the "empirical" pairwise ones (PBE-D2/ D3) improve the PBE value, but they are still smaller than those from the vdW-corrected functionals. This can be attributed to the fact that the empirical corrections never deform the charge densities (or equivalently wave functions) self-consistently as the total energies just decrease due to the empirical pair potentials. In other words, the charge densities of the pairwise vdW corrections (PBE-D2/D3) are the same as those of their underlying XC functional (PBE). In contrast, since the other pairwise vdW corrections (PBE-TS and -XDM) can take the local chemical environment into account, they significantly improve the θ s values compared to the PBE value, though they still underestimate slightly compared to the vdWcorrected functionals. As for the other properties, PBE + TS overestimates the adsorption energy (E ads ) and underestimates the bond length (R H 2 ), while PBE + XDM give E ads and R H 2 closer to the vdW-corrected functionals (R18). Considering the above findings, we may conclude that the self-consistent charge deformation is important, and pairwise vdW correction alone is insufficient for an accurate structure determination.

CONCLUSIONS
We performed DMC and various DFT simulations to evaluate the adsorption energies of a hydrogen molecule on a SiCNT with a (5,5) armchair structure. Recently developed XC functionals designed to reproduce vdW interactions (PBE + D3, vdW-DF2, and rVV10) and conventional XC functionals (PZ, PBE) were compared to DMC as a reference. Overall, all of the vdW-corrected XC functionals agree well with DMC, whereas PZ (PBE) over (under)-binds. The self-consistent nonlocal correlation functionals, vdW-DF2 and rVV10, give almost the same adsorption energies. Differences in the structural properties were found to closely correlate with differences in the charge density distribution. A higher charge density in the bonding region leads to a shorter distance between H 2 and SiCNT and larger surface angle. The magnitude of the vdW interaction was estimated to be ca. 1.2 kcal/mol, which corresponds to 9−29% of the ideal adsorption energy for hydrogen storage. This finding implies the importance of vdW corrections within the framework of DFT. We thus conclude that protocols based on vdWcorrected XC functionals will advance the computational investigation and exploration of storage materials in the near future. In addition, our results support previous findings that pristine SiCNT alone is unlikely to provide a sufficient binding for realizing ambient temperature hydrogen adsorption; thus, further research studies into doping and other surface modifications should be pursued to use it as the basis for a hydrogen storage system.