The total energy has contributions from the Coulomb interactions, the van der waals interactions, and the bonded hydroxyl interactions. The bonded terms include the bond stretch term and the angle bend energy term which is helpful in describing hydroxyl bend behavior of metal-hydroxyl. The van der waals (VDW) energy term is represented by Lennard– Jones 12– 6 (LJ12– 6) function

where the “ arithmetic mean” combining rule is used for the distance parameters R_{o, ij} and the “ geometric mean” combining rule for the energy parameters D_{o, ij}. Combination rules allow Lennard– Jones parameters for single atom species to be combined in one of several ways to avoid the explicit definition for every possible atom– atom interaction in a simulation

where R_{o, i}, R_{o, j} and D_{o, i}, D_{o, j} are the distance parameters and the energy parameters, respectively, which can be derived from the fitting of the observed structures.

The bond stretch term uses a harmonic function to describe hydroxyl group stretch energy

where k₁ is a force constant and r_o represents the equilibrium bond length, while r_ij is the distance between atoms i and j. There are no Coulomb interactions between atoms coupled by a bond stretch term, as it is assumed that this potential describes all components of the interactions between both atoms.

An angle bend term is included to improve the vibrational behavior of surface hydroxyl groups

where k₂ is a force constant and θ _o represents the equilibrium bend angle, while θ _ijk is the bend angle for metal– oxygen– hydrogen. This bend term is considered as an optional feature of CLAYFF and does not affect the total energy obviously, but it is computationally expensive to apply, ^[10] however, in the present paper, we still included it in modeling for more accuracy. Moreover, in order to describe the relationship of two splitting H sites in the OH-topaz well, a three-body function with the above harmonic form was also used, then θ _ijk presents the angle for hydroxyl hydrogen (site 1, H1)-hydroxyl oxygen-hydroxyl hydrogen (site 2, H2). All the potentials parameters are given in Table 1.

Table 1.

Table 1.

Table 1. Potentials parameters. LJ12– 6 potentials parameters and ions charges Atom type description i Do/eV Ro/Å Charge/e Hydroxyl hydrogen HOH 0 0 0.425 Hydroxyl oxygen OOH 6.7389 × 10− 3 3.5532 – 0.950 Oxygen bridging Si and Si or Si and Al OB 6.7389 × 10− 3 3.5532 – 1.050 Tetrahedral silicon Si 7.9814 × 10− 8 3.7064 2.100 Octahedral aluminum Al 5.7667 × 10− 8 4.7943 1.575 Fluorine substituting hydroxyl F 1.7000 × 10− 5 4.9321 – 0.525 Bond stretch potential parameters i j k1//eV· Å 2 r0/Å OOH HOH 48.1000 1.0 Angle bend potential parameters i j k k2/eV· Å 2 θ 0/(° ) Al OOH HOH 2.6020 109. 47 HOH OOH HOH 342.4925 109. 47

Table 1. Potentials parameters.

No parameters about F ions are included in the initial CLAYFF, so LJ12– 6 parameters for F ions were derived from the fitting of the observed structures. Aluminium trifluoride structure ^[18] was first used to fit for the initial F ions parameters, then F-topaz structure^[1] was supplemented to a further fitting, finally a relax fitting^[19] was needed to avoid the worse structure. When fitting the angle potential for H1-hydroxyl oxygen-H2 in OH-topaz, simple point charge (SPC) water model^[20] was used as initial parameters, then a conventional fitting and next a relax was processed. In order to keep the electrostatic equilibrium of the unit cell and just modify the initial force field slightly, the charge of F ions was assigned to be equal to hydroxyl, − 0.525e, where e is the charge of an electron. All fitting calculations were performed by means of the general utility lattice program (GULP), ^[19] and the Newton– Raphson method with Broyden– Fletcher– Goldfarb– Shanno (BFGS) updating of Hessian was used.^[21]

The cell parameters and atom coordinates for F- and OH-topaz ^[1] were used as the initial structures in the calculations. In OH-topaz, every hydroxyl hydrogen has two possible sites, H1 and H2, as shown in Fig. 1. Each H site is one-half occupied and the diagonal configuration of H atoms provides a geometrically reasonable model.^[22] Here, a mean field (or virtual crystal) approximation^[9] was adopted to describe this configuration and the space group of pbnm was still adopted.^[22] In order to model nature topazes with different F percentages, an attempt of the substitution of F for OH in OH-topaz was carried out, and the results would be used to compare with the experimental results.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Schematic diagram of H sites in OH-topaz viewed along the b axis.

