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.

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.

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.

In this project, the energies and lattice constants of the ferromagnetic and antiferromagnetic states of Co, MnF2, and MnBi were found and compared. To perform these calculations, a suitable energy cutoff and k-space sampling were selected to obtain a 0.01eV tolerance, and research was done in order to determine a proper starting c/a value of the lattice constant. From there, the the materials were placed into a ferromagnetic or antiferromagnetic state by aligning or misaligning their atomic spin polarizations and their energies were found at varying values of lattice constant a using CASTEP plane wave calculations [1].

These calculations were performed in Materials Studio using the GGA Perdew Burke Ernzerhof (PBE) functional and the pseudopotential was solved using the Koelling-Harmon atomic solver. In Materials Studio, spin polarizations are defined by the the total number of up electron spins minus down electron spins, so a spin polarization of +3, for example, implies that the atom has 3 more spin up than spin down electrons.

Cobalt

Fig. 1. Picture of cobalt cell used

For these calculations, cobalt was imported using the default Materials Studio structure, which gives a lattice constant of 2.507Å, a c/a value of 1.623, and forms the P63/MMC symmetry. From there, the spins in the 2 atom unit cell were polarized into a ferromagnetic +3, +3 state and energy calculations with varying energy cutoffs and k-samplings were performed until the results converged to a 0.01eV tolerance. The psuedopotential was calculating using the 1s2 2s2 2p6 3s2 3p6 electrons for the inner shells and 3d7 4s2 as the valence shells. Ultimately, a 420eV cutoff with a total k space sampling of 108 points were used.

Fig. 2. Energy cutoff calculations for Co

Fig. 3. k-space Sampling calculations for Co

From here, the lattice constant a was varied in 0.1Å and 0.005Å increments to find a suitable ferromagnetic lattice constant. These same calculations were repeated for the antiferromagnetic state, which were started with opposite +3 and -3 spin polarizations.

Fig. 4. Lattice constant calculations for ferromagnetic and antiferromagnetic states of Co

In the antiferromagnetic calculations, as the atoms approached one another, the calculations ceased to reach a local minimum in the antiferromagnetic state and began converging to non-spin polarized states. Additionally, some calculations failed to converge within 200 steps. Regardless, Co favors a ferromagnetic state, with a lattice constant of approximately 2.49Å, over the antiferromagnetic state, which appears to have a lattice constant between 2.465Å and 2.45Å (energy decreases past 2.465Å, but calculations converge to a non-spin-polarized state). This is in good agreement with the experimental value of 2.507Å in the ferromagnetic state.

Table 1. Results of the cobalt ferromagnetic and antiferromagnetic calculations

MnF2

Fig. 5. Picture of cell used

For these calculations, MnF2 was constructed using the analogous rutile TiO2 structure, which has the P42/MNM symmetry. The starting value of lattice constants c/a and a, 0.68 and 4.875Å respectively, used in these calculations were found by Dormann et. al [2]. From there, the spins in the atomic unit cell were polarized into a ferromagnetic state and energy calculations with varying energy cutoffs and k-samplings were performed until the results converged to a 0.01eV tolerance. The psuedopotential was calculating using the 1s2 electrons for the inner shells of fluorine and 2s2 2p5 as the valence shells. For manganese, the 1s2 2s2 2p6 electrons were used for the inner shells and 3s2 3p6 3d5 4s2 for the valence shells.

Fig. 6. Energy cutoff calculations for MnF2

Fig. 7. k-space Sampling calculations for MnF2

Calculations were made by shifting the lattice constant a in 0.1Å increments to get a rough idea of the lattice constant and 0.01Å measurements to take a finer sweep around the minimum energies of the two states. Ultimately, MnF2 converged to a lattice constant of approximately 4.930Å in the favored antiferromagnetic state and 4.945Å in the less favorable ferromagnetic state. This is relatively close to the experimental value of 4.875Å in the antiferromagnetic state.

Fig. 8. Graph of energies vs. lattice constants for the MnF2 lattice with ferromagnetic and antiferromagnetic spin polarization.

Table 2. Table of energies vs. lattice constants for MnF2.

MnBi

Fig. 9. Picture of cell used

For these calculations, MnBi was constructed using the analogous NiAs structure, which has the P63/MMC symmetry. The starting value of lattice constants c/a and a, 1.426 and 4.270Å respectively, used in these calculations were found by Koyama et. al [3]. From there, the spins in the atomic unit cell were polarized into a ferromagnetic state and energy calculations with varying energy cutoffs and k-samplings were performed until the results converged to a 0.01eV tolerance. The psuedopotential was calculating using the 1s2 2s2 2p6 3s2 3p6 3d10 4s2 4p6 4d10 4f14 5s2 5p6 electrons for the inner shells of bismuth and 5d10 6s2 6p3 as the valence shells. For manganese, the 1s2 2s2 2p6 electrons were used for the inner shells and 3s2 3p6 3d5 4s2 for the valence shells.

Fig. 10. Energy cutoff calculations for MnBi

Fig. 11. k-space Sampling calculations for MnBi

Next, the lattice constant was varied in 0.05Å increments to find an approximate energy minimum for the ferromagnetic and antiferromagnetic states. The lattice constant was then incremented in finer 0.01Å increments around the minimum value to obtain more accurate results. Ultimately, the lattice constant for both states came out to 4.23Å, with the ferromagnetic state being favored over the antiferromagnetic state. This is extremely close to the experimental value of 4.27Å in the ferromagnetic state.

Fig. 12. Graph of energy vs. lattice constant for the two ordered magnetic states of MnBi

Table 3. Table of energy vs. lattice constant for the two ordered states of MnBi

Discussion

