Collection of small fixes

.conda/meta.yaml

@@ -43,7 +43,6 @@ requirements:
     - lpsolve55
     - markupsafe
     - matplotlib >=1.5
-    - mock
     - mopac
     - mpmath
     - networkx

arkane/pdepTest.py

@@ -104,7 +104,9 @@ def test_pdep_job(self):
             reaction_list = read_reactions_block(chem, dictionary)
         rxn = reaction_list[0]
         self.assertIsInstance(rxn.kinetics, Chebyshev)
-        self.assertAlmostEquals(rxn.kinetics.get_rate_coefficient(1000.0, 1.0), 88.88253229631246)
+        # Accept a delta of 0.2, which could result from numerical discrepancies
+        # See RMG-Py #1682 on GitHub for discussion
+        self.assertAlmostEquals(rxn.kinetics.get_rate_coefficient(1000.0, 1.0), 88.88253229631246, delta=0.2)
 
         files = [f for f in os.listdir(os.path.join(self.directory, 'sensitivity', ''))
                  if os.path.isfile(os.path.join(self.directory, 'sensitivity', f))]
@@ -116,7 +118,7 @@ def test_pdep_job(self):
                 if '1000.0' in line:
                     break
         sa_coeff = line.split()[-2]
-        self.assertEquals(float(sa_coeff), -8.23e-6)
+        self.assertAlmostEquals(float(sa_coeff), -8.23e-6, delta=0.02e-6)
 
     @classmethod
     def tearDown(cls):

environment.yml

@@ -24,7 +24,6 @@ dependencies:
   - lpsolve55
   - markupsafe
   - matplotlib >=1.5
-  - mock
   - mopac
   - mpmath
   - networkx

rmgpy/chemkin.pyx

@@ -592,6 +592,7 @@ def _remove_line_breaks(comments):
                                 'This direction matched an entry in ', 'From training reaction',
                                 'This reaction matched rate rule', 'family: ', 'Warning:',
                                 'Chemkin file stated explicit reverse rate:', 'Ea raised from',
+                                'Fitted to', 'Reaction library',
                                 ]
     for indicator in new_statement_indicators:
         comments = comments.replace(' ' + indicator, '\n' + indicator, 1)

rmgpy/chemkinTest.py

@@ -27,9 +27,9 @@
 #                                                                             #
 ###############################################################################
 
-import unittest
-import mock
 import os
+import unittest
+from unittest import mock
 
 import rmgpy
 from rmgpy.chemkin import get_species_identifier, load_chemkin_file, load_transport_file, mark_duplicate_reactions, \

rmgpy/constraintsTest.py

@@ -32,7 +32,7 @@
 """
 
 import unittest
-import mock
+from unittest import mock
 
 from rmgpy.rmg.main import RMG
 from rmgpy.constraints import fails_species_constraints

rmgpy/data/kinetics/family.py

@@ -2758,63 +2758,6 @@ def add_atom_labels_for_reaction(self, reaction, output_with_resonance=True):
             for species in reaction.reactants + reaction.products:
                 species.generate_resonance_structures()
 
-    def get_w0(self, rxn):
-        """
-        calculates the w0 for Blower Masel kinetics by calculating wf (total bond energy of bonds formed)
-        and wb (total bond energy of bonds broken) with w0 = (wf+wb)/2
-        """
-        mol = None
-        a_dict = {}
-        for r in rxn.reactants:
-            m = r.molecule[0]
-            a_dict.update(m.get_all_labeled_atoms())
-            if mol:
-                mol = mol.merge(m)
-            else:
-                mol = m.copy(deep=True)
-
-        recipe = self.forward_recipe.actions
-
-        wb = 0.0
-        wf = 0.0
-        for act in recipe:
-
-            if act[0] == 'BREAK_BOND':
-                bd = mol.get_bond(a_dict[act[1]], a_dict[act[3]])
-                wb += bd.get_bde()
-            elif act[0] == 'FORM_BOND':
-                bd = Bond(a_dict[act[1]], a_dict[act[3]], act[2])
-                wf += bd.get_bde()
-            elif act[0] == 'CHANGE_BOND':
-                bd1 = mol.get_bond(a_dict[act[1]], a_dict[act[3]])
-
-                if act[2] + bd1.order == 0.5:
-                    mol2 = None
-                    for r in rxn.products:
-                        m = r.molecule[0]
-                        if mol2:
-                            mol2 = mol2.merge(m)
-                        else:
-                            mol2 = m.copy(deep=True)
-                    bd2 = mol2.get_bond(a_dict[act[1]], a_dict[act[3]])
-                else:
-                    bd2 = Bond(a_dict[act[1]], a_dict[act[3]], bd1.order + act[2])
-
-                if bd2.order == 0:
-                    bd2_bde = 0.0
-                else:
-                    bd2_bde = bd2.get_bde()
-                bde_diff = bd2_bde - bd1.get_bde()
-                if bde_diff > 0:
-                    wf += abs(bde_diff)
-                else:
-                    wb += abs(bde_diff)
-
-        return (wf + wb) / 2.0
-
-    def get_w0s(self, rxns):
-        return list(map(self.get_w0, rxns))
-
     def get_training_depository(self):
         """
         Returns the `training` depository from self.depositories

