## collecting internal energies and solvation terms
mut150 <- NULL
mut150$hept1 <- read.table("20-50ns_freeEnergy_1_1_energies.out_all_energies.dat", sep=",", header=TRUE)
mut150$hept2 <- read.table("20-50ns_freeEnergy_1_2_energies.out_all_energies.dat", sep=",", header=TRUE)
mut150$complex <- read.table("20-50ns_freeEnergy_1_3_energies.out_all_energies.dat", sep=",", header=TRUE)

mut151 <- NULL
mut151$hept1 <- read.table("20-50ns_freeEnergy_2_1_energies.out_all_energies.dat", sep=",", header=TRUE)
mut151$hept2 <- read.table("20-50ns_freeEnergy_2_2_energies.out_all_energies.dat", sep=",", header=TRUE)
mut151$complex <- read.table("20-50ns_freeEnergy_2_3_energies.out_all_energies.dat", sep=",", header=TRUE)

mut161 <- NULL
mut161$hept1 <- read.table("20-50ns_freeEnergy_3_1_energies.out_all_energies.dat", sep=",", header=TRUE)
mut161$hept2 <- read.table("20-50ns_freeEnergy_3_2_energies.out_all_energies.dat", sep=",", header=TRUE)
mut161$complex <- read.table("20-50ns_freeEnergy_3_3_energies.out_all_energies.dat", sep=",", header=TRUE)

mut162 <- NULL
mut162$hept1 <- read.table("20-50ns_freeEnergy_4_1_energies.out_all_energies.dat", sep=",", header=TRUE)
mut162$hept2 <- read.table("20-50ns_freeEnergy_4_2_energies.out_all_energies.dat", sep=",", header=TRUE)
mut162$complex <- read.table("20-50ns_freeEnergy_4_3_energies.out_all_energies.dat", sep=",", header=TRUE)

wt116 <- NULL
wt116$hept1 <- read.table("_20-50ns_freeEnergy_5_1_energies.out_all_energies.dat", sep=",", header=TRUE)
wt116$hept2 <- read.table("_20-50ns_freeEnergy_5_2_energies.out_all_energies.dat", sep=",", header=TRUE)
wt116$complex <- read.table("_20-50ns_freeEnergy_5_3_energies.out_all_energies.dat", sep=",", header=TRUE)

wt129 <- NULL
wt129$hept1 <- read.table("_20-50ns_freeEnergy_6_1_energies.out_all_energies.dat", sep=",", header=TRUE)
wt129$hept2 <- read.table("_20-50ns_freeEnergy_6_2_energies.out_all_energies.dat", sep=",", header=TRUE)
wt129$complex <- read.table("_20-50ns_freeEnergy_6_3_energies.out_all_energies.dat", sep=",", header=TRUE)




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

mut151$hept1$mean <- rowMeans( mut151$hept1 )
mut151$hept2$mean <- rowMeans( mut151$hept2 )
mut151$complex$mean <- rowMeans( mut151$complex )

mut161$hept1$mean <- rowMeans( mut161$hept1 )
mut161$hept2$mean <- rowMeans( mut161$hept2 )
mut161$complex$mean <- rowMeans( mut161$complex )

mut162$hept1$mean <- rowMeans( mut162$hept1 )
mut162$hept2$mean <- rowMeans( mut162$hept2 )
mut162$complex$mean <- rowMeans( mut162$complex )

wt116$hept1$mean <- rowMeans( wt116$hept1 )
wt116$hept2$mean <- rowMeans( wt116$hept2 )
wt116$complex$mean <- rowMeans( wt116$complex )

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$complex$g.gas <- wt116$complex$mean[4] + wt116$complex$mean[5]


## calculating g_solv
wt116$hept1$g.solv <- wt116$hept1$mean[6] + wt116$hept1$mean[9]
wt116$hept2$g.solv <- wt116$hept2$mean[6] + wt116$hept2$mean[9]
wt116$complex$g.solv <- wt116$complex$mean[6] + wt116$complex$mean[9]

## calculating totals
wt116$hept1$total <- wt116$hept1$g.gas + wt116$hept1$g.solv
wt116$hept2$total <- wt116$hept2$g.gas + wt116$hept2$g.solv
wt116$complex$total <- wt116$complex$g.gas + wt116$complex$g.solv




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