Before each calculation, the structure optimization was firstly performed by means of the Newton– Raphson minimization method with BFGS updating of Hessian, ^[21] the accuracy factor for lattice sums was default to 8.0 and the max gradient tolerance was 0.001. The Coulomb energy summation used the standard Ewald method, the short VDW energy used a cutoff radius of 8.5 Å . For preserving a continuous transition of potential surface, a polynomial was used to taper potentials to zero between 8.0 and 8.5 Å .

The defect energy calculations were performed by the Mott– Littleton method, ^[23] which divides the crystal into inner and outer regions, with the former being modeled by an explicit atomistic approach and the latter by quasi-continuum methods.^{[24, 25]} The calculations employed an inner region size of 6 Å and an outer region size of 12 Å . When performing defect calculations, the Newton– Raphson optimizer with the BFGS hessian update and the maximum gradient tolerance of 0.001 was used to get the minimum energy for the final defect.

In order to investigate the role of the hydrogen bonds, bond valence theory^{[26, 27]} was used to estimate the strength of hydrogen bonds. The bond valence theory is based on the point that the bond length can describe the bond strength, it suggests that the longer the bond, the weaker the bond, and the smaller the bond valence; the bond valence sum of all bonds to one atom is equal to the valence of this atom. The relationship between the bond length and the bond valence was deducted from the accumulated experimental data^[25]

where S_ij and r_ij are the bond valence and the distance between atoms i and j, R_o and N are experiential parameters relating to the atom type and the valence which can be found in “ Accumulated table of bond valence parameters” .^[28] Here, R_o is 0.087 nm and N is 2.2 for the bond valence between hydroxyl hydrogen and hydroxyl oxygen.

The vibrational frequencies of a molecule can be calculated from the Hessian matrix which is produced in the process of energy minimum.^[29] In calculation, the Hessian matrix in Cartesian coordinates was firstly transformed to the matrix in mass– weight coordinates, and then the corresponding eigen-equation was solved for its eigenvalue λ _i. The normal frequency ν _i could be determined by

The calculated frequencies can be used to compare with the experimental results and test the force field.^{[2, 5, 17]} In the present work, the k point used in the frequency calculation was default to the gamma point, and the direction of approach to the gamma point for the nonanalytic correction to the dynamical matrix was default to (1, 1, 1).

Calculated perfect lattice properties for F- and OH-topaz are shown in Table 2. It can be seen that a good agreement with experimental values is gained and all lattice parameters reproduced by CLAYFF have the difference of less than 1%, which are better than Jackson’ s results.^{[2, 5, 17]} In addition, the H– H distance between two diagonally non-equivalent sites calculated here is 2.5 Å which is close to the experimental result.^[1] Table 2 also gives the lattice properties for structures with different F/OH ratios, Al₂(SiO₄)[F_0.70(OH)_0.30]₂^[30] (here named as 0.70F-topaz), Al₂(SiO₄)[F_0.72(OH)_0.28]₂^[31] (0.72F-topaz), Al₂(SiO₄)[F_0.78(OH)_0.22]₂ (0.78F-topaz), and Al₂(SiO₄)[F_0.79(OH)_0.21]₂ (0.79F-topaz).^[32] Once it was noted that the F-topaz structure was used to derive LJ12– 6 parameters, then calculations for other complex topazes containing F ions might be used to test supplemented CLAYFF. All calculated lattice parameters agree well with the observed structures and have the difference of less than 1%. The small deviations further indicate that the improved CLAYFF can model topazes with different F/OH ratios well.

Table 2.

Table 2.

Table 2. The calculated lattice parameters for topazes with different F/OH ratios. Lattice parameters/Å Calculated Experimental Difference/% F-topaz Ref. [1] a 4.657 4.652 0.11 b 8.834 8.801 0.37 c 8.477 8.404 0.86 OH-topaz Ref. [1] a 4.767 4.720 0.98 b 8.954 8.920 0.37 c 8.494 8.418 0.90 0.70F-topaz Ref. [30] a 4.706 4.670 0.78 b 8.874 8.849 0.29 c 8.453 8.392 0.73 0.72F-topaz Ref. [31] a 4.704 4.665 0.83 b 8.871 8.838 0.37 c 8.454 8.398 0.66 0.78F-topaz Ref. [32] a 4.696 4.667 0.62 b 8.862 8.834 0.32 c 8.458 8.395 0.75 0.79F-topaz Ref. [32] a 4.695 4.657 0.82 b 8.860 8.838 0.25 c 8.458 8.426 0.38

Table 2. The calculated lattice parameters for topazes with different F/OH ratios.

