library(bio3d)
library(ncdf)


pdb1 <- read.pdb("1XCK_reference.pdb")
ca.inds <- atom.select(pdb1, "///9:515///CA/")


prefix116 <- "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/rmsd_monomer/"
prefix129 <- "/net/gulrotkake/slars/groel_md/1XCK/129_1XCK_MGATP/results/rmsd_monomer/"
prefix116b <- "/net/lutefisk/slars/groel_md/1XCK/116_1XCK_apo_temp350/results/rmsd_monomer/"
prefix129b <- "/net/lutefisk/slars/groel_md/1XCK/129_1XCK_MGATP_temp350/results/rmsd_monomer/"

prefix117 <- "/net/gulrotkake/slars/groel_md/1XCK/117_1XCK_hepta/results/rmsd_monomer/"
prefix128 <- "/net/gulrotkake/slars/groel_md/1XCK/128_1XCK_hepta_MGATP/results/rmsd_monomer/"


prefix171 <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/rmsd_monomer/"


sim116 <- NULL
sim129 <- NULL
sim117 <- NULL
sim128 <- NULL
sim171 <- NULL

for ( i in 1:14 ) {
  tmprmsd <- read.table(paste(prefix116, "50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim116$rmsd=rbind(sim116$rmsd, tmprmsd[,"V2"])
}

for ( i in 1:14 ) {
  tmprmsd <- read.table(paste(prefix129, "50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim129$rmsd=rbind(sim129$rmsd, tmprmsd[,"V2"])
}

for ( i in 1:7 ) {
  tmprmsd <- read.table(paste(prefix117, "50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim117$rmsd=rbind(sim117$rmsd, tmprmsd[,"V2"])
}
for ( i in 1:7 ) {
  tmprmsd <- read.table(paste(prefix128, "50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim128$rmsd=rbind(sim128$rmsd, tmprmsd[,"V2"])
}


for ( i in 1:14 ) {
  tmprmsd <- read.table(paste(prefix171, "50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim171$rmsd=rbind(sim171$rmsd, tmprmsd[,"V2"])
}


for ( i in 1:14 ) {
  tmprmsd <- read.table(paste(prefix116b, "50ns_5000frames_RMSD_", i , ".dat", sep=""))
  #sim116$rmsdB=rbind(sim116$rmsdB, tmprmsd[,"V2"])
}

for ( i in 1:14 ) {
  tmprmsd <- read.table(paste(prefix129b, "50ns_5000frames_RMSD_", i , ".dat", sep=""))
  #sim129$rmsdB=rbind(sim129$rmsdB, tmprmsd[,"V2"])
}

sim116$rmsd=cbind(sim116$rmsd, sim116$rmsdB)
sim129$rmsd=cbind(sim129$rmsd, sim129$rmsdB)


sim116$rmsd.mean <- colMeans(sim116$rmsd)

sim129$rmsd.mean$atp <- colMeans(sim129$rmsd[1:7,])
sim129$rmsd.mean$free <- colMeans(sim129$rmsd[8:14,])

sim171$rmsd.mean$atp <- colMeans(sim171$rmsd[1:7,])
sim171$rmsd.mean$free <- colMeans(sim171$rmsd[8:14,])

sim117$rmsd.mean <- colMeans(sim117$rmsd)
sim128$rmsd.mean <- colMeans(sim128$rmsd)






pdf("rmsd_monomers.pdf", w=10, h=5)
plot(sim116$rmsd[1,], type='l', col="grey80", ylim=c(0.5, 5.5),
     ylab="RMSD (Å)", xlab="Time (ns)", axes=FALSE )
     

#for ( i in 2:14 ) {
#  lines(sim116$rmsd[i,], lty=2, col="grey80")
#}

for ( i in 1:7 ) {
  lines(sim171$rmsd[i,], lty=2, col="grey80")
}






lines(sim116$rmsd.mean, lty=2, col="blue")
lines(sim129$rmsd.mean$atp, lty=2, col="red")
#lines(sim129$rmsd.mean$free, lty=2, col="orange")
#lines(sim171$rmsd.mean$free, lty=2, col="green")
#lines(sim171$rmsd.mean$atp, lty=2, col="black")

lines(sim117$rmsd.mean, lty=2, col="orange")
lines(sim128$rmsd.mean, lty=2, col="green")



axis(2)
axis(1, at=c(0,1000,2000,3000,4000,5000), labels=c(0,10,20,30,40,50))
box()
legend(0,5.5, legend=c("Individual monomers", "Means of ATP Bound", "Means of ATP free", "Means of APO"), col=c("grey80", "red", "orange", "blue"), lty=2)

dev.off()













sim116 <- NULL
sim129 <- NULL

for ( i in 1:14 ) {
  tmptraj <- read.ncdf(paste("../116_1XCK_apo/results/traj_monomer/500ns_noWAT_5000frames_chain_", i, ".nc", sep=""))
  sim116$traj = c(sim116$traj, tmptraj)
}

for ( i in 1:14 ) {
  tmptraj <- read.ncdf(paste("../129_1XCK_MGATP/results/traj_monomer/500ns_noWAT_5000frames_chain_", i, ".nc", sep=""))
  sim129$traj = c(sim129$traj, tmptraj)
}

sim116$rmsd <- rmsd( pdb1$xyz, sim116$traj, a.inds= ca.inds$xyz, b.inds=ca.inds$xyz, fit = TRUE )
sim129$rmsd <- rmsd( pdb1$xyz, sim129$traj, a.inds= ca.inds$xyz, b.inds=ca.inds$xyz, fit = TRUE )