Compute RMSD between two molecules

Let us compute the root-mean-square deviation (RMSD) of atomic positions between the atoms of two superimposed molecules.

from math import sqrt

def compute_rmsd(mol1_indexer, mol2_indexer):
    """Computes RMSD between two superimposed molecules and returns RMSD in angstrom"""

    # check if indexers are of the same size
    if len(mol1_indexer) != len(mol2_indexer):
        print('Molecules are not the same')
        return -1

    sqrd_sum = SBQuantity.squareLength(0.0)
    # go through all atoms in the indexer
    for i in range(len(mol1_indexer)):
        if mol1_indexer[i].elementType != mol2_indexer[i].elementType:
            print('Atoms are not the same')
            return -1
        # get positions of atoms from indexers
        pos1 = mol1_indexer[i].getPosition()
        pos2 = mol2_indexer[i].getPosition()
        # compute distance between atoms
        # and accumulate squared norm for distances
        sqrd_sum += (pos1 - pos2).norm2()

    # take a square root of the sum of squared distances divided by the number of atoms
    rmsd = SBQuantity.sqrt(sqrd_sum / len(mol1_indexer))

    # return RMSD in angstrom
    return rmsd.angstrom.value

def compute_rmsd_2mol(mol1, mol2):
    """Computes several RMSD between two superimposed molecules:
    1) between all atoms,
    2) between backbone atoms,
    3) between alpha-Carbon atoms.
    """
    # get indexers of all atoms
    mol1_all_indexer      = mol1.getNodes('node.type atom')
    mol2_all_indexer      = mol2.getNodes('n.t a')
    # get indexers of all atoms in the backbone
    mol1_backbone_indexer = mol1.getNodes('node.type atom in node.type backbone')
    mol2_backbone_indexer = mol2.getNodes('n.t a in n.t bb')
    # get indexers of all alpha Carbons in the backbone
    mol1_CA_indexer       = mol1.getNodes('"CA" in node.type backbone')
    mol2_CA_indexer       = mol2.getNodes('"CA" in n.t bb')

    rmsd_all      = compute_rmsd(mol1_all_indexer, mol2_all_indexer);
    rmsd_backbone = compute_rmsd(mol1_backbone_indexer, mol2_backbone_indexer);
    rmsd_CA       = compute_rmsd(mol1_CA_indexer, mol2_CA_indexer);

    if (rmsd_all != -1):      print("RMSD (all atoms): %e Angstrom" % rmsd_all)
    if (rmsd_backbone != -1): print("RMSD (backbone) : %e Angstrom" % rmsd_backbone)
    if (rmsd_CA != -1):       print("RMSD (C_alpha)  : %e Angstrom" % rmsd_CA)

    # return a list of RMSDs
    return [rmsd_all, rmsd_backbone, rmsd_CA]


# get all molecules in the active document
indexer = SAMSON.getNodes('node.type molecule')

if len(indexer) < 2:
    print("Less than two molecules are present in the data graph")
else:
    # we assume here that we need to compute RMSD between the first two molecules in the indexer
    mol1 = indexer[0]    # get the first molecule
    mol2 = indexer[1]    # get the second molecule
    # compute RMSD between the two superimposed molecules
    rmsd = compute_rmsd_2mol(mol1, mol2)