A new first principles approach to calculate phonon spectra of disordered alloys
Abstract
The lattice dynamics in substitutional disordered alloys with constituents having large size differences is driven by strong disorder in masses, interatomic force constants and local environments. In this letter, a new firstprinciples approach based on special quasi random structures and itinerant coherent potential approximation to compute the phonon spectra of such alloys is proposed and applied to NiPt alloy. The agreement between our results with the experiments is found to be much better than for previous models of disorder due to an accurate treatment of the interplay of interatomic forces among various pairs of chemical species. This new formalism serves as a potential solution to the longstanding problem of a proper microscopic understanding of lattice dynamical behavior of disordered alloys.
pacs:
63.20.dk, 63.50.Gh, 71.20.BeThe energy dispersion of lattice waves in ordered intermetallics and disordered alloys has recently gained significant importance due to its role in stabilizing a particular structural phase by contributing towards the system’s entropy and in pinpointing the fundamental mechanism behind the behavior of materials, which are important in technological applications. Examples include the connection between martensitic transformations in shapememory alloys and the anomalous behavior of their phonon spectra [1], anomalously low thermalexpansion coefficients in Febased Invar alloys and the unusual behavior of phonon frequencies under external pressure [2], and the dramatic alteration of relative stabilities of different precipitates in AlCu system due to vibrational entropy [3].
Theoretical investigation on lattice dynamics in disordered metallic alloys is one of the challenging areas in the state of the art materials science research. Most of the theoretical studies have been done on the alloy compounds although a wealth of neutronscattering data is available for various disordered alloys with arbitrary compositions. The dearth of theoretical works is due to the lack of a suitable selfconsistent, analytical methodology to carry out the averaging over the random configurations, taking into account the disorder due to the fluctuations in mass, force constant and local environment around an atom. The most widely used method is the singlesite coherent potential approximation (SCPA) [4], which being a singlesite theory, is unable to handle the nonlocal fluctuations. The recently developed itinerant coherent potential approximation (ICPA) [5], is a successful generalization of the SCPA in handling disordered alloys with large mass as well as force constant disorder. While the ICPA stands as a selfconsistent, analytic, sitetranslational invariant tool for computing the phonon spectra in disordered alloys, a missing key component in establishing it as an accurate and reliable formalism, is its integration with a tractable structural model for positional disorder which can closely mimic the fluctuations in the interatomic force constants in a random alloy. Moreover, due to the presence of shortrange order, the phonon spectra of such systems would be affected due to the dominant vibrations of unlike species pairs. Initially, attempts were made in the form of firstprinciples computation of interatomic force constants from selected ordered structures and use them in a random alloy environment [6, 7]. This is not a proper remedy as the force constants are not directly transferable across the environments [8, 9].
In this letter, we present a new firstprinciples based formalism to compute the lattice dynamics of substitutional disordered alloys which incorporates the effects of disorder in mass, force constant and environment. We demonstrate the formalism by computing the phonon dispersion spectra of NiPt alloy. This alloy is chosen as its constituents have large mass differences(Pt being 3 times heavier than Ni), large size differences ( 11) and large force constant differences (PtPt force constants being 55 larger than the NiNi ones) and thus significant effects of all three kinds of disorder on the lattice dynamics can be expected. Also, the neutronscattering experiments [10] show anomalous features in the dispersion curves in the form of resonance modes and splitting of dispersion branches. In previous works, neither SCPA [10] nor ICPA [5] with empirical set of force constants could have good quantitative agreement with the experiments, signifying that the interatomic interactions have not been modeled accurately.
Our formalism has three important steps: first we generate a structural model of substitutional disorder by special quasirandom structure (SQS) method, then the averaging of the force constant tensor is done to recover the original symmetry of the solid solution followed by the ICPA for performing the averaging over configurations.
(i) Structural model: First principles supercell calculations of disordered alloys involve the direct averaging of the properties obtained from a number of disordered configurations of the alloy. Current computational capabilities limit the size of the supercell necessary to describe each configuration, as well as the number of configurations sampled. Zunger et. al. developed a computationally tractable approach to model the disorder through the introduction of SQS method [11], an atom per cell periodic structure designed so that their distinct correlation functions best match the ensembleaveraged correlationfunctions of the random alloy. Here corresponds to the figure defined by the number of of atoms located on its vertices (=2, 3, 4…. are pairs, triangles, tetrahedra etc.) with being the order of neighbor distances separating them (=1, 2….are first,second neighbors etc.). As the properties like the equilibrium volume, local density of states and bond lengths around an atom are influenced by the local environment, SQS forms an adequate approximation. The biggest advantage of the SQS over a conventional supercell to model positional disorder is that the former uses the knowledge of the paircorrelation functions, a key property of the random alloys, to decide the positions of the atoms in the unit cell, instead of inserting them randomly as is done in the conventional supercell technique, and thus is guaranteed to provide a better description of the environments in an actual random alloy.
(ii) Averaging procedure: The force constant tensors between a given pair of atoms calculated using the SQS will have nontrivial offdiagonal elements due to low symmetry. Moreover, due to the atomic relaxations in a local environment, a distribution of bond distances for a pair of species gives rise to a distribution in their force constants. To extract the force constant tensors having the symmetry of the underlying crystal structure of the disordered alloy, we average the calculated SQS forceconstant matrices for a particular set of displacements of the atoms from their ideal positions in the lattice, related by symmetry that will transform the force constant matrices to one displaying the symmetry of the underlying crystal structure. These averaged forceconstant tensors are then used as the inputs to the ICPA for the calculation of the configurationaveraged quantities. Below, we demonstrate the idea for a FCC lattice.
In a FCC lattice, the distances between a given atom and its 12 nearest neighbors are specified by the vectors , and ; being the lattice parameter. The force constant matrix for two atoms separated by the vector is of the form
The other nearest neighbor force constant matrices are of the same form and are related by point group operations. However, when the atoms are allowed to relax, the vector separating a pair of nearest neighbor atoms are modified to and the matrix of force constants corresponding to a given pair of atoms reflects the loss of symmetry, taking the general form
(1) 
For a particular set of , and , one needs to perform the averaging such that we obtain the following force constant matrix:
(2) 
To achieve this, one needs to pick up the relevant ones out of the symmetry operations corresponding to the cubic group which transform the force constant matrices to one displaying the symmetry of the FCC structure. To elaborate on this, we provide the following example:
Suppose, two atoms are separated by the separation vector . The force constant matrix obtained from firstprinciples calculations on the SQS is given by Eq. 1. To recover the symmetry of the FCC force constants, we need to find out the transformations that transform the separation vector to a new separation vector and . These transformations, in this particular case, mimic the nearest neighbor separation for the unrelaxed FCC lattice. They are:
(3) 
Corresponding to each of these transformation, there is a transformation matrix . We transform the SQS forceconstant matrix to the FCC forceconstant matrix by performing the operation in each of the 4 cases and then adding them. This simple procedure produces , and . All other offdiagonal elements vanish due to this symmetrization and the FCC symmetry is recovered. This procedure is repeated for all pairs of force constants and an arithmetic average is finally computed.
(iii)The ICPA: The ICPA is a Green’s function based formalism that generalizes the SCPA by considering scattering from more than one site embedded in an effective medium within which the effect of this correlated disorder is built in. The medium is constructed in a selfconsistent way so that site translational invariance and analyticity of the Green’s function are ensured. The analyticity of the Green’s function is ensured by using the principles of the traveling cluster approximation [12] which showed that once there are nondiagonal terms in the Green’s function, the selfenergy must include itineration of the scatterer through the sample to preserve analyticity. This means that the physical observables would be site translationally invariant. This required translational invariance is ensured by expressing the operators associated with the physical observables in an extended Hilbert space which accounts for the statistical fluctuations in site occupancies due to disorder.
We employ firstprinciples plane wave projector augmented wave method within generalized gradient approximation (GGA) [13] as implemented in the VASP code [14] for the calculation of the HellmanFeynman forces on the atoms in a 64atom SQS cell (4x4x4 supercell of the primitive fcc unit cell) with the optimized lattice parameter of 3.93 Å. Our value of the lattice parameter is slightly larger than the experimental value of 3.78 Ådue to the use of GGA. The SQS generated structure has zero short range order in the first four neighboring shells indicating a homogeneous disorder. The cutoff energy for the electronic wave functions is 400 eV. 432 kpoints were used in the irreducible Brillouin zone for the calculations of the force constant matrix. We have tested the convergence of our results with respect to the relevant parameters and have concluded that our choice yields quite accurate results. To obtain the forceconstant matrix, first, the equilibrium geometry is obtained by relaxing the atomic positions in the SQS cell until the forces converged to 10 eV/Å. Then each atom in the SQS cell is moved by 0.01 Å from the equilibrium position along three cartesian axes and forces on the atoms are calculated. Due to the lack of any symmetry in the SQS cell, 64 different force constant matrices were generated by using the PHON code [15]. The symmetry averaging procedure is then used to obtain an effective 3x3 force constant matrix to be used in ICPA calculation. In ICPA, the disorder is considered in the nearestneighbor shell only as the further neighbor force constants have one order of magnitude less. The ICPA calculations are done with a mesh and 1000 energy points.
In Fig. 1, the variations of the number of nearest neighbors of a given type are plotted as a function of the bond distances. It is observed that the interatomic bond distances depend sensitively on the number of unlike atoms, which clearly proves the dependence of the bond distances on local environments. As a result of this, the interatomic force constants too undergo variations. The calculated average NiNi (), NiPt () and PtPt () bond distances are 2.64 Å, 2.66 Å and 2.73 Å respectively, indicating a significant dispersion among the three pairs of bonds. A comparison to the unrelaxed alloy bond distances reveal that the relaxed remains the same, is less than smaller and is larger. The force constants for three different models of disorder, viz., SQSaveraged, empirical [5] and SCPA are presented in Table 1. In Ref. [5], it was argued that the NiNi bonds in the alloy are softer compared to those in pure Ni, the PtPt bonds are stiffer compared to those in pure Pt and NiPt bonds are even softer than the NiNi ones because the NiPt bond distances were thought to be the largest. A comparison between the force constants obtained by the SQSaveraged and by this empirical scheme shows that the NiNi force constants computed by the SQSaveraging scheme are softer, the NiPt and the PtPt force constants are and harder on an average. More importantly, the NiPt force constants computed by the SQSaveraging scheme are harder than the NiNi ones, a result in contradiction with the empirical scheme. This difference occurs due to the proper inclusion of environmental disorder in our case. in the alloy as computed by the SQS are about larger than that in pure Ni and thus these bonds suffer most severe dilution in strength, reduces by about compared to that in pure Pt and thus they are harder by about only compared to those in pure Pt. Inspite of nearly same bond distances, the NiPt bonds computed by the SQS are stiffer than the NiNi bonds, due to the following reasons: in case of the NiPt pairs, the Ni atoms find much larger Pt atoms as their nearest neighbors, roughly with the same available space as they get in case of NiNi pairs. As a result, the smaller Ni atoms try to accommodate bigger Pt atoms within the same volume as that available for NiNi pairs, resulting in a hardening of NiPt nearest neighbor interactions compared to the NiNi ones.

