---
tags: SoftQC, self-assembly
---
# Stability of hard disks self-assembly

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

---

---

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

Frank, Giuseppe, Etienne

---

---

## Background

We have computed the infinite pressure phase diagram of binary (non-)additive hard disk mixtures in the size-ratio--composition plane.
This revealed several random-tiling quasicrystal regions.
https://arxiv.org/abs/2003.08889

We now would like to investigate the stability of the QC at finite pressures.

As a first step, we will use simple cell theory for free energy calculations of the crystal structures.

We will also perform large scale simulations of the mixtures at the points where QC is expected, to see whether we can equilibrate it over a large area, and compute some of its properties (diffraction diagram for instance).

## Plans

## To Do
:::success
green is for finished tasks
:::
### Simulate large systems of NAHD mixtures
:::success
* Particles overlap in large scale simulations : fix this bug !
* Run large scale simulations.
* Compute bond orientational parameters.
* Restart simulations $N_0=2000, N_1=1500$ density 1.65 and $N_0=2000, N_1=1000$ density 1.45.
* Compute scattering pattern.
* Scan packing fractions at $x_S=1/3$.
:::
* Run simulations with $q=0.47, \Delta\approx0.06725787749645651$ (in elongated box ?).
* Measure the number of squares and triangles of each orientation in the simulations.

### Compute finite pressure phase diagram at the size ratios of interest
:::success
* Learn about cell theory for free energy calculation.
* Compute the free energy of hexagonal crystals.
* Compute the free energy of the S1 crystal.
* Compute the free energy of the fluid.
:::
* Build the phase diagram.
* Include proper treatment of the random tiling.

### OPEN QUESTIONS

* What is the role of defects in random tiling QC ? Could we compute their equilibrium concentration.

## Large scale simulations

### First results

Simulations are still running. Here are snapshots of the systems as simulations go. When looking at the movies dynamically, regions of fluid appear clearly (particles move a lot). In particular, even though it is not very clear from snapshots, the systems with $x=15/35=3/7$ are composed of crystal-like fragments in a see of fluid.

|
:---------:|:---------:
*$N_S=2000$, $N_L=1000$, $t=0$*|*$N_S=2000$, $N_L=1000$, $t=250000$. HexL + fluid ?*

|
:---------:|:---------:
*$N_S=2000$, $N_L=1500$, $t=0$*|*$N_S=2000$, $N_L=1500$, $t=310000$. Random tiling + fluid ?*

|
:---------:|:---------:
*$N_S=2000$, $N_L=2000$, $t=0$*|*$N_S=2000$, $N_L=2000$, $t=183 000$. S1 + fluid ?*

>[name=Frank] Times are in MD time units? $\tau = \sqrt{\beta m \sigma_L^2}$? Might be nice to quantify the changes in the system via e.g. bond-orientational order. One way would be to define for particle $i$
>$$\phi_k (i)= \left|\frac{1}{N_b} \sum_{j}^{N_b} \exp(i k \theta_{ij})\right|,$$
>with $\theta_{ij}$ the angle between the bond between particles $i$ and $j$ and the x-axis, $k$ an integer, and $N_b$ the number of neighbors of $i$. I think if you consider only the large spheres in this calculation, particles in the square phase will have a high value of $\phi_4$, particles in the hex phase will have a high $\phi_6$, and particles in the random tiling will be mixed (and may have relatively high $\phi_5$ on average?). So plotting the average values of these quantities over time may let you see in what direction the system is moving.
>>[name=Etienne] Time is MD time yes. Good idea for the BOO, I will compute it !

### Boop calculation

Simulations ended at MD time $10^6$. I computed $q_4$, $q_5$ and $q_6$ only considering large particles and averaged it in each frame.

#### $N_0=2000, N1=1000$
* Density 1.45 ($\phi=0.836$):

$\tau=0$|$\tau=10^6$
:---:|:---:
|
|
*$q_5$ color coding.*|*$q_5$ color coding.*
|
*Number of neighbours.*|*Number of neighbours.*
|


Looking at the frame in sequence clearly shows channels of fluid between chunks of crysal. The crystals look like random tilings. From the boops plot, it is clear that the system has not reached equilibrium yet.
$q_5$ kind of catches the transition (second row), but is unable to discriminate clearly between crystal en fluid regions. The fluid regions are more easily seen as regions with lower number of neighbours (last row).
I will restart this run to try reaching equilibrium.

* Density 1.5 ($\phi=0.865$):


The system is frozen.

#### $N_0=2000, N_1=1500$

* Density 1.55 ($\phi=0.801$):

The system is fluid.

* Density 1.6 ($\phi=0.827$):

$\tau=0$|$\tau=10^6$
:---:|:---:
|
|
*Number of neighbours.*|*Number of neighbours*
|
*Bonds*|*Bonds*


Fluid with chunks of crystal floating around. The number of neighbours highlights crystal regions (2nd row). Crystal regions look like random tiling. In particular, dodecagonal wheels are observed (last row).
The system looks equilibrated. $q_5$ dominates.

* Density 1.65 ($\phi=0.853$):

$\tau=0$|$\tau=10^6$
:---:|:---:
|
|
*$q_4$*|*$q_4$*
|
*Number of neighbours*|*Number of neighbours*


The system is poly crystalline with pockets of liquid. The crystal is made of S1, HexL and random tiling chunks.
In contrast to previous cases, $q_4$ raises over time.
Finally, the system is still out of equilibrium : one crystal phase might take over the others.
I will restart this simuation as well.

#### $N_0=2000, N_1=2000$

* Density 1.75 ($\phi=0.826$):


