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

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


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

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

prefix <- "50ns_noWAT_5000frames_noH_chain_"
trj.inds <- seq(2001, 5000, by=15)


sim116 <- NULL
prefix116 <- "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/traj_monomer/"

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

sim129 <- NULL
prefix129 <- "/net/gulrotkake/slars/groel_md/1XCK/129_1XCK_MGATP/results/traj_monomer/"


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

for ( i in 8:14 ) {
  tmptrj <- read.ncdf(paste(prefix129, prefix, i , ".nc", sep=""))
  sim129$trj$free = rbind(sim129$trj$free, tmptrj[trj.inds,])
}


#sim116$dm <- dm.xyz.trj( sim116$trj, pdb,  threads=5 )
#sim129$dm$atp <- dm.xyz.trj( sim129$trj$atp, pdb, threads=5 )
#sim129$dm$free <- dm.xyz.trj( sim129$trj$free, pdb, threads=5 )

#save(sim116, sim129, file="sim116_129.RData")


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


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

pdf("116_129_dcm.pdf")
plot.dcm(dmat, pdb, xlim=c(0,524), ylim=c(0,524),
         sse=sse, sse.grid=TRUE,
         )
pointLabel(inds[,1], inds[,2], labels, cex=0.6, offset=10)
dev.off()



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


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


pdf("129free_116free_dcm.pdf")
plot.dcm(dmat, pdb, xlim=c(0,524), ylim=c(0,524),
         sse=sse, sse.grid=TRUE,
         )
pointLabel(inds[,1], inds[,2], labels, cex=0.6, offset=10)
dev.off()