Phase properties and derivatives

Homogeneous-state properties

torch_flash.properties

Homogeneous-state thermodynamic property calculations.

PhaseResponseDerivatives dataclass

PhaseResponseDerivatives(
    molar_volume: Tensor,
    isothermal_compressibility: Tensor,
    thermal_expansion_coefficient: Tensor,
    isothermal_compressibility_temperature_derivative: Tensor,
    thermal_expansion_temperature_derivative: Tensor,
)

Autodiff response functions used by the two derivative criteria.

Attributes:

Name	Type	Description
molar_volume	Tensor	Homogeneous physical molar volume in m3/mol.
isothermal_compressibility	Tensor	kappa = -(1/V)(dV/dP)_T in 1/Pa.
thermal_expansion_coefficient	Tensor	alpha = (1/V)(dV/dT)_P in 1/K.
isothermal_compressibility_temperature_derivative	Tensor	(d kappa/dT)_P in 1/(Pa K).
thermal_expansion_temperature_derivative	Tensor	(d alpha/dT)_P in 1/K^2.

Notes

The derivatives are assembled from first and second PyTorch autodiff derivatives of the model's explicit P(T, V, x) function. The supplied homogeneous root is evaluated once per state; no finite-difference step is used. Matching leading state batches are supported.

ThermodynamicDerivatives dataclass

ThermodynamicDerivatives(
    dfugacity_coefficient_dlogits: Tensor,
    dfugacity_coefficient_dindependent_composition: Tensor,
    dfugacity_coefficient_dtemperature: Tensor,
    dfugacity_coefficient_dpressure: Tensor,
    dfugacity_coefficient_dmoles: Tensor,
    dlog_fugacity_coefficient_dlogits: Tensor,
    dlog_fugacity_coefficient_dindependent_composition: Tensor,
    dlog_fugacity_coefficient_dtemperature: Tensor,
    dlog_fugacity_coefficient_dpressure: Tensor,
    dlog_fugacity_coefficient_dmoles: Tensor,
    dfugacity_dlogits: Tensor,
    dfugacity_dindependent_composition: Tensor,
    dfugacity_dtemperature: Tensor,
    dfugacity_dpressure: Tensor,
    dfugacity_dmoles: Tensor,
    dlog_fugacity_dlogits: Tensor,
    dlog_fugacity_dindependent_composition: Tensor,
    dlog_fugacity_dtemperature: Tensor,
    dlog_fugacity_dpressure: Tensor,
    dlog_fugacity_dmoles: Tensor,
    dchemical_potential_dlogits: Tensor,
    dchemical_potential_dindependent_composition: Tensor,
    dchemical_potential_dtemperature: Tensor,
    dchemical_potential_dpressure: Tensor,
    dchemical_potential_dmoles: Tensor,
    dreduced_chemical_potential_dlogits: Tensor,
    dreduced_chemical_potential_dindependent_composition: Tensor,
    dreduced_chemical_potential_dtemperature: Tensor,
    dreduced_chemical_potential_dpressure: Tensor,
    dreduced_chemical_potential_dmoles: Tensor,
    dmolar_volume_dlogits: Tensor,
    dmolar_volume_dindependent_composition: Tensor,
    dmolar_volume_dtemperature: Tensor,
    dmolar_volume_dpressure: Tensor,
    dmolar_volume_dmoles: Tensor,
    dgibbs_dtemperature: Tensor,
    dgibbs_dpressure: Tensor,
)

First homogeneous-state derivatives.

Component derivatives are provided in both unconstrained softmax-logit coordinates and n-1 independent mole fractions, where the last mole fraction is 1 - sum(x_independent). Temperature and pressure derivatives hold composition fixed.

Mole-number derivatives hold T and P fixed and are evaluated at n_i = x_i mol (a one-mole total basis). Because all returned properties are intensive, the corresponding derivative at another total amount is this value divided by that amount in mol.

Fugacity derivatives have Pa per coordinate units. dlog_fugacity_* differentiates the dimensionless ln(f_i / p_standard). Chemical-potential derivatives have J/mol per coordinate units; the reduced variants differentiate mu_i / (R*T). Molar-volume derivatives have m3/mol per coordinate units.

Attributes:

Name	Type	Description
dfugacity_coefficient_dlogits	Tensor	Jacobian of dimensionless phi_i with respect to unconstrained softmax logits.
dfugacity_coefficient_dindependent_composition	Tensor	Jacobian of phi_i with respect to the first n-1 mole fractions.
dfugacity_coefficient_dtemperature, dfugacity_coefficient_dpressure		Fixed-composition derivatives of phi_i in 1/K and 1/Pa.
dfugacity_coefficient_dmoles	Tensor	Unit-total-mole-basis derivative of phi_i with respect to moles.
dlog_fugacity_coefficient_dlogits	Tensor	Jacobian of ln(phi_i) with respect to logits.
dlog_fugacity_coefficient_dindependent_composition	Tensor	Jacobian of ln(phi_i) with respect to independent mole fractions.
dlog_fugacity_coefficient_dtemperature, dlog_fugacity_coefficient_dpressure		Fixed-composition derivatives of ln(phi_i).
dlog_fugacity_coefficient_dmoles	Tensor	Unit-basis mole-number derivative of ln(phi_i).
dfugacity_dlogits, dfugacity_dindependent_composition		Fugacity Jacobians in Pa with respect to the two composition coordinate systems.
dfugacity_dtemperature, dfugacity_dpressure, dfugacity_dmoles		Fugacity derivatives in Pa/K, dimensionless Pa/Pa, and Pa/mol on the documented one-mole basis.
dlog_fugacity_dlogits, dlog_fugacity_dindependent_composition		Jacobians of dimensionless ln(f_i / p_standard).
dlog_fugacity_dtemperature, dlog_fugacity_dpressure, dlog_fugacity_dmoles		Temperature, pressure, and unit-basis mole-number derivatives of dimensionless log fugacity.
dchemical_potential_dlogits, dchemical_potential_dindependent_composition		Chemical-potential Jacobians in J/mol per composition coordinate.
dchemical_potential_dtemperature, dchemical_potential_dpressure		Fixed-composition derivatives in J/(mol K) and J/(mol Pa).
dchemical_potential_dmoles	Tensor	Unit-basis mole-number Jacobian in J/mol2.
dreduced_chemical_potential_dlogits	Tensor	Jacobian of mu_i/(RT) with respect to logits.
dreduced_chemical_potential_dindependent_composition	Tensor	Jacobian of mu_i/(RT) with respect to independent mole fractions.
dreduced_chemical_potential_dtemperature	Tensor	Fixed-composition derivative of mu_i/(RT) in 1/K.
dreduced_chemical_potential_dpressure	Tensor	Fixed-composition derivative of mu_i/(RT) in 1/Pa.
dreduced_chemical_potential_dmoles	Tensor	Unit-basis mole-number Jacobian of mu_i/(RT).
dmolar_volume_dlogits, dmolar_volume_dindependent_composition		Molar-volume Jacobians in m3/mol per composition coordinate.
dmolar_volume_dtemperature, dmolar_volume_dpressure		Fixed-composition derivatives in m3/(mol K) and m3/(mol Pa).
dmolar_volume_dmoles	Tensor	Unit-basis mole-number derivative in m3/mol2.
dgibbs_dtemperature, dgibbs_dpressure		Fixed-composition total molar Gibbs derivatives in J/(mol K) and J/(mol Pa).

Notes

Composition Jacobians place the property/component output axis before the coordinate axis. Scalar temperature and pressure derivatives retain only the property/component output shape.

ThermalProperties dataclass

ThermalProperties(
    molar_enthalpy: Tensor,
    molar_internal_energy: Tensor,
    molar_entropy: Tensor,
    molar_helmholtz_energy: Tensor,
    molar_gibbs_energy: Tensor,
    reduced_helmholtz_energy: Tensor,
    reduced_gibbs_energy: Tensor,
    reduced_residual_helmholtz_energy: Tensor,
    reduced_residual_gibbs_energy: Tensor,
    isobaric_heat_capacity: Tensor,
    isochoric_heat_capacity: Tensor,
    joule_thomson_coefficient: Tensor,
    speed_of_sound: Tensor | None,
    residual_enthalpy: Tensor,
    residual_entropy: Tensor,
)

Caloric and response properties of one homogeneous phase.

All quantities are molar SI except the Joule-Thomson coefficient (K/Pa) and speed of sound (m/s). reduced_* free energies are dimensionless molar quantities divided by R*T.

Attributes:

Name	Type	Description
molar_enthalpy, molar_internal_energy		Total energies in J/mol.
molar_entropy	Tensor	Total entropy in J/(mol K).
molar_helmholtz_energy, molar_gibbs_energy		Total free energies in J/mol.
reduced_helmholtz_energy, reduced_gibbs_energy		Dimensionless total free energies divided by RT.
reduced_residual_helmholtz_energy, reduced_residual_gibbs_energy		Dimensionless EOS departure free energies.
isobaric_heat_capacity, isochoric_heat_capacity		C_p and C_v in J/(mol K).
joule_thomson_coefficient	Tensor	Isenthalpic pressure response in K/Pa.
speed_of_sound	Tensor | None	Speed in m/s, or None when molar mass is unavailable.
residual_enthalpy	Tensor	EOS departure enthalpy in J/mol.
residual_entropy	Tensor	EOS departure entropy in J/(mol K).

identify_flash_phases

identify_flash_phases(
    phases: tuple[PhaseProperties, ...],
    *,
    ambiguity_relative_tolerance: float = DEFAULT_AMBIGUITY_RELATIVE_TOLERANCE,
) -> tuple[PhaseProperties, ...]

Fill wholly unavailable multiphase identities by molar-volume ordering.

Parameters:

Name	Type	Description	Default
phases	tuple[PhaseProperties, ...]	Equilibrium homogeneous-phase properties for one flash result. Every molar volume must be a finite positive scalar when density ordering is required. Existing identifications produced by a named method are treated as authoritative and are not mixed with this fallback.	required
ambiguity_relative_tolerance	float	Non-negative relative separation required between the largest and second-largest molar volumes. If their ratio is no greater than 1 + ambiguity_relative_tolerance, all fallback identities are reported as unknown.	DEFAULT_AMBIGUITY_RELATIVE_TOLERANCE

Returns:

Type	Description
tuple[PhaseProperties, ...]	The original tuple when it is empty, contains one phase, or at least one phase already has an available named-method identification. Otherwise, copies of the supplied properties carry method="density-ordering" identifications. The unique largest-volume phase is vapor and the remaining phases are liquid when the volume separation is unambiguous.

Raises:

Type	Description
ValueError	If the ambiguity tolerance is negative or nonfinite, or if fallback ordering is required but any molar volume is batched, nonfinite, or nonpositive.

Notes

This is a conservative last-resort labeling heuristic for phases already returned by a converged flash. It does not alter compositions or phase fractions, solve equilibrium, distinguish multiple liquid types, or prove that the largest-volume phase is thermodynamically vapor.

Source code in src/torch_flash/properties/phase_identification.py

