RMG Surface Solver does not match Cantera Surface Solver

[sevyharris]

Bug Description

When I run equivalent surface simulations in RMG and Cantera, I get different results, and I can't figure out whether I'm doing something wrong or the solvers are just expected to behave differently.

This matters to me because I'm trying to use RMG's sensitivity calculations for the uncertainty tool, but I can't test any of it against Cantera if the simulated concentrations don't even match.

How To Reproduce

Here's a minimal RMG mechanism and Cantera mechanism of O2+CH4 reacting on Pt111:

short-gas.inp.txt
short-surface.inp.txt
short.yaml
species_dictionary.txt

and an example Jupyter Notebook to run them/see how I'm running them: https://github.com/comocheng/uncertainty_estimator/blob/main/uncertainty_tool_dev/surface_sensitivities/minimal_solver_difference.ipynb

Any and all advice is appreciated here! I would love it if the answer is I just screwed up the units somewhere.

Expected Behavior

I would expect RMG and Cantera to come up with the same results for a mechanism with only 9 reactions in a batch reactor, but maybe I've misunderstood the assumptions here.

The RMG SurfaceReactor claims to be isothermal, and constant volume in the class docstring. (The file docstring says isothermal and isobaric but that looks like it was accidentally left in while copying the SimpleReactor class).

You can get the results to match up for slightly simpler cases, like O2 as the only starting species adsorbing on the surface.

Installation Information

Describe your installation method and system information.

• OS: Ubuntu 24.04.4 LTS
• Installation method: source
• RMG version information:
  • RMG-Py: 9cc2cfd (latest main)
  • RMG-database: ac2cc758d00c59add46528d6ff38b815691b1726

[mjohnson541]

It's not impossible this is a numerics issue, but I think that's unlikely. These are pretty huge relative errors here. I would either compare with RMS or setup the ODE independently to verify which is the black sheep.

General thoughts:
Cantera's default solver is less robust than RMG's
Multi-phase simulations complicate the state vector and depending on how this is handled can cause simulation variables to have widely different values which has implications for ODE solve tolerances
RMG/Cantera might be loading different kinetic parameters

[rwest]

Check assumptions about total volume change, control volume, etc. ?

[sevyharris]

Thanks for these suggestions! I'll give them a try

[sevyharris]

Looks like it was the kinetics.

For the Ninth reaction:

RMG

From the chemkin file (mol/cm^2):

! Reaction index: Chemkin #12; RMG #4
! Library reaction: Surface/CPOX_Pt/Deutschmann2006_adjusted
! Flux pairs: CH4(2), CH3X(26); X(1), OHX(30); OX(24), OHX(30); 
X(1)+OX(24)+CH4(2)<=>OHX(30)+CH3X(26)               5.000000e+18 0.700     10.038   

RMG gets:
SurfaceArrhenius(A=(5e+18,'cm^4/(mol^2*s)'), n=0.7, Ea=(10.038,'kcal/mol'), T0=(1,'K'))
A=5.0e10 m^4/(mol^2*s)

Cantera

Meanwhile the Cantera file reads:

- equation: X(1) + OX(24) + CH4(2) <=> OHX(30) + CH3X(26)  # Reaction 9
  rate-constant: {A: 5.0e+18, b: 0.7, Ea: 10.038}
  note: |
    Reaction index: Chemkin #12; RMG #4
    Library reaction: Surface/CPOX_Pt/Deutschmann2006_adjusted
    Flux pairs: CH4(2), CH3X(26); X(1), OHX(30); OX(24), OHX(30); 

Giving A factor of 5e14 m^5 / kmol^2 / s
A=5.0e8 m^5/(mol^2*s)

If you set the reaction multiplier in Cantera to 100x, the plots match up:

Next Steps

Now I just need to look into why the Cantera file gets a different rate. My money's on the discrepancy between the units: Cantera's m^5/(mol^2*s) vs. RMG's m^4/(mol^2*s)

I think it's because they pick different reference species for calculating the units of k_fwd. But this looks like a ck2yaml problem.

[sevyharris]

Yes, it's the units. Here's my demo notebook on RMG vs. Cantera units for rate coefficients.

For second order reactions with both gas and surface species, there is ambiguity in whether the rate units are m^5/(mol^2 *s) or m^4/(mol^2*s).

In the examples I have, RMG seems to always choose m^4/(mol^2*s), but I probably need a bigger sample size or to figure out where RMG is setting these units.

[sevyharris]

I found a few more examples where RMG and Cantera disagree about the rate coefficient, but I think for different reasons than above. I think the reaction above is a special case where because it's a library reaction, it's using a different set of assumptions when being read/written to Cantera compared to family reactions.

I suspect the examples below (all Sticking Coefficients) are a result of a different implementation in Cantera vs RMG to go from sticking coefficient to rate coefficient (see this section of RMG):

RMG-Py/rmgpy/reaction.py

Line 934 in 9a5490a

if isinstance(self.kinetics, StickingCoefficient):

Just noting the pictures for myself for later: