- Email: [email protected]
Bond-order bond energy model for alloys
Acta Materialia 179 (2019) 406e413
Contents lists available at ScienceDirect
Acta Materialia journal homepage: www.elsevier.com/locate/actamat
Full length article
Bond-order bond energy model for alloys Christian Oberdorfer, Wolfgang Windl* Department of Materials Science and Engineering, The Ohio State University, Columbus, OH, 43210, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 6 February 2019 Received in revised form 26 July 2019 Accepted 25 August 2019 Available online 28 August 2019
We introduce a novel way to parameterize alloy energies in the form of a bond-order bond energy model. There, a bond order function models the transition between competing phases and switches their respective bond energies on and off. For the case of the Ni-Cr-Mo alloys investigated here, which assume face- or body-centered cubic structures, we propose a sigmoidal switching function fitted to the c/a ratio of “Bain-like” cells. With that, the model does an excellent job in describing the DFT-calculated alloy energies. We also show that the average atomic charge density can vary considerably as a function of composition, which can significantly modify alloy bond energies from simple expectation. The fitted bond energies can among others be used to determine phase diagrams, where we find excellent agreement with previous assessments in the solid range for Ni-Cr and Cr-Mo phase diagrams once the necessary entropy terms are added. They also allow quantitative, composition-dependent calculation of chemical potentials, which we use to determine vacancy formation energies in binary NiCr alloys for configurations that have been energy-minimized with Monte Carlo simulations. We show that the resulting regime of negative formation energies is a sign for thermodynamic instability of the underlying crystal and lies in the two-phase concentration range in the phase diagram, resulting in a holistic picture that unites defect and phase stability through the fully quantitative link between stoichiometry and chemical potentials enabled by the proposed bond energy model. © 2019 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.
Keywords: Alloy energetics and thermodynamics Phase diagrams First-principles calculations
1. Introduction and background Thermodynamic modeling of alloys requires the ability to determine the energy of the different atoms in the alloy, which traditionally has been a drawback of first-principles calculations where the energy is typically only known for the entire system, but not for the different atoms. For example, a fast parameterization of the total energy in terms of the atoms is desirable for Monte Carlo calculations such that many steps can be run quickly, or that chemical potentials (which are determined from the average energy of each species in the alloy) needed for acceptance criteria in e.g. semi-grandcanonical ensembles can be determined on the fly [1]. The chemical potentials are also needed for solidsolution theory to determine phase stability within e.g. the regular solution model [2,3], and calculation of the thermodynamics of point defects [4], non-stoichiometric compounds [5], or interface energetics [6]. One of the oldest approaches to determine such atomic energies
* Corresponding author. E-mail address: [email protected] (W. Windl). https://doi.org/10.1016/j.actamat.2019.08.048 1359-6454/© 2019 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.
is what is in its most general version called quasichemical solution theory (QST), where energies are assigned to the different bond types. As an example, in an A-B alloy, there are A-A, A-B, and B-B bonds with energies εAA, εAB, and εBB which contribute to the average atomic energy in fractions that depend on the overall composition as discussed below. However, QST has at times been labeled as a “naïve theory” [3] and not rigorous, and more complicated approaches with increased numbers of (fitting) parameters have been developed, most notably the cluster expansion [7,8], which has found a large number of very successful applications. While the cluster expansion works without question more successfully for a large number of alloy systems that stay within the same crystal structure and can be easily automated [9,10], it sacrifices the inherent simplicity and especially direct access to the bond energies that the QST offers, which are important input for assessment of surface processes that involve local adhesion, relevant for corrosion [11] and atom-probe tomography [12,13]. An important step to overcome these shortcomings is the energy density method, originally derived in the early 1990's by Martin et al. [14] and recently adopted to modern DFT approaches [15], which is elegant and in principle rather accurate, but as we found is much more involved and often hard to converge in comparison to a
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
simple bond energy fit within QST as we will propose here, and thus is not commonly used in alloy calculations. QST is in principle a pair-potential approach with fixed (average) bond length, whose drawbacks for both elemental and alloy systems are well known since many decades. To overcome these drawbacks, corrections to the pair-potential have been developed, historically first and overall most successful the bond-order potentials, which we understand here in the most general sense that the pair energies are corrected depending on the number (and often distance) of the neighbors. This concept has been successfully applied to both metals and non-metals. For metals, if the same number of electrons on an atom forms different numbers of bonds, the average number of electrons that “glue” the atoms together in a metallic bond changes, and with it the bond strength. In a very simplified way, one could say the more glue, the stronger the bond. Potentials that adjust the pair-potential bond strength in this way are for example embedded-atom method (EAM) [16] or glue [17] potentials. For covalently bonded solids, the hybridization can change with the number of bonds, leading to an adjustment of bond energy and also bond length. Examples include the Tersoff family of potentials [18,19] and Pettifor-style bond-order potentials [20]. Thus, while the bond-order idea has been applied for many decades to interatomic potentials, it has to date not been applied to bond energies within the QST. Here, we will show that introducing appropriate bond-order functions to establish a bond-order bond energy (BOBE) model can result in models that provide rather high quality energy fits, are simple, provide bond energies needed for assessment of thermodynamics and bond strengths, and even overcome one of the major drawbacks of cluster-expansion methods, which is the inability to parameterize alloys that span different crystal structures such as face-centered cubic (fcc) and body-centered cubic (bcc).
407
Fig. 1. Model crystal with two atom types A (red) and B (blue) to demonstrate the bond alloy energy model, where red lines are AA, blue lines BB, and purple lines AB bonds. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
2. Alloy energies e theoretical background The most traditional way to parameterize alloy energies, such as described by the early work of Becker [21], was to assign energies to different bonds, determine the energy of each atom by counting its bonds to the different types of neighbors, and add up the respective bond energies (which are equally distributed to the two atoms at the ends of the bond). As an example, the energy in a 2D model crystal as shown in Fig. 1 is
Eð1Þ ¼
2EAA þ 2EAB E þ 3EAB and Eð2Þ ¼ BB 2 2
(1)
for atoms 1 and 2, respectively. While it is easy to count bonds in ordered alloys within one phase, it becomes unclear what to do in a system that, depending on composition, can switch between different phases such as the fcc and bcc structures. Although the fcc structure has 12, and the bcc structure 8 nearest neighbors, both can be represented within the same unit cell, but with different c/a ratios (Fig. 2) as proposed by Bain for the martensitic transformation [22]. To deal with the transition, we propose here a model that determines the nature of the structure from the relaxed c/a ratio of a p 32-atom fcc supercell s, ffiffiffi which changes from sfcc ¼ 1 to sbcc ¼ 1 = 2 when the structure changes from fcc to bcc. In our model, we couple the number of neighbors used for the bond-energy fit to the c/a ratio, which allows to smoothly transition from fcc to bcc structures upon alloying, where the relaxed c/a ratio as a function of composition takes on the role of a bond-order parameter with a sigmoidal dependence on composition. Fig. 3 shows the examples of Ni-Cr and Ni-Mo binary alloys, whose bond order parameters can be well fitted by an error function to
Fig. 2. Bain cell with (a) c/a ¼ 1 (fcc) and (b) c/a ¼ 1/√2 (bcc). Nearest neighbors are indicated by sticks, 8 red ones and 4 blue ones (all equivalent) resulting in 12 for fcc, and 8 red ones for bcc. The thin blue dashed lines in (b) indicate the now broken former nearest neighbor bonds after transition from fcc to bcc. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
results of DFT calculations described below, with a resulting parameterization of
x 0:556 s þ sbcc þ fcc 0:118 2 2 sfcc sbcc x 0:597 sfcc þ sbcc þ erf fMo ðxÞ ¼ 0:133 2 2 fCr ðxÞ ¼
sfcc sbcc
erf
(2)
The center of the error function is at 56% for Cr and 60% for Mo, indicating that Cr is a stronger beta-stabilizer than Mo, which we therefore can quantify through this fit. Using the free-form modeling tool Eureqa [23], which can identify best-fitting
408
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
3. Results for bond energies Overall, 74 SQS structures were included, whose compositions are shown in Fig. 5. The calculated bonding or cohesive energies for the different structures, which were determined from subtracting the atomic energies in vacuum from the calculated total energies, were then fitted by a bond-energy expression analogous to what was discussed for the model structure in Fig. 1, where the 12 bonds shown in Fig. 2(a) were considered for fcc, and the 8 shown in Fig. 2(b) for bcc structures, with the ternary switching function f(xA,xB,xC) from Eq. (3) being used to activate either fcc or bcc energies, depending on the relaxed c/a ratio. The fitting equation used is therefore Fig. 3. c/a ratio for Ni-Cr and Ni-Mo alloys with 32-atom SQS cells after relaxation with DFT-PBE. The solid lines are fits with error functions, where solid lines are the singlesolute fits from Eq. (2) and dashed lines the double-solute fit from Eq. (3).
analytical equations and their parameters for multi-variable datasets, we found that both sigmoids can be described with nearly the same accuracy by a single sigmoid,
fCr;Mo ðxCr ;xMo Þ
sfcc sbcc 2
erf
x s þs x Cr bcc þ Mo 4:424 þ fcc 0:126 0:135 2 (3)
as illustrated by the comparison shown in Fig. 3. This allows extending our approach to the ternary NiCrMo system. For the calculation of the alloy energies, we have used DFT within VASP with PAW-PBE pseudo-potentials, relaxing both atomic positions and lattice vectors with a constant-pressure algorithm at zero pressure. We employed 6 6 6 k-points for the 32-atom supercells and a cutoff energy of 330 eV. Our structures were highly optimized special quasirandom structures (SQS) using the ATAT tool mcsqs [24]. This is demonstrated by the bond fractions x in an A1-xBx alloy displayed in Fig. 4 as determined from the various SQS cells, which follow very closely the ideal random quadratic distribution functions derived from statistical mechanics [21,25],
xAA ¼
nAA n n ¼ ð1 xÞ2 ; xAB ¼ AB ¼ 2xð1 xÞ; xBB ¼ BB ¼ x2 : ntot ntot ntot
(4)
Fig. 4. Bond fractions of the SQS cells used in the present work, where red squares, blue diamonds, and green triangles show AB, AA, and BB bond fractions in an A1-xBx alloy as a function of x. The solid lines represent the ideal fully random bond fractions shown in Eq. (4). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
f ðx ;x ;x Þ s A B C fcc bcc fcc fcc fcc fcc EðxA ;xB ;xC Þ ¼ nfcc AA EAA þ nAB E AB þ nBB EBB þ::: sfcc sbcc s f ðx ;x ;x Þ A B C bcc fcc bcc bcc bcc bcc þ nbcc AA EAA þ nAB EAB þ nBB E BB þ ::: sfcc sbcc (5) fcc fcc where nfcc AA ; nAB ; nBB are the number of the different bond types in bcc bcc the supercell in fcc, and nbcc AA ; nAB ; nBB those in bcc. The E's are the corresponding bond energies. The resulting fit is shown in Fig. 6(a), and the resulting bond energies are in Table 1. The fit can be done with very high fidelity, as is shown by the fact that all points in Fig. 6(a) are very closely following the x ¼ y relationship (the fitted equation with error bars determined from the 95% confidence intervals is Efit ¼ (0.00 ± 0.02) þ (0.999 ± 0.003)$ EDFT), and that in Fig. 6(b), the enthalpy of mixing can be reproduced within a few hundredths of an eV. Also, the use of the bond-order function obviously allows reproducing a concentration dependence of the enthalpy of mixing way outside of the parabolic dependence of the regular solution model. When looking at the bond energies in Table 1, the probably biggest surprise is that the Ni-Ni bond energy is found to be lower in bcc than in its native fcc structure, as indicated by the fact that 2/3 of the bcc bond energy (which should be indicative of the energy of the atom being distributed over 12 instead of 8 atoms) is 0.05 eV lower than the fcc bond energy of 0.80 eV. In order to understand that, we perform a Bader analysis and calculate the average number of valence electrons on Ni and Cr atoms as well as average Ni-Ni and Cr-Cr bond lengths in the Ni-Cr system, shown in Fig. 7 in
Fig. 5. Compositions for which SQS cells within the Cr-Mo-Ni ternary system were relaxed.
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
409
Fig. 6. (a) Energies from fitting Eq. (6) (y-axis) vs. the DFT-calculated energies used for the fit (x-axis). Perfect agreement means that a point is exactly on the y ¼ x line. (b) Enthalpy of mixing for Ni-Cr, Ni-Mo, and Cr-Mo binary alloys as a function of composition; dots are directly calculated values; lines are calculated with the fitted bond energy model. Dotted lines are negative values of entropy of mixing energy contributions, eTSmix, for 500, 1000, 1500, and 2000 K.
Table 1 Bond energies in eV fitted with Eq. (5) to DFT cohesive energies for 74 32-atom SQS cells and 4-atom fcc and L12 intermetallics (IM).
Ni-Ni Ni-Cr Ni-Mo Cr-Cr Cr-Mo Mo-Mo
S E SQ fcc
S E SQ bcc
S 2 =3,E SQ bcc
EIM fcc
E IM bcc
0.80 0.73 0.94 0.63 0.80 0.99
1.27 1.03 1.32 1.01 1.26 1.60
0.85 0.69 0.88 0.68 0.84 1.07
0.80 0.73 0.93 0.62 0.79 1.00
1.18 1.04 1.33 1.01 1.29 1.60
comparison to the corresponding values for elemental fcc and bcc Ni and Cr with relaxed lattice constants. For the elemental systems, we see a strong dependence of bond length on valence electron number, which makes sense since the Coulomb repulsion gives electrons a large “volume” in solids which changes the volume of the solid accordingly [26]. However, in the alloy, we find that Cr is shrinking to a similar degree as Ni is growing, and so the average homogeneous bond lengths around Ni and Cr atoms stay nearly constant at average values of 2.49 Å and 2.46 Å, respectively (Fig. 7). Thus, while the number of electrons increases on Ni and shrinks on Cr, their sizes stay nearly unchanged, and the alloy-bcc-Ni atom
ends up being different form the bulk-bcc-Ni atom. The effect of the changing valence electron count with constant volume on the energy difference between fcc and bcc is shown in Fig. 8. While for relaxed structures, fcc stays always lower in energy in the case of Ni, this situation is reversed for constant-volume calculations in fcc structures with the calculated average bond lengths of 2.49 Å for Ni-Ni and 2.46 Å for Cr-Cr. For Ni, once there are more than 10.145 electrons per atom, which happens at a Cr-concentration of ~25%, bcc starts to have a lower energy than fcc (Fig. 8(a)); thus Ni is never in a bcc configuration within the fcc stability range, and the lower bcc energy is a result of the “compressed state” where the Ni bond lengths are too small for the number of electrons on the Ni atom. Within the charge range of the alloy system (down to 5.55 electrons per atom), bcc stays always lower in energy than fcc for the case of Cr (Fig. 8(b)). What influence does the random alloy in total have on the bond energies, or in other words, how local are they? In order to answer this question, we have studied if the bond energies can be derived from simple intermetallic structures without going through larger SQS supercells. For fcc, we have chosen four-atom conventional fcc cells, which can represent both the fcc (A1) structure when all four positions are occupied by the same elemental species, or the L12 AB3 intermetallic structure where the corner position is occupied
Fig. 7. Average (a) Ni-Ni and (b) Cr-Cr bond lengths vs. respectively Ni and Cr valence electron count in Ni-Cr alloys over the entire composition range, in comparison to bond length vs. valence electron count in the same range for relaxed lattice constants in fcc and bcc elemental Ni (a) and Cr (b). The valence electron count is linear with composition, which is given on the top axis. The inset in (b) shows the average calculated Bader charges in e/atom on Ni and Cr relative to the neutral valence electron numbers of 10 and 6, respectively, as a function of Cr concentration.
410
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
Fig. 8. Energy difference between fcc and bcc for (a) Ni and (b) Cr over the calculated valence electron range in the alloy system. Blue: relaxed lattice vectors; red: fixed lattice constant. Dashed lines are at zero energy. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
by a different element. For bcc, we use the 2-atom bcc cell (A2 structure), which forms B2 when occupied by two different elements. As can be seen from the results in Table 1, for fcc structures, the intermetallic results are very close to the SQS results, with differences of no more than 0.01 eV. For bcc structures most of the bonds are rather close as well (last column Table 1), while the “compressed” state that Ni, the only native fcc-element in the mix, feels in random alloys as discussed above makes its alloy energy different from the bcc energy by 0.09 eV. Thus, while a quick exploration may be sensible from the simple intermetallics and can be performed in a few minutes even on a desktop system, for good accuracy, the alloy effects need to be included into the calculations, and larger supercells than simple intermetallics are needed. The 32-atom cells we have used here seem to be a good compromise between computational economy and accuracy. 4. Application of bond energies and discussion The fitted bond energies are useful by themselves, for example to analyze alloy effects in field ion emission and atom probe tomography or assess dissolution in corrosion processes. For the former, bond energies are used to assess surface adhesion energies, which are needed to calculate evaporation fields [12], since we have shown in the past that the order of the magnitude of surface adhesion energies can be predicted from bulk bond energies [13]. For the latter, a similar correlation is used to predict the surface adhesion of metals in alloys and predict how easily atoms dissolve from a surface in corrosive environment [11,27]. Beyond that, we demonstrate here the usefulness of the BOBE model for general thermodynamic questions with the help of two examples. We use it first as a tool to calculate alloy phase diagrams, and second for the calculation of the chemical potentials of the atoms in the alloy. For well-behaved phase diagrams, such as simple eutectics or systems with a miscibility gap, phase boundaries can be determined from the intersections of the enthalpy of mixing with the negative of the entropy-related energy. On the most accessible level, the entropy can just contain the mixing contribution, calculated from the Stirling formula, Emix ðx;TÞ ¼ kB Tðx ln x þ ð1 xÞlnð1 xÞÞ. Fig. 6(b)
mA ðxA ; xB ; xC Þ ¼
shows as example the entropy of mixing energy contributions for different temperatures in comparison to the enthalpy of mixing curves for Ni-Cr, Cr-Mo, and Ni-Mo binaries from the BOBE model. From the intersections of the curves, one can immediately tell that there is considerable solubility on the Ni side for both Ni-Cr and NiMo alloys, which is not the case on the bcc-metal ends. The corresponding phase boundaries in the solid range are shown for the examples of Ni-Cr and Cr-Mo in Fig. 9 in comparison to wellestablished phase diagrams from Refs. [28,29]. While the overall topology of the phase regions is already visible with just the mixing entropy added, the excellent energetic quality of the BOBE becomes clear once the vibrational contribution to the entropy is also added. We calculate the latter from phonon G-point calculations for our 32atom SQS cells as described in Ref. [30]. With that, the phase boundaries in both diagrams can be captured in rather close agreement, demonstrating the unique predictive power of the BOBE model, which is easily transferable to multicomponent alloys due to the natural pairwise separability of bond energies. For the second example, we determine quantitative chemical potentials for our model system. In the explicit definition, the chemical potential of an atom is given by the energy change of the system when the atom is added or removed. Since we have chosen free atoms as reference for our bond energies, adding or removing an atom from the system changes its energy exactly by the energy of the bonds that the atom forms when entering the ensemble. In thermodynamic average in a random alloy, the chemical potential can thus be simply determined from the perfect bond fractions in Eq. (4) and the corresponding bond energies in Table 1. This is very useful, since except for the tedious energy-density approach [14], determining quantitative chemical potentials is usually not possible in standard DFT approaches, where only the total energy for all atoms in the system is determined. With the chemical potentials, phase stability can be discussed and important quantities can be determined, such as defect formation or surface energies. Mathematically, the chemical potential of one atom of type A in the alloy can be determined from all bonds that connect to at least one A atom. Considering that all homo but only half of the hetero bonds belong to A,
NN fcc fcc fcc fcc fcc fcc f ðxA ; xB ; xC Þ sbcc 2xAA EAA þ xAB Efcc AB þ xAB EAB 4xA sfcc sbcc
NN bcc bcc bcc bcc bcc sfcc f ðxA ; xB ; xC Þ bcc 2xAA EAA þ xbcc þ AB EAB þ xAB EAB 4xA sfcc sbcc
(6)
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
411
Fig. 10. Chemical potentials in Cr-Ni and Ni-Mo binary alloys, calculated from Eq. (6) using the bond energies in Table 1 and the bond order parameter functions from Eq. (3). Blue, purple, and red lines are chemical potentials of Cr, Mo, and Ni as a function of composition, respectively. The green dots are the average energy of all atoms in the system. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
Fig. 9. (a) Ni-Cr phase diagram. Black lines are CALPHAD results, redrawn after Ref. [28]. (b) Cr-Mo phase diagram. Black lines are redrawn after Ref. [29]. In both panels, red lines are phase boundaries determined from sign change in free-energy of mixing, calculated from the sum of the enthalpy of mixing from the bond-energy model results (Fig. 6(b)) and entropy of mixing approximated by the Stirling formula. For the circles, the full vibrational entropy has been added, calculated from harmonic phonons as described in Ref. [30]. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
The resulting chemical potentials for Cr-Ni and Ni-Mo are shown in Fig. 10. The switching functions (Fig. 3) already contain most of the phase stability information, such as that fcc and bcc alone rule in either of the flat regions of the sigmoid. The non-flat regions of the sigmoids indicate the transition from one to the other, resulting in the undulation of the otherwise linear chemical potentials for close-to-equiatomic concentrations, and it is a priori not clear what meaning the transition region has thermodynamically. As suggested in Kirchheim's work on the extended Gibbs adsorption isotherm [31,32], where the necessary lowering of defect formation energies by solute segregation in alloys is shown and the potential meaning of negative formation energies is discussed, a useful way to probe the possible existence of the transition-region “martensite” consists in using the derived chemical potentials and atomic energies to calculate vacancy formation energies, which we do for the example of CrNi alloys. We show that our approach also offers a solution to one of the most confusing problems for multicomponent alloys, which is to determine the appropriate chemical potential needed to calculate the formation energy of a defect. The formation energy of a defect in a multicomponent (here, as an example, a vacancy in a binary) alloy is in general given by
Ef ðV; An Bm Þ ¼ EðV; An Bm Þ nmA nmB ;
(7)
where E(V, AnBm) is the energy of a cell with one vacancy, n A, and m
B atoms. In other words, the formation energy is the difference between the defective cell and a perfect cell with the same stoichiometry, i.e., the same number of n and m atoms, which of course is usually not known, since n and m, once a vacancy is introduced, do not add up to the stoichiometry of a perfect periodic cell, but are one atom short of what can be calculated. Therefore, usually the energy of either a perfect Anþ1Bm or AnBmþ1 cell is taken as the basis, and corrected for by the chemical potential of the extra atom, which is usually determined from the assumption that it is equal to the energy of an atom in elemental A or B material as needed, which then is labeled as A or B-rich case, with frequently no clear discussion which choice would be most physical, or what A or B rich even means for an alloy [33]. In our case, the procedure is rather clear, since from Eq. (6), the chemical potentials can simply be calculated for the exact stoichiometry needed in Eq. (7) with no ambiguity. We note that interstitials, or defects with any other type of stoichiometry (including multidimensional defects such as dislocations or grain boundaries) can be calculated in an analogous way. Also, in order to minimize the overall error, E(Anþ1Bm) or E(AnBmþ1) can and probably should be used in Eq. (7) and only the stoichiometry difference expressed by the chemical potential. The next question is, what configuration one should choose for the defective structure, since from an Anþ1Bm or AnBmþ1 SQS supercell, any of the A or B atoms respectively can be removed resulting in the same stoichiometry. We decided here to pick the lowest-energy configuration which we determined by running a Monte-Carlo energy minimization [34,35] with respect to the distribution of A and B atoms and the vacancy, as will be described in more detail elsewhere. This assumes that in the random alloy, the atomic mobility around a vacancy (which mediates transport in these metal alloys) is high enough that the local configuration reaches the thermodynamically most favorable state, while the long-range mobility is slow enough that the rest of the crystal remains a random alloy. This is a sensible, but not the only possible choice, and depends on the situation of interest. The resulting energies as a function of concentration are shown in Fig. 11, in comparison to the previously determined Ni-Cr switching function from Eq. (2). It can be seen that the formation energies from both ends decrease towards the center, as suggested by Kirchheim's work [31,32], and assume negative values in the transition region of the switching function, which has its inflexion point at 55.6% Cr. Since, once chemical potentials are linked to stoichiometry as is the case here, negative formation energies indicate an unstable bulk crystal for dilute vacancy concentrations [31e33], the vacancy formation energy within this approach can be
412
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
Fig. 11. (a) Formation energies for vacancies in Ni-Cr alloys, calculated from MCminimized alloy-vacancy configurations and ideal chemical potentials as defined in Eq. (6). (b) Ni-Cr bond switching function from Eq. (3).
used to estimate phase stability limits within alloy systems. Thus, the negative formation energies suggest that there should be a concentration gap between fcc and bcc phases (at a minimum at low temperatures), which the eutectic Ni-Cr phase diagram indeed shows in the form of the central two-phase region, whose width of course is determined by the full free energies rather than the zerotemperature formation energies discussed here. This then gives an interpretation to the transition region of the switching function, which thus would not indicate some “mixed intermediate” martensite-like structure between bcc and fcc, but instead a concentration regime of unstable single-phase bulk phases and thus a multiphase region of the phase diagram. This example thus illustrates that the chemical potentials from the bond-alloy approach enable exciting new angles for a more in-depth examination of the interdependencies between formation energies and phase stability. Since our approach is easily extendible to more-than-two component systems (and in fact here already includes the ternary Ni-CrMo system), it can be extended to multicomponent alloys in a straightforward way where point defects in many-component systems are an active field of research [36,37]. 5. Conclusions In summary, we have introduced a novel way to parameterize alloy energies in the form of a bond-order bond energy model. There, we have introduced a bond order function that models the switching between competing phases. The bond order function switches the respective bond energies for the different phases on and off, which allows to have a continuous model through the phase instability composition. For the case of the Ni-Cr-Mo alloys investigated here, the c/a ratios pffiffiffi calculated from 32-atom cubic SQS cells vary from 1 for fcc to 1 = 2 for bcc, and their function dependence suggests using a sigmoidal switching function, which we chose to describe by an error function. We have applied this model to the Ni-Cr-Mo alloy as a test system, since its stable crystal structures fluctuate between fcc and bcc. We find that our bond-order bond energy model describes the DFT energetics in this system surprisingly well and allows extracting well validated bond energies for both fcc and bcc phases that without concentration-dependent adjustment result in an excellent fit throughout the entire composition range.
We have also shown that for some elements in alloys, in our case specifically Ni, their average atomic charge density can vary considerably as a function of composition. This “charge compression” explains why Ni has lower bcc than fcc bond energies in the examined alloys. It is important to note that the higher charge density that lowers the bcc bond energy is an alloying effect, since the Ni atoms attract additional electrons from the other alloying elements while the alloy structure restricts their resulting volume expansion. In applying the fitted BOBE model, we have shown that, besides the obvious usefulness of bond energies for questions that relate bond energies to processes such as those involved in field ion emission or corrosion, bond energies can be used to determine phase diagrams which, once vibrational entropy is added, show surprisingly good agreement with experimentally derived phase diagrams, as well as quantitative, composition-dependent chemical potentials. From vacancy calculations for binary NiCr alloys for configurations that have been configurationally energy-minimized with Monte Carlo simulations, we could show that the regime of negative formation energies, a sign for thermodynamic instability of the underlying crystal, suggests the existence of a two-phase concentration range in the phase diagram. This suggests strongly that these concentration-dependent chemical potentials can be used to study thermodynamics, defects, and phase stability within a holistic approach since there is a fully quantitative link between stoichiometry and chemical potentials, which opens exciting possibilities for future work. Acknowledgements This work was funded by the Center for Performance and Design of Nuclear Waste Forms and Containers, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award # DE-SC0016584. C.O. also gratefully acknowledges support from the Alexander von Humboldt Foundation through the Feodor-Lynen research fellowship. Our work was supported by allocation of computing time from the Ohio Supercomputer Center under grant number PAS0072. References [1] D.A. Kofke, E.D. Glandt, Monte Carlo simulation of multicomponent equilibria in a semigrand canonical ensemble, Mol. Phys. 64 (6) (1988) 1105e1131. [2] D.R. Gaskell, Introduction to the Thermodynamics of Materials, Taylor and Francis, New York, 2008. [3] R. DeHoff, Thermodynamics in Materials Science, Chapman and Hall/CRC, 2006. [4] A. Agrawal, R. Mishra, L. Ward, K. Flores, M.,W. Windl, An embedded atom method potential of beryllium, Model. Simul. Mater. Sci. Eng. 21 (8) (2013), 085001. [5] R. Mishra, O.D. Restrepo, P.M. Woodward, W. Windl, First-principles study of defective and nonstoichiometric Sr2FeMoO6, Chem. Mater. 22 (22) (2010) 6092e6102. [6] D. Sen, W. Windl, Ab-initio study of energetics of the Si(001)-LaAlO3 interface, J. Comput. Theor. Nanosci. 4 (1) (2007) 57e64. [7] A. Zunger, First-principles statisitical mechanics of semiconductor alloys and intermetallic compounds, in: P.E.A. Turchi, A. Gonis (Eds.), Statistics and Dynamics of Alloy Phase Transformations, Springer, New York, 1994, pp. 361e420. [8] D. De Fontaine, Cluster Approach to order-disorder transformations in alloys, Solid State Phys. 47 (1994) 33e176. [9] A. van de Walle, M. Asta, G. Ceder, The alloy theoretic automated toolkit: a user guide, Calphad 26 (4) (2002) 539e553. [10] A. van de Walle, Multicomponent multisublattice alloys, nonconfigurational entropy and other additions to the Alloy Theoretic Automated Toolkit, Calphad 33 (2) (2009) 266e278. [11] P. Marcus, On some fundamental factors in the effect of alloying elements on passivation of alloys, Corros. Sci. 36 (12) (1994) 2155e2158. [12] L. Yao, T. Withrow, O.D. Restrepo, W. Windl, E.A. Marquis, Effects of the local structure dependence of evaporation fields on field evaporation behavior, Appl. Phys. Lett. 107 (24) (2015) 241602. [13] C. Oberdorfer, T. Withrow, L.J. Yu, K. Fisher, E.A. Marquis, W. Windl, Influence
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
[14] [15]
[16]
[17] [18] [19] [20] [21] [22] [23] [24]
[25] [26]
of surface relaxation on solute atoms positioning within atom probe tomography reconstructions, Mater. Char. 146 (2018) 324e335. N. Chetty, R.M. Martin, First-principles energy density and its applications to selected polar surfaces, Phys. Rev. B 45 (11) (1992) 6074e6088. B. Lee, D.R. Trinkle, Energetics of rutile TiO2 vicinal surfaces with 〈001〉 steps from the energy density method, J. Phys. Chem. C 119 (32) (2015) 18203e18209. M.S. Daw, M.I. Baskes, Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals, Phys. Rev. B 29 (12) (1984) 6443e6453. F. Ercolessi, M. Parrinello, E. Tosatti, Simulation of gold in the glue model, Philos. Mag. A 58 (1) (1988) 213e226. J. Tersoff, Empirical interatomic potential for carbon, with application to amorphous carbon, Phys. Rev. Lett. 61 (25) (1988) 2879e2882. D.W. Brenner, The art and science of an analytic potential, Phys. Status Solidi b 217 (1) (2000) 23e40. D.G. Pettifor, I.I. Oleinik, Analytic bond-order potentials beyond TersoffBrenner. I. Theory, Phys. Rev. B 59 (13) (1999) 8487e8499. R. Becker, Die Keimbildung bei der Ausscheidung in metallischen Mischkristallen, Ann. Phys. 424 (1e2) (1938) 128e140. E.C. Bain, The nature of martensite, Trans. AIME 70 (1924) 25e35. M. Schmidt, H. Lipson, Eureqa (software). www.nutonian.com, 2014. (Accessed 5 May 2018). A. van de Walle, P. Tiwary, M. de Jong, D.L. Olmsted, M. Asta, A. Dick, D. Shin, Y. Wang, L.Q. Chen, Z.K. Liu, Efficient stochastic generation of special quasirandom structures, Calphad 42 (2013) 13e18. T.L. Hill, Chapter 6. Distribution Functions and the Theory of the Liquid State, Statistical Mechanics, McGraw-Hill, New York, 1956. W. Windl, M.S. Daw, Predictive process simulation and Ab-Initio calculation of
[27]
[28] [29] [30] [31]
[32]
[33] [34] [35] [36]
[37]
413
the physical volume of electrons in silicon, in: Technical Proceedings of the 2002 International Conference on Modeling and Simulation of Microsystems: Quantum Effects, Quantum Devices and Spintronics, 2002, pp. 494e497. G.S. Frankel, J.D. Vienna, J. Lian, J.R. Scully, S. Gin, J.V. Ryan, J. Wang, S.H. Kim, W. Windl, J. Du, A comparative review of the aqueous corrosion of glasses, crystalline ceramics, and metals, npj Mater. Degrad. 2 (1) (2018) 15. P. Nash, The Cr-Ni (Chromium-Nickel) system, Bull. Alloy. Phase Diagrams 7 (5) (1986) 466. B. Predel, Cr-Mo (Chromium-Molybdenum), in: O. Madelung (Ed.), Cr-Cs e Cu-Zr, Springer Berlin Heidelberg, Berlin, Heidelberg, 1994, pp. 1e5. W. Luo, W. Windl, First principles study of the structure and stability of carbynes, Carbon 47 (2) (2009) 367e383. R. Kirchheim, Reducing grain boundary, dislocation line and vacancy formation energies by solute segregation. I. Theoretical background, Acta Mater. 55 (15) (2007) 5129e5138. R. Kirchheim, Reducing grain boundary, dislocation line and vacancy formation energies by solute segregation: II. Experimental evidence and consequences, Acta Mater. 55 (15) (2007) 5139e5148. R. Mishra, O. Restrepo, D.,A. Kumar, W. Windl, Native point defects in binary InP semiconductors, J. Mater. Sci. 47 (21) (2012) 7482e7497. C. Niu, W. Windl, M. Ghazisaeidi, Multi-Cell Monte Carlo Relaxation method for predicting phase stability of alloys, Scr. Mater. 132 (2017) 9e12. C. Niu, Y. Rao, W. Windl, M. Ghazisaeidi, Multi-Cell Monte Carlo Method for Phase Prediction, 2018, 04092 eprint arXiv:1811, arXiv:1811.04092. S.C. Middleburgh, D.M. King, G.R. Lumpkin, M. Cortie, L. Edwards, Segregation and migration of species in the CrCoFeNi high entropy alloy, J. Alloy. Comp. 599 (2014) 179e182. Z. Wang, C.T. Liu, P. Dou, Thermodynamics of vacancies and clusters in highentropy alloys, Phys. Rev. Mater. 1 (4) (2017), 043601.
Contents lists available at ScienceDirect
Acta Materialia journal homepage: www.elsevier.com/locate/actamat
Full length article
Bond-order bond energy model for alloys Christian Oberdorfer, Wolfgang Windl* Department of Materials Science and Engineering, The Ohio State University, Columbus, OH, 43210, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 6 February 2019 Received in revised form 26 July 2019 Accepted 25 August 2019 Available online 28 August 2019
We introduce a novel way to parameterize alloy energies in the form of a bond-order bond energy model. There, a bond order function models the transition between competing phases and switches their respective bond energies on and off. For the case of the Ni-Cr-Mo alloys investigated here, which assume face- or body-centered cubic structures, we propose a sigmoidal switching function fitted to the c/a ratio of “Bain-like” cells. With that, the model does an excellent job in describing the DFT-calculated alloy energies. We also show that the average atomic charge density can vary considerably as a function of composition, which can significantly modify alloy bond energies from simple expectation. The fitted bond energies can among others be used to determine phase diagrams, where we find excellent agreement with previous assessments in the solid range for Ni-Cr and Cr-Mo phase diagrams once the necessary entropy terms are added. They also allow quantitative, composition-dependent calculation of chemical potentials, which we use to determine vacancy formation energies in binary NiCr alloys for configurations that have been energy-minimized with Monte Carlo simulations. We show that the resulting regime of negative formation energies is a sign for thermodynamic instability of the underlying crystal and lies in the two-phase concentration range in the phase diagram, resulting in a holistic picture that unites defect and phase stability through the fully quantitative link between stoichiometry and chemical potentials enabled by the proposed bond energy model. © 2019 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.
Keywords: Alloy energetics and thermodynamics Phase diagrams First-principles calculations
1. Introduction and background Thermodynamic modeling of alloys requires the ability to determine the energy of the different atoms in the alloy, which traditionally has been a drawback of first-principles calculations where the energy is typically only known for the entire system, but not for the different atoms. For example, a fast parameterization of the total energy in terms of the atoms is desirable for Monte Carlo calculations such that many steps can be run quickly, or that chemical potentials (which are determined from the average energy of each species in the alloy) needed for acceptance criteria in e.g. semi-grandcanonical ensembles can be determined on the fly [1]. The chemical potentials are also needed for solidsolution theory to determine phase stability within e.g. the regular solution model [2,3], and calculation of the thermodynamics of point defects [4], non-stoichiometric compounds [5], or interface energetics [6]. One of the oldest approaches to determine such atomic energies
* Corresponding author. E-mail address: [email protected] (W. Windl). https://doi.org/10.1016/j.actamat.2019.08.048 1359-6454/© 2019 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.
is what is in its most general version called quasichemical solution theory (QST), where energies are assigned to the different bond types. As an example, in an A-B alloy, there are A-A, A-B, and B-B bonds with energies εAA, εAB, and εBB which contribute to the average atomic energy in fractions that depend on the overall composition as discussed below. However, QST has at times been labeled as a “naïve theory” [3] and not rigorous, and more complicated approaches with increased numbers of (fitting) parameters have been developed, most notably the cluster expansion [7,8], which has found a large number of very successful applications. While the cluster expansion works without question more successfully for a large number of alloy systems that stay within the same crystal structure and can be easily automated [9,10], it sacrifices the inherent simplicity and especially direct access to the bond energies that the QST offers, which are important input for assessment of surface processes that involve local adhesion, relevant for corrosion [11] and atom-probe tomography [12,13]. An important step to overcome these shortcomings is the energy density method, originally derived in the early 1990's by Martin et al. [14] and recently adopted to modern DFT approaches [15], which is elegant and in principle rather accurate, but as we found is much more involved and often hard to converge in comparison to a
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
simple bond energy fit within QST as we will propose here, and thus is not commonly used in alloy calculations. QST is in principle a pair-potential approach with fixed (average) bond length, whose drawbacks for both elemental and alloy systems are well known since many decades. To overcome these drawbacks, corrections to the pair-potential have been developed, historically first and overall most successful the bond-order potentials, which we understand here in the most general sense that the pair energies are corrected depending on the number (and often distance) of the neighbors. This concept has been successfully applied to both metals and non-metals. For metals, if the same number of electrons on an atom forms different numbers of bonds, the average number of electrons that “glue” the atoms together in a metallic bond changes, and with it the bond strength. In a very simplified way, one could say the more glue, the stronger the bond. Potentials that adjust the pair-potential bond strength in this way are for example embedded-atom method (EAM) [16] or glue [17] potentials. For covalently bonded solids, the hybridization can change with the number of bonds, leading to an adjustment of bond energy and also bond length. Examples include the Tersoff family of potentials [18,19] and Pettifor-style bond-order potentials [20]. Thus, while the bond-order idea has been applied for many decades to interatomic potentials, it has to date not been applied to bond energies within the QST. Here, we will show that introducing appropriate bond-order functions to establish a bond-order bond energy (BOBE) model can result in models that provide rather high quality energy fits, are simple, provide bond energies needed for assessment of thermodynamics and bond strengths, and even overcome one of the major drawbacks of cluster-expansion methods, which is the inability to parameterize alloys that span different crystal structures such as face-centered cubic (fcc) and body-centered cubic (bcc).
407
Fig. 1. Model crystal with two atom types A (red) and B (blue) to demonstrate the bond alloy energy model, where red lines are AA, blue lines BB, and purple lines AB bonds. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
2. Alloy energies e theoretical background The most traditional way to parameterize alloy energies, such as described by the early work of Becker [21], was to assign energies to different bonds, determine the energy of each atom by counting its bonds to the different types of neighbors, and add up the respective bond energies (which are equally distributed to the two atoms at the ends of the bond). As an example, the energy in a 2D model crystal as shown in Fig. 1 is
Eð1Þ ¼
2EAA þ 2EAB E þ 3EAB and Eð2Þ ¼ BB 2 2
(1)
for atoms 1 and 2, respectively. While it is easy to count bonds in ordered alloys within one phase, it becomes unclear what to do in a system that, depending on composition, can switch between different phases such as the fcc and bcc structures. Although the fcc structure has 12, and the bcc structure 8 nearest neighbors, both can be represented within the same unit cell, but with different c/a ratios (Fig. 2) as proposed by Bain for the martensitic transformation [22]. To deal with the transition, we propose here a model that determines the nature of the structure from the relaxed c/a ratio of a p 32-atom fcc supercell s, ffiffiffi which changes from sfcc ¼ 1 to sbcc ¼ 1 = 2 when the structure changes from fcc to bcc. In our model, we couple the number of neighbors used for the bond-energy fit to the c/a ratio, which allows to smoothly transition from fcc to bcc structures upon alloying, where the relaxed c/a ratio as a function of composition takes on the role of a bond-order parameter with a sigmoidal dependence on composition. Fig. 3 shows the examples of Ni-Cr and Ni-Mo binary alloys, whose bond order parameters can be well fitted by an error function to
Fig. 2. Bain cell with (a) c/a ¼ 1 (fcc) and (b) c/a ¼ 1/√2 (bcc). Nearest neighbors are indicated by sticks, 8 red ones and 4 blue ones (all equivalent) resulting in 12 for fcc, and 8 red ones for bcc. The thin blue dashed lines in (b) indicate the now broken former nearest neighbor bonds after transition from fcc to bcc. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
results of DFT calculations described below, with a resulting parameterization of
x 0:556 s þ sbcc þ fcc 0:118 2 2 sfcc sbcc x 0:597 sfcc þ sbcc þ erf fMo ðxÞ ¼ 0:133 2 2 fCr ðxÞ ¼
sfcc sbcc
erf
(2)
The center of the error function is at 56% for Cr and 60% for Mo, indicating that Cr is a stronger beta-stabilizer than Mo, which we therefore can quantify through this fit. Using the free-form modeling tool Eureqa [23], which can identify best-fitting
408
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
3. Results for bond energies Overall, 74 SQS structures were included, whose compositions are shown in Fig. 5. The calculated bonding or cohesive energies for the different structures, which were determined from subtracting the atomic energies in vacuum from the calculated total energies, were then fitted by a bond-energy expression analogous to what was discussed for the model structure in Fig. 1, where the 12 bonds shown in Fig. 2(a) were considered for fcc, and the 8 shown in Fig. 2(b) for bcc structures, with the ternary switching function f(xA,xB,xC) from Eq. (3) being used to activate either fcc or bcc energies, depending on the relaxed c/a ratio. The fitting equation used is therefore Fig. 3. c/a ratio for Ni-Cr and Ni-Mo alloys with 32-atom SQS cells after relaxation with DFT-PBE. The solid lines are fits with error functions, where solid lines are the singlesolute fits from Eq. (2) and dashed lines the double-solute fit from Eq. (3).
analytical equations and their parameters for multi-variable datasets, we found that both sigmoids can be described with nearly the same accuracy by a single sigmoid,
fCr;Mo ðxCr ;xMo Þ
sfcc sbcc 2
erf
x s þs x Cr bcc þ Mo 4:424 þ fcc 0:126 0:135 2 (3)
as illustrated by the comparison shown in Fig. 3. This allows extending our approach to the ternary NiCrMo system. For the calculation of the alloy energies, we have used DFT within VASP with PAW-PBE pseudo-potentials, relaxing both atomic positions and lattice vectors with a constant-pressure algorithm at zero pressure. We employed 6 6 6 k-points for the 32-atom supercells and a cutoff energy of 330 eV. Our structures were highly optimized special quasirandom structures (SQS) using the ATAT tool mcsqs [24]. This is demonstrated by the bond fractions x in an A1-xBx alloy displayed in Fig. 4 as determined from the various SQS cells, which follow very closely the ideal random quadratic distribution functions derived from statistical mechanics [21,25],
xAA ¼
nAA n n ¼ ð1 xÞ2 ; xAB ¼ AB ¼ 2xð1 xÞ; xBB ¼ BB ¼ x2 : ntot ntot ntot
(4)
Fig. 4. Bond fractions of the SQS cells used in the present work, where red squares, blue diamonds, and green triangles show AB, AA, and BB bond fractions in an A1-xBx alloy as a function of x. The solid lines represent the ideal fully random bond fractions shown in Eq. (4). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
f ðx ;x ;x Þ s A B C fcc bcc fcc fcc fcc fcc EðxA ;xB ;xC Þ ¼ nfcc AA EAA þ nAB E AB þ nBB EBB þ::: sfcc sbcc s f ðx ;x ;x Þ A B C bcc fcc bcc bcc bcc bcc þ nbcc AA EAA þ nAB EAB þ nBB E BB þ ::: sfcc sbcc (5) fcc fcc where nfcc AA ; nAB ; nBB are the number of the different bond types in bcc bcc the supercell in fcc, and nbcc AA ; nAB ; nBB those in bcc. The E's are the corresponding bond energies. The resulting fit is shown in Fig. 6(a), and the resulting bond energies are in Table 1. The fit can be done with very high fidelity, as is shown by the fact that all points in Fig. 6(a) are very closely following the x ¼ y relationship (the fitted equation with error bars determined from the 95% confidence intervals is Efit ¼ (0.00 ± 0.02) þ (0.999 ± 0.003)$ EDFT), and that in Fig. 6(b), the enthalpy of mixing can be reproduced within a few hundredths of an eV. Also, the use of the bond-order function obviously allows reproducing a concentration dependence of the enthalpy of mixing way outside of the parabolic dependence of the regular solution model. When looking at the bond energies in Table 1, the probably biggest surprise is that the Ni-Ni bond energy is found to be lower in bcc than in its native fcc structure, as indicated by the fact that 2/3 of the bcc bond energy (which should be indicative of the energy of the atom being distributed over 12 instead of 8 atoms) is 0.05 eV lower than the fcc bond energy of 0.80 eV. In order to understand that, we perform a Bader analysis and calculate the average number of valence electrons on Ni and Cr atoms as well as average Ni-Ni and Cr-Cr bond lengths in the Ni-Cr system, shown in Fig. 7 in
Fig. 5. Compositions for which SQS cells within the Cr-Mo-Ni ternary system were relaxed.
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
409
Fig. 6. (a) Energies from fitting Eq. (6) (y-axis) vs. the DFT-calculated energies used for the fit (x-axis). Perfect agreement means that a point is exactly on the y ¼ x line. (b) Enthalpy of mixing for Ni-Cr, Ni-Mo, and Cr-Mo binary alloys as a function of composition; dots are directly calculated values; lines are calculated with the fitted bond energy model. Dotted lines are negative values of entropy of mixing energy contributions, eTSmix, for 500, 1000, 1500, and 2000 K.
Table 1 Bond energies in eV fitted with Eq. (5) to DFT cohesive energies for 74 32-atom SQS cells and 4-atom fcc and L12 intermetallics (IM).
Ni-Ni Ni-Cr Ni-Mo Cr-Cr Cr-Mo Mo-Mo
S E SQ fcc
S E SQ bcc
S 2 =3,E SQ bcc
EIM fcc
E IM bcc
0.80 0.73 0.94 0.63 0.80 0.99
1.27 1.03 1.32 1.01 1.26 1.60
0.85 0.69 0.88 0.68 0.84 1.07
0.80 0.73 0.93 0.62 0.79 1.00
1.18 1.04 1.33 1.01 1.29 1.60
comparison to the corresponding values for elemental fcc and bcc Ni and Cr with relaxed lattice constants. For the elemental systems, we see a strong dependence of bond length on valence electron number, which makes sense since the Coulomb repulsion gives electrons a large “volume” in solids which changes the volume of the solid accordingly [26]. However, in the alloy, we find that Cr is shrinking to a similar degree as Ni is growing, and so the average homogeneous bond lengths around Ni and Cr atoms stay nearly constant at average values of 2.49 Å and 2.46 Å, respectively (Fig. 7). Thus, while the number of electrons increases on Ni and shrinks on Cr, their sizes stay nearly unchanged, and the alloy-bcc-Ni atom
ends up being different form the bulk-bcc-Ni atom. The effect of the changing valence electron count with constant volume on the energy difference between fcc and bcc is shown in Fig. 8. While for relaxed structures, fcc stays always lower in energy in the case of Ni, this situation is reversed for constant-volume calculations in fcc structures with the calculated average bond lengths of 2.49 Å for Ni-Ni and 2.46 Å for Cr-Cr. For Ni, once there are more than 10.145 electrons per atom, which happens at a Cr-concentration of ~25%, bcc starts to have a lower energy than fcc (Fig. 8(a)); thus Ni is never in a bcc configuration within the fcc stability range, and the lower bcc energy is a result of the “compressed state” where the Ni bond lengths are too small for the number of electrons on the Ni atom. Within the charge range of the alloy system (down to 5.55 electrons per atom), bcc stays always lower in energy than fcc for the case of Cr (Fig. 8(b)). What influence does the random alloy in total have on the bond energies, or in other words, how local are they? In order to answer this question, we have studied if the bond energies can be derived from simple intermetallic structures without going through larger SQS supercells. For fcc, we have chosen four-atom conventional fcc cells, which can represent both the fcc (A1) structure when all four positions are occupied by the same elemental species, or the L12 AB3 intermetallic structure where the corner position is occupied
Fig. 7. Average (a) Ni-Ni and (b) Cr-Cr bond lengths vs. respectively Ni and Cr valence electron count in Ni-Cr alloys over the entire composition range, in comparison to bond length vs. valence electron count in the same range for relaxed lattice constants in fcc and bcc elemental Ni (a) and Cr (b). The valence electron count is linear with composition, which is given on the top axis. The inset in (b) shows the average calculated Bader charges in e/atom on Ni and Cr relative to the neutral valence electron numbers of 10 and 6, respectively, as a function of Cr concentration.
410
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
Fig. 8. Energy difference between fcc and bcc for (a) Ni and (b) Cr over the calculated valence electron range in the alloy system. Blue: relaxed lattice vectors; red: fixed lattice constant. Dashed lines are at zero energy. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
by a different element. For bcc, we use the 2-atom bcc cell (A2 structure), which forms B2 when occupied by two different elements. As can be seen from the results in Table 1, for fcc structures, the intermetallic results are very close to the SQS results, with differences of no more than 0.01 eV. For bcc structures most of the bonds are rather close as well (last column Table 1), while the “compressed” state that Ni, the only native fcc-element in the mix, feels in random alloys as discussed above makes its alloy energy different from the bcc energy by 0.09 eV. Thus, while a quick exploration may be sensible from the simple intermetallics and can be performed in a few minutes even on a desktop system, for good accuracy, the alloy effects need to be included into the calculations, and larger supercells than simple intermetallics are needed. The 32-atom cells we have used here seem to be a good compromise between computational economy and accuracy. 4. Application of bond energies and discussion The fitted bond energies are useful by themselves, for example to analyze alloy effects in field ion emission and atom probe tomography or assess dissolution in corrosion processes. For the former, bond energies are used to assess surface adhesion energies, which are needed to calculate evaporation fields [12], since we have shown in the past that the order of the magnitude of surface adhesion energies can be predicted from bulk bond energies [13]. For the latter, a similar correlation is used to predict the surface adhesion of metals in alloys and predict how easily atoms dissolve from a surface in corrosive environment [11,27]. Beyond that, we demonstrate here the usefulness of the BOBE model for general thermodynamic questions with the help of two examples. We use it first as a tool to calculate alloy phase diagrams, and second for the calculation of the chemical potentials of the atoms in the alloy. For well-behaved phase diagrams, such as simple eutectics or systems with a miscibility gap, phase boundaries can be determined from the intersections of the enthalpy of mixing with the negative of the entropy-related energy. On the most accessible level, the entropy can just contain the mixing contribution, calculated from the Stirling formula, Emix ðx;TÞ ¼ kB Tðx ln x þ ð1 xÞlnð1 xÞÞ. Fig. 6(b)
mA ðxA ; xB ; xC Þ ¼
shows as example the entropy of mixing energy contributions for different temperatures in comparison to the enthalpy of mixing curves for Ni-Cr, Cr-Mo, and Ni-Mo binaries from the BOBE model. From the intersections of the curves, one can immediately tell that there is considerable solubility on the Ni side for both Ni-Cr and NiMo alloys, which is not the case on the bcc-metal ends. The corresponding phase boundaries in the solid range are shown for the examples of Ni-Cr and Cr-Mo in Fig. 9 in comparison to wellestablished phase diagrams from Refs. [28,29]. While the overall topology of the phase regions is already visible with just the mixing entropy added, the excellent energetic quality of the BOBE becomes clear once the vibrational contribution to the entropy is also added. We calculate the latter from phonon G-point calculations for our 32atom SQS cells as described in Ref. [30]. With that, the phase boundaries in both diagrams can be captured in rather close agreement, demonstrating the unique predictive power of the BOBE model, which is easily transferable to multicomponent alloys due to the natural pairwise separability of bond energies. For the second example, we determine quantitative chemical potentials for our model system. In the explicit definition, the chemical potential of an atom is given by the energy change of the system when the atom is added or removed. Since we have chosen free atoms as reference for our bond energies, adding or removing an atom from the system changes its energy exactly by the energy of the bonds that the atom forms when entering the ensemble. In thermodynamic average in a random alloy, the chemical potential can thus be simply determined from the perfect bond fractions in Eq. (4) and the corresponding bond energies in Table 1. This is very useful, since except for the tedious energy-density approach [14], determining quantitative chemical potentials is usually not possible in standard DFT approaches, where only the total energy for all atoms in the system is determined. With the chemical potentials, phase stability can be discussed and important quantities can be determined, such as defect formation or surface energies. Mathematically, the chemical potential of one atom of type A in the alloy can be determined from all bonds that connect to at least one A atom. Considering that all homo but only half of the hetero bonds belong to A,
NN fcc fcc fcc fcc fcc fcc f ðxA ; xB ; xC Þ sbcc 2xAA EAA þ xAB Efcc AB þ xAB EAB 4xA sfcc sbcc
NN bcc bcc bcc bcc bcc sfcc f ðxA ; xB ; xC Þ bcc 2xAA EAA þ xbcc þ AB EAB þ xAB EAB 4xA sfcc sbcc
(6)
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
411
Fig. 10. Chemical potentials in Cr-Ni and Ni-Mo binary alloys, calculated from Eq. (6) using the bond energies in Table 1 and the bond order parameter functions from Eq. (3). Blue, purple, and red lines are chemical potentials of Cr, Mo, and Ni as a function of composition, respectively. The green dots are the average energy of all atoms in the system. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
Fig. 9. (a) Ni-Cr phase diagram. Black lines are CALPHAD results, redrawn after Ref. [28]. (b) Cr-Mo phase diagram. Black lines are redrawn after Ref. [29]. In both panels, red lines are phase boundaries determined from sign change in free-energy of mixing, calculated from the sum of the enthalpy of mixing from the bond-energy model results (Fig. 6(b)) and entropy of mixing approximated by the Stirling formula. For the circles, the full vibrational entropy has been added, calculated from harmonic phonons as described in Ref. [30]. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
The resulting chemical potentials for Cr-Ni and Ni-Mo are shown in Fig. 10. The switching functions (Fig. 3) already contain most of the phase stability information, such as that fcc and bcc alone rule in either of the flat regions of the sigmoid. The non-flat regions of the sigmoids indicate the transition from one to the other, resulting in the undulation of the otherwise linear chemical potentials for close-to-equiatomic concentrations, and it is a priori not clear what meaning the transition region has thermodynamically. As suggested in Kirchheim's work on the extended Gibbs adsorption isotherm [31,32], where the necessary lowering of defect formation energies by solute segregation in alloys is shown and the potential meaning of negative formation energies is discussed, a useful way to probe the possible existence of the transition-region “martensite” consists in using the derived chemical potentials and atomic energies to calculate vacancy formation energies, which we do for the example of CrNi alloys. We show that our approach also offers a solution to one of the most confusing problems for multicomponent alloys, which is to determine the appropriate chemical potential needed to calculate the formation energy of a defect. The formation energy of a defect in a multicomponent (here, as an example, a vacancy in a binary) alloy is in general given by
Ef ðV; An Bm Þ ¼ EðV; An Bm Þ nmA nmB ;
(7)
where E(V, AnBm) is the energy of a cell with one vacancy, n A, and m
B atoms. In other words, the formation energy is the difference between the defective cell and a perfect cell with the same stoichiometry, i.e., the same number of n and m atoms, which of course is usually not known, since n and m, once a vacancy is introduced, do not add up to the stoichiometry of a perfect periodic cell, but are one atom short of what can be calculated. Therefore, usually the energy of either a perfect Anþ1Bm or AnBmþ1 cell is taken as the basis, and corrected for by the chemical potential of the extra atom, which is usually determined from the assumption that it is equal to the energy of an atom in elemental A or B material as needed, which then is labeled as A or B-rich case, with frequently no clear discussion which choice would be most physical, or what A or B rich even means for an alloy [33]. In our case, the procedure is rather clear, since from Eq. (6), the chemical potentials can simply be calculated for the exact stoichiometry needed in Eq. (7) with no ambiguity. We note that interstitials, or defects with any other type of stoichiometry (including multidimensional defects such as dislocations or grain boundaries) can be calculated in an analogous way. Also, in order to minimize the overall error, E(Anþ1Bm) or E(AnBmþ1) can and probably should be used in Eq. (7) and only the stoichiometry difference expressed by the chemical potential. The next question is, what configuration one should choose for the defective structure, since from an Anþ1Bm or AnBmþ1 SQS supercell, any of the A or B atoms respectively can be removed resulting in the same stoichiometry. We decided here to pick the lowest-energy configuration which we determined by running a Monte-Carlo energy minimization [34,35] with respect to the distribution of A and B atoms and the vacancy, as will be described in more detail elsewhere. This assumes that in the random alloy, the atomic mobility around a vacancy (which mediates transport in these metal alloys) is high enough that the local configuration reaches the thermodynamically most favorable state, while the long-range mobility is slow enough that the rest of the crystal remains a random alloy. This is a sensible, but not the only possible choice, and depends on the situation of interest. The resulting energies as a function of concentration are shown in Fig. 11, in comparison to the previously determined Ni-Cr switching function from Eq. (2). It can be seen that the formation energies from both ends decrease towards the center, as suggested by Kirchheim's work [31,32], and assume negative values in the transition region of the switching function, which has its inflexion point at 55.6% Cr. Since, once chemical potentials are linked to stoichiometry as is the case here, negative formation energies indicate an unstable bulk crystal for dilute vacancy concentrations [31e33], the vacancy formation energy within this approach can be
412
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
Fig. 11. (a) Formation energies for vacancies in Ni-Cr alloys, calculated from MCminimized alloy-vacancy configurations and ideal chemical potentials as defined in Eq. (6). (b) Ni-Cr bond switching function from Eq. (3).
used to estimate phase stability limits within alloy systems. Thus, the negative formation energies suggest that there should be a concentration gap between fcc and bcc phases (at a minimum at low temperatures), which the eutectic Ni-Cr phase diagram indeed shows in the form of the central two-phase region, whose width of course is determined by the full free energies rather than the zerotemperature formation energies discussed here. This then gives an interpretation to the transition region of the switching function, which thus would not indicate some “mixed intermediate” martensite-like structure between bcc and fcc, but instead a concentration regime of unstable single-phase bulk phases and thus a multiphase region of the phase diagram. This example thus illustrates that the chemical potentials from the bond-alloy approach enable exciting new angles for a more in-depth examination of the interdependencies between formation energies and phase stability. Since our approach is easily extendible to more-than-two component systems (and in fact here already includes the ternary Ni-CrMo system), it can be extended to multicomponent alloys in a straightforward way where point defects in many-component systems are an active field of research [36,37]. 5. Conclusions In summary, we have introduced a novel way to parameterize alloy energies in the form of a bond-order bond energy model. There, we have introduced a bond order function that models the switching between competing phases. The bond order function switches the respective bond energies for the different phases on and off, which allows to have a continuous model through the phase instability composition. For the case of the Ni-Cr-Mo alloys investigated here, the c/a ratios pffiffiffi calculated from 32-atom cubic SQS cells vary from 1 for fcc to 1 = 2 for bcc, and their function dependence suggests using a sigmoidal switching function, which we chose to describe by an error function. We have applied this model to the Ni-Cr-Mo alloy as a test system, since its stable crystal structures fluctuate between fcc and bcc. We find that our bond-order bond energy model describes the DFT energetics in this system surprisingly well and allows extracting well validated bond energies for both fcc and bcc phases that without concentration-dependent adjustment result in an excellent fit throughout the entire composition range.
We have also shown that for some elements in alloys, in our case specifically Ni, their average atomic charge density can vary considerably as a function of composition. This “charge compression” explains why Ni has lower bcc than fcc bond energies in the examined alloys. It is important to note that the higher charge density that lowers the bcc bond energy is an alloying effect, since the Ni atoms attract additional electrons from the other alloying elements while the alloy structure restricts their resulting volume expansion. In applying the fitted BOBE model, we have shown that, besides the obvious usefulness of bond energies for questions that relate bond energies to processes such as those involved in field ion emission or corrosion, bond energies can be used to determine phase diagrams which, once vibrational entropy is added, show surprisingly good agreement with experimentally derived phase diagrams, as well as quantitative, composition-dependent chemical potentials. From vacancy calculations for binary NiCr alloys for configurations that have been configurationally energy-minimized with Monte Carlo simulations, we could show that the regime of negative formation energies, a sign for thermodynamic instability of the underlying crystal, suggests the existence of a two-phase concentration range in the phase diagram. This suggests strongly that these concentration-dependent chemical potentials can be used to study thermodynamics, defects, and phase stability within a holistic approach since there is a fully quantitative link between stoichiometry and chemical potentials, which opens exciting possibilities for future work. Acknowledgements This work was funded by the Center for Performance and Design of Nuclear Waste Forms and Containers, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award # DE-SC0016584. C.O. also gratefully acknowledges support from the Alexander von Humboldt Foundation through the Feodor-Lynen research fellowship. Our work was supported by allocation of computing time from the Ohio Supercomputer Center under grant number PAS0072. References [1] D.A. Kofke, E.D. Glandt, Monte Carlo simulation of multicomponent equilibria in a semigrand canonical ensemble, Mol. Phys. 64 (6) (1988) 1105e1131. [2] D.R. Gaskell, Introduction to the Thermodynamics of Materials, Taylor and Francis, New York, 2008. [3] R. DeHoff, Thermodynamics in Materials Science, Chapman and Hall/CRC, 2006. [4] A. Agrawal, R. Mishra, L. Ward, K. Flores, M.,W. Windl, An embedded atom method potential of beryllium, Model. Simul. Mater. Sci. Eng. 21 (8) (2013), 085001. [5] R. Mishra, O.D. Restrepo, P.M. Woodward, W. Windl, First-principles study of defective and nonstoichiometric Sr2FeMoO6, Chem. Mater. 22 (22) (2010) 6092e6102. [6] D. Sen, W. Windl, Ab-initio study of energetics of the Si(001)-LaAlO3 interface, J. Comput. Theor. Nanosci. 4 (1) (2007) 57e64. [7] A. Zunger, First-principles statisitical mechanics of semiconductor alloys and intermetallic compounds, in: P.E.A. Turchi, A. Gonis (Eds.), Statistics and Dynamics of Alloy Phase Transformations, Springer, New York, 1994, pp. 361e420. [8] D. De Fontaine, Cluster Approach to order-disorder transformations in alloys, Solid State Phys. 47 (1994) 33e176. [9] A. van de Walle, M. Asta, G. Ceder, The alloy theoretic automated toolkit: a user guide, Calphad 26 (4) (2002) 539e553. [10] A. van de Walle, Multicomponent multisublattice alloys, nonconfigurational entropy and other additions to the Alloy Theoretic Automated Toolkit, Calphad 33 (2) (2009) 266e278. [11] P. Marcus, On some fundamental factors in the effect of alloying elements on passivation of alloys, Corros. Sci. 36 (12) (1994) 2155e2158. [12] L. Yao, T. Withrow, O.D. Restrepo, W. Windl, E.A. Marquis, Effects of the local structure dependence of evaporation fields on field evaporation behavior, Appl. Phys. Lett. 107 (24) (2015) 241602. [13] C. Oberdorfer, T. Withrow, L.J. Yu, K. Fisher, E.A. Marquis, W. Windl, Influence
C. Oberdorfer, W. Windl / Acta Materialia 179 (2019) 406e413
[14] [15]
[16]
[17] [18] [19] [20] [21] [22] [23] [24]
[25] [26]
of surface relaxation on solute atoms positioning within atom probe tomography reconstructions, Mater. Char. 146 (2018) 324e335. N. Chetty, R.M. Martin, First-principles energy density and its applications to selected polar surfaces, Phys. Rev. B 45 (11) (1992) 6074e6088. B. Lee, D.R. Trinkle, Energetics of rutile TiO2 vicinal surfaces with 〈001〉 steps from the energy density method, J. Phys. Chem. C 119 (32) (2015) 18203e18209. M.S. Daw, M.I. Baskes, Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals, Phys. Rev. B 29 (12) (1984) 6443e6453. F. Ercolessi, M. Parrinello, E. Tosatti, Simulation of gold in the glue model, Philos. Mag. A 58 (1) (1988) 213e226. J. Tersoff, Empirical interatomic potential for carbon, with application to amorphous carbon, Phys. Rev. Lett. 61 (25) (1988) 2879e2882. D.W. Brenner, The art and science of an analytic potential, Phys. Status Solidi b 217 (1) (2000) 23e40. D.G. Pettifor, I.I. Oleinik, Analytic bond-order potentials beyond TersoffBrenner. I. Theory, Phys. Rev. B 59 (13) (1999) 8487e8499. R. Becker, Die Keimbildung bei der Ausscheidung in metallischen Mischkristallen, Ann. Phys. 424 (1e2) (1938) 128e140. E.C. Bain, The nature of martensite, Trans. AIME 70 (1924) 25e35. M. Schmidt, H. Lipson, Eureqa (software). www.nutonian.com, 2014. (Accessed 5 May 2018). A. van de Walle, P. Tiwary, M. de Jong, D.L. Olmsted, M. Asta, A. Dick, D. Shin, Y. Wang, L.Q. Chen, Z.K. Liu, Efficient stochastic generation of special quasirandom structures, Calphad 42 (2013) 13e18. T.L. Hill, Chapter 6. Distribution Functions and the Theory of the Liquid State, Statistical Mechanics, McGraw-Hill, New York, 1956. W. Windl, M.S. Daw, Predictive process simulation and Ab-Initio calculation of
[27]
[28] [29] [30] [31]
[32]
[33] [34] [35] [36]
[37]
413
the physical volume of electrons in silicon, in: Technical Proceedings of the 2002 International Conference on Modeling and Simulation of Microsystems: Quantum Effects, Quantum Devices and Spintronics, 2002, pp. 494e497. G.S. Frankel, J.D. Vienna, J. Lian, J.R. Scully, S. Gin, J.V. Ryan, J. Wang, S.H. Kim, W. Windl, J. Du, A comparative review of the aqueous corrosion of glasses, crystalline ceramics, and metals, npj Mater. Degrad. 2 (1) (2018) 15. P. Nash, The Cr-Ni (Chromium-Nickel) system, Bull. Alloy. Phase Diagrams 7 (5) (1986) 466. B. Predel, Cr-Mo (Chromium-Molybdenum), in: O. Madelung (Ed.), Cr-Cs e Cu-Zr, Springer Berlin Heidelberg, Berlin, Heidelberg, 1994, pp. 1e5. W. Luo, W. Windl, First principles study of the structure and stability of carbynes, Carbon 47 (2) (2009) 367e383. R. Kirchheim, Reducing grain boundary, dislocation line and vacancy formation energies by solute segregation. I. Theoretical background, Acta Mater. 55 (15) (2007) 5129e5138. R. Kirchheim, Reducing grain boundary, dislocation line and vacancy formation energies by solute segregation: II. Experimental evidence and consequences, Acta Mater. 55 (15) (2007) 5139e5148. R. Mishra, O. Restrepo, D.,A. Kumar, W. Windl, Native point defects in binary InP semiconductors, J. Mater. Sci. 47 (21) (2012) 7482e7497. C. Niu, W. Windl, M. Ghazisaeidi, Multi-Cell Monte Carlo Relaxation method for predicting phase stability of alloys, Scr. Mater. 132 (2017) 9e12. C. Niu, Y. Rao, W. Windl, M. Ghazisaeidi, Multi-Cell Monte Carlo Method for Phase Prediction, 2018, 04092 eprint arXiv:1811, arXiv:1811.04092. S.C. Middleburgh, D.M. King, G.R. Lumpkin, M. Cortie, L. Edwards, Segregation and migration of species in the CrCoFeNi high entropy alloy, J. Alloy. Comp. 599 (2014) 179e182. Z. Wang, C.T. Liu, P. Dou, Thermodynamics of vacancies and clusters in highentropy alloys, Phys. Rev. Mater. 1 (4) (2017), 043601.