_{1-y-z}Ag

_{y}Cu

_{z}H

_{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.

**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):25; 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.

**Abstract**

The Pd-H system has attracted extensive attention. Pd can absorb considerable amount of H at room temperature, this ability is reversible, so it is suitable for multiple energy applications. Pd-Ag alloys possess higher H permeability, solubility and narrower miscibility gap with better mechanical properties than pure Pd, but sulfur poisoning remains an issue. Pd-Cu alloys have excellent resistance to sulfur and carbon monoxide poisoning and hydrogen embrittlement, good mechanical properties, and broader temperature working environments over pure Pd, but relatively lower hydrogen permeability and solubility than pure Pd and Pd-Ag alloys. This suggests that alloying Pd with Ag and Cu to create Pd-Ag-Cu ternary alloys can optimize the overall performance and substantially lowers the cost. Thus, in this paper, we provide the first embedded atom method potentials for the quaternary hydrides Pd_{1-y-z}Ag* _{y}*Cu

*H*

_{z}*. The fully analytical potentials are fitted utilizing the central atom method without performing time-consuming molecular dynamics simulations.*

_{x}**Keywords **

EAM; DFT; metal hydrides; hydrogen storage; miscibility Gap; Pd-Ag-Cu-H

**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_{2}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_{77}Ag_{23 }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_{2}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*and

_{Cu}*F*, 4 electron density functions,

_{H}*ρ*

*,*

_{Pd}*ρ*

*,*

_{Ag}*ρ*

*and*

_{Cu}*ρ*

*, and 10 pair functions,*

_{H}*φ*

*,*

_{Pd–Pd}*φ*

*,*

_{Ag–Ag}*φ*

*,*

_{Cu–Cu}*φ*

*,*

_{H–H}*φ*

*,*

_{Pd–Ag}*φ*

*,*

_{Pd–Cu}*φ*

*,*

_{Ag–Cu}*φ*

_{Pd–H}_{, }

*φ*

*and*

