Chemical potentials return `NaN` for vanishing components

[david-zwicker]

Compositions where one component has a fraction of identical zero are problematic in some methods. Apparently, the minimization works well (which I actually find surprising). However, calculating chemical potentials and free energies returns NaNs. For example

import flory
f = flory.FloryHuggins(2, 0)
f.free_energy_density([0, 1])

I think we can alleviate some of the problem by defining x * log(x) == 0 for x == 0, which needs some care. The associated chemical potentials would still diverge (but we should then return -inf and inf). These cases of infinity should probably treated separately in the equilibration_error method.

[qiangyicheng]

We can do this. One thing I'm not very sure about is whether these special fp numbers would propagate as we expected, since in these cases we might encounter arithmetics between inf and 0, which might potentially leads to nan. The main issue is the calculation of chemical potential:

f = self.free_energy_density(phis)
j = self.jacobian(phis)
ans = (
    np.atleast_1d(f)[..., None]
    - np.einsum("...i,...i->...", phis, j)[..., None]
    + j
)

Although I can fix jacobian and free_energy_density such that they returns 0 and inf correctly, the np.einsum("...i,...i->...", phis, j) becomes a bit tricky. The problem is that this function in FreeEnergyBase is ignorant of how jacobian and free_energy_density are calculated, so it cannot decide what 0 * inf should be. Maybe we need to provide a specialization for FloryHuggins.

The minimization indeed allows negative/zero volume fractions, since the minimization doesn't truly touch the chemical potential.

[david-zwicker]

I think we can choose either 0 * inf = 0 or 0 * inf = inf, because the result of the expression above will be inf (due to the last term j) in both cases. We could for instance simply return j for cases where np.isinf(j) is True. We just need to document this behavior clearly in the method.