Crystal structure prediction by ionic network analysis: the example of (p–T–X)-structure relationships in olivines

The method of Ionic Network Analysis (INA) is defined by reference to the known crystal structures of olivine minerals. It is based on a reversible transformation between two alternative representations of ionic crystal structures: (a) the crystallographic and (b) the interactional. Whereas the former encompasses unit-cell parameters and atomic coordinates, the latter consists of selected interaction vectors between ions. Since the lengths and orientations of these vary only slightly between crystal structures obtained under systematically varying (p, T, X) conditions, they may be used to predict the crystal structures at intermediate (p, T, X) values by interpolation. Two interactional networks are constructed, one for the anions and the other for cations. As both networks lead to independent calculated values of the unit-cell parameters, it is possible to exploit the known, continuous (p, T, X) variations of cell parameters as normative constraints for the prediction of atomic coordinates within a predictive structural refinement procedure. Continuously varying structurally based parameters such as the volumes of cation coordination polyhedra may likewise be used. The choice of olivines for developing the method has been guided by the availability of pressure, temperature and compositional structural data for them. However, the ideas are expounded sufficiently generally for the method to be applied to other minerals.

1. Introduction

Olivines continue to attract widespread attention from a variety of standpoints. From a geophysical perspective, the seismic discontinuity at a depth of 400 km within the Upper Mantle has long been attributed to a phase transition of (Mg,Fe)₂SiO₄ olivine to α- and β-spinels (Ringwood & Major, 1970). Olivines also serve as a model mineralogical system for investigating cation-ordering phenomena, since their two octahedral sites, termed M1 and M2, are divergent, i.e. symmetrically inequivalent. This has been the focus of many studies, for example, McCormick et al. (1987), Müller-Sommer et al. (1997), Redfern et al. (2000), Henderson et al. (2001), Heuer (2001), Rinaldi et al. (2005) and Heinemann et al. (2006).

As a contribution to the crystal chemistry of olivines, Lumpkin & Ribbe (1983) took the occupancies of M1 and M2 sites into account in calculating effective ionic radii for composite M1 and M2 cations. These were subsequently employed as the independent variables in a system of linear equations leading to cell parameters a, b, c and cell volume V, their coefficients having being determined by multiple linear regressions. With similar objectives, Della Giusta et al. (1990) developed an alternative system of linear equations relating site occupancies as independent variables to the 24 different interatomic distances, cell parameters and volumes in olivines. The coefficients in their equations were also dependent on ionization potentials and polarizabilities of the ions. Thus, the interatomic distances in (Mg, Fe, Ca, Mn, Co, Ni, Zn) silicate olivines could be predicted satisfactorily on this basis.

In the current work, no attempt is made to relate chemical composition or site occupancy explicitly to interatomic distances or cell parameters. Nor are ionic radii, ionization potentials and polarizabilities used. The independent parameters adopted are the crystallographic parameters relating to atomic coordinates and unit-cell constants. For a generalized olivine, (M1)(M2)TO₄, in space group Pbnm there are 14 such parameters, i.e. a, b, c, x(O1), y(O1), x(O2), y(O2), x(O3), y(O3), z(O3), x(M2), y(M2), x(T) and y(T). Ion M1 does not enter the parameterization, since it is located at a centre of symmetry, 000.

Some elements of the current work have been previously used (Thomas & Beitollahi, 1994; Thomas, 1996, 1998). Here the oxygen ion network was used to derive relationships between octahedral tilt angle and the volume ratio of AO₁₂ and BO₆ coordination polyhedra in perovskites, ABO₃. Linked ionic networks are likewise used here, with the possibility of incorporating the volumes of cation coordination polyhedra in the algorithm for structure prediction.

The article is arranged as follows. In §2 the interactional representation of olivines is defined for both the oxygen and the M2—T2 cation networks, which are defined independently of each other. Since this process of Ionic Network Analysis is crucial, it is proposed that the acronym INA be used to describe the whole method. In §3, values are listed of the calculated interactional parameters for olivines under varying (p, T, X) conditions. Subsequently the algorithm for crystal structure prediction at interpolated values of p, T or X is described and applied in §4. Finally in §5, the results and methodology are discussed.

2. The interactional representation of olivines

2.1. The oxygen ionic network

Continuous pathways of pairwise interactions between oxygen ions are to be defined, so that the unit-cell parameters may be recalculated. These interactions are generally octahedral stalks or tetrahedral edges, which are characterized by lengths and orientational angles to the crystal axes. In order to define the method, the structure of an (Mg,Fe)₂SiO₄ olivine at 293 K (Heinemann et al., 2006) is taken into consideration (Table 1).

Table 1
Cell parameters (Å) and oxygen substructure of an olivine in space group Pbnm

Cell parameters			O1 in 4c		O2 in 4c		O3 in 8d
a	b	c	x	y	x	y	x	y	z
4.7942	10.3500	6.0436	0.76627	0.09204	0.21521	0.45025	0.28278	0.16425	0.03493

The ten degrees of freedom here correspond to the crystallographic representation. Therefore, the equivalent interactional representation, which is to be based on links between oxygen ions, will also require ten variables. The mirror planes at in space group Pbnm impose severe restrictions on the oxygen substructure: O1 and O2 ions lie within these planes, with pairs of O3 ions equidistant from them. The intermediate planes at are characterized by centres of symmetry at (x, y) = (0, 0); (½, ½), at which the M1 ions are located (Fig. 1).

Figure 1 The olivine structure in the z = 0 plane. a, b: orthorhombic cell parameters; s1,xy: projection of O3—O3′ stalk of centrosymmetric (M1)O6 octahedron in plane; t1: tie-line between stalk projections; ϕ1, ϕ2: angles of stalk projection or tie-line to the y-axis; (+, −): deviation of oxygen ion z-coordinates out of plane; JK, LM: b glide planes at x = ¼,¾.

Variables s_1,xy, t₁, ϕ₁ and ϕ₂ are parameters in the interactional representation. Recalculation of cell parameters a and b can proceed from just these four parameters. In this connection, continuous pathways OABCDE and ABF may be used to formulate equations (1) and (3), respectively.

These terms correspond to the projections of lengths OA, AB, BC, CD and DE, respectively, on to the y-axis. Simplification of equation (1) leads to equation (2).

With respect to pathway ABF, the corresponding simplified relationship is as follows.

A similar structure applies to the plane z = ½, the only difference being a sign-reversal of the O3 out-of-plane deviations of Fig. 1. Although the green rods linking the oxygen ions in Fig. 1 have been introduced as projections of octahedral stalks, they could also be regarded as projections of planes parallel to the z-axis, which are shown in Fig. 2.

Figure 2 The olivine structure in clinographic projection. Green planes contain O3—O3′ octahedral stalks. Red circles represent O3 oxygen ions. Pink planes correspond to one of the three girths of each (M2)O6 octahedron. Each light-blue-coloured link is the edge of a TO4 tetrahedron.

The green planes in Fig. 2 have three functions: (a) to contain the O3—O3′ stalks of (M1)O₆ octahedra; (b) to provide fixed mountings for one of the three girths of (M2)O₆ octahedra; (c) to provide hinge-mountings for TO₄ tetrahedra (not shown).

