com <- read.table("_MMPBSA_complex_gb.mdout_all_energies.dat", sep=",", header=TRUE)
rec <- read.table("_MMPBSA_receptor_gb.mdout_all_energies.dat", sep=",", header=TRUE)
lig <- read.table("_MMPBSA_ligand_gb.mdout_all_energies.dat", sep=",", header=TRUE)





## collecting internal energies and solvation terms
com <- read.table("com.results_all_energies.dat", sep=",", header=TRUE)
rec <- read.table("rec.results.chain_1_all_energies.dat", sep=",", header=TRUE)
lig <- read.table("lig.results.chain_1_all_energies.dat", sep=",", header=TRUE)

com$means <- rowMeans(com[,1:8])
rec$means <- rowMeans(rec[,1:8])
lig$means <- rowMeans(lig[,1:8])

deltag <- com$means - rec$means - lig$means





## fetching H - only the last 600 frames due to entropy calculations
apo.h <- c()
apo.h <- rbind(apo.h, as.numeric(apo1["total",1401:2000]))
apo.h <- rbind(apo.h, as.numeric(apo2["total",1401:2000]))
apo.h <- rbind(apo.h, as.numeric(apo3["total",1401:2000]))
apo.h <- rbind(apo.h, as.numeric(apo4["total",1401:2000]))

holo.h <- c()
holo.h <- rbind(holo.h, as.numeric(holo1["total",1401:2000]))
holo.h <- rbind(holo.h, as.numeric(holo2["total",1401:2000]))
holo.h <- rbind(holo.h, as.numeric(holo3["total",1401:2000]))
holo.h <- rbind(holo.h, as.numeric(holo4["total",1401:2000]))



## calculating G = H - S
apo.g <- apo.h  - apo.s
holo.g <- holo.h - holo.s

apo.g.means <- rowMeans(apo.g)
holo.g.means <- rowMeans(holo.g)
##delta.g <- apo.g.means - holo.g.means

apo.h.means <- rowMeans(apo.h)
holo.h.means <- rowMeans(holo.h)
##delta.h <- apo.h.means - holo.h.means

t.test(apo.h.means, holo.h.means, "greater", var.equal=T)


plot(apo.g.means)



# argghh... no differences



pdf( "energies.pdf", height=8 )


par(mfrow=c(3,1))
plot(matrix(apo1["g_solv",]), type="l")
lines(matrix(holo1["g_solv",]), type="l", col="red")

plot(matrix(apo1["e_int",]), type="l")
lines(matrix(holo1["e_int",]), type="l", col="red")

plot(matrix(apo1["total",]), type="l")
lines(matrix(holo1["total",]), type="l", col="red")



dev.off()


apo.means = c()
apo.means = c(apo.means, mean(as.matrix(do.call(cbind, apo1["total",])[1,])))
apo.means = c(apo.means, mean(as.matrix(do.call(cbind, apo2["total",])[1,])))
apo.means = c(apo.means, mean(as.matrix(do.call(cbind, apo3["total",])[1,])))
apo.means = c(apo.means, mean(as.matrix(do.call(cbind, apo4["total",])[1,])))

holo.means = c()
holo.means = c(holo.means, mean(as.matrix(do.call(cbind, holo1["total",])[1,])))
holo.means = c(holo.means, mean(as.matrix(do.call(cbind, holo2["total",])[1,])))
holo.means = c(holo.means, mean(as.matrix(do.call(cbind, holo3["total",])[1,])))
holo.means = c(holo.means, mean(as.matrix(do.call(cbind, holo4["total",])[1,])))










mean(as.matrix(do.call(cbind, apo1["e_int",])[1,]))
mean(as.matrix(do.call(cbind, apo2["e_int",])[1,]))
mean(as.matrix(do.call(cbind, apo3["e_int",])[1,]))
mean(as.matrix(do.call(cbind, apo4["e_int",])[1,]))


mean(as.matrix(do.call(cbind, holo1["e_int",])[1,]))
mean(as.matrix(do.call(cbind, holo2["e_int",])[1,]))
mean(as.matrix(do.call(cbind, holo3["e_int",])[1,]))
mean(as.matrix(do.call(cbind, holo4["e_int",])[1,]))



mean(as.matrix(do.call(cbind, apo1["g_solv",])[1,]))
mean(as.matrix(do.call(cbind, apo2["g_solv",])[1,]))
mean(as.matrix(do.call(cbind, apo3["g_solv",])[1,]))
mean(as.matrix(do.call(cbind, apo4["g_solv",])[1,]))


mean(as.matrix(do.call(cbind, holo1["g_solv",])[1,]))
mean(as.matrix(do.call(cbind, holo2["g_solv",])[1,]))
mean(as.matrix(do.call(cbind, holo3["g_solv",])[1,]))
mean(as.matrix(do.call(cbind, holo4["g_solv",])[1,]))







apo.h <- c()
apo.h <- rbind(apo.h, as.numeric(apo1["total",1:1500]))
apo.h <- rbind(apo.h, as.numeric(apo2["total",1:1500]))
apo.h <- rbind(apo.h, as.numeric(apo3["total",1:1500]))
apo.h <- rbind(apo.h, as.numeric(apo4["total",1:1500]))

holo.h <- c()
holo.h <- rbind(holo.h, as.numeric(holo1["total",1:1500]))
holo.h <- rbind(holo.h, as.numeric(holo2["total",1:1500]))
holo.h <- rbind(holo.h, as.numeric(holo3["total",1:1500]))
holo.h <- rbind(holo.h, as.numeric(holo4["total",1:1500]))

apo.h.means <- rowMeans(apo.h)
holo.h.means <- rowMeans(holo.h)
##delta.h <- apo.h.means - holo.h.means

t.test(apo.h.means, holo.h.means, "greater", var.equal=T)




par(mfcol=c(4,2))

plot( as.numeric(apo1["e_int",]), as.numeric(apo1["g_solv",]))
plot( as.numeric(apo2["e_int",]), as.numeric(apo2["g_solv",]))
plot( as.numeric(apo3["e_int",]), as.numeric(apo3["g_solv",]))
plot( as.numeric(apo4["e_int",]), as.numeric(apo4["g_solv",]))

plot( as.numeric(holo1["e_int",]), as.numeric(holo1["g_solv",]))
plot( as.numeric(holo2["e_int",]), as.numeric(holo2["g_solv",]))
plot( as.numeric(holo3["e_int",]), as.numeric(holo3["g_solv",]))
plot( as.numeric(holo4["e_int",]), as.numeric(holo4["g_solv",]))