Generating an Input File¶

The first step for any kinetic model is to properly generate an input file. Input files are processed by a “parser” class of CatMAP. The job of the parser is to convert some standard input into a Python data structure compatible with the kinetic model. The default parser is the TableParser which has been designed to accept inputs in a tabular format. This tutorial is designed to explain how users can create their own input files for any model which uses the TableParser (currently the only “parser” class implemented). If users are interested in creating inputs in some other format then it is also possible to design custom “parser” classes, but this is advanced and will not be discussed further here.

The tutorial is broken into two parts: an overview of input file structure and an example of how to create an input file for methane synthesis on \({\rm{Rh}}(111)\) from DFT.

surface_name    site_name   species_name  formation_energy    bulk_structure  frequencies reference
None            gas         CH4             0                   None            []          By Definition
None            gas         H2O             0                   None            []          By Definition
None            gas         H2              0                   None            [4401]      By Definition
None            gas         CO              2.74                None            [2170]      Energy Environ. Sci., 3, 1311-1315 (2010)
Pt              211         CO              1.113               fcc             []          J. Phys. Chem. C, 113 (24), 10548-10553 (2009)

abinitio_energies = {
         'CO_gas': -626.611970497,
         'H2_gas': -32.9625308725,
         'CH4_gas': -231.60983421,
         'H2O_gas': -496.411394229,
         'CO_111': -115390.445596,
         'C_111': -114926.212205,
         'O_111': -115225.106527,
         'H_111': -114779.038569,
         'CH_111': -114943.455431,
         'OH_111': -115241.861661,
         'CH2_111': -114959.776961,
         'CH3_111': -114976.7397,
         'C-O_111': -115386.76440668429,
         'H-OH_111': -115257.78796158083,
         'H-C_111': -114942.25042955727,
         'slab_111': -114762.254842,
         }

ref_dict = {}
ref_dict['H'] = 0.5*abinitio_energies['H2_gas']
ref_dict['O'] = abinitio_energies['H2O_gas'] - 2*ref_dict['H']
ref_dict['C'] = abinitio_energies['CH4_gas'] - 4*ref_dict['H']
ref_dict['111'] = abinitio_energies['slab_111']

from ase.atoms import string2symbols

def get_formation_energies(energy_dict,ref_dict):
    formation_energies = {}
    for key in energy_dict.keys(): #iterate through keys
        E0 = energy_dict[key] #raw energy
        name,site = key.split('_') #split key into name/site
        if 'slab' not in name: #do not include empty site energy (0)
            if site == '111':
                E0 -= ref_dict[site] #subtract slab energy if adsorbed
            #remove - from transition-states
            formula = name.replace('-','')
            #get the composition as a list of atomic species
            composition = string2symbols(formula)
            #for each atomic species, subtract off the reference energy
            for atom in composition:
                E0 -= ref_dict[atom]
            #round to 3 decimals since this is the accuracy of DFT
            E0 = round(E0,3)
            formation_energies[key] = E0
    return formation_energies

formation_energies = get_formation_energies(abinitio_energies,ref_dict)

for key in formation_energies:
    print key, formation_energies[key]

>>
>> OH_111 0.323
>> H_111 -0.302
>> C_111 1.727
>> H2O_gas 0.0
>> CH_111 0.965
>> CO_111 0.943
>> H2_gas 0.0
>> C-O_111 4.624
>> CO_gas 2.522
>> O_111 0.597
>> CH3_111 0.644
>> CH4_gas 0.0
>> CH2_111 1.125
>> H-OH_111 0.878
>> H-C_111 2.17
>>

frequency_dict = {
                'CO_gas': [2170],
                'H2_gas': [4401],
                'CH4_gas':[2917,1534,1534,3019,3019,3019,1306,
                           1306,1306],
                'H2O_gas': [3657, 1595, 3756],
                'CO_111': [60.8, 230.9, 256.0, 302.9, 469.9, 1747.3],
                'C_111': [464.9, 490.0, 535.9],
                'O_111': [359.5, 393.3, 507.0],
                'H_111': [462.8, 715.9, 982.5],
                'CH_111': [413.3, 437.5, 487.6, 709.6, 735.1, 3045.0],
                'OH_111': [55, 340.9, 396.1, 670.3, 718.0, 3681.7],
                'CH2_111': [55, 305.5, 381.3, 468.0, 663.4, 790.2, 1356.1,
                            2737.7, 3003.9],
                'CH3_111': [55, 113.5, 167.4, 621.8, 686.0, 702.5, 1381.3,
                            1417.5, 1575.8, 3026.6, 3093.2, 3098.9],
                'C-O_111': [],
                'H-OH_111': [],
                'H-C_111': []
                }

def make_input_file(file_name,energy_dict,frequency_dict):

    #create a header
    header = '\t'.join(['surface_name','site_name',
                        'species_name','formation_energy',
                        'frequencies','reference'])

    lines = [] #list of lines in the output
    for key in energy_dict.keys(): #iterate through keys
        E = energy_dict[key] #raw energy
        name,site = key.split('_') #split key into name/site
        if 'slab' not in name: #do not include empty site energy (0)
            frequency = frequency_dict[key]
            if site == 'gas':
                surface = None
            else:
                surface = 'Rh'
            outline = [surface,site,name,E,frequency,'Input File Tutorial.']
            line = '\t'.join([str(w) for w in outline])
            lines.append(line)

    lines.sort() #The file is easier to read if sorted (optional)
    lines = [header] + lines #add header to top
    input_file = '\n'.join(lines) #Join the lines with a line break

    input = open(file_name,'w') #open the file name in write mode
    input.write(input_file) #write the text
    input.close() #close the file

    print 'Successfully created input file'

file_name = 'energies.txt'
make_input_file(file_name,formation_energies,frequency_dict)

>> Successfully created input file

surface_name    site_name   species_name    formation_energy    frequencies reference
None    gas CH4 0.0 [2917, 1534, 1534, 3019, 3019, 3019, 1306, 1306, 1306]  Input File Tutorial.
None    gas CO  2.522   [2170]  Input File Tutorial.
None    gas H2  0.0 [4401]  Input File Tutorial.
None    gas H2O 0.0 [3657, 1595, 3756]  Input File Tutorial.
Rh  111 C   1.727   [464.9, 490.0, 535.9]   Input File Tutorial.
Rh  111 C-O 4.624   []  Input File Tutorial.
Rh  111 CH  0.965   [413.3, 437.5, 487.6, 709.6, 735.1, 3045.0] Input File Tutorial.
Rh  111 CH2 1.125   [55, 305.5, 381.3, 468.0, 663.4, 790.2, 1356.1, 2737.7, 3003.9] Input File Tutorial.
Rh  111 CH3 0.644   [55, 113.5, 167.4, 621.8, 686.0, 702.5, 1381.3, 1417.5, 1575.8, 3026.6, 3093.2, 3098.9] Input File Tutorial.
Rh  111 CO  0.943   [60.8, 230.9, 256.0, 302.9, 469.9, 1747.3]  Input File Tutorial.
Rh  111 H   -0.302  [462.8, 715.9, 982.5]   Input File Tutorial.
Rh  111 H-C 2.17    []  Input File Tutorial.
Rh  111 H-OH    0.878   []  Input File Tutorial.
Rh  111 O   0.597   [359.5, 393.3, 507.0]   Input File Tutorial.
Rh  111 OH  0.323   [55, 340.9, 396.1, 670.3, 718.0, 3681.7]    Input File Tutorial.

#Test that input is parsed correctly
from catmap.model import ReactionModel
from catmap.parsers import TableParser
rxm = ReactionModel()
#The following lines are normally assigned by the setup_file
#and are thus not usually necessary.
rxm.surface_names = ['Rh']
rxm.adsorbate_names = ['CO','C','O','H','CH','OH','CH2','CH3']
rxm.transition_state_names = ['C-O','H-OH','H-C']
rxm.gas_names = ['CO_g','H2_g','CH4_g','H2O_g']
rxm.species_definitions = {'s':{'site_names':['111']}}
#Now we initialize a parser instance (also normally done by setup_file)
parser = TableParser(rxm)
parser.input_file = file_name
parser.parse()
#All structured data is stored in species_definitions; thus we can
#check that the parsing was successful by ensuring that all the
#data in the input file was collected in this dictionary.
for key in rxm.species_definitions:
    print key, rxm.species_definitions[key]

Creating a Microkinetic Model¶

This tutorial provides a walk-through of how to create a very basic micro-kinetic model for CO oxidation on transition-metal (111) surfaces. More advanced features are explored in Refining a Microkinetic Model and the Topics section.

surface_name    site_name   species_name    formation_energy    bulk_structure  frequencies other_parameters    reference
None    gas CO2 2.45    None    [1333,2349,667,667] []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
None    gas CO  2.74    None    [2170]  []  "Energy Environ. Sci., 3, 1311-1315 (2010)"
None    gas O2  5.42    None    [1580]  []  "Falsig et al (2012)"
Ru  111 O   -0.07   fcc []  []  "Falsig et al (2012)"
Ni  111 O   0.35    fcc []  []  "Falsig et al (2012)"
Rh  111 O   0.55    fcc []  []  "Falsig et al (2012)"
Cu  111 O   1.07    fcc []  []  "Falsig et al (2012)"
Pd  111 O   1.55    fcc []  []  "Falsig et al (2012)"
Pt  111 O   1.62    fcc []  []  "Falsig et al (2012)"
Ag  111 O   2.05    fcc []  []  "Falsig et al (2012)"
Au  111 O   2.61    fcc []  []  "Falsig et al (2012)"
Ru  111 CO  1.3 fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Rh  111 CO  1.34    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Pd  111 CO  1.55    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Ni  111 CO  1.63    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Pt  111 CO  1.7 fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Cu  111 CO  2.58    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Ag  111 CO  2.99    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Au  111 CO  3.04    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Ru  111 O-CO    2.53    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Rh  111 O-CO    3.1 fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Ni  111 O-CO    3.25    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Pt  111 O-CO    4.04    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Cu  111 O-CO    4.18    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Pd  111 O-CO    4.2 fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Ag  111 O-CO    5.05    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Au  111 O-CO    5.74    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Ag  111 O-O 5.98    fcc []  []  "Falsig et al (2012)"
Au  111 O-O 7.22    fcc []  []  "Falsig et al (2012)"
Cu  111 O-O 4.74    fcc []  []  "Falsig et al (2012)"
Pt  111 O-O 5.35    fcc []  []  "Falsig et al (2012)"
Rh  111 O-O 3.79    fcc []  []  "Falsig et al (2012)"
Ru  111 O-O 3.34    fcc []  []  "Falsig et al (2012)"
Pd  111 O-O 5.34    fcc []  []  "Falsig et al (2012)"

rxn_expressions = [

               '*_s + CO_g -> CO*',
               '2*_s + O2_g <-> O-O* + *_s -> 2O*',
               'CO* +  O* <-> O-CO* + * -> CO2_g + 2*',

                   ]

surface_names = ['Pt', 'Ag', 'Cu','Rh','Pd','Au','Ru','Ni']
#surfaces to include in scaling (need to have descriptors defined for each)

descriptor_names= ['O_s','CO_s'] #descriptor names

descriptor_ranges = [[-1,3],[-0.5,4]]

resolution = 15

temperature = 500

species_definitions = {}
species_definitions['CO_g'] = {'pressure':1.} #define the gas pressures
species_definitions['O2_g'] = {'pressure':1./3.}
species_definitions['CO2_g'] = {'pressure':0}

species_definitions['s'] = {'site_names': ['111'], 'total':1} #define the sites

data_file = 'CO_oxidation.pkl'

input_file = 'energies.txt'

gas_thermo_mode = "shomate_gas"
adsorbate_thermo_mode = "frozen_adsorbate"

scaling_constraint_dict = {
                           'O_s':['+',0,None],
                           'CO_s':[0,'+',None],
                           'O-CO_s':'initial_state',
                           'O-O_s':'final_state',
                           }

decimal_precision = 100 #precision of numbers involved

tolerance = 1e-50 #all d_theta/d_t's must be less than this at the solution

max_rootfinding_iterations = 100

max_bisections = 3

from catmap import ReactionModel

mkm_file = 'CO_oxidation.mkm'
model = ReactionModel(setup_file=mkm_file)
model.run()

>> mapper_iteration_0: status - 100 points do not have valid solution.
>> minresid_iteration_0: success - [ 3.00, 4.00] using coverages from [ 3.00, 4.00]
>> minresid_iteration_0: success - [ 3.00, 3.50] using coverages from [ 3.00, 3.50]
>> ...
>> ...
>> ...
>> minresid_iteration_0: success - [-1.00, 0.00] using coverages from [-1.00, 0.00]
>> minresid_iteration_0: success - [-1.00,-0.50] using coverages from [-1.00,-0.50]
>> mapper_iteration_1: status - 0 points do not have valid solution.

