Allow specifying species concentration range

docs/how_to.md

@@ -8,8 +8,9 @@ a chemical species. The following are possible keys and corresponding
 values for each species dictionary:
 
 - `label` (str, required): The species label.
-- `concentration` (float): concentration units are mole fraction for gas phase
+- `concentration` (Union[float, Tuple[float, float]]): Concentration units are mole fraction for gas phase
   and mol/cm3 for liquid phase. Defaults to `0`.
+  A concentration range can also be specified (a length-2 list/tuple).
 - `smiles` (str): The
   <a href="https://en.wikipedia.org/wiki/Simplified_molecular-input_line-entry_system" target="_blank">
   SMILES</a> representation.

examples/commented/input.yml

@@ -74,7 +74,7 @@ rmg:
   species:
     - label: ethane
       smiles: CC
-      concentration: 1
+      concentration: [1, 1.75]  # a concentration range can be defined (a length-2 list)
       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

t3/main.py

@@ -707,7 +707,8 @@ def determine_species_to_calculate(self) -> bool:
         if self.t3['options']['all_core_species']:
             for species in self.rmg_species:
                 if self.species_requires_refinement(species=species):
-                    species_keys.append(self.add_species(species=species, reasons=['All core species']))
+                    species_keys.append(self.add_species(species=species,
+                                                         reasons=[f'(i {self.iteration}) All core species']))
         else:
             # 1. SA observables
             sa_observables_exist = False
@@ -729,7 +730,8 @@ def determine_species_to_calculate(self) -> bool:
 
         species_keys = list(set([key for key in species_keys if key is not None]))
 
-        additional_calcs_required = bool(len(species_keys))
+        additional_calcs_required = bool(len(species_keys)) \
+            or any(spc['converged'] is None for spc in self.species.values())
         self.logger.info(f'Additional calculations required: {additional_calcs_required}\n')
         if additional_calcs_required:
             self.logger.log_species_to_calculate(species_keys, self.species)
@@ -1042,7 +1044,8 @@ def species_requires_refinement(self, species: Species) -> bool:
             bool: Whether the species thermochemical properties should be calculated. ``True`` if they should be.
         """
         thermo_comment = species.thermo.comment.split('Solvation')[0]
-        if self.get_species_key(species=species) is None \
+        if (self.get_species_key(species=species) is None
+            or self.species[self.get_species_key(species=species)]['converged'] is None) \
                 and ('group additivity' in thermo_comment or '+ radical(' in thermo_comment):
             return True
         return False
@@ -1139,7 +1142,7 @@ def add_species(self,
             key = len(list(self.species.keys()))
             qm_species = species.copy(deep=False)
             legalize_species_label(species=qm_species)
-            qm_species.label += f'_{key}'
+            qm_species.label = f's{key}_{qm_species.label}'
             self.species[key] = {'RMG label': species.label,
                                  'Chemkin label': species.to_chemkin(),
                                  'QM label': qm_species.label,

t3/schema.py

@@ -6,7 +6,7 @@
 """
 
 from enum import Enum
-from typing import Dict, List, Optional, Union
+from typing import Dict, List, Optional, Tuple, Union
 
 from pydantic import BaseModel, conint, confloat, constr, root_validator, validator
 
@@ -148,7 +148,7 @@ class RMGSpecies(BaseModel):
     A class for validating input.RMG.species arguments
     """
     label: str
-    concentration: confloat(ge=0) = 0
+    concentration: Union[confloat(ge=0), Tuple[confloat(ge=0), confloat(ge=0)]] = 0
     smiles: Optional[str] = None
     inchi: Optional[str] = None
     adjlist: Optional[str] = None
@@ -165,6 +165,37 @@ class RMGSpecies(BaseModel):
     class Config:
         extra = "forbid"
 
+    @validator('constant')
+    def check_ranged_concentration_not_constant(cls, value, values):
+        """RMGSpecies.constant validator"""
+        label = ' for ' + values['label'] if 'label' in values else ''
+        if value and isinstance(values['concentration'], tuple):
+            raise ValueError(f"A constant species cannot have a concentration range.\n"
+                             f"Got{label}: {values['concentration']}.")
+        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"""
+        label = ' for ' + values['label'] if 'label' in values else ''
+        if isinstance(value, tuple):
+            if value[0] == value[1]:
+                raise ValueError(f"A concentration range cannot contain to identical concentrations.\n"
+                                 f"Got{label}: {value}.")
+            if value[0] > value[1]:
+                value = (value[1], value[0])
+        return value
+
+    @validator('balance')
+    def check_concentration_of_balance_species(cls, value, values):
+        """Make sure the concentration of the balance species is defined, default to 1"""
+        if value and 'concentration' in values:
+            if not isinstance(values['concentration'], (int, float)):
+                raise ValueError(f"The balance species concentration cannot be defined as a range, "
+                                 f"got: {values['concentration']}.")
+            values['concentration'] = values['concentration'] or 1
+        return value
+
 
 class RMGReactor(BaseModel):
     """
@@ -492,7 +523,7 @@ def check_model(cls, value):
     @validator('pdep')
     def check_pdep_only_if_gas_phase(cls, value, values):
         """RMG.pdep validator"""
-        if value is not None and values['reactors'] is not None:
+        if value is not None and 'reactors' in values and values['reactors'] is not None:
             reactor_types = set([reactor.type for reactor in values['reactors']])
             if value is not None and not any(['gas' in reactor for reactor in reactor_types]):
                 raise ValueError(f'A pdep section can only be specified for gas phase reactors, got: {reactor_types}')

t3/utils/writer.py

@@ -94,7 +94,7 @@ def write_rmg_input_file(rmg: dict,
 
     # reactors
     reactors = rmg['reactors']
-    gas_batch_constant_t_p_template_template = """
+    gas_batch_constant_t_p_template = """
 simpleReactor(
     temperature=${temperature},
     pressure=${pressure},
@@ -104,7 +104,12 @@ def write_rmg_input_file(rmg: dict,
 )
 <%def name="concentrations()">
 % for spc in species_list:
+    % if isinstance(spc["concentration"], (int, float)):
         '${spc["label"]}': ${spc["concentration"]},
+    % endif
+    % if isinstance(spc["concentration"], (tuple, list)):
+        '${spc["label"]}': [${spc["concentration"][0]}, ${spc["concentration"][1]}],
+    % endif
 % endfor
 </%def>
 """
