Development of empirical potentials for sodium borosilicate glass systems

Journal of Non-Crystalline Solids 357 (2011) 3313–3321

Contents lists available at ScienceDirect

Journal of Non-Crystalline Solids j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j n o n c r y s o l

Development of empirical potentials for sodium borosilicate glass systems Le-Hai Kieu a,⁎, Jean-Marc Delaye a, Laurent Cormier b, Claude Stolz c a b c

Service d'Études et Comportement des Matériaux de Conditionnement, DEN/DTCD/SECM, CEA EA Marcoule, BP 17171, 30207 Bagnols sur Cèze, France Institut de Minéralogie et de Physique des Milieux Condensés (IMPMC), Université Paris 6, CNRS UMR7590, 140 rue de Lourmel, 75015 Paris, France Laboratoire de Mécanique des Solides, CNRS UMR7649, École Polytechnique, 91128 Palaiseau, France

a r t i c l e

i n f o

Article history: Received 16 February 2011 Received in revised form 18 May 2011 Available online 21 June 2011 Keywords: Glass; Alkali borosilicate; Empirical potentials; Elastic moduli

a b s t r a c t New parameter values are proposed for the empirical potentials used to describe SiO2–B2O3–Na2O alkali borosilicate glass systems. They are based on Buckingham potentials, but include dependence between the fitting parameters and the glass chemical composition to improve the representation of the complex environment around the boron atoms. In particular, the boron anomaly (observed when the [Na2O]/[B2O3] ratio varies) is correctly reproduced. The structural and mechanical properties of a wide range of glass compositions and of reedmergnerite crystals are correctly simulated: bond distances, mean angles, densities, elastic moduli. The deviations from the experimental values are small. © 2011 Elsevier B.V. All rights reserved.

1. Introduction Classical molecular dynamics simulation of alkali borosilicate glasses is complex because their structures vary in a nonlinear manner with the [Na2O]/[B2O3] ratio. As this ratio increases, these nonlinearities, known as the boron anomaly, result in an initial increase in certain properties (glass transition temperature, viscosity, etc.) related to the increasing boron coordination number, after which the values then decrease. Yun, Dell and Bray [1–3] proposed a model describing glass structural changes by two ratios, R = [Na2O]/[B2O3] and K = [SiO2]/[B2O3] to explain, based on interactions between sodium and boron atoms, why the boron coordination number increases to a maximum and then decreases. This model has been validated by nuclear magnetic resonance (NMR) experiments over a range of glass compositions with different R and K ratios. The model has recently been refined with the latest advances in nuclear magnetic resonance spectroscopy [4,5]. The difficulty of simulating SiO2–B2O3–Na2O systems is the coexistence of two types of environments for boron atoms. The boron environment is either tricoordinate or tetracoordinate, and the transition from one to the other is accompanied by elongation of the B\O bond from 1.37 Å to 1.47 Å and by a change in the electronic structure and therefore in the local charge. Classical molecular dynamics simulation of borosilicate glasses traditionally uses Born-Mayer-Huggins (BMH) pair potentials [6–8]. The earliest work used integer charges; in general the atomic

⁎ Corresponding author. E-mail address: [email protected] (L.-H. Kieu). 0022-3093/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2011.05.024

structures were in good agreement with experimental results, but the elastic moduli were far from the experimental values, with deviations of up to 70%–100% [9]. More recently, models have been fitted in which the parameters vary with the local atom coordination numbers (Takada-Catlow-Price model [10,11] modified by CormackPark [12], Huang-Kieffer model [13]), but they are relatively complex. To mitigate the current lack of empirical potentials for alkali borosilicate glasses over a wide range of compositions, we propose a new parameter configuration that is dependent on the glass chemical composition. It was necessary to take into account a dependency between the glass chemical composition and the fitting parameters to obtain a satisfactory description of the environments of boron atoms. The simplicity of this formalism will later allow these force fields to be used for modeling structures containing several tens or hundreds of thousands of atoms. This article is organized as follows. Section 2 deals with the simulation method and the potential-fitting procedure. The results obtained over the range of compositions considered are presented in Section 3, followed by a discussion and conclusions. 2. Simulation methods 2.1. Initial model The literature is rich in potentials configured for the simulation of silicate systems, for example the Vessal potential [14–18], the Matsui potential [19–21], or the Guillot–Sator (GS) potential [22]. None of these potentials take boron into account. Our objective is to reproduce the structures and elastic properties of the SiO2–Na2O–B2O3 system over a wide range of compositions. We started with the existing potentials for the basic SiO2–Na 2O

3314

L.-H. Kieu et al. / Journal of Non-Crystalline Solids 357 (2011) 3313–3321

Temperature (K)

Table 1 Parameter values of the Guillot-Sator potentials. Cij (eV·Å6)

ρij (Å)

Aij (eV)

Si–O Na–O O–O Charge

50 329.69 120 360.22 9027.03 qSi = 1.89; qO = − 0.945;

0.161 0.17 0.265 qNa = 0.4725

46.1395 0 85.0321

compound and attempted to supplement them with the force fields around the boron atoms. We chose to use the Guillot–Sator parameters because of their classical form and good overall performance to reproduce the structures and properties of solid glasses. The Guillot–Sator potentials propose a parameter set (Table 1) for the Si\O, Na\O and O\O bonds. These values have been validated for the thermodynamic properties of a set of geological silicate materials [22]. For the B–O interactions we opted for a Buckingham pair potential completed by the classical Coulomb term: !
 qi qj rij Cij ϕ rij = − 6 + Aij exp − rij ρij rij

