NGP_alpha_on_new_PMMA_and_PS.Rmd

---
title: "NGP and alpha on new PMMA and PS"
author: "Zijun Lu"
date: "October 30, 2017"
output:
  html_document:
    toc: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

We are tryng to use a **short time simulation** to get the glass transition temperature compatible with long time simulation. The new simulation on PMMA has 38528 atoms, 64 molecules, 40 monomers per molecule, which leads to 15 atoms per monomer with 2 hydrogen atom on each end of the PMMA molecule. 



```{r,echo=FALSE,include=FALSE}

library(foreach)
library(dplyr)
Data.import=function(Path="~/Dropbox/lammps" , filename = "MSD.colmean.matrix.1.txt", temp=c(320,345,370,395,420,445,470,500,550,600),polymer="PMMA"){
  
  # read in data into a matrix
 MSDcol=foreach(i=1:length(temp),.combine = cbind)%do% {
   # set correct path for the data file
   path=paste(Path,"/", polymer,"/atom",temp[i], sep='')
   path
   setwd(path)
   # read in the data
   MSDcol=read.table(file=filename,header=TRUE,sep=',')
   library(magrittr)
   MSDcol%>%dim
   MSDcol=MSDcol[[1]]
   return(MSDcol)
 }
 
 # add a time variable 
 time=1:dim(MSDcol)[1]
 # change data into data frame 
  MSDcol=cbind(time,MSDcol)%>%as.data.frame()
  
}

```


# Plot the g0,g1 and g2 MSD for PS (25800 atoms) and PMMA (38528 atoms)

To get the more details about the first few steps, I will need to look into the velocity behavior. 
## MSD g0 plot for PS

```{r}

#Run (optimizaed data import function.R) file first 
# temp=c(200,250,300,350,400,450,500,550,600)
# import g0 MSD colmean data for PS 
########################################################
library(dplyr)

MSD.PS.g0=Data.import(filename="MSD.g0.colmean.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS")

MSD.PS.g0=MSD.PS.g0%>%as.data.frame()
MSD.PS.g0=MSD.PS.g0%>%filter(time>1)
# add correct colname 
temp=c(200,250,300,350,400,450,500,550,600)
colnames(MSD.PS.g0)=c("time",paste("T",temp,sep=""))

#regular plot 
###########################################
library(reshape2)
MSD.PS.g0.melt=melt(MSD.PS.g0,id.vars="time")
library(dplyr)

library(ggplot2)
ggplot(data=MSD.PS.g0.melt,aes(x=time,y=value,colour=variable)) +geom_point(size=0.5)+
  ylab("MSD")+
  ggtitle("PS g0 MSD vs time plot")

# log10-log10 plot 
###########################################
library(reshape2)
MSD.PS.g0.melt=melt(MSD.PS.g0,id.vars="time")


library(ggplot2)
ggplot(data=MSD.PS.g0.melt,aes(x=log10(time),y=log10(value),colour=variable)) +geom_point(size=0.5)+
  ylab("log10(MSD)")+
  ggtitle("PS g0 MSD vs time log10-log10 plot")


```

Take a look at the first few timesteps in different temperatures. 
Discard the first data point in plotting and here. 

```{r}
MSD.PS.g0%>%head(10)
```




## MSD g0 plot for PMMA


```{r}

#Run (optimizaed data import function.R) file first 
# temp=c(200,250,300,350,400,450,500,550,600)
# import g0 MSD colmean data for PMMA_big 
########################################################
library(dplyr)

MSD.PMMA_big.g0=Data.import(filename="MSD.g0.colmean.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big")

MSD.PMMA_big.g0=MSD.PMMA_big.g0%>%as.data.frame()
MSD.PMMA_big.g0=MSD.PMMA_big.g0%>%filter(time>1)
# add correct colname 
temp=seq(300,600,by=20)
colnames(MSD.PMMA_big.g0)=c("time",paste("T",temp,sep=""))

#regular plot 
###########################################
library(reshape2)
MSD.PMMA_big.g0.melt=melt(MSD.PMMA_big.g0,id.vars="time")
library(dplyr)

library(ggplot2)
ggplot(data=MSD.PMMA_big.g0.melt,aes(x=time,y=value,colour=variable)) +geom_point(size=0.5)+
  ylab("MSD")+
  ggtitle("PMMA_big g0 MSD vs time plot")

# log10-log10 plot 
###########################################
library(reshape2)
MSD.PMMA_big.g0.melt=melt(MSD.PMMA_big.g0,id.vars="time")


library(ggplot2)
ggplot(data=MSD.PMMA_big.g0.melt,aes(x=log10(time),y=log10(value),colour=variable)) +geom_point(size=0.5)+
  ylab("log10(MSD)")+
  ggtitle("PMMA_big g0 MSD vs time log10-log10 plot")


```

