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..b6ea356fb7 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -42,6 +42,7 @@ from rmgpy.data.reference import Reference, Article, Book, Thesis from rmgpy.exceptions import DatabaseError, InvalidAdjacencyListError from rmgpy.kinetics.uncertainties import RateUncertainty +from rmgpy.kinetics.arrhenius import ArrheniusChargeTransfer, ArrheniusChargeTransferBM from rmgpy.molecule import Molecule, Group @@ -228,6 +229,8 @@ def load(self, path, local_context=None, global_context=None): local_context['shortDesc'] = self.short_desc local_context['longDesc'] = self.long_desc local_context['RateUncertainty'] = RateUncertainty + local_context['ArrheniusChargeTransfer'] = ArrheniusChargeTransfer + local_context['ArrheniusChargeTransferBM'] = ArrheniusChargeTransferBM local_context['metal'] = self.metal local_context['site'] = self.site local_context['facet'] = self.facet diff --git a/rmgpy/data/kinetics/database.py b/rmgpy/data/kinetics/database.py index d37e0de432..85c4794f9d 100644 --- a/rmgpy/data/kinetics/database.py +++ b/rmgpy/data/kinetics/database.py @@ -48,7 +48,7 @@ from rmgpy.molecule import Molecule, Group from rmgpy.reaction import Reaction, same_species_lists from rmgpy.species import Species - +from rmgpy.data.solvation import SoluteData ################################################################################ @@ -79,7 +79,8 @@ def __init__(self): 'SurfaceArrhenius': SurfaceArrhenius, 'SurfaceArrheniusBEP': SurfaceArrheniusBEP, 'R': constants.R, - 'ArrheniusBM': ArrheniusBM + 'ArrheniusBM': ArrheniusBM, + 'SoluteData': SoluteData, } self.global_context = {} diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index dc05acaad5..db7db5108a 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -55,7 +55,8 @@ from rmgpy.exceptions import ActionError, DatabaseError, InvalidActionError, KekulizationError, KineticsError, \ ForbiddenStructureException, UndeterminableKineticsError from rmgpy.kinetics import Arrhenius, SurfaceArrhenius, SurfaceArrheniusBEP, StickingCoefficient, \ - StickingCoefficientBEP, ArrheniusBM + StickingCoefficientBEP, ArrheniusBM, SurfaceChargeTransfer, ArrheniusChargeTransfer, \ + ArrheniusChargeTransferBM from rmgpy.kinetics.uncertainties import RateUncertainty, rank_accuracy_map from rmgpy.molecule import Bond, GroupBond, Group, Molecule from rmgpy.molecule.atomtype import ATOMTYPES @@ -308,7 +309,7 @@ def _apply(self, struct, forward, unique): if info < 1: raise InvalidActionError('Attempted to change a nonexistent bond.') # If we do not have a bond, it might be because we are trying to change a vdW bond - # Lets see if one of that atoms is a surface site, + # Lets see if one of that atoms is a surface site, # If we have a surface site, we will make a single bond, then change it by info - 1 is_vdW_bond = False for atom in (atom1, atom2): @@ -444,8 +445,8 @@ def apply_reverse(self, struct, unique=True): class KineticsFamily(Database): """ - A class for working with an RMG kinetics family: a set of reactions with - similar chemistry, and therefore similar reaction rates. The attributes + A class for working with an RMG kinetics family: a set of reactions with + similar chemistry, and therefore similar reaction rates. The attributes are: =================== =============================== ======================== @@ -699,11 +700,11 @@ def distribute_tree_distances(self): def load(self, path, local_context=None, global_context=None, depository_labels=None): """ Load a kinetics database from a file located at `path` on disk. - + If `depository_labels` is a list, eg. ['training','PrIMe'], then only those depositories are loaded, and they are searched in that order when generating kinetics. - + If depository_labels is None then load 'training' first then everything else. If depository_labels is not None then load in the order specified in depository_labels. """ @@ -796,7 +797,7 @@ def load(self, path, local_context=None, global_context=None, depository_labels= # depository and add them to the RMG rate rules by default: depository_labels = ['training'] if depository_labels: - # If there are depository labels, load them in the order specified, but + # If there are depository labels, load them in the order specified, but # append the training reactions unless the user specifically declares it not # to be included with a '!training' flag if '!training' not in depository_labels: @@ -835,7 +836,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. ' @@ -862,12 +864,12 @@ def save_training_reactions(self, reactions, reference=None, reference_type='', rank=3): """ This function takes a list of reactions appends it to the training reactions file. It ignores the existence of - duplicate reactions. - - The rank for each new reaction's kinetics is set to a default value of 3 unless the user specifies differently + duplicate reactions. + + The rank for each new reaction's kinetics is set to a default value of 3 unless the user specifies differently for those reactions. - - For each entry, the long description is imported from the kinetics comment. + + For each entry, the long description is imported from the kinetics comment. """ if not isinstance(reference, list): @@ -970,7 +972,7 @@ def save_training_reactions(self, reactions, reference=None, reference_type='', def save(self, path): """ - Save the current database to the file at location `path` on disk. + Save the current database to the file at location `path` on disk. """ self.save_groups(os.path.join(path, 'groups.py')) self.rules.save(os.path.join(path, 'rules.py')) @@ -987,7 +989,7 @@ def save_depository(self, depository, path): def save_groups(self, path): """ - Save the current database to the file at location `path` on disk. + Save the current database to the file at location `path` on disk. """ entries = self.groups.get_entries_to_save() @@ -1023,7 +1025,7 @@ def save_groups(self, path): f.write('reactantNum = {0}\n\n'.format(self.reactant_num)) if self.product_num is not None: f.write('productNum = {0}\n\n'.format(self.product_num)) - + if self.auto_generated is not None: f.write('autoGenerated = {0}\n\n'.format(self.auto_generated)) @@ -1147,14 +1149,14 @@ def generate_product_template(self, reactants0): def has_rate_rule(self, template): """ - Return ``True`` if a rate rule with the given `template` currently + Return ``True`` if a rate rule with the given `template` currently exists, or ``False`` otherwise. """ return self.rules.has_rule(template) def get_rate_rule(self, template): """ - Return the rate rule with the given `template`. Raises a + Return the rate rule with the given `template`. Raises a :class:`ValueError` if no corresponding entry exists. """ entry = self.rules.get_rule(template) @@ -1340,7 +1342,7 @@ def add_rules_from_training(self, thermo_database=None, train_indices=None): def get_root_template(self): """ Return the root template for the reaction family. Most of the time this - is the top-level nodes of the tree (as stored in the + is the top-level nodes of the tree (as stored in the :class:`KineticsGroups` object), but there are a few exceptions (e.g. R_Recombination). """ @@ -1493,7 +1495,10 @@ def apply_recipe(self, reactant_structures, forward=True, unique=True, relabel_a product_num = self.product_num or len(template.products) # Split product structure into multiple species if necessary - product_structures = product_structure.split() + if self.auto_generated and isinstance(reactant_structures[0],Group) and self.product_num == 1: + product_structures = [product_structure] + else: + product_structures = product_structure.split() # Make sure we've made the expected number of products if product_num != len(product_structures): @@ -1535,8 +1540,20 @@ 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)) @@ -1659,10 +1676,10 @@ def _generate_product_structures(self, reactant_structures, maps, forward, relab def is_molecule_forbidden(self, molecule): """ Return ``True`` if the molecule is forbidden in this family, or - ``False`` otherwise. + ``False`` otherwise. """ - # check family-specific forbidden structures + # check family-specific forbidden structures if self.forbidden is not None and self.forbidden.is_molecule_forbidden(molecule): return True @@ -1706,7 +1723,7 @@ def _create_reaction(self, reactants, products, is_forward): def _match_reactant_to_template(self, reactant, template_reactant): """ - Return a complete list of the mappings if the provided reactant + Return a complete list of the mappings if the provided reactant matches the provided template reactant, or an empty list if not. """ @@ -1882,8 +1899,8 @@ def calculate_degeneracy(self, reaction): For a `reaction` with `Molecule` or `Species` objects given in the direction in which the kinetics are defined, compute the reaction-path degeneracy. - This method by default adjusts for double counting of identical reactants. - This should only be adjusted once per reaction. To not adjust for + This method by default adjusts for double counting of identical reactants. + This should only be adjusted once per reaction. To not adjust for identical reactants (since you will be reducing them later in the algorithm), add `ignoreSameReactants= True` to this method. """ @@ -1957,7 +1974,7 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson will return an empty list. Each item in the list of reactants should be a list of :class:`Molecule` objects, each representing a resonance structure of the species of interest. - + This method returns all reactions, and degenerate reactions can then be found using `rmgpy.data.kinetics.common.find_degenerate_reactions`. @@ -1982,7 +1999,7 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson rxn_list = [] - # Wrap each reactant in a list if not already done (this is done to + # Wrap each reactant in a list if not already done (this is done to # allow for passing multiple resonance structures for each molecule) # This also makes a copy of the reactants list so we don't modify the # original @@ -2303,7 +2320,7 @@ def generate_products_and_reactions(order): if not forward and ('adsorption' in self.label.lower() or 'eleyrideal' in self.label.lower()): # Desorption should have desorbed something (else it was probably bidentate) # so delete reactions that don't make a gas-phase desorbed product - # Eley-Rideal reactions should have one gas-phase product in the reverse direction + # Eley-Rideal reactions should have one gas-phase product in the reverse direction # Determine how many surf reactants we expect based on the template n_surf_expected = len([r for r in self.forward_template.reactants if r.item.contains_surface_site()]) @@ -2373,7 +2390,7 @@ def get_reaction_pairs(self, reaction): """ pairs = [] if len(reaction.reactants) == 1 or len(reaction.products) == 1: - # When there is only one reactant (or one product), it is paired + # When there is only one reactant (or one product), it is paired # with each of the products (reactants) for reactant in reaction.reactants: for product in reaction.products: @@ -2542,7 +2559,7 @@ def get_kinetics_for_template(self, template, degeneracy=1, method='rate rules') `template` and reaction-path `degeneracy`. There are two possible methods to use: 'group additivity' (new possible RMG-Py behavior) and 'rate rules' (old RMG-Java behavior, and default RMG-Py behavior). - + Returns a tuple (kinetics, entry): If it's estimated via 'rate rules' and an exact match is found in the tree, then the entry is returned as the second element of the tuple. @@ -2560,7 +2577,7 @@ def get_kinetics_for_template(self, template, degeneracy=1, method='rate rules') def get_kinetics_from_depository(self, depository, reaction, template, degeneracy): """ Search the given `depository` in this kinetics family for kinetics - for the given `reaction`. Returns a list of all of the matching + for the given `reaction`. Returns a list of all of the matching kinetics, the corresponding entries, and ``True`` if the kinetics match the forward direction or ``False`` if they match the reverse direction. @@ -2632,7 +2649,7 @@ def get_kinetics(self, reaction, template_labels, degeneracy=1, estimator='', re for kinetics, entry, is_forward in kinetics_list0: kinetics_list.append([kinetics, depository, entry, is_forward]) - # If estimator type of rate rules or group additivity is given, retrieve the kinetics. + # If estimator type of rate rules or group additivity is given, retrieve the kinetics. if estimator: try: kinetics, entry = self.get_kinetics_for_template(template, degeneracy, method=estimator) @@ -2644,7 +2661,7 @@ def get_kinetics(self, reaction, template_labels, degeneracy=1, estimator='', re if not return_all_kinetics: return kinetics, estimator, entry, True kinetics_list.append([kinetics, estimator, entry, True]) - # If no estimation method was given, prioritize rate rule estimation. + # If no estimation method was given, prioritize rate rule estimation. # If returning all kinetics, add estimations from both rate rules and group additivity. else: try: @@ -2674,7 +2691,7 @@ def estimate_kinetics_using_group_additivity(self, template, degeneracy=1): """ Determine the appropriate kinetics for a reaction with the given `template` using group additivity. - + Returns just the kinetics, or None. """ warnings.warn("Group additivity is no longer supported and may be" @@ -2698,7 +2715,7 @@ def estimate_kinetics_using_rate_rules(self, template, degeneracy=1): """ Determine the appropriate kinetics for a reaction with the given `template` using rate rules. - + Returns a tuple (kinetics, entry) where `entry` is the database entry used to determine the kinetics only if it is an exact match, and is None if some averaging or use of a parent node took place. @@ -2709,8 +2726,8 @@ def estimate_kinetics_using_rate_rules(self, template, degeneracy=1): def get_reaction_template_labels(self, reaction): """ - Retrieve the template for the reaction and - return the corresponding labels for each of the + Retrieve the template for the reaction and + return the corresponding labels for each of the groups in the template. """ template = self.get_reaction_template(reaction) @@ -2723,8 +2740,8 @@ def get_reaction_template_labels(self, reaction): def retrieve_template(self, template_labels): """ - Reconstruct the groups associated with the - labels of the reaction template and + Reconstruct the groups associated with the + labels of the reaction template and return a list. """ template = [] @@ -2735,9 +2752,9 @@ def retrieve_template(self, template_labels): def get_labeled_reactants_and_products(self, reactants, products, relabel_atoms=True): """ - Given `reactants`, a list of :class:`Molecule` objects, and products, a list of - :class:`Molecule` objects, return two new lists of :class:`Molecule` objects with - atoms labeled: one for reactants, one for products. Returned molecules are totally + Given `reactants`, a list of :class:`Molecule` objects, and products, a list of + :class:`Molecule` objects, return two new lists of :class:`Molecule` objects with + atoms labeled: one for reactants, one for products. Returned molecules are totally new entities in memory so input molecules `reactants` and `products` won't be affected. If RMG cannot find appropriate labels, (None, None) will be returned. If ``relabel_atoms`` is ``True``, product atom labels of reversible families @@ -2916,7 +2933,7 @@ def add_entry(self, parent, grp, name): def _split_reactions(self, rxns, newgrp): """ divides the reactions in rxns between the new - group structure newgrp and the old structure with + group structure newgrp and the old structure with label oldlabel returns a list of reactions associated with the new group the list of reactions associated with the old group @@ -2941,14 +2958,14 @@ def _split_reactions(self, rxns, newgrp): comp.append(rxn) return new, comp, new_inds - + def reaction_matches(self, rxn, grp): rmol = rxn.reactants[0].molecule[0] for r in rxn.reactants[1:]: rmol = rmol.merge(r.molecule[0]) rmol.identify_ring_membership() return rmol.is_subgraph_isomorphic(grp, generate_initial_map=True, save_order=True) - + def eval_ext(self, parent, ext, extname, template_rxn_map, obj=None, T=1000.0): """ evaluates the objective function obj @@ -2972,15 +2989,15 @@ def get_extension_edge(self, parent, template_rxn_map, obj, T, iter_max=np.inf, finds the set of all extension groups to parent such that 1) the extension group divides the set of reactions under parent 2) No generalization of the extension group divides the set of reactions under parent - + We find this by generating all possible extensions of the initial group. Extensions that split reactions are added - to the list. All extensions that do not split reactions and do not create bonds are ignored + to the list. All extensions that do not split reactions and do not create bonds are ignored (although those that match every reaction are labeled so we don't search them twice). Those that match - all reactions and involve bond creation undergo this process again. - - Principle: Say you have two elementary changes to a group ext1 and ext2 if applying ext1 and ext2 results in a + all reactions and involve bond creation undergo this process again. + + Principle: Say you have two elementary changes to a group ext1 and ext2 if applying ext1 and ext2 results in a split at least one of ext1 and ext2 must result in a split - + Speed of this algorithm relies heavily on searching non bond creation dimensions once. """ out_exts = [[]] @@ -2991,7 +3008,7 @@ def get_extension_edge(self, parent, template_rxn_map, obj, T, iter_max=np.inf, n_splits = len(template_rxn_map[parent.label][0].reactants) iter = 0 - + while grps[iter] != []: grp = grps[iter][-1] @@ -3110,7 +3127,7 @@ def get_extension_edge(self, parent, template_rxn_map, obj, T, iter_max=np.inf, out_exts.append([]) grps[iter].pop() names.pop() - + for ind in ext_inds: # collect the groups to be expanded grpr, grpcr, namer, typr, indcr = exts[ind] if len(grps) == iter+1: @@ -3120,17 +3137,17 @@ def get_extension_edge(self, parent, template_rxn_map, obj, T, iter_max=np.inf, if first_time: first_time = False - + if grps[iter] == [] and len(grps) != iter+1 and (not (any([len(x)>0 for x in out_exts]) and iter+1 > iter_max)): iter += 1 if len(grps[iter]) > iter_item_cap: logging.error("Recursion item cap hit not splitting {0} reactions at iter {1} with {2} items".format(len(template_rxn_map[parent.label]),iter,len(grps[iter]))) iter -= 1 gave_up_split = True - + elif grps[iter] == [] and len(grps) != iter+1 and (any([len(x)>0 for x in out_exts]) and iter+1 > iter_max): logging.error("iter_max achieved terminating early") - + out = [] # compile all of the valid extensions together # may be some duplicates here, but I don't think it's currently worth identifying them @@ -3141,7 +3158,7 @@ def get_extension_edge(self, parent, template_rxn_map, obj, T, iter_max=np.inf, def extend_node(self, parent, template_rxn_map, obj=None, T=1000.0, iter_max=np.inf, iter_item_cap=np.inf): """ - Constructs an extension to the group parent based on evaluation + Constructs an extension to the group parent based on evaluation of the objective function obj """ exts, gave_up_split = self.get_extension_edge(parent, template_rxn_map, obj=obj, T=T, iter_max=iter_max, iter_item_cap=iter_item_cap) @@ -3197,10 +3214,10 @@ def extend_node(self, parent, template_rxn_map, obj=None, T=1000.0, iter_max=np. parent.item.clear_reg_dims() # this almost always solves the problem return True return False - + if gave_up_split: return False - + vals = [] for grp, grpc, name, typ, einds in exts: val, boo = self.eval_ext(parent, grp, name, template_rxn_map, obj, T) @@ -3296,7 +3313,7 @@ def extend_node(self, parent, template_rxn_map, obj=None, T=1000.0, iter_max=np. logging.error(prod.label) logging.error(prod.to_adjacency_list()) raise ValueError - + template_rxn_map[extname] = new_entries if complement: @@ -3312,21 +3329,21 @@ def generate_tree(self, rxns=None, obj=None, thermo_database=None, T=1000.0, npr """ Generate a tree by greedy optimization based on the objective function obj the optimization is done by iterating through every group and if the group has - more than one training reaction associated with it a set of potential more specific extensions - are generated and the extension that optimizing the objective function combination is chosen + more than one training reaction associated with it a set of potential more specific extensions + are generated and the extension that optimizing the objective function combination is chosen and the iteration starts over at the beginning - + additionally the tree structure is simplified on the fly by removing groups that have no kinetics data associated if their parent has no kinetics data associated and they either have only one child or have two children one of which has no kinetics data and no children (its parent becomes the parent of its only relevant child node) - + Args: rxns: List of reactions to generate tree from (if None pull the whole training set) obj: Object to expand tree from (if None uses top node) thermo_database: Thermodynamic database used for reversing training reactions T: Temperature the tree is optimized for - nprocs: Number of process for parallel tree generation + nprocs: Number of process for parallel tree generation min_splitable_entry_num: the minimum number of splitable reactions at a node in order to spawn a new process solving that node min_rxns_to_spawn: the minimum number of reactions at a node to spawn a new process solving that node @@ -3370,9 +3387,9 @@ def rxnkey(rxn): min_splitable_entry_num=min_splitable_entry_num, min_rxns_to_spawn=min_rxns_to_spawn, extension_iter_max=extension_iter_max, extension_iter_item_cap=extension_iter_item_cap) logging.error("built tree with {} nodes".format(len(list(self.groups.entries)))) - + self.auto_generated = True - + def get_rxn_batches(self, rxns, T=1000.0, max_batch_size=800, outlier_fraction=0.02, stratum_num=8): """ Breaks reactions into batches based on a modified stratified sampling scheme @@ -3475,7 +3492,7 @@ def make_tree_nodes(self, template_rxn_map=None, obj=None, T=1000.0, nprocs=0, d entries.remove(entry) else: psize = float(len(template_rxn_map[root.label])) - + logging.error(psize) mult_completed_nodes = [] # nodes containing multiple identical training reactions boo = True # if the for loop doesn't break becomes false and the while loop terminates @@ -3601,9 +3618,12 @@ def make_bm_rules_from_template_rxn_map(self, template_rxn_map, nprocs=1, Tref=1 inds = inds.tolist() revinds = [inds.index(x) for x in np.arange(len(inputs))] - pool = mp.Pool(nprocs) + if nprocs > 1: + pool = mp.Pool(nprocs) + kinetics_list = np.array(pool.map(_make_rule, inputs[inds])) + else: + kinetics_list = np.array(list(map(_make_rule, inputs[inds]))) - kinetics_list = np.array(pool.map(_make_rule, inputs[inds])) kinetics_list = kinetics_list[revinds] # fix order for i, kinetics in enumerate(kinetics_list): @@ -3628,11 +3648,11 @@ def make_bm_rules_from_template_rxn_map(self, template_rxn_map, nprocs=1, Tref=1 index += 1 - def cross_validate(self, folds=5, template_rxn_map=None, test_rxn_inds=None, T=1000.0, iters=0, random_state=1, ascend=False): + def cross_validate(self, folds=5, template_rxn_map=None, test_rxn_inds=None, T=1000.0, iters=0, random_state=1): """ Perform K-fold cross validation on an automatically generated tree at temperature T after finding an appropriate node for kinetics estimation it will move up the tree - iters times. + iters times. Returns a dictionary mapping {rxn:Ln(k_Est/k_Train)} """ @@ -3684,44 +3704,44 @@ def cross_validate(self, folds=5, template_rxn_map=None, test_rxn_inds=None, T=1 if entry.parent: entry = entry.parent + boo = True + + while boo: + if entry.parent is None: + break + kin = self.rules.entries[entry.label][0].data + kinparent = self.rules.entries[entry.parent.label][0].data + err_parent = abs(kinparent.uncertainty.data_mean + kinparent.uncertainty.mu - kin.uncertainty.data_mean) + np.sqrt(2.0*kinparent.uncertainty.var/np.pi) + err_entry = abs(kin.uncertainty.mu) + np.sqrt(2.0*kin.uncertainty.var/np.pi) + if err_entry <= err_parent: + break + else: + entry = entry.parent + uncertainties[rxn] = self.rules.entries[entry.label][0].data.uncertainty - - if not ascend: - L = list(set(template_rxn_map[entry.label]) - set(rxns_test)) - if L != []: + + L = list(set(template_rxn_map[entry.label]) - set(rxns_test)) + + if L != []: + if isinstance(L[0].kinetics, Arrhenius): kinetics = ArrheniusBM().fit_to_reactions(L, recipe=self.forward_recipe.actions) - kinetics = kinetics.to_arrhenius(rxn.get_enthalpy_of_reaction(T)) - k = kinetics.get_rate_coefficient(T) - errors[rxn] = np.log(k / krxn) + if kinetics.E0.value_si < 0.0 or len(L) == 1: + kinetics = average_kinetics([r.kinetics for r in L]) + else: + kinetics = kinetics.to_arrhenius(rxn.get_enthalpy_of_reaction(298.0)) else: - raise ValueError('only one piece of kinetics information in the tree?') - else: - boo = True - rlist = list(set(template_rxn_map[entry.label]) - set(rxns_test)) - kinetics = _make_rule((self.forward_recipe.actions,rlist,T,1.0e3,"",[rxn.rank for rxn in rlist])) - logging.error("determining fold rate") - c = 1 - while boo: - parent = entry.parent - if parent is None: - break - rlistparent = list(set(template_rxn_map[parent.label]) - set(rxns_test)) - kineticsparent = _make_rule((self.forward_recipe.actions,rlistparent,T,1.0e3,"",[rxn.rank for rxn in rlistparent])) - err_parent = abs(kineticsparent.uncertainty.data_mean + kineticsparent.uncertainty.mu - kinetics.uncertainty.data_mean) + np.sqrt(2.0*kineticsparent.uncertainty.var/np.pi) - err_entry = abs(kinetics.uncertainty.mu) + np.sqrt(2.0*kinetics.uncertainty.var/np.pi) - if err_entry > err_parent: - entry = entry.parent - kinetics = kineticsparent - logging.error("recursing {}".format(c)) - c += 1 + kinetics = ArrheniusChargeTransferBM().fit_to_reactions(L, recipe=self.forward_recipe.actions) + if kinetics.E0.value_si < 0.0 or len(L) == 1: + kinetics = average_kinetics([r.kinetics for r in L]) else: - boo = False - - kinetics = kinetics.to_arrhenius(rxn.get_enthalpy_of_reaction(T)) + kinetics = kinetics.to_arrhenius_charge_transfer(rxn.get_enthalpy_of_reaction(298.0)) + k = kinetics.get_rate_coefficient(T) errors[rxn] = np.log(k / krxn) - + else: + raise ValueError('only one piece of kinetics information in the tree?') + return errors, uncertainties def cross_validate_old(self, folds=5, T=1000.0, random_state=1, estimator='rate rules', thermo_database=None, get_reverse=False, uncertainties=True): @@ -3731,7 +3751,7 @@ def cross_validate_old(self, folds=5, T=1000.0, random_state=1, estimator='rate """ errors = {} uncs = {} - + kpu = KineticParameterUncertainty() rxns = np.array(self.get_training_set(remove_degeneracy=True,get_reverse=get_reverse)) @@ -3777,7 +3797,7 @@ def cross_validate_old(self, folds=5, T=1000.0, random_state=1, estimator='rate boo,source = self.extract_source_from_comments(testrxn) sdict = {"Rate Rules":source} uncs[rxn] = kpu.get_uncertainty_value(sdict) - + if uncertainties: return errors, uncs else: @@ -3787,19 +3807,19 @@ def simple_regularization(self, node, template_rxn_map, test=True): """ Simplest regularization algorithm All nodes are made as specific as their descendant reactions - Training reactions are assumed to not generalize + Training reactions are assumed to not generalize For example if an particular atom at a node is Oxygen for all of its descendent reactions a reaction where it is Sulfur will never hit that node - unless it is the top node even if the tree did not split on the identity + unless it is the top node even if the tree did not split on the identity of that atom - - The test option to this function determines whether or not the reactions - under a node match the extended group before adding an extension. - If the test fails the extension is skipped. - - In general test=True is needed if the cascade algorithm was used + + The test option to this function determines whether or not the reactions + under a node match the extended group before adding an extension. + If the test fails the extension is skipped. + + In general test=True is needed if the cascade algorithm was used to generate the tree and test=False is ok if the cascade algorithm - wasn't used. + wasn't used. """ for child in node.children: @@ -4005,7 +4025,7 @@ def clean_tree(self): def save_generated_tree(self, path=None): """ - clears the rules and saves the family to its + clears the rules and saves the family to its current location in database """ if path is None: @@ -4019,7 +4039,7 @@ def get_training_set(self, thermo_database=None, remove_degeneracy=False, estima """ retrieves all reactions in the training set, assigns thermo to the species objects reverses reactions as necessary so that all reactions are in the forward direction - and returns the resulting list of reactions in the forward direction with thermo + and returns the resulting list of reactions in the forward direction with thermo assigned """ @@ -4089,7 +4109,7 @@ def get_reactant_thermo(reactant,metal): root_labels = [x.label for x in root.atoms if x.label != ''] root_label_set = set(root_labels) - + for i, entry in enumerate(entries): if estimate_thermo: # parse out the metal to scale to @@ -4119,6 +4139,8 @@ def get_reactant_thermo(reactant,metal): else: mol = deepcopy(react.molecule[0]) + mol.update_atomtypes() + if fix_labels: for prod in rxns[i].products: fix_labels_mol(prod.molecule[0], root_labels) @@ -4137,12 +4159,12 @@ def get_reactant_thermo(reactant,metal): mol = mol.merge(react.molecule[0]) else: mol = deepcopy(react.molecule[0]) - + if fix_labels: mol_label_set = set([x.label for x in get_label_fixed_mol(mol, root_labels).atoms if x.label != '']) else: mol_label_set = set([x.label for x in mol.atoms if x.label != '']) - + if mol_label_set == root_label_set and ((mol.is_subgraph_isomorphic(root, generate_initial_map=True) or (not fix_labels and get_label_fixed_mol(mol, root_labels).is_subgraph_isomorphic(root, generate_initial_map=True)))): @@ -4204,6 +4226,8 @@ def get_reactant_thermo(reactant,metal): else: mol = deepcopy(react.molecule[0]) + mol.update_atomtypes() + if (mol.is_subgraph_isomorphic(root, generate_initial_map=True) or (not fix_labels and get_label_fixed_mol(mol, root_labels).is_subgraph_isomorphic(root, generate_initial_map=True))): # try product structures @@ -4232,7 +4256,7 @@ def get_reactant_thermo(reactant,metal): def get_reaction_matches(self, rxns=None, thermo_database=None, remove_degeneracy=False, estimate_thermo=True, fix_labels=False, exact_matches_only=False, get_reverse=False): """ - returns a dictionary mapping for each entry in the tree: + returns a dictionary mapping for each entry in the tree: (entry.label,entry.item) : list of all training reactions (or the list given) that match that entry """ if rxns is None: @@ -4325,10 +4349,10 @@ def retrieve_original_entry(self, template_label): """ Retrieves the original entry, be it a rule or training reaction, given the template label in the form 'group1;group2' or 'group1;group2;group3' - + Returns tuple in the form (RateRuleEntry, TrainingReactionEntry) - + Where the TrainingReactionEntry is only present if it comes from a training reaction """ template_labels = template_label.split()[-1].split(';') @@ -4345,13 +4369,13 @@ def get_sources_for_template(self, template): """ Returns the set of rate rules and training reactions used to average this `template`. Note that the tree must be averaged with verbose=True for this to work. - + Returns a tuple of rules, training - - where rules are a list of tuples containing + + where rules are a list of tuples containing the [(original_entry, weight_used_in_average), ... ] - + and training is a list of tuples containing the [(rate_rule_entry, training_reaction_entry, weight_used_in_average),...] """ @@ -4359,7 +4383,7 @@ def get_sources_for_template(self, template): def assign_weights_to_entries(entry_nested_list, weighted_entries, n=1): """ Assign weights to an average of average nested list. Where n is the - number of values being averaged recursively. + number of values being averaged recursively. """ n = len(entry_nested_list) * n for entry in entry_nested_list: @@ -4433,7 +4457,7 @@ def assign_weights_to_entries(entry_nested_list, weighted_entries, n=1): rules[rule_entry] += weight else: rules[rule_entry] = weight - # Each entry should now only appear once + # Each entry should now only appear once training = [(k[0], k[1], v) for k, v in training.items()] rules = list(rules.items()) @@ -4444,11 +4468,11 @@ def extract_source_from_comments(self, reaction): Returns the rate rule associated with the kinetics of a reaction by parsing the comments. Will return the template associated with the matched rate rule. Returns a tuple containing (Boolean_Is_Kinetics_From_Training_reaction, Source_Data) - + For a training reaction, the Source_Data returns:: [Family_Label, Training_Reaction_Entry, Kinetics_In_Reverse?] - + For a reaction from rate rules, the Source_Data is a tuple containing:: [Family_Label, {'template': originalTemplate, @@ -4483,7 +4507,7 @@ def extract_source_from_comments(self, reaction): 'but does not match the training reaction {1} from the ' '{2} family.'.format(reaction, training_reaction_index, self.label)) - # Sometimes the matched kinetics could be in the reverse direction..... + # Sometimes the matched kinetics could be in the reverse direction..... if reaction.is_isomorphic(training_entry.item, either_direction=False, save_order=self.save_order): reverse = False else: @@ -4497,7 +4521,7 @@ def extract_source_from_comments(self, reaction): elif line.startswith('Multiplied by'): degeneracy = float(line.split()[-1]) - # Extract the rate rule information + # Extract the rate rule information full_comment_string = reaction.kinetics.comment.replace('\n', ' ') # The rate rule string is right after the phrase 'for rate rule' @@ -4599,25 +4623,61 @@ def _make_rule(rr): rxns = np.array(rxns) data_mean = np.mean(np.log([r.kinetics.get_rate_coefficient(Tref) for r in rxns])) if n > 0: - kin = ArrheniusBM().fit_to_reactions(rxns, recipe=recipe) - if n == 1: - kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label) + if isinstance(rxns[0].kinetics, Arrhenius): + arr = ArrheniusBM else: - dlnks = np.array([ - np.log( - ArrheniusBM().fit_to_reactions(rxns[list(set(range(len(rxns))) - {i})], recipe=recipe) + arr = ArrheniusChargeTransferBM + if n > 1: + kin = arr().fit_to_reactions(rxns, recipe=recipe) + if n == 1 or kin.E0.value_si < 0.0: + kin = average_kinetics([r.kinetics for r in rxns]) + #kin.comment = "Only one reaction or Arrhenius BM fit bad. Instead averaged from {} reactions.".format(n) + if n == 1: + kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label) + else: + dlnks = np.array([ + np.log( + average_kinetics([r.kinetics for r in rxns[list(set(range(len(rxns))) - {i})]]).get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref) + ) for i, rxn in enumerate(rxns) + ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref + varis = (np.array([rank_accuracy_map[rxn.rank].value_si for rxn in rxns]) / (2.0 * 8.314 * Tref)) ** 2 + # weighted average calculations + ws = 1.0 / varis + V1 = ws.sum() + V2 = (ws ** 2).sum() + mu = np.dot(ws, dlnks) / V1 + s = np.sqrt(np.dot(ws, (dlnks - mu) ** 2) / (V1 - V2 / V1)) + kin.uncertainty = RateUncertainty(mu=mu, var=s ** 2, N=n, Tref=Tref, data_mean=data_mean, correlation=label) + return kin + else: + if n == 1: + kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label) + else: + if isinstance(rxns[0].kinetics, Arrhenius): + dlnks = np.array([ + np.log( + arr().fit_to_reactions(rxns[list(set(range(len(rxns))) - {i})], recipe=recipe) .to_arrhenius(rxn.get_enthalpy_of_reaction(Tref)) .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref) ) for i, rxn in enumerate(rxns) ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref - varis = (np.array([rank_accuracy_map[rxn.rank].value_si for rxn in rxns]) / (2.0 * 8.314 * Tref)) ** 2 - # weighted average calculations - ws = 1.0 / varis - V1 = ws.sum() - V2 = (ws ** 2).sum() - mu = np.dot(ws, dlnks) / V1 - s = np.sqrt(np.dot(ws, (dlnks - mu) ** 2) / (V1 - V2 / V1)) - kin.uncertainty = RateUncertainty(mu=mu, var=s ** 2, N=n, Tref=Tref, data_mean=data_mean, correlation=label) + else: + dlnks = np.array([ + np.log( + arr().fit_to_reactions(rxns[list(set(range(len(rxns))) - {i})], recipe=recipe) + .to_arrhenius_charge_transfer(rxn.get_enthalpy_of_reaction(Tref)) + .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref) + ) for i, rxn in enumerate(rxns) + ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref + varis = (np.array([rank_accuracy_map[rxn.rank].value_si for rxn in rxns]) / (2.0 * 8.314 * Tref)) ** 2 + # weighted average calculations + + ws = 1.0 / varis + V1 = ws.sum() + V2 = (ws ** 2).sum() + mu = np.dot(ws, dlnks) / V1 + s = np.sqrt(np.dot(ws, (dlnks - mu) ** 2) / (V1 - V2 / V1)) + kin.uncertainty = RateUncertainty(mu=mu, var=s ** 2, N=n, Tref=Tref, data_mean=data_mean, correlation=label) return kin else: return None @@ -4646,7 +4706,87 @@ def _child_make_tree_nodes(family, child_conn, template_rxn_map, obj, T, nprocs, family.groups.entries[root_label].parent = None family.make_tree_nodes(template_rxn_map=template_rxn_map, obj=obj, T=T, nprocs=nprocs, depth=depth + 1, - min_splitable_entry_num=min_splitable_entry_num, min_rxns_to_spawn=min_rxns_to_spawn, + min_splitable_entry_num=min_splitable_entry_num, min_rxns_to_spawn=min_rxns_to_spawn, extension_iter_max=extension_iter_max, extension_iter_item_cap=extension_iter_item_cap) child_conn.send(list(family.groups.entries.values())) + +def average_kinetics(kinetics_list): + """ + Based on averaging log k. + Hence we average n, Ea, arithmetically, but we + average log A (geometric average) + """ + logA = 0.0 + n = 0.0 + Ea = 0.0 + alpha = 0.5 + electrons = None + if isinstance(kinetics_list[0], SurfaceChargeTransfer) or isinstance(kinetics_list[0], ArrheniusChargeTransfer): + if electrons is None: + electrons = kinetics_list[0].electrons.value_si + assert all(np.abs(k.V0.value_si) < 0.0001 for k in kinetics_list), [k.V0.value_si for k in kinetics_list] + assert all(np.abs(k.alpha.value_si - 0.5) < 0.001 for k in kinetics_list), [k.alpha for k in kinetics_list] + V0 = 0.0 + count = 0 + for kinetics in kinetics_list: + count += 1 + logA += np.log10(kinetics.A.value_si) + n += kinetics.n.value_si + Ea += kinetics.Ea.value_si + + logA /= count + n /= count + alpha /= count + Ea /= count + Aunits = kinetics_list[0].A.units + if Aunits == 'cm^3/(mol*s)' or Aunits == 'cm^3/(molecule*s)' or Aunits == 'm^3/(molecule*s)': + Aunits = 'm^3/(mol*s)' + elif Aunits == 'cm^6/(mol^2*s)' or Aunits == 'cm^6/(molecule^2*s)' or Aunits == 'm^6/(molecule^2*s)': + Aunits = 'm^6/(mol^2*s)' + elif Aunits == 's^-1' or Aunits == 'm^3/(mol*s)' or Aunits == 'm^6/(mol^2*s)': + # they were already in SI + pass + elif Aunits in ['m^2/(mol*s)', 'cm^2/(mol*s)', 'm^2/(molecule*s)', 'cm^2/(molecule*s)']: + # surface: bimolecular (Langmuir-Hinshelwood) + Aunits = 'm^2/(mol*s)' + elif Aunits in ['m^5/(mol^2*s)', 'cm^5/(mol^2*s)', 'm^5/(molecule^2*s)', 'cm^5/(molecule^2*s)']: + # surface: dissociative adsorption + Aunits = 'm^5/(mol^2*s)' + elif Aunits == '': + # surface: sticking coefficient + pass + else: + raise Exception('Invalid units {0} for averaging kinetics.'.format(Aunits)) + + if type(kinetics) not in [Arrhenius,SurfaceChargeTransfer,ArrheniusChargeTransfer]: + raise Exception('Invalid kinetics type {0!r} for {1!r}.'.format(type(kinetics), self)) + + if isinstance(kinetics, SurfaceChargeTransfer): + averaged_kinetics = SurfaceChargeTransfer( + A=(10 ** logA, Aunits), + n=n, + electrons=electrons, + alpha=alpha, + V0=(V0,'V'), + Ea=(Ea * 0.001, "kJ/mol"), + comment="Only one reaction or BM with E0 < 0. Averaged from {} reactions.".format(len(kinetics_list)), + ) + elif isinstance(kinetics, ArrheniusChargeTransfer): + averaged_kinetics = ArrheniusChargeTransfer( + A=(10 ** logA, Aunits), + n=n, + electrons=electrons, + alpha=alpha, + V0=(V0,'V'), + Ea=(Ea * 0.001, "kJ/mol"), + comment="Only one reaction or BM with E0 < 0. Averaged from {} reactions.".format(len(kinetics_list)), + ) + else: + averaged_kinetics = Arrhenius( + A=(10 ** logA, Aunits), + n=n, + Ea=(Ea * 0.001, "kJ/mol"), + comment="Only one reaction or BM with E0 < 0. Averaged from {} reactions.".format(len(kinetics_list)), + ) + return averaged_kinetics diff --git a/rmgpy/kinetics/__init__.py b/rmgpy/kinetics/__init__.py index a9dc32e49b..7acf5ceb53 100644 --- a/rmgpy/kinetics/__init__.py +++ b/rmgpy/kinetics/__init__.py @@ -29,7 +29,8 @@ from rmgpy.kinetics.model import KineticsModel, PDepKineticsModel, TunnelingModel, \ get_rate_coefficient_units_from_reaction_order, get_reaction_order_from_rate_coefficient_units -from rmgpy.kinetics.arrhenius import Arrhenius, ArrheniusEP, PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, ArrheniusBM +from rmgpy.kinetics.arrhenius import Arrhenius, ArrheniusEP, PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, \ + ArrheniusBM, ArrheniusChargeTransfer, ArrheniusChargeTransferBM from rmgpy.kinetics.chebyshev import Chebyshev from rmgpy.kinetics.falloff import ThirdBody, Lindemann, Troe from rmgpy.kinetics.kineticsdata import KineticsData, PDepKineticsData diff --git a/rmgpy/kinetics/arrhenius.pxd b/rmgpy/kinetics/arrhenius.pxd index d3f2b78dc3..fda28a80f5 100644 --- a/rmgpy/kinetics/arrhenius.pxd +++ b/rmgpy/kinetics/arrhenius.pxd @@ -128,4 +128,55 @@ cdef class MultiPDepArrhenius(PDepKineticsModel): cpdef bint is_identical_to(self, KineticsModel other_kinetics) except -2 + + cpdef change_rate(self, double factor) + +################################################################################ +cdef class ArrheniusChargeTransfer(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 + + +################################################################################ + +cdef class ArrheniusChargeTransferBM(KineticsModel): + + cdef public ScalarQuantity _A + cdef public ScalarQuantity _n + cdef public ScalarQuantity _E0 + cdef public ScalarQuantity _w0 + 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_rate_coefficient_from_potential(self, double T, double V, double dGrxn) except -1 + + cpdef ArrheniusChargeTransfer to_arrhenius_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/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index a02806dfb0..db248462bd 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -37,6 +37,7 @@ import rmgpy.quantity as quantity from rmgpy.exceptions import KineticsError from rmgpy.kinetics.uncertainties import rank_accuracy_map from rmgpy.molecule.molecule import Bond +import logging # Prior to numpy 1.14, `numpy.linalg.lstsq` does not accept None as a value RCOND = -1 if int(np.__version__.split('.')[1]) < 14 else None @@ -58,15 +59,16 @@ cdef class Arrhenius(KineticsModel): `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 + `solute` Transition state solute data `comment` Information about the model (e.g. its source) =============== ============================================================= """ def __init__(self, A=None, n=0.0, Ea=None, T0=(1.0, "K"), Tmin=None, Tmax=None, Pmin=None, Pmax=None, - uncertainty=None, comment=''): + uncertainty=None, solute=None, comment=''): KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, uncertainty=uncertainty, - comment=comment) + solute=solute, comment=comment) self.A = A self.n = n self.Ea = Ea @@ -83,6 +85,7 @@ cdef class Arrhenius(KineticsModel): 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.solute: string += ', solute={0!r}'.format(self.solute) if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) string += ')' return string @@ -92,7 +95,7 @@ cdef class Arrhenius(KineticsModel): A helper function used when pickling an Arrhenius object. """ return (Arrhenius, (self.A, self.n, self.Ea, self.T0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, - self.uncertainty, self.comment)) + self.uncertainty, self.solute, self.comment)) property A: """The preexponential factor.""" @@ -196,6 +199,7 @@ cdef class Arrhenius(KineticsModel): self.T0 = (T0, "K") self.Tmin = (np.min(Tlist), "K") self.Tmax = (np.max(Tlist), "K") + self.solute = None 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])), @@ -295,6 +299,7 @@ cdef class Arrhenius(KineticsModel): Pmin=self.Pmin, Pmax=self.Pmax, uncertainty=self.uncertainty, + solute=self.solute, comment=self.comment) return aep ################################################################################ @@ -317,15 +322,16 @@ cdef class ArrheniusEP(KineticsModel): `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 + `solute` Transition state solute data `comment` Information about the model (e.g. its source) =============== ============================================================= """ def __init__(self, A=None, n=0.0, alpha=0.0, E0=None, Tmin=None, Tmax=None, Pmin=None, Pmax=None, uncertainty=None, - comment=''): + solute=None, comment=''): KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, uncertainty=uncertainty, - comment=comment) + solute=solute, comment=comment) self.A = A self.n = n self.alpha = alpha @@ -342,6 +348,7 @@ cdef class ArrheniusEP(KineticsModel): 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 is not None: string += ', uncertainty={0!r}'.format(self.uncertainty) + if self.solute is not None: string += ', solute={0!r}'.format(self.solute) if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) string += ')' return string @@ -351,7 +358,7 @@ cdef class ArrheniusEP(KineticsModel): A helper function used when pickling an ArrheniusEP object. """ return (ArrheniusEP, (self.A, self.n, self.alpha, self.E0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, - self.uncertainty, self.comment)) + self.uncertainty, self.solute, self.comment)) property A: """The preexponential factor.""" @@ -422,6 +429,7 @@ cdef class ArrheniusEP(KineticsModel): Pmin=self.Pmin, Pmax=self.Pmax, uncertainty=self.uncertainty, + solute=self.solute, comment=self.comment, ) @@ -475,15 +483,16 @@ cdef class ArrheniusBM(KineticsModel): `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 + `solute` Transition state solute data `comment` Information about the model (e.g. its source) =============== ============================================================= """ def __init__(self, A=None, n=0.0, w0=(0.0, 'J/mol'), E0=None, Tmin=None, Tmax=None, Pmin=None, Pmax=None, - uncertainty=None, comment=''): + uncertainty=None, solute=None, comment=''): KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, uncertainty=uncertainty, - comment=comment) + solute=solute, comment=comment) self.A = A self.n = n self.w0 = w0 @@ -500,6 +509,7 @@ cdef class ArrheniusBM(KineticsModel): 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 is not None: string += ', uncertainty={0!r}'.format(self.uncertainty) + if self.solute is not None: string += ', solute={0!r}'.format(self.solute) if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) string += ')' return string @@ -509,7 +519,7 @@ cdef class ArrheniusBM(KineticsModel): A helper function used when pickling an ArrheniusEP object. """ return (ArrheniusBM, (self.A, self.n, self.w0, self.E0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, - self.uncertainty, self.comment)) + self.uncertainty, self.solute, self.comment)) property A: """The preexponential factor.""" @@ -580,6 +590,7 @@ cdef class ArrheniusBM(KineticsModel): Tmin=self.Tmin, Tmax=self.Tmax, uncertainty=self.uncertainty, + solute=self.solute, comment=self.comment, ) @@ -600,26 +611,24 @@ cdef class ArrheniusBM(KineticsModel): if len(rxns) == 1: T = 1000.0 rxn = rxns[0] - dHrxn = rxn.get_enthalpy_of_reaction(T) + dHrxn = rxn.get_enthalpy_of_reaction(298.0) A = rxn.kinetics.A.value_si n = rxn.kinetics.n.value_si Ea = rxn.kinetics.Ea.value_si - + def kfcn(E0): Vp = 2 * w0 * (2 * w0 + 2 * E0) / (2 * w0 - 2 * E0) out = Ea - (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (2 * w0) * (2 * w0) + dHrxn * dHrxn) return out - if abs(dHrxn) > 4 * w0 / 10.0: - E0 = w0 / 10.0 - else: - E0 = fsolve(kfcn, w0 / 10.0)[0] + E0 = fsolve(kfcn, w0 / 10.0)[0] self.Tmin = rxn.kinetics.Tmin self.Tmax = rxn.kinetics.Tmax + self.solute = None self.comment = 'Fitted to {0} reaction at temperature: {1} K'.format(len(rxns), T) else: - # define optimization function + # define optimization function def kfcn(xs, lnA, n, E0): T = xs[:,0] dHrxn = xs[:,1] @@ -628,7 +637,7 @@ cdef class ArrheniusBM(KineticsModel): Ea = np.where(dHrxn< -4.0*E0, 0.0, Ea) Ea = np.where(dHrxn > 4.0*E0, dHrxn, Ea) return lnA + np.log(T ** n * np.exp(-Ea / (8.314 * T))) - + # get (T,dHrxn(T)) -> (Ln(k) mappings xdata = [] ydata = [] @@ -637,7 +646,7 @@ cdef class ArrheniusBM(KineticsModel): # approximately correct the overall uncertainties to std deviations s = rank_accuracy_map[rxn.rank].value_si/2.0 for T in Ts: - xdata.append([T, rxn.get_enthalpy_of_reaction(T)]) + xdata.append([T, rxn.get_enthalpy_of_reaction(298.0)]) ydata.append(np.log(rxn.get_rate_coefficient(T))) sigmas.append(s / (8.314 * T)) @@ -666,6 +675,7 @@ cdef class ArrheniusBM(KineticsModel): self.Tmin = (np.min(Ts), "K") self.Tmax = (np.max(Ts), "K") + self.solute = None self.comment = 'Fitted to {0} reactions at temperatures: {1}'.format(len(rxns), Ts) # fill in parameters @@ -1090,6 +1100,581 @@ cdef class MultiPDepArrhenius(PDepKineticsModel): for i, arr in enumerate(self.arrhenius): arr.set_cantera_kinetics(ct_reaction[i], species_list) +################################################################################ + +cdef class ArrheniusChargeTransfer(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 + `solute` The transition state solute data + `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, solute=None, uncertainty=None, comment=''): + + KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, solute=solute, 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 = 'ArrheniusChargeTransfer(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.solute: string += ', solute={0!r}'.format(self.solute) + 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 ArrheniusChargeTransfer object. + """ + return (ArrheniusChargeTransfer, (self.A, self.n, self.Ea, self.V0, self.alpha, self.electrons, self.T0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, + self.solute, 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.solute = None, + 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, ArrheniusChargeTransfer): + 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 + +cdef class ArrheniusChargeTransferBM(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 + `w0` The average of the bond dissociation energies of the bond formed and the bond broken + `E0` The activation energy for a thermoneutral reaction + `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 + `solute` Transition state solute data + `comment` Information about the model (e.g. its source) + =============== ============================================================= + + """ + + def __init__(self, A=None, n=0.0, w0=(0.0, 'J/mol'), E0=None, V0=(0.0,'V'), alpha=0.5, electrons=-1, Tmin=None, Tmax=None, + Pmin=None, Pmax=None, solute=None, uncertainty=None, comment=''): + + KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, solute=solute, uncertainty=uncertainty, + comment=comment) + + self.alpha = alpha + self.A = A + self.n = n + self.w0 = w0 + 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 = 'ArrheniusChargeTransferBM(A={0!r}, n={1!r}, w0={2!r}, E0={3!r}, V0={4!r}, alpha={5!r}, electrons={6!r}'.format( + self.A, self.n, self.w0, 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.solute: string += ', solute={0!r}'.format(self.solute) + 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 ArrheniusChargeTransfer object. + """ + return (ArrheniusChargeTransferBM, (self.A, self.n, self.w0, self.E0, self.V0, self.alpha, self.electrons, self.Tmin, self.Tmax, self.Pmin, self.Pmax, + self.solute, 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 w0: + """The average of the bond dissociation energies of the bond formed and the bond broken.""" + def __get__(self): + return self._w0 + def __set__(self, value): + self._w0 = quantity.Energy(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_rate_coefficient(self, double T, double dHrxn=0.0) except -1: + """ + Return the rate coefficient in the appropriate combination of m^3, + mol, and s at temperature `T` in K and enthalpy of reaction `dHrxn` + in J/mol. + """ + cdef double A, n, Ea + Ea = self.get_activation_energy(dHrxn) + A = self._A.value_si + n = self._n.value_si + return A * T ** n * exp(-Ea / (constants.R * T)) + + cpdef double get_activation_energy(self, double dHrxn) except -1: + """ + Return the activation energy in J/mol corresponding to the given + enthalpy of reaction `dHrxn` in J/mol. + """ + cdef double w0, E0 + E0 = self._E0.value_si + if dHrxn < -4 * self._E0.value_si: + return 0.0 + elif dHrxn > 4 * self._E0.value_si: + return dHrxn + else: + w0 = self._w0.value_si + Vp = 2 * w0 * (2 * w0 + 2 * E0) / (2 * w0 - 2 * E0) + return (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) ** 2 / (Vp ** 2 - (2 * w0) ** 2 + dHrxn ** 2) + + cpdef double get_rate_coefficient_from_potential(self, double T, double V, double dHrxn) 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 + heat of reaction `dHrxn` in J/mol. + """ + cdef double A, n, Ea + Ea = self.get_activation_energy_from_potential(V,dHrxn) + Ea -= self._alpha.value_si * self._electrons.value_si * constants.F * (V-self._V0.value_si) + A = self._A.value_si + n = self._n.value_si + return A * T ** n * exp(-Ea / (constants.R * T)) + + def fit_to_reactions(self, rxns, w0=None, recipe=None, Ts=None): + """ + Fit an ArrheniusBM model to a list of reactions at the given temperatures, + w0 must be either given or estimated using the family object + """ + assert w0 is not None or recipe is not None, 'either w0 or recipe must be specified' + + for rxn in rxns: + if rxn.kinetics._V0.value_si != 0.0: + rxn.kinetics.change_v0(0.0) + + if Ts is None: + Ts = [300.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1500.0] + if w0 is None: + #estimate w0 + w0s = get_w0s(recipe, rxns) + w0 = sum(w0s) / len(w0s) + + if len(rxns) == 1: + rxn = rxns[0] + dHrxn = rxn.get_enthalpy_of_reaction(298.0) + A = rxn.kinetics.A.value_si + n = rxn.kinetics.n.value_si + Ea = rxn.kinetics.Ea.value_si + + def kfcn(E0): + Vp = 2 * w0 * (2 * w0 + 2 * E0) / (2 * w0 - 2 * E0) + out = Ea - (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (2 * w0) * (2 * w0) + dHrxn * dHrxn) + return out + + E0 = fsolve(kfcn, w0 / 10.0)[0] + + self.Tmin = rxn.kinetics.Tmin + self.Tmax = rxn.kinetics.Tmax + self.comment = 'Fitted to 1 reaction' + else: + # define optimization function + def kfcn(xs, lnA, n, E0): + T = xs[:,0] + dHrxn = xs[:,1] + Vp = 2 * w0 * (2 * w0 + 2 * E0) / (2 * w0 - 2 * E0) + Ea = (w0 + dHrxn / 2.0) * (Vp - 2 * w0 + dHrxn) * (Vp - 2 * w0 + dHrxn) / (Vp * Vp - (2 * w0) * (2 * w0) + dHrxn * dHrxn) + Ea = np.where(dHrxn< -4.0*E0, 0.0, Ea) + Ea = np.where(dHrxn > 4.0*E0, dHrxn, Ea) + return lnA + np.log(T ** n * np.exp(-Ea / (8.314 * T))) + + # get (T,dHrxn(T)) -> (Ln(k) mappings + xdata = [] + ydata = [] + sigmas = [] + for rxn in rxns: + # approximately correct the overall uncertainties to std deviations + s = rank_accuracy_map[rxn.rank].value_si/2.0 + for T in Ts: + xdata.append([T, rxn.get_enthalpy_of_reaction(298.0)]) + ydata.append(np.log(rxn.get_rate_coefficient(T))) + + sigmas.append(s / (8.314 * T)) + + xdata = np.array(xdata) + ydata = np.array(ydata) + + # fit parameters + boo = True + xtol = 1e-8 + ftol = 1e-8 + while boo: + boo = False + try: + params = curve_fit(kfcn, xdata, ydata, sigma=sigmas, p0=[1.0, 1.0, w0 / 10.0], xtol=xtol, ftol=ftol) + except RuntimeError: + if xtol < 1.0: + boo = True + xtol *= 10.0 + ftol *= 10.0 + else: + raise ValueError("Could not fit BM arrhenius to reactions with xtol<1.0") + + lnA, n, E0 = params[0].tolist() + A = np.exp(lnA) + + self.Tmin = (np.min(Ts), "K") + self.Tmax = (np.max(Ts), "K") + self.comment = 'Fitted to {0} reactions at temperatures: {1}'.format(len(rxns), Ts) + + # fill in parameters + A_units = ['', 's^-1', 'm^3/(mol*s)', 'm^6/(mol^2*s)'] + order = len(rxns[0].reactants) + self.A = (A, A_units[order]) + + self.n = n + self.w0 = (w0, 'J/mol') + self.E0 = (E0, 'J/mol') + self._V0.value_si = 0.0 + self.electrons = rxns[0].electrons + + return self + + cpdef ArrheniusChargeTransfer to_arrhenius_charge_transfer(self, double dHrxn): + """ + Return an :class:`ArrheniusChargeTransfer` instance of the kinetics model using the + given heat of reaction `dHrxn` to determine the activation energy. + """ + return ArrheniusChargeTransfer( + A=self.A, + n=self.n, + electrons=self.electrons, + Ea=(self.get_activation_energy(dHrxn) * 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, + solute=self.solute, + 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, ArrheniusChargeTransferBM): + 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.w0.equals(other_kinetics.w0) + 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.') + def get_w0(actions, rxn): """ calculates the w0 for Blower Masel kinetics by calculating wf (total bond energy of bonds formed) diff --git a/rmgpy/kinetics/model.pxd b/rmgpy/kinetics/model.pxd index 2dcc0d0609..c640b4d94d 100644 --- a/rmgpy/kinetics/model.pxd +++ b/rmgpy/kinetics/model.pxd @@ -29,7 +29,7 @@ cimport numpy as np from rmgpy.quantity cimport ScalarQuantity, ArrayQuantity from rmgpy.kinetics.uncertainties cimport RateUncertainty - +from rmgpy.data.solvation import SoluteData ################################################################################ cpdef str get_rate_coefficient_units_from_reaction_order(n_gas, n_surf=?) @@ -43,6 +43,7 @@ cdef class KineticsModel: cdef public ScalarQuantity _Tmin, _Tmax cdef public ScalarQuantity _Pmin, _Pmax cdef public RateUncertainty uncertainty + cdef public object solute cdef public str comment diff --git a/rmgpy/kinetics/model.pyx b/rmgpy/kinetics/model.pyx index 5f38d38c99..f4a99f4cb0 100644 --- a/rmgpy/kinetics/model.pyx +++ b/rmgpy/kinetics/model.pyx @@ -117,16 +117,18 @@ cdef class KineticsModel: `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 + `solute` Solute data for the transition state `comment` Information about the model (e.g. its source) =============== ============================================================ """ - def __init__(self, Tmin=None, Tmax=None, Pmin=None, Pmax=None, uncertainty=None, comment=''): + def __init__(self, Tmin=None, Tmax=None, Pmin=None, Pmax=None, uncertainty=None, solute=None, comment=''): self.Tmin = Tmin self.Tmax = Tmax self.Pmin = Pmin self.Pmax = Pmax + self.solute = solute self.uncertainty = uncertainty self.comment = comment @@ -136,13 +138,13 @@ cdef class KineticsModel: KineticsModel object. """ return 'KineticsModel(Tmin={0!r}, Tmax={1!r}, Pmin={2!r}, Pmax={3!r}, uncertainty={4!r}, comment="""{5}""")'.format( - self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.uncertainty, self.comment) + self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.solute, self.uncertainty, self.comment) def __reduce__(self): """ A helper function used when pickling a KineticsModel object. """ - return (KineticsModel, (self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.uncertainty, self.comment)) + return (KineticsModel, (self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.solute, self.uncertainty, self.comment)) property Tmin: """The minimum temperature at which the model is valid, or ``None`` if not defined.""" diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index dfb598784f..f721c6d903 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -30,8 +30,8 @@ """ This module contains classes and functions for working with chemical reactions. -From the `IUPAC Compendium of Chemical Terminology -`_, a chemical reaction is "a process that +From the `IUPAC Compendium of Chemical Terminology +`_, a chemical reaction is "a process that results in the interconversion of chemical species". In RMG Py, a chemical reaction is represented in memory as a :class:`Reaction` @@ -53,7 +53,7 @@ from rmgpy.exceptions import ReactionError, KineticsError from rmgpy.kinetics import KineticsData, ArrheniusBM, ArrheniusEP, ThirdBody, Lindemann, Troe, Chebyshev, \ PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, get_rate_coefficient_units_from_reaction_order, \ - SurfaceArrheniusBEP, StickingCoefficientBEP + SurfaceArrheniusBEP, StickingCoefficientBEP, ArrheniusChargeTransfer, ArrheniusChargeTransferBM 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.diffusionLimited import diffusion_limiter @@ -68,7 +68,7 @@ class Reaction: """ A chemical reaction. The attributes are: - + =================== =========================== ============================ Attribute Type Description =================== =========================== ============================ @@ -92,7 +92,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, @@ -169,7 +169,7 @@ def __str__(self): def to_labeled_str(self, use_index=False): """ - the same as __str__ except that the labels are assumed to exist and used for reactant and products rather than + the same as __str__ except that the labels are assumed to exist and used for reactant and products rather than the labels plus the index in parentheses """ arrow = ' <=> ' if self.reversible else ' => ' @@ -234,7 +234,7 @@ def degeneracy(self, new): def to_chemkin(self, species_list=None, kinetics=True): """ Return the chemkin-formatted string for this reaction. - + If `kinetics` is set to True, the chemkin format kinetics will also be returned (requires the `species_list` to figure out third body colliders.) Otherwise, only the reaction string will be returned. @@ -321,9 +321,9 @@ def to_cantera(self, species_list=None, use_chemkin_identifier=False): if isinstance(ct_reaction, list): for rxn in ct_reaction: rxn.reversible = self.reversible - # Set the duplicate flag to true since this reaction comes from multiarrhenius or multipdeparrhenius + # Set the duplicate flag to true since this reaction comes from multiarrhenius or multipdeparrhenius rxn.duplicate = True - # Set the ID flag to the original rmg index + # Set the ID flag to the original rmg index rxn.ID = str(self.index) else: ct_reaction.reversible = self.reversible @@ -662,7 +662,7 @@ def get_rate_coefficient(self, T, P=0, surface_site_density=0): Return the overall rate coefficient for the forward reaction at temperature `T` in K and pressure `P` in Pa, including any reaction path degeneracies. - + If diffusion_limiter is enabled, the reaction is in the liquid phase and we use a diffusion limitation to correct the rate. If not, then use the intrinsic rate coefficient. @@ -672,7 +672,7 @@ def get_rate_coefficient(self, T, P=0, surface_site_density=0): """ if isinstance(self.kinetics, StickingCoefficient): if surface_site_density <= 0: - raise ValueError("Please provide a postive surface site density in mol/m^2 " + raise ValueError("Please provide a postive surface site density in mol/m^2 " f"for calculating the rate coefficient of {StickingCoefficient.__name__} kinetics") else: return self.get_surface_rate_coefficient(T, surface_site_density) @@ -757,7 +757,7 @@ def fix_barrier_height(self, force_positive=False): """ Turns the kinetics into Arrhenius (if they were ArrheniusEP) and ensures the activation energy is at least the endothermicity - for endothermic reactions, and is not negative only as a result + for endothermic reactions, and is not negative only as a result of using Evans Polanyi with an exothermic reaction. If `force_positive` is True, then all reactions are forced to have a non-negative barrier. @@ -766,47 +766,44 @@ def fix_barrier_height(self, force_positive=False): 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)): - 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)) + 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, ArrheniusChargeTransferBM): + Ea = self.kinetics.E0.value_si # temporarily using Ea to store the intrinsic barrier height E0 + self.kinetics = self.kinetics.to_arrhenius_charge_transfer(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, ArrheniusChargeTransfer)): # 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 def reverse_arrhenius_rate(self, k_forward, reverse_units, Tmin=None, Tmax=None): @@ -880,9 +877,60 @@ 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 +<<<<<<< HEAD +======= + 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=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 reverse_arrhenius_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(Tlist=np.ndarray, klist=np.ndarray, i=cython.int, V0=cython.double) + kf = k_forward + if not isinstance(kf, ArrheniusChargeTransfer): # Only reverse SurfaceChargeTransfer rates + raise TypeError(f'Expected a ArrheniusChargeTransfer 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 = ArrheniusChargeTransfer(alpha=kf.alpha.value, electrons=-1*self.electrons, V0=(V0,'V')) + kr.fit_to_data(Tlist, klist, reverse_units, kf.T0.value_si) + return kr + +>>>>>>> 3afc9cc91 (handle ChargeTransfer kinetics with in reaction.py) 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. @@ -905,6 +953,7 @@ def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, T Lindemann.__name__, Troe.__name__, StickingCoefficient.__name__, + ArrheniusChargeTransfer.__name__, ) # Get the units for the reverse rate coefficient @@ -912,7 +961,7 @@ def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, T surf_prods = [spcs for spcs in self.products if spcs.contains_surface_site()] except IndexError: surf_prods = [] - logging.warning(f"Species do not have an rmgpy.molecule.Molecule " + logging.warning(f"Species do not have an rmgpy.molecule.Molecule " "Cannot determine phases of species. We will assume gas" ) n_surf = len(surf_prods) @@ -920,7 +969,18 @@ 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 +<<<<<<< HEAD if isinstance(kf, KineticsData): +======= + + if isinstance(kf, SurfaceChargeTransfer): + return self.reverse_surface_charge_transfer_rate(kf, kunits, Tmin, Tmax) + + elif isinstance(kf, ArrheniusChargeTransfer): + return self.reverse_arrhenius_charge_transfer_rate(kf, kunits, Tmin, Tmax) + + elif isinstance(kf, KineticsData): +>>>>>>> 3afc9cc91 (handle ChargeTransfer kinetics with in reaction.py) Tlist = kf.Tdata.value_si klist = np.zeros_like(Tlist) @@ -939,7 +999,7 @@ def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, T elif isinstance(kf, StickingCoefficient): if surface_site_density <= 0: - raise ValueError("Please provide a postive surface site density in mol/m^2 " + raise ValueError("Please provide a postive surface site density in mol/m^2 " f"for calculating the rate coefficient of {StickingCoefficient.__name__} kinetics") else: return self.reverse_sticking_coeff_rate(kf, kunits, surface_site_density, Tmin, Tmax) @@ -1068,14 +1128,14 @@ def calculate_microcanonical_rate_coefficient(self, e_list, j_list, reac_dens_st reactant density of states is required; if the reaction is reversible, then both are required. This function will try to use the best method that it can based on the input data available: - + * If detailed information has been provided for the transition state (i.e. the molecular degrees of freedom), then RRKM theory will be used. - + * If the above is not possible but high-pressure limit kinetics - :math:`k_\\infty(T)` have been provided, then the inverse Laplace + :math:`k_\\infty(T)` have been provided, then the inverse Laplace transform method will be used. - + The density of states for the product `prod_dens_states` and the temperature of interest `T` in K can also be provided. For isomerization and association reactions `prod_dens_states` is required; for dissociation reactions it is @@ -1118,6 +1178,14 @@ def is_balanced(self): if reactant_elements[element] != product_elements[element]: return False +<<<<<<< HEAD +======= + if self.electrons < 0: + reactants_net_charge += self.electrons + elif self.electrons > 0: + products_net_charge -= self.electrons + +>>>>>>> 3afc9cc91 (handle ChargeTransfer kinetics with in reaction.py) return True def generate_pairs(self): @@ -1125,7 +1193,7 @@ def generate_pairs(self): Generate the reactant-product pairs to use for this reaction when performing flux analysis. The exact procedure for doing so depends on the reaction type: - + =================== =============== ======================================== Reaction type Template Resulting pairs =================== =============== ======================================== @@ -1134,8 +1202,8 @@ def generate_pairs(self): Association A + B -> C (A,C), (B,C) Bimolecular A + B -> C + D (A,C), (B,D) *or* (A,D), (B,C) =================== =============== ======================================== - - There are a number of ways of determining the correct pairing for + + There are a number of ways of determining the correct pairing for bimolecular reactions. Here we try a simple similarity analysis by comparing the number of heavy atoms. This should work most of the time, but a more rigorous algorithm may be needed for some cases. @@ -1195,9 +1263,9 @@ def _repr_png_(self): # Build the transition state geometry def generate_3d_ts(self, reactants, products): """ - Generate the 3D structure of the transition state. Called from + Generate the 3D structure of the transition state. Called from model.generate_kinetics(). - + self.reactants is a list of reactants self.products is a list of products """ @@ -1207,7 +1275,7 @@ def generate_3d_ts(self, reactants, products): atoms involved in the reaction. If a radical is involved, can find the atom with radical electrons. If a more reliable method can be found, would greatly improve the method. - + Repeat for the products """ for i in range(0, len(reactants)): diff --git a/rmgpy/yml.py b/rmgpy/yml.py index d1873805f8..3243042042 100644 --- a/rmgpy/yml.py +++ b/rmgpy/yml.py @@ -40,7 +40,7 @@ from rmgpy.reaction import Reaction from rmgpy.thermo.nasa import NASAPolynomial, NASA from rmgpy.thermo.wilhoit import Wilhoit -from rmgpy.kinetics.arrhenius import Arrhenius, PDepArrhenius, MultiArrhenius, MultiPDepArrhenius +from rmgpy.kinetics.arrhenius import Arrhenius, PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, ArrheniusChargeTransfer from rmgpy.kinetics.falloff import Troe, ThirdBody, Lindemann from rmgpy.kinetics.chebyshev import Chebyshev from rmgpy.data.solvation import SolventData @@ -136,6 +136,21 @@ 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, ArrheniusChargeTransfer): + obj.change_t0(1.0) + obj.change_v0(0.0) + result_dict["type"] = "Arrheniusq" + result_dict["A"] = obj.A.value_si + result_dict["Ea"] = obj.Ea.value_si + result_dict["n"] = obj.n.value_si + result_dict["q"] = obj._alpha.value_si*obj._electrons.value_si + elif isinstance(obj, SurfaceChargeTransfer): + obj.change_v0(0.0) + result_dict["type"] = "Arrheniusq" + result_dict["A"] = obj.A.value_si + result_dict["Ea"] = obj.Ea.value_si + result_dict["n"] = obj.n.value_si + result_dict["q"] = obj._alpha.value_si*obj._electrons.value_si elif isinstance(obj, StickingCoefficient): obj.change_t0(1.0) result_dict["type"] = "StickingCoefficient"