research.amnh.org/~debel ... codes/B83a/B83acts1.txt

R. Berman (PhD, U. British Columbia, 1983) CaO-MgO-Al2O3-SiO2 [CMAS] Liquid model
C codes for calculating G(mix) as of 26-Nov-2001, copyright Denton S. Ebel.

**  NOTE: Indexing is 1-based locally !!!
	  /* Binary parameters  (These have all been triple-checked!) */
		WG[1][1][1][2]=    318144.87 - t * 154.79;       /* "_WB[1][1]_"  1 for  C  C  C  M  */
		WG[1][1][1][3]=   -617537.61 - t * -76.83;       /* "_WB[1][2]_"  2 for  A  C  C  C  */
		WG[1][1][1][4]=   -960867.88 - t * -247.11;      /* "_WB[1][3]_"  3 for  S  C  C  C  */
		WG[2][2][2][3]=   -691193.17 - t * -227.94;      /* "_WB[1][5]_"  5 for  A  M  M  M  */
		WG[2][2][2][4]=   -610906.87 - t * -195.01;      /* "_WB[1][6]_"  6 for  S  M  M  M  */
		WG[3][3][3][4]=    258911.07 - t * 110.47;       /* "_WB[1][8]_"  8 for  S  A  A  A  */
		WG[1][2][2][2]=    114076.93 - t * 58.66;        /* "_WB[2][1]_" 11 for  C  M  M  M  */
		WG[1][3][3][3]=   -197743.47 - t * -1.11;        /* "_WB[2][2]_" 12 for  A  A  A  C  */
		WG[1][4][4][4]=    -25525.64 - t * 34.19;        /* "_WB[2][3]_" 13 for  S  S  S  C  */
		WG[2][3][3][3]=   -641890.87 - t * -224.47;      /* "_WB[2][5]_" 15 for  A  A  A  M  */
		WG[2][4][4][4]=     94145.91 - t * 52.77;        /* "_WB[2][6]_" 16 for  S  S  S  M  */
		WG[3][4][4][4]=   -161039.81 - t * -60.41;       /* "_WB[2][8]_" 18 for  S  S  S  A  */
		WG[1][1][2][2]=    590616.72 - t * 322.37;       /* "_WB[3][1]_" 21 for  C  C  M  M  */
		WG[1][1][3][3]=    -734020.1 - t * -251.94;      /* "_WB[3][2]_" 22 for  A  A  C  C  */
		WG[1][1][4][4]=   -341962.81 - t * 56.62;        /* "_WB[3][3]_" 23 for  S  S  C  C  */
		WG[2][2][3][3]=     727706.3 - t * 290.88;       /* "_WB[3][5]_" 25 for  A  A  M  M  */
		WG[2][2][4][4]=    -270581.7 - t * -59.98;       /* "_WB[3][6]_" 26 for  S  S  M  M  */
		WG[3][3][4][4]=   1803871.61 - t * 844.79;       /* "_WB[3][8]_" 28 for  S  S  A  A  */
	  /* Ternary parameters */
		WG[1][1][2][3]=  -2440837.52 - t * -526.78;      /* "_WT[1][1]_"  1 for  A  C  C  M  */
		WG[1][1][2][4]=   -2464803.7 - t * -669.0;       /* "_WT[1][2]_"  2 for  S  C  C  M  */
		WG[1][1][3][4]=     580678.7 - t * 526.17;       /* "_WT[1][4]_"  4 for  S  A  C  C  */
		WG[2][2][3][4]=    652384.49 - t * 397.38;       /* "_WT[1][7]_"  7 for  S  A  M  M  */
		 /* Note: WT[1][7] SAMM is in error in Berman's PhD and in diCapitani & Brown '87 */
		 /*       In Beckett's code (and Yoneda's), the correct numbers are printed */
		WG[1][2][2][3]=  -3334297.25 - t * -1148.1;      /* "_WT[2][1]_" 11 for  A  C  M  M  */
		WG[1][2][2][4]=   -2026666.9 - t * -555.03;      /* "_WT[2][2]_" 12 for  S  C  M  M  */
		WG[1][3][3][4]=  -2833471.13 - t * -976.8;       /* "_WT[2][4]_" 14 for  S  A  A  C  */
		WG[2][3][3][4]=  -3201173.35 - t * -1382.29;     /* "_WT[2][7]_" 17 for  S  A  A  M  */
		WG[1][2][3][3]=    343546.79 - t * 160.77;       /* "_WT[3][1]_" 21 for  A  A  C  M  */
		WG[1][2][4][4]=  -1143506.91 - t * -279.9;       /* "_WT[3][2]_" 22 for  S  S  C  M  */
		WG[1][3][4][4]=  -2685775.05 - t * -917.87;      /* "_WT[3][4]_" 24 for  S  S  A  C  */    
		WG[2][3][4][4]=  -1828080.99 - t * -693.44;      /* "_WT[3][7]_" 27 for  S  S  A  M  */
	  /* quaternary parameters */
		WG[1][2][3][4]=   2179011.74 - t * 1328.5;      /* "_WQ[1]_";    1 for  S  A  C  M  */


  /* get activities of endmember liquid components a[] */
    long double rtlngamma[NA+1], q, pp=-3.0; 	/* pp=(1-p) where p is polynomial degree=4 */
	/* Calculate RTln(gamma) as per Berman & Brown 1984 eqn. 22        */
	/*  nc=NA=4; (nc-1)=NR=3 All local variables are indexed to base 1 */
	/* note that if l==i, the W=0, so conditional under last sum not necessary */
	for (m=1; m<=NA; m++) {          /* CHECKED OK vs. independent calculations 2-July-97 */
		rtlngamma[m] = 0.0;	
		if (x[m]>0.0) {
		  for (i=1; i<=(NA-1); i++) {			
			for (j=i; j<=NA; j++) {			
				for (k=j; k<=NA; k++) {				
					for (l=k; l<=NA; l++) {
						  q = ((m==i)?1.0:0.0) + ((m==j)?1.0:0.0) 
								+ ((m==k)?1.0:0.0) + ((m==l)?1.0:0.0);
						   rtlngamma[m] += WG[i][j][k][l] * ( (q*x[i]*x[j]*x[k]*x[l])/x[m]
										+ pp*x[i]*x[j]*x[k]*x[l] );
					}
				}
			}
		  }
		}
	}
	for (i=1; i<=NA; i++) { 	/* note a[0:3] has 0-based index */
		a[i-1] = (rtlngamma[i] != 0.0) ? ( x[i]*exp(rtlngamma[i]/(rJmolk*t)) ) : 0.0;
	}
 /**************************************************************************/