Discard the first data point in plotting and here. 

```{r}
MSD.PMMA_big.g0%>%head()
```




## MSD g2 plot for PS

```{r}

#Run (optimizaed data import function.R) file first 
# import g2 MSD colmean data for PS 
########################################################
library(dplyr)

MSD.PS.g2=Data.import(filename="MSD.g2.colmean.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS")

MSD.PS.g2=MSD.PS.g2%>%as.data.frame()
MSD.PS.g2=MSD.PS.g2%>%filter(time>1)
# add correct colname 
temp=c(200,250,300,350,400,450,500,550,600)
colnames(MSD.PS.g2)=c("time",paste("T",temp,sep=""))

#regular plot 
###########################################
library(reshape2)
MSD.PS.g2.melt=melt(MSD.PS.g2,id.vars="time")
library(dplyr)

library(ggplot2)
ggplot(data=MSD.PS.g2.melt,aes(x=time,y=value,colour=variable)) +geom_point(size=0.5)+
  ylab("MSD")+
  ggtitle("PS g2 MSD vs time plot")

# log10-log10 plot 
###########################################
library(reshape2)
MSD.PS.g2.melt=melt(MSD.PS.g2,id.vars="time")


library(ggplot2)
ggplot(data=MSD.PS.g2.melt,aes(x=log10(time),y=log10(value),colour=variable)) +geom_point(size=0.5)+
  ylab("log10(MSD)")+
  ggtitle("PS g2 MSD vs time log10-log10 plot")


```


```{r}
MSD.PS.g2%>%head(10)
```


## MSD g2 plot for PMMA_big

```{r}

#Run (optimizaed data import function.R) file first 
# import g2 MSD colmean data for PMMA_big 
########################################################
library(dplyr)

MSD.PMMA_big.g2=Data.import(filename="MSD.g2.colmean.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big")

MSD.PMMA_big.g2=MSD.PMMA_big.g2%>%as.data.frame()
MSD.PMMA_big.g2=MSD.PMMA_big.g2%>%filter(time>1)
# add correct colname 
temp=seq(300,600,by=20)
colnames(MSD.PMMA_big.g2)=c("time",paste("T",temp,sep=""))

#regular plot 
###########################################
library(reshape2)
MSD.PMMA_big.g2.melt=melt(MSD.PMMA_big.g2,id.vars="time")
library(dplyr)

library(ggplot2)
ggplot(data=MSD.PMMA_big.g2.melt,aes(x=time,y=value,colour=variable)) +geom_point(size=0.5)+
  ylab("MSD")+
  ggtitle("PMMA_big g2 MSD vs time plot")

# log10-log10 plot 
###########################################
library(reshape2)
MSD.PMMA_big.g2.melt=melt(MSD.PMMA_big.g2,id.vars="time")


library(ggplot2)
ggplot(data=MSD.PMMA_big.g2.melt,aes(x=log10(time),y=log10(value),colour=variable)) +geom_point(size=0.5)+
  ylab("log10(MSD)")+
  ggtitle("PMMA_big g2 MSD vs time log10-log10 plot")


```

```{r}
MSD.PMMA_big.g2%>%head(10)
```



## MSD g1 plot for PMMA_big

```{r}

#Run (optimizaed data import function.R) file first 

# import g1 MSD colmean data for PMMA_big 
########################################################
library(dplyr)

MSD.PMMA_big.g1=Data.import(filename="MSD.g1.colmean.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big")

MSD.PMMA_big.g1=MSD.PMMA_big.g1%>%as.data.frame()
MSD.PMMA_big.g1=MSD.PMMA_big.g1%>%filter(time>1)
# add correct colname 
temp=seq(300,600,by=20)
colnames(MSD.PMMA_big.g1)=c("time",paste("T",temp,sep=""))

#regular plot 
###########################################
library(reshape2)
MSD.PMMA_big.g1.melt=melt(MSD.PMMA_big.g1,id.vars="time")
library(dplyr)

library(ggplot2)
ggplot(data=MSD.PMMA_big.g1.melt,aes(x=time,y=value,colour=variable)) +geom_point(size=0.5)+
  ylab("MSD")+
  ggtitle("PMMA_big g1 MSD vs time plot")

# log10-log10 plot 
###########################################
library(reshape2)
MSD.PMMA_big.g1.melt=melt(MSD.PMMA_big.g1,id.vars="time")


library(ggplot2)
ggplot(data=MSD.PMMA_big.g1.melt,aes(x=log10(time),y=log10(value),colour=variable)) +geom_point(size=0.5)+
  ylab("log10(MSD)")+
  ggtitle("PMMA_big g1 MSD vs time log10-log10 plot")


```