With OH-topaz as initial structures, lattice parameters with different F contents were also calculated by changing the F/OH ratio, the polynomial fit curves are shown in Fig. 2. The parameters of polynomial regression of a, b, and c, and the linear regression parameters of b are given in Table 3. For comparison, the experimental values of 0.70F-topaz, ^[30] 0.72F-topaz, ^[31] 0.78F-topaz, and 0.79F-topaz, ^[32] F- and OH-topaz^[1] are also given in Fig. 2, it can be seen that the trends of calculated values agree with the experimental values, especially for b and c, although the values of a are a little higher than the experimental results in the middle range. In addition, although the polynomial regression seems more appropriate than linear regression, the values of b do not deviate from the simple linear correlations with F contents apparently.^[33] This can be confirmed by the mean deviation between experiment and linear fit of b values, 0.0034, which is close to that between experiment and polynomial fit, 0.0032. The quasi-linear nature supports that the F content in topazes can be conveniently measured by the optical method because of the strong correlation between b and the optic axial angle.^{[33, 34]}

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Lattice constants (a (a), b (b), and c (c)) of different F contents and their fit curves, and some experimental values for comparison.

Table 3.

Table 3.

Table 3. Polynomial regression equations of lattice parameters of different F contents (standard errors for the regression coefficients are given in brackets). Polynomial regression equations R-square Standard deviation of the fit a = 4.76438 − 0.04258x − 0.06138x2 0.99765 0.0019     [0.00145] [0.00673] [0.00648] b = 8.9545 − 0.08976x − 0.0347x2 0.99827 0.00192     [0.00147] [0.00682] [0.00657] c = 8.49317 − 0.15831x + 0.14295x2 0.99715 8.47541 × 10− 4     [6.45702 × 10− 4] [0.003] [0.00289] Linear regression equation R-square Standard deviation of the fit b = 8.95971 − 0.12446x – 0.99612 0.00384     [0.00217] [0.00366]

Table 3. Polynomial regression equations of lattice parameters of different F contents (standard errors for the regression coefficients are given in brackets).

In order to investigate the role of the hydrogen bonds, bond valence theory^[27] was used to estimate the strength of hydrogen bonds. Here, the sum of hydrogen bond interaction was considered and derived by subtracting OH bond valence from the H atom bond valence sum. It is noted that the hydrogen bonds include not only H· · · O interaction but also H· · · F interaction. The curves of the bond valence sum of the hydrogen bond around each H site are given in Fig. 3. As can be seen, with the increase of the F content, both bond valence sums of hydrogen bonds gradually increase, this means that the more F there is in hydroxyl of the structure, the stronger the hydrogen bonds. This may be the reason why there is a uniform contraction of cation polyhedra in F-topaz relative to OH-topaz.^[1] In is also shown in Fig. 3 that the hydrogen bond sum around the H1 site is higher than the H2 site, which is consistent with the results^{[1, 35]} and suggests the different circumstances around them.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Bond valence sum of hydrogen bond in topazes with different F contents.

The formation energies for the basic intrinsic defects in F and OH-topaz are given in Table 4 where the Krö ger– Vink defect notation^[36] is used. All defect optimizations converge well, an exception is for Al_i in OH-topaz, but this could be overcome by using a unit diagonal matrix as the starting point for Hessian when the exact Hessian is ill-conditional; ^[25] in Jackson’ s results, ^{[2, 5, 17]} the calculation for Si_i did not converge and the corresponding Si_i and Si Frenkel formation energies were absent. For both F-topaz and OH-topaz, cation vacancies have the highest formation energies which indicate that V_Al and V_Si defects are difficult to form. Al_i, Si_i, and Si_Al in both topazes and O_F in F-topaz have the negative formation energies which means that these defects are much easier to form and more stable in the lattice than the other basic degects. The formation energies of above basic defects are defined as the energies calculated from the difference between the pure and the defective lattices. These energies were then used in the calculation of the Frenkel and Schottky energies, which are reported in Table 5. Here, the interstitial site is (0.89723, 0.3691, 0.75) for F-topaz, (0.90189, 0.3679, 0.75) for OH-topaz and the hydrogen site of the interstitial hydroxyl is (0.69839, 0.40087, 0.75), as shown in Fig. 4. For F-topaz, the exchange between O and F, and the exchange between Al and Si are the lowest energy processes, which will be favored over Schottky defect formation. For OH-topaz the lowest energy processes is exchange between Al and Si, which is favored over Schottky defect formation and OH Frenkel formation, this contradicts Jackson’ s results, ^{[2, 17]} in which OH Frenkel formation is the lowest energy process. It can be reasonably explained: the bigger volume of OH than F will bring more lattice distortion, therefore, more forming energy is required. Al– Si exchange is supported by the observation that the spectra of the luminescence is due to AlO₄, ^{[37– 39]} a centre formed when Al substitutes for Si and an extra hole is captured to neutralize the defect.