Here is a macro for g(mixing), but it is NOT what I actually use in the
condensation/evap. work:
/*   Start of macros generated 17-April-97 in MAPLE V.
**   w:=array([1..4][1..4][1..4][1..4]);
**   x:=array([1..4]);
*> gexp:=sum(sum(sum(sum(w[i,j,k,l]*x[i]*x[j]*x[k]*x[l],l=k..4),k=j..4),j=i..4),i=1..3)
*           -sum(w[i,i,i,i]*x[i]^4,i=1..4)
*/
#define M2_GMIX   WG[3][3][3][4]*CUBE(x[3])*x[4]+WG[1][1][2][2]*SQUARE(x[1])*SQUARE(x[2]) \
+WG[1][1][2][3]*SQUARE(x[1])*x[2]*x[3]+WG[1][1][2][4]*SQUARE(x[1])*               \
x[2]*x[4]+WG[1][1][3][3]*SQUARE(x[1])*SQUARE(x[3])+WG[1][1][3][4]*SQUARE(x[1])*   \
x[3]*x[4]+WG[1][1][4][4]*SQUARE(x[1])*SQUARE(x[4])+WG[1][2][2][3]*x[1]*           \
SQUARE(x[2])*x[3]+WG[1][1][1][4]*CUBE(x[1])*x[4]+WG[1][1][1][2]*CUBE(x[1])*       \
x[2]+WG[1][4][4][4]*x[1]*CUBE(x[4])+WG[1][3][4][4]*x[1]*x[3]*SQUARE(x[4])+        \
WG[1][3][3][4]*x[1]*SQUARE(x[3])*x[4]+WG[1][3][3][3]*x[1]*CUBE(x[3])+             \
WG[1][2][3][3]*x[1]*x[2]*SQUARE(x[3])+WG[1][2][3][4]*x[1]*x[2]*x[3]*x[4]+         \
WG[1][2][4][4]*x[1]*x[2]*SQUARE(x[4])+WG[1][2][2][4]*x[1]*SQUARE(x[2])*           \
x[4]+WG[1][2][2][2]*x[1]*CUBE(x[2])+WG[2][2][2][3]*CUBE(x[2])*x[3]+               \
WG[2][2][2][4]*CUBE(x[2])*x[4]+WG[1][1][1][3]*CUBE(x[1])*x[3]+WG[2][2][4][4]*     \
SQUARE(x[2])*SQUARE(x[4])+WG[2][2][3][3]*SQUARE(x[2])*SQUARE(x[3])+               \
WG[2][2][3][4]*SQUARE(x[2])*x[3]*x[4]+WG[2][4][4][4]*x[2]*CUBE(x[4])+WG[2][3][3][3]*  \
x[2]*CUBE(x[3])+WG[2][3][3][4]*x[2]*SQUARE(x[3])*x[4]+WG[2][3][4][4]*x[2]*x[3]*   \
SQUARE(x[4])+WG[3][3][4][4]*SQUARE(x[3])*SQUARE(x[4])+WG[3][4][4][4]*             \
x[3]*CUBE(x[4])