feos-core: binary critical point at fixed p, tp_flash warm start, newton (implicit AD), residual μ contributions, dew_point_line #401

the numbers use PC-SAFT with esper2023.json + rehner2023_binary.json.

1. critical_point_binary at fixed pressure never converges (state/critical_point.rs:573)

criticality_conditions_p evaluates ∂A/∂V at V = 1 (reduced), with mole fractions as the moles, so the EoS sees ρ = 1 Å⁻³ instead of the actual density:

let (_, da) = first_derivative(a, D::one());              // current
let (_, da) = first_derivative(a, density.recip());       // fix

Check: compute the critical point at T, then compute it again at the resulting p_c. For methane/propane, CO₂/n-butane and ethane/heptane at three temperatures each, all 9 runs fail with "Critical point did not converge", even when started at the exact solution. With the fix, all 9 return the original T, x and ρ.

The same call is used by Python State.critical_point_binary(eos, p) and by PhaseDiagram::binary_vle at fixed p when one component is supercritical. present since v0.9.0

2. tp_flash with initial_state returns the old vapor fraction when only the feed changes (phase_equilibria/tp_flash.rs:131-141)

At the same T and p, the old x and y already satisfy K·x = y, whatever the feed. So the first SS step passes the tolerance, and Rachford–Rice never runs for the new feed.

The result is the previous flash: same β, same phase compositions, with the new feed's total moles. No error is raised.

With the natural loop vle = tp_flash(T, p, z, initial_state=vle), the whole sweep freezes on the first result. Propane/n-butane/pentane at 300 K, 4 bar, z₁ = 0.40 → 0.70:

Warm-starting each flash from a cold one gives a one-step lag instead: 24 of 25 points wrong for a binary, 30 of 30 for a ternary.

No in-repo caller is affected; PhaseDiagram::lle changes T or p, and the ph/ps flashes reuse the same feed. Suggested fix (untested): run Rachford–Rice for the new feed before the first convergence check. For example, call rachford_rice with the old K in the initial_state branch, instead of reusing vapor_phase_fraction(). That fixes β, and the successive-substitution loop then moves x and y to the new tie line.

3. newton returns the state from before the final implicit-AD step (state/mod.rs:620-627)

In the convergence branch, state = s is taken from f(x0) before x0 -= fx / dfx, so the last dual update is never used. If the initial guess is already converged in its real part, the derivative is lost completely. This happens when you warm start from the previous solution.

Example: State::new_nph with PcSaftBinary seeded on k_ij (the parameters from px_flashes.rs, 5 bar, x = 0.3, 600 K vapor):

This affects new_nph, new_nps, new_nth, new_nts and new_nvu with dual D; f64 and Python are not affected.

Fix: split the new feed with the old K values before iterating. If that fails, fall back to the cold start as before.

if let Some(initial_state) = initial_state {
    let mut init = initial_state.clone();
    // the feed may have changed: split it with the old K values
    let k = initial_state
        .vapor()
        .molefracs
        .component_div(&initial_state.liquid().molefracs);
    if init.rachford_rice_inplace(self, &k).is_ok() {
        let vle = self.tp_flash_(init, options, non_volatile_components.clone());
        if vle.is_ok() {
            return vle;
        }
    }
}

With this change, the chained sweeps above match the cold flashes at every point: the ternary VLE, and LLE in methanol/hexane at 280 K, which without the fix is stuck at β = 0.6172.

4. residual_chemical_potential_contributions passes T as the molar volume (state/residual_properties.rs:562)

let v = Dual::from_re(self.temperature.into_reduced());          // current
let v = Dual::from_re(self.density.into_reduced().recip());      // fix (as in chemical_potential_contributions)

This only affects rust users.

5. dew_point_line repeats the point where the T sweep hands over to the p sweep (phase_equilibria/phase_envelope.rs:110-112)

The pressure sweep starts at min_pressure, which is the last point of the temperature sweep, so that dew point is computed twice. Changing n_p to n_p - 1 would remove the duplicate, but it would also make the result one point shorter than bubble_point_line. Starting the sweep one step above min_pressure keeps the count.