First principles modeling of the structural, electronic, and vibrational properties of Ni$_{40}$Pd$_{40}$P$_{20}$ bulk metallic glass

The structural, vibrational, and electronic properties of Ni$_{40}$Pd$_{40}$P$_{20}$ bulk metallic glass have been studied using ${\it ab\,initio}$ molecular-dynamics simulations and total-energy optimization. Structural analyses of the resulting ${\it ab\,initio}$ models show the presence of few to no P-P bonds and two main building blocks, consisting of tricapped trigonal prism (TTP) and capped square anti-prism (CSAP) with P as the center of these blocks. The computed Pd and Ni K-edge spectra of extended x-ray absorption fine structure (EXAFS) are found to be in good agreement with experimental data. The configurational-average static structure factor and the generalized vibrational density of states are also observed to be in good agreement with experimental data.


I. INTRODUCTION
Metallic glasses (MG) are a class of amorphous metallic alloys that are produced by rapid quenching from a molten state at a sufficiently fast cooling rate so that crystal nucleation and growth can be avoided. The ability of a metal alloy to form a glass -the so called glass forming ability (GFA)depends on the critical cooling rate, which is defined as the slowest cooling rate at which the melt can be quenched into a glassy state. The smaller the critical cooling rate, the higher the GFA of a metallic alloy. 1 Another measure for the GFA of an alloy is the Turnbull criterion. 2,3 The latter predicts that glass formation should occur if the reduced glass transition temperature, T rg , satisfies the condition T rg > 2/3, where T rg is defined as T rg = T g /T m , with T g and T m being the glass transition and melting temperatures of the alloy, respectively. The formation of metallic glasses is a difficult process in comparison with non-metallic glasses, such as silicates, ceramic, and polymeric glasses, because pure metals and ordinary metallic alloys readily crystallize when quenched from a molten state. The first metallic glass, Au 75 Si 25 , was reported to be obtained by rapid quenching from the liquid state at a very high cooling rate of 10 5 -10 6 K/s, with a critical thickness of ∼ 50 µm. 4 In general, metallic glasses formed via very high cooling rates (10 5 -10 10 K/s) have thicknesses in the sub-millimeter range; such metallic glasses are known as conventional metallic glasses (MG). By contrast, a bulk metallic glass (BMG) is defined to be one with a critical thickness of at least a millimeter -the word "bulk" being indicative of the millimeter length scale. 5 The first BMG was produced from a Pd-Cu-Si alloy by quenching from the melt at a relatively lower cooling rate of 10 3 K/s. In comparison to their crystalline metallic counterparts, BMGs possess unique properties such as higher strength, lower Young's modulus, improved wear resistance, good fatigue endurance, and excellent corrosion resistance, due mainly to the combination of metallic bonding and amorphous structure. 6 As a result, they are extensively employed in biomaterials research and they have potential uses for improved corrosion resistance, biocompatibility, strength, and longevity of biomedical implants, which make BMGs as promising candidates for cardiovascular applications. 1,5,[7][8][9] Bulk metallic glasses also find applications in the manufacture of electromagnetic instruments, cell-phone cases, optical mirrors, striking face plate in golf clubs, connecting parts for optical fibers, soft magnetic high-frequency power coils, etc. 10 Currently, a wide range of multi-component BMG alloys 5,11 exists; however, the industrial production of BMGs remains a difficult problem due to the requirement of high cooling rates to produce a metallic glass with a thickness on the macroscopic length scale.
Ni 40 Pd 40 P 20 was first prepared by Turnbull and coworkers, 12,13 who employed a critical cooling rate of about 1.4 K/s and the application of B 2 O 3 fluxing method to purify the melt and avoid heterogeneous nucleation. The study indicates that a value of T rg close to 2/3 was achieved in their experiments. In a later work, He et al. 14 produced Ni 40 Pd 40 P 20 BMG using a critical cooling rate of 10 −3 K/s. The computational modeling of BMGs presents a significant challenge due to high atomic coordination and the presence of strong chemical order in the system. This is further compounded by the heavy constituent atoms of BMGs, often involving d and f electrons, which exacerbate the problem by increasing the computational complexity of density-functional calculations. Furthermore, the quantum-mechanical computations of the energies and forces are hampered by the existence of finite electronic density of states at the Fermi level, which requires extra care for computation of the density matrix and a sophisticated mixing scheme for self-consistent field calculations of energies and forces.
In this paper, we address the structural, electronic, and vibrational properties of Ni 40 Pd 40 P 20 bulk metallic glass using ab initio molecular-dynamics simulations based on densityfunctional theory by employing a bigger system size and relatively longer cooling rates. Thus far, only a few ab initio studies on the electronic and structural properties of Ni 40 Pd 40 P 20 have been reported in literature. [15][16][17][18] The scope of the present study goes much further by addressing the vibrational properties of glassy Ni 40 Pd 40 P 20 , starting with a notably larger model and a relatively slower cooling rate than the ones that were used in earlier studies. This was addressed by examining the phonon/vibrational density of states of the constituent atoms and it was further validated by making comparisons with available experimental data. Since the vibrational prop-erties of a solid are more sensitive to the local structure and composition than their electronic counterpart, it is of crucial importance in the glass-structure determination. The degree of distortion of the local structural motif(s) can readily manifest in the vibrational excitation spectra, given the miniscule amount of energy (with a few tens of meV) needed to excite vibrational motion. More importantly, unlike earlier studies where the results are solely based on a single atomic configuration, the present study is the first of its kind where the effects of configuration averaging on the physical observables have been discussed at length. Given the importance of configuration fluctuations in small multi-component glassy/disordered systems, which are characterized by several distinctly different bonding environments, an averaging over independent configurations is essential to yield statistically reproducible results on structural, electronic, and vibrational properties of glassy systems. To this end, the results presented here were averaged over ten (10) independent structural configurations to minimize the effects of configuration fluctuations.
The layout of the paper is as follows. In sec. II, we present the computational methodology that involves density functional theory based ab initio molecular-dynamics simulations, followed by total-energy relaxation. The results and associated discussions from our simulations are presented in sec. III, with an emphasis on the electronic, structural, and vibrational properties of the glass. Finally the summary and conclusions of our work are presented in sec. IV.

