Modified Embedded Atom Method Potential for Modeling the Thermodynamic Properties of High Thermal Conductivity Beryllium Oxide.

Modified embedded atom method potential parameters of beryllium oxide (BeO) have been developed, which can well reproduce the thermodynamic properties of beryllium oxide. To accurately describe the interactions between the atoms in the BeO structure, the density functional theory is used to calculate the fundamental properties such as the lattice constant, bulk modulus, and elastic constant, which are used for the potential fitting. The properties such as the enthalpy and specific heat are used to test the validity of the potential parameters. The calculated results by the developed potential parameters are compared with the experimental and other theoretical data as a function of temperature. The good agreement between the calculated results by the new potential and the experimental data verifies the potential parameters. The developed potential parameters have also been used to predict the thermal conductivity of BeO as a function of temperature for further applications of beryllium oxide.

Introduction

Recently, beryllium oxide (BeO) has attracted much attention due to its outstanding properties, such as high hardness,[1] high electrical resistivity,[2] high thermal conductivity,[3] high melting point, and excellent radiation resistance.[4] Beryllium oxide is widely used in many high-performance semiconductor parts due to its good thermal conductivity and good electrical resistivity. It is also used as a structural ceramic for high-performance microwave devices, vacuum tubes, magnetrons, and gas lasers. The band gap of BeO is about 10.63 eV,[5,6] which makes it suitable for application in optoelectronic devices.[7] BeO also has a long history of being used in the nuclear area for both fission and fusion power generation due to its excellent mechanical, thermal, and chemical properties.[1,8] The excellent chemical compatibility between BeO and UO2 makes it possible to prepare UO2 and BeO composite fuel. By coupling with BeO, the poor thermal conductivity of UO2 could be effectively improved.[9,10] Moreover, BeO possesses good irradiation resistance and neutron scattering performance, which increase the preference of its application in nuclear reactors.[11,12]

In the past several decades, the electronic, optical, and mechanical properties of BeO have been reported in many experimental and theoretical works.[5,6,13−19] Those properties are attributed to the unique structure of BeO, which is the only alkaline-earth oxide crystallizing in the hexagonal wurtzite structure[20] at room temperature, while the other oxides in this family have the rock salt structure. Compared with the rock salt structures, the wurtzite structure has a stronger covalent bond, and the Be–O bond is shorter than that of the rock salt structures. This particular structure makes BeO have some unique physical properties. For the BeO system, it has a hexagonal phase of wurtzite (WZ), cubic phase of zinc blende (ZB), and cubic phase with the rock salt structure. Among these phases, WZ BeO is the ground state, and the other two phases of BeO are hard to obtain. The experimental samples keep the wurtzite structure at least up to 55GPa without any transition observed.[21] Using a first-principles pseudopotential method, Park et al.[22] theoretically studied the series of phase transitions, wurtzite → zinc blende → rock salt; the two corresponding transition pressures for the wurtzite-to-zinc blende transition and the zinc blende-to-rock salt transition are 91 and 147 GPa, respectively. Boettger et al.[23] found the same trend by using two all-electron, full-potential electronic structure techniques. In the following parts, we will mainly focus on the WZ phase since it is the most common stable sample in the experiment.

The studies of BeO from an atomic aspect are very important in understanding the mechanism of the properties mentioned above. However, it is comparatively difficult for the experimental research compared with theoretical studies. Although first-principles calculation is a very effective method to study BeO from an atomic aspect, it is challenging to study the structure with a large number of atoms especially for the system with defects, such as vacancy and doping. Besides the first-principles calculations, another method to study the materials from an atomic aspect is the molecule dynamics (MD) simulations. Compared with the DFT calculations, MD has some advantages when studying the properties of the material, especially for the microcosmic system and the system with defects, which are very difficult for the first-principles calculations. Besides, the MD simulation is a very useful tool to study the thermal conductivity and thermodynamic properties compared with the first-principles calculations. Although Malakkal et al.[24] studied the thermal conductivity using first-principles calculations combined with Boltzmann transport equation (BTE) and the Slack model, the studied system is a perfect crystal structure. To study the system with the defect by MD, the potential of BeO is needed.

