Merge pull request #1810 from ReactionMechanismGenerator/arkane_bugfixes

Arkane bugfixes

arkane/gaussian.py

@@ -360,7 +360,7 @@ def load_scan_energies(self):
             line = f.readline()
             while line != '':
                 # If the job contains a "freq" then we want to ignore the last energy
-                if ' freq ' in line:
+                if ' Freq' in line and ' Geom=' in line:
                     opt_freq = True
                 # if # scan is keyword instead of # opt, then this is a rigid scan job
                 # and parsing the energies is done a little differently
@@ -404,7 +404,7 @@ def load_scan_energies(self):
 
         return vlist, angle
 
-    def _load_scan_specs(self, letter_spec):
+    def _load_scan_specs(self, letter_spec, get_after_letter_spec=False):
         """
         This method reads the ouptput file for optional parameters
         sent to gaussian, and returns the list of optional parameters
@@ -413,6 +413,10 @@ def _load_scan_specs(self, letter_spec):
         `letter_spec` is a character used to identify whether a specification
         defines pivot atoms ('S'), frozen atoms ('F') or other attributes.
 
+        `get_after_letter_spec` is a boolean that, if True, will return the
+        parameters after letter_spec is found. If not specified or False, it will
+        return the preceeding letters, which are typically the atom numbers.
+
         More information about the syntax can be found http://gaussian.com/opt/
         """
         output = []
@@ -437,7 +441,10 @@ def _load_scan_specs(self, letter_spec):
                     if len(terms) > action_index:
                         # specified type explicitly
                         if terms[action_index] == letter_spec:
-                            output.append(terms[1:action_index])
+                            if get_after_letter_spec:
+                                output.append(terms[action_index+1:])
+                            else:
+                                output.append(terms[1:action_index])
                     else:
                         # no specific specification, assume freezing
                         if letter_spec == 'F':
@@ -466,6 +473,20 @@ def load_scan_frozen_atoms(self):
         """
         return self._load_scan_specs('F')
 
+    def _load_scan_angle(self):
+        """
+        Return the angle difference (degrees) for a gaussian scan.
+        """
+        output = self._load_scan_specs('S', get_after_letter_spec=True)
+        return float(output[0][1])
+
+    def _load_number_scans(self):
+        """
+        Return the number of scans for a gaussian scan specified in the input file
+        """
+        output = self._load_scan_specs('S', get_after_letter_spec=True)
+        return int(output[0][0])
+
     def load_negative_frequency(self):
         """
         Return the negative frequency from a transition state frequency

arkane/gaussianTest.py

@@ -179,6 +179,20 @@ def test_load_symmetry_and_optics(self):
                 found_rotor = True
         self.assertTrue(found_rotor)
 