These results all agree with the experimentally accepted lattice constant values for the three materials to a good degree of accuracy. In addition, they correctly determined that Co and MnBi are ferromagnetic, while MnF2 is antiferromagnetic. However, the effect of spin polarization in these materials is quite small – on the order of <1eV – suggesting that the energy change is relatively small even for materials considered to be ferromagnets and antiferromagnets.This provides some justification for ignoring spin polarization in our previous calculations, such as those involving copper, which is diamagnetic.

Another interesting note is that the antiferromagnetic calculations required breaking the initial lattice symmetry. This resulted in the antiferromagnetic calculations taking approximately 3-4 times longer than the ferromagnetic calculations; for example, the a = 4.18Å MnBi calculations  took 1701s for the P1 antiferromagnetic state versus 458s for the P63/MMC ferromagnetic state.

Though all the antiferromagnetic calculations in this post were performed with only P1 symmetry, the author realized afterwards that Materials Studio enables one to find a higher symmetry state factoring in the spins. By going into build->Symmetry->Find Symmetry->Options, one can force Materials Studio to consider the formal spins of each atom when searching for symmetries. Using the symmetry finder with these settings, one can reduce the  a = 4.18Å MnBi antiferromagnetic calculation time to 295s, a nearly 6-fold decrease in calculation time.

References

In this post, the phase diagram of lithium and lithium hydride is predicted using energy calculations from DFT and thermodynamics in the grand canonical ensemble. This is achieved by minizming the grand potential as a function of pressure and temperature. That is, the configuration (either Li or LiH) that has the lowest grand potential is the stable configuration of the system. The grand potential is defined as [1][5]

\begin{equation} \Omega=E-TS-\sum_i \mu_i N_i \end{equation}

where the sum is over the chemical species present in the system. To determine the equilibrium phase, we should compare \( \Omega_{Li} \) and \( \Omega_{LiH} \). However, assuming that the number of lithium atoms and the entropy does not change between the two phases, we can instead compare the following quantities:

\begin{equation} \Omega’_{Li} =E_{Li} \\ \Omega’_{LiH}=E_{LiH}-\mu_H N_H \end{equation}

We have in mind a system which is in chemical equilibrium with gaseous molecular hydrogen. Therefore, we can relate the chemical potential of the hydrogen in LiH to the chemical potential of the gaseous hydrogen (which will be treated as an ideal gas) [1]:

\begin{equation} \mu_H=\frac{1}{2} \mu_{H_2} =\frac{1}{2}\big(E_{H_2}+\tilde{\mu}_{H_2}(T,p_0) +kT \ln(p/p_0)\big)\end{equation}

where \( \tilde{\mu}_{H_2}(T,p_0) \) is the difference in chemical potential between the temperatures \( 0~K \) and \( T \) at the reference pressure \( p_0 \). The phase boundary occurs when the expressions in \( (2) \) are equal. This gives an expression for the phase boundary:

\begin{equation} p=p_0e^{\frac{1}{kT}\big[E_{LiH}-E_{Li}-E_{H_2}-\tilde{\mu}_{H_2}(T,p_0)\big] }\end{equation}

\( \tilde{\mu}_{H_2}(T,p_0) \) can be determined using data tabulated at [2] as well as the relation between the chemical potential and enthalpy (\(H\)) and entropy (\(S\)) per molecule:

\begin{equation} \tilde{\mu}=\big[H(T)-H(T_r)\big]-\big[H(0)-H(T_r)\big]-TS(T)  \end{equation}

The energies are found via DFT calculations.

Calculation of Energies

Energies of the Li (bcc crystal) LiH (NaCl structure) and gaseous hydrogen were calculated using DFT calculations in Materials Studio. The energies of the crystals (Li and LiH) were calculated using geometry optimization in CASTEP [3]. The functional used were GGA-PBE [3]. The energy cutoff was set to \( 600~eV \) with \( k \)-point set of \(10 \times 10 \times 10\). Full optimization of the unit cell, in each case, was used. Energy, max. force, max. stress, and max. displacement were \( 5~\mu eV/\mathrm{atom}\), \(.01~eV/Å\), \(.02~GPa\), and \(5\cdot 10^-4 Å\). The energy of gaseous hydrogen was calculated using DMol3 geometry optimization using the LDA-PWC functional. The energy, max. force, and max. displacement were \(10~\mu \mathrm{Ha}\), \(.002 ~\mathrm{Ha}/Å\), and \(.005~Å\). The results of each calculation are reported below:

\begin{align} E_{H_2}&=-31.6746~eV \\ E_{Li}&=-404.8750541~eV \\ E_{LiH}&=-438.267581~eV \end{align}

Phase Diagram

The phase diagram is calculated by setting the grand potentials equal to one another and using data from [2] as well as the energies calculated above. The result is shown below

Figure 1: Phase diagram of Li and LiH.

At atmospheric pressure, LiH is known to dissociate into Li and hydrogen gas at about 1200 K [6]. This fact appears to agree with the phase diagram above and provides a simple accuracy check of this calculation.

References

[1] D. Chandler Modern Statistical Mechanics (Oxford University Press, 1987)

[2] https://janaf.nist.gov/tables/H-050.html

[3]  Clark, et. al., First principles methods using CASTEP, (Z. Kristall, 2005), Vol. 220, P. 567-570

[4]  John P. Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized Gradient Approximation Made Simple”, Phys. Rev. Lett. 77, 3865 – Published 28 October 1996; Erratum Phys. Rev. Lett. 78, 1396 (1997)

[5] D. Sholl and J. Steckel, Density Functional Theory: A Practical Introduction. (Wiley 2009)

[6] Johnson, D.A. Metals and Chemical Change, Volume 1 (Royal Society of Chemistry) p. 167