Getting started

This guide takes a new user from installation to the calculations most often needed in a thermodynamic workflow:

1. configure PyTorch and construct a model;
2. evaluate a supplied homogeneous state and its derivatives;
3. solve a known two-phase state;
4. solve a fixed three-phase problem; and
5. check convergence, fugacity residuals, and material balance.

All public thermodynamic inputs and outputs use SI units. The examples use torch.float64, the reference precision for phase-equilibrium calculations.

Install torch-flash

torch-flash requires Python 3.11 or newer. Create an isolated Conda environment and install the package from Conda-Forge:

conda create --name torch-flash -c conda-forge python=3.11 torch-flash
conda activate torch-flash

To install into an existing Conda environment instead, use:

conda install -c conda-forge torch-flash

The default package is also available from PyPI:

python -m pip install torch-flash

The Conda-Forge package currently provides the default capability. Optional package capabilities are available as PyPI extras; these remain intentionally narrower than the Pixi development environments:

# Group-contribution helpers
python -m pip install "torch-flash[groups]"

# Optional sparse linear algebra on Linux or Windows
python -m pip install "torch-flash[intel]"

# Optional CUDA-oriented solver integrations on Linux or Windows
python -m pip install "torch-flash[gpu]"

The gpu extra does not install a GPU driver or make every calculation faster. Use a PyTorch build compatible with the local CUDA driver, keep tensors on the device, and batch enough work to amortize launch and transfer costs. See high-performance setup before selecting an accelerator.

For a source checkout, use the locked Pixi environment:

pixi install -e default
pixi run -e default python -c "import torch_flash; print(torch_flash.__version__)"

Configure the runtime first

Runtime dtype, device, thread counts, and deterministic behavior affect every model tensor created afterward. Configure them before constructing component sets or models:

import torch

from torch_flash import configure

runtime = configure(
    device="cpu",
    dtype=torch.float64,
    num_threads=4,
)
temperature = runtime.tensor(300.0)

device="auto" selects the first available CUDA, XPU, or MPS device that supports the chosen dtype, then falls back to CPU. device="gpu" requires an accelerator and raises if none supports that dtype. Float64 is not available on every accelerator; do not silently replace it with float32 near critical states or coalescing roots.

The model-state-calculation pattern

Most torch-flash calculations follow the same pattern:

import torch

from torch_flash import ChemicalState, component_set, configure, peng_robinson_1978

runtime = configure(device="cpu", dtype=torch.float64)
components = component_set(("methane", "n_butane"))
model = peng_robinson_1978(components)
state = ChemicalState(
    temperature=runtime.tensor(300.0),  # K
    pressure=runtime.tensor(5.0e6),  # Pa
    composition=runtime.tensor([0.70, 0.30]),  # mol/mol
)

The order of every composition vector is the order supplied to component_set. Component names are canonicalized at this API boundary. Named constructors use the versioned parameter databases documented in the parameter guide.

Compute properties and derivatives

phase_properties evaluates the supplied homogeneous state. It does not perform a flash. state_derivatives then differentiates the same state with respect to temperature, pressure, composition coordinates, and mole numbers.

"""Evaluate homogeneous phase properties and their TP derivatives."""

import torch

from torch_flash import (
    ChemicalState,
    component_set,
    configure,
    peng_robinson_1978,
    phase_properties,
    state_derivatives,
)

runtime = configure(device="cpu", dtype=torch.float64)
model = peng_robinson_1978(component_set(("methane", "n_butane")))
state = ChemicalState(
    temperature=runtime.tensor(300.0),  # K
    pressure=runtime.tensor(5.0e6),  # Pa
    composition=runtime.tensor([0.70, 0.30]),  # mol/mol
)

# This evaluates the supplied homogeneous state; it does not run a flash.
properties = phase_properties(model, state, phase="stable")
derivatives = state_derivatives(model, state, phase="stable")

print(f"Z = {float(properties.compressibility_factor):.6f}")
print(f"molar volume = {float(properties.molar_volume):.6e} m^3/mol")
print("fugacities [Pa] =", properties.fugacities)
print("d ln(f_i/p°) / dT [1/K] =", derivatives.dlog_fugacity_dtemperature)
print("dV_m / dP [m^3/(mol Pa)] =", derivatives.dmolar_volume_dpressure)

