Dear all,

I am trying to implement the basic multiphase D2Q9 free energy method (Swift, Orlandini, Yeomans). I am finding this one a lot tougher to implement than Shan-Chen and the color model, in particular I have run into a snag with the chemical potential, usually written \Delta \mu.

According to several papers, including the original Orlandini, Swift and Yeomans, this involves a term of the type

```
log((rho+phi)/(rho-phi))
```

where “rho” is the total density, and “phi” is the difference in density between fluids 1 and 2.

Now, I tried implementing starting with two pure phases, and found that there is a problem for pure phases, that is rho=phi or rho=-phi. In these cases the logarithm turns to plus or minus infinity, and the program well inevitably crash and burn producing "not a number"s galore over my matrices.

Does anyone know what Swift, Orlandini, and Yeomans had in mind when dealing with the pure phases?