In this equation, rij is the distance between 2 atoms i and j; qi and qj are the effective atom charges; Aij, ρij and Cij are the parameters describing repulsion and dispersion force between 2 atoms. The parameters of the first B–O potentials used [23] were renormalized with respect to the GS potentials to obtain a first estimate of the B–O force field. The initial parameter values were progressively improved to better reproduce the glass properties listed in Table 2. This table describes the compositions and some of the properties of the glass formulations considered. 2.2. Glass simulation The classical molecular dynamics code DL_POLY [27] was used to fabricate six glass formulations as follows (Fig. 1). First, a cubic cell containing 1000 randomly placed Si, B, O, and Na atoms was equilibrated for 100 000 time steps at 5000 K (1 step = 1 fs). The structure was then quenched in steps to room temperature (300 K). The temperature difference between two successive steps was 100 K, and the quenching rate was 5 × 10 12 K/s. After quenching, the resulting glass was relaxed at room temperature and pressure in the NPT ensemble for 20 000 time steps. Periodic conditions were systematically applied throughout the simulation. The pair potentials were applied with a cutoff radius of 11 Å. The full Ewald sum was used for the Coulomb interactions. A final relaxation of 5000 time steps in the NVE ensemble is applied. Table 2 Compositions, densities, Young's moduli, bulk moduli, and boron coordination numbers of six glasses and of reedmergnerite. ρ (g/cm3)

Glass composition (mol%)

a

SB SBN3a SBN10a SBN12 SBN14 SBN55 Reedb a b

SiO2

B2O3

Na2O

69.5 48 44.4 59.66 67.73 55.30 75

30.5 48.7 46.1 28.14 18.04 14.71 12.5

0 3.3 9.6 12.20 14.23 29.99 12.5

Reference [24]. References [25,26].

E (GPa)

NVT

5000

Kb (GPa)

CB

NVT

P=0

NVE

NPT

300

Time (ps) 100

0

1040

1060

1065

Fig. 1. Glass simulation diagram.

2.3. Potential fitting The parameters of the B–O potential were determined in two steps. A constant non-integer charge potential was first developed. A series of trial-and-error calculations yielded the parameters reproducing the structures and properties of the glasses as closely as possible. With this first version of the potentials, the structures and mechanical properties of SBN14 and SBN12 glasses are correctly reproduced. However, when the [Na2O]/[B2O3] ratio lies outside the range for SBN14 and SBN12 glass, problems appear in reproducing the structures. Specifically, with SB glass (not containing sodium atoms), experimentally the boron atoms are all at coordination number 3 whereas in the simulated structure about 20% of the boron atoms are found in a tetracoordinate environment. The boron coordination number in the simulated SBN55 glass is too high, and the nonlinearity of the boron coordination number with the R and K ratios is not correctly reproduced. To correct these deviations, a more complex potential was implemented in which the fitting parameters vary with the glass chemical composition. This new form of potential is based on the physics and chemistry of borosilicate glass. It has been demonstrated, especially by ab initio calculations, that the electronic structure around a boron atom varies with its coordination number. Analysis by the Mulliken method gives different charge values for [3]B and [4]B in BO3H3 and (BO4H4) − groups [28]. The use of charges dependent on the boron coordination number has already been tested in earlier potential models [12,13] but these potentials are too complex to model large systems. 1.75

1.7

qB'/qO'

Bond

K=6.09

1.65

K=4.67 1.6

K=3.76 K=2.12

1.55 2.042 2.069 2.181 2.37 2.45 2.54 2.78

34.31 35.65 45.63 71.8 82 69.4 110.4

23.74 23.47 28.17 42 45 48.6 68.7

3 3.07 3.21 3.43 3.73 3.62 4

K=1.25 1.5 0

1

2

3

R Fig. 2. Modeling the composition-dependence of the ratio between boron and oxygen charges (- - -): values calculated with the Yun and Bray model; (___): values fitted using Eq. (1).

L.-H. Kieu et al. / Journal of Non-Crystalline Solids 357 (2011) 3313–3321

determined by the following procedure. The mean coordination number of the boron atoms in several glass compositions was first calculated according to the rules determined by Yun, Dell and Bray. The relative tri- and tetracoordinate boron concentrations in the glass are designated fB3 and fB4. We also assumed that the charge ratio of boron atoms at coordination numbers 3 and 4 remains constant at qB4/qB3 = 1.14 [13,28], and that qB3/qO = − 1.5 and qB4/qO = − 1.71. The second assumption is that the ratio between the mean charges of the boron and oxygen atoms in any SiO2–B2O3–Na2O glass composition depends on the qB3/qO and qB4/qO ratios weighted by the relative concentrations of boron at coordination numbers 3 and 4.

Table 3 Parameters of model Si–O, Na–O, O–O, Si–Si, Si–B, B–B potentials. Bond

Aij (eV)

ρij (Å)

Cij (eV. Å6)

Si–O Na–O O–O Si–Si Si–B B–B

45 296.72 120 360.22 9027.03 834.40 337.70 121.10

0.161 0.17 0.265 0.29 0.29 0.35

46.1395 0 85.0321 0 0 0

In addition, the mean coordination number of the boron atoms depends on the glass composition. It is thus logical to introduce a dependence between the glass composition and the force field around a boron atom. We therefore use composition-dependent variable charges. In our model, the charge of the boron atom is adjusted depending on the R and K ratios. The analytical form of this dependence was

a

3315

  qB′    q′ 

O YB

    q  q  = fB3  B3  + fB4  B4  qO qO

b B-O

Si-O

SBN14

SBN14

SBN55 SBN12

g(r)

g(r)

SBN55 SBN12

SBN10

SBN10

SBN3

SBN3 SB

SB

1.45

1.55

1.75

1.65

1.3

1.4

r (Å)

1.5

1.6

r (Å)

c

d Na-O

B-B SBN14

SBN14

SBN55

g(r)

g(r)

SBN55

SBN12

SBN12

SBN10 SBN3

SBN10

SB

