DFT association: analytic branches use sigma of component 0 instead of the associating component #390

In crates/feos/src/association/dft.rs (v0.10.1, still on main), helmholtz_energy_density_ab_analytic (l. 297) and helmholtz_energy_density_cc_analytic (l. 333) call

self.model.association_strength_ij(temperature, 0, 0, par)

The component indices are already determined a few lines above (a.sites_a[0].assoc_comp, a.sites_b[0].assoc_comp, a.sites_c[0].assoc_comp) and used for the diameters, but not here. For PC-SAFT, association_strength_ij scales κ_AB with (σ_i σ_j)^1.5, so whenever the only associating component is not component 0, the functional uses σ_0³ instead of σ_i³. The EOS (association/mod.rs) passes comp_i, comp_j correctly, so functional and EOS disagree.

When it triggers: exactly one associating component with one A and one B site type (or one C site), and that component is not first, e.g. [nitrogen, water] with non-associating N2. If N2 also has a site (e.g. esper2023 N2 with nb), the general cross-association path is taken and everything is consistent.

Reproducer:

import feos
import si_units as si

water = feos.PureRecord(feos.Identifier(name="water"), 18.015, m=2.36948, sigma=2.15072,
                        epsilon_k=230.71557, association_sites=[
                            {"na": 1.0, "nb": 1.0, "kappa_ab": 0.35319, "epsilon_k_ab": 2195.10176}])
nitrogen = feos.PureRecord(feos.Identifier(name="nitrogen"), 28.01, m=1.1504, sigma=3.3848,
                           epsilon_k=91.4)

T, p = 323.15 * si.KELVIN, 40 * si.BAR
for records, feed in [([nitrogen, water], [0.1, 0.9]), ([water, nitrogen], [0.9, 0.1])]:
    params = feos.Parameters.new_binary(records, k_ij=0.0)
    iw = 0 if records[0] is water else 1
    for name, f in [("EOS", feos.EquationOfState.pcsaft(params)),
                    ("DFT", feos.HelmholtzEnergyFunctional.pcsaft(params))]:
        vle = feos.PhaseEquilibrium.tp_flash(f, T, p, feed)
        rho = vle.liquid.density / (si.KILO * si.MOL / si.METER**3)
        print(f"{records[0].identifier.name:>8} first, {name}: rho_liq = {rho:.4f} kmol/m3, "
              f"y_water = {vle.vapor.molefracs[iw]:.5f}")

Output (feos 0.10.1):

nitrogen first, EOS: rho_liq = 56.0585 kmol/m3, y_water = 0.00367
nitrogen first, DFT: rho_liq = 56.2089 kmol/m3, y_water = 0.00099
   water first, EOS: rho_liq = 56.0585 kmol/m3, y_water = 0.00367
   water first, DFT: rho_liq = 56.0585 kmol/m3, y_water = 0.00367

Here the association volume of water is inflated by (σ_N2/σ_water)³ ≈ 3.9. In my water–N2–quartz contact-angle calculations this shifted θ by 4–13°.

Suggested fix: association_strength_ij(temperature, i, j, par) in the AB branch and (temperature, i, i, par) in the CC branch.

Workaround: put the associating component first.