The Role of Cobalt Clusters (Con, n = 1–5) Supported on Defective γ–Graphyne for Efficient Hydrogen Adsorption: A First Principles Study

In this theoretical work, density functional theory calculations show the effect of small cobalt clusters (Con, n = 1–5) adsorbed on pristine γ‐graphyne (γ‐GY), and modify N‐doped γ‐GY monolayers (GYNs‐def). Different geometrical configurations are assessed with the adsorption energy, charge transfer, and density of states. The system with vacancy defects shows a large adsorption energy (19.96 eV) for the Co5 cluster. This behavior may be associated to the overlapping of the electronic state contributions between cobalt and carbon atoms in the valence states. This indicates that the Co5 cluster could be deposited on N‐doped γ‐GY monolayers (Con@GYNs‐def). The lowest‐energy systems are evaluated to estimate the strength of the interaction with hydrogen molecules (xH2, where x = 1–5). According to the adsorption energy values, the modified γ‐GY monolayers are allowed to be a suitable support material to capture H2 molecules via the small Con clusters. The hydrogen retention capacity for the supported cobalt atoms corresponding to the lowest‐energy configurations and larger systems are evaluated by using molecular dynamics simulations with the Born–Oppenheimer approximation. The role played by defects in the GYNs‐def monolayers is important, since the Con clusters remain attached to the vacancy with the absence of surface diffusion. This study may represent a guide to tailor novel nanostructures based on small cobalt clusters supported on graphyne monolayers modified to be applied in H2 adsorption.