II. COMPUTATIONAL METHODOLOGY
Ab initio molecular-dynamics (AIMD) simulations and subsequent geometry relaxations were performed using the density functional theory code SIESTA. 19 SIESTA employs local basis sets, and the electron-ion interactions were described using norm-conserving Troullier-Martins pseudopotentials, 20 factorized in the Kleinman-Bylander form. 21 The Perdew-Zunger parameterization 22,23 of the local density approximation (LDA) to the exchange-correlation energy was employed. Single-ζ (SZ) basis functions were used for the finite temperature molecular-dynamics simulation, while double-ζ with polarization (DZP) basis functions were used for geometry relaxations (i.e., total-energy minimizations). Only the Γ-point, k=0, was used to perform the Brillouin zone integration.
The Ni 40 Pd 40 P 20 simulations commenced with the generation of ten (10) independent initial random configurations, each comprising N = 300 atoms (120 Ni atoms, 120 Pd atoms, and 60 P atoms) placed in a cubic supercell at the experimental density 24 of 9.4 g/cm 3 . Each configuration was subjected to independent NVT runs by integrating the equations of motion using the velocity-Verlet algorithm, with a time step of ∆t = 1 fs. The temperature of the system was controlled by a chain of Nosé-Hoover thermostats. [25][26][27] The AIMD simulations commenced by equilibrating each system at 2400 K for 20 ps. After equilibration, each system was quenched to 300 K over a period of 200 ps; this corresponds to an average cooling rate of 10.5 K/ps (or 1.05×10 11 K/s). The geometry at the end of the 300 K dynamics was optimized by minimizing the total energy using a conjugate-gradient algorithm. The final relaxed geometry was used for structural analysis. The convergence criterion for the geometry optimization was set to a maximum force of 0.01 eV/Å on each atom.
The vibrational density of states of each of the ten (10) geometrically-relaxed configurations was computed within the harmonic approximation. Small atomic displacements of 0.005Å along six coordinate directions (±x, ±y, ±z) were applied to compute the total force on each atom. This corresponds to a total of 6N force computations, where N is the number of atoms. The dynamical matrix was constructed from the numerical derivatives of the forces with respect to the atomic displacement and diagonalized to obtain the vibrational eigenfrequencies and eigenvectors. The procedure was repeated for each of ten (10) configurations. All the properties presented in the sections to follow were averaged over the ten (10) independent configurations.

