fixed the hill formula to construct in alphabetical order.

Author: joshhorton

openforcefield/tests/test_molecule.py

@@ -538,6 +538,15 @@ def test_hill_formula(self):
         # make sure is not order dependent
         molecule_smiles_reverse = create_reversed_ethanol()
         assert molecule_smiles.hill_formula == molecule_smiles_reverse.hill_formula
+        # make sure single element names are put first
+        order_mol = Molecule.from_smiles('C(Br)CB')
+        assert order_mol.hill_formula == 'C2H6BBr'
+        # test molecule with no carbon
+        no_carb_mol = Molecule.from_smiles('OS(=O)(=O)O')
+        assert no_carb_mol.hill_formula == 'H2O4S'
+        # test no carbon and hydrogen
+        br_i = Molecule.from_smiles('BrI')
+        assert br_i.hill_formula == 'BrI'
         # make sure files and smiles match
         molecule_file = Molecule.from_file(get_data_file_path('molecules/ethanol.sdf'))
         assert molecule_smiles.hill_formula == molecule_file.hill_formula
@@ -549,6 +558,7 @@ def test_hill_formula(self):
         # make sure the networkx matches
         assert molecule_smiles.hill_formula == Molecule.to_hill_formula(molecule_smiles.to_networkx())
 
+
     def test_isomorphic_general(self):
         """Test the matching using different input types"""
         # check that hill formula fails are caught

openforcefield/tests/test_topology.py

@@ -429,7 +429,7 @@ def test_from_openmm_missing_conect(self):
 
         molecules = []
         molecules.append(Molecule.from_smiles('CCO'))
-        with pytest.raises(ValueError, match='No match found for molecule C1. This would be a '
+        with pytest.raises(ValueError, match='No match found for molecule C. This would be a '
                                              'very unusual molecule to try and parameterize, '
                                              'and it is likely that the data source it was '
                                              'read from does not contain connectivity '

openforcefield/topology/molecule.py

@@ -2864,6 +2864,7 @@ def to_hill_formula(molecule):
         """
 
         from openforcefield.topology import TopologyMolecule
+        import collections
         # check for networkx then assuming we have a Molecule or TopologyMolecule instance just try and
         # extract the info. Note we do not type check the TopologyMolecule due to cyclic dependencies
         if isinstance(molecule, nx.Graph):
@@ -2881,32 +2882,29 @@ def to_hill_formula(molecule):
                                           f'openforcefield.topology.topology.TopologyMolecule or networkx representaion '
                                           f'of the molecule.')
 
-        # Now sort the elements, method taken from topology._networkx_to_hill_formula
-        # Count the number of instances of each atomic number
-        at_num_to_counts = dict([(unq, atom_nums.count(unq)) for unq in atom_nums])
-
-        symbol_to_counts = {}
-        # Check for C and H first, to make a correct hill formula (remember dicts in python 3.6+ are ordered)
-        if 6 in at_num_to_counts:
-            symbol_to_counts['C'] = at_num_to_counts[6]
-            del at_num_to_counts[6]
-
-        if 1 in at_num_to_counts:
-            symbol_to_counts['H'] = at_num_to_counts[1]
-            del at_num_to_counts[1]
-
-        # Now count instances of all elements other than C and H, in order of ascending atomic number
-        sorted_atom_nums = sorted(at_num_to_counts.keys())
-        for atom_num in sorted_atom_nums:
-            symbol_to_counts[Element.getByAtomicNumber(atom_num).symbol] = at_num_to_counts[atom_num]
-
-        # Finally format the formula as string
-        formula = ''
-        for ele, count in symbol_to_counts.items():
-            if count == 1:
-                count = ''
-            formula += f'{ele}{count}'
-        return formula
+        # make a correct hill formula representation following this guide
+        # https://en.wikipedia.org/wiki/Chemical_formula#Hill_system
+
+        # create the counter dictionary using chemical symbols
+        atom_symbol_counts = collections.Counter(Element.getByAtomicNumber(atom_num).symbol for atom_num in atom_nums)
+
+        formula = []
+        # Check for C and H first, to make a correct hill formula
+        for el in ['C', 'H']:
+            if el in atom_symbol_counts:
+                count = atom_symbol_counts.pop(el)
+                formula.append(el)
+                if count > 1:
+                    formula.append(str(count))
+
+        # now get the rest of the elements in alphabetical ordering
+        for el in sorted(atom_symbol_counts.keys()):
+            count = atom_symbol_counts.pop(el)
+            formula.append(el)
+            if count > 1:
+                formula.append(str(count))
+
+        return "".join(formula)
 
     def chemical_environment_matches(self,
                                      query,

openforcefield/topology/topology.py

@@ -1653,7 +1653,7 @@ def from_openmm(cls, openmm_topology, unique_molecules=None):
             if match_found is False:
                 hill_formula = Molecule.to_hill_formula(omm_mol_G)
                 msg = f'No match found for molecule {hill_formula}. '
-                probably_missing_conect = ['C1', 'H1', 'O1', 'N1', 'P1', 'S1', 'F1', 'Cl1', 'Br1']
+                probably_missing_conect = ['C', 'H', 'O', 'N', 'P', 'S', 'F', 'Cl', 'Br']
                 if hill_formula in probably_missing_conect:
                     msg += ('This would be a very unusual molecule to try and parameterize, '
                             'and it is likely that the data source it was read from does not '