scratch page

x<-read.table("residence.txt",sep=",",stringsAsFactors = FALSE)

vector<-{}

for (i in seq (1,length(x$V1))) {vector<-c(vector,lengths(strsplit(x[i,3]," ")))}

write.table(vector,file="residence_vector.dat") //this takes about two days, save the results for later analysis.

edit residence_vector.dat, change the first row to "index" "residence"

awk 'NR%100==1 {print $0}' residence_vector.dat> residence_vector_short.dat

r<-read.table("residence_vector.dat",header=TRUE,sep=" ",stringsAsFactors = FALSE)

ggplot(r,aes(x=r$residence))+geom_histogram(binwidth=1,colour="black",fill="white")+geom_density(alpha=0.2,fill="#FF6666")

cd ./WATERBOX_amber

R

library(data.table)

library(ggplot2)

filelist<-c("./Pure1/residence_vector.dat","./NaCl_1/residence_vector.dat","./Na20Cl20_1/residence_vector.dat","./Mg20Cl40_1/residence_vector.dat","./Al20Cl60_1/residence_vector.dat")

filelist<-c("./Pure1/residence_vector_short.dat","./NaCl_1/residence_vector_short.dat","./Na20Cl20_1/residence_vector_short.dat","./Mg20Cl40_1/residence_vector_short.dat","./Al20Cl60_1/residence_vector_short.dat")

temp <- lapply(filelist, fread, sep=" ")

data<-melt(temp)

for (i in 1:nrow(data)) {if (data$L1[i] ==1) {data$L1[i]<-"Pure"} else if (data$L1[i] == 2) {data$L1[i]<-"NaCl"} else if (data$L1[i] ==3) {data$L1[i]<-"Na20Cl20"} else if (data$L1[i] == 4) {data$L1[i]<-"Mg20Cl40"} else if (data$L1[i] == 5) {data$L1[i]<-"Al20Cl60"}}

ggplot(data,aes(x=value,colour=L1,fill=L1))+geom_density(alpha=0.1)

hist(vector,col="red",freq=F,xlab="Occupancy of hbonds(in 10000 steps)",main="histogram of occupancy")+axis(side=1,at=seq(0,10000,1000))

As we can see ~60%(0.0012*500)hbonds spans around 500 steps, which is equal to 0.002ps*500=1ps.