The calculation of the mixing energy is extremely complicated and only approximate solutions have been obtained as far as I am aware. I shall attempt to outline one of these methods*.
All molecules are assumed to be spherical and of equal size and situated on a lattice with co-ordination number z . Each pair of neighbours AA of type A molecules contributes energy $w_{AA}$, similarly $w_{BB}$ for B molecules and energy $w_{AB}$ for AB pairs. The total energy due to these interactions
is
$$E_T=N_{AA}w_{AA} + N_{BB}w_{BB}+N_{AB}w_{AB}$$
where $N_{AA}$ is the number of AA pairs and so on. Counting the number of $A$ interactions and number of B ’s gives
$$zN_A=2N_{AA} + N_{AB} ~~~\text{and }~~~ zN_B=2N_{BB} + N_{AB} $$
and thus combining these equations gives
$$E_T=\frac{zN_Aw_{AA}}{2}+\frac{zN_Bw_{BB}}{2}+wN_{AB}=E^T_{AA}+E^T_{BB}+wN_{AB}$$
where $w=w_{AB}-w_{AA}/2-w_{BB}/2$ is the energy gain on mixing per pair of AB neighbours.
In statistical mechanics the partition function is the key to all other thermodynamic functions. In this case partition function is $$
Z= \sum_{N_{AB}}Z_AZ_Bg(N_A,N_B,N_{AB})e^{(-E^T_{AA}-E^T_{BB}-wN_{AB}/KT)}$$
and where $g(..) $ is the number of arrangements of $N_A$ A’s , $N_B$ B ’s and $N_{AB}$ pairs of molecules. Separating out the terms due to pure A and B molecules gives
$$Z=Z’_AZ’_B \sum_{N_{AB}}g(N_A,N_B,N_{AB})e^{-wN_{AB}/KT}$$
and as the Helmholtz energy is
$F=-kT\ln(Z)=F_A+F_B+\Delta F$
and for mixing only the $\Delta F$ term is important the others being the energy of each pure component.
Thus
$$\Delta F = -kT\ln\left(\sum_{N_{AB}}g(N_A,N_B,N_{AB})e^{-wN_{AB}/kT} \right) $$
The change in energy can be obtained from the Gibbs-Helmholtz equation $(E = F-T(dF/dT))_V$ and is
$$\Delta E= w\frac{\sum_{N_{AB}} N_{AB} g(N_A,N_B,N_{AB})e^{-wN_{AB}/kT} } { \sum_{N_{AB}}g(N_A,N_B,N_{AB})e^{-wN_{AB}/kT} } $$
and this is an equation that seems to have no exact solution and so this is where approximations have to be made.
The difficulty is in the exponential terms, as they contain $N_{AB}$ the aim is therefore to replace these with an average value, and thus the exponential terms become constant wrt $N_{AB}$. It can be shown that
$$\sum_{N_{AB} }g(N_A,N_B,N_{AB}) =\frac{(N_A + N_B)!}{N_A!N_B!}$$
and
$$\bar N_{AB}=w\frac{\sum_{N_{AB}}N_{AB}g(N_A,N_B,N_{AB})}{\sum_{N_{AB}}g(N_A,N_B,N_{AB})} = z\frac{N_AN_B}{N_A + N_B} $$
By replacing $N_{AB}$ by $\bar N_{AB}$ in the partition function produces
$$Z=Z’_AZ’_B \sum_{N_{AB}}g(N_A,N_B,N_{AB})e^{-w\bar N_{AB}/KT)}
= \frac{(N_A+N_B)!}{N_A!N_B!}\exp\left(-\frac{zw}{kT}\frac{N_AN_B}{N_A+N_B}\right)$$
the change in energy is then
$$\Delta E = zw\frac{N_AN_B}{N_A+N_B}= w\bar N_{AB}$$
and this can be re-written as molec fractions to produce
$$ \Delta E= z(N_A+N_B)x_1x_2w$$
- based on text in Rushbrooke ‘Introduction to Statistical Mechanics’