The O3 oxygen-ion structure also enables a continuous pathway to be formed in the c-direction. Thus, two further parameters are brought into the interactional set, i.e. s_1,z and e_z (Fig. 3).

Figure 3 Parameters s1z and ez together with continuous pathway ABCDE. Crosses: O3 ions.

The notation e_z stands for a tetrahedral edge in the z-direction. Analysis of pathway ABCDE leads to the relationship , which may be simplified to give equation (4).

The remaining four variables in the interactional set are associated with the z = ¼ plane (Fig. 4).

Figure 4 The z = ¼ plane, with projections of green planes (Fig. 2) and their tie-lines (Fig. 1) drawn in. (M2)O6 octahedra and TO4 tetrahedra are shown in pink and light blue, respectively. Contributions of this plane and the symmetry-related z = ¾ plane to an (M1)O6 octahedron at ½½½ are also shown. The definitive pathway ABCDE is made up alternately of tetrahedral edges exy and octahedral stalks s2. Parameters t2 and ft relate to the oxygen ion structure of TO4 tetrahedra.

The six interactional variables introduced so far relate to the O3 x,y,z coordinates and the three cell parameters a, b, c in the crystallographic representation. The remaining four interactional variables, e_xy, s₂, t₂ and f_t relate to the x and y crystallographic parameters of atoms O1 and O2, whose centres lie in the z = ¼ plane. The first of these parameters represents the length of TO₄ tetrahedral edges lying in the xy plane and the second the O1—O2 stalk-length of (M2)O₆ octahedra. It is these lengths that make up the pathway ABCDE in the y-direction (Fig. 4). Thus, the sum of the projections of these edges and stalks along y is constrained to be equal to the b parameter. Parameter t₂ is the length of the perpendicular tie-line linking the hinge-mounting of each tetrahedron to the opposite edge of length e_xy. Parameter f_t is the fraction of length e_xy in the positive x direction at which the tie-line meets the edge.

In summary, a set of ten interactional parameters has been defined to correspond to crystallographic parameters a, b, c, x(O1), y(O1), x(O2), y(O2), x(O3), y(O3) and z(O3). Values of the ten interactional parameters, s_1,xy, s_1,z, s₂, e_z, e_xy, t₁, t₂, ϕ₁, ϕ₂ and f_t, may be straightforwardly calculated for known crystal structures. The reverse calculation from the six interactional parameters s_1,xy, s_1,z, e_z, t₁, ϕ₁ and ϕ₂ to crystallographic parameters a, b, c, x(O3), y(O3) and z(O3) proceeds equally straightforwardly. Calculation of crystallographic parameters x(O1), y(O1), x(O2) and y(O2) from interactional parameters e_xy, s₂, t₂ and f_t requires the rotational angle about the tetrahedral hinge-mounting to be found that is consistent with the value of the b-parameter. This calculation is best performed computationally (Thomas, 2016).

2.2. The cation network

Just as for the oxygen ion network, continuous pathways between cations are defined so that the unit-cell parameters may be recalculated. In olivines, the location of both M2 and and T1 ions within z = ¼,¾ planes leads to a tessellated hexagonal structure (Fig. 5).

Figure 5 Interactional parameters L1, L2, L3, α and θ2 for the cation network in the z = ¼ plane (see also Fig. 4). White and black circles denote M2 and T ions, respectively. Values of parameters a, b, Δx and θ3 result, whereas parameter Δy remains unfixed, as signified by the double-headed arrow.

Each hexagon is associated with three unique side-lengths, L₁, L₂ and L₃, which correspond to interaction-lengths between M2 and T1 ions. Angles α and θ₂ are required in order to recalculate cell parameters a and b, in accordance with equations (5) and (6).

Angle θ₃ is merely used to simplify the form of equation (5), it being defined by the constraint . Equation (5) results from the continuous pathway BCD and equation (6) from continuous pathway ABC, whereby M2 ions A and C are related by the glide plane JK.

It is only possible to fix the network within the unit cell in the x-direction, as specified by parameter Δx, whose value results from the five definitive interactional variables, i.e. L₁, L₂, L₃, α and θ₂. By contrast, parameter Δy is free to float. Given an x-coordinate of −Δx for ion A, the x-coordinate of ion C is constrained by the glide plane to be ½ + Δx. It follows from pathway ABC that the value of parameter Δx is given by equation (7).

Since parameter cannot be fixed by the parameters of the cationic network, it is necessary to anchor the whole cationic network relative to the oxygen network. The bond-valence method may be used for this purpose, employing the functional form of Brown & Altermatt (1985) [equation (8)].

Here R_ij is the bond length and S_ij the bond valence. Parameters R₀ and B were taken from the so-called GRG (generalized reduced gradient) set of Gagné & Hawthorne (2015) in the current work. With respect to the 293 K structure of (Mg,Fe)₂SiO₄ (Heinemann et al., 2006), a bond-valence sum of 3.8870 arises for the coordination tetrahedra, which are fully occupied by Si ions. Similarly a bond-valence sum of 1.8979 applies to the M2 ions in their coordination octahedra, assuming the following occupancy factors, as specified by Heinemann et al.: Mg²⁺: 0.5291; Fe²⁺: 0.4225; Mn²⁺: 0.0402; Ca²⁺: 0.0082. Although these calculated bond-valence sums deviate from the idealized values of 4.0000 and 2.0000, this is not important here, since the objective is merely to fix the Δy value of the cation network. In this connection, it may be argued that the correct value of Δy for this structure is associated with bond-valence sums of 3.8870 and 1.8979 in the Si and M2 coordination polyhedra.

Since the absolute values of the bond-valence sums are of secondary importance, a further simplification results from the assumption that the M2 site is always fully occupied by Mg²⁺ ions. Under this assumption, the M2 bond valence for the 293 K structure of Heinemann et al. is calculated to be 1.7682. Similarly the tetrahedral sites may be assumed always to be fully occupied by Si⁴⁺ ions. This simplification eliminates the requirement of taking variable occupancy factors explicitly into account. Similarly these parameter values may be applied here to structures at elevated temperatures and non-ambient pressures, although bond-valence parameters generally apply to a temperature of 293 K and a pressure of 1 bar. Thus, the general validity of the functional form of equation (8) is assumed, in order to anchor the cationic network of Fig. 5. In the tables of the following section, the two bond-valence sums are denoted by the symbols ΣS(Si) and ΣS(Mg).

3. Values of interactional parameters for olivines under varying p, T, X conditions

