Quaternary Hydrides Pd_1-y-zAg_yCu_zH_x Embedded Atom Method Potentials for Hydrogen Energy Applications

Chaonan Zhang ^1,†, Robert Fuller ^1,†, Iyad Hijazi ^1,†,*

Weisberg Department of Mechanical Engineering, Marshall University, Huntington, United States

† These authors contributed equally to this work.

* Correspondence: Iyad Hijazi

Academic Editor: Alfonso Chinnici

Special Issue: Hydrogen Energy: Sustainable Production, Storage and Utilisation

Received: November 10, 2020 | Accepted: January 12, 2021 | Published: January 25, 2021

Journal of Energy and Power Technology 2021, Volume 3, Issue 1, doi:10.21926/jept.2101006

Recommended citation: Zhang C, Fuller R, Hijazi I. Quaternary Hydrides Pd1-y-zAgyCuzHx Embedded Atom Method Potentials for Hydrogen Energy Applications. Journal of Energy and Power Technology 2021; 3(1): 006; doi:10.21926/jept.2101006.

© 2021 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

The palladium hydride (Pd-H) system has attracted various studies [1,2,3,4,5,6,7,8]. Palladium (Pd) possesses the capacity to absorb a considerable amount of hydrogen (H) at room temperature and atmospheric pressure. It has excellent hydrogen selectivity, permeability and H diffusivity with high mobility [9,10,11,12]. Therefore, it has diverse applications such as fuel cells, hydrogen storage, refrigeration, hydrogen separation and purification, nuclear radiation adsorption and catalytic converter [9,10,11,12,13,14,15].

However, the coexistence of low H alpha (α) and high H beta (β) phases in Pd hydride causes a miscibility gap. The strain caused by the lattice mismatch between the two phases increases the possibility of mechanical failure [16,17,18]. In addition, Pd is vulnerable to hydrogen embrittlement after some cycles of α phase and β phase transformations in Pd membranes or absorption/desorption cycles in hydrogen storage [18,19]. Pd membranes are also prone to H₂S and CO poisoning [20,21,22]. Finally, palladium is an expensive noble metal, with a price currently higher than gold.

Pd can be alloyed with silver (Ag) and copper (Cu) to solve the above listed problems and reduce costs [23,24,25]. Although H is practically insoluble in Ag, Pd-Ag alloys possess higher H solubility than Pd [24]. As with Pd, Pd-Ag alloys are 100% selective when absorbing H, and therefore, are also ideal for separating H from mixture of gases [18]. Pd-Ag alloys also possess excellent H permeability, with the greatest permeability observed in the Pd₇₇Ag₂₃ structure [23,26]. Depending on Ag concentration, H diffusion rate can be also increased in Pd-Ag alloys as compared with pure Pd [5]. Pd-Ag alloys are reported to possess more steady mechanical properties than pure Pd, and better resistance to H embrittlement [18], however, strong poisoning by H₂S remains an issue [27].

Pd-Cu alloys are also known as one of the most effective materials in hydrogen separation field, as Cu is a much cheaper metal than Pd. Pd-Cu alloys also have better mechanical properties, high thermal stability [28], avoidance of hydrogen embrittlement at room temperature, improved sulfur poisoning resistance than Pd and Pd-Ag alloys [29], and more resistance to carbon monoxide poisoning than pure Pd [30]. In addition, Pd-Cu alloys are also completely selective for hydrogen [31], and H diffusion rate in body-centered-cubic (bcc) Pd-Cu alloys is faster than those in pure Pd and face-centered-cubic (fcc) Pd-Cu alloys [32]. However, H permeation rate in bcc Pd-Cu has been found to be slower than in pure Pd, but faster than in fcc Pd-Cu [33], and H solubility in Pd-Cu alloy is lower than pure Pd, and decreases strongly as Cu content increases [33].

Since Pd-Ag alloys have excellent H selectivity, permeability and solubility, and Pd-Cu alloys have excellent resistance to sulfur and carbon monoxide poisoning and hydrogen embrittlement, good mechanical properties, broader working conditions than pure Pd and relatively low cost, one way to optimize the overall performance and substantially lower the cost is alloying Pd with Ag and Cu [34,35,36,37,38]. To accelerate the search for Pd-Ag-Cu alloys with optimal hydrogen selectivity, permeability, diffusivity, absorption, thermal stability and resistance to sulfur and CO poisoning, an accurate and efficient modeling of the Pd-Ag-Cu-H system is needed. The most reliable simulation technique utilizes first principle (ab initio) calculations, but its high processing costs make it infeasible for simulations containing a large number of atoms. Molecular dynamics (MD) simulations using the embedded atom method (EAM) offers an efficient way to investigate alloys with larger atomic structures. EAM is well suitable to model binary and ternary hydrides with metallic crystal structures and interstitial H atoms. An EAM Pd-H potential which can predict the miscibility gap was formulated by Zhou et al. [4]. The Pd-H EAM potential was then expanded into Pd-Ag-H ternary potential by Hale et al. [5]. Their Pd-Ag-H potential predicted smooth changes in structures parameters, elastic properties and energy with increasing H concentrations, H sites occupation shift and the disappearance of the miscibility gap by the addition of Ag at 300 K. However, their model was built on Foiles et al.’s Pd potential, which is available in a tabular form but does not include a full explanation of the analytical form and parameters [39], and therefore impeding further improvement of the Pd-Ag-H ternary system.

