Compute IDT and IDT SA using Cantera

examples/commented/input.yml

@@ -35,7 +35,7 @@ t3:
     adapter: RMGConstantTP  # *required* (this is how SA is requested), can use any implemented simulation adapter
     atol: 1e-6  # optional, default: 1e-6
     rtol: 1e-4  # optional, default: 1e-4
-    global_observables: ['IgD', 'ESR', 'SL']  # optional, only implemented in the Cantera adapter, default: ``None``
+    global_observables: ['IDT', 'ESR', 'SL']  # optional, only implemented in the Cantera adapter, default: ``None``
     SA_threshold: 0.01  # optional, default: 0.01
     pdep_SA_threshold: 0.001  # optional, used to determine wells and reactions (to be implemented) to calculate
                               # thermo and rates for from a PES of a sensitive reaction, default: 0.001
@@ -89,6 +89,8 @@ rmg:
       reactive: true  # optional, default: ``True``
       xyz: [ethane.log]  # each entry could be a string/dictionary XYZ format or a file path to parse the information from
       seed_all_rads: ['radical', 'alkoxyl', 'peroxyl']  # radical derivatives that will be added the RMG input file
+      role: 'fuel'  # optional mark if this is a 'fuel', 'oxygen', or 'nitrogen' **instead** of specifying concentration, default: ``None``
+      equivalence_ratios: [0.5, 1.0, 1.5]  # optional, but must be given if the role is 'fuel', default: ``None``
     - label: OH
       smiles: '[OH]'
       observable: true  # optional, will be used as both SA and UA observable, default: ``False``
@@ -97,13 +99,15 @@ rmg:
     - label: O2
       smiles: '[O][O]'
       concentration: 3.5
+      role: 'oxygen'  # optional mark if this is a 'fuel', 'oxygen', or 'nitrogen' **instead** of specifying concentration, default: ``None``
     - label: N2
       adjlist: |
         1 N u0 p1 c0 {2,T}
         2 N u0 p1 c0 {1,T}
       concentration: 3.5 * 3.76
       constant: true  # optional, only allowed to be ``True`` in liquid phase simulations, default: ``False``
       balance: true  # optional, only allowed to be ``True`` in gas phase simulations, default: ``False``
+      role: 'nitrogen'  # optional mark if this is a 'fuel', 'oxygen', or 'nitrogen' **instead** of specifying concentration, default: ``None``
     - label: water
       smiles: O
       solvent: true  # optional, only allowed to be ``True`` in liquid phase simulations, required for one species for liquid phase simulations, default: ``False``

t3/common.py

@@ -5,7 +5,7 @@
 import datetime
 import os
 import string
-from typing import Dict, Optional, Tuple, Union
+from typing import Dict, List, Optional, Tuple, Union
 
 from rmgpy.species import Species
 
@@ -255,3 +255,72 @@ def get_chem_to_rmg_rxn_index_map(chem_annotated_path: str) -> Dict[int, int]:
                 splits = line.split()
                 rxn_map[int(splits[4].split('#')[1].split(';')[0])] = int(splits[-1].split('#')[1])
     return rxn_map
