diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 7a4095a206..1bcce1d4c5 100644 --- a/.github/workflows/CI.yml +++ b/.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} diff --git a/rmgpy/constants.pxd b/rmgpy/constants.pxd index 292ec667ca..491c98a708 100644 --- a/rmgpy/constants.pxd +++ b/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 diff --git a/rmgpy/constants.py b/rmgpy/constants.py index 9e7d0f47f5..212d5c3fca 100644 --- a/rmgpy/constants.py +++ b/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 }) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index 2b617a6b27..a381425389 100644 --- a/rmgpy/data/base.py +++ b/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 diff --git a/rmgpy/data/kinetics/database.py b/rmgpy/data/kinetics/database.py index d37e0de432..bccad66598 100644 --- a/rmgpy/data/kinetics/database.py +++ b/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 } diff --git a/rmgpy/data/kinetics/depository.py b/rmgpy/data/kinetics/depository.py index ad24c7db36..3dbd387c32 100644 --- a/rmgpy/data/kinetics/depository.py +++ b/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.' diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 36d024a1bb..e55b7d1143 100644 --- a/rmgpy/data/kinetics/family.py +++ b/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,13 +1630,17 @@ 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', @@ -1603,20 +1648,25 @@ def apply_recipe(self, reactant_structures, forward=True, unique=True, relabel_a 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: diff --git a/rmgpy/data/kinetics/familyTest.py b/rmgpy/data/kinetics/familyTest.py index af39d28659..6eb2841558 100644 --- a/rmgpy/data/kinetics/familyTest.py +++ b/rmgpy/data/kinetics/familyTest.py @@ -43,6 +43,7 @@ from rmgpy.data.thermo import ThermoDatabase from rmgpy.molecule import Molecule from rmgpy.species import Species +from rmgpy.reaction import Reaction ################################################### @@ -68,7 +69,8 @@ def setUpClass(cls): 'Intra_R_Add_Exo_scission', 'intra_substitutionS_isomerization', 'R_Addition_COm', - 'R_Recombination' + 'R_Recombination', + 'Surface_Proton_Electron_Reduction_Alpha', ], ) cls.family = cls.database.families['intra_H_migration'] @@ -623,6 +625,51 @@ def test_r_addition_com(self): self.assertEqual(len(products), 1) self.assertTrue(expected_products[0].is_isomorphic(products[0])) + def test_surface_proton_electron_reduction_alpha(self): + """ + Test that the Surface_Proton_Electron_Reduction_Alpha family can successfully match the reaction and returns properly product structures. + """ + family = self.database.families['Surface_Proton_Electron_Reduction_Alpha'] + m_proton = Molecule().from_smiles("[H+]") + m_x = Molecule().from_adjacency_list("1 X u0 p0") + m_ch2x = Molecule().from_adjacency_list( + """ + 1 C u0 p0 c0 {2,S} {3,S} {4,D} + 2 H u0 p0 c0 {1,S} + 3 H u0 p0 c0 {1,S} + 4 X u0 p0 c0 {1,D} + """ + ) + m_ch3x = Molecule().from_adjacency_list( + """ + 1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} + 2 H u0 p0 c0 {1,S} + 3 H u0 p0 c0 {1,S} + 4 H u0 p0 c0 {1,S} + 5 X u0 p0 c0 {1,S} + """ + ) + + reactants = [m_proton,m_ch2x] + expected_products = [m_ch3x] + + labeled_rxn = Reaction(reactants=reactants, products=expected_products) + family.add_atom_labels_for_reaction(labeled_rxn) + prods = family.apply_recipe([m.molecule[0] for m in labeled_rxn.reactants]) + self.assertTrue(expected_products[0].is_isomorphic(prods[0])) + + self.assertEqual(len(prods), 1) + self.assertTrue(expected_products[0].is_isomorphic(prods[0])) + reacts = family.apply_recipe(prods, forward=False) + self.assertEqual(len(reacts), 2) + + prods = [Species(molecule=[p]) for p in prods] + reacts = [Species(molecule=[r]) for r in reacts] + + fam_rxn = Reaction(reactants=reacts,products=prods) + + self.assertTrue(fam_rxn.is_isomorphic(labeled_rxn)) + def test_save_family(self): """ diff --git a/rmgpy/data/kinetics/rules.py b/rmgpy/data/kinetics/rules.py index 5a9f05bc4a..b8132b682d 100644 --- a/rmgpy/data/kinetics/rules.py +++ b/rmgpy/data/kinetics/rules.py @@ -44,7 +44,8 @@ from rmgpy.data.base import Database, Entry, get_all_combinations from rmgpy.data.kinetics.common import save_entry from rmgpy.exceptions import KineticsError, DatabaseError -from rmgpy.kinetics import ArrheniusEP, Arrhenius, StickingCoefficientBEP, SurfaceArrheniusBEP +from rmgpy.kinetics import ArrheniusEP, Arrhenius, StickingCoefficientBEP, SurfaceArrheniusBEP, \ + SurfaceChargeTransfer, SurfaceChargeTransferBEP from rmgpy.quantity import Quantity, ScalarQuantity from rmgpy.reaction import Reaction @@ -565,12 +566,23 @@ def _get_average_kinetics(self, kinetics_list): n = 0.0 E0 = 0.0 alpha = 0.0 - count = len(kinetics_list) + electrons = None + V0 = None + count = 0 for kinetics in kinetics_list: + if isinstance(kinetics, SurfaceChargeTransfer): + continue + count += 1 logA += math.log10(kinetics.A.value_si) n += kinetics.n.value_si alpha += kinetics.alpha.value_si E0 += kinetics.E0.value_si + if isinstance(kinetics, SurfaceChargeTransferBEP): + if electrons is None: + electrons = kinetics.electrons.value_si + if V0 is None: + V0 = kinetics.V0.value_si + logA /= count n /= count alpha /= count @@ -595,15 +607,25 @@ def _get_average_kinetics(self, kinetics_list): else: raise Exception('Invalid units {0} for averaging kinetics.'.format(Aunits)) - if type(kinetics) not in [ArrheniusEP, SurfaceArrheniusBEP, StickingCoefficientBEP]: + if type(kinetics) not in [ArrheniusEP, SurfaceArrheniusBEP, StickingCoefficientBEP, SurfaceChargeTransferBEP]: raise Exception('Invalid kinetics type {0!r} for {1!r}.'.format(type(kinetics), self)) - averaged_kinetics = type(kinetics)( - A=(10 ** logA, Aunits), - n=n, - alpha=alpha, - E0=(E0 * 0.001, "kJ/mol"), - ) + if isinstance(kinetics, SurfaceChargeTransferBEP): + averaged_kinetics = SurfaceChargeTransferBEP( + A=(10 ** logA, Aunits), + n=n, + electrons=electrons, + alpha=alpha, + V0=(V0,'V'), + E0=(E0 * 0.001, "kJ/mol"), + ) + else: + averaged_kinetics = type(kinetics)( + A=(10 ** logA, Aunits), + n=n, + alpha=alpha, + E0=(E0 * 0.001, "kJ/mol"), + ) return averaged_kinetics def estimate_kinetics(self, template, degeneracy=1): diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py index 9a981625ed..4dff418da9 100644 --- a/rmgpy/data/solvation.py +++ b/rmgpy/data/solvation.py @@ -1260,6 +1260,7 @@ def estimate_radical_solute_data_via_hbi(self, molecule, stable_solute_data_esti saturated_struct.remove_bond(bond) saturated_struct.remove_atom(H) atom.increment_radical() + saturated_struct.update_charge() # we need to update charges before updating lone pairs saturated_struct.update() try: self._add_group_solute_data(solute_data, self.groups['radical'], saturated_struct, {'*': atom}) diff --git a/rmgpy/kinetics/__init__.py b/rmgpy/kinetics/__init__.py index a9dc32e49b..2b895ca956 100644 --- a/rmgpy/kinetics/__init__.py +++ b/rmgpy/kinetics/__init__.py @@ -35,4 +35,5 @@ from rmgpy.kinetics.kineticsdata import KineticsData, PDepKineticsData from rmgpy.kinetics.tunneling import Wigner, Eckart from rmgpy.kinetics.surface import SurfaceArrhenius, SurfaceArrheniusBEP, \ - StickingCoefficient, StickingCoefficientBEP + StickingCoefficient, StickingCoefficientBEP, \ + SurfaceChargeTransfer, SurfaceChargeTransferBEP diff --git a/rmgpy/kinetics/surface.pxd b/rmgpy/kinetics/surface.pxd index 8e37ac54b7..9f60039e26 100644 --- a/rmgpy/kinetics/surface.pxd +++ b/rmgpy/kinetics/surface.pxd @@ -34,7 +34,7 @@ from rmgpy.quantity cimport ScalarQuantity, ArrayQuantity ################################################################################ cdef class StickingCoefficient(KineticsModel): - + cdef public ScalarQuantity _A cdef public ScalarQuantity _n cdef public ScalarQuantity _Ea @@ -48,7 +48,7 @@ cdef class StickingCoefficient(KineticsModel): cpdef fit_to_data(self, np.ndarray Tlist, np.ndarray klist, str kunits, double T0=?, np.ndarray weights=?, bint three_params=?) cpdef bint is_identical_to(self, KineticsModel other_kinetics) except -2 - + cpdef change_rate(self, double factor) cdef class StickingCoefficientBEP(KineticsModel): @@ -67,9 +67,62 @@ cdef class StickingCoefficientBEP(KineticsModel): ################################################################################ cdef class SurfaceArrhenius(Arrhenius): cdef public dict _coverage_dependence - pass + cpdef SurfaceChargeTransfer to_surface_charge_transfer(self, double V0, double electrons=?) + ################################################################################ cdef class SurfaceArrheniusBEP(ArrheniusEP): cdef public dict _coverage_dependence pass +################################################################################ +cdef class SurfaceChargeTransfer(KineticsModel): + + cdef public ScalarQuantity _A + cdef public ScalarQuantity _n + cdef public ScalarQuantity _Ea + cdef public ScalarQuantity _T0 + cdef public ScalarQuantity _V0 + cdef public ScalarQuantity _alpha + cdef public ScalarQuantity _electrons + + cpdef double get_activation_energy_from_potential(self, double V=?, bint non_negative=?) + + cpdef double get_rate_coefficient(self, double T, double V=?) except -1 + + cpdef change_rate(self, double factor) + + cpdef change_t0(self, double T0) + + cpdef change_v0(self, double V0) + + cpdef fit_to_data(self, np.ndarray Tlist, np.ndarray klist, str kunits, double T0=?, np.ndarray weights=?, bint three_params=?) + + cpdef bint is_identical_to(self, KineticsModel other_kinetics) except -2 + + cpdef SurfaceArrhenius to_surface_arrhenius(self) + + cpdef SurfaceChargeTransferBEP to_surface_charge_transfer_bep(self, double dGrxn, double V0=?) + +################################################################################ +cdef class SurfaceChargeTransferBEP(KineticsModel): + + cdef public ScalarQuantity _A + cdef public ScalarQuantity _n + cdef public ScalarQuantity _E0 + cdef public ScalarQuantity _V0 + cdef public ScalarQuantity _alpha + cdef public ScalarQuantity _electrons + + cpdef change_v0(self, double V0) + + cpdef double get_activation_energy(self, double dGrxn) except -1 + + cpdef double get_activation_energy_from_potential(self, double V, double dGrxn) except -1 + + cpdef double get_rate_coefficient_from_potential(self, double T, double V, double dGrxn) except -1 + + cpdef SurfaceChargeTransfer to_surface_charge_transfer(self, double dGrxn) + + cpdef bint is_identical_to(self, KineticsModel other_kinetics) except -2 + + cpdef change_rate(self, double factor) diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx index 2c108e9204..2d2e0c06bb 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -511,6 +511,25 @@ cdef class SurfaceArrhenius(Arrhenius): return (SurfaceArrhenius, (self.A, self.n, self.Ea, self.T0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.coverage_dependence, self.uncertainty, self.comment)) + cpdef SurfaceChargeTransfer to_surface_charge_transfer(self, double V0, double electrons=-1): + """ + Return an :class:`SurfaceChargeTransfer` instance of the kinetics model with reference + potential `V0` in Volts and electron stochiometric coeff `electrons` + """ + return SurfaceChargeTransfer( + A=self.A, + n=self.n, + electrons=electrons, + Ea=self.Ea, + V0=(V0,'V'), + T0=(1, "K"), + Tmin=self.Tmin, + Tmax=self.Tmax, + uncertainty = self.uncertainty, + comment=self.comment, + ) + + ################################################################################ cdef class SurfaceArrheniusBEP(ArrheniusEP): @@ -625,3 +644,504 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): coverage_dependence=self.coverage_dependence, comment=self.comment, ) + +################################################################################ + +cdef class SurfaceChargeTransfer(KineticsModel): + + """ + A kinetics model for surface charge transfer reactions + + It is very similar to the :class:`SurfaceArrhenius`, but the Ea is potential-dependent + + + The attributes are: + + =============== ============================================================= + Attribute Description + =============== ============================================================= + `A` The preexponential factor + `T0` The reference temperature + `n` The temperature exponent + `Ea` The activation energy + `electrons` The stochiometry coeff for electrons (negative if reactant, positive if product) + `V0` The reference potential + `alpha` The charge transfer coefficient + `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined + `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined + `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined + `comment` Information about the model (e.g. its source) + =============== ============================================================= + + """ + + def __init__(self, A=None, n=0.0, Ea=None, V0=None, alpha=0.5, electrons=-1, T0=(1.0, "K"), Tmin=None, Tmax=None, + Pmin=None, Pmax=None, uncertainty=None, comment=''): + + KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, uncertainty=uncertainty, + comment=comment) + + self.alpha = alpha + self.A = A + self.n = n + self.Ea = Ea + self.T0 = T0 + self.electrons = electrons + self.V0 = V0 + + def __repr__(self): + """ + Return a string representation that can be used to reconstruct the + Arrhenius object. + """ + string = 'SurfaceChargeTransfer(A={0!r}, n={1!r}, Ea={2!r}, V0={3!r}, alpha={4!r}, electrons={5!r}, T0={6!r}'.format( + self.A, self.n, self.Ea, self.V0, self.alpha, self.electrons, self.T0) + if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin) + if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) + if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin) + if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax) + if self.uncertainty: string += ', uncertainty={0!r}'.format(self.uncertainty) + if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) + string += ')' + return string + + def __reduce__(self): + """ + A helper function used when pickling a SurfaceChargeTransfer object. + """ + return (SurfaceChargeTransfer, (self.A, self.n, self.Ea, self.V0, self.alpha, self.electrons, self.T0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, + self.uncertainty, self.comment)) + + property A: + """The preexponential factor.""" + def __get__(self): + return self._A + def __set__(self, value): + self._A = quantity.SurfaceRateCoefficient(value) + + property n: + """The temperature exponent.""" + def __get__(self): + return self._n + def __set__(self, value): + self._n = quantity.Dimensionless(value) + + property Ea: + """The activation energy.""" + def __get__(self): + return self._Ea + def __set__(self, value): + self._Ea = quantity.Energy(value) + + property T0: + """The reference temperature.""" + def __get__(self): + return self._T0 + def __set__(self, value): + self._T0 = quantity.Temperature(value) + + property V0: + """The reference potential.""" + def __get__(self): + return self._V0 + def __set__(self, value): + self._V0 = quantity.Potential(value) + + property electrons: + """The number of electrons transferred.""" + def __get__(self): + return self._electrons + def __set__(self, value): + self._electrons = quantity.Dimensionless(value) + + property alpha: + """The charge transfer coefficient.""" + def __get__(self): + return self._alpha + def __set__(self, value): + self._alpha = quantity.Dimensionless(value) + + cpdef double get_activation_energy_from_potential(self, double V=0.0, bint non_negative=True): + """ + Return the effective activation energy (in J/mol) at specificed potential (in Volts). + """ + cdef double electrons, alpha, Ea, V0 + + electrons = self._electrons.value_si + alpha = self._alpha.value_si + Ea = self._Ea.value_si + V0 = self._V0.value_si + + Ea -= alpha * electrons * constants.F * (V-V0) + + if non_negative is True: + if Ea < 0: + Ea = 0.0 + + return Ea + + cpdef double get_rate_coefficient(self, double T, double V=0.0) except -1: + """ + Return the rate coefficient in the appropriate combination of m^2, + mol, and s at temperature `T` in K. + """ + cdef double A, n, V0, T0, Ea + + A = self._A.value_si + n = self._n.value_si + V0 = self._V0.value_si + T0 = self._T0.value_si + + if V != V0: + Ea = self.get_activation_energy_from_potential(V) + else: + Ea = self._Ea.value_si + + return A * (T / T0) ** n * exp(-Ea / (constants.R * T)) + + cpdef change_t0(self, double T0): + """ + Changes the reference temperature used in the exponent to `T0` in K, + and adjusts the preexponential factor accordingly. + """ + self._A.value_si /= (self._T0.value_si / T0) ** self._n.value_si + self._T0.value_si = T0 + + cpdef change_v0(self, double V0): + """ + Changes the reference potential to `V0` in volts, and adjusts the + activation energy `Ea` accordingly. + """ + + self._Ea.value_si = self.get_activation_energy_from_potential(V0) + self._V0.value_si = V0 + + cpdef fit_to_data(self, np.ndarray Tlist, np.ndarray klist, str kunits, double T0=1, + np.ndarray weights=None, bint three_params=False): + """ + Fit the Arrhenius parameters to a set of rate coefficient data `klist` + in units of `kunits` corresponding to a set of temperatures `Tlist` in + K. A linear least-squares fit is used, which guarantees that the + resulting parameters provide the best possible approximation to the + data. + """ + import scipy.stats + if not all(np.isfinite(klist)): + raise ValueError("Rates must all be finite, not inf or NaN") + if any(klist<0): + if not all(klist<0): + raise ValueError("Rates must all be positive or all be negative.") + rate_sign_multiplier = -1 + klist = -1 * klist + else: + rate_sign_multiplier = 1 + + assert len(Tlist) == len(klist), "length of temperatures and rates must be the same" + if len(Tlist) < 3 + three_params: + raise KineticsError('Not enough degrees of freedom to fit this Arrhenius expression') + if three_params: + A = np.zeros((len(Tlist), 3), np.float64) + A[:, 0] = np.ones_like(Tlist) + A[:, 1] = np.log(Tlist / T0) + A[:, 2] = -1.0 / constants.R / Tlist + else: + A = np.zeros((len(Tlist), 2), np.float64) + A[:, 0] = np.ones_like(Tlist) + A[:, 1] = -1.0 / constants.R / Tlist + b = np.log(klist) + if weights is not None: + for n in range(b.size): + A[n, :] *= weights[n] + b[n] *= weights[n] + x, residues, rank, s = np.linalg.lstsq(A, b, rcond=RCOND) + + # Determine covarianace matrix to obtain parameter uncertainties + count = klist.size + cov = residues[0] / (count - 3) * np.linalg.inv(np.dot(A.T, A)) + t = scipy.stats.t.ppf(0.975, count - 3) + + if not three_params: + x = np.array([x[0], 0, x[1]]) + cov = np.array([[cov[0, 0], 0, cov[0, 1]], [0, 0, 0], [cov[1, 0], 0, cov[1, 1]]]) + + self.A = (rate_sign_multiplier * exp(x[0]), kunits) + self.n = x[1] + self.Ea = (x[2] * 0.001, "kJ/mol") + self.T0 = (T0, "K") + self.Tmin = (np.min(Tlist), "K") + self.Tmax = (np.max(Tlist), "K") + self.comment = 'Fitted to {0:d} data points; dA = *|/ {1:g}, dn = +|- {2:g}, dEa = +|- {3:g} kJ/mol'.format( + len(Tlist), + exp(sqrt(cov[0, 0])), + sqrt(cov[1, 1]), + sqrt(cov[2, 2]) * 0.001, + ) + + return self + + cpdef bint is_identical_to(self, KineticsModel other_kinetics) except -2: + """ + Returns ``True`` if kinetics matches that of another kinetics model. Must match temperature + and pressure range of kinetics model, as well as parameters: A, n, Ea, T0. (Shouldn't have pressure + range if it's Arrhenius.) Otherwise returns ``False``. + """ + if not isinstance(other_kinetics, SurfaceChargeTransfer): + return False + if not KineticsModel.is_identical_to(self, other_kinetics): + return False + if (not self.A.equals(other_kinetics.A) or not self.n.equals(other_kinetics.n) + or not self.Ea.equals(other_kinetics.Ea) or not self.T0.equals(other_kinetics.T0) + or not self.alpha.equals(other_kinetics.alpha) or not self.electrons.equals(other_kinetics.electrons) + or not self.V0.equals(other_kinetics.V0)): + return False + + return True + + cpdef change_rate(self, double factor): + """ + Changes A factor in Arrhenius expression by multiplying it by a ``factor``. + """ + self._A.value_si *= factor + + cpdef SurfaceArrhenius to_surface_arrhenius(self): + """ + Return an :class:`SurfaceArrhenius` instance of the kinetics model + """ + return SurfaceArrhenius( + A=self.A, + n=self.n, + Ea=self.Ea, + T0=(1, "K"), + Tmin=self.Tmin, + Tmax=self.Tmax, + uncertainty = self.uncertainty, + comment=self.comment, + ) + + cpdef SurfaceChargeTransferBEP to_surface_charge_transfer_bep(self, double dGrxn, double V0=0.0): + """ + Converts an SurfaceChargeTransfer object to SurfaceChargeTransferBEP + """ + cdef double E0 + + self.change_t0(1) + self.change_v0(V0) + + E0 = self.Ea.value_si - self._alpha.value_si * dGrxn + if E0 < 0: + E0 = 0.0 + + aep = SurfaceChargeTransferBEP( + A=self.A, + electrons=self.electrons, + n=self.n, + alpha=self.alpha, + V0=self.V0, + E0=(E0, 'J/mol'), + Tmin=self.Tmin, + Tmax=self.Tmax, + Pmin=self.Pmin, + Pmax=self.Pmax, + uncertainty=self.uncertainty, + comment=self.comment) + return aep + +cdef class SurfaceChargeTransferBEP(KineticsModel): + """ + A kinetics model based on the (modified) Arrhenius equation, using the + Evans-Polanyi equation to determine the activation energy. The attributes + are: + + =============== ============================================================= + Attribute Description + =============== ============================================================= + `A` The preexponential factor + `n` The temperature exponent + `E0` The activation energy at equilibiurm + `electrons` The stochiometry coeff for electrons (negative if reactant, positive if product) + `V0` The reference potential + `alpha` The charge transfer coefficient + `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined + `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined + `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined + `comment` Information about the model (e.g. its source) + =============== ============================================================= + + """ + + def __init__(self, A=None, n=0.0, E0=None, V0=(0.0,'V'), alpha=0.5, electrons=-1, Tmin=None, Tmax=None, + Pmin=None, Pmax=None, uncertainty=None, comment=''): + + KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, uncertainty=uncertainty, + comment=comment) + + self.alpha = alpha + self.A = A + self.n = n + self.E0 = E0 + self.electrons = electrons + self.V0 = V0 + + def __repr__(self): + """ + Return a string representation that can be used to reconstruct the + Arrhenius object. + """ + string = 'SurfaceChargeTransferBEP(A={0!r}, n={1!r}, E0={2!r}, V0={3!r}, alpha={4!r}, electrons={5!r}'.format( + self.A, self.n, self.E0, self.V0, self.alpha, self.electrons) + if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin) + if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) + if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin) + if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax) + if self.uncertainty: string += ', uncertainty={0!r}'.format(self.uncertainty) + if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) + string += ')' + return string + + def __reduce__(self): + """ + A helper function used when pickling a SurfaceChargeTransfer object. + """ + return (SurfaceChargeTransferBEP, (self.A, self.n, self.E0, self.V0, self.alpha, self.electrons, self.Tmin, self.Tmax, self.Pmin, self.Pmax, + self.uncertainty, self.comment)) + + property A: + """The preexponential factor.""" + def __get__(self): + return self._A + def __set__(self, value): + self._A = quantity.SurfaceRateCoefficient(value) + + property n: + """The temperature exponent.""" + def __get__(self): + return self._n + def __set__(self, value): + self._n = quantity.Dimensionless(value) + + property E0: + """The activation energy.""" + def __get__(self): + return self._E0 + def __set__(self, value): + self._E0 = quantity.Energy(value) + + property V0: + """The reference potential.""" + def __get__(self): + return self._V0 + def __set__(self, value): + self._V0 = quantity.Potential(value) + + property electrons: + """The number of electrons transferred.""" + def __get__(self): + return self._electrons + def __set__(self, value): + self._electrons = quantity.Dimensionless(value) + + property alpha: + """The charge transfer coefficient.""" + def __get__(self): + return self._alpha + def __set__(self, value): + self._alpha = quantity.Dimensionless(value) + + cpdef change_v0(self, double V0): + """ + Changes the reference potential to `V0` in volts, and adjusts the + activation energy `E0` accordingly. + """ + + self._E0.value_si = self.get_activation_energy_from_potential(V0,0.0) + self._V0.value_si = V0 + + cpdef double get_activation_energy(self, double dGrxn) except -1: + """ + Return the activation energy in J/mol corresponding to the given + free energy of reaction `dGrxn` in J/mol at the reference potential. + """ + cdef double Ea + Ea = self._alpha.value_si * dGrxn + self._E0.value_si + + if Ea < 0.0: + Ea = 0.0 + elif dGrxn > 0.0 and Ea < dGrxn: + Ea = dGrxn + + return Ea + + cpdef double get_activation_energy_from_potential(self, double V, double dGrxn) except -1: + """ + Return the activation energy in J/mol corresponding to the given + free energy of reaction `dGrxn` in J/mol. + """ + cdef double Ea + Ea = self.get_activation_energy(dGrxn) + Ea -= self._alpha.value_si * self._electrons.value_si * constants.F * (V-self._V0.value_si) + + return Ea + + cpdef double get_rate_coefficient_from_potential(self, double T, double V, double dGrxn) except -1: + """ + Return the rate coefficient in the appropriate combination of m^3, + mol, and s at temperature `T` in K, potential `V` in volts, and + free of reaction `dGrxn` in J/mol. + """ + cdef double A, n, Ea + Ea = self.get_activation_energy_from_potential(V,dGrxn) + A = self._A.value_si + n = self._n.value_si + return A * T ** n * exp(-Ea / (constants.R * T)) + + cpdef SurfaceChargeTransfer to_surface_charge_transfer(self, double dGrxn): + """ + Return an :class:`SurfaceChargeTransfer` instance of the kinetics model using the + given free energy of reaction `dGrxn` to determine the activation energy. + """ + return SurfaceChargeTransfer( + A=self.A, + n=self.n, + electrons=self.electrons, + Ea=(self.get_activation_energy(dGrxn) * 0.001, "kJ/mol"), + V0=self.V0, + T0=(1, "K"), + Tmin=self.Tmin, + Tmax=self.Tmax, + Pmin=self.Pmin, + Pmax=self.Pmax, + uncertainty=self.uncertainty, + comment=self.comment, + ) + + cpdef bint is_identical_to(self, KineticsModel other_kinetics) except -2: + """ + Returns ``True`` if kinetics matches that of another kinetics model. Must match temperature + and pressure range of kinetics model, as well as parameters: A, n, Ea, T0. (Shouldn't have pressure + range if it's Arrhenius.) Otherwise returns ``False``. + """ + if not isinstance(other_kinetics, SurfaceChargeTransferBEP): + return False + if not KineticsModel.is_identical_to(self, other_kinetics): + return False + if (not self.A.equals(other_kinetics.A) or not self.n.equals(other_kinetics.n) + or not self.E0.equals(other_kinetics.E0) or not self.alpha.equals(other_kinetics.alpha) + or not self.electrons.equals(other_kinetics.electrons) or not self.V0.equals(other_kinetics.V0)): + return False + + return True + + cpdef change_rate(self, double factor): + """ + Changes A factor by multiplying it by a ``factor``. + """ + self._A.value_si *= factor + + def set_cantera_kinetics(self, ct_reaction, species_list): + """ + Sets a cantera ElementaryReaction() object with the modified Arrhenius object + converted to an Arrhenius form. + """ + raise NotImplementedError('set_cantera_kinetics() is not implemented for ArrheniusEP class kinetics.') diff --git a/rmgpy/kinetics/surfaceTest.py b/rmgpy/kinetics/surfaceTest.py index f5a4a13ee0..a4ddf4052a 100644 --- a/rmgpy/kinetics/surfaceTest.py +++ b/rmgpy/kinetics/surfaceTest.py @@ -35,10 +35,11 @@ import numpy as np -from rmgpy.kinetics.surface import StickingCoefficient, SurfaceArrhenius +from rmgpy.kinetics.surface import StickingCoefficient, SurfaceArrhenius, SurfaceChargeTransfer from rmgpy.species import Species from rmgpy.molecule import Molecule import rmgpy.quantity as quantity +import rmgpy.constants as constants ################################################################################ @@ -457,3 +458,317 @@ def test_is_identical_to(self): Test that the SurfaceArrhenius.is_identical_to method works on itself """ self.assertTrue(self.surfarr.is_identical_to(self.surfarr)) + + def test_to_surface_charge_transfer(self): + """ + Test that the SurfaceArrhenius.to_surface_charge_transfer method works + """ + + surface_charge_transfer = self.surfarr.to_surface_charge_transfer(2,-2) + self.assertIsInstance(surface_charge_transfer,SurfaceChargeTransfer) + surface_charge_transfer0 = SurfaceChargeTransfer( + A = self.surfarr.A, + n = self.surfarr.n, + Ea = self.surfarr.Ea, + T0 = self.surfarr.T0, + Tmin = self.surfarr.Tmin, + Tmax = self.surfarr.Tmax, + electrons = -2, + V0 = (2,'V') + ) + self.assertTrue(surface_charge_transfer.is_identical_to(surface_charge_transfer0)) + + +################################################################################ + + +class TestSurfaceChargeTransfer((unittest.TestCase)): + """ + Contains unit tests of the :class:`SurfaceChargeTransfer` class. + """ + + def setUp(self): + """ + A function run before each unit test in this class. + """ + self.A = 1.44e18 + self.n = -0.087 + self.Ea = 63.4 + self.T0 = 1. + self.Tmin = 300. + self.Tmax = 3000. + self.V0 = 1 + self.electrons = -1 + self.comment = 'CH3x + Hx <=> CH4 + x + x' + self.surfchargerxn_reduction = SurfaceChargeTransfer( + A=(self.A, 'm^2/(mol*s)'), + n=self.n, + electrons=self.electrons, + V0=(self.V0, "V"), + Ea=(self.Ea, "kJ/mol"), + T0=(self.T0, "K"), + Tmin=(self.Tmin, "K"), + Tmax=(self.Tmax, "K"), + comment=self.comment, + ) + + self.surfchargerxn_oxidation = SurfaceChargeTransfer( + A=(self.A, 'm^2/(mol*s)'), + n=self.n, + electrons=1, + V0=(self.V0, "V"), + Ea=(self.Ea, "kJ/mol"), + T0=(self.T0, "K"), + Tmin=(self.Tmin, "K"), + Tmax=(self.Tmax, "K"), + comment=self.comment, + ) + + def test_A(self): + """ + Test that the SurfaceChargeTransfer A property was properly set. + """ + self.assertAlmostEqual(self.surfchargerxn_reduction.A.value_si, self.A, delta=1e0) + + def test_n(self): + """ + Test that the SurfaceChargeTransfer n property was properly set. + """ + self.assertAlmostEqual(self.surfchargerxn_reduction.n.value_si, self.n, 6) + + def test_ne(self): + """ + Test that the SurfaceChargeTransfer electrons property was properly set. + """ + self.assertEqual(self.surfchargerxn_reduction.electrons.value_si, -1.0) + + def test_alpha(self): + """ + Test that the SurfaceChargeTransfer alpha property was properly set. + """ + self.assertEqual(self.surfchargerxn_reduction.alpha.value_si, 0.5) + + def test_Ea(self): + """ + Test that the SurfaceChargeTransfer Ea property was properly set. + """ + self.assertAlmostEqual(self.surfchargerxn_reduction.Ea.value_si * 0.001, self.Ea, 6) + + def test_T0(self): + """ + Test that the SurfaceChargeTransfer T0 property was properly set. + """ + self.assertAlmostEqual(self.surfchargerxn_reduction.T0.value_si, self.T0, 6) + + def test_Tmin(self): + """ + Test that the SurfaceChargeTransfer Tmin property was properly set. + """ + self.assertAlmostEqual(self.surfchargerxn_reduction.Tmin.value_si, self.Tmin, 6) + + def test_Tmax(self): + """ + Test that the SurfaceChargeTransfer Tmax property was properly set. + """ + self.assertAlmostEqual(self.surfchargerxn_reduction.Tmax.value_si, self.Tmax, 6) + + def test_V0(self): + """ + Test that the SurfaceChargeTransfer V0 property was properly set. + """ + self.assertAlmostEqual(self.surfchargerxn_reduction.V0.value_si, self.V0, 1) + + def test_Tmax(self): + """ + Test that the SurfaceChargeTransfer Tmax property was properly set. + """ + self.assertAlmostEqual(self.surfchargerxn_reduction.Tmax.value_si, self.Tmax, 6) + + def test_comment(self): + """ + Test that the SurfaceChargeTransfer comment property was properly set. + """ + self.assertEqual(self.surfchargerxn_reduction.comment, self.comment) + + def test_is_temperature_valid(self): + """ + Test the SurfaceChargeTransfer.is_temperature_valid() method. + """ + T_data = np.array([200, 400, 600, 800, 1000, 1200, 1400, 1600, 1800, 4000]) + valid_data = np.array([False, True, True, True, True, True, True, True, True, False], np.bool) + for T, valid in zip(T_data, valid_data): + valid0 = self.surfchargerxn_reduction.is_temperature_valid(T) + self.assertEqual(valid0, valid) + + def test_pickle(self): + """ + Test that an SurfaceChargeTransfer object can be pickled and unpickled with no loss + of information. + """ + import pickle + surfchargerxn_reduction = pickle.loads(pickle.dumps(self.surfchargerxn_reduction, -1)) + self.assertAlmostEqual(self.surfchargerxn_reduction.A.value, surfchargerxn_reduction.A.value, delta=1e0) + self.assertEqual(self.surfchargerxn_reduction.A.units, surfchargerxn_reduction.A.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.n.value, surfchargerxn_reduction.n.value, 4) + self.assertAlmostEqual(self.surfchargerxn_reduction.electrons.value, surfchargerxn_reduction.electrons.value, 4) + self.assertAlmostEqual(self.surfchargerxn_reduction.Ea.value, surfchargerxn_reduction.Ea.value, 4) + self.assertEqual(self.surfchargerxn_reduction.Ea.units, surfchargerxn_reduction.Ea.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.T0.value, surfchargerxn_reduction.T0.value, 4) + self.assertEqual(self.surfchargerxn_reduction.T0.units, surfchargerxn_reduction.T0.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.V0.value, surfchargerxn_reduction.V0.value, 4) + self.assertEqual(self.surfchargerxn_reduction.V0.units, surfchargerxn_reduction.V0.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.Tmin.value, surfchargerxn_reduction.Tmin.value, 4) + self.assertEqual(self.surfchargerxn_reduction.Tmin.units, surfchargerxn_reduction.Tmin.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.Tmax.value, surfchargerxn_reduction.Tmax.value, 4) + self.assertEqual(self.surfchargerxn_reduction.Tmax.units, surfchargerxn_reduction.Tmax.units) + self.assertEqual(self.surfchargerxn_reduction.comment, surfchargerxn_reduction.comment) + self.assertEqual(dir(self.surfchargerxn_reduction), dir(surfchargerxn_reduction)) + + def test_repr(self): + """ + Test that an SurfaceChargeTransfer object can be reconstructed from its repr() + output with no loss of information. + """ + namespace = {} + exec('surfchargerxn_reduction = {0!r}'.format(self.surfchargerxn_reduction), globals(), namespace) + self.assertIn('surfchargerxn_reduction', namespace) + surfchargerxn_reduction = namespace['surfchargerxn_reduction'] + self.assertAlmostEqual(self.surfchargerxn_reduction.A.value, surfchargerxn_reduction.A.value, delta=1e0) + self.assertEqual(self.surfchargerxn_reduction.A.units, surfchargerxn_reduction.A.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.n.value, surfchargerxn_reduction.n.value, 4) + self.assertAlmostEqual(self.surfchargerxn_reduction.Ea.value, surfchargerxn_reduction.Ea.value, 4) + self.assertEqual(self.surfchargerxn_reduction.Ea.units, surfchargerxn_reduction.Ea.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.T0.value, surfchargerxn_reduction.T0.value, 4) + self.assertEqual(self.surfchargerxn_reduction.T0.units, surfchargerxn_reduction.T0.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.Tmin.value, surfchargerxn_reduction.Tmin.value, 4) + self.assertEqual(self.surfchargerxn_reduction.Tmin.units, surfchargerxn_reduction.Tmin.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.Tmax.value, surfchargerxn_reduction.Tmax.value, 4) + self.assertEqual(self.surfchargerxn_reduction.Tmax.units, surfchargerxn_reduction.Tmax.units) + self.assertEqual(self.surfchargerxn_reduction.comment, surfchargerxn_reduction.comment) + self.assertEqual(dir(self.surfchargerxn_reduction), dir(surfchargerxn_reduction)) + + def test_copy(self): + """ + Test that an SurfaceChargeTransfer object can be copied with deepcopy + with no loss of information. + """ + import copy + surfchargerxn_reduction = copy.deepcopy(self.surfchargerxn_reduction) + self.assertAlmostEqual(self.surfchargerxn_reduction.A.value, surfchargerxn_reduction.A.value, delta=1e0) + self.assertEqual(self.surfchargerxn_reduction.A.units, surfchargerxn_reduction.A.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.n.value, surfchargerxn_reduction.n.value, 4) + self.assertAlmostEqual(self.surfchargerxn_reduction.Ea.value, surfchargerxn_reduction.Ea.value, 4) + self.assertEqual(self.surfchargerxn_reduction.Ea.units, surfchargerxn_reduction.Ea.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.T0.value, surfchargerxn_reduction.T0.value, 4) + self.assertEqual(self.surfchargerxn_reduction.T0.units, surfchargerxn_reduction.T0.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.Tmin.value, surfchargerxn_reduction.Tmin.value, 4) + self.assertEqual(self.surfchargerxn_reduction.Tmin.units, surfchargerxn_reduction.Tmin.units) + self.assertAlmostEqual(self.surfchargerxn_reduction.Tmax.value, surfchargerxn_reduction.Tmax.value, 4) + self.assertEqual(self.surfchargerxn_reduction.Tmax.units, surfchargerxn_reduction.Tmax.units) + self.assertEqual(self.surfchargerxn_reduction.comment, surfchargerxn_reduction.comment) + self.assertEqual(dir(self.surfchargerxn_reduction), dir(surfchargerxn_reduction)) + + def test_is_identical_to(self): + """ + Test that the SurfaceChargeTransfer.is_identical_to method works on itself + """ + self.assertTrue(self.surfchargerxn_reduction.is_identical_to(self.surfchargerxn_reduction)) + + def test_to_surface_arrhenius(self): + """ + Test that the SurfaceChargeTransfer.to_surface_arrhenius method works + """ + surface_arr = self.surfchargerxn_reduction.to_surface_arrhenius() + self.assertIsInstance(surface_arr,SurfaceArrhenius) + surface_arrhenius0 = SurfaceArrhenius( + A = self.surfchargerxn_reduction.A, + n = self.surfchargerxn_reduction.n, + Ea = self.surfchargerxn_reduction.Ea, + T0 = self.surfchargerxn_reduction.T0, + Tmin = self.surfchargerxn_reduction.Tmin, + Tmax = self.surfchargerxn_reduction.Tmax, + ) + + self.assertTrue(surface_arr.is_identical_to(surface_arrhenius0)) + + def test_get_activation_energy_from_potential(self): + """ + Test that the SurfaceChargeTransfer.get_activation_energy_from_potential method works + """ + + electrons_ox = self.surfchargerxn_oxidation.electrons.value_si + V0_ox = self.surfchargerxn_oxidation.V0.value_si + Ea0_ox = self.surfchargerxn_oxidation.Ea.value_si + alpha_ox = self.surfchargerxn_oxidation.alpha.value_si + + electrons_red = self.surfchargerxn_reduction.electrons.value_si + V0_red = self.surfchargerxn_reduction.V0.value_si + Ea0_red = self.surfchargerxn_reduction.Ea.value_si + alpha_red = self.surfchargerxn_reduction.alpha.value_si + + Potentials = (V0_ox + 1, V0_ox, V0_ox - 1) + + for V in Potentials: + Ea = self.surfchargerxn_oxidation.get_activation_energy_from_potential(V, False) + self.assertAlmostEqual(Ea0_ox - (alpha_ox * electrons_ox * constants.F * (V-V0_ox)), Ea, 6) + Ea = self.surfchargerxn_oxidation.get_activation_energy_from_potential(V, True) + self.assertTrue(Ea>=0) + Ea = self.surfchargerxn_reduction.get_activation_energy_from_potential(V, False) + self.assertAlmostEqual(Ea0_red - (alpha_red * electrons_red * constants.F * (V-V0_red)), Ea, 6) + + def test_get_rate_coefficient(self): + """ + Test that the SurfaceChargeTransfer.to_surface_arrhenius method works + """ + + A_ox = self.surfchargerxn_oxidation.A.value_si + A_red = self.surfchargerxn_reduction.A.value_si + electrons_ox = self.surfchargerxn_oxidation.electrons.value_si + electrons_red = self.surfchargerxn_reduction.electrons.value_si + n_ox = self.surfchargerxn_oxidation.n.value_si + n_red = self.surfchargerxn_reduction.n.value_si + Ea_ox = self.surfchargerxn_oxidation.Ea.value_si + Ea_red = self.surfchargerxn_reduction.Ea.value_si + V0_ox = self.surfchargerxn_oxidation.V0.value_si + V0_red = self.surfchargerxn_reduction.V0.value_si + T0_ox = self.surfchargerxn_oxidation.T0.value_si + T0_red = self.surfchargerxn_reduction.T0.value_si + alpha_ox = self.surfchargerxn_oxidation.alpha.value + alpha_red = self.surfchargerxn_reduction.alpha.value + + Potentials = (V0_ox + 1, V0_ox, V0_ox - 1) + for V in Potentials: + for T in np.linspace(300,3000,10): + Ea = Ea_ox - (alpha_ox * electrons_ox * constants.F * (V-V0_ox)) + k_oxidation = A_ox * (T / T0_ox) ** n_ox * np.exp(-Ea / (constants.R * T)) + Ea = Ea_red - (alpha_red * electrons_red * constants.F * (V-V0_red)) + k_reduction = A_red * (T / T0_red) ** n_red * np.exp(-Ea / (constants.R * T)) + self.assertAlmostEqual(k_oxidation,self.surfchargerxn_oxidation.get_rate_coefficient(T,V)) + self.assertAlmostEqual(k_reduction,self.surfchargerxn_reduction.get_rate_coefficient(T,V)) + + def test_change_v0(self): + + V0 = self.surfchargerxn_oxidation.V0.value_si + electrons = self.surfchargerxn_oxidation.electrons.value + alpha = self.surfchargerxn_oxidation.alpha.value + for V in (V0 + 1, V0, V0 - 1, V0): + delta = V - self.surfchargerxn_oxidation.V0.value_si + V_i = self.surfchargerxn_oxidation.V0.value_si + Ea_i = self.surfchargerxn_oxidation.Ea.value_si + self.surfchargerxn_oxidation.change_v0(V) + self.assertEqual(self.surfchargerxn_oxidation.V0.value_si, V_i + delta) + self.assertAlmostEqual(self.surfchargerxn_oxidation.Ea.value_si, Ea_i - (alpha *electrons * constants.F * delta), 6) + + V0 = self.surfchargerxn_reduction.V0.value_si + electrons = self.surfchargerxn_reduction.electrons.value + alpha = self.surfchargerxn_reduction.alpha.value + for V in (V0 + 1, V0, V0 - 1, V0): + delta = V - self.surfchargerxn_reduction.V0.value_si + V_i = self.surfchargerxn_reduction.V0.value_si + Ea_i = self.surfchargerxn_reduction.Ea.value_si + self.surfchargerxn_reduction.change_v0(V) + self.assertEqual(self.surfchargerxn_reduction.V0.value_si, V_i + delta) + self.assertAlmostEqual(self.surfchargerxn_reduction.Ea.value_si, Ea_i - (alpha *electrons * constants.F * delta), 6) + +if __name__ == '__main__': + unittest.main(testRunner=unittest.TextTestRunner(verbosity=2)) \ No newline at end of file diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index 327e1f2c12..c8d16116e3 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -89,7 +89,7 @@ def check_partial_charge(atom): the theoretical one: """ - if atom.symbol == 'X': + if atom.symbol in ('X','e','H+'): return # because we can't check it. valence = PeriodicSystem.valence_electrons[atom.symbol] diff --git a/rmgpy/molecule/atomtype.pxd b/rmgpy/molecule/atomtype.pxd index 3aa2081a63..105ae31891 100644 --- a/rmgpy/molecule/atomtype.pxd +++ b/rmgpy/molecule/atomtype.pxd @@ -39,6 +39,8 @@ cdef class AtomType: cdef public list decrement_radical cdef public list increment_lone_pair cdef public list decrement_lone_pair + cdef public list increment_charge + cdef public list decrement_charge cdef public list single cdef public list all_double diff --git a/rmgpy/molecule/atomtype.py b/rmgpy/molecule/atomtype.py index 06620fcc17..badc15ea11 100644 --- a/rmgpy/molecule/atomtype.py +++ b/rmgpy/molecule/atomtype.py @@ -65,6 +65,8 @@ class AtomType: `break_bond` ``list`` The atom type(s) that result when an existing single bond to this atom type is broken `increment_radical` ``list`` The atom type(s) that result when the number of radical electrons is incremented `decrement_radical` ``list`` The atom type(s) that result when the number of radical electrons is decremented + `increment_charge` ``list`` The atom type(s) that result when the number of radical electrons is decremented and charge is incremented + `decrement_charge` ``list`` The atom type(s) that result when the number of radical electrons is incremented and charge is decremented `increment_lone_pair` ``list`` The atom type(s) that result when the number of lone electron pairs is incremented `decrement_lone_pair` ``list`` The atom type(s) that result when the number of lone electron pairs is decremented @@ -106,6 +108,8 @@ def __init__(self, label='', self.break_bond = [] self.increment_radical = [] self.decrement_radical = [] + self.increment_charge = [] + self.decrement_charge = [] self.increment_lone_pair = [] self.decrement_lone_pair = [] self.single = single or [] @@ -136,6 +140,8 @@ def __reduce__(self): 'break_bond': self.break_bond, 'increment_radical': self.increment_radical, 'decrement_radical': self.decrement_radical, + 'increment_charge': self.increment_charge, + 'decrement_charge': self.decrement_charge, 'increment_lone_pair': self.increment_lone_pair, 'decrement_lone_pair': self.decrement_lone_pair, 'single': self.single, @@ -164,6 +170,8 @@ def __setstate__(self, d): self.break_bond = d['break_bond'] self.increment_radical = d['increment_radical'] self.decrement_radical = d['decrement_radical'] + self.increment_charge = d['increment_charge'] + self.decrement_charge = d['decrement_charge'] self.increment_lone_pair = d['increment_lone_pair'] self.decrement_lone_pair = d['decrement_lone_pair'] self.single = d['single'] @@ -178,7 +186,7 @@ def __setstate__(self, d): self.charge = d['charge'] def set_actions(self, increment_bond, decrement_bond, form_bond, break_bond, increment_radical, decrement_radical, - increment_lone_pair, decrement_lone_pair): + increment_lone_pair, decrement_lone_pair, increment_charge, decrement_charge): self.increment_bond = increment_bond self.decrement_bond = decrement_bond self.form_bond = form_bond @@ -187,6 +195,8 @@ def set_actions(self, increment_bond, decrement_bond, form_bond, break_bond, inc self.decrement_radical = decrement_radical self.increment_lone_pair = increment_lone_pair self.decrement_lone_pair = decrement_lone_pair + self.increment_charge = increment_charge + self.decrement_charge = decrement_charge def equivalent(self, other): """ @@ -202,7 +212,7 @@ def is_specific_case_of(self, other): atom type `atomType2` or ``False`` otherwise. """ return self is other or self in other.specific - + def get_features(self): """ Returns a list of the features that are checked to determine atomtype @@ -240,6 +250,9 @@ def get_features(self): ATOMTYPES = {} +# Electron +ATOMTYPES['e'] = AtomType(label='e', generic=[], specific=[], lone_pairs=[0], charge=[-1]) + # Surface sites: ATOMTYPES['X'] = AtomType(label='X', generic=[], specific=['Xv', 'Xo']) @@ -254,7 +267,7 @@ def get_features(self): # Non-surface atomTypes, R being the most generic: ATOMTYPES['R'] = AtomType(label='R', generic=[], specific=[ - 'H', + 'H','H0','H+', 'R!H', 'R!H!Val7', 'Val4','Val5','Val6','Val7', @@ -313,7 +326,12 @@ def get_features(self): 'I','I1s', 'F','F1s']) -ATOMTYPES['H'] = AtomType('H', generic=['R'], specific=[]) + +ATOMTYPES['H'] = AtomType('H', generic=['R'], specific=['H0','H+'], charge=[0,+1]) +ATOMTYPES['H0'] = AtomType('H0', generic=['R','H'], specific=[], single=[0,1], all_double=[0], r_double=[0], o_double=[0], s_double=[0], triple=[0], + quadruple=[0], benzene=[0], lone_pairs=[0], charge=[0]) +ATOMTYPES['H+'] = AtomType('H+', generic=['R','H'], specific=[], single=[0], all_double=[0], r_double=[0], o_double=[0], s_double=[0], triple=[0], + quadruple=[0], benzene=[0], lone_pairs=[0], charge=[+1]) ATOMTYPES['He'] = AtomType('He', generic=['R', 'R!H', 'R!H!Val7'], specific=[]) ATOMTYPES['Ne'] = AtomType('Ne', generic=['R', 'R!H', 'R!H!Val7'], specific=[]) @@ -631,158 +649,162 @@ def get_features(self): single=[0,1], all_double=[0], r_double=[], o_double=[], s_double=[], triple=[0], quadruple=[0], benzene=[0], lone_pairs=[3], charge=[0]) # examples for F1s: HF, [F], FO, CH3F, F2 -ATOMTYPES['X'].set_actions(increment_bond=['X'], decrement_bond=['X'], form_bond=['X'], break_bond=['X'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Xv'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Xo'], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Xo'].set_actions(increment_bond=['Xo'], decrement_bond=['Xo'], form_bond=[], break_bond=['Xv'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['R'].set_actions(increment_bond=['R'], decrement_bond=['R'], form_bond=['R'], break_bond=['R'], increment_radical=['R'], decrement_radical=['R'], increment_lone_pair=['R'], decrement_lone_pair=['R']) -ATOMTYPES['R!H'].set_actions(increment_bond=['R!H'], decrement_bond=['R!H'], form_bond=['R!H'], break_bond=['R!H'], increment_radical=['R!H'], decrement_radical=['R!H'], increment_lone_pair=['R!H'], decrement_lone_pair=['R!H']) -ATOMTYPES['R!H!Val7'].set_actions(increment_bond=['R!H!Val7'], decrement_bond=['R!H!Val7'], form_bond=['R!H!Val7'], break_bond=['R!H!Val7'], increment_radical=['R!H!Val7'], decrement_radical=['R!H!Val7'], increment_lone_pair=['R!H!Val7'], decrement_lone_pair=['R!H!Val7']) -ATOMTYPES['Val4'].set_actions(increment_bond=['Val4'], decrement_bond=['Val4'], form_bond=['Val4'], break_bond=['Val4'], increment_radical=['Val4'], decrement_radical=['Val4'], increment_lone_pair=['Val4'], decrement_lone_pair=['Val4']) -ATOMTYPES['Val5'].set_actions(increment_bond=['Val5'], decrement_bond=['Val5'], form_bond=['Val5'], break_bond=['Val5'], increment_radical=['Val5'], decrement_radical=['Val5'], increment_lone_pair=['Val5'], decrement_lone_pair=['Val5']) -ATOMTYPES['Val6'].set_actions(increment_bond=['Val6'], decrement_bond=['Val6'], form_bond=['Val6'], break_bond=['Val6'], increment_radical=['Val6'], decrement_radical=['Val6'], increment_lone_pair=['Val6'], decrement_lone_pair=['Val6']) -ATOMTYPES['Val7'].set_actions(increment_bond=['Val7'], decrement_bond=['Val7'], form_bond=['Val7'], break_bond=['Val7'], increment_radical=['Val7'], decrement_radical=['Val7'], increment_lone_pair=['Val7'], decrement_lone_pair=['Val7']) - -ATOMTYPES['H'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['H'], break_bond=['H'], increment_radical=['H'], decrement_radical=['H'], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['He'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=['He'], decrement_radical=['He'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Ne'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=['Ne'], decrement_radical=['Ne'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Ar'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['C'].set_actions(increment_bond=['C'], decrement_bond=['C'], form_bond=['C'], break_bond=['C'], increment_radical=['C'], decrement_radical=['C'], increment_lone_pair=['C'], decrement_lone_pair=['C']) -ATOMTYPES['Ca'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['C2s']) -ATOMTYPES['Cs'].set_actions(increment_bond=['Cd', 'CO', 'CS'], decrement_bond=[], form_bond=['Cs', 'Csc'], break_bond=['Cs'], increment_radical=['Cs'], decrement_radical=['Cs'], increment_lone_pair=['C2s'], decrement_lone_pair=['C2s']) -ATOMTYPES['Csc'].set_actions(increment_bond=['Cdc'], decrement_bond=[], form_bond=['Csc'], break_bond=['Csc', 'Cs'], increment_radical=['Csc'], decrement_radical=['Csc'], increment_lone_pair=['C2sc'], decrement_lone_pair=['C2sc']) -ATOMTYPES['Cd'].set_actions(increment_bond=['Cdd', 'Ct', 'C2tc'], decrement_bond=['Cs'], form_bond=['Cd', 'Cdc'], break_bond=['Cd'], increment_radical=['Cd'], decrement_radical=['Cd'], increment_lone_pair=['C2d'], decrement_lone_pair=[]) -ATOMTYPES['Cdc'].set_actions(increment_bond=[], decrement_bond=['Csc'], form_bond=['Cdc'], break_bond=['Cdc', 'Cd', 'CO', 'CS'], increment_radical=['Cdc'], decrement_radical=['Cdc'], increment_lone_pair=['C2dc'], decrement_lone_pair=[]) -ATOMTYPES['CO'].set_actions(increment_bond=['Cdd', 'C2tc'], decrement_bond=['Cs'], form_bond=['CO', 'Cdc'], break_bond=['CO'], increment_radical=['CO'], decrement_radical=['CO'], increment_lone_pair=['C2d'], decrement_lone_pair=[]) -ATOMTYPES['CS'].set_actions(increment_bond=['Cdd', 'C2tc'], decrement_bond=['Cs'], form_bond=['CS', 'Cdc'], break_bond=['CS'], increment_radical=['CS'], decrement_radical=['CS'], increment_lone_pair=['C2d'], decrement_lone_pair=[]) -ATOMTYPES['Cdd'].set_actions(increment_bond=[], decrement_bond=['Cd', 'CO', 'CS'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Ct'].set_actions(increment_bond=['Cq'], decrement_bond=['Cd', 'CO', 'CS'], form_bond=['Ct'], break_bond=['Ct'], increment_radical=['Ct'], decrement_radical=['Ct'], increment_lone_pair=['C2tc'], decrement_lone_pair=[]) -ATOMTYPES['Cb'].set_actions(increment_bond=['Cbf'], decrement_bond=[], form_bond=['Cb'], break_bond=['Cb'], increment_radical=['Cb'], decrement_radical=['Cb'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Cbf'].set_actions(increment_bond=[], decrement_bond=['Cb'], form_bond=[], break_bond=['Cb'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['C2s'].set_actions(increment_bond=['C2d'], decrement_bond=[], form_bond=['C2s'], break_bond=['C2s'], increment_radical=['C2s'], decrement_radical=['C2s'], increment_lone_pair=['Ca'], decrement_lone_pair=['Cs']) -ATOMTYPES['C2sc'].set_actions(increment_bond=['C2dc'], decrement_bond=[], form_bond=['C2sc'], break_bond=['C2sc'], increment_radical=['C2sc'], decrement_radical=['C2sc'], increment_lone_pair=[], decrement_lone_pair=['Cs']) -ATOMTYPES['C2d'].set_actions(increment_bond=['C2tc'], decrement_bond=['C2s'], form_bond=['C2dc'], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['Cd', 'CO', 'CS']) -ATOMTYPES['C2dc'].set_actions(increment_bond=[], decrement_bond=['C2sc'], form_bond=[], break_bond=['C2d'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['Cdc']) -ATOMTYPES['C2tc'].set_actions(increment_bond=[], decrement_bond=['C2d'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['Ct']) -ATOMTYPES['Cq'].set_actions(increment_bond=[], decrement_bond=['Ct'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['N'].set_actions(increment_bond=['N'], decrement_bond=['N'], form_bond=['N'], break_bond=['N'], increment_radical=['N'], decrement_radical=['N'], increment_lone_pair=['N'], decrement_lone_pair=['N']) -ATOMTYPES['N0sc'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['N0sc'], break_bond=['N0sc'], increment_radical=['N0sc'], decrement_radical=['N0sc'], increment_lone_pair=[], decrement_lone_pair=['N1s', 'N1sc']) -ATOMTYPES['N1s'].set_actions(increment_bond=['N1dc'], decrement_bond=[], form_bond=['N1s'], break_bond=['N1s'], increment_radical=['N1s'], decrement_radical=['N1s'], increment_lone_pair=['N0sc'], decrement_lone_pair=['N3s', 'N3sc']) -ATOMTYPES['N1sc'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['N1sc'], break_bond=['N1sc'], increment_radical=['N1sc'], decrement_radical=['N1sc'], increment_lone_pair=[], decrement_lone_pair=['N3s', 'N3sc']) -ATOMTYPES['N1dc'].set_actions(increment_bond=['N1dc'], decrement_bond=['N1s', 'N1dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['N3d']) -ATOMTYPES['N3s'].set_actions(increment_bond=['N3d'], decrement_bond=[], form_bond=['N3s'], break_bond=['N3s'], increment_radical=['N3s'], decrement_radical=['N3s'], increment_lone_pair=['N1s', 'N1sc'], decrement_lone_pair=['N5sc']) -ATOMTYPES['N3sc'].set_actions(increment_bond=['N3d'], decrement_bond=[], form_bond=['N3sc'], break_bond=['N3sc'], increment_radical=['N3sc'], decrement_radical=['N3sc'], increment_lone_pair=['N1s', 'N1sc'], decrement_lone_pair=['N5sc']) -ATOMTYPES['N3d'].set_actions(increment_bond=['N3t'], decrement_bond=['N3s', 'N3sc'], form_bond=['N3d'], break_bond=['N3d'], increment_radical=['N3d'], decrement_radical=['N3d'], increment_lone_pair=['N1dc'], decrement_lone_pair=['N5dc']) -ATOMTYPES['N3t'].set_actions(increment_bond=[], decrement_bond=['N3d'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['N3b'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['N5sc'].set_actions(increment_bond=['N5dc'], decrement_bond=[], form_bond=['N5sc'], break_bond=['N5sc'], increment_radical=['N5sc'], decrement_radical=['N5sc'], increment_lone_pair=['N3s', 'N3sc'], decrement_lone_pair=[]) -ATOMTYPES['N5dc'].set_actions(increment_bond=['N5ddc', 'N5tc'], decrement_bond=['N5sc'], form_bond=['N5dc'], break_bond=['N5dc'], increment_radical=['N5dc'], decrement_radical=['N5dc'], increment_lone_pair=['N3d'], decrement_lone_pair=[]) -ATOMTYPES['N5ddc'].set_actions(increment_bond=['N5dddc'], decrement_bond=['N5dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['N5dddc'].set_actions(increment_bond=[], decrement_bond=['N5ddc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['N5tc'].set_actions(increment_bond=[], decrement_bond=['N5dc'], form_bond=['N5tc'], break_bond=['N5tc'], increment_radical=['N5tc'], decrement_radical=['N5tc'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['N5b'].set_actions(increment_bond=['N5bd'], decrement_bond=[], form_bond=['N5b'], break_bond=['N5b'], increment_radical=['N5b'], decrement_radical=['N5b'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['N5bd'].set_actions(increment_bond=[], decrement_bond=['N5b'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['O'].set_actions(increment_bond=['O'], decrement_bond=['O'], form_bond=['O'], break_bond=['O'], increment_radical=['O'], decrement_radical=['O'], increment_lone_pair=['O'], decrement_lone_pair=['O']) -ATOMTYPES['Oa'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['O0sc'], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['O2s', 'O2sc']) -ATOMTYPES['O0sc'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['O0sc'], break_bond=['Oa', 'O0sc'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['O2s', 'O2sc']) -ATOMTYPES['O2s'].set_actions(increment_bond=['O2d'], decrement_bond=[], form_bond=['O2s', 'O2sc'], break_bond=['O2s'], increment_radical=['O2s'], decrement_radical=['O2s'], increment_lone_pair=['Oa', 'O0sc'], decrement_lone_pair=['O4sc']) -ATOMTYPES['O2sc'].set_actions(increment_bond=['O2d'], decrement_bond=[], form_bond=[], break_bond=['O2s'], increment_radical=['O2sc'], decrement_radical=['O2sc'], increment_lone_pair=[], decrement_lone_pair=['O4sc']) -ATOMTYPES['O2d'].set_actions(increment_bond=[], decrement_bond=['O2s', 'O2sc'], form_bond=[], break_bond=[], increment_radical=['O2d'], decrement_radical=['O2d'], increment_lone_pair=[], decrement_lone_pair=['O4dc', 'O4tc']) -ATOMTYPES['O4sc'].set_actions(increment_bond=['O4dc'], decrement_bond=[], form_bond=['O4sc'], break_bond=['O4sc'], increment_radical=['O4sc'], decrement_radical=['O4sc'], increment_lone_pair=['O2s', 'O2sc'], decrement_lone_pair=[]) -ATOMTYPES['O4dc'].set_actions(increment_bond=['O4tc'], decrement_bond=['O4sc'], form_bond=['O4dc'], break_bond=['O4dc'], increment_radical=['O4dc'], decrement_radical=['O4dc'], increment_lone_pair=['O2d'], decrement_lone_pair=[]) -ATOMTYPES['O4tc'].set_actions(increment_bond=[], decrement_bond=['O4dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['O4b'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['Ne'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=['Ne'], decrement_radical=['Ne'], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['Si'].set_actions(increment_bond=['Si'], decrement_bond=['Si'], form_bond=['Si'], break_bond=['Si'], increment_radical=['Si'], decrement_radical=['Si'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Sis'].set_actions(increment_bond=['Sid', 'SiO'], decrement_bond=[], form_bond=['Sis'], break_bond=['Sis'], increment_radical=['Sis'], decrement_radical=['Sis'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Sid'].set_actions(increment_bond=['Sidd', 'Sit'], decrement_bond=['Sis'], form_bond=['Sid'], break_bond=['Sid'], increment_radical=['Sid'], decrement_radical=['Sid'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Sidd'].set_actions(increment_bond=[], decrement_bond=['Sid', 'SiO'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Sit'].set_actions(increment_bond=['Siq'], decrement_bond=['Sid'], form_bond=['Sit'], break_bond=['Sit'], increment_radical=['Sit'], decrement_radical=['Sit'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['SiO'].set_actions(increment_bond=['Sidd'], decrement_bond=['Sis'], form_bond=['SiO'], break_bond=['SiO'], increment_radical=['SiO'], decrement_radical=['SiO'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Sib'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Sib'], break_bond=['Sib'], increment_radical=['Sib'], decrement_radical=['Sib'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Sibf'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Siq'].set_actions(increment_bond=[], decrement_bond=['Sit'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['P'].set_actions(increment_bond=['P'], decrement_bond=['P'], form_bond=['P'], break_bond=['P'], increment_radical=['P'], decrement_radical=['P'], increment_lone_pair=['P'], decrement_lone_pair=['P']) -ATOMTYPES['P0sc'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['P0sc'], break_bond=['P0sc'], increment_radical=['P0sc'], decrement_radical=['P0sc'], increment_lone_pair=[], decrement_lone_pair=['P1s', 'P1sc']) -ATOMTYPES['P1s'].set_actions(increment_bond=['P1dc'], decrement_bond=[], form_bond=['P1s'], break_bond=['P1s'], increment_radical=['P1s'], decrement_radical=['P1s'], increment_lone_pair=['P0sc'], decrement_lone_pair=['P3s']) -ATOMTYPES['P1sc'].set_actions(increment_bond=['P1dc'], decrement_bond=[], form_bond=['P1sc'], break_bond=['P1sc'], increment_radical=['P1sc'], decrement_radical=['P1sc'], increment_lone_pair=['P0sc'], decrement_lone_pair=['P3s']) -ATOMTYPES['P1dc'].set_actions(increment_bond=[], decrement_bond=['P1s'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['P3d']) -ATOMTYPES['P3s'].set_actions(increment_bond=['P3d'], decrement_bond=[], form_bond=['P3s'], break_bond=['P3s'], increment_radical=['P3s'], decrement_radical=['P3s'], increment_lone_pair=['P1s', 'P1sc'], decrement_lone_pair=['P5s', 'P5sc']) -ATOMTYPES['P3d'].set_actions(increment_bond=['P3t'], decrement_bond=['P3s'], form_bond=['P3d'], break_bond=['P3d'], increment_radical=['P3d'], decrement_radical=['P3d'], increment_lone_pair=['P1dc'], decrement_lone_pair=['P5d', 'P5dc']) -ATOMTYPES['P3t'].set_actions(increment_bond=[], decrement_bond=['P3d'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['P5t', 'P5tc']) -ATOMTYPES['P3b'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['P5s'].set_actions(increment_bond=['P5d', 'P5dc'], decrement_bond=[], form_bond=['P5s'], break_bond=['P5s'], increment_radical=['P5s'], decrement_radical=['P5s'], increment_lone_pair=['P3s'], decrement_lone_pair=[]) -ATOMTYPES['P5sc'].set_actions(increment_bond=['P5dc'], decrement_bond=[], form_bond=['P5sc'], break_bond=['P5sc'], increment_radical=['P5sc'], decrement_radical=['P5sc'], increment_lone_pair=['P3s'], decrement_lone_pair=[]) -ATOMTYPES['P5d'].set_actions(increment_bond=['P5dd', 'P5ddc', 'P5t', 'P5tc'], decrement_bond=['P5s'], form_bond=['P5d'], break_bond=['P5d'], increment_radical=['P5d'], decrement_radical=['P5d'], increment_lone_pair=['P3d'], decrement_lone_pair=[]) -ATOMTYPES['P5dd'].set_actions(increment_bond=['P5td'], decrement_bond=['P5d', 'P5dc'], form_bond=['P5dd'], break_bond=['P5dd'], increment_radical=['P5dd'], decrement_radical=['P5dd'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['P5dc'].set_actions(increment_bond=['P5dd', 'P5ddc', 'P5tc'], decrement_bond=['P5sc'], form_bond=['P5dc'], break_bond=['P5dc'], increment_radical=['P5dc'], decrement_radical=['P5dc'], increment_lone_pair=['P3d'], decrement_lone_pair=[]) -ATOMTYPES['P5ddc'].set_actions(increment_bond=[], decrement_bond=['P5dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['P5t'].set_actions(increment_bond=['P5td'], decrement_bond=['P5d'], form_bond=['P5t'], break_bond=['P5t'], increment_radical=['P5t'], decrement_radical=['P5t'], increment_lone_pair=['P3t'], decrement_lone_pair=[]) -ATOMTYPES['P5td'].set_actions(increment_bond=[], decrement_bond=['P5t', 'P5dd'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['P5tc'].set_actions(increment_bond=[], decrement_bond=['P5dc'], form_bond=['P5tc'], break_bond=['P5tc'], increment_radical=['P5tc'], decrement_radical=['P5tc'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['P5b'].set_actions(increment_bond=['P5bd'], decrement_bond=[], form_bond=['P5b'], break_bond=['P5b'], increment_radical=['P5b'], decrement_radical=['P5b'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['P5bd'].set_actions(increment_bond=[], decrement_bond=['P5b'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['S'].set_actions(increment_bond=['S'], decrement_bond=['S'], form_bond=['S'], break_bond=['S'], increment_radical=['S'], decrement_radical=['S'], increment_lone_pair=['S'], decrement_lone_pair=['S']) -ATOMTYPES['S0sc'].set_actions(increment_bond=['S0sc'], decrement_bond=['S0sc'], form_bond=['S0sc'], break_bond=['Sa', 'S0sc'], increment_radical=['S0sc'], decrement_radical=['S0sc'], increment_lone_pair=[], decrement_lone_pair=['S2s', 'S2sc', 'S2dc', 'S2tc']) -ATOMTYPES['Sa'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['S0sc'], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['S2s']) -ATOMTYPES['S2s'].set_actions(increment_bond=['S2d', 'S2dc'], decrement_bond=[], form_bond=['S2s', 'S2sc'], break_bond=['S2s'], increment_radical=['S2s'], decrement_radical=['S2s'], increment_lone_pair=['Sa', 'S0sc'], decrement_lone_pair=['S4s', 'S4sc']) -ATOMTYPES['S2sc'].set_actions(increment_bond=['S2dc'], decrement_bond=[], form_bond=['S2sc'], break_bond=['S2sc', 'S2s'], increment_radical=['S2sc'], decrement_radical=['S2sc'], increment_lone_pair=['S0sc'], decrement_lone_pair=['S4s', 'S4sc']) -ATOMTYPES['S2d'].set_actions(increment_bond=['S2tc'], decrement_bond=['S2s'], form_bond=['S2d'], break_bond=['S2d'], increment_radical=['S2d'], decrement_radical=['S2d'], increment_lone_pair=[], decrement_lone_pair=['S4dc', 'S4d']) -ATOMTYPES['S2dc'].set_actions(increment_bond=['S2tc', 'S2dc'], decrement_bond=['S2sc', 'S2s', 'S2dc'], form_bond=['S2dc'], break_bond=['S2dc'], increment_radical=['S2dc'], decrement_radical=['S2dc'], increment_lone_pair=['S0sc'], decrement_lone_pair=['S4d', 'S4dc']) -ATOMTYPES['S2tc'].set_actions(increment_bond=[], decrement_bond=['S2d', 'S2dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=['S0sc'], decrement_lone_pair=['S4t']) -ATOMTYPES['S4s'].set_actions(increment_bond=['S4d', 'S4dc'], decrement_bond=[], form_bond=['S4s'], break_bond=['S4s'], increment_radical=['S4s'], decrement_radical=['S4s'], increment_lone_pair=['S2s', 'S2sc'], decrement_lone_pair=['S6s']) -ATOMTYPES['S4sc'].set_actions(increment_bond=['S4d', 'S4dc'], decrement_bond=[], form_bond=['S4s', 'S4sc'], break_bond=['S4sc'], increment_radical=['S4sc'], decrement_radical=['S4sc'], increment_lone_pair=['S2s', 'S2sc'], decrement_lone_pair=['S6s']) -ATOMTYPES['S4d'].set_actions(increment_bond=['S4dd', 'S4dc', 'S4t', 'S4tdc'], decrement_bond=['S4s', 'S4sc'], form_bond=['S4dc', 'S4d'], break_bond=['S4d', 'S4dc'], increment_radical=['S4d'], decrement_radical=['S4d'], increment_lone_pair=['S2d', 'S2dc'], decrement_lone_pair=['S6d', 'S6dc']) -ATOMTYPES['S4dc'].set_actions(increment_bond=['S4dd', 'S4dc', 'S4tdc'], decrement_bond=['S4sc', 'S4dc'], form_bond=['S4d', 'S4dc'], break_bond=['S4d', 'S4dc'], increment_radical=['S4dc'], decrement_radical=['S4dc'], increment_lone_pair=['S2d', 'S2dc'], decrement_lone_pair=['S6d', 'S6dc']) -ATOMTYPES['S4b'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['S4dd'].set_actions(increment_bond=['S4dc'], decrement_bond=['S4dc', 'S4d'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['S6dd']) -ATOMTYPES['S4t'].set_actions(increment_bond=[], decrement_bond=['S4d'], form_bond=['S4t'], break_bond=['S4t'], increment_radical=['S4t'], decrement_radical=['S4t'], increment_lone_pair=['S2tc'], decrement_lone_pair=['S6t', 'S6tdc']) -ATOMTYPES['S4tdc'].set_actions(increment_bond=['S4tdc'], decrement_bond=['S4d', 'S4tdc'], form_bond=['S4tdc'], break_bond=['S4tdc'], increment_radical=['S4tdc'], decrement_radical=['S4tdc'], increment_lone_pair=['S6tdc'], decrement_lone_pair=['S6td', 'S6tdc']) -ATOMTYPES['S6s'].set_actions(increment_bond=['S6d', 'S6dc'], decrement_bond=[], form_bond=['S6s'], break_bond=['S6s'], increment_radical=['S6s'], decrement_radical=['S6s'], increment_lone_pair=['S4s', 'S4sc'], decrement_lone_pair=[]) -ATOMTYPES['S6sc'].set_actions(increment_bond=['S6dc'], decrement_bond=[], form_bond=['S6sc'], break_bond=['S6sc'], increment_radical=['S6sc'], decrement_radical=['S6sc'], increment_lone_pair=['S4s', 'S4sc'], decrement_lone_pair=[]) -ATOMTYPES['S6d'].set_actions(increment_bond=['S6dd', 'S6t', 'S6tdc'], decrement_bond=['S6s'], form_bond=['S6d', 'S6dc'], break_bond=['S6d', 'S6dc'], increment_radical=['S6d'], decrement_radical=['S6d'], increment_lone_pair=['S4d', 'S4dc'], decrement_lone_pair=[]) -ATOMTYPES['S6dc'].set_actions(increment_bond=['S6dd', 'S6ddd', 'S6dc', 'S6t', 'S6td', 'S6tdc'], decrement_bond=['S6sc', 'S6dc'], form_bond=['S6d', 'S6dc'], break_bond=['S6d', 'S6dc'], increment_radical=['S6dc'], decrement_radical=['S6dc'], increment_lone_pair=['S4d', 'S4dc'], decrement_lone_pair=[]) -ATOMTYPES['S6dd'].set_actions(increment_bond=['S6ddd', 'S6td'], decrement_bond=['S6d', 'S6dc'], form_bond=['S6dd', 'S6dc'], break_bond=['S6dd'], increment_radical=['S6dd'], decrement_radical=['S6dd'], increment_lone_pair=['S4dd'], decrement_lone_pair=[]) -ATOMTYPES['S6ddd'].set_actions(increment_bond=[], decrement_bond=['S6dd', 'S6dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['S6t'].set_actions(increment_bond=['S6td'], decrement_bond=['S6d', 'S6dc'], form_bond=['S6t'], break_bond=['S6t'], increment_radical=['S6t'], decrement_radical=['S6t'], increment_lone_pair=['S4t'], decrement_lone_pair=[]) -ATOMTYPES['S6td'].set_actions(increment_bond=['S6tt', 'S6tdc'], decrement_bond=['S6dc', 'S6t', 'S6dd', 'S6tdc'], form_bond=['S6td'], break_bond=['S6td'], increment_radical=['S6td'], decrement_radical=['S6td'], increment_lone_pair=['S4tdc'], decrement_lone_pair=[]) -ATOMTYPES['S6tt'].set_actions(increment_bond=[], decrement_bond=['S6td', 'S6tdc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['S6tdc'].set_actions(increment_bond=['S6td', 'S6tdc', 'S6tt'], decrement_bond=['S6dc', 'S6tdc'], form_bond=['S6tdc'], break_bond=['S6tdc'], increment_radical=['S6tdc'], decrement_radical=['S6tdc'], increment_lone_pair=['S4t', 'S4tdc'], decrement_lone_pair=[]) - -ATOMTYPES['Cl'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Cl'], break_bond=['Cl'], increment_radical=['Cl'], decrement_radical=['Cl'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Cl1s'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Cl1s'], break_bond=['Cl1s'], increment_radical=['Cl1s'], decrement_radical=['Cl1s'], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['Br'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Br'], break_bond=['Br'], increment_radical=['Br'], decrement_radical=['Br'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['Br1s'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Br1s'], break_bond=['Br1s'], increment_radical=['Br1s'], decrement_radical=['Br1s'], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['I'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['I'], break_bond=['I'], increment_radical=['I'], decrement_radical=['I'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['I1s'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['I1s'], break_bond=['I1s'], increment_radical=['I1s'], decrement_radical=['I1s'], increment_lone_pair=[], decrement_lone_pair=[]) - -ATOMTYPES['F'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['F'], break_bond=['F'], increment_radical=['F'], decrement_radical=['F'], increment_lone_pair=[], decrement_lone_pair=[]) -ATOMTYPES['F1s'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['F1s'], break_bond=['F1s'], increment_radical=['F1s'], decrement_radical=['F1s'], increment_lone_pair=[], decrement_lone_pair=[]) + +ATOMTYPES['e'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['X'].set_actions(increment_bond=['X'], decrement_bond=['X'], form_bond=['X'], break_bond=['X'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Xv'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Xo'], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Xo'].set_actions(increment_bond=['Xo'], decrement_bond=['Xo'], form_bond=[], break_bond=['Xv'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['R'].set_actions(increment_bond=['R'], decrement_bond=['R'], form_bond=['R'], break_bond=['R'], increment_radical=['R'], decrement_radical=['R'], increment_lone_pair=['R'], decrement_lone_pair=['R'], increment_charge=['R'], decrement_charge=['R']) +ATOMTYPES['R!H'].set_actions(increment_bond=['R!H'], decrement_bond=['R!H'], form_bond=['R!H'], break_bond=['R!H'], increment_radical=['R!H'], decrement_radical=['R!H'], increment_lone_pair=['R!H'], decrement_lone_pair=['R!H'], increment_charge=['R!H'], decrement_charge=['R!H']) +ATOMTYPES['R!H!Val7'].set_actions(increment_bond=['R!H!Val7'], decrement_bond=['R!H!Val7'], form_bond=['R!H!Val7'], break_bond=['R!H!Val7'], increment_radical=['R!H!Val7'], decrement_radical=['R!H!Val7'], increment_lone_pair=['R!H!Val7'], decrement_lone_pair=['R!H!Val7'], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Val4'].set_actions(increment_bond=['Val4'], decrement_bond=['Val4'], form_bond=['Val4'], break_bond=['Val4'], increment_radical=['Val4'], decrement_radical=['Val4'], increment_lone_pair=['Val4'], decrement_lone_pair=['Val4'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['Val5'].set_actions(increment_bond=['Val5'], decrement_bond=['Val5'], form_bond=['Val5'], break_bond=['Val5'], increment_radical=['Val5'], decrement_radical=['Val5'], increment_lone_pair=['Val5'], decrement_lone_pair=['Val5'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['Val6'].set_actions(increment_bond=['Val6'], decrement_bond=['Val6'], form_bond=['Val6'], break_bond=['Val6'], increment_radical=['Val6'], decrement_radical=['Val6'], increment_lone_pair=['Val6'], decrement_lone_pair=['Val6'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['Val7'].set_actions(increment_bond=['Val7'], decrement_bond=['Val7'], form_bond=['Val7'], break_bond=['Val7'], increment_radical=['Val7'], decrement_radical=['Val7'], increment_lone_pair=['Val7'], decrement_lone_pair=['Val7'],increment_charge=[], decrement_charge=[]) + +ATOMTYPES['H'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['H'], break_bond=['H'], increment_radical=['H'], decrement_radical=['H'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=['H'], decrement_charge=['H']) +ATOMTYPES['H0'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['H0'], break_bond=['H0'], increment_radical=['H0'], decrement_radical=['H0'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=['H+'], decrement_charge=[]) +ATOMTYPES['H+'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=['H0']) + +ATOMTYPES['He'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=['He'], decrement_radical=['He'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Ne'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=['Ne'], decrement_radical=['Ne'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Ar'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['C'].set_actions(increment_bond=['C'], decrement_bond=['C'], form_bond=['C'], break_bond=['C'], increment_radical=['C'], decrement_radical=['C'], increment_lone_pair=['C'], decrement_lone_pair=['C'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['Ca'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['C2s'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['Cs'].set_actions(increment_bond=['Cd', 'CO', 'CS'], decrement_bond=[], form_bond=['Cs', 'Csc'], break_bond=['Cs'], increment_radical=['Cs'], decrement_radical=['Cs'], increment_lone_pair=['C2s'], decrement_lone_pair=['C2s'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['Csc'].set_actions(increment_bond=['Cdc'], decrement_bond=[], form_bond=['Csc'], break_bond=['Csc', 'Cs'], increment_radical=['Csc'], decrement_radical=['Csc'], increment_lone_pair=['C2sc'], decrement_lone_pair=['C2sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['Cd'].set_actions(increment_bond=['Cdd', 'Ct', 'C2tc'], decrement_bond=['Cs'], form_bond=['Cd', 'Cdc'], break_bond=['Cd'], increment_radical=['Cd'], decrement_radical=['Cd'], increment_lone_pair=['C2d'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Cdc'].set_actions(increment_bond=[], decrement_bond=['Csc'], form_bond=['Cdc'], break_bond=['Cdc', 'Cd', 'CO', 'CS'], increment_radical=['Cdc'], decrement_radical=['Cdc'], increment_lone_pair=['C2dc'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['CO'].set_actions(increment_bond=['Cdd', 'C2tc'], decrement_bond=['Cs'], form_bond=['CO', 'Cdc'], break_bond=['CO'], increment_radical=['CO'], decrement_radical=['CO'], increment_lone_pair=['C2d'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['CS'].set_actions(increment_bond=['Cdd', 'C2tc'], decrement_bond=['Cs'], form_bond=['CS', 'Cdc'], break_bond=['CS'], increment_radical=['CS'], decrement_radical=['CS'], increment_lone_pair=['C2d'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Cdd'].set_actions(increment_bond=[], decrement_bond=['Cd', 'CO', 'CS'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Ct'].set_actions(increment_bond=['Cq'], decrement_bond=['Cd', 'CO', 'CS'], form_bond=['Ct'], break_bond=['Ct'], increment_radical=['Ct'], decrement_radical=['Ct'], increment_lone_pair=['C2tc'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Cb'].set_actions(increment_bond=['Cbf'], decrement_bond=[], form_bond=['Cb'], break_bond=['Cb'], increment_radical=['Cb'], decrement_radical=['Cb'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Cbf'].set_actions(increment_bond=[], decrement_bond=['Cb'], form_bond=[], break_bond=['Cb'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['C2s'].set_actions(increment_bond=['C2d'], decrement_bond=[], form_bond=['C2s'], break_bond=['C2s'], increment_radical=['C2s'], decrement_radical=['C2s'], increment_lone_pair=['Ca'], decrement_lone_pair=['Cs'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['C2sc'].set_actions(increment_bond=['C2dc'], decrement_bond=[], form_bond=['C2sc'], break_bond=['C2sc'], increment_radical=['C2sc'], decrement_radical=['C2sc'], increment_lone_pair=[], decrement_lone_pair=['Cs'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['C2d'].set_actions(increment_bond=['C2tc'], decrement_bond=['C2s'], form_bond=['C2dc'], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['Cd', 'CO', 'CS'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['C2dc'].set_actions(increment_bond=[], decrement_bond=['C2sc'], form_bond=[], break_bond=['C2d'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['Cdc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['C2tc'].set_actions(increment_bond=[], decrement_bond=['C2d'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['Ct'], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Cq'].set_actions(increment_bond=[], decrement_bond=['Ct'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['N'].set_actions(increment_bond=['N'], decrement_bond=['N'], form_bond=['N'], break_bond=['N'], increment_radical=['N'], decrement_radical=['N'], increment_lone_pair=['N'], decrement_lone_pair=['N'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['N0sc'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['N0sc'], break_bond=['N0sc'], increment_radical=['N0sc'], decrement_radical=['N0sc'], increment_lone_pair=[], decrement_lone_pair=['N1s', 'N1sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['N1s'].set_actions(increment_bond=['N1dc'], decrement_bond=[], form_bond=['N1s'], break_bond=['N1s'], increment_radical=['N1s'], decrement_radical=['N1s'], increment_lone_pair=['N0sc'], decrement_lone_pair=['N3s', 'N3sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['N1sc'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['N1sc'], break_bond=['N1sc'], increment_radical=['N1sc'], decrement_radical=['N1sc'], increment_lone_pair=[], decrement_lone_pair=['N3s', 'N3sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['N1dc'].set_actions(increment_bond=['N1dc'], decrement_bond=['N1s', 'N1dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['N3d'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['N3s'].set_actions(increment_bond=['N3d'], decrement_bond=[], form_bond=['N3s'], break_bond=['N3s'], increment_radical=['N3s'], decrement_radical=['N3s'], increment_lone_pair=['N1s', 'N1sc'], decrement_lone_pair=['N5sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['N3sc'].set_actions(increment_bond=['N3d'], decrement_bond=[], form_bond=['N3sc'], break_bond=['N3sc'], increment_radical=['N3sc'], decrement_radical=['N3sc'], increment_lone_pair=['N1s', 'N1sc'], decrement_lone_pair=['N5sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['N3d'].set_actions(increment_bond=['N3t'], decrement_bond=['N3s', 'N3sc'], form_bond=['N3d'], break_bond=['N3d'], increment_radical=['N3d'], decrement_radical=['N3d'], increment_lone_pair=['N1dc'], decrement_lone_pair=['N5dc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['N3t'].set_actions(increment_bond=[], decrement_bond=['N3d'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['N3b'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['N5sc'].set_actions(increment_bond=['N5dc'], decrement_bond=[], form_bond=['N5sc'], break_bond=['N5sc'], increment_radical=['N5sc'], decrement_radical=['N5sc'], increment_lone_pair=['N3s', 'N3sc'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['N5dc'].set_actions(increment_bond=['N5ddc', 'N5tc'], decrement_bond=['N5sc'], form_bond=['N5dc'], break_bond=['N5dc'], increment_radical=['N5dc'], decrement_radical=['N5dc'], increment_lone_pair=['N3d'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['N5ddc'].set_actions(increment_bond=['N5dddc'], decrement_bond=['N5dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['N5dddc'].set_actions(increment_bond=[], decrement_bond=['N5ddc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['N5tc'].set_actions(increment_bond=[], decrement_bond=['N5dc'], form_bond=['N5tc'], break_bond=['N5tc'], increment_radical=['N5tc'], decrement_radical=['N5tc'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['N5b'].set_actions(increment_bond=['N5bd'], decrement_bond=[], form_bond=['N5b'], break_bond=['N5b'], increment_radical=['N5b'], decrement_radical=['N5b'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['N5bd'].set_actions(increment_bond=[], decrement_bond=['N5b'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['O'].set_actions(increment_bond=['O'], decrement_bond=['O'], form_bond=['O'], break_bond=['O'], increment_radical=['O'], decrement_radical=['O'], increment_lone_pair=['O'], decrement_lone_pair=['O'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['Oa'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['O0sc'], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['O2s', 'O2sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['O0sc'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['O0sc'], break_bond=['Oa', 'O0sc'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['O2s', 'O2sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['O2s'].set_actions(increment_bond=['O2d'], decrement_bond=[], form_bond=['O2s', 'O2sc'], break_bond=['O2s'], increment_radical=['O2s'], decrement_radical=['O2s'], increment_lone_pair=['Oa', 'O0sc'], decrement_lone_pair=['O4sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['O2sc'].set_actions(increment_bond=['O2d'], decrement_bond=[], form_bond=[], break_bond=['O2s'], increment_radical=['O2sc'], decrement_radical=['O2sc'], increment_lone_pair=[], decrement_lone_pair=['O4sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['O2d'].set_actions(increment_bond=[], decrement_bond=['O2s', 'O2sc'], form_bond=[], break_bond=[], increment_radical=['O2d'], decrement_radical=['O2d'], increment_lone_pair=[], decrement_lone_pair=['O4dc', 'O4tc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['O4sc'].set_actions(increment_bond=['O4dc'], decrement_bond=[], form_bond=['O4sc'], break_bond=['O4sc'], increment_radical=['O4sc'], decrement_radical=['O4sc'], increment_lone_pair=['O2s', 'O2sc'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['O4dc'].set_actions(increment_bond=['O4tc'], decrement_bond=['O4sc'], form_bond=['O4dc'], break_bond=['O4dc'], increment_radical=['O4dc'], decrement_radical=['O4dc'], increment_lone_pair=['O2d'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['O4tc'].set_actions(increment_bond=[], decrement_bond=['O4dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['O4b'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['Ne'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=['Ne'], decrement_radical=['Ne'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['Si'].set_actions(increment_bond=['Si'], decrement_bond=['Si'], form_bond=['Si'], break_bond=['Si'], increment_radical=['Si'], decrement_radical=['Si'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Sis'].set_actions(increment_bond=['Sid', 'SiO'], decrement_bond=[], form_bond=['Sis'], break_bond=['Sis'], increment_radical=['Sis'], decrement_radical=['Sis'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Sid'].set_actions(increment_bond=['Sidd', 'Sit'], decrement_bond=['Sis'], form_bond=['Sid'], break_bond=['Sid'], increment_radical=['Sid'], decrement_radical=['Sid'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Sidd'].set_actions(increment_bond=[], decrement_bond=['Sid', 'SiO'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Sit'].set_actions(increment_bond=['Siq'], decrement_bond=['Sid'], form_bond=['Sit'], break_bond=['Sit'], increment_radical=['Sit'], decrement_radical=['Sit'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['SiO'].set_actions(increment_bond=['Sidd'], decrement_bond=['Sis'], form_bond=['SiO'], break_bond=['SiO'], increment_radical=['SiO'], decrement_radical=['SiO'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Sib'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Sib'], break_bond=['Sib'], increment_radical=['Sib'], decrement_radical=['Sib'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Sibf'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Siq'].set_actions(increment_bond=[], decrement_bond=['Sit'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['P'].set_actions(increment_bond=['P'], decrement_bond=['P'], form_bond=['P'], break_bond=['P'], increment_radical=['P'], decrement_radical=['P'], increment_lone_pair=['P'], decrement_lone_pair=['P'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['P0sc'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['P0sc'], break_bond=['P0sc'], increment_radical=['P0sc'], decrement_radical=['P0sc'], increment_lone_pair=[], decrement_lone_pair=['P1s', 'P1sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['P1s'].set_actions(increment_bond=['P1dc'], decrement_bond=[], form_bond=['P1s'], break_bond=['P1s'], increment_radical=['P1s'], decrement_radical=['P1s'], increment_lone_pair=['P0sc'], decrement_lone_pair=['P3s'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['P1sc'].set_actions(increment_bond=['P1dc'], decrement_bond=[], form_bond=['P1sc'], break_bond=['P1sc'], increment_radical=['P1sc'], decrement_radical=['P1sc'], increment_lone_pair=['P0sc'], decrement_lone_pair=['P3s'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['P1dc'].set_actions(increment_bond=[], decrement_bond=['P1s'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['P3d'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['P3s'].set_actions(increment_bond=['P3d'], decrement_bond=[], form_bond=['P3s'], break_bond=['P3s'], increment_radical=['P3s'], decrement_radical=['P3s'], increment_lone_pair=['P1s', 'P1sc'], decrement_lone_pair=['P5s', 'P5sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['P3d'].set_actions(increment_bond=['P3t'], decrement_bond=['P3s'], form_bond=['P3d'], break_bond=['P3d'], increment_radical=['P3d'], decrement_radical=['P3d'], increment_lone_pair=['P1dc'], decrement_lone_pair=['P5d', 'P5dc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['P3t'].set_actions(increment_bond=[], decrement_bond=['P3d'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['P5t', 'P5tc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['P3b'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['P5s'].set_actions(increment_bond=['P5d', 'P5dc'], decrement_bond=[], form_bond=['P5s'], break_bond=['P5s'], increment_radical=['P5s'], decrement_radical=['P5s'], increment_lone_pair=['P3s'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['P5sc'].set_actions(increment_bond=['P5dc'], decrement_bond=[], form_bond=['P5sc'], break_bond=['P5sc'], increment_radical=['P5sc'], decrement_radical=['P5sc'], increment_lone_pair=['P3s'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['P5d'].set_actions(increment_bond=['P5dd', 'P5ddc', 'P5t', 'P5tc'], decrement_bond=['P5s'], form_bond=['P5d'], break_bond=['P5d'], increment_radical=['P5d'], decrement_radical=['P5d'], increment_lone_pair=['P3d'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['P5dd'].set_actions(increment_bond=['P5td'], decrement_bond=['P5d', 'P5dc'], form_bond=['P5dd'], break_bond=['P5dd'], increment_radical=['P5dd'], decrement_radical=['P5dd'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['P5dc'].set_actions(increment_bond=['P5dd', 'P5ddc', 'P5tc'], decrement_bond=['P5sc'], form_bond=['P5dc'], break_bond=['P5dc'], increment_radical=['P5dc'], decrement_radical=['P5dc'], increment_lone_pair=['P3d'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['P5ddc'].set_actions(increment_bond=[], decrement_bond=['P5dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['P5t'].set_actions(increment_bond=['P5td'], decrement_bond=['P5d'], form_bond=['P5t'], break_bond=['P5t'], increment_radical=['P5t'], decrement_radical=['P5t'], increment_lone_pair=['P3t'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['P5td'].set_actions(increment_bond=[], decrement_bond=['P5t', 'P5dd'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['P5tc'].set_actions(increment_bond=[], decrement_bond=['P5dc'], form_bond=['P5tc'], break_bond=['P5tc'], increment_radical=['P5tc'], decrement_radical=['P5tc'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['P5b'].set_actions(increment_bond=['P5bd'], decrement_bond=[], form_bond=['P5b'], break_bond=['P5b'], increment_radical=['P5b'], decrement_radical=['P5b'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['P5bd'].set_actions(increment_bond=[], decrement_bond=['P5b'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['S'].set_actions(increment_bond=['S'], decrement_bond=['S'], form_bond=['S'], break_bond=['S'], increment_radical=['S'], decrement_radical=['S'], increment_lone_pair=['S'], decrement_lone_pair=['S'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S0sc'].set_actions(increment_bond=['S0sc'], decrement_bond=['S0sc'], form_bond=['S0sc'], break_bond=['Sa', 'S0sc'], increment_radical=['S0sc'], decrement_radical=['S0sc'], increment_lone_pair=[], decrement_lone_pair=['S2s', 'S2sc', 'S2dc', 'S2tc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['Sa'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['S0sc'], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['S2s'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S2s'].set_actions(increment_bond=['S2d', 'S2dc'], decrement_bond=[], form_bond=['S2s', 'S2sc'], break_bond=['S2s'], increment_radical=['S2s'], decrement_radical=['S2s'], increment_lone_pair=['Sa', 'S0sc'], decrement_lone_pair=['S4s', 'S4sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S2sc'].set_actions(increment_bond=['S2dc'], decrement_bond=[], form_bond=['S2sc'], break_bond=['S2sc', 'S2s'], increment_radical=['S2sc'], decrement_radical=['S2sc'], increment_lone_pair=['S0sc'], decrement_lone_pair=['S4s', 'S4sc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S2d'].set_actions(increment_bond=['S2tc'], decrement_bond=['S2s'], form_bond=['S2d'], break_bond=['S2d'], increment_radical=['S2d'], decrement_radical=['S2d'], increment_lone_pair=[], decrement_lone_pair=['S4dc', 'S4d'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S2dc'].set_actions(increment_bond=['S2tc', 'S2dc'], decrement_bond=['S2sc', 'S2s', 'S2dc'], form_bond=['S2dc'], break_bond=['S2dc'], increment_radical=['S2dc'], decrement_radical=['S2dc'], increment_lone_pair=['S0sc'], decrement_lone_pair=['S4d', 'S4dc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S2tc'].set_actions(increment_bond=[], decrement_bond=['S2d', 'S2dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=['S0sc'], decrement_lone_pair=['S4t'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S4s'].set_actions(increment_bond=['S4d', 'S4dc'], decrement_bond=[], form_bond=['S4s'], break_bond=['S4s'], increment_radical=['S4s'], decrement_radical=['S4s'], increment_lone_pair=['S2s', 'S2sc'], decrement_lone_pair=['S6s'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S4sc'].set_actions(increment_bond=['S4d', 'S4dc'], decrement_bond=[], form_bond=['S4s', 'S4sc'], break_bond=['S4sc'], increment_radical=['S4sc'], decrement_radical=['S4sc'], increment_lone_pair=['S2s', 'S2sc'], decrement_lone_pair=['S6s'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S4d'].set_actions(increment_bond=['S4dd', 'S4dc', 'S4t', 'S4tdc'], decrement_bond=['S4s', 'S4sc'], form_bond=['S4dc', 'S4d'], break_bond=['S4d', 'S4dc'], increment_radical=['S4d'], decrement_radical=['S4d'], increment_lone_pair=['S2d', 'S2dc'], decrement_lone_pair=['S6d', 'S6dc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S4dc'].set_actions(increment_bond=['S4dd', 'S4dc', 'S4tdc'], decrement_bond=['S4sc', 'S4dc'], form_bond=['S4d', 'S4dc'], break_bond=['S4d', 'S4dc'], increment_radical=['S4dc'], decrement_radical=['S4dc'], increment_lone_pair=['S2d', 'S2dc'], decrement_lone_pair=['S6d', 'S6dc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S4b'].set_actions(increment_bond=[], decrement_bond=[], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['S4dd'].set_actions(increment_bond=['S4dc'], decrement_bond=['S4dc', 'S4d'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=['S6dd'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S4t'].set_actions(increment_bond=[], decrement_bond=['S4d'], form_bond=['S4t'], break_bond=['S4t'], increment_radical=['S4t'], decrement_radical=['S4t'], increment_lone_pair=['S2tc'], decrement_lone_pair=['S6t', 'S6tdc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S4tdc'].set_actions(increment_bond=['S4tdc'], decrement_bond=['S4d', 'S4tdc'], form_bond=['S4tdc'], break_bond=['S4tdc'], increment_radical=['S4tdc'], decrement_radical=['S4tdc'], increment_lone_pair=['S6tdc'], decrement_lone_pair=['S6td', 'S6tdc'],increment_charge=[], decrement_charge=[]) +ATOMTYPES['S6s'].set_actions(increment_bond=['S6d', 'S6dc'], decrement_bond=[], form_bond=['S6s'], break_bond=['S6s'], increment_radical=['S6s'], decrement_radical=['S6s'], increment_lone_pair=['S4s', 'S4sc'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['S6sc'].set_actions(increment_bond=['S6dc'], decrement_bond=[], form_bond=['S6sc'], break_bond=['S6sc'], increment_radical=['S6sc'], decrement_radical=['S6sc'], increment_lone_pair=['S4s', 'S4sc'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['S6d'].set_actions(increment_bond=['S6dd', 'S6t', 'S6tdc'], decrement_bond=['S6s'], form_bond=['S6d', 'S6dc'], break_bond=['S6d', 'S6dc'], increment_radical=['S6d'], decrement_radical=['S6d'], increment_lone_pair=['S4d', 'S4dc'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['S6dc'].set_actions(increment_bond=['S6dd', 'S6ddd', 'S6dc', 'S6t', 'S6td', 'S6tdc'], decrement_bond=['S6sc', 'S6dc'], form_bond=['S6d', 'S6dc'], break_bond=['S6d', 'S6dc'], increment_radical=['S6dc'], decrement_radical=['S6dc'], increment_lone_pair=['S4d', 'S4dc'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['S6dd'].set_actions(increment_bond=['S6ddd', 'S6td'], decrement_bond=['S6d', 'S6dc'], form_bond=['S6dd', 'S6dc'], break_bond=['S6dd'], increment_radical=['S6dd'], decrement_radical=['S6dd'], increment_lone_pair=['S4dd'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['S6ddd'].set_actions(increment_bond=[], decrement_bond=['S6dd', 'S6dc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['S6t'].set_actions(increment_bond=['S6td'], decrement_bond=['S6d', 'S6dc'], form_bond=['S6t'], break_bond=['S6t'], increment_radical=['S6t'], decrement_radical=['S6t'], increment_lone_pair=['S4t'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['S6td'].set_actions(increment_bond=['S6tt', 'S6tdc'], decrement_bond=['S6dc', 'S6t', 'S6dd', 'S6tdc'], form_bond=['S6td'], break_bond=['S6td'], increment_radical=['S6td'], decrement_radical=['S6td'], increment_lone_pair=['S4tdc'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['S6tt'].set_actions(increment_bond=[], decrement_bond=['S6td', 'S6tdc'], form_bond=[], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['S6tdc'].set_actions(increment_bond=['S6td', 'S6tdc', 'S6tt'], decrement_bond=['S6dc', 'S6tdc'], form_bond=['S6tdc'], break_bond=['S6tdc'], increment_radical=['S6tdc'], decrement_radical=['S6tdc'], increment_lone_pair=['S4t', 'S4tdc'], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['Cl'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Cl'], break_bond=['Cl'], increment_radical=['Cl'], decrement_radical=['Cl'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Cl1s'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Cl1s'], break_bond=['Cl1s'], increment_radical=['Cl1s'], decrement_radical=['Cl1s'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['Br'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Br'], break_bond=['Br'], increment_radical=['Br'], decrement_radical=['Br'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['Br1s'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Br1s'], break_bond=['Br1s'], increment_radical=['Br1s'], decrement_radical=['Br1s'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['I'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['I'], break_bond=['I'], increment_radical=['I'], decrement_radical=['I'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['I1s'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['I1s'], break_bond=['I1s'], increment_radical=['I1s'], decrement_radical=['I1s'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) + +ATOMTYPES['F'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['F'], break_bond=['F'], increment_radical=['F'], decrement_radical=['F'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) +ATOMTYPES['F1s'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['F1s'], break_bond=['F1s'], increment_radical=['F1s'], decrement_radical=['F1s'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[]) # these are ordered in priority of picking if a more general atomtype is encountered -allElements = ['H', 'C', 'O', 'N', 'S', 'P', 'Si', 'F', 'Cl', 'Br', 'I', 'Ne', 'Ar', 'He', 'X'] +allElements = ['H', 'C', 'O', 'N', 'S', 'P', 'Si', 'F', 'Cl', 'Br', 'I', 'Ne', 'Ar', 'He', 'X', 'e'] # list of elements that do not have more specific atomTypes -nonSpecifics = ['H', 'He', 'Ne', 'Ar',] +nonSpecifics = ['He', 'Ne', 'Ar', 'e'] for atomtype in ATOMTYPES.values(): for items in [atomtype.generic, atomtype.specific, atomtype.increment_bond, atomtype.decrement_bond, atomtype.form_bond, atomtype.break_bond, atomtype.increment_radical, atomtype.decrement_radical, atomtype.increment_lone_pair, - atomtype.decrement_lone_pair]: + atomtype.decrement_lone_pair, atomtype.increment_charge, atomtype.decrement_charge]: for index in range(len(items)): items[index] = ATOMTYPES[items[index]] - def get_features(atom, bonds): """ Returns a list of features needed to determine atomtype for :class:'Atom' diff --git a/rmgpy/molecule/atomtypeTest.py b/rmgpy/molecule/atomtypeTest.py index 9fd6e5640d..ff97d231b0 100644 --- a/rmgpy/molecule/atomtypeTest.py +++ b/rmgpy/molecule/atomtypeTest.py @@ -83,6 +83,11 @@ def test_pickle(self): self.assertEqual(len(self.atomtype.decrement_radical), len(atom_type.decrement_radical)) for item1, item2 in zip(self.atomtype.decrement_radical, atom_type.decrement_radical): self.assertEqual(item1.label, item2.label) + for item1, item2 in zip(self.atomtype.increment_charge, atom_type.increment_charge): + self.assertEqual(item1.label, item2.label) + self.assertEqual(len(self.atomtype.decrement_charge), len(atom_type.decrement_charge)) + for item1, item2 in zip(self.atomtype.decrement_charge, atom_type.decrement_charge): + self.assertEqual(item1.label, item2.label) def test_output(self): """ @@ -120,13 +125,17 @@ def test_set_actions(self): self.atomtype.increment_radical, self.atomtype.decrement_radical, self.atomtype.increment_lone_pair, - self.atomtype.decrement_lone_pair) + self.atomtype.decrement_lone_pair, + self.atomtype.increment_charge, + self.atomtype.decrement_charge) self.assertEqual(self.atomtype.increment_bond, other.increment_bond) self.assertEqual(self.atomtype.decrement_bond, other.decrement_bond) self.assertEqual(self.atomtype.form_bond, other.form_bond) self.assertEqual(self.atomtype.break_bond, other.break_bond) self.assertEqual(self.atomtype.increment_radical, other.increment_radical) self.assertEqual(self.atomtype.decrement_radical, other.decrement_radical) + self.assertEqual(self.atomtype.increment_charge, other.increment_charge) + self.assertEqual(self.atomtype.decrement_charge, other.decrement_charge) ################################################################################ @@ -610,6 +619,9 @@ def setUp(self): 11 H u0 p0 {5,S} 12 H u0 p0 {6,S}''') + self.electron = Molecule().from_adjacency_list('''1 e u1 p0 c-1''') + self.proton = Molecule().from_adjacency_list('''1 H u0 p0 c+1''') + def atom_type(self, mol, atom_id): atom = mol.atoms[atom_id] atom_type = get_atomtype(atom, mol.get_bonds(atom)) @@ -622,7 +634,7 @@ def test_hydrogen_type(self): """ Test that get_atomtype() returns the hydrogen atom type. """ - self.assertEqual(self.atom_type(self.mol3, 0), 'H') + self.assertEqual(self.atom_type(self.mol3, 0), 'H0') def test_carbon_types(self): """ @@ -788,7 +800,7 @@ def test_occupied_surface_atom_type(self): """ Test that get_atomtype() works for occupied surface sites and for regular atoms in the complex. """ - self.assertEqual(self.atom_type(self.mol76, 0), 'H') + self.assertEqual(self.atom_type(self.mol76, 0), 'H0') self.assertEqual(self.atom_type(self.mol76, 1), 'Xo') def test_vacant_surface_site_atom_type(self): @@ -796,10 +808,22 @@ def test_vacant_surface_site_atom_type(self): Test that get_atomtype() works for vacant surface sites and for regular atoms in the complex. """ self.assertEqual(self.atom_type(self.mol77, 0), 'Cs') - self.assertEqual(self.atom_type(self.mol77, 1), 'H') + self.assertEqual(self.atom_type(self.mol77, 1), 'H0') self.assertEqual(self.atom_type(self.mol77, 3), 'Xv') self.assertEqual(self.atom_type(self.mol78, 0), 'Xv') + def test_electron(self): + """ + Test that get_atomtype() returns the electron (e) atom type. + """ + self.assertEqual(self.atom_type(self.electron, 0), 'e') + + def test_proton(self): + """ + Test that get_atomtype() returns the proton (H+) atom type. + """ + self.assertEqual(self.atom_type(self.proton, 0), 'H+') + ################################################################################ diff --git a/rmgpy/molecule/draw.py b/rmgpy/molecule/draw.py index e6b37e315c..f6b5135364 100644 --- a/rmgpy/molecule/draw.py +++ b/rmgpy/molecule/draw.py @@ -1228,7 +1228,7 @@ def _render_atom(self, symbol, atom, x0, y0, cr, heavy_first=True, draw_lone_pai heavy_atom = symbol[0] # Split label by atoms - labels = re.findall('[A-Z][a-z]*[0-9]*', symbol) + labels = re.findall('[A-Z]*[a-z]*[0-9]*', symbol) if not heavy_first: labels.reverse() if 'C' not in symbol and 'O' not in symbol and len(atoms) == 1: @@ -1333,6 +1333,8 @@ def _render_atom(self, symbol, atom, x0, y0, cr, heavy_first=True, draw_lone_pai cr.set_source_rgba(0.5, 0.0, 0.5, 1.0) elif heavy_atom == 'X': cr.set_source_rgba(0.5, 0.25, 0.5, 1.0) + elif heavy_atom == 'e': + cr.set_source_rgba(1.0, 0.0, 1.0, 1.0) else: cr.set_source_rgba(0.0, 0.0, 0.0, 1.0) diff --git a/rmgpy/molecule/element.py b/rmgpy/molecule/element.py index 88e2408e99..af67d16657 100644 --- a/rmgpy/molecule/element.py +++ b/rmgpy/molecule/element.py @@ -78,7 +78,7 @@ def __init__(self, number, symbol, name, mass, isotope=-1, chemkin_name=None): self.mass = mass self.isotope = isotope self.chemkin_name = chemkin_name or self.name - if symbol == 'X': + if symbol in ('X','e'): self.cov_radius = 0 else: try: @@ -122,11 +122,11 @@ class PeriodicSystem(object): https://sciencenotes.org/list-of-electronegativity-values-of-the-elements/ isotopes of the same element may have slight different electronegativities, which is not reflected below """ - valences = {'H': 1, 'He': 0, 'C': 4, 'N': 3, 'O': 2, 'F': 1, 'Ne': 0, + valences = {'H+':0, 'e': 0, 'H': 1, 'He': 0, 'C': 4, 'N': 3, 'O': 2, 'F': 1, 'Ne': 0, 'Si': 4, 'P': 3, 'S': 2, 'Cl': 1, 'Br': 1, 'Ar': 0, 'I': 1, 'X': 4} - valence_electrons = {'H': 1, 'He': 2, 'C': 4, 'N': 5, 'O': 6, 'F': 7, 'Ne': 8, + valence_electrons = {'H+':0, 'e': 1, 'H': 1, 'He': 2, 'C': 4, 'N': 5, 'O': 6, 'F': 7, 'Ne': 8, 'Si': 4, 'P': 5, 'S': 6, 'Cl': 7, 'Br': 7, 'Ar': 8, 'I': 7, 'X': 4} - lone_pairs = {'H': 0, 'He': 1, 'C': 0, 'N': 1, 'O': 2, 'F': 3, 'Ne': 4, + lone_pairs = {'H+':0, 'e': 0, 'H': 0, 'He': 1, 'C': 0, 'N': 1, 'O': 2, 'F': 3, 'Ne': 4, 'Si': 0, 'P': 1, 'S': 2, 'Cl': 3, 'Br': 3, 'Ar': 4, 'I': 3, 'X': 0} electronegativity = {'H': 2.20, 'D': 2.20, 'T': 2.20, 'C': 2.55, 'C13': 2.55, 'N': 3.04, 'O': 3.44, 'O18': 3.44, 'F': 3.98, 'Si': 1.90, 'P': 2.19, 'S': 2.58, 'Cl': 3.16, 'Br': 2.96, 'I': 2.66, 'X': 0.0} @@ -173,6 +173,8 @@ def get_element(value, isotope=-1): # Recommended IUPAC nomenclature is used throughout (including 'aluminium' and # 'caesium') +# electron +e = Element(-1, 'e', 'electron' , 5.486e-7) # Surface site X = Element(0, 'X', 'surface_site' , 0.0) @@ -310,6 +312,7 @@ def get_element(value, isotope=-1): # A list of the elements, sorted by increasing atomic number element_list = [ + e, X, H, D, T, He, Li, Be, B, C, C13, N, O, O18, F, Ne, diff --git a/rmgpy/molecule/group.pxd b/rmgpy/molecule/group.pxd index d75951f8e5..a1b138a914 100644 --- a/rmgpy/molecule/group.pxd +++ b/rmgpy/molecule/group.pxd @@ -71,6 +71,10 @@ cdef class GroupAtom(Vertex): cpdef bint is_bonded_to_surface(self) except -2 + cpdef bint is_proton(self) + + cpdef bint is_electron(self) + cpdef bint is_oxygen(self) cpdef bint is_sulfur(self) @@ -186,6 +190,10 @@ cdef class Group(Graph): cpdef bint is_surface_site(self) except -2 + cpdef bint is_proton(self) + + cpdef bint is_electron(self) + cpdef bint contains_surface_site(self) except -2 cpdef list get_surface_sites(self) diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py index 01443c1e4d..257b05b8a6 100644 --- a/rmgpy/molecule/group.py +++ b/rmgpy/molecule/group.py @@ -260,6 +260,54 @@ def _lose_radical(self, radical): # Set the new radical electron counts self.radical_electrons = radical_electrons + def _gain_charge(self, charge): + """ + Update the atom group as a result of applying a GAIN_CHARGE action, + where `charge` specifies the charge gained. + """ + atomtype = [] + + for atom in self.atomtype: + atomtype.extend(atom.increment_charge) + + if any([len(atom.increment_charge) == 0 for atom in self.atomtype]): + raise ActionError('Unable to update GroupAtom due to GAIN_CHARGE action: ' + 'Unknown atom type produced from set "{0}".'.format(self.atomtype)) + + if isinstance(self.charge,list): + charges = [] + for c in self.charge: + charges.append(c+charge) + self.charge = charges + else: + self.charge += 1 + + self.atomtype = list(set(atomtype)) + + def _lose_charge(self, charge): + """ + Update the atom group as a result of applying a LOSE_CHARGE action, + where `charge` specifies lost charge. + """ + atomtype = [] + + for atom in self.atomtype: + atomtype.extend(atom.decrement_charge) + + if any([len(atomtype.decrement_charge) == 0 for atomtype in self.atomtype]): + raise ActionError('Unable to update GroupAtom due to LOSE_CHARGE action: ' + 'Unknown atom type produced from set "{0}".'.format(self.atomtype)) + + if isinstance(self.charge,list): + charges = [] + for c in self.charge: + charges.append(c-charge) + self.charge = charges + else: + self.charge -= 1 + + self.atomtype = list(set(atomtype)) + def _gain_pair(self, pair): """ Update the atom group as a result of applying a GAIN_PAIR action, @@ -331,8 +379,12 @@ def apply_action(self, action): self._break_bond(action[2]) elif act == 'GAIN_RADICAL': self._gain_radical(action[2]) + elif act == 'GAIN_CHARGE': + self._gain_charge(action[2]) elif act == 'LOSE_RADICAL': self._lose_radical(action[2]) + elif act == 'LOSE_CHARGE': + self._lose_charge(action[2]) elif action[0].upper() == 'GAIN_PAIR': self._gain_pair(action[2]) elif action[0].upper() == 'LOSE_PAIR': @@ -492,6 +544,18 @@ def is_bonded_to_surface(self): return True return False + def is_electron(self): + """ + Return ``True`` if the atom represents a surface site or ``False`` if not. + """ + return self.atomtype[0] == ATOMTYPES['e'] + + def is_proton(self): + """ + Return ``True`` if the atom represents a surface site or ``False`` if not. + """ + return self.atomtype[0] == ATOMTYPES['H+'] + def is_oxygen(self): """ Return ``True`` if the atom represents an oxygen atom or ``False`` if not. @@ -628,6 +692,7 @@ def make_sample_atom(self): 'I': 3, 'Ar': 4, 'X': 0, + 'e': 0 } for element_label in allElements: @@ -1156,6 +1221,14 @@ def get_surface_sites(self): cython.declare(atom=GroupAtom) return [atom for atom in self.atoms if atom.is_surface_site()] + def is_proton(self): + """Returns ``True`` iff the group is a proton""" + return len(self.atoms) == 1 and self.atoms[0].is_proton() + + def is_electron(self): + """Returns ``True`` iff the group is an electron""" + return len(self.atoms) == 1 and self.atoms[0].is_electron() + def remove_atom(self, atom): """ Remove `atom` and all bonds associated with it from the graph. Does @@ -2722,12 +2795,14 @@ def make_sample_molecule(self): else: raise UnexpectedChargeError(graph=new_molecule) # check hardcoded atomtypes - positive_charged = ['Csc', 'Cdc', + positive_charged = ['H+', + 'Csc', 'Cdc', 'N3sc', 'N5sc', 'N5dc', 'N5ddc', 'N5tc', 'N5b', 'O2sc', 'O4sc', 'O4dc', 'O4tc', 'P5sc', 'P5dc', 'P5ddc', 'P5tc', 'P5b', 'S2sc', 'S4sc', 'S4dc', 'S4tdc', 'S6sc', 'S6dc', 'S6tdc'] - negative_charged = ['C2sc', 'C2dc', 'C2tc', + negative_charged = ['e', + 'C2sc', 'C2dc', 'C2tc', 'N0sc', 'N1sc', 'N1dc', 'N5dddc', 'O0sc', 'P0sc', 'P1sc', 'P1dc', 'P5sc', diff --git a/rmgpy/molecule/groupTest.py b/rmgpy/molecule/groupTest.py index 3b3e2b3f54..c77d14704d 100644 --- a/rmgpy/molecule/groupTest.py +++ b/rmgpy/molecule/groupTest.py @@ -189,6 +189,66 @@ def test_apply_action_lose_radical(self): atom1.apply_action(action) self.assertListEqual(atom1.radical_electrons, [0, 1, 2, 3]) + def test_apply_action_gain_charge(self): + """ + Test the GroupAtom.apply_action() method for a GAIN_CHARGE action. + """ + action = ['GAIN_CHARGE', '*1', 1] + for label, atomtype in ATOMTYPES.items(): + atom0 = GroupAtom(atomtype=[atomtype], radical_electrons=[0], charge=[0], label='*1', lone_pairs=[0]) + atom = atom0.copy() + try: + atom.apply_action(action) + self.assertEqual(len(atom.atomtype), len(atomtype.increment_charge)) + for a in atomtype.increment_charge: + self.assertTrue(a in atom.atomtype, + "GAIN_CHARGE on {0} gave {1} not {2}".format(atomtype, atom.atomtype, + atomtype.increment_charge)) + # self.assertEqual(atom0.radical_electrons, [r + 1 for r in atom.radical_electrons]) + self.assertEqual(atom0.charge, [c - 1 for c in atom.charge]) + self.assertEqual(atom0.label, atom.label) + self.assertEqual(atom0.lone_pairs, atom.lone_pairs) + except ActionError: + self.assertEqual(len(atomtype.increment_charge), 0) + + # test when radicals unspecified + group = Group().from_adjacency_list(""" + 1 R ux + """) # ux causes a wildcard for radicals + atom1 = group.atoms[0] + atom1.apply_action(action) + #self.assertListEqual(atom1.radical_electrons, [0, 1, 2, 3]) + + def test_apply_action_lose_charge(self): + """ + Test the GroupAtom.apply_action() method for a LOSE_CHARGE action. + """ + action = ['LOSE_CHARGE', '*1', 1] + for label, atomtype in ATOMTYPES.items(): + atom0 = GroupAtom(atomtype=[atomtype], radical_electrons=[1], charge=[0], label='*1', lone_pairs=[0]) + atom = atom0.copy() + try: + atom.apply_action(action) + self.assertEqual(len(atom.atomtype), len(atomtype.decrement_charge)) + for a in atomtype.decrement_charge: + self.assertTrue(a in atom.atomtype, + "LOSE_CHARGE on {0} gave {1} not {2}".format(atomtype, atom.atomtype, + atomtype.decrement_charge)) + # self.assertEqual(atom0.radical_electrons, [r - 1 for r in atom.radical_electrons]) + self.assertEqual(atom0.charge, [c + 1 for c in atom.charge]) + self.assertEqual(atom0.label, atom.label) + self.assertEqual(atom0.lone_pairs, atom.lone_pairs) + except ActionError: + self.assertEqual(len(atomtype.decrement_charge), 0) + + # test when radicals unspecified + group = Group().from_adjacency_list(""" + 1 R ux + """) # ux causes a wildcard for radicals + atom1 = group.atoms[0] + atom1.apply_action(action) + #self.assertListEqual(atom1.radical_electrons, [1, 2, 3, 4]) + def test_apply_action_gain_pair(self): """ Test the GroupAtom.apply_action() method for a GAIN_PAIR action when lone_pairs is either specified or not. @@ -409,6 +469,19 @@ def test_make_sample_atom(self): self.assertEquals(new_atom.charge, 0) self.assertEquals(new_atom.lone_pairs, 0) + def test_is_electron(self): + """ + Test the GroupAtom.is_electron() method. + """ + electron = GroupAtom(atomtype=[ATOMTYPES['e']]) + self.assertTrue(electron.is_electron()) + + def test_is_proton(self): + """ + Test the GroupAtom.is_proton() method. + """ + proton = GroupAtom(atomtype=[ATOMTYPES['H+']]) + self.assertTrue(proton.is_proton()) ################################################################################ @@ -586,6 +659,35 @@ def test_apply_action_lose_radical(self): except ActionError: pass + + def test_apply_action_gain_charge(self): + """ + Test the GroupBond.apply_action() method for a GAIN_RADICAL action. + """ + action = ['GAIN_CHARGE', '*1', 1] + for order0 in self.orderList: + bond0 = GroupBond(None, None, order=order0) + bond = bond0.copy() + try: + bond.apply_action(action) + self.fail('GroupBond.apply_action() unexpectedly processed a GAIN_CHARGE action.') + except ActionError: + pass + + def test_apply_action_lose_charge(self): + """ + Test the GroupBond.apply_action() method for a LOSE_CHARGE action. + """ + action = ['LOSE_CHARGE', '*1', 1] + for order0 in self.orderList: + bond0 = GroupBond(None, None, order=order0) + bond = bond0.copy() + try: + bond.apply_action(action) + self.fail('GroupBond.apply_action() unexpectedly processed a LOSE_CHARGE action.') + except ActionError: + pass + def test_equivalent(self): """ Test the GroupBond.equivalent() method. @@ -692,6 +794,22 @@ def test_is_surface_site(self): surface_site = Group().from_adjacency_list("1 *1 X u0") self.assertTrue(surface_site.is_surface_site()) + def test_is_electron(self): + """ + Test the Group.is_electron() method. + """ + self.assertFalse(self.group.is_electron()) + electron = Group().from_adjacency_list("""1 *1 e u1 p0 c-1""") + self.assertTrue(electron.is_electron()) + + def test_is_proton(self): + """ + Test the Group.is_proton() method. + """ + self.assertFalse(self.group.is_proton()) + proton = Group().from_adjacency_list("""1 *1 H+ u0 p0 c+1""") + self.assertTrue(proton.is_proton()) + def test_get_labeled_atom(self): """ Test the Group.get_labeled_atoms() method. diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index cb720ef84f..fe65b425a0 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -54,6 +54,10 @@ cdef class Atom(Vertex): cpdef Vertex copy(self) + cpdef bint is_electron(self) + + cpdef bint is_proton(self) + cpdef bint is_hydrogen(self) cpdef bint is_non_hydrogen(self) @@ -85,7 +89,11 @@ cdef class Atom(Vertex): cpdef increment_radical(self) cpdef decrement_radical(self) - + + cpdef increment_charge(self) + + cpdef decrement_charge(self) + cpdef set_lone_pairs(self, int lone_pairs) cpdef increment_lone_pairs(self) @@ -160,6 +168,10 @@ cdef class Molecule(Graph): cpdef bint has_bond(self, Atom atom1, Atom atom2) + cpdef bint is_electron(self) + + cpdef bint is_proton(self) + cpdef bint contains_surface_site(self) cpdef bint is_surface_site(self) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index d0a057bc2f..42042c9c49 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -321,6 +321,23 @@ def copy(self): a.props = deepcopy(self.props) return a + def is_electron(self): + """ + Return ``True`` if the atom represents an electron or ``False`` if + not. + """ + return self.element.number == -1 + + def is_proton(self): + """ + Return ``True`` if the atom represents a proton or ``False`` if + not. + """ + + if self.element.number == 1 and self.charge == 1: + return True + return False + def is_hydrogen(self): """ Return ``True`` if the atom represents a hydrogen atom or ``False`` if @@ -459,6 +476,18 @@ def decrement_radical(self): raise gr.ActionError('Unable to update Atom due to LOSE_RADICAL action: ' 'Invalid radical electron set "{0}".'.format(self.radical_electrons)) + def increment_charge(self): + """ + Update the atom pattern as a result of applying a GAIN_CHARGE action + """ + self.charge += 1 + + def decrement_charge(self): + """ + Update the atom pattern as a result of applying a LOSE_CHARGE action + """ + self.charge -= 1 + def set_lone_pairs(self, lone_pairs): """ Set the number of lone electron pairs. @@ -500,6 +529,10 @@ def update_charge(self): if self.is_surface_site(): self.charge = 0 return + if self.is_electron(): + self.charge = -1 + return + valence_electron = elements.PeriodicSystem.valence_electrons[self.symbol] order = self.get_total_bond_order() self.charge = valence_electron - order - self.radical_electrons - 2 * self.lone_pairs @@ -522,6 +555,10 @@ def apply_action(self, action): for i in range(action[2]): self.increment_radical() elif act == 'LOSE_RADICAL': for i in range(abs(action[2])): self.decrement_radical() + elif act == 'GAIN_CHARGE': + for i in range(action[2]): self.increment_charge() + elif act == 'LOSE_CHARGE': + for i in range(abs(action[2])): self.decrement_charge() elif action[0].upper() == 'GAIN_PAIR': for i in range(action[2]): self.increment_lone_pairs() elif action[0].upper() == 'LOSE_PAIR': @@ -1102,6 +1139,14 @@ def is_surface_site(self): """Returns ``True`` iff the molecule is nothing but a surface site 'X'.""" return len(self.atoms) == 1 and self.atoms[0].is_surface_site() + def is_electron(self): + """Returns ``True`` iff the molecule is nothing but an electron 'e'.""" + return len(self.atoms) == 1 and self.atoms[0].is_electron() + + def is_proton(self): + """Returns ``True`` iff the molecule is nothing but a proton 'H+'.""" + return len(self.atoms) == 1 and self.atoms[0].is_proton() + def remove_atom(self, atom): """ Remove `atom` and all bonds associated with it from the graph. Does @@ -1148,17 +1193,21 @@ def sort_atoms(self): for index, vertex in enumerate(self.vertices): vertex.sorting_label = index + def update_charge(self): + + for atom in self.atoms: + atom.update_charge() + def update(self, log_species=True, raise_atomtype_exception=True, sort_atoms=True): """ - Update the charge and atom types of atoms. + Update the lone_pairs, charge, and atom types of atoms. Update multiplicity, and sort atoms (if ``sort_atoms`` is ``True``) Does not necessarily update the connectivity values (which are used in isomorphism checks) If you need that, call update_connectivity_values() """ - for atom in self.atoms: - atom.update_charge() - + self.update_lone_pairs() + self.update_charge() self.update_atomtypes(log_species=log_species, raise_exception=raise_atomtype_exception) self.update_multiplicity() if sort_atoms: @@ -1720,7 +1769,7 @@ def from_smarts(self, smartsstr, raise_atomtype_exception=True): return self def from_adjacency_list(self, adjlist, saturate_h=False, raise_atomtype_exception=True, - raise_charge_exception=True): + raise_charge_exception=False): """ Convert a string adjacency list `adjlist` to a molecular structure. Skips the first line (assuming it's a label) unless `withLabel` is @@ -1899,7 +1948,7 @@ def find_h_bonds(self): ONinds = [n for n, a in enumerate(self.atoms) if a.is_oxygen() or a.is_nitrogen()] for i, atm1 in enumerate(self.atoms): - if atm1.atomtype.label == 'H': + if atm1.atomtype.label == 'H0': atm_covs = [q for q in atm1.bonds.keys()] if len(atm_covs) > 1: # H is already H bonded continue @@ -2159,10 +2208,12 @@ def is_aryl_radical(self, aromatic_rings=None): def generate_resonance_structures(self, keep_isomorphic=False, filter_structures=True, save_order=False): """Returns a list of resonance structures of the molecule.""" - return resonance.generate_resonance_structures(self, keep_isomorphic=keep_isomorphic, - filter_structures=filter_structures, - save_order=save_order, - ) + + try: + return resonance.generate_resonance_structures(self, keep_isomorphic=keep_isomorphic, filter_structures=filter_structures, save_order=save_order) + except: + logging.warning("Resonance structure generation failed for {}".format(self)) + return [self.copy(deep=True)] def get_url(self): """ @@ -2192,7 +2243,7 @@ def update_lone_pairs(self): """ cython.declare(atom1=Atom, atom2=Atom, bond12=Bond, order=float) for atom1 in self.vertices: - if atom1.is_hydrogen() or atom1.is_surface_site(): + if atom1.is_hydrogen() or atom1.is_surface_site() or atom1.is_electron(): atom1.lone_pairs = 0 else: order = atom1.get_total_bond_order() diff --git a/rmgpy/molecule/moleculeTest.py b/rmgpy/molecule/moleculeTest.py index 576d73c374..2c77504f00 100644 --- a/rmgpy/molecule/moleculeTest.py +++ b/rmgpy/molecule/moleculeTest.py @@ -109,6 +109,31 @@ def test_is_hydrogen(self): else: self.assertFalse(atom.is_hydrogen()) + def test_is_proton(self): + """ + Test the Atom.is_proton() method. + """ + for element in element_list: + atom = Atom(element=element, radical_electrons=0, charge=1, label='*1', lone_pairs=0) + if element.symbol == 'H': + self.assertTrue(atom.is_hydrogen()) + self.assertTrue(atom.is_proton()) + atom.charge = 0 + self.assertFalse(atom.is_proton()) + else: + self.assertFalse(atom.is_proton()) + + def test_is_electron(self): + """ + Test the Atom.is_electron() method. + """ + for element in element_list: + atom = Atom(element=element, radical_electrons=1, charge=-1, label='*1', lone_pairs=0) + if element.symbol == 'e': + self.assertTrue(atom.is_electron()) + else: + self.assertFalse(atom.is_electron()) + def test_is_non_hydrogen(self): """ Test the Atom.is_non_hydrogen() method. @@ -363,6 +388,34 @@ def test_apply_action_lose_radical(self): self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) + def test_apply_action_gain_charge(self): + """ + Test the Atom.apply_action() method for a GAIN_CHARGE action. + """ + action = ['GAIN_CHARGE', '*1', 1] + for element in element_list: + atom0 = Atom(element=element, radical_electrons=0, charge=0, label='*1', lone_pairs=0) + atom = atom0.copy() + atom.apply_action(action) + self.assertEqual(atom0.element, atom.element) + # self.assertEqual(atom0.radical_electrons, atom.radical_electrons + 1) + self.assertEqual(atom0.charge, atom.charge - 1) + self.assertEqual(atom0.label, atom.label) + + def test_apply_action_lose_charge(self): + """ + Test the Atom.apply_action() method for a LOSE_CHARGE action. + """ + action = ['LOSE_CHARGE', '*1', 1] + for element in element_list: + atom0 = Atom(element=element, radical_electrons=0, charge=0, label='*1', lone_pairs=0) + atom = atom0.copy() + atom.apply_action(action) + self.assertEqual(atom0.element, atom.element) + # self.assertEqual(atom0.radical_electrons, atom.radical_electrons - 1) + self.assertEqual(atom0.charge, atom.charge + 1) + self.assertEqual(atom0.label, atom.label) + def test_equivalent(self): """ Test the Atom.equivalent() method. @@ -873,6 +926,18 @@ def setUp(self): self.mol2 = Molecule(smiles='C') self.mol3 = Molecule(smiles='CC') + def test_is_proton(self): + """Test the Molecule `is_proton()` method""" + proton = Molecule().from_adjacency_list("""1 H u0 c+1""") + hydrogen = Molecule().from_adjacency_list("""1 H u1""") + self.assertTrue(proton.is_proton()) + self.assertFalse(hydrogen.is_proton()) + + def test_is_electron(self): + """Test the Molecule `is_electron()` method""" + electron = Molecule().from_adjacency_list("""1 e u1 c-1""") + self.assertTrue(electron.is_electron()) + def test_equality(self): """Test that we can perform equality comparison with Molecule objects""" self.assertEqual(self.mol1, self.mol1) diff --git a/rmgpy/molecule/translator.py b/rmgpy/molecule/translator.py index cc159a3be3..8d44a7fe4c 100644 --- a/rmgpy/molecule/translator.py +++ b/rmgpy/molecule/translator.py @@ -103,7 +103,17 @@ """ multiplicity 1 1 X u0 + """, + 'e': + """ + multiplicity 1 + 1 e u0 p0 c-1 + """, + '[H+]': """ + multiplicity 1 + 1 H u0 p0 c+1 + """, } #: This dictionary is used to shortcut lookups of a molecule's SMILES string from its chemical formula. @@ -128,6 +138,8 @@ 'ClH': 'Cl', 'I2': '[I][I]', 'HI': 'I', + 'H': 'H+', + 'e': 'e' } RADICAL_LOOKUPS = { @@ -155,7 +167,8 @@ 'I': '[I]', 'CF': '[C]F', 'CCl': '[C]Cl', - 'CBr': '[C]Br' + 'CBr': '[C]Br', + 'e': 'e' } @@ -506,7 +519,7 @@ def _read(mol, identifier, identifier_type, backend, raise_atomtype_exception=Tr if _lookup(mol, identifier, identifier_type) is not None: if _check_output(mol, identifier): - mol.update_atomtypes(log_species=True, raise_exception=raise_atomtype_exception) + mol.update(log_species=True, raise_atomtype_exception=raise_atomtype_exception, sort_atoms=False) return mol for option in _get_backend_list(backend): @@ -518,7 +531,7 @@ def _read(mol, identifier, identifier_type, backend, raise_atomtype_exception=Tr raise NotImplementedError("Unrecognized backend {0}".format(option)) if _check_output(mol, identifier): - mol.update_atomtypes(log_species=True, raise_exception=raise_atomtype_exception) + mol.update(log_species=True, raise_atomtype_exception=raise_atomtype_exception, sort_atoms=False) return mol else: logging.debug('Backend {0} is not able to parse identifier {1}'.format(option, identifier)) diff --git a/rmgpy/quantity.py b/rmgpy/quantity.py index c68ccd396e..ade8fadfa6 100644 --- a/rmgpy/quantity.py +++ b/rmgpy/quantity.py @@ -776,6 +776,8 @@ def __call__(self, *args, **kwargs): Momentum = UnitType('kg*m/s^2') +Potential = UnitType('V') + Power = UnitType('W') Pressure = UnitType('Pa', common_units=['bar', 'atm', 'torr', 'psi', 'mbar']) diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd index 586c4a6796..500317abae 100644 --- a/rmgpy/reaction.pxd +++ b/rmgpy/reaction.pxd @@ -31,7 +31,7 @@ from rmgpy.molecule.molecule cimport Atom, Molecule from rmgpy.molecule.element cimport Element from rmgpy.kinetics.model cimport KineticsModel from rmgpy.kinetics.arrhenius cimport Arrhenius -from rmgpy.kinetics.surface cimport SurfaceArrhenius, StickingCoefficient +from rmgpy.kinetics.surface cimport SurfaceArrhenius, StickingCoefficient, SurfaceChargeTransfer cimport numpy as np @@ -48,8 +48,11 @@ cdef class Reaction: cdef public KineticsModel kinetics cdef public Arrhenius network_kinetics cdef public SurfaceArrhenius + cdef public SurfaceChargeTransfer cdef public bint duplicate cdef public float _degeneracy + cdef public int electrons + cdef public int _protons cdef public list pairs cdef public bint allow_pdep_route cdef public bint elementary_high_p @@ -69,6 +72,10 @@ cdef class Reaction: cpdef bint is_surface_reaction(self) + cpdef bint is_charge_transfer_reaction(self) + + cpdef bint is_surface_charge_transfer_reaction(self) + cpdef bint has_template(self, list reactants, list products) cpdef bint matches_species(self, list reactants, list products=?) @@ -76,28 +83,36 @@ cdef class Reaction: cpdef bint is_isomorphic(self, Reaction other, bint either_direction=?, bint check_identical=?, bint check_only_label=?, bint check_template_rxn_products=?, bint generate_initial_map=?, bint strict=?, bint save_order=?) except -2 + + cpdef double _apply_CHE_model(self, double T) cpdef double get_enthalpy_of_reaction(self, double T) cpdef double get_entropy_of_reaction(self, double T) - cpdef double get_free_energy_of_reaction(self, double T) + cpdef double _get_free_energy_of_charge_transfer_reaction(self, double T, double potential=?) + + cpdef double get_free_energy_of_reaction(self, double T, double potential=?) - cpdef double get_equilibrium_constant(self, double T, str type=?, double surface_site_density=?) + cpdef double get_reversible_potential(self, double T) + + cpdef double set_reference_potential(self, double T) + + cpdef double get_equilibrium_constant(self, double T, double potential=?, str type=?, double surface_site_density=?) cpdef np.ndarray get_enthalpies_of_reaction(self, np.ndarray Tlist) cpdef np.ndarray get_entropies_of_reaction(self, np.ndarray Tlist) - cpdef np.ndarray get_free_energies_of_reaction(self, np.ndarray Tlist) + cpdef np.ndarray get_free_energies_of_reaction(self, np.ndarray Tlist, double potential=?) - cpdef np.ndarray get_equilibrium_constants(self, np.ndarray Tlist, str type=?) + cpdef np.ndarray get_equilibrium_constants(self, np.ndarray Tlist, double potential=?, str type=?) cpdef int get_stoichiometric_coefficient(self, Species spec) - cpdef double get_rate_coefficient(self, double T, double P=?, double surface_site_density=?) + cpdef double get_rate_coefficient(self, double T, double P=?, double surface_site_density=?, double potential=?) - cpdef double get_surface_rate_coefficient(self, double T, double surface_site_density) except -2 + cpdef double get_surface_rate_coefficient(self, double T, double surface_site_density, double potential=?) except -2 cpdef fix_barrier_height(self, bint force_positive=?) @@ -107,6 +122,8 @@ cdef class Reaction: cpdef reverse_sticking_coeff_rate(self, StickingCoefficient k_forward, str reverse_units, double surface_site_density, Tmin=?, Tmax=?) + cpdef reverse_surface_charge_transfer_rate(self, SurfaceChargeTransfer k_forward, str reverse_units, Tmin=?, Tmax=?) + cpdef generate_reverse_rate_coefficient(self, bint network_kinetics=?, Tmin=?, Tmax=?, double surface_site_density=?) cpdef np.ndarray calculate_tst_rate_coefficients(self, np.ndarray Tlist) diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index dfb598784f..2f98ef3ede 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -55,12 +55,13 @@ PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, get_rate_coefficient_units_from_reaction_order, \ SurfaceArrheniusBEP, StickingCoefficientBEP from rmgpy.kinetics.arrhenius import Arrhenius # Separate because we cimport from rmgpy.kinetics.arrhenius -from rmgpy.kinetics.surface import SurfaceArrhenius, StickingCoefficient # Separate because we cimport from rmgpy.kinetics.surface +from rmgpy.kinetics.surface import SurfaceArrhenius, StickingCoefficient, SurfaceChargeTransfer, SurfaceChargeTransferBEP # Separate because we cimport from rmgpy.kinetics.surface from rmgpy.kinetics.diffusionLimited import diffusion_limiter from rmgpy.molecule.element import Element, element_list from rmgpy.molecule.molecule import Molecule, Atom from rmgpy.pdep.reaction import calculate_microcanonical_rate_coefficient from rmgpy.species import Species +from rmgpy.thermo import ThermoData ################################################################################ @@ -92,7 +93,7 @@ class Reaction: `is_forward` ``bool`` Indicates if the reaction was generated in the forward (true) or reverse (false) `rank` ``int`` Integer indicating the accuracy of the kinetics for this reaction =================== =========================== ============================ - + """ def __init__(self, @@ -110,10 +111,11 @@ def __init__(self, pairs=None, allow_pdep_route=False, elementary_high_p=False, - allow_max_rate_violation=False, rank=None, + electrons=0, comment='', is_forward=None, + allow_max_rate_violation=False, ): self.index = index self.label = label @@ -129,11 +131,12 @@ def __init__(self, self.pairs = pairs self.allow_pdep_route = allow_pdep_route self.elementary_high_p = elementary_high_p + self.rank = rank + self.electrons = electrons self.comment = comment self.k_effective_cache = {} self.is_forward = is_forward self.allow_max_rate_violation = allow_max_rate_violation - self.rank = rank def __repr__(self): """ @@ -157,7 +160,8 @@ def __repr__(self): if self.elementary_high_p: string += 'elementary_high_p={0}, '.format(self.elementary_high_p) if self.comment != '': string += 'comment={0!r}, '.format(self.comment) if self.rank is not None: string += 'rank={0!r},'.format(self.rank) - string = string[:-2] + ')' + if self.electrons != 0: string += 'electrons={0:d},'.format(self.electrons) + string = string[:-1] + ')' return string def __str__(self): @@ -198,7 +202,8 @@ def __reduce__(self): self.allow_pdep_route, self.elementary_high_p, self.rank, - self.comment + self.electrons, + self.comment, )) @property @@ -231,6 +236,24 @@ def degeneracy(self, new): # set new degeneracy self._degeneracy = new + @property + def protons(self): + """ + The stochiometric coeff for protons in charge transfer reactions + """ + if self.is_charge_transfer_reaction(): + self._protons = 0 + for prod in self.products: + if prod.is_proton(): + self._protons += 1 + for react in self.reactants: + if react.is_proton(): + self._protons -= 1 + else: + self._protons = 0 + + return self._protons + def to_chemkin(self, species_list=None, kinetics=True): """ Return the chemkin-formatted string for this reaction. @@ -397,6 +420,22 @@ def is_surface_reaction(self): return True return False + def is_charge_transfer_reaction(self): + """ + Return ``True`` if one or more reactants or products are electrons + """ + if self.electrons != 0: + return True + return False + + def is_surface_charge_transfer_reaction(self): + """ + Return ``True`` if one or more reactants or products are electrons + """ + if self.is_surface_reaction() and self.is_charge_transfer_reaction(): + return True + return False + def has_template(self, reactants, products): """ Return ``True`` if the reaction matches the template of `reactants` @@ -464,6 +503,10 @@ def is_isomorphic(self, other, either_direction=True, check_identical=False, che strict=strict, save_order=save_order) + # compare stoichiometry of electrons in reaction + if self.electrons != other.electrons: + return False + # Compare reactants to reactants forward_reactants_match = same_species_lists(self.reactants, other.reactants, check_identical=check_identical, @@ -507,6 +550,33 @@ def is_isomorphic(self, other, either_direction=True, check_identical=False, che # should have already returned if it matches forwards, or we're not allowed to match backwards return reverse_reactants_match and reverse_products_match and collider_match + + def _apply_CHE_model(self, T): + """ + Apply the computational hydrogen electrode (CHE) model at temperature T (in 'K'). + + Returns the free energy (in J/mol) of 'N' proton/electron couple(s) in the reaction + using the Reversible Hydrogen Electrode (RHE) as referernce so that + N * deltaG(H+ + e-) = N * 1/2 deltaG(H2(g)) at 0V. + """ + + if not self.is_charge_transfer_reaction(): + raise ReactionError("CHE model is only applicable to charge transfer reactions!") + + if self.electrons != self.protons: + raise ReactionError("Number of electrons must equal number of protons! " + f"{self} has {self.electrons} protons and {self.electrons} electrons") + + + H2_thermo = ThermoData(Tdata=([300,400,500,600,800,1000,1500],'K'), + Cpdata=([6.895,6.975,6.994,7.009,7.081,7.219,7.72],'cal/(mol*K)'), + H298=(0,'kcal/mol'), S298=(31.233,'cal/(mol*K)','+|-',0.0007), + Cp0=(29.1007,'J/(mol*K)'), CpInf=(37.4151,'J/(mol*K)'), + label="""H2""", comment="""Thermo library: primaryThermoLibrary""") + # deltG_H+ + deltaG_e- -> 1/2 deltaG_H2 # only at 298K ??? + + return self.electrons * 0.5 * H2_thermo.get_free_energy(T) + def get_enthalpy_of_reaction(self, T): """ @@ -534,12 +604,40 @@ def get_entropy_of_reaction(self, T): dSrxn += product.get_entropy(T) return dSrxn - def get_free_energy_of_reaction(self, T): + def _get_free_energy_of_charge_transfer_reaction(self, T, potential=0.): + + cython.declare(dGrxn=cython.double, reactant=Species, product=Species) + + dGrxn = 0 + for reactant in self.reactants: + try: + dGrxn -= reactant.get_free_energy(T) + except Exception: + logging.error("Problem with reactant {!r} in reaction {!s}".format(reactant, self)) + raise + + for product in self.products: + try: + dGrxn += product.get_free_energy(T) + except Exception: + logging.error("Problem with product {!r} in reaction {!s}".format(reactant, self)) + raise + + if potential != 0.: + dGrxn -= self.electrons * constants.F * potential + + return dGrxn + + def get_free_energy_of_reaction(self, T, potential=0.): """ Return the Gibbs free energy of reaction in J/mol evaluated at - temperature `T` in K. + temperature `T` in K and potential in Volts (if applicable) """ cython.declare(dGrxn=cython.double, reactant=Species, product=Species) + + if self.is_charge_transfer_reaction(): + return self._get_free_energy_of_charge_transfer_reaction(T, potential=potential) + dGrxn = 0.0 for reactant in self.reactants: try: @@ -553,9 +651,33 @@ def get_free_energy_of_reaction(self, T): except Exception: logging.error("Problem with product {!r} in reaction {!s}".format(reactant, self)) raise + return dGrxn - def get_equilibrium_constant(self, T, type='Kc', surface_site_density=2.5e-05): + def get_reversible_potential(self, T): + """ + Get the Potential in `V` at T in 'K' at which the charge transfer reaction is at equilibrium + """ + cython.declare(deltaG=cython.double, V0=cython.double) + if not self.is_charge_transfer_reaction(): + raise KineticsError("Cannot get reversible potential for non charge transfer reactions") + + deltaG = self._get_free_energy_of_charge_transfer_reaction(T) #J/mol + V0 = deltaG / self.electrons / constants.F # V = deltaG / n / F + return V0 + + def set_reference_potential(self, T): + """ + Set the reference Potential of the `SurfaceChargeTransfer` kinetics model to the reversible potential + of the reaction + """ + if self.kinetics is None: + raise KineticsError("Cannot set reference potential for reactions with no kinetics attribute") + + if isinstance(self.kinetics, SurfaceChargeTransfer) and self.kinetics.V0 is None: + self.kinetics.V0 = (self.get_reversible_potential(T),'V') + + def get_equilibrium_constant(self, T, potential=0., type='Kc', surface_site_density=2.5e-05): """ Return the equilibrium constant for the reaction at the specified temperature `T` in K and reference `surface_site_density` @@ -564,11 +686,12 @@ def get_equilibrium_constant(self, T, type='Kc', surface_site_density=2.5e-05): ``Kc`` for concentrations (default), or ``Kp`` for pressures. This function assumes a reference pressure of 1e5 Pa for gas phases species and uses the ideal gas law to determine reference concentrations. For - surface species, the `surface_site_density` is the assumed reference. + surface species, the `surface_site_density` is the assumed reference. For protons (H+), + a reference concentration of 1000 mol/m^3 (1 mol/L) is assumed """ - cython.declare(dGrxn=cython.double, K=cython.double, C0=cython.double, P0=cython.double) + cython.declare(dN_gas=cython.int, dN_surf=cython.int, dGrxn=cython.double, K=cython.double, C0=cython.double, P0=cython.double) # Use free energy of reaction to calculate Ka - dGrxn = self.get_free_energy_of_reaction(T) + dGrxn = self.get_free_energy_of_reaction(T, potential) K = np.exp(-dGrxn / constants.R / T) # Convert Ka to Kc or Kp if specified # Assume a pressure of 1e5 Pa for gas phase species @@ -625,14 +748,14 @@ def get_entropies_of_reaction(self, Tlist): """ return np.array([self.get_entropy_of_reaction(T) for T in Tlist], np.float64) - def get_free_energies_of_reaction(self, Tlist): + def get_free_energies_of_reaction(self, Tlist, potential=0.): """ Return the Gibbs free energies of reaction in J/mol evaluated at temperatures `Tlist` in K. """ - return np.array([self.get_free_energy_of_reaction(T) for T in Tlist], np.float64) + return np.array([self.get_free_energy_of_reaction(T, potential=potential) for T in Tlist], np.float64) - def get_equilibrium_constants(self, Tlist, type='Kc'): + def get_equilibrium_constants(self, Tlist, potential=0., type='Kc'): """ Return the equilibrium constants for the reaction at the specified temperatures `Tlist` in K. The `type` parameter lets you specify the @@ -640,7 +763,7 @@ def get_equilibrium_constants(self, Tlist, type='Kc'): ``Kc`` for concentrations (default), or ``Kp`` for pressures. Note that this function currently assumes an ideal gas mixture. """ - return np.array([self.get_equilibrium_constant(T, type) for T in Tlist], np.float64) + return np.array([self.get_equilibrium_constant(T, potential=potential, type=type) for T in Tlist], np.float64) def get_stoichiometric_coefficient(self, spec): """ @@ -657,7 +780,7 @@ def get_stoichiometric_coefficient(self, spec): if product is spec: stoich += 1 return stoich - def get_rate_coefficient(self, T, P=0, surface_site_density=0): + def get_rate_coefficient(self, T, P=0, surface_site_density=0, potential=0.): """ Return the overall rate coefficient for the forward reaction at temperature `T` in K and pressure `P` in Pa, including any reaction @@ -670,7 +793,9 @@ def get_rate_coefficient(self, T, P=0, surface_site_density=0): If the reaction has sticking coefficient kinetics, a nonzero surface site density in `mol/m^2` must be provided """ - if isinstance(self.kinetics, StickingCoefficient): + if isinstance(self.kinetics,SurfaceChargeTransfer): + return self.get_surface_rate_coefficient(T, surface_site_density=surface_site_density, potential=potential) + elif isinstance(self.kinetics, StickingCoefficient): if surface_site_density <= 0: raise ValueError("Please provide a postive surface site density in mol/m^2 " f"for calculating the rate coefficient of {StickingCoefficient.__name__} kinetics") @@ -686,14 +811,15 @@ def get_rate_coefficient(self, T, P=0, surface_site_density=0): else: return self.kinetics.get_rate_coefficient(T, P) - def get_surface_rate_coefficient(self, T, surface_site_density): + def get_surface_rate_coefficient(self, T, surface_site_density, potential=0.): """ Return the overall surface rate coefficient for the forward reaction at temperature `T` in K with surface site density `surface_site_density` in mol/m2. Value is returned in combination of [m,mol,s] """ cython.declare(rateCoefficient=cython.double, - molecularWeight_kg=cython.double, ) + molecularWeight_kg=cython.double, + Ea=cython.double, deltaG=cython.double) if diffusion_limiter.enabled: raise NotImplementedError() @@ -727,6 +853,19 @@ def get_surface_rate_coefficient(self, T, surface_site_density): if isinstance(self.kinetics, SurfaceArrhenius): return self.kinetics.get_rate_coefficient(T, P=0) + if isinstance(self.kinetics, SurfaceChargeTransfer): + Ea = self.kinetics.get_activation_energy_from_potential(potential) + deltaG = self._get_free_energy_of_charge_transfer_reaction(298,potential) + if deltaG > 0 and Ea < deltaG: + corrected_kinetics = deepcopy(self.kinetics) + corrected_kinetics.V0.value_si = potential + corrected_kinetics.Ea.value_si = deltaG + logging.info("For reaction {0!s} Ea raised from {1:.1f} to {2:.1f} kJ/mol at {3:.2f} V".format( + self, self.kinetics.Ea.value_si / 1000., deltaG / 1000., potential)) + return corrected_kinetics.get_rate_coefficient(T, potential) + else: + return self.kinetics.get_rate_coefficient(T, potential) + raise NotImplementedError("Can't get_surface_rate_coefficient for kinetics type {!r}".format(type(self.kinetics))) def fix_diffusion_limited_a_factor(self, T): @@ -762,17 +901,17 @@ def fix_barrier_height(self, force_positive=False): If `force_positive` is True, then all reactions are forced to have a non-negative barrier. """ - cython.declare(H0=cython.double, H298=cython.double, Ea=cython.double) + cython.declare(H0=cython.double, H298=cython.double, Ea=cython.double, V0=cython.double, + deltaG=cython.double) if self.kinetics is None: raise KineticsError("Cannot fix barrier height for reactions with no kinetics attribute") - H298 = self.get_enthalpy_of_reaction(298) - H0 = sum([spec.get_thermo_data().E0.value_si for spec in self.products]) \ - - sum([spec.get_thermo_data().E0.value_si for spec in self.reactants]) - if isinstance(self.kinetics, (ArrheniusEP, SurfaceArrheniusBEP, StickingCoefficientBEP, ArrheniusBM)): + if isinstance(self.kinetics, SurfaceChargeTransferBEP): Ea = self.kinetics.E0.value_si # temporarily using Ea to store the intrinsic barrier height E0 - self.kinetics = self.kinetics.to_arrhenius(H298) + V0 = self.kinetics.V0.value_si + deltaG = self._get_free_energy_of_charge_transfer_reaction(298,V0) + self.kinetics = self.kinetics.to_surface_charge_transfer(deltaG) if self.kinetics.Ea.value_si < 0.0 and self.kinetics.Ea.value_si < Ea: # Calculated Ea (from Evans-Polanyi) is negative AND below than the intrinsic E0 Ea = min(0.0, Ea) # (the lowest we want it to be) @@ -781,33 +920,48 @@ def fix_barrier_height(self, force_positive=False): logging.info("For reaction {0!s} Ea raised from {1:.1f} to {2:.1f} kJ/mol.".format( self, self.kinetics.Ea.value_si / 1000., Ea / 1000.)) self.kinetics.Ea.value_si = Ea - if isinstance(self.kinetics, (Arrhenius, StickingCoefficient)): # SurfaceArrhenius is a subclass of Arrhenius - Ea = self.kinetics.Ea.value_si - if H0 >= 0 and Ea < H0: - self.kinetics.Ea.value_si = H0 - self.kinetics.comment += "\nEa raised from {0:.1f} to {1:.1f} kJ/mol to match endothermicity of " \ - "reaction.".format( Ea / 1000., H0 / 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., H0 / 1000., self)) - if force_positive and isinstance(self.kinetics, (Arrhenius, StickingCoefficient)) and self.kinetics.Ea.value_si < 0: - self.kinetics.comment += "\nEa raised from {0:.1f} to 0 kJ/mol.".format(self.kinetics.Ea.value_si / 1000.) - logging.info("For reaction {1!s} Ea raised from {0:.1f} to 0 kJ/mol.".format( - self.kinetics.Ea.value_si / 1000., self)) - self.kinetics.Ea.value_si = 0 - if self.kinetics.is_pressure_dependent() and self.network_kinetics is not None: - Ea = self.network_kinetics.Ea.value_si - if H0 >= 0 and Ea < H0: - self.network_kinetics.Ea.value_si = H0 - self.network_kinetics.comment += "\nEa raised from {0:.1f} to {1:.1f} kJ/mol to match endothermicity of" \ - " reaction.".format(Ea / 1000., H0 / 1000.) - logging.info("For reaction {2!s}, Ea of the high pressure limit kinetics raised from {0:.1f} to {1:.1f}" - " kJ/mol to match endothermicity of reaction.".format(Ea / 1000., H0 / 1000., self)) - if force_positive and isinstance(self.kinetics, Arrhenius) and self.kinetics.Ea.value_si < 0: - self.network_kinetics.comment += "\nEa raised from {0:.1f} to 0 kJ/mol.".format( - self.kinetics.Ea.value_si / 1000.) - logging.info("For reaction {1!s} Ea of the high pressure limit kinetics raised from {0:.1f} to 0" - " kJ/mol.".format(self.kinetics.Ea.value_si / 1000., self)) + else: + H298 = self.get_enthalpy_of_reaction(298) + H0 = sum([spec.get_thermo_data().E0.value_si for spec in self.products]) \ + - sum([spec.get_thermo_data().E0.value_si for spec in self.reactants]) + if isinstance(self.kinetics, (ArrheniusEP, SurfaceArrheniusBEP, StickingCoefficientBEP, ArrheniusBM)): + Ea = self.kinetics.E0.value_si # temporarily using Ea to store the intrinsic barrier height E0 + self.kinetics = self.kinetics.to_arrhenius(H298) + if self.kinetics.Ea.value_si < 0.0 and self.kinetics.Ea.value_si < Ea: + # Calculated Ea (from Evans-Polanyi) is negative AND below than the intrinsic E0 + Ea = min(0.0, Ea) # (the lowest we want it to be) + self.kinetics.comment += "\nEa raised from {0:.1f} to {1:.1f} kJ/mol.".format( + self.kinetics.Ea.value_si / 1000., Ea / 1000.) + logging.info("For reaction {0!s} Ea raised from {1:.1f} to {2:.1f} kJ/mol.".format( + self, self.kinetics.Ea.value_si / 1000., Ea / 1000.)) + self.kinetics.Ea.value_si = Ea + if isinstance(self.kinetics, (Arrhenius, StickingCoefficient)): # SurfaceArrhenius is a subclass of Arrhenius + Ea = self.kinetics.Ea.value_si + if H0 >= 0 and Ea < H0: + self.kinetics.Ea.value_si = H0 + self.kinetics.comment += "\nEa raised from {0:.1f} to {1:.1f} kJ/mol to match endothermicity of " \ + "reaction.".format( Ea / 1000., H0 / 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., H0 / 1000., self)) + if force_positive and isinstance(self.kinetics, (Arrhenius, StickingCoefficient)) and self.kinetics.Ea.value_si < 0: + self.kinetics.comment += "\nEa raised from {0:.1f} to 0 kJ/mol.".format(self.kinetics.Ea.value_si / 1000.) + logging.info("For reaction {1!s} Ea raised from {0:.1f} to 0 kJ/mol.".format( + self.kinetics.Ea.value_si / 1000., self)) self.kinetics.Ea.value_si = 0 + if self.kinetics.is_pressure_dependent() and self.network_kinetics is not None: + Ea = self.network_kinetics.Ea.value_si + if H0 >= 0 and Ea < H0: + self.network_kinetics.Ea.value_si = H0 + self.network_kinetics.comment += "\nEa raised from {0:.1f} to {1:.1f} kJ/mol to match endothermicity of" \ + " reaction.".format(Ea / 1000., H0 / 1000.) + logging.info("For reaction {2!s}, Ea of the high pressure limit kinetics raised from {0:.1f} to {1:.1f}" + " kJ/mol to match endothermicity of reaction.".format(Ea / 1000., H0 / 1000., self)) + if force_positive and isinstance(self.kinetics, Arrhenius) and self.kinetics.Ea.value_si < 0: + self.network_kinetics.comment += "\nEa raised from {0:.1f} to 0 kJ/mol.".format( + self.kinetics.Ea.value_si / 1000.) + logging.info("For reaction {1!s} Ea of the high pressure limit kinetics raised from {0:.1f} to 0" + " kJ/mol.".format(self.kinetics.Ea.value_si / 1000., self)) + self.kinetics.Ea.value_si = 0 def reverse_arrhenius_rate(self, k_forward, reverse_units, Tmin=None, Tmax=None): """ @@ -880,16 +1034,42 @@ def reverse_sticking_coeff_rate(self, k_forward, reverse_units, surface_site_den kr.fit_to_data(Tlist, klist, reverse_units, kf.T0.value_si) return kr + def reverse_surface_charge_transfer_rate(self, k_forward, reverse_units, Tmin=None, Tmax=None): + """ + Reverses the given k_forward, which must be a SurfaceChargeTransfer type. + You must supply the correct units for the reverse rate. + The equilibrium constant is evaluated from the current reaction instance (self). + """ + cython.declare(kf=SurfaceChargeTransfer, kr=SurfaceChargeTransfer) + cython.declare(Tlist=np.ndarray, klist=np.ndarray, i=cython.int, V0=cython.double) + kf = k_forward + self.set_reference_potential(298) + if not isinstance(kf, SurfaceChargeTransfer): # Only reverse SurfaceChargeTransfer rates + raise TypeError(f'Expected a SurfaceChargeTransfer object for k_forward but received {kf}') + if Tmin is not None and Tmax is not None: + Tlist = 1.0 / np.linspace(1.0 / Tmax.value, 1.0 / Tmin.value, 50) + else: + Tlist = np.linspace(298, 500, 30) + + V0 = self.kinetics.V0.value_si + klist = np.zeros_like(Tlist) + for i in range(len(Tlist)): + klist[i] = kf.get_rate_coefficient(Tlist[i],V0) / self.get_equilibrium_constant(Tlist[i],V0) + kr = SurfaceChargeTransfer(alpha=1-kf.alpha.value, electrons=-1*self.electrons, V0=(V0,'V')) + kr.fit_to_data(Tlist, klist, reverse_units, kf.T0.value_si) + return kr + def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, Tmax=None, surface_site_density=0): """ - Generate and return a rate coefficient model for the reverse reaction. + Generate and return a rate coefficient model for the reverse reaction. Currently this only works if the `kinetics` attribute is one of several (but not necessarily all) kinetics types. If the reaction kinetics model is Sticking Coefficient, please provide a nonzero surface site density in `mol/m^2` which is required to evaluate the rate coefficient. """ - cython.declare(Tlist=np.ndarray, Plist=np.ndarray, K=np.ndarray, + cython.declare(n_gas=cython.int, n_surf=cython.int, prod=Species, k_units=str, + Tlist=np.ndarray, Plist=np.ndarray, K=np.ndarray, rxn=Reaction, klist=np.ndarray, i=cython.size_t, Tindex=cython.size_t, Pindex=cython.size_t) @@ -897,6 +1077,7 @@ def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, T KineticsData.__name__, Arrhenius.__name__, SurfaceArrhenius.__name__, + SurfaceChargeTransfer.__name__, MultiArrhenius.__name__, PDepArrhenius.__name__, MultiPDepArrhenius.__name__, @@ -920,7 +1101,11 @@ def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, T kunits = get_rate_coefficient_units_from_reaction_order(n_gas, n_surf) kf = self.kinetics - if isinstance(kf, KineticsData): + + if isinstance(kf, SurfaceChargeTransfer): + return self.reverse_surface_charge_transfer_rate(kf, kunits, Tmin, Tmax) + + elif isinstance(kf, KineticsData): Tlist = kf.Tdata.value_si klist = np.zeros_like(Tlist) @@ -1090,10 +1275,13 @@ def is_balanced(self): Return ``True`` if the reaction has the same number of each atom on each side of the reaction equation, or ``False`` if not. """ - cython.declare(reactantElements=dict, productElements=dict, molecule=Molecule, atom=Atom, element=Element) + cython.declare(reactantElements=dict, productElements=dict, molecule=Molecule, atom=Atom, element=Element, + reactants_net_charge=cython.int, products_net_charge=cython.int) reactant_elements = {} product_elements = {} + reactants_net_charge = 0 + products_net_charge = 0 for element in element_list: reactant_elements[element] = 0 product_elements[element] = 0 @@ -1104,6 +1292,7 @@ def is_balanced(self): elif isinstance(reactant, Molecule): molecule = reactant for atom in molecule.atoms: + reactants_net_charge += atom.charge reactant_elements[atom.element] += 1 for product in self.products: @@ -1112,12 +1301,21 @@ def is_balanced(self): elif isinstance(product, Molecule): molecule = product for atom in molecule.atoms: + products_net_charge += atom.charge product_elements[atom.element] += 1 for element in element_list: if reactant_elements[element] != product_elements[element]: return False + if self.electrons < 0: + reactants_net_charge += self.electrons + elif self.electrons > 0: + products_net_charge -= self.electrons + + if reactants_net_charge != products_net_charge: + return False + return True def generate_pairs(self): diff --git a/rmgpy/reactionTest.py b/rmgpy/reactionTest.py index 64e8eb002a..47a9a63321 100644 --- a/rmgpy/reactionTest.py +++ b/rmgpy/reactionTest.py @@ -37,10 +37,11 @@ import numpy from external.wip import work_in_progress from copy import deepcopy +import math import rmgpy.constants as constants from rmgpy.kinetics import Arrhenius, ArrheniusEP, MultiArrhenius, PDepArrhenius, MultiPDepArrhenius, \ - ThirdBody, Troe, Lindemann, Chebyshev, SurfaceArrhenius, StickingCoefficient + ThirdBody, Troe, Lindemann, Chebyshev, SurfaceArrhenius, StickingCoefficient, SurfaceChargeTransfer from rmgpy.molecule import Molecule from rmgpy.quantity import Quantity from rmgpy.reaction import Reaction @@ -52,6 +53,10 @@ from rmgpy.statmech.vibration import HarmonicOscillator from rmgpy.thermo import Wilhoit, ThermoData, NASA, NASAPolynomial +################################################################################ + +def order_of_magnitude(number): + return math.floor(math.log(number, 10)) ################################################################################ @@ -241,6 +246,15 @@ def setUp(self): T0=(1.0, 'K'), coverage_dependence={Species().from_adjacency_list('1 X u0 p0 c0'): {'E': (0.1, 'J/mol'), 'm': -1.0, 'a': 1.0}, Species().from_adjacency_list('1 X u0 p0 c0 {2,S}\n2 H u0 p0 c0 {1,S}'): {'E': (0.2, 'J/mol'), 'm': -2.0, 'a': 2.0}})) + def test_electrons(self): + """Test electrons property""" + self.assertEquals(self.rxn1s.electrons,0) + self.assertEquals(self.rxn1s.electrons,0) + + def test_protons(self): + """Test protons property""" + self.assertEquals(self.rxn1s.protons,0) + self.assertEquals(self.rxn1s.protons,0) def test_is_surface_reaction_species(self): """Test is_surface_reaction for reaction based on Species """ @@ -313,6 +327,193 @@ def test_reverse_sticking_coeff_rate(self): krevrev = rxn_copy.get_rate_coefficient(T, P, surface_site_density=2.5e-5) self.assertAlmostEqual(korig / krevrev, 1.0, 0) +class TestChargeTransferReaction(unittest.TestCase): + """Test charge transfer reactions""" + + def setUp(self): + m_electron = Molecule().from_smiles("e") + m_proton = Molecule().from_smiles("[H+]") + m_x = Molecule().from_adjacency_list("1 X u0 p0") + m_ch2x = Molecule().from_adjacency_list( + """ + 1 C u0 p0 c0 {2,S} {3,S} {4,D} + 2 H u0 p0 c0 {1,S} + 3 H u0 p0 c0 {1,S} + 4 X u0 p0 c0 {1,D} + """ + ) + m_ch3x = Molecule().from_adjacency_list( + """ + 1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} + 2 H u0 p0 c0 {1,S} + 3 H u0 p0 c0 {1,S} + 4 H u0 p0 c0 {1,S} + 5 X u0 p0 c0 {1,S} + """ + ) + + s_electron = Species( + molecule=[m_electron], + thermo=ThermoData(Tdata=([300, 400, 500, 600, 800, 1000, 1500, 2000], + "K"), + Cpdata=([0., 0., 0., 0., 0., 0., 0., 0.], "cal/(mol*K)"), + H298=(0.0, "kcal/mol"), + S298=(0.0, "cal/(mol*K)"))) + + s_proton = Species( + molecule=[m_proton], + thermo = ThermoData( + Tdata=([300,400,500,600,800,1000,1500],'K'), + Cpdata=([3.4475,3.4875,3.497,3.5045,3.5405,3.6095,3.86],'cal/(mol*K)'), + H298=(0,'kcal/mol'), S298=(15.6165,'cal/(mol*K)','+|-',0.0007), + comment = '1/2 free energy of H2(g)')) + s_x = Species( + molecule=[m_x], + thermo=ThermoData(Tdata=([300, 400, 500, 600, 800, 1000, 1500, 2000], + "K"), + Cpdata=([0., 0., 0., 0., 0., 0., 0., 0.], "cal/(mol*K)"), + H298=(0.0, "kcal/mol"), + S298=(0.0, "cal/(mol*K)"))) + + s_ch2x = Species( + molecule=[m_ch2x], + thermo=ThermoData(Tdata=([300,400,500,600,800,1000,1500],'K'), + Cpdata=([28.4959,36.3588,42.0219,46.2428,52.3978,56.921,64.1119],'J/(mol*K)'), + H298=(0.654731,'kJ/mol'), S298=(19.8341,'J/(mol*K)'), Cp0=(0.01,'J/(mol*K)'), + CpInf=(99.7737,'J/(mol*K)'), + comment="""Thermo library: surfaceThermoPt111 Binding energy corrected by LSR (0.50C)""")) + + s_ch3x = Species( + molecule=[m_ch3x], + thermo=ThermoData(Tdata=([300,400,500,600,800,1000,1500],'K'), + Cpdata=([37.3325,44.9406,51.3613,56.8115,65.537,72.3287,83.3007],'J/(mol*K)'), + H298=(-45.8036,'kJ/mol'), S298=(57.7449,'J/(mol*K)'), Cp0=(0.01,'J/(mol*K)'), + CpInf=(124.717,'J/(mol*K)'), + comment="""Thermo library: surfaceThermoPt111 Binding energy corrected by LSR (0.25C)""")) + + # X=CH2 + H+ + e- <--> X-CH3 + rxn_reduction = Reaction(reactants=[s_proton, s_ch2x], + products=[s_ch3x], + electrons = -1, + kinetics=SurfaceChargeTransfer( + A = (2.483E21, 'cm^3/(mol*s)'), + V0 = (0, 'V'), + Ea = (10, 'kJ/mol'), + Tmin = (200, 'K'), + Tmax = (3000, 'K'), + electrons = -1, + )) + + rxn_oxidation = Reaction(products=[s_proton, s_ch2x], + reactants=[s_ch3x], + electrons = 1, + kinetics=SurfaceChargeTransfer( + A = (2.483E21, 'cm^3/(mol*s)'), + V0 = (0, 'V'), + Ea = (10, 'kJ/mol'), + Tmin = (200, 'K'), + Tmax = (3000, 'K'), + electrons = 1, + )) + + self.rxn_reduction = rxn_reduction + self.rxn_oxidation = rxn_oxidation + + def test_electrons(self): + """Test electrons property""" + self.assertEquals(self.rxn_reduction.electrons,-1) + self.assertEquals(self.rxn_oxidation.electrons,1) + + def test_protons(self): + """Test n_protons property""" + self.assertEquals(self.rxn_reduction.protons,-1) + self.assertEquals(self.rxn_oxidation.protons,1) + + def test_is_surface_reaction(self): + """Test is_surface_reaction() method""" + self.assertTrue(self.rxn_reduction.is_surface_reaction()) + self.assertTrue(self.rxn_oxidation.is_surface_reaction()) + + def test_is_charge_transfer_reaction(self): + """Test is_charge_transfer_reaction() method""" + self.assertTrue(self.rxn_reduction.is_charge_transfer_reaction()) + self.assertTrue(self.rxn_oxidation.is_charge_transfer_reaction()) + + def test_is_surface_charge_transfer(self): + """Test is_surface_charge_transfer() method""" + self.assertTrue(self.rxn_reduction.is_surface_charge_transfer_reaction()) + self.assertTrue(self.rxn_oxidation.is_surface_charge_transfer_reaction()) + + def test_get_reversible_potential(self): + """Test get_reversible_potential() method""" + V0_reduction = self.rxn_reduction.get_reversible_potential(298) + V0_oxidation = self.rxn_oxidation.get_reversible_potential(298) + + self.assertAlmostEqual(V0_reduction,V0_oxidation,6) + self.assertAlmostEqual(V0_reduction, 0.3967918, 6) + self.assertAlmostEqual(V0_oxidation, 0.3967918, 6) + + def test_get_rate_coeff(self): + """Test get_rate_coefficient() method""" + + # these should be the same + kf_1 = self.rxn_reduction.get_rate_coefficient(298,potential=0) + kf_2 = self.rxn_reduction.kinetics.get_rate_coefficient(298,0) + + self.assertAlmostEqual(kf_1, 43870506959779.0, 6) + self.assertAlmostEqual(kf_1, kf_2, 6) + + #kf_2 should be greater than kf_1 + kf_1 = self.rxn_oxidation.get_rate_coefficient(298,potential=0) + kf_2 = self.rxn_oxidation.get_rate_coefficient(298,potential=0.1) + self.assertGreater(kf_2,kf_1) + + def test_equilibrium_constant_surface_charge_transfer_kc(self): + """ + Test the equilibrium constants of type Kc of a surface charge transfer reaction. + """ + Tlist = numpy.arange(400.0, 1600.0, 200.0, numpy.float64) + Kclist0 = [1.39365463e+03, 1.78420988e+01, 2.10543835e+00, 6.07529099e-01, + 2.74458007e-01, 1.59985450e-01] #reduction + Kclist_reduction = self.rxn_reduction.get_equilibrium_constants(Tlist, type='Kc') + Kclist_oxidation = self.rxn_oxidation.get_equilibrium_constants(Tlist, type='Kc') + # Test a range of temperatures + for i in range(len(Tlist)): + self.assertAlmostEqual(Kclist_reduction[i] / Kclist0[i], 1.0, 6) + self.assertAlmostEqual(1 / Kclist_oxidation[i] / Kclist0[i], 1.0, 6) + + V0 = self.rxn_oxidation.get_reversible_potential(298) + Kc_reduction_equil = self.rxn_reduction.get_equilibrium_constant(298, V0) + Kc_oxidation_equil = self.rxn_oxidation.get_equilibrium_constant(298, V0) + self.assertAlmostEqual(Kc_oxidation_equil * Kc_reduction_equil, 1.0, 4) + C0 = 1e5 / constants.R / 298 + self.assertAlmostEqual(Kc_oxidation_equil, C0, 4) + self.assertAlmostEqual(Kc_reduction_equil, 1/C0, 4) + + @work_in_progress + def test_reverse_surface_charge_transfer_rate(self): + """ + Test the reverse_surface_charge_transfer_rate() method + """ + kf_reduction = self.rxn_reduction.kinetics + kf_oxidation = self.rxn_oxidation.kinetics + kr_oxidation = self.rxn_oxidation.generate_reverse_rate_coefficient() + kr_reduction = self.rxn_reduction.generate_reverse_rate_coefficient() + self.assertEquals(kr_reduction.A.units, 's^-1') + self.assertEquals(kr_oxidation.A.units, 'm^3/(mol*s)') + + Tlist = numpy.linspace(298., 500., 30.,numpy.float64) + for T in Tlist: + for V in (-0.25,0.,0.25): + kf = kf_reduction.get_rate_coefficient(T,V) + kr = kr_reduction.get_rate_coefficient(T,V) + K = self.rxn_reduction.get_equilibrium_constant(T,V) + self.assertEquals(order_of_magnitude(kf/kr),order_of_magnitude(K)) + kf = kf_oxidation.get_rate_coefficient(T,V) + kr = kr_oxidation.get_rate_coefficient(T,V) + K = self.rxn_oxidation.get_equilibrium_constant(T,V) + self.assertEquals(order_of_magnitude(kf/kr),order_of_magnitude(K)) + class TestReaction(unittest.TestCase): """ Contains unit tests of the Reaction class. @@ -1323,6 +1524,7 @@ def test_output(self): from its repr() output with no loss of information. """ namespace = {} + print('reaction = {0!r}'.format(self.reaction)) exec('reaction = {0!r}'.format(self.reaction), globals(), namespace) self.assertIn('reaction', namespace) reaction = namespace['reaction'] diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 23a58c7546..982fc6da15 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -44,6 +44,7 @@ from rmgpy.solver.mbSampled import MBSampledReactor from rmgpy.solver.simple import SimpleReactor from rmgpy.solver.surface import SurfaceReactor +from rmgpy.solver.electrode import ElectrodeReactor from rmgpy.util import as_list from rmgpy.data.surface import MetalDatabase @@ -521,6 +522,109 @@ def surface_reactor(temperature, rmg.reaction_systems.append(system) system.log_initial_conditions(number=len(rmg.reaction_systems)) +# Reaction systems +def electrode_reactor(potential, + temperature, + initialPressure, + initialGasMoleFractions, + initialSurfaceCoverages, + surfaceVolumeRatio, + nSims=4, + terminationConversion=None, + terminationTime=None, + terminationRateRatio=None, + sensitivity=None, + sensitivityThreshold=1e-3): + logging.debug('Found electrodeReactor reaction system') + + for value in list(initialGasMoleFractions.values()): + if value < 0: + raise InputError('Initial mole fractions cannot be negative.') + total_initial_moles = sum(initialGasMoleFractions.values()) + if total_initial_moles != 1: + logging.warning('Initial gas mole fractions do not sum to one; renormalizing.') + logging.debug('') + logging.debug('Original composition:') + for spec, molfrac in initialGasMoleFractions.items(): + logging.debug("{0} = {1}".format(spec, molfrac)) + for spec in initialGasMoleFractions: + initialGasMoleFractions[spec] /= total_initial_moles + logging.info('') + logging.debug('Normalized mole fractions:') + for spec, molfrac in initialGasMoleFractions.items(): + logging.debug("{0} = {1}".format(spec, molfrac)) + + if not isinstance(potential, list): + potential = Quantity(potential) + else: + if len(potential) != 2: + raise InputError('potential ranges can either be in the form of (number,units) or a list with 2 entries ' + 'of the same format') + potential = [Quantity(v) for v in potential] + + if not isinstance(temperature, list): + T = Quantity(temperature) + else: + if len(temperature) != 2: + raise InputError('Temperature ranges can either be in the form of (number,units) or a list with 2 entries ' + 'of the same format') + T = [Quantity(t) for t in temperature] + + if not isinstance(initialPressure, list): + P_initial = Quantity(initialPressure) + else: + if len(initialPressure) != 2: + raise InputError('Initial pressure ranges can either be in the form ''of (number,units) or a list with ' + '2 entries of the same format') + P_initial = [Quantity(p) for p in initialPressure] + + if not isinstance(temperature, list) and not isinstance(initialPressure, list): + nSims = 1 + if any([isinstance(x, list) for x in initialGasMoleFractions.values()]) or \ + any([isinstance(x, list) for x in initialSurfaceCoverages.values()]): + raise NotImplementedError("Can't do ranges on species concentrations for surface reactors yet.") + + termination = [] + if terminationConversion is not None: + for spec, conv in terminationConversion.items(): + termination.append(TerminationConversion(species_dict[spec], conv)) + if terminationTime is not None: + termination.append(TerminationTime(Quantity(terminationTime))) + if terminationRateRatio is not None: + termination.append(TerminationRateRatio(terminationRateRatio)) + if len(termination) == 0: + raise InputError('No termination conditions specified for reaction system #{0}.'.format(len(rmg.reaction_systems) + 2)) + + sensitive_species = [] + if sensitivity: + for spec in sensitivity: + sensitive_species.append(species_dict[spec]) + if not isinstance(T, list): + sensitivityTemperature = T + if not isinstance(initialPressure, list): + sensitivityPressure = initialPressure + sens_conditions = None + if sensitivity: + raise NotImplementedError("Can't currently do sensitivity with surface reactors.") + # The problem is inside base.pyx it reads the dictionary 'sensConditions' + # and guesses whether they're all concentrations (liquid reactor) or + # mole fractions (simple reactor). In fact, some may be surface coverages. + + system = ElectrodeReactor(potential=potential, + T=T, + P_initial=P_initial, + initial_gas_mole_fractions=initialGasMoleFractions, + initial_surface_coverages=initialSurfaceCoverages, + surface_volume_ratio=surfaceVolumeRatio, + surface_site_density=rmg.surface_site_density, + n_sims=nSims, + termination=termination, + sensitive_species=sensitive_species, + sensitivity_threshold=sensitivityThreshold, + sens_conditions=sens_conditions) + rmg.reaction_systems.append(system) + system.log_initial_conditions(number=len(rmg.reaction_systems)) + # Reaction systems def mb_sampled_reactor(temperature, @@ -967,6 +1071,7 @@ def read_input_file(path, rmg0): 'simpleReactor': simple_reactor, 'liquidReactor': liquid_reactor, 'surfaceReactor': surface_reactor, + 'electrodeReactor': electrode_reactor, 'mbsampledReactor': mb_sampled_reactor, 'simulator': simulator, 'solvation': solvation, diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index deaf99304a..ed9d663304 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -73,6 +73,7 @@ from rmgpy.rmg.settings import ModelSettings from rmgpy.solver.base import TerminationTime, TerminationConversion from rmgpy.solver.simple import SimpleReactor +from rmgpy.solver.electrode import ElectrodeReactor from rmgpy.stats import ExecutionStatsWriter from rmgpy.thermo.thermoengine import submit from rmgpy.tools.plot import plot_sensitivity @@ -166,6 +167,8 @@ def __init__(self, input_file=None, output_directory=None, profiler=None, stats_ self.Tmax = 0.0 self.Pmin = 0.0 self.Pmax = 0.0 + self.potential_min = 0.0 + self.potential_max = 0.0 self.database = None def clear(self): @@ -685,6 +688,8 @@ def execute(self, initialize=True, **kwargs): try: self.Pmin = min([x.Prange[0].value_si if hasattr(x, 'Prange') and x.Prange else x.P.value_si for x in self.reaction_systems]) self.Pmax = max([x.Prange[1].value_si if hasattr(x, 'Prange') and x.Prange else x.P.value_si for x in self.reaction_systems]) + self.potential_min = min([x.potential_range[0].value_si if x.potential_range else x.potential.value_si for x in self.reaction_systems]) + self.potential_max = max([x.potential_range[1].value_si if x.potential_range else x.potential.value_si for x in self.reaction_systems]) except AttributeError: pass @@ -698,6 +703,16 @@ def execute(self, initialize=True, **kwargs): self.rmg_memories.append(RMG_Memory(reaction_system, self.balance_species)) self.rmg_memories[index].generate_cond() log_conditions(self.rmg_memories, index) + conditions = self.rmg_memories[index].get_cond() + temperature = potential = None + if isinstance(reaction_system, ElectrodeReactor): + if conditions: + potential = conditions.get('potential') + temperature = conditions.get('T') + if not potential: + potential = reaction_system.potential.value_si + if not temperature: + temperature = reaction_system.T.value_si # Update react flags if self.filter_reactions: @@ -711,7 +726,7 @@ def execute(self, initialize=True, **kwargs): atol=self.simulator_settings_list[0].atol, rtol=self.simulator_settings_list[0].rtol, filter_reactions=True, - conditions=self.rmg_memories[index].get_cond(), + conditions=conditions, ) self.update_reaction_threshold_and_react_flags( @@ -732,7 +747,9 @@ def execute(self, initialize=True, **kwargs): self.reaction_model.enlarge(react_edge=True, unimolecular_react=self.unimolecular_react, bimolecular_react=self.bimolecular_react, - trimolecular_react=self.trimolecular_react) + trimolecular_react=self.trimolecular_react, + temperature=temperature, + potential=potential) if not np.isinf(self.model_settings_list[0].thermo_tol_keep_spc_in_edge): self.reaction_model.set_thermodynamic_filtering_parameters( @@ -793,6 +810,15 @@ def execute(self, initialize=True, **kwargs): objects_to_enlarge = [] conditions = self.rmg_memories[index].get_cond() + if isinstance(reaction_system, ElectrodeReactor): + if conditions: + potential = conditions.get('potential') + temperature = conditions.get('T') + if not potential: + potential = reaction_system.potential.value_si + if not temperature: + temperature = reaction_system.T.value_si + if conditions and self.solvent: T = conditions['T'] # Set solvent viscosity @@ -822,7 +848,7 @@ def execute(self, initialize=True, **kwargs): prune=prune, model_settings=model_settings, simulator_settings=simulator_settings, - conditions=self.rmg_memories[index].get_cond() + conditions=conditions ) except: logging.error("Model core reactions:") @@ -835,9 +861,9 @@ def execute(self, initialize=True, **kwargs): self.make_seed_mech() # Just in case the user wants to restart from this raise + log_conditions(self.rmg_memories, index) self.rmg_memories[index].add_t_conv_N(t, x, len(obj)) self.rmg_memories[index].generate_cond() - log_conditions(self.rmg_memories, index) reactor_done = self.reaction_model.add_new_surface_objects(obj, new_surface_species, new_surface_reactions, reaction_system) @@ -922,7 +948,9 @@ def execute(self, initialize=True, **kwargs): self.reaction_model.enlarge(react_edge=True, unimolecular_react=self.unimolecular_react, bimolecular_react=self.bimolecular_react, - trimolecular_react=self.trimolecular_react) + trimolecular_react=self.trimolecular_react, + potential=potential, + temperature=temperature) if old_edge_size != len(self.reaction_model.edge.reactions) or old_core_size != len( self.reaction_model.core.reactions): @@ -2212,6 +2240,9 @@ def __init__(self, reaction_system, bspc): if hasattr(reaction_system, 'Prange') and isinstance(reaction_system.Prange, list): Prange = reaction_system.Prange self.Ranges['P'] = [np.log(P.value_si) for P in Prange] + if hasattr(reaction_system, 'potential_range') and isinstance(reaction_system.potential_range, list): + potential_range = reaction_system.potential_range + self.Ranges['potential'] = [V.value_si for V in potential_range] if hasattr(reaction_system, 'initial_mole_fractions'): if bspc: self.initial_mole_fractions = deepcopy(reaction_system.initial_mole_fractions) @@ -2326,6 +2357,8 @@ def obj(y): enumerate(ykey)} if 'P' in list(new_cond.keys()): new_cond['P'] = np.exp(new_cond['P']) + if 'potential' in list(new_cond.keys()): + new_cond['potential'] = new_cond['potential'] if hasattr(self, 'initial_mole_fractions'): for key in self.initial_mole_fractions.keys(): @@ -2355,6 +2388,8 @@ def log_conditions(rmg_memories, index): s += 'T = {0} K, '.format(item) elif key == 'P': s += 'P = {0} bar, '.format(item / 1.0e5) + elif key == 'potential': + s += 'potential = {0} V, '.format(item) else: s += key.label + ' = {0}, '.format(item) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 840d7d1550..10e6381185 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -48,7 +48,7 @@ from rmgpy.data.rmg import get_db from rmgpy.display import display from rmgpy.exceptions import ForbiddenStructureException -from rmgpy.kinetics import KineticsData, Arrhenius +from rmgpy.kinetics import KineticsData, Arrhenius, SurfaceChargeTransfer from rmgpy.quantity import Quantity from rmgpy.reaction import Reaction from rmgpy.rmg.pdep import PDepReaction, PDepNetwork @@ -528,7 +528,8 @@ def make_new_pdep_reaction(self, forward): return forward def enlarge(self, new_object=None, react_edge=False, - unimolecular_react=None, bimolecular_react=None, trimolecular_react=None): + unimolecular_react=None, bimolecular_react=None, trimolecular_react=None, + temperature=None, potential=None): """ Enlarge a reaction model by processing the objects in the list `new_object`. If `new_object` is a @@ -660,9 +661,13 @@ def enlarge(self, new_object=None, react_edge=False, # self.new_reaction_list only contains *actually* new reactions, all in the forward direction. for reaction in self.new_reaction_list: # convert KineticsData to Arrhenius forms - if isinstance(reaction.kinetics, KineticsData): + if isinstance(reaction.kinetics, KineticsData) and not isinstance(reaction.kinetics, SurfaceChargeTransfer): reaction.kinetics = reaction.kinetics.to_arrhenius() - # correct barrier heights of estimated kinetics + + if isinstance(reaction.kinetics, SurfaceChargeTransfer): + reaction.set_reference_potential(temperature) + + # correct barrier heights of estimated kinetics if isinstance(reaction, TemplateReaction) or isinstance(reaction, DepositoryReaction): # i.e. not LibraryReaction reaction.fix_barrier_height() # also converts ArrheniusEP to Arrhenius. diff --git a/rmgpy/solver/__init__.py b/rmgpy/solver/__init__.py index b109055e7f..4d3a42e876 100644 --- a/rmgpy/solver/__init__.py +++ b/rmgpy/solver/__init__.py @@ -31,4 +31,5 @@ from rmgpy.solver.simple import SimpleReactor from rmgpy.solver.liquid import LiquidReactor from rmgpy.solver.surface import SurfaceReactor +from rmgpy.solver.electrode import ElectrodeReactor from rmgpy.solver.mbSampled import MBSampledReactor diff --git a/rmgpy/solver/base.pyx b/rmgpy/solver/base.pyx index 76277c69ef..430f1a78e3 100644 --- a/rmgpy/solver/base.pyx +++ b/rmgpy/solver/base.pyx @@ -213,6 +213,8 @@ cdef class ReactionSystem(DASx): self.T = Quantity(conditions['T'], 'K') if 'P' in keys and hasattr(self, 'P'): self.P = Quantity(conditions['P'], 'Pa') + if 'potential' in keys and hasattr(self, 'potential'): + self.potential = Quantity(conditions['potential'], 'V') for k in keys: if is_conc: if k in self.initial_concentrations.keys(): diff --git a/rmgpy/solver/electrode.pyx b/rmgpy/solver/electrode.pyx new file mode 100644 index 0000000000..6d32f4d133 --- /dev/null +++ b/rmgpy/solver/electrode.pyx @@ -0,0 +1,547 @@ +############################################################################### +# # +# RMG - Reaction Mechanism Generator # +# # +# Copyright (c) 2002-2020 Prof. William H. Green (whgreen@mit.edu), # +# Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu) # +# # +# Permission is hereby granted, free of charge, to any person obtaining a # +# copy of this software and associated documentation files (the 'Software'), # +# to deal in the Software without restriction, including without limitation # +# the rights to use, copy, modify, merge, publish, distribute, sublicense, # +# and/or sell copies of the Software, and to permit persons to whom the # +# Software is furnished to do so, subject to the following conditions: # +# # +# The above copyright notice and this permission notice shall be included in # +# all copies or substantial portions of the Software. # +# # +# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING # +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER # +# DEALINGS IN THE SOFTWARE. # +# # +############################################################################### + +""" +Contains the :class:`ElectrodeReactor` class, providing a reaction system +consisting of a homogeneous, isothermal, isobaric batch reactor. +""" + +import itertools +import logging + +cimport cython +import numpy as np +cimport numpy as np + +import rmgpy.constants as constants +cimport rmgpy.constants as constants +from rmgpy.quantity import Quantity +from rmgpy.quantity cimport ScalarQuantity +from rmgpy.solver.base cimport ReactionSystem + +cdef class ElectrodeReactor(ReactionSystem): + """ + A reaction system consisting of a heterogeneous, isothermal, constant volume batch + reactor with an applied potential + """ + + cdef public ScalarQuantity T + cdef public ScalarQuantity P_initial + cdef public ScalarQuantity potential + cdef public double V + cdef public ScalarQuantity I + cdef public bint constant_volume + + cdef public list Trange + cdef public list Prange + cdef public list potential_range + cdef public int n_sims + cdef public dict sens_conditions + + cdef public dict initial_gas_mole_fractions + cdef public dict initial_surface_coverages + cdef public ScalarQuantity surface_volume_ratio + cdef public ScalarQuantity surface_site_density + cdef public np.ndarray reactions_on_electrode # (electrode surface, not core/edge surface) + cdef public np.ndarray species_on_electrode # (electrode surface, not core/edge surface) + cdef public np.ndarray charge_transfer_reactions # (electrode surface, not core/edge surface) + + + def __init__(self, + potential, + T, + P_initial, + initial_gas_mole_fractions, + initial_surface_coverages, + surface_volume_ratio, + surface_site_density, + n_sims=None, + termination=None, + sensitive_species=None, + sensitivity_threshold=1e-3, + sens_conditions=None, + ): + ReactionSystem.__init__(self, + termination, + sensitive_species, + sensitivity_threshold) + + if isinstance(potential, list): + self.potential_range = [Quantity(v) for v in potential] + else: + self.potential = Quantity(potential) + + if isinstance(T, list): + self.Trange = [Quantity(t) for t in T] + else: + self.T = Quantity(T) + if isinstance(P_initial, list): + raise NotImplementedError("Can't do ranges of initial pressures for surface reactors yet") + else: + self.P_initial = Quantity(P_initial) + self.initial_gas_mole_fractions = initial_gas_mole_fractions + self.initial_surface_coverages = initial_surface_coverages + self.surface_volume_ratio = Quantity(surface_volume_ratio) + self.surface_site_density = Quantity(surface_site_density) + self.V = 0 # will be set from ideal gas law in initialize_model + self.constant_volume = True + self.sens_conditions = sens_conditions + self.n_sims = n_sims + + def convert_initial_keys_to_species_objects(self, species_dict): + """ + Convert the initial_gas_mole_fractions and initial_surface_coverages dictionaries + from species names into species objects, + using the given dictionary of species. + """ + initial_gas_mole_fractions = {} + for label, moleFrac in self.initial_gas_mole_fractions.items(): + initial_gas_mole_fractions[species_dict[label]] = moleFrac + self.initial_gas_mole_fractions = initial_gas_mole_fractions + initial_surface_coverages = {} + for label, surfaceCoverage in self.initial_surface_coverages.items(): + initial_surface_coverages[species_dict[label]] = surfaceCoverage + self.initial_surface_coverages = initial_surface_coverages + + def generate_reactant_product_indices(self, core_reactions, edge_reactions): + """ + Creates a matrix for the reactants and products. + """ + + self.reactant_indices = -np.ones((self.num_core_reactions + self.num_edge_reactions, 3), np.int) + self.product_indices = -np.ones_like(self.reactant_indices) + + for rxn in itertools.chain(core_reactions, edge_reactions): + j = self.reaction_index[rxn] + for l, spec in enumerate(rxn.reactants): + i = self.get_species_index(spec) + self.reactant_indices[j, l] = i + for l, spec in enumerate(rxn.products): + i = self.get_species_index(spec) + self.product_indices[j, l] = i + + cpdef initialize_model(self, + list core_species, + list core_reactions, + list edge_species, + list edge_reactions, + list surface_species=[], + list surface_reactions=[], + list pdep_networks=None, + atol=1e-16, + rtol=1e-8, + sensitivity=False, + sens_atol=1e-6, + sens_rtol=1e-4, + filter_reactions=False, + dict conditions=None, + ): + """ + Initialize a simulation of the simple reactor using the provided kinetic + model. + """ + + # First call the base class version of the method + # This initializes the attributes declared in the base class + ReactionSystem.initialize_model(self, + core_species=core_species, + core_reactions=core_reactions, + edge_species=edge_species, + edge_reactions=edge_reactions, + surface_species=surface_species, + surface_reactions=surface_reactions, + pdep_networks=pdep_networks, + atol=atol, + rtol=rtol, + sensitivity=sensitivity, + sens_atol=sens_atol, + sens_rtol=sens_rtol, + filter_reactions=filter_reactions, + conditions=conditions, + ) + cdef np.ndarray[np.int_t, ndim=1] species_on_electrode, reactions_on_electrode + cdef int index + #: 1 if it's on a surface, 0 if it's in the gas phase + reactions_on_electrode = np.zeros((self.num_core_reactions + self.num_edge_reactions), np.int) + charge_transfer_reactions = np.zeros((self.num_core_reactions + self.num_edge_reactions), np.int) + species_on_electrode = np.zeros((self.num_core_species), np.int) + for spec, index in self.species_index.items(): + if index >= self.num_core_species: + continue + if spec.contains_surface_site(): + species_on_electrode[index] = 1 + for rxn, index in self.reaction_index.items(): + if rxn.is_surface_reaction(): + reactions_on_electrode[index] = 1 + if rxn.is_charge_transfer_reaction(): + charge_transfer_reactions[index] = rxn.electrons + + self.species_on_electrode = species_on_electrode + self.reactions_on_electrode = reactions_on_electrode + self.charge_transfer_reactions = charge_transfer_reactions + + # Set initial conditions + self.set_initial_conditions() + + # Compute reaction thresholds if reaction filtering is turned on + if filter_reactions: + ReactionSystem.set_initial_reaction_thresholds(self) + + # Generate forward and reverse rate coefficients k(T,P) + self.generate_rate_coefficients(core_reactions, edge_reactions) + + ReactionSystem.compute_network_variables(self, pdep_networks) + + ReactionSystem.set_initial_derivative(self) + + # Initialize the model + ReactionSystem.initialize_solver(self) + + def generate_rate_coefficients(self, core_reactions, edge_reactions): + """ + Populates the kf, kb and equilibriumConstants + arrays with the values computed at the temperature, (effective) pressure, and potential of the + reaction system. + """ + + cdef double P, surface_volume_ratio_si, kb + + surface_volume_ratio_si = self.surface_volume_ratio.value_si + # ToDo: Pressure should come from ideal gas law? + P = self.P_initial.value_si + + + warned = False + for rxn in itertools.chain(core_reactions, edge_reactions): + j = self.reaction_index[rxn] + + # ToDo: get_rate_coefficient should also depend on surface coverages vector + + if rxn.is_surface_reaction(): + """ + Be careful! From here on kf and kb will now be in Volume units, + even for surface reactions (which you may expect to be in Area units). + This is to avoid repeatedly multiplying a bunch of things inside every + loop of the ODE solver. + """ + + self.kf[j] = (surface_volume_ratio_si * + rxn.get_surface_rate_coefficient(self.T.value_si, + self.surface_site_density.value_si, + self.potential.value_si + )) + + else: + if not warned and rxn.kinetics.is_pressure_dependent(): + logging.warning("Pressure may be varying, but using initial pressure to evaluate k(T,P) expressions!") + warned = True + self.kf[j] = rxn.get_rate_coefficient(self.T.value_si, P) + if rxn.reversible: + # ToDo: get_equilibrium_constant should be coverage dependent + self.Keq[j] = rxn.get_equilibrium_constant(self.T.value_si, self.potential.value_si) + kb = self.kf[j] / self.Keq[j] + self.kb[j] = kb + + def log_initial_conditions(self, number=None): + """ + Log to the console some information about this reaction system. + + Should correspond to the calculations done in set_initial_conditions. + """ + logging.info("\nSurface reaction system {}".format(number if number is not None else "")) + logging.info("Gas phase mole fractions:") + total_gas_moles = 0 + for spec, moleFrac in self.initial_gas_mole_fractions.items(): + logging.info(" {0:20s} {1:.5g}".format(spec, moleFrac)) + total_gas_moles += moleFrac + logging.info("Total gas phase: {:.3g} moles".format(total_gas_moles)) + logging.info("Pressure: {:.3g} Pa".format(self.P_initial.value_si)) + logging.info("Temperature: {} K".format(self.T.value_si)) + logging.info("Potential: {} V".format(self.potential.value_si)) + V = constants.R * self.T.value_si * total_gas_moles / self.P_initial.value_si + logging.info("Reactor volume: {:.3g} m3".format(V)) + surface_volume_ratio_si = self.surface_volume_ratio.value_si # 1/m + logging.info("Surface/volume ratio: {:.3g} m2/m3".format(surface_volume_ratio_si)) + logging.info("Surface site density: {:.3g} mol/m2".format(self.surface_site_density.value_si)) + total_surface_sites = V * surface_volume_ratio_si * self.surface_site_density.value_si # total surface sites in reactor + logging.info("Surface sites in reactor: {:.3g} moles".format(total_surface_sites)) + logging.info("Initial surface coverages (and amounts):") + for spec, coverage in self.initial_surface_coverages.items(): + logging.info(" {:18s} {:.5g} = {:.5g} moles".format(spec, coverage, total_surface_sites * coverage)) + + def set_initial_conditions(self): + """ + Sets the initial conditions of the rate equations that represent the + current reactor model. + + The volume is set to the value in m3 required to contain + one mole total of gas phase core species at start. + + The total surface sites are calculated from surface_volume_ratio and surface_site_density + allowing initial_surface_coverages to determine the number of moles of surface species. + The number of moles of gas phase species is taken from initial_gas_mole_fractions. + + The core_species_concentrations array is then determined, in mol/m3 for gas phase + and mol/m2 for surface species. + + The initial number of moles of a species j in the reactor is computed and stored in the + y0 instance attribute. + + """ + ### + ### WARNING -- When updating this method, please be sure + ### to also update the log_initial_conditions above + ### which unfortunately is maintained separately. + ### + cdef double V, P, surface_volume_ratio_si + + ReactionSystem.set_initial_conditions(self) + # self.y0 tracks number of moles of each core species + + # add only the gas phase species first + self.y0 *= 0. + for spec, moleFrac in self.initial_gas_mole_fractions.items(): + i = self.get_species_index(spec) + self.y0[i] = moleFrac # moles in reactor + + # Use ideal gas law to compute reactor volume + V = constants.R * self.T.value_si * np.sum(self.y0[:self.num_core_species]) / self.P_initial.value_si + self.V = V # volume in m^3 (assume reactor volume is gas phase volume, i.e electrode takes no space) + + surface_volume_ratio_si = self.surface_volume_ratio.value_si # 1/m + total_surface_sites = V * surface_volume_ratio_si * self.surface_site_density.value_si # total surface sites in reactor + + for spec, coverage in self.initial_surface_coverages.items(): + i = self.get_species_index(spec) + self.y0[i] = total_surface_sites * coverage # moles in reactor + + for j, isSurfaceSpecies in enumerate(self.species_on_electrode): # should only go up to core species + if isSurfaceSpecies: + self.core_species_concentrations[j] = self.y0[j] / V / surface_volume_ratio_si # moles per m2 of surface + else: + self.core_species_concentrations[j] = self.y0[j] / V # moles per m3 of gas + + def compute_network_variables(self, pdep_networks=None): + # ToDo: this should allow pressure to vary? + # for now, just call the base class version. + ReactionSystem.compute_network_variables(self, pdep_networks) + + def get_threshold_rate_constants(self, model_settings): + """ + Get the threshold rate constants for reaction filtering. + """ + raise NotImplementedError("filter_reactions=True for SurfaceReactor") + # Set the maximum unimolecular rate to be kB*T/h + unimolecular_threshold_rate_constant = 2.08366122e10 * self.T.value_si + # Set the maximum bi/trimolecular rate by using the user-defined rate constant threshold + bimolecular_threshold_rate_constant = model_settings.filter_threshold + # Maximum trimolecular rate constants are approximately three + # orders of magnitude smaller (accounting for the unit + # conversion from m^3/mol/s to m^6/mol^2/s) based on + # extending the Smoluchowski equation to three molecules + trimolecular_threshold_rate_constant = model_settings.filter_threshold / 1e3 + return (unimolecular_threshold_rate_constant, + bimolecular_threshold_rate_constant, + trimolecular_threshold_rate_constant) + + @cython.boundscheck(False) + def residual(self, + double t, + np.ndarray[np.float64_t, ndim=1] N, + np.ndarray[np.float64_t, ndim=1] dNdt, + np.ndarray[np.float64_t, ndim=1] senpar = np.zeros(1, np.float64) + ): + + """ + Return the residual function for the governing DAE system for the + simple reaction system. + """ + cdef np.ndarray[np.int_t, ndim=2] ir, ip, inet + cdef np.ndarray[np.int_t, ndim=1] reactions_on_electrode, species_on_electrode + cdef np.ndarray[np.float64_t, ndim=1] res, kf, kr, knet, delta, equilibrium_constants + cdef int num_core_species, num_core_reactions, num_edge_species, num_edge_reactions, num_pdep_networks + cdef int i, j, z, first, second, third + cdef double k, V, reaction_rate, surface_volume_ratio_si + cdef np.ndarray[np.float64_t, ndim=1] core_species_concentrations, core_species_rates, core_reaction_rates + cdef np.ndarray[np.float64_t, ndim=1] edge_species_rates, edge_reaction_rates, network_leak_rates + cdef np.ndarray[np.float64_t, ndim=1] C + cdef np.ndarray[np.float64_t, ndim=2] jacobian, dgdk + + ir = self.reactant_indices + ip = self.product_indices + equilibrium_constants = self.Keq + + kf = self.kf # are already 'per m3 of reactor' even for surface reactions + kr = self.kb # are already 'per m3 of reactor' even for surface reactions + + inet = self.network_indices + knet = self.network_leak_coefficients + + num_core_species = len(self.core_species_rates) + num_core_reactions = len(self.core_reaction_rates) + num_edge_species = len(self.edge_species_rates) + num_edge_reactions = len(self.edge_reaction_rates) + num_pdep_networks = len(self.network_leak_rates) + + res = np.zeros(num_core_species, np.float64) + + core_species_concentrations = np.zeros_like(self.core_species_concentrations) + core_species_rates = np.zeros_like(self.core_species_rates) + core_reaction_rates = np.zeros_like(self.core_reaction_rates) + edge_species_rates = np.zeros_like(self.edge_species_rates) + edge_reaction_rates = np.zeros_like(self.edge_reaction_rates) + network_leak_rates = np.zeros_like(self.network_leak_rates) + + reactions_on_electrode = self.reactions_on_electrode + species_on_electrode = self.species_on_electrode + surface_volume_ratio_si = self.surface_volume_ratio.value_si + + C = np.zeros_like(self.core_species_concentrations) + V = self.V # constant volume reactor + + for j in range(num_core_species): + if species_on_electrode[j]: + C[j] = (N[j] / V) / surface_volume_ratio_si + else: + C[j] = N[j] / V + #: surface species are in mol/m2, gas phase are in mol/m3 + core_species_concentrations[j] = C[j] + + for j in range(ir.shape[0]): + k = kf[j] + if ir[j, 0] >= num_core_species or ir[j, 1] >= num_core_species or ir[j, 2] >= num_core_species: + reaction_rate = 0.0 + elif ir[j, 1] == -1: # only one reactant + reaction_rate = k * C[ir[j, 0]] + elif ir[j, 2] == -1: # only two reactants + reaction_rate = k * C[ir[j, 0]] * C[ir[j, 1]] + else: # three reactants!! (really?) + reaction_rate = k * C[ir[j, 0]] * C[ir[j, 1]] * C[ir[j, 2]] + k = kr[j] + if ip[j, 0] >= num_core_species or ip[j, 1] >= num_core_species or ip[j, 2] >= num_core_species: + pass + elif ip[j, 1] == -1: # only one reactant + reaction_rate -= k * C[ip[j, 0]] + elif ip[j, 2] == -1: # only two reactants + reaction_rate -= k * C[ip[j, 0]] * C[ip[j, 1]] + else: # three reactants!! (really?) + reaction_rate -= k * C[ip[j, 0]] * C[ip[j, 1]] * C[ip[j, 2]] + + "reaction_rate is now in mol/m3/s" + # Set the reaction and species rates + if j < num_core_reactions: + # The reaction is a core reaction + core_reaction_rates[j] = reaction_rate + + # Add/subtract the total reaction rate from each species rate + # Since it's a core reaction we know that all of its reactants + # and products are core species + first = ir[j, 0] + core_species_rates[first] -= reaction_rate + second = ir[j, 1] + if second != -1: + core_species_rates[second] -= reaction_rate + third = ir[j, 2] + if third != -1: + core_species_rates[third] -= reaction_rate + first = ip[j, 0] + core_species_rates[first] += reaction_rate + second = ip[j, 1] + if second != -1: + core_species_rates[second] += reaction_rate + third = ip[j, 2] + if third != -1: + core_species_rates[third] += reaction_rate + + else: + # The reaction is an edge reaction + edge_reaction_rates[j - num_core_reactions] = reaction_rate + + # Add/substract the total reaction rate from each species rate + # Since it's an edge reaction its reactants and products could + # be either core or edge species + # We're only interested in the edge species + first = ir[j, 0] + if first >= num_core_species: edge_species_rates[first - num_core_species] -= reaction_rate + second = ir[j, 1] + if second != -1: + if second >= num_core_species: edge_species_rates[second - num_core_species] -= reaction_rate + third = ir[j, 2] + if third != -1: + if third >= num_core_species: edge_species_rates[third - num_core_species] -= reaction_rate + first = ip[j, 0] + if first >= num_core_species: edge_species_rates[first - num_core_species] += reaction_rate + second = ip[j, 1] + if second != -1: + if second >= num_core_species: edge_species_rates[second - num_core_species] += reaction_rate + third = ip[j, 2] + if third != -1: + if third >= num_core_species: edge_species_rates[third - num_core_species] += reaction_rate + + for j in range(inet.shape[0]): + if inet[j, 0] != -1: #all source species are in the core + k = knet[j] + if inet[j, 1] == -1: # only one reactant + reaction_rate = k * C[inet[j, 0]] + elif inet[j, 2] == -1: # only two reactants + reaction_rate = k * C[inet[j, 0]] * C[inet[j, 1]] + else: # three reactants + reaction_rate = k * C[inet[j, 0]] * C[inet[j, 1]] * C[inet[j, 2]] + network_leak_rates[j] = reaction_rate + else: + network_leak_rates[j] = 0.0 + + self.core_species_concentrations = core_species_concentrations + self.core_species_rates = core_species_rates + self.core_reaction_rates = core_reaction_rates + self.edge_species_rates = edge_species_rates + self.edge_reaction_rates = edge_reaction_rates + self.network_leak_rates = network_leak_rates + + res = core_species_rates * V + # mol/s + + if self.sensitivity and False: + delta = np.zeros(len(N), np.float64) + delta[:num_core_species] = res + if self.jacobian_matrix is None: + jacobian = self.jacobian(t, N, dNdt, 0, senpar) + else: + jacobian = self.jacobian_matrix + dgdk = ReactionSystem.compute_rate_derivative(self) + for j in range(num_core_reactions + num_core_species): + for i in range(num_core_species): + for z in range(num_core_species): + delta[(j + 1) * num_core_species + i] += jacobian[i, z] * N[(j + 1) * num_core_species + z] + delta[(j + 1) * num_core_species + i] += dgdk[i, j] + else: + delta = res + delta = delta - dNdt + + # Return DELTA, IRES. IRES is set to 1 in order to tell DASPK to evaluate the sensitivity residuals + return delta, 1 + diff --git a/rmgpy/species.pxd b/rmgpy/species.pxd index b949b24182..9d3a95448c 100644 --- a/rmgpy/species.pxd +++ b/rmgpy/species.pxd @@ -57,6 +57,8 @@ cdef class Species: cdef str _inchi cdef str _smiles + cpdef get_net_charge(self) + cpdef generate_resonance_structures(self, bint keep_isomorphic=?, bint filter_structures=?, bint save_order=?) cpdef bint is_isomorphic(self, other, bint generate_initial_map=?, bint save_order=?, bint strict=?) except -2 @@ -75,6 +77,10 @@ cdef class Species: cpdef bint is_surface_site(self) except -2 + cpdef bint is_electron(self) except -2 + + cpdef bint is_proton(self) except -2 + cpdef bint has_statmech(self) except -2 cpdef bint has_thermo(self) except -2 diff --git a/rmgpy/species.py b/rmgpy/species.py index c02f5b2ee3..7ddc7c42f8 100644 --- a/rmgpy/species.py +++ b/rmgpy/species.py @@ -263,6 +263,14 @@ def molecular_weight(self): def molecular_weight(self, value): self._molecular_weight = quantity.Mass(value) + def get_net_charge(self): + """ + Iterate through the atoms in the structure and calculate the net charge + on the overall molecule. + """ + + return self.molecule[0].get_net_charge() + def generate_resonance_structures(self, keep_isomorphic=True, filter_structures=True, save_order=False): """ Generate all of the resonance structures of this species. The isomers are @@ -343,7 +351,7 @@ def is_structure_in_list(self, species_list): ' should be a List of Species objects.'.format(species)) return False - def from_adjacency_list(self, adjlist, raise_atomtype_exception=True, raise_charge_exception=True): + def from_adjacency_list(self, adjlist, raise_atomtype_exception=True, raise_charge_exception=False): """ Load the structure of a species as a :class:`Molecule` object from the given adjacency list `adjlist` and store it as the first entry of a @@ -454,6 +462,22 @@ def is_surface_site(self): """Return ``True`` if the species is a vacant surface site.""" return self.molecule[0].is_surface_site() + def is_electron(self): + """Return ``True`` if the species is an electron""" + + if len(self.molecule) == 0: + return False + else: + return self.molecule[0].is_electron() + + def is_proton(self): + """Return ``True`` if the species is a proton""" + + if len(self.molecule) == 0: + return False + else: + return self.molecule[0].is_proton() + def get_partition_function(self, T): """ Return the partition function for the species at the specified diff --git a/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/groups.py b/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/groups.py new file mode 100644 index 0000000000..656846b7f4 --- /dev/null +++ b/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/groups.py @@ -0,0 +1,211 @@ +#!/usr/bin/env python +# encoding: utf-8 + +name = "Surface_Proton_Electron_Reduction_Alpha/groups" +shortDesc = u"" +longDesc = u""" + + *1 *1-*3H + || + *3H+ + *e- ----> | + ~*2~ ~*2~~ + +The rate, which should be in mol/m2/s, +will be given by k * (mol/m2) * (mol/m3) * 1 +so k should be in (m3/mol/s). +""" + +template(reactants=["Adsorbate", "Proton"], products=["Reduced"], ownReverse=False) + +reverse = "Surface_Proton_Electron_Oxidation_Alpha" + +reactantNum = 2 +productNum = 1 +allowChargedSpecies = True +electrons = -1 + +recipe(actions=[ + ['LOSE_CHARGE', '*3', 1], + ['CHANGE_BOND', '*1', -1, '*2'], + ['FORM_BOND', '*1', 1, '*3'], +]) + +entry( + index = 1, + label = "Adsorbate", + group = +""" +1 *1 R!H u0 {2,[D,T,Q]} +2 *2 X u0 {1,[D,T,Q]} +""", + kinetics = None, +) + +entry( + index = 2, + label = "Proton", + group = +""" +1 *3 H+ u0 p0 c+1 +""", + kinetics = None, +) + +entry( + index = 4, + label = "CX", + group = +""" +1 *1 C u0 {2,[D,T,Q]} +2 *2 X u0 {1,[D,T,Q]} +""", + kinetics = None, +) + +entry( + index = 5, + label = "CTX", + group = +""" +1 *1 C u0 {2,T} +2 *2 X u0 {1,T} +""", + kinetics = None, +) + +entry( + index = 6, + label = "HCX", + group = +""" +1 *1 C u0 {2,T} {3,S} +2 *2 X u0 {1,T} +3 H u0 {1,S} +""", + kinetics = None, +) + +entry( + index = 7, + label = "C=X", + group = +""" +1 *1 C u0 {2,D} +2 *2 X u0 {1,D} +""", + kinetics = None, +) + +entry( + index = 8, + label = "H2C=X", + group = +""" +1 *1 C u0 {2,D} {3,S} {4,S} +2 *2 X u0 {1,D} +3 H u0 {1,S} +4 H u0 {1,S} +""", + kinetics = None, +) + +entry( + index = 9, + label = "O=C=X", + group = +""" +1 *1 C u0 {2,D} {3,D} +2 *2 X u0 {1,D} +3 O2d u0 {1,D} +""", + kinetics = None, +) + + +entry( + index = 10, + label = "OX", + group = +""" +1 *1 O u0 {2,D} +2 *2 X u0 {1,D} +""", + kinetics = None, +) + + +entry( + index = 11, + label = "NX", + group = +""" +1 *1 N u0 {2,[D,T]} +2 *2 X u0 {1,[D,T]} +""", + kinetics = None, +) + +entry( + index = 12, + label = "NTX", + group = +""" +1 *1 N u0 {2,T} +2 *2 X u0 {1,T} +""", + kinetics = None, +) + +entry( + index = 13, + label = "N=X", + group = +""" +1 *1 N u0 {2,D} +2 *2 X u0 {1,D} +""", + kinetics = None, +) + +entry( + index = 14, + label = "HN=X", + group = +""" +1 *1 N u0 {2,D} {3,S} +2 *2 X u0 {1,D} +3 N u0 {1,S} +""", + kinetics = None, +) + +entry( + index = 15, + label = "N-N=X", + group = +""" +1 *1 N u0 {2,D} {3,S} +2 *2 X u0 {1,D} +3 N u0 {1,S} +""", + kinetics = None, +) + +tree( +""" +L1: Adsorbate + L2: CX + L3: CTX + L4: HCX + L3: C=X + L4: O=C=X + L4: H2C=X + L2: OX + L2: NX + L3: NTX + L3: N=X + L4: HN=X + L4: N-N=X + +L1: Proton +""" +) diff --git a/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/rules.py b/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/rules.py new file mode 100644 index 0000000000..67a06429f2 --- /dev/null +++ b/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/rules.py @@ -0,0 +1,26 @@ +#!/usr/bin/env python +# encoding: utf-8 + +name = "Surface_Proton_Electron_Reduction_Alpha/rules" +shortDesc = u"" +longDesc = u""" +Surface adsorption of a single radical forming a single bond to the surface site +""" + +# entry( +# index = 1, +# label = "Adsorbate;Proton;Electron", +# kinetics = SurfaceChargeTransfer( +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff, 0 default +# V0 = None, # Reference potential +# Ea = (15, 'kJ/mol'), # activation energy at the reversible potential +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# rank = 0, +# shortDesc = u"""Default""", +# longDesc = u"""https://doi.org/10.1021/jp4100608""" +# ) + diff --git a/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/training/dictionary.txt b/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/training/dictionary.txt new file mode 100644 index 0000000000..c4d8a3f3e8 --- /dev/null +++ b/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/training/dictionary.txt @@ -0,0 +1,64 @@ +H +1 *3 H u0 p0 c+1 + +NX +1 *1 N u0 p1 c0 {2,T} +2 *2 X u0 p0 c0 {1,T} + +HNX +1 *1 N u0 p1 c0 {2,D} {3,S} +2 *2 X u0 p0 c0 {1,D} +3 H u0 p0 c0 {1,S} + +HNX_p +1 *1 N u0 p1 c0 {2,D} {3,S} +2 *2 X u0 p0 c0 {1,D} +3 *3 H u0 p0 c0 {1,S} + +H2NX +1 *1 N u0 p1 c0 {2,S} {3,S} {4,S} +2 *2 X u0 p0 c0 {1,S} +3 *3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} + +CX +1 *1 C u0 p0 c0 {2,Q} +2 *2 X u0 p0 c0 {1,Q} + +CHX +1 *1 C u0 p0 c0 {2,T} {3,S} +2 *2 X u0 p0 c0 {1,T} +3 H u0 p0 c0 {1,S} + +CHX_p +1 *1 C u0 p0 c0 {2,T} {3,S} +2 *2 X u0 p0 c0 {1,T} +3 *3 H u0 p0 c0 {1,S} + +CH2X +1 *1 C u0 p0 c0 {2,D} {3,S} {4,S} +2 *2 X u0 p0 c0 {1,D} +3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} + +CH2X_p +1 *1 C u0 p0 c0 {2,D} {3,S} {4,S} +2 *2 X u0 p0 c0 {1,D} +3 *3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} + +CH3X +1 *1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} +2 *2 X u0 p0 c0 {1,S} +3 *3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} +5 H u0 p0 c0 {1,S} + +OX +1 *1 O u0 p2 c0 {2,D} +2 *2 X u0 p0 c0 {1,D} + +HOX +1 *1 O u0 p2 c0 {2,S} {3,S} +2 *2 X u0 p0 c0 {1,S} +3 *3 H u0 p0 c0 {1,S} \ No newline at end of file diff --git a/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/training/reactions.py b/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/training/reactions.py new file mode 100644 index 0000000000..3866fa74ac --- /dev/null +++ b/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/training/reactions.py @@ -0,0 +1,516 @@ +#!/usr/bin/env python +# encoding: utf-8 + +name = "Surface_Proton_Electron_Reduction_Alpha/training" +shortDesc = u"Reaction kinetics used to generate rate rules" +longDesc = u""" +Put kinetic parameters for specific reactions in this file to use as a +training set for generating rate rules to populate this kinetics family. +""" + +# entry( +# index = 1, +# label = "CX + H <=> CHX_p", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.20, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (0.24, 'V'), # reference potential +# Ea = (0.61, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Tafel""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +entry( + index = 2, + label = "CX + H <=> CHX_p", + degeneracy = 1, + kinetics = SurfaceChargeTransfer( + alpha = 0.20, # charge transfer coeff + A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ + n = 0, # temperature coeff + V0 = (-0.5, 'V'), # reference potential + Ea = (0.44, 'eV/molecule'), # activation energy + Tmin = (200, 'K'), + Tmax = (3000, 'K'), + electrons = -1, # electron stochiometric coeff + ), + shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", + longDesc = u"""Tafel""", + metal = "Pt", + facet = "111", + site = "", + rank = 5, +) + +# entry( +# index = 1, +# label = "CX + H <=> CHX_p", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.06, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (0.29, 'V'), # reference potential +# Ea = (0.19, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Heyrovsky""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +# entry( +# index = 1, +# label = "CX + H <=> CHX_p", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.06, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (-0.5, 'V'), # reference potential +# Ea = (0.13, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Heyrovsky""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +# entry( +# index = 3, +# label = "CHX + H <=> CH2X_p", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.31, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (0.32, 'V'), # reference potential +# Ea = (0.77, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Tafel""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +entry( + index = 3, + label = "CHX + H <=> CH2X_p", + degeneracy = 1, + kinetics = SurfaceChargeTransfer( + alpha = 0.31, # charge transfer coeff + A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ + n = 0, # temperature coeff + V0 = (-0.5, 'V'), # reference potential + Ea = (0.44, 'eV/molecule'), # activation energy + Tmin = (200, 'K'), + Tmax = (3000, 'K'), + electrons = -1, # electron stochiometric coeff + ), + shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", + longDesc = u"""Tafel""", + metal = "Pt", + facet = "111", + site = "", + rank = 5, +) + +# entry( +# index = 3, +# label = "CHX + H <=> CH2X_p", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.05, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (0.30, 'V'), # reference potential +# Ea = (0.59, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Heyrovsky""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +# entry( +# index = 4, +# label = "CHX + H <=> CH2X_p", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.05, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (-0.5, 'V'), # reference potential +# Ea = (0.53, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Heyrovsky""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +# entry( +# index = 6, +# label = "CH2X + H <=> CH3X", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.23, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (0.38, 'V'), # reference potential +# Ea = (0.62, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Tafel""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +entry( + index = 4, + label = "CH2X + H <=> CH3X", + degeneracy = 1, + kinetics = SurfaceChargeTransfer( + alpha = 0.23, # charge transfer coeff + A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ + n = 0, # temperature coeff + V0 = (-0.5, 'V'), # reference potential + Ea = (0.37, 'eV/molecule'), # activation energy + Tmin = (200, 'K'), + Tmax = (3000, 'K'), + electrons = -1, # electron stochiometric coeff + ), + shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", + longDesc = u"""Tafel""", + metal = "Pt", + facet = "111", + site = "", + rank = 5, +) + +# entry( +# index = 8, +# label = "NX + H <=> HNX", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.07, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (-0.12, 'V'), # reference potential +# Ea = (0.15, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2017.01.050""", +# longDesc = u""" +# """, +# metal = "Cu", +# facet = "111", +# site = "", +# rank = 5, +# ) + +# entry( +# index = 9, +# label = "NX + H <=> HNX_p", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.23, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (0.24, 'V'), # reference potential +# Ea = (0.78, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Tafel""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +entry( + index = 5, + label = "NX + H <=> HNX_p", + degeneracy = 1, + kinetics = SurfaceChargeTransfer( + alpha = 0.23, # charge transfer coeff + A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ + n = 0, # temperature coeff + V0 = (-0.5, 'V'), # reference potential + Ea = (0.59, 'eV/molecule'), # activation energy + Tmin = (200, 'K'), + Tmax = (3000, 'K'), + electrons = -1, # electron stochiometric coeff + ), + shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", + longDesc = u"""Tafel""", + metal = "Pt", + facet = "111", + site = "", + rank = 5, +) + +# entry( +# index = 10, +# label = "NX + H <=> HNX_p", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.07, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (0.3, 'V'), # reference potential +# Ea = (0.17, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Heyrovsky""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +# entry( +# index = 4, +# label = "NX + H <=> HNX_p", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.07, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (-0.5, 'V'), # reference potential +# Ea = (0.09, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Heyrovsky""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +# entry( +# index = 11, +# label = "HNX + H <=> H2NX", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.27, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (0.28, 'V'), # reference potential +# Ea = (1.20, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Tafel""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +entry( + index = 6, + label = "HNX + H <=> H2NX", + degeneracy = 1, + kinetics = SurfaceChargeTransfer( + alpha = 0.27, # charge transfer coeff + A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ + n = 0, # temperature coeff + V0 = (-0.5, 'V'), # reference potential + Ea = (0.97, 'eV/molecule'), # activation energy + Tmin = (200, 'K'), + Tmax = (3000, 'K'), + electrons = -1, # electron stochiometric coeff + ), + shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", + longDesc = u"""Tafel""", + metal = "Pt", + facet = "111", + site = "", + rank = 5, +) + +# entry( +# index = 11, +# label = "HNX + H <=> H2NX", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.64, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (0.31, 'V'), # reference potential +# Ea = (0.99, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Heyrovsky""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +# entry( +# index = 5, +# label = "HNX + H <=> H2NX", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.64, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (-0.5, 'V'), # reference potential +# Ea = (0.43, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Heyrovsky""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +# entry( +# index = 12, +# label = "OX + H <=> HOX", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.41, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (0.31, 'V'), # reference potential +# Ea = (0.87, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Tafel""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +entry( + index = 7, + label = "OX + H <=> HOX", + degeneracy = 1, + kinetics = SurfaceChargeTransfer( + alpha = 0.42, # charge transfer coeff + A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ + n = 0, # temperature coeff + V0 = (-0.5, 'V'), # reference potential + Ea = (0.48, 'eV/molecule'), # activation energy + Tmin = (200, 'K'), + Tmax = (3000, 'K'), + electrons = -1, # electron stochiometric coeff + ), + shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", + longDesc = u"""Tafel""", + metal = "Pt", + facet = "111", + site = "", + rank = 5, +) + +# entry( +# index = 12, +# label = "OX + H <=> HOX", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.02, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (0.57, 'V'), # reference potential +# Ea = (0.06, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Heyrovsky""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) + +# entry( +# index = 6, +# label = "OX + H <=> HOX", +# degeneracy = 1, +# kinetics = SurfaceChargeTransfer( +# alpha = 0.02, # charge transfer coeff +# A = (2.5e14, 'm^3/(mol*s)'), # pre-exponential factor estimate 10^11 s^-1 * 2.5e5 m^2/mol / 1000 m^3/mol H+ +# n = 0, # temperature coeff +# V0 = (-0.5, 'V'), # reference potential +# Ea = (0.03, 'eV/molecule'), # activation energy +# Tmin = (200, 'K'), +# Tmax = (3000, 'K'), +# electrons = -1, # electron stochiometric coeff +# ), +# shortDesc = u"""https://doi.org/10.1016/j.cattod.2018.03.048""", +# longDesc = u"""Heyrovsky""", +# metal = "Pt", +# facet = "111", +# site = "", +# rank = 5, +# ) diff --git a/rmgpy/yml.py b/rmgpy/yml.py index d1873805f8..1b10c10432 100644 --- a/rmgpy/yml.py +++ b/rmgpy/yml.py @@ -44,7 +44,7 @@ from rmgpy.kinetics.falloff import Troe, ThirdBody, Lindemann from rmgpy.kinetics.chebyshev import Chebyshev from rmgpy.data.solvation import SolventData -from rmgpy.kinetics.surface import StickingCoefficient +from rmgpy.kinetics.surface import StickingCoefficient, SurfaceChargeTransfer from rmgpy.util import make_output_subdirectory @@ -136,6 +136,14 @@ def obj_to_dict(obj, spcs, names=None, label="solvent"): result_dict["A"] = obj.A.value_si result_dict["Ea"] = obj.Ea.value_si result_dict["n"] = obj.n.value_si + elif isinstance(obj, SurfaceChargeTransfer): + result_dict["type"] = "SurfaceChargeTransfer" + result_dict["A"] = obj.A.value_si + result_dict["Ea"] = obj.Ea.value_si + result_dict["n"] = obj.n.value_si + result_dict["electrons"] = obj.electrons.value_si + result_dict["V0"] = obj.V0.value_si + result_dict["alpha"] = obj.alpha.value_si elif isinstance(obj, StickingCoefficient): obj.change_t0(1.0) result_dict["type"] = "StickingCoefficient" diff --git a/setup.py b/setup.py index c09e9d07c3..a040817c0c 100644 --- a/setup.py +++ b/setup.py @@ -127,6 +127,7 @@ Extension('rmgpy.solver.liquid', ['rmgpy/solver/liquid.pyx'], include_dirs=['.']), Extension('rmgpy.solver.mbSampled', ['rmgpy/solver/mbSampled.pyx'], include_dirs=['.']), Extension('rmgpy.solver.surface', ['rmgpy/solver/surface.pyx'], include_dirs=['.']), + Extension('rmgpy.solver.electrode', ['rmgpy/solver/electrode.pyx'], include_dirs=['.']), ] arkane_ext_modules = [