Results Analysis

Upon completion of a kinetic simulation run by solve_network (or through the end of a CRN exploration using explore_network, which also calls this function), Kinetica returns an ODESolveOutput object. This binds together the CRN that this simulation was run on (after modifications such as those described in Filtering CRNs and Removing Low-Rate Reactions), the parameters the simulation was run with, and the simulation's results. This data can be used for analysing the results of a kinetic simulation.

We will demonstrate this analysis using the results of the CRN we generated and simulated in Getting Started:

using Kinetica, Plots
res = load_output("../my_CRN_out/direct_network_final.bson")

Analysing the CRN

The output CRN is stored across two data structures: a SpeciesData for holding information about the species within the CRN, and a RxData for holding information about the reactions between those species. These can be accessed through ODESolveOutput.sd and ODESolveOutput.rd respectively. For more information on these objects, see CRN Representation.

sd, rd = res.sd, res.rd

By default, Kinetica stores information about these species and reactions in a compact form to maintain efficiency when working with large CRNs. For example, reactions are stored as Vectors of integer species IDs. To inspect individual reactions, we can use the print_rxn function:

# Print the first reaction in the CRN.
print_rxn(sd, rd, 1)

C --> [CH2] + [H][H]

The total number of reactions in a CRN is stored in RxData.nr, so printing all of the reactions in a CRN is as simple as this loop:

for i in 1:rd.nr
    print_rxn(sd, rd, i)
end

C --> [CH2] + [H][H]
[CH2] + [H][H] --> [CH3] + [H]
[CH2] + [H][H] --> C
[CH3] + [H] --> [CH2] + [H][H]
C --> [CH3] + [H]
[CH3] + [H] --> C
[H][H] --> 2 [H]
2 [H] --> [H][H]
2 [CH2] --> [CH]C
[CH]C + [H][H] --> CC
CC --> [CH2]C + [H]
[CH]C --> 2 [CH2]
CC --> [CH]C + [H][H]
[CH2]C + [H] --> CC
C + [CH3] --> CC + [H]
CC + [H] --> C + [CH3]
[CH2] + [CH3] --> [CH2]C
[CH2]C + [H] --> C=C + [H][H]
C=C --> [CH]=C + [H]
[CH2]C --> [CH2] + [CH3]
C=C + [H][H] --> [CH2]C + [H]
[CH]=C + [H] --> C=C
[CH3] --> [CH2] + [H]
[CH2] + [H] --> [CH3]
C + [CH2] --> [CH2]C + [H]
[CH2]C + [H] --> C + [CH2]
C + [CH2] --> CC
[CH2]C --> [CH]=C + [H][H]
CC --> C + [CH2]
[CH]=C + [H][H] --> [CH2]C

Other fields of RxData can also be useful for analysis, such as RxData.mapped_rxns. This field contains the atom-mapped reaction SMILES of all reactions in a CRN, and can be used for constructing multi-species reactant and product geometries with consistent atom indeces (vital for any reaction path techniques):

println("Atom-mapped reaction SMILES: $(rd.mapped_rxns[1])\n")

# Split reaction SMILES into individual SMILES for reactants/products.
am_reacs, am_prods = String.(split(rd.mapped_rxns[1], ">>"))

# Generate non-overlapping systems of molecules for reactants and products.
reac_species = [sd.toStr[sid] for sid in rd.id_reacs[1]]
reacsys = system_from_smiles(reac_species; dmin=3.0)
prod_species = [sd.toStr[sid] for sid in rd.id_prods[1]]
prodsys = system_from_smiles(prod_species; dmin=3.0)

# Remap atom indeces using atom-mapped SMILES and RDKit substructure matching.
mapped_reacs = atom_map_frame(am_reacs, reacsys)
mapped_prods = atom_map_frame(am_prods, prodsys)

# Add extra info to frames, convert to ExtXYZ and print.
mapped_reacs["info"]["type"] = "reactant"
mapped_reacs["info"]["SMILES"] = join(sort(reac_species), ".")
mapped_prods["info"]["type"] = "product"
mapped_prods["info"]["SMILES"] = join(sort(prod_species), ".")
println("Reaction XYZ:")
print(frame_to_xyz(mapped_reacs))
print(frame_to_xyz(mapped_prods))

Atom-mapped reaction SMILES: [C:1]([H:2])([H:3])([H:4])[H:5]>>[C:1]([H:4])[H:5].[H:2][H:3]

Reaction XYZ:
5
SMILES=C type=reactant comment=
C 0.9447905649425166 0.866894886438913 0.3396095028483221
H 0.5807205649425164 1.598654886438913 -0.384870497151678
H 0.5807205649425164 -0.12640511356108697 0.06812950284832209
H 2.0369905649425166 0.866894886438913 0.3396095028483221
H 0.5807305649425165 1.1284348864389129 1.3355795028483222
5
SMILES=[CH2].[H][H] type=product comment=
C -0.6586838128363594 9.912063656304642 -2.5841004766158817
H 20.337829224506123 -3.675939487784442 27.089117847582727
H 21.04582922450612 -3.675939487784442 27.089117847582727
H 0.21085618716364063 10.439513656304642 -2.900230476615882
H -1.5282838128363596 9.384703656304643 -2.268030476615882

As a result of the atom mapping, the hydrogen atoms which dissociate from the reactant ethane are now consistently indexed between reactant and product geometries.

Species-Reaction Graphs

Kinetica hooks in to the same graph plotting code as Catalyst.jl, which uses Graphviz to generate figures that show the relationship between species and reactions in a CRN. Graph(::SpeciesData, ::RxData) has been extended to allow for easy graph plotting directly from a Kinetica CRN:

g = Graph(sd, rd)

Graphs generated this way are compatible with Catalyst's graph methods and can be saved to file using Catalyst.savegraph.

In addition, Kinetica allows for passing Graphviz attributes directly through to the plotter for extra customisability. These take the form of the graph_attrs, species_attrs, rxn_attrs and edge_attrs keyword arguments to Graph(::SpeciesData, ::RxData), which can be supplied with Dict{Symbol, String}s to modify parameters of the overall graph, the species nodes and the reaction nodes, and the edges respectively:

g = Graph(sd, rd;
    graph_attrs=Dict(
        :layout => "sfdp",
        :overlap => "prism",
        :overlap_scaling => "-8"
    ), species_attrs=Dict(
        :shape => "hexagon",
        :color => "aquamarine1"
    ), rxn_attrs=Dict(
        :shape => "box",
        :color => "coral"
    ), edge_attrs=Dict(
        :color => "grey59"
    ))

The exploration level that each species and reaction was found in (see Iterative CRN Exploration) is additionally added as a Graphviz attribute to each respective node. While this is (probably) not useful within the context of Graphviz's layout engines, more advanced graphical interfaces such as Gephi can be used to create level-driven layouts:

Species Analysis

Kinetica makes use of many functions within both RDKit and Open Babel to assist with basic property prediction and conversion between geometry and SMILES. A comprehensive list of functions implemented are available in the Open Babel and RDKit API pages.

Kinetica supports cached calculation of properties such as species molecular weights and hard-sphere radii through Open Babel through the get_species_stats! function. This places these properties in the SpeciesData.cache, which is useful for storing per-species values during a calculation.

Kinetica also stores species geometries within SpeciesData.xyz. These are stored as ExtXYZ.jl frames and can be manipulated as such.

Analysing the Kinetic Simulation

Most of the information regarding the results of a kinetic simulation are stored in ODESolveOutput.sol, which is a DifferentialEquations ODESolution (in most cases, some extensions may apply depending on the simulation method but this can always be interacted as if it is an ODESolution). As such, the DiffEq documentation applies for directly analysing simulation results.

However, when plotting simulation results, it can often be useful to have information about the experimental conditions that the simulation took place under, as well as some of the other parameters used in the simulation which are not included within the ODESolution. This includes the names of each of the species, which are only referred to be an integer ID within Kinetica's solution methods. Kinetica therefore includes Plots.jl plot recipes which act on the ODESolveOutput object, as this has access to all of these values alongside the ODESolution. The plots created by these recipes can be modified in the same way as any other Plots.jl figure.

Concentration-Time Plots

This is the 'default' plot type for ODESolveOutput, with the call signature:

plot(::ODESolveOutput; tunit="s", label_above=0.1, ignore_species=nothing, ignore_below=nothing)

where

• tunit is the unit of time to display on the x-axis;
• label_above is a concentration above which species should contribute to the plot's legend and be plotted in a colour (species is grey otherwise);
• ignore_species can be a Vector of SMILES strings representing species which should not be plotted (nothing if this is not desired);
• ignore_below can be a concentration threshold where species with a maximum concentration below are not plotted (nothing if this is not desired).

For example, plotting the same results as in Getting Started, but excluding the radical species from the plot is simple:

is_radical(smi) = ('[' in smi) && !(smi == "[H][H]")
radical_species = [spec for spec in [sd.toStr[i] for i in 1:sd.n] if is_radical(spec)]

plot(res; label_above=0.01, ignore_species=radical_species, linewidth=5)

Conditions Plots

Kinetica defines a separate plot recipe for variable condition profiles, conditionsplot:

conditionsplot(::Union{::ODESolveOutput, ConditionSet}, ::Symbol; tunit="s")

This recipe can take either an ODESolveOutput or a ConditionSet, but it has a second mandatory argument - a Symbol representing the condition profile to plot. For example, if a variable temperature profile is required (and it exists in the ConditionSet being plotted), then this is usually available under :T. Again, tunit is the unit of time to display on the x-axis.

# Again, this is compatible with Plots.jl keyword arguments.
conditionsplot(res, :T; linecolor=:red)

Final Concentration Plots

It can often be useful to examine the final species concentrations at the end of a kinetic simulation, as these often correspond to experimental observables. Kinetica implements a specialised bar chart recipe for this, finalconcplot:

finalconcplot(::ODESolveOutput; 
              quantity=:conc, n_top=10, highlight_radicals=false, 
              ignore_species=nothing, xscale=:identity)

where

• quantity is a Symbol, either :conc or :percent. The former plots species concentrations, while the latter plots percentage concentrations with respect to the overall reaction mixture;
• n_top controls the number of species to include in the plot. These are sorted in decreasing order, with species outside n_top going into their own 'Others' bar at the end. For example, with n_top=10, the 10 species with the highest final concentrations will be plotted, followed by a bar for the combined final concentration of all other species;
• highlight_radicals colours radical species in red when true, as these are often not desired at the end of simulations;
• ignore_species accepts a Vector of SMILES strings for species to exclude from the plot. This can be useful when inert species are present with concentrations that are unimportant to the final state of the CRN (nothing if this is not desired);
• xscale mimics Plots.jl's xscale option for setting the scaling of the x-axis. In this case, only :identity (default) and :log10 are supported, as additional calculations are required to obtain the correct x-axis limits for correct bar chart rendering.

finalconcplot(res; quantity=:percent, n_top=5, highlight_radicals=true)