# For a binary mixture, the single independent coordinate is x_methane;
# x_n-butane = 1 - x_methane.
print(
    "d mu_i / d x_methane [J/mol] =",
    derivatives.dchemical_potential_dindependent_composition[:, 0],
)

The default chemical-potential convention uses a zero ideal-gas standard chemical potential at 1 bar. It is not a formation-property reference. Reference-independent residual energies are provided separately.

Composition derivatives require strictly positive mole fractions. For \(n\) components, the independent-composition Jacobian has \(n-1\) columns, with the final fraction defined by \(x_n=1-\sum_{i=1}^{n-1}x_i\). See Phase Properties and Their Derivatives for the complete coordinate and unit conventions.

Solve a two-phase flash

At fixed temperature, pressure, and overall composition, a two-phase TP flash finds the phase fraction and phase compositions that satisfy material balance and equality of component fugacities.

"""Solve a known two-phase methane/n-butane TP state."""

import torch

from torch_flash import (
    ChemicalState,
    component_set,
    configure,
    peng_robinson_1978,
    two_phase_flash,
)

runtime = configure(device="cpu", dtype=torch.float64)
model = peng_robinson_1978(component_set(("methane", "n_butane")))
state = ChemicalState(
    temperature=runtime.tensor(270.0),  # K
    pressure=runtime.tensor(3.0e6),  # Pa
    composition=runtime.tensor([0.50, 0.50]),  # overall z
)

# This example deliberately solves a state already known to be on a
# two-phase branch. Leave check_stability=True when the phase count is unknown.
result = two_phase_flash(
    model,
    state,
    check_stability=False,
    tolerance=1.0e-10,
    raise_on_failure=True,
)

liquid, vapor = result.phases
print("phase fractions [liquid, vapor] =", result.phase_fractions)
print("liquid x =", liquid.composition)
print("vapor y =", vapor.composition)
print("maximum log-fugacity residual =", float(result.residual_norm))

# Always check material balance and the numerical result before using it.
reconstructed_feed = (
    result.phase_fractions[0] * liquid.composition + result.phase_fractions[1] * vapor.composition
)
torch.testing.assert_close(reconstructed_feed, state.composition, atol=1.0e-9, rtol=0.0)
assert result.converged

This compact example sets check_stability=False because it deliberately solves a state already known to lie on a two-phase branch. Keep the default check_stability=True when the phase count is unknown. Disabling stability analysis is not a general shortcut: fugacity equality is necessary, but a lower-Gibbs-energy phase may still exist.

Do not use a result unless:

• result.converged is true;
• result.residual_norm is appropriate for the model and dtype;
• all phase fractions and compositions are physical;
• the reconstructed feed matches the specified feed; and
• the selected roots and phase identities make physical sense.

The Two-Phase Flash page derives the Rachford-Rice balance and documents the solver sequence.

Solve a fixed multiphase flash

multiphase_flash accepts one row of equilibrium ratios for every phase after the reference phase. Therefore, an initial \(K\) matrix with shape (2, n_components) requests a fixed three-phase solve.

The following example reproduces the model definition and rounded initialization from Tables 6.5-6.6 of Pedersen, Christensen, and Shaikh. The source values initialize the calculation; the converged result is checked by fugacity residual and material balance.

Complete runnable fixed-three-phase example

"""Solve a fixed three-phase flash initialized from a published example."""

import torch

from torch_flash import ChemicalState, ComponentSet, multiphase_flash, peng_robinson_1978

# Pedersen, Christensen, and Shaikh (2024), Tables 6.5-6.6. The published
# constants and phase compositions are rounded, so they are used here as a
# physically informed initializer rather than an exact numerical baseline.
names = (
    "nitrogen",
    "carbon_dioxide",
    "methane",
    "ethane",
    "propane",
    "isobutane",
    "n_butane",
    "isopentane",
    "n_pentane",
    "n_hexane",
    "n_heptane",
    "n_octane",
    "n_decane",
)
components = ComponentSet(
    names,
    (
        torch.tensor(
            [
                -147.0,
                31.1,
                -82.6,
                32.3,
                96.7,
                135.0,
                152.1,
                187.3,
                196.5,
                234.3,
                280.4,
                352.5,
                473.1,
            ],
            dtype=torch.float64,
        )
        + 273.15
    ),
    torch.tensor(
        [
            33.94,
            73.76,
            46.00,
            48.84,
            42.46,
            36.48,
            38.00,
            33.84,
            33.74,
            29.69,
            26.72,
            21.29,
            16.67,
        ],
        dtype=torch.float64,
    )
    * 1.0e5,
    torch.tensor(
        [
            0.040,
            0.225,
            0.008,
            0.098,
            0.152,
            0.176,
            0.193,
            0.227,
            0.251,
            0.296,
            0.373,
            0.518,
            0.803,
        ],
        dtype=torch.float64,
    ),
    torch.ones(13, dtype=torch.float64),
)

feed = torch.tensor(
    [0.08, 2.01, 82.51, 5.81, 2.88, 0.56, 1.24, 0.52, 0.60, 0.72, 1.66, 0.91, 0.49],
    dtype=torch.float64,
)
vapor_guess = torch.tensor(
    [
        0.18,
        1.08,
        96.45,
        1.86,
        0.33,
        0.03,
        0.05,
        0.01,
        0.01,
        1.0e-4,
        1.0e-4,
        1.0e-4,
        1.0e-4,
    ],
    dtype=torch.float64,
)
liquid_1_guess = torch.tensor(
    [0.08, 1.88, 87.95, 5.28, 2.17, 0.38, 0.76, 0.28, 0.29, 0.29, 0.48, 0.14, 0.01],
    dtype=torch.float64,
)
liquid_2_guess = torch.tensor(
    [0.05, 2.36, 75.66, 7.28, 4.00, 0.81, 1.83, 0.79, 0.93, 1.14, 2.73, 1.54, 0.87],
    dtype=torch.float64,
)

reference = liquid_1_guess / liquid_1_guess.sum()
initial_k = torch.stack(
    (
        vapor_guess / vapor_guess.sum() / reference,
        liquid_2_guess / liquid_2_guess.sum() / reference,
    )
)
state = ChemicalState(
    torch.tensor(201.15, dtype=torch.float64),
    torch.tensor(5.2e6, dtype=torch.float64),
    feed,
)
result = multiphase_flash(
    peng_robinson_1978(components),
    state,
    initial_k_values=initial_k,
    tolerance=1.0e-10,
    max_iterations=30,
)

print("phase fractions =", result.phase_fractions)
print("phase kinds =", result.phase_kinds)
print("maximum log-fugacity residual =", float(result.residual_norm))
assert result.converged
assert result.nphases == 3

reconstructed_feed = sum(
    fraction * phase.composition
    for fraction, phase in zip(result.phase_fractions, result.phases, strict=True)
)
torch.testing.assert_close(reconstructed_feed, state.composition, atol=2.0e-8, rtol=0.0)

The emitted ExperimentalModelWarning is intentional. The current API does not discover the correct phase count automatically, and a converged fixed-phase solution is not by itself a global-stability proof. Physical phase identification is also a diagnostic rather than an equilibrium equation; with the rounded constants in this example, the current heuristic labels the converged phases vapor, vapor, and liquid. Inspect densities, roots, stability, and the application context before assigning physical names.

See Multiphase Flash for the generalized Rachford-Rice equations, initialization contract, and limitations.

Differentiate a custom objective

Native model calculations remain in the PyTorch graph. For example, an optimization objective can depend directly on a fugacity coefficient:

import torch

from torch_flash import component_set, configure, peng_robinson_1978

runtime = configure(device="cpu", dtype=torch.float64)
model = peng_robinson_1978(component_set(("methane", "n_butane")))
temperature = runtime.tensor(300.0, requires_grad=True)
pressure = runtime.tensor(5.0e6)
composition = runtime.tensor([0.70, 0.30])

log_phi = model.log_fugacity_coefficients(
    temperature,
    pressure,
    composition,
    phase="stable",
)
objective = log_phi.square().sum()
objective.backward()
print(temperature.grad)

Avoid .item(), NumPy conversion, or detach() inside a differentiable calculation. Use them only at an intentional reporting or control-flow boundary.

Where to go next

• Chemical Equilibrium Overview explains equilibrium, fugacity equality, stability, and the tangent-plane distance.
• Activity Models covers NRTL, Wilson, original UNIFAC, and Huron-Vidal activity terms.
• Fugacity Models connects equation-of-state roots to phase fugacity.
• Characterization and Pseudo-Components covers plus-fraction splitting, property correlations, and lumping.
• High-performance setup covers batching, compilation, threading, accelerators, gradients, and trustworthy timing.
• Model scope and validation evidence define which scientific claims are currently supported.