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.
We want to read all data files and store them into a tidy tibble for ease of plotting.
tidyverse
and ggplot2
packagesGrt_files
Grt
that will store all
datafor
loop:
r.A
column containing the r values%>%
)
r.A
,
t.ps
and r2Grt
epskt
column containing the \(\varepsilon/k_BT\) valueGrt
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
## # … with 2,815,397 more rows
We will now plot this as a colored heatmap using
geom_raster
from ggplot2
:
Grt
only for times below 100 ps and
\(\varepsilon/k_BT=0.5\), so
recursively filter the data accordingly, then provide this to
ggplot
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.