+    def test_load_scan_angle(self):
+        """
+        Ensures proper scan angle found in Gaussian scan job
+        """
+        log = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'isobutanolQOOH_scan.log'))
+        self.assertAlmostEqual(log._load_scan_angle(), 10.0)
+
+    def test_load_number_scans(self):
+        """
+        Ensures proper scan angle found in Gaussian scan job
+        """
+        log = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'isobutanolQOOH_scan.log'))
+        self.assertAlmostEqual(log._load_number_scans(), 36)
+
     def test_determine_qm_software(self):
         """
         Ensures that determine_qm_software returns a GaussianLog object
@@ -197,6 +211,16 @@ def test_gaussian_log_error_termination(self):
         self.assertTrue(f'The Gaussian job in {file_path} did not converge.' in str(log_error.exception))
 
 
+    def test_load_scan_with_freq(self):
+        """
+        Ensures that the length of enegies with hr scans and freq calc is correct
+        """
+        log = GaussianLog(os.path.join(os.path.dirname(__file__), 'data', 'hr_scan_with_freq.log'))
+        self.assertAlmostEqual(log._load_number_scans(), 36)
+        self.assertAlmostEqual(log._load_scan_angle(), 10.0)
+        vlist, _ = log.load_scan_energies()
+        self.assertEqual(len(vlist), 37)
+
 ################################################################################
 
 if __name__ == '__main__':

arkane/input.py

@@ -290,12 +290,13 @@ def reaction(label, reactants, products, transitionState=None, kinetics=None, tu
     products = sorted([species_dict[spec] for spec in products])
     if transitionState:
         transitionState = transition_state_dict[transitionState]
-    if tunneling.lower() == 'wigner':
+    if transitionState and (tunneling == '' or tunneling is None):
+        transitionState.tunneling = None
+    elif tunneling.lower() == 'wigner':
         transitionState.tunneling = Wigner(frequency=None)
     elif tunneling.lower() == 'eckart':
         transitionState.tunneling = Eckart(frequency=None, E0_reac=None, E0_TS=None, E0_prod=None)
-    elif transitionState and (tunneling == '' or tunneling is None):
-        transitionState.tunneling = None
+
     elif transitionState and not isinstance(tunneling, TunnelingModel):
         raise ValueError('Unknown tunneling model {0!r}.'.format(tunneling))
     rxn = Reaction(label=label, reactants=reactants, products=products, transition_state=transitionState,

arkane/statmech.py

@@ -318,7 +318,7 @@ def load(self, pdep=False, plot=False):
             except KeyError:
                 raise InputError('Model chemistry {0!r} not found in from dictionary of energy values in species file '
                                  '{1!r}.'.format(self.modelChemistry, path))
-        if not os.path.isfile(energy.path):
+        if isinstance(energy, Log) and not os.path.isfile(energy.path):
             modified_energy_path = os.path.join(directory, energy.path)
             if not os.path.isfile(modified_energy_path):
                 raise InputError('Could not find single point energy file for species {0} '
@@ -639,7 +639,7 @@ def load(self, pdep=False, plot=False):
 
         # save supporting information for calculation
         self.supporting_info = [self.species.label]
-        symmetry_read, optical_isomers_read, point_group_read = statmech_log.get_symmetry_properties()
+        optical_isomers_read, symmetry_read, point_group_read = statmech_log.get_symmetry_properties()
         self.supporting_info.append(external_symmetry if external_symmetry else symmetry_read)
         self.supporting_info.append(optical_isomers if optical_isomers else optical_isomers_read)
         self.supporting_info.append(point_group_read)

rmgpy/pdep/configuration.pyx

@@ -200,9 +200,12 @@ cdef class Configuration(object):
         # Evaluate configuration integral
         tred = constants.R * T / epsilon
         omega22 = 1.16145 * tred ** (-0.14874) + 0.52487 * exp(-0.77320 * tred) + 2.16178 * exp(-2.43787 * tred)
-
         # Evaluate collision frequency
-        return omega22 * sqrt(8 * constants.kB * T / constants.pi / mu) * constants.pi * sigma*sigma * gas_concentration
+        collision_freq = omega22 * sqrt(8 * constants.kB * T / constants.pi / mu) * constants.pi * sigma*sigma * gas_concentration
+        logging.debug("Obtained the following collision parameters: sigma {0}, epsilon {1}, "
+                      "omega22 {2}, collision rate {3}".format(sigma, epsilon, omega22, collision_freq))
+
+        return collision_freq
 
     cpdef np.ndarray generate_collision_matrix(self, double T, np.ndarray dens_states, np.ndarray e_list,
                                                np.ndarray j_list=None):
@@ -292,9 +295,11 @@ cdef class Configuration(object):
                         raise AttributeError('Length of masses should be three for termolecular reactants. '
                                              'We got {0}.'.format(len(mass)))
                     mu = 1.0 / (1.0 / mass[0] + 1.0 / mass[1])
-                    modes.insert(0, IdealGasTranslation(mass=(mu / constants.amu, "amu")))
-                    mu2 = 1.0 / (1.0 / mass[0] + 1.0 / mass[2])
-                    modes.insert(0, IdealGasTranslation(mass=(mu2 / constants.amu, "amu")))
+                    mu2 = 1.0 / (1.0 / mass[1] + 1.0 / mass[2])
+                    mu3 = 1.0 / (1.0 / mass[0] + 1.0 / mass[2])
+                    mu_avg = (mu + mu2 + mu3) / 3 * 1.1447142
+                    modes.insert(0, IdealGasTranslation(mass=(mu_avg / constants.amu, "amu")))
+                    modes.insert(0, IdealGasTranslation(mass=(mu_avg / constants.amu, "amu")))
             if rmgmode:
                 # Compute the density of states by direct count
                 # This is currently faster than the method of steepest descents,

rmgpy/pdep/configurationTest.py

@@ -31,6 +31,7 @@
 This script contains unit tests of the :mod:`rmgpy.pdep.network` module.
 """
 