III. RESULTS AND DISCUSSION
A. Structural Properties  Table II). The configurationally averaged total and partial pair-correlation functions (PCF) are shown in Fig. 2.
While the total PCF exhibits a first minimum near 3.5Å, the partial PCFs show varying first shell minima, from 2.3Å for Ni-P to 3.8Å for Pd-P and P-P. A direct implication of this variation is that there are virtually no or very little P-P bonds in the network. This is evident from the pair-correlation data in Fig. 2(b); the P-P correlation begins from just under 3.0Å but a significant number of P atoms form a bond with Ni atoms within a distance of 3.0Å. Table I lists the average coordination number, N , of each atom, as well as the average number N Ni , N Pd , and N P of Ni, Pd, and P, respectively. An examination of Figs. 2(b) and 2(c) suggests that Pd and P has the highest and lowest coordination numbers, respectively. Of all the different bonds formed by Pd, the Pd-Pd bond is the most dominant, followed by Pd-Ni, and Pd-P. In the presence of virtually no P-P bonds, P atoms are segregated in the network from each other. This observation has been noted by Kumar et al. 16 Figure 3 shows the distribution of the average coordination numbers for each atomic species. A further breakdown of the nuances of the coordination numbers is presented in Table II. From Fig. 3, we see that the coordination numbers for Ni lies in the range from 9 to 15 with 40% of the Ni atoms are 11-fold coordinated and 31% are 12-fold coordinated. A breakdown of 11-fold and 12-fold coordinated atoms for Ni in Table II (top) suggests Ni mostly forms bonds with Pd (45%-47%) -  an observation that agrees with N Pd =5.2 being the largest coordination number for Ni in Table II. Likewise, the results for Pd in Fig. 3 indicates that the coordination numbers for Pd ranges from 10 to 17, with the dominant coordination numbers being 13 and 14. Pd atoms largely form bonds with another Pd atoms. A similar analysis confirms that nearly 50% of P atoms  7 8 9 10 11 12 13 14 15  are 9-fold coordinated, with the remainder being 8 and 10fold coordinated. An examination of these 8-, 9-, and 10-fold coordinated P atoms in Table II shows that P atoms prefers to bond with Pd and Ni atoms (approximately 50% and 46-49%, respectively). Table II for P atoms show only 0.3% P-P bonds have been realized in the network. In summary, the bonding trend between atom types observed in this study is consistent with the results from previous ab initio studies based on small models. 16,18 Figure 4 presents the atom-centered bond-angle distributions for Ni, Pd, and P. The distributions show a bimodal character with peaks in the vicinity of 60 • ±10 • and 120 • ±10 • . The peak positions for the P-centered bondangle distribution are particularly noteworthy because they bear close resemblance to the angular distribution of water molecules around +3 lanthanide and actinide metal ions. 28 In +3 lanthanide and actinide ions, the geometry of water molecules with an angular distribution similar to the Pcentered bond-angle distribution takes the form of tricapped trigonal prism (TTP) and capped square anti-prism (CSAP), both correspond to a metal coordination of 9. The bond-angle distributions of ideal TTP and CSAP geometries are shown in Fig. 4 as a palisade of δ-functions. It is evident that the peaks in the P-centered bond-angle distribution closely follow the center of the δ-functions from ideal TTP and CSAP structures, indicating the presence of distorted TTP and CSAP structural units in Ni 40 Pd 40 P 20 glass. To examine this further, we performed a search for approximate TTP and CSAP structural units, centered on P, in the network and found multiple  largely attributed to the broadening of the angular distributions in glassy Ni 40 Pd 40 P 20 environment but contributions from finite-temperature effects are also expected to play a part. The existence of such structural units was reported by Kumar et al. 16 in ab initio studies of Ni 40 Pd 40 P 20 glass and by Rainer et al. 29 and in phosphide-based metal clusters calculations. A characteristic feature of TTP and CSAP structures is that they can easily transform from one to other under the influence of temperature. It has been observed that an array of coordina-tion structures exists that are closely related to the TTP and CSAP structures. For example, an 8-fold coordinated atom is usually a TTP or a CSAP structure with a missing capping atom; a 10-fold coordinated atom could be a bicapped square anti-prism (CSAP structure plus an additional capping atom); a 12-fold coordinated atom is a quadruple-capped square anti-prism (CSAP structure plus three additional capping atoms). TTP, CSAP, and their low and high coordination geometric derivatives are expected be the building blocks of on the nanometer-scale structural order in this BMG. However, given the size of the supercell of about 8Å containing only 300 atoms, it is difficult for the system to form structural order on the nanometer length scale. Figure 6 depicts the simulated static structure factor, averaged over 10 configurations, along with the experimental structure-factor data from a recent differential scanning calorimetry (DSC) experiments, due to Lan et al. 30 The static structure factor, S(q), can be obtained directly from the position of the atoms in the network: where r jk = r k − r j is the vector from the position r j of the j-th atom to the position r k of the k-th atom (with periodic boundary conditions taken into account to mimic bulk effects), f j is the atomic form factor of the j-th atom, q is wave-vector transfer of magnitude q and symbol . indicates the rotational averaging of q of magnitude q over a solid angle of 4π 31 and averaging over all 10 configurations. The atomic form factor was computed using the Gaussian representation: where a j , b j , and c are atom-type specific empirical constants taken from Ref. 32. The match between the simulated structure factor with experiment in Fig. 6 is very good, further validating the accuracy of the structural models. The simulated data is also in good agreement with the experimental data by Egami et al. 33

