Abstract
We investigated the vibration modes of hydrogen atom when it is absorbed on a copper (111) surface at fcc site. The effect of zero-point energy is also discussed at the end of this paper.
Copper (111) Surface with add-on Hydrogen
Fig.1 Hydrogen absorbed at hcp site
Fig. 2 Hydrogen absorbed at fcc site
We build a slab containing 3 top layers of copper atoms on (111) surface. The hydrogen atom is added on hcp and fcc sites separately, as shown in figure 1 and figure 2. Then the geometry optimization is executed, during which the position of the hydrogen is fully adjusted while only the vertical positions of the top copper atoms are adjusted. The hcp and fcc sites are two local energy minima, we then apply harmonic approximation to find out the vibration modes near those two points. The energy is -5057.7236 eV (5057.7283 eV) for hcp (fcc) in our calculation.
Hessian Matrix and Vibrational modes
Close to the local energy minima, the energy can be approximated by harmonic oscillators. The vibration modes can hence be represented as the eigenvectors of the mass-weighted energy matrix, called Hessian matrix [2]:
\begin{equation}H_{ij}=\frac{\partial^2 E}{\partial x_i \partial x_j}\end{equation}
and
\begin{equation}A_{ij}=H_{ij}/m_i=\frac{1}{m_i} \frac{\partial^2 E}{\partial x_i \partial x_j}\end{equation}
The vibration frequencies can be obtained by diagonalizing the Hessian matrix:
\begin{equation} \nu_i=\sqrt{\lambda_i}/(2 \pi) \end{equation},
where 
We apply the following numerical derivative in our calculation, which has 2nd order accuracy [1],
\begin{equation}\frac{\partial^2 E(x)}{\partial x_i^2}=\frac{E(x+\delta)+(E(x-\delta)-2E(x))}{\delta^2},\end{equation}
and
\begin{equation}\frac{\partial^2 E(x,y)}{\partial x\partial y}=\frac{E(x+\delta,y+\delta)+E(x-\delta,y-\delta)-E(x+\delta,y-\delta)-E(x-\delta,y+\delta)}{4\delta^2}\end{equation}
In this report, due to the significant mass difference between copper and hydrogen, only the vibration of hydrogen is considered. We therefore calculated the energy of the system with hydrogen atom at certain displacements from its equilibrium. The energy data are listed in table 1 and table 2.
Table 1 Energy at different displacements for hcp site
Table 2 Energy at different displacements for fcc site
The frequencies for hcp and fcc sites based on this method are:
\begin{equation}\nu_1=8.738 \times 10^{13} Hz, \nu_2=8.206\times 10^{13} Hz, \nu_3=8.181\times 10^{13} Hz\end{equation}
for hcp, and
\begin{equation}\nu_1=8.964\times 10^{13} Hz, \nu_2=8.584\times 10^{13} Hz, \nu_3=8.511\times 10^{13} Hz\end{equation}
for fcc, respectively.
Zero-point Energy
According to quantum mechanics, each vibration mode contributes zero-point energy 
\begin{equation} E_0=\frac{1}{2} h (\nu_1+\nu_2+\nu_3)=0.5196 eV\end{equation}
and
\begin{equation}E_0’=\frac{1}{2} h (\nu_1+\nu_2+\nu_3)=0.5389 eV\end{equation}
Considering this energy correction, the ground state energy for two sites are
\begin{equation}E_{hcp}=-5057.2040 eV\end{equation}
and
\begin{equation}E_{fcc}=-5057.1894 eV\end{equation}.
The astonishing result of taking the zero-point energy into account is that the hcp has a lower energy, which is opposite to the classical result, showing that the former is more energetically favorable.
Summary
In this paper, we report a systematical method to evaluate the ground state energy involving the vibrations of atoms. Our result shows the significance of zero-point energy, especially for light atoms. However, due to the time limit, we did not investigate the effect of neighbouring hydrogen atoms (i.e. the hydrogen is also influenced by the hydrogen absorbed next it).
Acknowledgment
We thank Prof. Sofo for his advice on doing DFT calculations and Run for his discussion on the calculation. We also express our gratitude for software support from Penn State University.
Appendix: Parameters for DFT Calculation
The parameters used in the geometry optimization are listed here:
Lattice parameters: a=b=2.58Å, c=13.91Å; α=β=90°, γ=120°
Atomic and pseudo atomic structure for H: 1s1
Atomic structure for Cu: 1s2 2s2 2p6 3s2 3p6 3d10 4s1
pseudo atomic structure for Cu: 3d10 4s1
Functional: GGA PBE
Dipole correction: Self-consistent
Cut off energy: 408.2eV
k-point set: 11*11*2
Total energy / atom convergence: 1.0E-004eV
Max ionic |force| tolerance: 0.03eV/Å
Max |stress component| tolerance: 0.05GPa
Max ionic |displacement| tolerance: 0.001Å
The top layer of copper are allowed to relax on z-direction, the positions of hydrogen atoms are fully optimized. The bottom layers and the super lattice parameters are fixed during the relaxiation.
For energy calculation, we use the same parameters. The spatial step value is chosen such the energy difference is about 10~100 times bigger than the energy convergence.
Reference
[1] http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h07/undervisningsmateriale/kap7.pdf
[2]Sholl, David, and Janice A. Steckel. Density functional theory: a practical introduction. John Wiley & Sons, 2011.