New B2O3 Crystals Predicted from Concurrent Molecular Dynamics

13712

J. Phys. Chem. C 2007, 111, 13712-13720

New B2O3 Crystals Predicted from Concurrent Molecular Dynamics Simulations and First-Principles Calculations Liping Huang,†,‡,§ Murat Durandurdu,| and John Kieffer*,† Department of Materials Science and Engineering, UniVersity of Michigan, Ann Arbor, Michigan 48109-2136, Department of Chemical and Biomolecular Engineering and Department of Physics, North Carolina State UniVersity, Raleigh, North Carolina 27695, and Department of Physics, UniVersity of Texas at El Paso, El Paso, Texas 79968 ReceiVed: May 9, 2007; In Final Form: July 9, 2007

Molecular dynamics (MD) simulations, based on a new coordination-dependent charge-transfer potential, were used to study the behavior of crystalline B2O3 in response to various thermal and mechanical constraints. This interaction potential allows for the charges on atoms to redistribute upon the formation and rupture of chemical bonds and dynamically adjusts to multiple coordination states for a given species. Our MD simulations predict that upon isotropic expansion and compression, the B2O3-I crystal transforms into new low- and highdensity B2O3 crystals, the stability of which we have further verified using first-principles calculations. The low-density B2O3 crystals (B2O3-0) provide a key to understanding the anomalous thermomechanical behaviors of vitreous B2O3 and the crystallization anomaly of this compound. The high-density B2O3 crystal (B2O3-III), predicted from concurrent MD simulations and first-principles calculations, is different from the known highpressure phase of B2O3-II crystal, even though the bonding is the same in these two phases. B2O3-III is characterized by a higher energy than B2O3-II at low pressures, but upon further compression the energies of these two phases become indistinguishable. The transformation from B2O3-I to B2O3-III appears to be kinetically favored, especially at low temperatures. Our studies indicate that the phase diagram of B2O3 is much richer than previously known.

1. Introduction Two crystalline forms are generally known for boron oxide, that is, B2O3-I, which is stable at ambient pressures and characterized by an infinite chain or ribbon of BO3 triangles,1 and B2O3-II, which is stable at higher pressure and consists of a framework of linked BO4 tetrahedra.2 The pressure-temperature phase diagram of B2O3, including relative stability of these phases and the melting curve, has been rather scantly studied in experiments.3-5 A tentative phase diagram was constructed by Mackenzie and Claussen3 in 1961 up to 1100 °C and 90 000 atm. They found that crystalline B2O3 was readily obtained from the glassy phase at high pressures and temperatures. At atmospheric pressure, however, the transitions between B2O3I, B2O3-II, and liquid phases are extremely sluggish. In fact, as determined by high-temperature X-ray diffraction experiments, the B2O3-II to B2O3-I transformation cannot be observed at ambient pressure, even when providing ample thermal activation by heating the crystal as high as the melting temperature of B2O3-I. Only recently did Brazhkin et al.5 observe the B2O3-II to B2O3-I transition for the first time in situ during pressure release at high temperature, thereby confirming that B2O3-II is indeed thermodynamically stable at higher pressures than B2O3I. Furthermore, boron oxide has never been observed to crystallize from a dry melt at ambient pressure. Even if the melt * To whom correspondence should be addressed. E-mail: kieffer@ umich.edu. † University of Michigan. ‡ Department of Chemical and Biomolecular Engineering, North Carolina State University. § Department of Physics, North Carolina State University. | University of Texas at El Paso.

is seeded with crystals and the crystals are melted back a bit at temperatures above TM ≈ 450 °C, no crystal growth is observed at any imposed undercooling.4 This behavior has been termed the “B2O3 crystallization anomaly”. Aziz et al.6 studied the crystal growth kinetics of boron oxide under pressure. Using the published thermodynamic information on the B2O3 system and a crude free-energy model for the crystal and glass phase, they explained the B2O3 crystallization anomaly based on twodimensional nucleation of monolayers being the rate-limiting step; the frequency of two-dimensional nucleation is negligible at all temperatures and at pressures less than 1 GPa.6,7 However, the microscopic mechanisms for the transformations between B2O3-I, B2O3-II, and liquid phases still by and large remain puzzling. On the other hand, most of the theoretical studies, either using classical molecular dynamics simulations8,9 or first-principles calculations,10-13 have been focusing on the structural and electronic properties of equilibrium B2O3-I and B2O3-II crystals. Very little research has been carried out to study the structural transformations or the behavior of crystalline phases under various thermal and mechanical conditions. Employing a reduced cell volume with full relaxation of all internal coordinates using ab initio total-energy pseudopotential calculations, Takada et al.10 studied the structural transformation from the polymorph containing the BO3 triangular unit into that containing the BO4 tetrahedral unit. Their study showed that the BO3 triangular structural unit in the minimized structure for the 60% cell volume turns into a BO4 tetrahedron. This corresponds to a pressure-induced phase transition. Although the original cell is only isotropically compressed and the final structure is not completely the same as B2O3-II, it is consistent with the

