# Directionality in glassy systems
###### tags: `PhD projects`

One of the approaches to understand the relation between structure and dynamics in glassy systems has been the use of machine learning. Up and til now we have used rotational invariant quantities to capture both the structure and the dynamics, however it would be interesting to see whether we can understand more about how the structure influences the directions particles travel to and, especially, escape their cage. Doing so would possibly help us understand to what extend the direction particles travel to, matters when one wants to understand glassy dynamics.

## Theory
### Cage breaking
To know whether directionality matters in the cage breaking process, we need to have a measure that tells you whether or not a particle is escaped. For this we use the $D_\text{min}$ defined in Ref [1].

$$
D^2(i, t, \Delta t) = \min_{\epsilon_{\alpha\beta}}\left[\frac{1}{z}\sum_\alpha\sum_j\left((r^\alpha_j(t) - r^\alpha_i(t))- \sum_\beta(\epsilon_{\alpha\beta}+\delta_{\alpha\beta})(r^\beta_j(t- \Delta t)-r_i^\beta(t-\Delta t))\right)^2\right]
$$

where
$$
\epsilon_{\alpha\beta} = \sum_\gamma X_{\alpha\gamma}Y^{-1}_{\beta\gamma}-\delta_{\alpha\beta}
$$

With
$$
X_{\alpha\gamma} = \sum_j[r^\alpha_j(t)-r^\alpha_i(t)][r^\gamma_j(t-\Delta t)- r_i^\gamma(t-\Delta t)]
$$
and
$$
Y_{\beta\gamma}= \sum_j[r^\beta_j(t-\Delta t)- r^\beta_i(t-\Delta t)][r^\gamma_j(t-\Delta t) - r^\gamma_i(t-\Delta t)]
$$
In these equations, the subscripts *i* and *j* indicate different particles. The sum $\sum_j$ sums over all particles that are in a radius of $2.5\sigma_A$ of particle *i*. The total number of neighbours is defined as $z$. The superscripts $\alpha$ and $\beta$ indicate the different coordinates $x,y$ and $z$. Note that in the formula for $\epsilon$ the matrix $Y^{-1}$ is the inverse matrix of $Y$, note that $Y$ (and thus $Y^{-1}$) is a symmetric matrix.

The value of $D_\text{min}^2(i, t, \Delta t)$ computes the difference between the actual displacemt and the affine-displacement (which is the displacement that would occur when the particles where in a region of uniform stress).If $D_\text{min}^2(i, t, \Delta t)$ is larger than a certain treshold value $r_c$ we say that the partice has escaped.

::: warning
Note that in the paper of Andrea Liu (Fragility in Glassy Liquids: A Structural Approach Based on MachineLearning), the defenition is slightly different (a minus sign before the $\epsilon_{ij}$). I however think, that since we look at the difference between the actual displacement and the affine displacement, the definition by [1] is correct.

Moreover in the definition of $D^2(i, t, \Delta t)$ it says that we minimize the value of $D^2(i, t, \Delta t)$ w.r.t $\epsilon_{\alpha\beta}$. In the origininal paper [1], they state that this minimum is found when the $\epsilon_{\alpha\beta}$ is compute as written down above, so that is what we do

In the original paper they do not average over the number of neighbours, however I do. First of all, doing say makes sure that $D^2(i, t, \Delta t)$ does not depend on how many neighbours there are. Moreover, in the paper of Andrea Liu they indeed average over numer of neighbours.
:::

