# Glass simulations KA and harmonic potential
###### tags: `PhD projects`
## General
### KA
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 $(\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$.

The corresponding force is given by
$$\mathbf{F}_{ab}(\mathbf{r})= -\frac{dV_{ab}}{d\mathbf{r}}=\frac{24\epsilon_{ab}}{r^2}\left[2\left(\frac{\sigma_{ab}}{r}\right)^{12}-\left(\frac{\sigma_{ab}}{r}\right)^{6}\right]\mathbf{r]} $$

### Harmonic potential
The binary harmonic potential is given by \begin{equation}V_{ab}(r)= \epsilon_{ab}(1-\frac{r}{\sigma_{ab}})^2,\tag{2}\end{equation}
where $(\epsilon_{aa}, \epsilon_{ab}, \epsilon_{bb}) = (1, 1,1)$ and $(\sigma_{aa}, \sigma_{ab}, \sigma_{bb}) = (1, 1.2,1.4)$. With composition $N_A:N_B = 0.5:0.5$.

The corresponding force is given by
$$\mathbf{F}_{ab}(\mathbf{r})= -\frac{dV_{ab}}{d\mathbf{r}}= \frac{2\epsilon_{ab}}{\sigma_{ab}}(1-\frac{r}{\sigma_{ab}})\frac{\mathbf{r}}{r},$$

## LAMMPS simulation
### CHECKS
#### KA simulation
##### Monodisperse Lennard Jones system:
*Found in folder glass/CHECKS/Monodisperse_MC_LJ and glass/CHECKS/lammps_LJ on HERA*

First we check for a mono disperse Lennard Jones system. We compare the statepoints of a Lennard Jones system with different simulational models. The difference between 'Lammps' and' 'Lammps sf potential' is that the first one has a lj/cut potential which is shifted, while the second one has a lj/sf potential (i.e. with a tail correction and no shift).

The literature also uses a otential with tail correction, while the MC code truncate and shifts. Note also that MC measure the pressure in terms of $P\sigma^3/(k_BT\epsilon)$, meaning that to compare it with the other pressures we have to multiply it with $k_BT$. Pressure is computed making use of the Virial expansion: $P/k_BT= \rho_n+\frac{\beta}{3V}\langle\sum_{i<j}\mathbf{f}_{ij}\cdot \mathbf{r_{ij}}\rangle$

| | | $k_BT/\epsilon =0.75$, $\rho\sigma^3=0.05$ | $k_BT/\epsilon = 1.00$, $\rho\sigma^3=0.01$ | $k_BT/\epsilon = 1.50$, $\rho\sigma^3=0.1$ |
| --- | --- | -------------------------------- | ------------------------------------------- | --------------------------------------------
| Literature* |$U/\epsilon$| $-0.4535$| $-0.067272$| $-0.5241$|
| | $P\sigma^3/\epsilon$| $0.02627$| $0.009677$| $0.1349$|
| Lammps sf potential |$U/\epsilon$ |$-0.452999\pm0.0174796$|$-0.06725\pm0.00477289$|$-0.523795\pm0.0127503$ |
| | $P\sigma^3/\epsilon$|$0.0267354\pm0.0017442$| $0.00968903\pm0.000233983$ |$0.13484\pm0.00617984$|
| Lammps |$U/\epsilon$|$-2.63335\pm0.0681041$|$-0.0799801\pm0.00509921$|$-0.621197\pm0.0138591$
| | $P\sigma^3/\epsilon$ |$0.0110137\pm0.00437606$|$0.00957315\pm0.000253802$|$0.12668\pm0.00619963$|
| MC |$U/\epsilon$|$-2.40476\pm0.181936$**|$-0.0800966\pm0.00534845$|$-0.620373\pm0.0144752$|
| | $P\sigma^3/\epsilon$|$0.0160026\pm0.00520$**|$0.00958366\pm0.000150328$|$0.127176\pm0.00544109$|

*https://www.nist.gov/mml/csd/informatics/lammps-md-equation-state-pressure-vs-density-linear-force-shifted-potential-25s (Uses Tail correction)

**This point is very unstable (see data in the link), probably MC has a hard time equilibrating

##### Binary Lennard Jones
*Found in folder glass/CHECKS/Binary_MC_KA and glass/CHECKS/lammps_KA on HERA*

The second check we do is for a binary lennard Jones system, i.e. the 80:20 KA mixture (see * ). Note that this paper does truncate and shift the potential, so that is what we also do in Lammps.

| | | $k_BT/\epsilon =0.44$, $\rho\sigma^3=1.203$ | $k_BT/\epsilon =0.8$, $\rho\sigma^3=1.203$ | $k_BT/\epsilon = 1.0$, $\rho\sigma^3=1.203$ |
| --- | ------------------------- | -------------------------------------------- | ------------------------------------------- | -------------------------------------------- |
| Literature* | $U/\epsilon$ | $-7$ | $-6.3$ | $-6$ |
| | $P\sigma^3/\epsilon$ | $3$ | $8$ | $10$ |
| Lammps |$U/\epsilon$|$-7.031\pm0.0132904$|$-6.34942\pm0.0223181$|$-6.01863\pm0.0278269$|
| |$P\sigma^3/\epsilon$|$3.07877\pm0.0885674$|$7.76942\pm0.146618$|$10.0505\pm0.17761$|
|MC |$U/\epsilon$|$-7.01117\pm0.0129798$|$-6.3497\pm0.0222575$|$-6.01937\pm0.0278651$ |
| | $P\sigma^3/\epsilon$|$3.25876\pm0.0879145$|$7.77107\pm0.147007$|$10.0418\pm0.17906$ |