>> initial_evaluation: success - initial guess at point [-1.00,-0.50]
>> initial_evaluation: success - initial guess at point [-1.00, 0.00]
>> initial_evaluation: success - initial guess at point [-1.00, 0.50]
>> ...

Analyzing the Output¶

Accessing Output¶

If you look in the working directory you should see 5 files:

• energies.txt (input file)
• CO_oxidation.mkm (setup file)
• mkm_job.py (submission script)
• CO_oxidation.log (log file)
• CO_oxidation.pkl (data file)

The log file and the data file contain all information about the solved model. The log file is human-readable. If you open it up you will notice that is is actually a python script which contains many of the same things as are found in ‘CO_oxidation.mkm’, but also contains a number of new variable definitions. You will also see that it automatically reads in ‘CO_oxidation.pkl’ and stores the variables from this pickle file in the local namespace. Thus, the “data file” is actually just an extension of the log file which is stored in binary form (this saves a lot of time since the data is often so large). There are two interesting things you can do with this log file:

Load it in as a setup_file to a ReactionModel¶

Make a new file called “test.py” and enter the lines:

from catmap import ReactionModel

model = ReactionModel(setup_file='CO_oxidation.log')

print model.rxn_expressions
print model.coefficient_matrix

Notice that the rxn_expressions are identical to those from the setup file, but that the coefficient_matrix also exists even though we did not define it in the setup file. The coefficient_matrix was created by the scaler during the process of solving the model. The variable “model” in test.py is actually equivalent to the variable “model” in mkm_job.py right after the line with “model.run()”. This is a useful way to load in a model which is already solved for future analysis.

View output in interactive python mode¶

The file can be opened and viewed interactively by entering:

python -i CO_oxidation.log

in the command line. You will now have an interactive python prompt where you can view the various outputs and attributes of the solved model. For example we can look at the coverages or rates as a function of descriptor space:

>>> coverage_map
[[0.7777777777777777, 2.0], [mpf('1.553678172737489e-14'), mpf('0.99999999999788455')]], [[0.33333333333333326, 1.0], [mpf('0.75379752729405923'), mpf('0.246202464223187')]], [[1.2222222222222223, 3.5], [mpf('9.4245829753741903e-28'), mpf('0.99999999982984852')]], [[0.7777777777777777, 0.0], [mpf('1.0'), mpf('4.0785327108804121e-19')]], [[-0.11111111111111116, 3.5], [mpf('3.6157703851402e-41'), mpf('1.0')]], [[0.33333333333333326, 0.5], .... ]
>>> coverage_map[0]
[[0.7777777777777777, 2.0], [mpf('1.553678172737489e-14'), mpf('0.99999999999788455')]]
>>> rate_map[0]
[[0.7777777777777777, 2.0], [mpf('3.0626957315361884e-11'), mpf('1.5313478657680942e-11'), mpf('3.0626957315361884e-11')]]

The format of the “rate_map” and “coverage_map” is a list of lists where the first entry of each list is the point in descriptor space and the second is the rate/coverage. This is not particularly useful if you don’t know what each number in the output corresponds to. You can find out by checking the “output_labels” dictionary:

>>> output_labels['coverage']
('CO_s','O_s')
>>> output_labels['rate']
([['s', 'CO_g'], ['CO_s']], [['s', 's', 'O2_g'], ['O-O_s', 's'], ['O_s', 'O_s']], [['CO_s', 'O_s'], ['O-CO_s', 's'], ['CO2_g', 's', 's']])

In this case the model only outputs the rate and coverage. Information on how to get more outputs can be found in Refining a Microkinetic Model.

Visualizing Output¶

Unless you possess extraordinary skills in raw data visualization then reading the raw output probably doesn’t do you much good. Of course it is possible to use the raw data and write your own plotting scripts, but some tools exist within the micro-kinetics module to get a quick look at the outputs. We will explore some of these tools here.

Rate “Volcano” and Coverage Plots¶

Often the most interesting result from such an analysis is the so-called “volcano” plot of the reaction rate as a function of descriptor space. We can achieve this with the VectorMap plotting class (the “Vector” here refers to the fact that the rates are output as a 1-dimensional list/vector). First we instantiate the plotter using the model of interest by adding the following lines in mkm_job.py after model.run():

from catmap import analyze
vm = analyze.VectorMap(model)

Next we need to give the plotter some information on what to plot and how to plot it:

vm.plot_variable = 'rate' #tell the model which output to plot
vm.log_scale = True #rates should be plotted on a log-scale
vm.min = 1e-25 #minimum rate to plot
vm.max = 1e3 #maximum rate to plot

Most of these attributes are self-explanatory. Finally we create the plot:

vm.plot(save='rate.pdf') #draw the plot and save it as "rate.pdf"

The “save” keyword tells the plotter where to save the plot. You can set “save=False” in order to not save the plot. The plot() function returns the matplotlib figure object which can be further modified if necessary. If we run this script with “python mkm_job.py” we get the following plot:

This looks pretty similar to previously published results by Falsig et. al., with minor differences to be expected since the model and inputs used here are slightly different.

We notice that the rates are given for CO adsorption and oxygen adsorption, but that associative CO2 desorption is not included. This is because it is identical to the plot for CO adsorption (due to the steady-state condition). If we want to include it we can do:

vm.unique_only = False
vm.plot('all_rates.pdf')
vm.unique_only = True

(we turn it back to unique_only right afterwards since this is generally less cluttered)

which gives us a plot for each elementary step:

We might also be interested in the production rate of CO2 rather than the rates of elementary steps (it is trivial to see that they are equivalent here, but this is not always the case). If we want to analyze this we need to include the “production_rate” in the outputs, re-run the model, and re-plot.

model.output_variables += ['production_rate']
model.run()
vm.production_rate_map = model.production_rate_map #attach map
vm.threshold = 1e-30 #do not plot rates below this
vm.plot_variable = 'production_rate'
vm.plot(save='production_rate.pdf')

In the line commented “attach map” we point the VectorMap instance to the new output from the model. This line is not necessary if the output had been included in the original “output_variables”. We also note that the “threshold” variable will be discussed in the next tutorial.

Now we can see whats going on, but its not very pretty (the colorbar is cutoff). We can make a few aesthetic improvements fairly simply:

vm.descriptor_labels = ['CO reactivity [eV]', 'O reactivity [eV]']
vm.subplots_adjust_kwargs = {'left':0.2,'right':0.8,'bottom':0.15}
vm.plot(save='pretty_production_rate.pdf')

Ok, so its still not publishable, but its better. There are ways to control the finer details of the plots, but that will come in a later tutorial.

One more thing we might be interested in is the coverages of various intermediates. This can also be plotted with the VectorMap (since coverages are output as a 1-dimensional “vector”). However, we are going to want to make a few changes to the settings:

vm.plot_variable = 'coverage'
vm.log_scale = False
vm.min = 0
vm.max = 1
vm.plot(save='coverage.pdf')

Not the prettiest plot ever, but you get the point. We could re-adjust the subplots_adjust_kwargs to make this more readable, but that is left as an independent exercise.

Finally, we might not always be interested in seeing all of the coverages. If we only wanted to see the CO coverage we could specify this by:

vm.include_labels = ['CO_s']
vm.plot(save='CO_coverage.pdf')

Note that the strings to use in “include_labels” can be found by examining the “output_labels” dictionary from the log file; alternatively you can specify “include_indices = [0,1,...]” where the integers correspond to the indices of the plots to include.

Refining a Microkinetic Model¶

In this tutorial we will take the simple CO oxidation model presented in Creating a Microkinetic Model and refine it to be more complete. This tutorial should show some of the more powerful capabilities of CatMAP, and highlight the ability to dynamically make changes to the kinetic model with minimal programming. The tutorial will show several possibilities of ways to refine the model towards something that correctly represents the physical system:

• Adding elementary steps (refining the mechanism)
• Adding multiple sites (refining the active site structure)
• Sensitivity analyses (refining the inputs to the model)
• Refining numerical accuracy (resolution, tolerance, etc.)

These sections do not need to be followed sequentially. For each one we will start with the same setup file and input file from Creating a Microkinetic Model. We will use a slightly more efficient submission script:

from catmap import ReactionModel

mkm_file = 'CO_oxidation.mkm'
model = ReactionModel(setup_file=mkm_file)
model.output_variables += ['production_rate']
model.run()

from catmap import analyze
vm = analyze.VectorMap(model)
vm.plot_variable = 'production_rate' #tell the model which output to plot
vm.log_scale = True #rates should be plotted on a log-scale
vm.min = 1e-25 #minimum rate to plot
vm.max = 1e3 #maximum rate to plot

vm.descriptor_labels = ['CO reactivity [eV]', 'O reactivity [eV]']
vm.threshold = 1e-25 #anything below this is considered to be 0
vm.subplots_adjust_kwargs = {'left':0.2,'right':0.8,'bottom':0.15}
vm.plot(save='pretty_production_rate.pdf')

Adding Elementary Steps¶

One way of refining a model is to include additional steps in the mechanism. To demonstrate this we will add molecular adsorption of oxygen prior to dissociation. In order to do this we need to include the relevant energetics, so add the following lines to the energies.txt:

Ru  111 O2  3.15    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Rh  111 O2  3.63    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Ni  111 O2  3.76    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Pd  111 O2  4.29    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Cu  111 O2  4.52    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Pt  111 O2  4.56    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Ag  111 O2  5.1 fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"

Next, we just need to define the new elementary step in the setup file (CO_oxidation.mkm):

rxn_expressions = [

                '*_s + CO_g -> CO*',
#               '2*_s + O2_g <-> O-O* + *_s -> 2O*',
                '*_s + O2_g -> O2_s',
                '*_s + O2_s <-> O-O* + *_s -> 2O*',
                'CO* +  O* <-> O-CO* + * -> CO2_g + 2*',

                    ]

Now we can run mkm_job.py to get the output. If you run in a clean directory you should see something like:

mapper_iteration_0: status - 225 points do not have valid solution.
minresid_iteration_0: success - [ 3.00, 4.00] using coverages from [ 3.00, 4.00]
minresid_iteration_0: success - [ 3.00, 3.68] using coverages from [ 3.00, 3.68]
minresid_iteration_0: success - [ 3.00, 3.36] using coverages from [ 3.00, 3.36]
minresid_iteration_0: success - [ 3.00, 3.04] using coverages from [ 3.00, 3.04]
minresid_iteration_0: success - [ 3.00, 2.71] using coverages from [ 3.00, 2.71]
...
minresid_iteration_1: success - [-1.00, 0.14] using coverages from [-1.00, 0.46]
rootfinding_iteration_2: fail - stagnated or diverging (residual = 3.85907297979e-29)
minresid_iteration_0: fail - [-1.00,-0.18] using coverages from [-1.00,-0.18]; initial residual was 8.73508143601e-21 (residual = 3.85907297979e-29)
minresid_iteration_1: success - [-1.00,-0.18] using coverages from [-1.00, 0.14]
minresid_iteration_0: success - [-1.00,-0.50] using coverages from [-1.00,-0.50]
mapper_iteration_1: status - 0 points do not have valid solution.

However, if you run in the same directory that you used for Creating a Microkinetic Model, you will see slightly different output. Either way, the model should converge.

If you look at “pretty_production_rate.pdf” it should look like the following:

If we compare this to the figure from the previous tutorial we can see that there are a few differences, but the general conclusions are unchanged. If we wanted to be thorough we could continue refining the model by adding more elementary steps (CO2 molecular adsorption, O-O-CO transition state, etc.). For brevity these extensions are omitted.

Using previous results as initial guesses¶

If you ran mkm_job.py in the same directory as you had the CO_oxidation.pkl data file from Creating a Microkinetic Model, you might have noticed that instead of getting output about “minresid_iterations” you get something like:

Length of guess coverage vectors are shorter than the number of adsorbates. Assuming undefined coverages are 0
initial_evaluation: success - initial guess at point [ 3.00, 4.00]
Length of guess coverage vectors are shorter than the number of adsorbates. Assuming undefined coverages are 0
initial_evaluation: success - initial guess at point [ 3.00, 3.68]
Length of guess coverage vectors are shorter than the number of adsorbates. Assuming undefined coverages are 0
initial_evaluation: success - initial guess at point [ 3.00, 3.36]
...

This happens because the model detects the data file (CO_oxidation.pkl) and loads in the coverages to use as an initial guess. However, it notices that there is now more adsorbates than there are coverages since we added O2*. In order to make the best of this, it just assumes that the additional coverages are 0 and uses that as an initial guess. As you can see, it works out okay here. One thing worth noting, however, is that since the code does not know what the order of adsorbates in the previous model was, it cannot properly assign the coverages. Adsorbates are parsed in the order they appear in rxn_expressions, so in this model the order is:

adsorbate_names = ['CO_s','O2_s','O_s']

