## Calculate top structures based on kappa filtering
## Usage:
## Rscript /data/PNAI/yuba/scripts/analysis/get_stats_3D.R kappa_file N cafemol_pdb rmsd_pdb atom_type F_cutoff M_cutoff
## F_cutoff is the no. of frames to cutoff for frame filtering, by default 2000 frames are cut off for the analysis

args = commandArgs(trailingOnly=TRUE)
if (length(args) == 0) {
  cat("Calculate top structures based on kappa filtering and frame filtering \n")
  cat("Kappa filtering is done based on linear regression and outlier removal by Cooks Distance\n")
  cat("Frame filtering is based on sorting by Etot*CC^N or EstackEhbond*CC^N or EstackEhbondErepulsion*CC^N\n")
  cat("Usage:\n")
  cat("       Rscript get_gbStats_scan3D_v11.R scan3Dlist.txt N pdb_cafemol pdb_rmsd atom_type F_cutoff M_cutoff K_input N_kappa N_frame\n\n")
  cat("       scan3Dlist.txt file with list of 3D scan directories values\n")
  cat("       N              power to which afmCC is raised while multipying with Energy Terms\n")
  cat("       pdb_cafemol    pdb file listed in the cafemol_input file\n")
  cat("       pdb_rmsd       reference pdb file for calculating the rmsd against\n")
  cat("       atom_type      atom type from the list (sugar, phosphate, base, all) for rmsd calculation \n")
  cat("       F_cutoff       integer value representing the # of beginning frames to cutoff for frame filtering\n")
  cat("       M_cutoff       integer value representing the number of frames to include for calculation summary statistics for each kappa\n")
  cat("       K_input        input kappa value upon which the filtering of XY and Z pixel combinations is carried\n")
  cat("       N_kappa        integer value representing the number of top kappa runs to write out\n")
  cat("       N_frame        integer value representing the number of top PDB frames to write out under each kappa\n")
  cat("       Calculates mean and std dev summary statistics based on lowest 10 energy terms sorted by Etot - Afmfit \n")
  cat("       Output files are in the folder top_statsStHb11_PowN, frame filtering results are under subdirectory topPDB_Etot, topPDB_StckHb and topPDB_StckHbRepul\n") 
  cat("       And also under the output directory of each kappa/scan3D directory (Escore_ccPowN, FrameFilter_ccPowN) \n") 
  quit()
}


#install.packages("ggplot2")
#library(gdata)
library(dplyr)
library(stringr)
library(bio3d)
#install.packages("gridExtra")
library(gridExtra)
library(ggplot2)


kfile = args[1]
PCC = as.integer(args[2])
pdb.cafe = args[3]
pdb.rmsd = args[4]
F_cutoff = as.integer(args[6])
M_cutoff = as.integer(args[7])
K_input = args[8]
N_kappa = as.integer(args[9])
N_frame = as.integer(args[10])
conn = file(kfile,open="r")
kvalues = readLines(conn)
input.folders = kvalues
#ind = kvalues
#input.folders = c(unlist(strsplit(str_sub(args[1], 2, -2),",")))
#ind <- c(unlist(strsplit(str_sub(args[2], 2, -2),",")))

pdb1 = read.pdb(pdb.cafe)
pdb2 = read.pdb(pdb.rmsd)
atom.type = as.character(args[5])
if (atom.type == "sugar") {
   ca.inds <- atom.select(pdb1, elety = "S")
}   else if (atom.type == "phosphate") {
    ca.inds <- atom.select(pdb1, elety = "P")
}  else if (atom.type == "base") {
    ca.inds <- atom.select(pdb1, elety = c("Ab","Ub", "Gb","Cb"))
}  else if (atom.type == "all") {
    ca.inds <- atom.select(pdb1)
}

#setwd("/Users/bhandariyr/Documents/projects/afm/cobalamin/miscen/scripts/stats/contact10_5")
#top_percent = as.numeric(args[2])
# print(paste("top percent is", top_percent,sep = ": "))

