Merge pull request #385 from ReactionMechanismGenerator/choose_arrhen…

…ius_fit

Feature: Allow user specifying two- or three- parameter Arrhenius fit

arc/main.py

@@ -119,6 +119,9 @@ class ARC(object):
         compute_rates (bool, optional): Whether to compute rate coefficients for converged reactions.
         compute_transport (bool, optional): Whether to compute transport properties for converged species.
         statmech_adapter (str, optional): The statmech software to use.
+        three_params (bool, optional): Compute rate coefficients using the modified three-parameter Arrhenius equation
+                                       format (``True``, default) or classical two-parameter Arrhenius equation format
+                                       (``False``).
 
     Attributes:
         project (str): The project's name. Used for naming the working directory.
@@ -182,6 +185,8 @@ class ARC(object):
         statmech_adapter (str): The statmech software to use.
         fine_only (bool): If ``self.job_types['fine'] and not self.job_types['opt']`` ARC will not run optimization
                           jobs without fine=True
+        three_params (bool): Compute rate coefficients using the modified three-parameter Arrhenius equation
+                             format (``True``) or classical two-parameter Arrhenius equation format (``False``).
     """
 
     def __init__(self, input_dict=None, project=None, arc_species_list=None, arc_rxn_list=None, level_of_theory='',
@@ -192,7 +197,7 @@ def __init__(self, input_dict=None, project=None, arc_species_list=None, arc_rxn
                  job_memory=None, ess_settings=None, bath_gas=None, adaptive_levels=None, freq_scale_factor=None,
                  calc_freq_factor=True, n_confs=10, e_confs=5, dont_gen_confs=None, keep_checks=False,
                  solvation=None, compare_to_rmg=True, compute_thermo=True, compute_rates=True, compute_transport=True,
-                 specific_job_type='', statmech_adapter='Arkane'):
+                 specific_job_type='', statmech_adapter='Arkane', three_params=True):
         self.__version__ = VERSION
         self.verbose = verbose
         self.output = dict()
@@ -214,6 +219,7 @@ def __init__(self, input_dict=None, project=None, arc_species_list=None, arc_rxn
             self.project = project
             self.compute_thermo = compute_thermo
             self.compute_rates = compute_rates
+            self.three_params = three_params
             self.compute_transport = compute_transport
             self.statmech_adapter = statmech_adapter
             self.T_min = T_min
@@ -364,6 +370,8 @@ def as_dict(self) -> dict:
             restart_dict['compute_thermo'] = self.compute_thermo
         if not self.compute_rates:
             restart_dict['compute_rates'] = self.compute_rates
+        if not self.three_params:
+            restart_dict['three_params'] = self.three_params
         if not self.compute_transport:
             restart_dict['compute_transport'] = self.compute_transport
         restart_dict['statmech_adapter'] = self.statmech_adapter
@@ -443,6 +451,7 @@ def from_dict(self, input_dict, project=None, project_directory=None):
         self.execution_time = None
         self.compute_thermo = input_dict['compute_thermo'] if 'compute_thermo' in input_dict else True
         self.compute_rates = input_dict['compute_rates'] if 'compute_rates' in input_dict else True
+        self.three_params = input_dict['three_params'] if 'three_params' in input_dict else True
         self.compute_transport = input_dict['compute_transport'] if 'compute_transport' in input_dict else True
         self.statmech_adapter = input_dict['statmech_adapter'] if 'statmech_adapter' in input_dict else 'Arkane'
         self.verbose = input_dict['verbose'] if 'verbose' in input_dict else self.verbose
@@ -618,7 +627,9 @@ def execute(self) -> dict:
                             T_count=self.T_count or 50,
                             lib_long_desc=self.lib_long_desc,
                             rmg_database=self.rmg_database,
-                            compare_to_rmg=self.compare_to_rmg)
+                            compare_to_rmg=self.compare_to_rmg,
+                            three_params=self.three_params,
+                            )
 
         status_dict = self.summary()
         log_footer(execution_time=self.execution_time)

arc/mainTest.py

@@ -42,7 +42,8 @@ def test_as_dict(self):
         """Test the as_dict() method of ARC"""
         spc1 = ARCSpecies(label='spc1', smiles='CC', compute_thermo=False)
         arc0 = ARC(project='arc_test', job_types=self.job_types1, job_shortcut_keywords={'gaussian': 'scf=(NDump=30)'},
-                   arc_species_list=[spc1], level_of_theory='ccsd(t)-f12/cc-pvdz-f12//b3lyp/6-311+g(3df,2p)')
+                   arc_species_list=[spc1], level_of_theory='ccsd(t)-f12/cc-pvdz-f12//b3lyp/6-311+g(3df,2p)',
+                   three_params=False,)
         restart_dict = arc0.as_dict()
         long_thermo_description = restart_dict['species'][0]['long_thermo_description']
         self.assertIn('Bond corrections:', long_thermo_description)
@@ -114,6 +115,7 @@ def test_as_dict(self):
                                       't1': None}],
                          'specific_job_type': '',
                          'statmech_adapter': 'Arkane',
+                         'three_params': False,
                          }
         self.assertEqual(restart_dict, expected_dict)
 
@@ -152,7 +154,9 @@ def test_from_dict(self):
                                      'optical_isomers': 1,
                                      'rotors_dict': {},
                                      'xyzs': []}],
-                        'use_bac': True}
+                        'use_bac': True,
+                        'three_params': False,
+                        }
         arc1 = ARC(project='wrong', freq_scale_factor=0.95)
         self.assertEqual(arc1.freq_scale_factor, 0.95)  # user input
         project = 'arc_project_for_testing_delete_after_usage_test_from_dict'
@@ -168,6 +172,7 @@ def test_from_dict(self):
         self.assertEqual(arc1.arc_species_list[0].label, 'testing_spc1')
         self.assertFalse(arc1.arc_species_list[0].is_ts)
         self.assertEqual(arc1.arc_species_list[0].charge, 1)
+        self.assertFalse(arc1.three_params)
 
     def test_from_dict_specific_job(self):
         """Test the from_dict() method of ARC"""

arc/processor.py

@@ -49,6 +49,7 @@ def process_arc_project(statmech_adapter: str,
                         lib_long_desc: str = '',
                         rmg_database: Type[RMGDatabase] = None,
                         compare_to_rmg: bool = True,
+                        three_params: bool = True,
                         ) -> None:
     """
     Process an RMG project, generate thermo and rate coefficients using statistical mechanics (statmech).
