General Papers
Molecular Dynamics Simulation of the Behavior of Beryllium Diffusion in Corundum
2015 Volume 14 Issue 4 Pages 111-116
Details
2015 Volume 14 Issue 4 Pages 111-116
Molecular dynamics (MD) simulations are highly useful for analyzing atomic behavior during diffusion, especially in systems that are difficult to investigate experimentally. The focus of the present study was the diffusion behavior of Be in corundum, which was analyzed by MD calculations. First, we derived new potential parameter sets for O, Al, and Be. This parameter set was verified to well reproduce the structures and properties of corundum, bromellite, and chrysoberyl. Based on MD simulations of corundum containing Be as interstitial atoms, where the simulations were performed using the newly derived potential parameters, the diffusion coefficient was estimated to be approximately 10−7 cm2/s at around 2100 K. This is consistent with previously published experimental results, which confirms the validity of the MD simulation. The present calculations also reveal the detailed atomic movement, where Be atoms jump between Al sites and/or interstitial sites, and that the activation energy of this process is approximately 1.1 × 102 kJ/mol.
Diffusion is an important process with respect to understanding the properties of materials (e.g., Zhang and Charniak [1]). In many cases, the diffusion processes are investigated experimentally. However, cases exist where diffusion is difficult to analyze by experiments. For example, the case where the diffused atoms are light elements and are difficult to detect using conventional analytical methods such as X-ray fluorescence or even electron microprobe analysis. Another example is the case where the diffused atom is poisonous; in such cases, experimental analysis is best avoided to the extent possible. In such cases, molecular dynamics (MD) simulations are highly useful for understanding the diffusion behavior. Of course, MD simulation can also reveal atomic movement during diffusion. As an example of these cases, we focus on the diffusion of Be in corundum.
Transparent stones of corundum, Al2O3, are often used as gemstones (Hughes [2]). Red stones of corundum are named ruby, which contains a few percent of chromium, whereas stones that exhibit other colors are named sapphire. Sapphire exhibits a wide variety of colors because of slight differences in its impurity content, and its color can be artificially changed by introducing a certain element from the outside. Therefore, in the field of gemology, the ability to identify stones with natural color and to distinguish them from artificially colored corundum is significant. An example of a color change induced by "diffusion treatment" is the conversion of a colorless sapphire to a deep-blue sapphire by the diffusion of Ti. In this case, Ti is easily detected using conventional X-ray fluorescence analysis and this treatment can be easily identified. Conversely, in other cases, a color-change treatment is performed by diffusion of Be. In this treatment, Be results in an attractive pinkish-orange color in combination with other impurities such as Mg2+ (Emett et al. [3], Pisutha-Arnond et al. [4]). In this case, the Be can be detected only using expensive instrumental methods such as laser ablation inductively coupled plasma mass spectrometry; it is difficult to detect using methods conventionally used for gem testing. Furthermore, Be is poisonous, therefore experimental analysis should be avoided as much as possible. Consequently, understanding the physical properties of Be diffusion in corundum using computational techniques is important. Therefore, in the present work, we performed molecular dynamics (MD) simulations of the behavior of Be in corundum. To our knowledge, the present study represents the first application of a computer simulation to gemology.
MD calculations were performed using the MXDTRICL program [5] with the two-body central force interatomic potential
| (1) |
| Atom | Z[e] | A[Å] | B[Å] | C[kcal1/2 Å3 mol−1/2] |
| O | −1.500 | 1.755 | 0.160 | 20.000 |
| Al | 2.250 | 1.295 | 0.099 | 0.000 |
| Be | 1.500 | 0.955 | 0.080 | 0.000 |
We used the Verlet method for integration of equations of atomic motion and the Ewald method for summation of electrostatic interactions. Temperature and pressure were controlled by scaling atomic velocities and simulating cell parameters, respectively. The time step was 2 fs throughout our simulations. Periodic boundary conditions for corundum were imposed with a basic MD cell of 72 crystallographic units of corundum(aMD = 6a, cMD = 2c). In the case of pure Al2O3, this calculation cell contains 2160 atoms. For the simulation of bromellite and chrysoberyl, basic MD cells were composed of 196 unit cells (896 atoms) and 32 unit cells (896 atoms), respectively.
In Table 2, the calculated structural parameters and bulk modulus of corundum, bromellite, and chrysoberyl, calculated at 300 K and 1 atm, are compared with the experimentally obtained values. In these calculations, 10,000-step run (20 ps) at 300 K were performed after 10,000 steps of preliminary annealing of the initial structure reported in the literature (Hazen and Finger [7]). These results show that the parameter sets derived in the present study reproduce the experimental results very well.
| Al2O3 | BeO | |||
| MD | Exp.a,b | MD | Expc | |
| a [Å] | 4.767 | 4.758 | 2.712 | 2.696 |
| b [Å] | ||||
| c [Å] | 12.942 | 12.991 | 4.317 | 4.379 |
| Density [g/cm3] | 3.988 | 3.989 | 3.019 | 3.017 |
| Bulk Modulus [GPa] | 238 | 252 | 212 | 212 |
| Melting point [K] | 2350 | 2320 | ||
| BeAl2O4 | ||||
| MD | Expd | |||
| a [Å] | 4.439 | 4.424 | ||
| b [Å] | 9.415 | 9.369 | ||
| c [Å] | 5.520 | 5.471 | ||
| Density [g/cm3] | 3.656 | 3.75 | ||
| Bulk Modulus [GPa] | 169 | 242 | ||
a: Hughes, [2]
b: Singh et al. [6]
c: Hazen and Finger [7]
d: Hazen and Finger [8]
To further ensure the validity of the newly derived parameter set, we also calculated the temperature dependence of the corundum structure. The temperature was increased from 300 K at a constant pressure using the optimized structure at 300 K as the initial structure. Crystallographic properties were obtained from calculations over at least 10,000 steps (20 ps) after 10,000 annealing steps. Figure 1 shows the temperature dependence of the structure parameters; these results show that the calculation results obtained with the newly derived parameter set are consistent with the experimental results. Furthermore, the structure could not be maintained at 2400 K, which means that the experimentally observed melting point (2320 K) is also reproduced well. These results indicate that the newly derived parameter set is sufficiently precise to be useful for the calculation of the diffusion of Be in corundum.
Temperature dependence of cell parameter of corundum. MD calculated values are shown with solid circles, and experimental data (Fiquet et al. [9]) are open circles.
The crystal structure of corundum belongs to the so-called "haematite" group. The structures of minerals of this group are described on the basis of hexagonal close packing of oxygen atoms, with the cations in octahedral coordination between them (Hughes [2]). However, only 2/3 of the octahedral sites are occupied by Al3+ cations. Therefore, this structure contains relatively large spaces in the octahedral sites that are not occupied by Al3+; in the present study, these spaces are referred to as "interstitial sites" hereafter. Figure 2 shows the structure of corundum, as illustrated by the program VESTA (Momma and Izumi [10]).
Atomic arrangement of corundum structure: projection along the c axis (top) and the stacking manner of Al and O layers (bottom). Interstitial sites are shown with dotted circles. These are illustrated after Hughes [2]. (Color online)
For the calculations of Be diffusion, Be was added in the interstitial sites because an ion appears to be able to easily enter this space. In the present study, to maintain overall charge neutrality, two Be atoms were added and two Al atoms and one O atom were omitted (i.e., vacancies were created). Two Be atoms and three vacancies were set separately as far as possible in order to avoid the influence of the interaction between two Be atoms.
The behavior of Be diffusion was analyzed according to the method of Tsuchiyama and Kawamura [11]: The 50,000-step runs (100 ps) were repeated ten times after preliminary annealing of the initial structure for 10,000 steps (20 ps). As the initial state, we used the structure in which Be atoms and vacancies of Al and O are introduced into the relaxed corundum structure calculated at a given temperature. The diffusion coefficient D at each run was obtained from diffusion distance x, which was derived from average mean square displacement of Be, by using the following equation (Zhang, [12]);
| (2) |
In addition, the temperature dependence of diffusion coefficient D is expressed as
| (3) |
| (4) |
In Figure 3, the logarithm of D is plotted against 1/T, which is expected to exhibit a linear relationship. Therefore, a straight line was fitted to these data to obtain the activation energy E and pre-exponential factor D0 for Be diffusion from the line's slope and intercept, respectively, using equation [4]:E is 1.1 × 102 kJ/mol and D0 is 1.0 × 10−4 cm2/s.
Relationship between the calculated diffusion coefficients and temperature.
From this relationship, the diffusion coefficient at a certain temperate was estimated. Previously, Emmett et al. [3] performed heating experiments of corundum with chrysoberyl powder and derived an approximate diffusion coefficient of 10−7 cm2/s at 1800 °C. Based on the present calculations, the diffusion coefficient D at the same temperature was estimated to be 1.9 × 10−7 cm2/s, which is completely consistent with the experimental value. This result confirms the validity of the MD simulation and may also indicate that the rate-determining step is the diffusion of Be in corundum itself, and not one of the other processes such as decomposition of chrysoberyl.
3.3 Diffusion mechanismThe analysis of the movement of Be atoms during diffusion shows that a Be ion jumps between the Al sites and interstitial sites. A Be ion remains and vibrates around one of these sites for a certain period and then suddenly jumps to the next site; this process is repeated numerous times. Therefore, the velocity of diffusion depends on the length of time a Be ion remains at one site and on the frequency of jumps between these sites: As the staying time becomes shorter and the jump frequency becomes larger, the diffusion velocity increases. However, in many cases, an Al ion is present at the site to which a Be ion will jump, as shown in Figures 2 and 4a. Detailed observation of the trajectory of ions revealed that the diffusion of Be proceeds via the following mechanism. (1) An Al ion moves to another interstitial site where no ion is present. In many cases, Al moves along the c axis, because an interstitial site exists also above or below that Al ion along the c axis shown in Figures 2 and 4a. At the same that the Al ion moves, a Be ion jumps to the vacated Al site (Figure 4b). (2) This Be ion then jumps again to next site, and simultaneously the Al ion migrates back to its original site (Figure 4c). This process is schematically illustrated in Figures 4a-c and the trajectories of diffused Be atom and Al atoms moving associated with Be diffusion are shown in Figure 4d. Al atoms just vibrate at around the original position when Be is not located at its nearest site, and Al ion which is forced out to the next site by Be also moves back to the original position after Be diffuses.
(a)-(c) Schematic drawing of detailed atomic movement of Be diffusion in corundum. Be and Al ions move along solid and dotted arrows, respectively. (d) Trajectories of Be and Al atoms simulated with the calculations at 1900 K, are shown with purple and light blue lines, respectively. The position of atoms is the snapshot after the movement labeled as No. 2. It is shown that Be diffuses out of the area enclosed with yellow lines after the movement No. 3. (Color online)
Our result of the calculated diffusion coefficient being approximately the same as the experimental value indicates that this kind of atomic process may actually occur during diffusion of Be in corundum as the rate-determining step with an activation energy of approximately 1.1 × 102 kJ/mol. As mentioned above, an Al atom hardly diffuses in this process. However, when a vacancy of Al site exists at the nearest, an Al atom can move to that vacancy. Especially at high temperature, there is the case that Al vacancy diffuses to the site near Be atom, which may affect the Be diffusion. The examination of the effect of vacancies on the Be diffusion should be a future subject.
MD simulation was performed to understand the property of diffusion of Be in corundum. A potential parameter set for O, Al, and Be was newly derived empirically, and was ensured to reproduce very well the crystallographic parameters and bulk modulus at 300 K of corundum, bromellite, and chrysoberyl, and also temperature dependence of these parameters and melting point of corundum.
MD calculation of corundum containing Be atom shows that Be diffuses with a diffusion coefficient of approximately10−7 cm2/s at 2100 K. This result is consistent with the experimental results, which confirms the validity of MD simulation. During the diffusion, Be atoms jump between Al sites and/or interstitial sites with an activation energy of approximately1.1 × 102 kJ/mol.
The present study was conducted in part at the Institute of Gemology and Jewelry Art Yamanashi prefectural government. The authors thank Dr. Yasushi Takahashi and Dr. Takaya Nagai for their support.