@@ -117,7 +122,12 @@ def write_rmg_input_file(rmg: dict,
 )
 <%def name="concentrations()">
 % for spc in species_list:
+    % if isinstance(spc["concentration"], (int, float)):
         '${spc["label"]}': (${spc["concentration"]}, 'mol/cm^3'),
+    % endif
+    % if isinstance(spc["concentration"], (tuple, list)):
+        '${spc["label"]}': [(${spc["concentration"][0]}, 'mol/cm^3'), (${spc["concentration"][1]}, 'mol/cm^3')],
+    % endif
 % endfor
 </%def>
 """
@@ -129,8 +139,13 @@ def write_rmg_input_file(rmg: dict,
         else:
             raise ValueError(f"The reactor temperature must be a float or a list,\n"
                              f"got {reactor['T']} which is a {type(reactor['T'])}.")
-        species_list = [{'label': spc['label'], 'concentration': spc['concentration']} for spc in species]
-        species_list.sort(key=lambda spc: spc['concentration'], reverse=True)
+        species_list = [{'label': spc['label'], 'concentration': spc['concentration']} for spc in species
+                        if isinstance(spc['concentration'], (list, tuple))
+                        or spc['concentration'] > 0
+                        or spc['balance'] or not spc['reactive']]
+        species_list.sort(key=lambda spc: spc['concentration'][0] if isinstance(spc['concentration'], (tuple, list))
+                          else spc['concentration'], reverse=True)
+        print(species_list)
         termination = ''
         if reactor['termination_conversion'] is not None:
             termination += f"terminationConversion={reactor['termination_conversion']},"
@@ -161,7 +176,7 @@ def write_rmg_input_file(rmg: dict,
                 if spc['balance']:
                     balance = f"\n    balanceSpecies='{spc['label']}',"
                     break
-            rmg_input += Template(gas_batch_constant_t_p_template_template).render(
+            rmg_input += Template(gas_batch_constant_t_p_template).render(
                 temperature=temperature,
                 pressure=pressure,
                 species_list=species_list,

tests/test_main.py

@@ -519,7 +519,7 @@ def test_determine_species_to_calculate():
     assert len(list(t3.species.keys())) == 3
     assert all([species_dict['reasons'] == ['All core species'] for species_dict in t3.species.values()])
     assert all([species_dict['RMG label'] in ['OH', 'HO2', 'H2O2'] for species_dict in t3.species.values()])
-    assert all([species_dict['QM label'] in ['OH_0', 'HO2_1', 'H2O2_2'] for species_dict in t3.species.values()])
+    assert all([species_dict['QM label'] in ['0_OH', '1_HO2', '2_H2O2'] for species_dict in t3.species.values()])
 
     # 3. collision violators
     t3.iteration = 3
@@ -773,7 +773,7 @@ def test_add_species():
 
     found_h2 = False
     for qm_species in t3.qm['species']:
-        if qm_species.label == 'H2_4':
+        if qm_species.label == '4_H2':
             found_h2 = True
             assert isinstance(qm_species, ARCSpecies)
             assert qm_species.conformers == [{'symbols': ('H', 'H'),

tests/test_schema.py

@@ -241,7 +241,6 @@ def test_rmg_database_schema():
 def test_rmg_species_schema():
     """Test creating an instance of RMGSpecies"""
     rmg_species = RMGSpecies(label='N2',
-                             concentration=1,
                              smiles='N#N',
                              inchi='1S/N2/c1-2',
                              adjlist="""