_{Ag–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*

_{Cu–H}*E*of an EAM system is given by [

_{tot}**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*is the embedding energy,

_{i}*ϕ*is the pair potential between atom

_{ij}*i*and atom

*j*,

*r*represents the distance from atom

_{ij}*i*and to atom

*j*,

*f*is the electron density function of distance from the center of atom

_{j}*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*, Ω is the atomic volume and

_{c}/Ω*E*is the cohesive energy,

_{c}*r*is the equilibrium closest distance, and χ is a fitting parameter. The pair potential function is the modified potential created by Rose et al. [

_{e}**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 _{0}* are fitting parameters,

*r*

_{0 }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

*δ*, while the embedding function for H has the form:

_{H}$$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*are fitting parameters, and

_{H }*ε*= 0.0540638.

_{H}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 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*,

_{11}*C*, C

_{12}_{44}, and

*E*. Where

_{vf}*a*is the lattice constant,

_{o}*E*is the cohesive energy,

_{c}*C*,

_{11}*C*, and C

_{12}

_{44}_{ }are three elastic constants, and

*E*is the vacancy formation energy. Table 1 below lists the Pd, Ag and Cu fitting parameters.

_{vf}**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*, 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

_{b}*, 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. [*

_{r}**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* _{x}*Ag

_{1-x}and Pd

*Cu*

_{x}_{1-x}(0 ≤

*x*≤ 1) structures are almost identical with the experimental data [

**51**]. For the Pd

*Ag*

_{x}_{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* _{x}*Cu

_{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

*Cu*

_{x}_{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-}* _{x}*Ag

*and Pd*

_{x}*Cu*

_{x}_{1-x}structures

*a*and

_{o}*E*results obtained from MD simulations, experimental and DFT calculations [

_{c}**28**,

**55**,

**52**,

**51**].

In Figure 3(a) and (b) the values for the elastic constants *C _{11}*

_{, }

*C*from MD simulations for Pd

_{12}_{1-}

*Ag*

_{x}*and Pd*

_{x}_{1-}

*Cu*

_{x}*structures show consistent trend with the DFT calculations and the results from Hale et al. and Kart et al. [*

_{x}**54**,

**5**,

**28**]. The bulk modulus for Pd

_{1-}

*Ag*

_{x}*and Pd*

_{x}_{1-}

*Cu*

_{x}*structures obtained from MD simulations match the softening trends predicted by the DFT calculations as well [*

_{x}**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).

_{12}**Figure 3 **Pd_{1-}* _{x}*Ag

*and Pd*

_{x}*Cu*

_{x}_{1-x}alloys

*C*,

_{11}*C*elastic constant and B

_{12}_{m}results obtained from MD simulations, experimental and DFT calculations [

**28**,

**56**].

In Figure 4, the values for the elastic constants *C _{44}*

_{ }and

_{ }

*C’*from MD simulations for Pd

_{ }_{1-}

*Ag*

_{x}*show that our results are closer to the experimental data at the edge of the composition range than those of Hale et al. [*

_{x}**5**]. As with the Hale et al. EAM Morse model [

**5**], our Pd

_{1-}

*Ag*

_{x}*potential underestimates*

_{x}*C*relative to the DFT results and have an overall decreasing trend. For Pd

_{44}_{1-}

*Cu*

_{x}*, our*

_{x}*C*and

_{44}*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-x}Ag* _{x}* and Pd

*Cu*

_{x}_{1-x}alloys

*C*and

_{44}*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*

*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

*φ*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 [

_{Cu–H}**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.75}Ag_{0.25}H_{x}_{ }and Pd_{0.75}Cu_{0.25}H* _{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-}* _{x}*Ag

*DFT results from Figure 2(a) and (b) [*

_{x}**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.75}Ag

_{0.25}and Pd

_{0.75}Cu

_{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.75}Ag

_{0.25}H

*structures, which have been applied in a similar manner to the Pd*

_{x}_{0.75}Cu

_{0.25}H

*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*

_{x}_{0.75}Ag

_{0.25}H

*and Pd*

_{x}_{0.75}Cu

_{0.25}H

*structures are given in Table 4 and the shifting factors in Table 5.*

_{x}$$\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.75}Ag_{0.25}H* _{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.75}Cu

_{0.25}H

*shifted data have a similar trend to the Pd*

_{x}_{0.75}Ag

_{0.25}H

*shifted data, but with lower cohesive energy values.*

_{x}**Figure 6** Pd_{0.75}Ag_{0.25}H* _{x}* and Pd

_{0.75}Cu

_{0.25}H

*shifted cohesive energies compared with Hale et al. DFT results.*

_{x}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

*φ*pair potential functions during the fitting process. The

_{Cu–H}*φ*and

_{Ag–H}*φ*pair functions take the same generalized Morse potential mathematical form, as used previously in the Pd-H interaction [

_{Cu–H}**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.75}Ag_{0.25}H* _{x}* and Pd

_{0.75}Cu

_{0.25}H

*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).*

_{x}**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.00}H* _{x}*, Pd

_{0.75}Ag

_{0.25}H

*and Pd*

_{x}_{0.75}Cu

_{0.25}H

_{x}*a*and

_{o}*E*results obtained from MD simulations and fitting.

_{c}*3.2 Results*

*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.75}Ag_{0.25}H* _{x}* and Pd

_{0.75}Cu

_{0.25}H

*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.*

_{x}*3.3 Lattice Constant and Cohesive Energy Results*

*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.75}Ag_{0.25}H* _{x} *and Pd

_{0.75}Cu

_{0.25}H

*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*

_{x}_{1.00}H

*results in our previous work [*

_{x}**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.75}Ag

_{0.25}H

