R Exercises - G(r,t)

1 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 files can be downloaded here. 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.

2 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)
}

3 Plotting data

We should have a Grt dataset that looks like this:

Grt

## # A tibble: 2,815,407 x 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
## # … with 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)