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

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.