Added reaction.fixBarrierHeight method to convert from Evans-Polanyi …

…and raise Endothermic barriers.

This is similar to what I implemented in the RMG-Java.
The conversion from ArrheniusEP to Arrhenius is needed for the MEASURE
Pdep code to work. The correction of estimated barrier heights is useful
to fix things that were estimated to be negative when they probably
shouldn't be.

rmgpy/reaction.pxd

@@ -75,6 +75,8 @@ cdef class Reaction:
 
     cpdef double getRate(self, double T, double P, dict conc, double totalConc=?)
 
+    cpdef fixBarrierHeight(self)
+
     cpdef generateReverseRateCoefficient(self)
 
     cpdef numpy.ndarray calculateTSTRateCoefficients(self, numpy.ndarray Tlist, str tunneling=?)

rmgpy/reaction.py

@@ -42,10 +42,11 @@
 import cython
 import math
 import numpy
+import logging
 
 from quantity import constants
 from species import Species
-from kinetics import Arrhenius, KineticsData
+from kinetics import Arrhenius, KineticsData, ArrheniusEP
 
 ################################################################################
 
@@ -372,6 +373,29 @@ def getRate(self, T, P, conc, totalConc=-1.0):
         # Return rate
         return rateConstant * (forward - reverse / equilibriumConstant)
 
+    def fixBarrierHeight(self):
+        """
+        Turns the kinetics into Arrhenius (if they were ArrheniusEP)
+        and ensures the activation energy is at least the endothermicity
+        for endothermic reactions, and is not negative only as a result 
+        of using Evans Polanyi with an exothermic reaction.
+        """
+        cython.declare(H=cython.double, Ea=cython.double)
+        H = self.getEnthalpyOfReaction(298)
+        if isinstance(self.kinetics, ArrheniusEP):
+            Ea = self.kinetics.E0.value # temporarily using Ea to store the intrinsic barrier height E0
+            self.kinetics = self.kinetics.toArrhenius(H)
+            if Ea > 0 and self.kinetics.Ea.value < 0:
+                self.kinetics.comment += "Ea raised from {0:.1f} to 0 kJ/mol.".format(self.kinetics.Ea.value/1000)
+                logging.info("For reaction {1!s} Ea raised from {0:.1f} to 0 kJ/mol.".format(self.kinetics.Ea.value/1000, self))
+                self.kinetics.Ea.value = 0
+        if isinstance(self.kinetics, Arrhenius):
+            Ea = self.kinetics.Ea.value
+            if H > 0 and Ea < H:
+                self.kinetics.Ea.value = H
+                self.kinetics.comment += "Ea raised from {0:.1f} to {1:.1f} kJ/mol to match endothermicity of reaction.".format(Ea/1000,H/1000)
+                logging.info("For reaction {2!s}, Ea raised from {0:.1f} to {1:.1f} kJ/mol to match endothermicity of reaction.".format(Ea/1000, H/1000, self))
+
     def generateReverseRateCoefficient(self):
         """
         Generate and return a rate coefficient model for the reverse reaction. 

rmgpy/rmg/model.py

@@ -755,6 +755,8 @@ def makeNewReaction(self, forward, checkExisting=True):
         The forward direction is determined using the "is_reverse" attribute of the
         reaction's family.  If the reaction family is its own reverse, then it is
         made such that the forward reaction is exothermic at 298K.
+        
+        The forward reaction is appended to self.newReactionList if it is new.
         """
 
         # Determine the proper species objects for all reactants and products
@@ -910,7 +912,14 @@ def enlarge(self, newObject, database):
         logging.info('Generating thermodynamics for new species...')
         for spec in self.newSpeciesList:
             spec.generateThermoData(database)
-
+
+        # For new reactions, convert ArrheniusEP to Arrhenius, and fix barrier heights.
+        # self.newReactionList only contains *actually* new reactions, all in the forward direction.
+        for reaction in self.newReactionList:
+            if isinstance(reaction,TemplateReaction) or isinstance(reaction,DepositoryReaction): # i.e. not LibraryReaction
+                reaction.fixBarrierHeight() # also converts to Arrhenius.
+
+
         # Update unimolecular (pressure dependent) reaction networks
         if settings.pressureDependence:
             # Merge networks if necessary
@@ -960,7 +969,6 @@ def processNewReactions(self, newReactions, newSpecies, pdepNetwork=None):
         species or explored isomer `newSpecies` in network `pdepNetwork`.
         """
         for rxn in newReactions:
-
             rxn, isNew = self.makeNewReaction(rxn)
             if not isNew:
                 continue