README.rst

NICE

(Net-Event) Kinetic Monte Carlo and exact solvers for simultaneous equilibria.

Dependencies

The following programs/libraries are required to build NICE:

• GFortran (or another Fortran compiler, e.g. Flang, PGF, Intel)
• NumPy (including C headers)

The following libraries are required to run NICE:

• NumPy
• SciPy
• CMA-ES (optional: to use "cma" mode with exact solver)

The following programs/libraries are required to build the NICE documentation:

• Sphinx
• Read the Docs Sphinx Theme

Installation

To download NICE via git:

git clone https://github.com/quantumelephant/nice.git
cd nice

To build NICE in-place for use in cloned folder:

python setup.py build_ext --inplace

To install NICE for a user (i.e., locally):

python setup.py install --user

Building Documentation

To install the Read the Docs Sphinx theme via pip:

pip install sphinx-rtd-theme --user

Then, after installing NICE, run the following to build the HTML documentation:

cd doc && make html

Example

We have a set of simultaneous equilibrium equations

1 A        <-->  2 B    K_eq = 1.0
1 A + 2 B  <-->  2 C    K_eq = 0.1

and initial concentrations

[A] = 1.0 M,  [B] = 0.2 M,  [C] = 0.4 M

We set up the NEKMC solver:

import numpy as np
import nice

concs  = np.array([1.0, 0.2, 0.4])
stoich = np.array([[-0.5, 1.0, 0.0], [-0.5, -1.0, 1.0]])
keqs   = np.array([1.0, 0.1])
nekmc  = nice.NEKMCSolver(concs, stoich, keq_values=keqs)

Each row of stoich represents a reaction. Each column is the stoichiometric coefficient of a species; the value is negative for a reactant, positive for a product, and zero if the species does not participate. Be careful, as the equilibrium constants depend on the stoichiometry you select. I.e.,

1 A + 2 B <--> 2 C K_eq = 0.1

0.5 A + B <--> C K_eq = (0.1)^{1/2} = 0.32

It is also possible to initialize the solver by passing the forward and reverse rate constants directly:

concs  = np.array([1.0, 0.2, 0.4])
stoich = np.array([[-0.5, 1.0, 0.0], [-0.5, -1.0, 1.0]])
consts = np.array([[1.0, 0.1],  # forward rate constants
                   [0.1, 1.0]]) # reverse rate constants
nekmc  = nice.NEKMCSolver(concs, stoich, rate_consts=consts)

We run the solver with concentration step-size matching our desired precision for 50,000 iterations:

nekmc.iterate(step=1e-6, niter=50000)

Finally, we set up and run the exact solver, with initial guess generated from the NEKMC solver, to get a more precise result:

exact = nice.ExactSolver(concs, stoich, keq_values=keqs)
exact.optimize(guess=nekmc.compute_zeta(), tol=1e-9)

Now our system is fully converged to equilibrium. Final concentrations for each solver are located at nekmc.concs and exact.concs. See the code documentation for more detailed information on running NICE.