proteus.outgas.atmodeller¶

def calc_surface_pressures_atmodeller(dirs: dict, config: Config, hf_row: dict):
    """Compute volatile partitioning using atmodeller.

    Solves for chemical equilibrium between the magma ocean and
    atmosphere, accounting for real gas EOS, solubility laws, and
    optionally condensation. Updates hf_row with partial pressures,
    VMRs, dissolved masses, and total surface pressure.

    Parameters
    ----------
    dirs : dict
        Directory paths.
    config : Config
        PROTEUS configuration object.
    hf_row : dict
        Current helpfile row (modified in place).

    Raises
    ------
    ValueError
        If the kg total of a constrained element or a noble gas is not finite, or
        if ``planet.fO2_source = "from_O_budget"``, the O budget is at or below
        ``outgas.mass_thresh``, and another element has a mass constraint (H, C,
        N or S above the threshold, or a noble gas total above 0).
    RuntimeError
        If an element outside the solve has atmosphere or melt mass, or if an
        O2-only solve gives ``P_surf`` other than ``O2_bar``.

    Notes
    -----
    Under ``user_constant`` with no budget above the threshold, the O2-only network
    is solved against the buffer. Under ``from_O_budget`` with no budget above the
    threshold, no solve runs: every species column but the vmr is 0, ``P_surf`` is 0,
    each noble gas keeps its kg total with the matching mole total, the vmrs and
    ``atm_kg_per_mol`` keep their values as in ``run_desiccated``, and
    ``fO2_shift_IW_derived`` and ``O_res`` are NaN.
    """
    from atmodeller import (
        ChemicalSpecies,
        EquilibriumModel,
        Planet,
        SpeciesNetwork,
    )
    from atmodeller.solubility import get_solubility_models
    from atmodeller.thermodata import IronWustiteBuffer

    from proteus.utils.constants import (
        M_earth,
        condensate_list,
        element_mmw,
        gas_list,
        noble_gases,
        vol_element_list,
    )

    atm_config = config.outgas.atmodeller
    solubility_models = get_solubility_models()

    # Real gas EOS models (None = ideal gas)
    from atmodeller.eos import get_eos_models

    eos_models = get_eos_models()
    _eos_map = {
        'H2O': atm_config.eos_H2O,
        'CO2': atm_config.eos_CO2,
        'H2': atm_config.eos_H2,
        'CH4': atm_config.eos_CH4,
        'CO': atm_config.eos_CO,
    }
    _eos_map = {k: v for k, v in _eos_map.items() if v is not None}

    # Build species network
    species_list = []

    # Map PROTEUS gas species to atmodeller species with solubility
    _sol_map = {
        'H2O': atm_config.solubility_H2O,
        'CO2': atm_config.solubility_CO2,
        'H2': atm_config.solubility_H2,
        'N2': atm_config.solubility_N2,
        'S2': atm_config.solubility_S2,
        'CO': atm_config.solubility_CO,
        'CH4': atm_config.solubility_CH4,
    }
    # Remove None entries (solubility disabled for that species)
    _sol_map = {k: v for k, v in _sol_map.items() if v is not None}

    # Noble gases dissolve by the Jambon, Weill and Braun (1986) Henry's law,
    # the same calibration CALLIOPE uses, so the two backends agree on noble
    # gas solubility. The models ship with atmodeller.
    for gas in noble_gases:
        _sol_map[gas] = f'{gas}_basalt_jambon86'

    # Only include species whose constituent elements ALL have non-zero budgets.
    # Atmodeller's solver fails when the species network introduces elements
    # with no mass constraint (under-determined system -> no convergence).
    # E.g., CH4 has both H and C; if only H has a budget, including CH4
    # leaves C unconstrained.
    _species_elements = {
        'H2O': {'H'},
        'H2': {'H'},
        'CO2': {'C'},
        'CO': {'C'},
        'CH4': {'H', 'C'},
        'N2': {'N'},
        'NH3': {'H', 'N'},
        'S2': {'S'},
        'SO2': {'S'},
        'H2S': {'H', 'S'},
        'O2': set(),  # always included for fO2
    }
    # Each noble gas is its own element and species.
    for gas in noble_gases:
        _species_elements[gas] = {gas}

    # Determine which elements have budgets above threshold. Under
    # planet.fO2_source = "from_O_budget" the user O budget
    # is authoritative, so O joins the active set; under user_constant
    # the IW buffer constrains O and we exclude it from the mass
    # constraint here.
    fO2_source = config.planet.fO2_source
    constrained_elements = (
        ('H', 'C', 'N', 'S', 'O') if fO2_source == 'from_O_budget' else ('H', 'C', 'N', 'S')
    )
    for element in constrained_elements + tuple(noble_gases):
        total = float(hf_row.get(f'{element}_kg_total', 0.0))
        if not np.isfinite(total):
            raise ValueError(
                f'{element}_kg_total is not finite ({total}) at the atmodeller call'
            )
    active_elements = set()
    for element in constrained_elements:
        if float(hf_row.get(f'{element}_kg_total', 0.0)) > config.outgas.mass_thresh:
            active_elements.add(element)

    # Noble gases join the active set on any positive inventory. Their budgets
    # are intrinsically trace and sit below the major-volatile mass threshold,
    # so gating them on mass_thresh would drop a realistic noble inventory.
    for gas in noble_gases:
        if float(hf_row.get(f'{gas}_kg_total', 0.0)) > 0.0:
            active_elements.add(gas)

    # A species is included only if ALL its required elements have budgets
    active_species = set()
    for sp, required_elements in _species_elements.items():
        if required_elements.issubset(active_elements):
            active_species.add(sp)

    _atm_gas_species = {
        'H2O': 'H2O',
        'H2': 'H2',
        'CO2': 'CO2',
        'CO': 'CO',
        'CH4': 'CH4',
        'N2': 'N2',
        # atmodeller keys species by Hill-notation formula, so ammonia is
        # 'H3N' and sulfur dioxide is 'O2S' (not 'NH3' / 'SO2'). The wrapper
        # maps the PROTEUS name to atmodeller's name on both the create_gas
        # input and the output-dict / quick_look lookup; a mismatch silently
        # drops that species' atmospheric mass, which under-counts its element
        # (e.g. S in an SO2-dominated oxidising atmosphere).
        'NH3': 'H3N',
        'S2': 'S2',
        'SO2': 'O2S',
        'H2S': 'H2S',
        'O2': 'O2',
    }
    # Noble gases map to their own atmodeller species names.
    for gas in noble_gases:
        _atm_gas_species[gas] = gas

    for proteus_name, atm_name in _atm_gas_species.items():
        if proteus_name not in active_species:
            continue

        # Build kwargs for create_gas
        kwargs = {}

        # Solubility law (if configured)
        sol_key = _sol_map.get(proteus_name)
        if sol_key:
            if sol_key in solubility_models:
                kwargs['solubility'] = solubility_models[sol_key]
            else:
                log.warning(
                    'Solubility model %r not found for %s; using no solubility',
                    sol_key,
                    proteus_name,
                )

        # Real gas EOS (if configured)
        eos_key = _eos_map.get(proteus_name)
        if eos_key:
            if eos_key in eos_models:
                kwargs['activity'] = eos_models[eos_key]
            else:
                log.warning(
                    'EOS model %r not found for %s; using ideal gas', eos_key, proteus_name
                )

        created = ChemicalSpecies.create_gas(atm_name, **kwargs)
        # The wrapper looks up this species' output by f'{atm_name}_g'. atmodeller
        # keys its output by the species' Hill-notation name, so if the mapped
        # atm_name disagrees with the created species' name the output lookup
        # silently misses and this species' atmospheric mass is dropped (for
        # example the SO2 / 'O2S' Hill-name case). Surface a mismatch loudly at
        # the source rather than letting it drop mass downstream.
        expected_key = f'{atm_name}_g'
        if created.name != expected_key:
            log.warning(
                "atmodeller species name %r for %s does not match the wrapper's "
                'output key %r; its atmospheric/dissolved mass will be dropped. '
                'Update _atm_gas_species[%r] to %r.',
                created.name,
                proteus_name,
                expected_key,
                proteus_name,
                created.name.removesuffix('_g'),
            )
        species_list.append(created)

    log.debug(
        'Atmodeller species: %s (active elements: %s)',
        [s.name for s in species_list],
        active_elements,
    )

    # Build planet state
    M_planet = float(hf_row.get('M_planet', 1.0 * M_earth))
    R_int = float(hf_row.get('R_int', 6.371e6))
    T_magma = float(hf_row.get('T_magma', 3000.0))
    Phi_global = float(hf_row.get('Phi_global', 1.0))

    # Do not allow low temperatures: clamp to the configured floor, the
    # same semantics as the calliope entry point.
    if T_magma < config.outgas.T_floor:
        T_magma = float(config.outgas.T_floor)
        log.warning('Outgassing temperature clipped to %.1f K' % T_magma)

    # Add each supported condensate whose elements have budgets (O is buffered or
    # constrained, like O2) and whose thermodynamic data cover T_magma, so the
    # solver never extrapolates a condensed phase outside its fitted range.
    if atm_config.include_condensates:
        for formula in condensate_list:
            try:
                condensate = ChemicalSpecies.create_condensed(formula)
            except Exception as e:
                # Not every condensate is available in all atmodeller versions. Say
                # so, since the run otherwise proceeds without it.
                log.warning(
                    'Could not add condensate %s_cd to the atmodeller network; '
                    'continuing without it: %s',
                    formula,
                    e,
                )
                continue
            required_elements = set(condensate.data.elements) - {'O'}
            if not required_elements.issubset(active_elements):
                continue
            thermo = condensate.data.thermo
            if not min(thermo.T_min) <= T_magma <= max(thermo.T_max):
                continue
            species_list.append(condensate)

    # Core mass fraction from config
    cmf = config.interior_struct.core_frac

    planet = Planet(
        planet_mass=M_planet,
        core_mass_fraction=cmf,
        mantle_melt_fraction=Phi_global,
        surface_radius=R_int,
        temperature=T_magma,
        pressure=np.nan,  # Thin-atmosphere approximation
    )

    # Fugacity and mass constraints. The two paths split here:
    #
    # - user_constant: fO2 is buffered to the configured IW
    #   offset and atmodeller derives the resulting O distribution.
    #   Mass constraints cover H/C/N/S only.
    # - from_O_budget: the user O budget is authoritative and
    #   atmodeller back-solves for the IW offset. Mass constraints cover
    #   all five elements; no fugacity constraint on O2_g.
    fugacity_constraints: dict = {}
    mass_constraints: dict = {}
    for element in constrained_elements:
        key = f'{element}_kg_total'
        mass_kg = float(hf_row.get(key, 0.0))
        if mass_kg > config.outgas.mass_thresh:
            mass_constraints[element] = mass_kg

    # Noble gas mass constraints. A noble gas is included on any positive
    # inventory, matching the active-species selection above, so a trace noble
    # budget is not dropped by the major-volatile threshold.
    for gas in noble_gases:
        mass_kg = float(hf_row.get(f'{gas}_kg_total', 0.0))
        if mass_kg > 0.0:
            mass_constraints[gas] = mass_kg

    if fO2_source == 'user_constant':
        fugacity_constraints['O2_g'] = IronWustiteBuffer(config.outgas.fO2_shift_IW)

    # Stash the authoritative O target for the 'from_O_budget' source:
    # atmodeller's output
    # element_density / dissolved_mass arithmetic equals target['O'] only
    # up to the solver's element residual, and we must not let that
    # residual drift hf_row['O_kg_total'] across iterations (the escape
    # pipeline reads this column to compute each iteration's debit).
    target_O_kg = mass_constraints.get('O') if fO2_source == 'from_O_budget' else None

    if not mass_constraints and fO2_source == 'from_O_budget':
        log.warning('No volatile element budgets above threshold; skipping atmodeller')
        # No solve ran: the atmosphere is empty, and the derived buffer offset and the
        # O residual are undefined (the user pre-seed would misrecord a skipped solve).
        for s in _atm_gas_species:
            kg_total = float(hf_row.get(f'{s}_kg_total', 0.0))
            hf_row[f'{s}_bar'] = 0.0
            hf_row[f'{s}_kg_atm'] = hf_row[f'{s}_kg_liquid'] = hf_row[f'{s}_kg_solid'] = 0.0
            hf_row[f'{s}_mol_atm'] = hf_row[f'{s}_mol_liquid'] = hf_row[f'{s}_mol_solid'] = 0.0
            hf_row[f'{s}_kg_total'] = hf_row[f'{s}_mol_total'] = 0.0
            if s in noble_gases:  # its kg_total is the escape-owned element total
                hf_row[f'{s}_kg_total'] = kg_total
                hf_row[f'{s}_mol_total'] = kg_total / eval_gas_mmw(s)
        hf_row['P_surf'] = 0.0
        _populate_volatile_element_reservoirs(hf_row, element_mmw)
        # No condensate either
        for formula in condensate_list:
            hf_row[f'{formula}_cd_activity'] = 0.0
            hf_row[f'{formula}_cd_kg'] = 0.0
        for el in vol_element_list:
            hf_row[f'{el}_kg_condensed'] = 0.0
        hf_row['fO2_shift_IW_derived'] = float('nan')
        hf_row['O_res'] = float('nan')
        return

    if fO2_source == 'from_O_budget' and 'O' not in mass_constraints:
        raise ValueError(
            'planet.fO2_source = "from_O_budget" requires '
            'hf_row["O_kg_total"] > outgas.mass_thresh, but the O '
            'budget is below threshold. Increase planet.elements.O_budget '
            'or switch to fO2_source = "user_constant".'
        )

    log.debug(
        'Atmodeller solve: T=%.0f K, Phi=%.2f, elements=%s',
        T_magma,
        Phi_global,
        {k: f'{v:.2e}' for k, v in mass_constraints.items()},
    )

    # Build solver parameters from config
    from atmodeller.containers import SolverParameters

    solver_params = SolverParameters(
        atol=config.outgas.solver_atol,
        rtol=config.outgas.solver_rtol,
        max_steps=atm_config.solver_max_steps,
        multistart=atm_config.solver_multistart,
    )

    # Reuse one EquilibriumModel per (species network, fugacity source, solver
    # mode) so atmodeller's JIT solver (cached on the model as ._solver) is
    # compiled once, not rebuilt every solve. The latter two are fixed per run;
    # they are in the key so a reused model can never carry a solver compiled for
    # a different structure. Built here, after the skip/validation guards, so a
    # skipped or rejected call never constructs a model.
    model_signature = (
        tuple(s.name for s in species_list),
        fO2_source,
        atm_config.solver_mode,
    )
    model = _cached_model(
        model_signature,
        lambda: EquilibriumModel(SpeciesNetwork(tuple(species_list))),
    )

    # Solve equilibrium
    try:
        model.solve(
            state=planet,
            fugacity_constraints=fugacity_constraints,
            mass_constraints=mass_constraints,
            solver=atm_config.solver_mode,
            solver_parameters=solver_params,
        )
    except Exception as e:
        log.error('Atmodeller solve failed: %s', e)
        raise

    # Extract results
    output = model.output
    quick_look = output.quick_look()
    total_P = float(np.squeeze(output.total_pressure()))

    log.debug('Atmodeller result: P_total=%.2f bar', total_P)

    # atmodeller's per-species output dict carries the thermodynamically
    # consistent atmospheric (``gas_mass``) and dissolved (``dissolved_mass``)
    # masses. Read it once and reuse for both reservoirs below.
    try:
        output_dict = output.asdict()
    except Exception:
        output_dict = {}

    # An atmodeller release that does not raise on its own keeps the start guess as
    # the output when no attempt converges. Refuse it, so a non-solution never
    # reaches the row, where a held layer would keep its condensate.
    solver_info = output_dict.get('solver')
    status = solver_info.get('status') if isinstance(solver_info, dict) else None
    if status is None:
        log.warning('atmodeller output has no solver status; convergence not checked')
    elif not np.all(np.asarray(status, dtype=bool)):
        masses = ', '.join(f'{el} {kg:.3e}' for el, kg in mass_constraints.items())
        raise RuntimeError(
            f'atmodeller did not converge at T_magma = {T_magma:.0f} K, Phi = '
            f'{Phi_global:.3f} ({masses} kg): no attempt converged within '
            f'solver_multistart = {atm_config.solver_multistart} attempts of '
            f'solver_max_steps = {atm_config.solver_max_steps}. Raise '
            'outgas.atmodeller.solver_multistart or solver_max_steps.'
        )

    # Record condensate saturation on every solve so the helpfile and the log
    # show when each one forms. An activity of 1 means the condensate is stable;
    # below 1 it is undersaturated and its condensed mass is only the solver's
    # floor. Zero marks a solve without that condensate in the network.
    for formula in condensate_list:
        activity = mass = None
        if f'{formula}_cd' in output_dict:
            activity = _species_output_scalar(output_dict, f'{formula}_cd', 'activity')
            mass = _species_output_scalar(output_dict, f'{formula}_cd', 'total_mass')
            # A missing or non-finite value is None, which the float formats reject.
            log.info(
                'Condensate %s: activity=%s, condensed mass=%s kg',
                formula,
                'n/a' if activity is None else f'{activity:.3f}',
                'n/a' if mass is None else f'{mass:.2e}',
            )
        hf_row[f'{formula}_cd_activity'] = float(activity) if activity is not None else 0.0
        hf_row[f'{formula}_cd_kg'] = float(mass) if mass is not None else 0.0

    # Map atmodeller output back to hf_row
    _reverse_map = {v: k for k, v in _atm_gas_species.items()}

    gravity = float(hf_row.get('gravity', 9.81))
    area = 4.0 * np.pi * R_int**2

    # A species outside the solve reports 0, as in CALLIOPE; the loops below fill the others
    for proteus_name in _atm_gas_species:
        hf_row[f'{proteus_name}_bar'] = 0.0
        hf_row[f'{proteus_name}_kg_atm'] = 0.0

    P_total = 0.0
    for atm_name, p_bar in quick_look.items():
        proteus_name = _reverse_map.get(atm_name.replace('_g', ''))
        if proteus_name is None:
            log.debug('Atmodeller quick_look key %r not in species map; skipping', atm_name)
            continue
        if proteus_name not in gas_list:
            continue
        p_val = float(np.squeeze(p_bar))
        hf_row[f'{proteus_name}_bar'] = p_val
        P_total += p_val

        # Atmospheric mass: use atmodeller's internal per-species atmospheric
        # mass (``gas_mass``), consistent with the ``dissolved_mass`` read below
        # so the per-element atm+liquid split closes against the conserved
        # element budget. The thin-atmosphere relation p*1e5*A/g recovers
        # x_i * M_atm, the species' mole-fraction share of the total column,
        # which mis-weights each species by mu_mean/mu_i and breaks the
        # per-element closure; it is kept only as a fallback for the case where
        # the solver output omits gas_mass.
        gas_mass = _species_output_scalar(output_dict, atm_name, 'gas_mass')
        if gas_mass is not None:
            # A small negative solver residual is clamped to zero rather than
            # routed to the fallback (a near-zero species has a near-zero column
            # mass anyway); the fallback is reserved for a genuinely absent
            # gas_mass, i.e. an atmodeller output-schema change.
            hf_row[f'{proteus_name}_kg_atm'] = max(0.0, gas_mass)
        elif gravity > 0 and area > 0:
            log.warning(
                "atmodeller output missing 'gas_mass' for %s; falling back to the "
                'p*A/g thin-atmosphere reconstruction, which mis-weights the '
                'per-element split by mu_mean/mu_species.',
                atm_name,
            )
            hf_row[f'{proteus_name}_kg_atm'] = p_val * 1e5 * area / gravity
        else:
            hf_row[f'{proteus_name}_kg_atm'] = 0.0

    # Total surface pressure
    hf_row['P_surf'] = P_total

    # VMRs
    if P_total > 0:
        for s in gas_list:
            hf_row[f'{s}_vmr'] = float(hf_row.get(f'{s}_bar', 0.0)) / P_total
    else:
        for s in gas_list:
            hf_row[f'{s}_vmr'] = 0.0

    # Dissolved masses from atmodeller output (thermodynamically consistent).
    def _as_mass(value):
        """Squeeze an atmodeller mass or mole-count output to a float, or None on failure."""
        if value is None:
            return None
        try:
            return float(np.squeeze(value))
        except (TypeError, ValueError):
            return None

    # Snapshot the escape-authoritative noble gas element totals before the
    # output overwrites them. For a noble gas the species and the element are
    # the same key, and {gas}_kg_total is owned by the escape pipeline (which
    # debits it each step). The solver's total_mass equals this constraint only
    # to its residual, so the total is restored below rather than overwritten,
    # mirroring how the oxygen element total is preserved.
    noble_total_in = {g: float(hf_row.get(f'{g}_kg_total', 0.0)) for g in noble_gases}

    for proteus_name, atm_name in _atm_gas_species.items():
        if proteus_name not in gas_list:
            continue
        is_noble = proteus_name in noble_gases
        species_key = f'{atm_name}_g'
        species_data = output_dict.get(species_key, {})
        if not isinstance(species_data, dict):
            species_data = {}

        dissolved_kg = _as_mass(species_data.get('dissolved_mass'))
        total_kg = _as_mass(species_data.get('total_mass'))
        dissolved_mol = _as_mass(species_data.get('dissolved_number'))
        total_mol = _as_mass(species_data.get('total_number'))
        gas_mol = _as_mass(species_data.get('gas_number'))

        hf_row[f'{proteus_name}_kg_solid'] = 0.0
        hf_row[f'{proteus_name}_mol_solid'] = 0.0

        if total_kg is not None:
            # gas_mass is the atmospheric write, total - dissolved its closure check
            # and fallback; see the atmospheric mass note in
            # docs/Validation/outgas/atmodeller.md.
            liquid = min(max(0.0, dissolved_kg if dissolved_kg is not None else 0.0), total_kg)
            hf_row[f'{proteus_name}_kg_liquid'] = liquid
            gas_kg = _as_mass(species_data.get('gas_mass'))
            residual_kg = max(0.0, total_kg - liquid)
            if gas_kg is not None:
                hf_row[f'{proteus_name}_kg_atm'] = max(0.0, gas_kg)
                if not np.isclose(max(0.0, gas_kg), residual_kg, rtol=1e-6, atol=1.0):
                    log.warning(
                        'atmodeller gas-phase closure gap for %s: gas_mass=%.6e kg '
                        'but total_mass - dissolved_mass=%.6e kg. A persistent gap '
                        'means the species holds a reservoir this mapping does not '
                        'track (e.g. a solid phase); the solver gas_mass write is '
                        'kept.',
                        atm_name,
                        gas_kg,
                        residual_kg,
                    )
            else:
                hf_row[f'{proteus_name}_kg_atm'] = residual_kg
            # A noble gas element total is escape-owned; restore the pre-solve
            # value so the solver's residual cannot drift the tracked budget
            # across iterations. Its atmosphere and melt split are still taken
            # from the solve.
            hf_row[f'{proteus_name}_kg_total'] = (
                noble_total_in[proteus_name] if is_noble else total_kg
            )

            # Mole counts mirror the mass split above. They are read directly
            # from atmodeller's own number-of-moles output rather than derived
            # via a molar mass, since a noble-gas element total is restored
            # from the escape-owned kg budget above (not atmodeller's own
            # total_mass) and mol_total must stay consistent with that.
            if total_mol is not None:
                liquid_mol = min(
                    max(0.0, dissolved_mol if dissolved_mol is not None else 0.0), total_mol
                )
                hf_row[f'{proteus_name}_mol_liquid'] = liquid_mol
                hf_row[f'{proteus_name}_mol_atm'] = (
                    max(0.0, gas_mol)
                    if gas_mol is not None
                    else max(0.0, total_mol - liquid_mol)
                )
                hf_row[f'{proteus_name}_mol_total'] = (
                    hf_row[f'{proteus_name}_kg_total'] / eval_gas_mmw(proteus_name)
                    if is_noble
                    else total_mol
                )
            else:
                # atmodeller reported a mass but no mole count for this species;
                # derive moles from the kg values just written above, matching
                # the sibling branches' fallback convention.
                species_mmw = eval_gas_mmw(proteus_name)
                hf_row[f'{proteus_name}_mol_liquid'] = (
                    hf_row[f'{proteus_name}_kg_liquid'] / species_mmw
                )
                hf_row[f'{proteus_name}_mol_atm'] = (
                    hf_row[f'{proteus_name}_kg_atm'] / species_mmw
                )
                hf_row[f'{proteus_name}_mol_total'] = (
                    hf_row[f'{proteus_name}_kg_total'] / species_mmw
                )
        elif is_noble:
            # An inactive noble gas is not in the solve: its atmosphere and melt
            # columns are 0, and its escape-owned element total stays untouched.
            hf_row[f'{proteus_name}_kg_liquid'] = 0.0
            hf_row[f'{proteus_name}_kg_atm'] = 0.0
            hf_row[f'{proteus_name}_bar'] = 0.0
            hf_row[f'{proteus_name}_vmr'] = 0.0
            hf_row[f'{proteus_name}_mol_liquid'] = 0.0
            hf_row[f'{proteus_name}_mol_atm'] = 0.0
            hf_row[f'{proteus_name}_mol_total'] = float(
                hf_row.get(f'{proteus_name}_kg_total', 0.0)
            ) / eval_gas_mmw(proteus_name)
        else:
            # Reactive species without total_mass: kg_total = kg_atm + kg_liquid, with
            # kg_liquid 0 without dissolved_mass (both 0 outside the solve).
            hf_row[f'{proteus_name}_kg_liquid'] = (
                max(0.0, dissolved_kg) if dissolved_kg is not None else 0.0
            )
            hf_row[f'{proteus_name}_kg_total'] = float(
                hf_row.get(f'{proteus_name}_kg_atm', 0.0)
            ) + float(hf_row[f'{proteus_name}_kg_liquid'])

            # Moles from the kg values above with the species' own molar mass, as dummy.py does.
            species_mmw = eval_gas_mmw(proteus_name)
            hf_row[f'{proteus_name}_mol_liquid'] = (
                hf_row[f'{proteus_name}_kg_liquid'] / species_mmw
            )
            hf_row[f'{proteus_name}_mol_atm'] = (
                float(hf_row.get(f'{proteus_name}_kg_atm', 0.0)) / species_mmw
            )
            hf_row[f'{proteus_name}_mol_total'] = (
                hf_row[f'{proteus_name}_kg_total'] / species_mmw
            )

    # Mean molecular weight (approximate from VMRs). The noble gases are
    # gas_list members with their own VMRs, so they enter this sum directly;
    # their molar mass is the element molar mass since the species and the
    # element are the same monatomic entity.
    _mmw = {
        'H2O': 18.015e-3,
        'CO2': 44.01e-3,
        'H2': 2.016e-3,
        'CO': 28.01e-3,
        'CH4': 16.04e-3,
        'N2': 28.01e-3,
        'S2': 64.13e-3,
        'SO2': 64.07e-3,
        'O2': 32.0e-3,
        'H2S': 34.08e-3,
        'NH3': 17.03e-3,
        **{gas: element_mmw[gas] for gas in noble_gases},
    }
    mmw = sum(float(hf_row.get(f'{s}_vmr', 0.0)) * _mmw.get(s, 28.0e-3) for s in gas_list)
    if mmw > 0:
        hf_row['atm_kg_per_mol'] = mmw

    log.debug(
        'Atmodeller: P_surf=%.2f bar, MMW=%.3f g/mol',
        P_total,
        mmw * 1e3,
    )

    # Populate the per-element reservoir split ({e}_kg_atm/liquid/solid) from the
    # per-species inventory so the escape step (which reads {e}_kg_atm under
    # reservoir="outgas") keeps a live element budget under atmodeller, as
    # CALLIOPE writes natively. The per-element {e}_kg_total slot is left to its
    # owner (escape for the running budget; the authoritative O_kg_total write
    # below for oxygen), so escape's debited total is never overwritten here.
    _populate_volatile_element_reservoirs(hf_row, element_mmw)

    # Condensed-phase mass (graphite, sulfur, ... when include_condensates is on)
    # goes to its own per-element reservoir, kept apart from the solid mantle. The
    # gas-species inventory above does not carry condensates, so the split
    # atm + liquid + solid + condensed closes against the element budget only with
    # this term. atmodeller reports the condensed mass per element; it is zero for
    # the elements (and conditions) that do not condense, and the columns stay 0
    # when condensates are disabled.
    for el in vol_element_list:
        condensed = None
        if atm_config.include_condensates:
            condensed = _species_output_scalar(output_dict, f'element_{el}', 'condensed_mass')
        hf_row[f'{el}_kg_condensed'] = max(0.0, condensed) if condensed else 0.0

    for el in ('H', 'C', 'N', 'S'):
        kg = (hf_row.get(f'{el}_kg_atm', 0.0), hf_row.get(f'{el}_kg_liquid', 0.0))
        if el not in active_elements and any(kg):
            raise RuntimeError(
                f'{el} is outside the atmodeller solve but {el}_kg_atm, {el}_kg_liquid = '
                f'{kg[0]:.6e}, {kg[1]:.6e} kg'
            )
    if not mass_constraints and hf_row['P_surf'] != hf_row['O2_bar']:
        raise RuntimeError(
            f'O2-only solve: P_surf {hf_row["P_surf"]} != O2_bar {hf_row["O2_bar"]}'
        )

    # Total volatile oxygen (atmospheric + melt-dissolved), summed over the
    # O-bearing volatile species. Computed for every fO2 source so the
    # whole-planet O accounting (M_ele) includes atmospheric O whether or not
    # O was a mass constraint; mirrors CALLIOPE, which always writes O_kg_total.
    volatile_O_kg = _total_volatile_oxygen_kg(hf_row, element_mmw)
    # Oxygen held in the O-bearing condensates (water, sulfuric acid) this solve
    # formed, written to O_kg_condensed above; zero without them.
    condensed_O_kg = float(hf_row.get('O_kg_condensed', 0.0))

    # Plumb the derived IW-buffer offset only for the 'from_O_budget'
    # source. Under
    # user_constant the wrapper's pre-dispatch echo of
    # config.outgas.fO2_shift_IW IS the offset the chemistry equilibrated
    # to (the fugacity constraint set it), so overwriting with
    # atmodeller's back-computed log10dIW_1_bar would lose bit-for-bit
    # echo of the user input and introduce a small Hirschmann-buffer
    # 1-bar-vs-P reconstruction drift. The IW buffer atmodeller uses
    # for the 'from_O_budget' source is Hirschmann combined; CALLIOPE uses O'Neill & Eggins
    # 2002 (~0.95 dex offset at 3000 K). Cross-backend comparison work
    # later in the framework will quantify this; for now the column is
    # backend-faithful (each wrapper reports its own buffer).
    if fO2_source == 'from_O_budget':
        o2_dict = output_dict.get('O2_g', {})
        log10dIW = o2_dict.get('log10dIW_1_bar') if isinstance(o2_dict, dict) else None
        if log10dIW is None:
            log.warning(
                "atmodeller output missing 'O2_g.log10dIW_1_bar' for the "
                "'from_O_budget' source; "
                'fO2_shift_IW_derived left at pre-dispatch default '
                '(%.3f) which is meaningless when O was a mass constraint.',
                hf_row.get('fO2_shift_IW_derived', float('nan')),
            )
        else:
            try:
                hf_row['fO2_shift_IW_derived'] = float(np.squeeze(log10dIW))
            except (TypeError, ValueError) as e:
                log.warning(
                    'atmodeller log10dIW_1_bar JAX conversion failed (%s); '
                    'fO2_shift_IW_derived left at pre-dispatch default '
                    '(%.3f) which is meaningless when O was a mass constraint.',
                    e,
                    hf_row.get('fO2_shift_IW_derived', float('nan')),
                )
        # O mass-balance residual: the solved O (atmospheric + dissolved +
        # condensed, computed above) minus the authoritative user target, which
        # constrains all three. Reported in kg to match CALLIOPE's
        # H/C/N/S/O_res convention.
        hf_row['O_res'] = volatile_O_kg + condensed_O_kg - float(target_O_kg)

        # Restore the authoritative user O budget. atmodeller's solver may
        # have converged to atm_O+liq_O slightly off from target_O_kg (within
        # its atol/rtol); the per-iteration drift would otherwise accumulate
        # in the escape pipeline that reads hf_row['O_kg_total'] to compute
        # the next debit.
        hf_row['O_kg_total'] = float(target_O_kg)
    else:
        # Buffered-fO2 sources (e.g. user_constant): O is not a mass
        # constraint, so the equilibrium-derived total volatile O is the
        # tracked element budget. Without this, M_ele / M_planet omit
        # atmospheric O entirely (e.g. the O2 atmosphere at oxidising IW),
        # which breaks the M_atm <= M_planet accounting at high fO2. The
        # escape pipeline's between-step debit of O_kg_total is superseded
        # here each call, exactly as in CALLIOPE under user_constant: the
        # buffered equilibrium re-derives O from the current H/C/S
        # inventory, which itself falls as escape depletes those budgets.
        # The O in condensates belongs to the budget as well, so the total
        # closes against atm + liquid + solid + condensed.
        hf_row['O_kg_total'] = volatile_O_kg + condensed_O_kg
        # No authoritative O target under a buffered fO2 source, so the O
        # mass-balance residual is zero by construction. Written explicitly
        # to match CALLIOPE rather than relying on the pre-dispatch seed.
        hf_row['O_res'] = 0.0