#pdbcenter_npts_finder.R calculates the geometric center of a PDB file

#!/usr/bin/env Rscript      
args = commandArgs(trailingOnly=TRUE) 

if(length(args) < 1){
  stop("USE: Rscript pdbcenter_npts_finder.R <receptor_name> <session_dir>")
}

#Arguments definition 
receptor_name <- args[1]
session_dir <- args[2]

#Read PDB file 
archivo <- scan(args[1], what = character(), quiet = TRUE)

#Make a table where positions of carbon (C) atoms and their x, y, z coordenates will be placed
salto <- 0
for(i in 1:(length(archivo)-2)){
  if(archivo[i] == "ATOM" && archivo[i+2] == "C"){
    salto <- salto + 1
  }
}
tablaC <- matrix(1, nrow = salto, ncol = 4)

#Header. Note: first column: C atom position (will be filled in when the number of C will be known)
#Row labels
colnames(tablaC) <- c("C Position","X", "Y", "z")
rownames(tablaC) <- 1:salto

#Fill in tablaC with x, y,z coordenates of the C atoms as numeric
salto <- 0
posicionC <-  c()
for(i in 1:(length(archivo)-2)){
  if(archivo[i] == "ATOM" && archivo[i+2] == "C"){
    salto <- salto + 1
    posicionC <- c(posicionC, archivo[i+1])
    
    tablaC[salto,2] <- as.numeric(archivo[i+6])
    tablaC[salto,3] <- as.numeric(archivo[i+7])
    tablaC[salto,4] <- as.numeric(archivo[i+8])
 }
}

#Add C atom positions in the column 1
for(i in 1:salto){
  tablaC[i,1] <- as.numeric(posicionC[i]) 
}

#Geometric center calculation as grid center x, y, z to be introduced in the GPF file 
pto_x_medio <- (max(tablaC[,2]) + min(tablaC[,2]))/2
pto_y_medio <- (max(tablaC[,3]) + min(tablaC[,3]))/2
pto_z_medio <- (max(tablaC[,4]) + min(tablaC[,4]))/2

pto_x_medio_r <- round(pto_x_medio,digits=3)
pto_y_medio_r <- round(pto_y_medio,digits=3)
pto_z_medio_r <- round(pto_z_medio,digits=3)

#print(paste(pto_x_medio_r, pto_y_medio_r, pto_z_medio_r, sep = ";"))

#Calculate distances in each coordenate to have de number of points (npts) x, y, z to be introduced in the GPF file 
eu_x <- max(tablaC[,2])-min(tablaC[,2])
eu_y <- max(tablaC[,3])-min(tablaC[,3])
eu_z <- max(tablaC[,4])-min(tablaC[,4])

eu_x_r <- round(eu_x,digits=0)
eu_y_r <- round(eu_y,digits=0)
eu_z_r <- round(eu_z,digits=0)

#print(paste(eu_x, eu_y, eu_z, sep = ";"))
#print(paste(eu_x_r, eu_y_r, eu_z_r, sep = ";"))

if(eu_x > 126 || eu_y > 126 || eu_x > 126){
  print("At least one dimension in number of points is out of range (>126) according AutoDock software")
  break
}

#Write PDB geometric center and eu_x_r, eu_y_r, eu_z_r in a file
#write(paste(pto_x_medio_r, pto_y_medio_r, pto_z_medio_r, eu_x_r, eu_y_r, eu_z_r,sep= "\n"), file = "/home/galaxy/galaxy/tools/proteindocking/scripts/pdbcenter_npts.txt")
write(paste(pto_x_medio_r, pto_y_medio_r, pto_z_medio_r, eu_x_r, eu_y_r, eu_z_r,sep= "\n"), file = "pdbcenter_npts.txt")