B. Comparison of Simulated and Experimental XAFS Spectra
To further verify the credibility of the models presented in this work, we compare a key structural property, namely the extended X-ray absorption fine-structure spectroscopy (EX-AFS), of our simulated atomistic configurations to available experimental data. For a given configuration, 6Å clusters centered on the absorber (that is the Pd or K atom) are carved out used as input to the FEFF9 ab initio multiple scattering code, 34,35 which produces an K edge EXAFS spectrum for each cluster. The separate spectra for the individual clusters  6. (Color online) The static structure factors of Ni40Pd40P20 from AIMD simulations at 300 K and experiments. Experimental data are from differential scanning calorimetry measurements at 313 K due to Lan et al.. 30 are then averaged to yield EXAFS spectrum for that configuration. The EXAFS spectra for the ten (10) atomistic configurations was then averaged to obtain to representative EXAFS spectrum χ(k), of Pd and Ni K-edges. Depicted in Figures are the k 2 χ(k) EXAFS spectra and its corresponding Fourier transformχ(R) given by: where Ω(k) is the Hanning window function. In Figs. 7(a) and 7(b), the simulated K-edge k 2 χ(k) EX-AFS and the magnitude of its Fourier transform, |χ(R)| for Pd are depicted and compared recent experimental data due to Kumar et al. 16 In Fig. 7(a), we observe that the frequency of oscillations and amplitude of the simulated spectra matches fairly well with the experimental spectra. The Fourier-transformed data is displayed in the bottom Figure. The Figure does not include the correction for the phaseshift associated with the absorber-scatterer interactions and hence the displayed radial distances are lower than the true radial distances in the partial RDF for Pd. The plots for Ni K-edge in Fig. 7(c) shows amplitudes in the simulated data which match fairly well with experimental data. The Fouriertransformed data for the Ni K-edge in Fig. 7(d) also correlates well with the experimental data.