In previous work, Pd-H EAM potentials with fewer fitting parameters than Zhou et al. [4] were fitted by Hijazi et al. [40] that can predict many of the properties of the Pd-H structures accurately. In this paper, the previously developed Pd-H potentials are expanded into the Pd-Ag-Cu-H quaternary potentials. First principles simulations using SIESTA were carried out to obtain the fitting parameters. The fully analytical potentials are fitted utilizing the central atom method without performing time-consuming MD simulations. The Pd-Ag-Cu-H EAM potentials can capture the miscibility gap behavior; preferred H occupancy sites; and predicts the trends for the lattice constants, cohesive energies, and elastic properties during MD simulations.

2. The Atomic Potentials

A total of 18 functions are needed to create the Pd-Ag-Cu-H quaternary atomic EAM potentials. Which include 4 embedding energy functions, F_Pd, F_Ag, F_Cu and F_H, 4 electron density functions, ρ_Pd, ρ_Ag, ρ_Cu and ρ_H, and 10 pair functions, φ_Pd–Pd, φ_Ag–Ag, φ_Cu–Cu, φ_H–H, φ_Pd–Ag, φ_Pd–Cu, φ_Ag–Cu, φ_Pd–H, φ_Ag–H and φ_Cu–H. In EAM, each atom is embedded into a lattice that include all host atoms. The pair potential between atoms, and the energy related to embedding an atom inside the host lattice is modeled. The total energy E_tot of an EAM system is given by [41].

\[E_{tot}=\sum_{i=1}^NF_i(\rho_i)+\frac{1}{2}\sum_{i=1}^N\sum_{\scriptstyle j=1\atop \scriptstyle j\neq i}^N\varphi_{ij}(r_{ij})\tag{1} \]

\[\rho_i=\sum_{j\neq i}^Nf_j(r_{ij})\tag{2} \]

where ρ_i is the electron density for atom i, F_i is the embedding energy, ϕ_ij is the pair potential between atom i and atom j, r_ij represents the distance from atom i and to atom j, f_j is the electron density function of distance from the center of atom j. The EAM model by Hijazi et al. was used in fitting the Pd-Pd, Ag-Ag, Cu-Cu, Pd-Ag, Pd-Cu and Ag-Cu potentials [42,43,44,45,46]. The embedding function is represented by:

\[F(\rho)=F(\rho_e )\left[1-\eta\ ln⁡\left(\frac{\rho}{\rho_e}\right) \right] \left(\frac{\rho}{\rho_e}\right)^\eta\tag{3} \]

The host electron density is given by:

\[f=f_e\ e^{-\chi(r-r_e)}\tag{4} \]

where f_e is a scaling factor that can be obtained from f_e = E_c/Ω, Ω is the atomic volume and E_c is the cohesive energy, r_e is the equilibrium closest distance, and χ is a fitting parameter. The pair potential function is the modified potential created by Rose et al. [47] and has the form:

\[\varphi=-\varphi_e\left[1+\delta\left(\frac{r}{r_e} -1\right)\right] e^{-\beta(\frac{r}{r_e} -1) }\tag{5} \]

where ϕ_e, δ, and β are the 3 adjustable parameters. Therefore, for an fcc metal, there are 6 fitting parameters χ, ϕ_e, δ, β, η, and ρ_e.

The analytical expressions proposed by Zhou et al. [4,5] was used in fitting the H-H, Pd-H, Ag-H and Cu-H potentials. The pair potentials have the form of the generalized Morse potential and is given by:

\[\varphi(r)=D\left(\beta e^{-\alpha(r-r_0)}-\alpha e^{\beta(r-r_0 )}\right)\tag{6} \]

where D, α, β, and r₀ are fitting parameters, r₀ defines the interatomic spacing between two atoms, and D (β − α) is the binding energy. The H electron density function is given by:

\[\rho_H (r)=C_H e^{-\delta_H⋅r}\tag{7} \]

which has 2 fitting parameters C_H and δ_H, while the embedding function for H has the form:

\[F_{H,u} (\rho)=-C_H\cdot \left( \begin{aligned}&\frac{1}{2+d_H}\cdot(\rho+\varepsilon_H)^{2+d_H}-\frac{a_H+b_H}{1+d_H}\cdot(\rho+\varepsilon_H)^{1+d_H} \\ &+\frac{a_H\cdot b_H}{d_H}\cdot (\rho+\varepsilon_H)^{d_H} \end{aligned}\right)\tag{8} \]

where a_H, b_H, c_H, d_H are fitting parameters, and ε_H = 0.0540638.

The EAM total energy is a linear summation of the embedding energy and the pair potentials. A unique feature of the elemental EAM potential is that it will not change due to the transformation of the embedded energy functions. Thus, the embedding and pair potentials for Pd-H can be transformed utilizing the two equations below:

\[F_i^{Final} (\rho_i )=F_i^{initial} (\rho_i )+k\rho_i\tag{9} \]

\[\varphi_{ij}^{Final}(r_{ij})=\varphi_{ij}^{initial}(r_{ij})+2k\rho_i (r_{ij})\tag{10} \]

k represents an arbitrary constant. The embedding and pair potentials for Pd-H were thus converted in this pattern according to the method of Zhou et al. [4].