def identify_flash_phases(
    phases: tuple[PhaseProperties, ...],
    *,
    ambiguity_relative_tolerance: float = DEFAULT_AMBIGUITY_RELATIVE_TOLERANCE,
) -> tuple[PhaseProperties, ...]:
    """Fill wholly unavailable multiphase identities by molar-volume ordering.

    Parameters
    ----------
    phases
        Equilibrium homogeneous-phase properties for one flash result. Every
        molar volume must be a finite positive scalar when density ordering is
        required. Existing identifications produced by a named method are
        treated as authoritative and are not mixed with this fallback.
    ambiguity_relative_tolerance
        Non-negative relative separation required between the largest and
        second-largest molar volumes. If their ratio is no greater than
        ``1 + ambiguity_relative_tolerance``, all fallback identities are
        reported as unknown.

    Returns
    -------
    tuple[PhaseProperties, ...]
        The original tuple when it is empty, contains one phase, or at least
        one phase already has an available named-method identification.
        Otherwise, copies of the supplied properties carry
        ``method="density-ordering"`` identifications. The unique
        largest-volume phase is vapor and the remaining phases are liquid when
        the volume separation is unambiguous.

    Raises
    ------
    ValueError
        If the ambiguity tolerance is negative or nonfinite, or if fallback
        ordering is required but any molar volume is batched, nonfinite, or
        nonpositive.

    Notes
    -----
    This is a conservative last-resort labeling heuristic for phases already
    returned by a converged flash. It does not alter compositions or phase
    fractions, solve equilibrium, distinguish multiple liquid types, or prove
    that the largest-volume phase is thermodynamically vapor.
    """
    _validate_options(DEFAULT_VOLUME_TO_COVOLUME_THRESHOLD, ambiguity_relative_tolerance)
    if not phases:
        return phases
    has_named_identification = tuple(
        phase.phase_identification is not None
        and phase.phase_identification.method != "unavailable"
        for phase in phases
    )
    if any(has_named_identification):
        return phases
    if len(phases) == 1:
        return phases

    volumes = torch.stack(tuple(phase.molar_volume for phase in phases))
    if volumes.ndim != 1 or not bool(torch.isfinite(volumes).all()) or bool((volumes <= 0.0).any()):
        raise ValueError("phase molar volumes must be finite positive scalars")
    sorted_volumes, _ = torch.sort(volumes)
    vapor_volume = sorted_volumes[-1]
    next_volume = sorted_volumes[-2]
    separator = torch.sqrt(vapor_volume * next_volume)
    separation = vapor_volume / next_volume
    ambiguous = bool(separation <= 1.0 + ambiguity_relative_tolerance)

    identified = []
    for phase, volume in zip(phases, volumes, strict=True):
        if ambiguous:
            kind: PhaseIdentityKind = "unknown"
        else:
            kind = "vapor" if bool(volume > separator) else "liquid"
        identification = PhaseIdentification(
            kind=kind,
            method="density-ordering",
            criterion_value=volume,
            threshold=separator,
            ambiguous=ambiguous,
        )
        identified.append(replace(phase, phase_identification=identification))
    return tuple(identified)

identify_phase

identify_phase(
    model: object,
    state: ChemicalState,
    phase: PhaseKind = "stable",
    *,
    method: PhaseIdentificationCriterion = DEFAULT_PHASE_IDENTIFICATION_METHOD,
    threshold: float = DEFAULT_VOLUME_TO_COVOLUME_THRESHOLD,
    pseudo_critical_temperature_factor: float = DEFAULT_PSEUDO_CRITICAL_TEMPERATURE_FACTOR,
    ambiguity_relative_tolerance: float = DEFAULT_AMBIGUITY_RELATIVE_TOLERANCE,
    pip_denominator_relative_tolerance: float | None = None,
) -> PhaseIdentification

Identify one homogeneous state with a selected literature criterion.

Parameters:

Name	Type	Description	Default
model	object	Homogeneous-state model exposing the properties required by method. The derivative methods require differentiable pressure(T, V, x) and molar_volume(T, P, x, phase) methods.	required
state	ChemicalState	One scalar temperature in K, scalar pressure in Pa, and one normalized composition vector.	required
phase	PhaseKind	Homogeneous EoS root-selection request. It is separate from the returned physical identity and does not perform a flash.	'stable'
method	PhaseIdentificationCriterion	Literature criterion used for the label. The default remains Pedersen's V/b rule for backward compatibility. The Venkatarathnam-Oellrich method compares its dimensionless phase-identification parameter with one.	DEFAULT_PHASE_IDENTIFICATION_METHOD
threshold	float	Positive dimensionless V/b threshold used only by Pedersen's method.	DEFAULT_VOLUME_TO_COVOLUME_THRESHOLD
pseudo_critical_temperature_factor	float	Positive multiplier in Li's volume-weighted pseudo-critical temperature, used only by that method.	DEFAULT_PSEUDO_CRITICAL_TEMPERATURE_FACTOR
ambiguity_relative_tolerance	float	Non-negative relative ambiguity band around the selected separator. The returned liquid/vapor label still follows the exact inequality; this value changes only the ambiguous flag.	DEFAULT_AMBIGUITY_RELATIVE_TOLERANCE
pip_denominator_relative_tolerance	float | None	Optional scaled singularity guard forwarded to :func:phase_identification_parameter. It is ignored by the other methods. None selects the square root of machine epsilon.	None

Returns:

Type	Description
PhaseIdentification	Physical identity, method, native criterion tensor, decision threshold, and ambiguity flag. A method whose required model properties are unavailable returns kind="unknown". A singular or mechanically unstable PIP state also returns an ambiguous unknown result rather than an infinite criterion.

Raises:

Type	Description
ValueError	If the state is batched, a numerical option is invalid, the method name is unknown, or an evaluated non-PIP criterion is invalid.
TypeError	If a mandatory model operation for the selected available criterion is missing.

Notes

Every returned criterion tensor remains connected to the PyTorch graph; only the string label and ambiguity flag require a detached scalar decision. Phase identification is post-processing: it neither establishes stability nor changes phase compositions, fractions, or phase count.

Source code in src/torch_flash/properties/phase_identification.py

def identify_phase(
    model: object,
    state: ChemicalState,
    phase: PhaseKind = "stable",
    *,
    method: PhaseIdentificationCriterion = DEFAULT_PHASE_IDENTIFICATION_METHOD,
    threshold: float = DEFAULT_VOLUME_TO_COVOLUME_THRESHOLD,
    pseudo_critical_temperature_factor: float = DEFAULT_PSEUDO_CRITICAL_TEMPERATURE_FACTOR,
    ambiguity_relative_tolerance: float = DEFAULT_AMBIGUITY_RELATIVE_TOLERANCE,
    pip_denominator_relative_tolerance: float | None = None,
) -> PhaseIdentification:
    """Identify one homogeneous state with a selected literature criterion.

    Parameters
    ----------
    model
        Homogeneous-state model exposing the properties required by
        ``method``. The derivative methods require differentiable
        ``pressure(T, V, x)`` and ``molar_volume(T, P, x, phase)`` methods.
    state
        One scalar temperature in K, scalar pressure in Pa, and one normalized
        composition vector.
    phase
        Homogeneous EoS root-selection request. It is separate from the
        returned physical identity and does not perform a flash.
    method
        Literature criterion used for the label. The default remains
        Pedersen's ``V/b`` rule for backward compatibility. The
        Venkatarathnam-Oellrich method compares its dimensionless
        phase-identification parameter with one.
    threshold
        Positive dimensionless ``V/b`` threshold used only by Pedersen's
        method.
    pseudo_critical_temperature_factor
        Positive multiplier in Li's volume-weighted pseudo-critical
        temperature, used only by that method.
    ambiguity_relative_tolerance
        Non-negative relative ambiguity band around the selected separator.
        The returned liquid/vapor label still follows the exact inequality;
        this value changes only the ``ambiguous`` flag.
    pip_denominator_relative_tolerance
        Optional scaled singularity guard forwarded to
        :func:`phase_identification_parameter`. It is ignored by the other
        methods. ``None`` selects the square root of machine epsilon.

    Returns
    -------
    PhaseIdentification
        Physical identity, method, native criterion tensor, decision
        threshold, and ambiguity flag. A method whose required model
        properties are unavailable returns ``kind="unknown"``. A singular or
        mechanically unstable PIP state also returns an ambiguous unknown
        result rather than an infinite criterion.

    Raises
    ------
    ValueError
        If the state is batched, a numerical option is invalid, the method
        name is unknown, or an evaluated non-PIP criterion is invalid.
    TypeError
        If a mandatory model operation for the selected available criterion
        is missing.

    Notes
    -----
    Every returned criterion tensor remains connected to the PyTorch graph;
    only the string label and ambiguity flag require a detached scalar
    decision. Phase identification is post-processing: it neither establishes
    stability nor changes phase compositions, fractions, or phase count.
    """
    _validate_options(threshold, ambiguity_relative_tolerance)
    _validate_positive_factor(pseudo_critical_temperature_factor)
    _validate_pip_denominator_tolerance(pip_denominator_relative_tolerance)
    if state.temperature.ndim != 0 or state.pressure.ndim != 0 or state.composition.ndim != 1:
        raise ValueError("phase identification currently accepts one scalar T-P state")

    if method == "pedersen-volume-to-covolume":
        select_z: Any = getattr(model, "select_z", None)
        if not callable(select_z):
            raise TypeError("phase-identification model must implement select_z")
        compressibility_factor = cast(
            Tensor,
            select_z(
                state.temperature,
                state.pressure,
                state.composition,
                phase,
            ),
        )
        return _from_state_values(
            model,
            state,
            compressibility_factor,
            threshold=threshold,
            ambiguity_relative_tolerance=ambiguity_relative_tolerance,
        )

    if method == "li-pseudo-critical-temperature":
        if (
            _model_tensor(model, "critical_temperature") is None
            or _model_tensor(model, "critical_volume") is None
        ):
            return _unavailable_identification()
        pseudo_critical = li_pseudo_critical_temperature(
            model,
            state.composition,
            factor=pseudo_critical_temperature_factor,
        )
        ratio = state.temperature / pseudo_critical
        one = ratio.new_tensor(1.0)
        return _scalar_identification(
            criterion=ratio,
            threshold=one,
            method=method,
            vapor_if_positive=True,
            ambiguity_margin=torch.log(ratio),
            ambiguity_limit=torch.log1p(ratio.new_tensor(ambiguity_relative_tolerance)),
        )

    if method == "perschke-negative-flash":
        required = ("critical_temperature", "critical_pressure", "acentric_factor")
        if not all(_model_tensor(model, name) is not None for name in required):
            return _unavailable_identification()
        residual = negative_flash_residual(model, state)
        scale = _negative_flash_scale(model, state)
        zero = residual.new_zeros(())
        return _scalar_identification(
            criterion=residual,
            threshold=zero,
            method=method,
            vapor_if_positive=True,
            ambiguity_margin=residual,
            ambiguity_limit=scale * ambiguity_relative_tolerance,
        )

    if method == "venkatarathnam-oellrich-phase-identification-parameter":
        if not callable(getattr(model, "pressure", None)) or not callable(
            getattr(model, "molar_volume", None)
        ):
            return _unavailable_identification()
        one = state.temperature.new_tensor(1.0)
        try:
            parameter = phase_identification_parameter(
                model,
                state,
                phase,
                denominator_relative_tolerance=pip_denominator_relative_tolerance,
            )
        except _PhaseIdentificationParameterSingularityError:
            return PhaseIdentification(
                kind="unknown",
                method=method,
                criterion_value=None,
                threshold=one,
                ambiguous=True,
            )
        return _scalar_identification(
            criterion=parameter,
            threshold=one,
            method=method,
            vapor_if_positive=False,
            ambiguity_margin=parameter - one,
            ambiguity_limit=one * ambiguity_relative_tolerance,
        )

    if method in (
        "pasad-isothermal-compressibility-derivative",
        "bennett-thermal-expansion-derivative",
    ):
        if not callable(getattr(model, "pressure", None)) or not callable(
            getattr(model, "molar_volume", None)
        ):
            return _unavailable_identification()
        response = _phase_response_details(model, state, phase)
        if method == "pasad-isothermal-compressibility-derivative":
            criterion = response.values.isothermal_compressibility_temperature_derivative
            scale = response.compressibility_derivative_scale
        else:
            criterion = response.values.thermal_expansion_temperature_derivative
            scale = response.expansion_derivative_scale
        zero = criterion.new_zeros(())
        return _scalar_identification(
            criterion=criterion,
            threshold=zero,
            method=method,
            vapor_if_positive=False,
            ambiguity_margin=criterion,
            ambiguity_limit=scale * ambiguity_relative_tolerance,
        )

    raise ValueError(f"unknown phase-identification method {method!r}")

li_pseudo_critical_temperature

li_pseudo_critical_temperature(
    model: object,
    composition: Tensor,
    *,
    factor: float = DEFAULT_PSEUDO_CRITICAL_TEMPERATURE_FACTOR,
    critical_volume: Tensor | None = None,
) -> Tensor

Return Li's volume-weighted pseudo-critical temperature.

Parameters:

Name	Type	Description	Default
model	object	Model exposing per-component critical temperatures and, unless supplied explicitly, critical molar volumes.	required
composition	Tensor	Mole fractions with components on the final axis.	required
factor	float	Dimensionless tuning factor r1. Bennett and Schmidt use one.	DEFAULT_PSEUDO_CRITICAL_TEMPERATURE_FACTOR
critical_volume	Tensor | None	Optional per-component critical molar volumes in m3/mol.	None

