Merge pull request #96 from ReactionMechanismGenerator/xyz_2d

Fix to xyz/2D loading and atom reordering in ARCSpecies

arc/species/species.py

@@ -212,13 +212,12 @@ def __init__(self, is_ts=False, rmg_species=None, mol=None, label=None, xyz=None
                 self.multiplicity = multiplicity
                 self.charge = charge
 
-            if self.mol is None and not self.is_ts:
-                xyz = self.final_xyz if self.final_xyz else self.initial_xyz
-                _, mol = molecules_from_xyz(xyz)
-                self.number_of_atoms = len(mol.atoms)
-            elif not self.is_ts:
-                if self.initial_xyz is not None:
+            if not self.is_ts:
+                # Perceive molecule from xyz coordinates
+                # This also populates mol_list
+                if self.final_xyz or self.initial_xyz:
                     self.mol_from_xyz()
+                # Generate bond list for applying bond corrections
                 if not self.bond_corrections:
                     self.bond_corrections = self.mol.enumerate_bonds()
                     if self.bond_corrections:
@@ -232,6 +231,7 @@ def __init__(self, is_ts=False, rmg_species=None, mol=None, label=None, xyz=None
             elif not self.bond_corrections and self.generate_thermo:
                 logging.warn('Cannot determine bond additivity corrections (BAC) for species {0} based on xyz'
                              ' coordinates only. For better thermoproperties, provide bond corrections.')
+
             if self.initial_xyz is not None:
                 # consider the initial guess as one of the conformers if generating others.
                 # otherwise, just consider it as the conformer.
@@ -800,24 +800,33 @@ def mol_from_xyz(self, xyz=None):
         """
         Make sure atom order in self.mol corresponds to xyz
         Important for TS discovery and for identifying rotor indices
+
+        This works by generating a molecule from the xyz and using the
+        2D structure to confirm that the perceived molecule is correct.
+        Resonance structures are generated and saved to ``self.mol_list``.
         """
         if xyz is None:
-            xyz = self.initial_xyz
-        original_mol = self.mol.copy(deep=True)
-        _, self.mol = molecules_from_xyz(xyz)
-        if self.mol is not None and not check_isomorphism(original_mol, self.mol):
-            raise InputError('XYZ and the 2D graph representation of the Molecule are not isomorphic.\n'
-                             'Got xyz:\n{0}\n\nwhich corresponds to {1}\n{2}\n\nand: {3}\n{4}'.format(
-                               xyz, self.mol.toSMILES(), self.mol.toAdjacencyList(),
-                               original_mol.toSMILES(), original_mol.toAdjacencyList()))
+            xyz = self.final_xyz or self.initial_xyz
+
+        if self.mol is not None:
+            # self.mol should have come from another source, e.g. SMILES or yml
+            original_mol = self.mol
+            self.mol = molecules_from_xyz(xyz)[1]
+
+            if self.mol is not None and not check_isomorphism(original_mol, self.mol):
+                raise InputError('XYZ and the 2D graph representation of the Molecule are not isomorphic.\n'
+                                 'Got xyz:\n{0}\n\nwhich corresponds to {1}\n{2}\n\nand: {3}\n{4}'.format(
+                                   xyz, self.mol.toSMILES(), self.mol.toAdjacencyList(),
+                                   original_mol.toSMILES(), original_mol.toAdjacencyList()))
+        else:
+            self.mol = molecules_from_xyz(xyz)[1]
+
         if self.mol_list is None:
+            # Assign atom ids first, so they carry through to the resonance structures
             self.mol.assignAtomIDs()
-            self.mol_list = self.mol.generate_resonance_structures(keep_isomorphic=False, filter_structures=True)
-            # The generate_resonance_structures methods changes atom order, reorder atoms in self.mol:
-            id_mol = self.mol.copy(deep=True)
-            _, self.mol = molecules_from_xyz(xyz)
-            for i, atom in enumerate(self.mol.atoms):
-                atom.id = id_mol.atoms[i].id
+            # The generate_resonance_structures method changes atom order
+            # Make a copy so we don't disturb the original order from xyz
+            self.mol_list = self.mol.copy(deep=True).generate_resonance_structures(keep_isomorphic=False, filter_structures=True)
         order_atoms_in_mol_list(ref_mol=self.mol, mol_list=self.mol_list)
 
 

arc/species/speciesTest.py

@@ -70,6 +70,13 @@ def setUpClass(cls):
         xyz1 = os.path.join(arc_path, 'arc', 'testing', 'xyz', 'AIBN.gjf')
         cls.spc7 = ARCSpecies(label='AIBN', smiles=str('N#CC(C)(C)N=NC(C)(C)C#N'), xyz=xyz1)
 
+        hso3_xyz = str("""S      -0.12383700    0.10918200   -0.21334200
+        O       0.97332200   -0.98800100    0.31790100
+        O      -1.41608500   -0.43976300    0.14487300
+        O       0.32370100    1.42850400    0.21585900
+        H       1.84477700   -0.57224200    0.35517700""")
+        cls.spc8 = ARCSpecies(label=str('HSO3'), xyz=hso3_xyz, multiplicity=2, charge=0, smiles=str('O[S](=O)=O'))
+
     def test_conformers(self):
         """Test conformer generation"""
         self.spc1.generate_conformers()  # vinoxy has two res. structures, each is assigned two conformers (RDkit/ob)
@@ -454,6 +461,42 @@ def test_get_min_energy_conformer(self):
         min_xyz = get_min_energy_conformer(xyzs, energies)
         self.assertEqual(min_xyz, 'xyz2')
 
+    def test_mol_from_xyz_atom_id_1(self):
+        """Test that atom ids are saved properly when loading both xyz and smiles."""
+        mol = self.spc6.mol
+        mol_list = self.spc6.mol_list
+
+        self.assertEqual(len(mol_list), 1)
+        res = mol_list[0]
+
+        self.assertTrue(mol.atomIDValid())
+        self.assertTrue(res.atomIDValid())
+
+        self.assertTrue(mol.isIsomorphic(res))
+        self.assertTrue(mol.isIdentical(res))
+
+    def test_mol_from_xyz_atom_id_2(self):
+        """Test that atom ids are saved properly when loading both xyz and smiles."""
+        mol = self.spc8.mol
+        mol_list = self.spc8.mol_list
+
+        self.assertEqual(len(mol_list), 2)
+        res1, res2 = mol_list
+
+        self.assertTrue(mol.atomIDValid())
+        self.assertTrue(res1.atomIDValid())
+        self.assertTrue(res2.atomIDValid())
+
+        self.assertTrue(mol.isIsomorphic(res1))
+        self.assertTrue(mol.isIdentical(res1))
+
+        # Check that atom ordering is consistent, ignoring specific oxygen ordering
+        mol_ids = [(a.element.symbol, a.id) if a.element.symbol != 'O' else (a.element.symbol,) for a in mol.atoms]
+        res1_ids = [(a.element.symbol, a.id) if a.element.symbol != 'O' else (a.element.symbol,) for a in res1.atoms]
+        res2_ids = [(a.element.symbol, a.id) if a.element.symbol != 'O' else (a.element.symbol,) for a in res2.atoms]
+        self.assertEqual(mol_ids, res1_ids)
+        self.assertEqual(mol_ids, res2_ids)
+
 
 class TestTSGuess(unittest.TestCase):
     """