Computation of the RMSD

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

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