7
$\begingroup$

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?

$\endgroup$
5
  • 1
    $\begingroup$ Prof. Dattani, can't itertools be used? It should be able to count the number of elements and can be printed right? This is just my thoughts, and I agree that this doesn't need PySCF $\endgroup$ Jan 19, 2023 at 12:58
  • $\begingroup$ @HemanthHaridas there's been three answers so far and none of them used itertools. Can you write an answer that uses itertools? $\endgroup$ Jan 20, 2023 at 1:19
  • $\begingroup$ I meant the use of Counter function from collections module. itertools is used for generating all possible combinations from a list. Apologies for the confusion. $\endgroup$ Jan 20, 2023 at 6:34
  • $\begingroup$ I do it via Regex, assuming you're using python. $\endgroup$
    – Hendrix13
    Feb 14, 2023 at 16:34
  • $\begingroup$ @Hendrix13 I'm using Python, perhaps you could expand that comment into an answer? The answer that I accepted is far less compact than what I had in mind! $\endgroup$ Feb 14, 2023 at 16:42

3 Answers 3

3
$\begingroup$

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))
$\endgroup$
4
  • $\begingroup$ Or just sorted_keys = sorted(c, key=lambda element: (element not in 'CH', element)) $\endgroup$
    – Stef
    Jan 19, 2023 at 20:58
  • $\begingroup$ One can clean this a bit more: ` def parse(symbol_c): symbol, c = symbol_c if c == 1: return symbol return f"{symbol}{c}" return "".join(map(parse, sorted(c.items()))) ` no need to complicate things, a tuple sorting is item based, so only the first item will be unique any ways. $\endgroup$
    – nickpapior
    Jan 24, 2023 at 9:54
  • $\begingroup$ A slightly modified version of the code in this answer, is the one I'm using currently. It's still taking up more space in my script than I would like (meaning that my scripts are now so long that I can't see them all at once on my monitor, and I have to do a lot more scrolling), so any more ideas to compactify this script would be very welcome! $\endgroup$ Feb 6, 2023 at 5:24
  • $\begingroup$ By the way Kevin, since you wrote your own programming language, you might be interested in the new Programming Language Design SE which is anticipated to go live in a few months. $\endgroup$ Feb 6, 2023 at 5:25
7
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ +1 thanks very much! This "pseudocode" is essentially real Python, so I got confused by the use of elements{atom} instead of elements[atom] and the fact that the elements[atom] elements were never initialized. To do it this way, I ended up having to do two loops: for atom in mol.elements: ... elements[atom] = 0 then for atom in mol.elements: ... elements[atom] += 1. When printed this gives me: {'O': 1, 'H': 2}. $\endgroup$ Jan 19, 2023 at 20:36
  • 3
    $\begingroup$ There is no need to count the elements "manually" like that. This is python. Just do from collections import Counter; elements = Counter(mol.elements) $\endgroup$
    – Stef
    Jan 19, 2023 at 20:55
  • $\begingroup$ @Stef Yes! Obligatory xkcd: import antigravity. $\endgroup$
    – Lodinn
    Jan 20, 2023 at 19:43
2
$\begingroup$

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.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .