Density functional theory description of random Cu-Au alloys

Density functional alloy theory is used to accurately describe the three core effects controlling the thermodynamics of random Cu-Au alloys. These three core effects are exchange correlation (XC), ...


I. INTRODUCTION
Density functional theory (DFT) [1,2] is used in many areas of science owing to its accuracy coupled with good computational efficiency.Even the earliest exchangecorrelation (XC) approximation, the local density approximation (LDA) [3][4][5], is still being used today, but many modern and more advanced approximations, for example, metageneralized gradient approximations [6,7] and doublehybrid approximations [8], offer fairly consistent improvements on LDA.The predictive power of DFT is undeniable, but there exists a good number of situations and paradigms for which a direct application of DFT is still too difficult or time consuming.One of the difficult areas is short-range order (SRO) [9][10][11], because its accurate description typically requires large unit cells of low symmetry.For best results the local lattice relaxations (LLRs) should be taken into account, which for systems with large atomic size mismatches, such as Cu-Au, complicates the ab initio calculations even further.
Short-range order is a fundamental, but not yet fully understood, phenomenon in alloy theory.For example, SRO makes an important contribution to the physical properties of transition-metal-alloyed solid solutions [12][13][14][15][16][17], semiconductor alloys [18], and even high-entropy alloys [19,20].Experimental and theoretical studies of SRO, especially in the case of binary alloys, have a long history, and one of the earliest studies of SRO focused on the Cu-Au system [21], which is probably the most widely studied binary system in alloy theory [22,23].Short-range order and its effect on the * levamaki@kth.seorder-disorder transition temperature has mostly been studied using the cluster-expansion method [24][25][26], cluster-variation method [27], or effective-medium theory [28].Studies based on direct DFT calculations are much more scarce because compared to, e.g., the cluster-expansion method they are computationally much more demanding and there needs to be good approximations for the different free-energy components of each calculated structure.Nevertheless, the indisputable advantage of a direct DFT approach is that all three core effects, which are SRO, LLR, and XC, are all available on equal footing, without further approximations.
Such a unified DFT approach is now becoming possible for metallic alloys such as Cu-Au.In the past it was difficult for the generalized gradient approximation (GGA) [29][30][31][32][33] level XC functionals to produce acceptably accurate results for Cu-Au and other binary alloys [23].In our previous work we have shown that for metallic alloys the XC problem can be solved without abandoning the computationally tractable GGA level [34], by using the quasinonuniform XC approximation (QNA) scheme [35,36].QNA is an XC scheme designed for metallic alloy and compound calculations, which use element-specific optimal parameters for an improved GGA-level description.Later, we found that QNA in combination with LLR still cannot reproduce experimental findings on a quantitative level [37].The deviations were ascribed to the missing SRO effects, which should be considered on equal footing with the LLR and XC effects.
In this paper we include all relevant effects within one DFT scheme, putting forward a unified alloy theory for systems with sizable XC, LLR, and SRO effects.We use this direct DFT approach to predict the ordering properties of random Cu 1−x Au x alloys.The ability to consider all significant effects within one DFT program allows for a precise mapping of their relative impacts and draws conclusions for future alloy modeling.

A. Computational details
Calculations were performed with the grid-based projector augmented-wave DFT code GPAW [38,39], which is based on the PAW formalism and the atomic simulation environment (ASE) [40,41].A linear combination of atomic orbitals (LCAO) mode [42] with a p-valence basis [43] was used.The k points were generated using a 3 × 3 × 3 Monkhorst-Pack grid [44].The width of the Fermi-Dirac-type smearing of the occupation numbers was set to 0.01 eV.To describe the XC effects both the Perdew-Burke-Ernzerhof (PBE) [45] and QNA XC functionals were used.The QNA scheme was implemented in GPAW following the Becke fuzzy cell approach [46].QNA divides to total system into atom-centered, Voronoi-type regions and applies a slightly different XC approximation to each region.These region-specific XC approximations are optimized based on which type of element the Voronoi region belongs.The total QNA XC functional can then be expressed as a superposition or sum of the atomcentered optimized XC functionals [35,36].The advantage of the QNA scheme is that it represents a way to circumvent the accuracy limits of conventional GGA-level approximations when it comes to metallic alloy and compound calculations.
Atomic positions were relaxed using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm as implemented in ASE.Relaxations were deemed finished when all forces were below 0.3 eV Å −1 , which ensured that the mixing energy would not change appreciably upon any further relaxation.Ground-state unit cells were found using equation-of-state fitting with the constraints that for Cu 0.75 Au 0.25 and Cu 0.25 Au 0.75 structures the cells remain cubic, and for Cu 0.5 Au 0.5 structures the c/a ratio is also optimized.

