Dms oxy changes

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:

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"""
 

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]

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"""

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):
         """

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 <increasing> 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 <decreasing> the bond order,
             # atom1 gains a lone pair, hence cannot already have more than two lone pairs

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)

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

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

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

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,