In this section, calculated values of the interactional, or `INA' parameters for oxygen and cation networks are listed for known olivine structures. The solid solution system (Co_xMg_1 − x)₂SiO₄ (0 ≤ x ≤ 1) investigated by Müller-Sommer et al. (1997) is considered first of all in §3.1. The structural data of Heinemann et al. (2006) for (Mg,Fe)₂SiO₄ at temperatures between 293 and 1180 K are used to calculate the interactional parameter values listed in §3.2. Finally, the pressure-dependence of interactional parameters is considered in §3.3. To this end, the structural dataset of Kudoh & Takéuchi (1985) for forsterite (Mg₂SiO₄) between 31 and 149 kbar is considered in §3.3.1, augmented by the structural data of Lager et al. (1981) at atmospheric pressure (0.001 kbar) and Finkelstein et al. (2014) for a pressure of 453 kbar. In §3.3.2 the interactional parameters are listed corresponding to the structural data of Kudoh & Takéuchi (1985) for fayalite (Fe₂SiO₄) at pressures between 0.001 and 140 kbar.

There are two aims of the approach: (a) to predict the crystal structures at intermediate (p–T–X) values relative to known structural data; (b) to provide a tool for crystallographers in seeking optimum structural solutions from the results of diffraction experiments. A central concern is to identify the options for calculating INA parameters at intermediate (p–T–X) points from their values at experimentally determined (p–T–X) points, this being a process of interpolation. Thereafter, a reverse transformation from interpolated INA parameters to crystallographic parameters can take place, utilizing the equations in §2.

A well known option for performing the interpolations would be to use polynomial fits. This option is discussed in §5.3, in comparison to a novel approach, which is developed in §4. In this exploratory work, linear interpolations between experimental data-points are used. This simple technique has the advantage of giving equal a priori weighting to all experimental points. However, since linearity of parameter variation cannot be assumed to be correct, it is also necessary to define a means of shifting the parameter values away from the starting values obtained by linear interpolation. The guiding principle is to make these shifts as small as possible, consistent with the fulfilment of certain constraints. Three types of constraint are investigated here. First, that cell parameters a and b obtained either by equations (2) and (3) or (5) and (6) are in agreement with each other. Secondly, that cell parameters a, b and c are in agreement with quadratic curves applying to the particular system over the whole range of p, T or X. Thirdly, that the parameter values are in agreement with quadratic functions or narrow limits relating to the volumes of cation coordination polyhedra over the whole range of p, T or X for a particular system.

3.1. Parameter values for the system (Co_xMg_1 − x)₂SiO₄ (0 ≤ x ≤ 1)

Calculated values of the interactional parameters for the nine structures of Müller-Sommer et al. (1997) are given in Table 2. These had been determined by Rietveld refinements of X-ray powder diffraction data of samples prepared by solid-state synthesis. Standard errors quoted in Table 2 have been calculated from the standard errors in the crystallographic parameters, i.e. unit-cell parameters and atomic coordinates. For this purpose a numerical method was adopted, in which the normal distributions of each of these parameters were simultaneously sampled by 200 discrete points. For each sampling point the corresponding INA parameters were calculated. The sets of interactional parameters generated were subsequently used to calculate the standard errors of individual interactional parameters. Correlation errors were avoided by random sampling of the normal distributions of crystallographic parameters (Thomas, 2016).

Table 2
INA parameter values for (Co_xMg_1 − x)₂SiO₄ olivine calculated from the structural data of Müller-Sommer et al. (1997)

x	1.000	0.875	0.750	0.625	0.500	0.375	0.250	0.125	0.000
s1,xy (Å)	4.3388 (105)	4.3882 (80)	4.3165 (60)	4.2732 (77)	4.2942 (60)	4.2785 (40)	4.2782 (44)	4.2737 (56)	4.2564 (40)
s1,z (Å)	0.3987 (84)	0.4608 (96)	0.4522 (72)	0.4484 (84)	0.4350 (60)	0.4659 (48)	0.4419 (48)	0.4500 (48)	0.4500 (36)
s2 (Å)	4.3001 (101)	4.3546 (95)	4.2645 (68)	4.3283 (76)	4.2768 (63)	4.2744 (49)	4.2543 (49)	4.2494 (50)	4.2471 (41)
ez (Å)	2.6034 (84)	2.5395 (96)	2.5464 (72)	2.5487 (84)	2.5606 (60)	2.5283 (48)	2.5518 (48)	2.5423 (48)	2.5422 (36)
exy (Å)	2.6786 (81)	2.6778 (83)	2.6906 (60)	2.6580 (67)	2.6716 (58)	2.6931 (44)	2.7091 (44)	2.7123 (48)	2.7003 (37)
t1 (Å)	2.9924 (74)	2.9736 (49)	2.9780 (37)	2.9882 (49)	2.9790 (37)	2.9885 (25)	2.9745 (25)	2.9649 (37)	2.9707 (24)
t2 (Å)	1.8270 (73)	1.8334 (62)	1.8478 (46)	1.8374 (55)	1.8358 (44)	1.8276 (32)	1.8374 (33)	1.8292 (39)	1.8283 (29)
ϕ1 (°)	39.396 (124)	39.665 (103)	38.852 (78)	38.660 (98)	38.886 (79)	39.154 (53)	38.833 (60)	38.639 (71)	38.713 (53)
ϕ2 (°)	53.044 (189)	53.512 (128)	53.331 (95)	53.038 (126)	53.210 (95)	52.911 (63)	53.197 (63)	53.366 (95)	53.138 (63)
ft	0.4407 (31)	0.4573 (33)	0.4360 (24)	0.4201 (28)	0.4278 (23)	0.4265 (18)	0.4266 (18)	0.4266 (19)	0.4272 (14)
L1 (Å)	3.3123 (42)	3.2968 (37)	3.2869 (23)	3.2923 (37)	3.2761 (29)	3.2775 (23)	3.2655 (23)	3.2527 (37)	3.2463 (15)
L2 (Å)	2.7955 (38)	2.8252 (41)	2.8072 (29)	2.7935 (37)	2.8051 (30)	2.8050 (25)	2.7956 (25)	2.7981 (32)	2.7930 (20)
L3 (Å)	3.2776 (37)	3.2560 (42)	3.2789 (31)	3.2787 (37)	3.2726 (30)	3.2587 (26)	3.2747 (26)	3.2753 (31)	3.2804 (21)
θ2 (°)	41.763 (80)	41.476 (82)	41.938 (58)	41.900 (76)	41.898 (61)	41.665 (51)	42.036 (51)	42.143 (67)	42.310 (39)
α (°)	6.650 (60)	6.587 (76)	6.786 (59)	6.965 (65)	6.758 (53)	6.675 (47)	6.794 (47)	6.788 (48)	6.837 (42)
ΣS(Si)	4.1326	4.2718	4.1558	4.2329	4.1944	4.2322	4.1219	4.1561	4.1891
ΣS(Mg)	1.7850	1.7493	1.8047	1.7554	1.8035	1.7880	1.8217	1.8310	1.8288

In general, the parameter values are found to vary within narrow intervals, with the difference in values for a particular parameter between neighbouring x values generally greater than the sum of the two relevant standard errors.

In connection with the linear interpolation method to be used in connection with the novel approach, the piecewise linear variation of parameter values between adjacent x-points requires the use of equations such as equation (9). As an example, parameter e_xy has been taken, with its value at x = 0.800 calculated.

3.2. Parameter values for forsterite (Mg₂SiO₄) at different temperatures

Interactional parameters calculated from the structural data of Heinemann et al. (2006) for forsterite between 293 and 1180 K are given in Table 3. There is a sizeable gap in temperature values between 579 and 874 K. Furthermore, the data apply to two different single crystals, the first associated with the eight temperature columns to the left of the table, and the second with the four temperature columns to the right. The standard errors in the interactional parameter values are generally smaller than for the solid-solution system of §3.1, this probably being fundamentally due to the use of single-crystal rather than powder X-ray diffraction. A significantly greater tendency is observed for parameter values to vary monotonically over the whole temperature range, whereby it should be noted that the temperature of the ninth data column is lower than the temperature in the eighth.

Table 3
Parameter values in the Ionic Network Analysis corresponding to the 12 structures of Heinemann et al. (2006) for (Mg,Fe)₂SiO₄

The first eight columns refer to sample `Bo-10' and the final four to sample `Bo-2' (ibid).

Temperature (K)	293	377	475	579	874	929	974	1026	1021	1077	1125	1180
s1,xy (Å)	4.3476 (21)	4.3510 (18)	4.3554 (21)	4.3625 (21)	4.3765 (21)	4.3811 (21)	4.3848 (21)	4.3886 (21)	4.3901 (21)	4.3947 (21)	4.3991 (21)	4.4053 (21)
s1,z (Å)	0.4255 (12)	0.4270 (12)	0.4311 (12)	0.4329 (12)	0.4432 (24)	0.4472 (24)	0.4511 (24)	0.4540 (24)	0.4502 (24)	0.4531 (24)	0.4572 (24)	0.4576 (24)
s2 (Å)	4.2872 (16)	4.2917 (16)	4.2957 (17)	4.3024 (15)	4.3246 (17)	4.3253 (17)	4.3303 (23)	4.3341 (23)	4.3323 (17)	4.3363 (17)	4.3438 (17)	4.3495 (23)
ez (Å)	2.5963 (13)	2.5970 (13)	2.5961 (13)	2.5986 (13)	2.6010 (24)	2.5991 (25)	2.5970 (25)	2.5971 (25)	2.6002 (25)	2.5998 (25)	2.5992 (25)	2.6012 (25)
exy (Å)	2.7330 (15)	2.7347 (15)	2.7362 (19)	2.7365 (14)	2.7376 (19)	2.7390 (19)	2.7392 (22)	2.7410 (22)	2.7408 (19)	2.7404 (19)	2.7412 (19)	2.7412 (22)
t1 (Å)	2.9846 (13)	2.9859 (13)	2.9896 (13)	2.9915 (13)	3.0047 (13)	3.0065 (13)	3.0081 (13)	3.0118 (13)	3.0099 (13)	3.0119 (13)	3.0158 (13)	3.0167 (13)
t2 (Å)	1.8710 (13)	1.8701 (12)	1.8680 (13)	1.8699 (13)	1.8672 (13)	1.8672 (13)	1.8673 (15)	1.8652 (15)	1.8662 (13)	1.8659 (13)	1.8645 (13)	1.8651 (15)
ϕ1 (°)	38.619 (27)	38.634 (21)	38.686 (27)	38.671 (27)	38.770 (27)	38.796 (27)	38.811 (27)	38.857 (27)	38.845 (27)	38.878 (27)	38.913 (27)	38.911 (27)
ϕ2 (°)	53.433 (33)	53.428 (33)	53.391 (33)	53.416 (33)	53.299 (32)	53.296 (32)	53.292 (32)	53.260 (32)	53.286 (32)	53.290 (32)	53.246 (32)	53.274 (32)
ft	0.4316 (6)	0.4313 (5)	0.4323 (7)	0.4329 (6)	0.4327 (7)	0.4350 (7)	0.4344 (7)	0.4344 (7)	0.4345 (7)	0.4349 (7)	0.4351 (7)	0.4353 (7)
L1 (Å)	3.3053 (7)	3.3059 (7)	3.3081 (8)	3.3103 (7)	3.3181 (8)	3.3191 (12)	3.3204 (9)	3.3214 (9)	3.3226 (8)	3.3235 (9)	3.3259 (9)	3.3275 (10)
L2 (Å)	2.8399 (7)	2.8420 (7)	2.8444 (7)	2.8482 (6)	2.8573 (7)	2.8594 (9)	2.8611 (8)	2.8648 (8)	2.8630 (7)	2.8648 (8)	2.8679 (8)	2.8699 (9)
L3 (Å)	3.2810 (8)	3.2831 (8)	3.2868 (8)	3.2917 (7)	3.3071 (8)	3.3104 (9)	3.3129 (9)	3.3161 (9)	3.3150 (8)	3.3183 (8)	3.3223 (8)	3.3260 (9)
θ2 (°)	41.848 (14)	41.869 (14)	41.896 (14)	41.935 (13)	42.062 (15)	42.093 (20)	42.115 (16)	42.132 (16)	42.122 (15)	42.142 (15)	42.180 (16)	42.209 (18)
α (°)	7.048 (10)	7.034 (9)	7.033 (10)	7.013 (10)	6.956 (12)	6.956 (12)	6.942 (13)	6.921 (13)	6.915 (12)	6.902 (12)	6.885 (13)	6.868 (14)
ΣS(Si)	3.8870	3.8841	3.8902	3.8794	3.8804	3.8834	3.8868	3.8883	3.8803	3.8833	3.8873	3.8818
ΣS(Mg)	1.7682	1.7613	1.7511	1.7404	1.6995	1.6945	1.6870	1.6777	1.6823	1.6751	1.6627	1.6561

3.3. Parameter values at different pressures

3.3.1. Forsterite, Mg₂SiO₄

Interactional parameters calculated from the structural data of Lager et al. (1981), Kudoh & Takéuchi (1985) and Finkelstein et al. (2014) for forsterite at pressures between 0.001 and 453 kbar are given in Table 4. Standard deviations of parameter values for samples under pressure are generally significantly higher than in §§3.1 and 3.2, thus reflecting the experimental difficulties associated with using pressure cells. Crystallographic parameters had been determined universally from single-crystal data, in the case of the 453 kbar data of Finkelstein et al. utilizing a synchrotron beamline. Significant departures of parameter values from monotonic trends with pressure dictate the use of the piecewise interpolation method for intermediate compositions, as discussed in connection with equation (9).

Table 4
INA parameter values for forsterite calculated from the structural data of Lager et al. (1981) at 0.001 kbar, Kudoh & Takéuchi (1985) at 31–149 kbar and Finkelstein et al. (2014) at 453 kbar