10.1021/jp0735583 CCC: $37.00 © 2007 American Chemical Society Published on Web 08/28/2007

B2O3 Crystals from MD and First-Principles Calculations

J. Phys. Chem. C, Vol. 111, No. 37, 2007 13713

TABLE 1: Crystallographic Data of B2O3-I and B2O3-II from Literature space group

crystal a

a

b

c

fractional coordinates

(Å)

(Å)

(Å)

B2O3-I

P31

4.336

4.336

8.340

B2O3-Ib

P3121

4.336

4.336

8.340

B2O3-IIc

Ccm21

4.613

7.803

4.129

a

B1 B2 O1 O2 O3 B1 . O1 O2 B1 O1 O2

x

y

z

.2229 .8281 .5468 .1485 .0045 3954 .6009 .1607 .1606 .3698 .2475

.3926 .6031 .3972 .6004 .1608 .2299 .1477 .0000 .1646 .2911 .0000

.0198 .0921 .0000 .0775 .1291 .2244 .1282 .3333 .4335 .5802 .5000

Z)3 3 B2O3/cell

Z)3 3 B2O3/cell Z)4 4 B2O3/cell

Ref 1. bRef 26. cRef 2.

TABLE 2: Crystallographic Data of B2O3 Crystals at 0 K from Our First-Principles Calculations P (GPa)

a

c

fractional coordinates

space group

(Å)

(Å)

(Å)

B2O3-I

P3121

4.361

4.361

8.360

18

B2O3-III

P3121

4.107

4.107

6.955

-8

R-B2O3-0

C121

8.415

4.823

12.095

β) 91.7°

-12

β-B2O3-0

C121

8.804

5.050

12.225

β) 91.8°

0

crystal

b

observed phase diagram in that the 4-fold BO4 structural unit is more stable than the 3-fold BO3 structural unit at high pressure. However, the lattice vectors were not optimized in the study by Takada et al.,10 and the temperature effect was not taken into account either, which may be crucial for certain structural transformations to take place. In this study, we have first employed classical molecular dynamics (MD) simulation to investigate the behavior of crystalline B2O3 under various thermal and mechanical conditions. Our objective was to carry out virtual experiments, that is, subjecting the material to compressive or tensile stress, and observe the structural transitions that ensue without preconception of the resulting polymorphs. Concurrently, we verified the stability of the newly discovered phases using first-principles calculations. To simulate a system like B2O3 in which the structural building block can exhibit multiple coordination states, it is essential to use a potential model that allows for coordination changes during simulations. We have successfully designed such a coordination-dependent charge-transfer reactive force field for the description of mixed ionic-covalent interactions, as occur in silica and boron oxide. The three main features of

B (6c) O (6c) O (3b) Z ) 3, 3 B2O3/cell B (6c) O (6c) O (3b) Z ) 4, 3 B2O3/cell B (4c) B (4c) B (4c) O (4c) O (4c) O (4c) O (4c) O (2a) Z ) 3, 6 B2O3/cell B (4c) B (4c) B (4c) O (4c) O (4c) O (4c) O (4c) O (2b) Z ) 3, 6 B2O3/cell

x

y

z

.7620 .4441 .1585

.1615 .8437 .0000

.0592 .0418 .8333

.8364 .2257 .8567

.3063 .6213 .0000

.8952 .9194 .8333

.2667 .4842 .2491 .4056 .8508 .7473 .2499 .0000

.2402 .5091 .6370 .3839 .1580 .8545 .1840 .9929

.7789 .1132 .4458 .8136 .1486 .4775 .6670 .0000

.7432 .5020 .7544 .6208 .1264 .2514 .7344 .0000

.2630 .5439 .6850 .4097 .1706 .9214 .2025 .0627

.7204 .3858 .0538 .6699 .3346 .0095 .8320 .5000

