WIP: RMG-Electrocat

.github/workflows/CI.yml

@@ -12,7 +12,7 @@ jobs:
     strategy:
       max-parallel: 5
     env: # update this if needed to match a pull request on the RMG-database
-      RMG_DATABASE_BRANCH: master
+      RMG_DATABASE_BRANCH: electrocat
     defaults:
       run:
         shell: bash -l {0}

rmgpy/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

rmgpy/constants.py

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

rmgpy/data/base.py

@@ -1350,8 +1350,8 @@ def is_molecule_forbidden(self, molecule):
                 raise NotImplementedError('Checking is only implemented for forbidden Groups, Molecule, and Species.')
 
         # Until we have more thermodynamic data of molecular ions we will forbid them
-        if molecule.get_net_charge() != 0:
-            return True
+        # if molecule.get_net_charge() != 0:
+        #     return True
 
         return False
 

rmgpy/data/kinetics/database.py

@@ -44,7 +44,8 @@
 from rmgpy.kinetics import Arrhenius, ArrheniusEP, ThirdBody, Lindemann, Troe, \
                            PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, \
                            Chebyshev, KineticsData, StickingCoefficient, \
-                           StickingCoefficientBEP, SurfaceArrhenius, SurfaceArrheniusBEP, ArrheniusBM
+                           StickingCoefficientBEP, SurfaceArrhenius, SurfaceArrheniusBEP, \
+                           ArrheniusBM, SurfaceChargeTransfer
 from rmgpy.molecule import Molecule, Group
 from rmgpy.reaction import Reaction, same_species_lists
 from rmgpy.species import Species
@@ -78,6 +79,7 @@ def __init__(self):
             'StickingCoefficientBEP': StickingCoefficientBEP,
             'SurfaceArrhenius': SurfaceArrhenius,
             'SurfaceArrheniusBEP': SurfaceArrheniusBEP,
+            'SurfaceChargeTransfer': SurfaceChargeTransfer,
             'R': constants.R,
             'ArrheniusBM': ArrheniusBM
         }

rmgpy/data/kinetics/depository.py

@@ -35,6 +35,7 @@
 
 from rmgpy.data.base import Database, Entry, DatabaseError
 from rmgpy.data.kinetics.common import save_entry
+from rmgpy.kinetics import SurfaceChargeTransfer, SurfaceArrheniusBEP
 from rmgpy.reaction import Reaction
 
 
@@ -60,7 +61,8 @@ def __init__(self,
                  pairs=None,
                  depository=None,
                  family=None,
-                 entry=None
+                 entry=None,
+                 electrons=None,
                  ):
         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
@@ -187,6 +190,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.'

rmgpy/data/kinetics/family.py

@@ -55,7 +55,7 @@
 from rmgpy.exceptions import ActionError, DatabaseError, InvalidActionError, KekulizationError, KineticsError, \
                              ForbiddenStructureException, UndeterminableKineticsError
 from rmgpy.kinetics import Arrhenius, SurfaceArrhenius, SurfaceArrheniusBEP, StickingCoefficient, \
-                           StickingCoefficientBEP, ArrheniusBM
+                           StickingCoefficientBEP, ArrheniusBM, SurfaceChargeTransfer
 from rmgpy.kinetics.uncertainties import RateUncertainty, rank_accuracy_map
 from rmgpy.molecule import Bond, GroupBond, Group, Molecule
 from rmgpy.molecule.atomtype import ATOMTYPES
@@ -101,6 +101,7 @@ def __init__(self,
                  estimator=None,
                  reverse=None,
                  is_forward=None,
+                 electrons=0,
                  ):
         Reaction.__init__(self,
                           index=index,
@@ -114,6 +115,7 @@ def __init__(self,
                           degeneracy=degeneracy,
                           pairs=pairs,
                           is_forward=is_forward,
+                          electrons=electrons
                           )
         self.family = family
         self.template = template
@@ -139,7 +141,8 @@ def __reduce__(self):
                                    self.template,
                                    self.estimator,
                                    self.reverse,
-                                   self.is_forward
+                                   self.is_forward,
+                                   self.electrons
                                    ))
 
     def __repr__(self):
@@ -161,6 +164,7 @@ def __repr__(self):
         if self.pairs is not None: string += 'pairs={0}, '.format(self.pairs)
         if self.family: string += "family='{}', ".format(self.family)
         if self.template: string += "template={}, ".format(self.template)
+        if self.electrons: string += "electrons={}, ".format(self.electrons)
         if self.comment != '': string += 'comment={0!r}, '.format(self.comment)
         string = string[:-2] + ')'
         return string
@@ -194,6 +198,7 @@ def copy(self):
         other.transition_state = deepcopy(self.transition_state)
         other.duplicate = self.duplicate
         other.pairs = deepcopy(self.pairs)