but, before adding the new elementary step the order was of course different ([‘CO_s’,’O_s’]). Since there are so few adsorbates here it turned out to be a decent initial guess that the coverage of O2* was equal to the coverage of O* from the previous model, and that the coverage of O* was 0. In general, this will not be the case. If you want to use initial guesses from previous models it is best to add the new elementary steps after the old ones. Then the new adsorbates will be assumed to have 0 coverage at the initial guess, rather than scrambling all the coverages around. This is one of the best strategies for obtaining convergence in very complex kinetic models: start with a simple version of the system and slowly add more elementary steps, converging the model along the way and using coverages from the simpler model as an initial guess to the more complex one.

More examples of how to add elementary steps are given in the following section.

Adding multiple sites¶

Structure dependence is a common phenomenon is catalysis, so it is important to use the correct active site structure in order to obtain accurate kinetics. Here we will look at both the (111) and (211) facets for CO oxidation using the previously defined model.

The first thing we will need to do is include the energetic inputs for (211) sites:

Ir  211 CO  0.673   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Re  211 CO  0.753   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Ru  211 CO  0.983   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Rh  211 CO  1.073   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Pt  211 CO  1.113   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Pd  211 CO  1.223   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Ni  211 CO  1.253   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Co  211 CO  1.403   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Fe  211 CO  1.413   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Cu  211 CO  2.283   fcc []   []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Au  211 CO  2.573   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Ag  211 CO  2.873   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Ru  211 O-CO    2.351   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Rh  211 O-CO    2.559   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Co  211 O-CO    2.732   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Ni  211 O-CO    2.768   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Pt  211 O-CO    3.528   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Cu  211 O-CO    3.918   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Pd  211 O-CO    3.992   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Ag  211 O-CO    5.099   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Au  211 O-CO    5.448   fcc []  []  "J. Phys. Chem. C, 113 (24), 10548-10553 (2009)"
Ag  211 O-O 5.34    fcc []  []  Falsig et al (2012)
Au  211 O-O 6.18    fcc []  []  Falsig et al (2012)
Pt  211 O-O 4.9 fcc []  []  Falsig et al (2012)
Pd  211 O-O 4.6 fcc []  []  Falsig et al (2012)
Re  211 O   -1.5    fcc []  []  Falsig et al (2012)
Co  211 O   -0.15   fcc []  []  Falsig et al (2012)
Ru  211 O   -0.1    fcc []  []  Falsig et al (2012)
Ni  211 O   0.18    fcc []  []  Falsig et al (2012)
Rh  211 O   0.28    fcc []  []  Falsig et al (2012)
Cu  211 O   0.93    fcc [] []  Falsig et al (2012)
Pt  211 O   1.32    fcc []  []  Falsig et al (2012)
Pd  211 O   1.58    fcc []  []  Falsig et al (2012)
Ag  211 O   2.11    fcc []  []  Falsig et al (2012)
Au  211 O   2.61    fcc []  []  Falsig et al (2012)
Fe  211 O   -0.73   fcc []  []  "Phys. Rev. Lett. 99, 016105 (2007)"
Ir  211 O   -0.04   fcc []  []  "Phys. Rev. Lett. 99, 016105 (2007)"

We note that there is no data readily available for molecular O2 adsorption on the (211) facet, so we need to make sure we move back to the simpler model from Creating a Microkinetic Model for the (211) analysis:

rxn_expressions = [

               '*_s + CO_g -> CO*',
               '2*_s + O2_g <-> O-O* + *_s -> 2O*',
#               '*_s + O2_g -> O2_s',
#               '*_s + O2_s <-> O-O* + *_s -> 2O*',
               'CO* +  O* <-> O-CO* + * -> CO2_g + 2*',

                   ]

If we check the “pretty_production_rate.pdf” then we see the following:

which is not very pretty. The plot on the right is showing up because the plotter says it is not empty; however, as you can see it looks pretty empty. This is happening because of numerical issues - there are some very small (<1e-50 - the tolerance) positive values for production of CO at some points in descriptor space. The quick way to get rid of this is to set a “threshold” for the plotter, so that it counts very small values as 0:

vm.descriptor_labels = ['CO reactivity [eV]', 'O reactivity [eV]']
vm.threshold = 1e-25
vm.subplots_adjust_kwargs = {'left':0.2,'right':0.8,'bottom':0.15}
vm.plot(save='pretty_production_rate.pdf')

Now we get the following:

The same thing can also be achieved by tightening the numerical precision/tolerance, as discussed later. When we look at the plot we see the leg going out towards Ni/Ru/Rh which, based on the previous section, we can predict will be reduced if molecular oxygen adsorption is considered. We also notice that the maximum is moved towards the nobler metals, which is roughly consistent with the findings of Falsig et. al. who show that nobler metals are more active when undercoordinated clusters are examined.

Of course in a real catalyst, there will be both (111) and (211) facets (along with lots of others, but lets focus on these two for now). We can use CatMAP to examine both facets simultaneously by adding new sites. First, we need to define the mechanisms on both sites:

rxn_expressions = [

               '\*_s + CO_g -> CO*',
               '2*_s + O2_g <-> O-O* + \*_s -> 2O*',
#               '\*_s + O2_g -> O2_s',
#               '\*_s + O2_s <-> O-O* + \*_s -> 2O*',
               'CO* +  O* <-> O-CO* + * -> CO2_g + 2*',

               '\*_t + CO_g -> CO_t',
#               '2*_t + O2_g <-> O-O* + \*_t -> 2O*',
               '\*_t + O2_g -> O2_t',
               '\*_t + O2_t <-> O-O_t + \*_t -> 2O_t',
               'CO_t +  O_t <-> O-CO_t + \*_t -> CO2_g + 2*_t',

               '\*_t + CO_s -> CO_t + \*_s',
               '\*_t + O_s -> O_t + \*_s',


                   ]

Here we use _s (or just * which is equivalent to _s) to denote step sites, and _t to denote terrace sites. We have included molecular oxygen adsorption on the terrace, but not the step since we don’t have the energetics. Diffusion between the step and terrace sites are also included, and they have no activation barrier which implies that there should be equilibrium between CO* and O* on the step/terrace. In addition to the new elementary steps, we also need to include this new “terrace site” in the species definitions:

species_definitions['s'] = {'site_names': ['211'], 'total':0.05} #define the sites
species_definitions['t'] = {'site_names': ['111'], 'total':0.95}

We also need to decide whether we want to use the (111) or (211) adsorption energies as descriptors. The proper way to do this would be to check the quality of the scaling relations and see which shows a better correlation to the parameters. However, lets just stick with the (211) sites for now.

Here we have assumed that there are 5% step sites, and 95% terrace sites. Now we can run mkm_job.py, and after a lot of fussing the model should converge. The new output looks like:

which clearly shows Pt and Pd as the best CO oxidation catalysts (as we would expect). It is a little worrying that Ag is predicted to be better than Rh, but this could be due to neglecting some mechanism (e.g. O-O-CO), neglecting zero-point and free energy contributions for adsorbates, lack of adsorbate-adsorbate interactions, or issues with the DFT input energies.

Free Energy Diagrams¶

A common way to evaluate or diagnose simpler microkinetic models is to examine the free energy diagrams that went into its creation. In the setup file, we can define any number of reaction mechanisms like the following:

rxn_mechanisms = {  # these are 1-indexed
   "steps": [1, 1, 2, 3, 3],
   "terraces": [4, 4, 5, 6, 7, 7],
}

Here we have defined two reaction mechanisms that follow the catalytic cycle of CO oxidation on steps and terraces. The array for each mechanism is composed of the 1-indexed reaction numbers as described in rxn_expressions. You can use the reverse of a given elementary step by prepending the index with a negative sign.

To actually plot the free energy diagrams, we add the following lines to mkm_job.py:

...
ma = analyze.MechanismAnalysis(model)
ma.energy_type = 'free_energy' #can also be free_energy/potential_energy
ma.include_labels = False #way too messy with labels
ma.pressure_correction = False #assume all pressures are 1 bar (so that energies are the same as from DFT)
ma.include_labels = True
fig = ma.plot(plot_variants=['Pt'], save='FED.png')
print(ma.data_dict)  # contains [energies, barriers] for each rxn_mechanism defined
...

This uses CatMAP’s built-in automatic plotter to generate free energy diagrams for your defined reaction mechanisms on all surfaces by default. For clarity, we are choosing to only plot a subset of these surfaces with the plot_variants=['Pt'] keyword argument. For electrochemical systems using ThermodynamicScaler, plot_variants instead refers to an array of voltages at which to plot free energy diagrams. The resulting plot is fairly simplistic, but feel free to generate your own nicer-looking free energy diagrams using the dictionary provided in ma.data_dict, which stores the values of free energies and barriers for each defined reaction mechanism.

The resulting free energy diagram is a good way to quickly determine if the results of your microkinetic model match with your expectations from its free energy inputs.

Sensitivity Analyses¶

Of course the kinetic models we are building follow the golden rule of mathematical modeling: garbage in, garbage out (i.e. your model is only as good as its inputs). Even if you have the correct mechanism and active site configuration, the results will not make sense if the data in the energy tables is inaccurate. However, in order to refine these inputs it is often useful to know which ones are most important. This can be analyzed using sensitivity analyses.

Rate Control¶

The degree of rate control is a powerful concept in analyzing reaction pathways. Although many varieties exist, the version published by Stegelmann and Campbell is the most general and is implemented in the micro-kinetics module. In this definition we have:

\(X_{ij} = \frac{\mathrm{d} \log(r_i)}{\mathrm{d} (-G_j/kT)}\)

where \(X_{ij}\) is the degree of rate control matrix, \(r_i\) is the rate of production for product i, \(G_j\) is the free energy of species j, k is Boltzmann’s constant, and T is the temperature. A positive degree of rate control implies that the rate will increase by making the species more stable, while a negative degree of rate control implies the opposite.

In order to get the degree of rate control we need to add it as an output_variable in mkm_job.py:

...
mkm_file = 'CO_oxidation.mkm'
model = ReactionModel(setup_file=mkm_file)
model.output_variables += ['production_rate','rate_control']
model.run()
...

We also want to make a plot to visualize the degree of rate control:

mm = analyze.MatrixMap(model)
mm.plot_variable = 'rate_control'
mm.log_scale = False
mm.min = -2
mm.max = 2
mm.plot(save='rate_control.pdf')

The MatrixMap class is very similar to the VectorMap, except that it is designed to handle outputs which are 2-dimensional. This is true of the rate_control (and most other sensitivity analyses) since it will have a degree of rate control for each gas product/intermediate species pair. We set the min/max to -2/2 here since we know that degree of rate control is of order 1. In fact it is bounded by the number of times an intermediate appears on the same side of an elementary step. In this case that is 2, since O2* → 2O* (O* appears twice on the RHS). We could also just let the plotter decide the min/max automatically, but this is sometimes problematic due to numerical issues with rate control.

Now we can run the code. You should see that the initial guesses are proving successful for each point, but you will probably notice that the code is executing significantly slower (factor of ~16). The reason for this will be discussed later. Unlike rates/coverages, the rate control will not converge quicker with a previous solution as an initial guess. In this case it may be desirable to load in the results of a previous simulation directly like:

mkm_file = 'CO_oxidation.mkm'
#model = ReactionModel(setup_file=mkm_file)
#model.output_variables += ['production_rate','rate_control']
#model.run()

model = ReactionModel(setup_file=mkm_file.replace('mkm','log'))

In general this is a good way to re-load the results of a simulation without recalculating it. Regardless, the rate control plot looks like:

This shows us that the rate is decreased when O* or CO* are bound more strongly (depending on descriptor values). Conversely, the rate can be increased by lowering the energy of the O-CO transition state, or sometimes by binding O* more strongly at the (211) site. The effect of lowering the O-O transition-state varies depending on where the surface is in descriptor space.

While these types of analyses are useful, they should be used with caution. If the energies of other intermediates change considerably then it could result in those intermediates controlling the rate. Furthermore, as discussed by Nørskov et. al, there are underlying correlations beneath the parameters, so if one wants to optimize a catalyst these must also be considered.

Other Sensitivity Analyses¶

Similar to the degree of rate control, the degree of selectivity control can also be defined:

\(X^S_{ij} = \frac{\mathrm{d}\log(s_i)}{\mathrm{d}(-G_j/kT)}\)

where \(X^S_{ij}\) is the degree of selectivity control, and \(s_i\) is the selectivity towards species i. This can be included analogously to rate control by adding ‘selectivity_control’ to the output variables and analyzing with the MatrixMap class.

There is also the reaction order with respect to external pressures of various gasses, given mathematically by:

\(R_{ij} = \frac{\mathrm{d} \log(r_i)}{\mathrm{d} \log(pj)}\)

where \(p_j\) is the pressure of gas species j. This can also be included in the same way as rate_control and selectivity control by including “rxn_order” in the output variables.

Numerical Issues in Sensitivity Analyses¶

All sensitivity analyses implemented in the micro-kinetics module are calculated via numerical differentiation. This causes them to be very slow. Furthermore, the fact that numerical differentiation is notoriously sensitive to the “infinitesimal” number used to calculate the derivative, combined with the extreme stiffness of the sets of differential equations behind the kinetic model, can lead to issues. The two most common are:

