There are two ways I have approached this type of problem in the past. Which method is preferable depends on the types of structures you are trying to filter through.
Using the Kabsch Algorithm:
Roughly the way this method works is as follows:
- Rotate all molecules into a common frame determined by the current structure
- Take the difference of cartesian coordinates between all pairs of molecules
- Throw out any structure which has differences smaller than some threshold
The hard part of this approach is how to rotate all the molecules into a common frame. The simplest way to do this that I am aware of is using the Kabsch algorithm. The kabsch algorithm calculates the optimal rotation matrix to a point which minimizes the RMSD between two points. The rotation matrix is typically calculated via the singular-value decomposition. There is some interesting math behind why the SVD is the appropriate way to do this, but I think it's probably unnecessary here.
Here is a python implementation on github which is designed to rotate molecules into the same frame using the kabsch algorithm. I have used this code before and it works as expected.
So, the way you do this in total is to read all of your structures into a list of numpy arrays, loop through this list and rotate every molecule after the current on onto the axes of the current molecule. Then, take the difference between this molecule and all subsequent molecules. If the resulting matrix is appropriately close to all zeros (probably by maximum vector length, but whatever metric is probably fine), then you can remove all those molecules which meet your convergence criteria.
Keep doing this process until you reach the end of the list. As a side note, you should pre-process the molecules by shifting everything by its centroid.
If you had to do this with like millions of structures which are very large, then the algorithm would be pretty expensive as you have to do a bunch of matrix factorizations and the algorithms is worst-case $O(N^2)$ for the comparison step.
Graph-Based Approach
Another way to filter out duplicate molecules is by representing each molecule as a graph. Forming the graph is fairly simple as long as you have a good measure of when two atoms are connected. So, for instance, the edges of the graph are likely to be represented by covalent bonds. One can also represent the edges by hydrogen bonds if you're working with a van der Waal's cluster of some kind.
The way this method works is as follows:
- Build a graph representing each molecule based on some connectivity criteria
- Perform an isomorphism check among all pairs of graphs, keeping only one of each unique graph
In theory, the hardest part of this is doing the isomorphism check, but there are great software packages which can do this for you, such as networkx for a Python option.
If you use networkx, probably the easiest way to build the graphs is to determine the connectivity of all atoms in each molecule, and build an adjacency matrix. networkx can then make a graph object from this adjacency matrix. Then you do just as in the previous method and loop through the pairs of molecules removing any which turn out to be isomorphic to the current reference graph.
Now, this method seems pretty easy, but a graph representation of a molecule is not unique. For instance, all graphs representing the boat, chair, and planar conformations of cyclohexane result in identical graphs. In order to make the conformations distinguishable, you have to attach weights to either the edges or the nodes of the graph. Probably the easiest thing to do would be to attach a list of angles of each triplet of atoms. You would also need to label the handedness of each chiral center, as I don't think angles would be enough there.
The advantage of the graph approach is that it avoids problems with numerical precision, and it can be fairly fast I think. The drawback is having to make these weights for the nodes if you have to keep around minima which only differ by rotation of atoms in space.
Hopefully this is helpful!