In the MD simulations, the interaction between atoms is described by the potential parameters. It should be mentioned that the potential parameters play a crucial role in accurately calculating the properties of materials. Recently, Byggmästar et al.[25] have developed a new empirical bond order potential to study the nanotube structure and nanosheet BeO. However, their potential aims to describe the interaction of atoms in the nanostructure and underestimated the thermal conductivity for the bulk structure. In the following parts, we will show that our present potential parameters exhibit more improvements on those properties, especially for thermal expansion and heat transfer in the bulk structures.

In this work, we theoretically studied the thermodynamic properties of BeO by using the MD simulations. First, we developed the potential parameter of the BeO system with the potential form of modified embedded atom method (MEAM). To obtain the accurate atomic interaction, the parameters for this potential are carefully tested and compared with the results calculated using the DFT theory and experimental data. Using the new potential, the expansion is examined over a wide temperature range and compared with the experimental data. This property was not included in the fitting procedure, and as such it provides one test of the approach. Similarly, the temperature dependence of the expansion, enthalpy, and heat capacity is determined and compared with the experimental data. Finally, the thermal conductivity is studied using the new fitted potential parameters.

Results and Discussion

Potential Parameters

For the MEAM potential, there are 14 parameters for each kind of element, as shown in Table . The Ec is the cohesive energy of the reference structures, which is hexagonal close-packed (hcp) for Be and dimer for O in this work. re, α, and A are the distance between two atoms in the equilibrium structure, the parameter related to the bulk modulus and cohesive energy, and adjusted parameters, respectively. The four parameters β(0), β(1), β(2), and β(3), which are related to the electron density, are used to adjust to the decay rate in different directions. The other four parameters t(0), t(1), t(2), and t(3) are the angular contributions to the electron density. In the hexagonal structure, the values of c/a (the ratio of lattice constant) and elastic constant are significantly affected by t(3), and the value of t(3)should be less than zero to fulfill the fact that c/a < 1.[26] For the density scaling parameters, which are element-dependent, the values of 1.0 and 3.77 are chosen for Be and O, respectively. Moreover, the screening parameters play an essential role in the MEAM potential. In our fitting progress, we find that the parameters of the Be and O pair screened by the O atoms have a significant effect on the elastic constant. For the static calculations, the potential gives a better agreement of the elastic constant with the experimental results if these parameters are not included. Here, these parameters are not used in the static calculations at 0 K. However, the parameters are critical when the temperature increases. So, this parameter is carefully tested and chosen as 0.5 and 2.80 in the MD simulations. For Be, the parameters partially come from the data of ref (27), while there are also some differences due to the Be atoms located at a different environment in the Be or BeO structure. Similar to the Be atoms, the parameters of the O atom partially come from ref (28). First, the parameters of Be and O use the same values as in refs (27) and (28). Next, we start the fitting of potential parameter of the BeO binary interatomic potential. Since the hcp structure, rather than the fcc (face-centered cubic) or bcc (body-centered cubic) structure, shows better similarities to the structure of WZ BeO, the hcp structure is chosen as the reference structure in our fitting process. The Be–O dimer distance is first calculated by the DFT calculations. Based on this optimized distance, the value of re is carefully adjusted in the MD simulations under the temperature conditions. Simultaneously, the values of cohesive energy and α are also chosen based on the first-principles calculations. Finally, these parameters are slightly adjusted to make the value of elastic constant closer to the experimental data. The final parameters’ set derived by this fitting process is summarized in Tables and 2, which are for single elements and pair interactions, respectively. The structure parameters for WZ BeO are listed in Table , including lattice constant a0 (Å), bulk moduli B0 (GPa), and cohesive energy Ec (Ry/mol). The experimental results are also included in Table for comparison. The elastic constant for BeO is listed in Table . To evaluate our potential parameters, some reported data for the theoretical and experimental elastic constants are also included. It can be seen that our results calculated using the fitted parameters are close to the other theoretical data and experimental results.

