library(bio3d)
library(ncdf)

pdb.open <- read.pdb("/net/gulrotkake/slars/groel_nma_korgpil/1SVT/chainA/start.pdb")
pdb.closed <- read.pdb("1XCK_reference.pdb")
ca.inds <- atom.select(pdb.open, "calpha")

xyz <- fit.xyz(pdb.open$xyz, pdb.closed$xyz, fixed.inds = ca.inds$xyz, mobile.inds = ca.inds$xyz,  pdb.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
pdb.closed$xyz = xyz
write.pdb(pdb.closed, xyz=xyz, file="moved.pdb")

diff <- as.numeric(pdb.open$xyz[ca.inds$xyz] - pdb.closed$xyz[ca.inds$xyz])
diff=normalizedVector(diff)

pc <- scan("200ns_evecs.dat")
pc = matrix(d, ncol=51, byrow=F)


#trj0 <- read.ncdf("../traj/200ns_short_2000frames.nc")
trj0 <- read.ncdf("../traj/300ns_short_15000frames_noWAT.nc")

xyz0 <- fit.xyz(fixed = pdb.open$xyz[ca.inds$xyz], 
                mobile =  rbind(pdb.open$xyz[ca.inds$xyz],
                  pdb.closed$xyz[ca.inds$xyz], trj0[1:12500, ca.inds$xyz]), 
                pdb.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)



pc.proj <- (xyz0 - pc[,1]) %*% pc[,2]
pc.proj <- cbind(pc.proj, (xyz0 - pc[,1]) %*% pc[,3])
pc.proj <- cbind(pc.proj, (xyz0 - pc[,1]) %*% pc[,4])


diff.proj <- (xyz0 - pc[,1]) %*% (diff)





pdf("projection2.pdf", h=3, w=9)
par(mfrow=c(1,3), mar=c(4,4,2,2))
plot(pc.proj[,1], diff.proj, type="p", main = "", xlab="PC1", ylab="Conf. Diff.",
     col=densCols(pc.proj[,1], diff.proj), pch=20, cex=1.5 )
points(pc.proj[3:300,1], diff.proj[3:300], col = "green", pch = 16 )
points(pc.proj[1:2,1], diff.proj[1:2], col = c("red", "orange"), pch = 16, cex.axis=1.5 )
mtext("A", side=3, line=-1.7, at=0, cex=1.3)


plot(pc.proj[,1], pc.proj[,2], type="p", main = "", xlab="PC1", ylab="PC2",
     col=densCols(pc.proj[,1], pc.proj[,2]), pch=20 )
points(pc.proj[3:300,1], pc.proj[3:300,2], col = "green", pch = 16 )
points(pc.proj[1:2,1], pc.proj[1:2,2], col = c("red", "orange"), pch = 16, cex=1.5 )
mtext("B", side=3, line=-1.7, at=0, cex=1.3)


plot(pc.proj[,1], pc.proj[,3], type="p", main = "", xlab="PC1", ylab="PC3",
     col=densCols(pc.proj[,1], pc.proj[,3]), pch=20 )
points(pc.proj[3:300,1], pc.proj[3:300,3], col = "green", pch = 16 )
points(pc.proj[1:2,1], pc.proj[1:2,3], col = c("red", "orange"), pch = 16, cex=01.5 )
mtext("C", side=3, line=-1.7, at=0, cex=1.3)


dev.off()





     









ptraj.proj <- read.table("../projection/project.out")

pdf("tmp.pdf", w=12)
plot(ptraj.proj[,"Mode1"], ptraj.proj[,"Mode2"])
plot(pc.proj[,1], pc.proj[,2])
dev.off()



"normalizedVector" <-
  function(v) {
    return( v/sqrt(dotProduct(v,v)) )
  }