```{r}
MSD.PMMA_big.g1%>%head(10)
```


## MSD g1 plot for PS

```{r}

#Run (optimizaed data import function.R) file first 

# import g1 MSD colmean data for PS 
########################################################
library(dplyr)

MSD.PS.g1=Data.import(filename="MSD.g1.colmean.1.txt", temp=seq(200,600,by=50),polymer="PS")

MSD.PS.g1=MSD.PS.g1%>%as.data.frame()
MSD.PS.g1=MSD.PS.g1%>%filter(time>1)
# add correct colname 
temp=seq(200,600,by=50)
colnames(MSD.PS.g1)=c("time",paste("T",temp,sep=""))

#regular plot 
###########################################
library(reshape2)
MSD.PS.g1.melt=melt(MSD.PS.g1,id.vars="time")
library(dplyr)

library(ggplot2)
ggplot(data=MSD.PS.g1.melt,aes(x=time,y=value,colour=variable)) +geom_point(size=0.5)+
  ylab("MSD")+
  ggtitle("PS g1 MSD vs time plot")

# log10-log10 plot 
###########################################
library(reshape2)
MSD.PS.g1.melt=melt(MSD.PS.g1,id.vars="time")


library(ggplot2)
ggplot(data=MSD.PS.g1.melt,aes(x=log10(time),y=log10(value),colour=variable)) +geom_point(size=0.5)+
  ylab("log10(MSD)")+
  ggtitle("PS g1 MSD vs time log10-log10 plot")


```


```{r}
MSD.PS.g1%>%head(10)
```


# Curvature and NGP methods on PS and PMMA_big


Define a function to read the results, calculate curvature term, R-squared term and non-gaussian parameter. 

In the calculation, I change the time units to ps by multiply it by 0.01 and change the length unit to nm by dividing it by 100.