SBN3

2

2.3

2.6

2.9

2.3

2.4

2.5

r (Å)

2.8

2.9

3

3.1

f Si-B

SBN14

O-O

SBN14

SBN55

SBN55

SBN12

SBN12

g(r)

g(r)

2.7

r (Å)

e

SBN10

SBN10

SBN3

SBN3 SB

SB

2.4

2.6

2.6

2.8

r (Å)

3

3.2

2

2.25

2.5

2.75

r (Å)

Fig. 3. Radial distribution functions: (a) Si–O, (b) B–O, (c) Na–O, (d) B–B, (e) Si–B, (f) O–O.

3

3316

L.-H. Kieu et al. / Journal of Non-Crystalline Solids 357 (2011) 3313–3321

where q′B and q′O are the effective charges of the boron and oxygen atoms in the glass, i.e. the value taken into account to compute the Coulomb interactions. The Yun, Dell gives the fB3 and  and  Bray  model  qB3  qB4  fB4 values for each glass. Since the   and   ratios are fixed, it was q q

Consider a configuration containing NSi Si atoms, NO O atoms, NB B atoms, and NNa Na atoms. System neutrality requires:

then possible to plot the variation of q′B/q′O versus the R ratio. The resulting curves (Fig. 2) were fitted by a nonlinear relation of the following type:

The following values are used in the Guillot–Sator potential model: qSi = 1.89; qO = − 0.945; qB = 1.4175; qNa = 0.4725. The charge values that will be used in our model are different, but they continue to ensure overall system neutrality. They are designated by q′Si, q′O, q′Na, and are determined from the following equation:

O

!

5

qB′ = −qO′ C6 K + ∑ Ci R + C0 2

qSi ·NSi + qO ·NO + qB ·NB + qNa ·NNa = 0

O

i

ð1Þ

i=1

qi′ = qi −NB

The fit gives:

qB′ −qB NSi + NO + NNa

where i = Si; O; Na

ð2Þ

C0 = 1:49643; C1 = 0:29504; C2 = −0:2565; C3 = 0:08721; C4 = −0:01323; C5 = 0:00073; C6 = 0:00315ðif R N 0:55Þ; C6 = 0 ðif R≤0:55Þ: The parameter C6 is equal to 0 for R ≤ 0.55 because, according to the Yun and Bray's model, in this region, the Boron coordination is independent of the K parameter. In Fig. 2 the values from our model satisfactorily reproduce those calculated using the Yun and Bray model; the deviations are systematically less than 3%. The charges of the other atoms were calculated as follows:

a

This equation allows the charges of each ion to vary according to the charge of the boron atom, while maintaining overall neutrality. The set of Eqs. (1) and (2) has only one solution (assuming the oxygen atoms are negatively charged) and unambiguously determines the charge of each ion for a given composition. Having determined the ion charges, the other parameters of the pair potentials must be fitted to reproduce as accurately as possible the structures and mechanical properties of alkali borosilicate glass. At this stage we had to introduce additional flexibility in the B–O force field by varying the parameters of the repulsion term with the composition to reproduce the B–O distance satisfactorily. We

b Si-O-Si

O-Si-O SBN14 SBN55

Frequency

Frequency

SBN14 SBN55 SBN12

SBN12 SBN10

SBN10 SBN3

SBN3

SB

SB

80

100

120

140

100

120

Bond angle (deg)

140

160

Bond angle (deg)

c

d O-B-O

B-O-B

SBN14

SBN55

Frequency

Frequency

SBN14 SBN55 SBN12

SBN12

SBN10

SBN10 SBN3

SBN3 SB

80

SB

100

120

Bond angle (deg)

140

90

110

130

150

Bond angle (deg)

Fig. 4. Angular distributions: (a) O–Si–O, (b) Si–O–Si, (c) O–B–O, (d) B–O–B.

170

L.-H. Kieu et al. / Journal of Non-Crystalline Solids 357 (2011) 3313–3321

3317

Table 4 Structural properties of simulated glass formulations. Glass

dB–O (Å)

SB SBN3 SBN10 SBN12 SBN14 SBN55 Experiment

CB

Calc

Exp

Calc

Y. & B.

1.37 1.38 1.39 1.41 1.45 1.44

1.37 1.38 1.39 1.41 1.44 1.43

3.01 3.09 3.23 3.41 3.72 3.58

3.0 3.07 3.21 3.43 3.73 3.62

dNa–O (Å)

dSi–Si (Å)

Si–O–Si angle (°)

O–Si–O angle (°)

B–O–B angle (°)

O–B–O angle (°)

– 2.53 2.53 2.53 2.51 2.46 2.29–2.62

3.06 3.16 3.08 3.13 3.15 3.16 3.08

143.1 144.74 144.57 145.16 145.73 147.69 144–147

109.35 109.38 109.37 109.4 109.39 109.38 109.5–109.7

143.4 139.49 140.24 141.99 138.41 151.70 129.4–143

119.4 115.94 115.63 113.73 111.14 112.17 111.6–118.6

therefore presumed that the AB–O parameter of the B–O potential, which expresses the intensity of the repulsive force between boron and oxygen, varies with the R* ratio:   AB−O = f R The R* ratio is defined as follows: 8   K K > > Rb + 0:5 > < min R; 16 + 0:5 ; 4  R =   > K > > + 0:5; K + 2 0; R∈ : 4 R* thus depends on the R and K ratios. The limits introduced to define the R* parameter are taken from Yun and Bray's model. The B–O repulsive term is modulated according to the evolution of the Boron coordination. The relation between the AB–O parameter and the R* ratio was refined by trial and error, ultimately yielding the final model represented by the following equation: 5

ABO = ∑ ai R



i

i=1

+ a0

