Check for sticking coeff rxns > 1 in rules.py

[ChrisBNEU]

Motivation or Problem

Issue #1981: sometimes rule-generated sticking coefficients can exceed 1, because of the reaction path degeneracy calculations. For example, a rate rule of A = 0.5 applied to ethane losing a hydrogen, increases the sticking coefficient by 6 (A=3.0):

Description of Changes

changed the rules.py file to cap the sticking coefficient at 1.0, and to include a message in the reaction comments detailing why:

Testing

Unit tests run successfully. I'd be open to adding a test but I think it would require changing the setup for rmgpy/data/kineticsTest.py.

Reviewer Tips

Run a surface mechanism with it. I forced my example to have the same problem by changing the corresponding surface dissociation rule from 0.1 to 0.5.

[bjkreitz]

Thanks for putting this togehter. I haven't tested it yet, but I think we need to modify it. The sticking coefficient for non-activated adsorption reactions should never exceed 1, which agrees with the changes you made. However, this is not the case for activated adsorption. All dissociative adsorption steps are activated, whereas non-dissociative adsorption is usually not activated. In the case of activated dissociative adsorption, we do not provide the sticking coefficient in the rule but a parameter to calculate the sticking coefficient via an Arrhenius equation (See for example the Cantera documentation for the equation). So the value we provide can be greater than 1 but the calculated sticking coefficient from the Arrhenius equation should not exceed 1. That said, I think we should just exclude dissociative adsorption from this check.

[davidfarinajr]

Yup, I agree with @bjkreitz here. We want to use this method to check sticking coefficient > 1:

    cpdef double get_sticking_coefficient(self, double T) except -1:
        """
        Return the sticking coefficient (dimensionless) at temperature `T` in K. 
        """
        cdef double A, n, Ea, T0, stickingCoefficient
        A = self._A.value_si
        n = self._n.value_si
        Ea = self._Ea.value_si
        T0 = self._T0.value_si

        stickingCoefficient = A * (T / T0) ** n * exp(-Ea / (constants.R * T))
        if stickingCoefficient < 0:
            raise ValueError("Sticking coefficients cannot be negative, check your preexponential factor.")
        return min(stickingCoefficient, 1.0)

However, this returns min(stickingCoefficient, 1.0), so we have to modify this function or evaluate stickingCoefficient = A * (T / T0) ** n * exp(-Ea / (constants.R * T)) to see if its > 1

[ChrisBNEU]

Thanks for clarifying, David and Bjarne. For anyone else reading I had been conflating the "a" from the expression below with the true sticking coefficient (gamma):

If we have to evaluate the Sticking coefficient, since it is temperature dependent, we won’t be evaluating “sticking coefficient violators” consistently, since in theory you could have a reaction that has a sticking coefficient <1 in a lower temp reactor and >1 in a higher temp one. Maybe that’s fine but I am not sure. You might be able to access the maximum reactor temp that the user has specified in their input file, but that might be tricky.

[davidfarinajr]

Maybe we could report sticking coefficient that are > 1 at some T with the collision rate violators after RMG job? Sticking coefficient will never exceed 1 internally in RMG when it evaluates it since we take min(1, sticking coefficient), but I do not think this is true for Cantera, so we should log it to make users aware.

[bjkreitz]

Yes, I agree with you. That depends on the a factor and the energy barrier. Highest temperature in input would make sense. We can also set an arbitrary evaluation temperature of 1273 or 1500 K for example. There is not much catalysis going on on the surface if the temperature is too high. So it doesn't really matter if the sticking coefficient violates the boundary outside of this temperature range.

[ChrisBNEU]

From the above discussion and the RMG-Cat meeting, I will modify this PR so that we treat "sticking coefficient violators" similar to how we treat collision limit violators in the gas phase. I.e., we will check and report them at the end of the run.

[github-actions[bot]]

This pull request is being automatically marked as stale because it has not received any interaction in the last 90 days. Please leave a comment if this is still a relevant pull request, otherwise it will automatically be closed in 30 days.