ExoEOS

Differentiable equations of state and excess free-energy models for planetary atmospheres, fluids, and melts, built with JAX.

Residual equation-of-state models use reduced residual Helmholtz energy as their source of truth. Mole-fraction solution models use reduced molar excess Gibbs energy, from which logarithmic activity coefficients are obtained by automatic differentiation. Total Helmholtz derivatives provide caloric properties and response functions for fluids with a supplied ideal closure. A calorically perfect ideal-gas mixture is also available through the original temperature-pressure state interface.

Installation

python -m pip install "exoeos>=0.2.0"

Version 0.2.0 includes the residual Helmholtz, TP inversion, excess Gibbs, and total Helmholtz thermodynamics APIs documented below, alongside the caloric ideal-gas API from 0.1.0.

To install from a repository checkout:

python -m pip install .

Documentation

Start with the capability overview and model selection and usage guide. They distinguish package APIs, checkout examples, and the evidence supporting each model. The plot gallery shows how each capability's outputs change with temperature, pressure or composition, with reproducible code and explicit source conditions. The Helmholtz derivative guide covers heat capacities, sound speed and atmospheric/RCE use. The Japanese explanation is maintained separately; current English documentation lives in documents/.

Install the documentation dependencies and build the committed notebook-based tutorials with:

python -m pip install -e ".[docs]"
./update_doc.sh

The independent-reference notebooks additionally require the pinned reference backend:

python -m pip install -e ".[docs,reference]"

The executable notebooks are the editable sources. After changing one, run it with jupyter nbconvert --to notebook --execute --inplace <notebook>, then run python documents/tutorials/convert_notebooks.py to refresh its committed RST and image assets.

Residual Helmholtz API

Models implement

alphar(T, rho, x) = A_res / (n R T),

where rho is molar density in mol m-3. The initial IdealEOS is a zero residual placeholder behind this interface. A_res excludes the complete ideal-gas Helmholtz contribution.

import jax
import jax.numpy as jnp

from exoeos import IdealEOS, psir, state_trho


eos = IdealEOS()
T = 1000.0
rho = 12.0
x = jnp.array([0.85, 0.15])

alpha_r = eos.alphar(T, rho, x)
psi_r = psir(eos, T, rho * x)
state = state_trho(eos, T, rho, x)

state.P
state.Z
state.alphar
state.mu_res_RT
state.lnphi
state.gres_RT

batched_Z = jax.vmap(state_trho, in_axes=(None, 0, 0, 0))(
    eos,
    jnp.array([800.0, 1000.0]),
    jnp.array([10.0, 12.0]),
    jnp.array([[0.9, 0.1], [0.85, 0.15]]),
).Z

IdealEOS also implements pressure inversion, so it can be compared with other TP-capable residual models through the same entry point:

from exoeos import state_tp

ideal_state = state_tp(IdealEOS(), T=1000.0, P=1.0e5, x=x)

alphar and state_trho accept one state at a time: scalar T, scalar rho, and x with shape (K,). psir instead accepts the component molar density vector rho_vec with shape (K,). Use jax.vmap for batches. Compositions are neither clipped nor normalized.

Models that can invert pressure additionally implement TPHelmholtzEOS by providing molar_density(T, P, x, phase="vapor"). The common state_tp(eos, T, P, x, phase="vapor") entry point delegates density and root selection to that hook, then evaluates state_trho and returns the same TRhoState. This gives ExoGibbs both molar density and fugacity coefficients from its natural T, P, and x inputs.

Mass density follows from the same TP density hook when component molar masses are supplied in kg mol-1:

from exoeos import additive_volume_mass_density, mass_density_tp


rho = mass_density_tp(eos, T, P, x, molar_masses, phase="vapor")
rho_mixture = additive_volume_mass_density(mass_fractions, component_densities)

Both results use kg m-3. The additive-volume closure evaluates 1 / rho_mixture = sum_i(w_i / rho_i) and does not normalize its inputs.

from exoeos import SecondVirialEOS, state_tp