Returns:

Type	Description
Tensor	Pseudo-critical temperature in K with the composition batch shape.

Notes

This implements Tc = r1*sum(xj*Vcj*Tcj)/sum(xj*Vcj). A phase is labeled vapor when its actual temperature exceeds this value.

Source code in src/torch_flash/properties/phase_identification.py

def li_pseudo_critical_temperature(
    model: object,
    composition: Tensor,
    *,
    factor: float = DEFAULT_PSEUDO_CRITICAL_TEMPERATURE_FACTOR,
    critical_volume: Tensor | None = None,
) -> Tensor:
    """Return Li's volume-weighted pseudo-critical temperature.

    Parameters
    ----------
    model
        Model exposing per-component critical temperatures and, unless
        supplied explicitly, critical molar volumes.
    composition
        Mole fractions with components on the final axis.
    factor
        Dimensionless tuning factor ``r1``. Bennett and Schmidt use one.
    critical_volume
        Optional per-component critical molar volumes in m3/mol.

    Returns
    -------
    Tensor
        Pseudo-critical temperature in K with the composition batch shape.

    Notes
    -----
    This implements ``Tc = r1*sum(xj*Vcj*Tcj)/sum(xj*Vcj)``. A phase is
    labeled vapor when its actual temperature exceeds this value.
    """
    _validate_positive_factor(factor)
    critical_temperature = _model_tensor(model, "critical_temperature")
    if critical_temperature is None:
        raise TypeError("pseudo-critical-temperature model lacks critical temperatures")
    volume = _model_tensor(model, "critical_volume") if critical_volume is None else critical_volume
    if volume is None:
        raise TypeError("pseudo-critical-temperature method requires critical volumes")
    if (
        critical_temperature.ndim != 1
        or volume.ndim != 1
        or critical_temperature.shape != volume.shape
    ):
        raise ValueError("critical temperatures and volumes must be equal-length vectors")
    if not bool(torch.isfinite(critical_temperature).all().detach()) or bool(
        (critical_temperature <= 0.0).any().detach()
    ):
        raise ValueError("critical temperatures must be finite and positive")
    if not bool(torch.isfinite(volume).all().detach()) or bool((volume <= 0.0).any().detach()):
        raise ValueError("critical volumes must be finite and positive")
    x = normalize_composition(composition)
    if x.shape[-1] != critical_temperature.numel():
        raise ValueError("composition and critical-property component counts differ")
    weights = x * volume
    return factor * torch.sum(weights * critical_temperature, dim=-1) / weights.sum(dim=-1)

negative_flash_residual

negative_flash_residual(
    model: object, state: ChemicalState
) -> Tensor

Evaluate Perschke's dimensionless negative-flash residual.

Parameters:

Name	Type	Description	Default
model	object	Component parameter source exposing critical temperature in K, critical pressure in Pa, and acentric factor tensors. These are the quantities used by the Wilson equilibrium-ratio estimate.	required
state	ChemicalState	Temperature in K, pressure in Pa, and normalized composition. Matching leading state batches are supported; components occupy the final composition axis.	required

Returns:

Type	Description
Tensor	Perschke's dimensionless residual .. math:: G(0.5)=\sum_i \frac{x_i(K_i^{\mathrm{Wilson}}-1)} {1+0.5(K_i^{\mathrm{Wilson}}-1)} with the leading state shape. Positive values are vapor-like and negative values are liquid-like under this diagnostic. The tensor retains its PyTorch graph to differentiable state inputs and component parameters used by the Wilson estimate.

Raises:

Type	Description
TypeError	If any required critical-temperature, critical-pressure, or acentric-factor tensor is unavailable.

Notes

This function evaluates a homogeneous-state diagnostic only. It does not perform a negative flash, select an EoS root, establish phase stability, or determine the number of equilibrium phases.

Source code in src/torch_flash/properties/phase_identification.py

def negative_flash_residual(model: object, state: ChemicalState) -> Tensor:
    r"""Evaluate Perschke's dimensionless negative-flash residual.

    Parameters
    ----------
    model
        Component parameter source exposing critical temperature in K,
        critical pressure in Pa, and acentric factor tensors. These are the
        quantities used by the Wilson equilibrium-ratio estimate.
    state
        Temperature in K, pressure in Pa, and normalized composition.
        Matching leading state batches are supported; components occupy the
        final composition axis.

    Returns
    -------
    Tensor
        Perschke's dimensionless residual

        .. math::

            G(0.5)=\sum_i
            \frac{x_i(K_i^{\mathrm{Wilson}}-1)}
                 {1+0.5(K_i^{\mathrm{Wilson}}-1)}

        with the leading state shape. Positive values are vapor-like and
        negative values are liquid-like under this diagnostic. The tensor
        retains its PyTorch graph to differentiable state inputs and component
        parameters used by the Wilson estimate.

    Raises
    ------
    TypeError
        If any required critical-temperature, critical-pressure, or
        acentric-factor tensor is unavailable.

    Notes
    -----
    This function evaluates a homogeneous-state diagnostic only. It does not
    perform a negative flash, select an EoS root, establish phase stability,
    or determine the number of equilibrium phases.
    """
    required = ("critical_temperature", "critical_pressure", "acentric_factor")
    if not all(isinstance(getattr(model, name, None), Tensor) for name in required):
        raise TypeError("negative-flash method requires critical constants and acentric factors")
    components = cast(ComponentSet, model)
    k_values = wilson_k_values(components, state.temperature, state.pressure)
    terms = state.composition * (k_values - 1.0) / (1.0 + 0.5 * (k_values - 1.0))
    return terms.sum(dim=-1)

phase_identification_parameter

phase_identification_parameter(
    model: object,
    state: ChemicalState,
    phase: PhaseKind = "stable",
    *,
    denominator_relative_tolerance: float | None = None,
) -> Tensor

Evaluate the Venkatarathnam-Oellrich phase-identification parameter.

The dimensionless parameter is

.. math::

\Pi = V\left[
\frac{(\partial^2 P/\partial V\partial T)_x}

     {(\partial P/\partial T)_{V,x}}

\frac{(\partial^2 P/\partial V^2)_{T,x}}
     {(\partial P/\partial V)_{T,x}}
\right].

Venkatarathnam and Oellrich classify a mechanically stable homogeneous state as liquid or liquid-like when Pi > 1 and vapor-like when Pi <= 1. The method labels a supplied homogeneous state; it does not select an EoS root, perform a stability test, or determine the equilibrium phase count.

Parameters:

Name	Type	Description	Default
model	object	Homogeneous-state model implementing pressure(T, V, x) and molar_volume(T, P, x, phase) with PyTorch tensor operations. pressure must return pressure in Pa for physical molar volume in m3/mol.	required
state	ChemicalState	Temperature in K, pressure in Pa, and normalized composition. Matching leading state batches are supported; components occupy the final composition axis.	required
phase	PhaseKind	Homogeneous EoS root used only to obtain the evaluation volume. The default "stable" selects the admissible root of lowest residual Gibbs energy; it does not perform a flash.	'stable'
denominator_relative_tolerance	float | None	Optional non-negative singularity threshold. |P_T| is compared with this value times |P/T| and |P_V| with this value times |P/V|. When omitted, the threshold is the square root of machine epsilon for the state's dtype. Set zero only for controlled numerical experiments that require exact-zero guards.	None

Returns:

Type	Description
Tensor	Dimensionless Pi with the leading state shape. A scalar state returns a scalar tensor. The result remains connected through the first and second pressure derivatives, the selected molar volume, the state tensors, and any trainable model parameters.

Raises:

Type	Description
TypeError	If the model does not expose differentiable pressure and molar_volume methods.
ValueError	If the state shapes are inconsistent, any molar volume or pressure derivative is nonfinite, any volume is nonpositive, P_V >= 0 indicates a mechanically unstable root, a scaled denominator is singular, or the relative tolerance is invalid. For a batch, invalid state checks reject the complete batch so callers can isolate or subdivide it explicitly.

Notes

Nested torch.func.jvp calls evaluate the first and second directional derivatives of the batched explicit pressure. Because each output state depends only on the matching input state, all-one temperature or volume tangents recover the independent Jacobian diagonals without constructing a dense cross-state Jacobian. No finite-difference step or cubic-specific derivative formula is used.

The paper describes an additional high-temperature inversion correction for isolated superheated states. It also states that this correction is unnecessary when identifying saturated split phases after a converged flash. This function returns the local parameter itself and does not perform a remote temperature root search.

Source code in src/torch_flash/properties/phase_identification.py

def phase_identification_parameter(
    model: object,
    state: ChemicalState,
    phase: PhaseKind = "stable",
    *,
    denominator_relative_tolerance: float | None = None,
) -> Tensor:
    r"""Evaluate the Venkatarathnam-Oellrich phase-identification parameter.

    The dimensionless parameter is

    .. math::

        \Pi = V\left[
        \frac{(\partial^2 P/\partial V\partial T)_x}
             {(\partial P/\partial T)_{V,x}}
        -
        \frac{(\partial^2 P/\partial V^2)_{T,x}}
             {(\partial P/\partial V)_{T,x}}
        \right].

    Venkatarathnam and Oellrich classify a mechanically stable homogeneous
    state as liquid or liquid-like when ``Pi > 1`` and vapor-like when
    ``Pi <= 1``. The method labels a supplied homogeneous state; it does not
    select an EoS root, perform a stability test, or determine the equilibrium
    phase count.

    Parameters
    ----------
    model
        Homogeneous-state model implementing ``pressure(T, V, x)`` and
        ``molar_volume(T, P, x, phase)`` with PyTorch tensor operations.
        ``pressure`` must return pressure in Pa for physical molar volume in
        m3/mol.
    state
        Temperature in K, pressure in Pa, and normalized composition. Matching
        leading state batches are supported; components occupy the final
        composition axis.
    phase
        Homogeneous EoS root used only to obtain the evaluation volume. The
        default ``"stable"`` selects the admissible root of lowest residual
        Gibbs energy; it does not perform a flash.
    denominator_relative_tolerance
        Optional non-negative singularity threshold. ``|P_T|`` is compared
        with this value times ``|P/T|`` and ``|P_V|`` with this value times
        ``|P/V|``. When omitted, the threshold is the square root of machine
        epsilon for the state's dtype. Set zero only for controlled numerical
        experiments that require exact-zero guards.

    Returns
    -------
    Tensor
        Dimensionless ``Pi`` with the leading state shape. A scalar state
        returns a scalar tensor. The result remains connected through the
        first and second pressure derivatives, the selected molar volume, the
        state tensors, and any trainable model parameters.

    Raises
    ------
    TypeError
        If the model does not expose differentiable ``pressure`` and
        ``molar_volume`` methods.
    ValueError
        If the state shapes are inconsistent, any molar volume or pressure
        derivative is nonfinite, any volume is nonpositive, ``P_V >= 0``
        indicates a mechanically unstable root, a scaled denominator is
        singular, or the relative tolerance is invalid. For a batch, invalid
        state checks reject the complete batch so callers can isolate or
        subdivide it explicitly.

    Notes
    -----
    Nested ``torch.func.jvp`` calls evaluate the first and second directional
    derivatives of the batched explicit pressure. Because each output state
    depends only on the matching input state, all-one temperature or volume
    tangents recover the independent Jacobian diagonals without constructing
    a dense cross-state Jacobian. No finite-difference step or cubic-specific
    derivative formula is used.

    The paper describes an additional high-temperature inversion correction
    for isolated superheated states. It also states that this correction is
    unnecessary when identifying saturated split phases after a converged
    flash. This function returns the local parameter itself and does not
    perform a remote temperature root search.
    """
    _validate_pip_denominator_tolerance(denominator_relative_tolerance)
    state_shape = state.temperature.shape
    if state.pressure.shape != state_shape or state.composition.shape[:-1] != state_shape:
        raise ValueError("phase-identification parameter state shapes are inconsistent")
    pressure_function: Any = getattr(model, "pressure", None)
    volume_function: Any = getattr(model, "molar_volume", None)
    if not callable(pressure_function) or not callable(volume_function):
        raise TypeError("phase-identification parameter requires pressure and molar_volume methods")
    volume = cast(
        Tensor,
        volume_function(
            state.temperature,
            state.pressure,
            state.composition,
            phase,
        ),
    )
    if (
        volume.shape != state_shape
        or not bool(torch.isfinite(volume).all().detach())
        or bool((volume <= 0.0).any().detach())
    ):
        raise ValueError(
            "phase-identification molar volume must match the state shape and be finite positive"
        )
    composition = state.composition

    def pressure_at_temperature_volume(
        temperature: Tensor,
        molar_volume: Tensor,
    ) -> Tensor:
        return cast(
            Tensor,
            pressure_function(
                temperature,
                molar_volume,
                composition,
            ),
        )

    temperature_direction = torch.ones_like(state.temperature)
    volume_direction = torch.ones_like(volume)
    zero_temperature = torch.zeros_like(state.temperature)
    zero_volume = torch.zeros_like(volume)

    def pressure_volume_at_temperature_volume(
        temperature: Tensor,
        molar_volume: Tensor,
    ) -> Tensor:
        return cast(
            Tensor,
            torch.func.jvp(
                pressure_at_temperature_volume,
                (temperature, molar_volume),
                (zero_temperature, volume_direction),
            )[1],
        )

    pressure_temperature = torch.func.jvp(
        pressure_at_temperature_volume,
        (state.temperature, volume),
        (temperature_direction, zero_volume),
    )[1]
    pressure_volume = pressure_volume_at_temperature_volume(
        state.temperature,
        volume,
    )
    pressure_volume_temperature = torch.func.jvp(
        pressure_volume_at_temperature_volume,
        (state.temperature, volume),
        (temperature_direction, zero_volume),
    )[1]
    pressure_volume_second = torch.func.jvp(
        pressure_volume_at_temperature_volume,
        (state.temperature, volume),
        (zero_temperature, volume_direction),
    )[1]
    derivatives = torch.stack(
        (
            pressure_temperature,
            pressure_volume,
            pressure_volume_temperature,
            pressure_volume_second,
        )
    )
    if not bool(torch.isfinite(derivatives).all().detach()):
        raise ValueError("phase-identification pressure derivatives are nonfinite")
    if bool((pressure_volume >= 0.0).any().detach()):
        raise _PhaseIdentificationParameterSingularityError(
            "phase-identification parameter requires a mechanically stable P_V < 0 root"
        )

    relative_tolerance = (
        math.sqrt(torch.finfo(volume.dtype).eps)
        if denominator_relative_tolerance is None
        else denominator_relative_tolerance
    )
    pressure_temperature_scale = torch.abs(state.pressure / state.temperature)
    pressure_volume_scale = torch.abs(state.pressure / volume)
    if bool(
        (torch.abs(pressure_temperature) <= relative_tolerance * pressure_temperature_scale)
        .any()
        .detach()
    ) or bool(
        (torch.abs(pressure_volume) <= relative_tolerance * pressure_volume_scale).any().detach()
    ):
        raise _PhaseIdentificationParameterSingularityError(
            "phase-identification parameter has a scaled singular P_T or P_V denominator"
        )
    return cast(
        Tensor,
        volume
        * (
            pressure_volume_temperature / pressure_temperature
            - pressure_volume_second / pressure_volume
        ),
    )

