12
$\begingroup$

Does there exist a program that can compare a series of (ideally.xyz) files to each other to identify duplicate structures, and then identify the structures in some way?

I know that obabel has the -oconfabreport function. However it only identifies if a duplicate exists (I have used -xr for RMSD), but not the pair of structures themselves. Using this wouldn’t be bad, but across a large range of conformers it is difficult.

I am looking to compare conformers of the same molecule, so no changes in number of atoms or charge.

$\endgroup$
  • 1
    $\begingroup$ To be clear, you are trying to find a solution which does not just check pairwise using obebel? I suspect that pairwise checking is going to be close to optimal assuming you check for # of atoms etc before checking to save time. Having an idea of structure size and count would help signifigantly $\endgroup$ – Tristan Maxson Aug 25 at 1:22
  • 2
    $\begingroup$ @kskinnerx16 DM, I have something maybe. github.com/craldaz/prune_structures $\endgroup$ – Cody Aldaz Aug 25 at 2:12
  • 1
    $\begingroup$ Crest has this implemented, so maybe you can just feed it a 'movie' and let it do its magic. $\endgroup$ – Martin - マーチン Aug 26 at 14:06
  • 1
    $\begingroup$ Yes indeed, sorry for not including a link. $\endgroup$ – Martin - マーチン Aug 26 at 14:30
  • 3
    $\begingroup$ I would do this using the Kabsch algorithm, which rotates all molecules into a common frame by minimizing the RMSD of pairs of vectors. You need to make sure the atoms in each molecule are sorted properly. I didn't look too closely, but I think @CodyAldaz code is doing this. You then just throw out structures with an RMSD smaller than some tolerance. There are many implementations of the Kabsch algorithm online, and it's pretty easy to write yourself. $\endgroup$ – jheindel Sep 1 at 20:45
7
$\begingroup$

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!

| cite | improve this answer | |
$\endgroup$
5
$\begingroup$

I am looking to compare conformers of the same molecule, so no changes in number of atoms or charge.

TM-align to perform rigid-body structural alignment and to identify the $RMSD$ between two structures.

If I understand your problem correctly, you have $N$ conformers of the same molecule. And you need to identify structural similarity between them.

So, choose one of the conformers as a reference (an arbitrary one) and align all the other $(N-1)$ onto this reference using TM-align. TM-align will output the $RMSD$ values for each such alignment and also the corresponding spatial transform (and the TM-score as well).

Case 1: your conformers have different conformations (e.g., different values of dihedrals of the backbones), then the $RMSD$ produced by TM-align will reflect this difference.

Case 2: In case your conformers have the same conformations (i.e., same backbone dihedrals) and differ only in their spatial positions. Then you may record the spatial transforms produced by TM-align. After that, you may pipeline these spatial transforms to another software --

RigidRMSD performs calculation of the $RMSD$ between two structures without performing spatial transforms! (this is achieved by formulating the problem in terms of the inertia tensor -- see more details in the corresponding publication: Rapid determination of RMSDs corresponding to macromolecular rigid body motions )

Thus, you get the $RMSD$ values that correspond to spatial transforms.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Excellent answers. Do you know of a wiki or some examples online anywhere for RigidRMSD? $\endgroup$ – kskinnerx16 Sep 22 at 13:49
  • 1
    $\begingroup$ @kskinnerx16 I do not think that any wiki exists for RigidRMSD, besides what is given on their web-page and in the paper. But I can confirm that their method and implementation work perfectly. I have performed an independent implementation and extensive testing of their method. All results perfectly coincide. $\endgroup$ – u.heap_f3 Sep 22 at 17:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.