Boroxol rings from diffraction data on vitreous boron trioxide

There has been a considerable debate about the nature of the short range atomic order in vitreous B2O3. Some authorities state that it is not possible to build a model of glassy boron oxide of the correct density containing a large number of six-membered rings which also fits experimental diffraction data, but recent computer simulations appear to overrule that view. To discover which view is correct I use empirical potential structure refinement (EPSR) on existing neutron and x-ray diffraction data to build two models of vitreous B2O3. One of these consists only of single boron and oxygen atoms arranged in a network to reproduce the diffraction data as closely as possible. This model has less than 10% of boron atoms in boroxol rings. The second model is made up of an equimolar mixture of B3O3 hexagonal ring ‘molecules’ and BO3 triangular molecules, with no free boron or oxygen atoms. This second model therefore has 75% of the boron atoms in boroxol rings. It is found that both models give closely similar diffraction patterns, suggesting that the diffraction data in this case are not sensitive to the number of boroxol rings present in the structure. This reinforces recent Raman, ab initio, and NMR claims that the percentage of boroxol rings in this material may be as high as 75%. The findings of this study probably explain why some interpretations based on different simulation techniques only find a small fraction of boroxol rings. The results also highlight the power of EPSR for the extraction of accurate atomistic representations of amorphous structures, provided adequate additional, non-scattering data (such as Raman and NMR in this case) are available.


