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

##pdb <- read.pdb("1XCK_reference.pdb")
pdb <- read.pdb("1XCK_apo_dimer_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=30)


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

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


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


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


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


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


##a <- paste(s[inds[,1]], (inds[,1]%%524)+1, sep="")
##b <- paste(s[inds[,2]], (inds[,2]%%524)+1, 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("116_129_dcm.pdf")
plot.dcm(dmat, pdb, xlim=c(0,524), ylim=c(524,1048),
         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()