Post-processing Zacros Output

Part 3: Post-processing and visualising event statistics

During the course of the KMC simulation, Zacros aggregates the statistics of the occurrence of elementary events (adsorption/desorption, diffusion, different types of reactions) into procstat_output.txt. Post-processing and visualisation of this file can be done with Zacros-post, by selecting the appropriate option from the main window's menu, as shown below (notice also in the Windows Explorer the extracted adlayer configuration file state_input_ExtrConf4.dat which we generated in the previous part of this Tutorial).

After selecting this menu option, you should see a window like the one in the following snapshot. By default, all events in the KMC simulation are listed, even ones that were defined in the mechanism_input.dat file, but did not occur (i.e. were not sampled by the KMC algorithm during this simulation). Zacros-post groups the forward and reverse steps of reversible events, and for each of these events it calculates the occurrence frequency of the forward and the reverse step (i.e. the number of events occurred divided by the duration of the simulation). These are plotted as the bars labelled "fwd" and "rev", respectively, in a bar graph like the one below. In addition, Zacros-post calculates and plots the net rates of the reversible events, as the forward rate minus the reverse rate. If the net rate is positive, a "net-fwd" bar is plotted, otherwise a "net-rev" bar is plotted. Thus, while for a reversible event you may see both "fwd"- and "rev"-labelled bars, you will always see only one of the "net-fwd" and "net-rev" bars.

Moreover, notice the dashed vertical line around the occurrence frequency of 0.5: this line corresponds to the occurrence of a single event in the entire duration of the simulation. If you open general_output.txt and scroll down to the bottom (after the "Simulation stopped" message in that file), you will see that the simulation covered approximately 2 KMC time units; thus, the frequency of 1 event would be about 1/2 = 0.5. This frequency, with the corresponding dashed vertical line in the bar graph, provides important information about the "threshold of accuracy" of your simulation: clearly you should never see a bar below that line, since you cannot have "a fraction" of 1 event happening during your KMC simulation. Moreover, bars close to this threshold correspond to frequencies of events that were estimated with only a few event occurrences, so they are expected to be subject to relatively large errors. For instance, if, for one of the events in our simulation, Zacros had only sampled two events, we would have seen a bar at occurrence frequency of around 1, and this would be quite an inaccurate estimate: what if in a different realisation of the KMC simulation we had 3 events, or just 1 event? The relative error of this frequency of 2±1 events would be quite large indeed.

Certain post-processing and styling options can be set using the Edit menu, which is shown above. Thus, custom axes limits can be set from Edit > Axis limits; this works in a similar way as we have discussed for the lattice structure plot, so we won't present is again here. Moreover, the style of the plot can be customised from the Edit > Plot style option as shown above. The controls of the pertinent popup window are quite intuitive, and below you can see a customised plot in which we are plotting only the net rates, we have set the x-axis to display things in linear scale, and we have decreased the font size of the event name labels. We have also added some white margin on the right of the plot, using the second to last button of the Matplotlib toolbar at the bottom of the window (whose icon looks like three sliders).

Plotting the frequencies in linear scale can be useful when trying to draw conclusions about the active pathways and check whether the simulation has reached stationarity (roughly speaking, the equivalent of steady-state behaviour for stochastic systems). For instance, below, we see that the forward rates of CO adsorption, H₂O adsorption and H₂O decomposition, are almost equal to each other, and also equal to the reverse rate of H₂ adsorption (thus, the forward rate of H₂ desorption). This is good news, as it is in line with the stoichiometry of the water-gas shift reaction (CO + H₂O → CO₂ + H₂). What about the rate of CO2 production? Note that the overal rate at which CO₂ is elicited from the surface is the sum of the rates of CO oxidation (CO* + O* → CO₂) and CO-OH disproportionation (CO* + OH* → CO₂ + H*). In this simulation, the former event is represented as the reverse of the CO2 dissociative adsorption event. Indeed, you can easily see that the length of e.g. the H₂ dissociative adsorption bar, is approximately equal to the sum of the lengths of the following two bars: the light purple bar for CO₂ dissociative adsorption plus the length of the bar for the CO-OH disproportionation reaction.

In the preceding discussion, you may have noticed that the CO adsorption rate seems actually slightly larger than e.g. the H₂ dissociateive adsorption event. This is because in our sampling we are including initial configurations for which the system has not yet reached stationarity. This can be easily remedied in Zacros-post: notice the tollbar with the slider and the two input boxes labelled "Skip configurations" and "Use sliding horizon" (the latter is actually a check-box, which enables the correspnoding option). By changing the number of skipped configurations to say 25, and zooming into the plot (see below) we see that the frequency of occurrence of CO adsorption is equal to those of the other events (of course, there is always some randomness in the KMC simulation which results to deviations from "perfect equality"). Note also, in the plot below, that the time interval in the title of the plot has been updated, with the low bound being now 0.25.

In a similar way, one can use the "sliding horizon" option to restrict the processing of the event statistics to a "moving window" of ΔN configurations (entries in the procstat_output.txt file). Below you can see an example in which we process configurations N₁ = 130 to N₂ = N₁ + ΔN = 140, the latter being the one noted by the slider (we also skip the first 25 configurations, this does not affect the numbering). For this example simulation the sampling of the event statistics is done at regular time intervals of Δt = 0.01 units of time (see simulation_input.dat), and therefore these configurations correspond to time points (N₁ − 1)*Δt = 1.29 and (N₂ − 1)*Δt = 1.39, as noted in the title of the plot. Compare this plot with the corresponding one we showed first in this section; you will notice that the plots look quite similar (taking into account the stochasticity of the KMC simulation), which is positive news, as it means we are sampling the system in the stationary regime. However, there is a key difference between these two plots: in the one below, the dashed horizontal line denoting the frequency of the "single event" (also, the threshold of accuracy of our simulation), has now moved to higher values. This is expected: previously we were considering the aggregate statistics over an interval of about 2 KMC time units. Now, we are considering a much smaller interval, or merely 0.1 time units. Thus, the threshold of accuracy now is much higher: we cannot "measure" frequencies lower than 1 event per unit time.

There are two more features that we will discuss before closing this section. From the menu, by choosing Edit > Events displayed, you can easily switch on or off any of the events displayed in the bar graph, e.g. you may choose to only plot reaction events (no adsorption/desorption events), as shown below. In addition, by choosing Edit > Normalization, you can normalise the event rates with respect to the number of sites. You can choose one or more site types (if applicable); in the latter case, the sum of the number of sites for the chosen types will be used, and the occurrence frequencies displayed in the bar graph, will be calculated as: number of events occurred divided by the time interval, divided by the number of sites specified.

Once you are done with the events statistics' post-processing you can close the window and your custom settings will be saved into postprocessing.xml, as discussed for the the other apps.

In the next and final section of this tutorial, we will look into the postprocessing of the species numbers and associated information written in the specnum_output.txt file.