R Exercises - G(r,t)

A bit of context

In this exercise, we will see how to make a color plot of $4\pi r^2\times G(r,t)$ data obtained from Molecular Dynamics simulations. These auto-correlators $G(r,t)$ are obtained from MD trajectories of a fluid encapsulated in a porous medium by computing the probability that each molecule moves by a distance $r$ over a time $t$. The MD trajectories are performed with LAMMPS and the trajectories treated by a home-made C program producing these files.

In these simulations, we varied the depth of the interaction well. It is given in the name of the file such as Grt_xx.dat, where xx = 0.01, 0.1, 0.2, 0.3, 0.5, 0.8, 1 corresponds to $\varepsilon/k_BT$. $\varepsilon/k_BT$ is the depth of the interaction well in units of temperature.

The Grt_xx.dat files contain the $4\pi r^2\times G(r,t)$ data under the form of a matrix, with $r$ increasing with the rows and $t$ with the columns. $r$ goes from 0 to 20 Å and $t$ goes from 0 to 200 ps.

Reading and tidying the data

We want to read all data files and store them into a tidy tibble for ease of plotting.

Load the tidyverse and ggplot2 packages
Find all files starting by “Grt” and store them in a vector Grt_files
Initiate an empty tibble Grt that will store all data
Using a for loop:
- read the files
- add columns names as the times values
- add the r.A column containing the r values
- and recursively (using the pipe |>)
  - pivot the table to make it tidy, with the columns r.A, t.ps and r2Grt
  - make sure that the times are read as numeric values
  - add the epskt column containing the $\varepsilon/k_BT$ value
- store this in Grt

library(tidyverse)
library(ggplot2)
Grt_files <- list.files(pattern = "Grt_", path="Data")
Grt       <- tibble()
for(Grt_file in Grt_files){#Grt_file <- Grt_files[1]
    d <- read_table(file.path("Data",Grt_file), col_names=FALSE)
    names(d) <- as.character(seq(0,200, length=ncol(d)))
    d$r.A    <- seq(0,20, length=nrow(d))
    eps <- gsub(".dat","",gsub("Grt_","",Grt_file))
    d <- d |> 
        pivot_longer(cols=-r.A,
                     names_to="t.ps",
                     names_transform=list(t.ps = as.numeric),
                     values_to="r2Grt") |> 
        mutate(epskt=as.numeric(eps))
    Grt <- bind_rows(Grt, d)
}

Plotting data

We should have a Grt dataset that looks like this:

Grt

# A tibble: 2,815,407 × 4
     r.A  t.ps      r2Grt epskt
   <dbl> <dbl>      <dbl> <dbl>
 1     0   0   77.7        0.01
 2     0   0.1  0.000512   0.01
 3     0   0.2  0.000203   0.01
 4     0   0.3  0.000148   0.01
 5     0   0.4  0.0000979  0.01
 6     0   0.5  0.0000817  0.01
 7     0   0.6  0.0000714  0.01
 8     0   0.7  0.0000683  0.01
 9     0   0.8  0.0000687  0.01
10     0   0.9  0.0000552  0.01
# ℹ 2,815,397 more rows

We will now plot this as a colored heatmap using geom_raster from ggplot2:

We will first plot Grt only for times below 100 ps and $\varepsilon/k_BT=0.5$, so recursively filter the data accordingly, then provide this to ggplot
Use a nicer axis labels
Do you see anything on the plot?
In fact, we need to saturate values above a certain limit in order to see the fine evolution of the $G(r,t)$ auto-correlator. Add a mutation of the tibble attributing the limit value to all points above this limit, and play with the limit value in order to have a nicer plot.
Add a vertical dashed line in x=8.18Å marking the mean pore size for this structure.
Use a nicer color scheme with more colors in order to better see the variations.

