API Documentation

NICE (N-species Ice Table) module – provides simultaneous equilibria solvers:

• simulate Net-event kinetic Monte Carlo solver
• solve Nonlinear equations solver (Newton trust region method)

The ReactionSystem type is provided as an interface for these solvers.

Definition of the simultaneous equilibria reaction system. This type is the input to the solver functions simulate and solve, and is updated in-place.

The reaction system consists of

• $M$ Reactions involving $N$ species;
• a matrix $C = \left(c_{mn}\right)$ defining the stoichiometry of the system; rows correspond to species and columns correspond to reactions; each element is the stoichoimetric coefficient of the $n$th species in the $m$th reaction, and is positive if the species is a product, and negative if it is a reactant;
• the forward and reverse rate constants $k_f$ and $k_r$ for each reaction;
• initial concentrations or activities $\mathbf{a} = \left(a_n\right)$ of each species;

and can be written (through some abuse of notation) as

\[\begin{aligned} \sum_{\left. n \middle| c_{1n} < 0 \right.} &\left|c_{1n}\right| I_n &\xrightleftharpoons[k_r^{\left(1\right)}]{k_f^{\left(1\right)}} \sum_{\left. n \middle| c_{1n} > 0 \right.} &\left|c_{1n}\right| I_n \\ \vdots & ~ & ~ & \vdots \\ \sum_{\left. n \middle| c_{mn} < 0 \right.} &\left|c_{mn}\right| I_n &\xrightleftharpoons[k_r^{\left(M\right)}]{k_f^{\left(M\right)}} \sum_{\left. n \middle| c_{Mn} > 0 \right.} &\left|c_{Mn}\right| I_n \end{aligned}\]

ReactionSystem(stoich, concs, rev_rate_consts, fwd_rate_consts)

Instantiate a ReactionSystem instance using the forward and reverse rate constants for each reaction.

ReactionSystem(stoich, concs, K_eqs; φ)

Instantiate a ReactionSystem instance using the equilibrium constants $K_\text{eq}^{\left(m\right)}$ to derive the forward and reverse rate constants for each reaction:

\[\begin{aligned} k_r &= \frac{\phi}{K_\text{eq} - 1} \\ k_f &= \phi - k_r \end{aligned}\]

hybrid_solve(rxn_system, K_eqs; n_iter=Int(1e+8), chunk_iter=Int(1e+4), ε=1.0e-4, ε_scale=1.0, ε_concs=0.0, tol_ε=0.0, maxiters=1000, abstol=1.0e-9, reltol=0.0)

Combines stochastic simulation and deterministic solving to find equilibrium concentrations. First, it runs a Net-Event Kinetic Monte Carlo simulation via simulate to approach equilibrium, then refines the result using the deterministic solve method with provided $K_\text{eq}$ values. The function updates rxn_system.concs and returns the final nonlinear solution object.

jacobian_function!(J, rxn_system, K_eqs, xi)

Computes the Jacobian matrix J of the nonlinear system in-place.

This implements the analytical Jacobian for the system of equations defined by fr = K^(r) * ∏(ai^(-νi^(r))) - ∏(ai^(νi^(r))) = 0 where the Jacobian element Jrs = ∂fr/∂ξs is computed using:

Jrs = K^(r) * ∑{i=1}^N [∏{j=1,j≠i}^N (aj^(0) + ∑{r=1}^R νj^(r) ξ^(r))^(-νj^(r))] * [-νi^(r) * (ai^(0) + ∑{r=1}^R νi^(r) ξ^(r))^(-νi^(r)-1)] * (νi^(s)) - ∑{i=1}^N [∏{j=1,j≠i}^N (aj^(0) + ∑{r=1}^R νj^(r) ξ^(r))^(νj^(r))] * [νi^(r) * (ai^(0) + ∑{r=1}^R νi^(r) ξ^(r))^(νi^(r)-1)] * (ν_i^(s))

Arguments:

• J::Matrix{Float64}: The output Jacobian matrix to be filled in-place.
• rxn_system::ReactionSystem: The reaction system containing stoichiometry and initial concentrations.
• K_eqs::AbstractVector{Float64}: The equilibrium constants for each reaction, length R.
• xi::Vector{Float64}: The reaction progress vector ξ, length R.

objective_function(f, rxn_system, K_eqs, xi)

Computes the objective function fr(ξ) for each reaction r using the equation: fr = K^(r) * ∏(ai^(-νi^(r))) - ∏(ai^(νi^(r))) = 0 where the first product is over reactants (νi^(r) < 0) and the second over products (νi^(r) > 0). Returns a vector of length R, where f_r = 0 at equilibrium.

Arguments:

• f::Vector{Float64}: The output vector to store the computed objective function values.
• rxn_system::ReactionSystem: The reaction system with stoichiometry, initial concentrations, and equilibrium constants.
• K_eqs::AbstractVector{Float64}: The equilibrium constants for each reaction, length R.
• xi::AbstractVector{T}: The reaction progress vector ξ, length R, where T is a subtype of Real (e.g., Float64 or Dual).

simulate(rxn_system;
         n_iter=Int(1e+8), chunk_iter=Int(1e+4),
         ε=1.0e-4, ε_scale=1.0, ε_concs=0.0,
         tol_t=Inf, tol_ε=0.0, tol_concs=0.0)

Run a Net-Event Kinetic Monte Carlo (NEKMC) simulation to find the equilibrium concentrations of the reaction system.

solve(rxn_system::ReactionSystem; maxiters, abstol, reltol)

Solves for the equilibrium state of the given ReactionSystem.

This function sets up and solves the system of nonlinear equations for the reaction extents (ξ) that bring all reactions to equilibrium. It uses the NonlinearSolve.jl package with an analytically provided Jacobian for efficiency and robustness.

Upon successful convergence, the concs field of the rxn_system object is updated with the final equilibrium concentrations.

Arguments

• rxn_system::ReactionSystem: The reaction system to solve. All necessary parameters (stoich, concs_init, keq_vals) are taken from this object.

Keyword Arguments

• maxiters::Integer=100: Maximum number of iterations for the solver.
• abstol::Real=1.0e-9: Absolute tolerance for the solver.
• reltol::Real=0.0: Relative tolerance for the solver.

Returns

• solution: The full solution object from NonlinearSolve.solve, which contains the final reaction extents (solution.u), the return code (solution.retcode), and other diagnostic information.

solve(rxn_system, K_eqs; max, φ=1.0, rate_consts=:forward, maxiters=1000, abstol=1.0e-9, reltol=0.0)

Solve the system deterministically.

This method solves a linear least squares problem to get an initial guess for $\mathbf{\xi}$,

\[{\left(c_{mn}\right)}^T \left(\xi_m^{\left(0\right)}\right) = \left(a_n - a_n^{\left(0\right)}\right)\]

and then optimizes the equation for the simultaneous equilibria,

\[K_\text{eq}^{\left(m\right)} = \prod_{n=1}^N {\left( a_n^{\left(0\right)} + \sum_{m=1}^M c_{mn} \xi_m \right)}^{c_{mn}}\]

using Newton's method + Trust Region. Instead of specifying the $K_\text{eq}$ values directly, they are approximated here using the forward or reverse rate constants and a parameter $\phi$. If :forward is chosen, equilibrium constants are derived as $K_\text{eq} = (k_f - 2\phi)/(k_f - \phi)$; if :reverse, as $K_\text{eq} = (k_r + \phi)/k_r$. The function updates rxn_system.concs upon successful convergence and returns a solution object.