Forbidden structs in input

documentation/source/users/rmg/input.rst

@@ -230,8 +230,69 @@ The following is an example of a typical species item, based on methane using SM
 			"""
 	)
 
-.. _reactionsystem:
+.. _forbidden_structures:
+
+Forbidden Structures
+=====================
+
+RMG exlores a wide variety of structures as it builds a mechanism.
+Sometimes, RMG makes nonsensical, unphysical, or unstable species that make their way into the core.
+
+Luckily, RMG has a remedy to prevent these undesirable structures from making their way into your mechanism.
+RMG's ForbiddenStructures database contains structures that are forbidden from entering the core 
+(see :ref:`kineticsDatabase` for more information)
+
+While this database forbids many undesirable structures, it may not contain a species
+or structure you would like to exclude from your model. If you would like to forbid a structure from your model,
+you can do so by adding it to your input file in the forbidden block.
+
+The label and structure of each forbidden structure must be specified.
+The structure of the forbidden structure can be defined using either SMILES or
+:ref:`adjacencyList <rmgpy.molecule.adjlist>`.
+
+The following is an example of a forbidden structure item, based on bidentate CO2 using SMILES or adjacency list to define the structure::
 
+	forbidden(
+		label='CO2_bidentate',
+		structure=SMILES("O=C(*)O*"),
+	)
+
+	forbidden(
+		label='CO2_bidentate',
+		structure=adjacencyList(
+			"""
+			1 O u0 p2 c0 {2,D}
+			2 C u0 p0 c0 {1,D} {3,S} {4,S}
+			3 X u0 p0 c0 {2,S}
+			4 O u0 p2 c0 {2,S} {5,S}
+			5 X u0 p0 c0 {4,S}
+			"""
+		)
+	)
+
+If you would like to forbid a functional group from your mechanism, use the ``adjacencyListGroup`` syntax in the forbidden block.
+
+The following is an example of a forbidden structure using ``adjacencyListGroup``.
+The ``R-O-Cl`` forbidden structure forbids any species which contain an "R-O-Cl" group from your mechanism.
+::
+
+	forbidden(
+		label = "R-O-Cl,
+		structure=adjacencyListGroup(
+			"""
+			1 R  ux {2,S}
+			2 O  ux {1,S} {3,S}
+			3 Cl ux {2,S}
+			"""
+		)
+	)
+
+Note: If you would like to allow Hypochlorous acid (`ClOH`) in your model and forbid any other species with an "R-O-Cl" group,
+you can do so by changing the `R` atomtype to `R!H` in the ``R-O-Cl`` forbidden structure above.
+Alternatively, you could explicitly allow `ClOH` by creating a `ClOH` species in your input file (see :ref:`species_list`) and add "input species"
+to the `generatedSpeciesConstraints` block discussed here :ref:`miscellaneousoptions`.
+
+.. _reactionsystem:
 
 Reaction System
 ===============
@@ -922,10 +983,6 @@ from either the input file, seed mechanisms, or reaction libraries to bypass the
 Note that this should be done with caution, since the constraints will still apply to subsequent
 products that form.
 
-Note that under all circumstances all forbidden species will still be banned unless they are
-manually removed from the database.  See :ref:`kineticsDatabase` for more information on
-forbidden groups.
-
 By default, the ``allowSingletO2`` flag is set to ``False``.  See :ref:`representing_oxygen` for more information.
 
 

examples/rmg/catalysis/ch4_o2/input.py

@@ -1,7 +1,7 @@
 # Data sources
 database(
     thermoLibraries=['surfaceThermoPt111', 'primaryThermoLibrary', 'thermo_DFT_CCSDTF12_BAC','DFT_QCI_thermo'], # 'surfaceThermoPt' is the default. Thermo data is derived using bindingEnergies for other metals 
-    reactionLibraries = [('Surface/CPOX_Pt/Deutschmann2006', False)], # when Ni is used change the library to Surface/Deutschmann_Ni 
+    reactionLibraries = [('Surface/CPOX_Pt/Deutschmann2006_adjusted', False)], # when Ni is used change the library to Surface/Deutschmann_Ni 
     seedMechanisms = [],
     kineticsDepositories = ['training'],
     kineticsFamilies = ['surface','default'],
@@ -40,6 +40,14 @@
     reactive=True,
     structure=adjacencyList("1 X u0"),
 )
+
+# If you would like to forbid the bidentate form of absorbed CO2 from your model,
+# use the following `CO2_bidentate` forbidden structure
+# forbidden(
+#     label='CO2_bidentate',
+#     structure=SMILES("O=C(*)O*"),
+# )
+
 #----------
 # Reaction systems
 surfaceReactor(

examples/rmg/commented/input.py

@@ -78,6 +78,35 @@
     structure=SMILES("O=C=O")
 )
 
+# You can also list forbidden structures to forbid a certain molecule or set of molecules from your model
+# To prevent a single molecule from your model, use SMILES or adjacencyList to define the strucutre 
+# For example, if you do not want cyclopropyne in your model, you can forbid it using SMILES
+# forbidden(
+#     label='cyclopropyne',
+#     structure=SMILES("C1#CC1")
+# )
+# or adjacencyList
+# forbidden(
+#     label='cyclopropyne',
+#     structure=adjacencyList("""
+#     1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+#     2 C u0 p0 c0 {1,S} {3,T}
+#     3 C u0 p0 c0 {1,S} {2,T}
+#     4 H u0 p0 c0 {1,S}
+#     5 H u0 p0 c0 {1,S}
+#     """),
+# )
+# If you want to exclude not just cyclopropyne but all three member rings, use `adjacencyListGroup` 
+# to define the forbidden group structure
+# forbidden(
+#     label='Three-member Ring',
+#     structure=adjacencyListGroup("""
+#     1 R ux {2,[S,D,T]} {3,[S,D,T]}
+#     2 R ux {1,[S,D,T]} {3,[S,D,T]}
+#     3 R ux {1,[S,D,T]} {2,[S,D,T]}
+#     """),
+# )
+
 # Reaction systems
 # currently RMG models only constant temperature and pressure as homogeneous batch reactors.
 # two options are: simpleReactor for gas phase or liquidReactor for liquid phase