Introduction
The increasing demand for fossil fuels implies a rapid growing of air pollution and subsequent effects to be associated with human health issues, such as respiratory diseases. [1,2] Therefore, in the scientific community it is a top priority to find new alternatives to replace such fossil fuels. Hydrogen has been considered as a potential fuel to support sustainable and clean energy development, [3,4] due to the high energy density of 33.3 kW h kg −1 . [5] Even though hydrogen storage represents a safe and efficient methodology, it also represents a challenge for the development of novel materials. [6] An important feature to be considered is that the adsorption energy required to obtain a reversible hydrogen storage should be in the range from 0.2 to 0.7 eV per H 2 . [7][8][9][10] One of the best candidate materials to be used as a hydrogen storage material is the carbon-based nanostructures, such as fullerenes, [11] carbon nanotubes, [7,12] graphene, [13] and graphynes. [9,14] One the most attractive 2D carbonbased material, which has been widely assessed and scrutinized as a potential hydrogen storage media is the -graphyne ( -GY) structure. The -GY was proposed by Baughman et al., [15] and exhibits a geometry formed of a single atomic layer with an arrangement of carbon atoms with sp 1 -and sp 2 -hybridizations. A peculiar characteristic rises upon the hybridization of sp 2 +sp 1 +sp 1 +sp 2 bonding, resulting in the formation of acetylene linkages (-C≡C-). Such bonding is connected by six-membered rings. [15][16][17][18] In recent years, the -GY synthesis has been achieved via a mechanochemical reaction using CaC 2 and PhBr 6 as precursors. [19][20][21][22][23] The intrinsic cavity of -GY structure is very attractive to metal adatoms which prefer to interact in the plane of metal-decorated -GY. Therefore, metal atoms, including alkaline metals (AM), alkaline-earth metals (AEM), and transition metal (TM) atoms, have been explored by experimental and theoretical works. [8,9,[24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42] Currently, the -GY structure has been used as an ideal support to reduce the amount of TM atoms. Chen et al. [43,44] reported that pristine -GY structure is sufficiently stable to be used as a support of small Au and Ag clusters. They found that the TM clusters prefer to be adsorbed on the hollow site over the acetylenic ring, in which an exothermic process is involved. Meanwhile, some 2D carbon-based materials have been widely documented by theoretical calculations through vacancy defects to highlight the adsorption improvement of 3d-TM clusters for hydrogen storage applications. [45][46][47][48][49][50][51][52][53][54][55][56][57][58][59] Moreover, the vacancy defects play an important role in the adsorption stability of 3d-TM clusters. Recently, the small Co clusters have been explored for hydrogen storage, both free-standing clusters and supported on graphene. [56,[60][61][62][63] The combination of different sizes of Co clusters on the graphene surface favors the hydrogen storage capacity. [56,60] Several theoretical works have explored the capacity of the small free Co clusters to adsorb H 2 molecules. [61][62][63][64] Nevertheless, the interaction of small TM clusters (Co n ) on the -GY structure with vacancy defects has not been fully understood. Motivated by the recent experimental evidence reported by Cui et al., [22,23] the novel synthesis of the 2D carbon monolayer structures derived from the pristine -GY, such as -C≡C-vacancy (GY-def) and the nitrogendoped (GYN-def) systems have already been performed. However, the ability to adsorb TM atoms or small TM clusters on these monolayers has not yet been explored.
In this work, DFT calculations were perform to understand the role of the local defects at the 2D carbon monolayers on the adsorption of small cobalt clusters (Co n , n= 1-5). A systematic study of the hydrogen storage is also performed. The adsorption energy values, charge difference density (Δ ), and density of states (DOS) were analyzed. The capacity to retain hydrogen was also evaluated by Born-Oppenheimer molecular dynamics simulations. The theoretical results may aid to tailor novel materials with specific properties and intended to be applied as hydrogen storage systems.

Electronic Structure of GYN Monolayers
The atomic structure of -GY monolayer can be described by two unit cells (hexagonal cells). The molecular representation of the two optimized structures are depicted in Figure S1, Supporting Information. Both hexagonal structures are equivalent. Figure S1a, Supporting Information represents a six-membered ring geometry, while Figure S1b, Supporting Information shows a twelve-membered ring structure. Note that both contain 12 carbon atoms. The representation of the total density of states (DOS), partial DOS (PDOS), and electronic band structure along the Γ -K -M -Γ path in the hexagonal cell are depicted in Figure S1, Supporting Information. They maintain the same electronic behavior (band gap close to 0.47 eV at the Γ point). The Fermi level (E F ) is located at the half between the valence band maximum (VBM) and conduction band minimum (CBM). It is worth noting that the structural parameters did not show significant changes as shown in Table S1, Supporting Information. These results are consistent with previous theoretical works. [14][15][16]65] Therefore, the Perdew-Becke-Ernzerhof (PBE)-D2/DZP level of theory adequately describes the electronic and structural properties of the pristine -GY monolayer.
The description of the electronic structure of the pristine -GY monolayer from the fully optimized 2 × 2 × 1 supercell is depicted in Figure 1a. The characteristic structural parameters of pristine -GY are presented in Figure 1a. These values are similar to those obtained in previous works. [14][15][16]65] Then, we designed the -GY system with a vacancy defect by removing the central -C≡Cbonding atoms (GY-def), as it is shown in Figure 1b. Finally, we obtained an additional structure by doping with nitrogen in the benzene-rings (in the GY-def model) to obtain two pyridine-rings (GYN-def), as it is depicted in Figure 1c. The characteristic structural parameters are inserted in the molecular representation in Figure 1a-c. The vacancy generated in the GY-def and GYN-def structures (GYNs-def) slightly varies by 1.0 Å. The vacancy defects generated in these structures (GY-def and GYN-def) allow the retention of Co atoms in the center of geometry, as occurs in graphene-like monolayers. [42,66,67] The charge density contours in the xy plane were calculated for the -GY, GY-def, and GYN-def systems, as shown in Figure 1d-f, respectively. The vacancy defect does not cause structural distortion or local rearrangement in the carbon bonds. As a consequence, the 2D hierarchy is kept for the GY-def, and GYNs-def structures (GYNs), which could be optimal to be implemented as a support. The electronic structure properties of the GYN systems were compared to the pristine -GY structure. Figure S2, Supporting Information shows the PDOS and electronic band structure for all monolayers under study. The pristine -GY maintained the semiconductor character with a direct band gap of 0.40 eV at the Γ point. The same behavior is maintained in the GYN-def system, where the N orbitals contribution is located in deep states (close to −1.9 eV). However, the GY-def system showed the absence of a band gap, which suggests the rising of a metallic character (see Figure S2b,c, Supporting Information). Electron localized function (ELF) analysis was also used to construct map slices, and to identify the high density locations in the carbon bonding (see Figure 1g-i). Figure 1g-i show the cases of the GYN-def structures where the presence of an electronic density in the vacancy of the GY-def structure is weaker as compared to the GYN-def system. Therefore, the local electronic environment in all the proposed monolayers are completely different.