this potential model are the following: (1) dynamic charge redistribution, a charge-transfer term controls the extent of charge polarization in a covalent bond, as well as the amount of charge transferred between atoms upon rupture or formation of such a bond; (2) conditional three-body interactions, the directional character of the covalent bonding, modeled by means of three-body terms, is coupled to the degree of covalency in atomic interactions and vice versa; and (3) variable coordination number, both two- and three-body interactions depend on the effective number of nearest neighbors of an atom, which is evaluated dynamically based on the local environment of this atom. The detailed formulation of this potential and the values of potential parameters can be found in our previous publications.14,15 Using this model, we have significantly advanced the understanding of densification mechanisms in inorganic compounds, which are relevant in the context of materials design and of geological processes in the Earth’s interior.16 On the basis of this potential model, we have successfully reproduced the known anomalous thermomechanical behavior of B2O3, that is, the minimum in its mechanical modulus in the molten state, and have provided a mechanistic explanation.17 Furthermore,

13714 J. Phys. Chem. C, Vol. 111, No. 37, 2007

Huang et al.

Figure 1. Density evolution with pressure from MD simulations at 300, 500, and 700 K and from first-principles calculations at 0 K. Pressures are hydrostatic and negative magnitudes indicate isotropic tensile stress.

our previous MD simulations predicted a resurfacing of this anomaly under tensile stress and ambient temperatures in vitreous B2O3 and demonstrated that the underlying mechanisms are universal for all network-forming glasses.15 The transferability of our potential has been further tested on crystalline B2O3, as described in this paper. The concurrent use of firstprinciples calculations and MD simulations proved to be highly synergistic in terms of predicting previously unknown structural states providing validation of these predictions. In this paper, we present new insights concerning the phase diagram of B2O3 system and a better understanding of the transformation mechanisms between different phases of B2O3. 2. Computational Details We performed the first-principles calculations within the density functional theory formalism and the generalized-gradient approximation (GGA), using the Perdew-Burke-Ernzerhof functional18 for the exchange-correlation energy. The calculations were carried out with the SIESTA program19 using a linear combination of atomic orbitals as the basis set, and normconserving Troullier-Martins pseudopotentials.20 An optimized split-valence double-ζ plus polarized basis set was employed. A real space grid was equivalent to a plane wave cutoff energy of 120 Ry. Γ-point sampling for the supercell’s Brillouin zone integration was used, which is reasonable for a cell with 120 atoms. For each value of the applied stresses, the lattice vectors were optimized together with the atomic coordinates by using the conjugate-gradient technique. The original form of Parrinello-Rahman constant pressure algorithm21 and the modified form proposed by Wentzcovitch22 and Cleveland23 were used to maintain a specified constant pressure in our MD simulations. The dynamics in the latter constant pressure algorithm has a symmetry invariant with respect to the primitive cell. This is not the case in the original form of Parrinello-Rahman constant pressure algorithm. However, the behaviors observed using these two constant pressure algorithms were similar. Our MD simulations were carried out by using velocity scaling to control temperature. The equations of motion were integrated using the fifth-order Gear predictorcorrector algorithm24 with a time step of 2 fs. MD simulations were carried out on systems of 64 unit cells of B2O3-I (15 atoms per cell) with periodic boundary conditions. Hydrostatic pressure was applied to each system in increments of 1 GPa, continuously ramped at a rate of 0.05 GPa/ps, followed by a 20 ps equilibration

Figure 2. Comparison of the configurations at different pressures from first-principles calculations (left side) and MD simulations (right side), viewed in the [010] direction of B2O3-0.

period in NPT ensembles (isobaric-isothermal). Structural and dynamic properties were calculated based on the equilibrium structure at each state point. The bulk modulus was calculated directly from the equation of state according to B ) F(dP/dF). The symmetry and unit cell parameters of the simulated structures were determined by KPLOT (a program for identifying the space group of a crystal and visualizing the structure).25 3. Results and Discussion Gurr et al.1 determined that B2O3-I is of space group P31 with cell dimensions a ) 4.336, and c ) 8.340 Å. The cell contains three B2O3 units. The structural unit is an infinite chain or ribbon of BO3 triangles. These ribbons are interconnected, so that an O atom has four O neighbors at 2.3 to 2.43 Å distance. All boron atoms are trigonally coordinated by O at bond lengths between 1.34 and 1.40 Å. Some controversy exists in the literature for the space group classification of B2O3-I. Recently, Effenberger et al.,26 based on a refined analysis of previously published X-ray data,1 concluded that the B2O3-I structure belongs to space group P3121 (No. 152), instead of P31 (No. 144) as suggested earlier.1 In the P3121 structure, all boron atoms are equivalent and only two nonequivalent O atoms are present, which are coordinated



