library(bio3d)
library(ncdf)


pdb.allatom <- read.pdb("1XCK_reference.pdb")          
pdb.noH <- read.pdb('1XCK_chainA_noWAT_noH.pdb')
pdb1 <- read.pdb("1XCK_chainA_noWAT_noH.pdb")
pdb2 <- read.pdb("1SVT_start_noH.pdb.1")
pdb3 <- read.pdb("2C7E_A.pdb_flsq.pdb")

inds.aa <- atom.select(pdb.allatom, "///9:515///CA/")
inds.noH <- atom.select(pdb.noH, "///9:515///CA/")
inds <- atom.select(pdb1, "///9:515///CA/")
ca.inds.2c7e <- atom.select(pdb3, "///10:516///CA/")


prefix116 <- "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/"
prefix129 <- "/net/gulrotkake/slars/groel_md/1XCK/129_1XCK_MGATP/results/"
prefix117 <- "/net/gulrotkake/slars/groel_md/1XCK/117_1XCK_hepta/results/"
prefix128 <- "/net/gulrotkake/slars/groel_md/1XCK/128_1XCK_hepta_MGATP/results/"
prefix171 <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/"
prefix129b <- "/net/lutefisk/slars/groel_md/1XCK/129_1XCK_MGATP_temp350/results/"

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

t <- seq(1,5000,20)


# 128: 50ns_noWAT_1000frames_noH_chain_1.nc
# 129: 50ns_noWAT_5000frames_noH_chain_
## monomers
for ( i in 1:7 ) {
  tmptrj = read.ncdf(paste(prefix116, "traj_monomer/50ns_noWAT_5000frames_noH_chain_", i, ".nc", sep=""))
  sim116$trj = rbind(sim116$trj, tmptrj[t,])
}

for ( i in 1:7 ) {
  tmptrj = read.ncdf(paste(prefix129, "traj_monomer/50ns_noWAT_5000frames_noH_chain_", i, ".nc", sep=""))
  sim129$trj = rbind(sim129$trj, tmptrj[t,])
}

for ( i in 1:7 ) {
  tmptrj = read.ncdf(paste(prefix171, "traj_monomer/50ns_noWAT_5000frames_noH_chain_", i, ".nc", sep=""))
  sim171$trj = rbind(sim171$trj, tmptrj[t,])
}

t <- seq(1,4000,20)
for ( i in 1:7 ) {
  tmptrj = read.ncdf(paste(prefix129b, "traj_monomer/500ns_noWAT_5000frames_noATP_noH_chain_", i, ".nc", sep=""))
  sim129b$trj = rbind(sim129b$trj, tmptrj[t,])
}


t2 <- seq(1, length(t)*7, 1)