+        other.electrons = self.electrons
 
         # added for TemplateReaction information
         other.family = self.family
@@ -259,6 +264,10 @@ def get_reverse(self):
                 other.add_action(['GAIN_RADICAL', action[1], action[2]])
             elif action[0] == 'GAIN_RADICAL':
                 other.add_action(['LOSE_RADICAL', action[1], action[2]])
+            elif action[0] == 'GAIN_CHARGE':
+                other.add_action(['LOSE_CHARGE', action[1], action[2]])
+            elif action[0] == 'LOSE_CHARGE':
+                other.add_action(['GAIN_CHARGE', action[1], action[2]])
             elif action[0] == 'LOSE_PAIR':
                 other.add_action(['GAIN_PAIR', action[1], action[2]])
             elif action[0] == 'GAIN_PAIR':
@@ -382,7 +391,7 @@ def _apply(self, struct, forward, unique):
                     atom1.apply_action(['BREAK_BOND', label1, info, label2])
                     atom2.apply_action(['BREAK_BOND', label1, info, label2])
 
-            elif action[0] in ['LOSE_RADICAL', 'GAIN_RADICAL']:
+            elif action[0] in ['LOSE_RADICAL', 'GAIN_RADICAL', 'LOSE_CHARGE', 'GAIN_CHARGE']:
 
                 label, change = action[1:]
                 change = int(change)
@@ -400,6 +409,10 @@ def _apply(self, struct, forward, unique):
                             atom.apply_action(['GAIN_RADICAL', label, 1])
                         elif (action[0] == 'LOSE_RADICAL' and forward) or (action[0] == 'GAIN_RADICAL' and not forward):
                             atom.apply_action(['LOSE_RADICAL', label, 1])
+                        elif (action[0] == 'LOSE_CHARGE' and forward) or (action[0] == 'GAIN_CHARGE' and not forward):
+                            atom.apply_action(['LOSE_CHARGE', label, 1])
+                        elif (action[0] == 'GAIN_CHARGE' and forward) or (action[0] == 'LOSE_CHARGE' and not forward):
+                            atom.apply_action(['GAIN_CHARGE', label, 1])
 
             elif action[0] in ['LOSE_PAIR', 'GAIN_PAIR']:
 