Figure 3. Coordination number (a) and g(r) (b) of crystalline B2O3 phases at -10, 0 and 20 GPa.

by two B atoms. In the previously suggested less symmetric P31 structure,1 the B atoms can be separated into types B1 and B2, which are bonded to three nonequivalent O atoms and form two slightly different BO3 units (see Table 1). On the basis of the comparison between calculated and experimental structure parameters and the calculated lattice energies of the two structures at the density functional theory level, Islam et al.27 concluded that the ground-state structure of B2O3 has a P3121 symmetry. The initial configurations used in our MD simulations and first-principles calculations are generated from atomic coordinates in the unit cell of B2O3-I, which belongs to the trigonal space group P31.1 After a certain equilibration time in MD simulations (300 K and 0 GPa) or first-principles calculations (0 K and 0 GPa), the ending structures have a space group of P3121 as determined by KPLOT.25 Our concurrent MD simulations and first-principles calculations therefore confirm the conclusion by Islam et al.27 about the ground-state space group of B2O3-I. The unit cell information and the fractional coordinates of atoms in the asymmetric unit from our first-principles calculations are listed in Table 2. When subject to isotropic tensile stress up to -7 GPa and compressive stress up to 17 GPa in our MD simulations at 300 K, B2O3-I crystal transforms into distinct low- and high-density phases, as indicated by the abrupt changes in the density versus pressure curve in Figure 1. Similar transformations are observed in our first-principles calculations at 0 K. In our MD simulations, the transformation

pressure in the compressive regime decreases with increasing temperature (from 17 to 12 GPa when the temperature changes from 300 to 700 K), while that in the tensile regime stays almost the same at the three temperatures in our study. A sequence of configurations at different pressures from our first-principles calculations at 0 K and MD simulations at 300 K are shown side by side in Figure 2. The coordination number and pair correlation functions, g(r), of B2O3 phases at -10, 0 and 20 GPa are plotted in Figure 3. Figure 3a shows that in the lowdensity phase at -10 GPa, the B atom is three-coordinated by O atoms and the O atom has four nearest O neighbors, which are the same as in B2O3-I at 0 GPa. The B-O and O-O distances in the low-density phase are slightly larger than those in B2O3-I (Figure 3b). In the high-density phase at 20 GPa, the B atom is four-coordinated by O atoms and the O atom has nine nearest O neighbors as seen in Figure 3a. Figure 3b shows that the first peak corresponding to g(r)B-O at 20 GPa clearly splits: one B-O bond of short distance (∼1.4 Å) and three B-O bonds of long distance (∼1.5 Å) in each BO4 tetrahedron. Figure 3 shows that the observed low-density B2O3 phase is characterized by BO3 triangles, the same as in B2O3-I, while the high-density phase consists of a framework of linked BO4 tetrahedra, the same as in B2O3-II. Neither the low-density nor the high-density phases revealed by our simulations have been previously observed in simulations or experiments. In the following, we describe in detail the characteristics of these newly predicted B2O3 polymorphs.


Figure 4. Change in density (a) and bulk modulus (b) of B2O3 glass and B2O3-0 crystals under pressure. Stability ranges of R-B2O3-0 and β-B2O3-0 are delineated.

We find as a result of isotropic expansion of B2O3-I that the low-density phase has a monoclinic space group of C121 and a density of 1.810 g/cm3 under ambient conditions. We tentatively call this phase B2O3-0 in light of the succession it assumes in terms of density relative to B2O3-I and B2O3-II. To the best of our knowledge, Cole and Taylor28 were the only ones to ever report an experimental observation of a low-density crystalline B2O3 polymorph. The phase they described, which they obtained by vacuum dehydration of H3BO3, had a density of 1.805 g/cm3, that is, close to that of B2O3-0. However, our B2O3-0 phase is monoclinic with 6 B2O3 molecules per unit cell (see Table 2), whereas theirs is cubic with 16 B2O3 molecules per unit cell. That others have not been able to repeat the Cole and Taylor experiments remains puzzling. Perhaps the starting materials are important. In light of our predictions reported here, further experimental studies of low-density B2O3 states seem indicated. Upon releasing the tensile stress, B2O3-0 does not revert to B2O3-I but instead undergoes another structural transition at -2 GPa in MD simulations that results in the partial collapse of network rings but almost no detectable densification beyond what can be associated with elastic deformation (Figure 4). The configurations in Figure 2 show the low- and high-density modifications of B2O3-0 viewed along the [010] direction. The partial collapse of network rings is particularly obvious for the smaller ring seen from this perspective. In fact, the resemblance to the ring geometries in β- and R-cristobalite is remarkable,29

