RMG-Py v2.2.1 release

Author: mliu49

documentation/source/users/rmg/releaseNotes.rst

@@ -4,6 +4,22 @@
 Release Notes
 *************
 
+RMG-Py Version 2.2.1
+====================
+Date July 23, 2018
+
+This release is minor patch which fixes a number of issues discovered after 2.2.0.
+
+- Collision limit checking:
+    - RMG will now output a list of collision limit violations for the generated model
+
+- Fixes:
+    - Ambiguous chemical formulas in SMILES lookup leading to incorrect SMILES generation
+    - Fixed issue with reading geometries from QChem output files
+    - React flags for reaction filter were not properly updated on each iteration
+    - Fixed issue with inconsistent symmetry number calculation
+
+
 RMG-Py Version 2.2.0
 ====================
 Date: July 5, 2018

rmgpy/cantherm/qchem.py

@@ -33,6 +33,7 @@
 import logging
 import os.path
 import rmgpy.constants as constants
+from rmgpy.exceptions import InputError
 from rmgpy.cantherm.common import checkConformerEnergy
 from rmgpy.statmech import IdealGasTranslation, NonlinearRotor, LinearRotor, HarmonicOscillator, Conformer
 ################################################################################
@@ -115,35 +116,41 @@ def loadGeometry(self):
         """
         atom = []; coord = []; number = []; 
 
-        f = open(self.path, 'r')
-        line = f.readline()
-        while line != '':
-            if 'Final energy is' in line:
-                print 'found a sucessfully completed Qchem Geometry Optimization Job'
-                line = f.readline()
-                atom = []; coord = []
+
+        with open(self.path) as f:
+            log = f.read().splitlines()
+
+        #First check that the Qchem job file (not necessarily a geometry optimization)
+        #has successfully completed, if not an error is thrown
+        completed_job = False
+        for line in reversed(log):
+            if 'Total job time:' in line:
+                logging.debug('Found a sucessfully completed Qchem Job')
+                completed_job = True
                 break
-            line = f.readline()
-        found = 0           
-        while line != '':        
+
+        if not completed_job:
+            raise InputError('Could not find a successfully completed Qchem job in Qchem output file {0}'.format(self.path))
+
+        #Now look for the geometry.
+        #Will return the final geometry in the file under Standard Nuclear Orientation.
+        geometry_flag = False
+        for i in reversed(xrange(len(log))):
+            line = log[i]
             if 'Standard Nuclear Orientation' in line:
-                found += 1
-                for i in range(3): line = f.readline() # skip  lines
-                while '----------------------------------------------------' not in line:
-                    data = line.split()
-                    atom.append((data[1]))
-                    coord.append([float(data[2]), float(data[3]), float(data[4])])
-                    line = f.readline()
-                # Read the next line in the file    
-                line = f.readline()
-            # Read the next line in the file
-            line = f.readline()
-            if found ==1: break
-        line = f.readline()
-        #print coord
-        f.close()
+                atom, coord, number = [], [], []
+                for line in log[(i+3):]:
+                    if '------------' not in line:
+                        data = line.split()
+                        atom.append(data[1])
+                        coord.append([float(c) for c in data [2:]])
+                        geometry_flag = True
+                    else:
+                        break
+                if geometry_flag:
+                    break
+
         coord = numpy.array(coord, numpy.float64)
-        mass = numpy.array(coord, numpy.float64)
         # Assign appropriate mass to each atom in molecule
         # These values were taken from "Atomic Weights and Isotopic Compositions" v3.0 (July 2010) from NIST
 
@@ -173,7 +180,9 @@ def loadGeometry(self):
                 number.append('17')
             else:
                 raise NotImplementedError('Atomic atom {0:d} not yet supported in loadGeometry().'.format(atom[i]))
-        number = numpy.array(number, numpy.int)       
+        number = numpy.array(number, numpy.int)
+        if len(number) == 0 or len(coord) == 0 or len(mass) == 0:
+            raise InputError('Unable to read the numbers and types of atoms from Qchem output file {0}'.format(self.path))
         return coord, number, mass
 
     def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, symfromlog=None, label=''):
@@ -301,7 +310,7 @@ def loadEnergy(self,frequencyScaleFactor=1.):
         if E0 is not None:
             return E0
         else:
-            raise Exception('Unable to find energy in Qchem output file.')
+            raise InputError('Unable to find energy in Qchem output file.')
 
     def loadZeroPointEnergy(self,frequencyScaleFactor=1.):
         """