Table 4.

Table 4.

Table 4. Basic defects formation energy. Defect and energy symbol Formation energy in F-topaz/(eV/atom) Formation energy in OH-topaz/(eV/atom) VAl (EVAl) 20.837 19.222 VSi (EVSi) 40.505 38.379 VO1 (EVO1) 8.682 8.811 VO2 (EVO2) 8.402 8.663 VO3 (EVO3) 8.646 8.659 VF (EVF) 1.539 – VOH (EVOH) – 1.527 Ali (EAli) – 10.890 – 8.757a Sii (ESii) – 21.116 – 22.228 Oi (EOi) 1.956 1.556 Fi (EFi) 2.657 – OH1i (EOH1i) – 4.788 OH2i (EOH2i) – 4.398 AlSi (EAlSi) 15.761 15.025 SiAl (ESiAl) – 12.118 – 12.147 OF (EOF) – 3.360 – FO1 (EFO1) 5.776 – FO2 (EFO2) 5.595 – FO3 (EFO3) 5.720 – OOH (EOOH) – – 3.258 a Minimiser to switch to BFGS with unit hessian when gradient norm is less than 0.2.

Table 4. Basic defects formation energy.

Table 5.

Table 5.

Table 5. Calculated Frenkel, Schottky and ion exchange energies (U1 and U2 are the lattice energies of F-topaz and OH-topaz per unit chemical formula). Defect and energy symbol Energy expression Value for F-topaz/(eV/atom) Value for OH-topaz/(eV/atom) Al Frenkel (E1) E1 = EVAl + EAli 9.947 10.465 Si Frenkel (E2) E2 = EVSi + ESii 19.389 16.151 O1 Frenkel (E3) E3 = EVO1 + EOi 6.726 7.255 O2 Frenkel (E4) E4 = EVO2 + EOi 6.446 7.107 O3 Frenkel (E5) E5 = EVO3 + EOi 6.690 7.103 F Frenkel (E6) E6 = EVF + EFi 4.196 – OH1 Frenkel (E7) E7 = EVOH + EOH1i – 6.315 OH2 Frenkel (E8) E8 = EVOH + EOH2i – 5.925 Per-defect Schottky (E9) E9 = (2EVAl + EVSi + EVO1 + EVO2    + 2EVO2 + EVF + U1)/9 (F − topaz) 3.825 3.199 E9 = (2EVAl + EVSi + EVO1 + EVO2    + 2EVO2 + EVOH + U2)/9 (OH − topaz) Al– Si exchange (E10) E10 = EAlSi + ESiAl 3.643 2.878 O1– F exchange (E11) E11 = EOF + EFO1 2.416 – O2– F exchange (E12) E12 = EOF + EFO2 2.235 – O3– F exchange (E13) E13 = EOF + EFO3 2.360 –

Table 5. Calculated Frenkel, Schottky and ion exchange energies (U1 and U2 are the lattice energies of F-topaz and OH-topaz per unit chemical formula).

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. The considered interstitial site of hydroxyl group in OH-topaz, interstitial oxygen site is equal to that of hydroxyl oxygen. (a) viewed along c axis; (b) viewed along b axis. Hydroxyl group occurs in the interspaces formed by four octahedrons and a tetrahedron.

In this paper, the derivation of F potential parameters supplemented in CLAYFF was reported. With supplemented CLAYFF, the lattice properties of F-topaz, OH-topaz, and some other complex topazes were calculated and the results agree with the available experiment data well. The calculations for topaz with different F contents indicate that b does not deviate from the simple linear correlations to fluorine content obviously, which supports that the content of F in topaz may be measured by an optical method. The role of the hydrogen bond was also discussed in this paper; the results reveal that the more F in hydroxyl of the structure, the stronger the hydrogen bonds; the hydrogen bond sum around the H1 site is higher than that around the H2 site. Defect calculations provide the formation energies of all common defects and can be used to judge the ease of formation of them. The frequency calculations reproduce experimental IR and Raman vibration well, and the 1080-cm^{− 1} frequency occurring in natural OH-topaz samples can be attributed to Si– F stretching caused by the F substitution to OH and the Al– Si exchange.