The system is liquid. Note that it is not the same liquid as in the case of $N_0=2000, N_1=1500$ for density 1.55 : $q_6$ is lower (which is normal since there are more small particles to destroy order between large ones).

* Density 1.8 ($\phi=0.850$):

$\tau=0$|$\tau=10^6$
:---:|:---:
|
|
*$q_4$*|*$q_4$*


Very clear transition to S1, as expected for $x_S=0.5$ at infinite pressure.

* Density 1.85 ($\phi=0.874$):
*

System is jammed.

#### Discussion

Bond orientational order parameters are quite inefficient at detecting random tiling regions.
In square-triangle random tilings, local environments with 4, 5 and 6 neighbours. Therefore, none of the corresponding q's are specific to random tiling crystals.


To do : try using scattering pattern instead.

### Restarted runs

#### $\phi=0.836, N0=2000, N1=1000$
Restarted for an extra $5\times 10^5$ MD times.
Equilibrium is still not reached.
Boops : $q_5$ and $q_6$ and very much anti-correlated.
Reduced pressure drops. It looks like a significant relaxation event accurs aroud $\tau=8\times 10^5$, but I cannot spot anything special in the movie. More squares might temporarily be created as hinted by the small jump in $q_4$.

I computed structure factor for this run.

|
:---:|:---:

12-fold symmetry appears very clearly !

{%youtube 9SvLMkjeWv4 %}

*Last snapshot at $\tau=1.5\times 10^6$*

#### $\phi=0.853, N0=2000, N1=1500$
Restarted for an extra $5\times 10^5$ MD times.


|
:---:|:---:

{%youtube P2H1mpbSyqg %}

The system remains amorphous. It is polycristalline with small bits of HexL, S1 and random tiling attached to each other.

*Last snapshot at $\tau=1.5\times 10^6$*

### Spheres on a substrate

It might be experimentally relevent to consider pairs of size ratios $q$ and non-additivity parameters $\Delta$ corresponding to spheres lying on a flat surface.


From the figure :
$$(1+q)^2r_L^2=(1-q)^2r_L^2+d_{LS}^2\\
d_{LS}^2=r_L^2[(1+q)^2-(1-q)^2]\\
d_{LS}^2=4qr_L^2$$
So $d_{LS}=2r_L\sqrt{q}$
And by definition of $\Delta$ : $d_{LS}=(1+\Delta)(r_L+r_S)=r_L(1+\Delta)(1+q)$
Therefore $1+\Delta_{plan}=\frac{2\sqrt{q}}{(1+q)}$, ie.
$$\Delta_{plan}=\frac{2\sqrt{q}}{(1+q)}-1$$


We are interested in S1-HexL RT region. The following figure highlights the RT region at inifinite pressure along with the parameters corresponding to spheres on a planar substrate.


The large scale simulations above are performed at $q=0.45$ and $\Delta=0.05$, while $\Delta_{plan}(0.45) \approx 0.0747$.

### Density scan at $x_S=1/3, q=0.45, \Delta=0.05, N0=2000, N1=1000$

Packing fraction|Boops|Final structure factors|Pressure
:---:|:---:|:---:|:---:
$\phi=0.815$, fluid + HexL|||
$\phi=0.82$, fluid + RT (+HexL)|||
$\phi=0.825$, fluid + RT (+HexL)|||
$\phi=0.83$, fluid + RT|||
$\phi=0.836$, fluid + RT|||

The proportion of fluid decreases as the packing fraction increases. The fluid-crystal coexistence pressure is about $28 k_B T\sigma_L^2$.

### Density scan at $x_S=1/3, q=0.4, \Delta=0, N0=2000, N1=1000$

Packing fraction|Boops|Final structure factors|Pressure
:---:|:---:|:---:|:---:
$\phi=0.80$, fluid + HexL|||
$\phi=0.81$, fluid + polycrystalline HexL|||
$\phi=0.815$, fluid + polycrystalline HexL|||

For larger packing fractions, the system is jammed.
It seems that the coexistence of HexL with S1, predictec by cell theory, is not dynamically accessible.

**Next step:** use cell theory to try predicting a phase diagram. Since spheres on a planar substrate are experimentally relevent, focus on parameters that represent this situation.
We want:
* $\sigma_{LS}=(1-\Delta)(1+q)\sigma_L/2$ to be slightly smaller than $1/\sqrt{2}\sigma_L\approx0.707\sigma_L$ such that small particles can rattle inside squares,
* finite non-additivity $\Delta$ to promote mixing of small and large particles,
* the pair of parameters $(q,\Delta)$ to satisfy the planar constraint.

A sensible solution:
$$q=0.47\\ \Delta\approx0.06725787749645651$$

With this set of parameters, $\sigma_{LS}\approx0.6856$ close to $0.68875$, the value corresponding to $q=0.45, \Delta=0.05$ where we know that we can equilibrate random tilings.

### Density scan at $x_S=1/3, q=0.457, \Delta=-0.07, N0=2000, N1=1000$

Packing fraction|Boops|Final structure factors|Final frame|Pressure
:---:|:---:|:---:|:---:|:---:
$\phi=0.81$, fluid + HexL||||
$\phi=0.82$, fluid + HexL||||
$\phi=0.83$, fluid + RT||||
$\phi=0.84$, fluid + RT||||
$\phi=0.85$, RT,S1 ? (very slow dynamics)||||

### Density scan at $x_S=0.429, q=0.457, \Delta=-0.07, N0=2000, N1=1500$

Packing fraction|Boops|Final structure factors|Final frame|Pressure
:---:|:---:|:---:|:---:|:---:
$\phi=0.83$, fluid||||
$\phi=0.84$, RT, S1, fluid||||
$\phi=0.85$, RT, S1, fluid||||