@@ -75,6 +76,9 @@ def process_arc_project(statmech_adapter: str,
         rmg_database (RMGDatabase, optional): The RMG database object.
         compare_to_rmg (bool, optional): If ``True``, ARC's calculations will be compared against estimations
                                          from RMG's database.
+        three_params (bool, optional): Compute rate coefficients using the modified three-parameter Arrhenius equation
+                                       format (``True``, default) or classical two-parameter Arrhenius equation format
+                                       (``False``).
     """
     T_min = T_min or (300, 'K')
     T_max = T_max or (3000, 'K')
@@ -135,6 +139,7 @@ def process_arc_project(statmech_adapter: str,
                                                         T_min=T_min,
                                                         T_max=T_max,
                                                         T_count=T_count,
+                                                        three_params=three_params,
                                                         )
                     statmech_adapter.compute_high_p_rate_coefficient()
                     if reaction.kinetics is not None:

arc/statmech/arkane.py

@@ -51,6 +51,9 @@ class ArkaneAdapter(StatmechAdapter):
         T_min (tuple, optional): The minimum temperature for kinetics computations, e.g., (500, 'K').
         T_max (tuple, optional): The maximum temperature for kinetics computations, e.g., (3000, 'K').
         T_count (int, optional): The number of temperature points between t_min and t_max for kinetics computations.
+        three_params (bool, optional): Instruct Arkane to compute the high pressure kinetic rate coefficients in the
+                                       modified three-parameter Arrhenius equation format (``True``, default) or
+                                       classical two-parameter Arrhenius equation format (``False``).
     """
 
     def __init__(self,
@@ -64,7 +67,9 @@ def __init__(self,
                  species_dict: dict = None,
                  T_min: tuple = None,
                  T_max: tuple = None,
-                 T_count: int = 50):
+                 T_count: int = 50,
+                 three_params: bool = True,
+                 ):
         self.output_directory = output_directory
         self.output_dict = output_dict
         self.use_bac = use_bac
@@ -76,6 +81,7 @@ def __init__(self,
         self.T_min = T_min
         self.T_max = T_max
         self.T_count = T_count
+        self.three_params = three_params
 
         if not self.output_directory:
             raise InputError('A project directory was not provided.')
@@ -190,8 +196,13 @@ def compute_high_p_rate_coefficient(self) -> None:
                                              products=product_labels,
                                              transitionState=self.reaction.ts_label,
                                              tunneling='Eckart')
-                kinetics_job = KineticsJob(reaction=arkane_rxn, Tmin=self.T_min, Tmax=self.T_max, Tcount=self.T_count)
-                logger.info(f'Calculating rate for reaction {self.reaction.label}')
+                kinetics_job = KineticsJob(reaction=arkane_rxn, Tmin=self.T_min, Tmax=self.T_max, Tcount=self.T_count,
+                                           three_params=self.three_params)
+                if self.three_params:
+                    msg = 'using the modified three-parameter Arrhenius equation k = A * (T/T0)^n * exp(-Ea/RT)'
+                else:
+                    msg = 'using the classical two-parameter Arrhenius equation k = A * exp(-Ea/RT)'
+                logger.info(f'Calculating rate for reaction {self.reaction.label} {msg}.')
                 try:
                     kinetics_job.execute(output_directory=arkane_output_path, plot=True)
                 except (ValueError, OverflowError) as e:

arc/statmech/factory.py

@@ -45,6 +45,7 @@ def statmech_factory(statmech_adapter_label: str,  # add everything that goes in
                      T_min: tuple = None,
                      T_max: tuple = None,
                      T_count: int = 50,
+                     three_params: bool = True,
                      ) -> Type[StatmechAdapter]:
     """
     A factory generating a statmech adapter corresponding to ``statmech_adapter``.
@@ -64,6 +65,9 @@ def statmech_factory(statmech_adapter_label: str,  # add everything that goes in
         T_min (tuple, optional): The minimum temperature for kinetics computations, e.g., (500, 'K').
         T_max (tuple, optional): The maximum temperature for kinetics computations, e.g., (3000, 'K').
         T_count (int, optional): The number of temperature points between t_min and t_max for kinetics computations.
+        three_params (bool, optional): Compute rate coefficients using the modified three-parameter Arrhenius equation
+                                       format (``True``, default) or classical two-parameter Arrhenius equation format
+                                       (``False``).
 
     Returns:
         StatmechAdapter: The requested StatmechAdapter instance, initialized with the respective arguments,
@@ -79,5 +83,6 @@ def statmech_factory(statmech_adapter_label: str,  # add everything that goes in
                                                                                    T_min=T_min,
                                                                                    T_max=T_max,
                                                                                    T_count=T_count,
+                                                                                   three_params=three_params
                                                                                    )
     return statmech_adapter_class

docs/source/advanced.rst

@@ -838,5 +838,28 @@ one can supply an ARC job ID, and ALL jobs related to the project of the given j
 Note that either a ``-a``, a ``-p``, or a ``-j`` flag must be given.
 All flags can be combined with the optional ``-s`` flag.
 
+Choose modified/classical Arrhenius equation form for rate coefficient fitting
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+ARC uses statistical mechanics software packages (e.g., Arkane) to compute rate coefficients for chemical reactions from
+the results of quantum chemistry calculations. By default, ARC instructs the statmech programs to compute rate
+coefficients in the modified three-parameter Arrhenius equation format:
+
+  .. math :: k(T) = A \left( \frac{T}{T_0} \right)^n \exp \left( -\frac{E_a}{RT} \right)
+
+Alternatively, the user may request to compute the rate coefficients in the classical two-parameter Arrhenius format:
+
+  .. math :: k(T) = A \exp \left( -\frac{E_a}{RT} \right)
+
+by setting the ``three_params`` attribute to ``False`` (it is ``True`` by default). For example::
+
+    project: use_classical_arrhenius_eqn_for_rate_calc_demo
+    three_params: False
+
+instructs the relevant statmech program to compute rate coefficients in the classical two-parameter Arrhenius format for
+all reactions in the same ARC project.
+
+Advanced: to recompute the rate coefficient in the modified three-parameter Arrhenius equation format, simply change
+``three_params`` to ``True`` in the ARC project's restart.yml file, and then restart ARC.
 
 .. include:: links.txt