# Voronoi model

Github link: https://github.com/Syrocco/Voronoi

Forces derived from:
$E = \sum_i K_A(A_i-A_0)^2+K_p (P_i -P_0)^2$

Periodic boundary conditions

### Equilibrium

Some test at small temperature (overdamped with kind of a langevin bath). $dt = 0.1$, $A_0=1$, $K_A=1$ and $K_p=1$.

| what am i looking at | $p_0=3.6$ | $p_0=4.2$ |
| --- | --------------------------------------------------- | --------------------------------------------------- |
| Voronoi diagram |  $~$| $~~~~~~~~~~~~~~~~~~~~~~~~~$ |
| Position of selected particle (red) with time |  |  |

### Sheared

$p_0=3.6$, overdamped without langevin bath

| $\gamma$| $\sim 1/4$ | $\sim 1/2$ | $\sim 1$ |
| --- | ---------- | -------- | -------- |
| Usual time |  |  |  |
| Strobo | |||

$p_0=4.1$

| $\gamma$ | $\sim 1/4$ |
| -------- | -------- |
| Usual time |  |
|Strobo| |

Appears irreversible, probably still irreversible with dt and $\dot\gamma$ small enough.

________________________________________

What about bidisperse?

$E = \sum_i K_A(A_i-A_{0, i})^2+K_p (P_i - q_0\sqrt{A_{0, i}})^2$

Bidisperse ratio of area = 4/3, $q_0=3.65$ leads to a frustrated/glassy state.

Unfortunately it does not seem to help to reach an irreversible state:

| | Monodisperse | Bidisperse |
| --------------------------------------- | ---------------------------------- | --------------------------------------------------- |
| Unsheared |  |  |
| Sheared $\gamma=1$: absorbing state |  |$~~~~~~~~~~~~~~~~~~~~~~~$  |
| Sheared $\gamma = 0.5$: absorbing state | | |

### FIRE:

Now, we minimize the energy using FIRE.

| | Monodisperse | Bidisperse |
| --------------------------------------- | ---------------------------------- | --------------------------------------------------- |
| Sheared $\gamma=1$: absorbing state |  |$~~~~~~~~~~~~~~~~~~~~~~~$  |
| Sheared $\gamma = 0.5$: absorbing state | | |

____________________

| large $\gamma_{max}$ | small $\gamma_{max}$ |
| --------------------------------------------------- | -------------------- |
|  | 

#### Problem with Fire:

In the solid region at $f_{norm}~\sim 10^{-12}$, $f_{norm}$ does not decrease anymore, smaller dt does not seem to help:

$q_0=3.65$ (forces = green, velocities = orange)

In the liquid, at $f_{norm}~\sim 10^{-3}$, $f_{norm}$ does not decrease anymore, smaller dt does not seem to help:

$q_0=3.9$
Problem seems to be around T1: zoom 

Note that an overdamped dynamics is able to relax it below what FIRE can do and seem to always decrease the energy:

| Evolution energy| Evolution force |
| -------- | -------- |
|  |  |

#### Comparison simple shear:

$K_A=0$ monodisperse, $f_{norm}=10^{-7}$ (seems independent on $f_{norm}$):


Not the same...

Close to the solid same initial condition but different $\gamma$ step:

|  |  |  |
| -------- | -------- | -------- |

In between

|  |  |  |
| -------- | -------- | -------- |

Close to the liquid:

#### Choice of fnorm:

$\Delta \gamma = 0.007075$

|  |  |  |  |
| -------- | -------- | --- | -------- |

$\Delta\gamma = 0.001$ not really helping


#### $K_a=0$ vs $K_a=1$

| $K_a=0$ | $K_a = 1$|
| -------- | -------- |
|  | |

#### Problem with intermediate dt:

Here is a pure gradient descent with a small dt step:

| $dt=0.01$ | $dt=0.001$|
| -------- | -------- |
|  | |

#### FIRE vs Conjugate gradient descent.

As FIRE cannot resolve nicely the T1, we tried Conjugate Gradient descent. Which is just a fancy gradient descent without inertia.