## initialize data frames for each energy term
df_top_all = data.frame(kappa = character(), frame = integer(), cc = numeric(), rmsd = numeric(), diffEn.CCscaled = numeric())
df_etot_tot = data.frame(kappa = character(), mean = numeric(), sd = numeric())
df_local_tot = data.frame(kappa = character(), mean = numeric(), sd = numeric())
df_go_tot = data.frame(kappa = character(), mean = numeric(), sd = numeric())
df_afmfit_tot = data.frame(kappa = character(), mean = numeric(), sd = numeric())
df_stack_tot = data.frame(kappa = character(), mean = numeric(), sd = numeric())
df_hbond_tot = data.frame(kappa = character(), mean = numeric(), sd = numeric())
df_repul_tot = data.frame(kappa = character(), mean = numeric(), sd = numeric())
df_elect_tot = data.frame(kappa = character(), mean = numeric(), sd = numeric())
df_afmcc_tot = data.frame(kappa = character(), mean = numeric(), sd = numeric())
df_EleRepul_tot = data.frame(kappa = character(), mean = numeric(), sd = numeric())
df_StackHbond_tot = data.frame(kappa = character(), mean = numeric(), sd = numeric())

## make a top level directory to output the results
out.dir = paste("top_statsStHb11_Pow",PCC, sep = "")
out.pdbdir1 = paste(out.dir, "topPDB_Etot", sep = "/")
out.pdbdir2 = paste(out.dir, "topPDB_StckHb", sep = "/")
out.pdbdir3 = paste(out.dir, "topPDB_StckHbRepul", sep = "/")

out.pdblist1 = paste(out.dir, "topPDB_Etot.txt", sep = "/")
out.pdblist2 = paste(out.dir, "topPDB_StckHb.txt", sep = "/")
out.pdblist3 = paste(out.dir, "topPDB_StckHbRepul.txt", sep = "/")
if (dir.exists(out.dir)) {
  unlink(out.dir, recursive = TRUE)
}
dir.create(out.dir)
dir.create(out.pdbdir1)
dir.create(out.pdbdir2)
dir.create(out.pdbdir3)

## Loop through the given folders
# folders = list("output1_0","output1_1", "output1_2","output1_3","output1_4","output1_5",
#                "output1_6","output1_7","output1_8","output1_9","output1_10")

#path = "/data/PNAI/yxwang/CafeMol/cafe_runs/trRNaseP/experimental/f5511199d-10pc-1m-1A/"
## Define path if it exists
if (length(args) == 11) {
  path = args[11]
}