rmgpy/data/kinetics/familyTest.py

@@ -32,8 +32,8 @@
 import os.path
 import shutil
 import unittest
+from unittest import mock
 
-import mock
 import numpy as np
 
 from rmgpy import settings

rmgpy/data/thermoTest.py

@@ -82,19 +82,19 @@ class TestThermoDatabase(unittest.TestCase):
     """
 
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """A function that is run ONCE before all unit tests in this class."""
         global database
-        self.database = database.thermo
+        cls.database = database.thermo
 
-        self.databaseWithoutLibraries = ThermoDatabase()
-        self.databaseWithoutLibraries.load(os.path.join(settings['database.directory'], 'thermo'), libraries=[])
+        cls.databaseWithoutLibraries = ThermoDatabase()
+        cls.databaseWithoutLibraries.load(os.path.join(settings['database.directory'], 'thermo'), libraries=[])
 
         # Set up ML estimator
         models_path = os.path.join(settings['database.directory'], 'thermo', 'ml', 'main')
         hf298_path = os.path.join(models_path, 'hf298')
         s298_cp_path = os.path.join(models_path, 's298_cp')
-        self.ml_estimator = MLEstimator(hf298_path, s298_cp_path)
+        cls.ml_estimator = MLEstimator(hf298_path, s298_cp_path)
 
     def test_pickle(self):
         """
@@ -1640,10 +1640,10 @@ class TestThermoCentralDatabaseInterface(unittest.TestCase):
     """
 
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """A function that is run ONCE before all unit tests in this class."""
         global database
-        self.database = database.thermo
+        cls.database = database.thermo
 
     def connect_to_test_central_database(self):
         host, port, username, password = get_testing_tcd_authentication_info()

rmgpy/data/transportTest.py

@@ -126,13 +126,13 @@ class TestTransportDatabase(unittest.TestCase):
     """
 
     @classmethod
-    def setUpClass(self):
+    def setUpClass(cls):
         """A function that is run ONCE before all unit tests in this class."""
-        self.database = TransportDatabase()
-        self.database.load(os.path.join(settings['database.directory'], 'transport'),
-                           ['GRI-Mech', 'PrimaryTransportLibrary'])
+        cls.database = TransportDatabase()
+        cls.database.load(os.path.join(settings['database.directory'], 'transport'),
+                          ['GRI-Mech', 'PrimaryTransportLibrary'])
 
