added tests for hill formula, is_isomorphic_with and are_isomorphic, and remaping tests. Also tweaked the toolkits to get the atom mapping if found and attach it to the moelecule properties dict.

Author: joshhorton

openforcefield/tests/test_molecule.py

@@ -58,7 +58,7 @@ def assert_molecule_is_equal(molecule1, molecule2, msg):
         Message to include if molecules fail to match.
 
     """
-    if not(molecule1.is_isomorphic(molecule2)):
+    if not(molecule1.is_isomorphic_with(molecule2)):
         raise AssertionError(msg)
 
 
@@ -409,7 +409,7 @@ def test_to_from_iupac(self, molecule):
                 Molecule.from_iupac(iupac)
 
         molecule_copy = Molecule.from_iupac(iupac, allow_undefined_stereo=undefined_stereo)
-        assert molecule.is_isomorphic(molecule_copy, compare_atom_stereochemistry=not undefined_stereo)
+        assert molecule.is_isomorphic_with(molecule_copy)
 
     @pytest.mark.parametrize('molecule', mini_drug_bank())
     def test_to_from_topology(self, molecule):
@@ -529,6 +529,109 @@ def test_name(self):
         molecule.name = name
         assert molecule.name == name
 
+    def test_hill_formula(self):
+        """Test that making the hill formula is consistent between input methods and ordering"""
+        # make sure smiles match reference
+        molecule_smiles = Molecule.from_smiles('CCO')
+        assert molecule_smiles.hill_formula == 'C2H6O1'
+        # make sure is not order dependent
+        molecule_smiles_reverse = Molecule.from_smiles('OCC')
+        assert molecule_smiles.hill_formula == molecule_smiles_reverse.hill_formula
+        # make sure files and smiles match
+        molecule_file = Molecule.from_file(get_data_file_path('molecules/ethanol.mol2'))
+        assert molecule_smiles.hill_formula == molecule_file.hill_formula
+        # make sure the topology molecule gives the same formula
+        from openforcefield.topology.topology import TopologyMolecule, Topology
+        topology = Topology.from_molecules(molecule_smiles)
+        topmol = TopologyMolecule(molecule_smiles, topology)
+        assert molecule_smiles.hill_formula == Molecule.to_hill_formula(topmol)
+        # make sure the networkx matches
+        assert molecule_smiles.hill_formula == Molecule.to_hill_formula(molecule_smiles.to_networkx())
+
+    def test_isomorphic(self):
+        """Test the different levels of isomorphic matching"""
+        # check that hill formula fails are caught
+        ethanol = Molecule.from_smiles('CCO')
+        acetaldehyde = Molecule.from_smiles('CC=O')
+        assert ethanol.is_isomorphic_with(acetaldehyde) is False
+        assert acetaldehyde.is_isomorphic_with(ethanol) is False
+        # check that different orderings work with full matching
+        ethanol_reverse = Molecule.from_smiles('OCC')
+        assert ethanol.is_isomorphic_with(ethanol_reverse) is True
+        # check a reference mapping between ethanol and ethanol_reverse matches that calculated
+        ref_mapping = {0: 2, 1: 1, 2: 0, 3: 6, 4: 7, 5: 8, 6: 4, 7: 5, 8: 3}
+        assert Molecule.are_isomorphic(ethanol, ethanol_reverse, return_atom_map=True)[1] == ref_mapping
+        # check matching with nx.Graph atomic numbers and connectivity only
+        assert Molecule.are_isomorphic(ethanol, ethanol_reverse.to_networkx(), aromatic_matching=False,
+                                       formal_charge_matching=False, bond_order_matching=False,
+                                       atom_stereochemistry_matching=False,
+                                       bond_stereochemistry_matching=False)[0] is True
+        # check matching with nx.Graph with full matching
+        assert ethanol.is_isomorphic_with(ethanol_reverse.to_networkx()) is True
+        # check matching with a TopologyMolecule class
+        from openforcefield.topology.topology import TopologyMolecule, Topology
+        topology = Topology.from_molecules(ethanol)
+        topmol = TopologyMolecule(ethanol, topology)
+        assert Molecule.are_isomorphic(ethanol, topmol, aromatic_matching=False, formal_charge_matching=False,
+                                       bond_order_matching=False, atom_stereochemistry_matching=False,
+                                       bond_stereochemistry_matching=False)[0] is True
+        # test hill formula passes but isomorphic fails
+        mol1 = Molecule.from_smiles('Fc1ccc(F)cc1')
+        mol2 = Molecule.from_smiles('Fc1ccccc1F')
+        assert mol1.is_isomorphic_with(mol2) is False
+        assert mol2.is_isomorphic_with(mol1) is False
+
+    def test_remap(self):
+        """Test the remap function which should return a new molecule in the requested ordering"""
+        # the order here is CCO
+        ethanol = Molecule.from_file(get_data_file_path('molecules/ethanol.mol2'))
+        # get ethanol in reverse order OCC
+        ethanol_reverse = Molecule.from_smiles('OCC')
+        # get the mapping between the molecules
+        mapping = Molecule.are_isomorphic(ethanol, ethanol_reverse, True)[1]
+        ethanol.add_bond_charge_virtual_site([0, 1], 0.3 * unit.angstrom)
+        # make sure that molecules with virtual sites raises an error
+        with pytest.raises(Exception):
+            remapped = ethanol.remap(mapping, current_to_new=True)
+
+        # remake with no virtual site and remap to match the reversed ordering
+        ethanol = Molecule.from_file(get_data_file_path('molecules/ethanol.mol2'))
+
+        new_ethanol = ethanol.remap(mapping, current_to_new=True)
+
+        def assert_molecules_match_after_remap(mol1, mol2):
+            """Check all of the attributes in a molecule match after being remapped"""
+            for atoms in zip(mol1.atoms, mol2.atoms):
+                assert atoms[0].to_dict() == atoms[1].to_dict()
+            # bonds will not be in the same order in the molecule and the atom1 and atom2 indecies could be out of order
+            remaped_bonds = [bond.to_dict() for bond in mol2.bonds]
+            for bond in mol1.bonds:
+                assert bond.to_dict() in remaped_bonds
+            assert mol1.n_bonds == mol2.n_bonds
+            assert mol1.n_angles == mol2.n_angles
+            assert mol1.n_propers == mol2.n_propers
+            assert mol1.n_impropers == mol2.n_impropers
+            assert mol1.total_charge == mol2.total_charge
+
+        # check all of the properties match as well, torsions and impropers will be in a different order
+        # due to the bonds being out of order
+        assert_molecules_match_after_remap(new_ethanol, ethanol_reverse)
+
+        # catch mappings that are the wrong size
+        too_small_mapping = {0: 1}
+        with pytest.raises(Exception):
+            new_ethanol = ethanol.remap(too_small_mapping, current_to_new=True)
+
+        wrong_index_mapping = {(i + 10, new_id) for i, new_id in enumerate(mapping.values())}
+        with pytest.raises(Exception):
+            new_ethanol = ethanol.remap(wrong_index_mapping, current_to_new=True)
+
+        # test round trip (double remapping a molecule)
+        new_ethanol = ethanol.remap(mapping, current_to_new=True)
+        round_trip_mapping = Molecule.are_isomorphic(new_ethanol, ethanol, return_atom_map=True)[1]
+        round_trip_ethanol = new_ethanol.remap(round_trip_mapping, current_to_new=True)
+        assert_molecules_match_after_remap(round_trip_ethanol, ethanol)
+
     @pytest.mark.parametrize('molecule', mini_drug_bank())
     def test_n_particles(self, molecule):
         """Test n_particles property"""

openforcefield/topology/molecule.py

@@ -1789,7 +1789,7 @@ def __eq__(self, other):
            No effort is made to ensure that the atoms are in the same order or that any annotated properties are preserved.
 
         """
