library(plotrix)
library(Hmisc)



a<-read.table("20ns_apo_Total.dat")
b<-read.table("40ns_apo_Total.dat")
c<-read.table("60ns_apo_Total.dat")
d<-read.table("80ns_apo_Total.dat")
e<-read.table("100ns_apo_Total.dat")


f<-read.table("20ns_holo_Total.dat")
g<-read.table("40ns_holo_Total.dat")
h<-read.table("60ns_holo_Total.dat")
i<-read.table("80ns_holo_Total.dat")
j<-read.table("100ns_holo_Total.dat")

j=rbind(j, j[5,])


apo <- c()
apo = rbind(apo, a[,"V4"])
apo = rbind(apo, b[,"V4"])
apo = rbind(apo, c[,"V4"])
apo = rbind(apo, d[,"V4"])
apo = rbind(apo, e[,"V4"])
apo = apo*300/1000
#apo = cbind(apo, rowMeans(apo))
apo.means <- rowMeans(apo)


holo <- c()
holo = rbind(holo, f[,"V4"])
holo = rbind(holo, g[,"V4"])
holo = rbind(holo, h[,"V4"])
holo = rbind(holo, i[,"V4"])
holo = rbind(holo, j[,"V4"])
holo = holo*300/1000
#holo = cbind(holo, rowMeans(holo))
holo.means <- rowMeans(holo)




x<-c(20,40,60,80,100)

pdf("entropy_vs_simlength.pdf", height=5)
plot(x, apo.means, type="o", col="black", ylab="Entropy (kcal/mol)", xlab="Simulation length (ns)", main="Entropy vs simulation length (Subunit)", ylim=c(3150, 3375))
lines(x, holo.means, type="o", col="red" )


stde <- std.error(t(apo))
errbar( x, apo.means, apo.means+stde, apo.means-stde, add=T, col="black")

stde <- std.error(t(holo))
errbar( x, holo.means, holo.means+stde, holo.means-stde, add=T, col="red", lty=2)

legend("bottomright", legend=c("Unbound", "ATP Bound"), col=c("black", "red"), lty=1)
#text(80,2870,"Based on 4 apo and 4 holo independent simulations")
#text(78,2860,"No significant change in entropy as calculated from MD") 


dev.off()