### Density scan at $x_S=0.5, q=0.457, \Delta=-0.07, N0=2000, N1=2000$

Packing fraction|Boops|Final structure factors|Final frame|Pressure
:---:|:---:|:---:|:---:|:---:
$\phi=0.83$, fluid||||
$\phi=0.84$, S1+fluid|||
$\phi=0.85$, S1+fluid||||
$\phi=0.86$, S1+fluid||||
$\phi=0.87$, S1 (polycrystalline)||||
$\phi=0.88$, jammed||||

## Finite pressure phase diagram

### General strategy

We want to build the phase diagram of (NA)HD with fixed size ratio $q$ and non-additivity parameter $\Delta$, in the $p-x_S$ plane because only two phases can coexist at the same time in this representation.
However, it might be necessary to go to $\eta_S-\eta_L$ at some point since these are the parameters accessible to experimentalists.
According to Gibbs phase rule, if we fix $N, p, T$ and $x_S$, we should have coexistence regions in the phase diagram.

*Wild guess of the shape of the phase diagram, just to have ideas*

>[name=Frank] The phase transitions for the pure phases (so at x=0 and x=1) occur at the same scaled pressure. So for the fluid to Hex_L+F point, $\beta p \sigma_L^2$ is the same as $\beta p \sigma_S^2$ for the fluid to Hex_S + F point at x=1. Which means that the ratio between those pressures is just the size ratio squared. (With the small-particle one being the higher pressure)
What can be nice in the figure is to write the coexistences in the left-right order of the phases. So where you now have S1 + Hex_L, you could have the reverse, since Hex_L is on the left, etc.

We will construct this phase diagram horizontal slice by horizontal slice (at constat pressure).
At a given pressure, we need to compute, for each phase in the game $G(x_S)$, and then construct the common tangent that minimises the Gibbs free energy.

*Wild guess for a common tangent construction look, just to have ideas again. This doesn't even correspond to any slice of wild guess 1 !*

For the crystal phases (HexL, HexS and S1) the free energy can be determined for any value of the pressure with cell theory (see below). The fluid free energy could be obtained by thermodynamic integration. Equations of state at fixed $x_S$ could be obtained from NPT MC simulations (or NVT EDMD) and then interpolated to get $G_{fluid}(x_S)$.

A few questions/caveats:
* Computation of the fluid free energy as a function of composition seems to be a bit heavy. Is there a way to get it directly without computing many equations of state ?
* Shall we treat the random tiling as a phase on its own ? It is clear that simple HexL-S1 coexistence is not equivalent to RT, since it lacks the configurational entropy configuration. Maybe we could simply say $G_{RT}=G_{coex} - TS_{RT}(x_S)$. The coexistence term takes care of the vibrationnal entropy while the second is the RT entropy that we can get from (eg.) Kalugin's paper. Note that this ignores local environments with 5 neighbours, which are the most common in a RT. We could perform perfect RT tiling sambling with zipper moves and estimate the proportion of each type of site from there. At finite pressure though, vibrationnal entropy might change that distribution by favouring some environments.
* In cell theory, we assume that the crystal phases expand uniformly when pressure decreases. In reality, unit cell deformation are likely to occur. Can we quantify the error made in cell theory approximation, and eventually compare it with other sources of error like this one, or the fact that in first approximation, we might approximate real accessible volumes by simpler shapes ?
* According to Laura, cell theory gives reasonnable results when comparing crystals free energies. However, comparing cell crystals free energies from cell theory with fluid free energy obtained by other means (eg: thermodynamic integration) may lead to significant errors.

### Cell theory

Refs.
Historical paper, Lennard-Jones and Devonshire :
10.1098/rspa.1937.0210
More involved methods, Cottin&Monson :
10.1063/1.472832
10.1063/1.465560

#### General case

Consider a crystal in 2D of $N$ particles occupying $n$ types of non-equivalent sites. Let $N_i$ be the number of particles occupying a site of type $i$.
The partition function reads
$$Z=\frac{1}{N_1!\,\cdots\,N_n!}\frac{1}{h^{2N}} \int\prod_{i=1}^{N}d\vec{p_i}\int\prod_{i=1}^{N}d\vec{r_i}e^{-\beta H(\{\vec{p_i}\},\{\vec{r_i}\})} \\
=\frac{1}{N_1!\,\cdots\,N_n!}\frac{1}{\Lambda^{2N}}\int\prod_{i=1}^{N}d\vec{r_i}e^{-\beta U(\{\vec{r_i}\})}$$
with $H$ the full Hamiltonian of the system, $U$ the potential energy and $\Lambda=\sqrt{mk_BT/\hbar}$ the thermal wavelength (assuming all particles have the same mass $m$).

We now assume that each particle is trapped in a cage formed by its neirest neighbours sitting at their ideal lattice sites, and that the partition function can be factorised (mean field approximation).

$$Z\approx\frac{N!}{N_1!\,\cdots\,N_n!}\frac{1}{\Lambda^{2N}}\prod_{j=1}^{n}\left(\int d\vec{r}e^{-\beta U_j(\vec{r})}\right)^{N_j}\\
=\frac{N!}{N_1!\,\cdots\,N_n!}\frac{1}{\Lambda^{2N}}\prod_{j=1}^{n}A_j^{N_j},\\
\text{with} \quad A_j\doteq \int d\vec{r}e^{-\beta U_j(\vec{r})}.$$

$U_j$ is the potential felt by the central particle rattling around a site of type $j$, caged by its neirest neighbours.