for (i in 1:(length(input.folders))){

## check if path exists
if (exists("path")){
ext = paste(path, input.folders[i], sep = "/")
} else {
ext = input.folders[i]
}
ext1 = paste(ext,"output", "*.ts", sep = "/")
print(paste0("Now reading ", ext1))
file = Sys.glob(ext1)
#print(file)
data.ts.lines = readLines(file)[grep("#all", readLines(file))]
#print("Hi I am here")
data.ts =  data.ts.lines %>%
  strsplit("\\s+") %>%
  do.call(rbind, .) %>%
  data.frame()
dim(data.ts)
colnames(data.ts) = unlist(strsplit(readLines(file)[grep("#unit", readLines(file))], "\\s+"))
#head(data.ts)
colnames(data.ts)[12] = "stack"
# View(data.ts)
data.ts$frame = as.numeric(rownames(data.ts)) - 1

data.ts.sort <- data.ts %>% select("frame", everything())
## subset to remove first 2000 frames
data.ts.sort = data.ts.sort[(F_cutoff+1):nrow(data.ts.sort),]
# View(data.ts.sort)

## Check for the existence of unstable calculation by looking at last Etot line of timestep file
#print(data.ts.sort[nrow(data.ts.sort),c("etot","afmcc", "afmfit") ])
last.Enstatus = any(is.na(as.numeric(data.ts.sort$etot))) 
if (last.Enstatus) {
  print(paste0("found NA in kappa: ", input.folders[i]))
  print(paste("skipping", input.folders[i], "because of unstable cafemol calculation", sep = " "))
  next #go to next iteration
}
## print the energy corresponding to top X percent afmcc sturctures
#top_percent = as.numeric(args[2])

#print(colnames(data.ts.sort))
#print(head(data.ts.sort))
## calculate correlation between afmcc and nromCCxEtot
afmcc = as.numeric(data.ts.sort$afmcc) #Afmcc score
toten = as.numeric(data.ts.sort$etot)  #total energy
afmfit = as.numeric(data.ts.sort$afmfit) 

#data.ts.sort$En.score = (toten-afmfit)*(afmcc^PCC)
data.ts.sort$diffEn = (toten-afmfit)
data.ts.sort$diffEn.CCscaled = (toten-afmfit)*(afmcc^PCC)

data.ts.rank = data.ts.sort[order(data.ts.sort$diffEn.CCscaled, decreasing = FALSE), ]
#print(head(data.ts.rank))

## Write the En.score in output file under each kappa/3Dscan folder
Escore.ifile = paste(input.folders[i], "output", paste0("Escore_ccPow", PCC, ".txt"), sep = "/")
write.table(format(data.ts.rank), file = Escore.ifile, sep = "\t", row.names = FALSE, quote = FALSE)
#print(head(data.ts.rank))

## write the pdb file
#print("Now extracting the top PDB from dcd")
ext2 = paste(ext,"output", "*.dcd", sep = "/")
dcdfile = Sys.glob(ext2)
#dcd <- read.dcd(dcdfile)
id.kappa = input.folders[i]
id.frame = as.integer(data.ts.rank[1,"frame"])
#out.pdb = paste(out.pdbdir,"/",id.kappa,"_",id.frame,".pdb", sep="")
#if (id.frame !=0){
#  write.pdb(pdb1, dcd[id.frame, ], file = out.pdb)
#} else {
#    print("Frame id is 0, unable to extract PDB file")
#}

## write rmsd
#pdb3 = read.pdb(out.pdb)
#rmsd.val = rmsd(pdb2$xyz[, ca.inds$xyz], pdb3$xyz[, ca.inds$xyz], fit = TRUE)

df_top = data.frame(kappa = input.folders[i], frame = data.ts.rank[1,"frame"], cc = data.ts.rank[1,"afmcc"],  diffEn.CCscaled = data.ts.rank[1,"diffEn.CCscaled"])
df_top_all = rbind(df_top_all, df_top)

#cc_norm = (afmcc - min(afmcc))/(max(afmcc) - min(afmcc))
#ccNormxEtot = cc_norm*toten

## Now carry out summary stats based on each energy term
data.ts.sort$EleRepul = (as.numeric(data.ts.sort$elect) + as.numeric(data.ts.sort$repul))
data.ts.sort$StackHbond = (as.numeric(data.ts.sort$stack) + as.numeric(data.ts.sort$hbond)) 
data.ts.rank = data.ts.sort[order(data.ts.sort$diffEn, decreasing = F), ]
data.rank10 = data.ts.rank[1:M_cutoff, ]
#print(data.rank10)

df_etot = data.frame(kappa = input.folders[i], mean = mean(as.numeric(data.rank10$diffEn)), sd = sd(as.numeric(data.rank10$diffEn)))
df_etot_tot = rbind(df_etot_tot, df_etot)
#print(paste(mean(as.numeric(data.rank10$diffEn)), sd(as.numeric(data.rank10$diffEn)), sep = " "))

df_local = data.frame(kappa = input.folders[i], mean = mean(as.numeric(data.rank10$local)), sd = sd(as.numeric(data.rank10$local)))
df_local_tot = rbind(df_local_tot, df_local)

df_go = data.frame(kappa = input.folders[i], mean = mean(as.numeric(data.rank10$go)), sd = sd(as.numeric(data.rank10$go)))
df_go_tot = rbind(df_go_tot, df_go)

df_afmfit = data.frame(kappa = input.folders[i], mean = mean(as.numeric(data.rank10$afmfit)), sd = sd(as.numeric(data.rank10$afmfit)))
df_afmfit_tot = rbind(df_afmfit_tot, df_afmfit)

df_stack = data.frame(kappa = input.folders[i], mean = mean(as.numeric(data.rank10$stack)), sd = sd(as.numeric(data.rank10$stack)))
df_stack_tot = rbind(df_stack_tot, df_stack)
#print(paste(mean(as.numeric(data.rank10$stack)), sd(as.numeric(data.rank10$stack)), sep = " "))

df_hbond = data.frame(kappa = input.folders[i], mean = mean(as.numeric(data.rank10$hbond)), sd = sd(as.numeric(data.rank10$hbond)))
df_hbond_tot = rbind(df_hbond_tot, df_hbond)

df_repul = data.frame(kappa = input.folders[i], mean = mean(as.numeric(data.rank10$repul)), sd = sd(as.numeric(data.rank10$repul)))
df_repul_tot = rbind(df_repul_tot, df_repul)

df_elect = data.frame(kappa = input.folders[i], mean = mean(as.numeric(data.rank10$elect)), sd = sd(as.numeric(data.rank10$elect)))
df_elect_tot = rbind(df_elect_tot, df_elect)

df_afmcc = data.frame(kappa = input.folders[i], mean = mean(as.numeric(data.rank10$afmcc)), sd = sd(as.numeric(data.rank10$afmcc)))
df_afmcc_tot = rbind(df_afmcc_tot, df_afmcc)
#print(paste(mean(as.numeric(data.rank10$afmcc)), sd(as.numeric(data.rank10$afmcc)), sep = " "))

df_EleRepul = data.frame(kappa = input.folders[i], mean = mean(as.numeric(data.rank10$EleRepul)), sd = sd(as.numeric(data.rank10$EleRepul)))
df_EleRepul_tot = rbind(df_EleRepul_tot, df_EleRepul) 

df_StackHbond = data.frame(kappa = input.folders[i], mean = mean(as.numeric(data.rank10$StackHbond)), sd = sd(as.numeric(data.rank10$StackHbond)))
df_StackHbond_tot = rbind(df_StackHbond_tot, df_StackHbond)


#print(data.rank10)
#print(paste0("selected kappa is: ", input.folders[i]))
}