C. Electronic properties
Depicted in Fig. 8 are the total and angular-momentum resolved (for each atom type) electronic density of states (DOS) of the single particle Kohn-Sham energy eigenvalues. The Fermi level was set to E=0 (dashed vertical line). The region between E = −15 eV and E = −11.5 eV is dominated by the P 3s states. The P 3s hybridizes with the Ni 4s and Pd 5s electron states corresponding to the Ni and Pd atoms located at the vertices of the TTP and CSAP geometries. 15 The region between E = −9 eV and E = −5 eV is dominated by the P 3p states, with a small admixture by the Ni and Pd s and d states; these bands are responsible for the capping bonds between P and the capping Ni and Pd atoms in the TTP and CSAP structures. The strong hybridization between the P electrons states and the Ni and Pd states suggest the formation of P-Ni and P-Pd covalent bonds. The region between E = −5 eV and the Fermi level is dominated by the Ni 3d and Pd 4d bands, with a small contribution from the P 3d band. This energy region is responsible for the Pd-Pd, Pd-Ni, and Ni-Ni metallic bonding. Overall, the DOS analysis indicate that the P-Ni and P-Pd bonds are much stronger than metallic bonds, and hence the presence of P makes the BMG more cohesive.

D. Vibrational properties
The vibrational analysis of our models was based on the vibrational eigenfrequencies and eigenvectors obtained by diagonalizing the dynamical matrix. Specifically, we computed the atom-projected vibrational density of states (VDOS), total VDOS, generalized VDOS, and the inverse participation ratio (IPR).
The vibrational density of states, g α (ω), projected on atoms of type α at frequency ω is given by where α=Ni, Pd or P, {ω i : i = 1, 3N } are the phonon eigenfrequencies, N α is the total number of atoms of type α, and |χ ij | 2 is the contribution of atom j to the square of the amplitude of the normalized phonon eigenvector corresponding to the i th phonon eigenfrequency ω i . From the definition of g α (ω), the total vibrational density of states, g(ω), inverse participation ratio, IPR(ω), and generalized vibrational density of states (GVDOS) 36 , G(ω), at frequency ω are given by: where N is the total number of atoms, exp(−2W α ) is the Debye-Waller factor for the atom of type α, c α is the concen-tration of atoms of type α in the system, σ α is the total scattering cross-section of atom α and M α is the atomic mass of atom of type α. For this work, the weight factors w α (in equations 7 and 8) from the experimental work of Suck 36 , namely w Ni = 0.768, w Pd = 0.102, w P = 0.13, were employed. In Fig. 9(a) the atom-resolved and total vibrational density of states are shown. Clearly the low frequency states are dominated by Ni and Pd while the high-frequency states originate from the atomic vibrations of Ni and P. The IPR shown in Fig. 9(b) is measure of the amount eigenmode amplitude localized at the atomic sites. For a perfectly localized eigenmode conjugate to the eigenfrequency ω i , only a single atom gives rise to the vibrational amplitude and hence IPR(ω i )=1. For a perfectly delocalized eignemode corresponding to the eigenfrequency ω i , each atomic site makes an equal contribution to the amplitude of the eigenmode, and so IPR(ω i )=1/N . From the IPR plot in Fig. 9(b), an extended state-to-localized phonon eigenstate transition in the vicinity of ω = 300 cm −1 can be observed. An analysis of the atomic contribution to the amplitudes of vibration of the phonon eigenmodes indicates that the amplitudes of vibrations of majority of the modes below ω = 300 cm −1 are delocalized over all the atoms, while the amplitudes of the high-frequency eigenvectors were localized on P and about four or five neighboring Pd and/or Ni atoms. As an example, Fig.10 depicts the magnitude and direction of an atom's contribution to localized eigenmode amplitude (IPR=0.66, ω = 485 cm −1 ). The central P atom, shown in black color, makes an 80% contribution to the amplitude, while two Ni atoms (red) and one Pd atom (blue) make contributions of approximately 6% (resulting in a net contribution of approximately 98% by the four atoms). More importantly, the Ni and Pd atoms are within a radius of less than 2.3 A from the P atom. The nature of the mode in Fig. 10 appears to be an admixture of stretching and bending modes, centered on P, with the stretching modes being dominant. The other localized modes were found to behave in a fashion similar to the one in Fig. 10. The small localization radius (centered on P) is a characteristic feature of all the localized eigenmodes. In Fig. 11, the GVDOS is computed and compared with experimental data. Ni vibrations dominate the lower energy of the GVDOS, while P vibrations dominate the higher energy bands. Apart from the small deviations of the simulated GVDOS from the experimentally measured data, the overall match is good, further validating the reliability of the applied theoretical model.