CG:
- Starting from a given configuration, pick a direction for which the energy decreases (Hence: close to the gradient direction)
- Search an approximation for this minimum along this direction 
- Iterate

One problem is that in the search of an approximation for the minimum along the choosen direction, we can overshoot the local minimum...

Anyway.

FIRE does not manage to minimize correctly the energy around T1 except by pushing the particles close to a T1 using a large $dt$, this is problematic because it is arbitrary and most importantly, T1 should be induced by the modification of the potential landscape and not by the algorithm.

CG manages to minimize the energy up to a given point in which the configuration always end up with particles that would like to do a T1:


The forces obtained from the theory here, are wrong as the gradient descent performed along this force never decrease the energy. This is due to the fact that the forces computed analytically do not take into account the change of topology of the graph during a T1.

#### Comparison between CG and FIRE
Binary mixture

##### $q_0=3.7$

Note that for FIRE, all large drops seem to correspond to the unphysical resolution of T1:

##### $q_0=3.85$


##### Quick simulations for absorbing state in bidisperse ($K_A=1$):

FIRE (crashes in the active region):


Largest $\gamma_{max}$ with an absorbing state:

| $\Delta\gamma = 0.001$ | $\Delta\gamma = 0.01$ |
| -------- | -------- |
|  |   |

_________________________________________

CONJUGATE GRADIENT:

First active state:
| Strobo | Steady state |
| -------- | -------- |
| |  |

Last absorbing state:

| Strobo | Steady state |
| -------- | -------- |
|  |  |

The transition is a trivial one, we lose reversibility when it yields. It is plausible that the irreversibility we observe comes from the problems of the algorithm and not of the system itself.

### Testing more parameters

Conjugate Gradient:


#### Last absorbing state:
shear vs strain:

stroboscopic energy evolution:

#### First active state:
shear vs strain:
stroboscopic energy evolution:

_____________________________________

FIRE:

#### Last absorbing state:
shear vs strain:
stroboscopic energy evolution:

#### First active state:
shear vs strain:
stroboscopic energy evolution:

#### What the hell is wrong with the energy?

I'll be shoing configurations in which CG does not manage to find a minimum of energy. This happens around $f_{norm}\sim 10^{-7}$, the difference in energy between neighboring points will be around $10^{-14}$. This is however not the case at large $q_0$ for which, even at $f_{norm}\sim 10^{-2}$ we don't converge anymore. This is really due to a T1.

Monodisperse $K_A = 1$:
| $q_0=3.$ | $q_0=3.65$ | $q_0=4.$ |
| ------- | ---------- | ------- |
| | | WTF|
| || $f_{norm}\sim 10^{-4}$|

Monodisperse $K_A = 0$:
| $q_0=3.65$ | $q_0=3.8$ | $q_0=4.$ |
| --------------------------------------------------- | --- | --------------------------------------------------------------------------- |
|  |  | converges |
|  | |  $f_{norm}\sim 10^{-10}$ |

Bidisperse:

| $q_0=3.$ | $q_0=3.65$ | $q_0=4.$ |
| ------- | ---------- | ------- |
| |  |  |
|  | |  |

for $q_0=4$, the forces computed theoretically are the right one as I verified them numerically by finite difference.

### 64-bits, 80-bits and 128-bits

The problem lies into the machine precision. Here, i'm running different minimization and looking at the energy landscape in the direction of the force gradient when $f_{norm}\simeq 5\times 10^{-7}$. Comparing the exact same system necessitated to be very precise with the position saving, i was lazy.

Don't compare the energy, I had to do some shenanigans to plot them for 80b and 128b systems.

|precision: | double (hardware accelerated with SIMD, 64-bits) | long double (hardware accelerated, 80-bits) | quad (software emulated, 128-bits) |
|-------| -------- | -------- | -------- |
|energy landscape| ||  |
|$\epsilon$ precision|$2\times10^{-16}$|$1\times 10^{-19}$|$2\times 10^{-34}$|
|Minimal $f_{norm}$ we can obtain with energy based CG|$f_{norm}^{min}\simeq 5\times 10^{-7}$|$f_{norm}^{min}\simeq 10^{-9\sim 10}$|$f_{norm}^{min}\simeq 10^{-15\sim?}$|
|relative slowdown with respect to double precision|$1\times$slower|$1.5\times$slower|$5.5\times$slower|

