library(ncdf)
library(bio3d)
source("confplot_funs.R")
library(car)

load("pca_dimers.RData")
#load("pca_dimers_t-only.RData")


pdb.closed <- read.pdb("1XCK_apo_dimer_noH.pdb")
ca.inds <- atom.select(pdb.closed, "calpha")

trj.inds <- seq(1, 1000, by=10)

##adapted.core <- c(core$c1A.xyz-3, core$c1A.xyz-1575)
adapted.core <- core.md


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

for ( i in 1:6 ) {
  tmptrj <- read.ncdf(paste(prefix116, prefix, i , ".nc", sep=""))
  sim116$xyzfit[[i]] <- fit.xyz(fixed = pdbs$xyz[1,gaps.pos$f.inds],
                                mobile = tmptrj[trj.inds, ca.inds$xyz],
                                fixed.inds = adapted.core,
                                mobile.inds = adapted.core, 
                                full.pdbs = FALSE)
}

for ( i in 8:13 ) {
  tmptrj <- read.ncdf(paste(prefix116, prefix, i , ".nc", sep=""))
  sim116$xyzfit[[i]] <- fit.xyz(fixed = pdbs$xyz[1,gaps.pos$f.inds],
                                mobile = tmptrj[trj.inds, ca.inds$xyz],
                                fixed.inds = adapted.core,
                                mobile.inds = adapted.core,
                                full.pdbs = FALSE)
}

sim129 <- NULL
prefix129 <- "/net/gulrotkake/slars/groel_md/1XCK/129_1XCK_MGATP/results/traj_dimers/"
prefix <- "50ns_noWAT_1000frames_noH_chain_chain_"

for ( i in 1:6 ) {
  tmptrj <- read.ncdf(paste(prefix129, prefix, i , ".nc", sep=""))
  sim129$xyzfit[[i]] <- fit.xyz(fixed = pdbs$xyz[1,gaps.pos$f.inds],
                                mobile = tmptrj[trj.inds, ca.inds$xyz],
                                fixed.inds = adapted.core, 
                                mobile.inds = adapted.core,
                                full.pdbs = FALSE)
}

for ( i in 8:13 ) {
  tmptrj <- read.ncdf(paste(prefix129, prefix, i , ".nc", sep=""))
  sim129$xyzfit[[i]] <- fit.xyz(fixed = pdbs$xyz[1,gaps.pos$f.inds],
                                mobile = tmptrj[trj.inds, ca.inds$xyz],
                                fixed.inds = adapted.core,
                                mobile.inds = adapted.core,
                                full.pdbs = FALSE)
}



sim129b <- NULL
prefix129b <- "/net/lutefisk/slars/groel_md/1XCK/129_1XCK_MGATP_temp350/results/traj_dimers/"
prefix <- "50-90ns_noWAT_1000frames_noATP_noH_chain_"
trj.inds <- seq(1, 1000, by=1)


for ( i in 1:6 ) {
  tmptrj <- read.ncdf(paste(prefix129b, prefix, i , ".nc", sep=""))
  sim129b$xyzfit[[i]] <- fit.xyz(fixed = pdbs$xyz[1,gaps.pos$f.inds],
                                mobile = tmptrj[trj.inds, ca.inds$xyz],
                                fixed.inds = adapted.core,
                                mobile.inds = adapted.core,
                                full.pdbs = FALSE)
}

for ( i in 8:13 ) {
  tmptrj <- read.ncdf(paste(prefix129b, prefix, i , ".nc", sep=""))
  sim129b$xyzfit[[i]] <- fit.xyz(fixed = pdbs$xyz[1,gaps.pos$f.inds],
                                mobile = tmptrj[trj.inds, ca.inds$xyz],
                                fixed.inds = adapted.core,
                                mobile.inds = adapted.core,
                                full.pdbs = FALSE)
}


sim171 <- NULL
prefix171 <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/traj_dimers/"
prefix <- "50ns_noWAT_1000frames_noATP_noH_chain_"

for ( i in 1:6 ) {
  tmptrj <- read.ncdf(paste(prefix171, prefix, i , ".nc", sep=""))
  sim171$xyzfit[[i]] <- fit.xyz(fixed = pdbs$xyz[1,gaps.pos$f.inds],
                                mobile = tmptrj[trj.inds, ca.inds$xyz],
                                fixed.inds = adapted.core,
                                mobile.inds = adapted.core,
                                full.pdbs = FALSE)
}

for ( i in 8:13 ) {
  tmptrj <- read.ncdf(paste(prefix171, prefix, i , ".nc", sep=""))
  sim171$xyzfit[[i]] <- fit.xyz(fixed = pdbs$xyz[1,gaps.pos$f.inds],
                                mobile = tmptrj[trj.inds, ca.inds$xyz],
                                fixed.inds = adapted.core,
                                mobile.inds = adapted.core,
                                full.pdbs = FALSE)
}




tmptrj <- NULL
for ( i in 1:6 ) {
  tmptrj=rbind(tmptrj, sim129b$xyzfit[[i]])
}
sim129b$cis$proj <- pca.project(tmptrj, pc.xray)


tmptrj <- NULL
for ( i in 8:13 ) {
  tmptrj=rbind(tmptrj, sim129b$xyzfit[[i]])
}
sim129b$trans$proj <- pca.project(tmptrj, pc.xray)

tmptrj <- NULL
for ( i in 1:6 ) {
  tmptrj=rbind(tmptrj, sim171$xyzfit[[i]])
}
sim171$cis$proj <- pca.project(tmptrj, pc.xray)


tmptrj <- NULL
for ( i in 8:13 ) {
  tmptrj=rbind(tmptrj, sim171$xyzfit[[i]])
}
sim171$trans$proj <- pca.project(tmptrj, pc.xray)



tmptrj <- NULL
for ( i in 1:6 ) {
  tmptrj=rbind(tmptrj, sim129$xyzfit[[i]])
}
sim129$cis$proj <- pca.project(tmptrj, pc.xray)


tmptrj <- NULL
for ( i in 8:13 ) {
  tmptrj=rbind(tmptrj, sim129$xyzfit[[i]])
}
sim129$trans$proj <- pca.project(tmptrj, pc.xray)


tmptrj <- NULL
for ( i in 1:6 ) {
  tmptrj=rbind(tmptrj, sim116$xyzfit[[i]])
}
sim116$cis$proj <- pca.project(tmptrj, pc.xray)


tmptrj <- NULL
for ( i in 8:13 ) {
  tmptrj=rbind(tmptrj, sim116$xyzfit[[i]])
}
sim116$trans$proj <- pca.project(tmptrj, pc.xray)

sim116$xyzfit <- NULL
sim129$xyzfit <- NULL

#sim116b$xyzfit <- NULL
sim129b$xyzfit <- NULL
sim171$xyzfit <- NULL


save(sim116, sim129, sim171, file="projected_50nsForTrajPlots.RData")


d <- matrix(sim116$cis$proj[,1], nrow=6, byrow=TRUE)
d2 <- matrix(sim116$trans$proj[,1], nrow=6, byrow=TRUE)

e <- matrix(sim129$cis$proj[,1], nrow=6, byrow=TRUE)
f <- matrix(sim171$cis$proj[,1], nrow=6, byrow=TRUE)


pdf("projectionsVsTime.pdf", w=12, h=18)
par(mfrow=c(6,1))

for ( i in 1:6 ) {
	plot(d[i,] , type='l', ylim=c(-70,100))
	lines(d2[i,], lty=3, col='black')
	lines(e[i,], lty=2, col='red')
        lines(f[i,], lty=2, col='green')
}
dev.off()



## find a conformer of interest
if(FALSE){
pdb.closed <- read.pdb("1XCK_apo_dimer_noH.pdb")

prefix129 <- "/net/gulrotkake/slars/groel_md/1XCK/129_1XCK_MGATP/results/traj_dimers/"
prefix <- "20-50ns_noWAT_1500frames_noATP_noH_chain_"
i=1
trj <- read.ncdf(paste(prefix129, prefix, i , ".nc", sep=""))
 

for ( i in 2:6 ) {
  tmptrj <- read.ncdf(paste(prefix129, prefix, i , ".nc", sep=""))
  trj = rbind(trj, tmptrj) 
}

frame.ind <- which.max(sim129$cis$proj[,2])

write.pdb(pdb=pdb.closed, xyz=trj[frame.ind,])

}