FF till May 25

molecule/constants.pxd

@@ -25,4 +25,4 @@
 #                                                                             #
 ###############################################################################
 
-cdef double pi, Na, kB, R, h, hbar, c, e, m_e, m_p, m_n, amu, a0, E_h
+cdef double pi, Na, kB, R, h, hbar, c, e, m_e, m_p, m_n, amu, a0, E_h, F

molecule/constants.py

@@ -110,6 +110,12 @@
 #: :math:`\pi = 3.14159 \ldots`      
 pi = float(math.pi)
 
+#: Faradays Constant F in C/mol
+F = 96485.3321233100184
+
+#: Vacuum permittivity 
+epsilon_0 = 8.8541878128
+
 ################################################################################
 
 # Cython does not automatically place module-level variables into the module
@@ -130,4 +136,6 @@
     'm_n': m_n,
     'm_p': m_p,
     'pi': pi,
+    'F': F,
+    'epsilon_0': epsilon_0,
 })

molecule/data/base.py

@@ -42,6 +42,7 @@
 from molecule.data.reference import Reference, Article, Book, Thesis
 from molecule.exceptions import DatabaseError, InvalidAdjacencyListError
 from molecule.kinetics.uncertainties import RateUncertainty
+from molecule.kinetics.arrhenius import ArrheniusChargeTransfer, ArrheniusChargeTransferBM
 from molecule.molecule import Molecule, Group
 
 
@@ -228,6 +229,8 @@ def load(self, path, local_context=None, global_context=None):
         local_context['shortDesc'] = self.short_desc
         local_context['longDesc'] = self.long_desc
         local_context['RateUncertainty'] = RateUncertainty
+        local_context['ArrheniusChargeTransfer'] = ArrheniusChargeTransfer
+        local_context['ArrheniusChargeTransferBM'] = ArrheniusChargeTransferBM
         local_context['metal'] = self.metal
         local_context['site'] = self.site
         local_context['facet'] = self.facet
@@ -236,13 +239,17 @@ def load(self, path, local_context=None, global_context=None):
             local_context[key] = value
 
         # Process the file
-        f = open(path, 'r')
+        with open(path, 'r') as f:
+            content = f.read()
         try:
-            exec(f.read(), global_context, local_context)
-        except Exception:
-            logging.error('Error while reading database {0!r}.'.format(path))
+            exec(content, global_context, local_context)
+        except Exception as e:
+            logging.exception(f'Error while reading database file {path}.')
+            line_number = e.__traceback__.tb_next.tb_lineno
+            logging.error(f'Error occurred at or near line {line_number} of {path}.')
+            lines = content.splitlines()
+            logging.error(f'Line: {lines[line_number - 1]}')
             raise
-        f.close()
 
         # Extract the database metadata
         self.name = local_context['name']

molecule/data/kinetics/common.py

@@ -222,9 +222,55 @@ def independent_ids():
         for species in input_species:
             species.generate_resonance_structures(keep_isomorphic=True)
 
