library(gclus)
library(maptools)
library(bio3d)
source("confplot_funs.R")

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


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

pdb.e434k <- read.pdb("3FBH_chainAB_noWAT_noH.amb.pdb.1")
#pdb.e461k <- read.pdb("E461K_apo_chainA_noWAT_noH.pdb")
ca.inds.e434k <- atom.select(pdb.e434k, "calpha")
#ca.inds.e461k <- atom.select(pdb.e461k, "calpha")

trj.inds <- seq(1, 1500, by=1)
if(FALSE) {
sim161 <- NULL
prefix161 <- "/net/lutefisk/slars/groel_md/mutants/161_E434K_apo/results/traj_dimers/"
prefix <- "20-50ns_noWAT_1500frames_noH_chain_"

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

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

sim162 <- NULL
prefix162 <- "/net/lutefisk/slars/groel_md/mutants/162_E434K_MGATP/results/traj_dimers/"
prefix <- "20-50ns_noWAT_1500frames_noH_chain_"
trj.inds <- seq(1, 1500, by=1)

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

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


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

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

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


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

#sim161$xyzfit <- NULL
sim162$xyzfit <- NULL
save(sim162, file="projected_mutant_E434K.RData")