> [name=Frank] At this point, we have restricted each particle to a specific lattice site, for which we must have chosen a permutation of particles. Hence, the partition function now indeed picks up this extra factor $N!$.

Then the free energy reads
$$-\beta F= \ln Z = \ln N! + \sum_{j=1}^n \{ -\ln{N_j!}+N_j\ln{A_j} \} -2N\ln{\Lambda}\\
\approx\sum_{j=1}^n N_j\ln{\left(\frac{N A_j}{N_j\Lambda^2}\right)}.$$

Let $\phi_j = \lim_{N\rightarrow \infty ,N_j\rightarrow \infty} N_j/N$ the fraction of sites of type j in the crystal. Then

$$\frac{\beta F}{N} \approx\sum_{j=1}^n \phi_j \ln{\left(\frac{\phi_j\Lambda^2}{A_j}\right)}.$$

The pressure is given by
$$p = -\frac{\partial F}{\partial V}
= -\frac{N}{\beta} \sum_{j=1}^n \phi_j \ln \frac{\phi_j}{A_j}=Nk_BT\sum_{j=1}^n\frac{\phi_j}{A_j}\frac{\partial A_j}{\partial V}.$$

Sanity check : for a perfect gas, there is only one type of site and $A=V/N$. Then, $p=k_BTN/V$, *ie.* $pV=Nk_BT$ as expected.

> [name=Frank] ... Interesting. To be honest, I'm surprised that works, given that we have clearly made the assumption that each particle is tied to a lattice site. I still agree with your expresion, though.
Also useful: if we define $f = F/N$, then $p = \rho^2 \frac{\partial f}{\partial \rho}$, with $\rho = N/V$.
> [name=Etienne] There is a discussion of this in Lennard-Jones paper (Critical phenomena in gases, 1937). In fact, the perfect gas is the case for which the cell theory approximation is the worst. Indeed, the relative error in free energy scales with $N$. However, the error only depends on $N$ and $T$, and therefore the equation of state is not affected.

For hard spheres, $A_j$ is the accessible volume for the central particle in site $j$. The cell approximation is no longer valid when the spacing between particles of the cage is large enough to let the central particle escape.

#### Hexagonal compact packing at $x_S = 0$ or $1$

There is only one type of site in the hexagonal crystal.
*Accessible volume for the central particle is depicted in pink.*

Let $L$ be the lattice parameter and $\sigma$ the particles diameter. Then
$$
\rho\doteq\frac{N}{V}=\frac{2N}{NL^2\sqrt{3}}=\frac{2}{L^2\sqrt{3}}\\
L=\sqrt{\frac{2}{\sqrt{3}\rho}}
$$

The accessible volume for the central particle $A_{\text{hex}}$ reads
$$
A_{\text{hex}}(\rho)\approx \pi(L-\sigma)^2=\pi\left(\sqrt{\frac{2}{\sqrt{3}\rho}}-\sigma\right)^2
$$

The maximum density is achieved when all disks are in contact, *ie.* $L=\sigma$, so $\rho_{\text{max}}=2/(\sigma^2\sqrt{3})$.
When $L=2\sigma$, the spacing between particles is large enough to let the central particle escape, the cell theory approximation is not valid anymore. Therefore $\rho_{\text{min}}=1/(2\sqrt{3}\sigma^2)$.

For densities $\rho$ between $\rho_{\text{min}}$ and $\rho_{\text{max}}$, cell theory allows to compute an approximate free energy.

$$
\frac{\beta F(\rho)}{N}\approx\ln\left(\frac{\Lambda^2}{A_{\text{hex}}(\rho)}\right)=2\ln\left(\frac{\Lambda}{\sqrt{\frac{2}{\sqrt{3\rho}}}-\sigma}\right)-\ln\pi
$$

Since we are interested in Gibbs free energy $G(N,p,T,x_S)$, we need to compute the pressure.
$$
\frac{\beta p}{\rho}=\rho\frac{\partial}{\partial \rho}\left(\frac{\beta F}{N}\right)
=-2\rho\frac{\partial}{\partial \rho}\ln\left(\sqrt{\frac{2}{\sqrt{3}\rho}}-\sigma\right)
=\frac{1}{1-\sqrt{\frac{\sqrt{3}\rho}{2}}\sigma}
$$

This relation can be inverted to get $\rho$ as a function of $p$.
$$
\rho(p)=\beta p + \frac{\sqrt{3} \beta^{2} p^{2} \sigma^{2}}{4} -\frac{\beta^{\frac{3}{2}} p^{\frac{3}{2}} \sigma \sqrt{3 \beta p \sigma^{2} + 8 \sqrt{3}}}{4}
$$

Then, Gibbs free energy is given by
$$
\frac{\beta G}{N}=\frac{\beta F}{N} + \frac{\beta p}{\rho}=2\ln\left(\frac{\Lambda}{\sqrt{\frac{2}{\sqrt{3\rho}}}-\sigma}\right)-\ln\pi+\frac{1}{1-\sqrt{\frac{\sqrt{3}\rho}{2}}\sigma}
$$
One can express it a function of $p$ by mixing the two last equations. The resulting expression is pretty long.
To compute actual values, without loss of generality, we can take $\beta=1$ since hard disk systems are athermal and $\Lambda=1$ fixes particles mass, which does not affect the phase behaviour. In our binary systems, since we measure lengths in large disk diameters, we take $\sigma=1\sigma_{LL}$ for hexagonal crystal of large particles and $\sigma=q\sigma_{LL}$ for hexagonal crystal of small ones.