Huang et al. and we therefore refer to the low- and high-density modifications of B2O3-0 as β-like and R-like (configurations in Figure 2 at -10 and -2 GPa in MD, respectively). Similar to β-cristobalite, β-B2O3-0 not only has more symmetric rings but also the higher bulk modulus of the two modifications as can be seen in Figure 4. Similar to cristobalite silica, the transformation between R-B2O3-0 to β-B2O3-0 can not only be introduced mechanically, but also thermally. The thermally induced R-B2O3-0 to β-B2O3-0 transformation can be easily detected in the variation of density or bulk modulus with temperature shown in Figure 5a. Accordingly, the thermally induced transformation between R- and β-B2O3-0 shows a significantly more abrupt density change than that resulting from a pressure release. Correspondingly, the bulk modulus exhibits a pronounced minimum at the transformation temperature, reflecting a shear mode instability oftentimes associated with displacive transformations. In our simulations, this transformation occurs at 1000 K and melting at 4000 K. These temperatures are considerably elevated compared to what one might expect to observe in experiments. Such an overestimation of the transition temperatures by simulations can be anticipated given the high-heating rates and the absence of defects or surfaces in simulated structures. An interesting point to notice in Figure 5a is that the bulk modulus of β-B2O3-0 increases with temperature and finally becomes higher than that of R-B2O3-0, even though the R-phase has a higher density. It is the more symmetric and puffed-up ring structure in the β-phase, as can be seen in Figure 5b, that gives rise to the higher rigidity of the structure compared to the R-phase, which has lower symmetry, partially twisted rings that appear puckered in the [010] projection. The small symmetric rings of the β-phase shown in Figure 5b are disordered due to thermal motion. However, these variations merely represent local configurations that are symmetrically equivalent upon rotation about the (001) axis. Consequently, the mechanical properties of the thermally induced β-phase are comparable to those of the stress-induced β-phase at low temperature, shown in Figure 4. Further compression of R-B2O3-0 past 1 GPa results in the amorphization of the structure in MD simulations. Figure 4 shows that the density of B2O3-0 is much closer to that of B2O3 glass (1.80 g/cm3)30 at ambient pressure than that of B2O3-I (2.56 g/cm3).1 Vitreous B2O3 can be thought to contain a mixture of puffed-up β-like and puckered R-like rings. Accordingly, compression converts puffed-up rings into puckered rings, causing the structure to compact and soften. Conversely, upon expansion the network rings puff up and the structure stiffens. The bulk modulus of B2O3 glass increases gradually as it expands under mechanical tension or thermal expansion due to the increasing population of β-like rings in the system.15 Figure 6 a shows that upon heating, B2O3-I transforms into B2O3-0 at 2800 K with a dramatic density drop before it finally melts at about 4000 K. The thermally introduced B2O3-0 at ambient pressure has an open structure, similar to that of the pressure-introduced one at ambient temperature, but is much more disordered (from the top view in Figure 6b). The total structure factor of B2O3-I is plotted in Figure 7a, consisting of many peaks, typical of a crystalline structure. After the transformation, the total structure factor of B2O3-0 at high temperature is very close to that of B2O3 liquid and that of a B2O3 glass with 63% of the boron atoms incorporated in boroxol rings (Figure 7b). The glass was prepared by removing Cs2O and some BO3 triplets from a cesium eneaborate crystal and annealing it for specific amount of time, as described in a previous publication.15 The discovery of the low-density B2O3-0 polymorphs provides a key to uncovering the mystery of the



Figure 5. (a) Change in density and bulk modulus with temperature when heating up R-B2O3-0 crystal; (b) the corresponding structural configurations in the stability ranges delineated in (a) and viewed in the [010] direction of B2O3-0.

