Cu-Au, Ag-Au, Cu-Ag and Ni-Au intermetallics: First-principles study of

0 downloads 0 Views 393KB Size Report
Sep 14, 1997 - and Cu-Au to form compounds and Ni-Au and Cu-Ag to phase separate at T = 0 K. ..... act for long-period superlattices, but represents a choice.
Cu-Au, Ag-Au, Cu-Ag and Ni-Au intermetallics: First-principles study of phase diagrams and structures

arXiv:cond-mat/9710225v1 [cond-mat.mtrl-sci] 21 Oct 1997

V. Ozoli¸ nˇs, C. Wolverton, and Alex Zunger National Renewable Energy Laboratory, Golden, CO 80401 (September 14, 1997)

The classic metallurgical systems – noble metal alloys – that have formed the benchmark for various alloy theories, are revisited. First-principles fully relaxed general potential LAPW total energies of a few ordered structures are used as input to a mixed-space cluster expansion calculation to study the phase stability, thermodynamic properties and bond lengths in Cu-Au, Ag-Au, CuAg and Ni-Au alloys. (i) Our theoretical calculations correctly reproduce the tendencies of Ag-Au and Cu-Au to form compounds and Ni-Au and Cu-Ag to phase separate at T = 0 K. (ii) Of all possible structures, Cu3 Au (L12 ) and CuAu (L10 ) are found to be the most stable low-temperature phases of Cu1−x Aux with transition temperatures of 530 K and 660 K, respectively, compared to the experimental values 663 K and ≈ 670 K. The significant improvement over previous first-principles studies is attributed to the more accurate treatment of atomic relaxations in the present work. (iii) LAPW formation enthalpies demonstrate that L12 , the commonly assumed stable phase of CuAu3 , is not the ground state for Au-rich alloys, but rather that ordered h100i superlattices are stabilized. (iv) We extract the non-configurational (e.g., vibrational) entropies of formation and obtain large values for the size mismatched systems: 0.48 kB /atom in Ni0.5 Au0.5 (T = 1100 K), 0.37 kB /atom in Cu0.141 Ag0.859 (T = 1052 K), and 0.16 kB /atom in Cu0.5 Au0.5 (T = 800 K). (v) Using 8 atom/cell special quasirandom structures we study the bond lengths in disordered Cu-Au and Ni-Au alloys and obtain good qualitative agreement with recent EXAFS measurements. PACS numbers: 61.66.Dk, 71.20.Gj, 81.30.Bx

T = 0 ground state structures and finite-temperature thermodynamic quantities can be predicted without any empirical information, it is interesting to take a global look at the noble metal alloy family. Table I summarizes some of the salient features1–4,14,15,18,67–69 of the four binary systems Cu-Au, Ag-Au, Cu-Ag and Ni-Au. We included the relative lattice constant mismatch ∆a/a = 2 |aA − aB | / |aA + aB | between the consituents,67 the electronegativity difference ∆χ = χA −χB on the Pauling scale,68 the mixing enthalpy of the equiatomic alloy,2,18 the sign of the calculated nearest neighbor pair interaction J2 (present study), the structural identity of the lowtemperature phases1–4,67 and the order-disorder transition (or miscibility gap) temperatures2,69 Tc . Some interesting observations and trends which we will attempt to reproduce theoretically, are apparent from this general survey: (i) Despite a large (12%) size mismatch in Cu-Au, and a small (≈ 0%) size mismatch in Ag-Au, both systems form ordered compounds at low temperatures and have negative mixing enthalpies, suggesting attractive (“antiferromagnetic”) A–B interactions. Thus, when the difference in the electronegativity ∆χ of the constituent atoms is sufficiently large, as it is in CuAu and AgAu, size mismatch apparently does not determine ordering vs. phase separation tendencies. (ii) Despite a similar size mismatch (12%) in Cu-Au and Cu-Ag, the former orders while the latter phaseseparates. Thus, the existence of large electronegativity

I. INTRODUCTION: CHEMICAL TRENDS IN NOBLE METAL ALLOYS

Noble metal alloys are, experimentally, among the most studied intermetallic systems.1–24 In addition, the Cu-Au system has been considered the classic paradigm system for applying different theoretical techniques of phase diagram and phase stability calculations.25–63 Most notably, this system has been considered as the basic test case for the classic Ising-hamiltonian statistical-mechanics treatment of alloys.25–32 More recently, noble metal binary alloys have been treated theoretically via empirical fitting of the constants in Ising hamiltonians,25–34 semiempirical interatomic potentials,35–47 and via first-principles cluster expansions.48–55 The essential difference in philosophy between the classical application of Ising models to CuAu25–30,33 and more modern approaches based on the density functional formalism64 is that in the former approach the range and magnitudes of the interactions are postulated at the outset (e.g., first or second neighbor pair interactions), while the latter approaches make an effort to determine the interactions from an electronic structure theory. However, despite recent attempts,48–54 it is still not clear whether the noble metal alloys can be essentially characterized as systems with short-range pair interactions, or not. Now that first-principles cluster expansion approaches65,66 have advanced to the stage where both

1

TABLE I. Major physical properties of Ag-Au, Cu-Ag, Cu-Au and Ni-Au alloys. We give constituent size mismatches, ∆a/a = 2(aA −aB )/(aA +aB ), electronegativity differences on the Pauling scale,68 ∆χ, mixing enthalpies of the disordered alloys at the equiatomic composition, ∆Hmix (x = 21 ), signs of the nearest-neighbor Ising interaction, J2 , order-disorder transition temperatures (or miscibility gap temperatures for Cu-Ag and Ni-Au), Tc (x = 21 ), and excess entropies of solid solutions, form ∆Stot − ∆Sideal . All phases are fcc-based. System

∆a/aa

∆χb

Cu-Au Ag-Au Cu-Ag Ni-Au

12% 0% 12% 15%

0.64 0.61 0.03 0.63

∆Hmix (x = 1/2) (meV/atom) −91c −48d +80e +76f

J2

Low-T phasesg

>0 >0 0

L12 , L10 , L12 (?) L12 , L10 , L12 Phase sep. Phase sep.

Tc (x = 12 ) (K) 683g 115-168h > Tm 1083d

form ∆Stot − ∆Sideal g (kB /atom) +0.36 −0.17 +0.04 +0.35

a

Ref. 67. Ref. 68. c Refs. 15, 14, 2. d Ref. 2. e Theoretically calculated value from this work. f Refs. 2, 18. g Refs. 2, 4. h Ref. 69. b

and CuAu24 show distinct RAA 6= RAB 6= RBB bond lengths, which need to be explained. In this work we will analyze the above mentioned trends in terms of a first-principles mixed-space cluster expansion,65,66 based on modern local density approximation (LDA) total energy calculations. We reproduce the observed trends (i)-(vi) in ordering preferences, mixing enthalpies ∆Hmix , transition temperatures Tc and interatomic bond lengths. In addition, we predict new, hitherto unsuspected ordered phases in Au-rich Cu-Au alloys.

difference in Cu-Au (as opposed to the small difference in Cu-Ag), seems to induce ordering tendencies. (iii) Cu-Ag and Ni-Au both phase-separate (and have positive ∆Hmix ) as they have large size mismatches. Yet, Ni-Au having a large electronegativity difference, shows an ordering-type nearest-neighbor pair interaction (J2 > 0), just like the compound forming Cu-Au and AgAu, while Cu-Ag has a clustering-type nearest-neighbor interaction (J2 < 0). Thus, the sign of J2 does not reflect the low temperature ordering vs. phase separation. expt − ∆Sideal (iv) The amount ∆SXS = ∆Stot expt 2 by which the measured entropy ∆Stot deviates from the ideal configurational entropy ∆Sideal = kB [x log x + (1 − x) log(1 − x)], is unexpectedly large in Cu-Ag and Ni-Au, indicating a large non-configurational entropy of formation. Other interesting facts about the noble metal binary intermetallics include: (v) Despite numerous studies,1–4,7,8,10–12 the structure of the ordered phases in Au-rich Cu-Au is not well established yet. It is often assumed1–4 that the stable Au-rich low-temperature phase is CuAu3 in the L12 structure, but direct experiments7,8,10 below the orderdisorder transition temperature Tc (x = 34 ) ≈ 500 K are difficult because the diffusion rates are very low and even the best ordered samples contain significant disorder. Possible further thermodynamic transformations at lower temperatures may be kinetically inhibited. (vi) The trends in bond lengths vs. composition are non-trivial. Traditionally, all coherent-potentialapproximation based theories70–72 of intermetallic alloys have assumed that the nearest-neighbor bond lengths are equal, RAA = RAB = RBB , and proportional to the average lattice constant. Recent theories73–75 suggested, however, that bond lengths relax in the alloy to new values, and this has a significant effect on the electronic structure.53,76,77 Recent EXAFS experiments on NiAu23

II. BASIC IDEOLOGY AND METHODOLOGY

There are many problems in solid state physics that require knowledge of the total energy E(σ) of a lattice with N sites as a function of the occupation pattern σ of these sites by atoms of types A and B. This information is needed, for example, in the ground state search problem,72 where one seeks the configuration with the lowest energy at T = 0 K. {E(σ)} is also needed for calculating the temperature- and composition-dependent thermodynamic functions and phase diagrams of an A1−x Bx alloy. A direct, quantum-mechanical calculation of the toˆ tal energy Edirect (σ) = hΨ|H|Ψi/hΨ|Ψi (where Ψ is the ˆ is the manyelectronic ground state wave function and H electron Hamiltonian) is possible only for a limited set of configurations σ. This is so because (i) the computational effort to solve the Schr¨odinger equation for a single configuration scales as the cube of the number of atoms per unit cell, so that only small unit cells can be considered, (ii) there are 2N configurations, and (iii) for each configuration, one has to find the atomic relaxations δumin (σ) which minimize the total energy. Consequently, one searches for a “cluster expansion” (CE) that 2

Here Nf is the number of nonzero effective interactions and Πf (σ) are lattice averages of the spin products in configuration σ. Sanchez, Ducastelle and Gratias80 have shown that there is a set of composition-independent interactions for Eq. (3) which can exactly reproduce the directly calculated total energies of all configurations σ. This statement is strictly true if all possible clusters are included in Eq. (3), and should hold for the truncated series Eq. (3) if the cluster expansion is well converged. Several methods81,82 yield concentration-dependent effective interactions, providing in principle equally valid schemes for representing ∆Hdirect (σ) in terms of a cluster expansion. In the present work, we select compositionindependent interactions, since these can be directly compared to previous Ising model treatments25–34,48–55 of the noble metal alloy phase diagrams. A number of issues arise in trying to construct a cluster expansion that satisfies Eq. (2): (i) The number of interactions and their types (pair, multibody) cannot be decided arbitrarily, but must be constrained by a microscopic electronic-structure theory according to Eqs. (1) and (2). (ii) In most configurations σ, atoms move away from the ideal lattice sites, which not only lowers the total energies Edirect (σ), but also slows down the convergence85 of the expansion Eq. (3). The solution is to have a cluster expansion with many interaction terms NJ that can represent such situations. We accomplish this by using a reciprocal space formulation, formally equivalent to an infinite number of real-space pair interactions. (iii) Some cluster expansions78 require that the number of interactions NJ must equal the number of configurations Nσ whose total energies need to be evaluated via the direct electronic-structure method. The number of such calculations may be excessive in view of (ii). We thus introduce a method in which Nσ ≪ NJ . Furthermore, interactions that are not needed to satisfy Eq. (2) are automatically discarded. (iv) One has to deal with the macroscopic elastic strain leading to a k → 0 singularity in the Fourier transform of the pair interactions, X (6) Dj Jpair (Ri − Rj )e−ikRj , Jpair (k) =

