From 088cad7bec2c446e0a8112c1d0bed0f85bfe6b3e Mon Sep 17 00:00:00 2001 From: Ryan Gillis Date: Thu, 21 Mar 2019 11:43:15 -0400 Subject: [PATCH 1/5] When estimating thermochemistry from GAV prioritize resonance forms that do not contain charged atoms --- rmgpy/data/thermo.py | 2 ++ rmgpy/data/thermoTest.py | 10 ++++++++++ rmgpy/molecule/molecule.py | 6 ++++++ 3 files changed, 18 insertions(+) diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index 3929a9d2fa..2fcd66bf1f 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -1806,6 +1806,8 @@ def prioritize_thermo(self, species, thermo_data_list): # For noncyclics, default to original algorithm of ordering thermo based on the most stable enthalpy H298 = np.array([t.get_enthalpy(298.) for t in thermo_data_list]) indices = H298.argsort().tolist() + # Sort indices again by the Molecule.has_charge() + indices.sort(key=lambda index: species.molecule[index].has_charge(), reverse=False) # Sort indices again by the Molecule.reactive flag indices.sort(key=lambda index: species.molecule[index].reactive, reverse=True) else: diff --git a/rmgpy/data/thermoTest.py b/rmgpy/data/thermoTest.py index de41a9a7eb..c99ef1f3e9 100644 --- a/rmgpy/data/thermoTest.py +++ b/rmgpy/data/thermoTest.py @@ -573,6 +573,16 @@ def test_lowest_h298_for_resonance_structures(self): msg="Did not select the molecule with the lowest H298 " "as a the thermo entry for C=C[CH][O] / C=CC=O") + smiles = 'CS(=O)#S(=O)C' + # when using group additivity should instead select CS(=O)S(=O)C which has a lower enthalpy + spec = Species().from_smiles(smiles) + thermo_gav1 = self.database.get_thermo_data_from_groups(spec) + spec.generate_resonance_structures() + thermo_gav2 = self.database.get_thermo_data_from_groups(spec) + self.assertTrue(thermo_gav2.get_enthalpy(298) < thermo_gav1.get_enthalpy(298), + msg="Did not select the molecule with the lowest H298 " + "as a the thermo entry for CS(=O)S(=O)C / CS(=O)#S(=O)C") + def test_thermo_for_mixed_reactive_and_nonreactive_molecules(self): """Test that the thermo entry of nonreactive molecules isn't selected for a species, even if it's more stable""" diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index d5bb05ad72..de5fe4a4fc 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -2066,6 +2066,12 @@ def has_lone_pairs(self): if atom.lone_pairs > 0: return True return False + + def has_charge(self): + for atom in self.vertices: + if atom.charge != 0: + return True + return False def is_aryl_radical(self, aromatic_rings=None): """ From 42ab29cbce0c0e9170c1c7e8681af652e76e425c Mon Sep 17 00:00:00 2001 From: Ryan Gillis Date: Mon, 22 Jul 2019 14:52:43 -0400 Subject: [PATCH 2/5] Allows constant concentration species in simple gas phase reactors --- rmgpy/rmg/input.py | 21 ++++++++++++++------- rmgpy/solver/mbSampled.pyx | 4 +++- rmgpy/solver/simple.pyx | 23 ++++++++++++++++++++++- rmgpy/tools/simulate.py | 8 ++++---- 4 files changed, 43 insertions(+), 13 deletions(-) diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index f20c496ca8..ac848fd7f2 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -190,7 +190,7 @@ def simple_reactor(temperature, sensitivityTemperature=None, sensitivityPressure=None, sensitivityMoleFractions=None, - ): + constantSpecies=None): logging.debug('Found SimpleReactor reaction system') for key, value in initialMoleFractions.items(): @@ -249,6 +249,14 @@ def simple_reactor(temperature, else: sensitive_species.append('all') + #Check the constant species exist + if constantSpecies is not None: + logging.debug(' Generation with constant species:') + for const_spc in constantSpecies: + logging.debug(" {0}".format(const_spc)) + if const_spc not in species_dict: + raise InputError('Species {0} not found in the input file'.format(const_spc)) + if not isinstance(T, list): sensitivityTemperature = T if not isinstance(P, list): @@ -262,8 +270,7 @@ def simple_reactor(temperature, sens_conditions['T'] = Quantity(sensitivityTemperature).value_si sens_conditions['P'] = Quantity(sensitivityPressure).value_si - system = SimpleReactor(T, P, initialMoleFractions, nSims, termination, sensitive_species, sensitivityThreshold, - sens_conditions) + system = SimpleReactor(T, P, initialMoleFractions, nSims, termination, sensitive_species, sensitivityThreshold, sens_conditions, constantSpecies) rmg.reaction_systems.append(system) assert balanceSpecies is None or isinstance(balanceSpecies, str), 'balanceSpecies should be the string corresponding to a single species' @@ -340,10 +347,10 @@ def liquid_reactor(temperature, # chatelak: check the constant species exist if constantSpecies is not None: logging.debug(' Generation with constant species:') - for constantSpecie in constantSpecies: - logging.debug(" {0}".format(constantSpecie)) - if constantSpecie not in species_dict: - raise InputError('Species {0} not found in the input file'.format(constantSpecie)) + for const_spc in constantSpecies: + logging.debug(" {0}".format(const_spc)) + if const_spc not in species_dict: + raise InputError('Species {0} not found in the input file'.format(const_spc)) if not isinstance(T, list): sensitivityTemperature = T diff --git a/rmgpy/solver/mbSampled.pyx b/rmgpy/solver/mbSampled.pyx index a467944065..173c124b2b 100644 --- a/rmgpy/solver/mbSampled.pyx +++ b/rmgpy/solver/mbSampled.pyx @@ -64,6 +64,7 @@ cdef class MBSampledReactor(ReactionSystem): cdef public bint constant_volume cdef public dict initial_mole_fractions cdef public list constantSpeciesList + cdef public list const_spc_names # collider variables @@ -112,13 +113,14 @@ cdef class MBSampledReactor(ReactionSystem): cdef public np.ndarray pdep_specific_collider_reaction_indices def __init__(self, T, P, initial_mole_fractions, k_sampling, constantSpeciesList, termination, sensitive_species=None, - sensitivity_threshold=1e-3): + sensitivity_threshold=1e-3, const_spc_names=None): ReactionSystem.__init__(self, termination, sensitive_species, sensitivity_threshold) self.T = Quantity(T) self.P = Quantity(P) self.initial_mole_fractions = initial_mole_fractions self.k_sampling = RateCoefficient(k_sampling) self.constantSpeciesList = constantSpeciesList + self.const_spc_names = const_spc_names self.V = 0 # will be set in initialize_model self.constant_volume = False diff --git a/rmgpy/solver/simple.pyx b/rmgpy/solver/simple.pyx index 64833b1368..db19c47532 100644 --- a/rmgpy/solver/simple.pyx +++ b/rmgpy/solver/simple.pyx @@ -56,6 +56,8 @@ cdef class SimpleReactor(ReactionSystem): cdef public double V cdef public bint constant_volume cdef public dict initial_mole_fractions + cdef public list const_spc_names + cdef public list const_spc_indices # collider variables @@ -110,7 +112,7 @@ cdef class SimpleReactor(ReactionSystem): cdef public int n_sims def __init__(self, T, P, initial_mole_fractions, n_sims=1, termination=None, sensitive_species=None, - sensitivity_threshold=1e-3, sens_conditions=None): + sensitivity_threshold=1e-3, sens_conditions=None, const_spc_names=None): ReactionSystem.__init__(self, termination, sensitive_species, sensitivity_threshold) if type(T) != list: @@ -125,6 +127,10 @@ cdef class SimpleReactor(ReactionSystem): self.initial_mole_fractions = initial_mole_fractions + #Constant Species Properties + self.const_spc_indices = None + self.const_spc_names = const_spc_names #store index of constant species + self.V = 0 # will be set in initialize_model self.constant_volume = False @@ -163,6 +169,17 @@ cdef class SimpleReactor(ReactionSystem): conditions[species_dict[label]] = value self.sens_conditions = conditions + def get_const_spc_indices (self, core_species): + """ + Allow to identify constant Species position in solver + """ + for spc in self.const_spc_names: + if self.const_spc_indices is None: # initialize once the list if constant SPC declared + self.const_spc_indices = [] + for spc in core_species: #Need to identify the species object corresponding to the the string written in the input file + if spc.label == spc: + self.const_spc_indices.append(core_species.index(spc)) + cpdef initialize_model(self, list core_species, list core_reactions, list edge_species, list edge_reactions, list surface_species=None, list surface_reactions=None, list pdep_networks=None, atol=1e-16, rtol=1e-8, sensitivity=False, sens_atol=1e-6, sens_rtol=1e-4, @@ -504,6 +521,10 @@ cdef class SimpleReactor(ReactionSystem): reaction_rate = k * C[inet[j, 0]] * C[inet[j, 1]] * C[inet[j, 2]] network_leak_rates[j] = reaction_rate + if self.const_spc_indices is not None: + for spc_index in self.const_spc_indices: + core_species_rates[spc_index] = 0 + self.core_species_concentrations = core_species_concentrations self.core_species_rates = core_species_rates self.core_species_production_rates = core_species_production_rates diff --git a/rmgpy/tools/simulate.py b/rmgpy/tools/simulate.py index 943605f757..984af5feee 100644 --- a/rmgpy/tools/simulate.py +++ b/rmgpy/tools/simulate.py @@ -84,14 +84,14 @@ def simulate(rmg, diffusion_limited=True): rmg.load_database() solvent_data = rmg.database.solvation.get_solvent_data(rmg.solvent) diffusion_limiter.enable(solvent_data, rmg.database.solvation) - - # Store constant species indices - if reaction_system.const_spc_names is not None: - reaction_system.get_const_spc_indices(rmg.reaction_model.core.species) elif rmg.uncertainty is not None: rmg.verbose_comments = True rmg.load_database() + # Store constant species indices + if reaction_system.const_spc_names is not None: + reaction_system.get_const_spc_indices(rmg.reaction_model.core.species) + reaction_system.simulate( core_species=rmg.reaction_model.core.species, core_reactions=rmg.reaction_model.core.reactions, From 39873608a97a5831a87231e346477e8c1df3393d Mon Sep 17 00:00:00 2001 From: Ryan Gillis Date: Fri, 5 Apr 2019 10:35:09 -0400 Subject: [PATCH 3/5] Allows valence 12 sulfur species, i.e. DMSO2 --- rmgpy/molecule/filtration.py | 2 +- rmgpy/molecule/filtrationTest.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/molecule/filtration.py b/rmgpy/molecule/filtration.py index b697285123..8eedd12948 100644 --- a/rmgpy/molecule/filtration.py +++ b/rmgpy/molecule/filtration.py @@ -130,7 +130,7 @@ def get_octet_deviation(mol, allow_expanded_octet=True): # This results in O=S=O as a representative structure for SO2 rather than O=[:S+][:::O-], # and in CS(=O)C as a representative structure for DMSO rather than C[:S+]([:::O-])C. if atom.lone_pairs <= 1: - octet_deviation += min(abs(8 - val_electrons), abs(10 - val_electrons)) # octet/dectet on S p[0,1] + octet_deviation += min(abs(8 - val_electrons), abs(10 - val_electrons), abs(12 - val_electrons)) # octet/dectet on S p[0,1] # eg [O-][S+]=O, O[S]=O, OS([O])=O, O=S(=O)(O)O elif atom.lone_pairs >= 2: octet_deviation += abs(8 - val_electrons) # octet on S p[2,3] diff --git a/rmgpy/molecule/filtrationTest.py b/rmgpy/molecule/filtrationTest.py index 8e74834cd2..3a1f50634e 100644 --- a/rmgpy/molecule/filtrationTest.py +++ b/rmgpy/molecule/filtrationTest.py @@ -127,7 +127,7 @@ def penalty_for_s_triple_s_test(self): """ mol = Molecule().from_adjacency_list(adj) octet_deviation = get_octet_deviation(mol) - self.assertEqual(octet_deviation, 3.0) + self.assertEqual(octet_deviation, 1.0) def radical_site_test(self): """Test that a charged molecule isn't filtered if it introduces new radical site""" From 27a0392865adf2e686b68494ec7a5d76accbf0c9 Mon Sep 17 00:00:00 2001 From: Ryan Gillis Date: Tue, 8 Oct 2019 09:59:14 -0400 Subject: [PATCH 4/5] prevents resonance generation of S#S structures by find_lone_pair_multiplebond_paths and find_adj_lone_pair_multiple_bond_delocalization_paths --- rmgpy/molecule/pathfinder.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/rmgpy/molecule/pathfinder.py b/rmgpy/molecule/pathfinder.py index 36022935db..64db88c1cb 100644 --- a/rmgpy/molecule/pathfinder.py +++ b/rmgpy/molecule/pathfinder.py @@ -292,8 +292,8 @@ def find_lone_pair_multiple_bond_paths(atom1): paths = [] for atom2, bond12 in atom1.edges.items(): - # Bond must be capable of gaining an order - if bond12.is_single() or bond12.is_double(): + #If both atom1 and atom2 are sulfur then don't do this type of resonance. Also, the bond must be capable of gaining an order. + if (not atom1.is_sulfur() or not atom2.is_sulfur()) and (bond12.is_single() or bond12.is_double()): for atom3, bond23 in atom2.edges.items(): # Bond must be capable of losing an order without breaking, atom3 must be able to gain a lone pair if atom1 is not atom3 and (bond23.is_double() or bond23.is_triple()) \ @@ -381,8 +381,9 @@ def find_adj_lone_pair_multiple_bond_delocalization_paths(atom1): if atom2.is_non_hydrogen(): # don't bother with hydrogen atoms. # Find paths in the direction the bond order, # atom1 must posses at least one lone pair to loose it + # the final clause of this prevents S#S from forming by this resonance pathway if ((bond12.is_single() or bond12.is_double()) - and is_atom_able_to_lose_lone_pair(atom1)): + and is_atom_able_to_lose_lone_pair(atom1)) and not (atom1.is_sulfur() and atom2.is_sulfur() and bond12.is_double()): paths.append([atom1, atom2, bond12, 1]) # direction = 1 # Find paths in the direction the bond order, # atom1 gains a lone pair, hence cannot already have more than two lone pairs From f61752e19812dd908d611d6d967c344cfb458c87 Mon Sep 17 00:00:00 2001 From: Ryan Gillis Date: Tue, 8 Oct 2019 17:07:15 -0400 Subject: [PATCH 5/5] Adds a unit test to confirm that S#S is not formed by the lone pair multiple bond resonance pathways --- rmgpy/molecule/resonanceTest.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/rmgpy/molecule/resonanceTest.py b/rmgpy/molecule/resonanceTest.py index cd72f47370..c633735aef 100644 --- a/rmgpy/molecule/resonanceTest.py +++ b/rmgpy/molecule/resonanceTest.py @@ -1419,3 +1419,11 @@ def test_surface_o(self): 1 X u0 p0 c0 {2,D} 2 O u0 p2 c0 {1,D}""")) self.assertEquals(len(mol_list), 1) + + def test_sulfur_triple_bond(self): + """ + Test the prevention of S#S formation through the find_lone_pair_multiplebond_paths and + find_adj_lone_pair_multiple_bond_delocalization_paths + """ + mol_list = generate_resonance_structures(Molecule(smiles="S1SSS1"), filter_structures=False) + self.assertEqual(len(mol_list), 10)