B2O3 crystallization anomaly. Because the local structure of B2O3 liquid is closer to that of B2O3-0 crystals than that of B2O3-I, the former will probably nucleate first from the melts upon cooling. Further cooling, however, renders B2O3-0 unstable at ambient pressure. At the same time, the thermal energy needed for the structural rearrangements to convert to B2O3-I is withdrawn, and B2O3 is arrested in an amorphous state, perhaps containing residual fragments of B2O3-0. Only by applying pressure can the nucleation of B2O3-0 be avoided, which explains why a pressure is always needed for B2O3-I crystal to form from B2O3 liquid. If substantiated, the prediction of the crystalline B2O3-0 phase may have important ramifications for understanding the structure and properties of vitreous B2O3, particularly with respect to the existence of boroxol rings. Let us now turn to the high-pressure end of the range studied. Prewitt et al.2 obtained the known high-pressure phase, B2O3II, by quenching the compound from 6.5 GPa and 1200 °C. They determined the structure to be orthorhombic with the space group Ccm21 and cell dimensions of a ) 4.613, b ) 7.803, and c ) 4.129 Å. This structure consists of corner-linked BO4 tetrahedra. Within each tetrahedron, there are three long B-O distances of 1.507, 1.506, and 1.512 Å, and one short distance of 1.373 Å. As was already determined previously by Li et al.,13 bonding is between 60 and 70% ionic in both BO3 and BO4 units with marked covalent contributions through sp2 and sp3 hybridization, respectively. Upon conversion between BO3 and BO4 units, charge localization is little affected. Charge neutrality is maintained by virtue of the fact that two-thirds of the oxygen atoms are shared between three boron atoms and one-third between two. In our simulations, we isotropically compressed B2O3-I and obtained a new high-density B2O3 phase. In this

phase, B atoms are also four-coordinated by O atoms, but the tetrahedron is more asymmetric: the three long B-O distances measure 1.533, 1.510, and 1.480 Å, respectively, and the short distance measures 1.372 Å. In both phases, the oxygen associated with the short distance is coordinated by only two boron atoms, whereas the others are coordinated by three boron atoms. Thus, two-thirds of O atoms are three-coordinated, and onethird of O atoms are two-coordinated by B atoms. Even though the bonding is same, the space group of our high-density crystal is P3121 (see Table 2), that is, different from that of B2O3-II (Ccm21). The structures of these phases were compared by using CMPZ,31 a newly developed algorithm in KPLOT25 for the efficient comparison of periodic structures. CMPZ found no agreement between these two structures. We tentatively call our high-density phase B2O3-III. The X-ray diffraction patterns of all known and newly discovered crystalline phases are compared in Figure 8. In experiments, a combination of high pressure and high temperature is needed for synthesizing B2O3-II, while our concurrent MD simulations and first-principles calculations demonstrate that high pressure alone can result in the formation of B2O3-III by compressing B2O3-I at relatively low temperature. Once formed under high pressure, B2O3-III does not revert back to B2O3-I upon releasing pressure, and BO4 tetrahedra are stable at ambient conditions. This behavior is similar to that of B2O3II. High-pressure experiments at low temperature are needed to substantiate the existence of B2O3-III in the regime where the phase diagram of B2O3 has not been well studied. To understand why under compression B2O3-I transforms into B2O3-III instead of B2O3-II, we carried out first-principles calculations at 0 K to obtain the energy-volume relationships


Huang et al.

Figure 7. Total structure factor of B2O3-I crystal (upper pane) and of B2O3-0 crystal at high temperature compared to that of B2O3 liquid, resulting from melting B2O3-0, and glass containing 63% of boron atoms in boroxol rings (lower pane). Figure 6. (a) Change in density with temperature when heating up B2O3-I crystal; (b) the representative structures in the stability ranges delineated in (a), and the lower structure of B2O3-0 is the pressureinduced one at low temperature.

for these various phases. Employing the thermodynamic criterion of equal free energies allows one to identify the respective stability ranges for the different phases of B2O3. The crystalline B2O3 phases were equilibrated at several volumes and their energy-volume relationships were fit to the third-order BirchMurnaghan equation of state.32 The computed total energies as a function of volume are shown in Figure 9. Accordingly, the threefold coordinated B2O3-I structure is the most stable one at zero pressure. At their respective energy minima, the B2O3-0, B2O3-II, and B2O3-III phases are found to be about 16.7, 54, and 70.0 meV/atom higher in energy relative to B2O3-I, respectively. Remarkably, the energies of B2O3-II and B2O3III crystals essentially overlap after a certain amount of compression, that is, below about 7 Å3/atom. This behavior is compatible with a continuous phase transition between these structures, which is also clearly reflected in the enthalpy calculation as described below. Previous studies33-35 have shown that constant-pressure MD simulations based on the Parrinello-Rahman algorithm typically predicted transition pressures much higher than experimentally observed values. Hence, to obtain a better estimate for the equilibrium transition pressures between the different phases of B2O3, we calculated the static lattice enthalpy H ) Etot + PV at 0 K, where P ) dEtot/dV is obtained by direct differentiation of the fitted energy-volume curves. The computed enthalpy curve as a function of pressure of the B2O3 phases is plotted in

Figure 8. X-ray diffraction patterns of crystalline B2O3 phases.