where a 0 = 180 390.53; a 1 = 47 166.67; a 2 = − 43 827.65; a 3 = 210 268.55; a4 = −52 520.42; a5 = − 139 041.69 Table 3 contains the adjustable parameter values for the other atomic pairs. Compared to the original potentials from [22], the Si–O repulsive term has been slightly modified, and new repulsive terms for the Si–Si, Si–B and B–B pairs were added. We can now compare the simulation results with the experimental findings. 3. Results 3.1. Structural properties The six glass formulations in Table 2 were used to validate the empirical potentials. The radial distribution functions and the angular distributions of the six glasses are shown in Figs. 3 and 4.

Table 5 Concentration of nonbridging oxygen atoms (NBO) and tricoordinate oxygen atoms in the simulated glasses. Glass

R

K

SBN3 SBN10 SBN12 SBN14 SBN55

0.07 0.21 0.43 0.79 2.04

0.99 0.96 2.12 3.75 3.76

% [3]O

% NBO Simulation 0.97 1.99 2.64 2.8 23.17

[3]

Y&B

0 1 23

0.64 2.48 1.81 1.98 0.18

O

In the local environment of silicon atoms, the Si–O firstneighbor distance measured by the position of the first peak of the radial distribution function (Fig. 3a) conserves the same value for all six glasses (dSi–O = 1.61 Å) consistent with the experimental value (dSi–O = 1.61 Å [29], 1.60 Å [30], 1.62 Å [31]). The silicon coordination number remains equal to 4 for all six glasses: all the silicon atoms form SiO4 tetrahedra with four neighboring oxygen atoms. The simulated Si–Si distance (Table 4) varies slightly but remains near the experimental value (dSi–Si = 3.08 Å [29]). The tetrahedral environment of the Si atoms affects the shape of the O–Si–O angular distributions. In Fig. 4a the three O–Si–O distributions are virtually identical with a maximum at about 109.4° corresponding to the intra-tetrahedral O–Si–O angle (experimentally this angle is between 109.5° and 109.7° [29,32]). The distribution for the Si–O–Si angle (Fig. 4b) shows a broader spread, with a maximum position consistent with the experimental results (144–147°) [32–34]. The mean Si–O–Si angle of the simulated glass ranges from 144.8° to 147.9° (Table 4). The local environment of boron is complex in borosilicate glass. Some boron atoms adopt a tetrahedral environment similar to the SiO4 group with an alkali ion in a charge-compensating role (a sodium ion in this case). The other boron atoms adopt tricoordinate triangular environments without a charge compensator. On the radial distribution function for the B–O pair, the first peak is divided into two subpeaks for three glasses: SBN14, SBN55, SBN12 (Fig. 3b). The first subpeak corresponds to tricoordinate boron atoms and the second to tetracoordinate boron atoms. The experimental value of the B–O distance depends on the local boron coordination number, and ranges from 1.37 Å (for [3]B) to 1.47 Å (for [4]B) [35]. For comparison with the simulation, we estimated the experimental distance using the formula: dB–O = % [3]B·1.37 + % [4]B·1.47. The boron coordination numbers were estimated using the Yun and Bray model. In the simulated SBN14 glass, the boron coordination number is 3.72; it therefore contains a majority of tetracoordinate boron atoms, which accounts for the fact that the second subpeak in the B–O radial distribution function is higher than the first. Conversely, in SBN12 glass, the simulated boron coordination number is 3.41, hence the first subpeak is more intense than the second. In SBN55 glass (CB = 3.58), the intensity of both subpeaks is of the same order, reflecting the small difference between the concentrations of tricoordinate and tetracoordinate boron atoms. The second subpeak decreases for the glass with low mean coordination numbers around the boron atoms (SBN10: CB = 3.23 and SBN3 CB = 3.09), and disappears in the case of SB glass (CB = 3.01). Table 6 Concentration of nonbridging oxygen atoms (NBO) and [42].

[3]

O in SiO2–B2O3–BaO glasses [3]

Glass

R = [BaO]/[B2O3]

K = [SiO2]/[B2O3]

%NOB

%

BBS252 BBS352 BBS433

0.46 0.58 1.33

0.46 0.5 1.00

26–30.1% 5.1–6.8% 4.1–4.8%

5.6% 4.0% 1.3%

O

3318

L.-H. Kieu et al. / Journal of Non-Crystalline Solids 357 (2011) 3313–3321 Table 8 Structural and mechanical properties of reedmergnerite crystal.

1.5

S(q)

MD Exp

Properties

1.0 0.5

S(q)

1.5

RMC Exp

1.0 0.5 0.0

0

2

4

6

8

10

12

14

16

q(Å-1) Fig. 5. Structure factor of experimental SBN14 glass (Exp), simulated by molecular dynamics (MD), and corrected by the “reverse Monte Carlo” (RMC) method.

