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

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.

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.