With this notebook you can  calculate the following quantities for a BINARY (or one-component) mixture of STICKY Hard Spheres:

--> yHsBmcsl[i,j,nco,eeta] =  the contact value of the cavity function for HS given by the BMCSL formula.
-->yShsPy[i,j,nco,con,eeta,taus] = the contact value of the cavity function for SHS given by the PY approximation.
--> yHsPy[i,j,nco,eeta] = the contact value of the cavity function for HS given by the PY approximation.
--> yProCon[i,j,nco,con,eeta,taus] = the contact value of the cavity function for SHS given by equation (3.1) of the paper  A. Santos, S. B. Yuste and M. López de Haro, (J. Chem. Phys. vol 109, num. 16, 22 October 1998.
--> 1/ chim1HsBmcsl[nco,eeta] = the isothermal susceptibility given by the BMCSL formula for the HS mixture.
--> chiShsPy[nco,eeta,taus] = the isothermal susceptibility given by the PY approximation for the SHS mixture.
--> 1/ chim1HsPy[nco,eeta] = the isothermal susceptibility given by the PY approximation for the HS mixture.
 --> chiPro[nco,eeta,taus] = the isothermal susceptibility of the SHS mixture  given by the approximate formula (3.2) of the paper  A. Santos, S. B. Yuste and M. López de Haro, (J. Chem. Phys. vol 109, num. 16, 22 October 1998).

The mixture is defined by:
-> nco=number of species
-> eeta=volume packing fraction
-> con=list with the molar fractions of each species,e.g, {x_1,x_2}, 
 where x_1 is the molar fraction of the species with the largest diameter
-> ratio=list with the diameter ratios, e.g., {1,ratio_2} 
where ratio_2=sigma_2/sigma_1=sigma_2
-> taus=list with the tau parameters characterizing the adhesiveness of the intermolecular potential , e.g., {{tau_11,tau_12},{tau_21,tau_22}}.

How to use this notebook:
(1) Introduce the  nco, con  and ratio  (the values of eta and taus can be given later)  of your hard-sphere mixture in the corresponding cell (see below) , and 
(2) Evaluate the next outest cell.

This cell must be evaluated for each mixture  with different values of nco, con and/or ratio

Here you have to introduce the values of nco, eeta, con , ratio and  taus corresponding to your sticky-hard-sphere mixture

For example, these are the values for our Examples  1 and 2  (binary mixtures)

\(\(nco = 2; 
\ con = {1/2, 1/2}; 
ratio = {1, 1}; \ \ \)\)

For the  Example 3  (binary mixture) we used this values:

nco=2;
 con={2/3,1/3};
ratio={1,1/2};

For the  Example 4  (one-component mixture) we used this values:

nco=2;   
con={1/2,1/2};
ratio={1,1};

The core of the program...

Off[General::spell];Off[General::spell1]

Clear[solu,ecu,sol,sol1,vartest,grial,kLpy,AA,invdelA,mH0,mH1,inv1hgorro]

\(\(PI\  = \ N[Pi]; \ 
\ \ \ \ \ 
    Do[c[n]\  = \ con[\([n]\)], \ {n, \ 1, \ nco}]; \ 
\ \ \ \ \ \ \ 
    Do[sig[n, n] = \ ratio[\([n]\)], \ {n, \ 1, \ nco}]; \ 

    sig[i_, \ j_]\  := \ 
      1/2*\((ratio[\([i]\)]\  + ratio[\([j]\)])\)\  /; \ i\  != \ j; \ \ 
    \)\)

\(\(\ \ \ \ \ \ 
    rhototal[nco_, eeta_] := \ 
      \((6*eeta)\)/\((PI*Sum[c[i]*sig[i, i]^3, \ {i, \ 1, \ nco}])\); \ 

    \ \ \ \ rho[n_, nco_, eeta_] := \ rhototal[nco, eeta]*c[n]; \ \ \ \ \ 

    \ \ \ \ zeta[n_, \ nco_, eeta_]\  := \ 
      Sum[rho[i, nco, eeta]*sig[i, \ i]^n, \ {i, \ 1, \ nco}]; \ 
\ \ \ \ 
    rhot[\ nco_, eeta_]\  := \ zeta[0, nco, eeta]; \ 
\ \ \ \ 
    eta[nco_, eeta_]\  := \ 1/6*PI*zeta[3, nco, eeta]; \ 
\ \ \ \ 
    lam[nco_, eeta_]\  := \ \((2*PI)\)/\((1\  - \ eta[nco, eeta])\); \ 

    \ \ \ \ lamp[nco_, eeta_]\  := \ 
      \((PI^2*zeta[2, \ nco, eeta])\)/\((1\  - \ eta[nco, eeta])\)^2; \)\)

\(solu[nco_, con_, eeta_, taus_]\  := 
      \(solu[nco, con, eeta, taus]\  = 
        \((Clear[ecu, sol, sol1, grial, kk3]; 
		
          ecu := Table[
              12*taus[\([i, j]\)] kLpy[i, j]/sig[i, j] == 
                lam[nco, eeta]*sig[i, j] + 
                  1/2*lamp[nco, eeta]*sig[i, i]*sig[j, j] - 
                  1/2*lam[nco, eeta]*
                    Sum[rho[k, nco, eeta]*sig[k, k]*
                        \((kLpy[k, i]*sig[j, j] + kLpy[k, j]*sig[i, i])\), {
                        k, 1, nco}] + \ 
                  Sum[rho[k, nco, eeta]*kLpy[k, i]*kLpy[k, j], {k, 1, nco}], {
                i, nco}, {j, nco}]; 
		
          sol = NSolve[Flatten[ecu], 
              Flatten[Table[kLpy[i, j], {i, nco}, {j, nco}]]\ ]; 
			
          sol1 = Table[
              Table[sol[\([n, m, 2]\)], {m, 1, nco^2}], {n, 1, nco^2}]\ ; 
          						
			vartest = False; 
          grial = "Ninguno real positivo"; 

          Do[\ 
	\ \ \ 
            Do[
				\ \ \ \ \ \ \ 
              If[Head[sol1[\([i, j]\)]] === Real, 
							
                If[sol1[\([i, j]\)] > 0, 
				\ \ vartest = True, 

                  				\ \ vartest = False; Break[]\ ], vartest = False; 
                Break[]\ ], 
		{j, 1, nco^2}]; 
            If[vartest, grial = Table[sol1[\([i, jj]\)], {jj, 1, nco^2}]; 
              Break[], ], 
{i, 1, nco^2}]; 
		
          Table[grial[\([i]\)], {i, nco^2}]\ )\)\)\)

Lpy[i_, j_, 2, nco_, eeta_, taus_] := 
      solu[nco, con, eeta, taus][[i + j - 1]] /; i <= j

Lpy[i_, j_, 2, nco_, eeta_, taus_] := 
      Lpy[j, i, 2, nco, eeta, taus] /; i > j

yShsPy[i_, j_, nco_, con_, eeta_, taus_]  := 
      solu[nco, con, eeta, taus][[i + j - 1]] *
          (12*taus[[i, j]])/(2*PI*sig[i, j]^2) /; i <= j

yShsPy[i_, j_, nco_, con_, eeta_, taus_]  := 
      yShsPy[j, i, nco, con, eeta, taus]  /; i > j

\(\(\ \(Lpy[i_, \ j_, \ 0, \ nc_, eeta_, taus_]\  := \ 
\ \ \ \ \ \ 
      lam[nc, eeta]\  + \ lamp[nc, eeta]*sig[j, \ j]\  - \ 
\ \ \ \ \ \ \ 
        lam[nc, eeta]*
          Sum[rho[k, nc, eeta]*sig[k, \ k]*
              Lpy[k, \ j, \ 2, \ nc, eeta, taus], \ {k, \ 1, \ nco}]; \)\ 

    \)\)

\(\(\(Lpy[i_, \ j_, \ 1, \ nc_, eeta_, taus_]\  := \ 
\ \ \ 
      lam[nc, eeta]*sig[i, \ j]\  + \ 
        \((1*lamp[nc, eeta]*sig[i, \ i]*sig[j, \ j])\)/2\  - \ 
\ \ \ \ 
        \((1*lam[nc, eeta]*sig[i, \ i]*
              Sum[rho[k, nc, eeta]*sig[k, \ k]*
                  Lpy[k, \ j, \ 2, \ nc, eeta, taus], \ {k, \ 1, \ nco}])\)/
          2; \)\ \)\)

alfa = 0;   L[i_, j_, 3, nco_, eeta_, taus_] = 0;

L[i_,  j_,  2,  nc_, eeta_, taus_]  :=  
      Lpy[i,  j,  2,  nc, eeta, taus];

L[i_,  j_,  0,  nc_, eeta_, taus_]  :=  
      Lpy[i,  j,  0,  nc, eeta, taus];

L[i_, j_, 1, nc_, eeta_, taus_] := Lpy[i, j, 1, nc, eeta, taus];

\(\(A[i_, j_, nc_, eeta_, taus_] := 
      rho[i, nc, eeta]*

        \((\((1 - s\ sig[i, i] + 1/2\ s^2\ sig[i, i]^2 - xx[i])\)/s^3*

              L[i, j, 0, nc, eeta, taus] + 

            \((1 - s\ sig[i, i] - xx[i])\)/s^2*L[i, j, 1, nc, eeta, taus] + 

            \((1 - xx[i])\)/s*L[i, j, 2, nc, eeta, taus] - 
            xx[i]*L[i, j, 3, nc, eeta, taus])\); \)\)

\(\(
\(L[i_, j_, nc_, eeta_, taus_] := 
      L[i, j, 0, nc, eeta, taus] + L[i, j, 1, nc, eeta, taus]*s + 
        L[i, j, 2, nc, eeta, taus]*s^2 + L[i, j, 3, nc, eeta, taus]*s^3; 
    \)\)\)

AA[i_, j_, n_, nc_,eeta_,taus_] := AA[i,j,n,nc,eeta,taus]=
          Simplify[ (-1)^n*rho[i,nc,eeta]*
          (-(L[i ,j ,3 ,nc,eeta,taus]*sig[i,i]^n)/n!
           +(L[i, j, 2, nc,eeta,taus]*sig[i, i]^(n + 1))/(n + 1)! 
           -(L[i, j, 1, nc,eeta,taus]*sig[i, i]^(n + 2))/(n + 2)! 
           +(L[i, j, 0, nc,eeta,taus]*sig[i, i]^(n + 3))/(n + 3)! )];

delA[nc_,eeta_,taus_]:=IdentityMatrix[nc]-
Table[AA[k,j,0,nc,eeta,taus],{k,nc},{j,nc}]

invdelA[nc_,eeta_,taus_]:=invdelA[nc,eeta,taus]=Inverse[delA[nc,eeta,taus]]

BB0[i_,j_,nc_,eeta_,taus_]:=
1/(2PI)L[i,j,2,nc,eeta,taus]+
Sum[AA[k,j,2,nc,eeta,taus],{k,1,nc}]-
alfa sig[i,j]+
Sum[sig[i,k]AA[k,j,1,nc,eeta,taus],{k,1,nc}]-
1/2 sig[i,j]^2+
Sum[1/2 sig[i,k]^2 AA[k,j,0,nc,eeta,taus],{k,1,nc}]
mBB0[nc_,eeta_,taus_]:=Table[BB0[i,j,nc,eeta,taus],{i,nc},{j,nc}]

mH0[nc_,eeta_,taus_]:=mH0[nc,eeta,taus]=
mBB0[nc,eeta,taus].invdelA[nc,eeta,taus]//Simplify

H[i_,j_,0,nc_,eeta_,taus_]:=mH0[nc,eeta,taus][[i,j]]

BB1[i_,j_,nc_,eeta_,taus_]:=
1/(2PI)L[i,j,3,nc,eeta,taus]+
Sum[AA[k,j,3,nc,eeta,taus],{k,1,nc}]+
Sum[sig[i,k] AA[k,j,2,nc,eeta,taus],{k,1,nc}]-
(1/2 sig[i,j]^2+H[i,j,0,nc,eeta,taus])alfa+
Sum[(1/2sig[i,k]^2+H[i,k,0,nc,eeta,taus])*
    AA[k,j,1,nc,eeta,taus],{k,1,nc}]-
(1/6 sig[i,j]^3+sig[i,j]H[i,j,0,nc,eeta,taus])+
Sum[(1/6 sig[i,k]^3+sig[i,k]H[i,k,0,nc,eeta,taus])*
AA[k,j,0,nc,eeta,taus],{k,1,nc}]

mBB1[nc_,eeta_,taus_]:=Table[BB1[i,j,nc,eeta,taus],{i,nc},{j,nc}]

mH1[nc_,eeta_,taus_]:=mH1[nc,eeta,taus]=
mBB1[nc,eeta,taus].invdelA[nc,eeta,taus]//Simplify

H[i_,j_,1,nc_,eeta_,taus_]:=mH1[nc,eeta,taus][[i,j]]

hgorro[i_,j_,nc_,eeta_,taus_]:=-4*PI*Sqrt[rho[i,nc,eeta]*rho[j,nc,eeta]]*H[i,j,1,nc,eeta,taus]

inv1hgorro[nc_,eeta_,taus_]:=inv1hgorro[nc,eeta,taus]=
Inverse[IdentityMatrix[nc]+
Table[hgorro[i,j,nc,eeta,taus],{i,nc},{j,nc}]]

chiShsPy[nc_,eeta_,taus_]:=1/Sum[Sqrt[c[i]*c[j]]*inv1hgorro[nc,eeta,taus][[i,j]],{i,nc},{j,nc}]

\(yHsBmcsl[i_, j_, nc_, et_] := 
      1/\((2\ PI)\)*

        \((lam[nc, et] + 
            1/2\ lamp[nc, et] \((sig[i, i] sig[j, j])\)/sig[i, j] + 

            1/18\ lamp[nc, et]^2/lam[nc, et] 
              \((sig[i, i] sig[j, j]/sig[i, j])\)^2)\)\)

\(yHsPy[i_, j_, nc_, et_] := 
      1/\((2\ PI)\)*

        \((lam[nc, et] + 
            1/2\ lamp[nc, et] \((sig[i, i] sig[j, j])\)/sig[i, j])\)\)

yProCon[i_, j_, nc_, con_, et_, taus_] := 
      yHsBmcsl[i, j, nc, et] + yShsPy[i, j, nc, con, et, taus] - 
        yHsPy[i, j, nc, et]

\(chim1HsBmcsl[nc_, et_] := 
      \((1/rhototal[nc, et])\)*
        \((
rhototal[nc, et]/\((1 - eta[nc, et])\)^2 + 

            PI\ zeta[1, nc, et] zeta[2, nc, et]/\((1 - eta[nc, et])\)^3 + 

            PI^2/36\ zeta[2, nc, et]^3*
\ \ \ \ \ 
              \((9 - 4  eta[nc, et] + eta[nc, et]^2)\)/
                \((1 - eta[nc, et])\)^4)\)\)

\(chim1HsPy[nc_, et_] := 
      \((1/rhot[nc, et])\)*
        \((
rhot[nc, et]/\((1 - eta[nc, et])\)^2 + 

            PI\ zeta[1, nc, et] zeta[2, nc, et]/\((1 - eta[nc, et])\)^3 + 

            PI^2/4\ zeta[2, nc, et]^3/\((1 - eta[nc, et])\)^4)\)\)

Valor de chi=χ según fórmula (3.2)

chiPro[nco_, eeta_, taus_] := 
      1/(chim1HsBmcsl[nco, eeta] + 1/chiShsPy[nco, eeta, taus] - 
            chim1HsPy[nco, eeta])

EXAMPLES

Example 1. A binary mixture with
 nco=2;   
con={1/2,1/2};
ratio={1,1};
eeta=0.4;
taus={{1,10},{10,1}};

eeta=0.4;
taus={{1,10},{10,1}};

{chim1HsBmcsl[2, eeta], 1/chiShsPy[nco, eeta, taus], 
      chim1HsPy[nco, eeta], 1/chiPro[nco, eeta, taus]}

\({Table[yHsBmcsl[i, j, nco, eeta], {i, nco}, {j, nco}], 
	
        Table[yShsPy[i, j, nco, con, eeta, taus], {i, nco}, {j, nco}], 
	
        Table[yHsPy[i, j, nco, eeta], {i, nco}, {j, nco}], 
	
        Table[yProCon[i, j, nco, con, eeta, taus], {i, nco}, {j, nco}]} // 
      TableForm\)

Example 2 A binary mixture with
 nco=2;   
con={1/2,1/2};
ratio={1,1};
eeta=0.4;
taus={{2,10},{10,2}};

eeta=0.4;
taus={{2,10},{10,2}};

{chim1HsBmcsl[2, eeta], 1/chiShsPy[nco, eeta, taus], 
      chim1HsPy[nco, eeta], 1/chiPro[nco, eeta, taus]}

\({Table[yHsBmcsl[i, j, nco, eeta], {i, nco}, {j, nco}], 
	
        Table[yShsPy[i, j, nco, con, eeta, taus], {i, nco}, {j, nco}], 
	
        Table[yHsPy[i, j, nco, eeta], {i, nco}, {j, nco}], 
	
        Table[yProCon[i, j, nco, con, eeta, taus], {i, nco}, {j, nco}]} // 
      TableForm\)

Example 3. A binary mixture with
 nco=2;
 con={2/3,1/3};
ratio={1,1/2};
eeta=0.4;
taus={{1,10},{10,1}};

eeta=0.4;
taus={{1,10},{10,1}};

{chim1HsBmcsl[2, eeta], 1/chiShsPy[nco, eeta, taus], 
      chim1HsPy[nco, eeta], 1/chiPro[nco, eeta, taus]}

\({Table[yHsBmcsl[i, j, nco, eeta], {i, nco}, {j, nco}], 
	
        Table[yShsPy[i, j, nco, con, eeta, taus], {i, nco}, {j, nco}], 
	
        Table[yHsPy[i, j, nco, eeta], {i, nco}, {j, nco}], 
	
        Table[yProCon[i, j, nco, con, eeta, taus], {i, nco}, {j, nco}]} // 
      TableForm\)

Example 4. A one-component mixture with eeta=0.49 and tau=1. We evaluate this example as a binary mixture with
nco=2;
con={1/2,1/2};
ratio={1,1}; 
eeta=0.49;
taus={{1,1},{1,1}};

eeta=0.49;
taus={{1,1},{1,1}};

{chim1HsBmcsl[nco, eeta], 1/chiShsPy[nco, eeta, taus], 
      chim1HsPy[nco, eeta], 1/chiPro[nco, eeta, taus]}

\({Table[yHsBmcsl[i, j, nco, eeta], {i, nco}, {j, nco}], 
	
        Table[yShsPy[i, j, nco, con, eeta, taus], {i, nco}, {j, nco}], 
	
        Table[yHsPy[i, j, nco, eeta], {i, nco}, {j, nco}], 
	
        Table[yProCon[i, j, nco, con, eeta, taus], {i, nco}, {j, nco}]} // 
      TableForm\)