Add training reactions and enlarge the nodes in the tree for some surface reaction families

[Tingchenlee]

• Purpose:
  Adding more surface training data into kinetic families.

• Method:

  • Adding new libraries
  • Processing the new libraries by using the script in David's RMG-Py branch add_training_reactions (RMG-Py/ipython/kinetics_library_to_training.ipynb)
  • Adding the reactions to the belonged families.

• Notice:

  • The reactions will be marked with their catalytic types and facets.
  • Some pre-exponential factors were calculated from their rate constants
  • Some pre-exponential factor were calculated from the theoretical calculation method proposed by this paper
  • The units of pre-exponential will be converted depending on their metals (surface site density), and the calculation will be recorded in longDesc

• Paper: Kinetic Prefactors of Reactions on Solid Surfaces

  • The value of entropies are from the Molecule Search on the RMG website
  • The below are equations for estimating A factors
    

[Tingchenlee]

Some A factors in Lu_Ir111 library used the default value A =10e13. Maybe we could figure out another way to calculate A with imaginary frequency...?

[davidfarinajr]

Great work @Tingchenlee ! Thanks for these additions! Have you tested this branch by building a model with all of the families that you add training reactions? Some of this is not covered in unit testing, and bugs may only show up when we are adding rules from training to build models

[Tingchenlee]

I think I tested most of them individually after I finished one but I didn’t test all of them at one time. I am testing all of the libraries and families in one model now!

[Tingchenlee]

I had run a couple of models which included CHx, NHx, and O2 species at different temperatures (300K-1500K) and pressure (1 bar - 200 bar). All the models were okay. However, there was a secret space in the Schneider_Rh211 library (previously used "Schneider_Rh211 " rather than "Schneider_Rh211"). I just revised the name of that library and pushed it.

[davidfarinajr]

Thanks for these additions @Tingchenlee! Since the tests pass and we are able to build models with this branch, I approve of this PR. However, there are a few caveats and things we should address either as part of this PR or in the future

1. some of the kinetics trees are too small to leverage the kinetics data added here. For example, a lot of training reactions are added to the Surface_Adsorption_Single family, but the tree is very small:

tree(
"""
L1: Adsorbate
    L2: N
        L3: N=O

L1: VacantSite
"""

We need more nodes in the tree (for C and O for example) for the new reactions to descend to, otherwise we will have a log jam at the top of the tree and only the first reaction with the best rank will be used. So currently some of our trees are too small to properly handle the influx of training data with this PR.

2. averaging SurfaceArrhenius and StickingCoefficient kinetics. I'm not sure if this is an issue, but I suspect it may be. For some of the families (Surface_Adsorption_Single for example), we have StickingCoefficient kinetics in the forward direction and SurfaceArrhenius in the reverse direction. I don't think we properly average the kinetics when you are averaging both StickingCoefficient kinetics and SurfaceArrhenius kinetics. This is something we should investigate.

3. minor detail - I think we should use X for the species label for empty sites in training reactions and not Pt or Cu because it will get reused, and if you have a reaction on a different metal it can be confusing (reaction could be on 'Ni' but species label may be 'Pt')

[ChrisBNEU]

For point 1 on david's comment, @Tingchenlee discovered that you get an error loading the database with the mixed surfaceArrhenius and stickingCoefficient types in the training reactions. Specifically, the error occurs when trying to average these two in the surface adsorption vdw family:

1. X + N2 <=> N2X:

StickingCoefficientBEP(
     A=1.17e-06, 
     n=0, 
     alpha=0, 
     E0=(0,'J/mol'), 
     Tmin=(200,'K'), 
     Tmax=(3000,'K'), 
     comment="""Average of [From training reaction 28 used for N2H4;VacantSite]"""), 
[<Entry index=33 label="N-N">, <Entry index=2 label="VacantSite">]

2. X + H4N2 <=> H4N2X:

SurfaceArrheniusBEP(
     A=(3.464e+21,'cm^3/(mol*s)'), 
     n=0, 
     alpha=0, 
     E0=(4000,'J/mol'), 
     Tmin=(200,'K'), Tmax=(3000,'K'), 
     comment="""From training reaction 13 used for N#N;VacantSite"""), 
[<Entry index=37 label="N#N">, <Entry index=2 label="VacantSite">]

This is not caught by unit tests currently. the error that you get is because the units do not match, but the real problem is, as david said, we probably shouldn't be averaging these two reaction types. I think to determine a workaround, we need to determine how (or if) we should average these two reaction types. The same rates are reported differently by different papers so I think we should be able to average them, but I think that would mean calculating the A, b, and Ea for the sticking coefficient reaction:

I can try to program a work around for this, the actual averaging routine does not look very complicated, but I'd like to know what everyone else thinks of this solution. in the context of this PR we could probably also just convert the offending reactions to surface Arrhenius form.

[davidfarinajr]

I think for adsorption type reactions, the StickingCoeff kinetics model is more "physically meaningful" than SurfaceArrh. Therefore we should convert the SurfaceArrh kinetics for adsorption training reactions to StickingCoeff, then average the StickingCoeffs. What we are missing, however, is a method to convert SurfaceArrh to StickingCoeff. We would need to create this method, then use it to convert SurfaceArrh training reactions

[Tingchenlee]

Thanks, @davidfarinajr and @ChrisBNEU
I tried to add more nodes in the tree of the Surface_Adsorption_Single, Surface_Adsorption_vdW, Surface_Dissociation, and Surface_Dissociation_Beta families since we have more training reactions now. I also dealt with the issue David mentioned (point 3, change surface sites from Cu, Ni, or Pt to X). Then the unit converting errors of the Sticking coefficient and surface Arrhenius happened.
One thing I want to mention is the function of _get_average_kinetics in rmgpy/data/kinetics/rules.py . We don't have the unit of (m^4/mol^2/s) for Surface_Bidentate_Dissociation family, which was mentioned in rmgpy/quantity.py. Should we also add (m^4/mol^2/s) to the _get_average_kinetics function? We now have some training reactions in the Roldan_Cu111 and the Vlachos_Pt111 libraries that their A factors have the unit of (m^4/mol^2/s), e.g.N2H3_X + NH2_X + X <=> [Pt]NN[Pt] + NH3_X

[davidfarinajr]

name = "Surface_Bidentate_Dissociation/groups"
shortDesc = u""
longDesc = u"""
A bidentate species with a single bond dissociates. The reverse reaction is when two adsorbates, each with double-bonds or higher to the surface, come together and form a single bond between them (thereby reducing the bond orders to the surface).

 *1--*2             *1   +  *2
  |  |     ---->    ||       ||
~*3~~*4~~          ~*3~~ +  ~*4~~ 

The rate, which should be in mol/m2/s,
will be given by k * (mol/m2)
so k should be in (1/s)

looks like the Surface_Bidentate_Dissociation units are 1/s in forward and m^2/(mol*s) in reverse. It doesn't look like we currently have a family that can make that reaction N2H3_X + NH2_X + X <=> [Pt]NN[Pt] + NH3_X. If we do add a family for that, then you are right that we would need to add `'m^4/(mol^2s)' to that function (depending on what direction the family is in).

[davidfarinajr]

Opened a PR to address the mixed kinetics types issue in training reactions and was able to build model with this PR DB branch ReactionMechanismGenerator/RMG-Py#2148.

[Tingchenlee]

I have added another library Mhadeshwar_Pt111 since it contains some new NO* mechanisms which I had seen in another paper.

After adding lots of training data, I decided to re-number the entries and also took a look at our previous data. I found some of the reactions used the wrong units in the long description, for example, the calculation of index=16 in Surface_Abstraction_vdW training reactions.py is:

A factor from paper / surface site density of Cu 
1.675e12 m^4/(mol^2 * s) / 2.943e‐5 mol/m^2 = 5.691e16 m^2/(mol*s)

I think the unit of A factor in that paper is 1/s, and the equation becomes:

A factor from paper / surface site density of Cu 
1.675e12 1/s / 2.943e‐5 mol/m^2 = 5.691e16 m^2/(mol*s)

I revised these long descriptions and made some comments but I'm not sure if I'm wrong. Also corrected some of the short and long descriptions I made and fixed some typos in the dictionary.txt files.

some other minor revise:

1. two citations of the training reactions in the Surface_Dissociation family were wrong
2. the A factor of the reaction CH3O* + H* <=> CH3OH_2* + X_4 in the Surface_Dissociation_vdW family is wrong. I also made some comments about the calculation.

[mazeau]

Oops wrong pr

[rwest]

1E13

[Tingchenlee]

Thanks for pointing out this issue, I was having fun with the crashed models these two days...
Just fixed it now!

[rwest]

If they published this in both directions, we should use the forward Adsorption direction sticking coefficient.

[Tingchenlee]

Okay, I have removed all the SurfaceArrhenius reactions in the Surface_Adsorption_Dissociative, the Surface_Adsorption_Single, and the Surface_Adsorption_vdW families and imported the StickingCoefficient reactions.
If the paper only calculated the reaction rates in the reverse direction (SurfaceArrhenius), I keep them in the libraries only.

[Tingchenlee]

I opened #506 (adding surface libraries) and #507 (adding tree nodes) to separate this PR, will add another "adding training data" one after these two are merged.

[Tingchenlee]

Opened RMG-database PR#512 and closed this one