eos = SecondVirialEOS(
    jnp.array([[1.0e-4, 2.0e-5], [2.0e-5, 8.0e-5]])  # B_ij, m3 mol-1
)
state = state_tp(eos, T=700.0, P=2.0e5, x=jnp.array([0.4, 0.6]))

state.rho
state.Z
state.lnphi
state.gres_RT

SecondVirialEOS is the first non-ideal fluid model. It uses a constant, symmetric pair-coefficient matrix and

B_mix = sum_i sum_j x_i x_j B_ij,
alphar = rho B_mix,
Z = 1 + rho B_mix,
P = rho R T (1 + rho B_mix),
mu_res_i / (R T) = 2 rho sum_j B_ij x_j,
ln(phi_i) = mu_res_i / (R T) - ln(Z),
g_res / (R T) = 2 rho B_mix - ln(Z).

For state_tp, let rho_0 = P / (R T) and D = 1 + 4 B_mix rho_0. The vapor root is evaluated as rho = 2 rho_0 / (1 + sqrt(D)). Its domain requires T > 0, P > 0, normalized nonnegative x, Z > 0, and D > 0. Only phase="vapor" is supported. The second-virial truncation is a low-density model and should be used only where neglected higher virial terms are small; the constant coefficients also omit real temperature dependence.

PengRobinsonEOS adds a cubic EOS using critical temperatures in K, critical pressures in Pa, acentric factors, and optional binary interaction parameters:

from exoeos import PengRobinsonEOS


eos = PengRobinsonEOS(
    critical_temperatures=jnp.array([190.564]),
    critical_pressures=jnp.array([4_599_200.0]),
    acentric_factors=jnp.array([0.01142]),
)
vapor = state_tp(eos, T=150.0, P=1.0e6, x=jnp.array([1.0]))
liquid = state_tp(
    eos,
    T=150.0,
    P=1.0e6,
    x=jnp.array([1.0]),
    phase="liquid",
)

The implementation uses the original PR76 alpha correlation for every component and the exact critical-condition coefficients Omega_a = 0.45723552892138218938 and Omega_b = 0.077796073903888455972.

PengRobinsonEOS.second_virial_coefficients(T) returns the exact low-density, density-form coefficient matrix of that PR model at T. It can be passed directly to SecondVirialEOS when comparing PR with its consistent second-virial truncation.

A small curated critical-property table is available through get_critical_properties(formula). It contains CO, H2O, CO2, H2, CH4, N2, NH3, H2S, and SO2, with CoolProp source URLs retained in every record. These records can populate the existing PengRobinsonEOS constructor; they supply no fitted mixture interactions or validation at magma temperatures. HCN is not included in the curated table.

The optional binary_interaction_parameters matrix defaults to zero and uses a_ij = (1 - k_ij) sqrt(a_i a_j). phase="vapor" selects the largest physical real root of the Peng-Robinson cubic, while phase="liquid" selects the smallest. Both selectors return the same root in a one-root region. Density derivatives are defined away from multiple roots; exact critical and spinodal states are not differentiable root-selection points.

Use float64 when evaluating very-low-pressure liquid roots: reconstructing a small pressure from dense-liquid Helmholtz terms is ill-conditioned in float32 even when the selected density root is accurate.

Like state_trho, state_tp accepts one scalar state at a time. Use jax.vmap for batches. The phase string is a static selector: capture it in the transformed function or mark it static rather than mapping it.

The reduced fields alphar, mu_res_RT, and gres_RT are dimensionless. psir = A_res / (R T V) has units of mol m-3. This differs from ThermodynamicState.residual_gibbs, which is a molar energy in J mol-1.

ZhangDuanEOS implements the Zhang-Duan (2009) corresponding-states EOS for C-O-H fluids. The published species parameters and the fitted H2O-CO2 and H2O-CH4 interactions are available through from_species:

from exoeos import ZhangDuanEOS


species = ("CO", "H2O", "CO2", "H2")
eos = ZhangDuanEOS.from_species(species)
state = state_tp(
    eos,
    T=1000.0,
    P=1.0e9,
    x=jnp.array([0.4, 0.4, 0.1, 0.1]),
)

