## Process the CG model of PDB file generated by cafemol initial run
## Does coordiante centering based on cafemol input afm and pdb formatting 
## for preparing the pdb file for Cafemol AFM fit run
## Usage: Rscript process_pdb_cg.R input.pdb cafemol_input_file

# install.packages("bio3d")
library(bio3d)
library(dplyr)
library(stringr)

#!/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)==0) {
  stop("Please provide the initial PDB file", call.=FALSE)
}

## Read input pdb
input.pdb = args[1]
#input.pdb = "/Users/bhandariyr/Documents/projects/afm/cafemol/miscen/scripts/centering/4gma_amber_offset.pdb"
pdb = read.pdb2(input.pdb)

## Read cafemol input file
file = args[2]
#file = "/Users/bhandariyr/Documents/projects/afm/cafemol/miscen/scripts/centering/cmpt1_5A_resample_init_cafemol.inp"
data.lines = readLines(file)[grep("AFM_IMAGE", readLines(file))]
data.inp =  data.lines %>%
  strsplit("\\s+") %>%
  do.call(rbind, .) %>%
  data.frame()

# head(data.inp)
x.CM = sum(as.numeric(data.inp$X2))/nrow(data.inp)
y.CM = sum(as.numeric(data.inp$X3))/nrow(data.inp)

# class(pdb)
# dim(pdb$xyz)
#min(pdb$atom$x)
x.pdb.CM = sum(pdb$atom$x)/length(pdb$atom$x)
y.pdb.CM = sum(pdb$atom$y)/length(pdb$atom$y)
# x.CM; x.pdb.CM
# y.CM; y.pdb.CM

pdb$atom$x =  pdb$atom$x + (x.CM - x.pdb.CM)
pdb$atom$y =  pdb$atom$y + (y.CM - y.pdb.CM)
pdb$atom$z =  pdb$atom$z - min(pdb$atom$z)
#pdb$atom$z =  pdb$atom$z - min(pdb$atom$z) + 2
x.pdb.CM = sum(pdb$atom$x)/length(pdb$atom$x)
y.pdb.CM = sum(pdb$atom$y)/length(pdb$atom$y)
z.pdb.CM = sum(pdb$atom$z)/length(pdb$atom$z)
# x.pdb.CM; y.pdb.CM; z.pdb.CM

xyz.mat = cbind(pdb$atom$x, pdb$atom$y, pdb$atom$z )
# class(xyz.mat)
xyz.vec = as.vector(t(xyz.mat))
## set the spacing with respect to resid
pdb$atom$resid = paste0(pdb$atom$resid, " ")
# write the new pdb object to file
#write.pdb(pdb, xyz=xyz, file="/Users/bhandariyr/Documents/projects/afm/cafemol/dimer/run1/cbl_dimer_frame1000_offset.pdb")
out.pdb = paste0(tools::file_path_sans_ext(basename(input.pdb)),"_center.pdb")
write.pdb(pdb, xyz=xyz.vec, 
          file=out.pdb)

# pdb$atom[1:4,]
# scp -pr /Users/bhandariyr/Documents/projects/afm/cafemol/test_run/cbl_run1_offset.pdb bhandariyr@biowulf.nih.gov:/gpfs/gsfs12/users/bhandariyr/cafemol/cafemol/example/cbl_afm_fit/pdb





