The function RDF[i,j,q] gives the radial distribution function g_{i,j}(r), as prescribed by the RFA method [S.B. Yuste, A. Santos and M. López de Haro, J. Chem Phys, 108, 3683 (1998)], for a hard-sphere mixture defined by:
->nco=number of species
->eeta=packing fraction
->con=list with the molar fractions of each species,e.g, {x_1,x_2,x_3}, 
        x_1 being the corresponding to the species with the largest diameter
->ratio=list with the diameter ratios, e.g., {1,ratio_2,ratio_3} where 
ratio_n=sigma_n/sigma_1=sigma_n
->alfa=the α (alpha) parameter of the RFA method (it can be obtained by means of 
the function ALFA[nco,eeta,con,ratio])

How to use this notebook:
(1) Introduce the  nco, eeta, con , ratio and  alfa values of your hard-sphere mixture in the corresponding cell (see below) , and 
(2) Evaluate the next outest cell. 
After this, you can use the function RDF[i,j,r] to evaluate the RFA radial distribution function at distance r for species i and j of your mixture.

This cell should be evaluated for each mixture

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

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

For example, here they are the values of our Example 3  (ternary mixture)

\(nco = 3; \ \ 
eeta = 49/100; 
con = {1/102, 1/102, 100/102}; 

    ratio = {1, 3/10, 1/10}; \)

ALFA[nco,eeta,con,ratio] gives α= 0.018370493983291558

alfa = 0.018370493983291558;

For example, here they are the values of our Example 1  (binary mixture)

nco=2;   
eeta=49/100;
con={1/2,1/2};
ratio={1,3/10};

ALFA[nco,eeta,con,ratio] gives α=0.027838955619564247

alfa=0.027838955619564247;

numterm is a parameter used in the inversion of the Laplace transform G_ij(s).  A good value is numterm=100, but if you prefer faster calculations  and are not very concerned with the precision of the results, try a smaller number . For example, numterm=50 leads also to very good results and  halves the  computation times.

numterm=100;

Now  we evaluate the Laplace transform, G_ij(s), of   r*g_ij(r)  as prescribed by the RFA method
You must evaluate this cell again if you changes the mixture and/or the value of α.

Clear[rho,reglas,G]

PI = N[Pi]; 
     Do[c[n] = con[[n]], {n, 1, nco}]; 
    Do[sig[n, n] = ratio[[n]], {n, 1, nco}]; 
    rhototal = (6*eeta)/(PI*Sum[c[i]*sig[i, i]^3, {i, 1, nco}]); 
    rho[n_] = rhototal*c[n]; 
    sig[i_, j_] := 1/2*(sig[i, i] + sig[j, j]) /; i != j; 
    zeta[n_, ncomp_] := Sum[rho[i]*sig[i, i]^n, {i, 1, ncomp}]; 
    rhot[ncomp_] := zeta[0, ncomp]; 
    eta[ncomp_] := 1/6*PI*zeta[3, ncomp]; 
    lam[ncomp_] := (2*PI)/(1 - eta[ncomp]); 
    lamp[ncomp_] := (PI^2*zeta[2, ncomp])/(1 - eta[ncomp])^2; 
    g[i_, j_, nc_] := 
     (lam[nc] + (lamp[nc]*sig[i, i]*sig[j, j])/(2*sig[i, j]) + 
        (lamp[nc]^2*((sig[i, i]*sig[j, j])/sig[i, j])^2)/
         (18*lam[nc]))/(2*PI); 
    L[i_, j_, 2, nc_] := 2*PI*alfa*sig[i, j]*g[i, j, nc]; 
    L[i_, j_, 0, nc_] := 
      lam[nc] + lamp[nc]*sig[j, j] + 2*lamp[nc]*alfa - 
       lam[nc]*Sum[rho[k]*sig[k, k]*L[k, j, 2, nc], {k, 1, nc}]; 
    L[i_, j_, 1, nc_] := 
      lam[nc]*sig[i, j] + 1/2*lamp[nc]*sig[i, i]*sig[j, j] + 
       (lam[nc] + lamp[nc]*sig[i, i])*alfa - 
       1/2*lam[nc]*sig[i, i]*
        Sum[rho[k]*sig[k, k]*L[k, j, 2, nc], {k, 1, nc}]; 
        
        A[i_,j_,nc_]:=rho[i]*
((1-s sig[i,i]+1/2 s^2 sig[i,i]^2-xx[i])/s^3*
L[i,j,0,nc]+
(1-s sig[i,i]-xx[i])/s^2*L[i,j,1,nc]+
(1-xx[i])/s*L[i,j,2,nc]);

L[i_,j_,nc_]:=L[i,j,0,nc]+L[i,j,1,nc]*s+L[i,j,2,nc]*s^2;
ML[nc_]:=Table[L[i,j,nc],{i,nc},{j,nc}];
MA[nc_]:=Table[A[i,j,nc],{i,nc},{j,nc}];
in[nc_]:=Inverse[(1+alfa s)IdentityMatrix[nc]-MA[nc]];
MJota[nc_,eet_,co_,rati_]:=
1/(2Pi) ML[nc].in[nc];

jo=MJota[nco,eeta,con,ratio]//Together//Chop//Simplify;

reglas=Table[xx[n]->Exp[-sig[n,n]s],{n,1,nco}];
G[i_,j_,eet_,co_,rati_]:=
G[i,j,eet,co,rati]=jo[[i,j]]/s^2//.reglas;
fs[ss_,i_,j_]:=Re[G[i,j,eeta,con,ratio]/.s->ss//N]//Chop;

G_ij(s) is  inverted numerically.

Geuler[ii_,jj_,t_,ntr_,meuler_]:=Module[
{preci=16,
aa=30,hh=Pi/t,
suma,avgsu=0,avgsu1=0,fun,fun1,c,xx,su},
c[m_]:=c[m]=Binomial[meuler,m-1];

uu=Exp[aa/2]/t;
xx=aa/(2*t);
suma=fs[xx,ii,jj]/2+
     Sum[(-1)^nn*fs[xx+I*nn*hh,ii,jj],
         {nn,1,ntr}];
su[1]=suma;
Do[nn=ntr+k;
   su[k+1]=
   N[su[k]+
    (-1)^nn*fs[xx+I*nn*hh,ii,jj],preci],
{k,1,meuler+1}];
Do[avgsu=N[avgsu+c[j]*su[j],preci];
   avgsu1=N[avgsu1+c[j]*su[j+1],preci],
{j,1,meuler+1}];
fun=N[uu*avgsu/2^meuler,preci];
fun1=N[uu*avgsu1/2^meuler,preci];
fun1];
RDF[i_,j_,r_]:=Geuler[i,j,r-sig[i,j],numterm,30]/r/;r!=sig[i,j];
RDF[i_,j_,r_]:=g[i,j,nco]/;r==sig[i,j];

EXAMPLES

Example 1. A binary mixture with
  nco=2;   
eeta=49/100;
con={1/2,1/2};
ratio={1,3/10};
alfa=0.027838955619564247;

Here we have used numterm=100

Here we have used numterm=50

Example 2. A binary mixture with
nco=2;   
eeta=49/100;
con={1/16,15/16};
ratio={1,3/10};
alfa=0.018859574773073649;

Here we have used numterm=50

Example 3. A ternary mixture with
nco=3;  
eeta=49/100;
con={1/102,1/102,100/102};
ratio={1,3/10,1/10};
alfa=0.018370493983291558;

Here we have used numterm=100

Here we have used numterm=50