#print(df_etot_tot)

#print(df_top_all)

df_top_all = df_top_all[order(df_top_all$diffEn.CCscaled, decreasing = FALSE), ]
Escore.file = paste(out.dir, paste0("Escore_ccPow", PCC, ".txt"), sep = "/") 
write.table(format(df_top_all), file = Escore.file, sep = "\t", row.names = FALSE, quote = FALSE)
#print(head(df_top_all))


## write out all stats in a single dataframe and print to file
#print(df_etot_tot$mean)
#print(data.frame(kappa = df_etot_tot$kappa, diff.EtotAfm = df_etot_tot$mean))
df.comb = data.frame(kappa = df_etot_tot$kappa, diff.EtotAfm = df_etot_tot$mean, diff.EtotAfm.sd = df_etot_tot$sd,
                     local = df_local_tot$mean, local.sd = df_local_tot$sd, go = df_go_tot$mean, go.sd = df_go_tot$sd,
                     afmfit = df_afmfit_tot$mean, afmfit.sd = df_afmfit_tot$sd, stack = df_stack_tot$mean, stack.sd = df_stack_tot$sd,
                     hbond = df_hbond_tot$mean, hbond.sd = df_hbond_tot$sd, repul = df_repul_tot$mean, repul.sd = df_repul_tot$sd,
                     elect = df_elect_tot$mean, elect.sd = df_elect_tot$sd, afmcc = df_afmcc_tot$mean, afmcc.sd = df_afmcc_tot$sd, 
                     EleRepul = df_EleRepul_tot$mean, EleRepul.sd = df_EleRepul_tot$sd, StackHbond = df_StackHbond_tot$mean, StackHbond.sd = df_StackHbond_tot$sd)

