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/"

prefix150 <-"/net/lutefisk/slars/groel_md/mutants/150_E461K_apo/results/rmsd_monomer/"
prefix151 <- "/net/lutefisk/slars/groel_md/mutants/151_E461K_ATP/results/rmsd_monomer/"

prefix161 <- "/net/lutefisk/slars/groel_md/mutants/161_E434K_apo/results/rmsd_monomer/"
prefix162 <- "/net/lutefisk/slars/groel_md/mutants/162_E434K_MGATP/results/rmsd_monomer/"


sim116 <- NULL
sim129 <- NULL

sim150 <- NULL
sim151 <- NULL

sim161 <- NULL
sim162 <- 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:14 ) {
  tmprmsd <- read.table(paste(prefix150, "50ns_5000frames_RMSD_", i , ".dat", sep=""))
  sim150$rmsd=rbind(sim150$rmsd, tmprmsd[,"V2"])
}

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

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

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


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




pdf("rmsd_monomers_mutants.pdf", w=10, h=5)
plot(sim151$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(sim151$rmsd[i,], lty=2, col="grey80")
}

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

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

for ( i in 1:14 ) {
  lines(sim162$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(sim150$rmsd.mean, lty=2, col="green")

lines(sim151$rmsd.mean$free, lty=2, col="green")
lines(sim151$rmsd.mean$atp, lty=2, col="black")

lines(sim161$rmsd.mean, lty=2, col="green")
lines(sim162$rmsd.mean$free, lty=2, col="green")
lines(sim162$rmsd.mean$atp, lty=2, col="black")


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", "WT ATP", "WT ATP free", "E434K ATP", "E461K ATP", "Mut. ATP Free", "Apo"),
       col=c("grey80", "blue", "red", "black", "black", "green", "green"), 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 )