args=commandArgs(TRUE)
args
commandArgs(trailingOnly=TRUE)
args
print(commandArgs(TRUE)[1])
filename=as.character(commandArgs(TRUE)[1])

filename

x=filename

x
x=2000
x

filename=args[1]
filename

nbstep=as.numeric(args[2])
nbstep

prot1=read.table(filename)
#prot1=read.table("calc_polyala23_2.txt")


#**********************************************
#Separate les energies pr3-mb et pr3-eau

mb.prot1=prot1[,1][seq(1,nbstep-1,2)]
eau.prot1=prot1[,1][seq(2,nbstep,2)]
mb.prot1.Et=prot1[,2][seq(1,nbstep-1,2)]
eau.prot1.Et=prot1[,2][seq(2,nbstep,2)]



#***********************************************
#Calcul la moyenne et l'erreur standard de l'energie effective et de l'energie totale

stderr <- function(x) sqrt(var(x)/length(x))

mean(mb.prot1)
meanEtot=mean(mb.prot1)
stderr(mb.prot1)
stdEtot=stderr(mb.prot1)

mean(mb.prot1.Et)
stderr(mb.prot1.Et)



#*******************************************************
#Calcul l'energie de liason et sa moyenne sur les 500 dernieres ps

ener_start=nbstep/4
ener_start

ener_stop=nbstep/2
ener_stop


WWimm.prot1=mb.prot1-eau.prot1
#WWimm.prot1
mean(WWimm.prot1[ener_start:ener_stop])
meanWWimm=mean(WWimm.prot1[ener_start:ener_stop])
stderr(WWimm.prot1[ener_start:ener_stop])
stdWWimm=stderr(WWimm.prot1[ener_start:ener_stop])

out1=paste(filename,"ener.png", sep="_")
#out2=paste(filename,"binding_ener.ps", sep="_")


ener_all=paste("grep",filename, meanEtot,stdEtot,meanWWimm,stdWWimm)
#write(ener_a)
ener_all


jpeg(out1)

par(mfrow=c(2,1))

#postscript(out1)
plot(mb.prot1, xlab="Time(ps)",ylab="Energy",type="l", cex.lab=2, cex.main=2)
#dev.off()


#postscript(out2)
plot(WWimm.prot1, main="Binding energy", xlab="Time(ps)",ylab="WW IMM1-GC",type="l", cex.lab=2, cex.main=2)
dev.off()