R version 3.1.0 (2014-04-10) -- "Spring Dance"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-redhat-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> args=commandArgs(TRUE)
> args
[1] "helix_4_zwit.txt" "6000"            
> commandArgs(trailingOnly=TRUE)
[1] "helix_4_zwit.txt" "6000"            
> args
[1] "helix_4_zwit.txt" "6000"            
> print(commandArgs(TRUE)[1])
[1] "helix_4_zwit.txt"
> filename=as.character(commandArgs(TRUE)[1])
> 
> filename
[1] "helix_4_zwit.txt"
> 
> x=filename
> 
> x
[1] "helix_4_zwit.txt"
> x=2000
> x
[1] 2000
> 
> filename=args[1]
> filename
[1] "helix_4_zwit.txt"
> 
> nbstep=as.numeric(args[2])
> nbstep
[1] 6000
> 
> 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)
[1] -390.1991
> meanEtot=mean(mb.prot1)
> stderr(mb.prot1)
[1] 0.2189926
> stdEtot=stderr(mb.prot1)
> 
> mean(mb.prot1.Et)
[1] -10.16196
> stderr(mb.prot1.Et)
[1] 0.2349225
> 
> 
> 
> #*******************************************************
> #Calcul l'energie de liason et sa moyenne sur les 500 dernieres ps
> 
> ener_start=nbstep/4
> ener_start
[1] 1500
> 
> ener_stop=nbstep/2
> ener_stop
[1] 3000
> 
> 
> WWimm.prot1=mb.prot1-eau.prot1
> #WWimm.prot1
> mean(WWimm.prot1[ener_start:ener_stop])
[1] 8.965136
> meanWWimm=mean(WWimm.prot1[ener_start:ener_stop])
> stderr(WWimm.prot1[ener_start:ener_stop])
[1] 0.03266502
> 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
[1] "grep helix_4_zwit.txt -390.199064693333 0.218992611963469 8.96513566289141 0.032665022815651"
> 
> 
> 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()
null device 
          1 
> 
> 
> 
> 
> 
> 
>