Table 1

MEAM Potential Parameters for Be and Oa

element	Ec	re	α	A	β	β	β	β	t(0)	t(1)	t(2)	t(3)	Cmin	Cmax
Beb	3.43	2.24	3.85	1.00	0.54826	0.30416	0.33524	0.00001	1.00	4.556	21.956	–10.234	0.30	2.80
Oc	2.558	1.21	5.60	1.36	2.24826	1.00416	1.33524	1.0001	1.00	4.01	0.01	–0.0001	0.60	2.80

The units for Ec and re are eV and Å, respectively.

Ref (27).

Ref (28).

Table 2

Pair Interaction Parameters for Element Pairs (Be–O)a

X	Y	Ec (eV)	re (Å)	α	d	Cmin OBeBe	Cmax OBeBe	Cmin BeOO	Cmax BeOO	Cmin BeBeO	Cmax BeBeO	Cmin OBeO	Cmax OBeO
Be	O	4.19	1.842	4.80	0.005	0.40	0.82	0.50	2.26	0.65	1.80	0.50	2.80

Ec is the cohesive energy of the hcp structure, re is the equilibrium nearest neighbor distance, α is the exponential decay factor for the universal energy, d is the justified parameters, and Cmax and Cmin are the screening parameters (CABC denotes the screening parameter of the type A element to the type B and type C elements). The units for Ec and re are eV and Å, respectively.

Table 3

Structure Parameters for WZ BeOa

calculation	a0	c/a	u	B0	Ec
our work	2.702	1.57	0.377	177	0.21
GGA[22]	2.703	1.620	0.377	203	0.952
LDA[22]	2.650	1.624	0.378	224	1.058
other calc[23]	2.668	1.633		239	1.034
other calc[29]	2.639	1.629	0.3769	228
other calc[21]	2.775	1.58		186
other cal[25]	2.758	1.633	0.375	208	0.817
expt[1]	2.698	1.622	0.378	212
expt[6,30]					0.9

a0, lattice constant in Å; u, internal parameter; B0, bulk moduli in GPa; Ec, cohesive energy in Ry/mol. The experimental results are also listed.

Table 4

Elastic Constants for BeOa

elastic constant	C11	C12	C13	C33	C44	C66
our results	340	110	43	543	125	114
WZ[31]	450	102	75	471	144	174
GGA[16]	439	105	72	463	142	156
PIB[21]	366	113	90	361	132	126
ABOP[25]	463	98	62	499	164	183
experimental[14]	461	126	88	491	147	167
experimental[15]	470	168	119	494	153	152

To evaluate our potential parameters, some reported data for the theoretical and experimental elastic constants are also included. It can be seen that our results calculated using the fitted parameters are close to the other theoretical data and experimental results. The unit for elastic parameters Cxx is GPa.

Potential Fitting

The potential fitting started from the lattice constant, which is the fundamental property of materials. The structure of BeO, as shown in Figure , is fully relaxed, and the structure parameters such as the equilibrium lattice constant (a), cohesive energy (Ec), and bulk modulus (B) are listed in Table . The previous experimental results and other theoretical data are also given for comparison. In our calculations, the lattice constant is first optimized by the first-principles calculations. The calculated results show that the lattice constants optimized by DFT are a = b = 2.71 Å and c = 4.41 Å, which are slightly larger than the experimental values of a = b = 2.70 Å, and c = 4.38 Å. It is a common characteristic that the generalized gradient approximation Perdew–Burke–Ernzerhof (PBE) function[32] usually overestimated the lattice constant compared to the experimental data.

Figure 1

Crystal structure of wurtzite BeO from (a) the side view and (b) the corresponding unit cell; the blue and yellow spheres represent Be and O atoms, respectively.