Jacobian Errors¶

You may sometimes notice that the model will give output like:

initial_evaluation: success - initial guess at point [ 2.71, 3.36]
rootfinding_iteration_3: fail - stagnated or diverging (residual = 5.22501330063e-13)
jacobian_evaluation: fail - stagnated or diverging (residual = 5.22501330063e-13). Assuming Jacobian is 0.
initial_evaluation: success - initial guess at point [ 2.71, 3.04]

This implies that the coverages for the unperturbed parameters failed when used as an initial guess for the perturbed parameters. Given that the perturbation size is, by default, 1e-14, this should only happen if the system is extremely stiff. However, its not impossible. Usually you can figure out what you need to know even when you skip the points where the Jacobian fails, but in the case you really need it converged at every point, you can decrease the “perturbation_size” attribute of the reaction model. When specifying perturbations below 1e-14 it is probably a good idea to do this using the multiple-precision representation as:

model.perturbation_size = model._mpfloat(1e-16)

Diverging or Erroneous sensitivities¶

It is also not uncommon for the sensitivities to diverge to extremely large numbers, or just appear to be random numbers. This generally happens if the perturbation size is too small so that there is no measurable change in the values of the function. The best thing to do here is to tune the perturbation size to a slightly larger number and hope for convergence. Sometimes this does not work, in which case it might also be necessary to increase the precision and decrease the tolerance of the model by many orders of magnitude (see Refining Numerical Accuracy).

Refining Numerical Accuracy¶

A final way to refine a kinetic model is via changing the numerical parameters used for convergence, etc. A few of these parameters will be briefly discussed here:

resolution¶

The resolution determines the number of points between the min/max of the descriptor space. It can be a single number (same resolution in both directions) or a list the length of the number of descriptors. The latter case allows taking a higher resolution in one dimension vs. the other, which is useful if the descriptors have very different scales. It is also worth mentioning that a single-point calculation can be done by setting the resolution to 1. It is important to find a resolution that is fine enough to capture all the features in descriptor space, but of course higher resolution requires more time. It is also worth mentioning that if you want to refine the resolution it is good to pick a number like 2*old_resolution - 1 since this allows you to re-use all the points from the previous solution.

The CO oxidation volcano is shown below at a resolution of 29 (as opposed to 15):

It looks nicer, but doesn’t give much new insight.

decimal_precision¶

This parameter represents the numerical accuracy of the kinetic model. The solutions are found using a multiple-precision representation of numbers, so it is possible to check them to “arbitrary accuracy”. Of course the model will run slower as the decimal_precision is increased, but if the precision is not high enough then the results will not make sense. Generally a decimal_precision of ~100 is sufficient, but for complex models, or when the sensitivity analyses do not behave well, the decimal_precision sometimes needs to be increased upwards of 200-300 digits. If the solutions are correct it should be possible to increase this number arbitrarily and continue to quickly refine the precision of the solutions.

tolerance¶

The tolerance is the maximum rate which is considered 0 by the model. Thus the tolerance should be set to several orders of magnitude below the lowest rate which is relevant for the model. Usually something on the order of 1e-50 to 1e-35 is sufficient. However, when dealing with a model where the maximum rate is very low, or when trying to make sensitivity analyses more accurate, it may be necessary to decrease the tolerance to as low as 1e-decimal_precision. Similar to the decimal_precision, if the solutions are correct then it should be possible to arbitrarily decrease the tolerance (although it should never be lower than 1e-decimal_precision).

max_rootfinding_iterations¶

This determines the maximum number of times the Newton’s method rootfinding algorithm can iterate. It is generally safe to set this to a very high number, since if the algorithm begins to diverge (or even stops converging) then it will automatically exit. Usually a number around 100-300 is practical. This parameter does not affect the solutions of the model, just if/how long the model takes to converge.

max_bisections¶

This determines the number of times a distance between two points in descriptor space can be “bisected” when looking for a new solution. For example, if we know the solution at (0,0) and want the solution at (0,1) the first thing to try is using the (0,0) solution as an initial guess. If that fails, the line will be bisected, and the (0,0) solutions will be tried at (0,0.5). If this fails, then (0,0.25) is tried. This continues a maximum of max_bisections times before the module gives up. This is a “desperation” parameter since it is the best way to get a model to converge, but can be very slow. It is best to start with a value of 0-3, and then slowly increase until the algorithm can find a solution at all points. If the number goes above ~6 then it is an indication that there is something fundamentally wrong with the convergence critera (i.e. the solution oscillates) and that there is no steady-state solution.

Like max_rootfinding_iterations, max_bisections will not change the overall answers to the model, but will determine if/how long it takes to converge.

Code Overview¶

Descriptor based analysis is a powerful tool for understanding the trends across various catalysts. In general, the rate of a reaction over a given catalyst is a function of many parameters - reaction energies, activation barriers, thermodynamic conditions, etc. The high dimensionality of this problem makes it very difficult and expensive to solve completely and even a full solution would not give much insight into the rational design of new catalysts. The descriptor based approach seeks to determine a few “descriptors” upon which the other parameters are dependent. By doing this it is possible to reduce the dimensionality of the problem - preferably to 1 or 2 descriptors - thus greatly reducing computational efforts and simultaneously increasing the understanding of trends in catalysis.

The “catmap” Python module seeks to standardize and automate many of the mathematical routines necessary to move from “descriptor space” to reaction rates. The module is designed to be both flexible and powerful. A “reaction model” can be fully defined by a configuration file, thus no new programming is necessary to change the complexity or assumptions of a model. Furthermore, various steps in the process of moving from descriptors to reaction rates have been abstracted into separate Python classes, making it easy to change the methods used or add new functionality. The downside of this approach is that it makes the code more complex. The purpose of this guide is to explain the general structure of the code, as well as its specific functions and how to begin using it.

Useful Definitions: (Note symbols may not appear in Safari/IE)

• descriptor - variable used to describe reaction kinetics at a high level. Most commonly these are the binding energies of atomic constituents in the adsorbates (e.g. carbon/nitrogen/oxygen adsorption energies) or other intermediates. However, in the general sense other variables such as electronic structure parameters (e.g. d-band center) or thermodynamic parameters (e.g. temperature/pressure) could also be descriptors.
• descriptor space ( \(D\) ) - the space spanned by the descriptors. Usually \(\mathbb{R}^2\).
• parameter space ( \(P\) ) - the space spanned by the full reaction parameters for the reaction model. Usually \(\mathbb{R}^{2n}\) where n is the number of elementary steps (one reaction energy and one reaction barrier per elementary step).
• reaction rates ( \(r\) ) - an n-vector of rates corresponding to each of the n elementary reactions.
• descriptor map - a map of some variable (rates,coverages,etc.) as a function of descriptor space.
• reaction model - a set of elementary steps and conditions which define the physics of the kinetic system, along with the mathematical methods and assumptions used to move from “descriptor space” to reaction rates.
• setup file - a file used to define the reaction model
• input file - a file used to store other data which is likely common to many reaction models (e.g. energetics data)

Code Structure:

• The interface to the code is handled through the ReactionModel class. This class acts as a messenger class which broadcasts all its attributes to the other classes used in the kinetics module. The class also handles common functions such as printing reactions/adsorbates or comparing/reversing elementary steps. The attributes of ReactionModel are automatically synchronized with the parser, scaler, solver, and mapper so it is useful to think of ReactionModel as a “toolbox” where all the necessary information and common functions are stored.
• The Parser class serves to extend the “setup file” by reading in various quantities from an “input file”. Technically the use of a parser is optional, but in practice it is extremely helpful for reading in common data such as adsorption energies or vibrational frequencies rather than re-typing them for every reaction model.

The process of creating a “descriptor map” is abstracted into three general processes, which are handled by the following classes within the kinetics module:

Scaler: Projects descriptor space into parameter space: D → P

Solver: Maps parameter space into reaction rates: P → r

Mapper: Moves through descriptor space. This becomes important for practical reasons since it is often necessary to use a solution at one point in descriptor space as an initial guess for a nearby point.

The ThermoCorrections class is responsible for applying thermodynamic corrections to the electronic energies which are used as direct inputs. This includes the contributions of entropy/enthalpy and zero-point energy due to temperature, pressure, or other thermodynamic variables.

There are also a variety of analysis classes which allow for automated analysis and visualization of the reaction model. These include the VectorMap and MatrixMap classes which create plots of outputs as a function of descriptor space. The VectorMap class is designed for outputs in the form of vectors (rates, coverages, etc.) while the MatrixMap is designed for outputs in the form of matrices (rate control, sensitivity analyses, etc.). The analysis folder also contains MechanismAnalysis which creates free energy diagrams, and ScalingAnalysis which can give a visual representation of how well the scaler projects descriptor space to parameter space.

Using the code:

Some examples can be found in the Tutorials, and these should explain the syntax necessary and serve as a good starting point. The currently implemented features are also briefly described below in order to allow a better understanding of the demos and creating original reaction setup files.

Using the kinetics module to conduct a descriptor analysis requires (at least) 2 files: the “setup file” which defines the reaction model, as well as another script to initialize the ReactionModel class and conduct analyses. The “setup file” is generally a static file (i.e. it is not a program) while the submission script will be an actual python program. Setup files typically end in ”.mkm” for micro-kinetic model (although this is not required) while submission scripts end in .py since they are just python scripts. In addition it is very useful to also have an “input file” which contains the raw data about the energetics of the reaction model. An example of how to create an input file based on a table-like format is given in the Generating an Input File tutorial.

Each class described in the Code Structure section will require some specialized parameters. Some of these parameters are common to all variants of these classes, while others are specific to certain implementations. The possible inputs for the currently implemented variants of each class are listed below. Required attributes are +underlined+.

• ReactionModel:
  • rxn_expressions - These expressions determine the elementary reaction, and are the most important part of the model. They must be defined unless the elementary_rxns, adsorbate_names, transition_state_names, and gas_names are all explicitly defined since the rxn_expressions are parsed into these 3 attributes. It is much easier to just define rxn_expressions, although it is important to note the syntax. There must be spaces between all elements of each expression (i.e. C*+O* is not okay, but C* + O* is), and species ending with _g are gasses by default. Adsorbed species may end with * or _x where * designates adsorption at the “s” site (by default), while _x designates adsorption at the “x” site (note that “x” may be any letter except “g”, and that X* and X_s are equivalent). Transition-states should include a -, and reactions with a transition-state are specified by ‘IS <-> TS -> FS’ while reactions without a transition-state are defined as ‘IS -> FS’ (where IS,TS,FS are expressions for the Initial/Transition/Final State). When the model initializes it checks the expressions for mass/site balances, and if it finds that they are not balanced it will raise an exception. [list of strings]. Instead of specifying rxn_expressions the following attributes may instead be defined:
    • elementary_rxns - list version of rxn_expressions. These will be automatically populated if rxn_expressions are defined. [list of lists of lists]
    • adsorbate_names - list of adsorbate names included in the analysis. Automatically populated if rxn_expressions are defined. [list of strings]
    • transition_state_names - list of transition-state names included in the analysis. Automatically populated if rxn_expressions are defined. [list of strings]
    • gas_names - list of gas names included in the analysis. [list of strings]
  • surface_names - list of surface names to be included in the analysis. [list of strings]
  • species_definitions - This is a dictionary where all species-specific information is stored. The required information will vary depending on the scaler/thermo corrections/solver/mapper used, and the “parser” generally fills in most information. However, there are a few things which generally need to be supplied explicitly:
    • species_definitions[site][‘site_names’] (where *site* is each site name in the model) - A list of “site names” which correspond to *site*. If the TableParser (default) is being used then the “site names” must also match the designations in the “site_name” column. For example, if you want the “s” site to correspond to the energetics of an adsorbate at a (211) site, and (211) sites are designated by ‘211’ in the site_name column of the input_file, then this would be specified by: species_definitions[‘s’] = {‘site_names’:[‘211’]}. Similarly, if you wanted the ‘t’ site to correspond to ‘fcc’ or ‘bridge’ sites then you could specify: species_definitions[‘t’] = {‘site_names’:[‘fcc’,’bridge’]}.
    • species_definitions[site][‘total’] (where *site* is each site name in the model) - A number to which the total coverage of *site* must sum. For example, if you wanted to have a total coverage of 1 with 10% ‘s’ sites and 90% ‘t’ sites (with the same site definitions as above) you would specify: species_definitions[‘s’] = {‘site_names’:[‘211’],’total’:0.1} and species_definitions[‘t’] = {‘site_names’:[‘fcc’,’bridge’],’total:0.9}.
    • species_definitions[gas][‘pressure’] (where *gas* is each gas name in the model including the trailing _g) - The pressure of each gas species in bar. For example, if you wanted a carbon monoxide pressure of 10 bar and hydrogen pressure of 20 bar you would specify: species_definitions[‘CO_g’][‘pressure’] = 10 and species_definitions[‘H2_g’][‘pressure’] = 20. Note that for some situations you may instead need to specify a ‘concentration’,’approach_to_equilibrium’, or some other key, but in almost every situation some method for obtaining the gas pressures must be specified for each gas in the model.
  • temperature - temperature used for the analysis. May not be defined if ThermodynamicScaler is being used with temperature as a descriptor. [number in Kelvin]
  • descriptor_names - names of variables to be used as descriptors. [list of strings]
  • descriptor_ranges - Used for mapping through descriptors space. Specify the limits of the descriptor values. Should be a list equal in length to the number of descriptors where each entry is a list of 2 floats (min and max for that descriptor). [list of lists of floats].
  • resolution - Used for mapping through descriptor space. Resolution used when discretizing over descriptor_range. [int]
  • parser - name of class to use for solver. Defaults to TableParser. [string]
  • mapper - name of class to use as a mapper. Defaults to MinResidMapper. [string]
  • scaler - name of class to use for scaler. Defaults to GeneralizedLinearScaler. [string]
  • solver - name of class to use for solver. Defaults to SteadyStateSolver. [string]
  • thermodynamics - name of class to use for thermodynamic corrections. Defaults to ThermoCorrections. [string]
  • data_file - file where large outputs will be saved as binary pickle files. Defaults to ‘data.pkl’ [filepath string]
  • numerical_representation - determines how to store numbers as binary. Can be ‘mpmath’ for multiple precision or ‘numpy’ for normal floats. Note that ‘numpy’ rarely works. Defaults to ‘mpmath’. [string]
