diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 102c5f13ff..5d63574db1 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -214,6 +214,32 @@ def copy(self): return other + def check_if_spin_allowed(self): + # get the combined spin for reactants and products + reactants_combined_spin, products_combined_spin = self.calculate_combined_spin() + # check if there are any matches for combined spin between reactants and products + if reactants_combined_spin.intersection(products_combined_spin) != set([]): + return True + else: + logging.debug(f"Reactants combined spin is {reactants_combined_spin}, but the products combined spin is {products_combined_spin}") + return False + + def calculate_combined_spin(self): + if len(self.reactants) == 1: + reactant_combined_spin = {self.reactants[0].multiplicity} + elif len(self.reactants) == 2: + reactant_spin_string = "+".join(sorted([str(reactant.multiplicity) for reactant in self.reactants])) + reactant_combined_spin = allowed_spin[reactant_spin_string] + else: + return None + if len(self.products) == 1: + product_combined_spin = {self.products[0].multiplicity} + elif len(self.products) == 2: + product_spin_string = "+".join(sorted([str(product.multiplicity) for product in self.products])) + product_combined_spin = allowed_spin[product_spin_string] + else: + return None + return reactant_combined_spin, product_combined_spin def apply_solvent_correction(self, solvent): """ apply kinetic solvent correction in this case the parameters are dGTSsite instead of GTS @@ -1678,7 +1704,7 @@ def is_molecule_forbidden(self, molecule): return False - def _create_reaction(self, reactants, products, is_forward): + def _create_reaction(self, reactants, products, is_forward, check_spin = True): """ Create and return a new :class:`Reaction` object containing the provided `reactants` and `products` as lists of :class:`Molecule` @@ -1713,7 +1739,11 @@ def _create_reaction(self, reactants, products, is_forward): for key, species_list in zip(['reactants', 'products'], [reaction.reactants, reaction.products]): for species in species_list: reaction.labeled_atoms[key] = dict(reaction.labeled_atoms[key], **species.get_all_labeled_atoms()) - + if check_spin: + if not reaction.check_if_spin_allowed(): + logging.info("Did not create the following reaction, which violates conservation of spin...") + logging.info(str(reaction)) + return None return reaction def _match_reactant_to_template(self, reactant, template_reactant): @@ -1776,6 +1806,7 @@ def generate_reactions(self, reactants, products=None, prod_resonance=True, dele specified reactants and products within this family. Degenerate reactions are returned as separate reactions. """ + check_spin = True reaction_list = [] # Forward direction (the direction in which kinetics is defined) @@ -1963,7 +1994,7 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson specified reactants and products within this family. Degenerate reactions are returned as separate reactions. """ - + check_spin = True rxn_list = [] # Wrap each reactant in a list if not already done (this is done to @@ -2019,7 +2050,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + if self.label in allowed_spin_violation_families: + check_spin = False + rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin) if rxn: rxn_list.append(rxn) # Bimolecular reactants: A + B --> products @@ -2062,7 +2095,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + if self.label in allowed_spin_violation_families: + check_spin = False + rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin) if rxn: rxn_list.append(rxn) @@ -2086,8 +2121,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, - forward) + if self.label in allowed_spin_violation_families: + check_spin = False + rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin) if rxn: rxn_list.append(rxn) @@ -2140,7 +2176,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + if self.label in allowed_spin_violation_families: + check_spin = False + rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin) if rxn: rxn_list.append(rxn) else: @@ -2205,7 +2243,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + if self.label in allowed_spin_violation_families: + check_spin = False + rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin) if rxn: rxn_list.append(rxn) @@ -4902,3 +4942,20 @@ def get_site_solute_data(rxn): return site_data else: return None + +allowed_spin_violation_families =['1,4_Cyclic_birad_scission'] + +allowed_spin = { + "1+1": set([1]), + "1+2": set([2]), + "1+3": set([3]), + "1+4": set([4]), + "1+5": set([5]), + "2+2": set([1,3]), + "2+3": set([2,4]), + "2+4": set([3,5]), + "2+5": set([4,6]), + "3+3": set([1,3,5]), + "3+4": set([2,4,6]), + "3+5": set([3,5,7]), +} \ No newline at end of file diff --git a/rmgpy/molecule/fragment_utils.py b/rmgpy/molecule/fragment_utils.py index 25fe5e5cb5..201f63ec5f 100644 --- a/rmgpy/molecule/fragment_utils.py +++ b/rmgpy/molecule/fragment_utils.py @@ -98,7 +98,7 @@ def match_concentrations_with_same_sums(conc1, conc2, rtol=1e-6): seq1 = [tup[1] for tup in conc1] seq2 = [tup[1] for tup in conc2] - matches_seq = FragList.match_sequences(seq1, seq2, rtol) + matches_seq = match_sequences(seq1, seq2, rtol) matches_conc = [] for match_seq in matches_seq: @@ -184,7 +184,7 @@ def match_concentrations_with_different_sums(conc1, conc2): # let matches_conc match with remaining seq2 elif pin1 == len(seq1) and pin2 < len(seq2): remain_conc2 = [(labels2[pin2], val2)] + conc2[(pin2 + 1):] - matches_conc = FragList.match_concentrations_with_different_sums( + matches_conc = match_concentrations_with_different_sums( matches_conc, remain_conc2 ) @@ -216,7 +216,7 @@ def flatten(combo): return_list = [] for i in combo: if isinstance(i, tuple): - return_list.extend(FragList.flatten(i)) + return_list.extend(flatten(i)) else: return_list.append(i) return return_list @@ -278,10 +278,10 @@ def merge_frag_list(to_be_merged): frag2 = to_be_merged[-1].smiles # last fragment in list if 'R' in frag1 and 'L' in frag2: - newfrag = FragList.merge_frag_to_frag(frag1, frag2, 'L') + newfrag = merge_frag_to_frag(frag1, frag2, 'L') elif 'L' in frag1 and 'R' in frag2: - newfrag = FragList.merge_frag_to_frag(frag1, frag2, 'R') + newfrag = merge_frag_to_frag(frag1, frag2, 'R') # warn user if last two fragments in list cannot be merged (no R/L # combo to be made) @@ -290,7 +290,7 @@ def merge_frag_list(to_be_merged): frag1, frag2)) if 'L' in frag1 and 'L' in frag2: - newfrag = FragList.merge_frag_to_frag( + newfrag = merge_frag_to_frag( frag1.replace('L', 'R'), frag2, 'L') if len(to_be_merged) > 2: cut = len(to_be_merged) - 2 diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index d677755b07..bd4bca92ff 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -311,7 +311,6 @@ def make_new_species(self, object, label="", reactive=True, check_existing=True, try: mols = molecule.cut_molecule(cut_through=False) except AttributeError: - # it's Molecule object, change it to Fragment and then cut molecule = Fragment().from_adjacency_list(molecule.to_adjacency_list()) mols = molecule.cut_molecule(cut_through=False) if len(mols) == 1: @@ -864,7 +863,9 @@ def process_new_reactions(self, new_reactions, new_species, pdep_network=None, g pdep = rxn.generate_high_p_limit_kinetics() elif any([any([x.is_subgraph_isomorphic(q) for q in self.unrealgroups]) for y in rxn.reactants + rxn.products for x in y.molecule]): pdep = False - + elif not rxn.reversible: + pdep = False + # If pressure dependence is on, we only add reactions that are not unimolecular; # unimolecular reactions will be added after processing the associated networks if not pdep: @@ -2030,6 +2031,7 @@ def update_unimolecular_reaction_networks(self): # Move to the next core reaction index += 1 + def mark_chemkin_duplicates(self): """ Check that all reactions that will appear the chemkin output have been checked as duplicates. diff --git a/rmgpy/rmg/pdep.py b/rmgpy/rmg/pdep.py index 834f0f7243..508a1892d6 100644 --- a/rmgpy/rmg/pdep.py +++ b/rmgpy/rmg/pdep.py @@ -816,14 +816,6 @@ def update(self, reaction_model, pdep_settings, requires_rms=False): spec.conformer = Conformer(E0=spec.get_thermo_data().E0) E0.append(spec.conformer.E0.value_si) - # Use the average E0 as the reference energy (`energy_correction`) for the network - # The `energy_correction` will be added to the free energies and enthalpies for each - # configuration in the network. - energy_correction = -np.array(E0).mean() - for spec in self.reactants + self.products + self.isomers: - spec.energy_correction = energy_correction - self.energy_correction = energy_correction - # Determine transition state energies on potential energy surface # In the absence of any better information, we simply set it to # be the reactant ground-state energy + the activation energy @@ -851,9 +843,9 @@ def update(self, reaction_model, pdep_settings, requires_rms=False): 'type "{2!s}".'.format(rxn, self.index, rxn.kinetics.__class__)) rxn.fix_barrier_height(force_positive=True) if rxn.network_kinetics is None: - E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.kinetics.Ea.value_si + energy_correction + E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.kinetics.Ea.value_si else: - E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.network_kinetics.Ea.value_si + energy_correction + E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.network_kinetics.Ea.value_si rxn.transition_state = rmgpy.species.TransitionState(conformer=Conformer(E0=(E0 * 0.001, "kJ/mol"))) # Set collision model