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, 1500, by=1)

##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 <- "20-50ns_noWAT_1500frames_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 <- "20-50ns_noWAT_1500frames_noATP_noH_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)
}


sim116b <- NULL
prefix116b <- "/net/lutefisk/slars/groel_md/1XCK/116_1XCK_apo_temp350/results/traj_dimers/"
prefix <- "50-90ns_noWAT_1000frames_noATP_noH_chain_"
trj.inds <- seq(1, 955, by=1)


for ( i in 1:6 ) {
  tmptrj <- read.ncdf(paste(prefix116b, prefix, i , ".nc", sep=""))
  sim116b$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(prefix116b, prefix, i , ".nc", sep=""))
  sim116b$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, sim116b$xyzfit[[i]])
}
sim116b$cis$proj <- pca.project(tmptrj, pc.xray)


tmptrj <- NULL
for ( i in 8:13 ) {
  tmptrj=rbind(tmptrj, sim116b$xyzfit[[i]])
}
sim116b$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

save(sim116, sim129, sim116b, sim129b, file="projected.RData")





## 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,])

}