Computation of the RMSD

Computing RMSD

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

from math import sqrt
def compute_rmsd(mol1_indexer, mol2_indexer):
        Computes RMSD between two superimposed molecules
        if mol1_indexer.size != mol2_indexer.size:                # check if indexers are of the same size
                print('molecules are not the same')
                return -1
        rmsd = 0.0
        for i in range(mol1_indexer.size):                        # go through all atoms in the indexer
                if mol1_indexer[i].elementType != mol2_indexer[i].elementType:
                        print('atoms are not the same')
                        return -1
                pos1 = mol1_indexer[i].getposition()              # get position of an atom from an indexer
                pos2 = mol2_indexer[i].getposition()
                distance = pos1 - pos2                            # compute distance between atoms
                rmsd += distance.norm2().angstrom                 # accumulate squared norm for distances [in angstrom]
        rmsd = sqrt(rmsd / mol1_indexer.size)                     # take a square root of the sum of squared distances divided by the number of atoms
        return rmsd
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
        mol1_all_indexer      = mol1.getNodes('n.t a')            # get all atoms
        mol2_all_indexer      = mol2.getNodes('n.t a')
        mol1_backbone_indexer = mol1.getNodes('n.t a in n.t bb')  # get all atoms in the backbone
        mol2_backbone_indexer = mol2.getNodes('n.t a in n.t bb')
        mol1_CA_indexer       = mol1.getNodes('"CA" in n.t bb')   # get all alpha Carbons in the 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 [rmsd_all, rmsd_backbone, rmsd_CA]                 # return the list of different RMSD
indexer = SAMSON.getNodes('n.t m')                                # get all molecules
if indexer.size < 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
        rmsd = compute_rmsd_2mol(mol1, mol2)                      # compute RMSD between two superimposed molecules

Comments are closed.