2.1 Palladium Silver Copper Alloys

2.1.1 Ag and Cu Fitting Parameters

For the pure metals Ag and Cu, the EAM potentials were fitted the same way as previously done for Pd [40]. For each metal, the six fitting parameters included a_o, E_c, C₁₁, C₁₂, C₄₄, and E_vf. Where a_o is the lattice constant, E_c is the cohesive energy, C₁₁, C₁₂, and C₄₄ are three elastic constants, and E_vf is the vacancy formation energy. Table 1 below lists the Pd, Ag and Cu fitting parameters.

Table 1 Fitting parameters for Pd, Ag and Cu.

As can be seen from Table 2, the calculated fitting results were almost identical to the experimental values [42] and those obtained by Foiles et al. [48]. In addition, as with Pd, the plots of cohesive energy vs. lattice constant for Ag and Cu were also in very good consistency with the equation of state obtained from Rose et al. [47], as shown in Figure 1.

Table 2 Fitting results for Ag and Cu.

Figure 1 Ag and Cu Cohesive Energy and Rose et al. Equation of State [47].

To account for the pair potential interactions for Pd-Ag, Pd-Cu and Ag-Cu alloys, the mixing rule between a type-a and a type-b atom interaction introduced by Johnson [49] was applied, and is given by the equation:

\[ \phi_{ab}(r)=\frac{1}{2}\left[\frac{f_b (r)}{f_a (r)}\phi_{aa}(r)+ \frac{f_a (r)}{f_b (r)}\phi_{bb} (r)\right]\tag{11} \]

For each type in the alloy, the electron density parameter can be calculated from the equation f_e=S(E_c/Ω), where Ω is the atomic volume and S is a scaling factor, with S = 1 for pure metals. For type-a atom as a host (solvent) and type-b as impurity (solute), the unrelaxed dilute limit heat of solution can be determined by the five steps given below:

(a) Host removal: $F^b (\bar{\rho}^a)+\sum_{i\neq1} \phi^{ab}(r_{li}^a)$

(b) Add impurity:$- F^a (\bar{\rho}^a)- \sum_{i\neq1} \phi^{aa}(r_{li}^a)$

(c) Adjust neighbors:$- \sum_{i\neq1} F^a (\bar{\rho}^a)- \sum_{i\neq1} F_i^a(\bar{\rho}^a-f^a(r_{li}^a)+f^b(r_{li}^a))$

(d) Adjust cohesive energy:$-E_c^a+E_c^b$

(e) Relaxation energy:$E_r=[1.167(\Omega_b/\Omega_a-1)]^2$

where is $\bar{\rho}^a$ the expression of electron density for type-a atom, r^a is the distance to its closest neighbor and E_r is the drop in total energy caused by relaxation and is predominantly dependent on the unit cell volume mismatch.

The electron density scaling factors for type-a and type-b atoms, S_a and S_b, for the Pd-Ag, Pd-Cu and Ag-Cu pair potentials, obtained from fitting the experimental heat of solutions, are listed in Table 3 along with the calculated heat of solutions values for each metal. After applying the relaxation energy E_r, the values for the relaxed heat of solution are very consistent with experimental obtained data and overall better than those obtained by Foiles et al. [48] and Hijazi et al. [42].

Table 3 Pd-Ag, Pd-Cu and Ag-Cu heat of solution and scaling factors from fitting.

2.1.2 Validation

The fitted parameters for Pd-Pd, Ag-Ag, Cu-Cu, Pd-Ag, Pd-Cu and Ag-Cu have been applied to create a tabulated EAM potential file in DYNAMO setfl format for the ternary Pd-Ag-Cu system. Utilizing the tabulated EAM potential file, MD annealing simulations with a Nose-Hoover NPT ensemble from 500 K to 1 K in 100 ns were performed for Ag-Ag, Cu-Cu, Pd-Ag, Pd-Cu and Ag-Cu structures using a LAMMPS script code [50]. A molecular statics (MS) simulation utilizing the conjugate gradient (cg) minimization method was then applied after each MD simulation. The MD simulation results for Ag-Ag, Cu-Cu were in excellent agreement with the calculated fitting results, as can be seen from Table 2, and the MD results proved the reliability of the Pd-Ag, Pd-Cu and Ag-Cu EAM potentials as illustrated in Figure 2, Figure 3 and Figure 4.

Figure 2(a) and (b) show our lattice constant and cohesive energy results for the Pd_xAg_1-x and Pd_xCu_1-x (0 ≤ x ≤ 1) structures are almost identical with the experimental data [51]. For the Pd_xAg_1-x structures, the lattice constant results are closer to the experimental values than the results calculated using the Hale et al.’s EAM potential with Morse pair function [5]. On the other hand, the density functional theory (DFT) data collected from Løvvik et al. [52] reveal a similar trend but tend to overestimate the values for all compositions. Figure 2(b) shows that our cohesive energy values for Pd-Ag are very much in line with the values that are derived using the Hale et al. with Morse function [5]. The Hale et al. EAM potentials were obtained from the Interatomic Potential Repository [53]. However, their EAM potential with the hybrid model produced erratic results and was not included with the Figures.

