Adding surface training reactions

[Tingchenlee]

• Purpose:
  Adding more surface training data into kinetic families.
• Method:
  • 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 Mhadeshwar_Pt111 library is for diesel oxidation catalysts (DOC) system and has 52 training reactions (Ea vary depending on the reactions, most of the A factors are in the range of 10^11 to 10^13 (1/s)), not sure if we want to add these 52 reactions now (might affect carbon-related models?).
  • The Vlachos_Pt111 library is for methane oxidation and has 43 training reactions (Ea vary depending on the reactions, most of the A factors are in the range of 10^10 to 10^13 (1/s), sticking coefficient from 0.01 to 1), not sure if we want to add these 43 reactions now (might affect carbon-related models?).
  • Will add other libraries into training data but hold these two libraries for further discussions.

[Tingchenlee]

Run the example of methane oxidation with new training data and the models completed: (RMG-Py branch is master)

• RMG-database branch: add_surface_training
  The final model core has 60 species and 707 reactions
  The final model edge has 185 species and 1322 reactions

• RMG-database branch: master
  The final model core has 61 species and 689 reactions
  The final model edge has 189 species and 1375 reactions

database(
    thermoLibraries=['surfaceThermoPt111', 'primaryThermoLibrary', 'thermo_DFT_CCSDTF12_BAC','DFT_QCI_thermo'], 
    reactionLibraries = [('Surface/CPOX_Pt/Deutschmann2006_adjusted', False)], 
    seedMechanisms = [],
    kineticsDepositories = ['training'],
    kineticsFamilies = ['surface','default'],
    kineticsEstimator = 'rate rules',

)

generatedSpeciesConstraints(
    allowed=['input species','seed mechanisms','reaction libraries'],
    maximumCarbonAtoms=2,
    maximumOxygenAtoms=2,
    maximumNitrogenAtoms=2,
    maximumSurfaceSites=2,
    maximumRadicalElectrons=2,
)

catalystProperties(
    metal = 'Pt111'
)

species(
    label='CH4',
    reactive=True,
    structure=SMILES("[CH4]"),
)

species(
   label='O2',
   reactive=True,
   structure=adjacencyList(
       """
1 O u1 p2 c0 {2,S}
2 O u1 p2 c0 {1,S}
"""),
)

species(
    label='N2',
    reactive=False,
    structure=SMILES("N#N"),
)

species(
    label='vacantX',
    reactive=True,
    structure=adjacencyList("1 X u0"),
)
#----------
# Reaction systems
surfaceReactor(
    temperature=(800,'K'),
    initialPressure=(1.0, 'bar'),
    initialGasMoleFractions={
        "CH4": 0.1,
        "O2": 0.2,
        "N2": 0.7,
    },
    initialSurfaceCoverages={
        "vacantX": 1.0,
    },
    surfaceVolumeRatio=(1.e5, 'm^-1'),
    terminationConversion = { "CH4":0.99,},
    terminationTime=(0.1, 's'),
)

simulator(
    atol=1e-18,
    rtol=1e-12,
)

model(
    toleranceKeepInEdge=0.0,
    toleranceMoveToCore=1e-5,
    toleranceInterruptSimulation=0.1,
    maximumEdgeSpecies=100000,
)

options(
    units='si',
    generateOutputHTML=True,
    generatePlots=True,
    saveEdgeSpecies=True,
    saveSimulationProfiles=True,
)

And run the ammonia oxidation examples (RMG-Py: master):

• RMG-database branch: add_surface_training
  The final model core has 24 species and 88 reactions
  The final model edge has 71 species and 184 reactions

• RMG-database branch: master
  The final model core has 20 species and 39 reactions
  The final model edge has 56 species and 114 reactions

database(
    thermoLibraries=['surfaceThermoPt111', 'primaryThermoLibrary', 'thermo_DFT_CCSDTF12_BAC','DFT_QCI_thermo'], 
    reactionLibraries = [('Surface/Ammonia/Schneider_Pt111', False)], 
    seedMechanisms = [],
    kineticsDepositories = ['training'],
    kineticsFamilies = ['surface','default'],
    kineticsEstimator = 'rate rules',

)