+
+
+def determine_concentrations_by_equivalence_ratios(species: List[dict]):
+    """
+    Determine species concentrations based on the equivalence ratios if given.
+
+    Args:
+        species (list): Entries are species dictionaries following the schema format for RMG species.
+
+    Returns:
+        dict: A species dictionary with the concentrations as a list (not a range).
+    """
+    objects = {'fuel': None, 'oxygen': None, 'nitrogen': None}
+    for spc in species:
+        if spc['role'] == 'fuel':
+            objects['fuel'] = spc.copy()
+        elif spc['role'] == 'oxygen':
+            objects['oxygen'] = spc.copy()
+        elif spc['role'] == 'nitrogen':
+            objects['nitrogen'] = spc.copy()
+    if objects['fuel'] is not None and objects['oxygen'] is not None and objects['fuel']['equivalence_ratios'] is not None:
+        if objects['fuel']['concentration'] == 0:
+            objects['fuel']['concentration'] = 1
+        o2_stoichiometry = get_o2_stoichiometry(smiles=objects['fuel']['smiles'],
+                                                adjlist=objects['fuel']['adjlist'] if 'adjlist' in objects['fuel'].keys() else None,
+                                                inchi=objects['fuel']['inchi'] if 'inchi' in objects['fuel'].keys() else None
+                                                )
+        objects['oxygen']['concentration'] = [eq_ratio * o2_stoichiometry for eq_ratio in objects['fuel']['equivalence_ratios']]
+        if objects['nitrogen'] is not None and objects['nitrogen']['concentration'] == 0:
+            objects['nitrogen']['concentration'] = [o2 * 3.76 for o2 in objects['oxygen']['concentration']]
+    return objects
+
+
+def get_o2_stoichiometry(smiles: Optional[str] = None,
+                         adjlist: Optional[str] = None,
+                         inchi: Optional[str] = None,
+                         ) -> float:
+    """
+    Get the stoichiometry number of O2 for complete combustion of the ``fuel`` molecule.
+
+    Args:
+        smiles (Optional[str]): A SMILES string for the fuel molecule.
+        adjlist (Optional[str]): An adjacency list string for the fuel molecule.
+        inchi (Optional[str]): An InChI string for the fuel molecule.
+
+    Returns:
+        float: The stoichiometry of O2 for complete combustion.
+    """
+    if smiles is None and adjlist is None and inchi is None:
+        raise ValueError('Must provide either a SMILES, an adjacency list, or an InChI string for the fuel molecule.')
+    if adjlist is None:
+        fuel = Species(smiles=smiles, inchi=inchi)
+    else:
+        fuel = Species().from_adjacency_list(adjlist)
+    c, h, n, o, other = 0, 0, 0, 0, 0
+    for atom in fuel.molecule[0].atoms:
+        if atom.is_carbon():
+            c += 1
+        elif atom.is_hydrogen():
+            h += 1
+        elif atom.is_nitrogen():
+            n += 1
+        elif atom.is_oxygen():
+            o += 1
+        else:
+            other += 1
+    if other:
+        raise ValueError(f'Cannot calculate O2 stoichiometry for {fuel.label} with {other} atoms which are not C/H/N/O.')
+    return 0.5 * (2 * c + 0.5 * h - o)

t3/main.py

@@ -49,7 +49,13 @@
 from arc.species.species import ARCSpecies, check_label
 from arc.species.converter import check_xyz_dict
 
-from t3.common import PROJECTS_BASE_PATH, VALID_CHARS, delete_root_rmg_log, get_species_by_label, time_lapse
+from t3.common import (PROJECTS_BASE_PATH,
+                       VALID_CHARS,
+                       delete_root_rmg_log,
+                       determine_concentrations_by_equivalence_ratios,
+                       get_species_by_label,
+                       time_lapse,
+                       )
 from t3.logger import Logger
 from t3.runners.rmg_runner import rmg_runner
 from t3.schema import InputBase