For the Pd_xCu_1-x structures, the lattice constant values from MD simulations are almost identical with the experimental values and those obtained by Kart et al. [28] as shown in Figure 2(a). For the cohesive energies for Pd_xCu_1-x, our MD values have an increasing trend similar to the data obtained by Kart et al. [28], as can be seen in Figure 2(b).

Figure 2 Pd_1-xAg_x and Pd_xCu_1-x structures a_o and E_c results obtained from MD simulations, experimental and DFT calculations [28,55,52,51].

In Figure 3(a) and (b) the values for the elastic constants C_11, C₁₂ from MD simulations for Pd_1-xAg_x and Pd_1-xCu_x structures show consistent trend with the DFT calculations and the results from Hale et al. and Kart et al. [54,5,28]. The bulk modulus for Pd_1-xAg_x and Pd_1-xCu_x structures obtained from MD simulations match the softening trends predicted by the DFT calculations as well [54,52,28] and match the given experimental data at the edge of the composition range quite closely, as shown in Figure 3(c). It’s worth noticing that the Hale et al. [5] EAM potential overestimates C_11, C₁₂ and bulk modulus for the pure Pd metal, as can be seen at the left edges of Figure 3(a), (b) and (c).

Figure 3 Pd_1-xAg_x and Pd_xCu_1-x alloys C₁₁, C₁₂ elastic constant and B_m results obtained from MD simulations, experimental and DFT calculations [28,56].

In Figure 4, the values for the elastic constants C₄₄ and C’ from MD simulations for Pd_1-xAg_x show that our results are closer to the experimental data at the edge of the composition range than those of Hale et al. [5]. As with the Hale et al. EAM Morse model [5], our Pd_1-xAg_x potential underestimates C₄₄ relative to the DFT results and have an overall decreasing trend. For Pd_1-xCu_x, our C₄₄ and C’ values have a slightly increasing trend, with C’ for Pd being underestimated, but still more consistent with the experimental data than the results from Kart et al. [28].

Figure 4 Pd_1-xAg_x and Pd_xCu_1-x alloys C₄₄ and C’ elastic constant results obtained from MD simulations, experimental and DFT calculations [28,56].

3. Palladium Silver and Palladium Copper Hydrides

3.1 DFT Calculations

Since H is almost insoluble in Ag [24,57], and no experimental fitting data was found for Ag-H and Cu-H systems [40], therefore, the Pd-Ag-H and Pd-Cu-H properties were used as fitting data to fit the φ_Ag–H and φ_Cu–H pair functions. However, only a limited experimental and ab initio data were available to utilize a full H concentration in the fitting procedure. Hale et al. obtained their fitting data by utilizing DFT calculations for the Pd-Ag-H system [5], but failed to provide the lattice constant values, and only the cohesive energy values were given. Wei et al. performed DFT calculations on Pd-Cu-H phase stability, heat of formations and elastic property based on generalized gradient approximations (GGA) for the range of hydrogen concentration 0≤x≤0.5, but they failed to report the exact values for the lattice constants and the cohesive energies [58]. In this paper, the open source SIESTA software was used to perform ab initio simulations to get full fitting data for the Pd-Ag-H and Pd-Cu-H structures. The SIESTA pseudopotentials were obtained from the Abinit's Fritz-Haber-Institute (FHI) pseudo database [59]. The local density approximation (LDA) method with Ceperley–Alder exchange and correlation form using the norm-conserving Troullier–Martins scheme was utilized in the pseudopotentials. Valence states were described using double zeta-polarized (DZP) basis sets with split-valence scheme for multiple-zeta. The ab initio simulations were conducted with a dense 18×18×18 Monkhorst–Pack grid, a cutoff energy of 100 Ry, a 25 K electronic temperature, and electron spin polarization during the DFT calculations. For our Pd-Ag-H and Pd-Cu-H structures, the calculations utilized periodic boundary conditions with a unit cell of 3 Pd atoms, 1 Ag or Cu atom, and 1 to 4 H atoms at different locations.

During the DFT simulations, the Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x structures were constructed with five different H concentrations: x = 0, 0.25, 0.50, 0.75, and 1.00. In the Pd-Ag and Pd-Cu fcc lattice, H atoms were located in three different interstitial positions. As shown in Figure 5, these positions included the octahedral (O) positions, in which O1 represent a body center position and O2 an edge center position, and the tetrahedral (TE) positions.

Figure 5 H atoms (white) interstitial different positions in the Pd atoms (gray) and Ag or Cu atoms (black) lattice.

As with Løvvik et al. Pd_1-xAg_x DFT results from Figure 2(a) and (b) [52], the lattice constant and cohesive energy results from our DFT simulations for the Pd-Ag-H and Pd-Cu-H structures were also overestimated in comparison to the available experimentally obtained data for the Pd, Ag, Cu, PdH_0.50, PdH_1.00, Pd_0.75Ag_0.25 and Pd_0.75Cu_0.25 structures. The calculated DFT values can be shifted, if multiplied with a selected factor, to make it consistent with the experimental data [60]. Equations 12 to 15 describe the shifting procedure for the cohesive energies for Pd_0.75Ag_0.25H_x structures, which have been applied in a similar manner to the Pd_0.75Cu_0.25H_x case by replacing Ag atoms with Cu. The shifting procedure was also applied in a similar fashion to the lattice constants case. The shifting data for Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x structures are given in Table 4 and the shifting factors in Table 5.