generatedSpeciesConstraints(
    allowed=['input species','seed mechanisms','reaction libraries'],
    maximumOxygenAtoms=2,
    maximumNitrogenAtoms=2,
    maximumSurfaceSites=2,
    maximumRadicalElectrons=2,
)

catalystProperties(
    metal = 'Pt111'
)

species(
    label='NH3',
    reactive=True,
    structure=adjacencyList(
"""
1 N u0 p1 c0 {2,S} {3,S} {4,S}
2 H u0 p0 c0 {1,S}
3 H u0 p0 c0 {1,S}
4 H u0 p0 c0 {1,S}
"""),
)

species(
   label='O2',
   reactive=True,
   structure=adjacencyList(
       """
1 O u1 p2 c0 {2,S}
2 O u1 p2 c0 {1,S}
"""),
)

species(
    label='He',
    reactive=False,
    structure=adjacencyList(
"""
1 He u0 p1 c0
"""),
)

species(
    label='X',
    reactive=True,
    structure=adjacencyList("1 X u0"),
)

species(
    label='H2O',
    reactive=True,
    structure=SMILES("O"),
)

#----------
# Reaction systems
surfaceReactor(
    temperature=(800,'K'),
    initialPressure=(1.0, 'bar'),
    initialGasMoleFractions={
        "NH3": 0.002,
        "O2": 0.01,
        "He": 0.938,
        "H2O":0.05,
    },
    initialSurfaceCoverages={
        "X": 1.0,
    },
    surfaceVolumeRatio=(1.e5, 'm^-1'),
    terminationConversion = { "NH3":0.95,},
    terminationTime=(0.1, 's'),
)

simulator(
    atol=1e-18,
    rtol=1e-12,
)

model(
    toleranceKeepInEdge=0.01,
    toleranceMoveToCore=0.05,
    toleranceInterruptSimulation=0.1,
    maximumEdgeSpecies=100000,
)

options(
    units='si',
    generateOutputHTML=True,
    generatePlots=True,
    saveEdgeSpecies=True,
    saveSimulationProfiles=True,
)

The selectivity of the ammonia oxidation has been improved by adding the training reactions. Without the new training data, the model misses the product of N2O. However, the surface coverage wasn’t improved.

[davidfarinajr]

Thanks @Tingchenlee ! Great work! I was able to generate auto trees for a bunch of the families you added training data. This is a good test to see if there is anything wrong with the new training reactions. Also, you were able to generate models and add rules from training, which is another good test. One thing I did notice, however, is that I get slightly different Ea in J/mol when converting from eV:

For example, you have Ea = 1.49eV = 143770.1J/mol
and using RMG I get

RMG uses conversion factor of 96485.3423627019 to convert to J/mol. Since it is a small difference on the order of a few J/mol (which is much less than the uncertainty of the calculation), I am approving this PR. In the future, we should try to use the same units from the paper, and let RMG do the conversion

[ChrisBNEU]

what I did do

• check the description
• check that the orders of magnitudes for all of the values made sense
• check that I could access the sources with the information here
• found a spot where I need clarification for coverage dependence.

what I did NOT do

• redo any of the calculations (for A-factor or Ea). I don't think it is necessary.
• check every value against the value in the paper because it is too time consuming.
• verify the rank and degeneracy (I think we'd fail a unit test if these were wrong right?)

what I found:

• a few small typos
• a value that did not match the written description. The paper matches the written description so I think it needs to be changed. if somebody can sanity check me please do and I'll update it.
• some reactions don't have a facet because the paper used a single crystal reactor.

I can go back through and check the values against the papers, but that will take significantly more time and I don't think it is worth it. I am going to revise what is here after somebody sanity checks my comment above. I can take over this PR, I think Ting-Chen is probably not going to have time to be working on it anymore and it is almost ready.

[ChrisBNEU]

might be unnecessary, but we should probably specify they use a "single site" reactor model, which is why there is no facet info.

[mazeau]

this looks ok to me and passes all tests locally