+def check_for_same_reactants(reactants):
+    """
+    Given a list reactants, check if the reactants are the same.
+    If they refer to the same memory address, then make a deep copy so they can be manipulated independently.
+
+    Returns a tuple containing the modified reactants list, and an integer containing the number of identical reactants in the reactants list. 
+
+    """
+
+    same_reactants = 0
+    if len(reactants) == 2:
+        if reactants[0] is reactants[1]:
+            reactants[1] = reactants[1].copy(deep=True)
+            same_reactants = 2
+        elif reactants[0].is_isomorphic(reactants[1]):
+            same_reactants = 2
+    elif len(reactants) == 3:
+        same_01 = reactants[0] is reactants[1]
+        same_02 = reactants[0] is reactants[2]
+        if same_01 and same_02:
+            same_reactants = 3
+            reactants[1] = reactants[1].copy(deep=True)
+            reactants[2] = reactants[2].copy(deep=True)
+        elif same_01:
+            same_reactants = 2
+            reactants[1] = reactants[1].copy(deep=True)
+        elif same_02:
+            same_reactants = 2
+            reactants[2] = reactants[2].copy(deep=True)
+        elif reactants[1] is reactants[2]:
+            same_reactants = 2
+            reactants[2] = reactants[2].copy(deep=True)
+        else:
+            same_01 = reactants[0].is_isomorphic(reactants[1])
+            same_02 = reactants[0].is_isomorphic(reactants[2])
+            if same_01 and same_02:
+                same_reactants = 3
+            elif same_01 or same_02:
+                same_reactants = 2
+            elif reactants[1].is_isomorphic(reactants[2]):
+                same_reactants = 2
+    elif len(reactants) > 3:
+        raise ValueError('Cannot check for duplicate reactants if provided number of reactants is greater than 3. ' 
+                         'Got: {} reactants'.format(len(reactants))) 
+
+    return reactants, same_reactants
 
 def find_degenerate_reactions(rxn_list, same_reactants=None, template=None, kinetics_database=None,
-                              kinetics_family=None, save_order=False):
+                              kinetics_family=None, save_order=False, resonance=True):
     """
     Given a list of Reaction objects, this method combines degenerate
     reactions and increments the reaction degeneracy value. For multiple
@@ -250,6 +296,7 @@ def find_degenerate_reactions(rxn_list, same_reactants=None, template=None, kine
         kinetics_database (KineticsDatabase, optional): provide a KineticsDatabase instance for calculating degeneracy
         kinetics_family (KineticsFamily, optional):     provide a KineticsFamily instance for calculating degeneracy
         save_order (bool, optional):                    reset atom order after performing atom isomorphism
+        resonance (bool, optional):                     whether to consider resonance when computing degeneracy 
 
     Returns:
         Reaction list with degenerate reactions combined with proper degeneracy values
@@ -340,7 +387,7 @@ def find_degenerate_reactions(rxn_list, same_reactants=None, template=None, kine
                 from molecule.data.rmg import get_db
                 family = get_db('kinetics').families[rxn.family]
             if not family.own_reverse:
-                rxn.degeneracy = family.calculate_degeneracy(rxn)
+                rxn.degeneracy = family.calculate_degeneracy(rxn, resonance=resonance)
 
     return rxn_list
 

molecule/data/kinetics/database.py

@@ -37,18 +37,22 @@
 import molecule.constants as constants
 from molecule.data.base import LogicNode
 from molecule.data.kinetics.common import ensure_species, generate_molecule_combos, \
-                                       find_degenerate_reactions, ensure_independent_atom_ids
+                                       find_degenerate_reactions, ensure_independent_atom_ids, \
+                                       check_for_same_reactants
 from molecule.data.kinetics.family import KineticsFamily
 from molecule.data.kinetics.library import LibraryReaction, KineticsLibrary
 from molecule.exceptions import DatabaseError
 from molecule.kinetics import Arrhenius, ArrheniusEP, ThirdBody, Lindemann, Troe, \
                            PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, \
                            Chebyshev, KineticsData, StickingCoefficient, \
-                           StickingCoefficientBEP, SurfaceArrhenius, SurfaceArrheniusBEP, ArrheniusBM
+                           StickingCoefficientBEP, SurfaceArrhenius, SurfaceArrheniusBEP, \
+                           ArrheniusBM, SurfaceChargeTransfer, KineticsModel, Marcus, \
+                           ArrheniusChargeTransfer, ArrheniusChargeTransferBM
+from molecule.kinetics.uncertainties import RateUncertainty
 from molecule.molecule import Molecule, Group
 from molecule.reaction import Reaction, same_species_lists
 from molecule.species import Species
-
+from molecule.data.solvation import SoluteData, SoluteTSData, SoluteTSDiffData
 
 ################################################################################
 
@@ -61,12 +65,14 @@ def __init__(self):
         self.recommended_families = {}
         self.families = {}
         self.libraries = {}
+        self.external_library_labels = {}
         self.library_order = []  # a list of tuples in the format ('library_label', LibraryType),
                                  # where LibraryType is set to either 'Reaction Library' or 'Seed'.
         self.local_context = {
             'KineticsData': KineticsData,
             'Arrhenius': Arrhenius,
             'ArrheniusEP': ArrheniusEP,
+            'ArrheniusChargeTransfer': ArrheniusChargeTransfer,
             'MultiArrhenius': MultiArrhenius,
             'MultiPDepArrhenius': MultiPDepArrhenius,
             'PDepArrhenius': PDepArrhenius,
@@ -78,8 +84,16 @@ def __init__(self):
             'StickingCoefficientBEP': StickingCoefficientBEP,
             'SurfaceArrhenius': SurfaceArrhenius,
             'SurfaceArrheniusBEP': SurfaceArrheniusBEP,
+            'SurfaceChargeTransfer': SurfaceChargeTransfer,
             'R': constants.R,
-            'ArrheniusBM': ArrheniusBM
+            'ArrheniusBM': ArrheniusBM,
+            'ArrheniusChargeTransferBM': ArrheniusChargeTransferBM,
+            'SoluteData': SoluteData,
+            'SoluteTSData': SoluteTSData,
+            'SoluteTSDiffData': SoluteTSDiffData,
+            'KineticsModel': KineticsModel,
+            'Marcus': Marcus,
+            'RateUncertainty': RateUncertainty,
         }
         self.global_context = {}
 
@@ -226,17 +240,18 @@ def load_libraries(self, path, libraries=None):
         The `path` points to the folder of kinetics libraries in the database,
         and the libraries should be in files like :file:`<path>/<library>.py`.
         """
