search

ReactionMechanismGenerator / RMG-Py

Python version of the amazing Reaction Mechanism Generator (RMG).
http://reactionmechanismgenerator.github.io/RMG-Py/
Other
400 stars 228 forks source link

Barrier heights not fixed for linear scaling of libraries #2726

Open sevyharris opened 3 weeks ago

sevyharris commented 3 weeks ago

Bug Description

I tried regenerating the cpox mechanism on Rh and found that the O2 dissociation reaction has crazy fast kinetics compared to Emily's original model. Untitled Untitled

I think the barrier heights are the issue here. The reaction is defined in the library in the desorption direction: Untitled

And the barrier height is not raised to the endothermicity of the reaction in the model I generated. The three points plotted are 1. enthalpy of reactants (set to 0 for reference), 2. Ea for the reaction, 3. enthalpy of reaction/enthalpy of products relative to reactants Untitled-1

If you raise the barrier, you get more reasonable kinetics for O2 adsorption. So, I think this reaction's barrier was supposed to be fixed due to the scaling, but fix_barrier_height is intentionally not called for libraries in the rmgpy.rmg.model.CoreEdgeReactionModel class.

How To Reproduce

Run the latest RMG on the main branch using this input file (basically the same as Emily's cpox input file but now using Deutschmann2006_adjusted instead of Deutschmann2006: https://github.com/comocheng/uncertainty_estimator/blob/main/cpox_pt/cpox_rh_20241028/input.py

You can see the mechanism it generated here: https://github.com/comocheng/uncertainty_estimator/blob/main/cpox_pt/cpox_rh_20241028/chem_annotated-surface.inp

You can also just look at the values in Surface/CPOX_Pt/Deutschmann2006_adjusted using this code:

import numpy as np
import os
import copy

import rmgpy.chemkin

import matplotlib.pyplot as plt
%matplotlib inline

database = rmgpy.data.rmg.RMGDatabase()

thermoLibraries=['surfaceThermoPt111', 'primaryThermoLibrary', 'thermo_DFT_CCSDTF12_BAC', 'DFT_QCI_thermo']
reactionLibraries = ['Surface/CPOX_Pt/Deutschmann2006_adjusted', 'BurkeH2O2inArHe']

database.load(
    path = rmgpy.settings['database.directory'],
    thermo_libraries = thermoLibraries,
    transport_libraries = [],
    reaction_libraries = reactionLibraries,
    seed_mechanisms = [],
    kinetics_families = [],
    kinetics_depositories = ['training'],
    depository = False, # Don't bother loading the depository information, as we don't use it
)

my_rxn = database.kinetics.libraries['Surface/CPOX_Pt/Deutschmann2006_adjusted'].entries[10].item
display(my_rxn)
scale_to_metals = ['Pt111', 'Rh111']
linestyles = ['solid', 'dashed']

for a, scale_to in enumerate(scale_to_metals):
    for p in range(len(my_rxn.products)):

        if my_rxn.products[p].is_surface_site():
            thermos = database.thermo.get_thermo_data_from_library(my_rxn.products[p], database.thermo.libraries['surfaceThermoPt111'])
            my_rxn.products[p].thermo = database.thermo.correct_binding_energy(
                thermos[0],
                my_rxn.products[p],
                metal_to_scale_from='Pt111',
                metal_to_scale_to=scale_to
            )
            my_rxn.products[p].thermo = thermos[0]
        else:
            thermos = database.thermo.get_thermo_data_from_library(my_rxn.products[p], database.thermo.libraries['primaryThermoLibrary'])
        my_rxn.products[p].thermo = thermos[0]

    for r in range(len(my_rxn.reactants)):
        thermos = database.thermo.get_thermo_data_from_library(my_rxn.reactants[r], database.thermo.libraries['surfaceThermoPt111'])
        my_rxn.reactants[r].thermo = database.thermo.correct_binding_energy(
            thermos[0],
            my_rxn.reactants[r],
            metal_to_scale_from='Pt111',
            metal_to_scale_to=scale_to
        )
    my_rxn.kinetics = database.kinetics.libraries['Surface/CPOX_Pt/Deutschmann2006_adjusted'].entries[10].data

    plt.plot([1, 2, 3], [0, my_rxn.kinetics.Ea.value_si, my_rxn.get_enthalpy_of_reaction(1000)], linestyle=linestyles[a], label=scale_to)

my_rxn.fix_barrier_height()
plt.plot([1, 2, 3], [0, my_rxn.kinetics.Ea.value_si, my_rxn.get_enthalpy_of_reaction(1000)], linestyle='dotted', label='Rh111 raised barrier')

plt.legend()
plt.xlabel('Reaction Coordinate')
plt.ylabel('H (J/mol)')

Expected Behavior

I expect the O2 adsorption rate to be much slower than it was the last time I ran RMG.

Installation Information

Describe your installation method and system information.