\[ \begin{aligned}(Cohesive\ Ener&gy_{Pd_{0.75} Ag_{0.25} H_x})_{SIESTA\ shifted}=(C_{Pd_{0.75} Ag_{0.25} }+x\cdot C_H)\\ & \cdot (Cohesive\ Energy_{Pd_{0.75} Ag_{0.25} H_x })_{SIESTA}\end{aligned}\tag{12}\]

\[C_{Pd_{0.75} Ag_{0.25}}=\frac{(Cohesive\ Energy_{Pd_{0.75} Ag_{0.25}})_{Experimental}}{(Cohesive\ Energy_{Pd_{0.75} Ag_{0.25}})_{SIESTA}}\tag{13} \]

\[ C_H=2\cdot \bigg\{ \frac{(Cohesive\ Energy_{PdH_{0.50}})_{Experimental}}{(Cohesive\ Energy_{PdH_{0.50}})_{SIESTA}} -C_{Pd}\bigg\} \tag{14}\]

\[C_{Pd}= \frac{(Cohesive\ Energy_{Pd})_{Experimental}}{(Cohesive\ Energy_{Pd})_{SIESTA}} \tag{15} \]

Table 4 Experimental values used in shifting ab initio data.

Table 5 Shifting factors for Pd-Ag-H and Pd-Cu-H ab initio data.

The shifted cohesive energy values for the Pd_0.75Ag_0.25H_x structures are close to those obtained by Hale et al. DFT calculations [5] and have similar trend, as can be seen in Figure 6. Figure 6 also shows that the Pd_0.75Cu_0.25H_x shifted data have a similar trend to the Pd_0.75Ag_0.25H_x shifted data, but with lower cohesive energy values.

Figure 6 Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x shifted cohesive energies compared with Hale et al. DFT results.

The shifted values for the cohesive energy and lattice constant obtained from the DFT simulations are listed in Table S1 in the supplementary data. The shifted cohesive energies for 7 OC structures and 7 TE structures with 4 different H concentrations, were used in fitting the φ_Ag–H and φ_Cu–H pair potential functions during the fitting process. The φ_Ag–H and φ_Cu–H pair functions take the same generalized Morse potential mathematical form, as used previously in the Pd-H interaction [46]. Since a third atom type was added to the binary Pd-H structures to create the ternary Pd-Ag-H and Pd-Cu-H structures, all Pd-H potentials and property equations utilized in the fitting procedure were expanded by adding a central atom expression as a third type [40]. For the ternary system, the cohesive energy equation has an additional host term and is given by:

\[ E_c=\frac{1}{x+y+z}=\left[ \begin{aligned}&x\left( F_{a,i}(\rho_{a,i})+\frac{1}{2}\sum_{{j=1}\atop{j≠i}}^{N_a}\varphi_{a-a,ij}(r_{ij})+ \frac{1}{2}\sum_{{j=1}\atop{j≠i}}^{N_b}\varphi_{a-b,ij}(r_{ij})+\frac{1}{2}\sum_{{j=1}\atop{j≠i}}^{N_c}\varphi_{a-c,ij}(r_{ij})\right)+\\&y\left( F_{b,i}(\rho_{b,i})+\frac{1}{2}\sum_{{j=1}\atop{j≠i}}^{N_b}\varphi_{b-b,ij}(r_{ij})+ \frac{1}{2}\sum_{{j=1}\atop{j≠i}}^{N_a}\varphi_{b-a,ij}(r_{ij})+\frac{1}{2}\sum_{{j=1}\atop{j≠i}}^{N_c}\varphi_{b-c,ij}(r_{ij})\right)+\\&z\left( F_{c,i}(\rho_{c,i})+\frac{1}{2}\sum_{{j=1}\atop{j≠i}}^{N_c}\varphi_{c-c,ij}(r_{ij})+ \frac{1}{2}\sum_{{j=1}\atop{j≠i}}^{N_a}\varphi_{c-a,ij}(r_{ij})+\frac{1}{2}\sum_{{j=1}\atop{j≠i}}^{N_b}\varphi_{c-b,ij}(r_{ij})\right) \end{aligned} \right]\tag{16} \]

\[\rho_{a,i}=\rho_{a-a,i}+\rho_{a-b,i}+\rho_{a-c,i}=\sum_{{j=1}\atop{j≠i}}^{N_a}f_a(r_{ij} )+\sum_{{j=1}\atop{j≠i}}^{N_b}f_b(r_{ij} )+\sum_{{j=1}\atop{j≠i}}^{N_c}f_c(r_{ij} )\tag{17} \]

\[\rho_{b,i}=\rho_{b-a,i}+\rho_{b-b,i}+\rho_{b-c,i}=\sum_{{j=1}\atop{j≠i}}^{N_a}f_a(r_{ij} )+\sum_{{j=1}\atop{j≠i}}^{N_b}f_b(r_{ij} )+\sum_{{j=1}\atop{j≠i}}^{N_c}f_c(r_{ij} )\tag{18} \]

\[\rho_{c,i}=\rho_{c-a,i}+\rho_{c-b,i}+\rho_{c-c,i}=\sum_{{j=1}\atop{j≠i}}^{N_a}f_a(r_{ij} )+\sum_{{j=1}\atop{j≠i}}^{N_b}f_b(r_{ij} )+\sum_{{j=1}\atop{j≠i}}^{N_c}f_c(r_{ij} )\tag{19} \]