-
+        self.external_library_labels = dict()
         if libraries is not None:
             for library_name in libraries:
                 library_file = os.path.join(path, library_name, 'reactions.py')
                 if os.path.exists(library_name):
                     library_file = os.path.join(library_name, 'reactions.py')
-                    short_library_name = os.path.split(library_name)[-1]
+                    short_library_name = os.path.basename(library_name.rstrip(os.path.sep))
                     logging.info(f'Loading kinetics library {short_library_name} from {library_name}...')
                     library = KineticsLibrary(label=short_library_name)
                     library.load(library_file, self.local_context, self.global_context)
                     self.libraries[library.label] = library
+                    self.external_library_labels[library_name] = library.label
                 elif os.path.exists(library_file):
                     logging.info(f'Loading kinetics library {library_name} from {library_file}...')
                     library = KineticsLibrary(label=library_name)
@@ -251,10 +266,18 @@ def load_libraries(self, path, libraries=None):
             self.library_order = []
             for (root, dirs, files) in os.walk(os.path.join(path)):
                 for f in files:
-                    name, ext = os.path.splitext(f)
-                    if ext.lower() == '.py':
+                    if f.lower() == 'reactions.py':
                         library_file = os.path.join(root, f)
-                        label = os.path.dirname(library_file)[len(path) + 1:]
+                        dirname = os.path.dirname(library_file)
+                        if dirname == path:
+                            label = os.path.basename(dirname)
+                        else:
+                            label = os.path.relpath(dirname, path)
+
+                        if not label:
+                            logging.warning(f"Empty label for {library_file}. Using 'default'.")
+                            label = "default"
+
                         logging.info(f'Loading kinetics library {label} from {library_file}...')
                         library = KineticsLibrary(label=label)
                         try:
@@ -428,7 +451,7 @@ def generate_reactions(self, reactants, products=None, only_families=None, reson
         if only_families is None:
             reaction_list.extend(self.generate_reactions_from_libraries(reactants, products))
         reaction_list.extend(self.generate_reactions_from_families(reactants, products,
-                                                                   only_families=None, resonance=resonance))
+                                                                   only_families=only_families, resonance=resonance))
         return reaction_list
 
     def generate_reactions_from_libraries(self, reactants, products=None):
