library(bio3d)
#load("xrayPCA.RData")
source("confplot_funs.R")
library(car)

#load("xrayPCA_Tstates.RData")


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

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


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

for ( i in 1:7 ) {
  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 = core$c1A.xyz,
                                mobile.inds = core$c1A.xyz,
                                full.pdbs = FALSE)
}

for ( i in 8:14 ) {
  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 = core$c1A.xyz,
                                mobile.inds = core$c1A.xyz,
                                full.pdbs = FALSE)
}

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

for ( i in 1:7 ) {
  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 = core$c1A.xyz,
                                mobile.inds = core$c1A.xyz,
                                full.pdbs = FALSE)
}

for ( i in 8:14 ) {
  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 = core$c1A.xyz,
                                mobile.inds = core$c1A.xyz,
                                full.pdbs = FALSE)
}


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


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


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


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

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

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