Print chemical formula from list of atoms in Python?

Question

Lets consider the following code which uses the Mole object:

import pyscf
mol = pyscf.M(
    atom = '''
        O
        H  1  0.9576257
        H  1  0.9576257  2 104.51
    '''
)

Is there an easy way to print the chemical formula: H2O?

In the above-linked documentation page, I couldn't find any example that resembles this type of feature, and when I run a PySCF kernel() the printout doesn't ever say H2O so perhaps PySCF doesn't "know" that the molecule can be called H2O.

I searched on Google for "chemical formula from ZMAT" and some variants of that search phrase, without any success. I understand that things can get quite complicated when there's even just a few atoms, but for atoms, diatomics, and millions of polyatomics, it doesn't seem unreasonable to have a utility that can print or at least guess a chemical formula, especially if the utility were to use a database of known molecules to assist it.

Is there a simple way to print the formula from the ZMAT or XYZ file that would work with PySCF?

Answer 1

I like Geoff's proposal to use the Hill system. Here's an implementation.

import collections

def to_hill_notation(elements):
    c = collections.Counter(elements)
    if "C" in c:
        #C and H come first
        priority = lambda x: -1 if x in "CH" else 0
    else:
        #no special priority
        priority = lambda x: 0
    sorted_keys = sorted(c.keys(), key=lambda element: (priority(element), element))
    result = []
    for element in sorted_keys:
        amount = c[element]
        if amount == 1:
            result.append(element)
        else:
            result.append(element + str(amount))
    return "".join(result)

print(to_hill_notation(['O', 'H', 'H'])) #H20
print(to_hill_notation(['C', 'H', 'C', 'C', 'H', 'C', 'H', 'H', 'H', 'C', 'H', 'C'])) #C6H6
print(to_hill_notation(['O', 'H', 'O', 'H', 'H', 'C', 'O', 'C', 'C', 'H', 'H', 'O', 'C', 'H', 'H', 'O', 'H', 'O', 'O', 'C', 'H', 'H', 'H', 'C', 'H', 'C', 'H', 'C', 'O', 'C', 'H', 'H', 'C', 'H', 'C', 'H', 'H', 'H', 'O', 'O', 'H', 'O', 'H', 'C', 'H'])) #C12H22O11
print(to_hill_notation(['Cl', 'H', 'C', 'H', 'C', 'H'])) #C2H3Cl
print(to_hill_notation(['O', 'O', 'O', 'O', 'Fe', 'O', 'O', 'O', 'O', 'Fe', 'S', 'O', 'S', 'O', 'S', 'O', 'O'])) #Fe2O12S3

You can't pass a Mole object to the function directly, but you should be able to get a list of the Mole's elements using mol.elements. So the call would look something like:

print(to_hill_notation(mol.elements))

Answer 2

The tricky part comes if you have multiple fragments and want disconnected formulas. But in general, the Hill system for chemical formulas is used, and available in Open Babel, RDKit, etc.

As indicated in the comments, you'd just iterate through the atoms and count the number of each. For example:

# this is pseudocode
elements = {} # amend the end point if you don't want to count transuranic
for atom in mol.elements:
   elements{atom} += 1

# print in Hill order

You'll want to make sure the C and H elements come first if they exist, otherwise it's in sorted alphabetical order.

Answer 3

I think Kevin's answer is stellar (both in terms of understanding atoms and formulas), and the recommendation and hints from Geoff's answer helped a lot in putting things into perspective, so here is the answer from a naive perspective:

atoms_string = '''
        O
        H  1  0.9576257
        H  1  0.9576257  2 104.51
    '''

elements = {}
for line in atoms_string.strip().split("\n"):
    element = line.split()[0]
    if element in elements:
        elements[element] += 1
    else:
        elements[element] = 1

formula = ''.join(f'{element}{count}' for element, count in elements.items())
print(formula)

Which print O1H2. I think the OP mentioned they would prefer it output H2O in this case, but I don't know much about atoms and whatnot, compared to most people here...so this is just a naive implementation. It is however compact(?), but if you want better formatting for the output, Kevin's answer is a better fit.

Note that the above depends heavily on the formatting of its input, and it may not work on different inputs.

EDIT: Still don't know what I'm doing, but the following output H2O:

atoms_string = '''
        O
        H  1  0.9576257
        H  1  0.9576257  2 104.51
    '''

elements = {}
for line in atoms_string.strip().split("\n"):
    element = line.split()[0]
    if element in elements:
        elements[element] += 1
    else:
        elements[element] = 1

hill_order = sorted(elements.items(), key=lambda x: x[0])
hill_order = sorted(hill_order, key=lambda x: x[1], reverse=True)

formula = ''.join(f'{element}{count if count > 1 else ""}' for element, count in hill_order)
print(formula)

I tried to implement Hill ordering and remove/ignore count when they were equal to 1. I don't think this will work on other inputs, but for the example above, this works.