@@ -331,7 +340,7 @@ def loadZeroPointEnergy(self,frequencyScaleFactor=1.):
         if ZPE is not None:
             return ZPE
         else:
-            raise Exception('Unable to find zero-point energy in Qchem output file.')
+            raise InputError('Unable to find zero-point energy in Qchem output file.')
 
     def loadScanEnergies(self):
         """
@@ -394,4 +403,4 @@ def loadNegativeFrequency(self):
         if frequency < 0:
             return frequency
         else:
-            raise Exception('Unable to find imaginary frequency in QChem output file {0}'.format(self.path))
+            raise InputError('Unable to find imaginary frequency in QChem output file {0}'.format(self.path))

rmgpy/data/thermoTest.py

@@ -424,7 +424,6 @@ def testNewThermoGeneration(self):
                 self.assertAlmostEqual(Cp, thermoData.getHeatCapacity(T) / 4.184, places=1,
                                        msg="Cp{3} error for {0}. Expected {1} but calculated {2}.".format(smiles, Cp, thermoData.getHeatCapacity(T) / 4.184, T))
 
-    @work_in_progress
     def testSymmetryNumberGeneration(self):
         """
         Test we generate symmetry numbers correctly.

rmgpy/molecule/vf2.pyx

@@ -166,6 +166,14 @@ cdef class VF2:
             graph1.restore_vertex_order()
             graph2.restore_vertex_order()
 