Introduction
The nature of the short range order in glassy boron trioxide (B 2 O 3 ) has been debated for a long time right up until the present. In fact the answer to the question posed in the abstract above was answered, conceptually at least, in 1995 when Wright et al [1] proposed that a random network containing a large fraction of boroxol (B 3 O 6 ) rings could be built with the correct density provided one did this via interpenetrating networks, as opposed to the singly connected networks that had been built up until that time. Prior to this Hannon et al [2] showed, based on neutron diffraction data, that the short range structure of B 2 O 3 could be understood in terms of a high percentage of boron atoms in boroxol rings (f = 0.8) where f is the fraction of boron atoms which occur in boroxol rings. This in turn corroborated even earlier evidence from NMR [3] which also concluded the structure had a large fraction of boroxol rings, although the interpretation of the data in that case is difficult to follow in detail. Against these conclusions, however, a computer simulated random network model of vitreous B 2 O 3 by Verhoef et al [4] failed to find boroxol rings and assigned the 805 cm −1 Raman mode to a 'breathing mode' involving the three oxygens of a BO 3 triangle. It has to be said that the latter model did not reproduce the neutron scattering radial distribution function particularly well, which raises questions as to its validity. However a number of other diffraction and computer simulations of this material have seriously questioned the presence of a high degree of boroxol rings and a useful summary of this literature can be found in Huang and Kieffer [5].
The advent of reverse Monte Carlo (RMC) methods [6] drew another perspective on this problem by allowing three-dimensional large scale atomistic models of networks to be built which were deliberately constrained to reproduce the radiation (x-ray or neutron or both) diffraction patterns of the material in question. Prominent among these was the study by Swenson et al [7] where it was shown that models with significant numbers of boroxol rings (f = 0.5) gave less satisfactory representations of the diffraction results than models with lower numbers, and the conclusion was that the best fits were obtained with f < 0.2. The authors also conclude that if present in significant numbers then the boroxol rings would need to be heavily distorted to be consistent with the diffraction data. On the other hand Kohara and co-workers [8,9] built RMC models of vitreous B 2 O 3 based on x-ray diffraction data which apparently had a large fraction of boroxol rings, though the precise numbers appear not to be stated. At about the same time Zwanziger and co-workers [10,11], using a comparison of NMR data from the glass with that from crystals obtained values of f in the range 0.66-0.7.
In 2005 there ensued a controversy surrounding the ab initio molecular dynamics (MD) simulations of Umari and Pasquarello [12]. In this work the authors observed only a small fraction of boroxol rings (f = 0.09) in their simulation, but also noted the sharp mismatch between the simulated Raman spectrum compared to measurement. Using this mismatch they estimated the true fraction of boron atoms in boroxol rings in the real material to be closer to f = 0.74 and suggested the reason why such a high level of boroxol rings did not actually form in their simulation is because of the necessarily orders of magnitude higher quench rate in the simulation compared to experiment. However Swenson and Börjession [13] objected on the basis that their earlier work [7] had already established that a model with more than about 30% boroxol rings could not produce a satisfactory agreement with diffraction data, that the simulated Raman spectrum contained many other discrepancies with the Raman data other than the dominant ring 'breathing' mode at 808 cm −1 , and that an earlier analysis of Raman data had already concluded the boroxol percentage was as low as 10%. They also argued that the use of NMR chemical shift data to establish the extent of boroxol rings was unreliable, since this measures primarily the average of the B-O-B and O-B-O angles and, since these angles will occur in models with and without boroxol rings, this information should not be used to infer the presence of complete rings. In reply, [14], Umari and Pasquarello presented a further analysis of the NMR response, suggesting that a clear distinction between the in-ring B-O-B angle and the out-of-ring B-O-B angle is possible with this method, and so claiming, on the contrary, that NMR might be very sensitive to the boroxol ratio.
Another substantive contribution appeared in 2008 when Ferlat et al presented further ab initio computer simulation evidence in favour of a model with a high degree of boroxol content [15]. They produced two models. One, with a direct quench from the liquid, had a boroxol fraction of around 22%. The other based on a previous method using an empirical force field had a boroxol fraction around 75%. They claimed both models gave an adequate representation of the diffraction data.
More importantly their simulated Raman and NMR spectra, although still not producing an exact replica of the data, show significantly better agreement with experiment for the 75% boroxol ring model than for the 22% boroxol ring model. Note that, as for the Pasquarello et al work, the size of these simulation models was quite small, consisting of some 100 atoms for each model. Most recently Holland and co-workers have published new double-rotation (DOR) NMR data on this topic [16,17], and show that the separation of the ring B-O-B angle and the non-ring B-O-B angle is readily discernible by this technique. They conclude that the boroxol fraction is f = 0.73. These latter authors do caution of course that the results they present apply to the particular sample of vitreous B 2 O 3 they used, and that different results might have been obtained had the sample been prepared in a different manner.
Given these latter two, seemingly incontrovertible analyses, it seems difficult to maintain that vitreous B 2 O 3 can have anything other than a high degree of boroxol ring content. In that case however, why is it that some analyses, particularly those based on analysis of the diffraction data, came to the conclusion that the boroxol content is much lower? Why is there such a disparity between the different estimates? To try to answer this question, I tackle here once again the problem of building a structural model of amorphous B 2 O 3 which contains a high degree of boroxol ring content. The technique used, namely empirical potential structure refinement, EPSR [18,19], is analogous to, but distinct from RMC. Whereas RMC attempts to build the most disordered model that is compatible with a set of diffraction data and applied bonding constraints, EPSR attempts to achieve the same result by building in as much prior information into the problem as is possible. Thus atoms are moved around under the influence of forces, as in conventional atomistic computer simulation, and molecules are defined by harmonic forces, which are much more realistic than the bonding constraints that typically appear in conventional RMC. In order to fit the diffraction data the starting, or reference, interatomic potential is modified by a series of perturbations which are generated from the difference between the simulated and measured diffraction pattern. Equally, restrictions on local bonding are enforced by potential energy functions, not by bonding constraints. The following paper presents an account of this EPSR study of vitreous B 2 O 3 with a view to understanding why a reliable estimate of boroxol content from diffraction data has proved elusive.