In each case, the B–O distance increases with the boron coordination number (Table 4), which is logical since the oxygen atoms are farther away from a boron atom at coordination number 4 than at coordination number 3. The increase in the mean coordination number of the boron atoms results in shift of the first peak of the B–B and Si–B radial distribution functions (Fig. 3d and e). The B–B distance increases from 2.69 Å (SB glass) to 2.80 Å (SBN14 glass) and the Si–B distance from 2.80 Å to 2.91 Å. The local increase in the B–O distance results in a greater distance between second-neighbor cations. The O–B–O angular distributions in the six glasses (Fig. 4c) are characterized by a high-amplitude peak whose position depends on the boron coordination number. This can be explained logically by the mixture of triangular BO3 (with a mean O–B–O angle of 118.6° in the B2O3–I crystal [36]: all the boron atoms are tricoordinate) and tetrahedral BO4 (mean O–B–O angle = 111.6° in the B2O3-II crystal [37]: all the boron atoms are tetracoordinate). As for the Si–O–Si angle, the B–O–B angle exhibits a broader spread. The experimental value of the B–O–B angle is between 129.4° in the B2O3 crystal [36] and 143° in borate glass [38]. The mean angular values are indicated in Table 4. The environment around the sodium atoms is of a different nature because this ion is less bound to the polymerized network and may adopt either a charge-compensating or a network-modifying role. The Na–O distance thus depends on the role of sodium, ranging from 2.29 Å (network modifier) to 2.62 Å (charge compensator) [31,39]. In alkali borosilicate glass xA2O–B2O3–rSiO2 (0.01 b x b 2.0, r = 1 and 2, A = Na et K), Takahashi et al. [40] reported that sodium atoms are mainly charge compensators around BO4 groups when x b 0.55–0.6. When x exceeds this threshold the network-modifying sodium atoms become increasingly numerous, confirming the composition-dependence of the role of sodium atoms. In SBN3, SBN10, and SBN12 glass (x=0.07, 0.21, and 0.43), sodium atoms are mainly charge compensators and the Na–O distance is 2.53 Å. When x=0.78 (SBN14 glass), some sodium

SB SBN3 SBN10 SBN12 SBN14 SBN55

ρ (g/cm3)

Lattice angles α (°) β (°) γ (°) dSi–O (Å) dB–O (Å) dNa–O (Å) Angle O–B–O (°) Angle O–Si–O (°) Angle Si–O–Si (°) Angle Si–O–B (°) Densité (g/cm3) Kb (GPa) E (GPa)

Simulation

7.84 12.37 6.81

7.41 12.66 6.80

93.32 116.37 92.04 1.61 1.47 2.54 109.41 109.41 142.24 135.57 2.78 68.7 110.4

93.0 115.7 91.9 1.60 1.47 2.60 109.46 109.41 145.51 135.88 2.83 64.34 119.4

atoms become network modifiers and the Na–O distance diminishes slightly (dNa–O =2.51 Å). The network modifiers become the majority in SBN55 glass (x=2.04), shortening the Na–O distance to an even greater degree (dNa–O =2.46 Å). The O–O radial distribution functions show two subpeaks due to the presence of different local entities: tetrahedra (SiO4, BO4) and triangles (BO3). The first peak corresponds to the O–O distances in the BO3 and BO4 entities, and the second to the O–O distances in the SiO4 entities. In conclusion with respect to the structural properties, the new potential model performs well compared with what is routinely observed in the literature: the differences are systematically less than 5%. 3.2. Behavior of nonbridging and tricoordinate oxygen atoms With the composition-dependent potential, the number of nonbridging oxygen atoms (NBO) in SBN14 glass is equal to 2.8%. According to the Yun and Bray model, the number of nonbridging oxygens in SBN14 glass would be about 1%. The difference is thus small. Similarly, the Yun and Bray model accounts for the number of nonbridging oxygen atoms in SBN12 and SBN55 glass (Table 5). Note that the experimental coordination number of the boron atoms in SBN12 glass determined by NMR and XANES differs to a greater degree from the value given by the Yun and Bray model (Y&B: 3.43; NMR: 3.34; XANES: 3.37) [41]. According to Yun and Bray, all the sodium atoms in SBN12 glass are charge compensators. In the actual glass, however, some sodium atoms behave as network modifiers, which explains the lower experimental coordination number around the boron atoms. The NBO concentration in SBN12 glass is therefore

Table 9 Structural properties of glass formulations SBN12, SBN14 (original charges) and SBN12′, SBN14′ (composition-dependent charges).

Table 7 Mechanical properties of simulated glass formulations. Glass

Experiment

Lattice lengths a (Å) b (Å) c (Å)

Glass

Kb (GPa)

E (GPa)

Calc

Exp

Calc

Exp

Calc

Exp

1.98 2.10 2.24 2.373 2.457 2.57

2.04 2.07 2.18 2.37 2.45 2.54

19.0 20.2 24.9 33.9 40.6 40.5

23.7 23.5 28.2 42 45 48.6

30.0 35.9 45.6 65.4 76.0 57.7

34.3 35.7 45.6 71.8 82 69.4

CB

CSi

dB–O dSi–O (Å) (Å)

SBN12_exp 3.43 4.00 1.41 1.6– 1.62 SBN12 3.38 4.00 1.42 1.61 SBN12′ 3.41 4.00 1.41 1.61 SBN14_exp 3.73 4.00 1.44 1.6– 1.62 SBN14 3.62 4.00 1.46 1.61 SBN14′ 3.72 4.00 1.45 1.61

dNa–O Si–O–Si O–Si–O B–O–B O–B–O (Å) angle (°) angle (°) angle (°) angle (°) 2.3– 2.62 2.50 2.53 2.3– 2.62 2.48 2.51

144– 151 145.38 145.16 144– 151 146.30 145.73

109.5– 109.7 109.36 109.4 109.5– 109.7 109.39 109.39

120– 143 138.60 141.99 120– 143 134.11 138.41

111.6– 120.3 113.93 113.73 111.6– 120.3 111.76 111.14

L.-H. Kieu et al. / Journal of Non-Crystalline Solids 357 (2011) 3313–3321

not zero, as correctly reproduced by the potential model proposed here. Zhao et al. [42] studied SiO2–B2O3–BaO glass by MAS (magicangle spinning) and 3QMAS (triple quantum MAS) nuclear magnetic resonance to analyze the nonbridging and tricoordinate oxygen sites. They found that the tricoordinate oxygen concentration depends on the glass composition (Table 6) and on the R and K ratios. The potentials developed here simulate structures that also contain a small fraction of oxygen atoms at coordination number 3 consistent with the experimental results (Table 6). This agreement remains qualitative because we do not have the experimental resources to quantify the concentration of tricoordinate oxygen in the glass.