With the optimized lattice constant, the energy surface of BeO is then calculated by the DFT calculations, which is used for the fitting lattice constant by the GULP (general utility lattice program) code.[33]Figure shows the calculated energy as a function of volume, and both the DFT and fitted results are shown for comparison. It can be seen that the energy surface exhibits a parabolic curve as a function of volume, and the lowest energy represents the most stable structure corresponding to the optimized lattice constant. Compared to the DFT results, the MEAM data exhibit a left shift toward the small volume direction due to the slightly smaller lattice constant. This is reasonable because the DFT lattice constants are slightly larger than the experimental results, that is, 0.37% larger for a and b and 0.685% larger for c. Importantly, the fitted MEAM data of the energy surface exhibit a similar curve to the DFT one, which indicates that the fitted bulk moduli will have similar values compared with the DFT results. The bulk modulus can be obtained by the fitting of energy as a function of volume based on the EOS (empirical equations of state) equations.[34]

Figure 2

Energy of BeO as a function of volume calculated by the first-principles calculation and MEAM potential.

Besides the lattice constant, the bond length and angles between atoms are also important for the characterization of BeO structure. Compared with the Be–Be bond and O–O bond, the Be–O bond has the smallest bond length, and the charge density analysis indicates that the bond of Be–O is stronger than the other two kinds of bonds. In other words, the interaction between Be and O atoms is very important for the BeO structure. Thus, the Be–O bond will be the main focus in this work. It should be mentioned that the BeO structure is not a perfect hexagonal structure, and the value of c/a is 1.624 (calculated by first-principles calculation), which is not equal to the ideal value of 1.633. For the Be–O bond, the first-principles calculations show that there are two types of bond for Be–O in the unit cell. The bond length along the z direction (c direction) is 1.664 Å, while the length of the bond in the x–y plane (a–b plane) is 1.656 Å, as shown in Figure . For the MEAM potential, our fitting results show that the bond length along the z direction is 1.602 Å, while the length of the bond in the x–y plane is 1.645 Å.

To better compare the fitted results with the DFT results, the pair correlation function is also calculated as shown in Figure . Note that the MD result is in a dynamic state at room temperature, while the DFT data is in a static state at 0 K. Compared with the DFT results, the values of the three peaks for the MD results are different, and the location site as a function of distance shows a shift to the direction of the small value. This is reasonable since the atoms vibrate much harder at the equilibrium position with increasing temperature, which leads to the increase in the width of the peaks as seen in Figure . Moreover, the different location sites are due to the fact that the optimized lattice constant by the MEAM is smaller than that by the first-principles calculations.

Figure 3

Pair correlation function g(r) of BeO calculated by DFT at 0 K and MD simulation at 300 K.

Predictions of Thermodynamic Properties

Next, we consider the linear expansion of the BeO structure. The expansion can reveal the size change of materials with the temperature, which is very important for the application for materials. In the atomic point of view, it responds to the rate of change in kinetic energy of the atom with the temperature. The expansion coefficient α can be obtained by the formulawhere P is the pressure. Figure shows the increase in linear expansion as a function of temperature averaged over 10 randomly generated 6 × 4 × 4 supercell. The calculated results by the MEAM potential are scaled by a factor of 0.8 to account for the possible effects caused by the impurities in the experimental samples, the experimental errors, and the difference between the idealized atomic scale models and the bulk material of BeO in the whole temperature region. The experimental data[35,36] and other theoretical calculations[25] are also included. At the temperature range from 300 to 1200 K, a linear trend is observed. Then, the slope of increase changes at the temperature range between 1100 and 1800 K. After 1800 K, the slope changes again, and the experimental data also show some differences, which may be caused by the defect in the samples. The MEAM results with the abovementioned scaling factor show a good agreement with the experimental data.

Figure 4

Linear expansion coefficient of BeO as a function of temperature. Experimental results are also plotted for comparison.

Figure shows the enthalpy increase as a function of temperature. Similar to the lattice constant, the enthalpy increases linearly with the temperature increase at the low-temperature region. However, there is a sharp increase in the enthalpy at the temperature range between 500 and 2000 K. Using the Bunsen ice calorimeter and the drop method, Victor and Douglas[37] measured the enthalpy of BeO relative to 273 K. Besides, the enthalpy of BeO is measured by Greenbaum et al.[38] at the high-temperature region. To better compare our results with the experimental data, the enthalpy for the BeO is calculated relative to room temperature (300 K), as shown in Figure . The calculated values show a good agreement with the experimental results,[30] especially at the low- and high-temperature region. The fitted experimental data is the mean enthalpy of the sample based on the least squares. At the middle-temperature region, the calculated results and experimental data show a small difference. However, if we consider the density of the single crystal in the experimental samples and the experimental error, then it is a reason for the existing small difference.