*Gibbs free energy of hexagonal crystals of large small disks as a function of pressure. Shouldn't they be convex ??*
> [name=Frank] It does not have to be, no. $\left(\frac{\partial G}{\partial p}\right)_{N,T} = V$, so the slope of the plot should just be the inverse density... which decreases as the pressure goes up. Note that since you can't have a coexistence between states at different pressures, concavity here does not give you an instability (like it would if you had a concave Helmholtz free as a function of V). Also perhaps useful to consider is the fact that the Gibbs free energy of an ideal gas is just $\beta G_{id}/N = \log (\beta P \Lambda^3)$, which is clearly not convex.

#### S1 crystal at $x_S = 0.5$

There are two types of sites in the S1 crystal. They both appear with the same frequency $\phi_1=\phi_2=0.5$.

*Type 1 and 2 sites in S1 crystal, respectively hosting a small and large particle.*

Let $L$ be the lattice parameter. $\rho=2/L^2$ so $L=\sqrt{2/\rho}$.
The minimum density is achieved when the small disk can escape the cage of large ones, *ie.* $L=(1+q)\delta\sigma_{LL}$. So $\rho_{\text{min}}=2/((1+q)\delta\sigma_{LL})^2$.
The maximum density is achieved when $L=\sigma_{LL}$ (large disks in contact) or $L=(1+q)\delta\sigma_{LL}/\sqrt{2}$ (large and small disks touching along the square diagonal).
The first case occurs for $q<\sqrt{2}/\delta -1$ and the second for larger $q$'s. The discriminating value corresponds to a magic ratio of S1.
The maximum densities read respectively $\rho_{\text{max}}=2/\sigma_{LL}^2$ or $\rho_{\text{max}}=4/((1+q)\delta\sigma_{LL})^2$.

We now need to compute the contribution to the free energy and pressure of both of the sites.

##### Type 1 site (holding a small particle)

We approximate $A_1(\rho)$ by a circle of radius $L/\sqrt{2} - (1+q)\delta\sigma_{LL}/2=1/\sqrt{\rho} - (1+q)\delta\sigma_{LL}/2$.

Therefore
$$
A_1(\rho)\approx\pi\left(\sqrt{\frac{1}{\rho}}-(1+q)\frac{\delta\sigma_{LL}}{2}\right)^2
$$

The corresponding contribution to the free energy is
$$
\frac{\beta F^{(1)}}{N}\approx\frac{1}{2}\ln\left(\frac{\Lambda^2}{2A_1(\rho)}\right)=\ln\left(\frac{\Lambda}{\sqrt{\frac{1}{\rho}}-(1+q)\frac{\delta\sigma_{LL}}{2}}\right) -\ln(\sqrt{2\pi})
$$
and the associatied pressure reads
$$
\frac{\beta p^{(1)}}{\rho}=\rho\frac{\partial}{\partial \rho}\frac{\beta F^{(1)}}{N}=\frac{1}{2}\frac{1}{1-(1+q)\frac{\delta\sigma_{LL}}{2}\sqrt{\rho}}
$$

##### Type 2 site (holding a large particle)

The smallest distances the central particle can move from its ideal lattice site before colliding with a neighouring small or large particle are given respectively by
$$
d_S=\frac{L}{\sqrt{2}} - (1+q)\frac{\delta\sigma_{LL}}{2}=\sqrt{\frac{1}{\rho}}-(1+q)\frac{\delta\sigma_{LL}}{2}\\
d_L=L-\sigma_{LL} = \sqrt{\frac{2}{\rho}}-\sigma_{LL}
$$
We approximate the accessible volume for the central particle as $A_2(\rho)=\min\{A_2^S(\rho), A_2^L(\rho)\}=\min\{\pi d_S^2(\rho), \pi d_L^2(\rho)\}$.
$A_2^S(\rho) = A_1(\rho)$ so the free energy contribution from type-2 sites is equal to that of type-1 sites if $d_S(\rho)<d_L(\rho)$.
The free energy contribution from type-2 sites when $d_L(\rho)<d_S(\rho)$ is
$$
\frac{\beta F^{(2-L)}}{N} \approx \frac{1}{2}\ln\left(\frac{\Lambda^2}{2A_2^L(\rho)}\right) = \ln\left(\frac{\Lambda}{\sqrt{\frac{2}{\rho}}-\sigma_{LL}}\right)-\ln\sqrt{2\pi}
$$
The associated pressure reads
$$
\frac{\beta p^{(2-L)}}{\rho}=\rho\frac{\partial}{\partial \rho}\frac{\beta F^{(2-L)}}{N}=\frac{1}{2}\frac{1}{1-\frac{\sigma_{LL}}{\sqrt{2}}\sqrt{\rho}}
$$

##### S1 free energy

We need to discriminate between two cases:

* $d_S(\rho)<d_L(\rho)$:
$$\frac{\beta p}{\rho}=2\frac{\beta p^{(1)}}{\rho}=\frac{1}{1-(1+q)\frac{\delta\sigma_{LL}}{2}\sqrt{\rho}}$$
$$\frac{\beta G}{N}=2\frac{\beta F^{(1)}}{N}+\frac{\beta p}{\rho}=2\ln\left(\frac{\Lambda}{\sqrt{\frac{1}{\rho}}-(1+q)\frac{\delta\sigma_{LL}}{2}}\right) -\ln(2\pi)+\frac{\beta p}{\rho}$$
* $d_L(\rho)<d_S(\rho)$:
$$\frac{\beta p}{\rho}=\frac{\beta p^{(1)}}{\rho}+\frac{\beta p^{(2-L)}}{\rho}=\frac{1/2}{1-(1+q)\frac{\delta\sigma_{LL}}{2}\sqrt{\rho}}+\frac{1/2}{1-\frac{\sigma_{LL}}{\sqrt{2}}\sqrt{\rho}}$$
$$\frac{\beta G}{N}=\frac{\beta (F^{(1)}+F^{(2-L)})}{N}+ \frac{\beta p}{\rho}=\ln\left(\frac{\Lambda^2}{\left(\sqrt{\frac{1}{\rho}}-(1+q)\frac{\delta\sigma_{LL}}{2}\right)\left(\sqrt{\frac{2}{\rho}}-\sigma_{LL}\right)}\right) -\ln(2\pi)+\frac{\beta p}{\rho}$$