mean.file = paste(out.dir, "summary_stats.txt", sep = "/")
write.table(format(df.comb), file = mean.file , sep = "\t", row.names = FALSE, quote = FALSE)

## plot summary statistics
plot.file = paste(out.dir, "summary_stats.pdf", sep = "/")
file.lines = readLines(mean.file)

data =  file.lines %>%
  strsplit("\\s+") %>%
  do.call(rbind, .) %>%
  data.frame()

colnames(data) = as.character(data[1,])
data = data[2:nrow(data), ]

# data$diff.EtotAfm = as.numeric(data$diff.EtotAfm)
# data$diff.EtotAfm.sd = as.numeric(data$diff.EtotAfm.sd)
# data$local = as.numeric(data$local)
# data = cbind(data$kappa, as.numeric(data[, 2:ncol(data)]))
## convert into numeric columns
data = cbind(kappa = data$kappa, as.data.frame(sapply(data[, 2:ncol(data)], as.numeric)))
## lock in factor level order
data$kappa = factor(data$kappa, levels = data$kappa)
#class(data)
#class(data$diff.EtotAfm.sd)

pdf(plot.file)

p1 = ggplot(data, aes(x=kappa, y = diff.EtotAfm)) +
  geom_point() + 
  geom_errorbar(aes(ymin= diff.EtotAfm - diff.EtotAfm.sd, 
                    ymax = diff.EtotAfm + diff.EtotAfm.sd ), width = 0.2) +
  theme_classic(base_size = 10)+ ylab("Etot - Afmfit") +
  ggtitle("Etot - Afmfit") +
  theme_bw() +   theme(plot.title = element_text(hjust = 0.5)) +
  xlab("kappa or 3Dscan") + theme(axis.text.x = element_text(angle =90))

p2 = ggplot(data, aes(x=kappa, y = local)) +
  geom_point() + 
  geom_errorbar(aes(ymin= local - local.sd, 
                    ymax = local + local.sd ), width = 0.2) +
  theme_classic(base_size = 10)+ ylab("local") +
  ggtitle("Local Energy") +
  theme_bw() +   theme(plot.title = element_text(hjust = 0.5)) +
  xlab("kappa or 3Dscan") + theme(axis.text.x = element_text(angle =90))

p3 = ggplot(data, aes(x=kappa, y = go)) +
  geom_point() + 
  geom_errorbar(aes(ymin= go - go.sd, 
                    ymax = go + go.sd ), width = 0.2) +
  theme_classic(base_size = 10)+ ylab("go") +
  ggtitle("Go Energy") +
  theme_bw() +   theme(plot.title = element_text(hjust = 0.5)) +
  xlab("kappa or 3Dscan") + theme(axis.text.x = element_text(angle =90))

p4 = ggplot(data, aes(x=kappa, y = afmfit)) +
  geom_point() + 
  geom_errorbar(aes(ymin= afmfit - afmfit.sd, 
                    ymax = afmfit + afmfit.sd ), width = 0.2) +
  theme_classic(base_size = 10)+ ylab("afmfit") +
  ggtitle("Afmfit Energy") +
  theme_bw() +   theme(plot.title = element_text(hjust = 0.5)) +
  xlab("kappa or 3Dscan") + theme(axis.text.x = element_text(angle =90))