• Parser:
  • input_file - file where input data is stored. File must be in the correct format for the parser used. See Generating an Input File for more information.
• Scaler:
  • gas_thermo_mode - Approximation used for obtaining gas-phase free energy corrections. Defaults to ideal_gas. Other possibilities are: shomate_gas (use Shomate equation), zero_point_gas (zero-point corrections only), fixed_entropy_gas (include zero-point and assume entropy is 0.002 eV/K) , frozen_gas (no corrections), frozen_zero_point_gas (no zero-point and entropy is 0.002 eV/K). [string]
  • adsorbate_thermo_mode - Approximation used for obtaining adsorbate free energy corrections. Defaults to harmonic_adsorbate (use statistical mechanics+vibrational frequencies). Other possibilities are: zero_point_adsorbate (zero-point corrections only), frozen_gas (no corrections). [string]
• Solver:
• SteadyStateSolver:
  • decimal_precision - number of decimals to explicitly store. Calculation will be slightly slower with larger numbers, but will become completely unstable below some threshhold. Defaults to 50. [integer]
  • tolerance - all rates must be below this number before the system is considered to be at “steady state”. Defaults to 1e-50. [number]
  • max_rootfinding_iterations - maximum number of times to iterate the rootfinding algorithm (multi-dimensional Newtons method). Defaults to 50. [integer]
  • internally_constrain_coverages - ensure that coverages are greater than 0 and sum to less than the site total within the rootfinding algorithm. Slightly slower, but more stable. Defaults to True. [boolean]
  • residual_threshold - the residual must decrease by this proportion in order for the calculation to be considered “converging”. Must be less than 1. Defaults to 0.5. [number]
• Mapper:
• MinResidMapper:
  • search_directions - list of “directions” to search for existing solutions. Defaults to [ [0,0],[0,1],[1,0],[0,-1],[-1,0],[-1,1],[1,1],[1,-1],[-1,-1] ] which are the nearest points on the orthogonals and diagonals plus the current point. More directions increase the chances of findinga good solution, but slow the mapper down considerably. Note that the current point corresponds to an initial guess coverage provided by the solver (i.e. Boltzmann coverages) and should always be included unless some solutions are already known. [list of lists of integers]
  • max_bisections - maximum number of time to bisect descriptor space when moving from one point to the next. Note that this is actually the number of iterations per bisection so that a total of 2max_bisections<> points could be sampled between two points in descriptor space. Defaults to 3. [integer]
  • descriptor_decimal_precision - number of decimals to include when comparing two points in descriptor space. Defaults to 2. [integer]
• ThermoCorrections:
  • thermodynamic_corrections - corrections to apply. Defaults to [‘gas’,’adsorbate’]. [list of strings]
  • thermodynamic_variables - variables/attributes upon which thermo corrections depend. If these variables do not change the corrections will not be updated. Defaults to [‘temperatures’,’gas_pressures’]. [list of strings]
  • frequency_dict - used for specifying vibrational frequencies of gasses/adsorbates. Usually populated by the parser. Defaults to {}. [dictionary of string:list of numbers in eV]
  • ideal_gas_params - parameters used for ase.thermochemistry.IdealGasThermo. Defaults to catmap.data.ideal_gas_params. [dictionary of string:string/int]
  • fixed_entropy_dict - entropies to use in the static entropy approximation. Defaults to catmap.data.fixed_entropy_dict. [dictionary of string:float]
  • atoms_dict - dictionary of ASE atoms objects to use for ase.thermochemistry.IdealGasThermo. Defaults to ase.structure.molecule(gas_name). [dictionary of string:ase.atoms.Atoms]
  • force_recalculation - re-calculate thermodynamic corrections even if thermodynamic_variables do not change. Slows the code down considerably, but is useful for sensitivity analyses where thermodynamic variables might be perturbed by very small amounts. Defaults to False. [boolean]
• Analysis:
• MechanismAnalysis:
  • rxn_mechanisms - dictionary of lists of integers. Each integer corresponds to an elementary step. Elementary steps are indexed in the order that they are input with 1 being the first index. Negative integers are used to designate reverse reactions. [dictionary of string:list of integers]

Accessing and reformatting output¶

The purpose of this tutorial is to provide an explanation of how to access the raw output of CatMAP. The plotting functions included in CatMAP are not meant to provide publication-quality figures by default, but rather to provide a quick way of analyzing the output and determining whether or not the solution is correct. If you want to create more complex plots, or more beautiful plots, then you will likely want to reformat the raw data and plot it with software of your choice (MATLAB, etc.) or at least know how to access the matplotlib figure object if you are brave enough to directly edit the figure in python.

The simplest approach to post-processing is to convert the CatMAP output into a data table and use methods you are familiar with. We will use the CO oxidation example from tutorial 3 to provide some concrete context. If you have run Tutorial 3 you should have the file “CO_oxidation.log” in the tutorial 3 directory. Assuming you finished the tutorial, this file should have outputs for the coverage, rate, production_rate, and rate_control. We will look at how to read in each of these and convert them to more conventional tables.

First, lets take a look at how output is stored natively by CatMAP. Create a script with the following commands in the tutorial 3 directory:

from catmap.model import ReactionModel

model = ReactionModel(setup_file='CO_oxidation.log')

which reads the model (along with outputs) into the script. The model, after its run, will have “map” attributes for each output variable. The map is a python list structured like:

[
[[x_0,y_0],[out_0_0, out_0_1, ... , out_0_n]],
[[x_1,y_1],[out_1_0, out_1_1, ... , out_1_n]],
...,
[[x_m,y_m],[out_m_0, out_m_1, ... , out_m_n]]
]

where x and y represent descriptor values, and out_i_j represents the output vector at each point. For example, the coverage map on the CO oxidation model. Add the following to the script:

for pt, cvgs in model.coverage_map:
    print 'descriptors:', pt
    print 'coverages', cvgs

If you run this you will get a bunch of numbers like:

descriptors [1.9289202008928572, 2.232142857142857]
coverages [mpf('2.2800190488939763e-5'), mpf('0.020278026143493406'), mpf('2.2031908471388691e-7'), mpf('1.4713059814195941e-5'), mpf('0.18181715627634512')]
descriptors [0.1428571428571428, 2.553571428571429]
coverages [mpf('2.7893704760253176e-27'), mpf('0.050000000000000003'), mpf('2.8850535476689032e-25'), mpf('1.4466885932549904e-12'), mpf('0.94999999999855326')]
descriptors [3.0, 1.267857142857143]
coverages [mpf('0.049999987577682677'), mpf('4.4864510401468476e-25'), mpf('0.69752183650972651'), mpf('4.653706806446857e-10'), mpf('1.789519988373698e-17')]
descriptors [2.1830357142857144, 2.8750000000000004]
coverages [mpf('6.4052886383458481e-12'), mpf('0.024816133631270126'), mpf('3.4777261296001614e-13'), mpf('6.8031802914834851e-9'), mpf('0.32310023705530582')]
...

You should note that the coverages are “mpf” objects, which indicates that they are multiple precision. You probably also don’t know which intermediate each coverage corresponds to. Try the following:

labels = model.output_labels['coverage']
for pt,cvg in model.coverage_map:
    print 'descriptors',pt
    print 'intermediates',labels
    print 'coverages', [float(c) for c in cvg]

which gives the slightly more readable output:

descriptors [1.9289202008928572, 2.232142857142857]
intermediates ('CO_s', 'O_s', 'CO_t', 'O2_t', 'O_t')
coverages [2.2800190488939762e-05, 0.020278026143493406, 2.203190847138869e-07, 1.471305981419594e-05, 0.1818171562763451]
descriptors [0.1428571428571428, 2.553571428571429]
intermediates ('CO_s', 'O_s', 'CO_t', 'O2_t', 'O_t')
coverages [2.7893704760253176e-27, 0.049999999999999996, 2.885053547668903e-25, 1.4466885932549903e-12, 0.9499999999985532]
descriptors [3.0, 1.267857142857143]
intermediates ('CO_s', 'O_s', 'CO_t', 'O2_t', 'O_t')
coverages [0.049999987577682675, 4.486451040146847e-25, 0.6975218365097264, 4.653706806446857e-10, 1.7895199883736977e-17]
descriptors [2.1830357142857144, 2.8750000000000004]
intermediates ('CO_s', 'O_s', 'CO_t', 'O2_t', 'O_t')
coverages [6.405288638345848e-12, 0.024816133631270124, 3.477726129600161e-13, 6.803180291483484e-09, 0.3231002370553058]
...

Based on this, you can probably see how to create a text table containing coverage outputs. All the other outputs follow the same basic format; however, there are a few tricky situations when looking at other outputs. For example, the “labels” for reaction-specific quantities (rates, rate constants, etc.) are actually lists which need to be flattened into strings. Even more difficult are “matrix” outputs like rate control, where the output is a list of lists rather than a single list. To make life easier I have created the following script which should create a tab-separated text table from any output (.log) file created by CatMAP. Just place this script into the output directory, and run it with the name of the output of interest as its first argument.

from glob import glob
import sys
from catmap.model import ReactionModel

output_variable = sys.argv[1]
logfile = glob('*.log')
if len(logfile) > 1:
    raise InputError('Ambiguous logfile. Ensure that only one file ends with .log')
model = ReactionModel(setup_file=logfile[0])

if output_variable == 'rate_control':
    dim = 2
else:
    dim = 1

labels = model.output_labels[output_variable]

def flatten_2d(output):
    "Helper function for flattening rate_control output"
    flat = []
    for x in output:
        flat+= x
    return flat

#flatten rate_control labels
if output_variable == 'rate_control':
    flat_labels = []
    for i in labels[0]:
        for j in labels[1]:
            flat_labels.append('d'+i+'/d'+j)
    labels = flat_labels

#flatten elementary-step specific labels
if output_variable in ['rate','rate_constant','forward_rate_constant','reverse_rate_constant']:
    str_labels = []
    for label in labels:
        states = ['+'.join(s) for s in label]
        if len(states) == 2:
            new_label = '<->'.join(states)
        else:
            new_label = states[0]+'<->'+states[1]+'->'+states[2]
        str_labels.append(new_label)
    labels = str_labels

table = '\t'.join(list(['descriptor-'+d for d in model.descriptor_names])+list(labels))+'\n'

for pt, output in getattr(model,output_variable+'_map'):
    if dim == 2:
        output = flatten_2d(output)
    table += '\t'.join([str(float(i)) for i in pt+output])+'\n'

f = open(output_variable+'_table.txt','w')
f.write(table)
f.close()

This should give you the ability to import CatMAP output into pretty much any other analysis or plotting program. However, if you are a matplotlib loyalist you may want to try to edit the figure objects directly, or perhaps even exploit the plotting capabilities of CatMAP to plot some “map” other than those created by CatMAP. For example, lets say that for whatever reason we wanted to plot the coverage of CO* times the rate of CO2 formation. We can do this by creating a python script:

from catmap.model import ReactionModel
from catmap.analyze import VectorMap

log_file = 'CO_oxidation.log'
model = ReactionModel(setup_file=log_file)

CO_cvg_CO2_rate_map = []
CO_idx = model.output_labels['coverage'].index('CO_s')
CO2_idx = model.output_labels['production_rate'].index('CO2_g')

for i,pt_cvg in enumerate(model.coverage_map):
    pt_rate = model.production_rate_map[i]
    pt,cvg = pt_cvg
    pt_i,rate = pt_rate
    assert pt == pt_i #ensure that points are the same

    CO_cvg = cvg[CO_idx]
    CO2_rate = rate[CO2_idx]
    CO_cvg_CO2_rate_map.append([pt,[CO_cvg*CO2_rate]]) #multiply the two and store in new map

