Notes on the creation and manipulation of solid solution models

A large class of solid solution models are formulated on the premise that exchange of chemical species takes place on a finite number of unique crystallographic sites, and that the thermodynamic properties of the solution are a function of the proportions of species occupying each of the sites. These models are broadly classified as being of Bragg–Williams-type. They form an excellent first order approximation of non-ideal mixing and long-range order. In this article we present the mathematical framework common to all Bragg–Williams models, introducing necessary concepts from geometry, set theory and linear algebra. We combine this with a set of mathematical tools which we have found useful in building and using such models. We include several worked examples to illustrate key concepts and provide general expressions which can be used for all models. This paper is split into two parts. In the first part, we show how the valences of the species occupying each site and the total charge of the species involved in site exchange are sufficient to define the space of valid site occupancies of a solid solution, and to compute the endmembers bounding that space. We show that this space can be visualised as a polytope, i.e, an n-dimensional polyhedron, and we describe the relationship between site-occupancy space and compositional space. In the second part of the paper, we present the linear algebra required to transform descriptions of modified van Laar and subregular solution models from one independent endmember basis to another. The same algebra can also be used to derive macroscopic endmember interactions from microscopic site interactions. This algebra is useful both in the initial design of solution models, and when performing thermodynamic calculations in restricted chemical subsystems. A polytope description of solid solutions is used in the thermodynamic software packages Perple_X and burnman. The algorithms described in this paper are made available as python code.


Introduction
Solid solution models are an essential component in modelling geological, metallurgical and other chemical processes. The motivation behind this paper is to review some mathematical background defining the structure and properties of commonly used solution models, and to introduce tools that aid in their construction and manipulation. We concentrate on solution models which have a constant number of sites per formula unit, and where interactions are dependent only on the total proportions of species occupying each site. This type of solution was pioneered by Bragg and Williams as a way to model long range ordering in the 1930s braggspswil-liamssps1934 Williams 1934, 1935;Williams 1935). Such models do not consider local interactions, such as bonding between pairs or clusters of species (Bethe 1935;Inden et al. 2001;Kikuchi and Masuda-Jindo 2002).
Many readers will be familiar with the key petrological concepts described in this paper. After all, the importance and consequences of solution chemistry and order-disorder reactions on mineral properties were already understood in the 1970s (Grover and Orville 1969;Thompson and James 1969;Wood and Nicholls 1978). However, mathematical descriptions of the models underlying the petrological concepts are usually explained briefly, or for only specific cases. In this paper, we provide a more complete description. We particularly emphasise the difference between

A geometric visualization of solid solutions
The solid solution constraints given by Eqs. 2-4 are geometrically equivalent to a class of convex polytopes (the n-dimensional equivalent of a polyhedron). Every vertex of the polytope corresponds to an endmember of the solution with a different distribution of site-species.
The geometry of the polytope corresponding to a particular solution depends only on the number of sites, the number of site-species on each site, and the charge balance constraint.

Single-site solutions
Single-site solutions can be represented as simplexes, which are the n-dimensional equivalent of a triangle. A solution with two site-species can be represented by a line (a 1-simplex), a solution with three site-species by a triangle (a 2-simplex), and a solution with four site-species by a tetrahedron (a 3-simplex). Examples are shown in Fig. 1a, b. (2) [A] [C] [D] [A] [B] [C] a b c [A 0.5 (Fig. 1c).

Multisite solutions
Multisite solutions are geometrically equivalent to the Cartesian product of the individual site-simplexes. This mathematical jargon encapsulates the idea that there is an endmember vertex for every possible combination of site-species. Endmember vertices are connected by an edge if they differ by only one site-species. In Fig. 2a (Grover and Orville 1969) and [Cu,Ag] 10 [Fe,Zn] 2 [Sb,As] 4S13 fahlore Sack (2017) respectively.
As with single site solutions, if charge balance is not automatically satisfied by the site constraints, the additional constraint is equivalent to cutting the polytope with a hyperplane. For example, the triangle outlined by black lines in Fig. 2c Fig. 2 Polytopes corresponding to some 2-site (a, c) and 3-site (b, d) solution models. a The Cartesian product of a 2-simplex and 1-simplex (a triangular prism), such as [Ca,Fe,Mg][Fe,Mg]Si 2 O 6 -clinopyroxene. b The Cartesian product of three 1-simplexes, e.g. [Cu,Ag]  Some solution model formalisms such as the compound energy formalism (Hillert 2001) express the excess energy of a solid solution in terms involving all of its endmembers, whether or not they are neutrally charged. Others express the excess energy in terms of the energies of an independent set of endmembers (Helffrich and Wood 1989;Holland and Powell 2003), where an "independent" set comprises the minimum number of endmembers that can span the entire site-species occupancy space of the system. For example, for the halide solution [Na, K] written in terms of three endmembers. We shall return to the relationship between independent and dependent endmembers, and to different formalisms for the excess energies of mixing in the following sections.

An algebraic description of solid solutions
The graphical representations of the previous section are a useful introduction to the site-species occupancy spaces of solution models, but it is also necessary to have a mathematical description that can be used in computations. In the following sections we shall outline just such a description, using vectors and matrices which are summarised in Table 1.
First, we formalise the constraints in Eqs. 2-4 using set notation. The set X of valid site-species occupancies x for any solution can be described using an n constraints ⋅ n site-species matrix P and a vector b of length n constraints : Each row of the matrix P and element of vector b corresponds to an independent equality constraint on the site populations.
We now present two examples of P and b . First, consider the classic two site clinopyroxene in the CFMS system, given by the formula (Grover and Orville 1969): Each of the two sites must be fully occupied, and the solution must be charge-balanced, which gives three equality constraints, which are tabulated in Table 2. Inspection of this table reveals that the charge balance constraint is a linear combination of the site constraints (2 ⋅ [row1] + 2 ⋅ [row2] = [row3]) and it is therefore not independent from the site constraints.
Our second example is a one-site pyrope-majorite garnet: Si 2 O 6 .

Set of isochemical reactions between independent endmembers
Null(S ind ) q n Vector of isochemical reaction amounts Here we define an ordered endmember as any instance of a solid solution where each site is occupied by a single species. This differs from other studies, which define ordered compounds as those where some species reside on only a subset of their potential sites (Holland and Powell 1996).
The equality constraints for this solution are given in Table 3. Unlike the clinopyroxene example, the charge balance constraint is not a linear sum of the site constraints, and so it represents an independent constraint on the system. Once the matrix P and vector b for a solid solution have been specified, standard algorithms can be used to determine the characteristics of the solution polytope which we represented graphically in "A geometric visualization of solid solutions". Specifically, we can: -Determine the number of independent endmembers. -List the site-occupancies of all the endmembers.
-Determine an independent set of endmembers.
-Create a set of equality constraints ( P and b ) from an independent set of endmembers.

Determining the number of independent endmembers of a solution
Firstly, the number of independent endmembers can be obtained from the number of rows and columns of P . For a linear system with a certain number of linearly-independent constraints ( n constraints ) and unknowns n unknowns , there are n unknowns − n constraints degrees of freedom: In this specific case, n unknowns is the number of site-species, and n constraints is the number of sites, plus one if the chargebalance constraint is linearly independent from the site constraints. Each degree of freedom can be considered to correspond to the proportions of an independent endmember. The total proportion of independent endmembers must sum to one, so there is one fewer degree of freedom than the number of independent endmembers. We can then write the following expression by substitution into Eq. 6: where c = 1 if charge-balance constraint is linearly independent from the site constraints, and c = 0 otherwise. Therefore, the CFMS clinopyroxene example ( (6) n dof = n unknowns − n constraints .

Listing endmember site-occupancies
The site-occupancies of all the endmembers of a solution which do not have an independent charge-balance constraint ( c = 0 ) can be tabulated by iterating through all the possible combinations of site-species. The resultant list is known mathematically as a cartesian product. The total number of endmembers in such a solution is given by: where n species,S is the number of species on Site S. For example, [Ca, Fe, Mg][Fe, Mg]Si 2 O 6 clinopyroxene ( Table 2) has 3 ⋅ 2 = 6 endmembers, which can be tabulated in an endmember site-occupancy matrix (Table 4).
In solutions with charge-balance constraints, listing all of the valid endmembers is complicated by the fact that the charge-balance constraint can create disordered endmembers. For example, the one site pyrope-majorite (Table 3) has two valid endmembers, one of which is a disordered majorite, with both Mg and Si residing on the same site (Table 5).
In mathematical parlance, finding all the possible endmembers of a solution (determining E ) is known as vertexenumeration, and involves not only counting the vertices, but also determining their positions in space (Matheiss and Rubin 1980;Avis and Fukuda 1992;Lasserre 2004). In the python programs accompanying this paper, we generate   Table 3 Equality constraints for a one-site pyrope-majorite garnet (of multiplicity 2), tabulated as a matrix P and vector b (Eq. 5) The charge balance constraint is not independent from the site constraint matrices E from P and b (e.g. Tables 2, 3) using the package pycddlib (Fukuda and Prodon 1995). This package implements the double description method of Motzkin et al. (1953). This method takes as input a set of inequality constraints (equalities are represented as two inequality constraints), and uses these to compute all of the vertices of the polytope which is bounded by those inequalities. Mathematically-inclined readers are referred to the original paper for more information.
The solution polytope implementation in burnman (Appendix A) includes a function which generates a polytope from charge balance equalities. This function can be called in a single line, and generates an object with attributes such as the matrix E.

Finding an independent basis set of endmembers
Once all the endmembers of a solution have been obtained, a set of independent endmembers and their site-occupancies E ind can be computed in several ways. One way is to compute the row-reduced en-echelon form of E . Row-reduction (also known as Gaussian elimination) reduces a matrix to upper triangular (en-echelon) form through a combination of row swapping, scalar multiplication of rows, and addition/ subtraction of two rows. It is a standard technique to solve systems of linear equations. The endmembers corresponding to the non-zero rows of the en-echelon form of the matrix constitute an independent set. An implementation of this technique is provided in the python programs provided with this paper (Appendix A). The independent endmember set determined in this way depends on the initial ordering of the rows in E and the exact algorithm chosen for the row reduction (see Supplementary Information).
The relationship between independent endmember proportions p ind and site-species occupancies x is given by the equation: In multisite solutions, valid site-species distributions can be obtained even if some elements of p ind are negative.

Creating a set of equality constraints from an independent set of endmembers
Creators of thermodynamic models often start from a preferred basis set of independent endmembers for certain solid solutions. If one wishes to list all of the endmembers that can be described using this independent basis, they must first convert the basis into a set of equality and inequality constraints as given in Eq. 5 and then run those constraints through a vertex-enumeration routine. To see how to do this, let us consider a two-site bridgmanite solution, defined by the independent endmembers given in Table 6.
We are looking for a set of independent equality constraints which are satisfied by all possible instances of this solution. A partial set is provided by the right nullspace (also known as the kernel) of E ind . The right nullspace N(E ind ) corresponds to the set of changes in site-species occupancies which cannot be achieved by changing the amounts of independent endmembers; i.e. the site-species occupancies of the solution must satisfy the expression One potential independent endmember set is denoted by asterisks Pyrope py 0 1 0 Disordered majorite dmaj 0.5 0 0.5 Table 6 The endmember site-occupancy matrix E for a two-site FMAS bridgmanite For the bridgmanite solution given in Table 6, two potential basis vectors for the right nullspace are [ 0, 0, −1, 1, 0 ] and [ −1, −1, 0, 0, 1 ]. These two vectors define compositions such that there is an equal amount of Al on each site and that the total Si on the B site is equal to the sum of Mg and Fe on the A site. In practise, the construction of the nullspace can be achieved by row reduction of E ind followed by back-substitution (see supplementary information for a worked example). An additional equality constraint fixes the total number of moles of the solution to be equal to one: As we did in "Site-species and solution constraints", we can now tabulate these equality constraints. For the bridgmanite example, this construction leads to matrix P and b as given in Table 7.
We now check that the equality constraints represented by this table are equivalent to those produced by the procedure in "Site-species and solution constraints". The endmembers in Table 6 Table 8), and therefore they represent the same equality constraints.

The composition space of a solution model
"Simple" solution models are defined as models where every point in site-occupancy space corresponds to a unique bulk composition. This is true of [Na + , K + ][Cl − , F − , I − ] halide, for example. Not all solution models are "simple"; in many models it is possible to redistribute site-species over sites without changing the bulk composition. In such solutions, which we call order-disorder solutions (also sometimes called "complex" solutions), the composition and site-occupancy spaces of the solution are not equivalent. For example, the following two instances of CFMS clinopyroxene have exactly the same bulk composition: Site-species redistributions which satisfy the solution model constraints can be described by one or more isochemical reactions, bounded by the positivity constraints for each individual site-species (Eq. 2). A basis set of independent reactions can be found by first constructing the matrix of elemental compositions of each of the independent endmembers. We call this the stoichiometric matrix S ind . The nullspace of S ind corresponds to the set of isochemical reactions. We define R ind as any independent basis for this nullspace. The matrices S ind and R ind for CFMS clinopyroxene are given in Table 9. The reaction vector R 1 in Table 9 indicates that the bulk composition of the solid solution will not change if we substitute 2 moles of diopside and 1 mole of clinoferrosilite for 2 moles of hedenbergite and 1 mole of clinoenstatite. To see the redistribution of site-species implied by this reaction, we can take the dot product of ( Many petrological studies choose to consider the compositional space of the solid solution, rather than the siteoccupancy space, expressing the state of order via one or more order-parameters. Within the framework we have described, the compositional space can be obtained by projecting the site-species occupancy polytope onto a (hyper) plane perpendicular to the ordering vector(s). We illustrate this graphically for CFMS clinopyroxene in Fig. 3 site-occupancy space of this solution is a triangular prism (Fig. 3a), and the isochemical ordering vector [Fe −1 Mg 1 ] [Fe 1 Mg −1 ] is collinear with the cmf-cfm vector. Planes perpendicular to the ordering vector satisfy expressions of the form where q 1 is a scalar corresponding to the distance along the ordering vector. If we project the polytope onto any one of these planes, we obtain the classic "CFMS quadrilateral", which defines the composition space of the solid solution (Fig. 3b).
Because the isochemical ordering vector R ′ 1 is a linear combination of the independent endmembers di, hed, cen and cfs, it is possible to replace any one of those endmembers with the vector to make a new independent basis set (e.g. di, hed, cen, R ′ 1 ). If this is done, then one can describe any instance of the solid solution in terms of compositionally-independent endmember amounts and order parameter(s) (e.g., p di , p hed , p cen , q 1 ). The composition of the solution is determined by the endmember amounts, which must sum to one. However, the site occupancy space is only uniquely defined using all four parameters. For any bulk composition, the stable distribution of site-species is found by minimizing the Gibbs energy over all valid values of q 1 .
Within any order-disorder solution polytope, there is one surface of particular interest: the surface which maximises the configurational entropy q(p) . This surface is of interest for two reasons: it is continuous, differentiable and valid over all of composition space, and therefore represents a suitable starting guess for finding the stable configuration of site-species at fixed composition; and in the disordered (high-temperature) limit, this surface represents the stable configuration of site-species at all bulk compositions. In our clinopyroxene example, the maximum entropy surface Table 9 Composition and siteoccupancy matrices for two-site CFMS clinopyroxene First four lines: stoichiometric matrix ind , reaction basis vector R ind and independent endmember site occupancy matrix E ind . Last line: the isochemical reaction vector corresponding to R 1 written in terms of composition and site-species occupancies . 3 a The CFMS clinopyroxene polytope. There are six endmembers of this polytope: diopside (di), hedenbergite (hed), clinoenstatite (cen) clinoferrosilite (cfs) and two ordered endmembers, (cfm and cmf). The purple ordering vector represents the isochemical exchange reaction Fe B + Mg A ⟶ Fe A + Mg B . The blue plane is perpendicular to the ordering vector, and points on that plane satisfy the expres- corresponds to points within the solution space where the Fe:Mg ratio is equal on both sites. The values of q 1 corresponding to the maximum entropy surface are shown as contours in Fig. 3b. This surface is given by the formula q 1 = 2n Ca (n Fe − n Mg )∕((n Fe + n Mg )(n Fe + n Mg + n Ca )) , where n i is the number of atoms of element i in the solution.

An illustrative example of the energetics of order-disorder in Bragg-Williams-type solid solutions
In order-disorder solution models, the stable arrangement of site species at any given bulk composition is determined by minimizing the Gibbs energy with respect to the isochemical reaction vectors. Order-disorder solutions can be assigned to one of three groups: (a) those that are completely disordered at all temperatures, (b) those that become completely disordered at high temperatures (convergent ordering), and (c) those that approach, but do not reach, a state of complete disorder at high temperatures (non-convergent ordering) (Thompson and James 1969 where the Gibbs energy G of mixing is represented as a regular solution (Wohl 1946;Wood and Nicholls 1978;Powell and Holland 1993) between the endmembers clinoenstatite (cen), clinoferrosilite (cfs) and ordered Fe-Mg clinopyroxene (cfm). This energy of mixing is split into a non-configurational part ( G * ) and a configurational part ( TS conf ): T is the temperature, S conf is the configurational entropy and R is the gas constant. The expression for S conf ignores any contribution from short-range ordering. G * depends on interaction energies between endmembers ( W ij ) and the Gibbs energy required to form 1 mole of endmember cfm from a mechanical mixture of clinoenstatite and clinoferrosilite ( 0.5cen + 0.5cfs ⟶ cfm;G * cfm ): G * cfm accounts for the volume and vibrational entropy change of the reaction, but does not account for the change in configurational entropy, which in this case is zero because the cen, cfs and cfm endmembers all have only one species occupying each site.
Parameters which result in disordered, convergent and non-convergent models are provided in Table 10, and the corresponding non-configurational Gibbs energy surfaces are shown in Fig. 4. In the disordered model (Fig. 4a), the ordered endmembers cfm and cmf have higher non-configurational Gibbs energies than their disordered counterparts. The symmetry of the nonconfigurational Gibbs energy surface and the configurational entropy about the line of complete disorder means that the solution remains perfectly disordered at all temperatures. Maintaining the symmetry about the line of complete disorder, but making G * cfm negative means that ordering is favoured at low temperatures. As temperature increases, the configurational entropy term in Eq. 12 causes progressive disordering; the dashed and dotted lines in Fig. 4b correspond to the equilibrium distributions of species at 400 and 500 K respectively. At high temperatures, the configurational entropy term dominates, resulting in complete disorder.
In reality, the two sites in clinopyroxene are not identical (Grover and Orville 1969;Holland et al. 2018), and therefore the excess energy is not symmetric about the line of disorder (Fig. 4c). The degree of disorder increases with increasing temperature, but partial ordering remains at all finite temperatures.
In the examples given here, the progress from order to disorder is monotonic with temperature at all bulk compositions. Furthermore, the value of the order vector [Fe −1 Mg 1 ][Fe 1 Mg −1 ] corresponding to the most stable arrangement of site-species always lies on the same side of the maximum entropy surface (the most stable configuration of species always lies in the upper left triangles in Fig. 4, regardless of bulk composition). Not all solutions behave in this way. For ordering reactions in which G * is a function of pressure and temperature, changes in conditions can potentially lead to inversions in the sign of the ordering vector (Redfern et al. 2000), and large non-ideal interaction energies can lead to islands of stability on both sides of the line of complete disorder/maximum entropy.

Constructing thermodynamically consistent solution models
If a solid solution model has fewer independent endmembers than required by Eq. 7, then there is at least one constraint that is not related to charge-balance. Sometimes, these choices are made consciously, with the aim of improving computational efficiency by reducing the number of order parameters. However, reducing the number of endmembers also restricts the valid site-occupancy space of the solution, that in turn will tend to artificially destabilise the solution. We advocate always using the full site-occupancy space.
As an illustration of how difficult it can be to assess the validity of order-disorder solution models without using the mathematical tools described here, we present an independent endmember basis set for clinoamphibole in the NCKF-MASHTO system. This solution includes 18 site-species distributed over six sites (Green et al. 2016;Holland et al. 2018), with the general formula Using Eq. 7, a complete model in this system has 12 independent endmembers (Table 11, computed using the code accompanying this paper). The published model in this system (Green et al. 2016) includes only the first 11 endmembers in this set. Without using Eq. 7, it would be extremely difficult to know that the published solution model was incomplete. Even armed with the tools in this paper, it is difficult to quantify the effect that completing the basis set of endmembers has on the properties of the solution. Certainly, the smaller model fails to span the entirety of site-occupancy space; the full model has 436 dependent endmembers, while the published model has only 156. One difference is that in the subsystem the full model allows 0 ≤ x ≤ 1 , whereas the partial model only allows x = 0.5 . The partial model also excludes hydrogen-free, highly aluminous endmembers and hydrogen-free, ferric-iron bearing endmembers.
We mention this clinoamphibole model because it is particularly complex, but is hardly unique among published models. The KFMASHTO biotite model published in Tajčmanová et al. (2009) is another solution that contains one fewer independent endmember than the number suggested by Eq. 7. The general formula is whose site occupancy space can be spanned by seven endmembers (Table 12).
As with the clinoamphibole example, the published biotite model couples titanium substitution with deficits in hydrogen content. In the subsystem the full model allows 0 ≤ x ≤ 1 , whereas the partial model only accepts x = 0.5 . The partial model also excludes all 12 titanium-bearing endmembers in the subsystem  Site multiplicities are given in brackets. Endmember abbreviations are-tr: tremolite, ts: tschermaks-substituted camph, parg: pargasite, gl: glaucophane, cumm: cummingtonite, grun: gruenerite, a/b: two ordered Mg-Fe-bearing endmembers on the cummingtonite-gruenerite compositional join, mrb: Mg-riebeckite, kprg: K-pargasite, tts: Ti-endmember. The famph endmember is one of many potential endmembers that completes the published endmember basis (Green et al. 2016;Holland et al. 2018). Only the mixing sites are shown

Changing independent endmember bases
In this second part of this paper, we show how to convert endmember and interaction energies between sets (or bases) of independent endmembers. This conversion has two primary purposes: -We often wish to solve thermodynamic problems in restricted compositional spaces. of A ij corresponds to the number of moles of original endmember i contained within the new endmember j. The proportions of the original endmember set p i in terms of the new endmember set p ′ l are thus given by: The following subsections outline the mathematics to convert endmember bases within the "modified van Laar" (Holland and Powell 2003) and "subregular" (Helffrich and Wood 1989) mixing model formulations. "Ideal", "symmetric/regular" solution models can be viewed as special cases of these more general formalisms.

The modified van Laar model
The van Laar model (van Laar 1906) was reformulated by Holland and Powell (2003) for use as a generalised macroscopic asymmetric model. The general idea behind this formalism is that independent endmembers are associated with a parameter that skews non-ideal energies toward or away from the endmember. Greater values skew the energies increasingly toward the endmember. The non-configurational Gibbs energy at any given distribution of site species is a function of the proportions of the independent set of endmembers p i : where G mbr i is the Gibbs energy of endmember i, i is the van Laar (asymmetry) parameter for that endmember, the components of W binary ij are equal to 2W ij ∕( i + j ) , and W is an upper triangular matrix containing the binary endmember interaction parameters. Site multiplicities are given in brackets. Endmember abbreviations are-fbi: ferribiotite, tbi: Ti-bearing biotite, east: eastonite, ann: annite, phl: phlogopite, obi: ordered biotite on the annite-phlogopite compositional join. The ffbi endmember is one of many potential endmembers that completes the published endmember basis (Tajčmanová et al. 2009  Changing the set of independent endmembers for the asymmetric model has been described in Diener et al. (2007). Here, we provide an alternative derivation using Einstein summation convention. Repeated indices are summed over unless they appear on both sides of the equation.
To change the set of independent endmembers, first combine the endmember and interaction terms into a single matrix W C : This equation can be transformed into an expression involving the new endmember set by substituting p with p ′ (Eq. 15): We now seek to express this equation in the same form as Eq. 16. First, define new asymmetry parameters ′ and matrices B and C: such that Substituting these expressions into Eq. 21 yields The interaction matrix can now be transformed to yield an expression in the form of Eq. 16: The transformed endmember properties are now removed from W ′C . This is the reverse of the operation in Eq. 18: Finally, the binary can be converted back to upper triangular form by summing the upper and lower triangular components of D: where the components of W ′binary lm are equal to 2W � lm ∕( � l + � m ) .

The subregular model
The subregular model (Andersen and Lindsley 1981) is another popular form of asymmetric solution model. Instead of defining asymmetries with a small number of parameters assigned to the endmembers, the subregular model is simply an extension of the regular solution model to cubic terms in the independent endmembers.
In the subregular model, the non-configurational Gibbs energy includes unary G mbr , binary W binary and ternary W ternary contributions (Helffrich and Wood 1989 To convert Eq. 32 into an expression involving a new independent set of endmembers, we first combine the various parameters into a single term involving a interaction matrix W C : The conversion is non-unique; one way to construct W C ijk is as follows: Equation 33 can be transformed into an expression involving the new endmember proportions using Eq. 15: To convert this back into the form of Eq. 32, we first remove the transformed endmember excesses from W ′C : The transformed binary matrix W ′binary is then given by Removing the contribution of W ′binary from B leaves us with matrix C: The transformed ternary components W ′ternary can then be found by summing the six contributing terms in C: In multi-site subregular models, asymmetric and ternary terms are both contributors to asymmetry in the energetics of mixing. Moreover, changing the independent endmember set can result in the appearance of non-zero ternary terms where none were present in the original set of interactions.
The simplest example where this is the case is a two site model [A, B][C, D]. Let us arbitrarily consider AC, BC and BD to be the original endmember set, and replace BD with AD in the new endmember set, such that: Let us say that we have information only on the AC-BC and BC-BD binaries, and that the AC-BC binary appears (41) to be nearly ideal, while the BC-BD binary appears moderately asymmetric. We make the reasonable assumption that the two sites do not interact, so that W binary AC,BD = W binary BC,BD and W binary BD,AC = W binary BD,BC . We naively assume that we can ignore ternary terms, and parameterise the model as follows: Transforming into our new endmember set [AC, BC, AD] yields the following The emergence of a non-zero ternary term in the transformed solution implies that there was no real justification for setting the ternary term equal to zero in the original formulation. Experimental and theoretical arguments for incorporating non-zero ternary terms have been discussed previously (Helffrich and Wood 1989;Cheng and Ganguly 1994), and our observation reinforces these arguments. However, to our knowledge, the implications of imposing a value of zero on ternary terms in multisite solid solutions has not been discussed. We can visualise the implications for our simple model by contouring the excess energy as a function of composition (Fig. 5). In the left hand panel of Fig. 5, we show the excess energy of mixing from the subregular model described above. As expected, the binary along the bottom (BC-AC) is ideal, and the binary along the left hand side (BC-BD) is asymmetric, with the energy skewed toward (43) 5 Excess energy of mixing of the subregular solution models as described in "The subregular model". Parameters for the left hand figure are as given in Eq. 43. The model parameters for the right hand figure are the same, except that the ternary term W ternary AC,BC,BD has been changed from 0 to 2 kJ/ mol the BC endmember. However, the system as a whole is not symmetric. This might be surprising, given the assumptions we put into the model. The right hand panel of Fig. 5 illustrates the effect of increasing the ternary term W ternary AC,BC,BD from to 2 kJ/mol. This modification yields what we might have expected from the original model; mixing on the first site is ideal, while mixing on the second site is non-ideal and mildly asymmetric.
For models such as this, it may not seem obvious how to choose appropriate ternary terms given the microscopicallymotivated model assumptions. The solution to this problem is to express the Gibbs energy of mixing in terms of microscopic interactions, rather than interactions between endmembers. For the example above, we set the microscopic interaction parameters w binary CD,Site 2 = 2 kJ/mol and w binary DC,Site 2 = 4 kJ/mol, with all other interactions set to zero. The microscopic interaction parameters can then be converted into their macroscopic equivalents using the conversion procedure described in the following section.

Converting microscopic interactions into endmember interactions
The models presented so far are macroscopic models; they deal with the energetics of interaction between endmembers. Solution models can also be formulated in terms of atomic interactions (Sack andGhiorso 1991, 1994;Pingfang et al. 1994). Microscopic models can be described by interaction matrices of the same form as their macroscopic counterparts, with the dimensions of the interaction matrices equal to the number of site-species, rather than the number of independent endmembers. The elements of the matrices correspond to the site-species interactions.
There are benefits to describing the properties of a solid solution in terms of microscopic interactions. Microscopic descriptions provide a much more direct link to the physics of the interactions (for example, order-disorder, Al-Al avoidance). They therefore neatly sidestep the problems encountered in the previous section. In addition, Powell et al. (2014) argue that well-constrained values of microscopic interactions in one mineral can be used as informed guesses for other minerals. For example, Fe-Mg exchange interaction has a value of around 5 kJ per mole of sites in many silicate minerals. This concept, termed micro- , may be useful for constructing complicated models where there is insufficient data to fully constrain all interactions.
Despite the benefits of parameterising solutions using microscopic models, it is often practical to convert these to macroscopic models (Powell and Holland 1993). Macroscopic models require fewer parameters, automatically maintain charge balance, and the set of independent endmembers can be chosen to improve computational efficiency. The linear relationship between site-species proportions and endmember proportions (Eq. 9) has the same form as the relationship between sets of independent endmembers (Eq. 15). Therefore, the conversion from a microscopic to a macroscopic formalism requires only the mathematics of "Changing independent endmember bases". This is true for both the asymmetric and subregular models, although previous work on micro-focused only on symmetric interactions ).
There are a couple of substitutions that need to be made to the mathematics of "Changing independent endmember bases" to convert microscopic interactions to macroscopic interactions. First, the untransformed endmember proportions p in Eq. 15 must be replaced by site-species proportions x , and the matrix A must be replaced by E indT . The endmember energies G ′mbr and (macroscopic) interaction matrices calculated using the asymmetric formalism equations of "The modified van Laar model" must be normalised to one mole of endmembers by multiplying them by n sites . Finally, the subregular equations require only that the 1 i W binary jk term in Eq. 34 be divided by n sites . The following subsections describe how to populate the binary and ternary terms in the microscopic interaction matrices. For worked examples, the reader is referred to the microphi package (see Appendix A).

Binary interaction parameters
The symmetric, modified van Laar and subregular formalisms all involve parameters describing the interactions between pairs of site-species. Two types of interactions are considered: simple mixing on a single site (e.g. Mg 2+ and Al 3+ on Site X), and two-site combinations of species (e.g. Mg 2+ on Site X and Ca 2+ on Site Y). Each entry in the binary interaction matrix corresponds to a Gibbs energy of formation of a cluster of sites from their constituents: For single-site mixing, the reaction can be rewritten as Two-site (XY) energies ( [A] X [C] Y ) cannot be uniquely determined from experimental analyses. Instead, the energies corresponding to cross-site reactions are used to populate the matrix. For example, the reaction is associated with the cross-site interaction energy w ACBD,XY , that is a function of the two-site energies (e.g. [A] X [C] Y , abbreviated as AC ): In each XY block of the microscopic interaction matrix, one component can be set to zero. For the example above, let us arbitrarily choose that element to be [B] X [D] Y . All the remaining elements (such as [A] X [C] Y ) of that block are filled with values corresponding to reaction with the excluded component −w ACBD,XY . Note that interaction energies involving repeated site-species (e.g. w BCBD,XY and w ADBD,XY ) are also set equal to zero, because the corresponding reactions have the same products and reactants (Powell and Holland 1993). A fully-worked and corrected symmetric example corresponding to the model discussed in Powell et al. (2014) is provided in the python package accompanying this paper (Appendix A), along with notes in Appendix B. A subregular example corresponding to the model described at the end of "The subregular model" (and in Fig. 5) is also provided.

Ternary interactions
The subregular model has ternary interaction parameters in addition to binary parameters. To construct a microscopic ternary matrix w ternary , we consider each component of the matrix to correspond to the energy of formation of a cluster of sites Mixing of species on the same site can be rewritten: When mixing three site-species on two distinct sites, we have The energy of formation of two-and three-site complexes includes site-bonding terms that cannot be uniquely determined by experimental means. Similar to the case of binary two-site exchange, we can arbitrarily select a single "special" component in each XXY or XYZ block of the ternary matrix to be equal to zero (e.g.
Finally, components sharing no site-species with the special component are assigned values based on exchange reactions involving five clusters: The number of non-zero (and independent) components in the XYZ blocks of the ternary matrix is therefore mno − m − n − o + 2 , where m, n and o are the number of potential site-species on the X, Y and Z sites. The XXY blocks have fewer independent endmembers (because , but are otherwise constructed in the same way.

A rationalization for using asymmetric interaction terms in solid solution models
At the microscopic scale, asymmetry in the energetics of mixing in solid solutions implies that the excess bonding energy associated with dissolution of a small amount of endmember B into endmember A is lower than the excess bonding energy associated with dissolution of a small amount of A into B. Such energetic contributions are inherently local in nature, and local interactions are not explicitly treated by Bragg-Williams-type models.
So, to what extent can asymmetric interaction parameters approximate energetic effects due to local interactions? While short-range effects cannot be modelled effectively by Bragg-Williams models, we know that they can approximate long-range order by splitting the site on which ordering takes place ("An illustrative example of the energetics of orderdisorder in Bragg-Williams-type solid solutions"). It turns out that a solution model that has >2 identical sites and symmetric interaction parameters can be reduced to a subregular model with a single site at high temperatures.
Take/consider as an example [Mg,Ca] 3 Al 2 Si 3 O 12 garnet, a solution that is believed to exhibit larger excess energies of mixing at Ca-rich compositions (Ganguly et al. 1996;White et al. 2014). If we split the Mg,Ca site into three sites, we can represent this solution with the general formula This solution has eight distinct endmembers (the solution polytope is a cube). One could create an expression for the non-configurational Gibbs energy of mixing for this solution that only involves symmetric terms: where S X is the set of species {Ca, Mg} on site X. Making the assumption that all three sites have identical mixing properties means that W CaMgMg = W MgCaMg = W MgMgCa and W MgCaCa = W CaMgCa = W CaCaMg . At high temperature, this solution will become completely disordered, such that the free energy can be described in terms of the bulk proportions of Ca and Mg: Note that this parameterisation is equivalent to a subregular mixing model; the model is asymmetric when W MgCaCa is not equal to W CaMgMg . Using a one-site asymmetric model is therefore equally as valid as using the three-site convergent order-disorder model, as long as we are only interested in phase relations at temperatures above the disappearance of long range order ("An illustrative example of the energetics of order-disorder in Bragg-Williams-type solid solutions").

Other solution model formalisms
Although we have restricted the discussion in this paper to the simplest of the asymmetric Bragg-Williams-type solution models, similar procedures can be applied to other formalisms. For example, the energetics of mixing in the compound energy formalism (Hillert 2001) can also be reformulated as a polynomial in site-occupancy space, and therefore expressed as a macroscopic model.

Appendix A: Research data
Python implementations of the polytope algorithms described in this paper are implemented in the open-source burnman software (Cottaar et al. 2014). The latest version can be downloaded from https:// github. com/ geody namics/ burnm an. Scripts used to create the tables and Fig. 4 in this (55) Mg paper can be found in the contrib/solution_polytope directory. Examples of solution model polytope generation, endmember determination, independent endmember (basis set) generation and gridding of site occupancy space can be found in the example script example_solution_creation_ and_manipulation.py. A separate project implementing the endmember and microscopic-macroscopic basis conversions for the subregular and modified van Laar formalisms (Helffrich and Wood 1989;Diener et al. 2007;Powell et al. 2014) can be downloaded from https:// github. com/ bobmy hill/ micro phi. A set of illustrative examples are provided.

Appendix B: A matrix formulation of the chlorite microphi example from Powell et al. (2014)
In the paper introducing the microphi technique to convert microscopic interactions into a macroscopic solution model, Powell et al. (2014) created an example of a four site, ten site-species chlorite model with general formula The independent endmember set that they used is given in Table 13. In that paper, the relationships between microscopic and macroscopic interactions were provided as individual equations. In this appendix, we provide the same interactions in a form which is compatible with the tensor expressions in "Changing independent endmember bases".
The chlorite solution proposed in Powell et al. (2014) is a regular/symmetric model, and so the microscopic interaction matrix can be represented as an upper triangular block matrix (Table 14). Each block contains microscopic interactions between pairs of species residing on specific sites. The blocks down the diagonal contain binary interactions between species on the same site, while the off-diagonal blocks contain cross-site binary interactions. As explained in "Binary interaction parameters", the elements of one row and one column of each cross-site block are set to zero, marked in Table 14 Table 14 corresponds to crosssite interactions between the M1 and M4 sites. The bottom right element of this block (the one labelled −w cro ) corresponds to the energy associated with the reaction where each of the first three terms represent one of the dashed elements. Non-zero interactions proposed in Powell et al. (2014) are assigned symbols in Table 14. The same-site interactions on the M-sites are w MgFe = 4 kJ/mol and w MgAl = 10 kJ/mol, and w FeAl is assumed to be related to w MgAl by a scaling factor ( w FeAl = w MgAl , = 0.7). The tetrahedral site interaction between Al and Si ( w AlSi ) is chosen to be equal to 0 kJ/mol. Finally, two cross-site terms are assigned non-zero parameters: w crt = w AlAlMgSi,MT2 = 14 kJ/mol, and w cro = w AlAlMgMg,M1M4 = -28 kJ/mol. There are several cross-site reactions involving Mg-Fe exchange, which were assumed to be negligible by Powell et al. (2014). These have been marked with "0" in Table 14.
Armed with the transformation matrix (Table 13) and the microscopic interaction matrix (Table 14), the steps described in Sections 3.2 can be used to derive Mg M1 Al M4 + Al M1 Mg M4 ⟶ Mg M1 Mg M4 + Al M1 Al M4 the macroscopic interaction matrix (Table 15). The interaction parameters agree with all the algebraic equations in Powell et al. (2014) (although the ames-(afchl,daph,ochl1,chl4) interaction parameters differ from the published table).
Funding This work was supported by STFC Grant No. ST/R001332/1.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.