rmgpy/rmg/input.py

@@ -34,8 +34,10 @@
 import numpy as np
 
 from rmgpy import settings
+from rmgpy.data.base import Entry
 from rmgpy.exceptions import DatabaseError, InputError
 from rmgpy.molecule import Molecule
+from rmgpy.molecule.group import Group
 from rmgpy.quantity import Quantity, Energy, RateCoefficient, SurfaceConcentration
 from rmgpy.rmg.model import CoreEdgeReactionModel
 from rmgpy.rmg.settings import ModelSettings, SimulatorSettings
@@ -191,6 +193,11 @@ def species(label, structure, reactive=True):
     rmg.initial_species.append(spec)
     species_dict[label] = spec
 
+def forbidden(label, structure):
+
+    if '+' in label:
+        raise InputError('forbidden species {0} label cannot include a + sign'.format(label))
+    rmg.forbidden_structures.append(Entry(label=label,item=structure))
 
 def smarts(string):
     return Molecule().from_smarts(string)
@@ -207,6 +214,9 @@ def inchi(string):
 def adjacency_list(string):
     return Molecule().from_adjacency_list(string)
 
+def adjacency_list_group(string):
+    return Group().from_adjacency_list(string)
+
 def react(tups):
     if not isinstance(tups, list):
         raise InputError("React takes a list of tuples of species strings.")
@@ -1148,10 +1158,12 @@ def read_input_file(path, rmg0):
         'database': database,
         'catalystProperties': catalyst_properties,
         'species': species,
+        'forbidden': forbidden,
         'SMARTS': smarts,
         'SMILES': smiles,
         'InChI': inchi,
         'adjacencyList': adjacency_list,
+        'adjacencyListGroup': adjacency_list_group,
         'react': react,
         'simpleReactor': simple_reactor,
         'constantVIdealGasReactor' : constant_V_ideal_gas_reactor,

rmgpy/rmg/main.py

@@ -187,6 +187,7 @@ def clear(self):
         self.surface_site_density = None
         self.binding_energies = None
         self.coverage_dependence = False
+        self.forbidden_structures = []
 
         self.reaction_model = None
         self.reaction_systems = None
@@ -417,6 +418,17 @@ def load_database(self):
             global solvent
             solvent = self.solvent
 
+        # add any forbidden structures in the input file to the forbidden structures database
+        for forbidden_structure_entry in self.forbidden_structures:
+            label = forbidden_structure_entry.label
+            if label in self.database.forbidden_structures.entries:
+                raise InputError("""
+        Forbidden structure {0} label is already in the forbidden structure database.
+        Please choose a different label for this structure that is not already in
+        {1}""".format(label,os.path.join(self.database_directory,'forbiddenStructures.py')))
+            logging.info('Adding {0} to the forbidden structures database...'.format(label))
+            self.database.forbidden_structures.entries[label] = forbidden_structure_entry
+
         if self.kinetics_estimator == 'rate rules':
             if '!training' not in self.kinetics_depositories:
                 logging.info('Adding rate rules from training set in kinetics families...')
@@ -587,12 +599,12 @@ def initialize(self, **kwargs):
         for spec in self.initial_species:
             if self.database.forbidden_structures.is_molecule_forbidden(spec.molecule[0]):
                 if 'allowed' in self.species_constraints and 'input species' in self.species_constraints['allowed']:
-                    logging.warning('Input species {0} is globally forbidden.  It will behave as an inert unless found '
-                                    'in a seed mechanism or reaction library.'.format(spec.label))
+                    spec.explicitly_allowed = True
+                    logging.warning("Input species {0} is globally forbidden but will be explicitly allowed "
+                            " since input species are permitted by the user's species constraints".format(spec.label))
                 else:
                     raise ForbiddenStructureException("Input species {0} is globally forbidden. You may explicitly "
-                                                      "allow it, but it will remain inert unless found in a seed "
-                                                      "mechanism or reaction library.".format(spec.label))
+                                                      "allow it by adding 'input species' to the `generatedSpeciesConstraints` `allowed` list.".format(spec.label))
             if fails_species_constraints(spec):
                 if 'allowed' in self.species_constraints and 'input species' in self.species_constraints['allowed']:
                     self.species_constraints['explicitlyAllowedMolecules'].append(spec.molecule[0])