3.3. Structure factor To estimate the validity of the overall structure of the simulated glass, a neutron diffraction experiment was performed on SBN14 glass in the Léon Brillouin Laboratory (CEA, France) to obtain the structure factor. The simulated glass structure factor was determined using the Faber–Ziman formula [43]:

Sα β ðqÞ = 1 +

i 4πρ0 ∞ h ∫ r gα β (r) −1 sin (qr) dr q 0

The experimental curve corresponds to: S(q) = ∑ cα cβ bα bβ Sα β ( q) α;β

where cα = Nα/N, Nα is the atomic concentration of species α, N is the total number of atoms, bα and bβ are the neutron scattering lengths for elements α and β. The scattering length values are taken from the NIST database [44]. It can be seen in Fig. 5 that the structure factor of the simulated glass satisfactorily reproduces the experimental values, with a slightly larger discrepancy on the first peak. This peak may correspond to longer-range patterns that are difficult to reproduce using molecular simulation methods for reasons of size or computing time. The “reverse Monte Carlo” method was applied to resolve the discrepancies between the experimental and simulated structure factors. Fig. 5 shows that the atomic model can be corrected to eliminate the structure factor deviations. It is then possible to compare the “refined” structure with the initial atomic structure to detect changes in the radial distribution functions. The radial distribution functions before and after refinement do not reveal any major differences. Moreover, there is no major evolution concerning the local coordinations, the local angles, or the ring distribution. 3.4. Density and elastic moduli

where Sαβ(q) is the partial structure factor corresponding to the α-β pair, ρo is the density, q is the scattering vector, gαβ(r) is the radial distribution function of the α-β pair.

a

The mechanical properties of all the glasses are satisfactorily reproduced with the usual deviations (b20%) for this type of empirical potential (Table 7). The elastic moduli were calculated using GULP

b 50 45 40

3319

B-O SBN12'

Si-O SBN12'

B-O SBN12

Si-O SBN12

3.0 Na-O SBN12'

2.5 Na-O SBN12

35

2.0

g(r)

g(r)

30 25

1.5

20 1.0

15 10

0.5

5 0 1.3

1.4

1.5

1.6

0.0 2.08

1.7

2.58

r(Å)

c

3.08

3.58

r(Å)

d

60

3.0

50

B-O SBN14'

Si-O SBN14'

B-O SBN14

Si-O SBN14

Na-O SBN14'

2.5 Na-O SBN14

2.0

g(r)

g(r)

40 30

1.5

20

1.0

10

0.5

0 1.35

1.4

1.45

1.5

1.55

r(Å)

1.6

1.65

1.7

0.0 2.08

2.58

3.08

3.58

r(Å)

Fig. 6. Radial distribution functions: (a) Si–O, B–O in SBN12, SBN12′ glass; (b) Na–O in SBN12 and SBN12′ glass; (c) Si–O, B–O in SBN14, SBN14′ glass; (d) Na–O in SBN14 and SBN14′ glass.

L.-H. Kieu et al. / Journal of Non-Crystalline Solids 357 (2011) 3313–3321

a

b

8 7

Frequency

6 B-O-B SBN12'

Si-O-Si SBN12'

5

O-Si-O SBN12'

6

Si-O-Si SBN12

5

O-Si-O SBN12

4 3

Frequency

3320

O-B-O SBN12' B-O-B SBN12

4

O-B-O SBN12

3 2

2 1

1 0 85

0 95

95 105 115 125 135 145 155 165 175

105 115 125 135 145 155 165 175

Bond angle (deg)

c

Bond angle (deg)

d

6

6 Si-O-Si SBN14'

5

B-O-B SBN14'

5

O-Si-O SBN14'

O-B-O SBN14'

4

O-Si-O SBN14

3 2 1 0 85

B-O-B SBN14

Frequency

Frequency

Si-O-Si SBN14

4

O-B-O SBN14

3 2 1

95 105 115 125 135 145 155 165 175

0 90

100 110 120 130 140 150 160 170 180

Bond angle (deg)

Bond angle (deg)

Fig. 7. Angular distributions: (a) Si–O–Si and O–Si–O in SBN14 and SBN14′ glass, (b) B–O–B and O–B–O in SBN14 and SBN14′ glass, (c) Si–O–Si and O–Si–O in SBN14 and SBN14′ glass, (d) B–O–B and O–B–O in SBN14 and SBN14′ glass.

[45] in the minimum energy configuration. The method was first tested on silica glass prepared using the BKS potentials [46]. The bulk modulus, Kb (37.87 GPa), and Young's modulus, E (64.98 GPa), differed slightly from the experimental values (36.9 GPa and 72.3 GPa [47]). 3.5. Validating the model on a reedmergnerite crystal With the newly developed potential we calculated the macroscopic properties of the crystal lattice of reedmergnerite (NaBSi3O8 with only tetracoordinate boron atoms). The experimental lattice properties are indicated in Table 8 and compared with the simulated values. All the structural properties satisfactorily reproduce the experimental results [48]. The bulk modulus is 64.34 GPa, which is about 8–10% lower than the experimental value (68.7 GPa [25]). In the case of Young's modulus, since this crystal is anisotropic, it deforms differently along the different axes. When a uniform pressure is applied on the material, the deformation along the X-axis is 62% of the total deformation (ε1 / (ε1 + ε2 + ε3) = 0.62) [26]. Knowing the deformation of reedmergnerite under pressure [26], we can thus estimate the mean value of Young's modulus at about 119.4 GPa, which is near the calculated value. The simulation results also reproduce the anisotropy of the crystal with a lower value for Young's modulus along the X-axis (44.5) and higher values on the other two axes (146.4–165.1). The difference between the simulated and experimental density is 0.05 g/cm 3, or 2%. The potential proposed here thus yields satisfactory results for the reedmergnerite crystal. 4. Discussion The inability of potentials independent of the chemical composition to reproduce correctly the evolution of the SiO2–B2O3–Na2O glasses on a large set of compositions is interesting. In fact, the electronic structure,

and hence the charge, of the boron atoms changes with the local environment and evolves, like the coordination, in a non linear manner with the glass composition. The new empirical potentials presented herein attempt to reproduce this non linear evolution by introducing a chemical composition dependence. The necessity to take into account the modification of the Boron charge with the chemical composition is quite logical since abinitio calculations and experiments show large differences between the three coordinated and four coordinated Boron atoms [13,28,49]. The proposed empirical potentials are closer to the physical reality, even if the use of point charges remains an approximation. To analyze the impact of composition-dependent charges, SBN12 and SBN14 glasses with modified charges (designated SBN12′ and SBN14′) were compared with simulated SBN12 and SBN14 glasses having the initial charges (i.e. the values of the Guillot-Sator model). The structural properties of the four glasses are listed in Table 9. Except for the modified charges, the glasses were fabricated with the same potentials and under the same conditions. Concerning the environment around silicon atoms, a difference is apparent in the intensity of the first peak of the Si–O radial distribution function. The coordination numbers, distances, and angular distributions are practically identical (Figs. 6 and 7), but the environment around the boron atoms is more significantly modified. The broadening of the two subpeaks of the B–O radial distribution function for SBN12′

Table 10 Mechanical properties of glass formulations SBN12, SBN14 (original charges) and SBN12′, SBN14′ (composition-dependent charges).

ρ (g/cm3) Kb (GPa) E (GPa)

SBN12

SBN12′

SBN12_exp

SBN14

SBN14′

SBN14_exp

2.378 33.39 58.8

2.373 33.9 65.4

2.37 42 71.8

2.473 40.48 73.75

2.457 40.56 75.96

2.45 45 82

L.-H. Kieu et al. / Journal of Non-Crystalline Solids 357 (2011) 3313–3321

and SBN14′ glass amplifies the difference between the tricoordinate and tetracoordinate boron atoms. The tetracoordinate boron atoms are 1.47 Å (SBN14′) and 1.46 Å (SBN12′) from the oxygen atoms, compared with [3]B–O distances of 1.40 Å (SBN14′) and 1.39 Å (SBN12′). The mean distance is closer to the experimental value for the modified charges. Moreover, the use of composition-dependent charges avoids the formation of certain exotic local structures (especially 2-member rings with two [4]B and two [3]O). The peak with the highest intensity in the B–O–B angular distribution shows a higher order for the second-neighbor around the boron atom. In Fig. 6b an increase can be observed in the Na–O distance related to the larger number of charge-compensating sodium atoms in glass with composition-dependent charges. The mechanical properties of SBN12′ and SBN14′ glasses are also slightly better reproduced. A decrease in density and an increase in the elastic moduli can be observed when the composition-dependent potential is used (Table 10). The difference between the simulated and experimental values diminishes. This method introduces a correlation between the ionic charges and the composition. Thus the potential is more efficient for varying local environments. For instance, we can imagine fruitful endeavor would be to apply this method to aluminosilicate glasses because the Al ions can adopt different local coordination depending on the chemical composition [50]. In relation to the potential, some results concerning the structure modifications with compositions deserve to be underlined. In particular, the separation between threefold Boron atoms and fourfold Boron atoms appears clearly on Fig. 3b. The B–O distances around fourfold Boron atoms are larger. On the other hand, the decrease of the R ratio along the series SBN14, SBN12, SBN10, SBN3 and SB induces a shift to smaller values of the B–O distances around the threefold Boron atoms. It seems that the Na atoms additionally facilitate the formation of fourfold coordinated Boron atoms and modify the local environment of the threefold Boron atoms, enlarging the B–O distances. In relation to the oxygen environment, the simulated structures contain some threefold coordinated O atoms in small concentrations. This result is in agreement with previous numerical and experimental studies. The validity of threefold O atoms has been examined experimentally via NRM studies [51]. The simulated structures contain more non bridging O atoms than predicted by Yun and Bray's model when the concentration is low. This result is in agreement with recent NRM results. Moreover, recent XANES experiments on a SBN12 glass confirm also this fact. Experimentally, it has been shown that the concentration of threefold Boron atoms is equal to 63% while Yun and Bray's model predicts a value of 57% [41]. Na atoms not trapped by the B environments allow the formation of non bridging oxygen atoms. This process is qualitatively reproduced by the empirical potentials developed here. 5. Conclusion This study highlights the performance of composition-dependent empirical potentials for SiO2–B2O3–Na2O glasses. It was necessary to introduce composition-dependent potentials to correctly reproduce the boron anomaly, i.e. the nonlinear evolution of the boron coordination number with changes in composition. Both the charges and the repulsive component of the B–O potential become composition-dependent, allowing us to reproduce with good accuracy the structural and mechanical properties of a series of glass compositions covering a wide range of [Na2O]/[B2O3] ratios. For one of these glasses the overall structure was correctly reproduced, i.e. with a simulated structure factor near the actual experimental structure factor. More generally, the technique of introducing a dependency between the ion charges and the composition could improve the simulation of

3321

other compositions containing network formers which, like boron atoms, are capable of adopting several types of local environments. In addition, the simplicity of these potentials allows large systems to be simulated.

Acknowledgments The authors are grateful to the staff of the LAIN at the University of Montpellier 2 (France) for measuring the elastic moduli of SBN14, SBN12 and SBN55 glasses.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51]