+import numpy as np
 import unittest
 
 from rmgpy.pdep.collision import SingleExponentialDown
@@ -160,6 +161,14 @@ def test_str(self):
         for label in attributes:
             self.assertNotIn(label, output)
 
+    def test_no_nan_in_densStates(self):
+        """
+        This test asserts that there shouldn't be any NaN in the density of
+        states produced by calculateDensityofStates
+        """
+        elist = np.linspace(0, 1e5)
+        self.configuration.calculate_density_of_states(elist)
+        self.assertFalse(np.isnan(self.configuration.dens_states).any())
 
 ################################################################################
 

rmgpy/pdep/network.py

@@ -50,16 +50,16 @@ class Network(object):
     ======================= ====================================================
     Attribute               Description
     ======================= ====================================================
-    `isomers`               A list of the unimolecular isomers in the network
+    `isomers`               A list of the unimolecular isomers (Configuration objects) in the network
     `reactants`             A list of the bimolecular reactant channels (Configuration objects) in the network
     `products`              A list of the bimolecular product channels (Configuration objects) in the network
-    `path_reactions`        A list of "path" reaction objects that connect adjacent isomers (the high-pressure-limit)
+    `path_reactions`        A list of Reaction objects that connect isomers to their unimolecular and bimolecular products (the high-pressure-limit)
     `bath_gas`              A dictionary of the bath gas species (keys) and their mole fractions (values)
-    `net_reactions`         A list of "net" reaction objects that connect any pair of isomers
+    `net_reactions`         A list of Reaction objects that connect any pair of isomers (pressure dependent reactions)
     ----------------------- ----------------------------------------------------
     `T`                     The current temperature in K
-    `P`                     The current pressure in bar
-    `e_list`                The current array of energy grains in kJ/mol
+    `P`                     The current pressure in Pa
+    `e_list`                The current array of energy grains in J/mol
     `j_list`                The current array of total angular momentum quantum numbers
     ----------------------- ----------------------------------------------------
     `n_isom`                 The number of unimolecular isomers in the network
@@ -70,17 +70,23 @@ class Network(object):
     ----------------------- ----------------------------------------------------
     `grain_size`            Maximum size of separation between energies
     `grain_count`           Minimum number of descrete energies separated
-    `E0`                    A list of ground state energies of isomers, reactants, and products
+    `E0`                    A list of ground state energies of isomers, reactants, and products (J/mol)
     `active_k_rotor`        ``True`` if the K-rotor is treated as active, ``False`` if treated as adiabatic
     `active_j_rotor`        ``True`` if the J-rotor is treated as active, ``False`` if treated as adiabatic
     `rmgmode`               ``True`` if in RMG mode, ``False`` otherwise
     ----------------------- ----------------------------------------------------
     `eq_ratios`             An array containing concentration of each isomer and reactant channel present at equilibrium
-    `coll_freq`             An array of the frequency of collision between
+    `coll_freq`             An array of the frequency of collision between isomers and the bath gas
     `Mcoll`                 Matrix of first-order rate coefficients for collisional population transfer between grains for each isomer
     `dens_states`           3D np array of stable configurations, number of grains, and number of J
+    ----------------------- ----------------------------------------------------
+    `Kij`                   The microcanonical rates to go from isomer $j$ to isomer $i$. 4D array with indexes: i, j, energies, rotational energies
+    `Gnj`                   The microcanonical rates to go from isomer $j$ to reactant/product $n$. 4D array with indexes: n, j, energies, rotational energies
+    `Fim`                   The microcanonical rates to go from reactant $m$ to isomer $i$. 4D array with indexes: n, j, energies, rotational energies
+    ----------------------- ----------------------------------------------------
+    `K`                     2D Array of phenomenological rates at the specified T and P
+    `p0`                    Pseudo-steady state population distributions
     ======================= ====================================================