EPSR analysis
The underlying methodology of EPSR is described in [18,19]. Two models of vitreous B 2 O 3 are developed. Model 1 consists of 800 individual B atoms and 1200 individual O atoms placed initially at random in a box of dimension 29.4880Å, giving an atomic number density of 0.078 atomsÅ −3 . The reference interaction potential energy between pairs of these atoms was defined by means of a Lennard-Jones plus Coulomb potential, with an additional repulsive exponential potential to prevent atomic overlap below a set of specified minimum distances: where the Lennard-Jones parameters are represented as combinations of the individual atom values and q α is the charge on atom α, R αβ is the minimum allowed separation of the respective atom pair, and w is a width value, set to 0.05Å in the present work. The coefficient, C αβ is determined from the requirement that no two atoms of the appropriate type get closer than the distance specified by R αβ . Hence the value of this coefficient could fluctuate up and down in the course of the simulation, depending on the degree of overlap, so that only typical values can be given here. The values of ( , σ ) are set to (0.5 kJ mol −1 , 0.42Å) for boron atoms and (0.6 kJ mol −1 , 2.7Å) for oxygen atoms respectively. The unit of electrical charge was set to 0.26e for the present work, so the charge on the boron atoms was +0.78e and that on the oxygen atoms −0.52e. For model 2 an equimolar mixture of hexagonal ring molecules of B 3 O 3 , with alternating B and O atoms, and BO 3 triangular molecules is built. To distinguish between ring and triangle atoms, the atoms within the boroxol ring are labelled BR and OR respectively, while the atoms in the triangular molecules are labelled BT and OT respectively. The intramolecular BR-OR and BT-OT distances are both set to 1.36Å and the BR-OR-BR, OR-BR-OR, and OT-BT-OT included angles are all set to 120 • . For the ring molecules an invisible 'Q' atom is set at the centre of the hexagon to stabilize the structure with all the other atoms of the ring equidistant from it. This atom has no charge or Lennard-Jones parameters, but it has the same atomic mass as the boron atoms. Thus it only contributes to the intramolecular energy of the system. Equally it has no x-ray or neutron scattering properties and so cannot contribute to the x-ray or neutron diffraction pattern. The intramolecular potential energy is defined by a series of harmonic potential energy functions: where the width parameter r ij is the average separation of each pair of atoms in the molecule and M i is the atomic mass of atom i. The function (5) is an effective Debye-Waller factor for each atom pair since the latter is not known a priori. (Computer simulations such as those described in [12,15] might be able to shed light on the correct values of these Debye-Waller factors.) In the present case the overall intramolecular energy amplitude, E intra , is set to 150 kJ mol −1 to give sharp intramolecular peaks in the radial distribution function consistent with the diffraction data. All the parameters for the intermolecular potential are exactly the same as those described for Model 1. Initial simulations were run using the reference potentials (including of course the intramolecular potentials for Model 2) alone at a temperature T = 1000 K. Once equilibrium was established, the diffraction data were used to generate an empirical potential to perturb these reference potentials in the manner described in [18]. Neutron diffraction data were obtained from [2]. X-ray diffraction data were obtained from [9] and then renormalised to the single atom x-ray scattering as described in [20]. The purpose of this renormalization is that when x-ray scattering are put on an absolute scale of differential scattering cross section it is necessary to refer to the single atom scattering line. This line corresponds to what would be observed if there were no atom correlations in the material. Hence normalizing to this line in such circumstances would give a constant as a function Q which is what would be expected for a system with no correlations. As explained in [20], to normalize to the square of the average form factor as is traditionally done introduces a Q dependent bias that can prevent EPSR from obtaining the optimum fit to those data. This of course does not preclude the use of the square of the average form factor normalization for other purposes. The initial refinement was run for some 2500 iterations of the empirical potential (approximately 5 × 10 7 individual atom moves for Model 1, 1 × 10 7 molecule moves for Model 2) until the fit, energy and pressure of the simulation were stable. Then, in order to sharpen up the peaks in the real space distribution functions, the temperature was lowered in steps, allowing equilibration at each step, to reach a final temperature of 25 K, to simulate the quench that might occur when making the real material. Note that this final temperature is much lower than the actual temperature of the diffraction measurements (∼300 K), but it should be borne in mind that in a Monte Carlo simulation the absolute temperature is simply a scaling factor on the interatomic potential energy, so the same results could have been achieved by scaling the simulation temperature back to 300 K and the potential energy functions ( s, qs, and empirical potential coefficients) by an equivalent factor. Inspection of Model 1 at the end of the refinement suggested only about 6% of the boron atoms in this model were in B 3 O 3 rings, whereas for Model 2 this percentage was by definition 75%.
The neutron and x-ray structure functions are each related to the site-site radial distribution functions via a Fourier transform, weighted according to the corresponding neutron scattering lengths or x-ray form factors for each atom pair. Hence for neutrons the interference structure function is defined as where c α , b α are the atomic fraction and average scattering length of component α, the angle brackets represent nuclear spin and isotope averages, and the site-site partial structure factors, H αβ (Q) are Fourier transforms of the corresponding site-site radial distribution functions: with ρ the atomic number density of the material. As stated above, for x-rays the data were renormalized to the single atom scattering to remove that part of the Q dependence that derives from the Q-dependence of the form factors. Hence in this work the x-ray structure function is defined as: . (8) The Fourier inversion of these functions, the total radial distribution function, defined here as gives the radial dependence of the average correlation function, with the individual site-site terms weighted according to the product of scattering lengths as in equations (6) and (8). Based on (6), for neutrons the value of the constant C is but for x-rays, given the different Q dependences of the form factors for different elements, it does not have a well defined value. However if we assume f α (Q) ∝ Z α , where Z α is the atomic number of component α, then Assuming the 11 B isotope is used for the neutron scattering data on B 2 O 3 as per [2], then C n = 0.3773 barns/sr/atom, while for x-rays C x ≈ 0.9554, which is dimensionless. Table 1 lists the x-ray and neutron weights outside individual site-site terms in the total structure and radial distribution functions appropriate to the present work. It can be seen that for both types of radiation B-B correlations are relatively weakly weighted in the total functions. For x-rays the O-O are the more strongly weighted than B-O, while for neutrons the B-O correlation is most strongly weighted. Figure 1 shows the simulated structure functions for both models and compares them with the respective diffraction data. The mean square deviation per measured data point was 0.000 31 for Model 1 and 0.000 73 for Model 2. Although the fits are not perfect, they are generally very close to the data, and to each other. The fact that Model 1 gives a slightly better fit than Model 2 is unsurprising given that Model 1, with its individual atoms free to move throughout the simulation box, has more degrees of freedom than are available to Model 2. On the basis of this investigation one would therefore have difficulty concluding that the model with few if any boroxol rings, Model 1, was significantly superior to that where there are a substantial number of rings, Model 2. There are greater discrepancies between x-ray data and fit for Model 2 in the Q range 0-6Å −1 compared to Model 1, but these may not be large enough to be significant, especially as the neutron fit for Model 2 appears slightly better than Model 1 in this same region.