Adsorption of Co n Clusters on -GY, GY-def, and GYN-def Monolayers
The geometry of the Co n clusters generated by Da Silva et al. [68] were considered throughout this work. Figure 2a shows the energy stability of the Co n clusters, which is reported with the values of the cohesion energies (E coh ). The trend of the cohesion energy shows that as the size of the Co n clusters increases, the energy stability is favored. Therefore, the nature of the Co n clusters to agglomerate into larger nanostructures could be expected. The Co-Co bond length increases as the size of the Co n clusters increase (see Figure 2a). All relaxed Co n systems maintain the initial structural geometry according to Da Silva et al. [68] The lowest-energy configurations in the adsorption of Co n (n = 1-5) clusters onto the pristine -GY, GY-def, and GYN-def surfaces were studied by considering four different sites; namely, in the center of the 12-membered ring (H1), in the six-membered ring (H2), in the bridge of the -C≡C-bonds (B1), and at the bridge of the -C=C-bonds (B2). Such locations were evaluated at the top initial position, at 2.5 Å with respect to the carbon monolayers (see Figure 2b). Additionally, the Co n clusters present two possible adsorption configurations. That is, horizontal and vertical   positions. Therefore, these configurations were also considered in the geometry relaxation procedure.
The proposed configurations, Co n @ -GY and Co n @GYNs-def were fully optimized. The results are presented in Figures S3-S7, Supporting Information, which correspond to the Co n clusters (n = 1-5) supported on the monolayers under study. The relative energy values are shown for each case to identify the lowestenergy configurations. Figure 3 and 4 show the adsorption energy (E ads-Co/GYN ) values and charge density difference (Δ Co n /GYN ) for all configurations of lowest-energy. According to Figure 3, the cobalt clusters are adsorbed through an exothermic process. For the cluster adsorption onto GYN monolayers, E ads-Co/GYN increases linearly as the size of the Co n cluster increases. This behavior is absent in the pristine -GY system. Therefore, vacancy defects induce the adsorption of Co n clusters. Experimentally, a high E ads-Co/GYN value would suggest that the Co n clusters lack of diffusion on the carbon monolayers. Figure 4, depicts the isosurfaces of electronic density difference, which were obtained in accordance to Equation 5. For the lowest-energy system of the form Co n @ -GY, the Co atoms are preferably adsorbed on the hollow site of the 12membered ring. The Co n @GY-def and Co n @GYN-def systems revealed a similar trend, excluding the interaction with the Co 3 cluster. In this case, one of the Co atoms is selective to migrate to a 12-membered ring (Co 3 @GYN-def). The region around the Co n atoms is subjected to electron depletion for most of the cases, except the Co 4 @ -GY and Co 4 @GYN systems, due to the presence of the Co atoms (see Figure 4). Additionally, the Hirshfeld charge (Q) total values on the Co n clusters were also reported. Charge transfer increases when the size of the Co n clusters increase. The largest amount of charge transferred (Q = 0.73 e) is observed for the Co 5 @GY-def system. It could be attributed to the strong adsorption energy of 19.68 eV. In the lowest-energy configurations, the Co atoms prefer to interact in the red regions, where the electronic density is high (red color) and located in accordance to the sliced ELF map (see Figure 1g-i). As a consequence, the vacancy defect plays an important role in the adsorption of the Co n clusters.
The influence of the Co n clusters on the carbon-monolayer was assessed with the PDOS in Figure S8, Supporting Information. The presence of Co atoms contributes with electronic states above at the Fermi level (E F ). This indicates that one to five Co atoms are the responsible to tune the systems into a metallic electronic behavior. The high contributions of the Co-PDOS are located from −1.0 eV up to the Fermi level for the Co 1−4 systems. While for the Co 5 interactions, the Co-PDOS high contributions are located from 2.0 eV to the Fermi level. In these energy ranges, C-PDOS are present, which may indicate an overlap of these electronic states, associated with probable chemical interactions among the Co atoms and the carbon monolayers (see Figure S8, Supporting Information). This behavior can be confirmed with the Δ Co n /GYN isosurfaces (see Figure 4).