-    
     """
 
     def __init__(self, label='', isomers=None, reactants=None, products=None,
@@ -327,7 +333,7 @@ def set_conditions(self, T, P, ymB=None):
 
         n_isom = self.n_isom
         n_reac = self.n_reac
-        n_proc = self.n_prod
+        n_prod = self.n_prod
 
         E0 = self.E0
         grain_size = self.grain_size
@@ -391,12 +397,27 @@ def set_conditions(self, T, P, ymB=None):
 
                 # Rescale densities of states such that, when they are integrated
                 # using the Boltzmann factor as a weighting factor, the result is unity
-                for i in range(n_isom + n_reac):
+                # this converts the denisty of states into a population distribution
+                # as described in Allen 2012, equation 1
+                # check for numerical errors in Boltzman distribution
+                if np.exp(-self.e_list.max() / constants.R / self.T) == 0.0 or\
+                   np.exp(-self.e_list.min() / constants.R / self.T) * self.dens_states.max() == np.inf:
+                    logging.warning("The energies of range ({0}, {1}) J/mol result in numerical rounding errors in the Bolzman distribution. "
+                                    "Check that your energy corrections are accurate.".format(self.e_list.min(), self.e_list.max()))
+                self.dens_states_raw = self.dens_states.copy()
+                for i in range(n_isom + n_reac + n_prod):
                     Q = 0.0
                     for s in range(n_j):
                         Q += np.sum(
                             self.dens_states[i, :, s] * (2 * j_list[s] + 1) * np.exp(-e_list / constants.R / T))
-                    self.dens_states[i, :, :] /= Q
+                    if Q == 0.:
+                        logging.warning('No density of states found for structure {1} in network {0}. '
+                                        'possibly a product witout any thermo.'.format(self.label, i))
+                    else:
+                        self.dens_states[i, :, :] /= Q
+                if np.isnan(self.dens_states).any():
+                    logging.warning('Density of states for network {0} has NaN values after '
+                                    'rescaling densities. Double check issues with network.'.format(self.label))
 
             # Update parameters that depend on temperature and pressure if necessary
             if temperature_changed or pressure_changed:
@@ -581,10 +602,8 @@ def map_densities_of_states(self):
                 logging.debug('Mapping density of states for product channel "{0}"'.format(self.products[n]))
                 self.dens_states[n + n_isom + n_reac, :, :] = self.products[n].map_density_of_states(self.e_list, self.j_list)
 
-        # import pylab
-        # for i in range(n_isom + n_reac + n_prod):
-        #    pylab.semilogy(self.e_list*0.001, self.dens_states[i,:])
-        # pylab.show()
+        if np.isnan(self.dens_states).any():
+            raise Exception('Density of states has a NaN value.\n{0}'.format(self.dens_states))
 
     def calculate_microcanonical_rates(self):
         """

rmgpy/pdep/networkTest.py

@@ -266,6 +266,16 @@ def test_collision_matrix_memory_handling(self):
         except:
             pass
 
+    def test_get_all_species(self):
+        """
+        Ensures all species are in the get_species_list
+        """
+        species_list = self.network.get_all_species()
+        self.assertIn(self.nC4H10O, species_list)
+        self.assertIn(self.nC4H8, species_list)
+        self.assertIn(self.H2O, species_list)
+        self.assertIn(self.N2, species_list)
+
 
 ################################################################################
 

rmgpy/quantity.py

@@ -748,7 +748,7 @@ def __call__(self, *args, **kwargs):
     'J/mol',
     common_units=['kJ/mol', 'cal/mol', 'kcal/mol', 'eV/molecule'],
     extra_dimensionality={'K': constants.R,
-                         'cm^-1': constants.h * constants.c * 100
+                         'cm^-1': constants.h * constants.c * 100 * constants.Na
                          # the following hack also allows 'J' and 'kJ' etc. to be specified without /mol[ecule]
                          # so is not advisable (and fails unit tests)
                          # 'eV': constants.Na, # allow people to be lazy and neglect the "/molecule"

rmgpy/reaction.py

@@ -1264,7 +1264,7 @@ def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
 
     def calculate_coll_limit(self, temp, reverse=False):
         """
-        Calculate the collision limit rate for the given temperature
+        Calculate the collision limit rate in m3/mol-s for the given temperature
         implemented as recommended in Wang et al. doi 10.1016/j.combustflame.2017.08.005 (Eq. 1)
         """
         reduced_mass = self.get_reduced_mass(reverse)