##### A few plots

All plots are for a S1 crystal with $q=0.45$ and $\delta=0.95$ (same parameters as the scan above).

The two cases are indeed necessary: for the large particle, the closest neighbour switches from a small particle to a large one at $\rho\approx 1.8$.

The free energy is a continuous function of the accessible volume, so the two free energies cross as expected.

However, the derivatives don't match at the crossing point, hence there is an artificial jump in the pressure (and thus in the Gibbs free energy as well). The mismatch is actually quite big. What can we do ? It is clearly no physical.
The obvious alternative I see is doing a proper sampling of the accessible volume (with MC for instance), so we don't have to switch between two functional forms at some point. It is clearly heavier, but also more precise.
Another solution could be to stitch a second order polynomial on a small region to join the two free energies with matching derivatives (a bit like a spline).
Third solution would be to interpolate the whole curve with cubic splines to smooth the kink.
Fourth solution is to find another approximation for the accessible volume that has a continuous derivative as a function of density. **This is the one we choose.**

##### Type-2 sites, second take

We approximate the accessible volume to the central large disk by the Voronoi cell.
When the closest neighbour is a large disk, the accessible volume is a square of area
$$
\rho A_2^{(LL)}(\rho)=4(\sqrt{\rho}\sigma-\sqrt{2})^2.
$$
When a small disk is encountered before a large one, the Voronoi cell has the shape depicted above and its volume reads
$$
\rho A_2^{(LS)}(\rho)=4(\sqrt{\rho}\sigma-\sqrt{2})^2-[2+2\sqrt{2}(\sqrt{\rho}\sigma-\sqrt{2})-\delta\sqrt{\rho}\sigma (q+1)]^2.
$$
The two corresponding curves cross when $d_{LS}=\sqrt{2} d_{LL}$, which occurs at $\sigma^2\rho_c=8/(\sqrt{2}\delta(1+q)-4)^2$.

At the crossing point, derivatives match so there will be no jump in the pressure !

The contribution to the free energy is
$$
\frac{\beta F_2}{N}=\Theta(\rho_c-\rho)\frac{\beta F_2^{(LS)}}{N}+\Theta(\rho-\rho_c)\frac{\beta F_2^{(LL)}}{N}\\
=0.5\Theta(\rho_c-\rho)\log{\left(\frac{0.5 \Lambda^{2} \rho}{4 \left(\sqrt{\rho} \sigma - \sqrt{2}\right)^{2} - \left(- \delta \sqrt{\rho} \sigma \left(q + 1\right) + 2 \sqrt{2} \left(\sqrt{\rho} \sigma - \sqrt{2}\right) + 2\right)^{2}} \right)}+0.5\Theta(\rho-\rho_c)\log{\left(\frac{0.125 \Lambda^{2} \rho}{\left(\sqrt{\rho} \sigma - \sqrt{2}\right)^{2}} \right)}
$$
The associated pressure is
$$
\frac{\beta p_2}{\rho}=\rho\frac{\partial}{\partial \rho}\frac{\beta F_2}{N}\\
=\Theta(\rho_c-\rho)\frac{ \left(\sqrt{\rho} \sigma \left( \delta q + \delta\right) - 2\right)}{\sqrt{\rho} \sigma \left(4 \delta q + 4 \delta\right) + \rho \sigma^{2} \left( \delta^{2} q^{2} + 2 \delta^{2} q + \delta^{2} - 4 \sqrt{2} \delta q - 4 \sqrt{2} \delta + 4\right) - 4}+\Theta(\rho-\rho_c)\frac{0.5 \sqrt{2}}{\sqrt{2}-\sqrt{\rho}\sigma}
$$

*Contribution to the pressure of type-2 sites*
We can invert (numerically) this relation to get $\rho(p_2)$.
Then we can compute the Gibbs free energy for type-2 sites $\beta G_2/N=\beta F_2/N+\beta p_2/\rho$.

##### Type-1 sites, second take

To be consistent we recompute site-1 free energy with Voronoi cell as free volume approximation.

The accessible volume is given by
$$
\rho A_1(\rho)=(\delta\sqrt{\rho}\sigma(q+1)-2)^2.
$$
The free energy reads
$$
\frac{\beta F_1}{N}=0.5 \log{\left(\frac{\Lambda^{2} \rho}{2 \left(\delta \sqrt{\rho} \sigma \left(q + 1\right) - 2\right)^{2}} \right)}.
$$
And the associated pressure is
$$\frac{\beta p_1}{\rho}=\frac{1}{2-\delta \sqrt{\rho} \sigma \left(q + 1\right)}.$$
Finally, Gibbs free energy is
$$
\frac{\beta G_1}{N}=\frac{0.5 \left(\delta \sqrt{\rho} \sigma \left(q + 1\right) - 2\right) \log{\left(\frac{\Lambda^{2} \rho}{2 \left(\delta \sqrt{\rho} \sigma \left(q + 1\right) - 2\right)^{2}} \right)} - 1.0}{\delta \sqrt{\rho} \sigma \left(q + 1\right) - 2} .
$$