```{r}

# all even terms and MSD increase linearly with time.
###############################################################################################################
#Run ( data_import.R) file first 




# write a function to plot the alpha term for different polymers and different MSD calculation method. 
#Run (optimizaed data import function.R) file first 


alpha.even.ngp.plot=function(filename1="MSD.g0.colmean.1.txt",filename2="NGP.g0.1.txt", temp=seq(200,600,by=50),polymer="PS",method="linear"){
  
  
  # import colmean data and do linear regression for alpha term 
  ##############################################################
  MSD.PS.g0=Data.import(filename=filename1, temp=temp,polymer=polymer)
  
  MSD.PS.g0=MSD.PS.g0%>%as.data.frame()
  
  
  
  # Regression method 1 : using all even order and no intercept
  ################################################################
  
  # change the time units to ps
  time=MSD.PS.g0[,1]*0.01
  # change  the length unit to nm. 
  MSD.PS.g0=MSD.PS.g0/100
  library(foreach)
  
  
  ############################################################################
  
  if(method=="linear"){
    
    #MSD is increasing linearly with time. 
  alpha=foreach(i=2:(dim(MSD.PS.g0)[2]),.combine = rbind)%do%{
    (lm(MSD.PS.g0[,i]~-1+I(time^2)+I(time^4)+I(time^6)+I(time^8)+I(time^10)+I(time^12)+I(time^14)+I(time^16)+I(time^18)+I(time^20))%>%confint(level=0.95))[2,]
  }
  Rsquared=foreach(i=2:(dim(MSD.PS.g0)[2]),.combine = rbind)%do%{
     c(  ( lm(MSD.PS.g0[,i]~-1+I(time^2)+I(time^4)+I(time^6)+I(time^8)+I(time^10)+I(time^12)+I(time^14)+I(time^16)+I(time^18)+I(time^20)) )%>%summary%>%.$adj.rsquared,
          ( lm(MSD.PS.g0[,i]~-1+I(time^2)+I(time^4)+I(time^6)+I(time^8)+I(time^10)+I(time^12)+I(time^14)+I(time^16)+I(time^18)+I(time^20)) )%>%summary%>%.$r.squared
        )
}
    
  }else if(method=="frac"){
   
  
  #MSD is increasing linearly with fractional power of time. 
  alpha=foreach(i=2:(dim(MSD.PS.g0)[2]),.combine = rbind)%do%{
     (lm(MSD.PS.g0[,i]~-1+I(time^(3/2))+I(time^(7/2))+I(time^(11/2))+I(time^(15/2))+I(time^(19/2)))%>%confint(level=0.95))[2,]
  }
  
  Rsquared=foreach(i=2:(dim(MSD.PS.g0)[2]),.combine = rbind)%do%{
     c(  
       ( 
       lm(MSD.PS.g0[,i]~-1+I(time^(3/2))+I(time^(7/2))+I(time^(11/2))+I(time^(15/2))+I(time^(19/2)) ) 
          )%>%summary%>%.$adj.rsquared,
          ( 
            lm(MSD.PS.g0[,i]~-1+I(time^(3/2))+I(time^(7/2))+I(time^(11/2))+I(time^(15/2))+I(time^(19/2)) ) 
               )  %>%summary%>%.$r.squared
        )
  }
  }else if(method=="solomon"){
    
      #MSD is increasing linearly with time. 
  alpha=foreach(i=2:(dim(MSD.PS.g0)[2]),.combine = rbind)%do%{
    (lm(MSD.PS.g0[,i]/time^2~I(time^2)+I(time^4)+I(time^6)+I(time^8)+I(time^10)+I(time^12)+I(time^14)+I(time^16)+I(time^18)+I(time^20))%>%confint(level=0.95))[2,]
  }
  Rsquared=foreach(i=2:(dim(MSD.PS.g0)[2]),.combine = rbind)%do%{
     c(  ( lm(MSD.PS.g0[,i]/time^2~I(time^2)+I(time^4)+I(time^6)+I(time^8)+I(time^10)+I(time^12)+I(time^14)+I(time^16)+I(time^18)+I(time^20)) )%>%summary%>%.$adj.rsquared,
          ( lm(MSD.PS.g0[,i]/time^2~I(time^2)+I(time^4)+I(time^6)+I(time^8)+I(time^10)+I(time^12)+I(time^14)+I(time^16)+I(time^18)+I(time^20)) )%>%summary%>%.$r.squared
        )
}
    
  }else{
    print("Enter a valid method like linear or frac.")
  }
  
  
  
  ###############################################################
  
  #first transform the matrix then turn it into a vector 
  alpha=alpha%>%t%>%as.vector
  
 
  # combine alpha and temp into a data frame 
  temp=rep(temp,each=2)
  alpha=data.frame(alpha,temp,Rsquared,row.names = NULL)
  
  
  
  # calculate the average between two alpha values 
  alpha_average=alpha%>%group_by(temp) %>% summarize(alpha_average=mean(alpha))
  alpha_average=rbind(alpha_average,alpha_average)%>%arrange(temp)
  
  alpha=cbind(alpha,alpha_average["alpha_average"])
  
  
  # unscaled plot for two 95% confidence interval and average of the end points of the range 
  library(ggplot2)
  
  p1=ggplot(data=alpha, aes(y=alpha/temp,x=temp))+geom_point()+geom_line(aes(x=temp,y=alpha_average/temp),color="red")+
    #geom_smooth(method= "loess")+
    xlab("Temperature")+ylab("Curvature of MSD/Temperature ")
  p3=ggplot(data=alpha, aes(y=Rsquared,x=temp))+geom_point(color="orange")+
    xlab("Temperature")+ylab("R^2 and adj R^2")
  
  
  ####################################################################################################################
  
  #NGP 
  
  NGP=Data.import(filename=filename2, temp=temp,polymer=polymer)
  
  
  temp=temp
  # add column names
  colnames(NGP)=c("time",paste("T",temp,sep=""))
  # delete NA values
  NGP=NGP[complete.cases(NGP),]
  
  NGP=NGP%>%colMeans()
  
  NGP=NGP[-1]
  
  
  
  
  NGP=data.frame(NGP,temp)
  library(ggplot2)
  p2=ggplot(data=NGP, aes(y=NGP,x=temp))+geom_line()+
    xlab("Temperature")+ylab("NGP value")
    
  
  
  ## Plot NGP and alpha term together
  
  
  library(grid)
  grid.newpage()
  grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(p3), size = "first"))
  
  
  #return the plot. 
  
  
  
}  



```

Use two different methods to calculate the alpha term (curvature of MSD), the linear methods is Solomon's method, the frac method is using t^0.5 instead of t in the long term trend. 

On a side note, Solomon's method that divides the MSD by t^2 before doing regression could lead to a very small number causing calculation error.


## PS

### PS g0 alpha term all even orders (t^2 to t^20)

```{r}
alpha.even.ngp.plot(filename1="MSD.g0.colmean.1.txt",filename2="NGP.g0.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS",method="linear")

```
```{r}
alpha.even.ngp.plot(filename1="MSD.g0.colmean.1.txt",filename2="NGP.g0.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS",method="solomon")

```