accurately reproduces the results of a direct, quantummechanical (e.g., LDA) calculation ECE (σ) ∼ = Edirect (σ),

(1)

without the unfavorable scaling of the computational cost with the size of the unit cell. In designing a cluster expansion, there are few choices of independent parameters. For example, one could choose to obtain a cluster expansion for the volume(V ) dependent equation of state Edirect (σ, V ) [see, e.g., Refs. 52, 78, 79], or to find a cluster expansion for the energy at the volume Vmin (σ) that minimizes Edirect (σ, V ). We choose the latter possibility. Furthermore, for each configuration σ, we wish to reproduce the total energy corresponding to the fully relaxed cell shape and atomic positions {δumin(σ)}. In other words, we choose to represent ECE (σ) ∼ = Edirect [σ; δumin (σ); Vmin (σ)] ≡ Edirect(σ). (2) Note that by focusing on the equilibrium energy of each configuration, we give up the possibility of studying nonequilibrium geometries (e.g., bond lengths) and equations of state. Instead, for each occupation pattern σ, we can find the total energy E(σ) of the atomically relaxed and volume-optimized geometry. The best-known cluster expansion is the generalized Ising model in which the equilibrium total energy of an arbitrary configuration σ is expanded in a series of basis functions defined as pseudospin products on the crystal sites: X 1X Ji Si + E(σ) = J0 + Jij Si Sj 2 i i6=j

1 X Jijk Si Sj Sk + . . . , + 3!

(3)

i6=j6=k

where in binary A1−x Bx alloys Si = +1 or −1, depending on whether the site i is occupied by an atom of type A or B. Equation (3) is valid whether the lattice is relaxed ot not, as long as a one-to-one correspondence exists between the actual positions of atoms and the ideal fcc sites. The practical usefulness of the cluster expansion Eq. (3) rests on the assumption that the effective cluster interactions (ECI’s), Jij , Jijk , . . . , are rapidly decreasing functions of the number of sites and intersite separation, so that only a finite number of terms has to be kept in Eq. (3) for the desired accuracy. In this case, we can write the formation enthalpy of structure σ, ∆Hdirect (σ) = E(σ) − xEA − (1 − x)EB ,

j

where Dj is the number of {Ri , Rj } pairs per lattice site. As shown by Laks et al.65 (see also the discussion below), in size mismatched systems the correct Jpair (k) depends on direction kˆ in the long-wavelength limit k → 0. To solve this, we express Jpair (k) as a sum of two parts,

(4)

where EA and EB are total energies of the pure constituents A and B, as the following cluster expansion (CE): ∆HCE (σ) = J0 +

Nf X

Df Jf Πf (σ).

ˆ Jpair (k) = JSR (k) + JCS (k),

(7)

where JSR (k) is an analytic function of k and can be obtained from short-ranged real space pair interactions, ˆ contains the nonanalytic behavior around while JCS (k)

(5)

f

3

ˆ To exk = 0 and depends only on the direction k. plain this singularity, we consider the energy of a coherent An Bn superlattice, formed by a periodic stacking b By of n layers of A and n layers of B in direction G. introducing the structure factor, X (8) Sj e−ikRj , S(k, σ) =

ELDA [σ, δumin (σ); Vmin (σ)] ≡ ELDA (σ) at the atomically relaxed geometry and equilibrium volume of configuration σ. The function ECE (σ) we consider includes composition- and volume-independent interactions, so as to maintain maximum similarity with the classical Ising model. The number and type of interactions are not decided arbitrarily, but are constrained by the electronic structure theory used (here, the LDA). Relaxation is treated accurately by including long-range pair interactions in the reciprocal space representation. The k → 0 singularity, affecting both short and long-period superlattices, is dealt with explicitly. The above requirements are satisfied by the mixed space cluster expansion (MSCE): X 2 ∆HCE (σ) = Jpair (k) |S(k, σ)|

j

the total pair interaction energy in Eq. (3) can be expressed as a reciprocal space sum: X 2 Epair (σ) = Jpair (k) |S(k, σ)| . (9) k

An Bn superlattice has nonzero structure factor at 1 b G, and therefore its energy is determined by k = 2n 1 b Jpair ( 2n G). As n → ∞ its formation energy is given by a sum of the epitaxial deformation energies of pure constituents needed to bring them to a common lattice b Since the epiconstant in the plane perpendicular to G. taxial deformation energy of pure constituents is direction dependent (e.g., it is easier to stretch Cu in [100] planes than in [111] planes, see Sec. III B), the formation energy ∆H(A∞ B∞ ) is also direction dependent. Therefore, limk→0 Jpair (k) must depend on the direction of approach to the origin, proving that Jpair (k) is singular. Physically, the nonanaliticity of Jpair (k) is caused by long-range interactions via macroscopic elastic strain and cannot be reproduced using finite-ranged real-space pair interactions, but must be accounted for explicitly in reciprocal space. If the singularity is neglected, then as explained in Ref. 65, the cluster expansion fails not only for long-period (n → ∞) superlattices An Bn , but also for those short-period (n > 2) superlattices which have not been explicitly included in the constraint Eq. (2). We emphasize that although the contribution of Jsing (k) to the formation energy is nonzero only in size-mismatched systems, it is not related to the atomic relaxation energy for a particular structure σ in any simple way (except if σ itself is a long-period superlattice). The singularity in Jpair (k) is similar to the singularity in the dynamical matrix Dαβ (κκ′ |k) of polar crystals in the long-wavelength limit,83 caused by long-range electrostatic interactions via macroscopic electric field. In lattice dynamics, Dαβ (κκ′ |k) is expressed as a sum of sing regular and singular parts, Dαβ (κκ′ |k) = Dαβ (κκ′ |k) + reg reg ′ ′ Dαβ (κκ |k), where Dαβ (κκ |k) (analytic as k → 0) is due to short-range force constants. The singular part sing Dαβ (κκ′ |k) gives rise to LO/TO splitting of the zonecenter optical frequencies ωΓ in polar crystals, and also ˆ in uniaxial leads to a directional dependence of ωΓ (k) crystals (e.g., CuPt-type GaInP2 ). These phenomena cannot be reproduced by any set of finite-ranged microscopic force constants, but have to be calculated explicitly using the macroscopic Maxwell equations.84 In summary, we seek to find a function ECE (σ) which accurately reproduces the LDA total energies

k

+

MB X

Df Jf Πf (σ) + ∆ECS (σ).

(10)

f

We have separated out the so-called equilibrium constituent strain energy term, ∆ECS (σ), which accounts for the k → 0 singularity.65 In Eq. (10) we do not need to calculate ∆ECS (σ) for each configuration σ, but only for the directions b k of the wave vectors with nonzero S(k, σ). In fact, it is constructed to coincide with the elastic strain energy of coherent superlattices in the longperiod limit:65 X 2 ∆ECS (σ) = JCS (x, b k) |S(k, σ)| , (11) k

eq k) (x, b ∆ECS JCS (x, b k) = , 4x(1 − x)

(12)

where S(k, σ) is the structure factor from Eq. (8). The eq quantity ∆ECS (x, b k) depends only on the direction b k, and will be given in Sec. III B. Equation (11) is exact for long-period superlattices, but represents a choice for short-period superlattices and non-superlattice (e.g., L12 ) structures. It has been found65 that the choice Eq. (11) improves the cluster expansion predictions also for short-period superlattices. Equation (10) is a generalized Ising model description of the formation energy of any relaxed configuration σ, even if a direct LDA calculation for this σ is impractical. The cluster interaction energies {Jpair (k)} and {Jf } are obtained by fitting Eq. (10) to the LDA formation energies. An additional smoothness requirement is imposed on Jpair (k), which ensures that the pair interactions are optimally short-ranged in real space. Namely, we minimize the sum 1 X 2 ∆2rms = wσ [∆HCE (σ) − ∆HLDA (σ)] Nσ σ λ/2  t X Jpair (k), (13) Jpair (k) −∇2k + α k

4

TABLE II. Definition of the small-unit-cell ordered structures used in the LDA total energy calculations.

Simple superlattices Composition AB A2 B AB2 A3 B AB3 A2 B2

(001) L10 (CuAu) “β1” (MoSi2 ) “β2” (MoSi2 ) “Z1” “Z3” “Z2”

(011) L10 (CuAu) “γ1” (MoPt2 ) “γ2” (MoPt2 ) “Y1” “Y3” “Y2”

Name

Prototype

L12 L12 D7a D4 D7b

Cu3 Au Cu3 Au

Orientation (111) L11 (CuPt) “α1” (CdI2 ) “α2” (CdI2 ) “V1” “V3” “V2”

(311) L11 (CuPt) “γ1” (MoPt2 ) “γ2” (MoPt2 ) “W1” “W3” “W2”

(201) L10 (CuAu) “γ1”(MoPt2 ) “γ2”(MoPt2 ) D022 (TiAl3 ) D022 (TiAl3 ) “40”(CuFeS2 )

Period

Reference

Other structures Composition A3 B1 A1 B3 A7 B A4 B4 AB7 A8 B AB8 A6 B2 A6 B2 A4 B4 A4 B4 A6 B2 A2 B6

D023 LPS-3 SQS8a SQS8b SQS14a SQS14b

Superlattice direction none none none none none none none (401) (601) (311) (311) (201) (201)

Ni8 Nb Ni8 Nb Al3 Zr

where λ and t are free parameters and α is a normalization constant.65 Typically we choose λ = 4 and t = 1, but the fit is not sensitive to this choice. This approach solves the four problems indicated above in the sense that (i) the fitting process itself automatically selects the pair interactions that are essential to obtain a good fit (process still does not select multibody figures), (ii) the pair interactions can be of arbitrary long range, facilitating treatment of systems with large elastic relaxations, (iii) the number of pairs can be much larger than the number of ordered structures in the fit, and (iv) the directly calculated constituent strain energy ∆ECS contains the k → 0 singularity. Unlike all CPA-based methods,70,71 the present approach includes full account of atomic relaxation and local environment effects. Unlike the classical Ising descriptions,25,27–33 the interaction energies are determined by the electronic structure rather than being guessed. Finally, unlike the computational alchemy linear response approach,85 multibody terms are included here. Having written the expression for the total energy of arbitrary configuration, Eq. (10), we can evaluate its constants from a limited number of LDA calculations on small unit cell (Natoms < 10) ordered structures with fully relaxed atomic positions. Equation (10) can then be employed in simulated annealing and Monte Carlo calculations86,87 yielding T = 0 ground states and T > 0 statistical and thermodynamic properties. Further details of the method are given in Sec. III.

(5,1,1,1) (5,1,1,1) (2,3,2,1) (3,2,1,2) (6,2) (2,6)

52 52 52 52 52 95 95 95 87 76 76 73 73

III. DETAILS OF THE METHOD A. T = 0 energetics

The calculations of T = 0 total energies employ the full-potential linearized augmented plane wave method88 (FLAPW). The basis set consists of plane waves in the interstitial region, augmented in a continuous and differentiable way with the solutions of the radial Schr¨odinger equation inside the non-overlapping muffin-tin spheres. Non-spherical potential and electronic charge density terms are calculated in all space and included in the Hamiltonian matrix. Core states are treated fully relativistically and recalculated in each self-consistency iteration. The wave equation for the valence states includes all relativistic effects except the spin-orbit interaction, i.e., they are treated scalar relativistically. FLAPW is the most accurate all-electron method, superior to the methods employing overlapping atomic spheres (atomicspheres approximation – ASA) and/or shape approximations to the potential. We use the Wigner exchange-correlation functional.89 As a check, we have performed several calculations using the Perdew-Zunger90 parametrization of the CeperleyAlder91 functional and the generalized gradient approximation of Perdew and Wang.92 We find (see Sec. IV A 1) that the various exchange-correlation functionals change the enthalpies of formation of ordered Cu-Au compounds by a negligible amount (less than 2 meV/atom).

5

TABLE III. LDA calculated formation [Eq. (4)] enthalpies for fcc superstructures (defined in Table II) of Ag-Au, Cu-Ag, Cu-Au and Ni-Au. The numbers in parentheses represent errors of the cluster expansion fit. All energies in meV/atom.

Structure

Ag-Au

Cu-Ag

Cu-Au

Ni-Au

LDA ∆Hunrel

LDA ∆Hunrel

LDA ∆Hrel

LDA ∆Hunrel

LDA ∆Hrel

LDA ∆Hunrel

LDA ∆Hrel

0.0 (−0.4) 0.0 (−0.5)

0.0 0.0

0.0 (−0.1) 0.0 (+0.3)

0.0 0.0

0.0 (+0.2) 0.0 (−0.4)

0.0 0.0

0.0 (+0.4) 0.0 (−0.2)

(−0.8) (−0.1) (+0.1) (−0.1) (+0.7) (−0.3)

+107.6 +130.2 +112.0 +126.4 +96.8 +164.7

(+0.1) (−2.6) (+0.6) (+0.3) (+1.8) (−1.0)

+20.3 (−0.1)

+98.1 +207.8 +151.7 +221.7 +142.0 +286.7 +273.3 +355.5 +576.2

+76.1 (+1.4) +105.7 (−0.1) +38.3 (+0.1) +89.9 (−4.2) +32.4 (+4.0) +70.2 (+0.1) +57.1 (−0.8) +63.2 (+0.7) +30.8 (0.0)

+60.3 +123.0 +86.4 +136.1 +79.5 +170.6

+32.5 +61.4 +2.1 +78.6 +5.1 +52.2 +95.8

(−0.1) (−7.7) (+7.7) (+4.1) (+0.8) (−2.5) (+0.3)

+192.3 +288.5 +200.9 +290.8 +172.8 +335.8 +576.2

+166.8 +202.2 +100.9 +193.7 +83.0 +162.4 +173.8

(+1.4) (−6.4) (+6.4) (+4.1) (+4.0) (−4.1) (+1.3)

(−0.6) (+0.8) (+3.5) (−1.3) (−1.1) (−1.2)

−14.2 +1.7 +21.8 +19.4 +59.5

−18.4 −6.7 −1.3 −1.0 −4.2 +66.1

(+3.3) (−5.2) (+3.8) (+0.1) (−2.0) (+0.3)

+123.3 +126.3 +148.5 +104.1 +192.3 +576.2

+98.9 +102.6 +99.2 +78.7 +96.6 +117.7

(−3.8) (+3.8) (+7.8) (+1.1) (−4.5) (+1.6)

+94.2 +91.4 +104.7 +65.9

(−0.2) (+9.0) (−4.4) (−1.4)

+22.0 +21.1 +15.7

+7.0 +7.8 −20.9 +69.5

(+1.5) (+0.6) (−1.0) (+0.4)

+125.7 +144.2 +576.2

+120.8 +88.4 +93.6 +119.8

(+5.2) (+5.3) (−5.3) (+1.9)

+85.1 +76.4 +107.5 +67.3

(+1.3) (−0.5) (−0.4) (+1.6)

−32.7 −10.6 −19.0

−32.8 −11.8 −23.0 +53.4

(+0.3) (−1.8) (−0.6) (−0.4)

+75.0 +68.7 +93.5 +576.2

+75.0 +68.6 +84.8 +84.8

(+5.6) (+1.5) (−3.6) (−2.0)

Superlattice

Name

A B

fcc fcc

(001) Struct: A1 B1 A2 B1 A1 B2 A3 B1 A1 B3 A2 B2 A2 B3 A3 B3 A∞ B∞

L10 “β1” “β2” “Z1” “Z3” “Z2” “Z5” “Z6”

(111) Struct: A1 B1 A2 B1 A1 B2 A3 B1 A1 B3 A2 B2 A∞ B∞

L11 “α1” “α2” “V1” “V3” “V2”

−43.0 (−0.4) −30.2 (0.0) −30.8 (0.0) −21.3 (+0.3) −21.4 (+0.6) −22.9 (−0.4) 0.0 (0.0)

+134.8 +152.4 +124.9 +145.9 +106.8 +177.1

+129.8 +120.4 +95.0 +108.4 +73.6 +109.1 +86.3

(−1.1) (−2.9) (+2.9) (+0.4) (+1.5) (−1.0) (−1.0)

(011) Struct: A2 B1 A1 B2 A3 B1 A1 B3 A2 B2 A∞ B∞

γ1 γ2 “Y1” “Y3” “Y2”

−49.7 (−0.4) −46.9 (+0.4) −37.0 (0.0) −35.4 (+0.6) −44.1 (−0.3) 0.0 (0.0)

+106.4 +97.2 +105.1 +85.5 +136.0

+100.3 +92.5 +85.4 +75.2 +105.7 +75.3

(113) Struct: A3 B1 A1 B3 A2 B2 A∞ B∞

“W1” “W3” “W2”

−35.9 (+0.5) −34.4 (−0.2) −50.6 (−0.1) 0.0 (0.0)

+104.7 +98.6 +121.9

(201) Struct: A3 B1 A1 B3 A2 B2 A∞ B∞

D022 D022 CH, or “40”

−42.3 (−0.2) −41.0 (−0.3) −55.3 (+0.3) 0.0 (0.0)

+85.2 +76.8 +109.6

−59.7 −40.8 −40.0 −29.2 −27.9 −28.8

0.0 (0.0)

+100.5 +90.8 +75.0 +79.8 +56.9 +77.8

(+0.4) (−0.7) (+1.0) (+1.8) (−0.2) (+0.4)

−36.1 +51.0 +40.1 +76.5 +50.0 +136.4

+20.4 (0.0)

(401) Struct: A5 B1 A1 B1

D023

−33.3

(601) Struct: A5 B1 A1 B1

LPS−3

−34.1

Other Struct: A3 B1 A1 B3 A7 B1 A4 B4 A1 B7 A8 B1 A1 B8

L12 L12 D7 D4 D7b Ni8 Nb-type Ni8 Nb-type

Random: A4 B4 A4 B4 A3 B1 A1 B3

SQS8a SQS8b SQS14a SQS14b

−43.4 −44.0 −20.8 −42.9 −20.0

(+0.4) (+0.3) (+0.6) (+1.1) (−0.1)

−48.2 −3.8 −40.8 +10.6 −28.2 −6.7

−33.6 (0.0)

+84.8 +76.0 +61.9

+84.8 (−1.4) +76.0 (+1.8) +61.9 (−3.1)

−37.3 −17.3 +6.8

−37.3 (−0.1) −17.3 (−0.8) +6.8 (−8.3)

+77.5 +78.9 +82.9

+77.5 (−2.7) +78.9 (−0.2) +82.9 (−15.8)

+47.1 +63.7 +42.7

+47.1 (−3.3) +47.7 (+0.4) +36.4 (−1.7)

+12.9 +9.3 +30.9

+12.9 (+1.9) −9.1 (−4.5) +18.2 (+13.3)

+56.8

+56.8 (−0.7)

+183.2 +118.2

+122.6 (+1.2) +97.5 (−9.7) +96.8 (+15.3) +59.8 (−15.3)

−42.5 (+0.2) −43.6 (−0.2) +116.2 +92.2

+77.3 (+7.0) +69.7 (−7.0)

6

+56.5 +37.8

+12.9 −15.2 +5.5 −5.2

(+5.7) (−5.7) (+7.7) (−7.7)

Elastic deformation energy

The total energies of the ordered structures and endpoint constituents are obtained keeping all computational parameters exactly equal. Specifically, we always use the same basis sets (RKmax = 9), charge density cutoffs (RKmax = 19), muffin-tin radii RAu = 2.4a0 , RAg = RCu = RNi = 2.2a0 , maximum difference in the angular momenta in the nonspherical Hamiltonian terms (lmax = 4), maximum angular momenta in the nonspherical charge densities and potentials inside the muffin-tin spheres (lmax = 8), and equivalent k point sets93 in the evaluation of Brillouin zone integrals. When the unit cell vectors of the ordered compound permit, we choose a k mesh equivalent to the 60 special points 8 × 8 × 8 fcc mesh. Several structures (e.g., those of A2 B or AB2 stoichiometry) have reciprocal unit cell vectors which are incommensurate with the 8 × 8 × 8 mesh. In these cases we calculate the total energies of the compounds and fcc constituents with a finer k point grid. This procedure ensures that, due to systematic cancellation of errors, the formation enthalpies ∆H(σ), Eq. (4), converge much faster than the total energies. Indeed, the tests for CuAu described in Sec. IV A 1 show that with our choice of parameters ∆H(σ) are converged to within 2 meV/atom. The atomic positions are relaxed using quantum mechanical forces94 obtained at the end of the selfconsistency iterations. Minimization of the total energy with respect to the cell-external degrees of freedom is done by distorting the shape of the unit cell and tracing the decrease in the total energy. We estimate that the formation enthalpies are converged to at least 5 meV/atom with respect to all relaxational degrees of freedom. Table II and its caption defines our small-unit-cell ordered structures. Many are actually superlattices along (100), (110), (111), (201) and (311) directions. Table III gives the calculated LDA formation energies [Eq. (4)] for these Au-Ag, Cu-Au, Cu-Ag and Ni-Au compounds.

∆Ebulk(as)

aeq

∆Eepi(as)

Substrate lattice constant as FIG. 1. A schematic illustration of the concept of the epib given by the ratio of the bulk taxial softening function q(G), (upper curve) and epitaxial (lower curve) deformation enerb is the ratio of the gies. In the harmonic approximation q(G) curvatures of these curves at the equilibrium point.

energies required to deform bulk A and B to the epitaxial geometry with a planar lattice constant as : h eq b b = min x∆E epi (as , G) ∆ECS (x, G) A as i b . + (1 − x)∆EBepi (as , G) (14)

b is the strain energy of the material Here ∆E epi (as , G) epitaxially stretched to the lattice constant as in the dib and then allowed to relax along rection orthogonal to G, epi b b G. ∆E (as , G) is related to the bulk equation of state b ∆E bulk (as ) via the epitaxial softening function q(as , G): epi b b ≡ ∆E (as , G) , q(as , G) bulk ∆E (as )

B. The constituent strain energy

(15)

bulk where ∆EA (as ) is the energy required to hydrostatically deform a bulk solid to the lattice constant as . Figure 1 illustrates the concept of epitaxial softening:97 when the bulk solid is deformed hydrostatically from aeq to as 6= aeq , its energy rises. Energy can then be lowered if we keep ax = ay = as but relax the third lattice vector b measures the relative to its equilibrium value. q(as , G) energy lowering. Figure 2 shows the calculated LDA q’s for Cu, obtained by minimizing the total energy with respect to the latb for each value of the subtice constant c parallel to G strate lattice parameter as . As explained in Ref. 96, this treatment neglects the so-called shear strain terms, but is exact for the high symmetry directions (100), (111) and b is seen to be a nontrivial (110). The calculated qCu (as , G) function of the substrate lattice parameter as and direcb In contrast, the harmonic elasticity theory,97–100 tion G.

It is well known66 that real-space cluster expansions with finite-ranged interactions incorrectly predict zero formation enthalpies per atom for coherent long-period Ap Bq superlattices, while the correct answers are nonb The zero and depend on the superlattice direction G. constituent strain energy term ∆ECS (σ) in Eq. (10) is specifically designed to reproduce these superlattice energies, which are calculated directly from the LDA as follows. In the long-period limit pq → ∞ the interfacial energy becomes negligibly small (O(1/p)) in comparison with the elastic strain energy needed to deform the constituents to a common in-plane lattice constant as .55,96 Therefore, the formation energy per atom of A∞ B∞ sub with composition x is given by the conperlattice along G b defined as the equilibstituent strain energy ∆ECS (x, G), rium (eq) value of the composition-weighted sum of the 7

Epitaxial parameters: Cu

Epitaxial softening qCu

0.5 0.4 0.3

b = γ(as , G)

(110)

fcc Cu (111)

(201)

(110)

(111)

(100)

0.1

a

(100)

eq

0.0 3.4 3.6 3.8 4.0 Substrate lattice constant as (Å) b of fcc Cu for principle directions as functions FIG. 2. q(G) of the substrate lattice parameter as . Directly calculated LDA values are represented by open symbols, lines show the b in Kubic harmonics. fit using the expansion of γ(G) routinely used for semiconductor systems,97,100,101 gives q’s which do not depend on as : b =1− qharm (G)

B b C11 + ∆ γharm (G)

,

(16)

b is a geometric function of the spherical where γharm (G) b angles formed by G: γharm (φ, θ) = sin2 (2θ) + sin4 (θ) sin2 2 4√ = 4π[K0 (φ, θ) − √ K4 (φ, θ)], 5 21

l=0

b bl (as ) Kl (G),

(18)

which has the property that in the harmonic limit (as → a0 ) all expansion coefficients with angular momenta higher than 4 tend to zero, reproducing γharm from Eq. (17). Due to the cubic symmetry, only terms with l = 0, 4, 6, 8, 10, 12, . . . enter in this expansion. Detailed discussion of the nonlinear epitaxial strain properties of elemental metals will be given in a separate publication.96 eq b is calcuThe constituent strain energy ∆ECS (x, G) lated numerically from Eq. (14) using the direct LDA b for six principle directions. The values of ∆E epi (as , G) eq obtained ∆ECS for these directions are shown in Fig. 3, illustrating several properties of the constituent strain which cannot be reproduced by the harmonic theory.65 First, the curves in Fig. 3 are skewed to different sides, eq while the harmonic ∆ECS must be all skewed to the same eq side. Second, the calculated ∆ECS cross for different directions, a property not allowed by the harmonic functional form. These crossings lead to (201) as the softest direction below x ≈ 0.2, and (110) as the hardest for Au-rich superlattices, while the harmonic theory gives eq eq ∆ECS (111) as the highest and ∆ECS (100) as the lowest constituent strain for all compositions of the studied noeq ble metal alloys. The behavior of ∆ECS for (100) is particularly interesting, since the curves in Fig. 3 abruptly change slope around x ≈ 0.15 and have very low values for x > 14 . As we show in Ref. 96, this is a manifestation of the low energy cost of deforming fcc Cu into the bodycentered tetragonal structure along the epitaxial Baines path. Small constituent strain of (100) superlattices has profound influence on the predicted ground states of CuAu (see Sec. IV A 1). The constituent strain energy for arbitrary direction b is then obtained by interpolating between the prinG ciple directions using the following expansion in Kubic harmonics:

(201)

0.2

lX max

(17)

and Kl are the Kubic harmonics of angular momentum l. Figure 2 shows that the harmonic approximation manifestly breaks down for large epitaxial strains in metals since there are several important qualitative differences between the behavior in Fig. 2 and that predicted b strongly deby the harmonic elasticity. First, q(as , G) pends on the substrate lattice constant, while the harb does not. Second, the harmonic expresmonic qharm (G) b as a function of the sion gives a definite order of q(G) direction, i.e., either (100) is the softest and then (111) must be the hardest, or vice versa. This order does not hold for large deformations. For instance, (201) becomes the softest direction for as ≪ a0 and (110) is the hardest for as ≫ a0 in Cu. Finally, q(100) exhibits a particularly dramatic softening for as ≫ a0 , which has important consequences for the constituent strain energy and stability of superlattices along this direction.96 The above mentioned properties of qCu can be described by generalizing Eq. (17) for γ to higher Kubic harmonics and strain-dependent expansion coefficients:

b = ∆ECS (x, G)

lX max l=0

b cl (x) Kl (G).

(19)

We have taken lmax = 10, which gives five compositiondependent fitting coefficients determined from a fit to the directly calculated values [Eq. (14)] for six principle directions. The characteristic errors of this fit at the equiatomic composition are 1 − 2 meV/atom. Equation (19) is then used in Eqs. (11)–(12). C. Constructing the Cluster Expansion

Once we have a closed-form expresion for the equilibrium constituent strain energy ∆ECS (σ) and a set {∆H LDA (σ)} of T = 0 formation enthalpies, we determine the unknown cluster interactions of Eq. (10) in the following two-step process: 8

First, the total energies of all structures from Table III are used in the fit to investigate the behavior of the rootmean-square (rms) error ∆rms of the fit, Eq. (13), as a function of the number of real-space pair and multibody interactions. Reciprocal space CE allows to add pair interactions systematically in the order of increasing intersite separation, up to any number of near-neighbor shells. The k-space smoothness criterion in Eq. (13) automatically selects optimally short-ranged interactions and chooses physically important pair interactions which are essential to produce a good fit to the directly calculated LDA energies. The dependence of the rms error on the number of pair and multibody interactions is shown in Fig. 4. Figure 4(a) is obtained by fixing the number of multibody interactions, and varying the number of pair interactions. It shows that in all systems the cluster expansion is well converged using 10 to 20 pair interactions. The convergence rate is fastest for Ag-Au and slowest for Ni-Au, which we attribute to increasing size mismatch going from Ag-Au to Ni-Au, with Cu-Ag and Cu-Au exhibiting intermediate convergence rates. Selection of important multibody interactions is more delicate. The number of pair interactions is fixed to a converged value (20 or more), and a large set of 3- to 4body figures is tested as to whether it improves the rms error of the overall fit. It is retained in the CE only if ∆rms decreases considerably. During the fitting process, we also monitor the overall stability of the CE, as measured by a change in other multibody interactions upon the addition of a particular figure. Unstable behavior usually signals of linear dependencies in the chosen set of clusters and an ill-conditioned inverse problem, necessitating a different choice of {Jf }. Figure 4(b) shows the convergence of the CE with respect to the number of multibody interactions, keeping Npairs equal to their converged values. An important thing to notice is that the multibody interactions produce a decrease in the rms error which is of the same magnitude as that due to the pair interactions. Furthermore, the effect of multibody interactions is largest in Ni-Au, and decreases in order of decreasing size mismatch, becoming negligible in Ag-Au. In the second step we test the stability of the fit and its predictive power. Using the trial set of figures obtained in the previous step, we exclude several structures which are fit rather well (e.g., Z2, β2, and L12 in Ni-Au), and repeat the fit, obtaining new values of the effective cluster interactions. These values are used to predict the total energies of the structures excluded from the fit. If the change in ∆HCE (σ) is not acceptable (more than few meV/atom), we return to the first step to search for a better set of interactions. The most severe test is to exclude structures with the poorest fit to their formation enthalpies, e.g., SQS14a and SQS14b in Ni-Au. If the predicted formation energy does not change significantly, the chosen set of figures is considered to be stable and predictive. The final cluster expansion is produced by using this set of figures and all structures from Table III. Figure 5 shows the calculated pair interactions as func-

Constituent strain energies 100

CuAg

(111)

80

(221)

(110)

60 (311)

(201)

Equilibrium constituent strain energies ∆ECS (meV/atom)

40 (100) 20 0 0.0

0.2

0.4

0.6

0.8

1.0

Ag

Cu 200 (111)

NiAu

(221) 150 (311)

(110)

100 (201) 50

0 0.0

(100)

0.2

0.4

0.6

0.8

1.0

Au

Ni 120

CuAu

(111)

100

(221)

80 (311)

(110)

60 (201) 40 (100) 20 0 0.0

Cu

0.2

0.4

0.6

Composition x

0.8

1.0

Au

FIG. 3. Equilibrium constituent strain energies for Cu-Au, Ni-Au and Cu-Ag. The constituent strain energy of Ag-Au is negligibly small and therefore not shown.

9

RMS error of the CE (meV/atom)

Convergence of Cluster Expansion 15

15

10

10 NiAu

NiAu

5

5

CuAu AgAu CuAg

0 0

0

CuAu CuAg AgAu 0 1 2 3 4 5 6 7 Number of multibody interactions

10 20 30 40 50 60 Number of pair interactions

FIG. 4. Root-mean-square errors ∆rms of the cluster expansions for Ag-Au, Cu-Ag, Cu-Au and Ni-Au as functions of the number of pair and multibody interactions.

tion of the near-neighbor fcc shell. There are several noteworthy trends in the four alloy systems: (i) Only in Ag-Au and Cu-Au are the nearest-neighbor pair interactions dominant: in Cu-Ag the 1-st and 3-rd neighbor pair interactions are of similar magnitude, while the 3-rd neighbor interaction dominates in Ni-Au. (ii) The dominant interactions have signs consistent with the observed phase diagrams: Ag-Au and CuAu have positive (“antiferromagnetic”) nearest-neighbor pair interactions J2 , corresponding to the tendency towards complete miscibility and ordering at low temperatures. The behavior of Ni-Au, in spite of positive 1-st and 2-nd neighbor pair interactions, is dominated by the “ferromagnetic” 3-rd neighbor interaction L2 (which causes phase separation at low temperatures). Both dominant 1-st and 3-rd neighbor pair interactions in Cu-Ag are negative, implying a miscibility gap. The constituent strain eq energy ∆ECS is always positive and therefore increases the propensity for incoherent phase separation. (iii) Although the nearest-neighbor pair interaction is clearly dominant in Cu-Au, other pair interactions show a long-ranged oscillatory behavior extending over approximately 15 shells. As found in other systems,65,85 this is a direct consequence of the atomic relaxation caused by the constituent size mismatch between Cu and Au. The pair interactions are slowly decaying in Cu-Ag and Ni-Au, too. The calculated multibody interaction energies are shown in Figure 6. J1 is the point interaction, J3 , K3 , N3 , ..., are triplets and J4 , K4 , and L4 are four-point clusters in increasing order of interatomic separation (see Lu et al.54 for a full description of the clusters). Figure 6 illustrates the importance of the multibody terms in our Hamiltonian.

D. Finding the T = 0 ground states and T > 0 properties

Having parametrized the configurational energies in terms of the mixed-space cluster expansion Eq. (10), we can use it with established statistical methods to predict various structural properties: T = 0 ground states, order-disorder transition temperatures, configurational entropies, free energies, phase stabilities and atomic short-range order parameters. Due to the presence of both reciprocal and real space terms in the Hamiltonian (10), traditional techniques, e.g., the Cluster Variation Method, are not readily applicable. Monte Carlo simulations must be used instead to calculate statistical properties at finite temperatures. The basic computational algorithm is as follows. We adopt the Metropolis algorithm in the canonical ensemble (fixed composition). For each attempted spin flip, the change in the multiplet interaction energy is evaluated in the real space. To obtain the reciprocal space energy (constituent strain and pair interaction energies), the Fourier transform of the spin function S(Ri , σ) is needed. It can be calculated either with the help of the Fast Fourier Transform (FFT) or evaluated directly taking advantage of the special method described in Ref. 87, which is much more economical: if the total number of sites in the simulation box is N , a full FFT has to be done only once after approximately every √ N accepted spin flips, which makes the whole computational effort for this special method scale as N 1.5 . A simulation box of N = 4096 atoms (16 × 16 × 16) is used to calculate all thermodynamic properties presented in this paper. The transition temperatures are computed by cooling the system from high temperatures and monitoring the discontinuities in the average energy and peaks in heat capacity. To eliminate possible hysteresis effects, the resulting low-temperature configurations are gradually heated up past the transition point. The former

10

Pair Interactions

Multibody Interactions

60 40

40

AgAu

AFM

AgAu

20

20

J0 J1 J3 0

0

-20

-20

L3 K3

FM

-40

J4

-40

40

CuAg

AFM

40

Multibody interactions x Df (meV/atom)

Real-space pair interactions x Df (meV/atom)

60

20 0 -20

FM

-40 60 40

NiAu

AFM

20 0 -20 FM

-40

J0

CuAg

20 0

K3 L3

-20

Q3

-40

40

NiAu

P3 Q3

20

N3

J0

J4 K4

0

K3 -20

R3 L4

J1

-40

60

N3

J1 J3

J3

CuAu

40

AFM

40

20 20

0 0

-20

FM

CuAu

P3 Q3 J4

J1 J0

K3

-20

L4

-40 -40

0

5

10

15

20

J3

-60

Number of fcc shells

FIG. 6. Multibody interactions for the studied noble metal alloy systems.

FIG. 5. Real space pair interactions for the studied noble metal alloy systems.

11

process provides the lower bound on the transition temperature, T1 , while the latter gives the upper bound, T2 . The heating and cooling rates are such that T1 and T2 differ by no more than 20K, an insignificant uncertainty compared to the inaccuracies of the LDA calculations and the fit errors of the cluster expansion. 1000 flips/site and a temperature decrease of 2% for each Monte Carlo step are usually sufficient, although in a few cases the results are checked using 2000 flips/site and 0.5% temperature change. Zero temperature ground states are found by cooling the system to T = 0 and checking whether the energy of the final configuration lies on the convex hull. This process is repeated for several random number seeds and starting temperatures, always yielding configurations with similar (usually identical) energies. We explore many equally spaced compositions with an interval ∆x = 0.05. The number of possible configurations for N! each x is Nconf = (xN )!(N (1−x))! . Configurational entropy of the disordered alloys at finite T is computed from the energy vs. temperature curves obtained by cooling the system from very high (“T = ∞”) temperatures. The following thermodynamic formula gives the configurational entropy at temperature T: Z β ∆Sconf (T ) = ∆Sideal + E(T )/T − kB E(β ′ ) dβ ′ ,

Ground state structures 0

CuAu L12

-20

Formation energies ∆E (meV/atom)

Z3

-40

-60

0

0.2

0.4

0.6

0.8

Cu

1

Au

0

AgAu -20 -40 L12

L12 L10

0

0

Ag

where β = 1/kB T and ∆Sideal = kB [x log x + (1 − x) log(1 − x)] is the configurational entropy of an ideal solid solution.

β2 L10

-60

(20)

(100) SL

L12

0.2

0.4 0.6 0.8 Composition x

1

Au

FIG. 7. T = 0 K ground state lines for Cu-Au and Ag-Au obtained from simulated annealing calculations. L12 CuAu3 is not only above the ground state line, but also has a higher formation enthalpy than other structures at the same composition, e.g. LDA calculations place the formation enthalpy of Z3 below that of L12 . Plots for Cu-Ag and Ni-Au are not shown since these systems phase separate at T = 0 K.

IV. RESULTS A. T = 0 Ground States 1. Ground states of Cu-Au

calculated LDA enthalpy of formation of Z3 (which is a Cu1 Au3 (001) superlattice) is considerably lower than that of L12 CuAu3 . We carefully checked whether the predicted new LDA ground states for Au-rich Cu-Au alloys artifacts of some approximation in our LDA calculations or the fit error of the cluster expansion. The latter possibility was quickly dismissed, since the directly calculated LDA enthalpies of formation for L10 , β2, L12 and Z3 agreed with the values derived from the cluster expansion to better than 2 meV/atom (see Table III), while the new (100) SL ground state is 14 meV/atom below L12 . To address the former possibility, we performed careful convergence tests for L10 , β2, L12 and Z3 with respect to the plane wave cutoff and number of k points in the first Brillouin zone. The cutoff was increased from RKmax = 9 to RKmax = 11 and the density of the Brillouin zone

Figure 7 shows the calculated T = 0 ground state lines of Cu-Au and Ag-Au which were obtained from simulated annealing quenches of a 16 × 16 × 16 system. In Cu-Au, we find the L12 (Cu3 Au) and L10 (CuAu) structures as the stable ground states of Cu-rich alloys, in agreement with the existing phase diagram data.1–4 These data also list L12 as the stable low-temperature phase of CuAu3 . However, we find new, previously unsuspected ground states of Au-rich compounds, all belonging to the family of (001) superlattices. At x = 23 we find a stable β2 (CuAu2 ) phase (prototype MoSi2 ), which is a Cu1 Au2 superlattice along (001). At x = 43 , our cluster expansion predicts that a complex Cu1 Au4 Cu1 Au4 Cu1 Au2 Cu1 Au2 (001) superlattice falls on the convex hull, although its energy is less than 2 meV below the tieline connecting β2 (CuAu2 ) and Au. Furthermore, even the directly 12

mesh was doubled from 8 × 8 × 8 to 16 × 16 × 16, an eightfold increase in the total number of k points. These tests showed that the formation enthalpies of L10 , β2 and L12 were converged to within 1 meV/atom with respect to the size of the basis set and the number of k points. Further, we checked how the choice of muffintin radii affected ∆H. Varying RMT (Au) between 2.3a0 and 2.5a0 changed the formation enthalpies by at most 2 meV/atom and did not shift the relative stabilities of phases. Finally, we repeated these calculations using the Perdew-Zunger90 parametrization of the CeperleyAlder91 LDA functional, as well as the generalized gradient approximation (GGA) of Perdew and Wang,92 and found insignificant (about 2 meV/atom) changes in the formation enthalpies. Inclusion of the spin-orbit interaction in the second variation procedure104 changed the formation enthalpy of L10 (CuAu) by only 3.7 meV/atom (from −48.2 to −51.9), indicating that it is not important for the energetics of Cu-Au. This conclusion is in line with the findings of Ref. 105 that the spinorbit interaction influences the band structure but has little effect on equilibrium lattice properties. Therefore, we conclude that state-of-the-art first-principles density functional calculations do not predict L12 to be a stable T = 0 ground state of CuAu3 . It is possible that van der Waals interactions, omitted by the LDA and important for large, polarizable atoms such as Au, can affect the formation energies and hence the ground states of Cu-Au. We next analyze the possibility that the correct T = 0 ground state around x = 34 is not L12 as has been assumed in the literature before. Although most compilations1–4 of binary alloy phase diagrams give L12 as the stable structure of CuAu3 , the experimental evidence7,8,10 seems inconclusive because of the difficulties in obtaining equilibrated long-range ordered samples. X-ray studies8 have found superlattice peaks consistent with the cubic L12 structure, but only very broad loworder reflections have been observed. These superlattice lines could not be sharpened by any heat treatment.8 It is not clear to us if the X-ray reflections can be reindexed according to some other non-L12 phase. It is also possible that at elevated (T ≈ 500 K) temperatures L12 is stabilized by the entropy (configurational and vibrational), while another transformation to the low-energy structure should occur but is kinetically inhibited below 500 K. The biggest experimental obstacles to verifying our predictions seem to be low diffusion rates below the ordering temperature of CuAu3 , Tc ≈ 500 K. Next we discuss the experimental signatures of the new LDA ground state structures. MoSi2 -type β2 CuAu2 has a superlattice reflection at ( 23 00), but CuAu3 (100) superlattice has reflections at (100) and ( 13 00). These reflections also manifest themselves in the predicted atomic short-range order of the disordered alloys (for details see Ref. 103).

2. Ground states of Ag-Au, Cu-Ag and Ni-Au

The ground state line of Ag-Au is shown in Figure 7(b), exhibiting L12 (Ag3 Au), L10 (AgAu) and L12 (AgAu3 ) stable low-temperature phases. Experimentally, these alloys are known to be completely miscible,2–4 and there are several indications69 that they would order below 200 K if not for the very low diffusion rates. Theoretical transition temperatures and short-range order patterns, as well as a complete discussion are given by Lu and Zunger.54 The calculated ground states of Cu-Ag and Ni-Au are found to be phase separation, in agreement with the experimental enthalpy data.2 Neither alloy has a single ordered or disordered structure with negative enthalpy of formation and therefore there are no stable T = 0 ground states except the phase-separated alloy.

B. Mixing enthalpies

It is interesting to compare the calculated mixing enthalpies of disordered Cu-Au alloys with the available theoretical and experimental data. Table IV summarizes the values of ∆Hmix (x, T ) for the completely random (T = ∞), short-range ordered (T = 800 K) and completely ordered (T = 0 K) Cu-Au alloys at compositions x = 41 , 12 and 43 . Several important points are apparent from this table: (i) Studies50,48,62 which have completely neglected atomic relaxations predict a substantially positive enthalpy of formation for the completely random alloy. In our calculations, relaxations in the random alloy reduce ∆Hmix (T = ∞) by a large amount, bringing it down to essentially zero. (ii) Comparison of the present results for the T = ∞ random alloys with those of Wei et al.51 shows the influence of the number of structures included in the cluster expansion. Since Wei et al. used the same FLAPW method88 , but included a set of only five high-symmetry ordered structures [A1 (Cu), L12 (Cu3 Au), L10 (CuAu), L12 (CuAu3 ) and A1 (Au)], the atomic relaxation effects were included incompletely. Indeed, their treatment gives much larger mixing enthalpies of the random Cu-Au alloys than the present work employing approximately 30 low-symmetry structures with large relaxations. Therefore we conclude that the Connolly-Williams set of five ordered structures cannot correctly capture the large decrease of the mixing enthalpy of random Cu-Au alloys caused by the atomic relaxations. (iii) The good agreement between the relaxed (this study) and “unrelaxed” (Wei et al.51 ) values of ∆Hmix at T = 800 K suggests that the short-range order in CuAu tends to decrease the role of the atomic relaxations. This effect can be qualitatively explained on the basis of the ordering tendency towards high-symmetry structures

13

TABLE IV. Calculated mixing enthalpies of disordered Cu1−x Aux alloys compared with the values obtained by other studies and experimental measurements (in meV/atom). FLAPW is the full-potential linearized augmented plane wave method, LMTO – linearized muffin-tin-orbitals method, KKR – Korringa-Kohn-Rostoker multiple scattering method, ASA – atomic-sphere approximation, CPA – coherent potential approximation, CWM - Connolly-Williams cluster expansion, MSCE – mixed-space cluster expansion used in this study, “Rel.” – incorporating atomic relaxations, and “Unrel.” – neglecting atomic relaxations. Composition

Expt.f

This study FLAPW MSCE (Rel.)

Wei et al.a FLAPW CWM (Rel.)

Amador et al.b LMTO-ASA CWM (Unrel.)

Terakura et al.c ASW CWM (Unrel.)

Ruban et al.d LMTO-ASA CPA (Unrel.)

Weinberger et al.e KKR-ASA CPA (Unrel.)

+2.6 +1.6 +5.4

+46.3 +38.0 +18.6

+59 +61 +39

+26.9 +30.4 +20.4

+54.6 +44.3 +19.8

−27 −57 −31

−46g −53g −31g

−17.3 −19.3 −1.2

−16.9 −2.6

−6 −5 +8

−74 −91 −59

−37.3 −48.2 −17.3

−36.0 −62.9 −26.4

−65.0 −69.7 −34.0

−60.7 −83.4 −56.1

−54 −76 −47

∆Hmix (T = ∞) Cu0.75 Au0.25 Cu0.50 Au0.50 Cu0.25 Au0.75 ∆Hmix (T = 800 K) Cu0.75 Au0.25 Cu0.50 Au0.50 Cu0.25 Au0.75 ∆Hmix (T = 0 K) L12 Cu3 Au L10 CuAu L12 CuAu3 a

Ref. 51 using the Connolly-Williams structures (relaxation of L10 only). Ref. 50. c Ref. 48. d Ref. 62. e Ref. 58. f Ref. 2. g Values obtained at T = 720K. b

loys.

which have little or no relaxation energy (L12 and L10 in Cu-rich alloys). (iv) The mixing enthalpies of the random alloy calculated by Weinberger et al.58 using the coherent-potential approximation (CPA) differ strongly not only from those obtained using the cluster expansion methods,51,50,48 but also from the numbers given in the CPA work of Ruban, Abrikosov, and Skriver.62 Since the CPA of Weinberger et al.58 neglects the (a) atomic relaxation, (b) charge transfer and (c) short-range order, which all lower the formation energies, the negative values obtained by Weinberger et al.58 are very puzzling. (v) There are significant discrepancies between the best calculated and experimentally measured15,14,2 values of ∆Hmix at both T = 0 K and T = 800 K. At present these discrepancies are hard to explain since the available general potential LDA calculations51,52,57 of ∆H(L12 ) and ∆H(L10 ) agree with each other reasonably well. On the other hand, formation energies in Cu-Au are numerically very small and present a severe test for any firstprinciples model of electronic exchange-correlation. It is noteworthy that several less accurate first-principles calculations, using the atomic-sphere approximation (ASA), have achieved better agreement with the experimental enthalpies of formation than the state of the art general potential techniques. We consider this to be fortuitous. In all cases, LDA calculations correctly predict the relative magnitudes of ∆H for L12 and L10 , as well as reproduce measured asymmetry in formation enthalpies towards more negative values of ∆Hmix for Cu-rich al-

C. Order-disorder transition temperatures

Order-disorder transitions have been investigated at compositions (x = 14 , 21 , 23 and 34 ) using the Monte Carlo simulation technique described in Sec. III D. The resulting transition temperatures, Tc , are given in Table V. All transitions are found to be first order, involving discontinuities in the energy and correlation functions. At x = 41 we find a transition from the disordered state to long-range ordered L12 Cu3 Au at Tc = 530 K, which is only 130 K lower than the experimentally observed transition temperature. For the equiatomic alloy at x = 12 the calculated and experimental transition temperatures agree to a few degrees Kelvin. However, we do not find the CuAu II phase which exists in a narrow temperature range between 658 K and 683 K. This phase is stabilized by the free energy differences between L10 and long-period superstructures of L10 which are as small as 1 meV/atom56 and therefore beyond the accuracy of selfconsistent LDA calculations. For x = 43 we obtain a sequence of transformations, the first one occuring at T = 750 K from the disordered A1 phase to a coherent two-phase mixture of β2 and A1. Then a subsequent transition at T = 635 K takes CuAu3 into the long-range ordered (100) superlattice which is predicted to be the stable T = 0 ground state at that composition (see Sec. IV A 1). The calculated transition 14

form calc TABLE VI. The experimentally measured2 entropy of formation ∆Stot , the calculated configurationl entropy ∆Sconf and form the derived non-configurational entropy of formation, ∆Snon−conf . All values are given in units of kB /atom.

System

x

Cu-Au Ag-Au Cu-Ag Ni-Au a

form ∆Stot

T (K)

0.5 0.5 0.141 0.5

800 800 1052 1100

∆Sideal

0.73 0.52 0.77 1.04

0.69 0.69 0.41 0.69

calc ∆Sconf

0.57 0.62 0.40a 0.56

form ∆Snon−conf = form calc ∆Stot − ∆Sconf 0.16 −0.10 0.37 0.48

This value was obtained at T = 1136 K, since a coherent phase separation starts at lower temperatures.

TABLE V. Calculated order-disorder transition temperatures (in K) for Cu-Au. A1 denotes the configurationally disordered fcc phase, and n/a means that the transition has not been observed (either experimentally or in the Monte Carlo simulation). Composition Cu0.75 Au0.25 Cu0.50 Au0.50 Cu0.33 Au0.66 Cu0.25 Au0.75

Transition A1 → L12 A1 → L10 A1 → β2 A1 → L12 A1 → β2 + A1 β2 + A1 → (100)SL

Expt. 663 683/658a n/a ≈ 500 n/a n/a

form ∆Stot (A1−x Bx , T ) because all other contributions cancel out in Eq. (21). The non-configurational entropy of formation, form form ∆Snon−conf (A1−x Bx , T ) = ∆Stot (A1−x Bx , T ) − ∆Sconf (A1−x Bx , T ),

This study 530 660 735 n/a 750 680

contributes to such important quantities as mutual solubility limits and miscibility gap temperatures. Noble metal alloys are excellent cases to test the valform ues of ∆Snon−conf since accurate experimental data on form the entropies of formation, ∆Stot , are available, and the configurational entropy ∆Sconf can be calculated accurately using the thermodynamic integration technique described in Sec. III D. Table VI gives the measured entropies of formation for disordered solid solutions form A1−x Bx , ∆Stot (x, T ), the maximum attainable configurational entropy ∆Sideal , as well as the theoretically calcalc culated configurational entropy ∆Sconf , and the derived value for the non-configurational entropy of formation, form ∆Snon−conf . It shows that the size-mismatched noble form metal systems have large amounts of ∆Snon−conf in the disordered solid solution. Since it is unlikely that these form values of ∆Snon−conf are of electronic or magnetic origin, we suggest that the excess entropy in the disordered solid solutions of Ni-Au, Cu-Ag and Cu-Au is vibrational. It is possible that the atomic relaxations lead to a softening of lattice vibrations, although the physical mechanism of this softening is unclear at present. Sanchez et al.49 in their study of the Cu-Ag system noted that even a very crude model of the vibrational entropy markedly improved the agreement with the experimental solubility data. In the case of Ni-Au, which form exhibits the largest ∆Snon−conf , it is possible to reconcile the experimentally measured and theoretically calculated miscibility gap temperatures only by taking into account the non-configurational entropy of formation.117 form The fact that Cu-Au also has a positive ∆Snon−conf has little qualitative effect on the phase diagram since Cu and Au are completely miscible from total energy and configurational entropy considerations alone. Ag-Au is form calculated to have a negative ∆Snon−conf , but its value is close to the experimental uncertainty in the measurement of ∆S.

a

CuAu undergoes a transition to CuAu-II at 683 K, subsequently transforming into L10 CuAu-I at 658 K.

at x = 23 goes straight into the β2 phase at T = 735 K. Therefore, a two-phase β2 + A1 field is predicted to exist at temperatures somewhere between 635 K and 730 K and around x = 34 . These predictions reflect the LDA. As stated in Sec. IV A 1, corrections to the LDA might be significant. D. Non-configurational entropy

The effect of the non-configurational entropy (electronic, vibrational, etc.) on the alloy phase stability has recently attracted considerable interest.106–116 For instance, it has been suggested108–115 that there are large differences in the vibrational entropies of ordervib vib ing Sordered − Sdisord , which should manifest themselves in shifts of the order-disorder transition temperatures. There is another important class of thermodynamic properties where the vibrational entropy may play a role, and which has often been overlooked. Namely, it is the entropy of formation with respect to the pure constituents, defined in analogy with ∆H in Eq. (4): form ∆Stot (A1−x Bx , T ) = S(A1−x Bx , T ) − (1 − x)S(A, T ) − xS(B, T ),

(22)

(21)

where S(A, T ) is the total entropy of the pure constituent A at temperature T . It is often assumed that the configurational entropy is the dominant contribution to 15

have been fully relaxed to minimize the total energy. The results for Cu-Au and Ni-Au interatomic bond lengths are shown in Fig. 8. The main features are: (i) In spite of the different phase diagram properties (Ni-Au phase separates and Cu-Au orders at T = 0 K), the calculated behavior of bond lengths is very similar, which we attribute to the similar size mismatch in both systems (12% in Cu-Au and 15% in Ni-Au). (ii) Our calculations give three distinct bond lengths at all compositions, which is also observed experimentally.23,24 Probably the most interesting feature in Fig. 8 is the crossing of RBB (x) and RAB (x) curves at x = 43 in both systems. The measurements for Cu-Au24 and Ni-Au23 indicate that this may indeed be correct, since the deduced values around this composition are very close and have large error bars. (iii) Another important feature, observed experimentally and reproduced by our SQS results, is that A − A bonds change much more as x varies from 0 to 1 than B − B bonds when x varies from 1 to 0, suggesting that the compressed bonds become increasingly stiff and the expanded bonds weaken. This behavior can be explained by the asymmetry in the interatomic potential curves, which are rapidly hardening upon compression and softening upon expansion. However, our results for RAA at x = 43 and RBB at x = 14 are obtained from an average of only 4 minority bonds in the SQS14 structures, and perhaps are not representative of a wider statistical sample. (iv) It is interesting to note that the predicted bond lengths between unlike atoms RAB do not follow the linear relation RAB = RAA + x(RBB − RAA ).

Average SQS bond lengths 2.9

Au-Au

CuAu

2.8

Cu-Au

2.7 2.6

Cu-Cu Bond length (Å)

2.5 2.4

0

0.2

0.4

0.6

0.8

Cu 2.9

1

Au Au-Au

NiAu

2.8 2.7

Ni-Au

2.6

Ni-Ni

2.5 2.4 0

Ni

0.2

0.4

0.6

0.8

Composition x

1

V. SUMMARY

Au We have showed that accurate first-principles studies of alloys with large size mismatches are now feasible using the mixed-space cluster expansion method. This method has been applied to noble metal alloys where vast amounts of experimental data and many theoretical studies are available. (i) The mixed-space cluster expansion has been generalized to include the effects of nonlinear strain on the formation energies of long-period superlattices. We find that the elastic energy, required to lattice-match Cu and Ni to (100) surfaces of Au and Ag, is anomalously low, leading to a very low constituent strain energy of (100) superlattices. This effect is partly responsible for the stabilization of new LDA ground states of Au-rich Cu-Au alloys. (ii) In Au-rich Cu-Au, we predict new T = 0 K ground states. Our LDA results place L12 (CuAu3 ), previously thought of as the stable T = 0 state of CuAu3 , higher in energy than a family of superlattices along (100) direction. In particular, MoPt2 -type CuAu2 [Cu1 Au2 superlattice along (100)] and a compli-

FIG. 8. SQS bondlengths for Cu-Au and Ni-Au.

E. Bond lengths in random alloys

Since recent experimental measurements of the composition-dependence of interatomic bond lengths in Cu-Au24 and Ni-Au23 have found several unusual features, it is interesting to address these trends from firstprinciples LDA calculations. In the present work we model the atomic positions in the random alloys using special quasirandom structures118 (SQS). These periodic structures are designed to reproduced the pair and multibody correlation functions of the perfectly disordered configuration as closely as possible. It has been shown118 that even small unit cell SQS’s can give rather accurate representation of the properties of random alloys. We have performed LDA calculations for 8 atom/cell SQS’s at x = 14 (SQS14a), x = 12 (SQS8a , SQS8b ) and x = 43 (SQS14b ). The atomic positions and cell coordinates 16

6

cated Cu1 Au4 Cu1 Au4 Cu1 Au2 Cu1 Au2 (100) superlattice are found to be the LDA ground states. (iii) There are significant discrepancies (up to 50%) between the experimentally measured and calculated LDA mixing enthalpies for Cu-Au alloys. This is surprising since the experimental mixing enthalpies of Ni-Au and Ag-Au are reproduced very well.54,117 (iv) The calculated order-disorder transition temperatures are in an excellent agreement with experiment. For instance, Tccalc (x = 14 ) = 530 K and Tccalc (x = 1 expt (x = 41 ) = 663 K and 2 ) = 660 K, compared with Tc 1 expt Tc (x = 2 ) = 683/658 K. (v) From the experimentally measured entropies of form formation ∆Stot and the calculated configurational calc entropies ∆Sconf , we obtain large non-configurational (probably vibrational) entropies of formation in the sizeform form calc mismatched systems, ∆Snon−conf = ∆Stot − ∆Sconf . These entropies allow one to reconcile the experimental miscibility gap temperature and formation enthalpies of Ni-Au with the theoretical LDA values.117 (vi) Bond length distributions in Ni-Au and CuAu have been studied via supercell calculations employing the special quasirandom structure technique. The important qualitative features of recent EXAFS measurements23,24 are correctly reproduced: existence of distinct A − A, B − B and A − B bond lengths at all compositions, possible crossing of RAA (x) and RAB (x) around x = 34 (where x is the composition of the larger constituent), softening of the shorter bond as x → 1, and deviations of the bond length RAB (x) between unlike atoms from the linear Vergard’s law.

J. M. Cowley, J. Appl. Phys. 21, 24 (1950); S. Ogawa and D. Watanabe, J. Appl. Phys. 22, 1502 (1951). 8 B. W. Batterman, J. Appl. Phys. 28, 556 (1957). 9 G. C. Kuczynski, M. Doyama, and M. E. Fine, J. Appl. Phys. 27, 651 (1956). 10 R. Kubiak and J. Janczak, Journal of Alloys and Compounds 176, 133 (1991). 11 F. M. d’Heurle and P. Gordon, Acta Met. 9, 304 (1961). 12 R. L. Orr, J. Luciat-Labry, and R. Hultgren, Acta Met. 8, 431 (1960). 13 R. A. Oriani, Acta Met. 2, 608 (1954). 14 R. A. Oriani and W. K. Murphy, J. Phys. Chem 62, 327 (1958). 15 R. L. Orr, Acta Met. 8, 489 (1960). 16 M. Hirabayashi, S. Nagasaki, and H. K¯ ono, J. Appl. Phys. 28, 1070 (1957). 17 M. Bienzle, T. Oishi, F. Sommer, and K. Ono, Materials Transactions, JIM 33, 51 (1992). 18 M. Bienzle, T. Oishi, and F. Sommer, J. of Alloys and Compounds 220, 182 (1995). 19 T. Claeson and J. B. Boyce, Phys. Rev. B 29, 1551 (1984). 20 G. S. Sohal, C. Carbone, E. Kisker, S. Krummacher, A. Fattah, W. Uselhoff, R. C. Albers, and P. Weinberger, Z. Phys. B - Condensed Matter 78, 295 (1990). 21 T.-U. Nahm, K.-H. Park, S.-J. Oh, S.-M. Chung, and G. K. Wertheim, Phys. Rev. B 52, 16 466 (1995). 22 T. Shinohara, S. Saitoh, F. Wagatuma, and S. Yamaguchi, Phil. Mag. Lett. 74, 153 (1996). 23 G. Renaud, N. Motta, F. Lan¸con, and M. Belakhovsky, Phys. Rev. B 38, 5944 (1988). 24 A. I. Frenkel, E. A. Stern, A. Rubshtein, A. Voronel, and Yu. Rosenberg, submitted to J. de Phys. (1996). 25 W. Shockley, J. Chem. Phys. 6, 130 (1938). 26 J. M. Cowley, Phys. Rev. 77, 669 (1950). 27 C. M. van Baal, Physica 64, 571 (1973). 28 N. S. Golosov, L. E. Popov, and L. Ya. Pudan, J. Phys. Chem. Solids 34, 1149 and 1157 (1973). 29 R. Kikuchi, J. Chem. Phys. 60, 1071 (1974). 30 D. de Fontaine and R. Kikuchi, National Bureau of Standards SP-496, 999 (1977). 31 R. Kikuchi, J. M. Sanchez, D. de Fontaine, and H. Yamauchi, Acta Met. 28, 651 (1980). 32 C. Sigli and J. M. Sanchez, Acta Met. 33, 1097 (1985). 33 R. Kikuchi, in Noble Metal Alloys, (The Metallurgical Society, Warrendale, 1986). 34 W. A. Oates, P. J. Spencer, and S. G. Fries, CALPHAD 20, 481 (1996). 35 D.-H. Wu and R. A. Tahir-Kheli, J. Phys. Soc. Japan 31, 641 (1971). 36 B. Chakraborty and Z. Xi, Phys. Rev. Lett. 68, 2039 (1992). 37 A. Chiolero, Ph.D. thesis, University of Fribourg, Switzerland, unpublished (1995). 38 B. Chakraborty, Europhys. Lett. 30, 531 (1995). 39 G. J. Ackland and V. Vitek, Phys. Rev. B 41, 10 324 (1990). 40 H. M. Polatoglou and G. L. Bleris, Interface Science 2, 31 (1994). 41 G. Bozzolo, J. Ferrante, and J. R. Smith, Phys. Rev. B 7

ACKNOWLEDGMENTS

This work has been supported by the Office of Energy Research, Basic Energy Sciences, Materials Science Division, U.S. Department of Energy, under contract DEAC36-83CH10093.

1

Phase Diagrams of Binary Copper Alloys, eds. P. R. Subramanian, D. J. Chakrabarti, and D. E. Laughlin (ASM International, Materials Park, OH, 1994). 2 R. Hultgren, P. D. Desai, D. T. Hawkins, M. Gleiser, and K. Kelley, Selected Values of the Thermodynamic Properties of Binary Alloys (American Society for Metals, Metals Park, OH, 1973). 3 M. Hansen, Constitution of Binary Alloys (Genium, Schenectady, N.Y., 1985). 4 T. B. Massalski, Binary Alloy Phase Diagrams (ASM International, Materials Park, OH, 1990). 5 Noble Metal Alloys, edited by T. B. Massalski (The Metallurgical Society, Warrendale, 1986).

17

45, 493 (1992). G. Mazzone, V. Rosato, and M. Pintore, Phys. Rev. B 55, 837 (1997). 43 C. Seok and D. W. Oxtoby, J. Phys. Condens. Matter 9, 87 (1997). 44 L. G. Ferreira and M. A. Boselli, Solid State Commun. 96, 745 (1995). 45 L. G. Ferreira and M. A. Boselli, to appear in J. Phys. Cond. Matter (1997). 46 M. Ahlers, Z. Phys. B 99, 491 (1996). 47 F. Cleri and V. Rosato, Phil. Mag. Lett. 67, 369 (1993). 48 K. Terakura, T. Oguchi, T. Mohri, and K. Watanabe, Phys. Rev. B 35, 2169 (1987); T. Mohri, K. Terakura, S. Takizawa, and J. M. Sanchez, Acta Met. 39, 493 (1991). 49 J. M. Sanchez, J. P. Stark, and V. L. Moruzzi, Phys. Rev. B 44, 5411 (1991). 50 C. Amador and G. Bozzolo, Phys. Rev. B 49, 956 (1994). 51 S.-H. Wei, A. A. Mbaye, L. G. Ferreira, and Alex Zunger, Phys. Rev. B 36, 4163 (1987). 52 Z. W. Lu, S.-H. Wei, A. Zunger, S. Frota-Pessoa, and L. G. Ferreira, Phys. Rev. B 44, 512 (1991). 53 Z.-W. Lu and A. Zunger, Phys. Rev. B 50, 6626 (1994). 54 Z. W. Lu, B. M. Klein, and A. Zunger, Modelling Simul. Mater. Sci. Eng. 3, 753 (1995). 55 Z. W. Lu, B. M. Klein, and A. Zunger, Superlattices and Microstructures 18, 161 (1995). 56 A. T. Paxton and H. M. Polatoglou, Phys. Rev. Lett. 78, 270 (1996). 57 J. W. Davenport, R. E. Watson, and M. Weinert, Phys. Rev. B 37, 9985 (1988). 58 P. Weinberger, V. Drchal, L. Szunyogh, J. Fritscher, and B. I. Bennett, Phys. Rev. B 49, 13 366 (1994). 59 D. D. Johnson and F. J. Pinski, Phys. Rev. B 48, 11 553 (1993). 60 J. Kudrnovsk´ y, S. K. Bose, and O. K. Andersen, Phys. Rev. B 43, 4613 (1991). 61 I. A. Abrikosov, Yu. H. Vekilov, and A. V. Ruban, Phys. Lett. A 154, 407 (1991). 62 A. V. Ruban, I. A. Abrikosov, and H. L. Skriver, Phys. Rev. B 51, 12 958 (1995). 63 B. Ginatempo and J. B. Staunton, J. Phys. F: Met. Phys 18, 1827 (1988). 64 P. Hohenberg and W. Kohn, Phys. Rev. 136, 864 (1964); W. Kohn and L. J. Sham, Phys. Rev. A 136, 1133 (1965). 65 D. B. Laks, L. G. Ferreira, S. Froyen, and A. Zunger, Phys. Rev. B 46, 12 587 (1992). 66 A. Zunger, in NATO ASI on Statics and Dynamics of Alloy Phase Transformations, (Plenum Press, New York, 1994), p. 361. 67 W. B. Pearson, Handbook of Lattice Spacings and Structures of Metals and Alloys (Pergamon, Oxford, 1967). 68 L. Pauling, The Nature of the Chemical Bond, 3rd ed. (Cornell University Press, New York, 1960), p. 93. 69 B. Sch¨ onfeld, J. Traube, and G. Kostorz, Phys. Rev. B 45, 613 (1992); N. Norman and B. E. Warren, J. Appl. Phys. 22, 483 (1951); K. Ziesemer, PhD Thesis Universit¨ at Frankfurt, 1976. 70 B. Velick` y, S. Kirkpatrick, and H. Ehrehreich, Phys. Rev. 175, 747 (1968); P. Soven, ibid. 178, 1136 (1969). 71 For a review of the coherent potential approximation, see

B. L. Gyorffy, D. D. Johnson, F. J. Pinski, D. M. Nicholson, and G. M. Stocks, in Alloy Phase Stability, edited by G. M. Stocks and A. Gonis (Kluwer, Dordrecht, 1898), p. 293. 72 F. Ducastelle, Order and Phase Stability in Alloys (NorthHolland, New York, 1991). 73 Z. W. Lu, S.-H. Wei, and A. Zunger, Phys. Rev. B 45, 10 314 (1992). 74 N. Mousseau and M. F. Thorpe, Phys. Rev. B 45, 2015 (1992). 75 N. Papanikolaou, N. Stefanou, R. Zeller, and P. H. Dederichs, Phys. Rev. B, in press (1997); Comp. Mat. Sci., in press (1997); N. Papanikolaou, N. Stefanou, R. Zeller, and P. H. Dederichs, in Stability of Materials NATO ASI, eds. A. Gonis, P. E. A. Turchi, and J. Kudrnovsky (Plenum, New York, 1996). 76 Z. W. Lu, S.-H. Wei, and A. Zunger, Phys. Rev. B 44, 10 470 (1991). 77 Z. W. Lu, S.-H. Wei, and A. Zunger, Phys. Rev. Lett. 66, 1753 (1991); ibid. 68, 1961 (1992). 78 J. W. D. Connolly and A. R. Williams, Phys. Rev. B 27, 5169 (1983). 79 D. de Fontaine, in Solid State Physics, eds. H. Ehrenreich and D. Turnbull, Vol. 51, (Academic Press, New York, 1994). 80 J. M. Sanchez, F. Ducastelle, and D. Gratias, Physica 128A, 334 (1984). 81 M. Asta, C. Wolverton, D. de Fontaine, and H. Dreyss´e, Phys. Rev. B 44, 4907 (1991). 82 J. M. Sanchez, Phys. Rev. B 48, 14 013 (1993). 83 M. Born and K. Huang, Dynamical Theory of Crystal Lattices (Oxford Univ. Press, London and New York, 1961). 84 A. A. Maradudin, E. W. Montroll, G. H. Weiss, and I. P. Ipatova, Theory of Lattice Dynamics in the Harmonic Approximation, in Solid State Physics, eds. H. Ehrenreich and D. Turnbull, Suppl. 3, (Academic Press, New York, 1971). 85 S. de Gironcoli, P. Giannozzi, and S. Baroni, Phys. Rev. Lett. 66, 2116 (1991); M. Peressi and S. Baroni, Phys. Rev. B 49, 7490 (1994). 86 K. Binder, D. W. Heerman, Monte Carlo Simulation in Statistical Physics (Springer-Verlag, Berlin, 1988). 87 Z.-W. Lu, D. B. Laks, S.-H. Wei, and A. Zunger, Phys. Rev. B 50, 6642 (1994). 88 S.-H. Wei and H. Krakauer, Phys. Rev. Lett. 55, 1200 (1985), and references therein. 89 E. Wigner, Phys. Rev. 46, 1002 (1934). 90 J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). 91 D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). 92 J. P. Perdew and Y. Wang, Phys. Rev. B 33, 8800 (1986); 40, 3399 (1989) (E); J. P. Perdew, in Electronic Structure of Solids ’91, eds. P. Ziesche and H. Eschrig (Akademie Verlag, Berlin, 1991), p. 11. 93 S. Froyen, Phys. Rev. B 39, 3168 (1989). 94 R. Yu, D. Singh, and H. Krakauer, Phys. Rev. B 43, 6411 (1991). 95 P. Villars and L. D. Calvert, Pearson’s Handbook of Crystallographic Data for Intermetallic Phases (American Society for Metals, Metals Park, OH, 1985).

42

18

96

V. Ozoli¸ nˇs, C. Wolverton, and A. Zunger, unpublished (1997). 97 A. Zunger, in Hanbook of Crystal Growth, Vol. 3, edited by D. T. J. Hurle (Elsevier, Amsterdam, 1994), p. 997. 98 J. Hornstra and W. J. Bartels, J. Crystal Growth 44, 513 (1978). 99 E. Anastassakis, J. Appl. Phys. 68, 4561 (1990); J. of Crystal Growth 114, 647 (1991). 100 D. J. Bottomley and P. Fons, J. of Crystal Growth 160, 406 (1996). 101 D. M. Wood, J. Vac. Sci. Technol. B 10, 1675 (1992). 102 A. G. Khachaturyan, Theory of Structural Transformations in Solids (Wiley, New York, 1983). 103 C. Wolverton, V. Ozoli¸ nˇs, and A. Zunger, unpublished (1997). 104 A. H. MacDonald, W. E. Pickett, and D. D. Koelling, J. Phys. C 13, 2675 (1980); W. E. Pickett, A. J. Freeman, and D. D. Koelling, Phys. Rev. B 23, 1266 (1981). 105 S.-H. Wei and A. Zunger, Phys. Rev. B 37, 8958 (1988). 106 C. Wolverton and A. Zunger, Phys. Rev. B 52, 8813 (1995). 107 R. E. Watson and M. Weinert, Phys. Rev. B 30, 1641 (1984). 108 L. J. Nagel, L. Anthony, and B. Fultz, Phil. Mag. Lett. 72, 421 (1995). 109 L. Anthony, J. K. Okamoto, and B. Fultz, Phys. Rev. Lett. 70, 1128 (1993). 110 B. Fultz, L. Anthony, L. J. Nagel, R. M. Nicklow, and S. Spooner, Phys. Rev. B 52, 3315 (1995). 111 L. Anthony, L. J. Nagel, J. K. Okamoto, and B. Fultz, Phys. Rev. Lett. 73, 3034 (1994). 112 B. Fultz, L. Anthony, J. L. Robertson, R. M. Nicklow, S. Spooner, and M. Mostoller, Phys. Rev. B 52, 3280 (1995). 113 L. J. Nagel, B. Fultz, J. L. Robertson, and S. Spooner, Phys. Rev. B 55, 2903 (1997). 114 G. D. Garbulsky and G. Ceder, Phys. Rev. B 49, 6327 (1994); ibid. 53, 8993 (1996). 115 J. D. Althoff, D. Morgan, D. de Fontaine, M. Asta, S. M. Foiles, and D. D. Johnson, Phys. Rev. B 56, R5705 (1997). 116 J. Desplat, F. Bley, and F. Livet, Acta Mater. 44, 4961 (1996). 117 C. Wolverton and A. Zunger, Comp. Mater. Sci. 8, 107 (1997). 118 A. Zunger, S.-H. Wei, L. G. Ferreira, and J. E. Bernard, Phys. Rev. Lett. 65, 353 (1990).

19