Results
The total radial distribution functions (rdfs) derived from both data and simulation are shown in figure 2. (For comparison the structure functions of simulation and data have both been Fourier transformed over the same measured Q range of the relevant diffraction data.) These reveal a slight misfit in r-space around the position of the cross-ring B-O distance in the boroxol ring, if it exists. For a perfectly planar hexagonal ring structure this correlation should occur at twice the nearest neighbour B-O distance, i.e. 2.72Å, and indeed in the case of Model 2 there is a weak peak on the large r side of the main O-O peak at 2.37Å, as expected. Such a distinct peak is not obviously present in the Fourier transform of either the x-ray or neutron data. The discrepancy is more obvious in the x-ray fit for Model 2 compared to the neutron fit, where an intra-ring BR-OR correlation is specified from the definition of the boroxol ring at this distance. The fact that The lower dataset is the neutron diffraction data from [2], while the upper dataset (shifted upwards for clarity) is from [9]. the peak does not show up in the Fourier transforms of either dataset, which were measured completely independently of each other using different radiations, nor in any previous x-ray or neutron datasets, is a curiosity that goes beyond simply being an experimental artefact. Assuming this peak would be caused by B-O cross-ring correlations if it existed then the boroxol rings would need to be somewhat distorted to avoid this peak becoming visible in the total rdfs, as is found in the experiment.
In fact the small peak in the simulated total rdfs near 2.7Å arises not only because of an intra-ring B-O peak at 2.72Å as anticipated above, but also from an unexpected splitting or broadening of the B-B peak. Figure 3   in a structure dominated by boroxol rings. According to the recent NMR data [16] there should be no bond angles at this large value, which according to that data should be ∼135 • for non-ring boron atoms. Hence neither of these models reproduces a small but important and independent experimental observation.
To correct this feature and to give greater emphasis to the feature at ∼3.6Å in the total rdf which does not appear as pronounced in Model 2 as in the data or in Model 1, a revised version of Model 2 was developed. Additional terms to the BR-BT and BR-OT reference potentials were introduced to give attractive interactions at particular distances of the form with the values of E g = −40 kJ mol −1 , r g = 2.37Å for BR-BT, and E g = −10 kJ mol −1 , r g = 2.63Å for BR-OT. The value of w g was 0.2Å in both cases. In addition for the intermolecular BR-BR, BR-OR, OR-OR, BT-BT, BT-OT and OT-OT correlations were given a new minimum separation of 2.8Å and the OR-BT minimum separation was set to 2.9Å. Other minimum separations were left unchanged. One important benefit of these revised parameters was that the OT atoms on neighbouring BO 3 molecules were forced nearer the plane of the B 3 O 3 molecule to which they were bonded, as anticipated from an earlier analysis [2]. This gives a new, rather strong OR-BT peak near 3.0Å, and improves the fit to the total rdf in this region. Figure 4 shows the fits for the revised Model 2, while figure 5 shows the new total and site-site rdfs. The overall mean square deviation per data point for this revised model is 0.000 65, which is marginally better than the original value of 0.000 73. It is likely that further variations of the intermolecular potential could improve the fit still further, particularly if refinement of individual atom positions within molecules were allowed, but this is not attempted here because the point of the present exercise is to demonstrate that a model with 75% of the boron atoms in boroxol rings can be built which is compatible with the diffraction data and other, independent (NMR) data. It will also be apparent in the revised Model 2 that the BR-BT peak has moved closer to the OR-OR peak at 2.37Å compared to the original Model 2, figure 3, which in turn gives rise to a BR-OT-BT average angle (not shown here) near 135 • as required by the NMR data.
To get an idea of the differences in the local order presented by the two models of the diffraction data, figure 6 shows fragments of the structure from Model 1 and the revised Model 2. Note that for Model 2 the hexagonal and trigonal molecules are shown slightly distorted: this is to encompass the zero-point energy distortions that would occur in any real molecules and which are denoted by the finite width of the peaks in real space. However their average structure is planar and hexagonal or trigonal as specified. For Model 1 it is to be noted that the planar BO 3 molecules form readily in these simulations, while the hexagonal B 3 O 3 units are rare in this case. For Model 2 however 75% of the boron atoms are in B 3 O 3 hexagonal units.