The model also supports CH4, O2, and C2H6. It is intended for the homogeneous-fluid calibration ranges reported in Zhang and Duan (2009). The principal mixture range is 673--2573 K and 1 MPa--10 GPa; H2O-CH4 starts at 10 MPa, and the pure-species data ranges differ. Pressure inversion selects the mechanically stable root connected to the low-density branch and supports only phase="vapor". Fugacity coefficients are obtained by differentiating the residual Helmholtz energy rather than by transcribing the paper's mixture fugacity equation. The implementation uses the physical P V / (R T) compressibility; the scaled left-hand side printed in Equation 8 does not reproduce the paper's Table 6 values.

Excess Gibbs API

Solution models implement

gex_RT(T, P, x) = g_ex / (R T),

where P is absolute pressure in Pa. The extensive helper and its amount derivative are

G_ex / (R T) = total_gex_RT(model, T, P, n)
             = n_total gex_RT(T, P, n / n_total),
ln(gamma_i) = partial [G_ex / (R T)] / partial n_i.

import jax
import jax.numpy as jnp

from exoeos import IdealSolution, solution_state, total_gex_RT


model = IdealSolution()
T = 1600.0
P = 1.0e5
x = jnp.array([0.4, 0.6])

gex_RT = model.gex_RT(T, P, x)
total = total_gex_RT(model, T, P, x)
state = solution_state(model, T, P, x)

state.gex_RT
state.lngamma

batched_lngamma = jax.vmap(solution_state, in_axes=(None, 0, 0, 0))(
    model,
    jnp.array([1400.0, 1600.0]),
    jnp.array([1.0e5, 2.0e5]),
    jnp.array([[0.3, 0.7], [0.4, 0.6]]),
).lngamma

IdealSolution is the zero-excess placeholder: gex_RT and lngamma are zero. It does not add the ideal-mixing Gibbs energy. The API uses only the symmetric mole-fraction standard-state convention, a_i = x_i gamma_i, and satisfies gex_RT = sum_i x_i ln(gamma_i) for normalized compositions. The reference is the pure component, or a specified pure endmember, at the same T and P. Standard/endmember Gibbs energies, component-basis mapping, and phase equilibrium remain responsibilities of the calling application.

gex_RT and solution_state accept one state at a time: scalar T, scalar P, and normalized x with shape (K,). total_gex_RT instead accepts a component amount vector and forms x = n / sum(n). The kernel does not clip or numerically validate inputs; the extensive construction uses n / sum(n) by definition. Use jax.vmap for batches.

Fe-Si-O liquid activities

MaFeSiOLiquid implements a native JAX excess-energy model for ordered atomic mole fractions (Fe, Si, O). Its defaults complete Young (2023)'s printed alloy coefficients with the Ma Fe solvent term and use explicit formal endmember standards.

import jax.numpy as jnp
from exoeos import MaFeSiOLiquid, solution_state

model = MaFeSiOLiquid()
T, P = 2350.0, 1.0e5  # K, Pa
x = jnp.array([0.85, 0.10, 0.05])  # Fe, Si, O
model.validate_state(T, P, x)  # Validate eagerly, before JAX transformations.
state = solution_state(model, T, P, x)
shift_RT = model.standard_state_shift_RT(T)
lngamma_source = state.lngamma + shift_RT

The consumer must also transform its source standard potentials: mu0_formal_RT = mu0_source_RT + shift_RT, where potentials are divided by R*T. The model supplies the conversion, not absolute thermochemical data. interaction_K is a differentiable array of shape (3,), ordered (Si-Si, O-O, Si-O), with defaults (12.41*1873, -16500, -5*1873) K. Custom interactions have no established physical calibration.

The activity domain requires positive T, P, and x_Fe, nonnegative fractions, and normalized composition. Call validate_state explicitly; the evaluator does not invoke it. Pure Si/O have zero formal scalar energy, but their activity derivatives are unsupported. No calibrated T/P/composition box or high-pressure validity is established, and liquid stability must be assessed separately. This metal-only model supplies no silicate activities or equilibrium calculation. See the model and reference specification for equations, provenance, limits, and independent fixtures.