model.CO_cvg_CO2_rate_map = CO_cvg_CO2_rate_map #trick the model into thinking it has this output
model.output_labels['CO_cvg_CO2_rate'] = ['theta_CO*r_CO2']

vm = VectorMap(model)
vm.plot_variable = 'CO_cvg_CO2_rate' #tell the model to plot the output you just created
vm.log_scale = True #rates should be plotted on a log-scale
vm.min = 1e-25 #minimum rate to plot
vm.max = 1e3 #maximum rate to plot
vm.threshold = 1e-25 #anything below this is considered to be 0
vm.subplots_adjust_kwargs = {'left':0.2,'right':0.8,'bottom':0.15}
fig = vm.plot(save='CO_cvg_CO2_rate.pdf')

If you run this script you will have a CatMAP-style plot of the CO* coverage multiplied by the CO2 formation rate. If you want to make post-processing modifications to the plot, then you should note that the output of the VectorMap.plot function is actually a matplotlib.figure object. You can get the handles for each axis by iterating through the figure.axes attribute. Sometimes it is convenient to label each axis the first time through to know which one you are editing. For example, add the following lines to the script:

for j,ax in enumerate(fig.axes):
    ax.annotate(str(j), [0.05,0.9], color='w', xycoords='axes fraction')

fig.savefig('figure_with_axes_labels.pdf')

Now if you look at the plot you will see the main axis is labeled 0, while the colorbar is 1. You can then edit the axes properties using matplotlib. Manipulating matplotlib figures is beyond the scope of this tutorial, but there is plenty of good documentation at http://matplotlib.org/.

Using Thermodynamic Descriptors¶

The previous tutorials have explained how to create a kinetic model which uses binding energies of key intermediates as “descriptors” in order to create a “map” of catalytic activity across different types of surfaces. This type of analysis is useful for catalyst screening, but it is also often interesting to take a closer look at a single catalyst surface and explore the effect of reaction conditions. In order to do this in CatMAP we will use the same framework as discussed in the Code Overview. The only difference is that instead of using energies as descriptors we will use “thermodynamic variables” (temperature, pressure for now) as descriptors. Naturally this also means that we will need to use a different “scaler” which, instead of creating a linear map uses some pre-defined expressions to add in the proper free energy contributions.

To make this clearer we will continue with the CO oxidation example and examine CO oxidation on Pt(111) as a function of temperature and pressure. Again, the purpose is not to make a publishable analysis but to become familiar with the functionality of CatMAP. We will be able to use a “submission script” (mkm_job.py) very similar to the one from Refining a Microkinetic Model.

from catmap import ReactionModel

mkm_file = 'CO_oxidation.mkm'
model = ReactionModel(setup_file=mkm_file)
model.output_variables += ['production_rate'] model.run()

from catmap import analyze

vm = analyze.VectorMap(model)
vm.plot_variable = 'production_rate'

#tell the model which output to plot

vm.log_scale = True #rates should be plotted on a log-scale
vm.min = 1e-25 #minimum rate to plot
vm.max = 1e3 #maximum rate to plot
vm.threshold = 1e-25 #anything below this is considered to be 0
vm.subplots_adjust_kwargs = {'left':0.2,'right':0.8,'bottom':0.15}
vm.plot(save='production_rate.pdf')

For simplicity we will only look at the “production_rate” output variable in this tutorial. Other outputs can be calculated/analyzed as described in tutorials 2-3, Creating a Microkinetic Model and Refining a Microkinetic Model.

We will use a reaction model similar to the simple one used in Creating a Microkinetic Model since this will run faster than the more complicated models. In order to refine this to something more sophisticated you can follow the strategy outlined in Refining a Microkinetic Model. This also means that we can re-use the energies.txt file from Creating a Microkinetic Model, although we can make it even shorter since we only care about the Pt(111) energies:

surface_name        site_name       species_name    formation_energy        bulk_structure  frequencies     other_parameters        reference
None        gas     CO2     2.45    None    [1333,2349,667,667]     []      "Angew. Chem. Int. Ed., 47, 4835 (2008)"
None        gas     CO      2.74    None    [2170]  []      "Energy Environ. Sci., 3, 1311-1315 (2010)"
None        gas     O2      5.42    None    [1580]  []      Falsig et al (2012)
Pt  111     O       1.62    fcc     []      []      Falsig et al (2012)
Pt  111     CO      1.7     fcc     []      []      "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Pt  111     O-CO    4.04    fcc     []      []      "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Pt  111     O-O     5.35    fcc     []      []      Falsig et al (2012)

Note that we could also leave all the other energies in since the parser will just ignore any lines it doesn’t need. Finally, there is the “setup file” - CO_oxidation.mkm. This is where all of the necessary changes must be made in order to use “thermodynamic descriptors”. We will start with the same setup file from Creating a Microkinetic Model, but we need to make the following changes:

• Change the scaler Add the following line somewhere in the setup file (usually near the beginning): scaler = 'ThermodynamicScaler' As you can probably guess this tells CatMAP to use the “ThermodynamicScaler” class to move from “descriptor space” to “parameter space” (see Code Overview). This class will use whatever “thermo modes” are defined for gas phase/adsorbate species in order to add free energy contributions (see Creating a Microkinetic Model).

• Choose the relevant surface Next we need to tell CatMAP which surface to use. For this example we will look at Pt(111). To do this we just need to change the “surface_names” variable:

  surface_names = ['Pt']
  

  Note that the (111) surface is already selected due to the line:

  species_definitions['s'] = {'site_names': ['111'], 'total':1}
  

• Change the descriptors Now we need to tell the “ThermodynamicScaler” which two variables we will be using for descriptors, and we need to modify the ranges over which to vary these descriptors. Currently only temperature and pressure are implemented, although there is also an option to use log(pressure) as discussed later. For now lets look at temperatures from 400 - 1000 K and pressures from 1e-8 to 1000 bar:

  descriptor_names= ['temperature','pressure']
  descriptor_ranges = [[400,1000],[1e-8,1e3]]
  

• Modify temperature/pressure to be compatible - In Creating a Microkinetic Model) we used a model where temperature and pressure were explicitly specified. This doesn’t really make sense now, since we are varying these two variables. The temperature ends up not really mattering since it will be over-written as CatMAP moves through descriptor space; however, just to be unambiguous its good practice to delete the following line:

  temperature = 500 #Temperature of the reaction
  

Finally, we need to tell CatMAP how to handle the pressures. Previously we just defined “static pressures” for each gas-phase species, but that doesn’t make sense if the total pressure is varying. In order to get around this we instead specify “concentrations” for each gas-phase species:

species_definitions['CO_g'] = {'concentration':2./3.}
species_definitions['O2_g'] = {'concentration':1./3.}
species_definitions['CO2_g'] = {'concentration':0}

Note that this “concentration” is not normalized - the total pressure of a gas at any total pressure will be given by concentration*P where P is the total pressure. Thus, if the concentrations do not sum to 1 then the “pressure” axis will be incorrect.

After making these changes we can run the “submission script” with:

sh python mkm_job.py

which should give the usual kind of output. When it finishes you should see the following “production_rate.pdf” in the folder:

As expected, the temperature dependence is much more drastic than the pressure dependence. In many cases it makes more sense to look at pressure dependence on a log scale. This is easily achieved by changing the descriptor names/ranges:

descriptor_names= ['temperature','logPressure']
descriptor_ranges = [[400,1000],[-8,3]]

Note that the “log” in this notation refers to a base 10 logarithm so that the plot produced is the same as before, but with pressure on a log scale. If we now run the submission script we get the following output:

This looks a little nicer than the previous plot since the low pressure behavior has higher resolution.

We can see from this tutorial that it is fairly easy to move between a micro-kinetic model for a screening study and one for a “reaction condition” study (and vice-versa). Only a few lines of the “setup file” need to be changed. This is one advantage of CatMAP - once you setup a reaction model once you can re-use it for several different analyses.

Including Adsorbate-Adsorbate Interactions¶

In all of the previous tutorials we have made some assumptions about how the adsorbed species interact with the surface lattice and each other. Specifically, we assume that every species adsorbs at a single site and does not interact with its neighbors. In this tutorial we will examine these assumptions more closely and consider some strategies for making more advanced assumptions. Please note that while most of the other features of CatMAP are relatively well tested, the features outlined in this tutorial have had very little testing/error checking so please be careful if using them in your research (and naturally let the developers list know about any issues).

It is also worth mentioning that while “ideal” mean-field models like the ones shown previously are usually mathematically well-behaved and thus exhibit “existence and uniqueness” (there is one and only one steady-state solution), all bets are off once adsorbate interactions are included. Thus it is likely possible to construct models which have oscillating solutions or no solutions at all (within the physical bounds). Naturally a system with no solution at all is unphysical, so as long as you make good assumptions then you should be okay. However, oscillating solutions are slightly trickier and it can be very hard to determine if/when this is an issue since CatMAP uses root finding to get the steady-state solution. I have yet to come across an example where CatMAP cannot find a solution to a realistic model, but be aware that it could happen.

We will continue with the CO oxidation example to avoid having to create completely new files. The tutorials will seek to demonstrate the options available within CatMAP, even if these options are not physical for the CO oxidation system, so don’t take the outputs too seriously. The tutorial is divided into 2 parts which do not necessarily need to be followed sequentially:

• Multi-site adsorbates and maximum coverages
• Coverage dependent adsorption energies

For both sections we will use a submission script (mkm_job.py) similar to the one in the Using Thermodynamic Descriptors tutorial:

from catmap import ReactionModel

mkm_file = 'CO_oxidation.mkm'
model = ReactionModel(setup_file=mkm_file)
model.output_variables += ['production_rate']
model.run()

from catmap import analyze
vm = analyze.VectorMap(model)
vm.plot_variable = 'production_rate' #tell the model which output to plot
vm.log_scale = True #rates should be plotted on a log-scale
vm.min = 1e-25 #minimum rate to plot
vm.max = 1e8 #maximum rate to plot
vm.threshold = 1e-25 #anything below this is considered to be 0
vm.subplots_adjust_kwargs = {'left':0.2,'right':0.8,'bottom':0.15}
vm.plot(save='production_rate.pdf')

vm.plot_variable = 'coverage' #tell the model which output to plot
vm.log_scale = False #coverage should not be plotted on a log-scale
vm.min = 0 #minimum coverage
vm.max = 1 #maximum coverage
vm.subplots_adjust_kwargs = {'left':0.2,'right':0.8,'bottom':0.15,'wspace':0.6}
vm.plot(save='coverage.pdf')

We use the same input file (energies.txt) and start with the same setup file (CO_oxidation.mkm) from Tutorial 2 . Note that each section assumes you are starting with the “fresh” CO_oxidation.mkm file from Tutorial 2.

Multi-site adsorbates and maximum coverages¶

One very simple way of including adsorbate interactions is the “hard sphere exclusion” model. Actually, we have already assumed this (each site can have only 0 or 1 adsorbate), but we can also extend the assumption to allow adsorbates to occupy more than 1 site, or to set a “maximum coverage” for an adsorbate. Lets take a look at the multi-site adsorbates first:

Multi-site adsorption¶

CatMAP includes the ability to allow and adsorbate to adsorb to multiple sites of the same type. For instance, lets say that we want to force CO to adsorb to 2 sites. This is achieved by editing the “rxn_expressions”:

rxn_expressions = [

               '2*_s + CO_g -> CO*',
               '2*_s + O2_g <-> O-O* + *_s -> 2O*',
               'CO* +  O* <-> O-CO* + 2* -> CO2_g + 3*',

                   ]

Note that we had to edit the number of sites on the CO adsorption and desorption reactions in order to make everything consistent. The next thing we need to do is tell CatMAP that CO occupies 2 sites, so that it doesn’t get confused about the site balance:

species_definitions['CO_s'] = {'n_sites':2}

Now we can run the model and get the following coverages:

and rate:

If we compare these to Tutorial 2 then we can see that the CO* coverage is suppressed and there is more O* in the bottom left of the plot. This is what we would expect to happen when we require an adsorbate to have an extra free site to adsorb.

There is one thing worth noting about this approach. If coverage was defined as number of CO per number of surface sites we would expect the maximum CO coverage to be 0.5 since it occupies 2 sites. However, it is clear from the plot that the coverage goes to 1. That is because we have re-defined the number of “total sites” to be a factor of 2 less for CO so that the maximum coverage of an adsorbate is always 1. This is equivalent to assuming that the probability of an adsorbate which occupies 2 sites reacting with another adsorbate is a factor of 2 higher since the site it sits on is 2 times larger. Depending on the system this may be a poor assumption, but it is the only option currently implemented in CatMAP.

Maximum coverages¶

There may also be circumstances where we wish to constrain certain adsorbates to have a maximum coverage. This can easily be achieved by adding the line:

species_definitions['CO_s'] = {'max_coverage':0.5}

