In [3]:

    

import MDAnalysis
import matplotlib.pyplot as plt
import numpy as np
from MDAnalysis.analysis.align import *
from MDAnalysis.analysis.rms import rmsd
%matplotlib inline


    

In [13]:

    

def ligRMSD(u,ref):
    """
    This function produces RMSD data and plots for ligand. 
    :input 
        1) Universe of Trajectory
        2) reference universe
    :return
        1) matplot object
        2) array for RMSD data.
        
    """
    RMSD_lig = []
    ligand = u.select_atoms("(resid 142:146) and not name H*") ## include selection based on user description
    #current = u.select_atoms("segname BGLC and not name H*")
    reference = ref.select_atoms("(resid 142:146) and not name H*")
    for ts in u.trajectory:
        A = ligand.coordinates()
        B = reference.coordinates()
        C = rmsd(A,B)
        RMSD_lig.append((u.trajectory.frame, C))
    RMSD_lig = np.array(RMSD_lig)
    #print RMSD_lig
    import matplotlib.pyplot as plt
    ax = plt.subplot(111)
    ax.plot(RMSD_lig[:,0], RMSD_lig[:,1], 'r--', lw=2, label=r"$R_G$")
    ax.set_xlabel("Frame")
    ax.set_ylabel(r"RMSD of ligand ($\AA$)")
    #ax.figure.savefig("RMSD_ligand.pdf")
    #plt.draw()
    handles, labels = ax.get_legend_handles_labels()
    ax.legend(handles, labels, loc = 'lower left')
    return ax,  RMSD_lig


    

In [6]:

    

trj = '50_frame.dcd'
top = './41wl_ff.psf'


    

In [12]:

    

u = MDAnalysis.Universe('41wl_ff.psf','50_frame.dcd')
ref = MDAnalysis.Universe('41wl_ff.psf','50_frame.dcd') 
ligandRMSD = []
fig,ligandRMSD = ligRMSD(u,ref)
#print caRMSD
np.savetxt("SimAnaRep-allRMSD.data", ligandRMSD)
fig.figure.savefig("SimAnaRep-RMSD.pdf")


    

    




[[  0.00000000e+00   1.08735602e-06]
 [  1.00000000e+00   3.68662714e+00]
 [  2.00000000e+00   2.12019302e+00]
 [  3.00000000e+00   1.60456652e+00]
 [  4.00000000e+00   1.57198797e+00]
 [  5.00000000e+00   2.84287226e+00]
 [  6.00000000e+00   4.90751869e+00]
 [  7.00000000e+00   3.36932707e+00]
 [  8.00000000e+00   3.23793130e+00]
 [  9.00000000e+00   2.81030673e+00]
 [  1.00000000e+01   2.96390197e+00]
 [  1.10000000e+01   1.80263100e+00]
 [  1.20000000e+01   1.57591836e+00]
 [  1.30000000e+01   2.64336319e+00]
 [  1.40000000e+01   2.64373063e+00]
 [  1.50000000e+01   3.43870722e+00]
 [  1.60000000e+01   3.09839813e+00]
 [  1.70000000e+01   3.24417248e+00]
 [  1.80000000e+01   2.66968311e+00]
 [  1.90000000e+01   3.73360857e+00]
 [  2.00000000e+01   1.95119562e+00]
 [  2.10000000e+01   3.17663864e+00]
 [  2.20000000e+01   3.87500922e+00]
 [  2.30000000e+01   2.48765420e+00]
 [  2.40000000e+01   2.36807138e+00]
 [  2.50000000e+01   2.67178418e+00]
 [  2.60000000e+01   3.63333534e+00]
 [  2.70000000e+01   1.96125171e+00]
 [  2.80000000e+01   3.24832209e+00]
 [  2.90000000e+01   3.41356751e+00]
 [  3.00000000e+01   2.59487464e+00]
 [  3.10000000e+01   2.35808403e+00]
 [  3.20000000e+01   4.10424226e+00]
 [  3.30000000e+01   3.00861139e+00]
 [  3.40000000e+01   2.99031474e+00]
 [  3.50000000e+01   3.16022017e+00]
 [  3.60000000e+01   2.94211087e+00]
 [  3.70000000e+01   3.00621353e+00]
 [  3.80000000e+01   1.52727089e+00]
 [  3.90000000e+01   3.05931957e+00]
 [  4.00000000e+01   1.28726690e+00]
 [  4.10000000e+01   3.67078576e+00]
 [  4.20000000e+01   2.19999724e+00]
 [  4.30000000e+01   3.36467798e+00]
 [  4.40000000e+01   2.94593173e+00]
 [  4.50000000e+01   4.65756391e+00]
 [  4.60000000e+01   3.56851330e+00]
 [  4.70000000e+01   1.84376409e+00]
 [  4.80000000e+01   1.51971426e+00]
 [  4.90000000e+01   2.93048855e+00]]