Pair SQS SCPA Empirical NiNi 8231 19365 15587 NiPt 17868 ” 13855 ” PtPt 33494 ” 28993 ” NiNi 525 3255 436 NiPt 2820 ” 348 ” PtPt 6854 ” 7040 ” NiNi 9580 22679 19100 NiPt 20740 ” 15280 ” PtPt 39655 ” 30317 ”
In Fig. 2, we present the phonon dispersion curves of NiPt alloy computed using force constants obtained by three models of disorder. The configuration averaging in case of the SQS and the empirical scheme are done by the ICPA. SQSICPA shows the best agreement with the experiments as compared to the SCPA and the empiricalICPA [5].The disorderinduced widths computed by the SQSICPA method (the shaded region in Fig. 2) also agree reasonably well with the experiments. The comparisons with the experiments were possible only for the direction because experimental results were not available for other directions. A closer look at Fig. 2 reveals that for the lowfrequency branches, the empiricalICPA and the SQSICPA methods are in close agreement while there are significant differences for the high frequency branches. The massdisorder treated in SCPA, on the other hand, fails to reproduce the experimental features of the dispersion relations both qualitatively and quantitatively. The frequencies of the high frequency branches (both transverse and longitudinal) computed by the SCPA are severely overestimated while the lower frequency branches extend all the way to the zone boundary, thus displaying a splitband behavior in the phonon dispersions which is not observed experimentally.
All these observations can be understood in terms of the interatomic force constants (shown in Table I) which are influenced by the fluctuations in the local environments. In the SCPA, the fluctuations in the force constants are completely neglected and hence NiNi, NiPt and PtPt all have the same force constants. In earlier studies [10, 5], the force constants used in the SCPA were those of pure Ni ones obtained experimentally [16]. In this study we have used a more realistic set of force constants for the SCPA calculations which are where are the cartesian directions and is the concentration of Ni. are the SQSaveraged force constants where denote atomic species. Such a choice of force constants has been used to incorporate, in an average way, the effects of alloy environment. The results, nevertheless, suggest that unless the fluctuations in the force constants are incoporated, the results do not even agree qualitatively with the measurements. The frequencies of the upper (lower) branches of longitudinal and transverse modes computed by SCPA are too high (low) as compared to the experiments. NiNi vibrations dominate the higher frequency branches and the average force constants being too high pushes the frequencies away from the ones measured by neutronscattering while the lower frequency branches being solely due to the PtPt vibrations have lower frequencies due to the underestimation of the PtPt force constants. The empiricalICPA results on the other hand agree qualitatively with the experiments because there is no splitband like behavior as was seen in the SCPA. This is due to the consideration of the NiPt correlated vibrations which renormalize the spectral weights associated with the contributions from NiNi and PtPt pairs [5]. Moreover, the splitting of the vibrational branches found experimentally around as a signature of the existence of resonance mode, is reproduced and the frequencies of the highfrequency branches are better as compared to the SCPA. The NiNi(PtPt) force constants shown in Table I in the empirical model are softer(harder) as compared to the SCPA ones explaining the reason for better agreement of the phonon frequencies computed by the empiricalICPA model with the experimental results. However, the overestimation(underestimation) of the NiNi(PtPt) interactions and an incorrect qualitative estimation of the NiPt interactions as compared to the NiNi ones gives rise to significant discrepancies for the highfrequency longitudinal and optical branches. With SQSICPA, one can have a much better quantitative agreement between theory and experiments as seen in Fig. 2. The highfrequency branches for both longitudinal and transverse vibrations computed by the SQSICPA method agree substantially with the experimental results. The normal mode frequencies for these branches are dominated by the vibrations of the Ni pairs and thus a softening of the NiNi bonds as computed by the SQS pushes the frequencies downwards compared to the empirical model making a better agreement with the experiments. Similarly, the relaxations of the Pt atoms result in the stiffening of the PtPt bonds, thus pushing the frequencies slightly upwards. However, the highfrequency transverse branch computed by the SQSICPA model is still overestimated. In the neutronscattering measurements [10], there were some ambiguities in determining the peak positions of the line shapes for the highfrequency transverse modes and thus the experimental results for this branch had larger uncertainties. Given this fact, the agreement between the theory and the experiment can be considered to be fairly good.
In conclusion, we have developed a reliable firstprinciples based approach for the calculation of phonon spectra in substitutional disordered alloys to treat mass, force constant and environmental disorder on equal footing. We demonstrate in case of NiPt alloy, the importance of an accurate structural model of disorder taking into account the role of fluctuations in the local environment through atomic relaxations in interpreting the microscopic features of the lattice dynamics for this class of complex alloys. The accurate modeling of the environmental disorder made possible by the SQS paves the way for a reliable description of phonon spectra in alloys with shortrange order where the forceconstants between a pair of species is dominated by a particular configuration of the nearest neighbor environment around an atom. Our future aim is to calculate the thermodynamic quantities [17] to compare with the experiments. Finally, we conclude that a combination of reliable force constants obtained from abinitio and ICPA as a selfconsistent analytic method for configurationaveraging enables us to solve the longstanding problem of theoretical computation of lattice dynamics in disordered alloys.
References
References
 [1] Bungaro C and Rabe K M 2003 Phys. Rev. B 68 134104.
 [2] Noda Y and Endoh Y 1988 J. Phys. Soc. Jpn. 57 4225.
 [3] Wolverton C and Ozolins V 2001 Phys. Rev. Lett. 86 5518.
 [4] Taylor D W 1967 Phys. Rev. 156 1017.
 [5] Ghosh S, Leath P L and Cohen M H 2002 Phys. Rev. B 66 214206.
 [6] Ghosh S, Neaton J B, Antons A H and Cohen M H 2004 Phys. Rev. B 70 024206.
 [7] Alam A, Ghosh S and Mookerjee A 2007Phys. Rev. B 75 134202.
 [8] Sluiter M H F, Weinert M and Kawazoe Y 1999 Phys. Rev. B 59 4100.
 [9] van de Walle A and Ceder G 2002 Rev. Mod. Phys. 74 11.
 [10] Tsunoda Y, Kunitomi N, Wakabayashi N, Nicklow R M and Smith H G 1979 Phys. Rev. B 19 2876.
 [11] Zunger A, Wei S H, Ferreira L G and Bernard J E 1990 Phys. Rev. Lett. 65 353; Lu Z W, Wei S H and Zunger A 1991 Phys. Rev. B 44 10470.
 [12] Mills R and Ratanavararaksha P 1978 Phys. Rev. B 18 5291.
 [13] Perdew J P, Burke K and Ernzerhof M 1996 Phys. Rev. Lett. 77 3865.
 [14] Kresse G and Hafner J 1993 Phys. Rev. B 47 RC558; Kresse G and Furthmüller J 1996 Phys. Rev. B 54 11169.
 [15] Alfé D 2009 Comp. Phys. Comm. 180 2622.
 [16] Dutton D H 1972 Can. J. Phys. 50 2915.
 [17] Grabowski B, Hickel T and Neugebauer J 2007 Phys. Rev. B 76 024309.