Long double is hardware accelerated by x87 instructions. SIMD not available, this is getting optimized, most loops are too:
Hot parts are linked list backed by arrays and cannot be SIMDed so the perf gain is not too dramatic

### Problem with minimizer

[Mirages in the Energy Landscape of Soft Sphere Packings](https://arxiv.org/pdf/2409.12113). Fast minimizers are not accuratly staying in the same local minima.

#### Fire:
##### PRO:
- Fast
- Force based, machine precision not too problematic
##### CONS:
- Possible use of large $dt$ to make convergence fast can make the system fly over minima.
- Convergence problem around T1s when particles are close.
- Large $dt$ (or maybe arbitrarly small ones?), can induce fake T1s when particles are very close to each others. Said differently, I think a T1 should be resolved at the first initial steps of the minimizations (due to an energy landscape change), not in the middle of the minimzation where a T1 can probbly only be induced by a (problematic) potential barrier jump.

____________________

#### Conjugate Gradient:
Select a direction in the energy landscape (usually close to the gradient descent direction) and perform a line search along this direction.
##### Energy based minimizations
BACKWARD SEARCH:
##### Pros:
- Less problem than FIRE to converge (deal with particles close to each others a bit bette)
##### Cons:
- 
- ~$5$ times slower than FIRE.
- Double precision is not enough for accurate energy minimizations (factor: $^2$ compared to force based minimizations)
____________________________________________________________

FORWARD SEARCH:

Qualitatively different than backjward search (compare purple (forward) vs green and red (backward)):

##### Pro
- Less problem than FIRE to converge (deal with particles close to each others a bit bette)
- Converge a bit faster than backward
- Probably jump less over minima than backward search

##### Cons
- Double not enough for energy minimizations

#### Force based minimizations:
Based on a Fortran routine in DYNAMO, [reimplemented in LAMMPS ](https://github.com/lammps/lammps/blob/develop/src/min_linesearch.cpp)
Here, we try to reduce the FORCE, not the energy.
##### Pro:
- Better convergence without double
##### Cons:
- Horribly slow (have to compute force, which are costly compared to energy)

### Single shear comparison energy minimization method

RKF45:

We solve the dynamics:

$$\Gamma d r_i/dt = F_i ~~~~\text{ along }~~~~ d\gamma/dt = \dot{\gamma}$$

by integrating these equations and advancing time with a small $dt$.
The solver is adaptative, we solve at order 4 and check accuracy at order 5. If order 4 - order 5 is above some small threshold, we take a smaller dt
##### Bidisperse, $q_0=3.5$


Individual realization:

| Fire dt=0.1 | fire dt=0.01 | Conjugate gradient | RK45 adaptative |
| -------- | -------- | --- | -------- |
|| ||  (don't mind the weird zoom) |
|||||

With RK45, we can see what is the optimal dt for a given error, here with an error between the 4th and 5th order term below $10^{-6}$, we obtain:

##### Monodisperse, $q_0=3.9$ and $K_A=0$


Individual realization:

| Fire dt=0.1 | fire dt=0.01 | RK45 $\dot\gamma = 0.0003$ |RK45 $\dot\gamma = 0.003$ |
| --------------------------------------------------- | --------------------------------------------------- | --------------------------------------------------- | --- |
|  |  | ||
|  |  |  | |

Concerning the optimal dt with rk4:

| Column 1 | Column 2 |
| -------- | -------- |
|  |  |

singular force to perform a T1:

| snapshot | zoom | stress | dt |
| ---------------------------------- | ---------------------------------- | --- | --- |
|  |  |||

Local rearrangment after the singular force is gotten rid off by the shear:

(individual images here: https://drive.google.com/drive/folders/14r4J-b-pgYFYbe-Y9EuGzi7kLaRzLFkB)
They arange as a long defect:

## Adding soft cores

Runge kutta soft cores:


Still the same problem. No overlaping anymore but still ill defined lattice.