library(gclus)
library(maptools)
library(bio3d)
source("confplot_funs.R")
#load("xrayPCA.RData")
load("xrayPCA_Tstates.RData")


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


      

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

sim171 <- NULL
prefix171 <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/traj_monomer/"

prefix <- "100ns_noWAT_2000frames_noATP_noH_chain_"

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

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



### THIS IS FOR 20-50 ns only ###
sim171b <- NULL
prefix171b <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/traj_monomer/"

prefix <- "20-50ns_noWAT_1500frames_noATP_noH_chain_"

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

for ( i in 1:7 ) {
  tmptrj <- read.ncdf(paste(prefix171b, prefix, i , ".nc", sep=""))
  sim171b$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(prefix171b, prefix, i , ".nc", sep=""))
  sim171b$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)
}

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

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

## 1XCK RHOD 20-50
tmptrj <- NULL
for ( i in 1:7 ) {
  tmptrj=rbind(tmptrj, sim171b$xyzfit[[i]])
}
sim171b$cis$proj <- pca.project(tmptrj, pc.xray)

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


sim171$xyzfit <- NULL
sim171b$xyzfit <- NULL
save(sim171, sim171b, file="projected_1XCK-rhod.RData")