Figure 5

Enthalpy of BeO as a function of temperature relative to room temperature.

The specific heat capacity of a material is defined as the heat required to raise the temperature of unit mass of a given substance by 1°. It should be mentioned that the temperature dependence of specific heat capacity is very important for BeO to be applied as a nuclear fuel additive. The heat capacity Cp is then calculated using the first derivative of enthalpy H as a function of temperature at constant pressure with the following formulawhere n is the number of mole and T is the temperature. Figure shows the calculated heat capacity as a function of temperature. Around room temperature, there are some differences between our results and the experimental data, and the capacities at the low-temperature region calculated by the MD are overestimated than the experimental data. The difference may come from the derivation process since the small difference in enthalpy may lead to a big difference in the heat capacity data. However, at the temperature range from 600 to 2800 K, our results calculated using the fitted potential show a good agreement with the experimental[37,39] and theoretical data.[20]

Figure 6

Heat capacity of BeO as a function of temperature.

Heat transfer is one of the most important properties for BeO. The experimental[3,40,41] and theoretical studies[24] show that the thermal conductivity of BeO is about 300 W/(m K) at room temperature, which can be comparable to that of the metal materials, such as aluminum and silver. Unlike the metal materials of which heat is mainly transported by electrons, the BeO system transports heat mainly by phonons since BeO is an insulating material with a large band gap of about 10 eV.[5,6] After the calculation of the lattice constant, elastic constant, expansion, enthalpy, and heat capacity, the reliability of the potential has been verified. In this section, the potential is used for the calculation of thermal conductivity. In the MD simulations, the lattice thermal conductivity can be obtained by the EMD method with the following relationshipwhere V is the system volume, kB is the Boltzmann constant, T is the system temperature, J is the heat current, τm is the integration time, and the angular brackets denote an average over time in the MD simulation. Figure shows the calculated thermal conductivity as a function of temperature. For a better comparison, the experimental data and other theoretical results are also plotted in Figure . In the whole-temperature region, the thermal conductivity decreases with increasing temperature, which is not like the metal material. This is reasonable because, in the BeO system, the phonon transport is dominant and the u process plays an important role in heat transport. The calculated results by the MEAM potential are scaled by a factor of 0.8 to account for the possible effects caused by the impurities in the experimental samples, the experimental errors, and the difference between the idealized atomic scale models and the bulk material of BeO in the whole-temperature region. We can see that our results show the same trend as the experimental data. The MEAM results with the abovementioned scaling factor show a good agreement with the experimental data. Compared with the experimental data, our results are reasonable since the experimental sample usually has some defects such as vacancy, impurity, and grain boundary, which block the heat transfer.

Figure 7

Thermal conductivity of BeO as a function of temperature.

Conclusions

The MEAM potential parameters are developed to describe the interactions of atoms for the BeO system. The first-principles calculations and experimental data are used for the determination of the potential parameters. The fitted potential shows good performance on the lattice constant and elastic constant. Besides, the temperature dependence of dynamic properties such as expansion, enthalpy, and heat capacity calculated by this potential is also tested and compared with the reported experimental and theoretical data. The good agreements between the calculated results and experimental data further validate the potential models. Finally, the potential has been used to predict the thermal conductivity of BeO, which is very important for its heat transfer applications in the energy field. The calculated results show a good agreement with the experimental data as a function of temperature. Furthermore, the potential can be used to study large systems with defects and help to understand the dynamic progress from the atomic aspect.

Methods

