# plot the MSD of a simulation and compute the diffusion coefficient
# The MSD data of all simulations must be in file all_msd.dat
#
# use from within R: > source("msd.R")
# use from the bash prompt: $ Rscript msd.R
y<-read.table("all_msd.dat")
n<-length(y[,1])
x<-c(1:n)
dt<-10

diffusion <- function(slope,dt) { return(0.192946*slope/dt) }

# compute the slope of a function in x[fitrange]
slope<-function(y){z<-line(x[fitrange],y[fitrange]); return(z[[2]][[2]])}

fitrange=20000:n
z<-line(x[fitrange],y[fitrange,1])
# compute all diffusion coefficients
s<-apply(y,2,slope)
D<-diffusion(s,dt)
cat("Diffusion coefficients (cm^2/s)\n")
print(D)
# save diffusion coeffs in Dcoeffs.dat
write(D,file="Dcoeffs.dat",ncolumns=1)
meanD<-mean(D)
cat("mean(D) = ",meanD,"\n")
sdD<-sd(D)
cat("sd(D) = ",sdD,"\n")

# compute average MSD
a<-rowMeans(y)
# fit linear function in fitrange
z<-line(x[fitrange],a[fitrange])
# compute diffusion coefficient
savg<-z[[2]][[2]]
Davg<-diffusion(savg,dt)
cat("Davg = ",Davg," cm^2/s\n")

# plot all msd
#plot(x,y[,1],type="l",xlab="time (step)",ylab="MSD (Bohr^2)")
#apply(y,2,lines,type="l")

# plot the average msd
#lines(x,a,type="l",lwd="3",col="red")