p5 = ggplot(data, aes(x=kappa, y = stack)) +
  geom_point() + 
  geom_errorbar(aes(ymin= stack - stack.sd, 
                    ymax = stack + stack.sd ), width = 0.2) +
  theme_classic(base_size = 10)+ ylab("stack") +
  ggtitle("Stack Energy") +
  theme_bw() +   theme(plot.title = element_text(hjust = 0.5)) +
  xlab("kappa or 3Dscan") + theme(axis.text.x = element_text(angle =90))

p6 = ggplot(data, aes(x=kappa, y = hbond)) +
  geom_point() + 
  geom_errorbar(aes(ymin= hbond - hbond.sd, 
                    ymax = hbond + hbond.sd ), width = 0.2) +
  theme_classic(base_size = 10)+ ylab("hbond") +
  ggtitle("hbond Energy") +
  theme_bw() +   theme(plot.title = element_text(hjust = 0.5)) +
  xlab("kappa or 3Dscan") + theme(axis.text.x = element_text(angle =90))

p7 = ggplot(data, aes(x=kappa, y = repul)) +
  geom_point() + 
  geom_errorbar(aes(ymin= repul - repul.sd, 
                    ymax = repul + repul.sd ), width = 0.2) +
  theme_classic(base_size = 10)+ ylab("repul") +
  ggtitle("Repulsive Energy") +
  theme_bw() +   theme(plot.title = element_text(hjust = 0.5)) +
  xlab("kappa or 3Dscan") + theme(axis.text.x = element_text(angle =90))

p8 = ggplot(data, aes(x=kappa, y = elect)) +
  geom_point() + 
  geom_errorbar(aes(ymin= elect - elect.sd, 
                    ymax = elect + elect.sd ), width = 0.2) +
  theme_classic(base_size = 10)+ ylab("elect") +
  ggtitle("Electrostatic Energy") +
  theme_bw() +   theme(plot.title = element_text(hjust = 0.5)) +
  xlab("kappa or 3Dscan") + theme(axis.text.x = element_text(angle =90))

p9 = ggplot(data, aes(x=kappa, y = EleRepul)) +
  geom_point() + 
  geom_errorbar(aes(ymin= EleRepul - EleRepul.sd, 
                    ymax = EleRepul + EleRepul.sd ), width = 0.2) +
  theme_classic(base_size = 10)+ ylab("EleRepul") +
  ggtitle("Electrostatic + Repulsive Energy") +
  theme_bw() +   theme(plot.title = element_text(hjust = 0.5)) +
  xlab("kappa or 3Dscan") + theme(axis.text.x = element_text(angle =90))

p1;p2;p3;p4;p5; p6;p7;p8;p9
dev.off()

## Now carry out the kappa filtering using linear modeling
cat("\n###########################################################\n")
cat("Now carrying out Kappa Filtering\n\n")
df.stats = read.table(mean.file, header = T)
## Here are the list of energy terms to filter against
#min_num = K_cutoff ## minimum number of structures to keep befor exiting the kappa filtering loop

pixel.gp =  unique(sapply(strsplit(df.stats$kappa, "_k"), function(y){return(y[1])}))
kappa.gp = unique(sapply(strsplit(df.stats$kappa, "_k"), function(y){return(y[2])}))
if (!(K_input %in% kappa.gp)){
  print("The kappa you provided is not in the calculation list")
  quit()
}
context.gp = paste(pixel.gp, K_input, sep = "_k")
print(context.gp) 


cat("\n###########################################################\n")
cat("Now carrying out Frame Filtering\n\n")

