# Granular Project V1.0

###### tags: `Granular Simulations` `SoftQC` `self-assembly`

---

---

### Owners (the only one with the permission to edit the main test)

Etienne, Andrea Jr., Andrea Sr., Frank, Giuseppe

---

---

## Background

The idea is to try to see if the *thermal* crystal structure discussed by Etienne in:

[https://arxiv.org/pdf/2003.08889.pdf](https://arxiv.org/pdf/2003.08889.pdf)

can be observed in 2D vibrofludized granular mixtures.
To this end, Andrea Jr and Andrea Sr will help us adapting their 3D model as discussed here:

[https://arxiv.org/pdf/2003.08662.pdf](https://arxiv.org/pdf/2003.08662.pdf)

to 2D vibrating plane with sphere of different size deposited on the top.

## Plans

Try to understand if we can observe in granular systems the same crystal that we see in thermal systems (colloids).

## Notations

* $\sigma_i$: diameter of disks of type $i$. $i=L, S$ for large and small disks respectively.
* $q=\sigma_S/\sigma_L$: size ratio of the particles.
* $N_i$: number of disks of type $i$.
* $x_S=N_S/(N_S+N_L)$: number fraction of small disks.
* $\Delta$: non-additivity parameter. It is defined by $\sigma_{LS}=(\sigma_L+\sigma_S)/2 \times (1+\Delta)$.
* $V$: area of the 2D box (we often call it *volume* even though it is technically an *area*).
* $\rho=(N_L+N_S)/V$: particles density.
* $\phi=\pi (N_S \sigma_S^2+N_L \sigma_L^2)/4V$: *packing fraction*, surface fraction occupied by the disks.

## Sample Preparation

We list here available starting configurations for granular simulations. There are three types of configurations:
* very dense jammed initial state configurations, that need to be expanded to the desired density
* crystallised configurations that are produced by MD simulations. They come at the density at which they self-assembled (not jammed)
* ideal S1 crystal.

In each subsection, we provide a link to download the correponding data-file. (One in LAMMPS-like format, and one with 4 columns: particle ID, particle type, xs, ys). Should you need another format, please ask !

Data-files contain **scaled coordinates**, *ie.* the real coordinates divided by the length of the box.
Box parameters can be found in the LAMMPS-style files.

### Initial state configurations

These configurations are obtained from a rapid compression of the mixtures in MC simulations. They are jammed and of very high density. We decrease their density to the desired one, and use them as starting configurations for simulations. To reach the desired density, we simply choose the corresponding box size (scaled coordinates don't change).

* $N_L=500, N_S=500, q=0.457, \Delta=-0.07, \rho=1.9073668228800094, \phi=0.9054535283662086$ [MC_q0.457_delta0.93_N0500_N1500_aspect1](https://catdrop.drycat.fr/r/5uYHkkh_#GoM6roMYLqM+PFHJi2OucJ4E8qjgyLyJctLhcbKPBfU=)
* $N_L=683, N_S=317, q=0.457, \Delta=-0.07, \rho=1.4938327116499532, \phi=0.8790075206502072$ [MC_q0.457_delta0.93_N0683_N1317_aspect1](https://catdrop.drycat.fr/r/HflBdkCX#g1vqvpbXsZMjFekSx6AeE6BLq69F3Kjw2TNazT5nYSo=)

### Crystallised configurations

* S1 $N_L=500, N_S=500, q=0.457, \Delta=-0.07, \rho=1.9073668228800094, \phi=0.855$ [EDMD_q0.457_delta0.93_N0500_N1500_packfrac0.855_aspect1](https://drop.infini.fr/r/dAAEMFT0G2#g1LnYuGqcncGEvrkckzQyzTmE9p2XzmlhCYljObBJnQ=)

* Random tiling $N_L=687, N_S=317, q=0.457, \Delta=-0.07, \rho=1.41904356463743436923, \phi=0.835$ [EDMD_q0.457_delta0.93_N0687_N1317_packfrac0.835_aspect1](https://drop.infini.fr/r/mbV0Lu1a5Q#S4KrfbQ42FinZBCtwn8/7H2OoYyXDiQpVX8dHecA414=)
>[name=Etienne] It seems to be harder to equilibrate a random tiling in a small system. This one is far from perfect, even though it exhibits the 12-fold structure factor.

### Ideal S1 crystal
Just an ideal S1 crystal with 1000 particles. Can be inflated as well to reach a specific density.

* $N_L=500, N_S=500, q=0.457, \Delta=-0.07, \rho=1.89035916824196597353, \phi=0.8973797584355784$ [idealS1_q0.457_delta0.93_N0500_N1500_aspect1](https://catdrop.drycat.fr/r/UJoRGGed#NZEp1BlENGc/ilS92zcpEh4a5n+C7GgmZI2812bY098=)

## Granular Simulations

* First Granular simulations starting with both initial and final configurations
* We can play a bit with the density for the kinetics
* The physica parameters will be the same as in Andreas'paper linked above (steel balls on a steel plane)

### Numerical setup

The numerical setup consists on a box of size $L\times L\times h$ with periodic boundaries on the $xy$ sides. $L$ is found by fixing $\sigma_s=0.001$ m, $\sigma_L=\sigma_s/q$: $L=\sqrt{\pi(N_s\sigma_s^2+N_L\sigma_L^2)/4\phi}$ where $\phi$ is the desired packin fraction. With this procedure we can convert all the configuration from MC and EDMD to the the scale of millimiter. In order to realize a pseudo-2D setup, the height of the box is taken as $h=\frac{6}{5}\sigma_L$. The system is subjected to the gravity acceleration $g$. The bottom of the box oscillates with a sinusoidal signal with amplitude $A$ and frequency $f$ $\Delta z(t)=A\sin(2\pi f t)$ while the upper side does not move. We use the Hertz-Mindlin model to describe both the particle-particle interactions and the particle-wall ones (see SM of https://arxiv.org/pdf/1907.02389.pdf). We note that the rotational degrees of freedom of the grains and tangetial forces are crucial for trasmission of energy from the vertical direction to the horizontal one.

### First results

For this preliminary study, we have choosen values for the material properties that have already demostrated to be coherent with a realistic behavior of vibrated granular media. Also the values of $A$ and $f$ are choosen with the same motivation. Up to now we have results from the evolution of the configuration generated by MC with $N_s=N_L=500$ (MCconf) and the crystalized one obtained after EDMD evolution (S1conf). We have fixed $A=2.6\times 10^{-5}$ m so we are in the limit $A\ll\sigma_s$ in which the kind of agitation induced on the grains is more 'thermal'. We have performed simulation for three values of $f=[200,500,800]$ Hz and the relative $\Gamma=(2\pi f)^2A/g$ values are $\Gamma=[4.2,26.2,67.0]$. In all the cases (but with more evidence at high driving frequencies), starting both from MCconf and S1conf, we observe size segregation i.e. clusters of particles with the same size are formed. This phenomenon is different from the brazil-nut effect observed in vertical setup where there is a real size separation in space (small particles on the bottom and large particles on top).
Changing $f$ seems to affect only the velocity with which the clusters are formed. Looking at the following snapshots we can also see that some S1 domains persist also until the end of the simulation signaling that this spatial structure is (in some sense) resistent.

* Starting from **MCconf**. $L=0.1$ m, $\sigma_s = 0.001$ m, $N_s=N_L=500$, $q = 0.457$, $\phi = 0.855$, $f = 200$ Hz, $A=2.6\times 10^{-5}$ m. Snapshots of the system for three different times.
 0 s  3373 s 6740 s

* Starting from **S1conf**. $L=0.1$ m, $\sigma_s = 0.001$ m, $N_s=N_L=500$, $q = 0.457$, $\phi = 0.855$, $f=200$ Hz $A=2.6\times 10^{-5}$ m. Snapshots of the system for three different times.
 0 s  3373 s 6740 s
The system has not yet reached equilibrium, but it should eventually look like the MCconf system that is equilibrated.

* Starting from **S1conf**. $L=0.1$ m, $\sigma_s = 0.001$ m, $N_s=N_L=500$, $q = 0.457$, $\phi = 0.855$, $f=800$ Hz $A=2.6\times 10^{-5}$ m. Snapshots of the initial and final configurations.  0 s  67 s

* Starting from **MCconf**. $L=0.1$ m, $\sigma_s = 0.001$ m, $N_s=N_L=500$, $q = 0.457$, $\phi = 0.855$, $f=200$ Hz $A=1.3\times 10^{-5}$ m. Snapshot of the final configuration.  6740 s
The small grains are apparently not piling up. It seems that with lower A the S1 domains has growth. Check with longer simulations. Check with the orientational bond parameters.

q4 is not clearly increasing, even though we can see grains of S1 forming in the snapshots. Much longer simulations are required (That we can run during the summer on our cluster). This simulation and the one starting from S1 should eventually reach the same state, and both barely evolve. The kinetics seems much slower with a smaller amplitude (which makes sense).

* Starting from **S1conf**. $L=0.1$ m, $\sigma_s = 0.001$ m, $N_s=N_L=500$, $q = 0.457$, $\phi = 0.855$, $f=200$ Hz $A=1.3\times 10^{-5}$ m. Snapshot of the final configuration.  6740 s
The defects don't nucleate but migrate in the crystal that remains stable (see the movie-> https://we.tl/t-AJzwLsZiEW). Check with longer simulations. Check with the orientational bond parameters.
Very little change along the simulation (q4 seems to decrease very slowly). In our simulations, we usually see a S1 crystal with a pocket of fluid moving around. It is encouraging to see the same behaviour here, although the moving liquid pocket seem to disrupt a bit the crystal.

* Starting from **S1conf**. $L=0.1$ m, $\sigma_s = 0.001$ m, $N_s=N_L=500$, $q = 0.457$, $\phi = 0.855$, $f=250$ Hz $A=1.3\times 10^{-5}$ m.
S1 (meta)stable

* Starting from **MCconf**. $L=0.1$ m, $\sigma_s = 0.001$ m, $N_s=N_L=500$, $q = 0.457$, $\phi = 0.855$, $f=250$ Hz $A=1.3\times 10^{-5}$ m.
Small increase in q4, let's try running longer !

* Starting from **MCconf**. $L=0.1$ m, $\sigma_s = 0.001$ m, $N_s=N_L=500$, $q = 0.457$, $\phi = 0.855$, $f=200$ Hz $A=1.3\times 10^{-5}$ m.
Slow increase of q4.

* Starting from **S1conf**. $L=0.1$ m, $\sigma_s = 0.001$ m, $N_s=N_L=500$, $q = 0.457$, $\phi = 0.855$, $f=200$ Hz $A=1.3\times 10^{-5}$ m.

## Hunting for S1

### Frequency-Packing fraction scan
Preliminary simulations by Andrea hinted the possibility to self-assemble S1. The dynamics appeared to be extremely slow, so I decided to perform some scans of the main parameters to try and find a region where S1 equilibration could be faster. The results of the scan are presented below in big pictures, first with last snapshot of the simulation, and then with boop analysis. Open the image on a new tab for closer inspection.
All simulations were performed with amplitude $A=1.5\times 10^{-5} m$.
Scanned parameters: $\phi \in [0.820,0.850]$, $fr \in [200 Hz, 600 Hz]$.
For frequencies above 400 Hz, small particles pile up. The only simulation with increasing $q_4$ is $(\phi=0.850,fr=200 Hz)$, but the trend is still weaker than in Andrea's simulation at $A=1.3\times 10^{-5} m$, $\phi=0.855$ and $fr=200 Hz$. It seems indeed that low frequency, low amplitude and large packing fractions are favorable. This is a bit surprising since, intuitively, I would expect slower dynamics in this region of the parameter space.
To do:
* Restart Andrea system for an extra 1e11 steps (on the way).
:::success
* Scan low amplitude-low frequency region of the parameters space (on the way: $A\in [0.5e-5, 1.3e-5]$, $fr \in [100,400]$ at $\phi=0.855$)
:::



{%youtube HZ5tGORziWQ %}

The systems at high packing fraction looks very close to jamming. Small particles remain trapped in their cages for a very long time. We need to speed up the dynamics.
To do:
:::success
* Try changing the timestep. So far, it is fixed at $0.2\times t_R$. Maybe this is unnecessarily small. Run simulations with different timesteps, and monitor $g(r)$ and velocity distribution to see when the simulation goes wild (on the way).
:::
* Compute MSD to assess the speed of the dynamics. The system drifts a lot during simulations, this drift must be substracted before MSD is computed.

### Amplitude - Frequency scan
Fixed parameters: $\phi=0.855$.
Scanned parameters: $A \in [0.5e-5m,1.3e-5m]$, $fr \in [100 Hz, 400 Hz]$.
For low frequencies and amplitudes, nothing moves. Small particles vibrate up and down without affecting the large ones.



The state point $A=1.3e-5\,m$,$fr=200\,Hz$ is the same as the MCconf8 run by Andrea. However, results are significantly different: $q_4$ is not growing. This behaviour is also observed in the restarted run (see next section). Is there a hidden parameter missmatch between our simulations, or are we juste seeing fluctuations ?

### Restarted promising run

I restarted the run (MCconf8) that exhibited a nice increse in $q4$ for an extra $1e11$ timesteps. The simulation is still running, but here are the preliminary results.
Parameters are: $A=1.3e-5\,m$, $fr=200\,Hz$, $\phi=0.855$, $dt=6.74645e-06\,s$.
The $1e10$ first timesteps:

The $3.5e10$ extra timesteps run so far:

$q_4$ drops, as in the simulation at the same state point in the frequency-amplitude scan (see previous section).

### Timestep optimisation

To speed up the simulations, one option is to increase the timestep $dt$.
The timestep is set as the Rayleigh time multiplied by a factor. Previous simulations were run with a factor 0.2, wich corresponds to a timestep of $dt=6.74645e-06\,s$.
The Rayleigh time only depends on material properties so I first tried to increase the timestep factor in dillute systems ($\phi = 0.3$), while monitoring the radial distribution function and the velocity distribution. Plots of the radial distribution functions are offset by 0.1 for clarity.


The radial distribution function is unaffected by the change in dtfactor, so we can assume that the simulations with larger timestep essentially have the same structural properties than the ealier ones.
The velocity distributions differ slightly: as dtfactor increases, the distribution is more peaked. Averaging the velocity distribution over early steps for the faster simulations does not change this behaviour, so this is likely a difference in the actual dynamics rather than the faster simulations being better equilibrated.

For larger dtfactors (4 and larger), particles are lost during the simulation, so there is no doubt that the timestep is too large at this point.
With a dtfactor of 3.5, the simulations are sped up by a factor 17.5 with respect to the ones with dtfactor of 0.2.

In order to check that such a promising acceleration could be achieved with high density systems, I ran simulations at $\phi=0.855$, $A=1.3e-5\,m$ and $fr=200\,Hz$ (same parameters as the restarted run).
Here, particles are lost when dtfactor is larger than 1.1.
To monitor the stability of the simulation, $q_4$ is computed, and plotted against the actual simulations time in seconds.
The simulation with dtfactor of 1.1 seems to behave slightly differently than the others. This suggests that we can safely increase dtfactor to $\sim 1$. This would accelerate the simulations by a factor 5, not as huge as the previous 17.5, but still decent.

### Model parameters 28/09/2020

Andrea's and LPS simulations at the same state point looked quite different, suggesting that there might be more than just fluctuations. I compared values of the model's parameters from 1/Andrea's paper (https://arxiv.org/pdf/1907.02389.pdf p.7), 2/ the script we used for our simulations and 3/ the exact script that Andrea Jr used for his simulations.
Discrepancies are found for some parameters:

|Parameter|Paper's value|Andrea's script|Provided script|
|:---:|:---:|:---:|:---:|
|$\gamma_n^{aa}$|$3\times10^4\,s^{-1}m^{-1}$|$1.08\times10^7\,s^{-1}m^{-1}$|$1.08\times10^7\,s^{-1}m^{-1}$|
|$\gamma_t^{aa}$|$1\times10^4\,s^{-1}m^{-1}$|$3.26\times10^7\,s^{-1}m^{-1}$|$3.26\times10^6\,s^{-1}m^{-1}$|
|$\gamma_n^{ab}$|$3\times10^7\,s^{-1}m^{-1}$|$3.2\times10^6\,s^{-1}m^{-1}$|$3.2\times10^7\,s^{-1}m^{-1}$|
|$\gamma_t^{ab}$|$1\times10^5\,s^{-1}m^{-1}$|$2.6\times10^4\,s^{-1}m^{-1}$|$2.6\times10^5\,s^{-1}m^{-1}$|

The discrepancy in $\gamma's$ is due to a inversion of the values of *Gamma_t_Pref* and *Gamma_t_12_Pref* (they have respective values of 3 and 0.3 in Andrea's script and 0.3 and 3 in ours). This is likely a mistake that occured during the translation of the script from the old granular package to the new one.
The discrepancy with respect to paper's values for $\gamma_t^{aa}$ and $\gamma_t^{ab}$ is more severe.

I ran simulations with corrected *Gamma_t_Pref* to see the inpact on the self assembly.

### "Corrected" parameter

I restarted Andrea's simulation at $\phi=0.855$, $A=1.3e-5\,m$, $fr=200\,Hz$ with dtfactors of 0.2 (same as Andrea), 0.5 and 1.0.
S1 self-assembles consistently !

|  |  |  |
|:------------------------------------:|:------------------------------------:|:------------------------------------:|
| dtfactor = 0.2 | dtfactor=0.5 | dtfactor=1.0 |


Increasing dtfactor seems to speed-up the simulation without affecting the physics.

Performed a new scan of frequency and packing fraction with the corrected paramters:


For frequency of 100 Hz, not enough energy is pumped into the system and nothing moves.
Above 400 Hz (included), small particles pile up. At 500 Hz, the 2D effective density is lowered enough to allow large particles to crystallise into a hexagonal packing.
At 200 Hz, S1 grows above $\phi = 0.845$. The growth rate decreases at high packing fractions because dynamics slow down.
At 300 Hz, the growth is much faster. A few small particles sometime pile up in the fluid pockets. At $\phi=0.855$ no fluid remains and the grown S1 crystal is essentially perfect. System at $\phi=0.860$ follows the same trend, but slower.
For reference , in molecular simulations, we observed the fluid crystal coexistence for packing fractions in $[0.84,0.87]$.

$q=0.457,\,x_S=0.5,\,A=1.3\times10^{-5}m,\,fr=300\,Hz,\, \phi=0.855$. There is only one defect in this crystal.

Large particles velocity distribution (no crystallisation)
Small particles velocity distribution (no crystallisation)

Large particles velocity distribution (crystallisation)
Small particles velocity distribution (crystallisation)

Large particles average velocities (no S1)
Large particles velocities std (no S1)

Large particles average velocities (S1)
Large particles velocities std (S1)

**Plan**
* Discuss with Andrea's to fix material properties before simulating the larger systems (with hard walls, and possibly no top plate) that could be used in experiments.

## Still to do

* Compute MSD to assess diffusion. Substract the drifts.
* Try running simulations with larger q.
* Run extra simulations around the promising point with the larger timestep.
* Restart simulations taking S1 crystal as initial configuration. (Ongoing)
* Try starting with an S1 cluster + fluid.

### OPEN QUESTIONS

* Will the system segregate?
* Does the kinetic allow rearrangments at the chosen density?
* Do we need to play with the physical paramenters?