Discussion and conclusion
It is worth mentioning that another model consisting of a mixture of planar B 3 O 6 molecules and free boron atoms was attempted, but abandoned early on when it became clear the free boron atoms were bonding to the ring oxygens in a manner that is not expected. Other models (not shown here) consisting of a mixture of B 3 O 3 molecules and free B and O atoms in the appropriate numbers were also attempted and these gave very acceptable fits to the diffraction data, but did not give such a clear demarcation between B 3 O 3 hexagons and BO 3 triangles as the present Model 2. Equally the idea of developing a model involving interpenetrating networks as proposed in [1] was felt to be too complicated and unnecessary for the purposes of the present discussion, but that concept remains a possibility for future work.
Based on the present analysis and preceding work [15] it is certainly feasible to build large scale models of vitreous B 2 O 3 which have the correct density and contain a large fraction of boroxol rings. It must be stressed that none of the structural models of vitreous B 2 O 3 presented in this paper produce an exact fit to the diffraction data (although the fits shown here are not significantly better or worse than previous RMC fits to these data, [7,9]). It is perfectly possible that there are small inconsistencies between x-ray and neutron data since they are obtained on different samples, plus the x-ray atomic form factors are never precisely known, particularly in cases such as the present one where there may be strong localization of electrons around chemical bonds. Such discrepancies would prevent such a perfect fit being obtained. Moreover by further adjusting of the potential parameters, bond lengths and Debye-Waller factors, it is conceivable better fits could have been obtained. Equally it must be emphasized that neither model given here is necessarily unique-other parameterizations might have given yet further possible structures. The revised Model 2 is preferred here to Model 1 because it fits better with independent (non-diffraction data), but this does not exclude the possibility that an even better model could have been found. The purpose here was not necessarily to obtain an exact fit to the diffraction data. Instead the aim was to establish whether a range of structural models with different local topologies would produce structure functions that are sufficiently different from one another to reliably distinguish between them, given the likely uncertainties in the diffraction data and the techniques used to model them. The answer to this question based on the preceding analysis appears to be no: the diffraction data are certainly sensitive constraints on any models of the structure that are proposed, but by themselves, at least in this particular case, they do not determine the structure uniquely. The extent to which this is true more widely requires each material to be considered on its own merit, as there are probably no general rules, but in the case of B 2 O 3 it appears that a range of structures can be generated which are consistent with the diffraction data. This almost certainly arises from the planar nature of the rings and triangles that apparently occur in this material. Such rings and triangles would not normally be expected to form spontaneously in a three-dimensional random network subject to stochastic forces and so therefore do not form in a computer simulation unless specifically forced to do so as in the present Model 2.
Why therefore have there been such conflicting reports on the percentage of boroxol rings? Almost certainly the answer relates to the results of the preceding analysis. As seen above the force fields that are used to build a particular structural model will have a strong impact on the local structural topology that is formed. A model can be formed from individual atoms (such as Model 1) or with some form of initial local topology (such as the 'molecules' used here in Model 2). How the atoms are constrained by those forces, both between atoms and molecules and between atoms within a given molecule, can have a profound effect on the local structure that emerges. In the case of vitreous B 2 O 3 it seems that a range of models with different local topologies can be made which are adequately consistent with the diffraction data. This happens in this case most likely because the B-B correlation is weakly weighted in the both x-ray and neutron scattering data so this correlation is not strongly constrained by the data themselves, even though variations in the B-B correlation can have profoundly different consequences for the overall structure. All models built by empirical potential computer simulation, RMC modelling, or EPSR modelling come within a class of models which mostly avoid many-body interactions, unless such interactions are introduced explicitly from the outset, as here in Model 2 in the form of 'molecules', or via an explicit three-body interaction such as in [5]. In that situation the only way to address the number of boroxol rings in vitreous B 2 O 3 is via other, non-scattering, experimental techniques.
Finally the results presented in this paper demonstrate how the empirical potential structure refinement (EPSR) method can be used as a sensitive and accurate probe of structure in amorphous and disordered systems. In recent work [21] EXAFS data were included in the EPSR analysis to gauge the local coordination of a dilute species in aqueous solution, while in the present example non-diffraction structural information from previous NMR data has been included in the structure refinement.