H 2 Adsorption on Isolated Co n Clusters
The capacity of hydrogen storage by the Co n clusters was evaluated with the adsorption of xH 2 molecules (x = 1-6) on free Co n clusters. That is, one to six H 2 molecules were added to interact with each of the Co n clusters. It is worth to denote that these configurations were fully optimized. The local minima are presented in Figures S9-S13, Supporting Information. The initial position of the H 2 molecules was located at a distance about 2.50 Å from each of the free Co n clusters. For all Co n clusters, adsorption and dissociation of the H 2 molecules are evident. As depicted in Figure 5, the number of adsorbed H 2 molecules is directly proportional to the adsorption energy per molecule. As a consequence, the free Co n clusters have the capacity to chemisorb xH 2 molecules. This behavior is in agreement with DFT calculations previously performed, in which the ability of Co n clusters to adsorb and dissociate H 2 has been predicted. [60][61][62][63] The adsorption energy (E ads-H 2 ∕Co n ) values are relatively similar (with an average of 1.94 eV per H 2 ) when the Co 2-5 clusters interact with up to 6H 2 molecules. The adsorption and dissociation of the H 2 molecules are present around the free Co n clusters (see Fig-ures S9-S13, Supporting Information). This behavior may predict the tendency to use a smaller number of Co atoms to obtain comparable adsorption energies.

Interaction of xH 2 Molecules (x = 1-5) with Co n @ -GY, Co n @GY-def, and Co n @GYN-def Systems
The effect of the carbon monolayers ( -GY, GY-def, and GYN-def systems) on the ability to adsorb hydrogen molecules was studied by evaluating the interaction of the H 2 molecules (x = 1-5) with Co n @ -GY, Co n @GY-def, and Co n @GYN-def systems. The hydrogen saturation was considered as hydrogen storage. Several configurations were proposed. The initial configuration corresponds to a H 2 molecule located at a distance about 2.50 Å from each of the Co atoms exposed on the surface. All model systems were fully optimized for a number of xH 2 molecules (x = 1-5) at the DFT/PBE-D2/DZP level of theory.
The adsorption energies (given in eV per H 2 ) for each number of H 2 molecules are depicted in Figure 6. The adsorption energies (E ads-H 2 ∕Co n @GYN ) of the xH 2 molecules on the Co n @--GY systems showed an increasing behavior when H 2 , 2H 2 , and 4H 2 molecules were adsorbed, excluding the Co 4 @ -GY system (see Figure 6a). The magnitude of the interaction with the 3H 2 molecules increases as the size of the Co 3 @ -GY systems Figure 4. Isosurface of the charge density difference for the lowest-energy systems: a-c) Co 1 @GYNs, d-f) Co 2 @GYNs, g-i) Co 3 @GYNs, j-l) Co 4 @GYNs, m-o) and Co 5 @GYNs. The values under each cobalt structure stand for the Hirshfeld charges (Q). The yellow and blue isosurfaces represent the accumulation and depletion of electron density, respectively. The isosurface level is set to be 0.0025 e Å −3 .
increase. Particularly, the xH 2 -Co n @GY-def systems exhibited adsorption energy values smaller than 0.8 eV per H 2 . This behavior in adsorption energy can be attributed to the position of the Co atoms exposed on a smaller surface area. In Figure 6b, a repulsive behavior was evidenced when the 3H 2 , and 4H 2 molecules interacted on the Co 2 @GY-def and Co 3 @GY-def systems. This could be due to a decrease in E ads-H 2 ∕Co n @GYN values. Small adsorption energy values ranging from 0.30 to 0.65 eV per H 2 for the Co 4 @GY-def and Co 5 @GY-def systems were evaluated. Figure 6c shows the adsorption energy behavior of the xH 2 molecules adsorbed on the Co n @GYN-def systems. A smaller range from 0.30 to 0.60 eV per H 2 was obtained when 3H 2 , 4H 2 , and 5H 2 molecules interact with the Co 4 @GYN-def and Co 5 @GYN-def systems. Based on the adsorption energy reported in Figure 6, the carbon monolayers act as ideal supports that effectively improve the hydrogen storage capacity of the Co n clusters. The molecular structure in each inset corresponds to the systems that effectively retained the H 2 molecules. Figure 7 shows the 5H 2 -Co 5 @ -GY, 5H 2 -Co 5 @GY-def, and 5H 2 -Co 5 @GYN-def systems from top and side views. The adsorption energy values www.advancedsciencenews.com www.advtheorysimul.com are within the ideal range to obtain a reversible hydrogen storage (0.2 to 0.70 eV per H 2 ). [10,42,69,70] For all cases, dissociation and adsorption of the H 2 molecule are present on the deposited Co 5 cluster. Particularly, the 5H 2 -Co 5 @GY-def system shows two H 2 molecules interacting close to the Co 5 cluster. To corroborate the interaction of the two H 2 molecules, the electronic density is shown in Figure S14, Supporting Information. According the computed isosurface of the 5H 2 -Co 5 @GY-def system, the electronic density of the H 2 molecules slightly interacts with the supported Co 5 cluster. As a consequence, the resulting interaction distance (H-Co ≈ 2.77 Å) could be associated with an electrostatic type interaction. The difference of adsorption energy between the PBE-D2 method and VDW functional corroborated this assumption (the difference amounts to 0.07 eV).