*[*

_{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.75}Ag_{0.25}H* _{x}* and Pd

_{0.75}Cu

_{0.25}H

*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*

_{x}_{1.00}H

*plot [*

_{x}**40**]. For the Pd

_{0.75}Ag

_{0.25}H

*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 [*

_{x}**5**].

*3.4 Elastic Constants and Bulk Modulus*

*3.4 Elastic Constants and Bulk Modulus*

Using the relaxed Pd_{0.75}Ag_{0.25}H* _{x}* and Pd

_{0.75}Cu

_{0.25}H

*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 [*

_{x}**50**], the elastic constants and bulk modulus values were estimated. Figure 10 shows that the elastic constants

*C*and

_{11}*C*and bulk modulus for Pd

_{12}_{0.75}Ag

_{0.25}H

*, Pd*

_{x}_{0.75}Cu

_{0.25}H

*and Pd*

_{x}_{1.00}H

*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*

_{x}_{0.75}Ag

_{0.25}H

*structures and in a good agreement with our previously obtained results for Pd*

_{x}_{1.00}H

*[*

_{x}**5**,

**40**]. The bulk modulus values from our MD simulations for Pd

_{0.75}Cu

_{0.25}H

*also have a similar decreasing trend with those obtained by Wei et al. DFT simulations for fcc Pd-Cu-H (O1) [*

_{x}**58**]. Other researchers have also documented this softening behavior with increasing H concentrations [

**61**,

**62**,

**63**]. The results for Pd

_{0.75}Ag

_{0.25}H

*structures obtained from Hale et al. Hybrid potential yielded an unstable trend but still had an overall similar softening trend [*

_{x}**5**,

**53**]. As for the Pd

_{0.75}Cu

_{0.25}H

*compositions, our simulation results for*

_{x}*C*and

_{11}*C*also have a similar smooth overall decreasing trend to those obtained by Wei et al. [

_{12}**58**].

**Figure 10 **Pd_{1.00}H_{x}, Pd_{0.75}Ag_{0.25}H_{x} and Pd_{0.75}Cu_{0.25}H_{x}*C _{11}*,

*C*elastic constant and B

_{12}_{m}results from MD simulations.

Figure 11 shows the elastic constants* C _{44}* and

*C’*for Pd

_{1.00}H

*, Pd*

_{x}_{0.75}Ag

_{0.25}H

*and Pd*

_{x}_{0.75}Cu

_{0.25}H

*alloys. As previously obtained for Pd*

_{x }_{1.00}H

*[*

_{x }**40**], the plot values for

*C*for Pd

_{44 }_{0.75}Ag

_{0.25}H

*and Pd*

_{x}_{0.75}Cu

_{0.25}H

*decrease while the shear elastic constants*

_{x}*C’*increase with increasing H composition.

_{ }**Figure 11 **Pd_{1.00}H* _{x}*, Pd

_{0.75}Ag

_{0.25}H

*and Pd*

_{x}_{0.75}Cu

_{0.25}H

*elastic constant and*

_{x}C_{44}*C’*shear elastic constant results from MD simulations.

In addition, our elastic constants values for the various Pd_{1.00}H* _{x}*, Pd

_{0.75}Ag

_{0.25}H

*and Pd*

_{x}_{0.75}Cu

_{0.25}H

*alloys shown in Figure 10 and Figure 11 satisfy the theory of strain energy for cubic structures [*

_{x}**64**]. According to strain energy theory, the following formulas can be applied to a mechanically stable cubic:

*C*-

_{11 }*C*> 0,

_{12}*C*+ 2

_{11 }*C*> 0 and

_{12}*C*> 0. From Figure 10 and Figure 11, it can also be seen that the Pd

_{44 }_{0.75}Ag

_{0.25}H

*and Pd*

_{x}_{0.75}Cu

_{0.25}H

*structures have smaller*

_{x}*C*and bigger

_{12}*C*than the Pd

_{44}_{1.00}H

*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 [*

_{x}**58**].

*3.5 Additional Compositions*

*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.50}Ag_{0.50}H* _{x}*, Pd

_{0.50}Cu

_{0.50}H

*and Pd*

_{x}_{0.50}Ag

_{0.25}Cu

_{0.25}H

*hydrides. The Pd*

_{x}_{50}Cu

_{50}structure was chosen on purpose because of its relative similarity to the Pd

_{52.5}Cu

_{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.00}H

*, Pd*

_{x}_{0.75}Ag

_{0.25}H

*and Pd*

_{x}_{0.75}Cu

_{0.25}H

*hydrides.*

_{x }**Figure 12 **Pd_{1.00}H* _{x}*, Pd

_{0.75}Ag

_{0.25}H

*, Pd*

_{x}_{0.50}Ag

_{0.50}H

*, Pd*

_{x}_{0.75}Cu

_{0.25}H

*, Pd*

_{x}_{0.50}Cu

_{0.50}H

*and Pd*

_{x}_{0.50}Ag

_{0.25}Cu

_{0.25}H

_{x}*a*and

_{o}*E*results from MD simulations.

_{c}In Figure 13(a), (b) and (c) the MD simulation results for the elastic constants and bulk modulus for Pd_{0.50}Ag_{0.50}H_{x}, Pd_{0.50}Cu_{0.50}H* _{x}* and Pd

_{0.50}Ag

_{0.25}Cu

_{0.25}H

*are plotted. It can be seen that they also have a consistent trend with our previous results for the Pd*

_{x}_{1.00}H

*, Pd*

_{x}_{0.75}Ag

_{0.25}H

*and Pd*

