## collecting internal energies and solvation terms
wt129 <- NULL
wt129$hept1 <- read.table("_SHORT_20-50ns_freeEnergy_5_1_energies.out_all_energies.dat", sep=",", header=TRUE)
wt129$hept2 <- read.table("_SHORT_20-50ns_freeEnergy_5_2_energies.out_all_energies.dat", sep=",", header=TRUE)
wt129$complex <- read.table("_SHORT_20-50ns_freeEnergy_5_3_energies.out_all_energies.dat", sep=",", header=TRUE)



# fetching H, extracting averages over trajectory
wt129$hept1$mean <- rowMeans( wt129$hept1 )
wt129$hept2$mean <- rowMeans( wt129$hept2 )
wt129$complex$mean <- rowMeans( wt129$complex )


## calculating g_gas
wt116$hept1$g.gas <- wt116$hept1$mean[4] + wt116$hept1$mean[5]
wt116$hept2$g.gas <- wt116$hept2$mean[4] + wt116$hept2$mean[5]
wt116$comlex$g.gas <- wt116$complex$mean[4] + wt116$complex$mean[5]


wt116$hept1$total <- wt116$hept1$g.gas + wt116$hept1$mean[11]
wt116$hept2$total <- wt116$hept2$g.gas + wt116$hept2$mean[11]
wt116$comlex$total <- wt116$comlex$g.gas + wt116$complex$mean[11]


wt129$hept1$g.gas <- wt129$hept1$mean[4] + wt129$hept1$mean[5]
wt129$hept2$g.gas <- wt129$hept2$mean[4] + wt129$hept2$mean[5]
wt129$comlex$g.gas <- wt129$complex$mean[4] + wt129$complex$mean[5]


wt129$hept1$total <- wt129$hept1$g.gas + wt129$hept1$mean[11]
wt129$hept2$total <- wt129$hept2$g.gas + wt129$hept2$mean[11]
wt129$comlex$total <- wt129$comlex$g.gas + wt129$complex$mean[11]



## calculating delta G
wt116$deltaG <- wt116$complex$total - wt116$hept1$total - wt116$hept2$total
wt129$deltaG <- wt129$complex$total - wt129$hept1$total - wt129$hept2$total













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",]))