Figure 10. The intersection of two enthalpy curves indicates the pressure at which the phase transition between these two phases is expected to occur. The enthalpy curves of the B2O3-I and B2O3-II phases intersect around 2.0 GPa. This predicted transition pressure is in excellent agreement with the experimentally established phase boundary between B2O3-I and B2O3-II3,5 and with the result of previous density functional calculations,12 keeping in mind that a direct transformation from B2O3-I to B2O3-II has not yet been observed in either experiments or simulations. Figure 10 also shows a possible phase transition between B2O3-I and B2O3-III at about 7.0 GPa, which is consistent with our observations in dynamical simulations


Figure 9. Calculated total energy as a function of volume for crystalline B2O3 phases. The energy of B2O3-I at zero pressure is taken as the origin of energy.

Figure 10. Enthalpy versus pressure curves of crystalline B2O3 phases. Inset is the magnification of the pressure range at which B2O3-I transforms into B2O3-0, B2O3-II, and B2O3-III.

using the Parrinello-Rahman method, in both classical MD and first-principles studies. The enthalpies of B2O3-II and B2O3-III phases, on the other hand, have practically the same values above 30.0 GPa, and it is impossible to distinguish which structure is more stable. This again would indicate a gradual phase transition between these two structures. In the tensile regime, our enthalpy calculation confirms a phase transformation from B2O3-I to B2O3-0 around -2.0 GPa. All these results strongly validate our constant pressure MD simulations at finite temperatures. According to our dynamical simulations, the phase transformation between B2O3-I and B2O3-II is kinetically impeded and a transformation from B2O3-I directly into B2O3-III is favored. On the basis of the locations in Figure 9 where common tangents would apply to the energy curves, the volume change associated with the transformation between B2O3-I and B2O3-III is somewhat smaller than that associated with the B2O3-I to B2O3II transformation. However, more likely the activation energy associated with the former transformation is smaller than that for the latter, probably because a more direct pathway connects the atomic positions of the beginning and ending configurations in the B2O3-I and B2O3-III transformation. Because either transformation involves a coordination change, the activation barriers are significant, and upon imparting sufficient elastic

