Critical temperature-thermodynamics-simulation

Hey fellows :)
How do you do these lovely days?
I'm at the end of twelfth grade righy now and I have to finish a work I am working on since about 2 years. It will count as an extra subject on my final diploma then it's really important!
I need your help guys, I'll be gratefull to get some...
(My question is actually pure thermodynamics therefore I post it here. Hope that's the place.)

So, here is the thing:
I wont make you all tired by telling you about my work, I'll just jump in to the problematic part - During the work I wrote down the free energy of binary system (such as water-oil system) as a function of one of the components concentration. By simple differentiation I got a function of the transition temperature (between one mixed phase or two separated phases) for every given concentration. My teacher confirmed it.
Simultaneously, I used Metropolis algorithm to create a simulation of the system for various temperatures, interaction values and concentration (which also passed my teacher's test).
I wanted to compare between the analytical results and the umerical results.
For such purpose I had to provide a way to determine when the system in the simulation is mixed or not... Something better than just "how it's looks".
It was decided (by me of course cuse I'm so cool :P ) to count how many neighbors any particle has from the same tipe as itself. Lets use the sign η for this parameter.
As we get closer to 4, we are getting closer to phases separation.
That is what I say: The mean field theory tells us that at the best mixed state we axpect too see η=4(c1^2+c2^2) -where c1,c2 are the concentrations- and at the opposite state we are axpecting to η=4, which in fact is almost the same as the results from the simulation (up to +/- 0.1).
My question to you is what value of μ should we get around the transition temperature? with naive thingking as mine I thought to see the average of 4 & 4(c1^2+c2^2) but the reality is quite different:
When I insert 0.5 as the concentrations, -300 for same tipe body interaction, 300 for different tipe body interaction and temperature of 1200 Kelvin, I get μ=2.55. One thing I forgot to mention is that the code had been written when Boltzmann constant is affixed to be 1. That's why the transition temperature should be simply J_T (1200 in our case).
Unfortunately, μ=2.55 is a bad result. From the average idea we want η to be 3 in this system, a value that comes out when I set the temperature to be 780 Kelvin - clearly not 1200.
Now, under the assumption that my temperatures's function and my simulation's code are both fine, what's the problem? The problem has to be, I think, with this stupid average idea.
I'm here cuse I need you, enthusiastic Young physicists (advice for life: don't waste your all life with The Big Bang Theory as I did), to come up with a clever solution and an explanation for that mismatch.

Thanks you very very much!!

 PhysOrg.com physics news on PhysOrg.com >> Promising doped zirconia>> New X-ray method shows how frog embryos could help thwart disease>> Bringing life into focus