Metadata-Version: 2.4
Name: exoeos
Version: 0.2.0
Summary: Differentiable equations of state and excess free-energy models for planetary atmospheres, fluids, and melts
Author-email: Hajime Kawahara <divrot@gmail.com>
Maintainer-email: Hajime Kawahara <divrot@gmail.com>
License-Expression: GPL-3.0-or-later
Project-URL: Homepage, https://github.com/HajimeKawahara/exoeos
Project-URL: Repository, https://github.com/HajimeKawahara/exoeos
Project-URL: Issues, https://github.com/HajimeKawahara/exoeos/issues
Project-URL: Changelog, https://github.com/HajimeKawahara/exoeos/releases
Keywords: activity coefficients,equations of state,excess Gibbs energy,exoplanets,jax,thermodynamics
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Operating System :: OS Independent
Classifier: Topic :: Scientific/Engineering :: Atmospheric Science
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: jax>=0.4.30
Requires-Dist: numpy>=1.22
Provides-Extra: docs
Requires-Dist: ipykernel>=6; extra == "docs"
Requires-Dist: matplotlib>=3.7; extra == "docs"
Requires-Dist: nbconvert==7.17.1; extra == "docs"
Requires-Dist: nbformat>=5; extra == "docs"
Requires-Dist: pypandoc-binary==1.17; extra == "docs"
Requires-Dist: sphinx>=5; extra == "docs"
Requires-Dist: sphinx-rtd-theme>=1; extra == "docs"
Provides-Extra: reference
Requires-Dist: teqp==0.23.2; extra == "reference"
Provides-Extra: test
Requires-Dist: pytest>=7; extra == "test"
Dynamic: license-file

# 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

```bash
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:

```bash
python -m pip install .
```

## Documentation

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

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

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

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

```bash
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

```text
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.

```python
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:

```python
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`:

```python
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.

```python
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

```text
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:

```python
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`:

```python
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)](https://doi.org/10.1016/j.gca.2009.01.021).
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

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

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

```text
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.
```

```python
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.

```python
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](documents/fe_si_o_reference.rst) 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`.

```python
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](documents/fe_si_o_h_reference.rst)
for its limits and runnable finite-H numerical check.

Independent [MELTS silicate references](documents/melts_silicate_reference.rst)
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](documents/melts_liquid_evaluator.rst)
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](documents/cns_provider_scope.rst)
for remaining standard alignment, alloy, dissolution, and calibration work.

## Fixed-composition tabulated H/He API

`ChabrierDebrasEOS` loads the published
[Chabrier-Debras (2021)](https://doi.org/10.3847/1538-4357/abfc48) H/He TP
and T-rho table pair for one fixed composition. The table loader downloads
the [official data archive](https://perso.ens-lyon.fr/gilles.chabrier/DirEOS/)
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:

```python
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()``:

```python
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](documents/potential_tables.rst).

## Composition-dependent tabulated silicate-hydrogen API

`MarcumSilicateHydrogenEOS` interpolates the published
[Marcum, Stixrude, and Young (2026)](https://arxiv.org/abs/2608.27401)
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)`.

```python
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](https://github.com/s-marcum/MgSiO3-H-EOS),
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.

```python
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](https://github.com/HajimeKawahara/exoeos/blob/main/documents/thermodynamic_state_contract.md).

## 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`:

```python
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](documents/thermodynamic_derivatives.rst) for equations,
units and a batched atmospheric example. Existing residual states and their
fugacity fields remain available through `state_trho` and `state_tp`.

## Development

```bash
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.