### PS g0 alpha term fractional orders (t^(3/2) to t^(19/2)) 

```{r}
alpha.even.ngp.plot(filename1="MSD.g0.colmean.1.txt",filename2="NGP.g0.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS",method="frac")

```


### PS g2 alpha term all even orders (t^2 to t^20)

```{r}
alpha.even.ngp.plot(filename1="MSD.g2.colmean.1.txt",filename2="NGP.g2.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS",method="linear")

```

```{r}
alpha.even.ngp.plot(filename1="MSD.g2.colmean.1.txt",filename2="NGP.g2.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS",method="solomon")

```

### PS g2 alpha term fractional orders (t^(3/2) to t^(19/2)) 

```{r}
alpha.even.ngp.plot(filename1="MSD.g2.colmean.1.txt",filename2="NGP.g2.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS",method="frac")

```

## PS

### PS g1 alpha term all even orders (t^2 to t^20)

```{r}
alpha.even.ngp.plot(filename1="MSD.g1.colmean.1.txt",filename2="NGP.g1.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS",method="linear")

```

```{r}
alpha.even.ngp.plot(filename1="MSD.g1.colmean.1.txt",filename2="NGP.g1.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS",method="solomon")

```

### PS g1 alpha term fractional orders (t^(3/2) to t^(19/2)) 

```{r}
alpha.even.ngp.plot(filename1="MSD.g1.colmean.1.txt",filename2="NGP.g1.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS",method="frac")

```

It seems that NGP method is more appropriate for smaller molecules while curvature method is more appropriate for large molecules. 



## PMMA_big 

### PMMA_big g0 alpha term even orders (t^2 to t^20)

```{r}
alpha.even.ngp.plot(filename1="MSD.g0.colmean.1.txt",filename2="NGP.g0.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big",method="linear")

```
```{r}
alpha.even.ngp.plot(filename1="MSD.g0.colmean.1.txt",filename2="NGP.g0.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big",method="solomon")

```

### PMMA_big g0 alpha term fractional orders (t^(3/2) to t^(19/2))

```{r}
alpha.even.ngp.plot(filename1="MSD.g0.colmean.1.txt",filename2="NGP.g0.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big",method="frac")

```

### PMMA_big g2 alpha term even orders (t^2 to t^20)

```{r}
alpha.even.ngp.plot(filename1="MSD.g2.colmean.1.txt",filename2="NGP.g2.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big",method="linear")

```

```{r}
alpha.even.ngp.plot(filename1="MSD.g2.colmean.1.txt",filename2="NGP.g2.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big",method="solomon")

```
### PMMA_big g2 alpha term fractional orders (t^(3/2) to t^(19/2))

```{r}
alpha.even.ngp.plot(filename1="MSD.g2.colmean.1.txt",filename2="NGP.g2.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big",method="frac")

```


### PMMA_big g1 alpha term even orders (t^2 to t^20)

```{r}
alpha.even.ngp.plot(filename1="MSD.g1.colmean.1.txt",filename2="NGP.g1.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big",method="linear")

```
```{r}
alpha.even.ngp.plot(filename1="MSD.g1.colmean.1.txt",filename2="NGP.g1.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big",method="solomon")

```

### PMMA_big g1 alpha term fractional orders (t^(3/2) to t^(19/2))

```{r}
alpha.even.ngp.plot(filename1="MSD.g1.colmean.1.txt",filename2="NGP.g1.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big",method="frac")

```

The result on PMMA_big confirms my assumption on the use of NGP on atom averaged result and use of curvature methods on molecule averaged result.


# Conclusion


## PS

Use the curvature method on g0 MSD and NGP on g0 MSD. Glass transition Temperature: 400K~450K

Tg value is about 100 °C (373K)

```{r}
alpha.even.ngp.plot(filename1="MSD.g0.colmean.1.txt",filename2="NGP.g0.1.txt", temp=c(200,250,300,350,400,450,500,550,600),polymer="PS",method="frac")

```


##PMMA_big

Use the curvature method on g0 MSD and NGP on g0 MSD. Glass transition Temperature: 440K~460K

Tg values of commercial grades of PMMA range from 85 to 165 °C (358K to 438 °F) from wiki (https://en.wikipedia.org/wiki/Poly(methyl_methacrylate))
```{r}
alpha.even.ngp.plot(filename1="MSD.g0.colmean.1.txt",filename2="NGP.g0.1.txt", temp=seq(300,600,by=20),polymer="PMMA_big",method="frac")
```