A Novel First-Principles Molecular Dynamics Approach for Testing the Attraction Forces in Hydrogen Storage on a 2D Material at Different Temperatures

1. Introduction

Scientific interest in finding alternatives to the worldwide problems stemming from population growth and limited fossil fuel resources has increased in recent years. The focus has been on producing green energy and using materials to store it. Hydrogen storage is a desirable option for this purpose.

The development of materials for hydrogen storage remains a significant challenge. The US Department of Energy (US-DOE) highlights the importance of finding a material that can store hydrogen with high gravimetric and volumetric density while being environmentally friendly [1]. To achieve optimal adsorption and desorption, the adsorption energies for storing H₂ molecules should fall within the range of 0.1 eV to 0.2 eV [2,3]. Various methods exist to generate renewable energy sources for hydrogen production. Electrolytic hydrogen production offers an alternative to hydrogen production for internal combustion engines and fuel cells [4]. The hydrogen evolution reaction (HER) is a critical electrocatalytic reaction, with hydrogen adsorption on the electrode playing a crucial role. In the past, Pt-based materials were considered the best HER catalysts. However, it is essential to find cheaper and more abundant alternatives to Pt materials [5]. Molybdenum sulfide was the first reported alternative to Pt-based materials for catalyzing the HER in acidic aqueous solutions [6]. Since then, Mo-based materials, including molybdenum phosphides and transition metal phosphides (TMPs), have gained prominence as low-cost HER electrocatalysis materials [7,8,9].

We investigated the direct interaction between a MoP₂ surface and H2 molecules for hydrogen storage using DFT calculations. Additionally, we explored the influence of Mo vacancies on the gravimetric capacity of hydrogen adsorption on a MoP₂ surface. Our study involved FPMD calculations at 300 K, 77 K, 4 K, and 3 K to assess the system's stability under varying conditions.

2. Materials and Methods

We performed ab initio DFT calculations to study hydrogen storage on the 2D material MoP₂. We considered the pristine material and the same material with 12.5% Mo vacancies. We employed the generalized gradient approximation (GGA) with the exchange-correlation energy given by Perdew-Burke-Ernzerhof double-zeta polarized basis sets and norm-conserving pseudopotentials [10,11]. Our calculations used the SIESTA code [12] and the spin-polarized density functional formalism. We followed the Monkhorst–Pack scheme [13] and considered a Brillouin zone sampling of 24x24x1 k-points. We took an energy cutoff of 180 Ry for numerical integrations. The geometries of the studied systems converged until forces between atoms were smaller than 0.01 eV/Å. To consider dispersion forces, we included the Van der Waals interactions using Grimme's semiempirical approach [14].

We carried out our simulations of these systems with a supercell 2x2, with a separation between layers of 30 Å. In the case of the pristine 2DMoP₂ system, the supercell consists of four molybdenum (Mo) atoms and eight phosphorus (P) atoms (Figure 1a). In comparison, the 2D MoP₂ with a Mo vacancy has three molybdenum (Mo) particles and eight phosphorus (P) atoms (Figure 1b). In this manner, the Mo vacancy in the cell of twelve atoms corresponds to 12.5% of vacancies. In the system's initial configuration (pristine or with 12.5% of vacancies), we added hydrogen molecules one by one above a Mo or P atom to optimize the system, and the original distance of the hydrogen molecules from the surface was 3 Å.

Every structure was relaxed to determine the most stable configuration before the saturation of H₂ molecules.

We calculated the adsorption energy using equation (1).

$${\mathrm{E}}_{\mathrm{ads}}={\mathrm{E}}_{\mathrm{final}\mathrm{configuration}}–\left(n{\mathrm{E}}_{\mathrm{H}}+{\mathrm{E}}_{\mathrm{surface}}\right)$$

(1)

Where E_ads is the adsorption energy, n is the number of hydrogen molecules we considered, E_H is the total energy of one free H₂ molecule, and E_surface is the total energy of the surface unit cell. We obtained by optimization an average adsorption energy of ‒ 0.13 eV per molecule.

We calculated the gravimetric capacity Wt% using equation (2).

$$\mathrm{Wt}\%=100\mathrm{n}{\mathrm{W}}_{\mathrm{H}}/\left(n{\mathrm{W}\mathrm{E}}_{\mathrm{H}}+{\mathrm{W}}_{\mathrm{surface}}\right)$$

(2)

W_H is the molecular hydrogen mass, n is the same as before, and W_suface is the mass in the unit cell surface.

We established a new criterion to test the attraction forces in hydrogen storage at a given temperature. Thus, we took a first-principles molecular dynamics approach at different temperatures. So, if we want to know how many molecules remained attached to the surface at temperature T₁, we proceed as follows. We performed an FPMD calculation at temperature T₁, taking the optimization as the initial system configuration. Then, we perform a second FPMD calculation at a temperature T₂ ( with T₂ ‹‹ T₁), bringing the stable configuration of the first FPMD calculation as the initial configuration. We identified as adsorbed H₂ molecules at temperature T₁, only those that moved back toward the surface at temperature T₂.

With the optimization of the system and the adsorbed hydrogen molecules as the starting point, we used the Quantum Espresso Code [15] to perform the first-principles molecular dynamics calculations at several temperatures. We employed XCrysden [16] to visualize the system’s configurations.

3. Results

3.1. Hydrogen Storage on Pristine 2D MoP₂

3.1.1. Optimization

The pristine monolayer can adsorb up to 16 H₂ molecules—eight molecules per side, as shown in Fig. 2, which presents the final configurations for one side. We calculated the adsorption energy using equation (1).

$${\mathrm{E}}_{\mathrm{ads}}={\mathrm{E}}_{\mathrm{final}\mathrm{configuration}}–\left(n{\mathrm{E}}_{\mathrm{H}}+{\mathrm{E}}_{\mathrm{surface}}\right)$$

(1)

Where E_ads is the adsorption energy, n is the number of hydrogen molecules we considered, E_H is the total energy of one free H₂ molecule, and E_surface is the total energy of the surface unit cell. We obtained by optimization an average adsorption energy of ‒ 0.13 eV per molecule.

We calculated the gravimetric capacity Wt% using equation (2).

$$\mathrm{Wt}\%=100\mathrm{n}{\mathrm{W}}_{\mathrm{H}}/\left(n{\mathrm{W}\mathrm{E}}_{\mathrm{H}}+{\mathrm{W}}_{\mathrm{surface}}\right)$$

(2)

W_H is the molecular hydrogen mass, n is the same as before, and W_suface is the mass in the unit cell surface.

We obtained a gravimetric capacity for hydrogen storage of 5.72%, and the average adsorption energy is ‒ 0.13 eV/molecule.

3.1.2. First-Principles Molecular Dynamics Calculations

We followed the criterion described in the Materials and Methods section to calculate the gravimetric capacity, determining how many molecules adsorbed the surface at a given temperature. Using the optimization of the system and the adsorbed hydrogen molecules as the starting point, we performed the first-principles molecular dynamics calculations at 300K.

We found that at 300 K, all the adsorbed H₂ molecules in the initial optimization moved away from the surface, as Figure 3 shows.

We then reduced the temperature to 3 K and performed a second FPMD calculation to determine how many remained bound to the surface. Thus, we considered the stable configuration obtained at 300 K as the initial state for the second FPMD calculation, as shown in Figure 4. We identified as adsorbed H₂ molecules only those that moved back to the surface. In this way, we obtained that at 300 K, the pristine surface retained only two H₂ molecules per side, which implies 1.28% of gravimetric storage capacity.

Again, at 77 K, we considered the optimization configuration the initial state and performed an FPMD calculation, as shown in Figure 5, and used the same criterion. Afterward, we reduced the temperature to 3 K and performed another FPMD calculation to determine how many remained bound to the surface (see Figure 6). Thus, we identified as adsorbed H2 molecules only those that moved back to the surface. There were only three H₂ molecules per side of the surface, with a gravimetric capacity for molecular hydrogen storage of 1.92%.

In the case of 4 K, Figure 7 shows the FPMD calculation where the optimization is the initial configuration. Notice that all the molecules remained close to the surface.

In the following FPMD calculation at 3 K, all the hydrogen molecules moved back to the surface, indicating that they adsorbed on the pristine 2D MoP₂ surface, as Figure 8 shows. The gravimetric capacity for molecular hydrogen storage is 5.1% for 16 adsorbed hydrogen molecules (eight per side).

In summary, in the optimization, the pristine 2D MoP₂ surface adsorbed up to 16 hydrogen molecules (8 per side), with a gravimetric molecular hydrogen storage capacity of 5.1% and an average adsorption energy of ‒ 0.13 eV/molecule. However, we performed an FPMD calculation at 300 K and applied our already mentioned criterion to determine how many molecules adsorbed the surface at a given temperature. At that temperature, we found that only two (per side) hydrogen molecules adsorbed on the pristine 2D MoP₂ surface, with a gravimetric capacity of 1.28%. When we considered a temperature of 77 K and applied the same criterion, we found a gravimetric capacity for molecular hydrogen storage of 1.92%. For 4 K, the gravimetric capacity increases to 5.1%.

3.2. Hydrogen Storage on 2DMoP₂ with Mo Vacancies

3.2.1. Optimization