where a, b, and c are three different type of atoms, and x, y, and z are the concentrations for each type of atom in the structure respectively. A constrained nonlinear optimization MATLAB code was used during the fitting procedure to obtain the fitting parameters. Table 6 lists the parameters for the Ag-H and Cu-H from fitting, and the Pd-H parameters from our previous work [40].

Table 6 Fitting data used for Pd-H, Cu-H and Ag-H.

The previously obtained fitting data for the H-H potential are also included in Table 7 [40]. The two body potential functions used in our Pd-Ag-Cu-H model are plotted in Figure 7 and Figure 8. The calculated cohesive energies for Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x structures from fitting are consistent with the fitting data for most of the H concentrations, but the results start to diverge from the fitting data at high H concentrations as can be seen in Figure 9(b).

Table 7 Fitting data for H-H.

Figure 7 H-H, Pd-H, Ag-H and Cu-H two body potential functions.

Figure 8 Pd-Pd, Ag-Ag, Cu-Cu, Pd-Ag, Pd-Cu and Ag-Cu two body potential functions.

Figure 9 Pd_1.00H_x, Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x a_o and E_c results obtained from MD simulations and fitting.

3.2 Results

To test the reliability of the Pd-Ag-Cu-H potentials, a tabulated potential file in DYNAMO setfl format was generated utilizing the final fitting parameters. Utilizing LAMMPS and the tabulated potential file, MD annealing simulations with a Nose-Hoover NPT ensemble from 500 K to 1 K in 100 ns with random hydrogen atom positions were performed for the Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x structures. A molecular statics (MS) simulation utilizing the conjugate gradient (cg) minimization method was then applied after each MD simulation. Ten sets of data were generated for each H composition and their average values were taken to ensure accuracy.

3.3 Lattice Constant and Cohesive Energy Results

The stress triggered by variation in the lattice constants in regions with different H concentrations at equilibrium is of important consideration [4]. Therefore, the influence of H concentration on the equilibrium lattice constant was investigated. As can be seen in Figure 9(a), the lattice constant values obtained from MD simulations for Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x structures are almost identical with the DFT results used in the fitting process. The lattice constant plots show an increasing trend with increasing H concentrations, similar with the DFT calculated data and the Pd_1.00H_x results in our previous work [40]. The increasing trend from our plots is also consistent with the results using the Hale et al. EAM potential with the Morse function for Pd_0.75Ag_0.25H_x [5], and the DFT simulation results for the fcc Pd-Cu-H (O1) from Wei et al. [58].

For the cohesive energy, the results for Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x structures shown in Figure 9(b), were in good agreement with our fitting results and DFT data, and also have a similar trend to the Pd_1.00H_x plot [40]. For the Pd_0.75Ag_0.25H_x structures, our MD results are closer to our DFT data than the Hale et al. EAM Morse potential MD results with respect to their own DFT fitting data [5].

3.4 Elastic Constants and Bulk Modulus

Using the relaxed Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x structures obtained from our MD + MS simulations and utilizing an adaptation of the general-purpose LAMMPS script written by Aiden P. Thompson at Sandia National Laboratories [50], the elastic constants and bulk modulus values were estimated. Figure 10 shows that the elastic constants C₁₁ and C₁₂ and bulk modulus for Pd_0.75Ag_0.25H_x, Pd_0.75Cu_0.25H_x and Pd_1.00H_x display smooth curves with similar trends. Although, as stated previously, that the Hale et al. EAM potential overestimates the bulk modulus for pure Pd, our overall decreasing trend matches well with the results obtained by Hale et al. EAM Morse potential for the Pd_0.75Ag_0.25H_x structures and in a good agreement with our previously obtained results for Pd_1.00H_x [5,40]. The bulk modulus values from our MD simulations for Pd_0.75Cu_0.25H_x also have a similar decreasing trend with those obtained by Wei et al. DFT simulations for fcc Pd-Cu-H (O1) [58]. Other researchers have also documented this softening behavior with increasing H concentrations [61,62,63]. The results for Pd_0.75Ag_0.25H_x structures obtained from Hale et al. Hybrid potential yielded an unstable trend but still had an overall similar softening trend [5,53]. As for the Pd_0.75Cu_0.25H_x compositions, our simulation results for C₁₁ and C₁₂ also have a similar smooth overall decreasing trend to those obtained by Wei et al. [58].

Figure 10 Pd_1.00H_x, Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x C₁₁, C₁₂ elastic constant and B_m results from MD simulations.

Figure 11 shows the elastic constants C₄₄ and C’ for Pd_1.00H_x, Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x alloys. As previously obtained for Pd_1.00H_x [40], the plot values for C₄₄ for Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x decrease while the shear elastic constants C’ increase with increasing H composition.

Figure 11 Pd_1.00H_x, Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x C₄₄ elastic constant and C’ shear elastic constant results from MD simulations.