_{x}_{0.75}Cu

_{0.25}H

*hydrides. The plots for bulk modulus for Pd*

_{x}_{0.50}Ag

_{0.50}H

*, Pd*

_{x}_{0.50}Cu

_{0.50}H

*and Pd*

_{x}_{0.50}Ag

_{0.25}Cu

_{0.25}H

*have a similar softening trend to our previous results for Pd*

_{x}_{1.00}H

*, Pd*

_{x}_{0.75}Ag

_{0.25}H

*and Pd*

_{x}_{0.75}Cu

_{0.25}H

*. For Pd*

_{x}_{0.50}Ag

_{0.25}Cu

_{0.25}H

*, the bulk modulus values*

_{x}_{ }are higher than those for Pd

_{0.50}Ag

_{0.50}H

*and closer to Pd*

_{x}_{0.50}Cu

_{0.50}H

*, indicating that adding Cu has a strengthening impact on the solid solution while adding Ag has a softening effect.*

_{x}The elastic constant values for *C _{11}* and

*C*for the various Pd

_{12}_{0.50}Ag

_{0.50}H

*, Pd*

_{x}_{0.50}Cu

_{0.50}H

*and Pd*

_{x}_{0.50}Ag

_{0.25}Cu

_{0.25}H

*structures shown in Figure 13 also satisfy the strain energy theory for cubic structures [*

_{x}**58**], implying these structures also possess mechanical stability.

**Figure 13 **Pd_{1.00}H* _{x}*, Pd

_{0.75}Ag

_{0.25}H

*, Pd*

_{x}_{0.50}Ag

_{0.50}H

*, Pd*

_{x}_{0.75}Cu

_{0.25}H

*, Pd*

_{x}_{0.50}Cu

_{0.50}H

*and Pd*

_{x}_{0.50}Ag

_{0.25}Cu

_{0.25}H

_{x}*C*,

_{11}*C*elastic constant and B

_{12}_{m}results from MD simulations.

*3.6 Dynamic Stability*

*3.6 Dynamic Stability*

In the Pd_{1.00}H* _{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.75}Ag

_{0.25}H

*and Pd*

_{x }_{0.75}Cu

_{0.25}H

*structures are highly energetically favorable to H atoms; this behavior was also observed and reported by other researchers [*

_{x}**5**,

**52**]. In order to verify the stability for the Pd

_{1-y}Ag

*H*

_{y}*, Pd*

_{x}_{1-y}Cu

*H*

_{y}*and Pd*

_{x}_{1-y-z}Ag

*Cu*

_{y}*H*

_{z}*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*

_{x}_{0.50}Ag

_{0.25}Cu

_{0.25}H

*is shown in Figure 14(b).*

_{x}**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*

*3.7 Miscibility Gap and Gibbs Free Energy of Mixing*

For the previously studied Pd_{1.00}H* _{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-y}Ag

*H*

_{y}*, Pd*

_{x}_{1-y}Cu

*H*

_{y}

_{x }_{}and Pd

_{1-y-z}Ag

*Cu*

_{y}*H*

_{z}

_{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.00}H* _{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.00}H

*at 300 K are*

_{x}*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-y}Ag

*H*

_{y}*at 300 K show that when Ag concentration increases, the values become more negative relative to the Pd*

_{x}_{1.00}H

*system, indicating more favorable mixing, and the miscibility gap become narrower. No miscibility gap observed when*

_{x}*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-y}Ag

*[*

_{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-y}Cu

*H*

_{y}*compositions, experimental data at 303 K indicated that by increasing Cu concentration, the*

_{x}*α*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.00}H

*structures. At*

_{x}*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.50}Ag

_{0.25}Cu

_{0.25}H

*compositions, the Gibbs free energy plot has a similar trend with those obtained from the Pd*

_{x}_{1.00}H

*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.*

_{x}**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.75}Ag_{0.25}H* _{x}* and 14 Pd

_{0.75}Cu

_{0.25}H

*structures, obtained from ab initio SIESTA simulations. The MD simulations utilizing LAMMPS demonstrated that our lattice constant and cohesive energy results for Pd*

_{x}_{0.75}Ag

_{0.25}H

*and Pd*

_{x}_{0.75}Cu

_{0.25}H

*structures were consistent with the ab initio fitting data for most of the H concentrations. The MD results for the Pd*

_{x}_{1-y}Ag

*H*

_{y}*, Pd*

_{x}_{1-y}Cu