lim    <- 0.1
nbreak <- 6
colors <- colorRampPalette(c("white","royalblue","seagreen","orange","red","brown"))(50)
Grt |> 
    filter(t.ps<=100, epskt==0.5) |> 
    mutate(r2Grt=ifelse(r2Grt>lim, lim, r2Grt)) |> 
    ggplot(aes(x=r.A, y=t.ps, fill=r2Grt))+
        geom_raster()+
        geom_vline(aes(xintercept=8.18), lty=2)+
        scale_fill_gradientn(
            colors = colors, 
            name   = '',
            breaks = seq(0, lim, length=nbreak), 
            labels = c(seq(0, lim, length=nbreak)[1:(nbreak-1)], paste(">",lim)) 
        ) +
        labs(x = "r [Å]", y="Time [ps]") +
        scale_x_continuous(expand = c(0,0)) + 
        scale_y_continuous(expand = c(0,0))

Now we want to see the evolution of $G(r,t)$ as a function of $\varepsilon/k_BT$. Using the procedure above, plot $G(r,t)$ for all $\varepsilon/k_BT$ on a grid.

lim    <- 0.1
nbreak <- 6
colors <- colorRampPalette(c("white","royalblue","seagreen","orange","red","brown"))(50)
Grt |> 
    filter(t.ps<=100) |> 
    mutate(r2Grt=ifelse(r2Grt>lim, lim, r2Grt)) |> 
    ggplot(aes(x=r.A, y=t.ps, fill=r2Grt))+
        geom_raster()+
        geom_vline(aes(xintercept=8.18), lty=2)+
        scale_fill_gradientn(
            colors = colors, 
            name   = '',
            breaks = seq(0, lim, length=nbreak), 
            labels = c(seq(0, lim, length=nbreak)[1:(nbreak-1)], paste(">",lim)) 
        ) +
        labs(x = "r [Å]", y="Time [ps]") +
        scale_x_continuous(expand = c(0,0)) + 
        scale_y_continuous(expand = c(0,0)) +
        facet_wrap(~epskt)