@@ -182,6 +188,7 @@ def __init__(self,
         self.rmg = self.schema['rmg']
         self.qm = self.schema['qm']
         self.verbose = self.schema['verbose']
+        self.update_species_concentrations()
 
         if clean_dir and os.path.isdir(self.project_directory):
             self.cleanup()
@@ -294,21 +301,35 @@ def execute(self):
                     for species in self.rmg['species']:
                         if species['observable'] or species['SA_observable']:
                             self.sa_observables.append(species['label'])
-
-                simulate_adapter = simulate_factory(simulate_method=self.t3['sensitivity']['adapter'],
-                                                    t3=self.t3,
-                                                    rmg=self.rmg,
-                                                    paths=self.paths,
-                                                    logger=self.logger,
-                                                    atol=self.rmg['model']['atol'],
-                                                    rtol=self.rmg['model']['rtol'],
-                                                    observable_list=self.sa_observables,
-                                                    sa_atol=self.t3['sensitivity']['atol'],
-                                                    sa_rtol=self.t3['sensitivity']['rtol'],
-                                                    global_observables=None,
-                                                    )
-                simulate_adapter.simulate()
-                self.sa_dict = simulate_adapter.get_sa_coefficients()
+                if self.sa_observables:
+                    simulate_adapter = simulate_factory(simulate_method=self.t3['sensitivity']['adapter'],
+                                                        t3=self.t3,
+                                                        rmg=self.rmg,
+                                                        paths=self.paths,
+                                                        logger=self.logger,
+                                                        atol=self.rmg['model']['atol'],
+                                                        rtol=self.rmg['model']['rtol'],
+                                                        observable_list=self.sa_observables,
+                                                        sa_atol=self.t3['sensitivity']['atol'],
+                                                        sa_rtol=self.t3['sensitivity']['rtol'],
+                                                        )
+                    simulate_adapter.simulate()
+                    self.sa_dict = simulate_adapter.get_sa_coefficients()
+                    save_yaml_file(path=self.paths['SA dict'], content=self.sa_dict)
+                if self.t3['sensitivity']['global_observables'] == 'IDT':
+                    simulate_adapter = simulate_factory(simulate_method='CanteraIDT',
+                                                        t3=self.t3,
+                                                        rmg=self.rmg,
+                                                        paths=self.paths,
+                                                        logger=self.logger,
+                                                        atol=self.rmg['model']['atol'],
+                                                        rtol=self.rmg['model']['rtol'],
+                                                        sa_atol=self.t3['sensitivity']['atol'],
+                                                        sa_rtol=self.t3['sensitivity']['rtol'],
+                                                        )
+                    simulate_adapter.simulate()
+                    self.sa_dict_idt = simulate_adapter.get_sa_coefficients()
+                    save_yaml_file(path=self.paths['SA idt dict'], content=self.sa_dict_idt)
 
             additional_calcs_required = self.determine_species_and_reactions_to_calculate()
 
@@ -366,12 +387,15 @@ def set_paths(self,
             'RMG job log': os.path.join(iteration_path, 'RMG', 'job.log'),
             'RMG coll vio': os.path.join(iteration_path, 'RMG', 'collision_rate_violators.log'),
             'RMS': os.path.join(iteration_path, 'RMG', 'rms'),
-            'cantera annotated': os.path.join(iteration_path, 'RMG', 'cantera', 'chem_annotated.cti'),
+            'cantera annotated': os.path.join(iteration_path, 'RMG', 'cantera', 'chem_annotated.yaml'),
             'chem annotated': os.path.join(iteration_path, 'RMG', 'chemkin', 'chem_annotated.inp'),
             'species dict': os.path.join(iteration_path, 'RMG', 'chemkin', 'species_dictionary.txt'),
+            'figs': os.path.join(iteration_path, 'Figures'),
             'SA': os.path.join(iteration_path, 'SA'),
             'SA solver': os.path.join(iteration_path, 'SA', 'solver'),
             'SA input': os.path.join(iteration_path, 'SA', 'input.py'),
+            'SA dict': os.path.join(iteration_path, 'SA', 'sa.yaml'),
+            'SA IDT dict': os.path.join(iteration_path, 'SA', 'sa_idt.yaml'),
             'PDep SA': os.path.join(iteration_path, 'PDep_SA'),
             'ARC': os.path.join(iteration_path, 'ARC'),
             'ARC input': os.path.join(iteration_path, 'ARC', 'input.yml'),
@@ -1500,6 +1524,18 @@ def load_species_and_reactions(self):
                                                   )
                 self.reactions[key] = mod_rxn_dict
 
+    def update_species_concentrations(self):
+        """
+        Update the species concentrations based on species roles and equivalence ratios.
+        """
+        objects = determine_concentrations_by_equivalence_ratios(species=self.rmg['species'])
+        for spc in self.rmg['species']:
+            if spc['role'] == 'fuel' and spc['concentration'] == 0 and objects['fuel'] is not None:
+                spc['concentration'] = objects['fuel']['concentration']
+            elif (spc['role'] == 'oxygen' or (spc['role'] == 'nitrogen' and spc['concentration'] == 0)) \
+                    and objects[spc['role']] is not None:
+                spc['concentration'] = [min(objects[spc['role']]['concentration']), max(objects[spc['role']]['concentration'])]
+
     def check_overtime(self) -> bool:
         """
         Check that the timer hasn't run out.

t3/runners/rmg_runner.py

@@ -12,6 +12,7 @@
 from arc.job.local import _determine_job_id, change_mode, execute_command, parse_running_jobs_ids, submit_job
 
 from t3.imports import local_t3_path, settings, submit_scripts
+from t3.utils.fix_cantera import fix_cantera
 
 if TYPE_CHECKING:
     from t3.logger import Logger
@@ -257,6 +258,7 @@ def rmg_runner(rmg_input_file_path: str,
                t3_project_name: Optional[str] = None,
                rmg_execution_type: Optional[str] = None,
                restart_rmg: bool = False,
+               fix_cantera_model: bool = True,
                ) -> bool:
     """
     Run an RMG job as a subprocess under the rmg_env.
@@ -271,6 +273,7 @@ def rmg_runner(rmg_input_file_path: str,
         t3_project_name (str, optional): The T3 project name, used for setting a job name on the server for the RMG run.
         rmg_execution_type (str, optional): The RMG execution type (incore or local). Also set via settings.py.
         restart_rmg (bool, optional): Whether to restart RMG from seed.
+        fix_cantera_model (bool, optional): Whether to fix the Cantera model before running the simulation.
 
     Returns:
         bool: Whether an exception was raised.
@@ -321,9 +324,11 @@ def rmg_runner(rmg_input_file_path: str,
                       and not(len(rmg_errors) >= 2 and error is not None and error == rmg_errors[-2])
             restart_rmg = False if error is not None and 'Could not find one or more of the required files/directories ' \
                                                          'for restarting from a seed mechanism' in error else True
+        if fix_cantera_model:
+            fix_cantera_model_files(rmg_path=os.path.dirname(rmg_input_file_path))
         return not converged
     else:
-        logger.warning(f'Expected wither "incore" or "local" execution type for RMG, got {rmg_execution_type}.\n'
+        logger.warning(f'Expected either "incore" or "local" execution type for RMG, got {rmg_execution_type}.\n'
                        f'Not executing RMG.')
         return True
 
@@ -376,6 +381,17 @@ def get_new_memory_for_an_rmg_run(job_log_path: str,
     return new_mem
 
 
+def fix_cantera_model_files(rmg_path: str) -> None:
+    """
+    Fix a Cantera model file that has undeclared duplicate reactions.
+
+    Args:
+        rmg_path (str): The path to the RMG folder.
+    """
+    fix_cantera(model_path=os.path.join(rmg_path, 'cantera', 'chem_annotated.yaml'))
+    fix_cantera(model_path=os.path.join(rmg_path, 'cantera', 'chem.yaml'))
+
+
 # def get_names_by_sub_folders(pwd: str) -> List[str]:
 #     """
 #     Get the names of the runs.

t3/schema.py

@@ -71,7 +71,7 @@ class T3Sensitivity(BaseModel):
     adapter: constr(max_length=255) = 'RMGConstantTP'
     atol: confloat(gt=0, lt=1e-1) = 1e-6
     rtol: confloat(gt=0, lt=1e-1) = 1e-4
-    global_observables: Optional[List[constr(min_length=2, max_length=3)]] = None
+    global_observables: Optional[List[constr(min_length=2, max_length=3)]] = None  # ['IDT', 'ESR', 'SL']
     SA_threshold: confloat(gt=0, lt=0.5) = 0.01
     pdep_SA_threshold: Optional[confloat(gt=0, lt=0.5)] = 0.001
     ME_methods: List[constr(min_length=2, max_length=3)] = ['CSE', 'MSC']
@@ -98,8 +98,8 @@ def check_global_observables(cls, value):
         """T3Sensitivity.global_observables validator"""
         if value is not None:
             for i, entry in enumerate(value):
-                if entry.lower() not in ['igd', 'esr', 'sl']:
-                    raise ValueError(f'The global observables list must contain a combination of "IgD", "ESR", and "SL", '
+                if entry.lower() not in ['idt', 'esr', 'sl']:
+                    raise ValueError(f'The global observables list must contain a combination of "IDT", "ESR", and "SL", '
                                      f'Got {entry} in {value}')
                 if entry.lower() in [value[j].lower() for j in range(i)]:
                     raise ValueError(f'The global observables list must not contain repetitions, got {value}')
@@ -170,6 +170,8 @@ class RMGSpecies(BaseModel):
     """
     label: str
     concentration: Union[confloat(ge=0), Tuple[confloat(ge=0), confloat(ge=0)]] = 0
+    equivalence_ratios: Optional[List[confloat(gt=0)]] = None
+    role: Optional[str] = None
     smiles: Optional[str] = None
     inchi: Optional[str] = None
     adjlist: Optional[str] = None
@@ -195,6 +197,15 @@ def check_ranged_concentration_not_constant(cls, value, values):
                              f"Got{label}: {values['concentration']}.")
         return value
 
+    @validator('role')
+    def check_species_role(cls, value, values):
+        """RMGSpecies.role validator"""
+        if value not in ['fuel', 'oxygen', 'nitrogen', None]:
+            raise ValueError(f'The species role must be either "fuel", "oxygen", or "nitrogen".\nGot: {value}')
+        if value == 'fuel' and value['equivalence_ratios'] is None:
+            raise ValueError(f'If the species role is "fuel", then the equivalence ratios must be specified.')
+        return value
+
     @validator('concentration')
     def check_concentration_range_order(cls, value, values):
         """Make sure the concentration range is ordered from the smallest to the largest"""

t3/simulate/adapter.py

@@ -35,15 +35,3 @@ def get_sa_coefficients(self):
              sa_dict (dict): a SA dictionary, whose structure is given in the docstring for T3/t3/main.py
         """
         pass
-
-    @abstractmethod
-    def get_idt_by_T(self):
-        """
-        Finds the ignition point by approximating dT/dt as a first order forward difference
-        and then finds the point of maximum slope.
-
-        Returns:
-            idt_dict (dict): Dictionary whose keys include 'idt' and 'idt_index' and whose values are lists of
-                             the ignition delay time in seconds and index at which idt occurs respectively.
-        """
-        pass