@@ -720,6 +733,8 @@ def load(self, path, local_context=None, global_context=None, depository_labels=
         local_context['reactantNum'] = None
         local_context['productNum'] = None
         local_context['autoGenerated'] = False
+        local_context['allowChargedSpecies'] = False
+        local_context['electrons'] = 0
         self.groups = KineticsGroups(label='{0}/groups'.format(self.label))
         logging.debug("Loading kinetics family groups from {0}".format(os.path.join(path, 'groups.py')))
         Database.load(self.groups, os.path.join(path, 'groups.py'), local_context, global_context)
@@ -732,6 +747,8 @@ def load(self, path, local_context=None, global_context=None, depository_labels=
         self.product_num = local_context.get('productNum', None)
 
         self.auto_generated = local_context.get('autoGenerated', False)
+        self.allow_charged_species = local_context.get('allowChargedSpecies', False)
+        self.electrons = local_context.get('electrons', 0)
 
         if self.reactant_num:
             self.groups.reactant_num = self.reactant_num
@@ -835,7 +852,8 @@ def load_recipe(self, actions):
         for action in actions:
             action[0] = action[0].upper()
             valid_actions = [
-                'CHANGE_BOND', 'FORM_BOND', 'BREAK_BOND', 'GAIN_RADICAL', 'LOSE_RADICAL', 'GAIN_PAIR', 'LOSE_PAIR'
+                'CHANGE_BOND', 'FORM_BOND', 'BREAK_BOND', 'GAIN_RADICAL', 'LOSE_RADICAL', 
+                'GAIN_CHARGE', 'LOSE_CHARGE', 'GAIN_PAIR', 'LOSE_PAIR'
             ]
             if action[0] not in valid_actions:
                 raise InvalidActionError('Action {0} is not a recognized action. '
@@ -1022,6 +1040,12 @@ def save_groups(self, path):
         if self.auto_generated is not None:
             f.write('autoGenerated = {0}\n\n'.format(self.auto_generated))
 
+        if self.allow_charged_species:
+            f.write('allowChargedSpecies = {0}\n\n'.format(self.allow_charged_species))
+
+        if self.electrons != 0:
+            f.write('electrons = {0}\n\n'.format(self.electrons))
+
         # Write the recipe
         f.write('recipe(actions=[\n')
         for action in self.forward_recipe.actions:
@@ -1226,6 +1250,21 @@ def add_rules_from_training(self, thermo_database=None, train_indices=None):
                     Tmax=deepcopy(data.Tmax),
                     coverage_dependence=deepcopy(data.coverage_dependence),
                 )
+            elif isinstance(data, SurfaceChargeTransfer):
+                for reactant in entry.item.reactants:
+                    # Clear atom labels to avoid effects on thermo generation, ok because this is a deepcopy
+                    reactant_copy = reactant.copy(deep=True)
+                    reactant_copy.molecule[0].clear_labeled_atoms()
+                    reactant_copy.generate_resonance_structures()
+                    reactant.thermo = thermo_database.get_thermo_data(reactant_copy, training_set=True)
+                for product in entry.item.products:
+                    product_copy = product.copy(deep=True)
+                    product_copy.molecule[0].clear_labeled_atoms()
+                    product_copy.generate_resonance_structures()
+                    product.thermo = thermo_database.get_thermo_data(product_copy, training_set=True)
+                V = data.V0.value_si
+                dGrxn = entry.item._get_free_energy_of_charge_transfer_reaction(298,V)
+                data = data.to_surface_charge_transfer_bep(dGrxn,0.0)
             else:
                 raise NotImplementedError("Unexpected training kinetics type {} for {}".format(type(data), entry))
 
@@ -1254,7 +1293,9 @@ def add_rules_from_training(self, thermo_database=None, train_indices=None):
         for entry in reverse_entries:
             tentries[entry.index].item.is_forward = False
 
-            assert isinstance(entry.data, Arrhenius)
+            if not isinstance(entry.data, Arrhenius):
+                print(self.label)
+                assert False
             data = deepcopy(entry.data)
             data.change_t0(1)
             # Estimate the thermo for the reactants and products
@@ -1589,34 +1630,43 @@ def apply_recipe(self, reactant_structures, forward=True, unique=True, relabel_a
                 struc.update()
             reactant_net_charge += struc.get_net_charge()
 
+
+        is_molecule = True
         for struct in product_structures:
             # If product structures are Molecule objects, update their atom types
             # If product structures are Group objects and the reaction is in certain families
             # (families with charged substances), the charge of structures will be updated
             if isinstance(struct, Molecule):
+                struct.update_charge()
                 struct.update(sort_atoms=not self.save_order)
             elif isinstance(struct, Group):
+                is_molecule = False
                 struct.reset_ring_membership()
                 if label in ['1,2_insertion_co', 'r_addition_com', 'co_disproportionation',
                              'intra_no2_ono_conversion', 'lone_electron_pair_bond',
                              '1,2_nh3_elimination', '1,3_nh3_elimination']:
                     struct.update_charge()
             else:
                 raise TypeError('Expecting Molecule or Group object, not {0}'.format(struct.__class__.__name__))
-            product_net_charge += struc.get_net_charge()
-        if reactant_net_charge != product_net_charge:
+            product_net_charge += struct.get_net_charge()
+
+
+        if self.electrons < 0:
+            if forward:
+                reactant_net_charge += self.electrons
+            else: 
+                product_net_charge += self.electrons
+        elif self.electrons > 0:
+            if forward:
+                product_net_charge -= self.electrons
+            else: 
+                reactant_net_charge -= self.electrons
+
+        if reactant_net_charge != product_net_charge and is_molecule:
             logging.debug(
                 'The net charge of the reactants {0} differs from the net charge of the products {1} in reaction '
                 'family {2}. Not generating this reaction.'.format(reactant_net_charge, product_net_charge, self.label))
             return None
-        # The following check should be removed once RMG can process charged species
-        # This is applied only for :class:Molecule (not for :class:Group which is allowed to have a nonzero net charge)
-        if any([structure.get_net_charge() for structure in reactant_structures + product_structures]) \
-                and isinstance(struc, Molecule):
-            logging.debug(
-                'A net charged species was formed when reacting {0} to form {1} in reaction family {2}. Not generating '
-                'this reaction.'.format(reactant_net_charge, product_net_charge, self.label))
-            return None
 
         # If there are two product structures, place the one containing '*1' first
         if len(product_structures) == 2:
@@ -1762,8 +1812,17 @@ def _create_reaction(self, reactants, products, is_forward):
             reversible=self.reversible,
             family=self.label,
             is_forward=is_forward,
+            electrons = self.electrons
         )
 
+        if not self.allow_charged_species:
+            charged_species = [spc for spc in (reaction.reactants + reaction.products) if spc.get_net_charge() != 0]
+            if charged_species:
+                return None
+
+        if not reaction.is_balanced():
+            return None
+
         # Store the labeled atoms so we can recover them later
         # (e.g. for generating reaction pairs and templates)
         for key, species_list in zip(['reactants', 'products'], [reaction.reactants, reaction.products]):
@@ -1958,6 +2017,9 @@ def calculate_degeneracy(self, reaction):
         # 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 reaction.is_charge_transfer_reaction():
+            # Not implemented yet for charge transfer reactions
+            return 1
         reactants = reaction.reactants
         same_reactants = 0
         if len(reactants) == 2: