Introduction

When H atoms are adsorbed on Cu(100), the coverages of these atoms are determined by the chemical potential of hydrogens. In order to gain better understanding of such behavior, we will carry out a thermodynamic ab initio study in this project. We will firstly start from the calculations of the energy of bulk Cu as well as gas-phase hydrogen atom. The energy of the adsorption system with different coverages of H atoms will then be calculated. Convergence tests with respect to energy cutoff and k-point sampling will be performed for each calculation.

The Vienna Ab initio Simulation Package (VASP) is used to perform all the periodic density functional theory (DFT) calculations,^1-3employing the projected augmented-wave (PAW) pseudopotentialsand the generalized gradient approximation with exchange-correlation functional by Perdew, Burke, and Ernzerhof (PBE).^4-6The conjugated gradient algorithm is chosen to perform the geometric optimization in VASP. For all the calculations in this project, the dipole correction along z-axis is included and the dispersion correction method of Grimme is introduced to better describe the interaction between H atoms and Cu surfaces.⁷

 

Bulk Cu and Gas-Phase Hydrogen

The unit cell with one Cu atom in it is constructed to calculate the energy of bulk Cu (the energy of one Cu atom in bulk phase) as shown below in Figure 1. The lattice vectors are (0.5, 0.5, 0.0), (0.0, 0.5, 0.5), and (0.5, 0.0, 0.5). The lattice constant is chosen to be the experimental value of 3.615 Å.

Figure 1: Unit cell of bulk Cu

 

According to the previous calculation on fcc bulk metal using the same unit cell, the k-point sampling of 15x15x15 will give the accuracy within 0.01 eV and here we use the sampling of 21x21x21 to test the energy cutoff. The results for the convergence tests are shown below in Table 1.

Energy Cutoff (eV)	400	450	500	550	600	650
Bulk Energy (eV)	-4.1428	-4.1410	-4.1414	-4.1417	-4.1416	-4.1417

Table 1. Convergence of energy cutoff

 

Here we can see that, the energy cutoff of 500.00 eV can give us the accuracy of 0.0003 eV and this energy cutoff will be used for the following calculations (the energy cutoff for H should be much lower than Cu, so the 500.00 eV energy cutoff is appropriate for H).

The energy of gas-phase hydrogen atom is calculated using a cubic box with side length of 25.0 Å and the energy is -6.77 eV and the energy for one H atom is thus -3.385 eV.

 

Adsorption of H atoms on Cu(100)

In order to calculate the adsorption of H on Cu(100), a six layer slab model is used to model the Cu surface. The bottom three layers are fixed at bulk position with lattice constant of 3.615 Å and the upper three layers along with the H atoms are fully relaxed during the geometry optimization. For the 1/9 ML coverage, a supercell of 3×3 is used with one H atom on it, while the 2×2 supercell is used for the 0.25 ML and 0.50 ML coverages with one and two H atoms on it, respectively. The 1×1 unit cell is used for the 1ML coverage.

Before we carry out the calculations of these adsorption system, we have to firstly determine the adsorption site of H atoms on Cu(100). Here We choose the Cu(100)-3×3 supercell with three different initial adsorption sites, including atop, bridge, and hollow, to study the energy-favored adsorption site. The energy cutoff of 500.00 eV and sampling of 6x6x1 are used here (whose convergence has not been teste but reliable for comparisons between several rough geometry optimizations). The optimized surfaces are shown below in Figure 2.

Figure 2: Binding conformation of H atoms at different sites, from left to right: atop, bridge and hollow, only the upmost layer is shown here.

 

The system energies are -213.132, -213.603, and -213.736 eV for the atop, bridge, and hollow cases, respectively. So we can see, the hollow site is most energy favored and will be chosen as the initial sites for the later calculations.

The convergence tests with respect to the k-point sampling on these three supercells (1×1, 2×2, and 3×3) are summarized below in Table 2.

1×2	17x17x1	19x19x1	21x21x1	23x23x23
Energy (eV)	-26.861	-26.864	-26.864	-26.866
2×2	7x7x1	8x8x1	9x9x1	10x10x1	11x11x1	13x13x1
Energy (eV)	-96.941	-96.998	-96.963	-96.987	-96.981	-96.985
3×3	5x5x1	6x6x1	7x7x1	8x8x1	9x9x1	11x11x1
Energy (eV)	-213.760	-213.736	-213.729	-213.711	-213.710	-213.711

Table2: Summary of convergence tests on k-point samplings

 

The k-point sampling for these three supercells are chosen to be 19x19x1, 11x11x1, and 9x9x1 to provide an accuracy of 0.006 eV. The system energies used for the later ab initio calculations are -26.864, -96.981, and -213.710, respectively.

 

Thermodynamic ab initio calculations

In order to find the equilibrium points (with respect to the chemical potential of H atom) between different coverages (In this project, 0.00 ML, 0.11 ML, 0.25 ML, 0.50 ML and 1.00 ML are considered), we define the normalized grand potential for each coverage as:

\begin{equation} \Omega_{i}(\mu_{H}) = \frac{E_{i}(N_{Cu, i}, N_{H, i})-\mu_{H}\times N_{H, i}-\mu_{Cu}\times N_{Cu,i}}{N_{Cu,i}} \end{equation}

Thus, we can have the grand potential diagram to show the potential-dependent coverages of H atoms on Cu(100). The grand potential diagram are shown below in Figure 3. The x-axis is the chemical potential difference while the DFT energy for H atom is used as the reference. The range of x-axis is [-2.7, -2.2], as the bare Cu(100) will be preferred at lower chemical potential out of this range and the 1.00 coverage case will be preferred at higher chemical potential out of this range. The four equilibrium points are -3.6074, -3.5983, -3.5301, and -3.4721 eV.

Figure 3: Grand potential vs chemical potential of H atoms

 

To have a clearer understand of such grand potential diagram, we further reduce the x-range to [-2.6, -2.3] and only keep lines which have the minimal grand potential at each different chemical potential of H atoms as shown below in Figure 4, which allows us to learn the potential-dependent coverages of H atoms on Cu(100).

Figure 4: Minimal grand potential at different chemical potential of H atoms

To get the phase diagram, we calculate the  according to the equation (2) as listed below, while in that equation is obtained from equation (3). All the data are acquired from NIST.⁸

\begin{equation}\ln{p/p^o}=\frac{1}{k_B T} (\mu_{H_2}-E_{H_2}^{DFT}-\tilde{\mu}_{H_2}(T, p^o))\end{equation}

\begin{equation}\tilde{\mu}_{H_2}(T, p^o)=[H(T)-H^o(T_r)]-TS(T)-[H(0)-H^o(T_r)]\end{equation}

The phase diagram is shown below in Figure 5.

 

Reference

1. Kresse, G. and J. Hafner,Ab initio molecular dynamics for liquid metals.Physical Review B, 1993. 47(1): p. 558.
2. Kresse, G. and J. Hafner,Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium.Physical Review B, 1994. 49(20): p. 14251.
3. Kresse, G. and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.Physical review B, 1996. 54(16): p. 11169.
4. Kresse, G. and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method.Physical Review B, 1999. 59(3): p. 1758.
5. Blöchl, P.E., Projector augmented-wave method.Physical review B, 1994. 50(24): p. 17953.
6. Perdew, J.P., K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple.Physical review letters, 1996. 77(18): p. 3865.
7. Grimme, S., Semiempirical GGA-type density functional constructed with a long-range dispersion correction. Journal of computational chemistry, 2006. 27(15): p. 1787-1799.
8. Stull, Daniel Richard, and Harold Prophet. JANAF thermochemical tables. No. NSRDS-NBS-37. National Standard Reference Data System, 1971.