In addition, our elastic constants values for the various Pd_1.00H_x, Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x alloys shown in Figure 10 and Figure 11 satisfy the theory of strain energy for cubic structures [64]. According to strain energy theory, the following formulas can be applied to a mechanically stable cubic: C₁₁ - C₁₂ > 0, C₁₁ + 2C₁₂ > 0 and C₄₄ > 0. From Figure 10 and Figure 11, it can also be seen that the Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x structures have smaller C₁₂ and bigger C₄₄ than the Pd_1.00H_x structures, implying that alloying Pd with Cu or Ag has a significant impact on the elastic constant properties for Pd-Cu-H and Pd-Ag-H phases [58].

3.5 Additional Compositions

To demonstrate the validity of our EAM potentials beyond the Pd, Ag and Cu concentrations utilized during the fitting process, Figure 12 shows the lattice constants and cohesive energies for the Pd_0.50Ag_0.50H_x, Pd_0.50Cu_0.50H_x and Pd_0.50Ag_0.25Cu_0.25H_x hydrides. The Pd₅₀Cu₅₀ structure was chosen on purpose because of its relative similarity to the Pd_52.5Cu_47.5 structure which proved to have the highest H permeability by experimental findings [58]. As can be seen from Figure 12, the lattice constant and cohesive energy results for these additional compositions display a similar trend consistent with our previously obtained results for Pd_1.00H_x, Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x hydrides.

Figure 12 Pd_1.00H_x, Pd_0.75Ag_0.25H_x, Pd_0.50Ag_0.50H_x, Pd_0.75Cu_0.25H_x, Pd_0.50Cu_0.50H_x and Pd_0.50Ag_0.25Cu_0.25H_x a_o and E_c results from MD simulations.

In Figure 13(a), (b) and (c) the MD simulation results for the elastic constants and bulk modulus for Pd_0.50Ag_0.50H_x, Pd_0.50Cu_0.50H_x and Pd_0.50Ag_0.25Cu_0.25H_x are plotted. It can be seen that they also have a consistent trend with our previous results for the Pd_1.00H_x, Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x hydrides. The plots for bulk modulus for Pd_0.50Ag_0.50H_x, Pd_0.50Cu_0.50H_x and Pd_0.50Ag_0.25Cu_0.25H_x have a similar softening trend to our previous results for Pd_1.00H_x, Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x. For Pd_0.50Ag_0.25Cu_0.25H_x, the bulk modulus values are higher than those for Pd_0.50Ag_0.50H_x and closer to Pd_0.50Cu_0.50H_x, indicating that adding Cu has a strengthening impact on the solid solution while adding Ag has a softening effect.

The elastic constant values for C₁₁ and C₁₂ for the various Pd_0.50Ag_0.50H_x, Pd_0.50Cu_0.50H_x and Pd_0.50Ag_0.25Cu_0.25H_x structures shown in Figure 13 also satisfy the strain energy theory for cubic structures [58], implying these structures also possess mechanical stability.

Figure 13 Pd_1.00H_x, Pd_0.75Ag_0.25H_x, Pd_0.50Ag_0.50H_x, Pd_0.75Cu_0.25H_x, Pd_0.50Cu_0.50H_x and Pd_0.50Ag_0.25Cu_0.25H_x C₁₁, C₁₂ elastic constant and B_m results from MD simulations.

3.6 Dynamic Stability

In the Pd_1.00H_x hydride, H atoms tend to take the OC sites in the Pd fcc lattice [65]. The DFT calculation results show that the OC sites in Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x structures are highly energetically favorable to H atoms; this behavior was also observed and reported by other researchers [5,52]. In order to verify the stability for the Pd_1-yAg_yH_x, Pd_1-yCu_yH_x and Pd_1-y-zAg_yCu_zH_x structures using our EAM potentials, structures with TE sites occupied by H atoms were created using LAMMPS, as shown in Figure 14(a). MD simulations were carried out with an NPT ensemble, each TE structure was annealed from 500 K to 1 K for 100 ns, and then followed by cg energy minimization. After each MD simulation, the H atoms moved to lower energy OC sites, as was reported. The resulting structure for a Pd_0.50Ag_0.25Cu_0.25H_x is shown in Figure 14(b).

Figure 14 (a) H atoms (red) in TE sites before simulation. (b) H atoms moved to OC sites after simulation.

3.7 Miscibility Gap and Gibbs Free Energy of Mixing

For the previously studied Pd_1.00H_x hydrides, the miscibility gap was predicted based on the Gibbs free energy of mixing as a function of H concentration [40]. Following the method of Hale et al. [5], the mixing enthalpy term was modified to obtain the Gibbs free energy functions for Pd_1-yAg_yH_x, Pd_1-yCu_yH_{x ­}and Pd_1-y-zAg_yCu_zH_x hydrides as follows:

\[\Delta G^\text{mix}=\Delta H^\text{mix} -\Delta S^\text{mix} \cdot T\tag{20} \]

\[\Delta H^\text{mix}=E_{Pd_{1-y}Ag_yH_x}-2X\cdot E_{Pd_{1-y}Ag_yH_{1.0}} -(1-2X)\cdot E_{Pd_{1-y}Ag_y}\tag{21} \]

where the cohesive energies $E_{Pd_{1-y}Ag_yH_x}, E_{Pd_{1-y}Ag_yH_{1.0}}, and E_{Pd_{1-y}Ag_y}$ were applied, and X is the mole fraction which X= x/(1+x).

\[\Delta S_t=-k_B\cdot \left[X\cdot \ln \left[\frac{X}{(1-X)}\right]+(1-2\cdot X)\ln \left[\frac{1-2X}{(1-X)}\right] \right]\tag{22} \]