*H*

_{y}*and Pd*

_{x}_{1-y-z}Ag

*Cu*

_{y}*H*

_{z}*structures also demonstrated a consistent trend with our previously obtained values for the Pd*

_{x}_{1.00}H

*hydride. The elastic constants trend was as expected, with the bulk modulus decreasing with increasing H concentration. As with Pd*

_{x}_{1.00}H

*, dynamic stability testing for the Pd*

_{x}_{1-y-z}Ag

*Cu*

_{y}*H*

_{z}*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*

_{x}_{1-y-z}Ag

*Cu*

_{y}*H*

_{z}*and predicted it to narrow and disappear when Ag and Cu concentration increases as was predicted by the experimental findings.*

_{x}**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.

**References**

- Nelin G. A neutron diffraction study of palladium hydride. Phys Status Solidi B. 1971; 45: 527-536. [CrossRef]
- Caputo R, Alavi A. Where do the H atoms reside in PdHx systems? Mol Phys. 2003; 101: 1781-1787. [CrossRef]
- Manchester FD, San-Martin A, Pitre JM. The H-Pd (hydrogen-palladium) system. J Phase Equilib. 1994; 15: 62-83. [CrossRef]
- Zhou XW, Zimmerman JA, Wong BM, Hoyt JJ. An embedded-atom method interatomic potential for Pd–H alloys. J Mater Res. 2008; 23: 704-718. [CrossRef]
- Hale LM, Wong BM, Zimmerman JA, Zhou XW. Atomistic potentials for palladium–silver hydrides. Modell Simul Mater Sci Eng. 2013; 21: 045005. [CrossRef]
- Hsu DK, Leisure RG. Elastic constants of palladium and β-phase palladium hydride between 4 and 300 K. Phys Rev B. 1979; 20: 1339. [CrossRef]
- Dong W, Ledentu V, Sautet P, Eichler A, Hafner J. Hydrogen adsorption on palladium: A comparative theoretical study of different surfaces. Surf Sci. 1998; 411: 123-136. [CrossRef]
- Goltsova MV, Artemenko YA, Zaitsev VI. Kinetics and morphology of the reverse β→ α hydride transformation in thermodynamically open Pd–H system. J Alloys Compd. 1999; 293-295: 379-384. [CrossRef]
- Alefeld G, Völkl J. Hydrogen in Metals II: Application-Oriented Properties. Berlin Heidelberg: Springer-Verlag; 1978.
- Muetterties EL. Transition metal--hydrogen interaction. Transition Met Hydrides. 1971; 4: 11-13.
- Siegel B, Libowitz GG, Mueller WM, Blackledge JP. Metal Hydrides. New York: Academic Press; 1968.
- Fukai Y. The Metal-Hydrogen System: Basic Bulk Properties. Berlin Heidelberg: Springer-Verlag; 1993. [CrossRef]
- Lässer R. Tritium and Helium-3 in Metals. Berlin Heidelberg: Springer-Verlag; 1989. [CrossRef]
- Alefeld G, Völkl J. Hydrogen in Metals I: Basic Properties. Berlin Heidelberg: Springer-Verlag; 1978.
- Povel R, Feucht K, Gelse W, Withalm G. Hydrogen fuel for motorcars. Interdiscip Sci Rev. 1989; 14: 365-373. [CrossRef]
- Knapton AG. Palladium alloys for hydrogen diffusion membranes. Platinum Met Rev. 1977; 21: 44-50.
- Adams BD, Chen A. The role of palladium in a hydrogen economy. Mater Today. 2011; 14: 282-289. [CrossRef]
- Jimenez G, Dillon E, Dahlmeyer J, Garrison T, Garrison T, Darkey S, et al. A comparative assessment of hydrogen embrittlement: Palladium and palladium-silver (25 weight% silver) subjected to hydrogen absorption/desorption cycling. Adv Chem Eng Sci. 2016; 6: 246. [CrossRef]
- Foletto EL, Wirbitzki Da Silveira JV, Jahn SL. Preparation of palladium-silver alloy membranes for hydrogen permeation. Lat Am Appl Res. 2008; 38: 79-84.
- Chen CH, Ma YH. The effect of H2S on the performance of Pd and Pd/Au composite membrane. J Membr Sci. 2010; 362: 535-544. [CrossRef]
- O’Brien CP, Gellman AJ, Morreale BD, Miller JB. The hydrogen permeability of Pd4S. J Membr Sci. 2011; 371: 263-267. [CrossRef]
- Gabitto J, Tsouris C. Modeling sulfur poisoning of palladium membranes used for hydrogen separation. Int J Chem Eng. 2019; 2019: 9825280. [CrossRef]
- Amandusson H, Ekedahl LG, Dannetun H. Hydrogen permeation through surface modified Pd and PdAg membranes. J Membr Sci. 2001; 193: 35-47. [CrossRef]
- Semidey-Flecha L, Sholl DS. Combining density functional theory and cluster expansion methods to predict H2 permeance through Pd-based binary alloy membranes. J Chem Phys. 2008; 128: 144701. [CrossRef]
- Qin L, Jiang C. First-principles based modeling of hydrogen permeation through Pd–Cu alloys. Int J Hydrog Energy. 2012; 37: 12760-12764. [CrossRef]
- Adhikari S, Fernando S. Hydrogen membrane separation techniques. Ind Eng Chem Res. 2006; 45:875-881. [CrossRef]
- Martin MH, Galipaud J, Tranchot A, Roué L, Guay D. Measurements of hydrogen solubility in CuxPd100− x thin films. Electrochim Acta. 2013; 90: 615-622. [CrossRef]
- Kart SÖ, Erbay A, Kılıç H, Cagin T, Tomak M. Molecular dynamics study of Cu-Pd ordered alloys. J Achiev Mater Manuf Eng. 2008; 31: 41-46.
- McKinley DL. Method for hydrogen separation and purification. Washington, DC: U.S. Patent and Trademark Office; 1969; 3,439,474.
- Illas F, López N, Ricart JM, Clotet A, Conesa JC, Fernández-García M. Interaction of CO and NO with PdCu (111) surfaces. J Phys Chem B. 1998; 102: 8017-8023. [CrossRef]
- O’Brien CP, Lee IC. The interaction of CO with PdCu hydrogen separation membranes: An operando infrared spectroscopy study. Catal Today. 2019; 336: 216-222. [CrossRef]
- Guerreiro BH, Martin MH, Roué L, Guay D. Hydrogen permeability of PdCuAu membranes prepared from mechanically-alloyed powders. J Membr Sci. 2016; 509: 68-82. [CrossRef]
- Liu LC, Wang JW, He YH, Gong HR. Solubility, diffusivity, and permeability of hydrogen at PdCu phases. J Membr Sci. 2017; 542: 24-30. [CrossRef]
- Tarditi AM, Braun F, Cornaglia LM. Novel PdAgCu ternary alloy: Hydrogen permeation and surface properties. Appl Surf Sci. 2011; 257: 6626-6635. [CrossRef]
- Semidey-Flecha L, Ling C, Sholl DS. Detailed first-principles models of hydrogen permeation through PdCu-based ternary alloys. J Membr Sci. 2010; 362: 384-392. [CrossRef]
- Ling C, Semidey-Flecha L, Sholl DS. First-principles screening of PdCuAg ternary alloys as H2 purification membranes. J Membr Sci. 2011; 371: 189-196. [CrossRef]
- Yuan L, Goldbach A, Xu H. Real-time monitoring of metal deposition and segregation phenomena during preparation of PdCu membranes. J Membr Sci. 2008; 322: 39-45. [CrossRef]
- Zhao L, Goldbach A, Xu H. Tailoring palladium alloy membranes for hydrogen separation from sulfur contaminated gas streams. J Membr Sci. 2016; 507: 55-62. [CrossRef]
- Foiles SM, Hoyt JJ. Computer simulation of bubble growth in metals due to He. Sandia Natl Lab. 2001; SAND2001-0661. [CrossRef]
- Hijazi I, Zhang Y, Fuller R. A simple embedded atom potential for Pd-H alloys. Mol Simul. 2018; 44: 1371-1379. [CrossRef]
- Daw MS, Baskes MI. Semiempirical, quantum mechanical calculation of hydrogen embrittlement in metals. Phys Rev Lett. 1983; 50: 1285. [CrossRef]
- Hijazi IA, Park YH. Consistent analytic embedded atom potential for face-centered cubic metals and alloys. J Mater Sci Technol. 2009; 25: 835-846.
- Park YH, Hijazi IA. Simple analytic embedded atom potential for FCC materials. Int J Microstruct. Mater Prop. 2011; 6: 378-396. [CrossRef]
- Hijazi IA, Park YH. Structure of pure metallic nanoclusters: Monte carlo simulation and ab initio study. Eur Phys J D. 2010; 59: 215-221. [CrossRef]
- Park YH, Hijazi IA. Critical size of transitional copper clusters for ground state structure determination: Empirical and ab initio study. Mol Simul. 2012; 38:241-247. [CrossRef]
- Hijazi IA. Structural, electronic and magnetic properties of pure metallic and bimetallic nanoclusters: Empirical and density functional studies. Ann Arbor: ProQuest Dissertations; 2010.
- Rose JH, Smith JR, Guinea F, Ferrante J. Universal features of the equation of state of metals. Phys Rev B. 1984; 29: 2963-2969. [CrossRef]
- Foiles SM, Baskes MI, Daw MS. Embedded-atom-method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys. Phys Rev B. 1986; 33: 7983. [CrossRef]
- Johnson RA. Alloy models with the embedded-atom method. Phys Rev B. 1989; 39:12554-12559. [CrossRef]
- Thompson A. In.elastic. Sandia National Laboratory. Available from: https://github.com/lammps/lammps/tree/master/examples/ELASTIC.
- Cai J, Ye YY. Simple analytical embedded-atom-potential model including a long-range force for fcc metals and their alloys. Phys Rev B. 1996; 54:8398-8410. [CrossRef]
- Løvvik OM, Olsen RA. Density functional calculations on hydrogen in palladium–silver alloys. J Alloys Compd. 2002; 330-332: 332-337. [CrossRef]
- Hale LM, Wong BM, Zimmerman JA, Zhou XW. Atomistic potentials for palladium–silver hydrides. Model Simul Mat Sci Eng. 2013; 21: 045005. [CrossRef]
- Delczeg-Czirjak EK, Delczeg L, Ropo M, Kokko K, Punkkinen MP, Johansson B, et al. Ab initio study of the elastic anomalies in Pd-Ag alloys. Phys Rev B. 2009; 79: 085107. [CrossRef]
- Uemiya S, Matsuda T, Kikuchi E. Hydrogen permeable palladium-silver alloy membrane supported on porous ceramics. J Membr Sci. 1991; 56: 315-325. [CrossRef]
- Burch R. On the role of silver atoms in the absorption of hydrogen by palladium-silver alloys. Solid State Commun. 1969; 7: 1313-1317. [CrossRef]
- Karakaya I, Thompson WT. The Ag− Pd (silver-palladium) system. Bull Alloy Phase Diagrams. 1988; 9:237-243. [CrossRef]
- Wei C, Kong FT, Gong HR. Phase stability and elastic property of PdH and PdCuH phases. Int J Hydrog Energy. 2013; 38: 16485-16494. [CrossRef]
- ABINIT. Accessed Dec. 31, 2020. [Online]. Available from: https://www.abinit.org/downloads/psp-links/pseudopotentials.
- Fuller R. Improved embedded atom method potentials for metal hydride systems. Huntington: Marshall University; 2018. Available from: https://mds.marshall.edu/etd/1124.
- Ilawe NV, Zimmerman JA, Wong BM. Breaking badly: DFT-D2 gives sizeable errors for tensile strengths in palladium-hydride solids. J Chem Theory Comput. 2015; 11: 5426-5435. [CrossRef]
- Schwarz RB, Bach HT, Harms U, Tuggle D. Elastic properties of Pd–hydrogen, Pd–deuterium, and Pd–tritium single crystals. Acta Mater. 2005; 53: 569-580. [CrossRef]
- Zhong W, Cai Y, Tomanek D. Computer simulation of hydrogen embrittlement in metals. Nature. 1993; 362: 435-437. [CrossRef]
- Mouhat F, Coudert FX. Necessary and sufficient elastic stability conditions in various crystal systems. Phys Rev B. 2014; 90: 224104. [CrossRef]
- Antonov VE. Phase transformations, crystal and magnetic structures of high-pressure hydrides of d-metals. J Alloys Compd. 2002; 330: 110-116. [CrossRef]
- Flanagan TB, Luo S, Clewley JD. Calorimetric enthalpies for the reaction of H2 with Pd–Cu alloys at 303 K. J Alloys Compd. 2003; 356: 13-16. [CrossRef]
- Huang W, Opalka SM, Wang D, Flanagan TB. Thermodynamic modelling of the Cu–Pd–H system. Calphad. 2007; 31: 315-329. [CrossRef]