phase_response_derivatives

phase_response_derivatives(
    model: object,
    state: ChemicalState,
    phase: PhaseKind = "stable",
) -> PhaseResponseDerivatives

Evaluate batched homogeneous response functions with PyTorch autodiff.

Parameters:

Name	Type	Description	Default
model	object	Homogeneous-state model implementing differentiable pressure(T, V, x) and molar_volume(T, P, x, phase) methods. Pressure is in Pa and physical molar volume is in m3/mol.	required
state	ChemicalState	Temperature in K, pressure in Pa, and normalized composition. Temperature and pressure share the leading state shape; components occupy the final composition axis. Scalar and leading-batch states are both supported.	required
phase	PhaseKind	Homogeneous EoS root request used to obtain the evaluation volume. Root selection remains separate from the returned physical response.	'stable'

Returns:

Type	Description
PhaseResponseDerivatives	Molar volume, isothermal compressibility, thermal-expansion coefficient, and their constant-pressure temperature derivatives. Every field has the leading state shape and remains connected to state tensors and trainable model parameters.

Raises:

Type	Description
TypeError	If the model lacks the required pressure or molar-volume operation.
ValueError	If state or volume shapes are inconsistent, any volume or derivative is nonfinite, any volume is nonpositive, or a selected root has P_V >= 0 and is therefore mechanically unstable.

Notes

Nested :func:torch.func.jvp calls evaluate the independent first and second pressure-response derivatives. For a leading batch, each output state depends only on the matching input state, so all-one temperature or volume tangents recover the Jacobian diagonals without constructing a dense cross-state Jacobian. The implementation is mathematically equivalent to differentiating each scalar response with :func:torch.func.grad and :func:torch.func.jacrev, and is tested against that reverse-mode scalar oracle.

Source code in src/torch_flash/properties/phase_identification.py

def phase_response_derivatives(
    model: object,
    state: ChemicalState,
    phase: PhaseKind = "stable",
) -> PhaseResponseDerivatives:
    r"""Evaluate batched homogeneous response functions with PyTorch autodiff.

    Parameters
    ----------
    model
        Homogeneous-state model implementing differentiable
        ``pressure(T, V, x)`` and ``molar_volume(T, P, x, phase)`` methods.
        Pressure is in Pa and physical molar volume is in m3/mol.
    state
        Temperature in K, pressure in Pa, and normalized composition.
        Temperature and pressure share the leading state shape; components
        occupy the final composition axis. Scalar and leading-batch states are
        both supported.
    phase
        Homogeneous EoS root request used to obtain the evaluation volume.
        Root selection remains separate from the returned physical response.

    Returns
    -------
    PhaseResponseDerivatives
        Molar volume, isothermal compressibility, thermal-expansion
        coefficient, and their constant-pressure temperature derivatives.
        Every field has the leading state shape and remains connected to state
        tensors and trainable model parameters.

    Raises
    ------
    TypeError
        If the model lacks the required pressure or molar-volume operation.
    ValueError
        If state or volume shapes are inconsistent, any volume or derivative
        is nonfinite, any volume is nonpositive, or a selected root has
        ``P_V >= 0`` and is therefore mechanically unstable.

    Notes
    -----
    Nested :func:`torch.func.jvp` calls evaluate the independent first and
    second pressure-response derivatives. For a leading batch, each output
    state depends only on the matching input state, so all-one temperature or
    volume tangents recover the Jacobian diagonals without constructing a
    dense cross-state Jacobian. The implementation is mathematically
    equivalent to differentiating each scalar response with
    :func:`torch.func.grad` and :func:`torch.func.jacrev`, and is tested
    against that reverse-mode scalar oracle.
    """
    return _phase_response_details(model, state, phase).values

volume_to_covolume_ratio

volume_to_covolume_ratio(
    model: object,
    state: ChemicalState,
    phase: PhaseKind = "stable",
) -> Tensor

Return the differentiable cubic-family V/b phase diagnostic.

Parameters:

Name	Type	Description	Default
model	object	Model implementing select_z and cubic mixture_parameters.	required
state	ChemicalState	Temperature in K, pressure in Pa, and mole fractions.	required
phase	PhaseKind	Homogeneous root-selection request.	'stable'

Returns:

Type	Description
Tensor	Positive dimensionless EOS volume-to-covolume ratio.

Raises:

Type	Description
TypeError	If the model lacks root selection or a mixture covolume.
ValueError	If the evaluated covolume or ratio is nonpositive/nonfinite.

Notes

Leading batch dimensions are supported. Cubic volume translations are excluded because Pedersen's criterion uses the volume entering the cubic repulsive term.

Source code in src/torch_flash/properties/phase_identification.py

def volume_to_covolume_ratio(
    model: object,
    state: ChemicalState,
    phase: PhaseKind = "stable",
) -> Tensor:
    """Return the differentiable cubic-family ``V/b`` phase diagnostic.

    Parameters
    ----------
    model
        Model implementing ``select_z`` and cubic ``mixture_parameters``.
    state
        Temperature in K, pressure in Pa, and mole fractions.
    phase
        Homogeneous root-selection request.

    Returns
    -------
    Tensor
        Positive dimensionless EOS volume-to-covolume ratio.

    Raises
    ------
    TypeError
        If the model lacks root selection or a mixture covolume.
    ValueError
        If the evaluated covolume or ratio is nonpositive/nonfinite.

    Notes
    -----
    Leading batch dimensions are supported. Cubic volume translations are
    excluded because Pedersen's criterion uses the volume entering the cubic
    repulsive term.
    """
    select_z: Any = getattr(model, "select_z", None)
    if not callable(select_z):
        raise TypeError("phase-identification model must implement select_z")
    covolume = _covolume(model, state)
    if covolume is None:
        raise TypeError("phase-identification model does not expose a mixture covolume")
    compressibility_factor = cast(
        Tensor,
        select_z(
            state.temperature,
            state.pressure,
            state.composition,
            phase,
        ),
    )
    equation_of_state_volume = compressibility_factor * R * state.temperature / state.pressure
    ratio = equation_of_state_volume / covolume
    if not bool(torch.isfinite(ratio).all().detach()) or bool((ratio <= 0.0).any().detach()):
        raise ValueError("volume-to-covolume ratio must be finite and positive")
    return ratio

fugacities_tv

fugacities_tv(
    model: HelmholtzStateModel,
    temperature: Tensor,
    volume: Tensor,
    moles: Tensor,
) -> Tensor

Return component fugacities from an explicit T-V-n state.

Parameters:

Name	Type	Description	Default
model	HelmholtzStateModel	Model exposing an extensive residual Helmholtz energy.	required
temperature	Tensor	Scalar positive temperature in K.	required
volume	Tensor	Scalar positive total physical volume in m3.	required
moles	Tensor	Strictly positive component amounts in mol with shape (ncomponents,).	required

Returns:

Type	Description
Tensor	Component fugacities in Pa, obtained by exponentiating :func:log_fugacities_tv.

Raises:

Type	Description
ValueError	If the state is not scalar/vector shaped or contains nonpositive or nonfinite inputs.

Source code in src/torch_flash/properties/state.py

def fugacities_tv(
    model: HelmholtzStateModel,
    temperature: Tensor,
    volume: Tensor,
    moles: Tensor,
) -> Tensor:
    """Return component fugacities from an explicit ``T-V-n`` state.

    Parameters
    ----------
    model
        Model exposing an extensive residual Helmholtz energy.
    temperature
        Scalar positive temperature in K.
    volume
        Scalar positive total physical volume in m3.
    moles
        Strictly positive component amounts in mol with shape
        ``(ncomponents,)``.

    Returns
    -------
    Tensor
        Component fugacities in Pa, obtained by exponentiating
        :func:`log_fugacities_tv`.

    Raises
    ------
    ValueError
        If the state is not scalar/vector shaped or contains nonpositive or
        nonfinite inputs.
    """
    return STANDARD_PRESSURE * torch.exp(log_fugacities_tv(model, temperature, volume, moles))

log_fugacities_tv

log_fugacities_tv(
    model: HelmholtzStateModel,
    temperature: Tensor,
    volume: Tensor,
    moles: Tensor,
) -> Tensor

Return ln(f_i/p_standard) from an explicit T-V-n state.

Parameters:

Name	Type	Description	Default
model	HelmholtzStateModel	Model exposing an extensive residual Helmholtz energy.	required
temperature	Tensor	Scalar positive temperature in K.	required
volume	Tensor	Scalar positive total physical volume in m3.	required
moles	Tensor	Strictly positive component amounts in mol.	required

Returns:

Type	Description
Tensor	Dimensionless component ln(f_i / 1 bar) values.

Raises:

Type	Description
ValueError	If state shapes or positivity/finite-value requirements are violated.

Notes

volume is the total physical volume in m3 and moles are component amounts in mol. The identity

ln(f_i/p_standard) = ln(n_i*R*T/(V*p_standard)) + d(A^R/RT)/dn_i

provides an independent Helmholtz-route check of the usual TP fugacity calculation. The model's residual Helmholtz energy must be extensive and referenced to an ideal gas at the same physical volume.

Source code in src/torch_flash/properties/state.py

