# Roadmap paper glassy systems
###### tags: `PhD projects`

For the glass roadmap paper we want to apply our caging methodology on two different systems (3D KA from Shiba and 2D KA of Jung). Note that to measure how well the predicted propensities correspond to the propensities, we use the Pearson correlation.

## 3D KA
*Data can be found on ODIN /home/users/5670772/glass/Road_map_paper/data_shiba_3DKA/*

The Shiba dataset can be found here: https://github.com/h3-Open-BDEC/pyg_botan and the paper can be found here: https://arxiv.org/pdf/2206.14024.pdf.

The binary Kob Anders potential is given by \begin{equation}V_{ab}(r)= 4\epsilon_{ab}\left[\left(\frac{\sigma_{ab}}{r}\right)^{12}-\left(\frac{\sigma_{ab}}{r}\right)^{6}\right],\tag{1}\end{equation}

Where in this case the potential is truncated and shifted such that

$$
U(r) = u(r)- u(r_c)-(r-r_c)\frac{du(r)}{dr}|_{r= r_c},\\
U(r) = 4\epsilon_{ab}\left[\left(\frac{\sigma_{ab}}{r}\right)^{12}-\left(\frac{\sigma_{ab}}{r}\right)^{6}- \left(\frac{1}{2.5}\right)^{12}+\left(\frac{1}{2.5}\right)^{6}\right] \\+ \frac{24\epsilon_{ab}(r- 2.5\sigma_{ab})}{2.5\sigma_{ab}}\left[2\left(\frac{1}{2.5}\right)^{12}-\left(\frac{1}{2.5}\right)^{6}\right]
$$

where $r_c = 2.5\sigma_{\alpha\beta}$ and where $(\epsilon_{aa}, \epsilon_{ab}, \epsilon_{bb}) = (1, 1.5,0.5)$ and $(\sigma_{aa}, \sigma_{ab}, \sigma_{bb}) = (1, 0.80,0.88)$. With composition $N_A:N_B = 0.8:0.2$ and $n=4096$. Time is measured in units of $\tau = \sqrt{\frac{m\sigma_{AA}}{\epsilon_{AA}}}$. The system has $N/V\sigma_{AA}^3 = 1.2$ and $T/(k_BT) = 0.44$, such that $\tau_\alpha/\tau = 4120$.

Note that this is a slightly different potential than the we used before (previously we did not have $U(r_\text{cut}) = 0$). Below we plot all the potentials and corresponding forces $AA$ particle interactions. . The red line corresponds to the bare potential, the blue line to the shifted and truncated potential we used before and the yellow line to the potential used by Shiba. Note that Shiba's potential is less negative, and therefore it is expected that both the potential energy and pressure of the system are slightly higher.
| $U/\epsilon$ | $P\sigma^3/\epsilon$ |
| ------------------------------------ | ------------------------------------ |
|||

:::info
We checked that this potential, as well as its derivative is zero at $r = 2.5\sigma_{ab}$
:::

### Data Shiba
*Found on ODIN /home/users/5670772/glass/Road_map_paper/data_shiba_3DKA/code_data_Shiba_to_ours/.*

The propensity, that is measured in the isoconfigurational ensemble, is measured by averaging over 32 runs. There are 500 initial configurations, 400 of those are used as training and 100 for testing. We checked whether when we convert Shiba's data to ours, we get back something that looks like a proper propensity (see graph below).

| Propensity | |
| -------- | -------- |
|| |
### MC simulation to measure cage state
*Checks from this code can largely be found in the data of the previous paper, see ODIN /home/users/5670772/glass/Analysis/code_Monte_Carlo/CHECKS/KA.*

We choose a restriction radius of $r/\sigma_{AA}= 1.25$, since we saw for the KA system that that worked the best (see SI of the caging paper).

#### Checks Monte Carlo
*Found on Odin /home/users/5670772/glass/Road_map_paper/data_shiba_3DKA/code_Monte_Carlo/CHECKS/*

For the MC code, we use the earlier used MC code for finding the cage state. In this code we measure the the pressure and the potential energy. For our previous system this was equal to respectively $-7.031\pm0.0132904$ and to $3.07877\pm0.0885674$ (see hackMD on Glass simulations KA and harmonic), however for reasons discussed before we expect both te pressure and potential to be slightly higher. Below we plot the results. Note that these agree with the potential energy found by Jung.
| $U/\epsilon$ | $P\sigma^3/\epsilon$ |
| ------------------------------------ | ------------------------------------ |
| | |

### BOPS
*See ODIN /home/users/5670772/glass/Road_map_paper/data_shiba_3DKA/code_BOPS/. Pictures are made with the Mathematica notebook that can be found in the folder CHECKS.*

We use the same parameters to capture the structure as we used for the KA system in the caging paper, i.e. :
* 366 in total
* 198 radial functions with $s\in{A,B}$:
* 60 equally spaced in the interval $r/\sigma_A\in(0.5,2.0]$ with $\delta=0.025$.
* 20 equally spaced in the interval $r/\sigma_A\in(2.0,3.0]$ with $\delta=0.050$.
* 19 equally spaced in the interval $r/\sigma_A\in(3.0,4.9]$ with $\delta=0.100$.
* 168 Angular fuctions
* 14 equally spaced in the interval $r/\sigma_A\in[1.0, 2.3]$ with $\delta=0.1$ and $l\in[1,12]$

See results below for one snapshot (one particle and average over all particles). This is the same code that we normally use for the KA BOPS.
| | Radial functions | BOPS |
| --- | --------------------------------------------- | --------------------------------------------- |
| Particle 1 | |  |
| Average | |  |

### Results
*The left picture and analysis can be found on ODIN /home/users/5670772/glass/Road_map_paper/data_shiba_3DKA/Linear_regression/CHECKS/ in the notebook Pearson correlation data Shiba and our data.nb. The right picture and analysis can be found on ODIN /home/users/5670772/glass/Road_map_paper/data_shiba_3DKA/Linear_regression/Analysis/ in the noteboek ResultsShiba.nb*

We train the data on 100 snapshots from the train set of Shiba and then test it on the 100 test sets. Below se show the results.

| Comparsion our data | Results |
| -------- | -------- |
| | |

#### CHECKS
*See Pearson correlation data Shiba and our data.nb on ODIN: /home/users/5670772/glass/Road_map_paper/data_shiba_3DKA/Linear_regression/CHECKS/*
We checked whether the non-normalized propensities give the same results as the normalized propensities.
## 2D KA

The second model is the 2D KA model used by Jung. The paper can be found here (https://arxiv.org/pdf/2210.16623.pdf) and the SI here (https://cloud.coulomb.umontpellier.fr/index.php/s/GZFLCi7BckiMzsC?dir=undefined&openfile=10671355). The data set can be found here (https://cloud.coulomb.umontpellier.fr/index.php/s/GZFLCi7BckiMzsC).

They consider a 2D KA mixture with 3 types, $N_1 = 600$, $N_2 = 330$ and $N_3 = 360$ (i.e. $N = 1290$) in a box of $L= 32.896 \sigma_{11}$ at probably $k_BT/\epsilon_{11}=0.3$. They use the following potential

$$
V_{\alpha\beta}(R_{ij}) = \begin{cases}4\epsilon_{ab}\left[\left(\frac{\sigma_{ab}}{R_{ij}}\right)^{12}-\left(\frac{\sigma_{ab}}{R_{ij}}\right)^{6}+C_0+C_2\left(\frac{R_{ij}}{\sigma_{\alpha\beta}}\right)^2 + C_4\left(\frac{R_{ij}}{\sigma_{\alpha\beta}}\right)^4\right]&R_{ij}<R^\text{cut}_{\alpha\beta}\\
0 &\text{otherwise}\end{cases}\\
$$
Here $\alpha\beta = \{1,2,3\}$ and $\epsilon_{11} =1.0$, $\epsilon_{12} = 1.5$, $\epsilon_{22} = 0.5$, $\sigma_{11} = 1.0$, $\sigma_{12} = 0.8$, $\sigma_{22} = 0.88$ (i.e. normal KA mixture). Then we also consider $\epsilon_{13}= 0.75$, $\epsilon_{23}= 1.5$, $\epsilon_{33}= 0.75$, $\sigma_{13} = 0.9$, $\sigma_{23} = 0.8$ and $\sigma_{33} = 0.94$ and $R^\text{cut}_{\alpha\beta} = 2.5$.

Furthermore we set $C_0 = 0.04049023795$, $C_2 = −0.00970155098$ and $C_4 = 0.00062012616$, such that the potential is continuous up to second derviative. This means that the potential does not have to be truncated/shifted.

Below we plot the potential for $AA$ particle interactions.

The propensity that is measured in the isoconfigurational ensemble is measured by averaging over 20 runs. They consider three temperatures, of which we will consider $k_BT/\epsilon_{AA} = 0.23$ and $k_BT/\epsilon_{AA} = 0.30$. Predictions are made for particle type 1.

### Data Jung
*Found on ODIN /home/users/5670772/glass/Road_map_paper/data_jung_2DKA/code_data_Jung_to_ours/.*

We check whether when we convert Jung's data format to ours, we indeed get back something that lookds like a propensity the propensity.

| Propensity |
| -------- |
||
### MC simulation to measure cage state
*Found on ODIN /home/users/5670772/glass/Road_map_paper/data_jung_2DKA/code_Monte_Carlo/*

For the MC simulation, we use the same MC code as before, however we changed it to 2D and implemented it for 3 particle types

#### Checks Monte Carlo
*Found on ODIN /home/users/5670772/glass/Road_map_paper/data_jung_2DKA/code_Monte_Carlo/CHECKS/, the MC_test_Upot_KA2D.c code contains the MC code, rewritten such that it computes the energy of every initial snapshot of reach run of the trainset.*

From Gerhard Jung we know that the average potential energy per particle in the model at $T=0.3$ is $-2.594\pm0.001$ (and -2.886+/-0.003 in the inherent states) for all configurations. If we take into account all the initial snapshots of the training set using our MC simulation (where we only compute our initia energy), we indeed find:

| $U/\epsilon = -2.59373 \pm 0.00933833$ | $P\sigma^3/\epsilon = 5.3249\pm5.3249$ |
| ------------------------------------ | ------------------------------------ |
|| |

### BOPS and radial function
(https://journals.aps.org/pre/abstract/10.1103/PhysRevE.67.031608)
The BOPS only work in 3D, so we turn to the 2D variant, which are defined as
$$
\Phi^l_m = \frac{1}{N_b}\sum_{n = 1}^{N_b}e^{li\theta_{mn}}
$$.
Where $\theta_{mn}$ is the angle between some fixed axis (e.g., $x$- or $y$-axis) and the bond joining the mth particle with another nth neighboring particle.
Note that we consider $l\in\{1,12\}$. In order to make them, as much as possible, comparable to our 3D BOPS, we slightly alter the to:

$$
\phi^l_i (r, \delta) = \frac{1}{Z}\sum_{i\neq j}e^{\frac{-(r_{ij}-r)^2}{2\delta^2}}e^{li\theta_{ij}},
$$

with
$$
Z = \sum_{i\neq j}e^{\frac{-(r_{ij}-r)^2}{2\delta^2}}
$$
And then we compute
$$
\Phi^l_i(r, \delta) = \sqrt{\phi^l_i\cdot (\phi^l_i)^*}
$$

Of course we also consider the radial density, defined as
$$
G_m(r, \delta, s) =\sum_{i\neq j, s_j = s} e^{\frac{-(r_{ij}-r)^2}{2\delta^2}},
$$

Obtaining higher order parameters, we just use the same averaging method as we normally do.

### Paircorrelation function
*Found on ODIN /home/users/5670772/glass/Road_map_paper/data_jung_2DKA/code_paircorrelation/*

In order to know to which distance we want to include paramaters, we measure the paircorrelation function for respectively 500 snapshots ($T=0.30$) and 150 snapshots ($T=0.23$) of the trainingdata.

| $T=0.23$| $T=0.30$ |
| -------- | -------- |
| ||

Since the $g(r)$ looks very similar for both temperatures, we can include parameters up to the same distance for both temperatures.

### BOPS
*Found on ODIN /home/users/5670772/glass/Road_map_paper/data_jung_2DKA/code_BOPS_2D/*
As usual we take into account particles up to the second minimum (BOPS), and the fifth minimum (radial parameters). Note furthermore that for the radial parameters we know have three values per distance $r$, since we consider three different particle types.

Leading to the following parameters
* 462 in total
* 294 radial functions with $s\in{A,B}$:
* 60 equally spaced in the interval $r/\sigma_A\in(0.5,2.0]$ with $\delta=0.025$.
* 20 equally spaced in the interval $r/\sigma_A\in(2.0,3.0]$ with $\delta=0.050$.
* 18 equally spaced in the interval $r/\sigma_A\in(3.0,4.8]$ with $\delta=0.100$.
* 168 Angular fuctions
* 14 equally spaced in the interval $r/\sigma_A\in[1.0, 2.3]$ with $\delta=0.1$ and $l\in[1,12]$

See results below for one snapshot (one particle and average over all particles)
| | Radial functions | BOPS |
| --- | --------------------------------------------- | --------------------------------------------- |
| Particle 1 |  |  |
| Average |  |  |

The data in these first two figures is evaluated for a type 1 particle. This means that we expect the first peaks at slightly higher than $1.0$ (11), $0.8$ (12) and $0.9$ (13), since $\sigma_{11} = 1.0$, $\sigma_{12} = 0.8$ and $\sigma_{13} = 0.9$.

### Results
*Found on /home/users/5670772/glass/Road_map_paper/data_jung_2DKA/Linear_regression/Analysis/ in the notebook ResultsJung.nb*

First of all we looked at different restriction radii ($r_c/\sigma_{AA} \in\{ 1.00, 1.25, 1.50\}$) and concluded that just like in the 3D case $r_c/\sigma_{AA}= 1.25$ yields the best results.


We train models based on different sets of parameters:
* $X^\text{init}$ are the parameters based on the initial configuration.
* $X^\text{IS}$ are the parameters based on the inherent state.
* $X^\text{CS}$ are the parameters based on the case state.
* $\Delta r^\text{IS}$ is the absolute difference between particles in the initial and the inherent state.
* $\Delta r^\text{CS}$ is the absolute difference between particles in the initial and the cage state.
* $\text{All}$ conclude respectively for the left and right picture the set $\{X^\text{init}, X^\text{IS},\Delta r^\text{IS}\}$ and $\{X^\text{init}, X^\text{CS},\Delta r^\text{CS}\}$.

Below we show the results. The first row shows the correlation obtained for $k_bT/\epsilon_{AA} = 0.23$, whereas the second row shows the correlations obtained for $k_bT/\epsilon_{AA} = 0.30$.

| | Inherent state|Cage state|
| -------- | -------- |-------- |
| $T = 0.23$|||
| $T = 0.30$| | |

As we can see, for both temperatures the cage state is a better predictor for the propensity than the initial state. Moreover, it turns out that the temperature has a significant influence on how well we can capture the cage state. For $k_bT/\epsilon_{AA} = 0.23$ we are able to predict the propensity with a lot more accuracy (note that the same is true for how good of a predictor the IS is for the propensity).

### Lower correlation 2D system
*See notebook Predictedpropvsrealprop.nb on ODIN /home/users/5670772/glass/Road_map_paper/data_jung_2DKA/Linear_regression/Analysis/*
Possible causes:
* 2D instead of 3D
* More particle types, i.e. less well defined cage center. We looked at correlation between which particletypes surround the particle and the difference between measured and predicted propensity in the caging regime. However this does correlation is not significant
* Less Pruns (only 20), so maybe the cage center is less well defined.

## TO DO
Shiba
- [ ] Run system for different temperatures and look at transferibility

Jung
- [ ] Do unsupervised learning (look at 2D BOPS for one shell of nearest neighbours and averaged over also one shell of neughbours) and then do PCA