prometheus_equilibrium.equilibrium.element_matrix

ElementMatrix — stoichiometric matrix relating species to elements.

The element matrix A is the core linear-algebra object shared by all three equilibrium solvers. Entry A[i, k] gives the number of atoms of element k in species i (from the species’ elements dict).

Mathematical role

The equilibrium constraints are:

A^T · n = b₀

where n is the vector of species mole amounts and b₀ is the vector of total element abundances (fixed by the reactant composition). The element residual

Δb = b₀ − A^T · n

appears in the right-hand side of the element-balance equations in all three solvers.

Basis selection (Browne’s method — used by PEP and Hybrid solvers)

A basis is a subset of S species (S = number of elements) whose S×S composition sub-matrix B is non-singular. The optimised basis (Browne, 1960) picks the S species with the largest mole amounts that satisfy this condition, found by sorting species descending by mole amount and accepting candidates via the Gram-Schmidt orthogonality test.

The reaction-coefficient matrix ν = C · B⁻¹ expresses every non-basis species as a stoichiometric combination of basis species. Together B and ν are the foundation of the Villars reaction-adjustment iteration and the Hybrid solver’s compressed Newton system.

Charge balance

The electron pseudo-element "e-" is included as an ordinary column so that ionic species are balanced automatically.

Singular-system handling (Gordon-McBride solver)

If the reactants do not supply every element represented in the product species list, the corresponding column of A has no matching b₀ entry and the Jacobian will be singular. independent_elements() identifies the linearly independent subset of elements via QR decomposition, allowing the G-McB solver to drop redundant constraints (RP-1311 §3.2).

Classes

ElementMatrix(species, elements)

Stoichiometric matrix A mapping species to element abundances.

class prometheus_equilibrium.equilibrium.element_matrix.ElementMatrix(species: List[Species], elements: List[str])

Bases: object

Stoichiometric matrix A mapping species to element abundances.

Parameters

specieslist of Species

All species in the mixture (gas first, condensed second).

elementslist of str

Ordered list of element symbols that define the columns of A. Typically derived from the union of all elements present in species.

classmethod from_mixture(mixture: Mixture) → ElementMatrix

Build from a Mixture.

The element list is the union of all elements across all species, sorted alphabetically (with "e-" placed last if present).

property matrix: ndarray

The raw stoichiometric matrix A, shape (n_species, n_elements).

property elements: List[str]

Ordered element symbol list (defines the column ordering of A).

property species: List[Species]

Species list (defines the row ordering of A).

property n_species: int

property n_elements: int

element_abundances(moles: ndarray) → ndarray

Compute element abundances b = A^T · n [mol of element k].

Parameters

molesarray-like, shape (n_species,)

Species mole amounts nⱼ.

Returns

np.ndarray, shape (n_elements,)

Total moles of each element in the mixture.

element_residuals(moles: ndarray, b0: ndarray) → ndarray

Compute the element-balance residual Δb = b₀ − A^T · n.

This vector must be zero at equilibrium. It forms the right-hand side of the element-balance equations in all three solvers.

Parameters:

• moles – Current species mole amounts, shape (n_species,).

• b0 – Target element abundances fixed by the reactant composition, shape (n_elements,).

Returns:

Residual vector b₀ − A^T·n, shape (n_elements,).

select_basis(moles: ndarray) → Tuple[List[int], List[int]]

Select the optimised basis using Browne’s method (Gram-Schmidt).

The optimised basis is the set of S = n_elements species with the largest mole amounts whose S×S composition sub-matrix B is non-singular. It is found by sorting species in descending mole- amount order and accepting each candidate via the Gram-Schmidt orthogonality test (Cruise, NWC TP 6037, eq. 5).

At the start of each outer iteration the basis is re-optimised (since the composition has changed). This automatic re-selection serves two purposes:

1. It keeps the basis well-conditioned (dominant species are always in the basis, so B⁻¹ is numerically stable).

2. It relieves the user from choosing a basis manually.

Parameters

molesarray-like, shape (n_species,)

Current mole amounts nⱼ. Species with nⱼ ≤ 0 are skipped.

Returns

basis_indiceslist of int

Indices (into self.species) of the S basis species, in the order they were accepted.

nonbasis_indiceslist of int

All remaining species indices, in original species order.

Raises

RuntimeError

If fewer than S linearly independent species are found (the system is under-determined — reactant elements are not all representable).

Notes

The Gram-Schmidt test replaces the more expensive determinant test. A candidate row cₘ is orthogonalised against the current incomplete basis rows; if the result is the zero vector, the candidate is linearly dependent and rejected (Cruise eq. 5). The PEP paper also describes updating the basis via a linear-programming tableau pivot (Smith & Missen, 1968) to avoid recomputing B⁻¹ from scratch each iteration; this optimisation can be added once the basic loop works.

basis_matrix(basis_indices: List[int]) → ndarray

Return B = A[basis_indices, :], the S×S basis composition matrix.

B[j, k] = atoms of element k in the j-th basis species.

Parameters

basis_indiceslist of int, length S

Row indices of the basis species in A.

Returns

np.ndarray, shape (S, S) where S = n_elements

reaction_coefficients(basis_indices: List[int]) → ndarray

Compute ν = C · B⁻¹, the reaction-coefficient matrix.

ν[i, j] is the stoichiometric coefficient of basis species j consumed when one mole of (non-basis) species i is formed from the basis. The chemical reaction for species i is:

Σⱼ ν[i,j] · basis_species[j]  →  species[i]

Used to:

• Compute equilibrium constants (PEP solver): ln Kᵢ = −Σⱼ ν[i,j] · gⱼ°/RT

• Correct basis mole amounts after each stoichiometric adjustment (PEP eq. 10).

Parameters

basis_indices : list of int, length S

Returns

np.ndarray, shape (n_species, n_elements)

Full ν matrix (rows for all species, including basis species themselves — their rows will be identity-like by construction).

rank() → int

Numerical rank of A (number of independent element constraints).

Computed via singular-value decomposition. A rank less than n_elements indicates redundant elements that will make the G-McB Jacobian singular if included naively.

independent_elements() → List[str]

Return the subset of element names forming a full-rank system.

Uses QR decomposition with column pivoting on Aᵀ to identify a maximal linearly independent set of element columns. Used by the G-McB solver to drop redundant element rows from the Jacobian (RP-1311 §3.2).

Returns

list of str

Element symbols in a linearly independent ordering (most influential elements first, as determined by QR pivoting).

reduced(active_elements: List[str] | None = None) → ElementMatrix

Return a new ElementMatrix restricted to active_elements.

If active_elements is None, independent_elements() is called automatically. Used by the G-McB solver to remove redundant constraints before building the Newton system.

Parameters

active_elementslist of str, optional

Element symbols to keep. Must be a subset of self.elements.

Returns

ElementMatrix

A new instance with the same species but only the specified element columns.

gas_rows() → ndarray

Sub-matrix of A for gas-phase species only, shape (n_gas, n_elements).

condensed_rows() → ndarray

Sub-matrix of A for condensed-phase species, shape (n_cnd, n_elements).