library(bio3d)
library(ncdf)
pdb.open <- read.pdb("/net/gulrotkake/slars/groel_nma_korgpil/1SVT/chainA/start.pdb")
pdb.closed <- read.pdb("1XCK_reference.pdb")
pdb <- pdb.closed
ca.inds <- atom.select(pdb.open, "calpha")
ca.inds.closed <- atom.select(pdb.closed, "calpha")

trj0 <- read.ncdf("../traj/300ns_short_15000frames_noWAT.nc")
trj1 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/123_1SVT_chainA_apo_seed1/results/20ns_2000frames_noWAT.nc")
trj2 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/123_1SVT_chainA_apo_seed2/results/20ns_2000frames_noWAT.nc")
trj3 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/123_1SVT_chainA_apo_seed3/results/20ns_2000frames_noWAT.nc")
trj4 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/123_1SVT_chainA_apo_seed4/results/20ns_2000frames_noWAT.nc")
trj5 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/123_1SVT_chainA_apo_seed5/results/20ns_2000frames_noWAT.nc")



rmsd0 <- rmsd( pdb.closed$xyz, trj0[1:10000,], a.inds=inds$xyz, b.inds=inds$xyz, fit = TRUE )
rmsd1 <- rmsd( pdb.closed$xyz, trj1, a.inds=inds$xyz, b.inds=inds$xyz, fit = TRUE )
rmsd2 <- rmsd( pdb.closed$xyz, trj2, a.inds=inds$xyz, b.inds=inds$xyz, fit = TRUE )
rmsd3 <- rmsd( pdb.closed$xyz, trj3, a.inds=inds$xyz, b.inds=inds$xyz, fit = TRUE )
rmsd4 <- rmsd( pdb.closed$xyz, trj4, a.inds=inds$xyz, b.inds=inds$xyz, fit = TRUE )
rmsd5 <- rmsd( pdb.closed$xyz, trj5, a.inds=inds$xyz, b.inds=inds$xyz, fit = TRUE )



pdf("rmsd_all.pdf", h=10)
par(mfrow=c(6,1), mar=c(2,4,0.5,1))
plot(rmsd0,  col="gray", type="o", pch = 16, cex = 0.4, ylab="RMSD", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(0,2500,5000,7500,10000), label=c(0,50,100,150,200))
box()
mtext("A", side=3, line=-1.7, at=0, cex=1.3)
plot(rmsd1,  col="gray", type="o", pch = 16, cex = 0.4, ylab="RMSD", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(0,5,10,15,20))
box()
mtext("B", side=3, line=-1.7, at=0, cex=1.3)
plot(rmsd2,  col="gray", type="o", pch = 16, cex = 0.4, ylab="RMSD", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(0,5,10,15,20))
box()
mtext("C", side=3, line=-1.7, at=0, cex=1.3)
plot(rmsd3,  col="gray", type="o", pch = 16, cex = 0.4, ylab="RMSD", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(0,5,10,15,20))
box()
mtext("D", side=3, line=-1.7, at=0, cex=1.3)
plot(rmsd4,  col="gray", type="o", pch = 16, cex = 0.4, ylab="RMSD", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(0,5,10,15,20))
box()
mtext("E", side=3, line=-1.7, at=0, cex=1.3)
plot(rmsd5,  col="gray", type="o", pch = 16, cex = 0.4, ylab="RMSD", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(0,5,10,15,20))
box()
mtext("F", side=3, line=-1.7, at=0, cex=1.3)
dev.off()