J. Phys. Chem. C, Vol. 111, No. 37, 2007 13719 energy during the compression of B2O3-I chances are that the point of the common tangent to the energy curve of B2O3-III is reached. Fitting the calculated energy data with the Birch-Murnaghan equation of state also yields the bulk moduli and their pressure derivatives for these various phases. Our calculated bulk modulus of 56.0 GPa for B2O3-I is in good agreement with the experimental and theoretical result of 55 ( 15 GPa5 and 59 GPa,12 respectively. The pressure derivative of our calculated modulus is 3.02. For B2O3-II, our calculated bulk modulus is 177.1 GPa, which is an overestimation compared to the experimental value of 90 ( 15 GPa5 and to the theoretically predicted 144 GPa,12 but by being larger than that of B2O3-I, it follows the correct trend. For the B2O3-III phase, the bulk modulus and its pressure derivative from our calculations is 186.47 GPa and 4.11. The B2O3-0 phase has the smallest bulk modulus, 34.6 GPa, as expected. 4. Conclusions On the basis of a newly developed coordination-dependent charge-transfer potential, we studied the polymorphic transitions in crystalline B2O3 under various thermomechanical conditions by using concurrent MD simulations and first-principles calculations. On the basis of our simulations, we predict the existence of previously unknown low- and high-density B2O3 polymorphs. The low-density B2O3 phases (B2O3-0) not only provide twostate structural motifs for vitreous B2O3 that can explain its thermomechanical anomalies, but also shed new light on the mystery of the B2O3 crystallization anomaly. The predicted highdensity B2O3-III crystal is different from the known highpressure crystalline phase of B2O3-II, even though the bonding is the same in these two phases. B2O3-III is characterized by an initially higher energy than B2O3-II, but when compressed to higher densities, the energies of these two phases become indistinguishable. The transformation from B2O3-I to B2O3-III appears to be kinetically favored, especially at low temperatures. Further experiments are needed to verify the existence of these predicted phases. If substantiated, the phase diagram of B2O3 will be much richer than previously known. The existence of these phases will also have important ramifications for the understanding of the structure and properties of vitreous B2O3, especially under high pressure. Acknowledgment. The authors acknowledge the support of National Institute of Standards and Technology, under Grant 60NANB9D0102 and of the National Science Foundation under Grant DMR-0230662. Most of the work was done on supercomputers of the Center for Advanced Computing (CAC) at the University of Michigan. L.H. is grateful to Dr. Rudolf Hundt at the University of Bonn, Germany, for helpful discussion on determining the symmetries of the simulated structures and comparing two periodic structures. The atomic structures in this paper were rendered with SimReplay.36 References and Notes (1) Gurr, G. E.; Montgomery, P. W.; Knutson, C. D.; Gorres, B. T. Acta Crystallogr., Sect. B. 1970, B26, 906. (2) Prewitt, C. T.; Shannon, R. D. Acta Crystallogr., Sect. B 1968, B24, 869. (3) Mackenzie, J. D.; Claussen, W. F. J. Am. Ceram. Soc. 1961, 44, 79. (4) Uhlmann, D. R.; Hays, J. F.; Turnbull, D. Phys. Chem. Glasses 1967, 8, 1. (5) Brazhkin, V. V.; Katayama, Y.; Inamura, Y.; Kondrin, M. V.; Lyapin, A. G.; Popova, S. V.; Voloshin, R. N. JETP Lett. 2003, 78, 393.

13720 J. Phys. Chem. C, Vol. 111, No. 37, 2007 (6) Aziz, M. J.; Nygren, E.; Hays, J. F.; Turnbull, D. J. Appl. Phys. 1985, 57, 2233. (7) Turnbull, D. J. Non-Cryst. Solids 1985, 75, 197. (8) Takada, A.; Catlow, C. R. A.; Price, G. D. J. Phys.: Condens. Matter 1995, 7, 8659. (9) Takada, A.; Catlow, C. R. A.; Price, G. D. J. Phys.: Condens. Matter 1995, 7, 8693. (10) Takada, A.; Catlow, C. R. A.; Lin, J. S.; Price, G. D.; Lee, M. H.; Milman, V.; Payne, M. C. Phys. ReV. B 1995, 51, 1447. (11) Takada, A.; Catlow, C. R. A.; Price, G. D.; Hayward, C. L. Phys. Chem. Miner. 1997, 24, 423. (12) Engberg, U. Phys. ReV. B 1997, 55, 2824. (13) Li, D.; Ching, W. Y. Phys. ReV. B 1996, 54, 13616. (14) Huang, L. P.; Kieffer, J. J. Chem. Phys. 2003, 118, 1487. (15) Huang, L. P.; Kieffer, J. Phys. ReV. B 2006, 74, 224107. (16) Huang, L. P.; Durandurdu, M.; Kieffer, J. Nat. Mater. 2006, 5, 977. (17) Youngman, R. E.; Kieffer, J.; Bass, J. D.; Duffre`ne, L. J. NonCryst. Solids 1997, 222, 190. (18) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (19) Ordejo´n, P.; Artacho, E.; Soler, J. M. Phys. ReV. B 1996, 53, 10441. (20) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993.

Huang et al. (21) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182. (22) Wentzcovitch, R. M. Phys. ReV. B 1991, 44, 2358. (23) Cleveland, C. L. J. Chem. Phys. 1988, 89, 4987. (24) Gear, C. W. Numerical Initial Value Problems in Ordinary Differential Equations; Prentice Hall: New York, 1971. (25) Hundt, R.; Scho¨n, J. C.; Hannemann, A.; Jansen, M. J. Appl. Crystallogr. 1999, 32, 413. (26) Effenberger, H.; Lengauer, C. L.; Parthe, E. Monatsh. Chem. 2001, 132, 1515. (27) Islam, M. M.; Bredow, T.; Minot, C. Chem. Phys. Lett. 2006, 418, 565. (28) Cole, S. S.; Taylor, N. W. J. Am. Chem. Soc. 1935, 18, 55. (29) Huang, L. P.; Kieffer, J. Phys. ReV. Lett. 2005, 95, 215901. (30) Johnson, P. A. V.; Wright, A. C.; Singlair, R. N. J. Non-Cryst. Solids 1982, 50, 281. (31) Hundt, R.; Schon, J. C.; Jansen, M. J. Appl. Crystallogr. 2006, 39, 6. (32) Birch, F. J. Geophys. Res. 1986, 91, 4949. (33) Mizushima, K.; Yip, S.; Kaxiras, E. Phys. ReV. B 1994, 50, 14952. (34) Lee, I. H.; Jeong, J. W.; Chang, K. J. Phys. ReV. B 1997, 55, 5689. (35) Durandurdu, M.; Drabold, D. A. Phys. ReV. B 2001, 6401. (36) Shi, Y.; Huang, L. A visualization tool for computer simulations; www.simreplay.org.

New B2O3 Crystals Predicted from Concurrent Molecular Dynamics

Recommend Documents