Fe-Si-O-H dilution control

MaFeSiOHLiquid extends the dry alloy in atomic (Fe, Si, O, H) order: gex_RT = (1 - x_H) * dry_model.gex_RT(T, P, x[:3] / (1 - x_H)). The four-component ideal activities use x; the dry excess term uses the normalized dry composition. Differentiating this scalar preserves the dry activity coefficients and gives ln(gamma_H) = 0.

import jax.numpy as jnp
from exoeos import MaFeSiOHLiquid, solution_state

model = MaFeSiOHLiquid()
x = jnp.array([0.765, 0.09, 0.045, 0.10])
model.validate_state(2350.0, 1.0e5, x)
state = solution_state(model, 2350.0, 1.0e5, x)
shift_RT = model.standard_state_shift_RT(2350.0)

This is a formal control with no H excess interactions or H partition calibration. It requires positive dry amount and the dry Fe-rich domain; pure H is unsupported. The returned H standard-state shift is zero and preserves the consumer's H standard, whose absolute potential must be supplied separately. See the control specification for its limits and runnable finite-H numerical check.

Independent MELTS silicate references provide three partially crystallized basalt states from a pinned external alphaMELTS release. They record phase masses, liquid endmember activities, chemical potentials, and the explicit oxide-basis conversion. The fixture and its offline consistency tests require no MELTS installation; regenerating it uses a separate runtime. These references add no native silicate model or coupled melt-metal equilibrium calculation.

The optional supplied-composition evaluator under examples/melts_liquid_evaluator.py calls the same pinned backend at requested K, Pa, and liquid endmember amounts. Each evaluation uses a fresh process and checks that the returned composition matches the request with oxygen buffering disabled. It returns full chemical potentials, phase Gibbs energy, standards, activities, and provenance, including explicit conversion to the consumer's mu/(R*T) convention. Its --validate command checks the three saved liquids, Fe/Si/O perturbations, amount scaling, and numerical derivatives. This property evaluator supplies no phase-stability result or alloy/gas standard-state alignment.

An explicit calculation_mode=4 selects a separate rhyolite-MELTS 1.2.0 carbon property path. It retains 19 independent input components and maps the backend's dependent CaCO3 species back to that basis. Positive CO2 requires this mode; SO3 is unsupported and N is absent from the basis. A separate carbon fixture checks finite-difference potentials, amount scaling, and exact zero C on the same model. See the C/N/S provider scope for remaining standard alignment, alloy, dissolution, and calibration work.

Fixed-composition tabulated H/He API

ChabrierDebrasEOS loads the published Chabrier-Debras (2021) H/He TP and T-rho table pair for one fixed composition. The table loader downloads the official data archive when the selected pair is absent, verifies SHA-256 checksums, and caches only the required files. Select one of the published Y0275, Y0292, or Y0297 variants:

from exoeos import ChabrierDebrasTableLoader


loader = ChabrierDebrasTableLoader(
    variant="Y0275",
    # cache_directory="/optional/custom/cache/DirEOS2021",
)
eos = loader.load()
tp_state = eos.state_tp(T=1.0e4, P=1.0e11)
trho_state = eos.state_trho(T=1.0e4, mass_density=1.0e3)

tp_state.rho
tp_state.u
tp_state.s
tp_state.nabla_ad

The default cache is $XDG_CACHE_HOME/exoeos/DirEOS2021, falling back to ~/.cache/exoeos/DirEOS2021. Loader metadata is available through expected_filenames, variant, checksum, checksums, citation, table_domain, and cache_directory. Existing local tables can still be opened directly with ChabrierDebrasEOS.from_directory(...).