--- title : "R Exercises - G(r,t)" format: html: toc : true theme : cosmo page-layout : full linkcolor : "#3d54cb" code-link : true code-tools : true code-copy : true code-block-border-left: true code-block-bg : true self-contained : true params: solution: true execute: warning: false message: false cache: false --- ```{r echo=FALSE, warning=FALSE, message=FALSE} xfun::embed_file("./Archive.zip", text = "Download data files") ``` --- ```{r setup, include=FALSE, warning = FALSE, cache=FALSE} library(plotly) library(ggplot2) library(tidyverse) my_theme <- theme_bw()+ theme(axis.text = element_text(size = 14,colour="black"), text = element_text(size = 14), axis.ticks = element_blank(), plot.title = element_text(size=14), legend.text = element_text(size = 14,colour="black"), legend.title = element_text(size = 14,colour="black"), panel.border = element_rect(colour = "black", fill=NA, size=1), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.key.height=unit(1,"cm"), legend.justification = "center", strip.background = element_blank(), strip.text=element_text(size = 14), strip.text.y = element_text(angle=0), panel.spacing = unit(1, "lines") ) theme_set(my_theme) ``` # A bit of context In this exercise, we will see how to make a color plot of $4\pi r^2\times G(r,t)$ data obtained from Molecular Dynamics simulations. These auto-correlators $G(r,t)$ are obtained from MD trajectories of a fluid encapsulated in a porous medium by computing the probability that each molecule moves by a distance $r$ over a time $t$. The MD trajectories are performed with [LAMMPS](https://lammps.sandia.gov/) and the trajectories treated by a home-made C program producing these files. In these simulations, we varied the depth of the interaction well. It is given in the name of the file such as `Grt_xx.dat`, where xx = 0.01, 0.1, 0.2, 0.3, 0.5, 0.8, 1 corresponds to $\varepsilon/k_BT$. $\varepsilon/k_BT$ is the depth of the interaction well in units of temperature. The `Grt_xx.dat` files contain the $4\pi r^2\times G(r,t)$ data under the form of a matrix, with $r$ increasing with the rows and $t$ with the columns. $r$ goes from 0 to 20 Å and $t$ goes from 0 to 200 ps. # Reading and tidying the data We want to read all data files and store them into a tidy tibble for ease of plotting. - Load the `tidyverse` and `ggplot2` packages - Find all files starting by "Grt" and store them in a vector `Grt_files` - Initiate an empty tibble `Grt` that will store all data - Using a `for` loop: - read the files - add columns names as the times values - add the `r.A` column containing the r values - and recursively (using the pipe `|>`) - pivot the table to make it tidy, with the columns `r.A`, `t.ps` and `r2Grt` - make sure that the times are read as numeric values - add the `epskt` column containing the $\varepsilon/k_BT$ value - store this in `Grt` ```{r include=params$solution} library(tidyverse) library(ggplot2) Grt_files <- list.files(pattern = "Grt_", path="Data") Grt <- tibble() for(Grt_file in Grt_files){#Grt_file <- Grt_files[1] d <- read_table(file.path("Data",Grt_file), col_names=FALSE) names(d) <- as.character(seq(0,200, length=ncol(d))) d$r.A <- seq(0,20, length=nrow(d)) eps <- gsub(".dat","",gsub("Grt_","",Grt_file)) d <- d |> pivot_longer(cols=-r.A, names_to="t.ps", names_transform=list(t.ps = as.numeric), values_to="r2Grt") |> mutate(epskt=as.numeric(eps)) Grt <- bind_rows(Grt, d) } ``` # Plotting data We should have a `Grt` dataset that looks like this: ```{r include=TRUE} Grt ``` We will now plot this as a colored heatmap using `geom_raster` from `ggplot2`: - We will first plot `Grt` only for times below 100 ps and $\varepsilon/k_BT=0.5$, so recursively filter the data accordingly, then provide this to `ggplot` + Use a nicer axis labels + Do you see anything on the plot? + In fact, we need to saturate values above a certain limit in order to see the fine evolution of the $G(r,t)$ auto-correlator. Add a mutation of the tibble attributing the limit value to all points above this limit, and play with the limit value in order to have a nicer plot. - Add a vertical dashed line in x=8.18Å marking the mean pore size for this structure. - Use a nicer color scheme with more colors in order to better see the variations. ```{r include=params$solution} lim <- 0.1 nbreak <- 6 colors <- colorRampPalette(c("white","royalblue","seagreen","orange","red","brown"))(50) Grt |> filter(t.ps<=100, epskt==0.5) |> mutate(r2Grt=ifelse(r2Grt>lim, lim, r2Grt)) |> ggplot(aes(x=r.A, y=t.ps, fill=r2Grt))+ geom_raster()+ geom_vline(aes(xintercept=8.18), lty=2)+ scale_fill_gradientn( colors = colors, name = '', breaks = seq(0, lim, length=nbreak), labels = c(seq(0, lim, length=nbreak)[1:(nbreak-1)], paste(">",lim)) ) + labs(x = "r [Å]", y="Time [ps]") + scale_x_continuous(expand = c(0,0)) + scale_y_continuous(expand = c(0,0)) ``` Now we want to see the evolution of $G(r,t)$ as a function of $\varepsilon/k_BT$. Using the procedure above, plot $G(r,t)$ for all $\varepsilon/k_BT$ on a grid. ```{r include=params$solution} lim <- 0.1 nbreak <- 6 colors <- colorRampPalette(c("white","royalblue","seagreen","orange","red","brown"))(50) Grt |> filter(t.ps<=100) |> mutate(r2Grt=ifelse(r2Grt>lim, lim, r2Grt)) |> ggplot(aes(x=r.A, y=t.ps, fill=r2Grt))+ geom_raster()+ geom_vline(aes(xintercept=8.18), lty=2)+ scale_fill_gradientn( colors = colors, name = '', breaks = seq(0, lim, length=nbreak), labels = c(seq(0, lim, length=nbreak)[1:(nbreak-1)], paste(">",lim)) ) + labs(x = "r [Å]", y="Time [ps]") + scale_x_continuous(expand = c(0,0)) + scale_y_continuous(expand = c(0,0)) + facet_wrap(~epskt) ```