to CO_oxidation.mkm. However, when you run the submission script you will notice that after a lot of complaining CatMAP will give the following:

mapper_iteration_3: fail - no solution at 99 points.

This is the first time we have encountered a model that will not converge. Normally we would try to get convergence by increasing “max_bisections” or other parameters as discussed in Tutorial 3. However, in this case it is hopeless. This is probably because there is no solution within the bounds we have defined (which means they are not physical). This isn’t too surprising since we just made the constraint up. We can still take a look at the points that did converge in coverages.pdf:

This is pretty consistent with what we might expect. The model converges everywhere that CO coverage is less than 0.5 in the unconstrained solution, but starts to break down when the constraint limits the CO coverage to less than what is found in the unconstrained solution. Although this approach does not really make physical sense here, there could be systems where it does. In these cases CatMAP should be able to find a valid solution. Note that the “max_coverage” only pertains to one adsorbate, and does not inhibit competitive adsorption (i.e. you could have CO coverage of 0.5 and O coverage of 0.5).

Coverage dependent adsorption eneriges¶

A more powerful method for including adsorbate-adsorbate interactions is to allow adsorption energies to depend on the coverages the adsorbates. This is still relatively crude compared to an explicit lattice method like kinetic Monte Carlo, but it should provide a good picture of the first-order effects of coverage . Of course there are many ways to parameterize such a model, but there is currently only one option implemented in CatMAP - the “first order adsorption energy” model. We will first introduce the model, then look at how to use it in CatMAP, and finally show an example of how to apply it to the CO oxidation example.

First order adsorption energy model¶

In this model we assume that adsorption energies follow the following relationship:

\[\begin{split}E_{i} =& E_{i}^{i} + \sum_{j} {\cal{F} (|\theta|_{j}) \varepsilon_{ij} \theta_{j}} \\ |\theta|_j =& \sum_{{\mathrm{site}_k}={\mathrm{site}_j}} \theta_{k}\end{split}\]

where \(E_{i}\) is the generalized formation energy for species \(i\), |θ|j is the total coverage of occupied sites for the site on which adsorbate \(j\) is adsorbed, \(\varepsilon_{ij}\) is the “interaction matrix”, and F is the “interaction response function” which is usually some smoothed piecewise linear function and will be discussed later. When computing the Jacobian matrix for the system we will also need the derivative of the energy with respect to coverages. This is given by:

\[\frac{\partial E_{i}}{\partial \theta_l} = \sum_{j} \varepsilon_{ij} \left( \frac{\rm{d} \cal{F} \left(|\theta|_j\right)}{\rm{d}|\theta|_j} \frac{\rm{d}|\theta|_j}{\rm{d}\theta_{l}} \theta_{j} + \cal{F} \left(|\theta|_j\right)\delta_{jl} \right)\]

The model is called “first order” since it includes only one term of coverage dependence, and this term is first order in the coverage (and \(\cal{F}\)).

We see that in order to calculate adsorption energies we need the function \(\cal{F}\), and the matrix \(\varepsilon\). We will also end up needing the derivative of the function \(\cal{F}\) w.r.t. \(|θ|_j\) . These two quantities will be discussed below.

Interaction response function¶

The “interaction response function” determines how much the adsorption energy changes as a function of the total coverage at a site. This is necessary because adsorption energies often follow non-linear behavior as a function of coverage. Some examples of possible response functions are shown below:

The “linear”, “piecewise_linear”, and “smooth_piecewise_linear” are implemented in CatMAP, while the “linear_step” is a hypothetical model which could be implemented. Depending on the complexity of the interaction response function it will require some parameters. The parameters are “site specific”, so that if you have a model with step sites and terrace sites you could use different “cutoffs” for the piecewise linear response function. However, the parameters do not vary by adsorbate which limits the complexity of the model.

Interaction matrix¶

The other key input for the “first order” interaction model is the “interaction matrix”, εij. There are two types of terms in this matrix - “self interaction” terms (\(\varepsilon_{ii}\)) and “cross interaction” terms (\(\varepsilon_{ij} (i\ne j)\)). As the name suggests the “self interaction” terms tell how much an adsorbate interacts with itself, while “cross interactions” tell how much it interacts with other adsorbates. The interaction matrix is symmetric (\(\varepsilon_{ij} = \varepsilon_{ji}\)). The values for the matrix are determined by fitting to data. If the differential binding energies are available at various coverages then the fitting is very straight-forward. However, in most cases density functional theory (DFT) will be used to calculate binding energies. Due to the discreet nature of coverages in DFT, it is impossible to calculate differential binding energies. Instead, average binding energies are calculated and used to obtain the interaction parameters. The definition of average binding energy is:

\[\begin{split}\bar{E}_i =& \frac{\int_{0}^{\theta_i} E_{i} \rm{d} \theta_{i}}{|\theta|_i} \\ =& \frac{\int_0^{\theta_i} \left(E_i^0 + \sum_k \cal{F}(|\theta|_k)\varepsilon_{ik}\theta_k \right) \rm{d}\theta_i}{|\theta|_i}\end{split}\]

from this we can solve for the self-interaction parameters:

These equations look nasty at first site, but the form of \(\cal{F}\) is usually simple enough that they aren’t so intimidating. CatMAP also includes the ability to fit the self-interaction functions automatically, as discussed later.

The cross interaction terms are very costly to calculate, since they require many DFT calculations (two per adsorbate per adsorbate, or Nadsorbates2). For this reason it is common to use some approximations. The most common approximations are:

\[\begin{split}{\mathrm{geometric\ mean:}}\ \varepsilon_{ij} =& \sqrt{|\varepsilon_{ii} \varepsilon_{jj}|} \\ {\mathrm{arithmetic\ mean:}}\ \varepsilon_{ij} =& (\varepsilon_{ii} + \varepsilon_{j})/2\\ {\mathrm{ neglect:}}\ \varepsilon_{ij} =& 0\end{split}\]

In CatMAP the cross interaction terms are between adsorbate-adsorbate and adsorbate-transition_states. This means that the interaction matrix is actually (Nadsorbate+Ntransition-state)2. Both self and cross interactions between transition states are neglected since by definition their coverage will always be negligible. However, cross interactions between adsorbates and transition-states is not negligible. Since we don’t have any self-interaction parameters for transition-states, we need some method of estimating them. This can be done by:

• transition-state scaling: transition-state scaling is used to estimate the cross parameters so that the transition-state scaling relation holds (best approximation if available)
• initial state: use the cross interaction parameters corresponding to the initial state (forward barrier is static)
• final state: use the cross interaction parameters corresponding to the final state (reverse barrier is static)
• intermediate state: use some weighted average of the initial and final state interactions (usually 0.5).
• neglect: assume to be 0 (all barriers decrease)

Implementation in CatMAP¶

The implementation of adsorbate-interactions requires modifications at many levels of CatMAP - specifically, the solver, scaler, and parser have been modified for the first order interaction model. However, the place that the “interactions” fit most logically into the design of CatMAP is in the “thermodynamics” since technically this is a modification of the assumption of non-interacting adsorbates. For this reason, most of the implementation has been abstracted into a class which is in the thermodynamics directory. If you are not developing then this is not a concern, but just be aware that in order to use the “first order” interaction model (or others in the future) you need to ensure that the parser, scaler, and solver are compatible. Currently the default parser (TableParser), scaler (GeneralizedLinearScaler), and solver (SteadyStateSolver) are the only ones compatible with interaction models.

Relevant Attributes¶

The implementation relies on the following attributes of the reaction model:

• adsorbate_interaction_model: Determines which model to use. Currently can be ‘ideal’ (default) or ‘first_order’

• interaction_response_function: The function \(F\) from above. Can be ‘linear’, ‘piecewise_linear’, or ‘smooth_piecewise_linear’. Can also be a callable function which takes the total coverage of a site as its first argument and the “interaction_response_parameters” dictionary as a **kwargs argument. Must return the value and derivative of the function at the specified total coverage.

• interaction_response_parameters: This is a dictionary of argument names/values to be used in the “interaction_response_function”. The “interaction_response_parameters” can be specified as an attribute of the ReactionModel (use the same parameters for all sites) or as a key/value in the “species_definitions” dictionary for different sites (use different parameters for different sites).

• self_interaction_parameters: These are the self interaction parameters for a given adsorbate. They should be specified as a key/value in the species definition entry for the adsorbate. The key should be “self_interaction_parameters” and the value should be a list of the same length as “surface_names”. The parameters should be entered for each surface in the same order that the surfaces appear in “surface_names”. If an interaction parameter is not available for a surface then None should be entered.

• cross_interaction_parameters: Cross-interaction parameters can be input as a key/value pair in the species_definitions entry for one of the two adsorbates. The key should be “cross_interaction_parameters” and the value should be a dictionary where the key is the other adsorbate of the cross interaction pair and the value is a list of the same length as “surface names” where the parameters are input similar to the “self_interaction_parameters”. The following is an example of how this might appear in the setup file for neglecting CO-O cross interactions on Pt, Pd, and Rh:

  ...
  surface_names = ['Pt','Pd','Rh']
  ...
  species_definitions['CO_s'] = {'cross_interaction_parameters':{'O_s':[0,0,0]}}
  ..
  

  Explicitly specifying cross interaction parameters is optional. Any parameters that are not explicitly specified will be estimated as specified by “cross_interaction_mode”. Note that one parameter must be specified for each surface, and None can be used if a value is unknown. The “surface_names” in the example above is the same “surface_names” which defines the surfaces in the entire model, and thus should only be defined once.

• max_self_interaction: Practically it is sometimes found that the self interaction parameter should not be larger than some cutoff. This can be specified by setting the “max_self_interaction” key in the species_definitions dictionary for the adsorbate to either the numerical value or a name of one of the “surface_names” to automatically bound the interaction parameter at the value for that surface. The attribute can also be added directly to the ReactionModel in order to bound all self interaction parameters (for this, especially, using a “surface name” as a bound is recommended).

• cross_interaction_mode: The cross interaction mode tells CatMAP how to approximate cross interaction parameters that are not specified explicitly. The values can be: ‘geometric_mean’ (default), ‘arithmetic_mean’ or ‘neglect’ as described above.

• transition_state_cross_interaction_mode: Similar to “cross_interaction_mode” but for transition-states. Can be ‘transition_state_scaling’, ‘initial_state’, ‘final_state’, ‘intermediate_state’, or ‘neglect’ as described above. Using ‘intermediate_state’ will assume a weight of 0.5, or you can specify ‘intermediate_state(X)’ to set a weight of X.

• interaction_scaling_constraint_dict: The equivalent of “scaling_constraint_dict” but for interaction parameters. By default, “scaling_constraint_dict” will be used, but constraints which force slopes to be positive/negative will be removed since sign changes are expected between the “adsorbate scaling” coefficient and the interaction parameter scaling coefficient. Any parameter which does not have scaling constraints defined will be set to the “default_interaction_constraints” attribute ([None,None,None] by default).. Cross interaction parameter names are defined by ‘A&B’ and can appear in either order. For example, to constrain the cross parameter between O* and CO* to scale only with the first descriptor we could do:

  interaction_scaling_constraint_dict['O_s&CO_s'] = [None,0,None]
  

  Defining the constraint for ‘CO_s&O_s’ would be equivalent. See Tutorial 2 for a refresher on the syntax of constraint definitions.

• non_interacting_site_pairs: Pairs of site names which are not interacting. All cross interactions between adsorbates on these sites will be set to 0. For example, to prevent cross interactions between adsorbates on the ‘s’ and ‘t’ site:

  non_interacting_site_pairs = [['s','t']]
  

  The order of adsorbates does not matter since the interaction matrix is symmetric.

• interaction_strength: All interaction parameters will be multiplied by this. Should be floatable. Defaults to 1. Useful for getting model to converge.

• interaction_fitting_mode: Determines how to construct fits to raw data. Can be None (default), ‘average_self’. None implies that CatMAP should not try to automatically do any fitting because the parameters are explicitly specified. Using “average_self” will fit the self interaction parameters assuming that there are coverage-dependent average adsorption energies in the input file.

In addition, the interaction matrix can be included as an output for error-checking (this is recommended since the interaction model is still relatively new). Simply include “interaction_matrix” in the “output_variables” and analyze the output as described in Tutorial 2 <../tutorials/creating_a_microkinetic_model>.

CO Oxidation Example¶

Including coverage-dependent interactions¶

First, lets assume that we already know the self-interaction parameters and want to include coverage dependent adsorbate interactions on top of the model discussed in Tutorial 2. In order to do this we need to add the following to the CO_oxidation.mkm setup file:

adsorbate_interaction_model = 'first_order' #use "first order" interaction model
interaction_response_function = 'smooth_piecewise_linear' #use "smooth piecewise linear" interactions
species_definitions['s']['interaction_response_parameters'] = {'cutoff':0.25,'smoothing':0.01}
#input the interaction paramters
#surface_names = ['Pt', 'Ag', 'Cu','Rh','Pd','Au','Ru','Ni'] #surface order reminder
species_definitions['CO_s'] = {'self_interaction_parameter':[3.248, 0.965, 3.289, 3.209, 3.68, None, None, None]}
species_definitions['O_s'] = {'self_interaction_parameter':[3.405, 5.252, 6.396, 2.708, 3.87, None, None, None]}
max_self_interaction = 'Pd' #self interaction parameters cannot be higher than the parameter for Pd
transition_state_cross_interaction_mode = 'transition_state_scaling' #use TS scaling for TS interaction
cross_interaction_mode = 'geometric_mean' #use geometric mean for cross parameters

If we use the same submission script as before we should get the following outputs for coverage and rate:

We can see that the coverages change much more gradually, as expected. The rate volcano is a little worrying since it now predicts Pt and Pd to be some of the worst catalysts. However, we recall that the reaction mechanism here is very simplistic, and that we are only looking at the (111) surfaces. A more realistic analysis would reveal that Pt and Pd are still the optimal catalysts, as shown by Grabow et. al..

Including scaled cross interactions¶

In the previous section we used the “geometric mean” approximation to get the cross-interaction terms from the self-interaction terms. While this is a good first approximation, it is sometimes not sufficiently accurate. In order to account for this it is possible to also include some cross-interaction terms as scaled parameters. For a very unphysical example, we will neglect cross-interactions between adsorbed O and CO, and between adsorbed CO and the O-O transition-state. This can be done by adding the following to the species definition for adsorbed CO:

species_definitions['CO_s'] = {'self_interaction_parameter':[3.248, 0.965, 3.289, 3.209, 3.68, None, None, None],
                    'cross_interaction_parameters':{'O_s':[0,0,0,0,0,0,0,0],'O-O_s':[0,0,0,0,0,0,0,0]}}

We note that the cross interactions could have equivalently been defined in the species definitions for adsorbed O and the O-O transition-states (where CO_s would be the key of the cross_interaction_parameters dictionary) but it is easier to group them both into the CO_s definition. If we now run the submission script we get the following outputs:

These results are not physical because there is no reason to expect that CO does not interact with O or O-O, but they do illustrate the syntax for specifying arbitrary cross interaction parameters. Note that the vector of zeroes here is the same length as the number of surfaces. Much like the self interaction parameters, the values of these cross interactions must be in the same order as the order of the surface names, with any unknown parameters given as None. If actual parameters were input instead of zeroes, then they would also be estimated using scaling relations in the same way the self interaction parameters are.

Using CatMAP to fit self interactions¶

In many cases the interaction parameters will not be available and they must be determined from some set of coverage dependent raw data. In this situation it is very convenient to have the interaction matrix automatically fit to this raw data to avoid typos and round off error in the interaction parameters. CatMAP is capable of automatically constructing this fit for the self-interaction parameters of the “first order” model described above. Fitting the second order parameters is more complicated, and should be done manually. In order to create the automatic fit it is necessary to have the energies as a function of coverage. For example, we can use the following input file with some soon to be published data for coverage dependent O and CO adsorption, along with transition-state energies from previous examples. Note that there is now a new “coverage” column:

surface_name    site_name   species_name    coverage    formation_energy    bulk_structure  frequencies other_parameters    reference
None    gas CO2 0   2.46    None    [1333,2349,667,667] []  "NIST"
None    gas CO  0   2.77    None    [2170]  []  "Energy Environ. Sci., 3, 1311-1315 (2010)"^M
None    gas O2  0   5.42    None    [1580]  []  NIST^M
Rh  111 O   0.25    0.54    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pt  111 O   0.25    1.62    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pd  111 O   0.25    1.55    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Cu  111 O   0.25    1.08    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Ag  111 O   0.25    2.04    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Au  111 O   0.25    2.75    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Rh  111 O   0.50    0.76    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pt  111 O   0.50    1.9 fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pd  111 O   0.50    1.88    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Cu  111 O   0.50    1.755   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Ag  111 O   0.50    2.585   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Au  111 O   0.50    3.065   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Rh  111 O   0.75    1.043   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pt  111 O   0.75    2.243   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pd  111 O   0.75    2.237   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Cu  111 O   0.75    2.423   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Ag  111 O   0.75    3.147   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Au  111 O   0.75    3.5 fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Rh  111 O   1.00    1.31    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pt  111 O   1.00    2.592   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pd  111 O   1.00    2.665   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Cu  111 O   1.00    2.925   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Ag  111 O   1.00    3.55    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Au  111 O   1.00    3.797   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Rh  111 CO  0.25    1.25    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pt  111 CO  0.25    1.49    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pd  111 CO  0.25    1.3 fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Cu  111 CO  0.25    2.53    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Ag  111 CO  0.25    2.96    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Rh  111 CO  0.50    1.58    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pt  111 CO  0.50    1.915   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Ag  111 CO  0.50    3.07    fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Rh  111 CO  1.00    2.193   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pt  111 CO  1.00    2.473   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Pd  111 CO  1.00    2.335   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Cu  111 CO  1.00    3.455   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Ag  111 CO  1.00    3.247   fcc []  []  Khan et. al. Parameterization of an interaction model for adsorbate-adsorbate interactions
Rh  111 O-CO    0.25    3.1 fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Pt  111 O-CO    0.25    4.04    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Pd  111 O-CO    0.25    4.2 fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Cu  111 O-CO    0.25    4.18    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Ag  111 O-CO    0.25    5.05    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Au  111 O-CO    0.25    5.74    fcc []  []  "Angew. Chem. Int. Ed., 47, 4835 (2008)"
Rh  111 O-O 0.25    3.79    fcc []  []  Falsig et al (2012)
Pt  111 O-O 0.25    5.35    fcc []  []  Falsig et al (2012)
Pd  111 O-O 0.25    5.34    fcc []  []  Falsig et al (2012)
Cu  111 O-O 0.25    4.74    fcc []  []  Falsig et al (2012)
Ag  111 O-O 0.25    5.98    fcc []  []  Falsig et al (2012)
Au  111 O-O 0.25    7.22    fcc []  []  Falsig et al (2012)

Naturally the transition-states only need to be computed at a single coverage, since they do not have self interaction parameters. It is also worth noting that even if not all metals have coverage dependent data, they can still be included in the analysis (their interaction parameters will be estimated from scaling).

You can find the above data table as coverage_energies.txt in the folder for this tutorial. If you make the following changes to CO_oxidation.mkm then the parameters will be determined automatically:

input_file = 'coverage_energies.txt'
interaction_fitting_mode = 'average_self'

The “average_self” fitting mode refers to the fact that the energies in the input file are average adsorption energies, and that only the self interaction parameters will be fit. The only other option is “differential_self” which assumes that the inputs are differential adsorption energies and fits self interaction parameters.

Now, if you run mkm_job.py then you will get the same output as when the self interaction parameters were input manually (because the parameters were pre-determined by this procedure). If you want to view the parameters then you can do so by looking at the “self_interaction_parameter_dict” in the CO_oxidation.log. You should notice that they match the parameters that were input manually earlier. The advantage of the automatic fitting procedure is that any changes in the “interaction response function” will automatically be compensated for in the fit (i.e. if the smoothing value is decreased, cutoff is changed, etc.). It also makes it easier to generalize the model to inputs coming from different calculation methods, functionals, etc.

Electrochemistry¶

Electrochemistry in CatMAP is implemented as a set of thermodynamic corrections in the enthalpy_entropy module. To run models containing electrochemical reactions in CatMAP, the user needs to provide the following in addition to all of the typical components of a CatMAP model described in the tutorials:

• Some sort of voltage-dependent species. Currently, pe and ele are supported as fictitious gas molecules that represents the free energy of a proton-electron pair at 0V vs RHE or an electron at 0V vs SHE, respectively. As part of the Computational Hydrogen Electrode, the pe gas species should have the same free energy of 1/2 of a molecule of H2 gas. This species’s energy (usually 0 if an H2 reference is used) can be defined manually in the input file, but defaults to half the free energy of the H2 gas molecule in your system (if using pe) or 0 if using “ele”. These are special electrochemistry-specific gas-phase species that should ignore your choice of free energy correction scheme.
• In the reaction definitions, any potential-dependent reaction should have the pe_g gas molecule or ele_g in the initial (if reductive) or final (if oxidative) state. In addition, in the reaction definitions and input file, transition states for these reactions should contain pe or ele in the species name. For example, the reduction of adsorbed oxygen to adsorbed OH can be written as 'O* + pe_g <-> O-pe* -> OH*'.
• An alternative formulation for electrochemical transition states is available via a special notation shown in the ORR tutorials. Instead of explicitly defining a transition state species in your input file and using it in a reaction expression, you may instead write something like ^0.26eV_a where the 0.26eV represents the free energy barrier of elementary step at that step’s equilibrium/limiting potential on site a.
• A voltage (vs RHE) defined as a parameter in the .mkm file e.g. voltage = -0.2. This is also a global thermodynamic variable, and can be manipulated as a descriptor as shown in the ORR_thermo tutorial.
• A transfer coefficient defined as a parameter in the .mkm file e.g. beta = 0.5. This value can be defined for each elementary step through a reaction-specific flag like ;prefactor=None, beta=0.45 at the end of the corresponding reaction expression string. See the ORR tutorials for an example of this notation.
• In the .mkm file, electrochemical_thermo_mode can be specified to one of three possible values (defaults to simple_electrochemical): simple_electrochemical, hbond_electrochemical, and hbond_with_estimates_electrochemical. “simple” only adds free energy corrections to adsorbates and transition states to account for voltage and beta. “hbond” will take a default or user-provided hbond_dict to correct each species for hydrogen bonding stabilization (see catmap/data/hbond_data.py for the default hbond_dict). “hbond_with_estimates” attempts to estimate a hydrogen bonding correction for a given species based on its chemical formula. This is a very crude process that is described in catmap/thermodynamics/enthalpy_entropy.py.

There are a few provided examples of using electrochemistry in CatMAP in the tutorials/electrochemistry directory. The HER example shows the hydrogen evolution reaction in the low-coverage limit, and the provided README file explains the details of the simplest of electrochemical reactions as done in CatMAP. The ORR example is meant to reproduce the results discussed in Hansen et al. (2014) DOI: 10.1021/jp4100608, which used a home-spun microkinetic model. The README in that directory goes into depth on how you can replicate some of the major features of this work with the provided scripts. The electron example shows how to equivalently use CatMAP for an SHE potential versus an RHE one.

Output Variables¶

There are a large variety of possible output variables catmap can produce as demonstrated by the “output_variables” tutorial, which provides a comprehensive list of them. You can verify that they all work by running the tutorial and inspecting the plots generated or inspect the raw data in the .pkl file. A brief description of some of the more common ones are presented here below.

• free_energy: contains the free energy of all species in the reaction model as calculated by the scaler. This, in conjunction with the MechanismAnalysis module can be useful in diagnosing issues in the model where the free energies of intermediates are unexpected in some way. Along with the related scaler-originated variables rxn_parameter, frequency, electronic_energy, zero_point_energy, enthalpy, and entropy, these are not typically plotted on a contour plot as they are in the tutorial without manual adjustment to the plotting parameters.
• rate_control and selectivity_control: see Sensitivity analyses
• production_rate: The rate of the production of each gas phase species. If the gas is being consumed, consumption_rate or turnover_frequency (which encapsulates both production and consumption rates) may be more useful.
• coverage: Coverage of each adsorbed species. Very useful for quick sanity checks of “does this model output make sense”.
• rate: net rate (forward - backward) for each elementary step. The derived quantities forward_rate and reverse_rate are the absolute values of rate corresponding to the overall reaction direction. Specifically for the forward and reverse rates before they are summed, use directional_rates.

Developer Info¶

To get started developing you obviously need git installed. You can make changes in your local directory, and if you make a change that you think is substantial and general enough that others would be interested you can “push” the change using git (see the git tutorial). If you are going to submit a change, please update the version number in catmap/init.py (see version == x.x.x). The first number in the version is used for major changes, and a 0 represents code which is still not completely stable. The second number represents moderate changes (significant new functionality added), and the final number corresponds to the git submission number. The final number is the most annoying to update, since you need to change it every time before you submit. In order to make this slightly easier you can use the following script:

from subprocess import Popen, PIPE

output = Popen(['git', 'rev-list', 'HEAD','--count'],stdout=PIPE)
rev = output.stdout.read().strip()
init = open('catmap/__init__.py')
init_txt = init.read()
init.close()
new_txt = []
for li in init_txt.split('\n'):
    if '__version__' in li:
        oldversion = li.rsplit('.',1)[0]
        newversion = oldversion + '.' + rev + '"'
        new_txt.append(newversion)
    else:
        new_txt.append(li)

new_init = '\n'.join(new_txt)
init = open('catmap/__init__.py','w')
init.write(new_init)
init.close()

message = raw_input('Submit message:')
os.system('git commit -am "'+message+'"')
os.system('git push')

Place the script in the base directory of the project and run it with python. If anyone knows a better way of keeping the version number synchronized I am open to suggestions.