Y.H. Yun, P.J. Bray, J. Non-Cryst. Solids 27 (1978) 363. W.J. Dell, P.J. Bray, S.Z. Xiao, J. Non-Cryst. Solids 58 (1983) 1. J. Zhong, P.J. Bray, J. Non-Cryst. Solids 111 (1989) 67. G. El-Damrawi, W. Müller-Warmuth, H. Doweidar, I.A. Gohar, Phys. Chem. Glasses 34 (1993) 52. L.-S. Du, J.F. Stebbins, J. Non-Cryst. Solids 351 (2005) 3508. B. Vessal, M. Amini, C.R.A. Catlow, J. Non-Cryst. Solids 159 (1993) 184. F. Gou, G.N. Greaves, W. Smith, R. Winter, J. Non-Cryst. Solids 293–295 (2001) 539. W. Soppe, H.W. den Hartog, J. Non-Crys. Solids 108 (1989) 260. N. Deladerriere, J.M. Delaye, F. Augereau, G. Despaux, S. Peuget, J. Nucl. Mater. 375 (2008) 120. A. Takada, C.R.A. Catlow, G.D. Price, J. Phys-Condens Mat. 7 (1995) 8659. A. Takada, C.R.A. Catlow, G.D. Price, J. Phys-Condens. Mat. 7 (1995) 8693. A.N. Cormack, B. Park, Phys. Chem. Glasses 41 (2000) 272. L. Huang, J. Kieffer, Phys. Rev. B 74 (2006) 224107. W. Smith, G.N. Greaves, J. Chem. Phys. 103 (1995) 3091. W. Smith, T.R. Forester, G.N. Greaves, S. Hayter, M.J. Gillan, J. Mater. Chem. 7 (1997) 331. B. Park, H. Li, L.R. Corrales, J. Non-Cryst. Solids 297 (2002) 220. G. Lusbardi, G. Malavasi, L. Menabue, M.C. Menziani, J. Phys. Chem. B 106 (2002) 9753. L. Linati, G. Lusvardi, G. Malavasi, L. Menabue, M.C. Menziani, P. Mustarelli, U. Segre, J. Phys. Chem. B 109 (2005) 4989. M. Matsui, Geophys. Res. Lett. 23 (1996) 395. M. Matsui, Phys. Chem. Miner. 23 (1996) 345. T. Taniguchi, S. Ito, Phys. Chem. Glasses 43 (2002) 493. B. Guillot, N. Sator, Geochim. Cosmochim. Acta 71 (2007) 1249. J.-M. Delaye, V. Louis-Achille, D. Ghaleb, J. Non-Cryst. Solids 210 (1997) 232. C.A. Maynen, G.A. Saunders, S. Scholes, J. Non-Cryst. Solids 12 (1973) 271. T.P. Hackwell, R.J. Angel, Eur. J. Mineral. 4 (1992) 1221. I. Parsons, Feldspars and their actions, Kluwer Academic Publishers, 1994. W. Smith, The DL_POLY Molecular Simulation Package (http://www.cse.scitech. ac.uk/ccg/ software/DL_POLY/). V.V. Maslyuk, T. Bredow, H. Pfnur, Eur. Phys. J. B 41 (2004) 281. D.I. Grimley, A.C. Wright, R.N. Sinclair, J. Non-Cryst. Solids 119 (1990) 49. V. Petkov, S.J.L. Billinge, S.D. Shastri, B. Himmel, Phys. Rev. Lett. 85 (2000) 3436. M. Fabian, P. Jovari, E. Svab, G. Meszaros, T. Proffen, E. Veress, J. Phys-Condens. Mat. 19 (2007) 335209. H.F. Poulsen, J. Neuefeind, H.B. Neumann, J.R. Schneider, M.D. Zeidler, J. Non-Cryst. Solids 188 (1995) 63. R.L. Mozzi, B.E. Warren, J. Appl. Crystallogr. 2 (1969) 164. F.L. Galeener, Philos. Mag. B 51 (1985) L1. R.L. Mozzi, B.E. Warren, J. Appl. Crystallogr. 3 (1970) 251. G.E. Gurr, P.W. Montgomery, C.D. Knutson, B.T. Gorres, Acta Cryst. B: Struct. Sci. 26 (1970) 906. C.T. Prewitt, R.D. Shannon, Acta Cryst. B: Struct. Sci. 24 (1968) 869. W. Soppe, C. van der Marel, W.F. van Gunsteren, H.W. den Hartog, J. Non-Cryst. Solids 103 (1988) 201. D.A. McKeown, G.A. Waychunas, G.E. Brown, J. Non-Cryst. Solids 74 (1985) 325. K. Takahashi, A. Osaka, R. Furuno, J. Non-Cryst. Solids 55 (1983) 15. G. Bureau, Ph. D. thesis, Université Paris VI (2008). P. Zhao, S. Kroeker, J.F. Stebbins, J. Non-Cryst. Solids 276 (2000) 122. X. Zhang, D.A. Drabold, Phys. Rev. B 62 (2000) 15695. Neutron News 3 (1992) 29 (http://www.ncnr.nist.gov/resources/n-lengths/). J.D. Gale, General Utility Lattice, Program (https://www.ivec.org/gulp/). B.W.H. Van Beest, G.J. Kramer, R.A. van Santen, Phys. Rev. Lett. 64 (1990) 1955. A. Polian, D. Vo-Thanh, P. Richet, Europhys. Lett. 57 (2002) 375–381. R.T. Downs, H. Yang, R.M. Hazen, L.W. Finger, C.T. Prewitt, Am. Mineral. 84 (1999) 333. C.H. Hsieh, H. Jain, A.C. Miller, E.I. Kamitsos, J. Non-Cryst. Solids 168 (1994) 247. J.F. Stebbins, S. Kroeker, S.K. Lee, T.J. Kiczenski, J. Non-Cryst. Solids 275 (2000). J.F. Stebbins, J.V. Oglesby, S. Kroeker, Am. Mineral. 86 (2001) 1307.