@@ -258,7 +257,7 @@ def test_rmg_species_schema():
                              seed_all_rads=['radical', 'alkoxyl', 'peroxyl'],
                              )
     assert rmg_species.label == 'N2'
-    assert rmg_species.concentration == 1
+    assert rmg_species.concentration == 1  # automatically set if not specified
     assert rmg_species.smiles == 'N#N'
     assert rmg_species.inchi == '1S/N2/c1-2'
     # adjlist must be in the same column as adjlist above
@@ -275,6 +274,12 @@ def test_rmg_species_schema():
     assert rmg_species.solvent is False
     assert rmg_species.seed_all_rads == ['radical', 'alkoxyl', 'peroxyl']
 
+    rmg_species = RMGSpecies(label='H2O',
+                             concentration=[0.203, 0.502],
+                             smiles='O',
+                             )
+    assert rmg_species.concentration == (0.203, 0.502)
+
     with pytest.raises(ValidationError):
         # check that concentration is constrained to >= 0
         RMGSpecies(concentration=-1)
@@ -299,6 +304,50 @@ def test_rmg_species_schema():
                    seed_all_rads=['radical', 'non-supported'],
                    )
 
+    with pytest.raises(ValidationError):
+        # check that the concentration of a balance species is not a list
+        RMGSpecies(label='H2O',
+                   concentration=[0.203, 0.502],
+                   smiles='O',
+                   balance=True,
+                   )
+
+    with pytest.raises(ValidationError):
+        # check that concentration cannot be a 3-length tuple
+        RMGSpecies(label='H2O',
+                   concentration=[0.203, 0.502, 0.809],
+                   smiles='O',
+                   )
+
+    with pytest.raises(ValidationError):
+        # check that concentration cannot be negative
+        RMGSpecies(label='H2O',
+                   concentration=-0.203,
+                   smiles='O',
+                   )
+
+    with pytest.raises(ValidationError):
+        # check that a concentration range cannot include a negative number
+        RMGSpecies(label='H2O',
+                   concentration=[0.203, -0.502],
+                   smiles='O',
+                   )
+
+    with pytest.raises(ValidationError):
+        # check that a concentration range does not contain two equal boundaries
+        RMGSpecies(label='H2O',
+                   concentration=[0.203, 0.203],
+                   smiles='O',
+                   )
+
+    with pytest.raises(ValidationError):
+        # check that species defined with a concentration range cannot be constant
+        RMGSpecies(label='H2O',
+                   concentration=[0.203, 0.502],
+                   smiles='O',
+                   constant=True,
+                   )
+
 
 def test_rmg_reactors_schema():
     """Test creating an instance of RMGReactor"""

tests/test_writer.py

@@ -197,7 +197,7 @@ def test_write_rmg_input_file_liquid():
                         'concentration': 0.0124},
                        {'label': 'water',
                         'smiles': 'O',
-                        'concentration': 0.0278,
+                        'concentration': [0.0278, 0.0502],
                         'solvent': True},
                        {'label': 'O2',
                         'smiles': '[O][O]',
@@ -251,7 +251,7 @@ def test_write_rmg_input_file_liquid():
                  "liquidReactor(\n",
                  "    temperature=[(293.0, 'K'), (393.0, 'K')],\n",
                  "    initialConcentrations={\n",
-                 "        'water': (0.0278, 'mol/cm^3'),\n",
+                 "        'water': [(0.0278, 'mol/cm^3'), (0.0502, 'mol/cm^3')],\n",
                  "        'AIBN': (4.9e-06, 'mol/cm^3'),\n",
                  "        'O2': (2.73e-07, 'mol/cm^3'),\n",
                  "        'cyanoisopropylOO': (0, 'mol/cm^3'),\n",
@@ -281,7 +281,7 @@ def test_write_rmg_input_file_seed_all_radicals():
                         'seed_all_rads': ['radical', 'alkoxyl', 'peroxyl']},
                        {'label': 'O2',
                         'smiles': '[O][O]',
-                        'concentration': 2},
+                        'concentration': [2, 2.5]},
                        {'label': 'N2',
                         'smiles': 'N#N',
                         'constant': True,
@@ -311,7 +311,8 @@ def test_write_rmg_input_file_seed_all_radicals():
 
     with open(file_path, 'r') as f:
         lines = f.readlines()
-    for line in ["    thermoLibraries=['BurkeH2O2'],\n",
+    for line in ["        'O2': [2, 2.5],\n",
+                 "    thermoLibraries=['BurkeH2O2'],\n",
                  "    label='methylethylester',\n",
                  "    label='methylethylester_radical_0',\n",
                  "    label='methylethylester_alkoxyl_0',\n",