This is a simple Python3 script that is used to find the chemically equivalent atoms in a molecule. The molecule is defined with any geometry file that can be read by MDAnalysis.
A pair of atoms is defined as equivalent if:
- they are the same element
- removing one or the other (along with all its bond) from the molecular graph generate omomorphic graphs
❌ bond orders menaing that all bonds are just equal (this is needed also because very few formats do support specification for bond order)
❌ enantiotopic atoms and enantiomerism in general since everything is encoded in a graph without any knowledge about threedimensional structure.
-
The molecule file is loaded. At this stage element information are guessed (if needed) as well as bond information (for this we use the algorithm described in "A rule-based algorithm for automatic bond type perception"); a check on overconnected atoms is also performed based on chemical valence rules.
-
The loaded and modified molecule is converted in a networkx grah. The final graph has one node for every atom labeled with the element and one edge for every bond without any label.
-
For each atom, the equivalence with any other atoms not already in an equivalence chain is checked. The following sequence of checks is done, if one fails, the subsequent are skipped (in order to improve efficiency)
a. The element of the two atom is checked,
b. The number and elements of all the first neighbour atoms of the two is checked,
c. If the two atoms have just one first neighbour (they are probably hydrogens) check for the second neighbour of the two,
d. Create two graphs removing the two atoms and their bonds from the original molecular graph and check the isomorphism of those two.
-
A list of equivalent atoms is printed or saved.
The easies way of using this script is with uv.
Clone the repository, enter it and then:
uv run main.py
For more details about its usage take a look at the help
uv run main.py --help
It could be useful especially in molecular mechanics contextes. The script was originally written to facilitaty parameters optimization performed with poltype2 which is particularly slow in doing this step. But you can use it for whatever you like!