*https://journals.aps.org/pre/pdf/10.1103/PhysRevE.51.4626 (Note that literature points are read from a graph, i.e. the corresponding values are only approximate).

#### Harmonic system simulation
*Found in folder glass/CHECKS/Binary_MC_har and glass/CHECKS/lammps_harmonic on HERA*

All statepoints are chosen such that $\tau_\alpha \approx 10^4$ (measured in units of $\tau = \sqrt{\frac{m\sigma_A^2}{\epsilon}}$).

Note that the potential in Lammps is defined ambiguously, namely: $v= k(r_c-r)^2$, with *r* in length units and *k* in energy units. We therefore use two different potentials, one where $k=\epsilon$ as defined in equation $(2)$ and one where $k = \epsilon /r_c^2$. As we can see, the latter yields results that agree with the MC simulation.

| | | $k_BT/\epsilon =0.00065$, $\rho\sigma^3=0.65$| $k_BT/\epsilon =0.0027$, $\rho\sigma^3=0.72$ | $k_BT/\epsilon =0.0045$, $\rho\sigma^3=0.82$|
| --- | -------------------- | -------------------------------------------------------------------- | ------------------------------------------------------------------- | ------------------------------------------------------------------- |
|Lammps $k=\epsilon$|$U/\epsilon$|$0.000538002\pm0.0000145836$|$0.00615876\pm0.0000902119$|$0.0251965\pm0.000145381$|
||$P\sigma^3/\epsilon$|$0.0119541\pm0.000212097$|$0.0629642\pm0.000517721$|$0.173745\pm0.000539723$|
|Lammps $k=\frac{\epsilon}{r_c^2}$|$U/\epsilon$|$0.000591555\pm0.0000163006$|$0.00578502\pm0.0000919411$|$0.0186804\pm0.000215572$|
||$P\sigma^3/\epsilon$|$0.0104195\pm0.000181439$|$0.0478648\pm0.000415494$|$0.113462\pm0.000634303$|
| MC| $U/\epsilon$|$0.000586436\pm0.0000151196$|$0.0057897\pm0.0000910177$|$0.0188086\pm0.000150003$|
||$P\sigma^3/\epsilon$|$0.0103451\pm0.000165723$|$0.0478933\pm0.000396463$|$0.113851\pm0.000465065$|

### Collecting data

We consider the following statepoints

| | $\rho\sigma^3$ | $k_BT/\epsilon$ | $\tau_\alpha/\tau$* | Mixture |
| -------- | --------------- | --------------- | ------------------ | --- |
| KA |1.203|0.44|4636|(80:20)|
| Harmonic |0.82|0.0045|10000|(50:50)|

*Note that $\tau = \sqrt{\frac{m\sigma_A^2}{\epsilon}}$
1. We start 100 simulations. In each simulation 2000 particles are placed randomly.
2. First an energy minimization is executed in order to avoid overlap.
3. Afterwards we perform an NVT simulation of $t = 10\tau_\alpha$, where $dt = 0.002$ and with a Nose Hoover thermostat with damping $dT=200*dt$.
4. This NVT simulation is followed by an NVE simulation that again runs for $t=10\tau_\alpha$.
5. Using this protocol provides us with 100 equilibrated snapshots. For each of these snapshots we subsequently perform 50 simulations in an NVE ensemble, with the following characteristics.
- The initial velocity of particles is drawn from a Gaussian distribution at the desired $k_BT/\epsilon$. The velocities of the particles is then shifted such that the overal linear momentum is zero.
- We output data at logarithmic intervals (i.e. $t_n = 1, 2, 3, 5, 10, 11, 12, etc.$). At each timestep a coordinate file is written out, which contains the unwrapped positions of particles (i.e. not the position of the particles projected back in the box). We save the unwrapped coordinates in order to keep track of the CofM drifting.
- The simulation runs for $10\tau_\alpha$
6. Afterwards the propensity is computed by averaging over the meandisplacemt of particles in these 50 runs.
## Propensity
### CHECKS
*Found in folder glass/Analysis/code_Propensity/CHECKS/ on HERA*
We perform two checks:
1. Whether the Center of Mass stays constant during the runs. This turns out to be the case (see mathematica notebook )
2. Whether the propensity of the systems looks normal. Time is measured in units of $\tau = \sqrt{\frac{m\sigma_A^2}{k_BT}}$
* Harmonic
<img src="https://i.imgur.com/qIgRSSD.png
" alt="drawing" width="364" />
* KA
<img src="https://i.imgur.com/C5PpXjR.png
" alt="drawing" width="300"/><img src="https://i.imgur.com/eAtKkH2.png
" alt="drawing" width="360" />

The left picture shows the propensity as found in the BAPST paper (*Unveiling the predictive power of static structure in glassy systems*). The darkest blue line corresponds to *A* particles in a system at the statepoint we are looking at. The right picture contains the propensity as we measure it in our simultion. As we can see the two agree with each other (small difference in time scale are explained by the fact that Bapst measures time in units of $\tau = \sqrt{\frac{m\sigma_A^2}{\epsilon}}$ and not in units of $\tau = \sqrt{\frac{m\sigma_A^2}{k_BT}}$)

Conversion of $t$: In units of \epsilon we have $\tau = \tau^*\tau_\epsilon$, where $\tau^*$ is the dimensionless time and $\tau_\epsilon$ are the units in which we measure, where $\tau_\epsilon= \ \sqrt{\frac{m\sigma^2}{\epsilon}}$ we then convert it to $\beta$ units by computing $\tau = \tau^*\cdot\sqrt{\frac{k_BT}{\epsilon}}\sqrt{\frac{\epsilon}{k_BT}}\tau_\epsilon = \sqrt{\frac{k_BT}{\epsilon}}\tau_\beta$, where $\tau_\beta = \sqrt{\frac{m\sigma^2}{k_BT}}$. I.e. to convert to $\beta$ units, we need to multiply the $\tau^*$ from the $\epsilon$ units with $\sqrt{\frac{k_BT}{\epsilon}}$.

## Restricted Monte Carlo simulation

This simulation is used to find the ensemble averaged coordinates of particles restricted to a certain spherical volume.
### Checks
#### Hard Spheres

We check whether the newly written MC code gives particle positions that agree with the HS simulation that was used earlier. In the Figure below, we see that this is indeed the case 

#### KA and harmonic
*Found in folder glass/Monte_Carlo/Checks/KA and glass/Analaysis/code_Monte_Carlo/Checks/Harmonic on HERA*
To check this, we perform a simulation with the MC restricted simulation code, where the restriction radius is set to more than half the bos size. If this simulation yields the same energy and pressure, we know the simulation is correct.

| | Statepoint | $U/\epsilon$ checked MC code | $U/\epsilon$ restricted MC code | $P\sigma^3/\epsilon$ checked MC code |$P\sigma^3/\epsilon$ restricted MC code |
| -------- | --------------------------------------- | ---------------------------- | ------------------------------- | ----------------------------------- | --- |
| KA | $k_BT/\epsilon = 0.44$, $\rho_N=1.203$ |$−7.01117\pm0.0129798$|$-7.02055\pm0.0119456$|$3.25876\pm0.0879145$|$3.19079\pm0.0835622$
| Harmonic | $k_BT/\epsilon = 0.0045$, $\rho_N=0.82$ |$0.0188086\pm0.000150003$|$0.0188015\pm0.000149946$|$0.113851\pm0.000465065$|$0.113764\pm0.000454576$|

### Correlation different restriction radii and propensity
For 1 snapshot we measure the correlation between the propensity and $|\Delta(\mathbf{r}_i^\text{MC}-\mathbf{r}_i^\text{init snap})|$, where the particles in the MC simulation are restricted to different radii. For the additive systems (HS and har), this restricted radius is simply a numer times the radius of the particle itself (where the radius of a harmonic particle is said to be respectively 1 for $A$ particles and 1.4 for $B$ particles). The average values are measured for $5\cdot 10^5$ init cycles and $ 10^6$ measure cycles

In the case of KA particles, we say that the radius of the $A$ particles is $r_A = \sigma_{AA}/2$ and $r_B = \sigma_{BB}/2$.
#### HS
<img src="https://i.imgur.com/li9v0rG.png
" alt="drawing" width="364" /><img src="https://i.imgur.com/dZEtQfh.png
" alt="drawing" width="364"/>

Note that it makes sense that compared to the figures in my thesis, the graphs for higher restriction radii are somewhat shifted to the left. In my thesis we ran the system for less initialization and less measure cycles, meaning that particles were less likely to escape their cage, even if they were restricted to a radius larger than te cage. If we run the system longer, more and more particles will escape, meaning that the peaks will shift to later times, i.e. to the right.

#### Har
<img src="https://i.imgur.com/G1YypEc.png
" alt="drawing" width="364"/><img src="https://i.imgur.com/svH4hdL.png
" alt="drawing" width="364" />
#### KA
<img src="https://i.imgur.com/AbEnA3h.png
" alt="drawing" width="364" /><img src="https://i.imgur.com/n4sLctZ.png
" alt="drawing" width="364"/>
The fact that we clearly see three peaks in the correlation between MC and propensity has probably to do with the three different regimes (ballistic, cage and diffusive), i.e.
* First peak: we see that particles are send to a specific direction due to forces of other particles.
* Second peak: caging regime
* Third peak : after caging how much do particles move.
### Correlation average position in voronoi cell and propensity
*Found in folder HERA/glass/Analysis/code_Voronoi/CHECKS/*

For 1 snapshot we measure the correlation between the propensity and $|\Delta(\mathbf{r}_i^\text{MC}-\mathbf{r}_i^\text{init snap})|$, where the particles in the MC simulation are restricted toit's initial voronoi cell. To determinde the voronoi cell, we first determine the nearest neighbours with the SANN algorithm. A particle is said to be in it's initial voronoi cell at time $t$ if
$$
\forall j \in \text{{nearest neighbours particle i}: }\frac{|\mathbf{r}_i(t>0)-\mathbf{r}_i(t=0)|}{|\mathbf{r}_i(t>0)-\mathbf{r}_j(t=0)|}\cdot\frac{\sigma_j}{\sigma_i} < 1
$$
#### HS
<img src="https://i.imgur.com/eDZBzC3.png
" alt="drawing" width="450"/>
#### Har

<img src="https://i.imgur.com/4YlZyZR.png
" alt="drawing" width="450"/>
#### KA
<img src="https://i.imgur.com/CyViPcA.png
" alt="drawing" width="450"/>

As we can see using the voronoi cell works very well.

## FIRE algorithm simulation
This simulation is used to find the inherent state of the system.
### CHECKS
*Found in the folder glass/Analaysis/code_Monte_Carlo/Checks/Analysis_checks*

#### HS
Since the new HS FIRE code has a particle dependent cutoff, the energy of the origional code and the new code will not match up. However, we checked that when in the new code the cutoff is set to $2.1/\sigma_A$ for all particle types, the energies match up. In the new code the cutoff is set to $1.35/\sigma_A$ (slightly bigger than after the first peak).

Below, we show the potential energy per particle and the kinectic energy of the total system as a function of time.

<img src="https://i.imgur.com/RhUmBA8.png
" alt="drawing" width="300"/><img src="https://i.imgur.com/9w0zLW4.png
" alt="drawing" width="300" />
<img src="https://i.imgur.com/sxvQjPj.png" alt="centered image" width="300" />

Displacement of particles between initial and quenched snapshot (dark green means small displacements, yellow indicates large displacement).

**NB!!!!** It turned out that the HS quench did not fully quench because of the cutoff. In the paper of Francesco Arceri * and Eric I. Corwin (http://corwinlab.uoregon.edu/publications/2020/Arceri-VibrationalPropertiesHard-2020.pdf) they say that one should only take into account the first shell of nearest neighbours, i.e. the potential should be cutoff at $\approx 1.35 \sigma_A$. However at this point the force is non-zero. As a result many particles end up just outside this cutoff. In order to fix this we shifted the force and with that the potential.

The force shift is such that at $r\sigma_A = 1.35$, $F=0$. Since $U = \epsilon \log(r-\sigma_{ij})$, the force is given by $\mathbf{F} = \frac{1}{r-\sigma_{ij}}\frac{\mathbf{r}}{r}$, i.e. $\mathbf{F}_\text{cut} = \frac{1}{r_c-\sigma_{ij}}\frac{\mathbf{r}}{r}$. And $\mathbf{F}_\text{res} = \frac{1}{r-\sigma_{ij}}\frac{\mathbf{r}}{r}-\mathbf{F}_\text{cut}$.

As a result the potential must be shifted with
$$
U_\text{res} = U + \mathbf{F}_\text{cut}\cdot \mathbf{r}=U + \frac{1}{r_c-\sigma_{ij}}r
$$ such that
$$
-\nabla U_\text{res} = \mathbf{F} - \frac{1}{r_c-\sigma_{ij}}\frac{\mathbf{r}}{r} = \mathbf{F}_\text{res}
$$
#### Harmonic
The initial potential energy of the FIRE simulations matches with the average potential energy of the system $U_\text{pot}/\epsilon(t=0)=0.018062$ (*see the folder initial_energy*).

Below, we show the potential energy per particle and the kinectic energy of the total system as a function of time for the snapshot of the the first dataset and first run (N.B. this is a different snapshot than the snapshot for which the potential energy was checked).

<img src="https://i.imgur.com/BpaU5tp.png
" alt="drawing" width="340"/><img src="https://i.imgur.com/pGcTMkB.png
" alt="drawing" width="300" />

<img src="https://i.imgur.com/zkEqowu.png" alt="centered image" width="300" />

Displacement of particles between initial and quenched snapshot (dark green means small displacements, yellow indicates large displacement).

#### KA
The initial potential energy of the FIRE simulations matches with the average potential energy of the system $U_\text{pot}/\epsilon(t=0)=-7.02584$ (*see the folder initial_energy*).

Below, we show the potential energy per particle and the kinectic energy of the total system as a function of time for the snapshot of the the first dataset and first run (N.B. this is a different snapshot than the snapshot for which the potential energy was checked).

<img src="https://i.imgur.com/vuqPezu.png
" alt="drawing" width="340"/><img src="https://i.imgur.com/sMVNLtH.png
" alt="drawing" width="300" />
<img src="https://i.imgur.com/djqv7Ls.png" alt="centered image" width="300" />

Displacement of particles between initial and quenched snapshot (dark green means small displacements, yellow indicates large displacement).

### Correlation quench and propensity for different maxtimes
*Found in folder HERA/glass/Analysis/code_FIRE/CHECKS/*

We compute the correlation between the absolute distance between the initial and the quenched snapshot for multipled maxtimes. As we can see for smaller maxtimes, the particle has a higher correlation with the propensity in the caging reigme. For larger maxtimes, the correlation increases for later times.

#### HS
<img src="https://i.imgur.com/WeGqMQ4.png" alt="drawing" width="300" /><img src="https://i.imgur.com/57mOD2A.png
" alt="drawing" width="350"/>

#### Harmonic
<img src="https://i.imgur.com/VusLPkr.png" alt="drawing" width="300" /><img src="https://i.imgur.com/HASvvBf.png
" alt="drawing" width="350"/>

If we look at the distribution of the distance between initial snapshot and inherent state, we also see that the quenched snapshot really overshoots. From this we can also understand that the quenched snapshot sometimes has a higher correlation for later times: There will be some correlation with the energetic minimum and the direction particles go to in a dynamical system. However, as we can conclude from the density analysis, density is still the driving factor.

<img src="https://i.imgur.com/RewOoli.png
" alt="drawing" width="400"/>

#### KA
<img src="https://i.imgur.com/Kb8ToEM.png" alt="drawing" width="300" /><img src="https://i.imgur.com/NO78d18.png
" alt="drawing" width="350"/>

## MC restriction vs MC voronoi vs
*Found in folder glass/Analysis/code_Monte_Carlo/correlation_quench_MCrest_MCvoronoi on HERA*

If we use a restriction of 1.25 times the particle radius, we find dat Voronoi and radius restriction lead to almost exactly the same results. Note by the way that these are not the restriction radii that lead to the best correlations.

#### HS
:warning: (old quench coordinates for HS, the correct pictures are in the paper)
<img src="https://i.imgur.com/JYSQCeH.png
" alt="drawing" height="200" /><img src="https://i.imgur.com/IrgwY4b.png
" alt="drawing" height="200" />
#### Har
<img src="https://i.imgur.com/QnaWbv9.png
" alt="drawing" height="190"/><img src="https://i.imgur.com/b8k8iBf.png" alt="drawing" height="190" />
#### KA
<img src="https://i.imgur.com/G0YPM6U.png
" alt="drawing" height="200"/><img src="https://i.imgur.com/0jKNWpU.png
" alt="drawing" height="200" />
We found that for A particles the correlation with the voronoi cell is the highest when all particles are restricted to a cell of 1 times. While for B particles it works the best when all particles are restricted to a cell of 1.25 times the particle radius.

**Question: Why would this be the case?**

**QUESTION: why would restricting the particles to a radius of 1.5 lead to the highest correlation? Does that have to do something with how exactly they relax?**
<img src="https://i.imgur.com/einHvKX.png
" alt="drawing" height="300" />
*Figure that shows the MC restriction with all particles restricted to 1 times the radius (left), all particles to 1.25 radius (middle) and voronoi cell (right). Note that in the triangles of the voronoi diagram no particles are allowed.*
### # of MC cycles needed
In order to make the simulations as short as possible, we want to know how long we need to run the system to get a good enough accuracy. To do so, we do simulation for a different number of initialization cycles ($I$), measurement cycles ($M$) and write intervals ($\Delta W$).

To have a good comparison, we perform the run for $I = 5\cdot 10^5$, $M= 10^6$, $\Delta W= 10^2$ twice and look how much they differ. Afterwards we then compare this difference to the difference between runs with less cycles and the run with $I = 5\cdot 10^5$, $M= 10^6$, $\Delta W= 10^2$. Below we show the results. The first row computes the average distances between the same particle in the run with less cycles and the reference run $\Delta r$. For the second row we first compute the difference in distance between the initial snapshot and the position in the MC snapshot in the reference run $\Delta r_{MC}$ and then computes the percentual difference $\Delta r/\Delta r_{MC}$.

| | HS: $I = 10^5$, $M= 5\cdot10^5$, $\Delta W= 50$ | HS: $I = 5\cdot 10^5$, $M= 10^6$, $\Delta W= 10^2$ |Har: $I = 10^5$, $M= 5\cdot10^5$, $\Delta W= 50$ |Har: $I = 5\cdot 10^5$, $M= 10^6$, $\Delta W= 10^2$ |KA: $I = 10^5$, $M= 5\cdot10^5$, $\Delta W= 50$ |KA: $I = 5\cdot 10^5$, $M= 10^6$, $\Delta W= 10^2$ |
| --------------------------------------------------- | ----------------------------------------------- | ------------------------------------------------- | ----------------------------------------------- | ------------------------------------------------- | --- | --- |
| Average difference per particle $\Delta r/\sigma_A$ | 0.0130657| 0.0104697|0.00987157|0.00934684|0.00731098|0.00555999|
| Percentual difference$\Delta r/\Delta r_{MC}$ | 0.128433|0.102915|0.077159|0.0730575|0.0687608|0.0522925|

N.B. Note that the results for $I = 5\cdot 10^5$, $M= 10^6$, $\Delta W= 10^2$ are computed by measuring two different runs, both with $I = 5\cdot 10^5$, $M= 10^6$, $\Delta W= 10^2$ and then comparing the difference.

From this we see that we probably get away with doing $6\cdot10^5$ cycles instead of $15\cdot10^6$, which already speads up the simulation by a factor of 3.
## BOPS

### Parameter settings
We use (in agreement with *Averaging Local Structure to Predict the Dynamic Propensity in Supercooled Liquids - E. Boatini*) the following parameters.

Radial functions are up to 5 shells of particles, angular functions up to 2. See analysis below
(*Found in the folder glass/Analysis/code_paircorrelation/*).

| | $g(r)$ (zoomed in) | second minimum | fifth minimum |
| --- | ------ | -------------- | ------------- |
| HS | | $2.1\sigma_A$| $4.5\sigma_A$ |
| har | |$2.7\sigma_A$ | $5.6\sigma_A$ |
| KA | |$2.3\sigma_A$ | $4.9\sigma_A$ |

* HS (306 in total)
* 162 radial functions with $s\in{A,B}$:
* 46 equally spaced in the interval $r/\sigma_A\in(0.85,2.0]$ with $\delta=0.025$.
* 20 equally spaced in the interval $r/\sigma_A\in(2.0,3.0]$ with $\delta=0.050$.
* 15 equally spaced in the interval $r/\sigma_A\in(3.0,4.5]$ with $\delta=0.100$.
* 144 Angular fuctions
* 12 equally spaced in the interval $r/\sigma_A\in[1.0, 2.1]$ with $\delta=0.1$ and $l\in[1,12]$
* harmonic (428 in total)
* 212 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$.
* 26 equally spaced in the interval $r/\sigma_A\in(3.0,5.6]$ with $\delta=0.100$.
* 216 Angular fuctions
* 18 equally spaced in the interval $r/\sigma_A\in[1.0, 2.7]$ with $\delta=0.1$ and $l\in[1,12]$
* KA (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]$

These are then averaged using the following formula.
$$
\mathbf{X}_i^{(n)} = \frac{1}{C}\sum_{j:r_{ij}<r_c}e^{-r_{ij}/r_c}\mathbf{X}_j^{(n-1)}
$$
We use the following cutoffs (placed on the secon minimum)
* HS*: $r_c/\sigma_A=2.1$
* KA*: $r_c/\sigma_A=2.3$
* HA**: $r_c/\sigma_A=2.7$

### CHECKS
#### HS
*Found in folder glass/Analaysis/code_BOPS/Checks/HS on HERA*

We rewrote the earlier code, such that it can easily be applied to multiple systems with different parameters. To check whether the new code works we check it with earlier results found for the HS system.
<img src="https://i.imgur.com/aeHajGR.png" alt="centered image" width="400" />

In the end we use slightly different cutoff values. We thus look whether for these new values the parameters also look sensible.
*Found in folder glass/Analaysis/code_BOPS/Checks/HS_newnumberparameters on HERA*
<img src="https://i.imgur.com/O1wGe9I.png" alt="drawing" width="350"/><img src="https://i.imgur.com/sPQi89c.png" alt="drawing" width="350"/>
#### harmonic
*Found in folder glass/Analaysis/code_BOPS/Checks/har on HERA*
We check harmonic by looking whether the found parameters look sensible. Below we show the radial distribution of A particles and the BOPS, both for particle 1 in run 1 and dataset 1 of the KA system. As we can see, the results look sensible.
<img src="https://i.imgur.com/0snfeo5.png" alt="drawing" width="350"/><img src="https://i.imgur.com/xy5s7pZ.png" alt="drawing" width="350"/>

#### KA
*Found in folder glass/Analaysis/code_BOPS/Checks/KA on HERA*
We check KA by looking whether the found parameters look sensible. Below we show the radial distribution of A particles and the BOPS, both for particle 1 in run 1 and dataset 1 of the KA system. As we can see, the results look sensible.

<img src="https://i.imgur.com/OAnu6vr.png" alt="drawing" width="350"/><img src="https://i.imgur.com/TT65IM9.png" alt="drawing" width="350"/>

*Averaging Local Structure to Predict the Dynamic Propensity in Supercooled Liquids - E. Boatini

## Pair correlation function different systems

We look at the $g_A(r)$ and $g_B(r)$ of the different systems we consider.
*Found on HERA in the fodler glass/Analysis/code_density/paircorrelation_particlesize/*
### HS
<img src="https://i.imgur.com/T4oU8ba.png" alt="drawing" width="350"/><img src="https://i.imgur.com/n0kTTEZ.png" alt="drawing" width="350"/>

### Harmonic
No results yet
### KA

From the KA analysis we can see that on average large particles have a diameter of 1.1 and small particles have a diamter of 0.9, i.e. for the density analysis we set $r_A/\sigma_A \approx 0.55$ and $r_B/\sigma_A \approx 0.45$
<img src="https://i.imgur.com/1ckMuFa.png" alt="drawing" width="350"/><img src="https://i.imgur.com/IKbsWiK.png" alt="drawing" width="350"/>

As we can see, quench and MC lead to narrower peaks. This is even more the case for MC than for quench: i.e. MC just pushes particles as far away from each other as possible.

## Correlation weighted density and propensity
*Found on HERA in the folder glass/Analysis/code_density/Correlation_density_propensity/*

We found that with the obtained cage centers we can also improved the propensity prediction at the long term. Not only by evaluating the structural parameters, but also by looking at the weighted packing fraction. This parameter is defined as
$$
\rho_i^\text{weighted}=\sum_{i\neq j}\left(\frac{\sigma_j}{2}\right)^3\exp[-(r^\text{MC}_{ij}-\sigma_i/2)]
$$
### HS
Below we show the correlation between this weighted density of both initial, quench and MC coordinates for the *A* particles in a HS system, where the density is computed with in a radius of 4.5 $\sigma_A$.

<img src="https://i.imgur.com/AwJzbWi.png" alt="drawing" width="500"/>

As we can see the correlation is very high, especiallly for times close to $\tau_\alpha$. Since the absolute distance between initial and quench/MC coordinates is a very good predictor for the propensity in the caging regime, this means that we can now do a pretty decent prediction of the propensity in both the caging and the diffusive regime with only two coordinates (instead of the normally used $2000$-ish).

Note that this correlation is only very high when we take into account particles from far away. If we take into account different shells to compute the correlation we see that the correlation increases when we include more shells
<img src="https://i.imgur.com/rSIPaHV.png" alt="drawing" width="500"/>

Moreover the high correlation around $\tau_\alpha$ also tells us that the mobility of particles is mainly density driven. However, to compute this density you need to get rid of thermal noise (which is what the quench and MC method do).

We also try to include higher order generations of density. This mainly increases the correlation at times in between the caging and the diffusive regime, however the increase is not very large

<img src="https://i.imgur.com/QjV0tNg.png" alt="drawing" width="500"/>

Finally we compute the correlation between propensity and both the distance between initial snapshot and MC coordinates, as well as the MC density.

<img src="https://i.imgur.com/2TXB5vv.png" alt="drawing" width="500"/>

### Har
cutoff is set to half the boxsize
<img src="https://i.imgur.com/daNMlmV.png" alt="drawing" width="500"/>

Looking at the result of the harmonic system, we see that appparently even though the quench does not really find the exact caging spot, there is still a lot of information in the density of the quenched system. This is quite unexpected.

### KA
<img src="https://i.imgur.com/TWIkDRV.png" alt="drawing" width="500"/>

We obtain the particle size of these particles using the pair correlation function: $r_A\sigma_A \approx 0.55$ and $r_B\sigma_A \approx 0.45$.

Note that the difference between init, mc and quenchd is very small (compared to HS and MC).
However this difference becomes larger if you decrease the size of the *B* particles to $r_B/\sigma_A$.
**QUESTION: WHY??**

<img src="https://i.imgur.com/vsobUyu.png" alt="drawing" width="500"/>

## Data structure
The different potentials will have the following data structure:
* <potentialname>_<dataset>, wheredataset is 1 or 2
* run_<runname>, where runname $\in [1001,1050]$.
* snap.sph (contains initial snapshot)
* snap_quenched.sph (contains initial positions in the quenched system)
* snap_MC.sph (contains average positions in a restricted MC simulation)
* snap_MC_quenched.sph (contains quenched positions of particles in the MC snapshot)
* density_init.dat (contains weighted density of particles in the initial snapshot)
* density_quenched.dat (contains weighted density of particles in the initial snapshot)
* density_MC.dat (contains weighted density of particles in the quenched snapshot)
* density_init.dat (contains weighted density of particles in the MC snapshot)
* parameters.dat (contains Bond Order parameters for each particle in the initial snapshot)
* parameters_quenched.dat (contains Bond Order parameters for each particle in the quenched snapshot)
* (contains Bond Order parameters for each particle in the MC snapshot)
* meandistdata.dat (contains propensity at logmarithmic time intervals)
* writetimes.dat (contains timestepsin terms of $t/\tau$).

| | HS |Har | KA |
| -------- | ---- | ---- | -------- |
|snap.sph| :heavy_check_mark: |:heavy_check_mark: |:heavy_check_mark:|
|snap_quenched.sph|:heavy_check_mark:| :heavy_check_mark:|:heavy_check_mark:|
|snap_MC_vor.sph|:heavy_check_mark: |:heavy_check_mark:| :heavy_check_mark: |
|parameters.dat |:heavy_check_mark:|:heavy_check_mark:| :heavy_check_mark: |
|parameters_quenched.dat|:heavy_check_mark:|:heavy_check_mark:| :heavy_check_mark:|
|parameters_MC_vor.dat |:heavy_check_mark:|:heavy_check_mark:|:heavy_check_mark: |
|meandistdata.dat|:heavy_check_mark:| :heavy_check_mark:|:heavy_check_mark:|
|writetimes.dat| :heavy_check_mark:|:heavy_check_mark:|:heavy_check_mark: |

## No Linear regression analysis
On *Odin /home/users/5670772/glass/correlation_initial_CS_IS/* one can find the correlation between the distance between IS and CS and the propensity where no LR is used.
## Linear regression analysis
### Checks
*Found in glass/Linear_regression/CHECKS/ on HERA*
We check our new LR code, by comparing the new found correlations for third generation order parameters (based on the initial snapshot), match with the correlations found with the earlier code.

<img src="https://i.imgur.com/Ububk0c.png" alt="drawing" width="500"/>

We do the following input analysis. Note that we always input three generations of structural parameters.
| | HS | KA | Har |
| -------- | ---- | ---- | -------- |
| Par. | | | |
| Par. quenched | | | |
| Par. MC | | | |
| $\Delta r_\text{quenched}$ | | | |
|$\Delta r_\text{MC}$ | | | |
| Par., Par. quenched, $\Delta r_\text{quenched}$| | | |
| Par., Par. MC, $\Delta r_\text{MC}$| | | |

## Comparison KA results with other papers
When comparing the earlier obtained results from other papers (Bapst (2020), Boattini (2021), Shiba (2022), Landes (2022), Jung (2023))with our results we run into a few problems.

First of all the time units. None of the other papers specify their time units. However, since the most intuitive choice while working with a KA system seems to be lennard-jones units I assume they measure time in units of $\sqrt{m\sigma_A^2/\epsilon}$. Since we measure time in $\sqrt{m\sigma_A^2/k_BT}$, we have to multiply their points in time by $\sqrt{k_BT/\epsilon}$ in order to be able to compare them with our results (we ).

>[name=Frank] Actually, the Bapst paper does specify their units: "On the basis of the Lennard–Jones potential, we define dimensionless units and measure distances in units of σAA, time in units of $\tau = \sqrt{m\sigma_{AA}^2/\epsilon_{AA}}$ (with m the atomic mass), ..." (in the Methods section). So that agrees with your assumption. Also the Jung paper does specify that they work in LJ units, and specify the time unit as well (again agreeing with your assumption).
>[name=Rinske] Oops sorry, then I didn't look carefulll enough. But then at least we know multiplyling with $\sqrt{k_BT/\epsilon}$ is the correct thing to do.

Below I will discuss per paper the choices I made and the problems I encouter.

* Jung (2023). In Figure 5 they compare the results of the GNN of Bapst with their MLP. The problem is that their GNN results do not match with the results of Bapst (Bapst considers slightly different points in time and has slightly higher correlations). It could be that Jung *et al* wrote their own GNN code, however I'm not sure about this.
>[name=Frank] My guess would be that they were slightly sloppy in extracting the data from the Bapst paper.
>[name=Rinske] But they even have an extra datapoint... That to me seems kind of weird.

* Landes (2022) uses the same weird time convention as Bapst (i.e. they consider stepwise intervals of $S(q,t)$). Since their time axis is unreadable. What I do now is use the times of Bapst paper (from that paper we have the actual data) and match them with the correlations found by Landes.
>[name=Frank] I don't really understand what you mean by the time axis being unreadable? Each point seems clearly labeled in Fig. 5? But the data is the data from Bapst, so it should agree.
>[name=Rinske] Sorry, maybe I wasn't clear enough. What I meant was that Landes does not use a linear or a log time scale. And therefore it is hard to read of the times (we had the same problem with the data from Bapst.)
>[name=Laura] The times are expliitly written on the x-axis-- you do not need to extract them since they are marked independently. You only need to extract the y-coordinates.

* My data does not match the data of Emanuele. This difference could maybe explained by the fact that we use different datasets (as far as I know Emanuele uses the dataset of Bapst that consists 400 configurations with 30 isoconfigurational ensemble runs, while we have 100 configurations with 50 isoconfigurational ensemble runs). Moreover we have a slightly different set of parameters.
>[name=Frank] Those both seem good explanations, and should be mentioned in the discussion.


## To do later
- [ ] Send references freezing paper to Laura en Frank
- [ ] Start new project on
- [ ] Compute correlation cage centers unfrozen and frozen for larger system

## Removed from paper
These sections where removed from the paper, but they might be useful for later reference:
- Comparing the correlations found by Bapst *et al* with our obtained correlations, we see that including the CS especially improves the predictions during the caging regime. Although the CS and even the IS are able to probe the positions of the cage centers, apparently the GNN is not able to. The same results was found by Alkemade *et al*, who showed in Ref. that also in the case of Hard-spheres the GNN has hard time predicting the propensities during the caging regime. We believe that there are at least two reasons why even complex methods like GNN's have a hard time dissecting the positions of the cage centers from the initial positions. % First of all, the caging process is a process in which directional information matters: On average the movement of particles in this regime can be described by a vector pointing from the initial position of particles to the associated cage center. However the GNN of Bapst *et al* and Alkemade *et al*, are implement such that they are rotationally invariant, meaning that all information about directionality is lost. It is therefore not surprising that these methods have a hard time dissecting quantities like $\Delta r_\text{IS/CS}$, which are quantities that are closely related to vectors. Note that in the recent paper of Pezzicoli *et al*, the authors have developed a equivariant GNN that is able to deal with vectors. Comparing their obtained correlations with the correlations of Bapst *et al*, we see that using directional information indeed improves the propensity predictions.

The second reason that even more complex machine learning methods will have a hard time dissecting the positions of the cage centers, is because it turns out that finding the location of the cage centers is a highly collective process. To show this, we compute the average position of particles in their cage while part of the system is frozen at its initial position. Concretely, this means that for each particle $i$ in the system, we freeze all particles at a distance further than $r_f$. Thereafter, we let the particles that have $|\mathbf{r}_i-\mathbf{r}_j|<r_f$ move around in their cage, while measuring the average position of particle $i$. We perform this simulation for various values of $r_f$, after which we compute the correlation between the found ensemble average positions for each particle and the associated propensity. In Figure \ref{fig:partiallyfrozen} we plot the results for a system of 2000 particles. As we can see, for lower values of $r_f$ the correlation peak is significantly lower and shifted to the left. Keeping in mind that the box to which the particles are restricted has a boxlength of approximately $5.5\sigma_A$, we see that to find correlations that match with the unfrozen system, practically the entire box needs to be allowed to move. Finding the positions of the cage centers thur truly is a collective phenomenon. Again, this explains, why even the more complex machine learning methods have a hard time making good propensity predictions during the caging regime.
## Dates
### Lammps simulations
* **2022-08-30** NVT simulations for both KA and Harmonic are started (100 snapshots for each system).
* **2022-09-01** NVE equlibration simulations for KA are started.
* **2022-09-07** NVE equlibration simulations for har are started.
* **2022-09-07** Measurements runs for KA are started with a crontab job on the cluster
* **2022-10-10** Measurements runs for harmonic are started with a crontab job on the cluster