def log_fugacities_tv(
    model: HelmholtzStateModel,
    temperature: Tensor,
    volume: Tensor,
    moles: Tensor,
) -> Tensor:
    """Return ``ln(f_i/p_standard)`` from an explicit ``T-V-n`` state.

    Parameters
    ----------
    model
        Model exposing an extensive residual Helmholtz energy.
    temperature
        Scalar positive temperature in K.
    volume
        Scalar positive total physical volume in m3.
    moles
        Strictly positive component amounts in mol.

    Returns
    -------
    Tensor
        Dimensionless component ``ln(f_i / 1 bar)`` values.

    Raises
    ------
    ValueError
        If state shapes or positivity/finite-value requirements are violated.

    Notes
    -----
    ``volume`` is the total physical volume in m3 and ``moles`` are component
    amounts in mol. The identity

    ``ln(f_i/p_standard) = ln(n_i*R*T/(V*p_standard)) + d(A^R/RT)/dn_i``

    provides an independent Helmholtz-route check of the usual TP fugacity
    calculation. The model's residual Helmholtz energy must be extensive and
    referenced to an ideal gas at the same physical volume.
    """
    if temperature.ndim != 0 or volume.ndim != 0:
        raise ValueError("log_fugacities_tv currently accepts scalar temperature and volume")
    if moles.ndim != 1:
        raise ValueError("log_fugacities_tv currently accepts one mole-number vector")
    if bool((temperature <= 0.0) | (volume <= 0.0)):
        raise ValueError("temperature and volume must be positive")
    if bool((~torch.isfinite(moles)).any() | (moles <= 0.0).any()):
        raise ValueError("moles must be finite and strictly positive")

    residual_mu_rt: Tensor = torch.func.jacrev(
        lambda current_moles: model.residual_helmholtz_rt(
            temperature,
            volume,
            current_moles,
        )
    )(moles)
    ideal_log_fugacity = torch.log(moles * R * temperature / (volume * STANDARD_PRESSURE))
    return ideal_log_fugacity + residual_mu_rt

phase_properties

phase_properties(
    model: StateModel,
    state: ChemicalState,
    phase: PhaseKind = "stable",
    *,
    caloric: bool = True,
    standard_state: StandardState | None = None,
) -> PhaseProperties

Evaluate a homogeneous state, without any equilibrium calculation.

Parameters:

Name	Type	Description	Default
model	StateModel	Homogeneous-state model.	required
state	ChemicalState	Temperature in K, pressure in Pa, and mole fractions.	required
phase	PhaseKind	Liquid, vapor, or stable-root request.	'stable'
caloric	bool	Evaluate scalar-state residual enthalpy and entropy when possible.	True
standard_state	StandardState | None	Optional component standard chemical potentials.	None

Returns:

Type	Description
PhaseProperties	Homogeneous properties, reference-aware chemical potentials, residual free energies, and physical phase-identification diagnostics.

Notes

Fugacities are returned in Pa and logarithmic fugacities are ln(f_i / p_standard) with p_standard = 1 bar. Chemical potentials and total free energies use a zero ideal-gas standard chemical potential at that pressure unless standard_state is supplied. The corresponding dimensionless chemical potential is mu_i / (R*T); it is the thermodynamically meaningful alternative to ln(mu_i), which is not generally defined for a dimensional, reference-dependent quantity that may be negative.

Absolute caloric reference terms are deliberately not invented; returned enthalpy and entropy are residual quantities.

The reduced molar free energies are reduced_gibbs_energy = g/(R*T) and reduced_helmholtz_energy = a/(R*T), with a = g - P*v. The reference-independent departures follow g^R/(R*T) = sum(x_i*ln(phi_i)) and a^R/(R*T) = g^R/(R*T) - Z + 1 + ln(Z), where Z = P*v/(R*T).

At an exactly zero mole fraction, the logarithmic component properties use the smallest positive value of the tensor dtype as a finite trace limit. state_derivatives instead requires a strictly positive, interior composition because logarithmic composition derivatives are singular on the simplex boundary.

Source code in src/torch_flash/properties/state.py

def phase_properties(
    model: StateModel,
    state: ChemicalState,
    phase: PhaseKind = "stable",
    *,
    caloric: bool = True,
    standard_state: StandardState | None = None,
) -> PhaseProperties:
    """Evaluate a homogeneous state, without any equilibrium calculation.

    Parameters
    ----------
    model
        Homogeneous-state model.
    state
        Temperature in K, pressure in Pa, and mole fractions.
    phase
        Liquid, vapor, or stable-root request.
    caloric
        Evaluate scalar-state residual enthalpy and entropy when possible.
    standard_state
        Optional component standard chemical potentials.

    Returns
    -------
    PhaseProperties
        Homogeneous properties, reference-aware chemical potentials, residual
        free energies, and physical phase-identification diagnostics.

    Notes
    -----
    Fugacities are returned in Pa and logarithmic fugacities are
    ``ln(f_i / p_standard)`` with ``p_standard = 1 bar``. Chemical potentials
    and total free energies use a zero ideal-gas standard chemical potential
    at that pressure unless ``standard_state`` is supplied. The corresponding
    dimensionless chemical potential is ``mu_i / (R*T)``; it is the
    thermodynamically meaningful alternative to ``ln(mu_i)``, which is not
    generally defined for a dimensional, reference-dependent quantity that
    may be negative.

    Absolute caloric reference terms are deliberately not invented; returned
    enthalpy and entropy are residual quantities.

    The reduced molar free energies are
    ``reduced_gibbs_energy = g/(R*T)`` and
    ``reduced_helmholtz_energy = a/(R*T)``, with ``a = g - P*v``. The
    reference-independent departures follow
    ``g^R/(R*T) = sum(x_i*ln(phi_i))`` and
    ``a^R/(R*T) = g^R/(R*T) - Z + 1 + ln(Z)``, where
    ``Z = P*v/(R*T)``.

    At an exactly zero mole fraction, the logarithmic component properties
    use the smallest positive value of the tensor dtype as a finite trace
    limit. ``state_derivatives`` instead requires a strictly positive,
    interior composition because logarithmic composition derivatives are
    singular on the simplex boundary.
    """
    temperature = state.temperature
    pressure = state.pressure
    composition = state.composition
    z_factor = model.select_z(temperature, pressure, composition, phase)
    volume = model.molar_volume(temperature, pressure, composition, phase)
    log_phi = model.log_fugacity_coefficients(temperature, pressure, composition, phase)
    log_fugacity = _log_fugacities_from_coefficients(log_phi, pressure, composition)
    fugacity = STANDARD_PRESSURE * torch.exp(log_fugacity)
    chemical_potentials = _chemical_potentials_from_log_fugacities(
        log_fugacity,
        temperature,
        standard_state,
    )
    reduced_chemical_potentials = chemical_potentials / (R * temperature[..., None])
    gibbs = torch.sum(composition * chemical_potentials, dim=-1)
    pressure_volume_over_rt = pressure * volume / (R * temperature)
    helmholtz = gibbs - pressure * volume
    reduced_gibbs = gibbs / (R * temperature)
    reduced_helmholtz = helmholtz / (R * temperature)
    reduced_residual_gibbs = torch.sum(composition * log_phi, dim=-1)
    reduced_residual_helmholtz = (
        reduced_residual_gibbs - pressure_volume_over_rt + 1.0 + torch.log(pressure_volume_over_rt)
    )

    residual_enthalpy: Tensor | None = None
    residual_entropy: Tensor | None = None
    if caloric and temperature.ndim == 0 and composition.ndim == 1:

        def residual_gibbs_over_t(current_temperature: Tensor) -> Tensor:
            current_log_phi = model.log_fugacity_coefficients(
                current_temperature, pressure, composition, phase
            )
            return R * torch.sum(composition * current_log_phi)

        derivative = torch.func.grad(residual_gibbs_over_t)(temperature)
        residual_enthalpy = -temperature.square() * derivative
        residual_gibbs = R * temperature * torch.sum(composition * log_phi)
        residual_entropy = (residual_enthalpy - residual_gibbs) / temperature

    properties = PhaseProperties(
        kind=phase,
        composition=composition,
        compressibility_factor=z_factor,
        molar_volume=volume,
        log_fugacity_coefficients=log_phi,
        fugacities=fugacity,
        log_fugacities=log_fugacity,
        chemical_potentials=chemical_potentials,
        reduced_chemical_potentials=reduced_chemical_potentials,
        molar_gibbs_energy=gibbs,
        molar_helmholtz_energy=helmholtz,
        reduced_gibbs_energy=reduced_gibbs,
        reduced_helmholtz_energy=reduced_helmholtz,
        reduced_residual_gibbs_energy=reduced_residual_gibbs,
        reduced_residual_helmholtz_energy=reduced_residual_helmholtz,
        residual_enthalpy=residual_enthalpy,
        residual_entropy=residual_entropy,
    )
    return replace(
        properties,
        phase_identification=identify_phase_from_properties(model, state, properties),
    )

state_derivatives

state_derivatives(
    model: StateModel,
    state: ChemicalState,
    phase: PhaseKind = "stable",
    *,
    standard_state: StandardState | None = None,
) -> ThermodynamicDerivatives

Differentiate homogeneous TP properties with PyTorch autodiff.

Parameters:

Name	Type	Description	Default
model	StateModel	Homogeneous-state model to differentiate.	required
state	ChemicalState	One scalar temperature-pressure state. Temperature is in K, pressure in Pa, and composition is a strictly positive mole-fraction vector.	required
phase	PhaseKind	Root-selection request held fixed throughout every derivative.	'stable'
standard_state	StandardState | None	Optional component standard chemical potentials. When omitted, chemical potentials contain only the ideal-mixing/fugacity departure relative to 1 bar.	None

Returns:

Type	Description
ThermodynamicDerivatives	Jacobians with respect to unconstrained logits, n-1 independent mole fractions, unit-basis mole numbers, temperature, and pressure. Field units and coordinate conventions are documented on the result class.

Raises:

Type	Description
ValueError	If temperature or pressure is not scalar, composition is not a vector, or a mole fraction is zero or negative.

Notes

The phase root is not differentiated through a root-switching event. Derivatives near coalescing roots or the composition-simplex boundary can be ill-conditioned even when the underlying property value is finite.

Source code in src/torch_flash/properties/state.py

