[yaml_cantera1] handle negative A factors (also fixes to_cantera)

We add handlers to the to_cantera methods to set the 
allow_negative_pre_exponential_factor property of the Cantera 
object, which Cantera then correctly passes to the input_data
dictionary.

We also have a safety-net in the yaml writing stage, but maybe this 
is unnecessary and slows us a bit. It doesn't seem to get triggered
by the unit tests.

Author: rwest

rmgpy/kinetics/arrhenius.pyx

@@ -272,7 +272,10 @@ cdef class Arrhenius(KineticsModel):
         if arrhenius_class:
             return ct.Arrhenius(A, b, E)
         else:
-            return ct.ArrheniusRate(A, b, E)
+            rate = ct.ArrheniusRate(A, b, E)
+            if A < 0:
+                rate.allow_negative_pre_exponential_factor = True
+            return rate
 
     def set_cantera_kinetics(self, ct_reaction, species_list):
         """
@@ -784,7 +787,10 @@ cdef class ArrheniusBM(KineticsModel):
         Ea = self._E0.value_si * 1000  # convert from J/mol to J/kmol
         w = self._w0.value_si * 1000  # convert from J/mol to J/kmol
 
-        return ct.BlowersMaselRate(A, b, Ea, w)
+        rate = ct.BlowersMaselRate(A, b, Ea, w)
+        if A < 0:
+            rate.allow_negative_pre_exponential_factor = True
+        return rate
 
     def set_cantera_kinetics(self, ct_reaction, species_list):
         """

rmgpy/kinetics/falloff.pyx

@@ -236,7 +236,12 @@ cdef class Lindemann(PDepKineticsModel):
 
         high_rate = self.arrheniusHigh.to_cantera_kinetics(arrhenius_class=True)
         low_rate = self.arrheniusLow.to_cantera_kinetics(arrhenius_class=True)
-        return ct.LindemannRate(low=low_rate, high=high_rate)
+        rate = ct.LindemannRate()
+        if high_rate.pre_exponential_factor < 0 or low_rate.pre_exponential_factor < 0:
+            rate.allow_negative_pre_exponential_factor = True
+        rate.high_rate = high_rate
+        rate.low_rate = low_rate
+        return rate
 
 
 ################################################################################
@@ -415,4 +420,10 @@ cdef class Troe(PDepKineticsModel):
 
         high = self.arrheniusHigh.to_cantera_kinetics(arrhenius_class=True)
         low = self.arrheniusLow.to_cantera_kinetics(arrhenius_class=True)
-        return ct.TroeRate(high=high, low=low, falloff_coeffs=falloff)
+        rate = ct.TroeRate()
+        if high.pre_exponential_factor < 0 or low.pre_exponential_factor < 0:
+            rate.allow_negative_pre_exponential_factor = True
+        rate.high_rate = high
+        rate.low_rate = low
+        rate.falloff_coeffs = falloff
+        return rate

rmgpy/kinetics/surface.pyx

@@ -299,7 +299,10 @@ cdef class StickingCoefficient(KineticsModel):
         b = self._n.value_si
         E = self._Ea.value_si * 1000  # convert from J/mol to J/kmol
 
-        return ct.StickingArrheniusRate(A, b, E)
+        rate = ct.StickingArrheniusRate(A, b, E)
+        if A < 0:
+            rate.allow_negative_pre_exponential_factor = True
+        return rate
 
 
     def set_cantera_kinetics(self, ct_reaction, species_list):
@@ -654,7 +657,10 @@ cdef class SurfaceArrhenius(Arrhenius):
 
         b = self._n.value_si
         E = self._Ea.value_si * 1000  # convert from J/mol to J/kmol
-        return ct.InterfaceArrheniusRate(A, b, E)
+        rate = ct.InterfaceArrheniusRate(A, b, E)
+        if A < 0:
+            rate.allow_negative_pre_exponential_factor = True
+        return rate
 
     def set_cantera_kinetics(self, ct_reaction, species_list):
         """

rmgpy/yaml_cantera1.py

@@ -459,6 +459,30 @@ def reaction_to_dicts(obj, spcs):
             # spurious 'efficiencies' entry in input_data.
             # see https://github.com/Cantera/cantera/issues/2115
             reaction_data.pop("efficiencies", None)
+
+        # This next block is hopefully unnecessary, and probably slows us a little
+        # so maybe should be removed. For now it is here for safety.
+        for rate_key in (
+            "rate-constant",
+            "high-P-rate-constant",
+            "low-P-rate-constant",
+            "sticking-coefficient",
+        ):
+            if reaction_data.get(rate_key, {}).get("A", 0) < 0:
+                if not reaction_data.get("negative-A", False):
+                    logging.warning(
+                        "%s\n%s",
+                        "Reaction did not have the negative-A flag set to True.",
+                        yaml.dump(
+                            reaction_data,
+                            Dumper=Dumper,
+                            default_flow_style=False,
+                            sort_keys=False,
+                        ).rstrip(),
+                    )
+                reaction_data["negative-A"] = True
+                break
+
         # Convert any AnyMap objects to regular dicts before appending
         reaction_data = _convert_anymap_to_dict(reaction_data)
 

test/rmgpy/yaml_cantera1Test.py

@@ -43,6 +43,7 @@
 from rmgpy.transport import TransportData
 from rmgpy.kinetics import (
     Arrhenius,
+    Lindemann,
     PDepArrhenius,
     Troe,
     ThirdBody,
@@ -164,6 +165,13 @@ def test_reaction_to_dicts_arrhenius_equation_and_rate(self):
         assert np.isclose(d["rate-constant"]["b"], 0.5)
         assert np.isclose(d["rate-constant"]["Ea"], 10e6)  # 10 kJ/mol → 1e7 J/kmol
 
+    def test_reaction_to_dicts_negative_a_arrhenius(self):
+        """Negative Arrhenius A factors are marked for Cantera."""
+        kin = Arrhenius(A=(-1e13, "s^-1"), n=0.5, Ea=(10, "kJ/mol"), T0=(1, "K"))
+        rxn = Reaction(reactants=[self.h2], products=[self.h, self.h], kinetics=kin)
+        entries = reaction_to_dicts(rxn, self.all_gas)
+        assert entries[0]["negative-A"] is True
+
     def test_reaction_to_dicts_thirdbody(self):
         """ThirdBody: type is three-body, equation uses M, efficiencies map present."""
         kin = ThirdBody(
@@ -234,6 +242,24 @@ def test_reaction_to_dicts_troe(self):
         assert "efficiencies" in d
         assert np.isclose(d["efficiencies"]["Ar(3)"], 2.0)
 
+    def test_reaction_to_dicts_negative_a_falloff(self):
+        """Falloff rates with negative high- and low-pressure A factors are marked for Cantera."""
+        kin = Lindemann(
+            arrheniusHigh=Arrhenius(
+                A=(-1e14, "s^-1"), n=0, Ea=(10, "kJ/mol"), T0=(1, "K")
+            ),
+            arrheniusLow=Arrhenius(
+                A=(-1e20, "cm^3/(mol*s)"), n=0, Ea=(10, "kJ/mol"), T0=(1, "K")
+            ),
+            efficiencies={self.ar.molecule[0]: 2.0},
+        )
+        rxn = Reaction(reactants=[self.h], products=[self.h], kinetics=kin)
+        entries = reaction_to_dicts(rxn, self.all_gas)
+        d = entries[0]
+        assert d["negative-A"] is True
+        assert d["low-P-rate-constant"]["A"] < 0
+        assert np.isclose(d["efficiencies"]["Ar(3)"], 2.0)
+
     # ------------------------------------------------------------------
     # reaction_to_dicts — surface kinetics
     # ------------------------------------------------------------------
@@ -271,6 +297,19 @@ def test_reaction_to_dicts_sticking_coefficient(self):
         assert np.isclose(d["sticking-coefficient"]["A"], 0.1)
         assert np.isclose(d["sticking-coefficient"]["Ea"], 0.0)
 
+    def test_reaction_to_dicts_negative_a_sticking_coefficient(self):
+        """Negative sticking-coefficient A factors are marked for Cantera."""
+        kin = StickingCoefficient(
+            A=(-0.1, ""), n=0, Ea=(0, "kJ/mol"), T0=(1, "K")
+        )
+        rxn = Reaction(
+            reactants=[self.h2, self.x, self.x],
+            products=[self.hx, self.hx],
+            kinetics=kin,
+        )
+        entries = reaction_to_dicts(rxn, self.all_surface)
+        assert entries[0]["negative-A"] is True
+
     def test_reaction_to_dicts_coverage_dependence(self):
         """Coverage-dependent kinetics: coverage-dependencies block present with correct units.