Cacher is the code snippet organizer for pro developers

# Better to run this script on Jupyter and print the result in every steps
# It would minimize the error on the final data

import matplotlib
from matplotlib import pyplot as plt
import pytraj as pt
import numpy as np

# Define the moving average function
def movingaverage(interval, window_size):
  window = np.ones(int(window_size))/float(window_size)
  return np.convolve(interval, window, 'same')


# Loading trajectory and Reference
traj = pt.iterload('1gh0-cyc-200ns.nc', '1gh0_cyc_nobox.prmtop') # the trajectory file loaded with 'iterload' because the RAM doesn't fit for 50 GB of data; if you have sufficient RAM size, use 'load' instead
ref = pt.load('ekuil5-nobox.rst7','1gh0_cyc_nobox.prmtop') # reference taken from the last step of equilibration

# Compute the RMSD for all protein and all ligand
data_rmsd = pt.pmap(pt.rmsd, traj, ref=ref, mask=":1-1011@CA,C,N", n_cores=4) # I used 4 cores because I run it on Core i5 that has 4 core 4 threads
data_rmsd_cyc_all = pt.pmap(pt.rmsd, traj, ref=ref, mask=":CYC", n_cores=4) # Just for the ligand

# Let's plot it
ax = plt.subplot(111)
x1 = 0.002*np.arange(data_rmsd['RMSD_00001'].size)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
major_ticks = np.arange(0, 201, 40)
ax.set_xticks(major_ticks)
av_all= movingaverage(data_rmsd['RMSD_00001'], 50)
av_all_cyc = movingaverage(data_rmsd_cyc_all['RMSD_00001'], 50)
ax.plot(x1, av_all, linewidth='0.6')
ax.plot(x1, av_all_cyc, linewidth='0.6')
plt.axis([0, 200, 0, 12])
plt.xlabel('Timestep (ns)', fontsize=18)
plt.ylabel('RMSD (Angstrom)', fontsize=18)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.savefig('RMSD-all.png', dpi = 1000, bbox_inches = 'tight')

# Save it to csv file to plot it on another plotting software
bx = np.array(av_all)
bx_cyc = np.array(av_all_cyc)
zipped_bx = zip(x1, bx)
zipped_bx_cyc = zip (x1, bx_cyc)
np.savetxt("RMSD_all.csv", zipped_bx, delimiter=',')
np.savetxt("RMSD_cyc_all.csv", zipped_bx_cyc, delimiter=',')

# Compute RMSD for protein and ligand on chain A only
data_rmsd_chaina = pt.pmap(pt.rmsd, traj, ref=ref, mask=":1-163@CA,C,N", n_cores=4)
data_rmsd_chaina_cyc = pt.pmap(pt.rmsd, traj, ref=ref, mask=":163", n_cores=4)

# And again, plot it
ax = plt.subplot(111)
xa = 0.002*np.arange(data_rmsd_chaina['RMSD_00001'].size)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
major_ticks = np.arange(0, 201, 40)
ax.set_xticks(major_ticks)
av_chaina = movingaverage(data_rmsd_chaina['RMSD_00001'], 50)
av_chaina_cyc = movingaverage(data_rmsd_chaina_cyc['RMSD_00001'], 50)
ax.plot(xa, av_chaina, linewidth='0.6')
ax.plot(xa, av_chaina_cyc, linewidth='0.6')
plt.axis([0, 200, 0, 6])
plt.xlabel('Timestep (ns)', fontsize=18)
plt.ylabel('RMSD (Angstrom)', fontsize=18)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.savefig('RMSD-chainA.png', dpi = 1000, bbox_inches = 'tight')

# And just repeat the steps for other chain