sim116$rmsd.closed <- rmsd(pdb1$xyz, sim116$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
sim116$rmsd.open <- rmsd(pdb2$xyz, sim116$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
sim116$rmsd.r <- rmsd(pdb3$xyz, sim116$trj[t2,], a.inds=ca.inds.2c7e$xyz, b.inds=inds$xyz, fit=TRUE)

sim129$rmsd.closed <- rmsd(pdb1$xyz, sim129$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
sim129$rmsd.open <- rmsd(pdb2$xyz, sim129$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
sim129$rmsd.r <- rmsd(pdb3$xyz, sim129$trj[t2,], a.inds=ca.inds.2c7e$xyz, b.inds=inds$xyz, fit=TRUE)

sim129b$rmsd.closed <- rmsd(pdb1$xyz, sim129b$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
sim129b$rmsd.open <- rmsd(pdb2$xyz, sim129b$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
sim129b$rmsd.r <- rmsd(pdb3$xyz, sim129b$trj, a.inds=ca.inds.2c7e$xyz, b.inds=inds$xyz, fit=TRUE)


sim128$rmsd.r <- rmsd(pdb3$xyz, sim128$trj, a.inds=ca.inds.2c7e$xyz, b.inds=inds$xyz, fit=TRUE)

sim171$rmsd.closed <- rmsd(pdb1$xyz, sim171$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
sim171$rmsd.open <- rmsd(pdb2$xyz, sim171$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
sim171$rmsd.r <- rmsd(pdb3$xyz, sim171$trj[t2,], a.inds=ca.inds.2c7e$xyz, b.inds=inds$xyz, fit=TRUE)

pdf("rmsd_R_histogram.pdf", w=15, h=5)
par(mfrow=c(1,3))
hist(sim116$rmsd.r, xlim=c(3.5, 6.5), ylim=c(0,200))
hist(sim129$rmsd.r, xlim=c(3.5, 6.5), ylim=c(0,200))
hist(sim171$rmsd.r, xlim=c(3.5, 6.5), ylim=c(0,200))
dev.off()

w <- wilcox.test(sim116$rmsd.r, sim129$rmsd.r, paired = FALSE, alternative="greater")
t <- t.test(sim116$rmsd.r, sim129$rmsd.r, paired = FALSE, alternative="greater")

t.test(sim116$rmsd.r, sim129$rmsd.r)


sim116$trj=NULL
sim129$trj=NULL
sim171$trj=NULL

#save(sim116,sim129,sim171, file="rmsd.RData")
load("rmsd.RData")

rmsd(pdb1$xyz, pdb2$xyz, a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
rmsd(pdb1$xyz, pdb3$xyz, a.inds=inds$xyz, b.inds=ca.inds.2c7e$xyz, fit=TRUE)

max(sim116$rmsd.closed)
max(sim129$rmsd.closed)
max(sim171$rmsd.closed)

min(sim116$rmsd.r)
min(sim129$rmsd.r)
min(sim171$rmsd.r)

min(sim116$rmsd.open)
min(sim129$rmsd.open)
min(sim171$rmsd.open)



mean(sim116$rmsd.closed)
mean(sim129$rmsd.closed)
mean(sim171$rmsd.closed)

var(sim116$rmsd.r)
var(sim129$rmsd.r)
var(sim171$rmsd.r)

mean(sim116$rmsd.r)
mean(sim129$rmsd.r)
mean(sim171$rmsd.r)

mean(sim116$rmsd.open)
mean(sim129$rmsd.open)
mean(sim171$rmsd.open)


tmp <- matrix(sim116$rmsd.closed, nrow=7, byrow=T)
mean(tmp[,100:500])
var(as.vector(tmp[,100:500]))

tmp <- matrix(sim129$rmsd.closed, nrow=7, byrow=T)
mean(tmp[,100:500])
var(as.vector(tmp[,100:500]))


tmp <- matrix(sim171$rmsd.closed, nrow=7, byrow=T)
mean(tmp[,100:500])
var(as.vector(tmp[,100:500]))




pdf("rmsd_toAll.pdf", w=12, h=12)
par(mfcol=c(3,1))

tmp <- matrix(sim116$rmsd.closed, nrow=7, byrow=T)
plot(colMeans(tmp), type='l', col="blue", ylim=c(1, 3.2))
tmp <- matrix(sim129$rmsd.closed, nrow=7, byrow=T)
lines(colMeans(tmp), col="red")
tmp <- matrix(sim171$rmsd.closed, nrow=7, byrow=T)
lines(colMeans(tmp), col="black")


tmp <- matrix(sim116$rmsd.r, nrow=7, byrow=T)
plot(colMeans(tmp), type='l', col="blue", ylim=c(4.5, 6))
tmp <- matrix(sim129$rmsd.r, nrow=7, byrow=T)
lines(colMeans(tmp), col="red")
tmp <- matrix(sim171$rmsd.r, nrow=7, byrow=T)
lines(colMeans(tmp), col="black")



tmp <- matrix(sim116$rmsd.open, nrow=7, byrow=T)
plot(colMeans(tmp), type='l', col="blue", ylim=c(10.5, 12.5))
tmp <- matrix(sim129$rmsd.open, nrow=7, byrow=T)
lines(colMeans(tmp), col="red")
tmp <- matrix(sim171$rmsd.open, nrow=7, byrow=T)
lines(colMeans(tmp), col="black")



dev.off()

pdf("rmsd_toR.pdf", w=12, h=6)
par(mfcol=c(1,1))
tmp <- matrix(sim116$rmsd.r, nrow=7, byrow=T)
plot(colMeans(tmp), type='l', col="blue", ylim=c(3.5, 6))
for ( i in 1:nrow(tmp) ) {
	lines(tmp[i,], col="grey50", lty=2)
}

tmp <- matrix(sim129$rmsd.r, nrow=7, byrow=T)
lines(colMeans(tmp), col="red")
for ( i in 1:nrow(tmp) ) {
        lines(tmp[i,], col="grey50", lty=3)
}

tmp <- matrix(sim171$rmsd.r, nrow=7, byrow=T)
lines(colMeans(tmp), col="black")

tmp <- matrix(sim129b$rmsd.r, nrow=7, byrow=T)
lines(colMeans(tmp), col="green")

dev.off()