+        # We're done, so clear the mappings to prevent downstream effects
+        for vertex1 in graph1.vertices:
+            vertex1.mapping = None
+            vertex1.terminal = False
+        for vertex2 in graph2.vertices:
+            vertex2.mapping = None
+            vertex2.terminal = False
+
     cdef bint match(self, int callDepth) except -2:
         """
         Recursively search for pairs of vertices to match, until all vertices

rmgpy/molecule/vf2Test.py

@@ -81,6 +81,15 @@ def test_feasible(self):
                     self.assertTrue(self.vf2.feasible(atom1, atom2))
                 else: # different connectivity values should return false
                     self.assertFalse(self.vf2.feasible(atom1, atom2))
+
+    def test_clear_mapping(self):
+        """Test that vertex mapping is cleared after isomorphism."""
+        self.vf2.isIsomorphic(self.mol, self.mol2, None)
+
+        for atom in self.mol.atoms:
+            self.assertIsNone(atom.mapping)
+            self.assertFalse(atom.terminal)
+
 ################################################################################
 
 if __name__ == '__main__':

rmgpy/reaction.pxd

@@ -113,4 +113,12 @@ cdef class Reaction:
 
     cpdef ensure_species(self, bint reactant_resonance=?, bint product_resonance=?)
 
+    cpdef list check_collision_limit_violation(self, float t_min, float t_max, float p_min, float p_max)
+
+    cpdef calculate_coll_limit(self, float temp, bint reverse=?)
+
+    cpdef get_reduced_mass(self, bint reverse=?)
+
+    cpdef get_mean_sigma_and_epsilon(self, bint reverse=?)
+
 cpdef bint isomorphic_species_lists(list list1, list list2, bint check_identical=?, bint only_check_label=?)

rmgpy/reaction.py

@@ -1146,6 +1146,106 @@ def ensure_species(self, reactant_resonance=False, product_resonance=True):
         except AttributeError:
             pass
 
+    def check_collision_limit_violation(self, t_min, t_max, p_min, p_max):
+        """
+        Warn if a core reaction violates the collision limit rate in either the forward or reverse direction
+        at the relevant extreme T/P conditions. Assuming a monotonic behaviour of the kinetics.
+        Returns a list with the reaction object and the direction in which the violation was detected.
+        """
+        conditions = [[t_min, p_min]]
+        if t_min != t_max:
+            conditions.append([t_max, p_min])
+        if self.kinetics.isPressureDependent() and p_max != p_min:
+            conditions.append([t_min, p_max])
+            if t_min != t_max:
+                conditions.append([t_max, p_max])
+        logging.debug("Checking whether reaction {0} violates the collision rate limit...".format(self))
+        violator_list = []
+        kf_list = []
+        kr_list = []
+        collision_limit_f = []
+        collision_limit_r = []
+        for condition in conditions:
+            if len(self.reactants) >= 2:
+                try:
+                    collision_limit_f.append(self.calculate_coll_limit(temp=condition[0], reverse=False))
+                except ValueError:
+                    continue
+                else:
+                    kf_list.append(self.getRateCoefficient(condition[0], condition[1]))
+            if len(self.products) >= 2:
+                try:
+                    collision_limit_r.append(self.calculate_coll_limit(temp=condition[0], reverse=True))
+                except ValueError:
+                    continue
+                else:
+                    kr_list.append(self.generateReverseRateCoefficient().getRateCoefficient(condition[0], condition[1]))
+        if len(self.reactants) >= 2:
+            for i, k in enumerate(kf_list):
+                if k > collision_limit_f[i]:
+                    ratio = k / collision_limit_f[i]
+                    condition = '{0} K, {1:.1f} bar'.format(conditions[i][0], conditions[i][1]/1e5)
+                    violator_list.append([self, 'forward', ratio, condition])
+        if len(self.products) >= 2:
+            for i, k in enumerate(kr_list):
+                if k > collision_limit_r[i]:
+                    ratio = k / collision_limit_r[i]
+                    condition = '{0} K, {1:.1f} bar'.format(conditions[i][0], conditions[i][1]/1e5)
+                    violator_list.append([self, 'reverse', ratio, condition])
+        return violator_list
+
+    def calculate_coll_limit(self, temp, reverse=False):
+        """
+        Calculate the collision limit rate 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)
+        sigma, epsilon = self.get_mean_sigma_and_epsilon(reverse)
+        Tr = temp * constants.kB * constants.Na / epsilon
+        reduced_coll_integral = 1.16145 * Tr ** (-0.14874) + 0.52487 * math.exp(-0.7732 * Tr) + 2.16178 * math.exp(
+            -2.437887 * Tr)
+        k_coll = (math.sqrt(8 * math.pi * constants.kB * temp * constants.Na / reduced_mass) * sigma ** 2
+                  * reduced_coll_integral * constants.Na)
+        return k_coll
+
+    def get_reduced_mass(self, reverse=False):
+        """
+        Returns the reduced mass of the reactants if reverse is ``False``
+        Returns the reduced mass of the products if reverse is ``True``
+        """
+        if reverse:
+            mass_list = [spc.molecule[0].getMolecularWeight() for spc in self.products]
+        else:
+            mass_list = [spc.molecule[0].getMolecularWeight() for spc in self.reactants]
+        reduced_mass = reduce((lambda x, y: x * y), mass_list) / sum(mass_list)
+        return reduced_mass
+
+    def get_mean_sigma_and_epsilon(self, reverse=False):
+        """
+        Calculates the collision diameter (sigma) using an arithmetic mean
+        Calculates the well depth (epsilon) using a geometric mean
+        If reverse is ``False`` the above is calculated for the reactants, otherwise for the products
+        """
+        sigmas = []
+        epsilons = []
+        if reverse:
+            for spc in self.products:
+                trans = spc.getTransportData()
+                sigmas.append(trans.sigma.value_si)
+                epsilons.append(trans.epsilon.value_si)
+            num_of_spcs = len(self.products)
+        else:
+            for spc in self.reactants:
+                trans = spc.getTransportData()
+                sigmas.append(trans.sigma.value_si)
+                epsilons.append(trans.epsilon.value_si)
+            num_of_spcs = len(self.reactants)
+        if any([x == 0 for x in sigmas + epsilons]):
+            raise ValueError
+        mean_sigmas = sum(sigmas) / num_of_spcs
+        mean_epsilons = reduce((lambda x, y: x * y), epsilons) ** (1 / len(epsilons))
+        return mean_sigmas, mean_epsilons
+
 def isomorphic_species_lists(list1, list2, check_identical=False, only_check_label=False):
     """
     This method compares whether lists of species or molecules are isomorphic