Molecular Dynamics Study
The structural and thermal stability of the systems under study were also investigated by using Born-Oppenheimer molecular dynamics (BOMD) calculations with the same periodic boundary conditions. The selected systems were those that showed largest adsorption energies, such as 5H 2 -Co 5 @ -GY, 5H 2 -Co 5 @GY-def, and 5H 2 -Co 5 @GYN-def systems at a temperature of 1000 K with a simulation time >1000 fs (1 ps). Figure 8 shows the evolution of these systems with the profile of the molecular dynamics simulation over time (BOMD simulations starting from the relaxed geometries). The 5H 2 -Co 5 @ -GY system exhibit variations of the potential energy in an interval of ≈3 eV. After 1 ps, the hydrogen atoms remain in the Co 5 cluster. The 5H 2 -Co 5 @GY-def system presents a potential energy variation of ≈8 eV, where an instability in the Co-H bonds is observed. Consequently, at the time step of 1 ps, the desorption of six-H 2 molecules could be identified. The 5H 2 -Co 5 @GYN-def system exhibits variations at an energy interval of ≈6 eV, and a bond elongation in the GYN-def monolayer. After 1 ps of time step, the Co 5 cluster maintained the adsorption with eight hydrogen atoms. Therefore, after completing the BOMD simulations, the behavior of the systems under study shows 100%, 60%, and 80% of hydrogen retention for Co 5 @ -GY, Co 5 @GY-def, and Co 5 @GYN-def systems, respectively. Mul- Figure 6. Adsorption energies (E ads-H 2 ∕Co n @GYN ) of xH 2 molecules (x = 1-5) for systems: a) Co n @ -GY, b) Co n @GY-def, c) Co n @GYN-def. The molecular structure inset in each plot corresponds to the interaction enclosed in a black circle. All numerical values are reported in Tables S2-S4,  Supporting Information. timedia material (video.mov) was also included in the Supporting Information, corresponding to the BOMD simulations. Therefore, the high retention capacity of H 2 molecules of the Co 5 @ -GY system after the BOMD simulation may be associated with the adsorption energy (0.74 eV per H 2 ). While the difference in retention capacity for Co 5 @GY-def (60%) and Co 5 @GYN-def (80%) systems may be associated to the geometrical disposition in which the Co atoms are exposed on the surface area (see Figure 7b,c). Likewise, a structural stability of the Co 5 cluster on the surfaces is observed. This is strongly correlated with the high adsorption energy values according to Figure 3. This feature may aid to reduce the instability issue in small Co n clusters [71,72] by depositing Co atoms in modified carbon monolayers.

Conclusions
By using DFT-D2 calculations, the capacity to adsorb Co n clusters (n = 1-5) on pristine -GY, GY-def, and GYN-def surfaces was studied. Electronic structure properties and charge transfer in the systems were analyzed. According to the E ads values, the Co n clusters are selective to the adsorption on surfaces with vacancy defects. Large E ads values may be associated to the large charge transfer by the Co atoms toward the carbon surfaces, according to the isosurface of charge density difference and Hirshfeld charge values. The structural stability of the supported Co 5 clusters is evident by keeping its trigonal bipyramidal geometry. The behavior of isolated and supported Co n clusters was evaluated by adsorption of xH 2 molecules (x = 1-6). The analysis of adsorption energies per H 2 molecule showed that the Co n @ -GY, Co n @GY-def, and Co n @GYN-def systems exhibited an improvement in the ideal energy range for efficient reversible hydrogen storage. Additionally, BOMD calculations were performed on the structurally stable systems 5H 2 -Co 5 @ -GY, 5H 2 -Co 5 @GYdef, and 5H 2 -Co 5 @GYN-def. The temperature effect plays an important role to reduce the amount of interacting H 2 molecules in our systems. These systems are capable to store hydrogen with 60% of efficiency. Furthermore, we contributed to understand the molecular interaction among H 2 and Co n @GYN-def systems. These results could be implemented as a tool to guide the design of novel nanomaterials for hydrogen storage systems.

