Quaternary Hydrides Pd 1- y-z Ag y Cu z H x Embedded Atom Method Potentials for Hydrogen Energy Applications

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 z H x . The fully analytical potentials are fitted utilizing the central atom method without performing time-consuming molecular dynamics 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 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.75 Ag 0.25 H x and Pd 0.75 Cu 0.25 H x structures were consistent with the ab initio fitting data for most of the H concentrations. The MD results for the Pd 1- y Ag y H x , Pd 1- y Cu y H x and Pd 1- y-z Ag y Cu z H x structures also demonstrated a consistent trend with our previously obtained values for the Pd 1.00 H x hydride. The elastic constants trend was as expected, with the bulk modulus decreasing with increasing H concentration. As with Pd 1.00 H x , dynamic stability testing for the Pd 1- y-z Ag y Cu z H 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-z Ag y Cu z H x and predicted it to narrow and disappear when Ag and Cu concentration increases as was predicted by the experimental findings.

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.

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, FPd, FAg, FCu and FH, 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 Etot of an EAM system is given by [41].
where ρi is the electron density for atom i, Fi is the embedding energy, ϕij is the pair potential between atom i and atom j, rij represents the distance from atom i and to atom j, fj 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: The host electron density is given by: where fe is a scaling factor that can be obtained from fe = Ec/Ω, Ω is the atomic volume and Ec is the cohesive energy, re 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: 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: where D, α, β, and r0 are fitting parameters, r0 defines the interatomic spacing between two atoms, and D (β − α) is the binding energy. The H electron density function is given by: which has 2 fitting parameters CH and δH, while the embedding function for H has the form: where aH, bH, cH, dH 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: 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].

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 ao, Ec, C11, C12, C44, and Evf. Where ao is the lattice constant, Ec is the cohesive energy, C11, C12, and C44 are three elastic constants, and Evf is the vacancy formation energy. Table 1 below lists the Pd, Ag and Cu fitting parameters. 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. 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: For each type in the alloy, the electron density parameter can be calculated from the equation fe=S(Ec/Ω), 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: where ̅ is the expression of electron density for type-a atom, r a is the distance to its closest neighbor and Er 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, Sa and Sb, 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 Er, 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]. 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. For the PdxCu1-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 PdxCu1-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 [54,5,28]. The bulk modulus for Pd1-xAgx and Pd1-xCux 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 C11, C12 and bulk modulus for the pure Pd metal, as can be seen at the left edges of Figure 3 [28,56].
In Figure 4, the values for the elastic constants C44 and C' from MD simulations for Pd1-xAgx 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 Pd1-xAgx potential underestimates C44 relative to the DFT results and have an overall decreasing trend. For Pd1-xCux, our C44 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].

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 Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx 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 [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, PdH0.50, PdH1.00, Pd0.75Ag0.25 and Pd0.75Cu0.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 Pd0.75Ag0.25Hx structures, which have been applied in a similar manner to the Pd0.75Cu0.25Hx 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 Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx structures are given in Table 4 and the shifting factors in Table 5.    [5] and have similar trend, as can be seen in Figure 6. Figure 6 also shows that the Pd0.75Cu0.25Hx shifted data have a similar trend to the Pd0.75Ag0.25Hx shifted data, but with lower cohesive energy values. 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: 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]. 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 Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx 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).

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 Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx 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.

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 Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx 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 Pd1.00Hx 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 Pd0.75Ag0.25Hx [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 Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx 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 Pd1.00Hx plot [40]. For the Pd0.75Ag0.25Hx 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].

Elastic Constants and Bulk Modulus
Using the relaxed Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx 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 C11 and C12 and bulk modulus for Pd0.75Ag0.25Hx, Pd0.75Cu0.25Hx and Pd1.00Hx 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 Pd0.75Ag0.25Hx structures and in a good agreement with our previously obtained results for Pd1.00Hx [5,40]. The bulk modulus values from our MD simulations for Pd0.75Cu0.25Hx 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 Pd0.75Ag0.25Hx 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 Pd0.75Cu0.25Hx compositions, our simulation results for C11 and C12 also have a similar smooth overall decreasing trend to those obtained by Wei et al. [58].  Figure 11 shows the elastic constants C44 and C' for Pd1.00Hx, Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx alloys. As previously obtained for Pd1.00Hx [40], the plot values for C44 for Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx decrease while the shear elastic constants C' increase with increasing H composition. In addition, our elastic constants values for the various Pd1.00Hx, Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx 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: C11 -C12 > 0, C11 + 2C12 > 0 and C44 > 0. From Figure 10 and Figure 11, it can also be seen that the Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx structures have smaller C12 and bigger C44 than the Pd1.00Hx 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].

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 Pd0.50Ag0.50Hx, Pd0.50Cu0.50Hx and Pd0.50Ag0.25Cu0.25Hx hydrides. The Pd50Cu50 structure was chosen on purpose because of its relative similarity to the Pd52.5Cu47.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 Pd1.00Hx, Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx hydrides.

Dynamic Stability
In the Pd1.00Hx 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 Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx 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 Pd1-yAgyHx, Pd1-yCuyHx and Pd1-y-zAgyCuzHx 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 Pd0.50Ag0.25Cu0.25Hx is shown in Figure 14

Miscibility Gap and Gibbs Free Energy of Mixing
For the previously studied Pd1.00Hx 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 Pd1-yAgyHx, Pd1-yCuyHx and Pd1-y-zAgyCuzHx hydrides as follows: where the cohesive energies Pd 1− , Pd 1− 1.0 , and Pd 1− were applied, and X is the mole fraction which X= x/(1+x).
where kB 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 Pd1.00Hx 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 Pd1.00Hx 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 Pd1-yAgyHx at 300 K show that when Ag concentration increases, the values become more negative relative to the Pd1.00Hx 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 Pd1-yAgy [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 Pd1-yCuyHx 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 Pd1.00Hx 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 Pd0.50Ag0.25Cu0.25Hx compositions, the Gibbs free energy plot has a similar trend with those obtained from the Pd1.00Hx 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.

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 Pd0.75Ag0.25Hx and 14 Pd0.75Cu0.25Hx structures, obtained from ab initio SIESTA simulations. The MD simulations utilizing LAMMPS demonstrated that our lattice constant and cohesive energy results for Pd0.75Ag0.25Hx and Pd0.75Cu0.25Hx structures were consistent with the ab initio fitting data for most of the H concentrations. The MD results for the Pd1-yAgyHx, Pd1-yCuyHx and Pd1-y-zAgyCuzHx structures also demonstrated a consistent trend with our previously obtained values for the Pd1.00Hx hydride. The elastic constants trend was as expected, with the bulk modulus decreasing with increasing H concentration. As with Pd1.00Hx, dynamic stability testing for the Pd1-y-zAgyCuzHx 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 Pd1-y-zAgyCuzHx 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.