The structure optimization of BeO is calculated by using a plane wave pseudopotential formulation[42−44] within the framework of the DFT. The code is implemented in the Vienna ab initio simulation package (VASP). The exchange–correlation energy is in the form of PBE function[32] with the generalized gradient approximation (GGA). The cutoff energy for the plane wave expansion is 400 eV. The k points are sampled on a uniform grid of 15 × 15 × 7. The system is fully relaxed until the magnitude of the forces acting on all the atoms becomes less than 0.001 eV/Å. For the heat transport, the lattice thermal conductivity is predicted using the equilibrium molecular dynamics (EMD) simulations combined with the Green–Kubo autocorrelation decay method.[45] The MEAM potential[46] is used to describe the interatomic interactions, and the time step is set to 0.5 fs. After a constant temperature simulation of 500,000 steps and a constant energy simulation of 400,000 steps, the heat current data are collected to calculate the heat current autocorrelation function.

To describe the atomic interaction in a BeO system, the MEAM potential is applied to describe the atomic interaction. Here, we applied the potential with the form as in ref (46), which is successfully used to describe the metal[47,48] and covalent compounds.[49] In the MEAM potential, the total energy of the system can be expressed as[48,50−52]where the sums are over the atoms i and j, F is the embedding function, ρ is the background electron density of atom i, S is the screen parameter between atom i and atom j with the distance R, and φ is the atomic interaction between two atoms. Thus, the energy of a single atom i can be expressed aswhere Zis the number of nearest neighbors in the structure. It should be mentioned that the function F has many different forms, such as power law, Johnson, VBO, and Baskes. However, for the pair interactions φ, there is no mathematical form. To obtain the φ between two atoms, the reference structure was constructed with the known energy. Using the energy expression (eq ), the energy of the reference structure can be expressed aswhere E can also be obtained by the zero-temperature universal equation of state by Rose et al.[53] as a function of R.where d is an adjustable parameter andandwhere B is the bulk modulus, and Ω is the equilibrium atomic volume. In our work, we also used the first-principles calculations to calculate the value of E, which could provide a reliable accuracy in determining the parameter α and re.

Combined with eq , the pair interaction φ can be derived asThen, we can obtain thatTherefore, the energy of the system can be rewritten as a function of the density of the background electron. In this work, the embedding function of Baskes form is usedwhere A is a fitting parameter and Ec is the cohesive energy of the compounds. The background electron density ρ̅ is composed of a spherically symmetric partial electron density ρ0 and angular contributions ρ1, ρ2, and ρ3. ρ0 has no angular dependency of electron densitywhere r is the distance between atom i and other atoms. The angular contributions can be expressed as followswhere α, β, and γ represent three different directions. The above formulas show that ρ has an exponential relationship with r. Then, the electron density can be expressed as the following formwhere β is the attenuation coefficient. Here, we assume that the contribution to the total electron density coming from the angle item is small. Under this assumption, the total background electron density can be expressed aswhere G(Γ) can be described aswhere Γ can be expressed aswhere t is the fitted parameter related to the electron density. In the MEAM, a screening function S is introduced to describe the cutoff for the many-body interactions. For arbitrary two atoms i and j in the compounds, the screening function S can be expressed aswhere S is the screening factor that denotes the influence of neighbor atoms k to the atoms i and j. To obtain the value of S, the following geometric construction is adopted from ref (54). Imagine an ellipse on the x and y plane, passing through atoms i, k, and j with the x axis of the ellipse determined by atoms i and j. The equation of the ellipse is given bywhere the value of the parameter C is related to the coordinate positions of the three atoms i, j, and k, as followswhere X = (R/R)2 and X = (R/R)2. The screening factor S is defined as a function of C as followswhere Cmin and Cmax are the limiting values of C determining the extent of screening, and the smooth cutoff function isIn the following part, the physical meaning of the screening factor is explained. Cmax and Cmin determine the maximum and minimum area of the ellipse, respectively. If the atom k is outside the maximum area, then the atom k has no effect on the interaction between atoms i and j (S = 1). If the atom k is inside the minimum area, then atoms i and j have no interactions at all (S = 0). If the atom k is located at the position between the minimum and maximum area, then the interactions between atoms i and j change gradually with the coordinate position of atom k.