write(paste("kappa", "frame", "rmsd", "afmcc", sep = "\t"), file = out.pdblist1, append=T)
write(paste("kappa", "frame", "rmsd", "afmcc", sep = "\t"), file = out.pdblist2, append=T)
write(paste("kappa", "frame", "rmsd", "afmcc", sep = "\t"), file = out.pdblist3, append=T)
for (item in df.stats$kappa[1:N_kappa]){
  file.path = paste(item, "output", paste0("Escore_ccPow", PCC, ".txt"), sep = "/")
  df.En = read.delim(file.path, header = TRUE, sep = "\t", dec = ".")
  ## first order by frames
  df.En = df.En[order(df.En$frame, decreasing = FALSE), ]
  #print(dim(df.En))
  ## skip first 4000 frames
  ##df.En = df.En[(F_cutoff+1):nrow(df.En),]
  #print(dim(df.En))
  afmcc = as.numeric(df.En$afmcc)
  df.En$Etot.CCpowN = (as.numeric(df.En$etot))*(afmcc^PCC)
  ## Sorting based on total energy
  df.En.sort = df.En[order(df.En$Etot.CCpowN, decreasing = FALSE), ]
  #print(head(df.En.sort))
  dcdfile = Sys.glob(paste(item, "output","*.dcd", sep = "/"))
  dcd = read.dcd(dcdfile)
  for (frame in df.En.sort$frame[1:N_frame]){
      #print(as.integer(frame))
      afmcc = df.En.sort[df.En.sort$frame == frame, "afmcc"]
      out.pdb = paste(out.pdbdir1,"/",item,"_",frame,".pdb", sep="")
      write.pdb(pdb1, dcd[as.integer(frame), ], file = out.pdb)
      ## write rmsd
      pdb3 = read.pdb(out.pdb)
      rmsd.val = rmsd(pdb2$xyz[, ca.inds$xyz], pdb3$xyz[, ca.inds$xyz], fit = TRUE)
      write(paste(item, frame, rmsd.val, afmcc, sep = "\t"), file = out.pdblist1, append=T)
      }

  stack = as.numeric(df.En$stack)
  hbond = as.numeric(df.En$hbond)
  repul = as.numeric(df.En$repul)
  ## Sorting based on stack and hbond energy
  df.En$EsumStackHbond.CCpowN = (stack+hbond)*(afmcc^PCC)
  df.En.sort = df.En[order(df.En$EsumStackHbond.CCpowN, decreasing = FALSE), ]
  #print(head(df.En.sort))
  for (frame in df.En.sort$frame[1:N_frame]){
      afmcc = df.En.sort[df.En.sort$frame == frame, "afmcc"]
      out.pdb = paste(out.pdbdir2,"/",item,"_",frame,".pdb", sep="")
      write.pdb(pdb1, dcd[as.integer(frame), ], file = out.pdb)
      ## write rmsd
      pdb3 = read.pdb(out.pdb)
      rmsd.val = rmsd(pdb2$xyz[, ca.inds$xyz], pdb3$xyz[, ca.inds$xyz], fit = TRUE)
      write(paste(item, frame, rmsd.val, afmcc, sep = "\t"), file = out.pdblist2, append=T)
      }
  ## Sorting based on stack, hbond, and repulsion energy
  df.En$EsumStackHbondRepul.CCpowN = (stack+hbond+repul)*(afmcc^PCC)
  df.En.sort = df.En[order(df.En$EsumStackHbondRepul.CCpowN, decreasing = FALSE), ]
  #print(head(df.En.sort))
  for (frame in df.En.sort$frame[1:N_frame]){
      afmcc = df.En.sort[df.En.sort$frame == frame, "afmcc"]
      out.pdb = paste(out.pdbdir3,"/",item,"_",frame,".pdb", sep="")
      write.pdb(pdb1, dcd[as.integer(frame), ], file = out.pdb)
      ## write rmsd
      pdb3 = read.pdb(out.pdb)
      rmsd.val = rmsd(pdb2$xyz[, ca.inds$xyz], pdb3$xyz[, ca.inds$xyz], fit = TRUE)
      write(paste(item, frame, rmsd.val, afmcc, sep = "\t"), file = out.pdblist3, append=T)
      }
  Escore.file = paste(item, "output", paste0("FrameFilter_ccPow", PCC, ".txt"), sep = "/")
  write.table(format(df.En), file = Escore.file, sep = "\t", row.names = FALSE, quote = FALSE)
}