##### S1 crystal

The total pressure is simply $\beta p/\rho=\beta p_1/\rho+\beta p_2/\rho$.

Likewise, the total Gibbs free energy reads $\beta G/N=\beta G_1/N+\beta G_2/N$.
By numerically inverting the pressure relation, we can express $G$ as a function of pressure.

### Crystals-only phase diagram

Using free energies of the two hexagonal crystals and that of S1, obtained above, we can perform a common tangent construction for various pressures and sketch a phase diagram with those crystals only, forgetting about the fluid for now.

The result is not super exciting. Lets see if the fluid brings in some more fun.

### Fluid free energy

We will compute the fluid free energy from thermodynamic integration of an equation of state for binary mixtures of disks.
In the additive case, we use the following formula provided by Santos *et al.* in [10.1080/00268979909482932].
$$
Z_{\text{eSHY}}(\eta)=\frac{\beta p}{\rho}=\frac{<\sigma>^2/<\sigma^2>}{1-2\eta+(2\eta_0-1)\eta^2/\eta_0^2}+\frac{1}{1-\eta}\left(1-\frac{<\sigma>^2}{<\sigma^2>}\right)
$$
$\eta$ is the packing fraction, $\eta_0=\sqrt{3}\pi/6$ is the packing fraction of the hexagonal close packing of monodisperse HD. $<\sigma^m>=\sum_{i=1}^{\text{nb of species}}x_i \sigma_i$ with $\sigma_i$ the diameter of disks of specie $i$ and $x_i$ their number fraction.

For a binary mixture with $q=0.4$ and $x_S=1/3$, we get the following curve.

