Skip to content

Merging Anna's changes from 'Pdep spin 39' to the main RMG branch #2837

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 6 commits into
base: main
Choose a base branch
from
Draft

Merging Anna's changes from 'Pdep spin 39' to the main RMG branch #2837

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b972628
fixed functions in fragment_utils.py
donerancl Feb 10, 2025
650e116
moved ammonia project PR changes to new branch after update of main t…
donerancl Jun 2, 2025
15cdb8c
pdep spin conservation updates
donerancl Jun 23, 2025
edc03c4
debugging pdep
donerancl Jun 23, 2025
aac9a56
revert update_unimolecular_reaction_network
donerancl Jun 25, 2025
960dfb4
revrt def merge
donerancl Jun 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
75 changes: 66 additions & 9 deletions rmgpy/data/kinetics/family.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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([]):
Copy link

Copilot AI Jul 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use set() instead of set([]) for creating an empty set, or use the more Pythonic if reactants_combined_spin.intersection(products_combined_spin):.

Suggested change
if reactants_combined_spin.intersection(products_combined_spin) != set([]):
if reactants_combined_spin.intersection(products_combined_spin):

Copilot uses AI. Check for mistakes.
return True
else:
logging.debug(f"Reactants combined spin is {reactants_combined_spin}, but the products combined spin is {products_combined_spin}")
Copy link

Copilot AI Jul 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The log message uses 'but' which implies a problem, but this might be called in valid cases. Consider using more neutral language like 'and' instead of 'but'.

Suggested change
logging.debug(f"Reactants combined spin is {reactants_combined_spin}, but the products combined spin is {products_combined_spin}")
logging.debug(f"Reactants combined spin is {reactants_combined_spin}, and the products combined spin is {products_combined_spin}")

Copilot uses AI. Check for mistakes.
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
Comment on lines +228 to +241
Copy link

Copilot AI Jul 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This method returns None for reactions with more than 2 reactants/products, but the caller may not handle this case properly. Consider raising an exception or documenting this limitation clearly.

Suggested change
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
"""
Calculate the combined spin multiplicities for reactants and products.
This method supports reactions with up to two reactants and two products.
If there are more than two reactants or products, a ValueError is raised.
Returns:
tuple: A tuple containing sets of combined spin multiplicities for reactants and products.
Raises:
ValueError: If there are more than two reactants or products.
"""
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:
raise ValueError("calculate_combined_spin only supports up to two reactants.")
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:
raise ValueError("calculate_combined_spin only supports up to two products.")

Copilot uses AI. Check for mistakes.
Comment on lines +228 to +241
Copy link

Copilot AI Jul 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This method returns None for reactions with more than 2 reactants/products, but the caller may not handle this case properly. Consider raising an exception or documenting this limitation clearly.

Suggested change
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
"""
Calculate the combined spin states for reactants and products.
This method supports reactions with up to two reactants and two products.
If there are more than two reactants or products, a ValueError is raised.
Returns:
tuple: A tuple containing the combined spin states for reactants and products.
Raises:
ValueError: If there are more than two reactants or products.
"""
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:
raise ValueError("calculate_combined_spin only supports up to two reactants.")
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:
raise ValueError("calculate_combined_spin only supports up to two products.")

Copilot uses AI. Check for mistakes.
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
Expand Down Expand Up @@ -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`
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand All @@ -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)

Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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]),
}
12 changes: 6 additions & 6 deletions rmgpy/molecule/fragment_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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
)

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down
6 changes: 4 additions & 2 deletions rmgpy/rmg/model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Copy link

Copilot AI Jul 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The comment explaining why this conversion is needed was removed. Consider restoring it as it provides important context for this non-obvious operation.

Suggested change
except AttributeError:
except AttributeError:
# Convert the Molecule to a Fragment using its adjacency list.
# This is necessary because the `cut_molecule` method is only available for Fragment objects.

Copilot uses AI. Check for mistakes.
# 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:
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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.
Expand Down
12 changes: 2 additions & 10 deletions rmgpy/rmg/pdep.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Copy link

Copilot AI Jul 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Trailing whitespace should be removed from this line.

Suggested change
E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.kinetics.Ea.value_si
E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.kinetics.Ea.value_si

Copilot uses AI. Check for mistakes.
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
Comment on lines +846 to +848
Copy link

Copilot AI Jul 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Trailing whitespace should be removed from this line.

Copilot uses AI. Check for mistakes.
rxn.transition_state = rmgpy.species.TransitionState(conformer=Conformer(E0=(E0 * 0.001, "kJ/mol")))

# Set collision model
Expand Down
Loading