Species Limit Expressions for Order-Disorder Solution Models
_____________________________________________________________________________________________

NOTE: an excellent (!) python script to convert endmember site fractions to Perple_X species 
limit expressions is available at https://github.com/bobmyhill/solution_limits
courtesy of Bob Myhill (Bristol). 
_____________________________________________________________________________________________

To compute the equilibrium speciation of order-disorder models (typically solution model types
6 and 8) efficiently, Perple_X requires limits on the physically realizable fractions of the 
ordered species. In this context, physically realizable means that the site fractions of every 
atom possible on a particular lattice site are greater than or equal to zero. The limits
are entered as a list within a section of the solution model demarked by "begin_limits" and 
"end_limits" keywords. 

The fully disordered speciation of an order-disorder solution in a c-component composition 
space with s-ordered species is uniquely defined by the site populations and stoichiometry
of the c+s endmembers chosen to represent the solution. The fractions of the species for the 
disordered composition will be designated p0(1..c+s), while the equilibrium speciation will
be designated p(1..c+s). The bulk composition of each of the s-ordered species can be 
described in terms a stoichiometric relationship between c-disordered species (disordered 
is a little misleading in this context because it may refer to species in which no 
configurational entropy is possible) of the form 

ordered_endmember(j) = sum(c(j,i) * disordered_endmember(i), i=1..c),  (1)

where ..._endmember(k) denotes the bulk composition of endmember k, and c(j,i) is the 
stoichiometric coefficient of the disordered_endmember(i) "in" the ordered_endmember(j). 
From (1), for the disordered species (1..c)

p(i) = p0(i) - sum[c(j,i)*[p(j,i) - p0(j,i)], j=c+1..c+s).  (2)

Given the known atomic site fractions of every endmember species, the site fraction of 
atom k on site l can be expressed 

z(k,l) = sum(p(i)*z(k,l,i), i=1..c+s),    (3)

where z(k,l,i) is the site fraction of atom k on site l of endmember i. Substituting
(2) into (3) and collecting the stoichiometric factors c(j,i) and z(k,l,i) on a 
particular fraction into the coefficient d(k,l,i) yields

z(k,l) = sum(d(k,l,i)*p0(i), i=1..c)  + sum(d(k,l,i)*p0(i), i=c+1..c+s) 
                                      + sum(z(k,l,i)*p(i), i=c+1..c+s)     (4)

Solving z(k,l) = 0 and z(k,l) = 1 for p(c+1) gives upper and lower, not necessarily 
respectively, physically realizable bounds on p(c+1). In Perple_X type 6 and 8 solution
models, the lower bound on the ordered endmember p(c+1) from atom k on site l is written

p(c+1) = f(k,l,c+1) 
         - d(k,l,1)/z(k,l,c+1)*p0(1) ....  - d(k,l,c)/z(k,l,c+1)*p0(c) 
         - d(k,l,c+1)/z(k,l,c+1) * p0(c+1)
         - d(k,l,c+2)/z(k,l,c+1) * p(c+2)  - d(k,l,c+2)/z(k,l,c+1) * p0(c+2)
         ....
         - d(k,l,c+s)/z(k,l,c+1) * p(c+s)  - d(k,l,c+s)/z(k,l,c+1) * p0(c+s),   (5)

where f(k,l,c+1) is a collection of stoichiometric constants, it being assumed that: 

1) if the endmember is disordered (i < c + 1), the endmember name represents p0(i).
2) if the expression is for ordered species m, its name occurs only once and represents
   p0(m)
3) if the endmember is ordered (i > c) and not the mth ordered species, its name occurs 
   twice, the first occurence corresponds to p(i) and the second to p0(i).

Because the site fractions are linear in p(i..c+s) the upper and lower abounds differ 
by a constant, delta(k,l). Only the expression for the lower bound and the value
for delta(k,l) are specified in Perple_X solution model input. For example for the
Chl(W) (chlorite) solution model with the following endmembers and site occupancies:

                  M1      M2+M3     M4     T2
                 ________________________________

  Mutliplicity    1         4       1      2 
                 ________________________________
  
    daph          Fe        Fe      Al     Al_Si   disordered
    f3clin        Mg        Mg      Fe3+   Al_Si   disordered   
    ames          Al        Mg      Al     Al_Al   disordered
    afchl         Mg        Mg      Mg     Si_Si   disordered
    clin          Mg        Mg      Al     Al_Si   disordered 
    mnchl         Mn        Mn      Al     Al_Si   disordered   
                  _______________________________ 

    och1          Mg        Fe      Fe     Si_Si   ordered
    och2          Fe        Mg      Mg     Si_Si   ordered


The site fraction of Mg on lattice site M1 is:

         z(Mg,M1) = p(f3clin)  + p(afchl)  + p(clin) + p(och1) 
                  = p0(f3clin) + p0(afchl) + p0(clin)+ p(och1) - 4/5*p(och2) + 4/5*p0(och2) 

and solving Z(Mg,M1) = 1 for p(och2) yields a lower bound of p(och2) as

         p(och2) = -5/4 + 5/4*p0(f3clin) + 5/4*p0(afchl) + 5/4*p0(clin) + 5/4*p(och1) + p0(och2)

and in Perple_X input this bound (and the corresponding upper bound) is written 

 och2 = -5/4 + 5/4 afchl  + 5/4 f3clin + 5/4 clin + 5/4 och1 + 1 och2       delta = 5/4

For additional information on the formulation of order-disorder/speciation/compound formation 
models refer to chapter 11 of http://www.perplex.ethz.ch/thermo_course/thermo_course.pdf