## Usage: Rscript --vanilla write_dcd2pdb.R pdbfile dcdfile time_stepfile top_percent
## the output will be a text file of timestep and cc, as well as 
## top percent of the pdbs written out

#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
#print(args[1])
#length(args)
## test if there is at least one argument: if not, return an error
if (length(args)!=4) {
  stop("Please provide the pdb, dcd, timestep file, and top_percent", call.=FALSE)
}
  
library(dplyr)
library(stringr)
library(bio3d)

#file = "/Users/bhandariyr/Documents/projects/afm/cafemol/dimer/twochain/run2_rescaled/output1/cbl_afm.ts"
file = args[3]
data.ts.lines = readLines(file)[grep("#all", readLines(file))]
data.ts =  data.ts.lines %>%
              strsplit("\\s+") %>%
              do.call(rbind, .) %>%
              data.frame()
dim(data.ts)
#View(data.ts)
top.percent = as.numeric(args[4])
afmcc.file = paste("top_", top.percent, "percent_cc.txt", sep = "")
top.frames = round(top.percent*nrow(data.ts)/100) 
data.ts.sort = data.ts[order(data.ts$X15, decreasing = TRUE), ]
data.ts.sort$frame = as.numeric(rownames(data.ts.sort)) - 1
data.ts.sort = data.ts.sort[1:top.frames,c("frame", "X2", "X15")]
colnames(data.ts.sort) = c("frame","time.step", "afmcc")
rownames(data.ts.sort) = 1:nrow(data.ts.sort)
write.table(format(data.ts.sort), file = afmcc.file, sep = "\t", row.names = FALSE, quote = FALSE)

# dcdfile = "~/Documents/projects/afm/cafemol/dimer/twochain/run2_rescaled/output1/cbl_afm.dcd"
# pdbfile = "~/Documents/projects/afm/cafemol/dimer/twochain/run1/pdb/cbl_dimer_offset.pdb"
pdbfile = args[1]
dcdfile = args[2]
dcd <- read.dcd(dcdfile)
pdb <- read.pdb(pdbfile)
ca.inds <- atom.select(pdb)
xyz <- fit.xyz(fixed=pdb$xyz, mobile=dcd,
               fixed.inds=ca.inds$xyz,
               mobile.inds=ca.inds$xyz)
dim(xyz) == dim(dcd)
for (val in data.ts.sort$frame){
  if (val != 0) {
  pdb.name = paste("best_cc_frame_",val, ".pdb", sep = "")
  write.pdb(pdb, xyz[val, ], file=pdb.name)
  }
}