-        self.speciesList = [
+        cls.speciesList = [
             Species().from_smiles('C'),
             Species().from_smiles('CCCC'),
             Species().from_smiles('O'),

rmgpy/kinetics/arrhenius.pyx

@@ -1081,10 +1081,10 @@ def get_w0(actions, rxn):
     and wb (total bond energy of bonds broken) with w0 = (wf+wb)/2
     """
     mol = None
-    aDict = {}
+    a_dict = {}
     for r in rxn.reactants:
         m = r.molecule[0]
-        aDict.update(m.get_all_labeled_atoms())
+        a_dict.update(m.get_all_labeled_atoms())
         if mol:
             mol = mol.merge(m)
         else:
@@ -1097,13 +1097,13 @@ def get_w0(actions, rxn):
     for act in recipe:
 
         if act[0] == 'BREAK_BOND':
-            bd = mol.get_bond(aDict[act[1]], aDict[act[3]])
+            bd = mol.get_bond(a_dict[act[1]], a_dict[act[3]])
             wb += bd.get_bde()
         elif act[0] == 'FORM_BOND':
-            bd = Bond(aDict[act[1]], aDict[act[3]], act[2])
+            bd = Bond(a_dict[act[1]], a_dict[act[3]], act[2])
             wf += bd.get_bde()
         elif act[0] == 'CHANGE_BOND':
-            bd1 = mol.get_bond(aDict[act[1]], aDict[act[3]])
+            bd1 = mol.get_bond(a_dict[act[1]], a_dict[act[3]])
 
             if act[2] + bd1.order == 0.5:
                 mol2 = None
@@ -1113,19 +1113,19 @@ def get_w0(actions, rxn):
                         mol2 = mol2.merge(m)
                     else:
                         mol2 = m.copy(deep=True)
-                bd2 = mol2.get_bond(aDict[act[1]], aDict[act[3]])
+                bd2 = mol2.get_bond(a_dict[act[1]], a_dict[act[3]])
             else:
-                bd2 = Bond(aDict[act[1]], aDict[act[3]], bd1.order + act[2])
+                bd2 = Bond(a_dict[act[1]], a_dict[act[3]], bd1.order + act[2])
 
             if bd2.order == 0:
-                bd2bde = 0.0
+                bd2_bde = 0.0
             else:
-                bd2bde = bd2.get_bde()
-            bdediff = bd2bde - bd1.get_bde()
-            if bdediff > 0:
-                wf += abs(bdediff)
+                bd2_bde = bd2.get_bde()
+            bde_diff = bd2_bde - bd1.get_bde()
+            if bde_diff > 0:
+                wf += abs(bde_diff)
             else:
-                wb += abs(bdediff)
+                wb += abs(bde_diff)
     return (wf + wb) / 2.0
 
 def get_w0s(actions, rxns):

rmgpy/molecule/converter.py

@@ -226,7 +226,10 @@ def to_ob_mol(mol, return_mapping=False):
     obmol = openbabel.OBMol()
     for atom in atoms:
         a = obmol.NewAtom()
-        a.SetAtomicNum(atom.number)
+        if atom.element.symbol == 'X':
+            a.SetAtomicNum(78)  # not sure how to do this with linear scaling when this might not be Pt
+        else:
+            a.SetAtomicNum(atom.number)
         if atom.element.isotope != -1:
             a.SetIsotope(atom.element.isotope)
         a.SetFormalCharge(atom.charge)

rmgpy/molecule/translator.py

@@ -531,6 +531,8 @@ def _write(mol, identifier_type, backend):
             logging.debug('Backend {0} is not able to generate {1} for this molecule:\n'
                           '{2}'.format(option, identifier_type, mol.to_adjacency_list()))
 
+    logging.error('Unable to generate identifier for this molecule:\n{0}'.format(mol.to_adjacency_list()))
+
     raise ValueError("Unable to generate identifier type {0} with backend {1}.".format(identifier_type, backend))
 
 

rmgpy/molecule/translatorTest.py

@@ -34,6 +34,7 @@
 import re
 import unittest
 from external.wip import work_in_progress
+from unittest.mock import patch
 
 from rmgpy.molecule.adjlist import ConsistencyChecker
 from rmgpy.molecule.atomtype import ATOMTYPES
@@ -52,6 +53,18 @@ def test_empty_molecule(self):
         self.assertEqual(mol.to_smiles(), '')
         self.assertEqual(mol.to_inchi(), '')
 
+    @patch('rmgpy.molecule.translator.logging')
+    def test_failure_message(self, mock_logging):
+        """Test that we log the molecule adjlist upon failure."""
+        mol = Molecule(smiles='[CH2-][N+]#N')
+
+        with self.assertRaisesRegex(ValueError, 'Unable to generate identifier type'):
+            to_inchi(mol, backend='rdkit')
+
+        mock_logging.error.assert_called_with(
+            "Unable to generate identifier for this molecule:\n{0}".format(mol.to_adjacency_list())
+        )
+
 
 class InChIGenerationTest(unittest.TestCase):
     def compare(self, adjlist, aug_inchi):
@@ -414,6 +427,32 @@ def test_isotopic_molecule_2(self):
 
         self.assertEqual(mol.to_inchi(), inchi)
 
+    def test_surface_molecule_rdkit(self):
+        """Test InChI generation for a surface molecule using RDKit"""
+        mol = Molecule().from_adjacency_list("""
+1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+2 H u0 p0 c0 {1,S}
+3 H u0 p0 c0 {1,S}
+4 H u0 p0 c0 {1,S}
+5 X u0 p0 c0 {1,S}
+""")
+        inchi = 'InChI=1S/CH3.Pt/h1H3;'
+
+        self.assertEqual(to_inchi(mol, backend='rdkit'), inchi)
+
+    def test_surface_molecule_ob(self):
+        """Test InChI generation for a surface molecule using OpenBabel"""
+        mol = Molecule().from_adjacency_list("""
+1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+2 H u0 p0 c0 {1,S}
+3 H u0 p0 c0 {1,S}
+4 H u0 p0 c0 {1,S}
+5 X u0 p0 c0 {1,S}
+""")
+        inchi = 'InChI=1S/CH3.Pt/h1H3;'
+
+        self.assertEqual(to_inchi(mol, backend='openbabel'), inchi)
+
 
 class SMILESGenerationTest(unittest.TestCase):
     def compare(self, adjlist, smiles):
@@ -741,6 +780,32 @@ def test_aromatics(self):
         self.assertNotEqual(smiles2, smiles3)
         self.assertNotEqual(smiles1, smiles3)
 
+    def test_surface_molecule_rdkit(self):
+        """Test InChI generation for a surface molecule using RDKit"""
+        mol = Molecule().from_adjacency_list("""
+1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+2 H u0 p0 c0 {1,S}
+3 H u0 p0 c0 {1,S}
+4 H u0 p0 c0 {1,S}
+5 X u0 p0 c0 {1,S}
+""")
+        smiles = 'C[Pt]'
+
+        self.assertEqual(to_smiles(mol, backend='rdkit'), smiles)
+
+    def test_surface_molecule_ob(self):
+        """Test InChI generation for a surface molecule using OpenBabel"""
+        mol = Molecule().from_adjacency_list("""
+1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+2 H u0 p0 c0 {1,S}
+3 H u0 p0 c0 {1,S}
+4 H u0 p0 c0 {1,S}
+5 X u0 p0 c0 {1,S}
+""")
+        smiles = 'C[Pt]'
+
+        self.assertEqual(to_smiles(mol, backend='openbabel'), smiles)
+
 
 class ParsingTest(unittest.TestCase):
     def setUp(self):