IV. CONCLUSIONS
The structural, electronic, and vibrational properties of Ni 40 Pd 40 P 20 bulk metallic glass have been studied using ab initio molecular-dynamics simulations and total-energy minimizations. Since multi-component BMGs are characterized by a variety of atomic environments, an ensemble of ten (10) independent 300-atom structural configurations were simulated and studied in an effort to take into account the effect of (Color online) An example of a localized mode in Ni40Pd40P20 with IPR = 0.66 and ω = 485 cm −1 . The red arrow represents an 80% contribution by P (black) to the vibrational amplitude while the yellow arrows represent approximately 6% contribution each from two Ni (red) atoms and one Pd (blue) atom. The average bond length between the central P atom and the neighboring three atoms is about 2.2Å, which is indicative of the vibrational localization length of the mode. structure factor obtained from our simulations is found to be in excellent agreement with experimental data. Likewise, the simulated Ni and Pd K-edge EXAFS spectra, averaged over 10 independent configurations, are also observed to be in good agreement with experimental measurements. The structural analysis of the resulting 10 configurations shows the absence of P-P bonds, with the segregated P atoms acting as the centers for distorted tricapped trigonal prisms and capped square anti-prism geometries. An analysis of the electronic density of states indicates that P atoms form strong covalent bonds with Ni and P at low energies, while the bonding states near the (Color online) Comparison of the calculated generalized vibrational densities of states (GVDOS) with inelastic neutronscattering experimental data, due to Suck. 36 indicates that the low-frequency vibrational modes are dominated by the vibrations of Ni and Pd atoms, while the highfrequency modes primarily arise from the atomic vibrations of P atoms. The amplitudes of the high-frequency eigenmodes are observed to be well-localized on P atoms, which comprise an admixture of P-Ni and P-Pd stretch modes and bending modes. The generalized vibrational density of states is also found to agree fairly closely with the experimental data from inelastic neutron-scattering measurements.
We conclude this section with the following observation. A question of considerable importance that the current study could not address is the role of medium-range order on the nanometer length scale, which has been postulated as one of the qualities that give rise to the unique properties of BMGs. Although the presence of a few structural motifs in the glassy Ni 40 Pd 40 P 20 environment is indicative of the medium-range order, the limited size of the models does not permit us to address this issue directly from our simulations. Nonetheless, the approach presented here can be employed profitably to simulate larger models on the nanometer length scale by using small high-quality structural configurations as the basic building blocks, which has the correct short-range order and the local chemistry built into it. In future, we shall address this building-block approach by combining the small structural units developed here to form a super structure/unit and annealing the resultant the super structure through a series of temperatures. Such an approach would be computationally efficient and the resulting models should provide an accurate structural basis to address the nanoscale properties of bulk metallic glasses.