We can then compute the Helmoltz free energy of the fluid as
$$
\frac{\beta F(\rho)}{N} = \frac{\beta F_{\text{id}}(\rho)}{N} + \int_{0}^\rho d\rho' \frac{\partial}{\partial \rho}\left(\frac{\beta F_{\text{fluid}}(\rho)}{N}-\frac{\beta F_{\text{id}}(\rho')}{N}\right)\\
=\ln(\rho \Lambda^2)+x_S\ln(x_S)+(1-x_S)\ln(1-x_S)-1+\int_0^{\rho}\left(\frac{\beta p}{\rho'}-1\right)\frac{d\rho'}{\rho'}
$$


Then, we use the equation of state again to compute Gibbs free energy $\beta G/N=\beta F/N + \beta p/\rho$. We can numerically invert the equation of state to obtain $G$ as a function of pressure or as a function of the composition $x_S$ for fixed pressure. The latter is used in common tangent construction.

### Adding fluid to the phase diagram

We can now perform common tangent constructions taking the fluid into account.
At very high pressure, we indeed recover the infinite pressure regime, where only hexagonal crystals, S1, and their coexistences, are found stable.

At low pressure, only the fluid is stable.
*Smile!*

Coexistences appear inbetween.

Scanning pressures, one gets the following phase diagram.

Large disks crystallise at $\beta p\sigma_{LL}^2\approx 13$ while small disks crystallise at $\beta p\sigma_{LL}^2\approx 81$, that corresponds to $\beta p \sigma_{SS}^2=\beta p\sigma_{LL}^2q^2\approx 13$.

At compositions of $x_S=1/3$, HexL+S1 coexistence becomes stable at $\beta p /\sigma_{LL}^2\approx 24$ which is about 15% smaller than the RT+fluid coexistence pressure of $28k_B T\sigma_{LL}^2$ found in simulation scan performed with slightly larger disks and a bit of non additivity (different system, and no proper treatment of RT, but that's the best comparison I have for the moment : scan of additive system demixes into HexL+fluid so far).

### Naive treatment of the RTQC

We want to distinguish the random tiling from the HexL+S1 coexistence. As a first, very rough approximation, let's assume that
$$
\frac{\beta G_{RT}}{N}=\frac{\beta G_{Hex_L+S1}}{N}-\frac{\beta TS_{RT}}{N}
$$
Pavel's paper [10.1088/0305-4470/27/11/010] provides $\sigma_v \doteq S_{RT}/N_v$ as a function of $\alpha_t$ the area fraction covered by the triangles. $\alpha_t$ is directly related to the composition $x_S$ in our system and $N_v$, the number of vertices, is the number of large disks. Therefore, $S_{RT}(\alpha_t)/N=\sigma_v(\alpha_t)(1-x_S)$.
Moreover, the entropy is concave and $G_{Hex_L+S1}$ is a straight line. Therefore only the point of maximal entropy has a chance to be stable. Thus, we treat the RTQC as a crystal point with composition $x_S=\sqrt{3}/(2+2\sqrt{3})\approx 0.317$ and Gibbs free energy $\beta G_{Hex_L+S1}(p)/N-\sigma_v^{max}(1-x_S)$.
The maximum value of the entropy per vertex is provided in Pavel's paper $\sigma_v^{max}/k_B=\ln(108)-2\sqrt{3}\ln{(2+\sqrt{3})}\approx 0.120$.

The phase diagram is modified as follows.

### Non-additive fluid

We use Santos equation of state (10.1063/1.1832591) for thermodynamic integration. It reads
$$
Z^{\text{SHY}}(\rho)=1+\frac{\eta}{1-\eta}\frac{b_3 \langle\sigma^2\rangle \overline{B_2} - b_2 \overline{B_3}}{(b_3-b_2)\langle \sigma^2\rangle^2} + (Z_{\text{pure}}(\eta)-1)\frac{\overline{B_3}-\langle \sigma^2\rangle\overline{B_2}}{(b_3-b_2)\langle\sigma^2\rangle^2}
$$
with $\eta$ the packing fraction of the mixture, $\langle \sigma^2\rangle=x_S \sigma_{SS}^2 + (1-x_S) \sigma_{LL}^2$, $b_2=2$ and $b_3=16/3-4\sqrt{3}/\pi$ the reduced virial coefficients of the monodisperse fluid of HD.
$\overline{B_n}$ are the reduced virial coefficients of the mixture. They can be expressed in terms of the composition-independant virial coefficients as
$$
\overline{B_2} = \overline{B_{SS}}x_S^2 + 2\overline{B_{SL}}x_S(1-x_S) + \overline{B_{LL}}(1-x_S)^2\\
\overline{B_3} = \overline{B_{SSS}}x_S^3 + 3\overline{B_{SSL}}x_S^2(1-x_S) + 3\overline{B_{SLL}}x_S(1-x_S)^2 +\overline{B_{LLL}}(1-x_S)^3
$$
The composition-independant virial coefficients are given by
$$
\overline{B_{ij}}=2\sigma_{ij}^2\\
\overline{B_{SSS}}=b_3\sigma_{SS}^4\\
\overline{B_{SSL}}=b_3\sigma_{SS}^4B_{\text{app}}\left(\frac{\sigma_{SL}}{\sigma_{SS}}\right)
$$
and permutation of S's and L's.
$B_{\text{app}}(s)=1/(3b_3)[4(s-1)^2+b_3(4s^2-1)]$ for all $s>=1/2$ (which is always valid for the small $\Delta$ that we consider here).
$Z_{\text{pure}}$ is a freely chosen equation of state for monodisperse hard disks. Here, we choose Henderson's one, wich is both simple and accurate.
$$
Z_{\text{pure}}=\frac{1+\eta^2/8}{(1-\eta)^2}
$$
In the expression above, it is evaluated at the packing fraction of the *mixture*.
For $\Delta=0$, it reproduces very well the additive equation of state.
Thermodynamic integration gives Helmoltz free energy.
Gibbs free energy is obtained by Legendre transform.

### Sanity check: additive case with non-additive expression

We check that the expression for non-additive disks reproduces the correct phase diagram in the additive case.

### Non-additive phase diagram

We can now compute the phase diagram for non-additive mixtures.
A close-up view of the transitions with the fluid
Note that for $q=0.4$ and $\Delta=-0.05$, we should have included other phases (T1). This is only a check of the procedure and tools.

#### Inside the RT region

A few points inside the RT region, on the line corresponding to spheres on a planar subtrate.





## Square-triangle random tilings

This section is devoted to analysis performed on square-triangle random tilings more specifically.
In the tilings, squares and triangles are colored according to their orientation. For squares, the orientation is defined as the orientation of any side, modulo $\pi/2$. For triangles, the orientation is defined as the orientation of the vector going from the middle of any side to the opposite vertex, modulo $2\pi/3$.
Appart from the colored tiling, I plotted the orientation distributions of squares and triangles in the last snapshot of the simulations.
Main directions and number of tiles in each orientation are determined using the k-mean clustering algorithm. The evolution of the number of tiles in each orientation over time, and the evolution of the main orientations aver time are also plotted.

In an square-triangle random tiling quasicrystal, the number of squares divided by the number of triangles is $N_{sq}/N_{tr}=\sqrt{3}/4\approx 0.433$. Moreover, there are 3 different orientations of squares, which appear in the same proportions, and 4 different orientations of triangles which also appear in the same proportions.

Note that when a specific orientation is too close to zero and the correponding groupe is wrapped around the angle interval, all orientations are shifted by a small angle so that each group gets assigned a proper color.

### First analysis of the equilibrated RT

* $q=0.45, \Delta=-0.05, N_L=2000, N_S=1000, \phi=0.836$ !
$N_{sq}/N_{tr}=820/2079\approx0.394$. There is an excess of triangles in this tiling compared to ideal one.
There seem to be two kinds of regions in the tiling: one with a big clusters of triangles separated by squares grain boundaries (eg. centre of the tiling) and the other in which squares and triangles are much more mixed, with a lot of sigma units (eg. below and right of the previous region).
The first kind reminds *triangular foams* corresponding to tilings with a lot of triangles. The second looks much closer to a RTQC.
Could we actually have a coexistence of random tilings ?

* $q=0.457, \Delta=-0.07, N_L=2000, N_S=1000, \phi=0.83$ 
$N_{sq}/N_{tr}\approx0.370$. The tiling differs from RTQC both by the ratio of tiles (excess of triangles) and by the orientation ditribution (excess of one pair of triangles). The excess of triangles is due to coexisting fluid capturing a lot of small particles, hence preventing the formation of squares. At larger packing fraction, the fluid part is reduced and this effect vanishes (see next tiling).

* $q=0.457, \Delta=-0.07, N_L=2000, N_S=1000, \phi=0.83$ 
$N_{sq}/N_{tr}\approx 0.429$ is pretty close the the RTQC value (0.433). However angular distributions are off : one pair of triangles dominates the other and squares are not evenly distributed either.
The total number of squares is slowly increasing and seems to not have reached equilibrium at the end of the simulation. Conversely the total number of triangles decreases very slowly.
The orientations of the various tiles are almost unchanged since the beginning of the simulation, and the orientations are the same as in the previous case, which used the same initial configuration at lower packing fraction. The RT is not actually formed from the melt, but grows from seeds already present in the initial configuration. To reach equilibrium, the fluid then has to come in and reorganise the tiling, which takes a lot of time.