-        return FrozenMolecule.are_isomorphic(self, other, return_atom_map=False)
+        return FrozenMolecule.are_isomorphic(self, other, return_atom_map=False)[0]
 
     def to_smiles(self, toolkit_registry=GLOBAL_TOOLKIT_REGISTRY):
         """
@@ -1900,7 +1900,6 @@ def are_isomorphic(
             bond_stereochemistry_matching=True,
     ):
         """
-        #TODO write tests for imorphic checks and mapping.
         Determines whether the two molecules are isomorphic by comparing their graphs.
 
         Parameters
@@ -1932,8 +1931,6 @@ def are_isomorphic(
         # Do a quick hill formula check first
         if FrozenMolecule.to_hill_formula(mol1) != FrozenMolecule.to_hill_formula(mol2):
             return False, None
-        else:
-            print(FrozenMolecule.to_hill_formula(mol1), FrozenMolecule.to_hill_formula(mol2))
 
         # Build the user defined matching functions
         def node_match_func(x, y):
@@ -1976,10 +1973,15 @@ def find_data_type(data):
             """Find the data type and return the networkx graph"""
 
             try:
+                # Molecule class instance
                 return data.to_networkx()
             except AttributeError:
-                if isinstance(data, nx.Graph):
-                    return data
+                try:
+                    # TopologyMolecule class instance
+                    return data.reference_molecule.to_networkx()
+                except AttributeError:
+                    if isinstance(data, nx.Graph):
+                        return data
 
         mol1_netx = find_data_type(mol1)
         mol2_netx = find_data_type(mol2)
@@ -2033,7 +2035,7 @@ def is_isomorphic_with(self, other):
 
         # what level of matching do we want here?
         # should we expose some options as well?
-        return FrozenMolecule.are_isomorphic(self, other, return_atom_map=False)
+        return FrozenMolecule.are_isomorphic(self, other, return_atom_map=False)[0]
 
     def generate_conformers(self,
                             toolkit_registry=GLOBAL_TOOLKIT_REGISTRY,
@@ -3374,44 +3376,37 @@ def from_mapped_smiles(cls, mapped_smiles):
         :return: openforcefield.topology.molecule.Molecule
         """
 
-        # extract the atomic mapping
-        mapping = []
-        for atom in mapped_smiles.split(':'):
-            try:
-                mapping.append(int(atom[:2]))
-            except ValueError:
-                try:
-                    mapping.append(int(atom[:1]))
-                except ValueError:
-                    pass
+        # create the molecule from the smiles and check we have the right number of indexes
+        # in the mapped SMILES
+        offmol = cls.from_smiles(mapped_smiles, hydrogens_are_explicit=True)
 
         # check we found some mapping
-        if len(mapping) == 0:
+        try:
+            mapping = offmol._properties['atom_map']
+        except KeyError:
             raise SmilesParsingError('The given SMILES has no indexing, please generate a valid explicit hydrogen '
                                      'mapped SMILES using cmiles.')
-        # create the molecule from the smiles and check we have the right number of indexes
-        # in the mapped SMILES
-        offmol = cls.from_smiles(mapped_smiles, hydrogens_are_explicit=True)
 
         if len(mapping) != offmol.n_atoms:
             raise SmilesParsingError('The mapped smiles does not contain enough indexes to remap the molecule.')
 
-        # Make a valid atom map dict Dict[new_index: old_index] and remap
-        # use -1 here as the cmiles maps from 1 not 0
-        atom_mapping = dict((i, index - 1) for i, index in enumerate(mapping))
+        # remap the molecule using the atom map found in the smiles
+        # the order is mapping = Dict[current_index: new_index]
 
-        return offmol.remap(atom_mapping, new_to_old=False)
+        return offmol.remap(mapping, new_to_current=False)
 
     @classmethod
-    def from_qcarchive(cls, qca_json, client_instance=None, allow_undefined_stereo=False):
+    def from_qcschema(cls, qca_dict, client=None, allow_undefined_stereo=False):
         """
         Create a Molecule from  a QCArchive entry based on the cmiles information.
 
-        If we also have a client instance we can go and attach the starting geometry.
+        If we also have a client instance/address we can go and attach the starting geometry.
 
         Parameters
         ----------
-        qca_mol : a QCArchive json representation
+        qca_dict : a QCArchive dict with json encoding
+        #TODO should this also accept the entry record instance that we can type check?
+        client : a qcportal.FractalClient instance or addess of an archive that we can search for the geometry.
 
         Returns
         -------
@@ -3420,10 +3415,13 @@ def from_qcarchive(cls, qca_json, client_instance=None, allow_undefined_stereo=F
 
         Examples
         --------
+        >>> import qcportal as ptl
+        >>> client = ptl.FractalClient()
+        >>>
 
         """
 
-        if 'canonical_isomeric_explicit_hydrogen_mapped_smiles' in qca_json['attributes'].keys():
+        if 'canonical_isomeric_explicit_hydrogen_mapped_smiles' in qca_dict['attributes'].keys():
             # make a new molecule that has been reordered to match the cmiles mapping
             offmol = cls.from_mapped_smiles(qca_json['attributes']['canonical_isomeric_explicit_hydrogen_mapped_smiles'])
             if client_instance is not None:
@@ -3492,44 +3490,54 @@ def from_pdb(cls, file_path, smiles, allow_undefined_stereo=False):
         else:
             raise InvalidConformerError('The PDB and SMILES structures do not match.')
 
-    def remap(self, mapping_dict, new_to_old=True):
+    def remap(self, mapping_dict, current_to_new=True):
         """
         Remap all of the indexes in the molecule to match the given mapping dict
-        :param mapping_dict: A dictionary of the mapping between in the indexes
-        :param new_to_old: The dict is {new_index: old_index} if True else {old_index: new_index}
-        :return:
+        :param mapping_dict: A dictionary of the mapping between in the indexes, this should start from 0.
+        :param current_to_new: The dict is {current_index: new_index} if True else {new_index: current_index}
+        :return: a new openforcefield.topology.molecule.Molecule instance with all attributes transferred
         """
 
         if self.n_virtual_sites != 0:
             raise NotImplementedError('We can not remap virtual sites yet!')
 
+        # make sure the size of the mapping matches the current molecule
+        if len(mapping_dict) != self.n_atoms:
+            raise ValueError(f'There are too many mapping indices({len(mapping_dict)}) for the amount of atoms in this '
+                             f'molecule({self.n_atoms})')
+
         # make two mapping dicts we need new to old for atoms
         # and old to new for bonds
-        if new_to_old:
-            n_to_o = mapping_dict
-            o_to_n = dict((v, u) for u, v in mapping_dict.items())
+        if current_to_new:
+            cur_to_new = mapping_dict
+            new_to_cur = dict((v, u) for u, v in mapping_dict.items())
         else:
-            o_to_n = mapping_dict
-            n_to_o = dict((v, u) for u, v in mapping_dict.items())
+            new_to_cur = mapping_dict
+            cur_to_new = dict((v, u) for u, v in mapping_dict.items())
 
         new_molecule = Molecule()
         new_molecule.name = self.name
 
-        # add the atoms list
-        for i in range(self.n_atoms):
-            # get the old atom info
-            old_atom = self._atoms[n_to_o[i]]
-            new_molecule.add_atom(atomic_number=old_atom.atomic_number,
-                                  formal_charge=old_atom.formal_charge,
-                                  is_aromatic=old_atom.is_aromatic,
-                                  stereochemistry=old_atom.stereochemistry,
-                                  name=old_atom.name)
-        # add the bonds
+        try:
+            # add the atoms list
+            for i in range(self.n_atoms):
+                # get the old atom info
+                old_atom = self._atoms[new_to_cur[i]]
+                new_molecule.add_atom(atomic_number=old_atom.atomic_number,
+                                      formal_charge=old_atom.formal_charge,
+                                      is_aromatic=old_atom.is_aromatic,
+                                      stereochemistry=old_atom.stereochemistry,
+                                      name=old_atom.name)
+        # this is the first time we access the mapping catch an index error here corresponding to mapping that starts
+        # from 0 or higher
+        except IndexError:
+            raise IndexError(f'The mapping supplied is missing a relation corresponding to atom({i})')
+
+            # add the bonds but with atom indexes in a sorted ascending order
         for bond in self._bonds:
-            a1 = o_to_n[bond.atom1_index]
-            a2 = o_to_n[bond.atom2_index]
-            new_molecule.add_bond(atom1=a1,
-                                  atom2=a2,
+            atoms = sorted([cur_to_new[bond.atom1_index], cur_to_new[bond.atom2_index]])
+            new_molecule.add_bond(atom1=atoms[0],
+                                  atom2=atoms[1],
                                   bond_order=bond.bond_order,
                                   is_aromatic=bond.is_aromatic,
                                   stereochemistry=bond.stereochemistry,
@@ -3538,15 +3546,15 @@ def remap(self, mapping_dict, new_to_old=True):
         # remap the charges
         new_charges = np.zeros(self.n_atoms)
         for i in range(self.n_atoms):
-            new_charges[i] = self.partial_charges[n_to_o[i]].value_in_unit(unit.elementary_charge)
+            new_charges[i] = self.partial_charges[new_to_cur[i]].value_in_unit(unit.elementary_charge)
         new_molecule.partial_charges = new_charges * unit.elementary_charge
 
         # remap the conformers there can be more than one
         if self.conformers is not None:
             for conformer in self.conformers:
                 new_conformer = np.zeros((self.n_atoms, 3))
                 for i in range(self.n_atoms):
-                    new_conformer[i] = conformer[n_to_o[i]].value_in_unit(unit.angstrom)
+                    new_conformer[i] = conformer[new_to_cur[i]].value_in_unit(unit.angstrom)
                 new_molecule.add_conformer(new_conformer * unit.angstrom)
 
         # move any properties across