# Ziff-Gulari-Barshad model: Steady State Conditions.¶

Note

To follow this tutorial, either:

Download

`SteadyState.py`

(run as`$AMSBIN/amspython SteadyState.py`

).Download

`SteadyState.ipynb`

(see also: how to install Jupyterlab)

The mechanism being studied frequently involves several intermediary
species, which concentrations (partial pressure for species in the gas
phase and coverage fraction for species that are adsorbed) change in a
particular manner as the simulation progresses. However, after a
particular simulation time, the concentration of all intermediate
species remains constant. At that point, we say the system has reached a
steady state configuration. Sadly, we are unable to anticipate the
simulation time at which that stage will be reached. Therefore, the more
logical approach to take is to: 1) run the simulation for a set period
of time (initial guess), 2) determine if the system has reached a steady
state, and 3) if so, consider the calculation to have converged and
stop; otherwise, increase the simulation time and resume repeating all
steps until convergence. With **pyZacros**, and much more directly with
**Zacros**, this is laborious to complete manually. Thus, this tutorial
aims to demonstrate how to accomplish this goal using the class
`ZacrosSteadyStateJob`

from the **pyZacros** extended components.

The first step is to import all packages we need:

```
import numpy
import scm.pyzacros as pz
import scm.pyzacros.models
```

Then, we initialize the **pyZacros** environment.

```
scm.pyzacros.init()
```

```
PLAMS working folder: /home/user/pyzacros/examples/ZiffGulariBarshad/plams_workdir
```

Notice this command created the directory where all **Zacros** input and
output files will be stored if they are needed for future reference
(`plams_workdir`

by default). Typically, the user doesn’t need to use
these files.

First, we define the physical system to study. Here we used the Ziff-Gulari-Bashard model. This simple model describes the catalytic processes of carbon monoxide oxidation to carbon dioxide (\(\text{CO}+\frac{1}{2}\text{O}_2\longrightarrow \text{CO}_2\)) and accurately captures the interesting property of the phase transition between two surface poisoned states (either CO- or O-poisoned) and a steady state in between. It is named after Robert M. Ziff, Erdogan Gulari, and Yoav Barshad’s pioneering work in 1986. Check the API documentation for more details about its implementation in pyZacros.

```
zgb = pz.models.ZiffGulariBarshad()
```

Then, we set up the Zacros calculation. All parameters are set using a
`Setting`

object. To begin, we define the physical parameters: the
molar fractions of the gas species (`CO`

and `O2`

), the temperature
(in K), and the pressure (in bar). The calculation parameters are then
set: `species numbers`

(in s) determines how frequently information
about the number of gas and surface species will be stored, `max time`

(in s) specifies the maximum allowed simulated time, and “random seed”
specifies the random seed to make the calculation precisely
reproducible. Keep in mind that `max time`

defines the calculation’s
stopping criterion, and it is the parameter we will control below to
achieve the steady-state configuration. Finally, we create the
`ZacrosJob`

, which uses the parameters we just defined as well as the
Ziff-Gulari-Bashard model’s lattice, mechanism, and cluster expansion.
Notice we do not run this job, we use it as a reference for the
steady-state calculation described below.

```
z_sett = pz.Settings()
z_sett.molar_fraction.CO = 0.42
z_sett.molar_fraction.O2 = 1.0 - z_sett.molar_fraction.CO
z_sett.temperature = 500.0
z_sett.pressure = 1.0
z_sett.species_numbers = ('time', 0.1)
z_sett.max_time = 100.0*0.1
z_sett.random_seed = 953129
job = pz.ZacrosJob( settings=z_sett,
lattice=zgb.lattice,
mechanism=zgb.mechanism,
cluster_expansion=zgb.cluster_expansion )
```

It is now time to set up the steady state calculation. It also needs a
`Settings`

object to set its parameters, as shown in the first block
of code below. To begin, we define the parameters for calculating the
turn-over frequency (TOF), which is the property that will be monitored
to determine convergence as the steady state is reached. In a nutshell,
for a given simulation time (that will be increased systematically), the
simulation is divided into an `turnover_frequency.nbatch`

ensemble of
contiguous batches (20 for this case) where the TOFs are calculated for
each one. If the estimated confidence level for these TOFs is higher
than the confidence level specified by the
`turnover frequency.confidence`

parameter (96% for this case),
convergence is then considered to have been achieved. The
`turnover frequency.nreplicas`

parameter allows several simulations to
run in parallel to speed up the calculation at the expense of more
computational power. For the time being, we will leave it at 1, but we
will return to it later.

In the second block of code, the `ZacrosSteadyStateJob.Parameters()`

class allows us to specify the grid in `max time`

, which in this case
ranges from 20 to 1000 every 100 seconds. Take note that the convergence
is verified for each point on this grid, and if it has not converged,
the calculation is resumed up to the next point in `max time`

.

Finally, we create `ZacrosSteadyStateJob`

, which references the
`ZacrosJob`

defined above as well as the `Settings`

object and
parameters we just defined:

```
ss_sett = pz.Settings()
ss_sett.turnover_frequency.nbatch = 20
ss_sett.turnover_frequency.confidence = 0.96
ss_sett.turnover_frequency.nreplicas = 1
parameters = pz.ZacrosSteadyStateJob.Parameters()
parameters.add( 'max_time', 'restart.max_time', numpy.arange(20.0, 1000.0, 100) )
ss_job = pz.ZacrosSteadyStateJob( settings=ss_sett, reference=job, parameters=parameters )
```