To simulate the adsorption of H₂ molecules on the surface of 2DMoP₂ with vacancies, we considered a supercell 2x2x1, as Fig. 1 (b) shows. In this manner, we made a Mo vacancy in the cell of twelve atoms, corresponding to 12.5% of vacancies. Again, in the system's initial configuration, we added hydrogen molecules one by one above a Mo or P atom to optimize the system, and the original distance of the hydrogen molecules from the surface was 3 Å. We obtained that 2D MoP₂ with Mo vacancies can adsorb up to 16 H2 molecules (eight molecules per side), as Figure 9 shows, and a gravimetric capacity of 6.0%. We calculated average adsorption energy of around ‒ 0.14 eV/molecule.

3.2.2. First-Principles Molecular Dynamics Calculations

We followed the procedure in Section 2.1 to calculate the gravimetric capacity at 300 k. Figure 10 shows the FPMD calculation at that temperature. As expected, all the molecules displaced away from the surface.

Figure 11 shows the subsequent FPMD calculation at 3 K, starting with the stable configuration at 300 K. Notice that seven hydrogen molecules (per side) moved towards the surface when the temperature decreased. Thus, the gravimetric capacity for hydrogen storage on the surface is 5.27%.

In the case of 77 K, we again considered the optimization configuration as the initial state and performed an FPMD calculation at that temperature, as Figure 12 shows. Again, all the molecules moved away from the surface.

In Figure 13, we present the subsequent FPMD calculation at 3 K, starting with the stable configuration at 77 K. Notice that all the eight hydrogen molecules (per side) moved towards the surface when the temperature decreased, showing a gravimetric capacity of 6.02%.

Figure 14 shows the FPMD calculation at 4 K, taking the optimization as the initial configuration. Again, we can notice that all the hydrogen molecules moved away from the surface. The elapsed time was two picoseconds.

In Figure 15, we present the subsequent FPMD calculation at 3 K, starting with the stable configuration at 4 K. Notice that all eight hydrogen molecules (per side) moved towards the surface when the temperature decreased, showing a gravimetric capacity of 6.02%.

In summary, when we considered the 2D MoP₂ with 12.5% Mo vacancies and optimized the system, the surface adsorbed up to 16 hydrogen molecules (8 per side). This corresponds to a gravimetric molecular hydrogen storage capacity of 6.02%, and the average adsorption energy is ‒ 0.14 eV/molecule. However, when we applied our already described approach to determine how many molecules adsorbed the surface at a given temperature and performed an FPMD calculation at 300 K, we found that only seven (per side) hydrogen molecules adsorbed on the surface. The gravimetric capacity is 5.27%. Using the same methodology at 77 K, we found a gravimetric capacity for molecular hydrogen storage of 6.02%. As we expected, for 4 K, the gravimetric capacity is 6.02%, too.

4. Discussion

We performed ab initio DFT calculations to study hydrogen storage on the 2D material MoP₂. We considered the pristine material and with Mo vacancies.

We found by optimization of the system that the pristine 2D MoP₂ surface adsorbed a maximum of 16 hydrogen molecules (8 per side), with an average adsorption energy of ‒ 0.13 eV/molecule. Furthermore, we considered the same surface with 12.5% Mo vacancies and obtained by optimization a gravimetric capacity of 6.02% and an average adsorption energy of ‒ 0.14 eV/molecule. We had to consider dispersion forces. Thus, we included the Van der Waals interactions using Grimme's semiempirical approach [14].

We established a new criterion to test the attraction forces in hydrogen storage at a given temperature described in Section 2. Thus, we took a first-principles molecular dynamics approach at different temperatures. So, if we want to know how many molecules remained attached to the surface at temperature T_1, we proceed as follows. We performed an FPMD calculation at temperature T₁, taking the optimization as the initial system configuration. Then, we perform a second FPMD calculation at a temperature T₂ ( with T₂ ‹‹ T₁), bringing the stable configuration of the first FPMD calculation as the initial configuration. We identified as adsorbed H₂ molecules at temperature T_1, only those that moved back toward the surface at temperature T₂.

We followed the new criterion described above to calculate the gravimetric storage capacity for molecular hydrogen storage and report the results described in Section 3 for 300K, 77 K, and 4 K.

Our results suggest that using only optimization may lead to unreliable predictions for hydrogen storage. Applying our proposed criterion, we conclude that the 2D MoP₂ surface with Mo vacancies was adequate for the same purpose. At 300 K, the defective 2D MoP₂ with 12.5% Mo vacancies had a 5.27% gravimetric capacity (the US-DOE target is 5.5% for automotive-grade molecular hydrogen storage). At 77 K and below this temperature, this system reached a gravimetric capacity of 6.02%, above the US-DOE target. Here, the disadvantage is maintaining the material at that temperature during the application for hydrogen storage in automotive-grade molecular hydrogen storage. However, our results indicate that further investigation of this class of 2D materials is needed to achieve successful molecular hydrogen storage at room temperature.