Pressure (kbar)	0.001	31	47	53	79	86	111	149	453
s1,xy (Å)	4.2453 (17)	4.2208 (134)	4.2008 (95)	4.1724 (94)	4.1660 (78)	4.1703 (78)	4.1991 (193)	4.1977 (468)	3.9426 (112)
s1,z (Å)	0.3960 (12)	0.3827 (594)	0.3883 (165)	0.3762 (212)	0.3939 (234)	0.3928 (222)	0.3758 (233)	0.4744 (505)	0.2585 (67)
s2 (Å)	4.2233 (13)	4.2022 (173)	4.1682 (97)	4.1587 (98)	4.1306 (102)	4.1155 (96)	4.1006 (211)	4.1433 (569)	3.8608 (118)
ez (Å)	2.5947 (12)	2.5883 (594)	2.5622 (165)	2.5718 (212)	2.5366 (234)	2.5297 (222)	2.5422 (233)	2.3975 (505)	2.5149 (67)
exy (Å)	2.7472 (10)	2.6753 (148)	2.7602 (93)	2.7077 (97)	2.7037 (94)	2.6894 (90)	2.6505 (201)	2.5176 (483)	2.6305 (95)
t1 (Å)	2.9649 (12)	2.9309 (84)	2.9289 (59)	2.9311 (60)	2.9041 (47)	2.8903 (46)	2.8796 (115)	2.7934 (260)	2.7131 (72)
t2 (Å)	1.8607 (11)	1.8631 (112)	1.8373 (69)	1.8303 (69)	1.8188 (65)	1.8345 (64)	1.7994 (145)	1.8959 (356)	1.8239 (84)
ϕ1 (°)	38.404 (20)	38.500 (177)	38.706 (124)	38.634 (125)	38.682 (106)	38.507 (105)	39.517 (260)	37.338 (637)	37.250 (142)
ϕ2 (°)	53.284 (32)	53.697 (222)	53.617 (158)	53.445 (157)	53.818 (127)	54.141 (127)	54.146 (318)	56.357 (801)	56.469 (229)
ft	0.4293 (4)	0.4194 (61)	0.4243 (34)	0.4195 (34)	0.4196 (35)	0.4178 (34)	0.4131 (76)	0.4114 (198)	0.4104 (41)
L1 (Å)	3.2524 (14)	3.2814 (71)	3.1961 (43)	3.1945 (43)	3.1844 (50)	3.1726 (42)	3.1724 (84)	3.1386 (188)	3.0160 (72)
L2 (Å)	2.7882 (11)	2.7247 (69)	2.7551 (45)	2.7511 (45)	2.7169 (50)	2.7326 (45)	2.7027 (99)	2.6951 (223)	2.5863 (52)
L3 (Å)	3.2728 (10)	3.2036 (68)	3.2424 (46)	3.2294 (46)	3.2121 (50)	3.1954 (45)	3.1695 (102)	3.1585 (230)	3.0205 (48)
θ2 (°)	42.133 (24)	40.834 (146)	41.977 (93)	41.804 (93)	41.628 (106)	41.402 (93)	40.979 (202)	41.006 (457)	39.610 (132)
α (°)	6.883 (12)	7.012 (117)	6.924 (84)	6.934 (84)	6.849 (91)	6.946 (85)	6.243 (196)	6.783 (451)	7.801 (63)
ΣS(Si)	3.8810	4.0318	4.0705	4.0898	4.2127	4.2071	4.3919	4.6960	4.3861
ΣS(Mg)	1.8702	1.9451	1.9907	2.0017	2.0638	2.0995	2.1368	2.2349	2.7937

3.3.2. Fayalite, Fe₂SiO₄

INA parameters calculated from the structural data of Kudoh & Takeda (1986) for fayalite at pressures between 0.001 and 140 kbar are given in Table 5. Standard deviations of parameter values are generally high, together with the general absence of monotonic trends. Both factors necessitate the use of the piecewise interpolation method.

Table 5
INA parameter values for fayalite calculated from the single-crystal structural data of Kudoh & Takeda (1986) at pressures up to 140 kbar

Pressure (kbar)	0.001	49	67	93	117	140
s1,xy (Å)	4.4367 (65)	4.4215 (221)	4.3299 (337)	4.3356 (416)	4.8509 (755)	4.6504 (1059)
s1,z (Å)	0.4631 (73)	0.4156 (181)	0.4857 (157)	0.4549 (203)	0.5522 (251)	0.5441 (323)
s2 (Å)	4.3597 (68)	4.3190 (257)	4.2585 (356)	4.1792 (416)	4.1684 (868)	4.0868 (1168)
ez (Å)	2.5839 (74)	2.6049 (183)	2.5273 (158)	2.5376 (205)	2.4358 (252)	2.4454 (323)
exy (Å)	2.7201 (68)	2.7100 (209)	2.6786 (258)	2.6924 (309)	2.7037 (611)	2.6716 (741)
t1 (Å)	2.9999 (38)	2.9211 (144)	2.9820 (247)	2.9269 (297)	2.6504 (401)	2.6917 (619)
t2 (Å)	1.8784 (46)	1.9217 (169)	1.8665 (249)	1.8565 (297)	1.9438 (587)	1.9684 (831)
ϕ1 (°)	38.774 (85)	37.988 (257)	39.563 (386)	38.970 (466)	36.998 (715)	37.841 (1055)
ϕ2 (°)	53.484 (97)	55.126 (401)	53.160 (634)	54.354 (810)	63.796 (1762)	61.679 (2450)
ft	0.4287 (24)	0.4153 (82)	0.4296 (112)	0.4288 (130)	0.4540 (286)	0.3755 (392)
L1 (Å)	3.3513 (25)	3.2785 (111)	3.2708 (148)	3.2257 (172)	3.2345 (286)	3.1466 (358)
L2 (Å)	2.8743 (23)	2.8542 (84)	2.8348 (100)	2.8182 (123)	2.8036 (212)	2.8120 (260)
L3 (Å)	3.2956 (23)	3.2783 (81)	3.2532 (92)	3.2584 (114)	3.2162 (198)	3.1907 (241)
θ2 (°)	41.838 (46)	41.877 (178)	41.605 (217)	41.908 (267)	40.970 (469)	40.768 (585)
α (°)	6.934 (37)	7.595 (105)	7.370 (113)	8.028 (146)	7.579 (264)	8.348 (298)
ΣS(Si)	3.9207	3.7859	4.1768	4.1763	4.2153	4.3871
ΣS(Mg)	1.6719	1.8069	1.7857	1.9151	2.1788	2.1776

4. Crystal structure prediction at interpolated values of X, T or p

4.1. The system (Co_xMg_1 − x)₂SiO₄ (0 ≤ x ≤ 1)

Table 2 may be used as a basis for predicting crystal structures at x values interpolated between the experimentally determined crystal structures. For this purpose it was found expedient to use Microsoft^® Excel 2010. The parameter values in Table 2 were pasted into an Excel spreadsheet, in order to act as a lookup table. Calculation of interpolated interactional parameters for x = 0.800 is shown, for example, in Fig. 6.

Figure 6 Screen dump from the Excel spreadsheet (German version) used to calculate starting values of INA parameters by linear interpolation. The associated crystallographic parameters, calculated by reverse transformation, are also shown.

The value 0.800 is entered for x in cell B1. This value calls up the two columns of parameter values in the lookup-table nearest to x = 0.800, i.e. x = 0.750 and x = 0.875 in columns B and C. Interpolated values for x = 0.800 are calculated in column D by means of equation (9). The crystallographic parameters associated with the interpolated parameters in column D are calculated by Excel formulae and output in the box occupying rows 6–21 and columns G–J. Different values are observed for cell parameters a and b in cells H10 and H11 compared with H19 and H20, depending on whether the oxygen ionic or the cationic network is used for the calculation. In the former case, equations (1) and (3) apply, and in the latter, equations (5) and (6).