@@ -488,43 +511,10 @@ def generate_reactions_from_families(self, reactants, products=None, only_famili
         Returns:
             List of reactions containing Species objects with the specified reactants and products.
         """
-        # Check if the reactants are the same
-        # If they refer to the same memory address, then make a deep copy so
-        # they can be manipulated independently
         if isinstance(reactants, tuple):
             reactants = list(reactants)
-        same_reactants = 0
-        if len(reactants) == 2:
-            if reactants[0] is reactants[1]:
-                reactants[1] = reactants[1].copy(deep=True)
-                same_reactants = 2
-            elif reactants[0].is_isomorphic(reactants[1]):
-                same_reactants = 2
-        elif len(reactants) == 3:
-            same_01 = reactants[0] is reactants[1]
-            same_02 = reactants[0] is reactants[2]
-            if same_01 and same_02:
-                same_reactants = 3
-                reactants[1] = reactants[1].copy(deep=True)
-                reactants[2] = reactants[2].copy(deep=True)
-            elif same_01:
-                same_reactants = 2
-                reactants[1] = reactants[1].copy(deep=True)
-            elif same_02:
-                same_reactants = 2
-                reactants[2] = reactants[2].copy(deep=True)
-            elif reactants[1] is reactants[2]:
-                same_reactants = 2
-                reactants[2] = reactants[2].copy(deep=True)
-            else:
-                same_01 = reactants[0].is_isomorphic(reactants[1])
-                same_02 = reactants[0].is_isomorphic(reactants[2])
-                if same_01 and same_02:
-                    same_reactants = 3
-                elif same_01 or same_02:
-                    same_reactants = 2
-                elif reactants[1].is_isomorphic(reactants[2]):
-                    same_reactants = 2
+
+        reactants, same_reactants = check_for_same_reactants(reactants)
 
         # Label reactant atoms for proper degeneracy calculation (cannot be in tuple)
         ensure_independent_atom_ids(reactants, resonance=resonance)
@@ -537,7 +527,8 @@ def generate_reactions_from_families(self, reactants, products=None, only_famili
                                                       prod_resonance=resonance))
 
         # Calculate reaction degeneracy
-        reaction_list = find_degenerate_reactions(reaction_list, same_reactants, kinetics_database=self)
+        reaction_list = find_degenerate_reactions(reaction_list, same_reactants, kinetics_database=self,
+                                                  resonance=resonance)
         # Add reverse attribute to families with ownReverse
         to_delete = []
         for i, rxn in enumerate(reaction_list):
@@ -658,7 +649,7 @@ def get_forward_reaction_for_family_entry(self, entry, family, thermo_database):
             elif len(reverse) == 1 and len(forward) == 0:
                 # The reaction is in the reverse direction
                 # First fit Arrhenius kinetics in that direction
-                T_data = 1000.0 / np.arange(0.5, 3.301, 0.1, np.float64)
+                T_data = 1000.0 / np.arange(0.5, 3.301, 0.1, float)
                 k_data = np.zeros_like(T_data)
                 for i in range(T_data.shape[0]):
                     k_data[i] = entry.data.get_rate_coefficient(T_data[i]) / reaction.get_equilibrium_constant(T_data[i])

molecule/data/kinetics/depository.py

@@ -35,6 +35,7 @@
 
 from molecule.data.base import Database, Entry, DatabaseError
 from molecule.data.kinetics.common import save_entry
+from molecule.kinetics import SurfaceChargeTransfer, SurfaceArrheniusBEP
 from molecule.reaction import Reaction
 
 
@@ -60,7 +61,8 @@ def __init__(self,
                  pairs=None,
                  depository=None,
                  family=None,
-                 entry=None
+                 entry=None,
+                 electrons=0,
                  ):
         Reaction.__init__(self,
                           index=index,
@@ -72,7 +74,8 @@ def __init__(self,
                           transition_state=transition_state,
                           duplicate=duplicate,
                           degeneracy=degeneracy,
-                          pairs=pairs
+                          pairs=pairs,
+                          electrons=electrons,
                           )
         self.depository = depository
         self.family = family
@@ -104,12 +107,34 @@ def get_source(self):
         """
         return self.depository.label
 
+    def apply_solvent_correction(self, solvent):
+        """
+        apply kinetic solvent correction
+        """
+        from molecule.data.rmg import get_db
+        solvation_database = get_db('solvation')
+        solvent_data = solvation_database.get_solvent_data(solvent)
+        solute_data = self.kinetics.solute
+        correction = solvation_database.get_solvation_correction(solute_data, solvent_data)
+        dHR = 0.0
+        dSR = 0.0
+        for spc in self.reactants:
+            spc_solute_data = solvation_database.get_solute_data(spc)
+            spc_correction = solvation_database.get_solvation_correction(spc_solute_data, solvent_data)
+            dHR += spc_correction.enthalpy
+            dSR += spc_correction.entropy
+
+        dH = correction.enthalpy-dHR
+        dA = np.exp((correction.entropy-dSR)/constants.R)
+        self.kinetics.Ea.value_si += dH
+        self.kinetics.A.value_si *= dA
+        self.kinetics.comment += "\nsolvation correction raised barrier by {0} kcal/mol and prefactor by factor of {1}".format(dH/4184.0,dA)
 
 ################################################################################
 
 class KineticsDepository(Database):
     """
-    A class for working with an RMG kinetics depository. Each depository 
+    A class for working with an RMG kinetics depository. Each depository
     corresponds to a reaction family (a :class:`KineticsFamily` object). Each
     entry in a kinetics depository involves a reaction defined either by a
     real reactant and product species (as in a kinetics library).
@@ -187,6 +212,9 @@ def load(self, path, local_context=None, global_context=None):
                                         ''.format(product, self.label))
                 # Same comment about molecule vs species objects as above.
                 rxn.products.append(species_dict[product])
+
+            if isinstance(entry.data, (SurfaceChargeTransfer, SurfaceArrheniusBEP)):
+                rxn.electrons = entry.data.electrons.value
 
             if not rxn.is_balanced():
                 raise DatabaseError('Reaction {0} in kinetics depository {1} was not balanced! Please reformulate.'