- Email: [email protected]
Interatomic potentials for Zirconium Diboride and Hafnium Diboride
Computational Materials Science 50 (2011) 2828–2835
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Interatomic potentials for Zirconium Diboride and Hafnium Diboride Murray S. Daw a,⇑, John W. Lawson b, Charles W. Bauschlicher Jr. c a
Dept of Physics & Astronomy, Clemson University, Clemson, SC 29634, USA MS 234-1, NASA Ames Research Center, Moffett Field, CA 94035, USA c Space Technology Division, MS 230-3, NASA Ames Research Center, Moffett Field, CA 94035, USA b
a r t i c l e
i n f o
Article history: Received 4 November 2010 Received in revised form 25 April 2011 Accepted 26 April 2011 Available online 31 May 2011 Keywords: Interatomic potential Tersoff potential Zirconium Diboride Hafnium Diboride
a b s t r a c t We report on the first interatomic potentials for Zirconium Diboride and Hafnium Diboride. The potentials are of the Tersoff form, and are obtained by fitting to a first-principles database of basic properties of elemental Zr, Hf, B, and the compounds ZrB2 and HfB2. Two variants of the Zr–B potentials have been obtained, and one for Hf–B. The potentials have been tested against a variety of properties of the compound, with the conclusion that they are stable and provide a reasonable representation of the desired properties of the two diborides. Ó 2011 Elsevier B.V. All rights reserved.
1. Introduction and motivation Ceramic borides, carbides, and nitrides (referred to generically here as Ultra High Temperature Ceramics, or UHTCs) comprise a family of materials characterized by high melting point, chemical inertness, good strength and good oxidation resistance [1,2]. They are candidate materials for both atmospheric reentry and hypersonic vehicles. In particular, high temperature materials are needed for sharp leading edges of hypersonic aircraft. Such edges will enable development of safer, more maneuverable aircraft. For these applications, the materials must maintain good properties and operate with no or limited oxidation or ablation in high temperature, reactive environments. Industrial applications for UHTCs are also of interest and include use in foundries, molten metal crucibles, and nuclear reactors. While initial studies of boride compounds such as UHTCs began in the 1960s [3], progress is still needed on a range of issues related to processing, properties, and performance of these materials. For example, thermal and chemical stability provide challenges for all these issues and also presents an obstacle for development of economical manufacturing. Better understanding of the effect of processing on microstructure and ultimately on properties are needed as well as a fuller elucidation of the basic mechanisms responsible for the response of these materials to extreme environments. Further innovations in the materials will require thorough knowledge of their properties and behavior under diverse conditions of temperature, environment and stress. ⇑ Corresponding author. E-mail address: [email protected] (M.S. Daw). 0927-0256/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2011.04.038
Computational modeling of these borides has been limited in part by the absence of appropriate interatomic potentials. While potentials exist for some carbides and nitrides, and there has been wide experience with their application, potentials for borides have been relatively less explored, especially metal borides. To our knowledge, until this work, there have been no interatomic potentials for ZrB2, HfB2 or any similar material. In the development of interatomic potentials, care must be exercised in the choice of semi-empirical form. A common form for carbides and nitrides is the Tersoff potential, by which one seeks to represent something of the angular dependence of whatever covalent bonding is present. Tersoff-style potentials have been developed for SiC, which is a common component of UHTCs. In this paper, we report on our development of the first interatomic potentials for ZrB2 and HfB2. We have chosen to develop Tersoff-style potentials for the diborides, based upon two considerations. First, the crystal structure (C32) exhibits closed-packed sheets of the transition metals, reminiscent of the HCP structure of the metals, interleaved with hexagonal B sheets. It seems quite possible that the bonding within the B sheets may show some degree of covalency. Second, UHTCs are frequently formed as composites of metal-diborides and SiC, and the latter already has well-established Tersoff-style potentials. The ultimate goal of this line of work would be to investigate some of the behavior of composites. For this reason, it would be desirable to combine older existing SiC potentials with the new potentials, and this is made most convenient if the new potentials already share the format of the older potentials. The task of the present work is to develop Tersoff-style interatomic potentials for ZrB2 and HfB2 which will be suitable for
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835
furthering the understanding of its mechanical and thermal properties. Development of potentials for compounds is a task which involves many intermediate steps and decisions, and it is appropriate to keep in mind that our aim is to study Zirconium and Hafnium Diborides in the context of UHTC applications. The Tersoff potentials of compounds contain a large number of adjustable parameters, and progress is easiest if we subdivide the large parameter space. We have therefore pursued a strategy of isolating subsets of the parameters in order to obtain some control of the fitting process. Our general strategy is to obtain first a Tersoff potential for elemental Zr and Hf, then a Tersoff potential for elemental B, and then determine the cross-type parameters for the compounds. Target properties for the elements include lattice constants, cohesive energies, elastic constants, and vacancy energies in a variety of crystal structures. For the compound we target the lattice constants and cohesive energy. The transition metals present no special challenge in this matter. There exist already several potentials for Zr — for example, those by Willaime and Massobrio [4,5] — and there are good reasons to have some level of confidence in these potentials. It is also convenient that the form of the potential by Willaime and Massobrio can be easily converted to Tersoff-form, as will be described below. We have also obtained a new Zr potential for this work by optimizing to the current database, as well as a new Hf potential. However, the story for the trivalent metalloid (B) is quite different [6–10]. There are several polymorphs established for elemental boron, and there are two phases which compete for stability at ambient conditions. Both of the two low-energy structures have large unit cells and complex structure. The subtle bonding variations seem to occur because of competition between metal and insulating behavior, with some covalent bonding likely in other cases. It seems unlikely that interatomic potentials will be able to capture such a broad range of behavior, which may also account for the dirth in the literature of reports on interatomic potentials for B. We have therefore focused at this stage on obtaining potentials which capture the gross behavior over a range of coordination, as will be described in more detail in the next section. It is especially important, then, to remember the context of this work: we are trying to develop potentials which are to be used in a relatively narrow range of circumstance (within the diboride) and only for a narrow range of properties (mechanical and thermal). The assertion is then that the behavior of B in this narrowly defined context might be represented with sufficient success by interatomic potentials. It is anticipated that further refinements of the potential may be necessary, with this first cut serving as a basis for evaluating possible weaknesses. Once the elemental Zr, Hf and B potentials are established, we use the properties of the compound ZrB2 (HfB2) to determine the cross-terms in the Tersoff potential. Because we have two, separate Zr potentials (one from the literature and one original to this work), we do this determination twice, so that we have two separate potentials for ZrB2. Only one potential for HfB2 is obtained here. We then test the two interatomic potential sets by comparing to the results of first principles calculations for a variety of properties that were not targets for the fitting procedure. The tests include elastic constants, point defect energies, phonon dispersion curves, and stability of two naturally occuring boron structures. Overall the results are reasonable, showing that the potentials are stable and generally reliable. Our conclusion is that these first potentials are ready to be used in more sophisticated calculations of mechanical and thermal properties. This paper is organized as follows. In the next section, we give the details of the determination of the interatomic potential, in terms of the specifics of the fitting process, including properties that were used as targets for the fitting. In the next section we
2829
describe the independent testing that we subjected the potentials to, as well as the results. Finally, we discuss the qualities of the potentials, as well as suggestions for the direction of future refinements. 2. Method As mentioned in the Introduction, our general strategy for determining a Tersoff potential for ZrB2 (HfB2) has been first to obtain functions for element Zr (Hf) and B, and then to obtain the cross-element parameters from properties of the compound. Our approach has been to perform first-principles calculations of specific properties to serve as targets for the fitting. In this section, we describe the generation of targets for fitting, then details of the fitting for elemental Zr and Hf, elemental B, and finally the compounds. 2.1. Generation of targets for fitting Target values for fitting were generated from ab initio calculations based on density functional theory (DFT). Results from a more comprehensive set of ab initio calculations for ZrB2 and HfB2 will appear in a separate publication [11] as will further computational details. We summarize briefly some of those details here. Computations were performed with the plane wave DFT codes VASP [12,13] and ABINIT [14,15]. While there is considerable overlap between the two codes, each also offers unique features. Using different codes also allows us to check our results. FHI norm-conserving pseudopotentials were used with ABINIT and the projected augmented wave (PAW) potentials with VASP. The functional of Perdew, Burke, and Ernzerhof (PBE) [16] was used throughout. For pure Zr (Hf), cohesive energies, lattice constants, elastic constants and vacancy energies were calculated for HCP, FCC and BCC crystal structures. For pure B, the only stable structures known have quite complex unit cells. Rather than deal with the difficulties of fitting to properties of structures with large unit cells, we have focused on the properties of simple structures (SC, FCC, hexagonal sheet, and triangular sheet). For those, we have calculated cohesive energies, lattice constants, and bulk moduli. Note that the full set of elastic constants for B was not calculated because for these simple structures, B is unstable to shear. For the diborides, the cohesive energies and lattice constants were obtained as targets for fitting. Lattice constants were obtained by energy minimization. Cohesive energies were derived by evaluating the total energy at the optimized geometry and subtracting the atomic energies. Elastic constants were obtained by the linear response method [17,18] as implemented in ABINIT and the finite difference method as implemented in VASP. We expect the linear response method to be slightly more accurate than finite differences since it does not rely on convergence with respect to supercell. All quantities were converged individually with respect to k-space sampling and plane wave cutoff. 2.2. Definition of the potentials Because some variations of the Tersoff potential for compounds have emerged in the literature [19–24], we make explicit in the following the forms we are using for this work. The original Tersoffstyle form for multi-component systems [23] was augmented by Nord et al. [25] to include the c factors in the calculation of the bond order. These factors alter the way in which a third atom enhances or weakens the bond between a pair. For example, the factor cBBZr alters the way that Zr affects the B–B bond, making it
2830
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835
different than the way a B atom would affect the bond. These c factors were found particularly helpful in fitting to the properties of ZrB2. The total energy for a configuration of atoms in a multi-component system is given by
E¼
X
Ei ¼
i
1X V ij 2 i–j
ð1Þ
where the pair term is
h i ½ij ½ij ½ij V ij ¼ fC ðdij Þ fR ðdij Þ þ bij fA ðdij Þ
ð2Þ
The repulsive piece is given by ½ij
fR ðdÞ ¼ Aij expðkij dÞ
ð3Þ
and the attractive piece by ½ij
fA ðdÞ ¼ Bij expðlij dÞ
ð4Þ
These functions are restricted to short-range by this rather abrupt cut-off function
8 d < Rij > < 1; ½ij fC ðdÞ ¼ 12 þ 12 cos½pðd Rij Þ=ðSij Rij Þ; Rij < d < Sij > : 0; Sij < d
ð5Þ
The bond-order function is n
n
bij ¼ ð1 þ bi i fiji Þ1=2ni
ð6Þ
which uses the parameter
fij ¼
X
fC ðr ik Þcijk g i ðhijk Þ expf½k3i ðdij dik Þmi g ½ij
ð7Þ
k–i;j
The angular dependence is supplied by
h i 2 2 g i ðhÞ ¼ 1 þ c2i =di c2i = di þ ðhi cos hÞ2
ð8Þ
Note that parameters such as Aij depend only on the types of atoms i and j (not on positions). The parameters that depend only on a single type are b, n, c, d, h, m, and k3. The parameters that depend on two types are A, B, k, l, R, and S. It is often taken that the cross-type parameters (e.g., AZrB and kZrB) are restricted to be means of the same-type parameters. For example, kij = (ki + kj)/2 and Aij = (AiAj)1/2. However, we do not make that choice here. Instead, we allow the cross-type parameters to be flexible in the fit. 2.3. Potentials for Zr, Hf Tersoff-style potentials for Zr [26] and Hf [27] have been developed recently as part of a study of their oxides. We have developed here an independent set of Tersoff potentials for Zr and Hf. It is usually argued that metallic bonding in highly coordinated environments does not display significant angular dependence. With that in mind, Embedded Atom Method/Finnis–Sinclair (EAM/FS) [28–30] forms were successfully obtained for Zr by several authors [31,4,5,32]. The FS potentials due to Willaime and Massobrio (WM) [4,5] are of such form that they can be easily converted to Tersoff form (see note by Brenner [33]). We therefore have converted the second WM potential (‘‘WM2’’) to Tersoff form as the first of our candidate Zr potentials which we will refer to as ‘‘Zr-1’’. We will present the parameters and the calculated properties below. We have also used the Zr WM2 potential as a starting point for a re-optimization to the current targets. The second Zr potential is somewhat altered from the WM2 potential and will be referred to here as ‘‘Zr-2’’.
Finally, we used the reoptimized Zr potential as the starting point for the new Hf potential. The target properties we chose for Zr and Hf included a full range of properties related to the close-packed structures FCC and HCP as well as the more open BCC structure which is expected to be similar in energy. To simplify the fitting process, we focused on the FCC form of Zr, which should be quite similar in its properties to those of HCP. We therefore chose to fit to these properties in FCC-Zr: lattice constant, cohesive energy, elastic constants, vacancy energy, as well as the pressure derivative of the bulk modulus. In addition we fitted to the cohesive energies of HCP and BCC phases, so ensure a good representation of their relative stabilities. As has been noted in many places (see, for example, Willaime and Massobrio), the FCC-BCC energy difference is sensitive to cutoff. We also have chosen to include 2nn in the FCC structure in the current potentials. The new set of Zr and Hf parameters were obtained by doing a least squares fit to the target data. We chose initial values for the Zr parameters to be those of the WM2 potential. For Hf, we used the newly obtained Zr potential as the starting point. Then the fitting proceeded in stages. The first adjustments were to A, B, k, l, and k3. The second stage allowed b and n to flex. In the third stage, all parameters were allowed to vary. We attempted to increase the value of c, which would have introduced some bond-angle dependence to the functions, but the fit to the target data was always better with this parameter at zero. In Tables 1 and 2 we show the fidelity of the Zr and Hf potentials in representing the chosen target properties. All potentials are reasonably good. They are all relatively weak in representing the derivative of the bulk modulus, which might be expected from a 2nn potential. This weakness should not be fatal in the present application. All three of the parameter sets (two for Zr and one for Hf) were subjected to simple stability tests, such as relaxing around a vacancy, adding significant random distortions to the FCC structure and relaxing, and running molecular dynamics at elevated temperature. All sets gave satisfactory results. The resulting parameters for the potentials are displayed in Table 3. We note that we also tried several other Zr forms, with shorter cut-offs and also allowing for angular dependence, but these EAM/FS-like potentials produced the most satisfactory results.
2.4. Potential for B One consideration in developing interatomic potentials is to ensure stability, at least within the range of environments that the potential is likely to see service. Boron is challenging for the usual fitting procedure, because all of the simple structures (hexagonal
Table 1 The target properties and the Tersoff values for Zr calculated using the parameters WM2 [4,5] (‘‘Zr-1’’) and the parameters derived in this work (‘‘Zr-2’’). B is the bulk modulus and B0 here is the derivative of the bulk modulus w.r.t. lattice constant at equilibrium. Vf is the vacancy formation energy. Property (units)
Target
‘‘Zr-1’’ (WM2)
‘‘Zr-2’’ (new)
a0(FCC) (Å) E0(FCC) (eV) B(FCC) (eV/Å3) B0 (FCC) (eV/Å4) C11(FCC) (eV/Å3) C12(FCC) (eV/Å3) C44(FCC) (eV/Å3) Vf(FCC) (eV) a0(HCP) (Å) E0(HCP) (eV) E0(BCC) (eV)
4.530 6.160 0.578 0.86 0.774 0.481 0.356 2.500 3.230 6.180 6.050
4.532 6.127 0.573 1.95 0.712 0.503 0.347 2.211 3.256 6.142 6.124
4.547 6.174 0.573 1.93 0.664 0.527 0.347 2.224 3.215 6.174 6.006
2831
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835 Table 2 The target properties and the Tersoff values for Hf calculated using the new parameters derived in this work. (See the caption for Table 1 for a description of the properties.) Property (units)
Target
New pot
a0(FCC) (Å) E0(FCC) (eV) B(FCC) (eV/Å3) B0 (FCC) (eV/Å4) C11(FCC) (eV/Å3) C12(FCC) (eV/Å3) C44(FCC) (eV/Å3) Vf(FCC) (eV) a0(HCP) (Å) E0(HCP) (eV) E0(BCC) (eV)
4.470 6.410 0.649 0.92 0.911 0.518 0.462 2.580 3.190 6.480 6.300
4.570 6.469 0.606 1.92 0.766 0.536 0.347 2.296 3.231 6.469 6.296
Table 3 Parameters describing the two different potentials for Zr and one potential for Hf. The ‘‘Zr-1’’ is from the work of Willaime and Massobrio. ‘‘Zr-2’’ and ‘‘Hf’’ are derived for this work. See also Brenner’s note on the relationship of various interatomic potential forms [33]. Parameter
‘‘Zr-1’’ (WM2)
‘‘Zr-2’’ (new)
‘‘Hf’’ (new)
A (eV) B (eV) k (Å1) l (Å1)
3.65374 103 3.68286 101 2.92969 0.661542
3.58787 103 3.71035 101 2.91078 0.659036
3.15839 103 3.96250 101 2.83165 0.65895
b n c d h
1 1 0 1 0
1 1 0 1 0
1 1 0 1 0
R (Å) S (Å)
4.88 5.18
4.6 5.0
4.6 5.0
m k3 (Å1)
1 1.32308
1 1.32308
1 1.32308
and triangular sheets, as well as SC, BCC, and FCC) are unstable to shear. While it is certainly possible to fit to unstable structures, one must surrender the assurance given by the usual stability tests like we did for Zr in the previous section. In our present approach, we attempted to separate the properties and parameters according to whether they relate to any possible angular dependence of B bonding. This can be achieved in the way suggested by the recent work of Liang, Phillpot, and Sinnott [34]. They noted that the 1nn Tersoff form can effectively factor into two parts, with the separation being afforded by the combined factor of the repulsive coefficient B and the bond order parameter b. Taking this product, bB, and noting that b depends only on the bond angles present in a given crystal, it is then possible to separate the bond-angle and bond-length dependence. This is carried out by fitting to the lattice constant, cohesive energy, and bond stiffness for four simple structures: hexagonal sheet, triangular sheet, SC, and FCC. This is especially simple if we assume 1nn form, as we have here. The fitting determines the parameters A, k, l, and the value of bB in each of the four structures. The result of this fit is shown in Table 4, and the same results are plotted Figs. 1–3, where one can see that the quality of the fit is reasonable. In Fig. 1, one can see that the cohesive energies are reproduced to within 0.1 eV and the bondlengths to within 0.1 Å. In Figs. 2 and 3, we can see that the trends of bond strength and bond stiffness are correctly captured. To complete the fit, the resulting bond order parameters must be related to the Tersoff form in Eq. (6), which incorporates some angular dependence as well as coordination dependence. However, the resulting bond order parameters show a simple dependence on
Table 4 The target properties and the fitted values for boron, used to establish the (bB) factors for the listed structures. Lattice constants are in Å, cohesive energies in eV, and bond stiffnesses in eV/Å2. [Note that these are calculated using 1nn only in each structure. Using the Tersoff potential with cut-off fixed by distance — as would ordinarily be done — may yield somewhat different answers.]. Structure
Property
Target
Fit
Hexagonal sheet
a0 E0 E000
2.91 5.15 11.35
2.89 5.08 7.98
Triangular sheet
a0 E0 E000
1.70 5.71 21.73
1.81 5.75 27.06
SC
a0 E0 E000
1.88 5.33 24.50
1.84 5.21 24.51
FCC
a0 E0 E000
2.86 5.07 21.85
2.84 5.22 12.28
Fig. 1. Comparison of fit to target for cohesive energy vs. bondlength for B in four structures. The target data are represented by circles and the fits by squares. The four structures are (in order of increasing bondlength): hexagonal sheet, triangular sheet, SC, and FCC.
1
Fig. 2. Comparison of fit to target for bondstrength vs. bondlength for B in four structures. The target data are represented by circles and fits by squares. The four structures are (in order of increasing bondlength): hexagonal sheet, triangular sheet, SC, and FCC.
coordination and no evidence of angular bonding, as shown when we plot the bond order parameters vs. coordination in Fig. 4. We therefore set c = 0 to eliminate the angular dependence and fitted the remaining coordination dependence using b and n. The quality of the fit is illustrated in the figure. The resulting parameters for B are shown in Table 5, which are very simple. As noted above, it is not possible to subject these B parameters to the same sort of stability tests as Zr.
2832
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835 Table 6 The target properties and the Tersoff values for ZrB2 (C32 structure) calculated using the two alloy sets obtained here. Property (units)
Target
Pot 1
Pot 2
a0(C32) (Å) c0(C32) (Å) E0(C32) (eV)
3.170 3.550 21.70
3.140 3.547 21.54
3.131 3.484 21.35
The procedure was repeated for HfB2. The targets and the fits are displayed in Table 7. The resulting cross-type Tersoff parameters for the alloys are presented in Table 8. Fig. 3. Comparison of fit to target for bond stiffness vs. bondlength for B in four structures. The target data are represented by circles and fits by squares. The four structures are (in order of increasing bondlength): hexagonal sheet, triangular sheet, SC, and FCC.
3. Tests for the compound potentials
Parameter
Value
A (eV) B (eV) k (Å1) l (Å1)
1.0619848 103 3.6767392 101 3.955278 1.190167
In order to test independently of the target properties, we have evaluated — in ZrB2 and HfB2 — point defects energies, elastic constants and phonon spectra, as well as energy vs. lattice constant for large strains. All three sets of interatomic potentials were tested. First, we compared certain point defect energies as predicted by the interatomic potentials to the results of first-principles. In computing energies of defects with varying stoichiometry, it is necessary to be somewhat cautious. The interatomic potentials are by nature constructed with respect to isolated atoms as the reference of energy. The first-principles calculations can be defined with respect to the same reference, but that involves a calculation of the atomic energies for which there is some uncertainty. The most secure way around this ambiguity is to compare sets of point defects that preserve overall stoichiometry. We have thus computed the energy of a pair of anti-sites (interchanging distant atoms of different type) as well as a trio of vacancies (removing one M and two B from three distant sites), and comparing this to the cohesive energy per unit cell of the compound. The results for the point defects are presented in Tables 9 and 10. In calculating the concentrations of various point defects in the stoichiometric compound, certain combinations of the point defect energies occur. Those combinations which have thermodynamic significance for the stoichiometric compounds we present in Tables 11 and 12. The defect energies are themselves quite large, and the differences between interatomic potentials and first-principles are also significant. However, it is not expected that the energies of point defects such as these are likely to be important in the properties for which these potentials are intended. With that in mind, it seems that the point defect energies are acceptable. The next test was to compare elastic constants and phonon spectra obtained from the potentials versus ab initio computations. As mentioned previously, the ab initio results will be discussed in more detail in a forthcoming publication [11]. In brief, for the ab initio calculations, both the elastic constants and the phonon spectra were calculated using linear response theory as implemented in ABINIT. As with the target calculations, FHI pseudopotentials and the PBE functional were used. All quantities were individually converged with respect to k-space sampling and plane wave cutoff. Elastic constants were computed both with and without internal relaxation; however, there is only a very small amount of relaxation.
b n c d h
0.541 21.885575 0 1 0
Table 7 The target properties and the Tersoff values for HfB2 (C32 structure) calculated using the alloy set obtained here.
R (Å) S (Å)
2.2 2.5
m k3 (Å1)
1 0.
Fig. 4. Bond order parameter vs. coordination for B. The points are the established values (see text) for four structures, in order of increasing coordination: hexagonal sheet, triangular sheet, SC, and FCC. The curve is our fit using the Tersoff form.
2.5. Potentials for Zr–B and Hf–B Given the two candidate Zr potentials and the new B potential described in the preceding sections, we then moved onto the determination of the alloy parameters for ZrB2. Our approach was to keep fixed those parameters determined to this point and to determine the cross-type parameters from the properties of ZrB2. Three simple properties were chosen for the fitting: the two lattice constants and the cohesive energy of ZrB2 in the C32 structure (see Table 6). These properties were fit easily by varying primarily the c factors, though the other cross-type parameters were also adjusted.
Table 5 The parameters of the Tersoff form from the new fit to boron described in the text.
Property (units)
Target
Current
a0(C32) (Å) c0(C32) (Å) E0(C32) (eV)
3.15 3.52 22.47
3.13 3.42 22.39
2833
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835 Table 8 The cross-type parameters of two Tersoff forms for the Zr–B and Hf–B alloys. Two potentials are shown. Zr–B ‘‘pot 1’’ combines with the ‘‘Zr-1’’ from Table 3, Zr–B ‘‘pot 2’’ with the ‘‘Zr-2’’ parameters, and Hf–B with the Hf parameters. All three diboride potentials use the same B potential (Table 5). Parameter
Zr–B (pot 1)
Zr–B (pot 2)
Hf–B
AMB (eV) BMB (eV) kMB (Å1) lMB (Å1)
1.9008612 103 4.0127474 101 3.3688840 0.9258545
1.7659542 103 4.0023041 101 3.360803944 0.9220569
1.6967221 103 4.2730462 101 3.330828218 0.92455101
RMB (Å) SMB (Å)
3.28 3.60
3.18 3.54
3.18 3.54
cMMM cMMB cMBM cMBB cBMM cBMB cBBM cBBB
1.0 0.1905300 9.5016209 0.9751164 0.9999995 1.0000000 0.1721612 1.0
1.0 0.2930705 9.0588265 3.5295546 0.94069107 0.82400007 0.17523123 1.0
1.0 0.25267505 8.0112269 3.7926378 0.9934950 0.81400090 0.1620252 1.0
Table 13 Calculated elastic constants for ZrB2, comparing first-principles results and the predictions of the two interatomic potentials (‘‘pot 1’’ and ‘‘pot 2’’) developed here. Elastic constants are in GPa. Method
c11
c12
c13
c33
c44
B
G
A(=c33/c11)
ABINIT Zr–B (pot 1) Zr–B (pot 2)
543 422 575
56 156 213
113 171 301
419 320 817
234 119 246
230 240 400
223 118 211
0.77 0.76 1.42
Table 14 Calculated elastic constants for HfB2, comparing first-principles results and the predictions of the interatomic potentials developed here. Elastic constants are in GPa. Method
c11
c12
c13
c33
c44
B
G
A(=c33/c11)
ABINIT Hf–B
578 603
71 224
118 303
430 866
238 259
244 415
231 224
0.74 1.43
K
M
900
Method
E(ZrB)
E(BZr)
E(VZr)
E(VB)
ABINIT pot 1 pot 2
11.2 7.2 8.5
8.6 5.3 5.2
11.2 9.6 9.3
9.3 8.2 8.1
800
Frequency (cm-1)
Table 9 Some (relaxed) point defect energies for ZrB2, comparing two potentials and first principles. Energies are in eV.
E(ZrB)
E(BZr)
E(VZr)
E(VB)
ABINIT new
11.6 8.6
8.4 5.9
10.6 10.2
9.0 8.3
E(VZr) + 2E(VB) + Ecoh
ABINIT pot 1 pot 2
19.80 12.54 13.72
7.80 4.57 4.15
Table 12 Stoichiometric combinations of point defect energies for HfB2, comparing the new potentials and first principles. Energies are in eV. Method
E(HfB) + E(BHf)
E(VHf) + 2E(VB) + Ecoh
ABINIT new
20.01 14.45
6.53 4.46
400 300
0
A
L
H
A
Fig. 5. Normal mode frequencies for ZrB2 calculated using interatomic potentials from ‘‘pot 1’’ (solid lines), compared to those for ABINIT (dashed lines).
900 800 700
Frequency (cm-1)
E(ZrB) + E(BZr)
500
100
Table 11 Stoichiometric combinations of point defect energies for ZrB2, comparing two potentials and first principles. Energies are in eV. Method
600
200
Table 10 Some (relaxed) point defect energies for HfB2, comparing the potential and first principles. Energies are in eV. Method
700
600 500 400 300 200 100 0
K
M
A
L
H
A
Fig. 6. Normal mode frequencies for ZrB2 calculated using interatomic potentials from ‘‘pot 2’’ (solid lines), compared to those for ABINIT (dashed lines).
To calculate the elastic constants for the potentials, we generated a set of five symmetry independent deformations of the crystal parameterized by a strain parameter a [35]. Curves were generated of energy versus a. Elastic constants were obtained by fitting Taylor’s expansions of the energy with respect to the strain parameter. To calculate the phonon spectra for the potentials, the ‘‘frozen phonon’’ method implemented in the code PHON [36] was used. In this method, symmetry independent atomic displacements were
generated and the corresponding forces were calculated using the potentials in conjunction with the LAMMPS [37] code. The dynamical matrices were then constructed, diagonalized and the phonon spectra obtained. Supercells up to 8 8 8 cells gave well-converged values. The calculated elastic constants are presented in Tables 13 and 14. Generally, the results are acceptable. Specifically, the worst agreement is with some of the shear constants. It is reasonable
2834
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835
900
harmonic part is very good for significant strains (up to 20%). For larger strains, the interatomic potentials become somewhat too stiff. We note that the overly stiff behavior under tension (a > a0) is unavoidable with interatomic potentials restricted to first or second nearest-neighbors. This would seem to be related to the over-estimate of the magnitude of B0 for the elemental metals (see Tables 1 and 2.) Similar results are shown for HfB2 in Fig. 9. Overall the results show that these first interatomic potentials for ZrB2 and HfB2 provide a reasonable representation of the properties of the compounds.
Frequency (cm-1)
800 700 600 500 400 300 200
4. Conclusions
100 0
K
M
A
L
H
A
Fig. 7. Normal mode frequencies for HfB2 calculated using interatomic potentials derived in this work, compared to those for ABINIT (dashed lines).
Fig. 8. Energy vs. lattice constant for ZrB2 for large strain, comparing the two interatomic potentials and first-principles results.
We have presented here our derivation of the first interatomic potentials for ZrB2 and HfB2, which are obtained in the Tersoff format. We have determined two variant sets of parameters by fitting to targets generated by first-principles for a variety of properties, including lattice constants, cohesive energies, vacancy energies, elastic constants, and others, for Zr, Hf, B, and the compounds ZrB2 and HfB2. The agreement with the targets is reasonable. We have also determined that the functions appear to provide a reasonable representation of additional properties of the compounds outside of the target sets, including elastic constants, phonon dispersion, and point defects. The potentials have so far passed all tests of stability. If there is an obvious weakness of these first potentials for ZrB2 and HfB2, it would seem to be the complete lack of angular bonding among the B atoms. This result was determined by our database of structures, which included but did not emphasize the hexagonal sheet. However, it is not obvious at this point how the lack of angular dependence should be expressed in the various predicted properties of the compound. It is therefore judged that these potentials represent a reasonable first version. From here, further testing will be carried out which may suggest directions for future refinements. Acknowledgement MSD was supported under a NASA prime contract to ELORET Corporation. References
Fig. 9. Energy vs. lattice constant for HfB2 for large strain, comparing the two interatomic potentials and first-principles results.
to expect that including angular bonding in the boron potentials would significantly improve the shear moduli. However, there is some differences in the quality even between our two Zr–B potentials, neither of which contains angular dependence for the boron. This would seem to be the best place to begin improvements. The calculated phonon spectra for ZrB2 are presented for in Figs. 5 and 6 (for ‘‘Zr-1’’ and ‘‘Zr-2’’, respectively) and for HfB2 in Fig. 7. Generally the results are reasonable. Finally, we test the anharmonic part of the potentials by applying a large hydrostatic strain. For ZrB2, Fig. 8 shows that the
[1] M.J. Gasch, D.T. Ellerby, S.M. Johnson, Ultra high temperature ceramic composites, in: N.P. Bansal (Ed.), Handbook of Ceramic Composites, Kluwer Academic, 2005, pp. 197–224. [2] R. Loehman, E. Corral, H.P. Dumm, P. Kotula, R. Tandon, Ultra high temperature ceramics for hypersonic vehicle applications, Technical Report SAND20062925, Sandia National Labs, Albuquerque, NM, 2006. [3] L. Kaufman, E.V. Clougherty, Investigation of boride compounds for very high temperature applications, Technical Report RTD-TRD-N63-4096 Part III, ManLabs, Inc., Cambridge, MA, 1966. [4] F. Willaime, C. Massobrio, Phys. Rev. Lett. 63 (1989) 2244. [5] F. Willaime, C. Massobrio, Phys. Rev. B 43 (14) (1991) 11653. [6] M.W. Chase Jr., NIST-JANAF Thermochemical Tables, 1998 (issued as J. Phys. Chem. Ref. Data monograph no. 9). [7] B.E. Douglas, S.-M. Ho, Structure and Chemistry of Crystalline Solids, Springer, 2006. [8] M.J. van Setten, M.A. Uijttewaal, G.A. de Wijs, R.A. de Groot, J. Am. Chem. Soc. 129 (9) (2007) 2458. [9] H. Tang, S. Ismail-Beigi, Phys. Rev. Lett. 99 (2007) 115501. [10] A.R. Organov, J. Chen, C. Gatti, Y. Ma, Y. Ma, C.W. Glass, Z. Liu, T. Yu, O.O. Kurakevych, V.L. Solozhenko, Nature 457 (2009) 863. and references therein. [11] J.W. Lawson, C.W. Bauschlicher, M.S. Daw, J. Am. Ceram. Soc. (2011), in press. [12] G. Kresse, J. Furthmller, Phys. Rev. B 54 (1996) 11169. [13] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758. [14] M. Fuchs, M. Scheffler, Comput. Mater. Sci. 119 (1999) 67. [15] X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, Ph. Ghosez, J.-Y. Raty, D.C. Allan, Comput. Mater. Sci. 25 (2002) 478. [16] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [17] X. Gonze, Phys. Rev. B 55 (1997) 10337. [18] X. Gonze, C. Lee, Phys. Rev. B 55 (1997) 10355.
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835 [19] [20] [21] [22] [23] [24] [25] [26] [27]
J. Tersoff, Phys. Rev. Lett. 56 (1986) 632. J. Tersoff, Phys. Rev. B 37 (1988) 6991. J. Tersoff, Phys. Rev. B 38 (1988) 9902. J. Tersoff, Phys. Rev. Lett. 61 (1988) 2879. J. Tersoff, Phys. Rev. B 39 (1989) 5566. J. Tersoff, Phys. Rev. B 49 (1994) 16349. J. Nord, K. Albe, P. Erhart, K. Nordlund, J. Phys.: Cond. Matter 15 (2003) 5649. T. Iwasaki, J. Mater. Res. 20 (2005) 1300. T.-R. Shan, B.D. Devin, T.W. Kemper, S.B. Sinnott, S.R. Phillpot, Phys. Rev. B 81 (2010) 125328. [28] M.S. Daw, M.I. Baskes, Phys. Rev. Lett. 50 (1983) 1285.
[29] [30] [31] [32] [33] [34] [35] [36] [37]
2835
M.S. Daw, M.I. Baskes, Phys. Rev. B 29 (1984) 6443. M.W. Finnis, J.E. Sinclair, Philos. Mag. 50 (1984) 45. D.J. Oh, R.A. Johnson, J. Mater. Res. 3 (33) (1988) 471. R.C. Pasianot, A.M. Monti, J. Nucl. Mater. 264 (1–2) (1999) 198–205. D.W. Brenner, Phys. Rev. Lett. 63 (1989) 1022. T. Liang, S.R. Phillpot, S.B. Sinnott, Phys. Rev. B 79 (2009) 245110. L. Fast, J.M. Wills, B. Johansson, O. Eriksson, Phys. Rev. B 15 (1995) 17431. D. Alfe, Comput. Phys. Commun. 180 (2009) 2522. S.J. Plimpton, J. Comput. Phys. 117 (1995) 1–19. further info at lammps.sandia.gov.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Interatomic potentials for Zirconium Diboride and Hafnium Diboride Murray S. Daw a,⇑, John W. Lawson b, Charles W. Bauschlicher Jr. c a
Dept of Physics & Astronomy, Clemson University, Clemson, SC 29634, USA MS 234-1, NASA Ames Research Center, Moffett Field, CA 94035, USA c Space Technology Division, MS 230-3, NASA Ames Research Center, Moffett Field, CA 94035, USA b
a r t i c l e
i n f o
Article history: Received 4 November 2010 Received in revised form 25 April 2011 Accepted 26 April 2011 Available online 31 May 2011 Keywords: Interatomic potential Tersoff potential Zirconium Diboride Hafnium Diboride
a b s t r a c t We report on the first interatomic potentials for Zirconium Diboride and Hafnium Diboride. The potentials are of the Tersoff form, and are obtained by fitting to a first-principles database of basic properties of elemental Zr, Hf, B, and the compounds ZrB2 and HfB2. Two variants of the Zr–B potentials have been obtained, and one for Hf–B. The potentials have been tested against a variety of properties of the compound, with the conclusion that they are stable and provide a reasonable representation of the desired properties of the two diborides. Ó 2011 Elsevier B.V. All rights reserved.
1. Introduction and motivation Ceramic borides, carbides, and nitrides (referred to generically here as Ultra High Temperature Ceramics, or UHTCs) comprise a family of materials characterized by high melting point, chemical inertness, good strength and good oxidation resistance [1,2]. They are candidate materials for both atmospheric reentry and hypersonic vehicles. In particular, high temperature materials are needed for sharp leading edges of hypersonic aircraft. Such edges will enable development of safer, more maneuverable aircraft. For these applications, the materials must maintain good properties and operate with no or limited oxidation or ablation in high temperature, reactive environments. Industrial applications for UHTCs are also of interest and include use in foundries, molten metal crucibles, and nuclear reactors. While initial studies of boride compounds such as UHTCs began in the 1960s [3], progress is still needed on a range of issues related to processing, properties, and performance of these materials. For example, thermal and chemical stability provide challenges for all these issues and also presents an obstacle for development of economical manufacturing. Better understanding of the effect of processing on microstructure and ultimately on properties are needed as well as a fuller elucidation of the basic mechanisms responsible for the response of these materials to extreme environments. Further innovations in the materials will require thorough knowledge of their properties and behavior under diverse conditions of temperature, environment and stress. ⇑ Corresponding author. E-mail address: [email protected] (M.S. Daw). 0927-0256/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2011.04.038
Computational modeling of these borides has been limited in part by the absence of appropriate interatomic potentials. While potentials exist for some carbides and nitrides, and there has been wide experience with their application, potentials for borides have been relatively less explored, especially metal borides. To our knowledge, until this work, there have been no interatomic potentials for ZrB2, HfB2 or any similar material. In the development of interatomic potentials, care must be exercised in the choice of semi-empirical form. A common form for carbides and nitrides is the Tersoff potential, by which one seeks to represent something of the angular dependence of whatever covalent bonding is present. Tersoff-style potentials have been developed for SiC, which is a common component of UHTCs. In this paper, we report on our development of the first interatomic potentials for ZrB2 and HfB2. We have chosen to develop Tersoff-style potentials for the diborides, based upon two considerations. First, the crystal structure (C32) exhibits closed-packed sheets of the transition metals, reminiscent of the HCP structure of the metals, interleaved with hexagonal B sheets. It seems quite possible that the bonding within the B sheets may show some degree of covalency. Second, UHTCs are frequently formed as composites of metal-diborides and SiC, and the latter already has well-established Tersoff-style potentials. The ultimate goal of this line of work would be to investigate some of the behavior of composites. For this reason, it would be desirable to combine older existing SiC potentials with the new potentials, and this is made most convenient if the new potentials already share the format of the older potentials. The task of the present work is to develop Tersoff-style interatomic potentials for ZrB2 and HfB2 which will be suitable for
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835
furthering the understanding of its mechanical and thermal properties. Development of potentials for compounds is a task which involves many intermediate steps and decisions, and it is appropriate to keep in mind that our aim is to study Zirconium and Hafnium Diborides in the context of UHTC applications. The Tersoff potentials of compounds contain a large number of adjustable parameters, and progress is easiest if we subdivide the large parameter space. We have therefore pursued a strategy of isolating subsets of the parameters in order to obtain some control of the fitting process. Our general strategy is to obtain first a Tersoff potential for elemental Zr and Hf, then a Tersoff potential for elemental B, and then determine the cross-type parameters for the compounds. Target properties for the elements include lattice constants, cohesive energies, elastic constants, and vacancy energies in a variety of crystal structures. For the compound we target the lattice constants and cohesive energy. The transition metals present no special challenge in this matter. There exist already several potentials for Zr — for example, those by Willaime and Massobrio [4,5] — and there are good reasons to have some level of confidence in these potentials. It is also convenient that the form of the potential by Willaime and Massobrio can be easily converted to Tersoff-form, as will be described below. We have also obtained a new Zr potential for this work by optimizing to the current database, as well as a new Hf potential. However, the story for the trivalent metalloid (B) is quite different [6–10]. There are several polymorphs established for elemental boron, and there are two phases which compete for stability at ambient conditions. Both of the two low-energy structures have large unit cells and complex structure. The subtle bonding variations seem to occur because of competition between metal and insulating behavior, with some covalent bonding likely in other cases. It seems unlikely that interatomic potentials will be able to capture such a broad range of behavior, which may also account for the dirth in the literature of reports on interatomic potentials for B. We have therefore focused at this stage on obtaining potentials which capture the gross behavior over a range of coordination, as will be described in more detail in the next section. It is especially important, then, to remember the context of this work: we are trying to develop potentials which are to be used in a relatively narrow range of circumstance (within the diboride) and only for a narrow range of properties (mechanical and thermal). The assertion is then that the behavior of B in this narrowly defined context might be represented with sufficient success by interatomic potentials. It is anticipated that further refinements of the potential may be necessary, with this first cut serving as a basis for evaluating possible weaknesses. Once the elemental Zr, Hf and B potentials are established, we use the properties of the compound ZrB2 (HfB2) to determine the cross-terms in the Tersoff potential. Because we have two, separate Zr potentials (one from the literature and one original to this work), we do this determination twice, so that we have two separate potentials for ZrB2. Only one potential for HfB2 is obtained here. We then test the two interatomic potential sets by comparing to the results of first principles calculations for a variety of properties that were not targets for the fitting procedure. The tests include elastic constants, point defect energies, phonon dispersion curves, and stability of two naturally occuring boron structures. Overall the results are reasonable, showing that the potentials are stable and generally reliable. Our conclusion is that these first potentials are ready to be used in more sophisticated calculations of mechanical and thermal properties. This paper is organized as follows. In the next section, we give the details of the determination of the interatomic potential, in terms of the specifics of the fitting process, including properties that were used as targets for the fitting. In the next section we
2829
describe the independent testing that we subjected the potentials to, as well as the results. Finally, we discuss the qualities of the potentials, as well as suggestions for the direction of future refinements. 2. Method As mentioned in the Introduction, our general strategy for determining a Tersoff potential for ZrB2 (HfB2) has been first to obtain functions for element Zr (Hf) and B, and then to obtain the cross-element parameters from properties of the compound. Our approach has been to perform first-principles calculations of specific properties to serve as targets for the fitting. In this section, we describe the generation of targets for fitting, then details of the fitting for elemental Zr and Hf, elemental B, and finally the compounds. 2.1. Generation of targets for fitting Target values for fitting were generated from ab initio calculations based on density functional theory (DFT). Results from a more comprehensive set of ab initio calculations for ZrB2 and HfB2 will appear in a separate publication [11] as will further computational details. We summarize briefly some of those details here. Computations were performed with the plane wave DFT codes VASP [12,13] and ABINIT [14,15]. While there is considerable overlap between the two codes, each also offers unique features. Using different codes also allows us to check our results. FHI norm-conserving pseudopotentials were used with ABINIT and the projected augmented wave (PAW) potentials with VASP. The functional of Perdew, Burke, and Ernzerhof (PBE) [16] was used throughout. For pure Zr (Hf), cohesive energies, lattice constants, elastic constants and vacancy energies were calculated for HCP, FCC and BCC crystal structures. For pure B, the only stable structures known have quite complex unit cells. Rather than deal with the difficulties of fitting to properties of structures with large unit cells, we have focused on the properties of simple structures (SC, FCC, hexagonal sheet, and triangular sheet). For those, we have calculated cohesive energies, lattice constants, and bulk moduli. Note that the full set of elastic constants for B was not calculated because for these simple structures, B is unstable to shear. For the diborides, the cohesive energies and lattice constants were obtained as targets for fitting. Lattice constants were obtained by energy minimization. Cohesive energies were derived by evaluating the total energy at the optimized geometry and subtracting the atomic energies. Elastic constants were obtained by the linear response method [17,18] as implemented in ABINIT and the finite difference method as implemented in VASP. We expect the linear response method to be slightly more accurate than finite differences since it does not rely on convergence with respect to supercell. All quantities were converged individually with respect to k-space sampling and plane wave cutoff. 2.2. Definition of the potentials Because some variations of the Tersoff potential for compounds have emerged in the literature [19–24], we make explicit in the following the forms we are using for this work. The original Tersoffstyle form for multi-component systems [23] was augmented by Nord et al. [25] to include the c factors in the calculation of the bond order. These factors alter the way in which a third atom enhances or weakens the bond between a pair. For example, the factor cBBZr alters the way that Zr affects the B–B bond, making it
2830
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835
different than the way a B atom would affect the bond. These c factors were found particularly helpful in fitting to the properties of ZrB2. The total energy for a configuration of atoms in a multi-component system is given by
E¼
X
Ei ¼
i
1X V ij 2 i–j
ð1Þ
where the pair term is
h i ½ij ½ij ½ij V ij ¼ fC ðdij Þ fR ðdij Þ þ bij fA ðdij Þ
ð2Þ
The repulsive piece is given by ½ij
fR ðdÞ ¼ Aij expðkij dÞ
ð3Þ
and the attractive piece by ½ij
fA ðdÞ ¼ Bij expðlij dÞ
ð4Þ
These functions are restricted to short-range by this rather abrupt cut-off function
8 d < Rij > < 1; ½ij fC ðdÞ ¼ 12 þ 12 cos½pðd Rij Þ=ðSij Rij Þ; Rij < d < Sij > : 0; Sij < d
ð5Þ
The bond-order function is n
n
bij ¼ ð1 þ bi i fiji Þ1=2ni
ð6Þ
which uses the parameter
fij ¼
X
fC ðr ik Þcijk g i ðhijk Þ expf½k3i ðdij dik Þmi g ½ij
ð7Þ
k–i;j
The angular dependence is supplied by
h i 2 2 g i ðhÞ ¼ 1 þ c2i =di c2i = di þ ðhi cos hÞ2
ð8Þ
Note that parameters such as Aij depend only on the types of atoms i and j (not on positions). The parameters that depend only on a single type are b, n, c, d, h, m, and k3. The parameters that depend on two types are A, B, k, l, R, and S. It is often taken that the cross-type parameters (e.g., AZrB and kZrB) are restricted to be means of the same-type parameters. For example, kij = (ki + kj)/2 and Aij = (AiAj)1/2. However, we do not make that choice here. Instead, we allow the cross-type parameters to be flexible in the fit. 2.3. Potentials for Zr, Hf Tersoff-style potentials for Zr [26] and Hf [27] have been developed recently as part of a study of their oxides. We have developed here an independent set of Tersoff potentials for Zr and Hf. It is usually argued that metallic bonding in highly coordinated environments does not display significant angular dependence. With that in mind, Embedded Atom Method/Finnis–Sinclair (EAM/FS) [28–30] forms were successfully obtained for Zr by several authors [31,4,5,32]. The FS potentials due to Willaime and Massobrio (WM) [4,5] are of such form that they can be easily converted to Tersoff form (see note by Brenner [33]). We therefore have converted the second WM potential (‘‘WM2’’) to Tersoff form as the first of our candidate Zr potentials which we will refer to as ‘‘Zr-1’’. We will present the parameters and the calculated properties below. We have also used the Zr WM2 potential as a starting point for a re-optimization to the current targets. The second Zr potential is somewhat altered from the WM2 potential and will be referred to here as ‘‘Zr-2’’.
Finally, we used the reoptimized Zr potential as the starting point for the new Hf potential. The target properties we chose for Zr and Hf included a full range of properties related to the close-packed structures FCC and HCP as well as the more open BCC structure which is expected to be similar in energy. To simplify the fitting process, we focused on the FCC form of Zr, which should be quite similar in its properties to those of HCP. We therefore chose to fit to these properties in FCC-Zr: lattice constant, cohesive energy, elastic constants, vacancy energy, as well as the pressure derivative of the bulk modulus. In addition we fitted to the cohesive energies of HCP and BCC phases, so ensure a good representation of their relative stabilities. As has been noted in many places (see, for example, Willaime and Massobrio), the FCC-BCC energy difference is sensitive to cutoff. We also have chosen to include 2nn in the FCC structure in the current potentials. The new set of Zr and Hf parameters were obtained by doing a least squares fit to the target data. We chose initial values for the Zr parameters to be those of the WM2 potential. For Hf, we used the newly obtained Zr potential as the starting point. Then the fitting proceeded in stages. The first adjustments were to A, B, k, l, and k3. The second stage allowed b and n to flex. In the third stage, all parameters were allowed to vary. We attempted to increase the value of c, which would have introduced some bond-angle dependence to the functions, but the fit to the target data was always better with this parameter at zero. In Tables 1 and 2 we show the fidelity of the Zr and Hf potentials in representing the chosen target properties. All potentials are reasonably good. They are all relatively weak in representing the derivative of the bulk modulus, which might be expected from a 2nn potential. This weakness should not be fatal in the present application. All three of the parameter sets (two for Zr and one for Hf) were subjected to simple stability tests, such as relaxing around a vacancy, adding significant random distortions to the FCC structure and relaxing, and running molecular dynamics at elevated temperature. All sets gave satisfactory results. The resulting parameters for the potentials are displayed in Table 3. We note that we also tried several other Zr forms, with shorter cut-offs and also allowing for angular dependence, but these EAM/FS-like potentials produced the most satisfactory results.
2.4. Potential for B One consideration in developing interatomic potentials is to ensure stability, at least within the range of environments that the potential is likely to see service. Boron is challenging for the usual fitting procedure, because all of the simple structures (hexagonal
Table 1 The target properties and the Tersoff values for Zr calculated using the parameters WM2 [4,5] (‘‘Zr-1’’) and the parameters derived in this work (‘‘Zr-2’’). B is the bulk modulus and B0 here is the derivative of the bulk modulus w.r.t. lattice constant at equilibrium. Vf is the vacancy formation energy. Property (units)
Target
‘‘Zr-1’’ (WM2)
‘‘Zr-2’’ (new)
a0(FCC) (Å) E0(FCC) (eV) B(FCC) (eV/Å3) B0 (FCC) (eV/Å4) C11(FCC) (eV/Å3) C12(FCC) (eV/Å3) C44(FCC) (eV/Å3) Vf(FCC) (eV) a0(HCP) (Å) E0(HCP) (eV) E0(BCC) (eV)
4.530 6.160 0.578 0.86 0.774 0.481 0.356 2.500 3.230 6.180 6.050
4.532 6.127 0.573 1.95 0.712 0.503 0.347 2.211 3.256 6.142 6.124
4.547 6.174 0.573 1.93 0.664 0.527 0.347 2.224 3.215 6.174 6.006
2831
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835 Table 2 The target properties and the Tersoff values for Hf calculated using the new parameters derived in this work. (See the caption for Table 1 for a description of the properties.) Property (units)
Target
New pot
a0(FCC) (Å) E0(FCC) (eV) B(FCC) (eV/Å3) B0 (FCC) (eV/Å4) C11(FCC) (eV/Å3) C12(FCC) (eV/Å3) C44(FCC) (eV/Å3) Vf(FCC) (eV) a0(HCP) (Å) E0(HCP) (eV) E0(BCC) (eV)
4.470 6.410 0.649 0.92 0.911 0.518 0.462 2.580 3.190 6.480 6.300
4.570 6.469 0.606 1.92 0.766 0.536 0.347 2.296 3.231 6.469 6.296
Table 3 Parameters describing the two different potentials for Zr and one potential for Hf. The ‘‘Zr-1’’ is from the work of Willaime and Massobrio. ‘‘Zr-2’’ and ‘‘Hf’’ are derived for this work. See also Brenner’s note on the relationship of various interatomic potential forms [33]. Parameter
‘‘Zr-1’’ (WM2)
‘‘Zr-2’’ (new)
‘‘Hf’’ (new)
A (eV) B (eV) k (Å1) l (Å1)
3.65374 103 3.68286 101 2.92969 0.661542
3.58787 103 3.71035 101 2.91078 0.659036
3.15839 103 3.96250 101 2.83165 0.65895
b n c d h
1 1 0 1 0
1 1 0 1 0
1 1 0 1 0
R (Å) S (Å)
4.88 5.18
4.6 5.0
4.6 5.0
m k3 (Å1)
1 1.32308
1 1.32308
1 1.32308
and triangular sheets, as well as SC, BCC, and FCC) are unstable to shear. While it is certainly possible to fit to unstable structures, one must surrender the assurance given by the usual stability tests like we did for Zr in the previous section. In our present approach, we attempted to separate the properties and parameters according to whether they relate to any possible angular dependence of B bonding. This can be achieved in the way suggested by the recent work of Liang, Phillpot, and Sinnott [34]. They noted that the 1nn Tersoff form can effectively factor into two parts, with the separation being afforded by the combined factor of the repulsive coefficient B and the bond order parameter b. Taking this product, bB, and noting that b depends only on the bond angles present in a given crystal, it is then possible to separate the bond-angle and bond-length dependence. This is carried out by fitting to the lattice constant, cohesive energy, and bond stiffness for four simple structures: hexagonal sheet, triangular sheet, SC, and FCC. This is especially simple if we assume 1nn form, as we have here. The fitting determines the parameters A, k, l, and the value of bB in each of the four structures. The result of this fit is shown in Table 4, and the same results are plotted Figs. 1–3, where one can see that the quality of the fit is reasonable. In Fig. 1, one can see that the cohesive energies are reproduced to within 0.1 eV and the bondlengths to within 0.1 Å. In Figs. 2 and 3, we can see that the trends of bond strength and bond stiffness are correctly captured. To complete the fit, the resulting bond order parameters must be related to the Tersoff form in Eq. (6), which incorporates some angular dependence as well as coordination dependence. However, the resulting bond order parameters show a simple dependence on
Table 4 The target properties and the fitted values for boron, used to establish the (bB) factors for the listed structures. Lattice constants are in Å, cohesive energies in eV, and bond stiffnesses in eV/Å2. [Note that these are calculated using 1nn only in each structure. Using the Tersoff potential with cut-off fixed by distance — as would ordinarily be done — may yield somewhat different answers.]. Structure
Property
Target
Fit
Hexagonal sheet
a0 E0 E000
2.91 5.15 11.35
2.89 5.08 7.98
Triangular sheet
a0 E0 E000
1.70 5.71 21.73
1.81 5.75 27.06
SC
a0 E0 E000
1.88 5.33 24.50
1.84 5.21 24.51
FCC
a0 E0 E000
2.86 5.07 21.85
2.84 5.22 12.28
Fig. 1. Comparison of fit to target for cohesive energy vs. bondlength for B in four structures. The target data are represented by circles and the fits by squares. The four structures are (in order of increasing bondlength): hexagonal sheet, triangular sheet, SC, and FCC.
1
Fig. 2. Comparison of fit to target for bondstrength vs. bondlength for B in four structures. The target data are represented by circles and fits by squares. The four structures are (in order of increasing bondlength): hexagonal sheet, triangular sheet, SC, and FCC.
coordination and no evidence of angular bonding, as shown when we plot the bond order parameters vs. coordination in Fig. 4. We therefore set c = 0 to eliminate the angular dependence and fitted the remaining coordination dependence using b and n. The quality of the fit is illustrated in the figure. The resulting parameters for B are shown in Table 5, which are very simple. As noted above, it is not possible to subject these B parameters to the same sort of stability tests as Zr.
2832
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835 Table 6 The target properties and the Tersoff values for ZrB2 (C32 structure) calculated using the two alloy sets obtained here. Property (units)
Target
Pot 1
Pot 2
a0(C32) (Å) c0(C32) (Å) E0(C32) (eV)
3.170 3.550 21.70
3.140 3.547 21.54
3.131 3.484 21.35
The procedure was repeated for HfB2. The targets and the fits are displayed in Table 7. The resulting cross-type Tersoff parameters for the alloys are presented in Table 8. Fig. 3. Comparison of fit to target for bond stiffness vs. bondlength for B in four structures. The target data are represented by circles and fits by squares. The four structures are (in order of increasing bondlength): hexagonal sheet, triangular sheet, SC, and FCC.
3. Tests for the compound potentials
Parameter
Value
A (eV) B (eV) k (Å1) l (Å1)
1.0619848 103 3.6767392 101 3.955278 1.190167
In order to test independently of the target properties, we have evaluated — in ZrB2 and HfB2 — point defects energies, elastic constants and phonon spectra, as well as energy vs. lattice constant for large strains. All three sets of interatomic potentials were tested. First, we compared certain point defect energies as predicted by the interatomic potentials to the results of first-principles. In computing energies of defects with varying stoichiometry, it is necessary to be somewhat cautious. The interatomic potentials are by nature constructed with respect to isolated atoms as the reference of energy. The first-principles calculations can be defined with respect to the same reference, but that involves a calculation of the atomic energies for which there is some uncertainty. The most secure way around this ambiguity is to compare sets of point defects that preserve overall stoichiometry. We have thus computed the energy of a pair of anti-sites (interchanging distant atoms of different type) as well as a trio of vacancies (removing one M and two B from three distant sites), and comparing this to the cohesive energy per unit cell of the compound. The results for the point defects are presented in Tables 9 and 10. In calculating the concentrations of various point defects in the stoichiometric compound, certain combinations of the point defect energies occur. Those combinations which have thermodynamic significance for the stoichiometric compounds we present in Tables 11 and 12. The defect energies are themselves quite large, and the differences between interatomic potentials and first-principles are also significant. However, it is not expected that the energies of point defects such as these are likely to be important in the properties for which these potentials are intended. With that in mind, it seems that the point defect energies are acceptable. The next test was to compare elastic constants and phonon spectra obtained from the potentials versus ab initio computations. As mentioned previously, the ab initio results will be discussed in more detail in a forthcoming publication [11]. In brief, for the ab initio calculations, both the elastic constants and the phonon spectra were calculated using linear response theory as implemented in ABINIT. As with the target calculations, FHI pseudopotentials and the PBE functional were used. All quantities were individually converged with respect to k-space sampling and plane wave cutoff. Elastic constants were computed both with and without internal relaxation; however, there is only a very small amount of relaxation.
b n c d h
0.541 21.885575 0 1 0
Table 7 The target properties and the Tersoff values for HfB2 (C32 structure) calculated using the alloy set obtained here.
R (Å) S (Å)
2.2 2.5
m k3 (Å1)
1 0.
Fig. 4. Bond order parameter vs. coordination for B. The points are the established values (see text) for four structures, in order of increasing coordination: hexagonal sheet, triangular sheet, SC, and FCC. The curve is our fit using the Tersoff form.
2.5. Potentials for Zr–B and Hf–B Given the two candidate Zr potentials and the new B potential described in the preceding sections, we then moved onto the determination of the alloy parameters for ZrB2. Our approach was to keep fixed those parameters determined to this point and to determine the cross-type parameters from the properties of ZrB2. Three simple properties were chosen for the fitting: the two lattice constants and the cohesive energy of ZrB2 in the C32 structure (see Table 6). These properties were fit easily by varying primarily the c factors, though the other cross-type parameters were also adjusted.
Table 5 The parameters of the Tersoff form from the new fit to boron described in the text.
Property (units)
Target
Current
a0(C32) (Å) c0(C32) (Å) E0(C32) (eV)
3.15 3.52 22.47
3.13 3.42 22.39
2833
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835 Table 8 The cross-type parameters of two Tersoff forms for the Zr–B and Hf–B alloys. Two potentials are shown. Zr–B ‘‘pot 1’’ combines with the ‘‘Zr-1’’ from Table 3, Zr–B ‘‘pot 2’’ with the ‘‘Zr-2’’ parameters, and Hf–B with the Hf parameters. All three diboride potentials use the same B potential (Table 5). Parameter
Zr–B (pot 1)
Zr–B (pot 2)
Hf–B
AMB (eV) BMB (eV) kMB (Å1) lMB (Å1)
1.9008612 103 4.0127474 101 3.3688840 0.9258545
1.7659542 103 4.0023041 101 3.360803944 0.9220569
1.6967221 103 4.2730462 101 3.330828218 0.92455101
RMB (Å) SMB (Å)
3.28 3.60
3.18 3.54
3.18 3.54
cMMM cMMB cMBM cMBB cBMM cBMB cBBM cBBB
1.0 0.1905300 9.5016209 0.9751164 0.9999995 1.0000000 0.1721612 1.0
1.0 0.2930705 9.0588265 3.5295546 0.94069107 0.82400007 0.17523123 1.0
1.0 0.25267505 8.0112269 3.7926378 0.9934950 0.81400090 0.1620252 1.0
Table 13 Calculated elastic constants for ZrB2, comparing first-principles results and the predictions of the two interatomic potentials (‘‘pot 1’’ and ‘‘pot 2’’) developed here. Elastic constants are in GPa. Method
c11
c12
c13
c33
c44
B
G
A(=c33/c11)
ABINIT Zr–B (pot 1) Zr–B (pot 2)
543 422 575
56 156 213
113 171 301
419 320 817
234 119 246
230 240 400
223 118 211
0.77 0.76 1.42
Table 14 Calculated elastic constants for HfB2, comparing first-principles results and the predictions of the interatomic potentials developed here. Elastic constants are in GPa. Method
c11
c12
c13
c33
c44
B
G
A(=c33/c11)
ABINIT Hf–B
578 603
71 224
118 303
430 866
238 259
244 415
231 224
0.74 1.43
K
M
900
Method
E(ZrB)
E(BZr)
E(VZr)
E(VB)
ABINIT pot 1 pot 2
11.2 7.2 8.5
8.6 5.3 5.2
11.2 9.6 9.3
9.3 8.2 8.1
800
Frequency (cm-1)
Table 9 Some (relaxed) point defect energies for ZrB2, comparing two potentials and first principles. Energies are in eV.
E(ZrB)
E(BZr)
E(VZr)
E(VB)
ABINIT new
11.6 8.6
8.4 5.9
10.6 10.2
9.0 8.3
E(VZr) + 2E(VB) + Ecoh
ABINIT pot 1 pot 2
19.80 12.54 13.72
7.80 4.57 4.15
Table 12 Stoichiometric combinations of point defect energies for HfB2, comparing the new potentials and first principles. Energies are in eV. Method
E(HfB) + E(BHf)
E(VHf) + 2E(VB) + Ecoh
ABINIT new
20.01 14.45
6.53 4.46
400 300
0
A
L
H
A
Fig. 5. Normal mode frequencies for ZrB2 calculated using interatomic potentials from ‘‘pot 1’’ (solid lines), compared to those for ABINIT (dashed lines).
900 800 700
Frequency (cm-1)
E(ZrB) + E(BZr)
500
100
Table 11 Stoichiometric combinations of point defect energies for ZrB2, comparing two potentials and first principles. Energies are in eV. Method
600
200
Table 10 Some (relaxed) point defect energies for HfB2, comparing the potential and first principles. Energies are in eV. Method
700
600 500 400 300 200 100 0
K
M
A
L
H
A
Fig. 6. Normal mode frequencies for ZrB2 calculated using interatomic potentials from ‘‘pot 2’’ (solid lines), compared to those for ABINIT (dashed lines).
To calculate the elastic constants for the potentials, we generated a set of five symmetry independent deformations of the crystal parameterized by a strain parameter a [35]. Curves were generated of energy versus a. Elastic constants were obtained by fitting Taylor’s expansions of the energy with respect to the strain parameter. To calculate the phonon spectra for the potentials, the ‘‘frozen phonon’’ method implemented in the code PHON [36] was used. In this method, symmetry independent atomic displacements were
generated and the corresponding forces were calculated using the potentials in conjunction with the LAMMPS [37] code. The dynamical matrices were then constructed, diagonalized and the phonon spectra obtained. Supercells up to 8 8 8 cells gave well-converged values. The calculated elastic constants are presented in Tables 13 and 14. Generally, the results are acceptable. Specifically, the worst agreement is with some of the shear constants. It is reasonable
2834
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835
900
harmonic part is very good for significant strains (up to 20%). For larger strains, the interatomic potentials become somewhat too stiff. We note that the overly stiff behavior under tension (a > a0) is unavoidable with interatomic potentials restricted to first or second nearest-neighbors. This would seem to be related to the over-estimate of the magnitude of B0 for the elemental metals (see Tables 1 and 2.) Similar results are shown for HfB2 in Fig. 9. Overall the results show that these first interatomic potentials for ZrB2 and HfB2 provide a reasonable representation of the properties of the compounds.
Frequency (cm-1)
800 700 600 500 400 300 200
4. Conclusions
100 0
K
M
A
L
H
A
Fig. 7. Normal mode frequencies for HfB2 calculated using interatomic potentials derived in this work, compared to those for ABINIT (dashed lines).
Fig. 8. Energy vs. lattice constant for ZrB2 for large strain, comparing the two interatomic potentials and first-principles results.
We have presented here our derivation of the first interatomic potentials for ZrB2 and HfB2, which are obtained in the Tersoff format. We have determined two variant sets of parameters by fitting to targets generated by first-principles for a variety of properties, including lattice constants, cohesive energies, vacancy energies, elastic constants, and others, for Zr, Hf, B, and the compounds ZrB2 and HfB2. The agreement with the targets is reasonable. We have also determined that the functions appear to provide a reasonable representation of additional properties of the compounds outside of the target sets, including elastic constants, phonon dispersion, and point defects. The potentials have so far passed all tests of stability. If there is an obvious weakness of these first potentials for ZrB2 and HfB2, it would seem to be the complete lack of angular bonding among the B atoms. This result was determined by our database of structures, which included but did not emphasize the hexagonal sheet. However, it is not obvious at this point how the lack of angular dependence should be expressed in the various predicted properties of the compound. It is therefore judged that these potentials represent a reasonable first version. From here, further testing will be carried out which may suggest directions for future refinements. Acknowledgement MSD was supported under a NASA prime contract to ELORET Corporation. References
Fig. 9. Energy vs. lattice constant for HfB2 for large strain, comparing the two interatomic potentials and first-principles results.
to expect that including angular bonding in the boron potentials would significantly improve the shear moduli. However, there is some differences in the quality even between our two Zr–B potentials, neither of which contains angular dependence for the boron. This would seem to be the best place to begin improvements. The calculated phonon spectra for ZrB2 are presented for in Figs. 5 and 6 (for ‘‘Zr-1’’ and ‘‘Zr-2’’, respectively) and for HfB2 in Fig. 7. Generally the results are reasonable. Finally, we test the anharmonic part of the potentials by applying a large hydrostatic strain. For ZrB2, Fig. 8 shows that the
[1] M.J. Gasch, D.T. Ellerby, S.M. Johnson, Ultra high temperature ceramic composites, in: N.P. Bansal (Ed.), Handbook of Ceramic Composites, Kluwer Academic, 2005, pp. 197–224. [2] R. Loehman, E. Corral, H.P. Dumm, P. Kotula, R. Tandon, Ultra high temperature ceramics for hypersonic vehicle applications, Technical Report SAND20062925, Sandia National Labs, Albuquerque, NM, 2006. [3] L. Kaufman, E.V. Clougherty, Investigation of boride compounds for very high temperature applications, Technical Report RTD-TRD-N63-4096 Part III, ManLabs, Inc., Cambridge, MA, 1966. [4] F. Willaime, C. Massobrio, Phys. Rev. Lett. 63 (1989) 2244. [5] F. Willaime, C. Massobrio, Phys. Rev. B 43 (14) (1991) 11653. [6] M.W. Chase Jr., NIST-JANAF Thermochemical Tables, 1998 (issued as J. Phys. Chem. Ref. Data monograph no. 9). [7] B.E. Douglas, S.-M. Ho, Structure and Chemistry of Crystalline Solids, Springer, 2006. [8] M.J. van Setten, M.A. Uijttewaal, G.A. de Wijs, R.A. de Groot, J. Am. Chem. Soc. 129 (9) (2007) 2458. [9] H. Tang, S. Ismail-Beigi, Phys. Rev. Lett. 99 (2007) 115501. [10] A.R. Organov, J. Chen, C. Gatti, Y. Ma, Y. Ma, C.W. Glass, Z. Liu, T. Yu, O.O. Kurakevych, V.L. Solozhenko, Nature 457 (2009) 863. and references therein. [11] J.W. Lawson, C.W. Bauschlicher, M.S. Daw, J. Am. Ceram. Soc. (2011), in press. [12] G. Kresse, J. Furthmller, Phys. Rev. B 54 (1996) 11169. [13] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758. [14] M. Fuchs, M. Scheffler, Comput. Mater. Sci. 119 (1999) 67. [15] X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, Ph. Ghosez, J.-Y. Raty, D.C. Allan, Comput. Mater. Sci. 25 (2002) 478. [16] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [17] X. Gonze, Phys. Rev. B 55 (1997) 10337. [18] X. Gonze, C. Lee, Phys. Rev. B 55 (1997) 10355.
M.S. Daw et al. / Computational Materials Science 50 (2011) 2828–2835 [19] [20] [21] [22] [23] [24] [25] [26] [27]
J. Tersoff, Phys. Rev. Lett. 56 (1986) 632. J. Tersoff, Phys. Rev. B 37 (1988) 6991. J. Tersoff, Phys. Rev. B 38 (1988) 9902. J. Tersoff, Phys. Rev. Lett. 61 (1988) 2879. J. Tersoff, Phys. Rev. B 39 (1989) 5566. J. Tersoff, Phys. Rev. B 49 (1994) 16349. J. Nord, K. Albe, P. Erhart, K. Nordlund, J. Phys.: Cond. Matter 15 (2003) 5649. T. Iwasaki, J. Mater. Res. 20 (2005) 1300. T.-R. Shan, B.D. Devin, T.W. Kemper, S.B. Sinnott, S.R. Phillpot, Phys. Rev. B 81 (2010) 125328. [28] M.S. Daw, M.I. Baskes, Phys. Rev. Lett. 50 (1983) 1285.
[29] [30] [31] [32] [33] [34] [35] [36] [37]
2835
M.S. Daw, M.I. Baskes, Phys. Rev. B 29 (1984) 6443. M.W. Finnis, J.E. Sinclair, Philos. Mag. 50 (1984) 45. D.J. Oh, R.A. Johnson, J. Mater. Res. 3 (33) (1988) 471. R.C. Pasianot, A.M. Monti, J. Nucl. Mater. 264 (1–2) (1999) 198–205. D.W. Brenner, Phys. Rev. Lett. 63 (1989) 1022. T. Liang, S.R. Phillpot, S.B. Sinnott, Phys. Rev. B 79 (2009) 245110. L. Fast, J.M. Wills, B. Johansson, O. Eriksson, Phys. Rev. B 15 (1995) 17431. D. Alfe, Comput. Phys. Commun. 180 (2009) 2522. S.J. Plimpton, J. Comput. Phys. 117 (1995) 1–19. further info at lammps.sandia.gov.