The two sets of values for a and b may be brought into agreement with each other by use of the Solver add-in tool within Excel. In order to do this, the values of cells in the box of Fig. 6 are copied to a second box of crystallographic parameters, as shown in Fig. 7(a). The interactional parameters corresponding to these crystallographic parameters are automatically back-calculated and listed in column E of Fig. 6. A refinement is then carried out, whereby the r.m.s. shift (Fig. 6; cell F23) between starting parameter values (Fig. 6; column D) and back-calculated parameters (Fig. 6; column E) is minimized. In so doing, the Solver is empowered to change the values in the cells representing variable atomic coordinates and cell parameters (Fig. 7; columns N–P), subject to the constraint that both sets of a, b cell parameters be equal. The extent of the changes carried out is limited by minimizing the r.m.s. shift. Thus, the starting values of the INA parameters, obtained by linear interpolation, should correspond to a crystal structure close to the optimal structure, but must not generate this directly. Employing the GRG non-linear algorithm within the Solver, unified cell constants a = 4.77892  and b = 10.28521 Å were obtained (Fig. 7b). Small differences in some atomic coordinates between Figs. 7(a) and (b) are observed. The r.m.s. shift for the refinement, which arises from discrepancies between starting and refined values of interactional parameters, was calculated to be 0.001%.

Figure 7 Calculated crystallographic parameters: (a) prior to refinement; (b) following refinement under the condition that both sets of a and b parameters be equal to each other. (c) following refinement under the condition that all three sets of cell parameters be equal to one other; (d) following setting of Δy for the cationic network through the closest fit to interpolated bond-valence values.

The flexibility of this approach allows additional, normative conditions to be applied to the refinement. Thus, a second refinement was carried out, this employing a third set of unit-cell parameters generated by fitting quadratic functions to the unit cell data of Müller-Sommer et al. [equation (10)].

The calculated values for x = 0.800 are given in cells E1 to E3 of Fig. 6. Upon use of the condition that all three sets of cell parameters, i.e. oxygen ionic network, cationic network and normative, be equal to one other, the results given in Fig. 7(c) were obtained, these being associated with an r.m.s. shift in INA parameter values of 0.005%.