def state_derivatives(
    model: StateModel,
    state: ChemicalState,
    phase: PhaseKind = "stable",
    *,
    standard_state: StandardState | None = None,
) -> ThermodynamicDerivatives:
    """Differentiate homogeneous TP properties with PyTorch autodiff.

    Parameters
    ----------
    model
        Homogeneous-state model to differentiate.
    state
        One scalar temperature-pressure state. Temperature is in K, pressure
        in Pa, and composition is a strictly positive mole-fraction vector.
    phase
        Root-selection request held fixed throughout every derivative.
    standard_state
        Optional component standard chemical potentials. When omitted,
        chemical potentials contain only the ideal-mixing/fugacity departure
        relative to 1 bar.

    Returns
    -------
    ThermodynamicDerivatives
        Jacobians with respect to unconstrained logits, ``n-1`` independent
        mole fractions, unit-basis mole numbers, temperature, and pressure.
        Field units and coordinate conventions are documented on the result
        class.

    Raises
    ------
    ValueError
        If temperature or pressure is not scalar, composition is not a vector,
        or a mole fraction is zero or negative.

    Notes
    -----
    The phase root is not differentiated through a root-switching event.
    Derivatives near coalescing roots or the composition-simplex boundary can
    be ill-conditioned even when the underlying property value is finite.
    """
    if state.temperature.ndim != 0 or state.pressure.ndim != 0:
        raise ValueError("state_derivatives currently accepts one scalar T-P state")
    if state.composition.ndim != 1:
        raise ValueError("state_derivatives currently accepts one composition vector")
    if bool((state.composition <= 0.0).any()):
        raise ValueError("state_derivatives requires strictly positive mole fractions")

    def composition_from_logits(logit_values: Tensor) -> Tensor:
        return torch.softmax(logit_values, dim=-1)

    def composition_from_independent(independent: Tensor) -> Tensor:
        return torch.cat((independent, (1.0 - independent.sum()).reshape(1)))

    def composition_from_moles(moles: Tensor) -> Tensor:
        return moles / moles.sum()

    def log_phi_at(composition: Tensor, temperature: Tensor, pressure: Tensor) -> Tensor:
        return model.log_fugacity_coefficients(
            temperature,
            pressure,
            composition,
            phase,
        )

    def log_fugacity_at(composition: Tensor, temperature: Tensor, pressure: Tensor) -> Tensor:
        return _log_fugacities(
            model,
            temperature,
            pressure,
            composition,
            phase,
        )

    def mu_at(composition: Tensor, temperature: Tensor, pressure: Tensor) -> Tensor:
        return _chemical_potentials(
            model,
            temperature,
            pressure,
            composition,
            phase,
            standard_state,
        )

    def volume_at(composition: Tensor, temperature: Tensor, pressure: Tensor) -> Tensor:
        return model.molar_volume(temperature, pressure, composition, phase)

    logits = torch.log(torch.clamp_min(state.composition, 1.0e-300))
    independent = state.composition[:-1]
    unit_moles = state.composition

    PropertyFunction = Callable[[Tensor, Tensor, Tensor], Tensor]

    def apply_logits(function: PropertyFunction, current: Tensor) -> Tensor:
        return function(
            composition_from_logits(current),
            state.temperature,
            state.pressure,
        )

    def apply_independent(function: PropertyFunction, current: Tensor) -> Tensor:
        return function(
            composition_from_independent(current),
            state.temperature,
            state.pressure,
        )

    def apply_moles(function: PropertyFunction, current: Tensor) -> Tensor:
        return function(
            composition_from_moles(current),
            state.temperature,
            state.pressure,
        )

    dlog_phi_dlogits = torch.func.jacrev(lambda value: apply_logits(log_phi_at, value))(logits)
    dlog_phi_dindependent = torch.func.jacrev(lambda value: apply_independent(log_phi_at, value))(
        independent
    )
    dlog_phi_dmoles = torch.func.jacrev(lambda value: apply_moles(log_phi_at, value))(unit_moles)
    dlog_phi_dt = torch.func.jacrev(lambda t: log_phi_at(state.composition, t, state.pressure))(
        state.temperature
    )
    dlog_phi_dp = torch.func.jacrev(lambda p: log_phi_at(state.composition, state.temperature, p))(
        state.pressure
    )

    dlog_fugacity_dlogits = torch.func.jacrev(lambda value: apply_logits(log_fugacity_at, value))(
        logits
    )
    dlog_fugacity_dindependent = torch.func.jacrev(
        lambda value: apply_independent(log_fugacity_at, value)
    )(independent)
    dlog_fugacity_dmoles = torch.func.jacrev(lambda value: apply_moles(log_fugacity_at, value))(
        unit_moles
    )
    dlog_fugacity_dt = torch.func.jacrev(
        lambda t: log_fugacity_at(state.composition, t, state.pressure)
    )(state.temperature)
    dlog_fugacity_dp = torch.func.jacrev(
        lambda p: log_fugacity_at(state.composition, state.temperature, p)
    )(state.pressure)

    dmu_dlogits = torch.func.jacrev(lambda value: apply_logits(mu_at, value))(logits)
    dmu_dindependent = torch.func.jacrev(lambda value: apply_independent(mu_at, value))(independent)
    dmu_dmoles = torch.func.jacrev(lambda value: apply_moles(mu_at, value))(unit_moles)
    dmu_dt = torch.func.jacrev(lambda t: mu_at(state.composition, t, state.pressure))(
        state.temperature
    )
    dmu_dp = torch.func.jacrev(lambda p: mu_at(state.composition, state.temperature, p))(
        state.pressure
    )

    dvolume_dlogits = torch.func.jacrev(lambda value: apply_logits(volume_at, value))(logits)
    dvolume_dindependent = torch.func.jacrev(lambda value: apply_independent(volume_at, value))(
        independent
    )
    dvolume_dmoles = torch.func.jacrev(lambda value: apply_moles(volume_at, value))(unit_moles)
    dvolume_dt = torch.func.jacrev(lambda t: volume_at(state.composition, t, state.pressure))(
        state.temperature
    )
    dvolume_dp = torch.func.jacrev(lambda p: volume_at(state.composition, state.temperature, p))(
        state.pressure
    )

    def gibbs(t: Tensor, p: Tensor) -> Tensor:
        mu = mu_at(state.composition, t, p)
        return torch.sum(state.composition * mu)

    dg_dt = torch.func.grad(lambda t: gibbs(t, state.pressure))(state.temperature)
    dg_dp = torch.func.grad(lambda p: gibbs(state.temperature, p))(state.pressure)

    log_fugacity = _log_fugacities(
        model,
        state.temperature,
        state.pressure,
        state.composition,
        phase,
    )
    fugacity = STANDARD_PRESSURE * torch.exp(log_fugacity)
    log_phi = log_phi_at(state.composition, state.temperature, state.pressure)
    phi = torch.exp(log_phi)
    dphi_dlogits = phi[:, None] * dlog_phi_dlogits
    dphi_dindependent = phi[:, None] * dlog_phi_dindependent
    dphi_dmoles = phi[:, None] * dlog_phi_dmoles
    dphi_dt = phi * dlog_phi_dt
    dphi_dp = phi * dlog_phi_dp
    dfugacity_dlogits = fugacity[:, None] * dlog_fugacity_dlogits
    dfugacity_dindependent = fugacity[:, None] * dlog_fugacity_dindependent
    dfugacity_dmoles = fugacity[:, None] * dlog_fugacity_dmoles
    dfugacity_dt = fugacity * dlog_fugacity_dt
    dfugacity_dp = fugacity * dlog_fugacity_dp

    chemical_potential = _chemical_potentials(
        model,
        state.temperature,
        state.pressure,
        state.composition,
        phase,
        standard_state,
    )
    rt = R * state.temperature
    dreduced_mu_dlogits = dmu_dlogits / rt
    dreduced_mu_dindependent = dmu_dindependent / rt
    dreduced_mu_dmoles = dmu_dmoles / rt
    dreduced_mu_dt = dmu_dt / rt - chemical_potential / (R * state.temperature.square())
    dreduced_mu_dp = dmu_dp / rt

    return ThermodynamicDerivatives(
        dfugacity_coefficient_dlogits=dphi_dlogits,
        dfugacity_coefficient_dindependent_composition=dphi_dindependent,
        dfugacity_coefficient_dtemperature=dphi_dt,
        dfugacity_coefficient_dpressure=dphi_dp,
        dfugacity_coefficient_dmoles=dphi_dmoles,
        dlog_fugacity_coefficient_dlogits=dlog_phi_dlogits,
        dlog_fugacity_coefficient_dindependent_composition=dlog_phi_dindependent,
        dlog_fugacity_coefficient_dtemperature=dlog_phi_dt,
        dlog_fugacity_coefficient_dpressure=dlog_phi_dp,
        dlog_fugacity_coefficient_dmoles=dlog_phi_dmoles,
        dfugacity_dlogits=dfugacity_dlogits,
        dfugacity_dindependent_composition=dfugacity_dindependent,
        dfugacity_dtemperature=dfugacity_dt,
        dfugacity_dpressure=dfugacity_dp,
        dfugacity_dmoles=dfugacity_dmoles,
        dlog_fugacity_dlogits=dlog_fugacity_dlogits,
        dlog_fugacity_dindependent_composition=dlog_fugacity_dindependent,
        dlog_fugacity_dtemperature=dlog_fugacity_dt,
        dlog_fugacity_dpressure=dlog_fugacity_dp,
        dlog_fugacity_dmoles=dlog_fugacity_dmoles,
        dchemical_potential_dlogits=dmu_dlogits,
        dchemical_potential_dindependent_composition=dmu_dindependent,
        dchemical_potential_dtemperature=dmu_dt,
        dchemical_potential_dpressure=dmu_dp,
        dchemical_potential_dmoles=dmu_dmoles,
        dreduced_chemical_potential_dlogits=dreduced_mu_dlogits,
        dreduced_chemical_potential_dindependent_composition=dreduced_mu_dindependent,
        dreduced_chemical_potential_dtemperature=dreduced_mu_dt,
        dreduced_chemical_potential_dpressure=dreduced_mu_dp,
        dreduced_chemical_potential_dmoles=dreduced_mu_dmoles,
        dmolar_volume_dlogits=dvolume_dlogits,
        dmolar_volume_dindependent_composition=dvolume_dindependent,
        dmolar_volume_dtemperature=dvolume_dt,
        dmolar_volume_dpressure=dvolume_dp,
        dmolar_volume_dmoles=dvolume_dmoles,
        dgibbs_dtemperature=dg_dt,
        dgibbs_dpressure=dg_dp,
    )

molar_enthalpy_of_mixing

molar_enthalpy_of_mixing(
    model: StateModel,
    temperature: Tensor,
    pressure: Tensor,
    composition: Tensor,
    phase: PhaseKind = "stable",
) -> Tensor

Return the isothermal-isobaric molar enthalpy change on mixing.

Parameters:

Name	Type	Description	Default
model	StateModel	Homogeneous-state model used for the mixture and pure-component reference states.	required
temperature	Tensor	Finite positive scalar temperature in K and pressure in Pa.	required
pressure	Tensor	Finite positive scalar temperature in K and pressure in Pa.	required
composition	Tensor	Nonnegative mole fractions with components on the final axis. Leading dimensions are evaluated as a batch of independent compositions at the supplied temperature and pressure.	required
phase	PhaseKind	Root-selection request applied independently to the mixture and every pure-component reference state.	'stable'

Returns:

Type	Description
Tensor	Molar enthalpy of mixing in J/mol. The result has the leading shape of composition and is scalar for a single composition, .. math:: h^\mathrm{M}(T,P,\mathbf{x}) = h(T,P,\mathbf{x}) - \sum_i x_i h_i(T,P).

Raises:

Type	Description
ValueError	If temperature, pressure, composition, or root selection is invalid.

Notes

Ideal-gas component enthalpies cancel exactly between the mixture and its unmixed pure-component reference states. The implementation therefore evaluates only residual enthalpies from

.. math::

h^\mathrm{R}
= -R T^2
  \left(\frac{\partial \sum_i x_i \ln\phi_i}
  {\partial T}\right)_{P,\mathbf{x}},

with the temperature derivative obtained by PyTorch automatic differentiation. This preserves derivatives with respect to equation-of- state parameters and avoids requiring an arbitrary ideal-gas caloric reference.

The definition matches the fixed-temperature, fixed-pressure enthalpy change on mixing used to parameterize E-PPR78 by Xu, Lasala, Privat, and Jaubert, International Journal of Greenhouse Gas Control 56 (2017) 126-154, doi:10.1016/j.ijggc.2016.11.015. The requested homogeneous root must represent the experimental single-phase state; this function does not perform a flash or silently average phase enthalpies inside a two-phase region.

Source code in src/torch_flash/properties/thermal.py

