library(maptools)
library(rparallel)
library(bio3d)
library(ncdf)
source("dmat_funs.R")

##pdb <- read.pdb("1XCK_reference.pdb")
pdb <- read.pdb("trimer_noH.pdb.1")
pdb.ref <- pdb


seq <- seq.pdb(pdb)
s <- array(seq)
#sse <- dssp(pdb)
sse <- stride(pdb)


ca.inds <- atom.select(pdb, "calpha")
dim <- length(ca.inds$atom)

prefix <- "50ns_noWAT_1000frames_noADP_noH_chain_"
trj.inds <- seq(801, 1000, by=5)

sim180 <- NULL
prefix180 <- "/net/gulrotkake/slars/groel_md/1SX4/180_1SX4/results/traj_trimer/"

for ( i in 2:7 ) {
  tmptrj <- read.ncdf(paste(prefix180, prefix, i , ".nc", sep=""))
  sim180$trj = rbind(sim180$trj, tmptrj[trj.inds,])
}

sim181 <- NULL
prefix181 <- "/net/gulrotkake/slars/groel_md/1SX4/181_1SX4_atp/results/traj_trimer/"

for ( i in 2:7 ) {
  tmptrj <- read.ncdf(paste(prefix181, prefix, i , ".nc", sep=""))
  sim181$trj = rbind(sim181$trj, tmptrj[trj.inds,])
}

#sim180$rmsd <- rmsd(pdb$xyz, b=sim180$trj, a.inds=ca.inds$xyz, b.inds=ca.inds$xyz, fit=TRUE)
#sim181$rmsd <- rmsd(pdb$xyz, b=sim181$trj, a.inds=ca.inds$xyz, b.inds=ca.inds$xyz, fit=TRUE)
#plot(sim180$rmsd)
#lines(sim181$rmsd)
#dev.off()

sim180$dm <- dm.xyz.trj3( sim180$trj, pdb,  threads=6 )
sim181$dm <- dm.xyz.trj3( sim181$trj, pdb, threads=6 )


#save(sim180, sim181, file="sim180_181.RData")
load("sim180_181.RData")


dmat <- dcm(sim180$dm, sim181$dm, occupancy=0.5, diff.cut=0)
inds <- which(dmat != 0, arr.ind=T)


#a <- paste(s[inds[,1]], ((inds[,1]-1)%%524)+2, sep="")
#b <- paste(s[inds[,2]], ((inds[,2]-1)%%524)+2, sep="")
a <- paste(s[inds[,1]], (inds[,1])+1, sep="")
b <- paste(s[inds[,2]], (inds[,2])+1, sep="")
labels <- paste(a, b, sep="-")


pdf("180_181_dcm.pdf")
plot.dcm(dmat, pdb, xlim=c(1,1572), ylim=c(1,1572),
         sse=sse, sse.grid=TRUE,
         )
pointLabel(inds[,1], inds[,2], labels, cex=0.3, offset=10)
dev.off()

rows=c(1:5)
tmp1=tail(diff.table1, n=1050)
pointLabel(tmp1[,"row"], tmp1[,"col"],
           paste(tmp1[,"res1"],tmp1[,"res2"],sep="-"),
           cex=0.3, offset=10)
tmp2=head(diff.table1, n=0)
pointLabel(tmp2[,"row"], tmp2[,"col"],
           paste(tmp2[,"res1"],tmp2[,"res2"],sep="-"),
           cex=0.3, offset=10)

abline(v=524, lty=2)
abline(h=524, lty=2)


dev.off()





md.avg.diff <- sim116$dm$dmat.avg - sim129$dm$dmat.avg


diff <- which((md.avg.diff!=0) & (dmat!=0), arr.ind=T )
#diff <- which((md.avg.diff!=0), arr.ind=T )
diff = cbind(diff, md.avg.diff[diff])
diff = data.frame(diff)
diff = cbind(diff,  paste( s[diff[,"row"]], ((diff[,"row"]-1)%%524)+2, sep="" ) )
diff = cbind(diff,  paste( s[diff[,"col"]], ((diff[,"col"]-1)%%524)+2, sep="" ) )
colnames(diff)<-c("row", "col", "avg.diff", "res1", "res2")
diff.table1 <- diff[abs(order(diff$avg.diff)),]

head(diff.table1)


t <- which(diff.table1[,"row"]<525 & diff.table1[,"col"]>524)
u <- which(diff.table1[,"col"]<525 & diff.table1[,"row"]>524)
t=c(t,u)