where k_B is Boltzmann’s constant.

At 0 K, the Gibbs free energy values in Figure 15, as expected, are all above zero for all structures and various H concentrations, indicating that the average attractive interactions between different atom types are weaker than those between the same atom types. At 300 K, Figure 16 shows that the Gibbs free energy plot for Pd_1.00H_x has two minima at x = 0.034 and 0.95, corresponding to the mole fraction of X = 0.033 and 0.49 and represent the α and β phases, respectively. They describe a miscibility gap region in an alloy, where two phases are more stable than a single one. The Hale et al. EAM Morse model predicts the α and β phases to be X = 0.0 and 0.47 [5]. Experimentally obtained phase boundaries for Pd_1.00H_x at 300 K are x = 0.03 and 0.6, corresponding to mole fraction of X = 0.029 and 0.375. Therefore, our model is in better agreement in predicting the α phase but the β phase is slightly overestimates compared to Hale et al. EAM Morse potential. In Figure 16, our MD results for Pd_1-yAg_yH_x at 300 K show that when Ag concentration increases, the values become more negative relative to the Pd_1.00H_x system, indicating more favorable mixing, and the miscibility gap become narrower. No miscibility gap observed when y = 0.5. At 300 K, the experimental values indicate that the α phase and β phase cease to be distinct at y = 0.25–0.30 for Pd_1-yAg_y [56]. This shows that our EAM potentials are able to detect the miscibility gap, and are consistent with the experimental results regarding the miscibility gap overall behavior as Ag concentration increases. For the Pd_1-yCu_yH_x compositions, experimental data at 303 K indicated that by increasing Cu concentration, the α phase and β phase shift to the right and to the left respectively, the miscibility gap narrows, and finally disappears at y = 0.29 [66,67]. Our values from MD simulations at 300 K in Figure 16 also indicate that adding Cu causes the Gibbs free energy to increase for all H concentrations in comparison to the Pd_1.00H_x structures. At y = 0.25, all calculated energies are positive indicating unfavorable mixing, no miscibility gap observed at y = 0.5. This shows that our model predicts the miscibility gap to disappear at high Cu concentration. For Pd_0.50Ag_0.25Cu_0.25H_x compositions, the Gibbs free energy plot has a similar trend with those obtained from the Pd_1.00H_x structures, indicating that the addition of copper and silver with equal concentration seems to have an opposite effect on the Gibbs free energy and tend to offset each other.

Figure 15 Gibbs free energy plot for different structures at 0 K.

Figure 16 Gibbs free energy plot for different structures at 300K.

4. Conclusion

In this paper, the central atom method was used to fit fully analytical Pd-Ag-Cu-H EAM potentials without utilizing the time-intensive MD simulations during the fitting process. The potentials were efficient in minimizing the objective functions during the fitting calculations, and the number of fitting parameters were reduced compared to previously developed EAM potentials. There were 6 fitting parameters for each Pd-Pd, Ag-Ag and Cu-Cu EAM potential, 2 scaling factors calculated by a mixing rule for each Pd-Ag, Pd-Cu and Ag-Cu pair interaction, 10 fitting parameters for H-H, and 4 for each Pd-H, Ag-H and Cu-H EAM potential. Our MD simulation results validated that these EAM potentials can be applied accurately in further simulations.

The experimentally obtained heat of solutions values were used in fitting the Pd-Ag, Pd-Cu and Ag-Cu pair potentials. The Ag-H and Cu-H EAM potentials were fitted to the cohesive energies for 14 Pd_0.75Ag_0.25H_x and 14 Pd_0.75Cu_0.25H_x structures, obtained from ab initio SIESTA simulations. The MD simulations utilizing LAMMPS demonstrated that our lattice constant and cohesive energy results for Pd_0.75Ag_0.25H_x and Pd_0.75Cu_0.25H_x structures were consistent with the ab initio fitting data for most of the H concentrations. The MD results for the Pd_1-yAg_yH_x, Pd_1-yCu_yH_x and Pd_1-y-zAg_yCu_zH_x structures also demonstrated a consistent trend with our previously obtained values for the Pd_1.00H_x hydride. The elastic constants trend was as expected, with the bulk modulus decreasing with increasing H concentration. As with Pd_1.00H_x, dynamic stability testing for the Pd_1-y-zAg_yCu_zH_x quaternary structures also predicted H atoms transferring from higher energy TE sites to lower energy OC sites. Our EAM potentials also captured the existence of a miscibility gap for the Pd_1-y-zAg_yCu_zH_x and predicted it to narrow and disappear when Ag and Cu concentration increases as was predicted by the experimental findings.

Additional Materials

The following additional materials are uploaded at the page of this paper:

1. Table S1: PdAgH and PdCuH ab initio data, fitting results, and MD results.

Author Contributions

Robert Fuller and Iyad Hijazi contributed to the work on the PdAgH EAM potentials and related calculations. Chaonan Zhang and Iyad Hijazi contributed to the work on the PdCuH and PdAgCuH EAM potentials and related calculations.

Funding

West Virginia Space Grant Consortium (Grant No. 217164, Funder ID. 10.13039/100005768).

Competing Interests

The authors have declared that no competing interests exist.