def molar_enthalpy_of_mixing(
    model: StateModel,
    temperature: Tensor,
    pressure: Tensor,
    composition: Tensor,
    phase: PhaseKind = "stable",
) -> Tensor:
    r"""Return the isothermal-isobaric molar enthalpy change on mixing.

    Parameters
    ----------
    model
        Homogeneous-state model used for the mixture and pure-component
        reference states.
    temperature, pressure
        Finite positive scalar temperature in K and pressure in Pa.
    composition
        Nonnegative mole fractions with components on the final axis. Leading
        dimensions are evaluated as a batch of independent compositions at
        the supplied temperature and pressure.
    phase
        Root-selection request applied independently to the mixture and every
        pure-component reference state.

    Returns
    -------
    Tensor
        Molar enthalpy of mixing in J/mol. The result has the leading shape of
        ``composition`` and is scalar for a single composition,

        .. math::

            h^\mathrm{M}(T,P,\mathbf{x})
            = h(T,P,\mathbf{x})
            - \sum_i x_i h_i(T,P).

    Raises
    ------
    ValueError
        If temperature, pressure, composition, or root selection is invalid.

    Notes
    -----
    Ideal-gas component enthalpies cancel exactly between the mixture and its
    unmixed pure-component reference states. The implementation therefore
    evaluates only residual enthalpies from

    .. math::

        h^\mathrm{R}
        = -R T^2
          \left(\frac{\partial \sum_i x_i \ln\phi_i}
          {\partial T}\right)_{P,\mathbf{x}},

    with the temperature derivative obtained by PyTorch automatic
    differentiation. This preserves derivatives with respect to equation-of-
    state parameters and avoids requiring an arbitrary ideal-gas caloric
    reference.

    The definition matches the fixed-temperature, fixed-pressure enthalpy
    change on mixing used to parameterize E-PPR78 by Xu, Lasala, Privat, and
    Jaubert, *International Journal of Greenhouse Gas Control* 56 (2017)
    126-154, doi:10.1016/j.ijggc.2016.11.015. The requested homogeneous root
    must represent the experimental single-phase state; this function does
    not perform a flash or silently average phase enthalpies inside a
    two-phase region.
    """
    if phase not in ("liquid", "vapor", "stable"):
        raise ValueError("enthalpy of mixing phase must be 'liquid', 'vapor', or 'stable'")
    if temperature.ndim != 0 or not bool(torch.isfinite(temperature) & (temperature > 0.0)):
        raise ValueError("enthalpy of mixing requires one finite positive temperature")
    if pressure.ndim != 0 or not bool(torch.isfinite(pressure) & (pressure > 0.0)):
        raise ValueError("enthalpy of mixing requires one finite positive pressure")
    if not composition.is_floating_point() or composition.ndim < 1:
        raise ValueError(
            "enthalpy of mixing composition must be a floating tensor with a final component axis"
        )
    if not bool(torch.isfinite(composition).all() & (composition >= 0.0).all()):
        raise ValueError("enthalpy of mixing composition must be finite and nonnegative")

    normalized = normalize_composition(composition)
    component_count = normalized.shape[-1]
    if component_count < 2:
        raise ValueError("enthalpy of mixing requires at least two components")
    solve_temperature = temperature.to(dtype=normalized.dtype, device=normalized.device)
    solve_pressure = pressure.to(dtype=normalized.dtype, device=normalized.device)
    mixture_enthalpy, _ = _residual_terms(
        model,
        solve_temperature,
        solve_pressure,
        normalized,
        phase,
    )
    pure_compositions = torch.eye(
        component_count,
        dtype=normalized.dtype,
        device=normalized.device,
    )
    pure_enthalpies, _ = _residual_terms(
        model,
        solve_temperature,
        solve_pressure,
        pure_compositions,
        phase,
    )
    return mixture_enthalpy - torch.sum(normalized * pure_enthalpies, dim=-1)

thermal_properties

thermal_properties(
    model: StateModel,
    state: ChemicalState,
    standard_state: CaloricStandardState,
    phase: PhaseKind = "stable",
    *,
    reference_pressure: float = STANDARD_PRESSURE,
    molar_mass: Tensor | None = None,
) -> ThermalProperties

Evaluate Pedersen Chapter 8 properties at a specified state.

Parameters:

Name	Type	Description	Default
model	StateModel	Homogeneous-state model.	required
state	ChemicalState	One scalar TP state with temperature in K, pressure in Pa, and one composition vector.	required
standard_state	CaloricStandardState	Component ideal-gas heat capacity, enthalpy, and entropy functions.	required
phase	PhaseKind	Root-selection request held fixed during differentiation.	'stable'
reference_pressure	float	Positive ideal-gas entropy reference pressure in Pa.	STANDARD_PRESSURE
molar_mass	Tensor | None	Optional scalar mixture or component molar masses in kg/mol. When omitted, the model's molar_mass attribute is used if available.	None

Returns:

Type	Description
ThermalProperties	Caloric energies, response functions, residual terms, and optional speed of sound in SI units.

Raises:

Type	Description
ValueError	If state shapes, reference pressure, standard-state component count, or molar-mass data are invalid.
InvalidStateError	If response functions violate positive heat-capacity or mechanical stability requirements.

Notes

No phase-equilibrium calculation is performed. The requested phase root remains fixed while differentiating. A volume-translated cubic receives Pedersen's P*DeltaV enthalpy correction (Eq. 8.12).

The entropy pressure term in Eq. 8.17 is interpreted as -R*ln(P/Pref); the glyph in the printed 2024 edition can resemble T, but dimensional consistency and the ideal-gas identity require the gas constant. Helmholtz energy follows a = g - P*v. Reduced total free energies use the supplied caloric reference, while reduced residual free energies are reference-independent EoS departures.

Source code in src/torch_flash/properties/thermal.py

def thermal_properties(
    model: StateModel,
    state: ChemicalState,
    standard_state: CaloricStandardState,
    phase: PhaseKind = "stable",
    *,
    reference_pressure: float = STANDARD_PRESSURE,
    molar_mass: Tensor | None = None,
) -> ThermalProperties:
    """Evaluate Pedersen Chapter 8 properties at a specified state.

    Parameters
    ----------
    model
        Homogeneous-state model.
    state
        One scalar TP state with temperature in K, pressure in Pa, and one
        composition vector.
    standard_state
        Component ideal-gas heat capacity, enthalpy, and entropy functions.
    phase
        Root-selection request held fixed during differentiation.
    reference_pressure
        Positive ideal-gas entropy reference pressure in Pa.
    molar_mass
        Optional scalar mixture or component molar masses in kg/mol. When
        omitted, the model's ``molar_mass`` attribute is used if available.

    Returns
    -------
    ThermalProperties
        Caloric energies, response functions, residual terms, and optional
        speed of sound in SI units.

    Raises
    ------
    ValueError
        If state shapes, reference pressure, standard-state component count,
        or molar-mass data are invalid.
    InvalidStateError
        If response functions violate positive heat-capacity or mechanical
        stability requirements.

    Notes
    -----
    No phase-equilibrium calculation is performed. The requested phase root
    remains fixed while differentiating. A volume-translated cubic receives
    Pedersen's ``P*DeltaV`` enthalpy correction (Eq. 8.12).

    The entropy pressure term in Eq. 8.17 is interpreted as
    ``-R*ln(P/Pref)``; the glyph in the printed 2024 edition can resemble
    ``T``, but dimensional consistency and the ideal-gas identity require
    the gas constant. Helmholtz energy follows ``a = g - P*v``. Reduced total
    free energies use the supplied caloric reference, while reduced residual
    free energies are reference-independent EoS departures.
    """
    if state.temperature.ndim != 0 or state.pressure.ndim != 0:
        raise ValueError("thermal_properties currently accepts one scalar T-P state")
    if state.composition.ndim != 1:
        raise ValueError("thermal_properties currently accepts one composition vector")
    if reference_pressure <= 0.0:
        raise ValueError("reference_pressure must be positive")

    temperature = state.temperature
    pressure = state.pressure
    composition = state.composition
    ideal_enthalpy = standard_state.enthalpy(temperature)
    ideal_entropy = standard_state.entropy(temperature)
    if ideal_enthalpy.shape != composition.shape or ideal_entropy.shape != composition.shape:
        raise ValueError("standard-state component count must match the composition")

    residual_enthalpy, residual_entropy = _residual_terms(
        model,
        temperature,
        pressure,
        composition,
        phase,
    )
    enthalpy = torch.sum(composition * ideal_enthalpy) + residual_enthalpy
    safe_composition = torch.clamp_min(composition, torch.finfo(composition.dtype).tiny)
    entropy = (
        torch.sum(composition * ideal_entropy)
        - R * torch.log(pressure / reference_pressure)
        - R * torch.sum(composition * torch.log(safe_composition))
        + residual_entropy
    )
    volume = model.molar_volume(temperature, pressure, composition, phase)
    internal_energy = enthalpy - pressure * volume
    gibbs = enthalpy - temperature * entropy
    helmholtz = gibbs - pressure * volume
    reduced_gibbs = gibbs / (R * temperature)
    reduced_helmholtz = helmholtz / (R * temperature)
    pressure_volume_over_rt = pressure * volume / (R * temperature)
    reduced_residual_gibbs = _weighted_log_fugacity(
        model,
        temperature,
        pressure,
        composition,
        phase,
    )
    reduced_residual_helmholtz = (
        reduced_residual_gibbs - pressure_volume_over_rt + 1.0 + torch.log(pressure_volume_over_rt)
    )

    enthalpy_at_temperature = lambda current_temperature: _enthalpy(  # noqa: E731
        model,
        standard_state,
        current_temperature,
        pressure,
        composition,
        phase,
    )
    enthalpy_at_pressure = lambda current_pressure: _enthalpy(  # noqa: E731
        model,
        standard_state,
        temperature,
        current_pressure,
        composition,
        phase,
    )
    cp = torch.func.grad(enthalpy_at_temperature)(temperature)
    dh_dp = torch.func.grad(enthalpy_at_pressure)(pressure)
    joule_thomson = -dh_dp / cp

    volume_at_temperature = lambda current_temperature: model.molar_volume(  # noqa: E731
        current_temperature,
        pressure,
        composition,
        phase,
    )
    volume_at_pressure = lambda current_pressure: model.molar_volume(  # noqa: E731
        temperature,
        current_pressure,
        composition,
        phase,
    )
    dv_dt = torch.func.grad(volume_at_temperature)(temperature)
    dv_dp = torch.func.grad(volume_at_pressure)(pressure)
    cv = cp + temperature * dv_dt.square() / dv_dp

    response_values = torch.stack((cp, cv, dv_dp))
    if bool(
        (
            (~torch.isfinite(response_values)).any() | (cp <= 0.0) | (cv <= 0.0) | (dv_dp >= 0.0)
        ).detach()
    ):
        raise InvalidStateError(
            "thermal response functions require positive Cp/Cv and "
            "negative isothermal compressibility"
        )

    mixture_molar_mass = _mixture_molar_mass(model, composition, molar_mass)
    speed_of_sound = None
    if mixture_molar_mass is not None:
        speed_squared = -volume.square() * cp / (mixture_molar_mass * cv * dv_dp)
        if bool((~torch.isfinite(speed_squared) | (speed_squared <= 0.0)).detach()):
            raise InvalidStateError("speed-of-sound expression is not positive and finite")
        speed_of_sound = torch.sqrt(speed_squared)

    return ThermalProperties(
        enthalpy,
        internal_energy,
        entropy,
        helmholtz,
        gibbs,
        reduced_helmholtz,
        reduced_gibbs,
        reduced_residual_helmholtz,
        reduced_residual_gibbs,
        cp,
        cv,
        joule_thomson,
        speed_of_sound,
        residual_enthalpy,
        residual_entropy,
    )

Ideal-gas and chemical-potential standard states

torch_flash.standard_state

Ideal-gas standard-state contributions and caloric reference choices.

The polynomial integrations implement Pedersen, Phase Behavior of Petroleum Reservoir Fluids, 3rd ed. (2024), Eqs. 8.2-8.4, doi:10.1201/9780429457418. The named coefficients are frozen from the Poling data bank distributed by chemicals 1.5.2 and trace back to Poling et al., The Properties of Gases and Liquids, 5th ed. (2001).

StandardState

Bases: Protocol

Protocol for component standard chemical potentials.

Implementations return one value per modeled component, retain leading temperature batch dimensions, and use the same dtype and device as the input temperature.

Methods:

Name	Description
chemical_potential	Return component standard chemical potentials in J/mol.

chemical_potential

chemical_potential(temperature: Tensor) -> Tensor

Return component standard chemical potentials in J/mol.

Source code in src/torch_flash/standard_state.py

def chemical_potential(self, temperature: Tensor) -> Tensor:
    """Return component standard chemical potentials in J/mol."""

IdealGasPolynomial

IdealGasPolynomial(
    heat_capacity_coefficients: Tensor,
    reference_enthalpy: Tensor,
    reference_entropy: Tensor,
    *,
    reference_temperature: float = 298.15,
    trainable: bool = False,
)

Bases: Module

Polynomial ideal-gas heat-capacity standard state.

For m columns, the heat capacity is sum(c[j]*T**j, j=0..m-1) in J/(mol K). Four columns recover Pedersen's Eq. 8.4; five columns support the Poling data bank. Reference enthalpies and entropies make the otherwise arbitrary caloric datum explicit and allow the coefficients to be fitted with PyTorch.

Parameters:

Name	Type	Description	Default
heat_capacity_coefficients	Tensor	Matrix of coefficients in ascending temperature power, with one row per component and SI units that make the evaluated heat capacity J/(mol K).	required
reference_enthalpy	Tensor	Component enthalpies in J/mol at reference_temperature.	required
reference_entropy	Tensor	Component standard-pressure entropies in J/(mol K) at the reference.	required
reference_temperature	float	Positive caloric reference temperature in K.	298.15
trainable	bool	Register heat-capacity coefficients as trainable parameters.	False