The final stage in the structural prediction is to set the parameter Δy of the cationic network. This is carried out by clicking the button `Determine Delta Y' in the Excel spreadsheet (Fig. 7c). Upon so doing, an Excel–VBA–Macro is activated to calculate the value of Δy that gives the closest agreement with interpolated bond-valence sums of 4.20220 for ΣS(Si) and 1.78254 for ΣS(Mg) (see Fig. 6; cells D22 and D23). A value for Δy of 0.27083 was obtained, as observed in Fig. 7(d). The box in Fig. 7(d) contains the final predicted crystallographic parameters for the olivine (Co_0.8Mg_0.2)₂SiO₄.

4.2. Forsterite, Mg₂SiO₄, at variable temperatures

Since Heinemann et al. (2006) did not solve the structures of forsterite at temperatures between 579 and 874 K, an appropriate exercise is to predict the crystal structures at temperatures of 673 and 773 K, corresponding to 400 and 500°C, respectively. Accordingly, the following normative quadratic functions were derived from the cell parameters of the 12 structures, in order to link temperature T in Kelvin with a, b and c [equation (11)].

By applying procedures in Excel similar to those in §4.1, the crystallographic parameters listed in Table 6 for temperatures of 673 and 773 K were obtained. The r.m.s. shifts quoted are the mean shifts of the INA parameters from their starting values in the two refinements.

Table 6
Predicted crystal structures of forsterite at 673 and 773 K employing cell parameters as normative constraints

	T = 673 K			T = 773 K
Unit-cell parameters (Å)	a	b	c	a	b	c
	4.80804	10.38781	6.07017	4.81273	10.40055	6.07853
Parameter	x	y	z	x	y	z
O1	0.76494	0.09233	1/4	0.76464	0.09236	1/4
O2	0.28457	−0.04902	1/4	0.28460	−0.04878	1/4
O3	0.28393	0.16400	0.03593	0.28417	0.16390	0.03616
Si	0.07156	0.59575	1/4	0.07156	0.59575	1/4
M2	−0.01233	0.27929	1/4	−0.01208	0.27941	1/4
R.m.s. shift (%)	0.011			0.012

In order to provide a further constraint, the volumes of Si, M1 and M2 coordination polyhedra were calculated for the 12 structures. Although the volumes of the SiO₄ tetrahedra were found to lie within a narrow band between 2.2115 and 2.2167 ų, no monotonic variation with temperature was found. In contrast, the volumes of (M1)O₆ and (M2)O₆ octahedra lie in the following broader ranges: 12.3201 ≤ V(M1) ≤ 12.7347 ų; 12.8078 ≤ V(M2) ≤ 13.2904 ų. Further, their variation with temperature could be satisfactorily modelled by the quadratic functions of (12).

The three constraints above were built into the refinement procedure: changes in the volumes of the Si, M1 and M2 coordination octahedra occur during the refinement when the Solver alters the variable parameters for the O1, O2 and O3 ions (see Figs. 7a–d). Accordingly, three functions were written in Excel–VBA to calculate the SiO₄, (M1)O₆ and (M2)O₆ volumes from these parameters based on the formula for the volume of a pyramid of basal area b and height h, either as a complete or partial polyhedron (Thomas, 1991). The predicted crystal structures, again at 673 and 773 K, are given in Table 7, whereby the changes in atomic coordinates at 673 K compared with Table 6 are small. The changes in oxygen ion coordinates at 773 K between Tables 6 and 7 are slightly larger, this being consistent with the higher r.m.s. shift obtained when polyhedral constraints are applied.

Table 7
Predicted crystal structures of forsterite at 673 and 773 K employing cell parameters and polyhedral volumes as normative constraints

	T = 673 K			T = 773 K
Unit-cell parameters (Å)	a	b	c	a	b	c
	4.80804	10.38781	6.07017	4.81273	10.40055	6.07853
Parameter	x	y	z	x	y	z
O1	0.76492	0.09237	1/4	0.76458	0.09251	1/4
O2	0.28466	−0.04900	1/4	0.28495	−0.04874	1/4
O3	0.28393	0.16401	0.03593	0.28417	0.16392	0.03616
Si	0.07156	0.59596	1/4	0.07156	0.59621	1/4
M2	−0.01233	0.27949	1/4	−0.01208	0.27987	1/4
R.m.s shift (%)	0.012			0.020
V(Si) (Å^3)	2.2154			2.2144
V(M1) (Å^3)	12.4565			12.5036
V(M2) (Å^3)	12.9747			13.0288

4.3. Forsterite and fayalite at variable pressures

4.3.1. Forsterite, Mg₂SiO₄

Just as for the temperature dependence of (Mg,Fe)₂SiO₄ olivine in §4.2, the structural data for forsterite at variable pressure allow the formation of normative relationships for the unit-cell parameters and the two octahedral volumes [equation (13)]. The pressure, p, is given in kbar here.

These allow the stable prediction of crystal structures at the two intermediate pressures of 250 and 350 kbar, for example (Table 8).

Table 8
Predicted crystal structures of forsterite at pressures of 250 and 350 kbar, employing cell parameters and polyhedral volumes as normative constraints

	p = 250 kbar			p = 350 kbar
Unit-cell parameters (Å)	a	b	c	a	b	c
	4.59190	9.53590	5.65453	4.55147	9.37764	5.58301
Parameter	x	y	z	x	y	z
O1	0.76331	0.08409	1/4	0.76478	0.08712	1/4
O2	0.28518	−0.05456	1/4	0.28049	−0.05696	1/4
O3	0.27093	0.17067	0.03558	0.26754	0.17020	0.02968
Si	0.07249	0.59659	1/4	0.07374	0.60181	1/4
M2	−0.01052	0.27670	1/4	−0.01282	0.28086	1/4
R.m.s. shift (%)	0.4181			0.4179
V(Si) (Å^3)	1.9407			1.9580
V(M1) (Å^3)	10.3215			9.9299
V(M2) (Å^3)	10.4120			9.9499

The higher values of r.m.s. shift of ca 0.42% in Table 8 compared with Tables 6 and 7 indicate that the Solver has applied larger shifts from starting parameter values in reaching the normative values of equations (13).

4.3.2. Fayalite, Fe₂SiO₄

The structural data of Kudoh & Takeda (1986) do not allow the use of monotonic normative relationships for unit-cell constants and polyhedral volumes over the whole pressure range to 140 kbar. Attempts to use normative relationships over narrower pressure ranges met with no success with regard to the prediction of the crystal structure at an example pressure of 25 kbar. This indicates that such normative relationships impose stringent conditions on the refinement, and that they must be essentially correct in order to achieve a structural prediction. Consequently, the only condition applied in order to predict a structure at 25 kbar was to require that the a and b unit-cell parameters arising from oxygen ion and cation networks are equal to one another. The results are given in Table 9.

Table 9
Crystal structure of fayalite predicted at 25 kbar without normative constraints

Unit-cell parameters (Å)	a	b	c
	4.79275	10.31482	6.04084
Parameter	x	y	z
O1	0.77345	0.08753	1/4
O2	0.28751	−0.04680	1/4
O3	0.28431	0.16910	0.03439
Si	0.07349	0.59432	1/4
M2	−0.01697	0.27911	1/4
R.m.s. shift (%)	0.024

5. Discussion

Three questions arise naturally from the method advocated: firstly, why the transformation from crystal structural to interactional parameters is desirable; secondly, why a structural refinement is necessary; and thirdly, why a piecewise linear interpolation method has been used, instead of polynomial fitting, in order to determine the starting values for the structural refinements. In order to address these issues, the structural data of Heinemann et al. (2006) are taken as a basis, as analysed in §3.2 and §4.2. This is primarily because this dataset is associated with the lowest standard errors in the parameter values. The three questions are considered in separate subsections.

5.1. Motivation for the INA transformation

Fig. 8 is a graphical representation of the data in Table 3, in which the parameters have been divided into four groups. Since the curves here have been drawn by linear interpolation between adjacent points, they represent starting values of the parameters in the structural refinements. The curves in Fig. 8(a), which are associated with cationic parameters, show clear, monotonic trends. This is also the case for most of the other parameters in Figs. 8(b) and (c), although jagged behaviour is a feature of several, notably e_xy, ϕ₂, t₂, e_z and f_t. Since the three length parameters here, i.e. e_xy, t₂ and e_z, have small left-to-right Δ values of 0.008, −0.006 and 0.005 Å, respectively, this jaggedness indicates variations within very narrow bandwidths. Similar considerations apply to the fractional parameter f_t, which has a Δ value of just 0.003. The relatively long error bars observed for parameters e_xy, t₂, e_z and f_t are thus an artefact of their Δ values being small: the errors in these parameters are not fundamentally larger than the errors in the other parameters.

Figure 8 Variation of the 15 length- and angle-based interactional parameters of forsterite with temperature (data of Table 3): (a) cationic parameters; (b) anionic parameters with an ascending trend; (c) anionic parameters with a descending trend; (d) approximately constant anionic parameters. The black crosses correspond to solutions of the Excel Solver for temperatures of 673 and 773 K, together with temperatures between 998 and 1123 K in 25 K intervals. Unit-cell and polyhedral volume constraints have been applied.

The black crosses at temperatures of 673 and 773 K indicate the INA values associated with the structural predictions of Table 7. The action of the Solver for structural predictions at higher temperatures between 998 and 1123 K is also indicated. This temperature region was chosen for scrutiny in the diagrams because of the relative unevenness of the curves. The r.m.s. shifts obtained were as follows: 998 K: 0.007%; 1023 K: 0.012%; 1048 K: 0.009%; 1073 K: 0.018%; 1098 K: 0.006%; 1123 K: 0.011%. A deviation of a cross from the corresponding curve is indicative of a departure of the end-point of the refinement from the starting value.

The cell parameters are generally fixed in the refinements by applying normative constraints. Therefore, it is logical to question whether the variations of untransformed x,y,z atomic parameters with temperature could also have been used for interpolations, since these are ultimately the only parameters free to vary. The behaviour of these parameters is shown in Fig. 9.

Figure 9 Variation of x,y,z fractional coordinates of the ions with temperature, with straight lines linking neighbouring points. The data and the error bars have been taken directly from Heinemann et al. (2006). (a) Cationic parameters; (b) anionic parameters with an ascending trend; (c) anionic parameters with a descending trend; (d) approximately constant anionic parameters.

A comparison of Figs. 8(a) and 9(a) points directly to the need for the INA transformation. Although both figures are based on the identical atomic parameters, i.e. x(Si), y(Si), x(M2) and y(M2), the INA parameters of Fig. 8(a) give rise to monotonic trends, whereas parameters x(Si) and y(Si) of Fig. 9(a) do not. Thus it is possible to generate reliable interpolated values from the INA parameters (as starting values for the refinements), but not from the x,y,z fractional coordinates. The structural reason for this is to be found in Fig. 5: INA parameters L₁, L₂, L₃, α and θ₂ relate directly to a two-dimensional planar network, which is constrained to respond in a correlated way to changes in temperature.

Similar considerations apply to the planar network of Fig. 1, which is associated with atomic parameters x(O3), y(O3) and INA parameters s_1,xy, ϕ₁, t₁, ϕ₂. Parameters s_1,xy and t₁ vary monotonically with temperature, whereas parameters ϕ₁ and ϕ₂ show monotony apart from the third point from the left, corresponding to a temperature of 475 K. The upward kink in the ϕ₁ curve and the downward kink in the ϕ₂ curve are correlated with each other. Given the error bars in the y(O3) parameter at this temperature (Fig. 9c), it is likely that the observed kinks are an artefact of the structural solution of Heinemann et al. at 475 K. Insights of this kind, provided by trends in the network parameters, would be of value to crystallographers whilst seeking to optimize structural refinements from diffraction data.

The four INA parameters already identified as having values within narrow bandwidths, i.e. e_xy, t₂, e_z and f_t, are all associated with the oxygen ions of the [SiO₄]⁴⁻ tetrahedra. Thus, this behaviour indicates strong chemical bonding restraints over the whole temperature range. Given this tetrahedral rigidity, it follows that flexibility in the value of the s₂ parameter is required: this parameter is associated with the (M2)O₆ octahedra, and its variation allows the length of chain ABCDE to adapt to changes in the cell parameter b (Fig. 4). Of all the INA length parameters, s₂ has the largest Δ value of 0.063 Å (Fig. 8). Adaptation to changes in cell parameter a is enabled by this parameter acting in combination with the rotational hinge (Fig. 4). The remaining INA parameter, s_1,z, is unconstrained. Acting together with the rigid parameter e_z [equation (4); Fig. 3], its variation leads to changes in the cell parameter c.

In the case of a rigid, i.e. narrow bandwidth INA parameter, linear interpolation does not give rise to significant errors. For systematically varying INA parameters, linear interpolations likewise lead to reliable starting values for the structural refinements. By comparison, the variations in untransformed atomic parameters (Fig. 9) do not allow a categorization into systematically varying and rigid parameters. Whereas jaggedness/conspicuous error bars are generally associated with rigid parameters in the INA representation, this feature in the curves of Fig. 9 cannot be associated with specific structural elements. Thus, the INA parameters are superior both for structural prediction and in providing information on structural adaptation to varying (p–T–X) conditions.

5.2. Reasons for adopting a structural refinement

There are several reasons for adopting a structural refinement strategy. First, the INA method generates, for intermediate (p, T, X) values, two independent sets of values for the cell parameters, one from anionic parameters and the other from cationic parameters.¹ These two sets of values need to be brought into agreement with each other. Secondly, although linear interpolation is a simple and robust method of predicting INA parameter values at intermediate (p, T, X) conditions, the reliability of prediction can be improved by applying normative constraints as part of such a refinement. In this work, two types of constraint have been applied on an exploratory basis, one based on cell parameters, and the other on the volumes of cation coordination polyhedra [see equations (11) and (12) with respect to forsterite at variable temperatures]. Additional types of constraint are conceivable. For example, one might wish to fix the value of f_t at 0.433. Alternatively, one might wish to build in constraints on interaction lengths obtained from calculated lattice energies, or indeed incorporate lattice energy calculations into the refinement procedure directly. It is maintained that the Excel Solver Add-In is a simple, but powerful front-end with its GRG non-linear refinement algorithm.

Having made the decision for a structural refinement strategy, a further feature of the INA parameterization compared with conventional crystallographic parameters comes to light. INA has brought the cell parameters and the atomic coordinates into a single unified framework. This is not the case with the untransformed alternative, since changes in structure resulting from varying (p, T, X) conditions are manifested by simultaneous changes in both cell parameters and atomic coordinates. Thus, the frame of reference for the structures changes, as well as the structures within that frame of reference.

5.3. The exclusion of polynomial fitting as an option

The proposed approach of linear interpolation in order to provide starting parameter values for a subsequent structural refinement has several advantages over polynomial fitting. Apart from its simplicity, it can be applied to all datasets, irrespective of whether monotonic/smooth trends in parameter values arise, or not. It also allows an equal a priori weighting to be given to all experimentally determined data-points. This principle of equal weighting, which is considered to be important, cannot be upheld by a polynomial fitting method.

In order to illustrate this point, a polynomial fitting of the variation with temperature of one of the INA parameters, s_1,z, is carried out here, as an example. This parameter shows a clear ascending trend, although it is not monotonic at the high-temperature end (Fig. 10). The maximum polynomial order consistent with 12 experimental temperatures would be 11 [equation (14)].

Here i represents the index of an INA parameter, in this case s_1,z. x represents reduced temperature, i.e. (T  (K) - 293 )/ (1180 - 293), with a_i,j the polynomial coefficients to be determined by fitting.

Figure 10 Comparison of fitting to the 12 s1,z values (black crosses) by linear interpolation (red lines), a polynomial of order 11 (blue curve), a polynomial of order 10 (green curve) and a polynomial of order 9 (pink curve).

Whereas the linear interpolations necessarily link all experimental points, this does not apply to the three polynomial fits. Although the 11th-order polynomial is more successful than the two lower-order polynomials in passing closer to all experimental points, it suffers from an undershoot at low reduced temperature and is characterized by significant curvature, in particular in the reduced temperature range between ca 0.4 and 0.6, this having no physical basis. Although the undershoot and the accentuated curvature can be removed by reducing the polynomial order from 11 to 9, such a procedure is arbitrary. It is observed that the trajectories of the three polynomials smooth out the details of the experimental points at the high temperature end. Whilst these polynomial fits might be closer to physical reality than the linear interpolations between experimental points, the sole requirement in this work is to find reliable starting points for structural refinements. The normative constraints applied in the refinements serve to shift the parameter values away from any physically unrealistic starting points.

The difficulties identified for polynomial fitting to the s_1,z parameter would also apply to the other INA parameters. In the case of non-smooth INA parameters in Fig. 8, such as t₂ and e_z, the best option would probably be to fit straight lines. The option of direct polynomial fitting to atomic parameters, i.e. without carrying out an INA transformation, would also be problematical: the jaggedness of many of the curves in Fig. 9 does not lend itself to polynomial representation. It is for this reason that methods such as parametric Rietveld refinement, as described by Stinton & Evans (2007), are based on the selective use of polynomials. In their approach, a continuous variation of parameter values with (p–T–X), as commonly described by polynomials, can only be assumed for some of the crystallographic parameters.

The idea of using a single polynomial for each INA parameter would also covertly assume that the responses of the parameters to the varying (p–T–X) conditions can be defined for the whole range of p, T or X. Although more than one polynomial per parameter could be employed to avoid this difficulty, each with its own (p–T–X)-range, this practice would introduce an unwelcome level of complexity into the method.

5.4. The context of the current work

In this initial work, the principles of INA have been described by exclusive reference to olivines. Therefore, a generalization to all minerals is called for. The two-dimensional nature of the cationic network (Fig. 5) and the z = 0 INA parameters (Fig. 1) is a peculiarity of olivines resulting from the dominant structural influence of the mirror planes in space group Pbnm. A three-dimensional cationic network is anticipated in the majority of mineral structures, thereby imposing tighter constraints on the development of INA parameter values under varying (p–T–X) conditions.

The use of piecewise linear interpolations, as described here, allows data with relatively high standard errors in atomic coordinates to be analysed. The observation of questionable irregularities in the variation of INA parameters with p, T or X could be exploited by crystallographers in order to arrive at better structural models.

The observed correlations between the temperature variations of cationic INA parameters are of particular interest (Fig. 8a). A guiding hypothesis for oxide structures is that cation–cation non-bonded interactions are more important than oxygen–oxygen non-bonded interactions in determining the sequences of phase transitions (Keeffe & Hyde, 1981a,b). If this is correct, it follows that a parameterization of cationic networks, along the lines of INA, could provide a robust, but flexible quantitative basis for modelling and/or rationalizing oxide mineral phase sequences.

Footnotes

¹In the particular case of olivines, this generalization applies only to the a and b parameters, since the c parameter value is generated by the anionic network only.