[1] *Dynamics of viscoplastic deformation in amorphous solid, M. L. Falk and J. S. Langer (https://journals.aps.org/pre/pdf/10.1103/PhysRevE.57.7192)*

#### Finding the critical $D_\text{min}$
*Can be found on ODIN in /home/users/5670772/glass/Analysis/code_Dmin/CHECKS/Critical_Dmin.nb. Note that the data that we use for the analysis can be found in /home/users/5670772/glass/Analysis/KA_Dmin/KA_1/run_1001/ runs $t_0$ to $t_{885}$*

In her paper, Andrea Liu defines the critical $D_\text{min}$ as the $D_\text{min}$ for which $1\%$ of the particles escapes in a time $\Delta t$. We can however not use her $D_\text{min}$ because she looks at the harmonic system.

In our case we have $\Delta t/\tau = 15$. To find the corresponding critical $D_\text{min}$ for the KA system, we run 59 simulations with $\Delta t= 15$ and plot the $D_\text{min}$ histogram (blue) together with a fitter distribution (black) together with the two found $D_\text{min}$'s (two vertical lines). Below we will explain how we found this critical $D_\text{min}$.


By computing
$$
\frac{\int_0^{D_\text{min}}f(x)dx}{\int_0^{\infty}f(x)dx}\approx 0.99
$$
we find $D_\text{min} \approx 0.165$. Note however that the blacak

However by computing
$$
\frac{\sum_{n=0}^{n_{Dmin}}N(n)}{\sum_{n=0}^{\infty}N(n)}\approx 0.99
$$
where we sum over the histogram bins and $N(n)$ is the number of particles that as associated with bin $n$. Computing this we find $D_\text{min} \approx 0.28$

Finally we use one other method to find the critical $D_\text{min}$. Instead of looking at the PDF, we look at the CDF (note that this essentially the same, however fitting the function is easier). We make an interpolation of the data of the CDF, and then compute for which $D_\text{min}$ we reach 0.99. Again this turns out to be $D_\text{min}\approx0.28$.

### Escape direction
We do not only want to know whether or not a particles escapes, we also want to know in what direction. For this we use the vector propensity, which is defined as
$$
\langle\Delta r^\alpha_i(t)\rangle = \sum_n r_i^\alpha(t)-r_i^\alpha(0),
$$
where we consider particle $i$ and $\alpha\in{x,y,z}$.

However, since we ideally want to direction in which the particle escapes w.r.t. the cage center and not the initial position, we instead propose the following new definition
$$
\langle\Delta r^\alpha_i(t)\rangle = \sum_n r_i^\alpha(t)-r_i^{\alpha,\text{ cage}},
$$
where we find the cage positions using either the inherent state and the cage state.

The reason that we use an altered form of the propensity (i.e. a configurational ensemble average) to investigate the cage breaking direction, is because cagebreaking is a stochastic process: Considering one run solely will not tell us whether there is a preferred escape direction. Using an ensemble average however has the downside that although for the first escaping particles it is possible to assign an escape direction, after the first cage escape events the difference between different runs will be too big. As a result we can only use propensity at short timescales as a measure for which direction particles escape their cage.

### Problem with defining escape direction
*Found on ODIN /home/users/5670772/glass/Analysis/code_Dmin/CHECKS/KA_Dmin_check/direction_propensity.nb*
We use the already obtained data for the KA system to investigate how the exact cage center influences the cage escape results.
First of all we plot the propensity w.r.t. either the cage state (obtained with MC simulation) or the inherent state (obtained by quenching the system).

From these results we see that the propensity w.r.t the cage state during the caging regime is somewhat lower than the propensity w.r.t. to the inherent state. This is something we expect because the cage state is better in capturing the positions of particles during the caging regime.

This difference between the cage state and the inherent state however, turns out to lead to different results when we plot the vector propensity with respect to the two different states. Below we show on the right side the vector propensity w.r.t. the inehrent state and on the right the vector propensity w.r.t to the cage state (both at $t/\tau=20$, note that vectors are increased by a factor of 5 to make them more visible). Note that both pictures display the same region of the system. As we can see, in the left picture particles all seem to move to the left, while on the right particles seem to move to the right.

<img src="https://i.imgur.com/p7D9bld.png
" alt="drawing" width="300"/> <img src="https://i.imgur.com/BccvgKz.png
" alt="drawing" width="300" />

From this we can see that directionality, up to a certain degree, might be an ill-defined concept in glassy system. For particles that are not located in a very well defined cage (which these particles are, judging on the very large difference between cage state and inherent state), whether or not a particle escapes to a certain direction might be highly dependend on where we choose the reference position of the particle.

## CHECKS
### Check DMIN calculation
*Found on ODIN in the folder /home/users/5670772/glass/Analysis/code_Dmin/CHECKS/KA_Dmin_check the analysis below can be found in the notebook Dmin.nb*

We start with an initial equilibrated configuration. Then we let the system evolve 50 times with different initial velocities for $t/\tau= 15$ (in lj units), which corresponds to the peak time of the caging regime. For each of these runs we calculate $D_\text{min}(i, t =15, \Delta = 15)$ of all the particles. We consider a Kob-Anderson mixture of 2000 particles at $k_BT/\epsilon =0.44$ and $\rho/\sigma^3 = 1.203$.

First of all we look at the ensemble average $\langle D_\text{min}(i, t =15, \Delta = 15)\rangle$ (i.e. averaged over all the different runs), which is plotted in the figure below for the $A$ particles.


Then we compute the correlation between the propensities of these particles at $t/\tau = 15$ and $\langle D_\text{min}(i, t =15, \Delta = 15)\rangle$, these turn out to be for the $A$ particles

$\rho(\Delta r(t/\tau = 15), \langle D_\text{min}(i, t =15, \Delta = 15)\rangle ) = 0.90$

And for the $B$ particles

$\rho(\Delta r(t/\tau = 15), \langle D_\text{min}(i, t =15, \Delta = 15)\rangle ) = 0.94$.

Since this there is such a high correlation between the propensity and $D_\text{min}$ (which we expect to be highly coupled), we conclude that the $D_\text{min}$ calculation is correct.

## Research proposal
Given the preliminary results, we propose the following:

We compute the vector propensity up to $t/\tau = 15$ (using 50 different runs starting from the same initial position).

Using the definition of $D_{min}$ we compute per run which particles have escaped. This gives us an escape probility per particle. For the particles that have escaped more than a certain treshold value (f.e. more than 25 times out of the 50 runs), we first qualitatively look at the escape directions.

With this qualitative analysis we want to answer the following questions:
1. Are there preferred escape directions for the early escaping particles (this question can also be answered by comparing the propensity and the vector propensity with each other. If particles always escape in the same direction the norm of the vector propensity should be close to the normal propensity. If particles escape in random directions, the vetor propensity will be small, while the normal propensity is high).

:::info
Below we show to which directions the three particles move that have escaped in run 1001 (:warning: I lost the notebook with which I made these picture, so do not base any real analysis on them). These particles have escaped in more than half of the runs. The dots that are show are the escaped directions with respect to of the escaped Pruns with respect to the IS (yellow) and the CS (red).


If we plot the points w.r.t. the particle positions we obtain the following plots (again yellow is IS, red is CS and black is the initial state).
:::

2. Does it matter w.r.t. to which reference point you define the escape direction (i.e. initial state, cage state and inherent state).

**Follow up steps**

If this analysis leads to the conclusion that there is indeed a good defined escape direction, the following step would be trying to predict this escape direction bases on parameters.

::: warning
This leads to the old question what kind of parameters are needed when one needs to make directional predictions. On the one hand we might use the normal BOPS and then use some exponential function in certain directions. We might also turn to more complex machine learning methods like the directional GNN (I read through the paper SE(3)-equivariant Graph Neural Networks for Learning Glassy Liquids Representations of Francesco Saverio Pezzicoli, Guillaume Charpiat, François P. Landes). Or design completely new parameters (maybe using SOAP parameters).
:::

Moreover, in order to train models and make predictions, we need significantly more data. The idea is to obtain this data by, after $t/\tau = 15$, taking one of the final configurations as a new initial snapshot. This way we obtain more escape data, and moreover we can track the dynamical evolution of the system.

## Results escaping particles
*Results can be found on ODIN /home/users/5670772/glass/Analysis/KA_Dmin/KA_1/run_1001/. The different t directories contain the data for consecutive runs.*

*Codes that generate the random vector data and the data itself can be found on Odin in the folder /home/users/5670772/glass/Analysis/code_direction_propensity/randomvector/*

*The propensity data for $t/\tau=15$ can be found on Odin in the folder /home/users/5670772/glass/Analysis/KA/KA_1/run_1001/*

*The mathematica notebook can be foun on Odin in the folder /home/users/5670772/glass/Analysis/KA_Dmin/Escaped particles.nb*

For 60 consecutive snapshots (obtained via a LAMMPS simulations that writes out snapshots every $t/\tau =15$) we perform 50 simulations (from now on called *pruns*) starting from the same snapshot, but with different velocities, over a time period of $t/\tau = 15$. Per snapshot we classify which particles have escaped. To do so, we calculate $D_\text{min}$ for each particle in each run and then classify a particle as being 'escaped' when it has a $D_\text{min}>D^c_\text{min}$ in more than 50\% of the Pruns.

:::warning
Note that when using $D^c_\text{min}=0.28$ , almost no particles escaped (only $\approx 100$ out of 96000). This makes sense because the escape criterium is very strict. For the rest of this analysis we therefore slightly alter the escape criterium: $D_\text{min}>0.20$ in more than 50\% of the runs. Note that this is still a strict criterium: We still only select $\approx 400$ particles. This means that we still only consider the most mobile particles.
:::

### Escape direction
We first investigate how much the escaped particles move in the same direction in the differen Pruns. For this we consider the inproduct between the vectors associated with the escaped particles in the different Pruns. I.e. we consider the following algorithm

1. For each escaped particle we select the Pruns where the particle has $D_\text{min}>0.20$.
2. Assume there are $N_E$ Pruns where the particle escaped. For these Pruns we compute the unit vector pointing in the direction between the IS and the position at $t/\tau=15$ (so for particle *i* we have a collection of unit vectors $\{\hat{n}_j^i\}$ where $j\in N_E^i$ and $N_E^i$ is the number of Pruns where particle *i* escapes).
::: danger
Note that I use the inherent state instead of the cage state to define the escape direction on. We know that in principle the cage state does a better job in finding the cage centers. However, since we perform the analysis for a lot of different snapshots, doing the MC simulation is a very expensive is computationally more demanding.

Moreover, it is the question whether, in the case of almost escaping particles, the CS is better than the IS. In the CS particles are constrained to their voronoi cell. For particles that are on the verge of escaping, we could expect the voronoi cell to be less well defined. For these particles the IS thus might be a better point of reference.

:question: What would be a good way to investigate this? The problem is that computing the CS is computationally so expensive. One possibility, I thought, might be to consider the difference between IS and CS (or initial state and cage state) for one run and then see how strongly this correlates with the associated $D_\text{min}$. If there is a strong correlation it shows that for escaping particles the cage center might be ill defined. Note sure though about this. Do you have better ideas?:question:
:::
4. Then we compute the following list for each particle *i*
$$
I^i= \{\hat{n}^i_s\cdot \hat{n}^i_t\} \text{ for } s \in [1, N_e] \text{ and } t \in[s+1, N_e]
$$

Note that $I^i$ gives a distribution of the inproduct between the different vectors associated with different Pruns. If particles always escape in the same direction we expect the inproduct to be always close to 1.

Below we plot the inproduct distribution (both normal [left] and cumulative [right]) for all escaped particles, together with the distribution you would get when the particles escaped in random directions.
<img src="https://i.imgur.com/hoMGxl8.png
" alt="drawing" width="300" /><img src="https://i.imgur.com/yOUhqRK.png
" alt="drawing" width="320" />

As we can see the distribution for the escaped particles is peaks around one, and is significantly different from the distribution you obtain when you look at random particles.

The same conclusion can be drawn when we look at the length of the average escape vector of each particle, i.e.
$$
v^i_M = \frac{|\sum_{s = 0}^{N_e}\hat{n}^i_s|}{N_e}
$$
Below we plot the histogram of the $v_M$ distribution of the escaped particles. The black line shows the value that we would expect when the vectors are distributed randomly.

From both analysis we can conclude that particles on average move in certain directions.
::: warning
The question is what to do with this? We clearly see that there is a preference in the escape direction of particles, which I think is an interesting finding. However, since it will be very hard to collect enough data to train models on, I think it will be not possible to make any real claims about the relation between structure and escape direction.

Maybe it would be possible to predict which particles are likely to soon escape, however in that case you are almost predicting the softness of Andrea Liu...

What do you think?
:::

### Escape direction w.r.t. system
Finally we look at where the fast moving particles lie with respect to what later will be the fast and slow moving area's. Below you see in the most left panel the IS configuration of the system at $t/\tau =0$, where the particles are coloured by the propensity at $t/\tau\approx \tau_\alpha$. On the right we color the particles that have escaped. Their color indicates at what time they have escaped (ranging from $t/\tau=0$ to $t/\tau=900$, see bar). Note that the position of the escaped particles is given by the IS at the associated times. Some particles that have escaped multiple times are therefore displayed multiple times. Note furthermore that we only consider large particles in the escape analysis.

We also look at the mean escape vector $\hat{v}^i_M = \frac{\sum_{s = 0}^{N_e}\hat{n}^i_s}{N_e}$ w.r.t. to the IS of the particles. The escape vectors with respect to the IS is displayed on the left (note that the vectors are made twice as long to increase visibility). Finally we also plot everyting together o the right.

From these figures we can conclude multiple things:
- Even particles that escape at very early times are located in the fast moving area. Note that this is quite interesting: We base this analysis on essentially a single run (the consecutive snapshots come from one continuous simulation). This means that the slow areas are not just on average slow, instead the slow particles truly never move.
- We see that escaped particles of the same colour are grouped together, i.e. we see cluster moves. This can also be seen if we look at individual time steps, f.e. $t/\tau=885$ as shown below.
<img src="https://i.imgur.com/RcPmRJS.png
" alt="drawing" width="364" />

I find it hard to conclude something qualitativ about the direction particles go w.r.t the slow or fast moving area's

## Results difference IS and CS
As mentioned above, we are also interested in which state (CS or IS)is better in capturing the cage center for the particles that nearly escape. Here we make a careful start of that analysis by looking at the correlation between $D_\text{min}$ and the difference between initial and CS/IS.

#### Initial and IS
We look at both the correlation between $|\mathbf{r}^{IS}-\mathbf{r}^{init}|$ and the average $D_\text{min}$ and the median $D_\text{min}$ value of each particle:

$\rho(|\mathbf{r}^{IS}-\mathbf{r}^{init}|, \langle D^\text{average}_\text{min}(i, t =15, \Delta = 15)\rangle ) = 0.56$
$\rho(|\mathbf{r}^{IS}-\mathbf{r}^{init}|, \langle D^\text{median}_\text{min}(i, t =15, \Delta = 15)\rangle ) = 0.57$
#### Initial and CS
$\rho(|\mathbf{r}^{CS}-\mathbf{r}^{init}|, \langle D^\text{average}_\text{min}(i, t =15, \Delta = 15)\rangle ) = 0.71$
$\rho(|\mathbf{r}^{CS}-\mathbf{r}^{init}|, \langle D^\text{median}_\text{min}(i, t =15, \Delta = 15)\rangle ) = 0.75$
#### CS and IS
$\rho(|\mathbf{r}^{CS}-\mathbf{r}^{IS}|, \langle D^\text{average}_\text{min}(i, t =15, \Delta = 15)\rangle ) = 0.53$
$\rho(|\mathbf{r}^{CS}-\mathbf{r}^{IS}|, \langle D^\text{median}_\text{min}(i, t =15, \Delta = 15)\rangle ) = 0.51$

#### Comparison
As we can see there is a clear correlation between $|\mathbf{r}^{CS}-\mathbf{r}^{init}|$ and both the average and median value of $D_\text{min}$. That this correlation is the strongest is something we expect: $D_\text{min}$ is a good predictor of the dynamics during the caging regime, and therefore we also expect it to capture the non-affine movement.

For this analysis we are interested more in the extreme points, i.e. the points associated with a very large $D_\text{min}$. Although it is impossible to base any real statistics on only a few point, we can make their behaviour more insightfull by plotting the initial state (black), IS (blue), CS (red)and the average position after escape (yellow)
<img src="https://i.imgur.com/tYsuEqB.png
" alt="drawing" width="200" /><img src="https://i.imgur.com/CbXZfXP.png
" alt="drawing" width="200" /><img src="https://i.imgur.com/vbZJhhy.png
" alt="drawing" width="200" />

As we can see, for both the first and the third of these escapes, the CS is already close to being escaped. This means that we can carefully conclude that the CS is probably does not capture the cage centre very well for the nearly escaped particles. Moroever, as we can see, the escape direction is in general influenced by what you call your cage centre. In that sense it is the question whether we can really talk about escape directions.
::: info
Note that with this I do not mean that particles will not move in a certain direction when they escape. I mean that it matters w.r.t. what origin you define the escape direction.
:::