Source code in src/torch_flash/standard_state.py

def __init__(
    self,
    heat_capacity_coefficients: Tensor,
    reference_enthalpy: Tensor,
    reference_entropy: Tensor,
    *,
    reference_temperature: float = 298.15,
    trainable: bool = False,
) -> None:
    super().__init__()
    if heat_capacity_coefficients.ndim != 2 or heat_capacity_coefficients.shape[1] < 1:
        raise ValueError("heat-capacity coefficients must be a nonempty matrix")
    expected = (heat_capacity_coefficients.shape[0],)
    if reference_enthalpy.shape != expected or reference_entropy.shape != expected:
        raise ValueError("one reference enthalpy and entropy are required per component")
    if reference_temperature <= 0.0:
        raise ValueError("reference temperature must be positive")
    self.reference_temperature = float(reference_temperature)
    if trainable:
        self.heat_capacity_coefficients = nn.Parameter(heat_capacity_coefficients.clone())
    else:
        self.register_buffer(
            "heat_capacity_coefficients",
            heat_capacity_coefficients.clone(),
        )
    self.register_buffer("reference_enthalpy", reference_enthalpy.clone())
    self.register_buffer("reference_entropy", reference_entropy.clone())

heat_capacity

heat_capacity(temperature: Tensor) -> Tensor

Evaluate component ideal-gas heat capacities.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Temperature in K with arbitrary leading batch dimensions.	required

Returns:

Type	Description
Tensor	Heat capacities in J/(mol K), with a final component axis.

Source code in src/torch_flash/standard_state.py

def heat_capacity(self, temperature: Tensor) -> Tensor:
    """Evaluate component ideal-gas heat capacities.

    Parameters
    ----------
    temperature
        Temperature in K with arbitrary leading batch dimensions.

    Returns
    -------
    Tensor
        Heat capacities in J/(mol K), with a final component axis.
    """
    powers = torch.arange(
        self.heat_capacity_coefficients.shape[1],
        dtype=temperature.dtype,
        device=temperature.device,
    )
    return torch.sum(
        self.heat_capacity_coefficients * temperature[..., None, None].pow(powers),
        dim=-1,
    )

enthalpy

enthalpy(temperature: Tensor) -> Tensor

Integrate component ideal-gas enthalpies from the reference state.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Temperature in K.	required

Returns:

Type	Description
Tensor	Component enthalpies in J/mol with a final component axis.

Source code in src/torch_flash/standard_state.py

def enthalpy(self, temperature: Tensor) -> Tensor:
    """Integrate component ideal-gas enthalpies from the reference state.

    Parameters
    ----------
    temperature
        Temperature in K.

    Returns
    -------
    Tensor
        Component enthalpies in J/mol with a final component axis.
    """
    powers = torch.arange(
        1,
        self.heat_capacity_coefficients.shape[1] + 1,
        dtype=temperature.dtype,
        device=temperature.device,
    )
    reference = temperature.new_tensor(self.reference_temperature)
    integral = (
        self.heat_capacity_coefficients
        * (temperature[..., None, None].pow(powers) - reference.pow(powers))
        / powers
    )
    return self.reference_enthalpy + torch.sum(integral, dim=-1)

entropy

entropy(temperature: Tensor) -> Tensor

Integrate component ideal-gas entropies at standard pressure.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Positive temperature in K.	required

Returns:

Type	Description
Tensor	Component standard-state entropies in J/(mol K).

Source code in src/torch_flash/standard_state.py

def entropy(self, temperature: Tensor) -> Tensor:
    """Integrate component ideal-gas entropies at standard pressure.

    Parameters
    ----------
    temperature
        Positive temperature in K.

    Returns
    -------
    Tensor
        Component standard-state entropies in J/(mol K).
    """
    reference = temperature.new_tensor(self.reference_temperature)
    coefficients = self.heat_capacity_coefficients
    leading = coefficients[:, 0] * torch.log(temperature[..., None] / reference)
    powers = torch.arange(
        1,
        coefficients.shape[1],
        dtype=temperature.dtype,
        device=temperature.device,
    )
    remaining = (
        coefficients[:, 1:]
        * (temperature[..., None, None].pow(powers) - reference.pow(powers))
        / powers
    )
    return self.reference_entropy + leading + torch.sum(remaining, dim=-1)

chemical_potential

chemical_potential(temperature: Tensor) -> Tensor

Return ideal-gas standard chemical potentials.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Temperature in K.	required

Returns:

Type	Description
Tensor	h_i - T s_i in J/mol with a final component axis.

Source code in src/torch_flash/standard_state.py

def chemical_potential(self, temperature: Tensor) -> Tensor:
    """Return ideal-gas standard chemical potentials.

    Parameters
    ----------
    temperature
        Temperature in K.

    Returns
    -------
    Tensor
        ``h_i - T s_i`` in J/mol with a final component axis.
    """
    return self.enthalpy(temperature) - temperature[..., None] * self.entropy(temperature)

ideal_gas_polynomial

ideal_gas_polynomial(
    names: tuple[str, ...] | list[str],
    parameter_set: ParameterSource,
    *,
    dtype: dtype | None = None,
    device: device | str | None = None,
    reference_temperature: float | None = None,
    trainable: bool = False,
) -> IdealGasPolynomial

Construct an ideal-gas polynomial from versioned parameters.

Parameters:

Name	Type	Description	Default
names	tuple[str, ...] | list[str]	Component names in the desired output order. Aliases are canonicalized through the component database.	required
parameter_set	ParameterSource	Bundled identifier, YAML path, mapping, or loaded standard-state parameter set.	required
dtype	dtype | None	Tensor placement. Omitted values use the configured runtime defaults.	None
device	dtype | None	Tensor placement. Omitted values use the configured runtime defaults.	None
reference_temperature	float | None	Caloric reference temperature in K. When omitted, the parameter set's default is required.	None
trainable	bool	Register heat-capacity coefficients as trainable parameters. Reference enthalpies and entropies remain buffers.	False

Returns:

Type	Description
IdealGasPolynomial	Standard-state module ordered according to names.

Raises:

Type	Description
ParameterDatabaseError	If the source has the wrong model kind or malformed/missing entries.
KeyError	If a requested component is absent from the parameter set.

Source code in src/torch_flash/standard_state.py

def ideal_gas_polynomial(
    names: tuple[str, ...] | list[str],
    parameter_set: ParameterSource,
    *,
    dtype: torch.dtype | None = None,
    device: torch.device | str | None = None,
    reference_temperature: float | None = None,
    trainable: bool = False,
) -> IdealGasPolynomial:
    """Construct an ideal-gas polynomial from versioned parameters.

    Parameters
    ----------
    names
        Component names in the desired output order. Aliases are canonicalized
        through the component database.
    parameter_set
        Bundled identifier, YAML path, mapping, or loaded standard-state
        parameter set.
    dtype, device
        Tensor placement. Omitted values use the configured runtime defaults.
    reference_temperature
        Caloric reference temperature in K. When omitted, the parameter set's
        default is required.
    trainable
        Register heat-capacity coefficients as trainable parameters. Reference
        enthalpies and entropies remain buffers.

    Returns
    -------
    IdealGasPolynomial
        Standard-state module ordered according to ``names``.

    Raises
    ------
    ParameterDatabaseError
        If the source has the wrong model kind or malformed/missing entries.
    KeyError
        If a requested component is absent from the parameter set.
    """
    dtype, device = resolve_tensor_options(dtype, device)
    loaded = load_model_parameters(parameter_set)
    if loaded.model_kind != "standard_state":
        raise ParameterDatabaseError(
            f"{loaded.identifier!r} is {loaded.model_kind!r}, not 'standard_state'"
        )
    records = loaded.parameters.get("components")
    if not isinstance(records, Mapping):
        raise ParameterDatabaseError(f"{loaded.identifier!r} requires a components mapping")
    canonical = tuple(component(name).name for name in names)
    coefficients: list[list[float]] = []
    reference_enthalpy: list[float] = []
    reference_entropy: list[float] = []
    for name in canonical:
        record = records.get(name)
        if not isinstance(record, Mapping):
            raise KeyError(f"no {loaded.identifier} heat-capacity coefficients for: {name}")
        values = record.get("coefficients")
        if (
            not isinstance(values, Sequence)
            or isinstance(values, str)
            or not values
            or any(not isinstance(value, int | float) for value in values)
        ):
            raise ParameterDatabaseError(
                f"{loaded.identifier!r} coefficients for {name!r} must be numeric"
            )
        coefficients.append([float(value) for value in values])
        enthalpy = record.get("reference_enthalpy", 0.0)
        entropy = record.get("reference_entropy", 0.0)
        if not isinstance(enthalpy, int | float) or not isinstance(entropy, int | float):
            raise ParameterDatabaseError(
                f"{loaded.identifier!r} reference properties for {name!r} must be numeric"
            )
        reference_enthalpy.append(float(enthalpy))
        reference_entropy.append(float(entropy))
    if len({len(values) for values in coefficients}) != 1:
        raise ParameterDatabaseError("ideal-gas polynomial orders must match")
    resolved_reference = reference_temperature
    if resolved_reference is None:
        stored_reference = loaded.parameters.get("default_reference_temperature")
        if not isinstance(stored_reference, int | float):
            raise ParameterDatabaseError(
                f"{loaded.identifier!r} requires default_reference_temperature"
            )
        resolved_reference = float(stored_reference)
    dimensionless = torch.tensor(coefficients, dtype=dtype, device=device)
    return IdealGasPolynomial(
        R * dimensionless,
        torch.tensor(reference_enthalpy, dtype=dtype, device=device),
        torch.tensor(reference_entropy, dtype=dtype, device=device),
        reference_temperature=resolved_reference,
        trainable=trainable,
    )

poling_ideal_gas

poling_ideal_gas(
    names: tuple[str, ...] | list[str],
    *,
    dtype: dtype | None = None,
    device: device | str | None = None,
    reference_temperature: float = 273.15,
    trainable: bool = False,
) -> IdealGasPolynomial

Construct a named Poling ideal-gas standard state.

Parameters:

Name	Type	Description	Default
names	tuple[str, ...] | list[str]	Component names in desired tensor order.	required
dtype	dtype | None	Tensor placement.	None
device	dtype | None	Tensor placement.	None
reference_temperature	float	Caloric reference temperature in K.	273.15
trainable	bool	Register heat-capacity coefficients as trainable parameters.	False

Returns:

Type	Description
IdealGasPolynomial	Poling coefficient model for the selected components.

Raises:

Type	Description
KeyError	If a selected component has no frozen Poling coefficient record.

Notes

The tabulated fits cover 50-1000 K for gases through propane and 200-1000 K for C4-C10. Argon and helium are the monatomic ideal-gas limits. This function deliberately does not extrapolate or clip temperatures; callers remain responsible for checking those source ranges in their application.

Source code in src/torch_flash/standard_state.py

def poling_ideal_gas(
    names: tuple[str, ...] | list[str],
    *,
    dtype: torch.dtype | None = None,
    device: torch.device | str | None = None,
    reference_temperature: float = 273.15,
    trainable: bool = False,
) -> IdealGasPolynomial:
    """Construct a named Poling ideal-gas standard state.

    Parameters
    ----------
    names
        Component names in desired tensor order.
    dtype, device
        Tensor placement.
    reference_temperature
        Caloric reference temperature in K.
    trainable
        Register heat-capacity coefficients as trainable parameters.

    Returns
    -------
    IdealGasPolynomial
        Poling coefficient model for the selected components.

    Raises
    ------
    KeyError
        If a selected component has no frozen Poling coefficient record.

    Notes
    -----
    The tabulated fits cover 50-1000 K for gases through propane and
    200-1000 K for C4-C10. Argon and helium are the monatomic ideal-gas
    limits. This function deliberately does not extrapolate or clip
    temperatures; callers remain responsible for checking those source
    ranges in their application.
    """
    try:
        return ideal_gas_polynomial(
            names,
            "standard-state.poling-2001",
            dtype=dtype,
            device=device,
            reference_temperature=reference_temperature,
            trainable=trainable,
        )
    except KeyError as exc:
        canonical = tuple(component(name).name for name in names)
        raise KeyError(
            f"no frozen Poling heat-capacity coefficients for: {', '.join(canonical)}"
        ) from exc