```
[27.01|09:11:22] JOB plamsjob Steady State Convergence: Using nbatch=20,confidence=0.96,ignore_nbatch=1,nreplicas=1
```

The steady-state calculation setup is ready. Therefore, we can start it
by invoking the function `run()`

, which will provide access to the
results via the `results`

variable after it has been completed. The
sentence involving the method `ok()`

, verifies that the calculation
was successfully executed, and waits for the completion of every
executed thread in case of parallel execution.

```
results = ss_job.run()
if not ss_job.ok():
print('Something went wrong!')
```

```
[27.01|09:11:22] JOB plamsjob STARTED
[27.01|09:11:22] JOB plamsjob RUNNING
[27.01|09:11:22] JOB plamsjob/ss_iter000 Steady State: NEW
[27.01|09:11:22] JOB plamsjob/ss_iter000 STARTED
[27.01|09:11:22] JOB plamsjob/ss_iter000 RUNNING
[27.01|09:11:22] JOB plamsjob/ss_iter000 FINISHED
[27.01|09:11:23] JOB plamsjob/ss_iter000 SUCCESSFUL
[27.01|09:11:23] species TOF error ratio conv?
[27.01|09:11:23] CO -0.57600 0.08915 0.15478 False
[27.01|09:11:23] O2 -0.29120 0.07190 0.24691 False
[27.01|09:11:23] CO2 0.57667 0.08573 0.14866 False
[27.01|09:11:23] JOB plamsjob Steady State Convergence: NO CONVERGENCE REACHED YET
[27.01|09:11:23] JOB plamsjob/ss_iter001 Steady State: NEW (dep=plamsjob/ss_iter000)
...
[27.01|09:11:44] JOB plamsjob/ss_iter007 FINISHED
[27.01|09:11:45] JOB plamsjob/ss_iter007 SUCCESSFUL
[27.01|09:11:46] species TOF error ratio conv?
[27.01|09:11:46] CO -0.60170 0.02183 0.03628 True
[27.01|09:11:46] O2 -0.30081 0.01095 0.03639 True
[27.01|09:11:46] CO2 0.60170 0.02183 0.03628 True
[27.01|09:11:46] JOB plamsjob Steady State Convergence: CONVERGENCE REACHED. DONE!
[27.01|09:11:46] JOB plamsjob FINISHED
[27.01|09:11:47] JOB plamsjob SUCCESSFUL
```

If the execution got up to this point, everything worked as expected. Hooray!

Now, in the following lines, we just nicely print the results in a
table. See the API documentation to learn more about how the `results`

object is structured. Here we show the history of the simulation and see
how it progresses as the `max_time`

is increased. We print the TOF for
CO2 (in mol/s/site), its error, and whether the calculation converged.
Notice that the calculation should have been converged at 720 s of
`max_time`

.

```
print(60*'-')
fline = "{0:>8s}{1:>10s}{2:>15s}{3:>12s}{4:>10s}"
print( fline.format('iter', 'max_time', 'TOF_CO2', 'error', 'conv?') )
print(60*'-')
for i,step in enumerate(results.history()):
fline = "{0:8d}{1:>10.2f}{2:15.5f}{3:>12.5f}{4:>10s}"
print( fline.format(i,
step['max_time'],
step['turnover_frequency']['CO2'],
step['turnover_frequency_error']['CO2'],
str(all(step['converged'].values()))) )
```

```
------------------------------------------------------------
iter max_time TOF_CO2 error conv?
------------------------------------------------------------
0 20.00 0.57667 0.08573 False
1 120.00 0.59301 0.03157 False
2 220.00 0.59360 0.03311 False
3 320.00 0.59293 0.03062 False
4 420.00 0.59345 0.02729 False
5 520.00 0.61241 0.02564 False
6 620.00 0.60345 0.02697 False
7 720.00 0.60170 0.02183 True
```

Now that all calculations are done, we can close the pyZacros environment:

```
scm.pyzacros.finish()
```

```
[27.01|09:11:47] PLAMS run finished. Goodbye
```

Additionally, you can see the aforementioned results visually if you have installed the package matplotlib. Please review the code below. In particular, pay close attention to how to obtain the “children results” (each thread executed) which are identified by their replica and iteration number. In the figure, each iteration is represented by a different color, and the stopping points on “max time” are represented by vertical gray dashed lines.

```
import matplotlib.pyplot as plt
fig = plt.figure()
ax = plt.axes()
ax.set_xlabel('Time (s)', fontsize=14)
ax.set_ylabel("CO$_2$ Production (mol/site)", fontsize=14)
colors = 'bgrcmykb'
for i in range(results.niterations()):
for j in range(results.nreplicas()):
molecule_numbers = results.children_results(i,j).molecule_numbers(['CO2'], normalize_per_site=True)
ax.plot( molecule_numbers['Time'], molecule_numbers['CO2'], lw=3, color=colors[i], zorder=-i )
ax.vlines( max(molecule_numbers['Time']) , 0, max(molecule_numbers['CO2']), colors='0.8', linestyles='--',)
plt.show()
```

Finally, if you run the entire script, replacing `ss_sett.nreplicas = 1`

with `ss_sett.nreplicas = 4`

, you should get the following result:

The calculation now converged in two iterations rather than the previous eight. When using replicas, each replica uses the same parameters as the reference job but with different random seeds, and the corresponding TOFs are evaluated as an average over the entire replica set. This accelerates convergence by increasing the accessible space’s sampling efficiency.