Inputs and returned quantities use SI units. The logarithmic derivative fields are dimensionless. This original-table backend interpolates each published field independently; its response columns can differ from automatic derivatives of the interpolated density or entropy, and first derivatives can jump at cell boundaries. The variants are separate fixed-composition datasets; ExoEOS does not interpolate in helium mass fraction. Queries outside the nominal rectangular grids return nan, and the tables do not provide a mask for unphysical states inside those rectangles. If conversion to the selected floating-point dtype makes any returned field non-finite, the complete state is returned as nan. Both evaluators accept one state at a time; use jax.vmap for batches. Y0292 and Y0297 are the effective-abundance variants defined by the authors.

For a separate potential-consistent reconstruction, enable JAX 64-bit mode before loading and call original.to_helmholtz():

import jax
from exoeos import ChabrierDebrasTableLoader

jax.config.update("jax_enable_x64", True)
original = ChabrierDebrasTableLoader(variant="Y0275").load()
potential = original.to_helmholtz()
state = potential.state_trho(T=1.0e4, mass_density=1.0e3)
residuals = original.helmholtz_residuals(potential, 1.0e4, 1.0e3)

This opt-in backend reconstructs a = u - T*s, interpolates a local C2 Helmholtz potential, and derives pressure, entropy and responses from it. It retains mass-specific SI units. residuals exposes signed differences from every original T-rho field, and state.is_stable checks local thermal and mechanical stability. Consistency does not guarantee source accuracy or stability; see the method, API and numerical comparison.

Composition-dependent tabulated silicate-hydrogen API

MarcumSilicateHydrogenEOS interpolates the published Marcum, Stixrude, and Young (2026) MgSiO3-H lookup table at temperature, pressure, and composition. The ordered endmember mole-fraction vector is (MgSiO3, MgSiO3H4), so X = x[1] and the corresponding hydrogen mass fraction is 4 X M_H / (M_MgSiO3 + 4 X M_H).

import jax.numpy as jnp

from exoeos import MarcumSilicateHydrogenTableLoader


eos = MarcumSilicateHydrogenTableLoader().load()
state = eos.state_tp(
    T=6000.0,
    P=1.0e11,
    x=jnp.array([0.5, 0.5]),
)

state.rho
state.h
state.s
state.cp
state.Ks
state.nabla_ad

The loader downloads the single CSV from the authors' table repository, verifies its SHA-256 checksum, and caches it. An existing file can be opened with MarcumSilicateHydrogenEOS.from_file(path). The default cache is $XDG_CACHE_HOME/exoeos/MgSiO3-H-EOS, falling back to ~/.cache/exoeos/MgSiO3-H-EOS. Provenance and domain metadata are exposed by expected_filename, checksum, commit, table_url, citation, and table_domain.

Inputs and outputs use SI units. The table coordinates are 3000--10000 K, 1--800 GPa, and 2.5e-5 <= X <= 1; the 1--4 GPa slices stop at 6000 K. Queries that require a missing cell or lie outside the table return an all-nan state; there is no clipping or extrapolation. The evaluator accepts one scalar state and a normalized, nonnegative composition of shape (2,); it never clips or renormalizes inputs. Use jax.vmap for batches. Its molar_masses and mass_density_tp method also satisfy the MassDensityProvider contract directly.

The table extends beyond the directly simulated calibration range of roughly 4000--8000 K and 4.9--615.83 GPa. Intermediate compositions are the authors' ideal Gibbs mixture of the dry and fully hydrogenated endmembers, not an additive-volume mixture with a separate pure-H2 EOS.

Composite density providers

The density-provider layer combines heterogeneous EOS backends without assigning workflow-specific species mappings to ExoEOS. Its public types are MassDensityProvider, DensityComponent, TPHelmholtzDensityProvider, FixedCompositionDensityProvider, and AdditiveVolumeCompositeDensityProvider. A composite maps an ordered global species tuple to components, converts mole fractions to component mass fractions, and applies the additive-volume law.

Component species must form an exact, non-overlapping partition of the global species tuple. TPHelmholtzDensityProvider adapts a TP Helmholtz EOS, while FixedCompositionDensityProvider adapts a fixed-composition backend such as ChabrierDebrasEOS. If the supplied within-group mass fractions do not match that backend's configured composition within its declared tolerance, the density result is nan.

