---
# System prepended metadata

title: Granular Simulation V3.0
tags: [granular, self-assembly, Granular Simulations]

---

# Granular Simulation V3.0

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



---


---


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

AP, EF, FS, GF



---


---


## Background

The summary of the results obtained in the first weeks is contained in the presentation: https://docs.google.com/presentation/d/1R3FBbRSL6oeysVHRWDwHkMgaEMKsXJ1ZfiKB_gNY96s/edit?usp=sharing


## Plans



## To Do
* START STRUCTURE A PAPER WITH FIGURE 
* QC8 realistic setup (Ntot=12000 + hard walls). Analize things
* S1 with PBC and realistic setup for different $\phi$ (Kinetic energy VS $\phi$)
* MSD of the particles in different moment of the self-assembly process (Understand cage size and cage escape times)
* Energy map and space distribution of particles contributing to the fat tails

:::success
green is for finished tasks
:::

## Direct Coexistence LAMMPS

### Analysis 2024

#### 03/04/2024

Profiles (density and temperature) averaged over 5 independent realizations for $h=0.0795$

![FIG_AveProfile3phi_AveOverRand_DCper_Ntot3600_h0.0795_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE](https://hackmd.io/_uploads/BymZjT9k0.png)

![FIG_AveProfile3phi_TempBig_AveOverRand_DCper_Ntot3600_h0.0795_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr10](https://hackmd.io/_uploads/Hk4zspckC.png)

![FIG_AveProfile3phi_TempSmall_AveOverRand_DCper_Ntot3600_h0.0795_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE](https://hackmd.io/_uploads/B1fmoTq1C.png)

![FIG_AveProfile3phi_TempTot_AveOverRand_DCper_Ntot3600_h0.0795_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE](https://hackmd.io/_uploads/Hy-Lo6q10.png)

![FIG_AveProfile3phi_TempBigXYZ_AveOverRand_DCper_Ntot3600_h0.0795_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE](https://hackmd.io/_uploads/BJ9Uja9kA.png)

![FIG_AveProfile3phi_TempSmallXYZ_AveOverRand_DCper_Ntot3600_h0.0795_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE](https://hackmd.io/_uploads/H1LPs6c1R.png)

![FIG_AveProfile3phi_TempTotXYZ_AveOverRand_DCper_Ntot3600_h0.0795_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE](https://hackmd.io/_uploads/HyTDsa9yA.png)



#### 16/02/2024

Here are the **density** profiles for two independent sets of simulations. Each line is averaged over the three good $h$s previously identified.

nr 12664
![FIG_AveProfileAll_AveOverh_DCper_Ntot3600_Nrand12664_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr10](https://hackmd.io/_uploads/S14vffaoT.png)

nr 42424
![FIG_AveProfileAll_AveOverh_DCper_Ntot3600_Nrand42424_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr10](https://hackmd.io/_uploads/rkJcff6oT.png)

Here are the **temperature** profiles for one set of simulations. Each line is averaged over the three good $h$s previously identified.

nr 12664

XY

![FIG_AveProfileAll_TempBig_AveOverh_DCper_Ntot3600_Nrand12664_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr10](https://hackmd.io/_uploads/HkFtepmn6.png)

![FIG_AveProfileAll_TempSmall_AveOverh_DCper_Ntot3600_Nrand12664_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr10](https://hackmd.io/_uploads/HkW9eT7n6.png)

![FIG_AveProfileAll_TempTot_AveOverh_DCper_Ntot3600_Nrand12664_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr10](https://hackmd.io/_uploads/BJacxaQnp.png)

XYZ

![FIG_AveProfileAll_TempBigXYZ_AveOverh_DCper_Ntot3600_Nrand12664_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr10](https://hackmd.io/_uploads/rk7Rx6Q26.png)

![FIG_AveProfileAll_TempSmallXYZ_AveOverh_DCper_Ntot3600_Nrand12664_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr10](https://hackmd.io/_uploads/r1iClpQ26.png)

![FIG_AveProfileAll_TempTotXYZ_AveOverh_DCper_Ntot3600_Nrand12664_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr10](https://hackmd.io/_uploads/r1QJ-6Q2p.png)




#### 14/02/2024

In order to write the long paper 'Understanding self-assembly' I restart the analysis on direct coexistence. At the beginning I'm focusing only on old data. 
I take again the global phi [0.83,0.835,0.84,0.845,0.85]. 

Other things fixed:
DCper_Ntot3600_Nrand12664_h0.0783_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr10

##### Lattice spacing analysis

Check on lattice spacing for all the phi in which S1 appears. I check that for the the three good h I identified $a_x/a_y$ is close to 1. I do averages over the last 10 snapshots (nr 42424) 
![FIG_lattSpVariohL10F_DCper_Ntot3600_AllPhi_Nrand42424_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr](https://hackmd.io/_uploads/BygDOOcip.png)
For all the $\phi$ in the coexistence I observe $a_x/a_y\sim 1$ for the three good $h$ identified. Is it strange that the more crystalline one is non-monotonic?
I cjecked that the simulation at h=0.0796 is a full crystal while the pther two are with a bit of liquid.


##### Stationarity analysis
Here I check stationarity for different $\phi$. For each $\phi$ I plot data for the identified three $h$ with the same color. 
![FIG_aveq4_DCper_Ntot3600_AllPhi_Nrand42424_fr350_amp0.000009_gaa1_gab1_dtf0.7_dumpf100000000_logf100000000_verTHE_npr10_allh](https://hackmd.io/_uploads/HyP_v_qip.png)
Maybe $\phi=0.855$ is not really a good crystal. 0.86 is now running.



### Simulation setup

#### Initial configuration S1
I fix the total number of grains to be four times a perfect square $N=4n^2$. In this way I can always have an initial crystal with $n$ rows and $2n$ columns. I fix the y-lenght of the box $h$ and the packing fraction $\phi$. Then, the x-length of the box $L$ and the distance between particles $2r_{hi}+d$ in the initial crystal is chosen accordingly. All the particles are subjected to a very small random kick at the beginning of the simulation (this allows to have random initial condition).
![](https://i.imgur.com/IsWm6hr.png)
![](https://i.imgur.com/obfX2yZ.png)

#### Test for parallelization
![](https://i.imgur.com/VzTDCbb.png)
It slows down around 6 than it increases again... This is a bit strange to me. Could it be related to the anisotropic (i.e. rectangular) shape of the box?
![](https://i.imgur.com/kluQR5i.png)

### Final run (in progress)

For now varying 6 $\phi$s each for three h

Chekc q4 for stability (realization Nr 42424)
![](https://i.imgur.com/eBkreXy.png)

$\phi=0.83$ is fully fluid. $\phi=0.855$ fully crystal (still running).


Density profiles mediated over 150 snapshots

All packing fraction togheter
![](https://i.imgur.com/AcdaPlX.png)
Local density of crystal and liquid phase as a function of packing fraction (inside coexistence region)
![](https://i.imgur.com/2AgxeJX.png)

Single Packing fractions
![](https://i.imgur.com/gRMkBjr.png)
![](https://i.imgur.com/rrjKyzY.png)
![](https://i.imgur.com/L42rPEy.png)
![](https://i.imgur.com/0vuaBO2.png)
![](https://i.imgur.com/CD9lPdt.png)

Analysis lattice spacing with final run
![](https://i.imgur.com/KsJBZ6C.png)
This plots confirm that the 3 $h$ we have chosen are good because they produce a crystal with a symmetric lattice spacing





### Preliminary run

#### longer runs ($>2$ bln dt)
Longer simulations (between 3 and 5.6 bilions timesteps). 
**xy lattice spacing**
For each last snapshots I take particles belonging to s1 enviroment and I compute the ratio between the horizontal lattice spacing and the vertical one in the last snapshot. So in this case there are points that come from longer simulations and other from shorte ones.
Average over particles in the last snapshot
![](https://i.imgur.com/wyofsB7.png)
![](https://i.imgur.com/StOeKoL.png)
Average over last 5 snapshots
![](https://i.imgur.com/pRL8wv2.png)
![](https://i.imgur.com/qTHAnGT.png)

COMMENT: strange behaviour until $h=$ 0.0785 then it is better and for both packing fractions the point at $h=0.0795$ is very close to $a_x/a_y=1$.

I show the last snapshot of the "best" condition for future runs namely $\phi=0.84$ and $h=0.0795$
![](https://i.imgur.com/sl19oSC.jpg)
Here are the q4 and the phi profile of the last snapshot
![](https://i.imgur.com/N2OIWIK.png)





#### longer runs (2 bln dt)
Longer (2 bilion timesteps) runs but probably still not enough to have full stationarity.

Last frame $\phi=0.84$ $h=0.0775$
![](https://i.imgur.com/JLFzL9U.jpg)
Last frame $\phi=0.84$ $h=0.0785$
![](https://i.imgur.com/OSd48sX.jpg)

I plot q4 over time to check stationarity
$\phi=0.84$ varying h
![](https://i.imgur.com/AtzUdsP.png)

Last frame $\phi=0.85$ $h=0.0775$
![](https://i.imgur.com/9g5SX0D.jpg)
Last frame $\phi=0.85$ $h=0.0785$
![](https://i.imgur.com/wurYnQ2.jpg)


$\phi=0.85$ varying h
![](https://i.imgur.com/AJj96xx.png)

Comparison between the two packing fractions
![](https://i.imgur.com/dqO62jV.png)


COOMENT: $\phi=0.84$ seems to have reached stationarity while $\phi=0.85$ is clearly slower and we observe more variability changing h. Probably we would define eroorbars for fixed $\phi$ at different $h$?

**xy lattice spacing**
For each last snapshots I take particles belonging to s1 enviroment and I compute the ratio between the horizontal lattice spacing and the vertical one in the last snapshot.
![](https://i.imgur.com/5Jqihiq.png)
![](https://i.imgur.com/Lk3mbYM.png)
Strange behaviour.. Probably the correct h-window in which we have to perform this analysis starts around 0.0783. From the discussion we had before christmas I remember that we would expect a monotonous decreasing of this ratio as a function of h. Indeed, from 0.0783 it seems to be monotonously decreasing and approaching 1. Simulations up to h=0.0795 are running.


#### short runs
I've tried short runs (20 mln timesteps) with two different packing fraction varying $h$. I show some cases. Q4 as a function of time id plotted to check stationarity.

$\phi=0.85$ $h=0.0775$
![](https://i.imgur.com/oo0eEPv.gif)
![](https://i.imgur.com/yI4DPGq.png)

$\phi=0.85$ $h=0.0785$
![](https://i.imgur.com/vAT0utn.gif)
![](https://i.imgur.com/raiay5n.png)

$\phi=0.84$ $h=0.0775$
![](https://i.imgur.com/ePNvSoL.gif)
![](https://i.imgur.com/l625TkN.png)

$\phi=0.84$ $h=0.0785$
![](https://i.imgur.com/kr9JlhF.gif)
![](https://i.imgur.com/ligyZbG.png)

COMMENT:
These simulations are clearly too short to draw conclusions. Anyway it seems that at $\phi0.85$ the "final" amount of crystal slightly depends on h but maybe it is only a transient.


Simulations with 20 bilions of timesteps are running

## Numerical granular QC8 

### Periodic boundaries (feb/mar 2023)

#### N-phi scan (d=2.5 mm)

![](https://i.imgur.com/bKfxjdW.jpg)
![](https://i.imgur.com/sEH7bln.jpg)
![](https://i.imgur.com/2G1JEla.jpg)

#### Other state points (d=2.5 mm)


| Ntot5000_q0.5_xL0.3_phi0.862_dhi0.0025_Nrand67542_fr350_amp0.000009_gaa1_gab1_dtf0.7_verD | Ntot5000_q0.5_xL0.3_phi0.864_dhi0.0025_Nrand67542_fr350_amp0.000009_gaa1_gab1_dtf0.7_verD |
| -------- | -------- |
|  ![](https://i.imgur.com/Y240pLH.png) |    ![](https://i.imgur.com/js3jcol.png)|
|     ![](https://i.imgur.com/c98TUSx.png) |      ![](https://i.imgur.com/H9XyG53.png)  |
| ![](https://i.imgur.com/LucrtfZ.png) | ![](https://i.imgur.com/imRSjxf.png) |

| Ntot5000_q0.5_xL0.325_phi0.86_dhi0.0025_Nrand67542_fr350_amp0.000009_gaa1_gab1_dtf0.7_verD | Ntot5000_q0.5_xL0.35_phi0.86_dhi0.0025_Nrand67542_fr350_amp0.000009_gaa1_gab1_dtf0.7_verD |
| -------- | -------- |
| ![](https://i.imgur.com/znAkgQv.png) |  ![](https://i.imgur.com/jvzAVYw.png) |
| ![](https://i.imgur.com/Gc1tcRI.png) |  ![](https://i.imgur.com/yMFQG59.png) |
| ![](https://i.imgur.com/IzlJkV6.png) |  ![](https://i.imgur.com/FBpFcR4.png) |

### Hard boundaries
last configurations (still running)
**N=5000**
![](https://i.imgur.com/zO2ePti.jpg)
![](https://i.imgur.com/T4g5Eb2.jpg)

**N=12000**
![](https://i.imgur.com/GCk6ujT.jpg)
![](https://i.imgur.com/qdj8ZRy.jpg)




### Periodic boundaries (old)

I put here now (12/2022) the results obtained in spring.
![](https://i.imgur.com/QYZbQmd.png)
Best results in the central row
![](https://i.imgur.com/OiSBB5g.gif)

Tile analysis 
![](https://i.imgur.com/niUjvau.png)
![](https://i.imgur.com/5Uv8jwf.png)
![](https://i.imgur.com/MfahXr4.png)
![](https://i.imgur.com/hDQyem8.png)
![](https://i.imgur.com/3Hwt4EF.png)
![](https://i.imgur.com/CWsM2xi.png)

With imposed alignment
![](https://i.imgur.com/csmDw4x.png)
Over time
![](https://i.imgur.com/iL9Y1LP.png)
![](https://i.imgur.com/B7aLVyH.png)
![](https://i.imgur.com/5KAlrVl.png)
In this case it seems we reached a plateau


## Analisi Gravity
Comparison between no tilted plate and slightly tilted.
NO GRAVITY
![](https://i.imgur.com/oCPEwqd.gif)
GRAVITY
![](https://i.imgur.com/Nn7HxLr.gif)




## Looking for QC12
Finished simulations with better initial condition, nothing really good probably need more time. 
Promising to restart are these state point:
q0.5_xL0.65_phi0.84
q0.5_xL0.6_phi0.84

You should do simulation N5000 with q0.5_xL0.65_phi0.84

Some state points N=2000
![](https://i.imgur.com/7lYiYa6.jpg)

Different state points for three different initial conditions

$x_L=0.55$ $\phi=0.85$
![](https://i.imgur.com/VomlH67.jpg)

$x_L=0.65$ $\phi=0.85$
![](https://i.imgur.com/Hf7gx6L.jpg)

$x_L=0.7$ $\phi=0.84$
![](https://i.imgur.com/40EctvB.jpg)

$x_L=0.6$ $\phi=0.84$
![](https://i.imgur.com/EKlWsG2.jpg)

Maybe this is the bast one $x_L=0.6$ $\phi=0.84$ random init2
![](https://i.imgur.com/aSGlmOI.png)
![](https://i.imgur.com/9W7xSa7.png)

![](https://i.imgur.com/SJlM5nM.png)



## Looking for Rapha phase with dimers in LAMMP
**q=0.38 (Best snapshots)**
![](https://i.imgur.com/kOQE8pK.png)
![](https://i.imgur.com/zYHtLVL.png)
These last two are from the same stat point
do longer/larger and more random realizations
![](https://i.imgur.com/g1b8SBD.png)
![](https://i.imgur.com/9DCqGVv.png)
These last two are from the same stat point. Good out of 8 realizations. Try longer/larger.
![](https://i.imgur.com/EFsy1XJ.png)
![](https://i.imgur.com/P8OKyVU.png)
These last two are from the same stat point. Good out of 8 realizations. Try longer/larger.
![](https://i.imgur.com/6SHuyo4.png)
is it good?

Longer simulations with N=2000 are running

## Crystal phases as a function of $\phi$

###  S1 with PBC and realistic setup for different $\phi$ (Kinetic energy VS $\phi$)
Simulation with xL=0.5, q=0.476, $A=9$ $\mu m$, $f=350$ Hz
![](https://i.imgur.com/M7kL6eq.png)
Energy in the xy plane: Small particles are more energetic than the large ones. Energy increases with the packing fraction (effect of more collisions transferring energy from z to xy), then seems to saturate but as S1 starts to form there is a jump.
![](https://i.imgur.com/eXIQ9BJ.png)
Total energy and z energy show that the global contrbute of collisions is dissipative indeed the trend up to $\phi \sim 0.83$ is decreasing. Then (see the ZOOM below) we still observe a jump in the total kinetic energy. This tells that particles in S1 are arranged in such a way to perform less collisions or less dissipative collisions (check with MSD coming soon). So there is not "only" a more efficient transfer of energy from z to xy at the transition but there is a real optimization of the energy trasfer from the plate to the granular system.
![](https://i.imgur.com/7lMeT9v.png)
ZOOM of the total kinetic energy near the transition.
![](https://i.imgur.com/O6mhPOT.png)
This is just to check that the jump in energy corresponds to the jump in the q4. 

#### Comparison with hexagonal crystal
This is what happens when we perform the same analisys in a monodispersed system around the liquid-solid transition. As we expect from previous studies, the transition occurs at a smaller packing fraction.
![](https://i.imgur.com/srArWTO.png)
In the xy plane we see that the energy increases with $\phi$ but it not present a marked jump at the transition as we observed in the bidispersed case. The last decrease of the energy is probably due to the fact that dissipation in collisions becomes dominant with respect z $\to$ xy energy transfer.
![](https://i.imgur.com/okNNpPl.png)
The total and z kinetic energy are globally decreasing with a little jump near the transition. Nevertheless, this jump is not so evident as in the bidispersed case. 

* To be done: Simulation of larger systems
* To be done: MSD to see cage size

## QC8 realistic setup (Ntot=12000 + hard walls)
* 4 simulations for some good values of $q$,$\phi$,$x_L$ are running
* Find the reuslts on the cluster

## MSD of the particles in different moment of the self-assembly process 
* In progress on QC8
* Will be applies on simulations now running. 
##  Energy map and space distribution of particles contributing to the fat tails

:::success
* Energy Map: Each particle is colored according to its kinetic energy. Not really interesting by eyes (see figures).  Maybe one can see some correlation between being in the fluid defects and not being yellow but probably the correlation between kin erg and q4 would give much bettere evidence

:::
EnergyMap_n76_GR_rLL0.00125_rSS0.000595_fr350_amp0.000009_NL10000_NS10000_packfrac0.850_dtfactor1.0_roof_hw_wiggle.dump **Only Large**
![](https://i.imgur.com/fmfF78p.png)
EnergyMap_n76_GR_rLL0.00125_rSS0.000595_fr300_amp0.000010_NL10000_NS10000_packfrac0.850_dtfactor1.0_roof_hw_wiggle.dump **Only Large**
![](https://i.imgur.com/zgF11cM.png)
EnergyMap_All_n76_GR_rLL0.00125_rSS0.000595_fr300_amp0.000010_NL10000_NS10000_packfrac0.850_dtfactor1.0_roof_hw_wiggle.dump **All with Fat tail Particles (red)**
![](https://imgur.com/gXtY4WM.png)
EnergyMap_All_n76_GR_rLL0.00125_rSS0.000595_fr350_amp0.000009_NL10000_NS10000_packfrac0.850_dtfactor1.0_roof_hw_wiggle.dump **All with Fat tail Particles (red)**
![](https://imgur.com/vrKnTbx.png)
Particles contributing to fat tails do not clusterize
* Correlation between q4 and kinetic energy 
 
### OPEN QUESTIONS

*  ?