Experimental Section
All calculations in this work were performed by using DFT as implemented in the SIESTA computational code. [73] In order to describe the exchange-correlation term, the generalized gradient approximation was applied with the exchange-correlation functional proposed by PBE. [74,75] The long-range van der Waals forces were described by the semiempirical correction in the Grimme's scheme (D2). [76] The Troullier-Martins [77] norm-conserving pseudopotentials were used and a double-plus basis set was considered, including a polarization function (DZP) with an energy shift of 0.2 eV for all calculations. The real-space mesh cut-off energy was set to 250 Ry (3400 eV). The convergence threshold of the density matrix, the total energy and Hellmann-Feynman forces were set to 10 −4 eV, 10 −5 eV, and 0.04 eV Å −1 , respectively. A Monkhorst-Pack grid [78] of 3 × 3 × 1 k-points was used for the Brillouin zone sampling. The DOS were calculated with 6 × 6 × 1 k-points. The optimized structures were modeled using the 2 × 2 × 1 supercell (48 carbon atoms) of the pristine -GY. A 20 Å of vacuum was left in the z-direction to separate the neighboring slabs, and avoid spurious interactions. The cohesive energy (E coh ) of isolated Co n clusters were obtained according to the following equation where E Co n corresponds to the total energy of the free cobalt cluster. The term E isolated-Co is the free energy of a single Co-atom, and n is the total number of atoms at the cluster. To evaluate the stability of the adsorbed Co n clusters (n = 1-5) onto the -GY, GY-def, and GYN-def monolayers (GYNs), several possible configurations with different relative positions of the M atoms were proposed. For the search of the lowest-energy configuration, the adsorption energy (E ads-Co/GYN ) was calculated as in which E t (Co n @GYN monolayer ) is the total energy of the proposed GYN monolayers, interacting with the Co atoms. E t (GYN monolayer ), and E t (Co n ) stand for the energy of the given GYN monolayers, and the isolated Co cluster adsorbates, respectively. In Equation (2), positive E ads-Co/GYN values indicate a stable exothermic adsorption. The interaction of the H 2 molecules with the bare Co n clusters was also assessed (E ads-H 2 /Co n ), in accordance to the following adsorption energy relationship In Equation (3), the term E t (Co n + xH 2 ) represents the total energy of the H 2 molecules that were grafted on the bare Co n clusters. Additionally, the terms E t (Co n ) and E t (H 2 ) are the total energy of the bare Co n clusters and the isolated H 2 molecule. The x indicates the number of H 2 molecules that were grafted to the Co n clusters.
in which x is the number of hydrogen molecules, E t (Co n @GYN monolayer + xH 2 ), is the optimized total energy, E t (Co n @GYN monolayer ) and E t (xH 2 ), corresponded to the Co n @GYN configuration energy, and the hydrogen molecule energy with no interaction, respectively. Since it is expected that the adsorption of the Co n clusters on the GYN monolayers, and the adsorption of the H 2 molecules on the composite Co n @GYN-monolayers, represent an exothermic reaction, a minus sign was imposed for Equations (2)(3)(4). This is considered throughout the discussion of the present work. The charge density difference (Δ Co n ∕GYN ) corresponding to the electron redistribution caused by the presence of the Co n clusters on the GYN monolayers, is given by the following equation Δ Co n ∕GYN = Co n @GYN monolayer − GYN monolayer − Co n In Equation (5), Co n @GYN monolayer , GY monolayer , and Co n corresponded to the electronic density of the given carbon-monolayers interacting with the Co n clusters, the isolated carbon-monolayer, and the isolated Co n clusters, respectively.
The temperature effect was also assessed, since it represents a relevant factor in hydrogen storage devices. [79] BOMD simulations were performed using the Quantum ESPRESSO (QE) code. [80] BOMD calculations were carried out from the relaxed structures, with an equilibrium time-step of 1.69 fs. Moreover, 600 steps were considered, and a constant temperature of 1000 K (NVT Nosé-Hoover thermostat) during 1000 fs (1 ps). This temperature was chosen to understand the structural stability of the supported Co clusters.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.