mass_density_tp(T, P, x) evaluates one state: T and P are scalars in K and Pa, x is a one-dimensional mole-fraction vector in the declared species order, molar masses use kg mol-1, and the result uses kg m-3. Use external jax.vmap for batches. Numerical normalization and positivity are caller contracts. MELTYQ-specific species aliases and EOS assignments, and conversions to or from ExoGibbs or ExoJAX units, remain at the calling workflow or example boundary.

Caloric ideal-gas API

All public quantities use SI units. Component molar masses are in kg mol-1, component molar heat capacities are in J mol-1 K-1, temperature is in K, and pressure is in Pa.

import jax
import jax.numpy as jnp

from exoeos import IdealGas


eos = IdealGas(
    molar_masses=jnp.array([2.01588e-3, 4.002602e-3]),
    molar_heat_capacities=jnp.array([28.84, 20.786]),
)

x = jnp.array([0.85, 0.15])
state = eos.state(T=1000.0, P=1.0e5, x=x)

state.Z
state.mass_density
state.number_density
state.log_fugacity_coefficients
state.residual_gibbs
state.residual_enthalpy
state.cp
state.cv
state.thermal_expansion
state.adiabatic_gradient

jitted_density = jax.jit(
    lambda temperature: eos.state(temperature, 1.0e5, x).mass_density
)
density_gradient = jax.grad(jitted_density)(1000.0)

x contains normalized mole fractions on its last axis. ExoEOS does not renormalize composition. Component molar masses and heat capacities are required because T, P, and x alone do not determine mass density or caloric properties. Reference enthalpies and entropies default to zero; pass physical reference data when absolute values are needed.

The complete units, shape, reference-state, and transformation contract is in the thermodynamic-state contract.

Total Helmholtz thermodynamics

HelmholtzThermodynamics(residual, ideal) combines an existing residual EOS with an ideal free energy and evaluates its first and second T/rho derivatives. It provides molar enthalpy, entropy, internal/Gibbs/Helmholtz energies, cp, cv, sound speed, adiabatic gradient, compressibilities and thermal expansion in one HelmholtzThermodynamicState:

import jax.numpy as jnp
from exoeos import HelmholtzThermodynamics, IdealGas, SecondVirialEOS

ideal = IdealGas(molar_masses=[0.028], molar_heat_capacities=[29.1])
eos = HelmholtzThermodynamics(SecondVirialEOS([[1.0e-5]]), ideal)
state = eos.state_tp(T=500.0, P=1.0e5, x=jnp.array([1.0]))
cp_mass = state.cp / state.mean_molar_mass  # J/(kg K), useful for RCE
sound_speed = state.sound_speed            # m/s
adiabatic_gradient = state.nabla_ad

state_trho(T, rho, x) accepts molar density in mol/m3; state aliases state_tp. The generic thermodynamic_state_trho(model, T, rho, x) also accepts a complete custom MolarHelmholtzEOS implementing molar_helmholtz(T, rho, x) in J/mol and molar_masses in kg/mol. IdealGas supplies a constant-cp ideal closure; custom ideal closures can represent temperature-dependent heat capacities. A residual EOS alone does not specify total caloric properties.

Calls evaluate a single state and support jax.jit and external jax.vmap. Responses hold composition fixed in a homogeneous phase; chemical-equilibrium and latent-heat responses require an additional closure. Values are not clipped at unstable or singular states. For ExoJAX pressure inputs in bar, convert to Pa with P_bar * 1.0e5. See the derivative guide for equations, units and a batched atmospheric example. Existing residual states and their fugacity fields remain available through state_trho and state_tp.

Development

python -m pip install -e ".[test]"
pytest tests/unittests

SecondVirialEOS, PengRobinsonEOS and ZhangDuanEOS are the non-ideal fluid backends. Additional fluid EOS and nonzero Gibbs-excess models can be added behind the separate TPHelmholtzEOS and GibbsExcessModel contracts.