B. Free-energy model
Free energy F ({α i }, T ) at a finite temperature T is approximated as where E 0 is the 0-K ground-state DFT energy and S conf is the configurational entropy, which both depend on the so-called Warren-Cowley SRO parameters {α i } [21,[47][48][49][50]. Warren-Cowley parameters are very often used in binary A 1−x B x alloy research to quantify the degree of SRO, and experimentally the presence of SRO can be measured by x-ray diffraction or scattering technologies.The SRO parameters describe the probability of finding certain types of atoms on the ith coordination shell that surrounds a given reference atom [51].Mathematically, they can be represented as [52] where q = 2x − 1 and i is the (thermally) averaged pair correlation function between lattice sites and their ith coordi-FIG. 1.The α 1 , α 2 , and α 3 values of the generated SRO structures (middle and bottom panels).The top panel shows the computed configurational entropy of the generated structures as a function of α 1 .Blue "+" symbols denote all generated structures and black/orange those that correspond to maximal/minimal configurational entropy.Red symbols represent experimentally measured α values of Cu 3 Au near the order-disorder transition temperature.nation shell.For perfectly random structures all α i are zero, whereas for the fully ordered L1 2 (Cu 3 Au, CuAu 3 ) and L1 0 (CuAu) structures α i = −1/3 for i odd and α i = 1 for i even.We neglect vibrational contributions to the free energy, not only because they can be time consuming to calculate but also because their contribution compared to the configurational entropy is expected to be small [53][54][55].A straightforward numerical scheme was conceived in order to obtain S conf as a function of α 1 , α 2 , and α 3 , where dependency on α 4 and higher has been neglected.As Eq. ( 1) shows, S conf (per atom) is needed as a function SRO parameters α i .In other words, for each macrostate of a unique combination of α i values, we want to estimate the number of microstates that corresponds to that macrostate, in order to get the S conf (α i ).We are implicitly assuming that a unique combination of α i values corresponds to a unique structure, or if not, then at least that every structure with some given combination of α i values gives (up to some small variation) the same DFT total energy.Let us have a binary alloy with two atom types, A and B. Our starting point is the fully ordered 3 × 3 × 3 108-atom L1 2 or L1 0 structure whose S conf = k B ln /N = 0. Using bigger cell sizes would improve the estimation, but here we use 3 × 3 × 3 cells, because we want the entropy to be consistent with the structures that will be used in DFT calculations, which are 3 × 3 × 3. Now we consider randomly switching one pair of atoms.Starting from the fully ordered structure, it is easy to see that the multiplicity becomes where N A and N B are the total number of type A and B atoms, respectively.For the L1 0 structure N A = N/2 and N B = N/2, where N is the total number of atoms.For the L1 2 structure N A = 3N/4 and N B = N/4.The variable p denotes the number of pairs that have been switched.The entropy after one pair switch becomes which can be evaluated numerically.The Warren-Cowley parameters α i can be routinely computed for any structure.We can thus compute the α i after the pair switch to obtain S conf as a function of α i .We now consider switching two pairs of atoms and the multiplicity becomes For a large cell and small number of pair switches p it is very likely that the "disturbed" nearest-neighbor (NN) clusters are all isolated, which means that for a small p the α i change linearly.For example, the α 1 of L1 0 grows linearly as a function of p as The configurational entropy as a function of α i towards the fully ordered end is thus known to a high degree of accuracy.
The general formula for the entropy after p pair switches is It should be noted that we only need to consider switching unique pairs of atoms, because any structure can be reached by switching unique pairs of atoms.This means that switching some atom two or more times is not needed, since it does not provide any new states or any new useful information.It is easy to see that when we randomly change p = N A /2 (L1 0 ) or p = 3N B /4 (L1 2 ) of all the unique pairs we get fully random structures, The ideal fully random entropies are where c A/B is N A/B /N.Our goal is to use the above pair switching strategy to map S conf as a function of α i parameters.For each p in the range 0 p N A /2 (L1 0 ) or 0 p 3N B /4 (L1 2 ) we randomly generate a large number, say, N iters = 30 000, of such structures.We then define an "effective" whose purpose is to keep track of the multiplicity of a given unique α i combination across all p values, and to normalize it according to the number of iterations N iters .For example, with p = 0 the correct multiplicity ˜ ( The effective multiplicity gives the entropy as a function of To decide which α i combinations are unique and which ones are not calls for some kind of binning of α i values.We have defined the bins to be the range −1/3 α 1 0 divided by 1000 for α 1 , and α 2 and α 3 are binned similarly.As a final step, due to the small size of the 108-atom unit cell, we scale the numerically computed entropies so that the maximum value becomes S random .

III. RESULTS AND DISCUSSION
The Cu-Au system is characterized by ordered structures at lower temperatures and random face-centered-cubic (fcc)  [25,28].CuAu has an additional, more complicated ordered CuAu-II phase, which is only stable between 656 and 683 K [25,28], but this phase is not part of our simulations.
Using the scheme of Sec.II B, we generate a large group of 108-atom SRO structures whose size is 3 × 3 × 3 times the conventional fcc unit cell.Figure 1 and Table I show the α 1 , α 2 , and α 3 values of the generated structures and the corresponding configurational entropies for L1 2 -based structures (Cu 0.75 Au 0.25 and Cu 0.25 Au 0.75 ).Similar results were obtained for the L1 0 structure as well, and those α values can be seen in Table I.For maximal information one would have to compute the DFT energy of every generated structure.This is, however, computationally far too prohibitive, so here we choose to decrease the number of calculated structures to a manageable number by parametrizing α 2 and α 3 , thus FIG. 2. Configurational entropy per atom of L1 2 and L1 0 binary systems as a function of α 1 .Polynomial fits to the data are drawn using black lines.The errors of the fits in terms of the R-squared value are also shown.a References [25,28].
leaving only α 1 as a free parameter.How α 2 and α 3 can be parametrized (as a function of α 1 ) is based on the following rationale.
In Fig. 1 the blue "+" symbols represent the α 1/2/3 distribution of all generated structures and the black dots are structures with such α 2/3 values that they yield maximum S conf for a given α 1 .We can see that the structures that have the maximum S conf are those structures that maximize disorder, i.e., the magnitudes of α 2 and α 3 are minimized.Conversely, the minimum entropy structures (orange symbols) are those that maximize order (magnitudes of α 2 and α 3 are maximized).The red symbols represent experimentally measured α 1/2/3 values of Cu 3 Au around the T c of 663 K.There is a convincing agreement in the α 1/2/3 values between the calculated maximum S conf structures and the experimental ones.This makes sense given that the S conf term stabilizes disordered structures at high temperatures most, and therefore stable structures at such temperatures are those with the highest S conf .The same trend, though not shown, is also true for CuAu and CuAu 3 .We note that our numerical results can also explain the experimentally observed anomalously small value of α 3 [57]; the value of α 3 is minimized in order to maximize configurational entropy at high temperatures.It is interesting to note that by normalizing the α values of Ref. [57] to α 0 = 1 (Norm.Moss 1964 in Fig. 1 and Table I), as has been done in Ref. [56], they fall nicely in line with other measurements, as well as our calculated maximum entropy α values.
Based on the above reasoning we choose to calculate the four most representative structures along the maximum entropy path, which in Fig. 1 are shown by the blue "×" symbols.By doing so we are able to fix the values of α 2 and α 3 using the maximum entropy condition and represent the free energy as a function of α 1 only.In order to make the free energy a continuous function of α 1 , we fit a seventh-order polynomial to the maximum entropy curve of the top panel of Fig. 1.By imposing the boundary conditions S conf (α 1 = −1/3) = 0, S conf (α 1 = 0) = S random , and ∂S conf /∂α 1 (α 1 = 0) = 0, we obtain where S random is the ideal fully random configurational entropy, and A, B, C, D, and E are the fitting coefficients whose values are listed in Table II for x = 0.25/0.75(L1 2 ) and x = 0.5 (L1 0 ).The graphical results of the fits using the fitting coefficients of Table II are shown in Fig. 2. The entropy curves of Fig. 2 resemble quite closely those of Ref. [58], which were calculated within the cluster-variation method.Both our curves and those of Ref. [58] are seen to have a "shoulder" around α 1 ≈ −0.15, which corresponds to the experimental α 1 value near the T c The calculated 0-K mixing energies are shown in Fig. 3, where the solid curves are a spline fit through the calculated DFT mixing energies.Previous studies have shown that the magnitudes of PBE mixing energies of CuAu alloys are FIG.3. The top panels show the mixing energies of Cu-Au alloys as a function of the Warren-Cowley α 1 parameter.Experimental mixing energies are from Ref. [48] (ordered at 320 K and random at 720 K) and in the top panels they have been coupled with experimental α 1 values from Refs.[52,56].The middle and bottom panels show temperature-dependent mixing energies of Cu-Au alloys as functions of temperature and the Warren-Cowley α 1 parameter.The black star symbols indicate the minima of the curve at a given temperature.underestimated quite significantly [23,34].Here, we report that this underestimation happens at all degrees of SRO.It has also been shown previously that the QNA scheme can produce vastly improved mixing energies while managing to stay on the same GGA level as PBE [34].At temperatures around the T c there is a good agreement between QNA and experimental α 1 versus mixing energy data points, which indicates that QNA predictions about the stability and the degree of SRO around the T c are in good agreement with experiments.
The middle and bottom rows of Fig. 3 show the free mixing energy F at various temperatures.By finding the energy minimum (black stars in Fig. 3) for a dense mesh of temperatures we can compute how α 1 evolves as a function of temperature T .These temperature-dependent α 1 curves are plotted in Fig. 4. The PBE and QNA curves are very similar and there is practically only a shift toward a higher T c from PBE to QNA, which means that there is no significant difference in the predicted transition α 1 values between PBE and QNA.We extract the T c from the data of Fig. 4 and list the resulting T c in Table III.It turns out that our DFT model, together with the numerically estimated S conf , is able to reproduce the experimental trend T c (x = 0.25) ≈ T c (x = 0.5) > T c (x = 0.75), which seems to have been previously possible only by cluster-expansion [26] or effective-medium theory methods [28].
Having access to the unified DFT model, we can now measure the impact of various core effects on the critical temperature.Table III lists the T c with LLR, SRO, or both either enabled ( ) or disabled ( ).The SRO-free T c is calculated by considering only the fully ordered (α 1 = −1/3) and fully random (α 1 = 0) structures and then calculating their energy difference divided by the ideal fully random entropy S random .Table III shows that the effect of SRO is to decrease the predicted T c by a significant amount.By repeating our calculations without relaxing atomic positions we can probe the effect LLR on T c .For the CuAu system, with its large atomic size mismatches, the LLR effect turns out to be the most crucial.Without LLR the predicted T c are highly overestimated, with and without taking XC or SRO effects into account.This is not surprising given the large size mismatch between Cu and Au atoms.When it comes to XC effects, despite predicting inaccurate mixing energies, PBE gives very accurate estimates for T c of Cu 0.75 Au 0.25 and Cu 0.5 Au 0.5 , while T c of Cu 0.25 Au 0.75 is heavily underestimated.The surprisingly good performance of PBE for T c can be attributed to the fact that even though the PBE mixing energies in the top row of Fig. 3 are highly underestimated, this large error cancels because the slopes of the PBE curves are not bad, and T c depends on the slope of the curve and not on its position on the mixing energy axis.QNA slopes are a bit larger than those of PBE, and QNA therefore somewhat overestimates T c of Cu 0.75 Au 0.25 and Cu 0.5 Au 0.5 and only slightly underestimates it for Cu 0.25 Au 0.75 .Overall, QNA gives balanced estimates for all three alloy concentrations, which are actually very similar to those of SRO-free PBE, and according to Table III it has the smallest maximum absolute relative error across all three concentrations.

IV. CONCLUSIONS
In summary, using a unified DFT approach, we studied the thermodynamics of Cu-Au alloys.The LLR effects are the strongest, but accurate descriptions of XC and SRO are also needed in order to quantitatively predict the order-disorder transition temperatures.Coupled with numerical estimates for the configurational entropy, our direct DFT approach can predict the transition temperatures with an accuracy that previously only specialized cluster-expansion or effective-medium theory methods could.One great advantage of the present model is that it can be extended to multicomponent systems, such as high-entropy alloys.

FIG. 4 .
FIG.4.Evolution of short-range order parameter α 1 as a function of temperature T .Dashed lines mark the order-disorder transition temperatures.

TABLE I .
Experimental and present calculated α values.Most of the Cu 0.75 Au 0.25 values can also be seen in Fig. 1.

TABLE II .
Polynomial fitting coefficient of the maximal configurational entropy.disorder transition temperature T c .Experimentally, the stable ordered structures of Cu 3 Au, CuAu, and CuAu 3 at low temperatures are L1 2 , L1 0 , and L1 2 , respectively.The experimental T c of Cu 3 Au, CuAu, and CuAu 3 are 663, 656, and ≈500 K, respectively

TABLE III .
Order-disorder transition temperatures T c (in K) and their relative errors (in %) compared to experimental data.Solid (open) tick boxes indicate which effects have been included (excluded) when calculating the transition temperatures.