---
tags: SoftQC, self-assembly
---
# Finite pressure phase diagram for binary NAHD mixtures

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

---

---

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

Frank, Giuseppe, Etienne

---

---

## Background

We computed phase diagrams at infinite pressure for binary NAHD mixtures, and found regions where dodecagonal random tiling quasicrystals are expected to form.

Using the cell approximation, and thermodynamic integration of the SHY equation of state for mixtures of NAHD, we computed finite pressure phase diagrams for those systems.


We now want to compute a more precise phase diagram by performing exact free energy calculations.

## Plans

The questions:
* Can we build a more precise phase diagram at finite pressure ? (We certainly can)
* Is the perfect square-triangle QC more stable than its random tiling counterparts ?
* Can we identify modes that contribute the most to the RTQC stability ?

To answer these question, precise free-energy calculations are required.

The SHY equation of state for NAHD fluids is only an approximation. We want to compute the fluid free energy from an equation of state obtained via NPT MC simulations instead.

The solids free energy will be computed exactly using Frenkel-Ladd method.

Since those calculations are a bit heavy, we focus on one specific size ratio and non-additivity parameter. We chose $q = 0.46$ and $\Delta = -0.0709$ that corresponds to spheres lying on a substrate.
We choose a size ratio on the lower end of the RTQC stability region, to favour the self-assembly (small particles escape squares more easily).


## Points to keep in mind

* Square-triangle tiling requires a very complicated fractal-shape acceptance region in perpendicular space. Simpler regions introduce other tiles in the tiling, including shields. To differentiate between a perfect QC with a finite amount of shields and a defective square-triangle RTQC, we should look at the diffraction pattern and see whether the Bragg peaks are broadened or not.
* Box equilibrium shape might change as pressure decreases. Remember to take that into account.
* Some phases that were not present in the infinite pressure phase diagram might be stable at finite pressure. Don't forget them.

## Fluid free energy

### Equation of state fit

[2021/02/22]

#### Simulated equation of state

We use the composition $x_S = 0.5$ to benchmark the tools we will then use for other compositions.

We measure the equation of state in NPT MC simulations with 500 particles of each type. In monodisperse systems, HD crystallise in simulation at $\beta p \sigma^2 \sim 9$. Since crystallisation of monodisperse systems of large and small disks must occur at the same reduced pressure, we run the simulations up to pressure $9 / q^2 \sim 45 \beta \sigma^2$.
The simulations are run in two steps. First the system is randomly initialised at low density, and $10^6$ MC sweeps are performed at the target pressure, while MC step amplitudes are adapted every $10^3$ MC sweeps to match target rejection rates. Then, the steps amplitudes are frozen, and $10^6$ MC sweeps are performed, during which the system's density is measured. Simulations even at the highest pressures show no sign of drift in the density during the production run.

We compare the simulated equation of state with the analytical expression provided by Santos et al. in J. Chem. Phys. **122**, 24514 (2005).


The SHY approximation is good at low densities, but quite bad above $\rho \sigma^2 \sim 1$, as could be expected.

We do not observe S1 crystallisation in the simulations even at the highest pressures. This could be because we only observe "supercompressed" liquids in our short, quickly compressed simulations.

#### Thermodynamic integration

To compute the free energy of the fluid, we integrate

$$
I(\rho) = \frac{1}{\rho} \left( \frac{\beta p}{\rho} - 1 \right) = \sum_{k=2}^\infty B_k \rho^{k-2}
$$

From the measured $\rho(p)$, the integrand $I(\rho)$ is computed and fitted with polynomials of higher and higher degrees.
To improve fitting at low density, I compute exactly the second Virial coefficient using expressions in Santos et al. paper.
$$
B_2 = \frac{\pi}{2} \left[ \sigma_{LL}^2 (1 - x_S)^2 + 2 \sigma_{LS}^2 x_S(1-x_S) + \sigma_{SS}^2 x_S^2 \right]
$$
$\sigma_{ij}$ is the minimal approach distance for disks of types $i$ and $j$. $x_S$ is the number fraction of small disks (here 0.5).

Santos et al. also provide an expression for the exact third coefficient, which is a bit more cumbersome, but also used in the fit procedure.
$$
B_3 = \left(\frac{\pi}{4}\right)^2 \left[ \overline{B}_{LLL} (1-x_S)^3 + 3 \overline{B}_{LSS}x_S^2(1-x_S) + 3 \overline{B}_{LLS} x_S (1-x_S)^2 + \overline{B}_{SSS} x_S^3 \right]
$$
with
$$
b_3 = \frac{16}{3} - \frac{4\sqrt{3}}{\pi} \\
\overline{B}_{iii} = b_3 \sigma_i^4 \\
\overline{B}_{iij} = b_3 \sigma_i^4 \mathcal{B}\left(\frac{\sigma_{ij}}{\sigma_{ii}}\right) \\
\mathcal{B}(s) = \frac{4}{3 \pi b_3} \left[ 4 \pi s^4 - 8 s^2 (s-1) \arccos\left(\frac{1}{2s}\right) - (2s^2 + 1) \sqrt{4s^2 - 1} \right] \quad \text{for} \quad s \geq \frac{1}{2}
$$

#### Least square fit

##### Polynomial fit degree

Fitting with up to the 20th Virial coefficient (17 free parameters) provides a good representation of the measured equation of state at high density.
However, since the sampling of the equation of state is performed with a fixed pressure step, the data points density is much higher at high particles density. The larger y-uncertainty at high density (taken into account by the LS fit) is not enough to compensate for this oversampling, which results in a poor fit a low density.

>[name=Frank] Agree. One thing to keep in mind for later is that it can be helpful later to reduce the regime of data you include. If it turns out later that the freezing point at this composition is at a density of e.g. $\rho \sigma^2 = 1.50$, then dropping the data beyond, say, $\rho \sigma^2 = 1.60$ is quite reasonable. This is for two reasons: 1) you may have trouble properly equilibrating the systems at densities much beyond freezing (glassy dynamics, risk of forming crystal regions), and 2) since the values of the function-to-be-fitted are high at high densities, this means that the fitting procedure is also inclined to spend a lot of effort fitting them accurately (i.e. a 10% error at high density "costs" more than one at low density in terms of the squared error which is minimized). Hence, it can be better not to include extreme values unless necessary.


##### Sampling accuracy

The following figure investigates the effect of the number of points in the data set. One point every $n$ is taken into account ($n \in [1, 5]$)
The more points are removed, the more oscillations propagate to higher densities.
No improvement in the fit is gained by pruning the data set uniformly.

>[name=Frank] This is indeed a nice test. (In principle, using a bigger pressure step at higher pressures would also be reasonable.)

> [name=gfoffi] I am slightly confused by the picture below. Are you not using analytic B2 and B3? In this case, I would expect the low density limit to be the same for all the curves no?


> [name=Etienne] The poor fitting at low density also surprised me. The second and third Virial coefficient only dictate the constant and linear term in the integrand. When taking a closer look, it seems that all curves indeed have the same slope at the origin, but then quickly deviate, probably because more weight is put on the high density part.


[04/03/21]
In order to avoid giving an unreasonable weight to the high density part of the integrand, I pruned the data set to impose a minimal $\rho$-spacing $\delta \rho = 0.01\sigma^{-2}$ between consecutive data points.
The resulting fit (up to 20th Virial coefficient) displays no weird oscillations at low density, and is very close to the fit to the full data set at high density.


This solution however in not very satisfying since it introduces an extra parameter $\delta \rho$ to be tuned. More importantly, it does not take advantage of the whole information present in the simulated data.

##### Fitting domain

Cutting out points with pressure higher than $35 k_B T \sigma^{-2}$ yields much better fits.

However, we do not know in advance the maximum relevant pressure. Therefore, it would be nice to have a fitting method that would be robust enough to compute a first reasonable estimate of the phase diagram, from which maximum pressures at every composition could be inferred. If necessary, this procedure could be repeated interatively until convergence.

#### Orthogonal distance regression fit

[04/03/21]

Least square (LS) fitting minimises the distance between the points and the fitted curve in the $y$ direction. Therefore, it can take $y$-uncertainties into account only.
In our case, we measure $\rho$ from simulations, which gives an uncertainty in both $x$ and $y$ directions when we want to fit the integrand.

At high density, the uncertainty on $\rho$ is actually much larger than the uncertainty on the computed integrand. Therefore weighting the data points with $y$-uncertainty only does not balance the much finer sampling at high density.


Orthogonal distance regression (ODR) is a fitting method that takes both $x$ and $y$ uncertainties into account. It yields better results than LS fit in our case.

When the full data set is used, the ODR fit is more robust at low densities than the LS fit. As a consequence, fitting virial coefficients up to $B_{15}$ already gives better results at low density than the LS fit with virial coefficients up to $B_{20}$.


When the data set is restricted to pressures below crystallisation, the two fitting methods yield similar results.

ODR appears more robust than LS. It allows to reduce the number of fitted parameters, increasing stability.
Therefore, I will use ODR to fit the equation of state up to $B_{15}$ (12 free parameters) in the following.

>[name=Frank] This seems like a good plan. From a physics point of view, it is slightly arbitrary since now the choice of units assigns a scale ratio to the horizontal and vertical axes that specifies how much we 'mind' an error in that direction. Since the units here result in reasonable numbers that seems fine, though.

### Free energy calculation

The Gibbs free-energy per particle is computed from the measured equation of state as

$$
\frac{\beta G(N, x_S, p)}{N} = \ln\left(\rho (p) \Lambda_T^2\right) - 1+ x_S \ln(x_S) + (1 - x_S) \ln(1 - x_S)+ \frac{\beta p}{\rho} + \int_{0}^{\rho(p)}\left(\frac{\beta p}{\rho'} - 1\right)\frac{1}{\rho'} d\rho'
$$

We can choose the particles masses such that $\Lambda_T = 1$ without affecting the stability properties.

The mixing entropy term $x_S \ln(x_S) + (1 - x_S) \ln(1 - x_S)$ has infinite derivatives at $x_S = 0$ or $1$ and needs to be accurately computed for a reliable estimate of the common tangent. Therefore, it is better to add it after fits or interpolations. This is not very complicated since it only depends on composition.

Using the measured equation of state for fluids with 500 particles of each type, we obtain the following free energy:

##### Finite size effects

To assess the magnitude of finite size effects, I did the same calculation on an equation of state simulated with 100, 250 and 1000 particles of each type.

The following figure displays the free energy difference with respect to the system with 1000 particles of each type.


There is no clear convergence towards the "large" system's free energy. In particular, the free energy computed from systems of 250 particles of each type is closer to the line obtained with 1000 particles of each type than that of the systems with 500 particles of each type.
I have no explanation for the oscillation around $p = 1$.

> [name=frank] My guess is that you mostly get some variation in actual pressure values for high pressures ($>30$). Your differences in between are pretty flat, meaning that between pressures of 5 and 30 your data is quite consistent... The wobbles at low pressure could be the result of the fit trying to deal with strong behavior at high pressures and neglecting to pay attention to the accuracy at low pressures, where values are smaller. Perhaps trying to fit only up to pressure 30 and seeing if it improves is worthwile. (as noted above, it is not unlikely that the fluid is no longer stable from some point on...).

> [name=gfoffi] I guess that you could try to vary the upper value of the interval on which you fit the pressure. Maybe this help to see if the higher part of the pressure is "dragging" your fit. As you said below, the appearance of the S1 phase is probably reducing the accuracy of the fit.

[02/03/2021]
Ignoring points above pressure $35 k_B T \sigma^2$ does not change much the behaviour of the curves.


[08/03/2021]
I added more points at low density and fitted with a polynomial of order 9 instead of 13.

The sampling in $p$ is the following:
* for $\beta p \sigma^2 \in [0.05, 1]$, $\beta \delta p \sigma^2 = 0.05$
* for $\beta p \sigma^2 \in [1, 3]$, $\beta \delta p \sigma^2 = 0.1$
* for $\beta p \sigma^2 > 3$, $\beta \delta p \sigma^2 = 0.2$


The free energy finite size scaling is more consistent.

Finite size error on the Gibbs free-energy is of the order of $10^{-2} k_B T$.

##### Fit degree

The following figure displays the free energy difference obtained on the same simulations data (500 particles of each type, improved sampling at low density) for various degrees of the fitting polynomial.


We choose to use a polynomial of order 9. The uncertainty in the Gibbs free-energy stemming from the fit order choice is of the order of $10^{-3}k_BT$.

##### Simulation length

To make sure that simulations were long enough, I measured again the equation of state of systems with 1000 particles of each type, with 10 times longer simulations.

With more time, S1 solids start appearing in the simulations. A shy kink appears in the equation of state at pressures of about $35 \beta \sigma^2$

>[name=Frank] Again a good test. This also means that you should not fit the data for pressures higher than where you see spontaneous crystallization, especially for the longer simulations. You will not need the fluid free energy there (since it is unstable), and fitting through a kink will always lead to errors.

The following figure displays the Gibbs free-energy difference with respect to the longer simulations.

The larger difference at high pressures is due to S1 nucleation in the longer simulations.
At the pressure of the fluid, the free-energy difference is of the order of $10^{-3} k_B T$.

[02/03/2021]
Excluding pressures above $35 k_B T \sigma^2$ from the fit further improves the improves the agreement between the two curves.

[08/03/2021]
Same as above, with a 9th order polynomial fit and more points at low densities.

I think I'll keep the shorter simulations and simulate systems with 500 particles of each types to compute the free energy at many different compositions.
In this setting, each simulation takes about 3 to 5 minutes. At each composition 261 data points are taken at various pressures.

>[name=Frank] Seems very reasonable!

> [name=gfoffi] Agree! Nice work indeed!

[27/09/2021]

The custom sampling scheme presented above works fine, but the domains boundaries as well as the sampling density within each domain is quite arbitrary. I've tried logarithmic sampling instead.

The fit with a polynomial of degree 9 looks good:


Taking the free energy computed from the EOS sampled with the previous scheme and $10^7$ MC sweeps in the production run as reference (blue curve), the next figure shows that the logarithmic sampling of the EOS with $10^6$ MC sweeps (gree) yields results similar to the previous sampling with $10^6$ MC sweeps (orange).


### Fluid Gibb-free-energy
[29/09/2021]

I computed the equation of state of mixtures of 2000 NAHD for various compositions and pressures, using NPT MC.
The equilibration and production runs are both $10^6$ MC sweeps long.
250 simulations are run for pressures geometrically distributed in $[0.1, p_{max}]$. The maximum pressure depends on the composition. Monodisperse hard disks freeze at a reduced pressure of about 9. Thus our systems will freeze at $\beta p \sigma_{LL} ~ 9$ when $x_S = 0$ and at $\beta p \sigma_{LL} ~ 45$ when $x_S = 1$. At intermediate composition, the freezing pressure is unknown, so we set the maximum pressure at $\beta p_{max} \sigma_{LL} = 35 x_S + 15$ which linearly interpolates between two points located above the monodisperse freezing pressures.
Compositions are sampled with a step size of $\delta x_S = 100/2000 = 0.05$.

The simulated equations of state at each composition are fitted with 12 Virial-like coefficients\footnote{This is just a way to report the number of fitting parameters. The values of the fitted coefficient do not correspond to real Virial coefficient, since the fit parameters can compensate each other in many ways.}, inverted and integrated to compute Gibbs free energies.

The following figure displays the Gibbs free-energy of the NAHD mixture.


Finally, the free-energies at constant pressure, are interpolated to obtain the free-energy as a function of composition that will be used in common tangent construction.
The following figure shows the result of this interpolation at reduced pressure of 5, before and after the inclusion of the mixing entropy term.


The curve is smooth, so the sampling in $x_S$ seems good for now.

## Frenkel-Ladd method

To compute the free energies of the solid phases, we will use Frenkel-Ladd method [D Frenkel and AJC Ladd, J. Chem. Phys. **81**, 3188 (1984)] with the Einstein molecule approach [C Vega and EG Noya, J. Chem. Phys. **127**, 154113 (2007)] [C Vega, E Sanz, JLF Abascal and EG Noya, J. Phys. Condens. Matter **20**, 153101 (2008)].

The general idea is to start from a reference system whose free energy is known, and then gradually transform it into the system of interest by a combination of explicit free-energy differences calculations and Hamiltonian integration.

The reference system is the *ideal Einstein molecule* in which all particles but one are tight to their lattice site by a harmonic spring. The particle that is not attached to a spring (say particle 1, without loss of generality) is called the carrier and defines the lattice.
$$
U_\text{ein-mol-id}(\vec{r_1}, \vec{r_2}, \ldots, \vec{r_N}) = \gamma \sum_{i=2}^N (\vec{r_i} - \vec{r_{i_0}}(\vec{r_1}))^2
$$
Notice that the sum starts at 2, so the carrier can move freely. The lattice sites postions are functions of the position of the carrier since the lattice moves with it. $\gamma$ is the stiffness of the springs.

The ideal Einstein molecule is first transformed into an interacting Einstein molecule with hard repulsion between particles added on top of the springs. Then, the springs are gradually turned off to transform the system into a pure hard disk mixture.

When the springs stiffness goes to zero, the whole solid can drift, leading to a quasi-divergence of the intergrand in the Hamiltonian integration. This can be avoided by fixing the position of one of the particles in the system. In the Einstein molecule, this is easily done by fixing the position of the carrier, which is the reason for choosing such a strange reference system to start with. Fixing a particle's position comes with a free energy contribution that needs to be accounted for, but allows to keep the numerical method stable.

Overall, the transformation path is the following:
* Start from an *ideal Einstein molecule* at the desired packing fraction. Its free energy $F_\text{ein-mol-id}$ can be computed analytically.
* Transform it into an *ideal Einstein molecule with one particle fixed*. The free-energy difference due to fixing the carrier can be computed analytically, and turns out to not play a role. The star superscript indicates that one particle is fixed in the system.
* Add the hard core interaction on top of the springs to get an *interacting Einstein molecule*. The resulting free energy $F_\text{ein-sol}^*$ (the sol subscript stands for *solid*) can be obtained from a average calculation in a simulation.
* Gradually turn off the springs and perform a Hamiltonian integration to get the free energy $F_\text{sol}^*$ of the solid with one particle fixed.
* Free the fixed particle to get the free energy $F_\text{sol}$ of system of interest.

After this procedure, we get the Helmoltz free-energy of the solid at a given packing fraction.
The Gibbs free-energy can be computed from the Helmoltz one by adding $\beta p / \rho$ that could be obtained from a simulated equation of state.

>[name=Etienne] Maybe it is possible to adapt the Frenkel-Ladd method to Gibbs free energies by using NpT simulations instead of NVT. I will see if I can find equations for this.
>> The Einstein molecule pressure is ill-defined since particles don't interact. Therefore, running NpT simulations on such system is going to be troublesome. Simply computing the equation of state seems a much better solution.

In the following, we derive the equations for the free-energy calculation of a binary mixture with $N_1$ particles of type 1 and $N_2$ particles of type 2. We assume that the carrier is a particle of type 1.

### Free energy of the ideal Einstein molecule

After integration over the momentum degrees of freedom, the partition function of the ideal Einstein molecule reads

$$
Z_\text{ein-mol-id} = \frac{1}{N_1!} \frac{1}{N_2!} \frac{1}{\Lambda_T^{2N}} \int{ d\vec{r_1} \ldots d\vec{r_N} e^{-\beta U_\text{ein-mol-id}}}\\
= \frac{1}{N_1!} \frac{1}{N_2!} \frac{1}{\Lambda_T^{2N}} \int{ d\vec{r_1}} \int{ d\vec{r_2}\ldots d\vec{r_N} e^{-\beta U_\text{ein-mol-id}}}
$$

The second integral is independant of the value of $\vec{r_1}$, ie. of the position of the lattice, because the integral runs over the whole space anyway. Therefore, the first integral over $\vec{r_1}$ can be performed and yields a volume factor.

$$
Z_\text{ein-mol-id} = \frac{1}{N_1!} \frac{1}{N_2!} \frac{1}{\Lambda_T^{2N}} V \int{ d\vec{r_2}\ldots d\vec{r_N} e^{-\beta U_\text{ein-mol-id}}}
$$

For the remaining configurational integral, one needs to take into account all possible permutations of the particles on the lattice sites.

$$
Z_\text{ein-mol-id} = \frac{(N_1 - 1)!}{N_1!} \frac{N_2!}{N_2!} \frac{1}{\Lambda_T^{2N}} V \int_\text{one perm}{d\vec{r_2} \ldots d\vec{r_N} e^{-\beta \gamma \sum_{i=2}^N (\vec{r_i} - \vec{r_{i_0}})^2}}\\
= \frac{1}{\Lambda_T^{2N}} \frac{V}{N_1} \left(\int{dr e^{-\beta \gamma r^2}}\right)^{2(N-1)}\\
= \frac{V}{\Lambda_T^{2}} \frac{1}{N_1}\left(\frac{\pi}{\beta \gamma \Lambda_T^2}\right)^{N-1}
$$

The corresponding Helmoltz free energy is
$$
\frac{\beta F_\text{ein-mol-id}}{N} = -\frac{1}{N} \ln(Z_\text{ein-mol-id})\\
= \frac{1}{N}\ln\left(\frac{\Lambda_T^2 N_1}{V}\right) + \left(1 - \frac{1}{N}\right)\ln\left(\frac{\Lambda_T^2\beta\gamma}{\pi}\right)
$$

### Fixing one particle

To avoid quasi-divergence in the Hamiltonian integration, we need to prevent drifts of the whole solid. This can be achieved by fixing the position of one of the particles in the system.

Consider a **general** mixture of two types of particles. After integration over the momentum degrees of freedom, the partition function reads

$$
Z = \frac{1}{N_1!} \frac{1}{N_2!} \frac{1}{\Lambda_T^{2N}} \int{ d\vec{r_1} \ldots d\vec{r_N} e^{-\beta U(\vec{r_1}, \vec{r_2},\ldots, \vec{r_N})}}
$$

If the potential energy is invariant under the stranslation of all the particles, we can write

$$
Z = \frac{1}{N_1!} \frac{1}{N_2!} \frac{1}{\Lambda_T^{2N}} \int{ d\vec{r_1}} \int{ d\vec{r_2} \ldots d\vec{r_N} e^{-\beta U(\vec{0}, \vec{r_2}-\vec{r_1},\ldots, \vec{r_N}-\vec{r_1})}}\\
= \frac{1}{N_1!} \frac{1}{N_2!} \frac{1}{\Lambda_T^{2N}} \int{ d\vec{r_1}} \int{ d\vec{r_2'} \ldots d\vec{r_N'} e^{-\beta U(\vec{0}, \vec{r_2'},\ldots, \vec{r_N'})}}
$$

The second integral does not depend on $\vec{r_1}$ anymore, so the first integral can be performed and simply yields a volume factor.

$$
Z = \frac{1}{N_1!} \frac{1}{N_2!} \frac{V}{\Lambda_T^{2N}} \int{ d\vec{r_2'} \ldots d\vec{r_N'} e^{-\beta U(\vec{0}, \vec{r_2'},\ldots, \vec{r_N'})}} = \frac{V}{\Lambda_T^2} Z^* \\
Z^* = \frac{1}{\Lambda_T^{2(N-1)}} \frac{1}{N_1!}\frac{1}{N_2!}\int{ d\vec{r_2} \ldots d\vec{r_N} e^{-\beta U(\vec{0}, \vec{r_2},\ldots, \vec{r_N})}}
$$

$Z^*$ is the partition function of the system with particle 1 fixed.
The free-energy can then be written
$$
\frac{\beta F}{N} = \frac{1}{N} \ln\frac{\Lambda_T^2}{V} - \frac{1}{N} \ln Z^*\\
= \frac{1}{N} \ln \frac{\Lambda_T^2}{V} + \frac{\beta F^*}{N}
$$

It is quite miraculous that the free-energy of a system with one particle fixed and the free-energy of the full system only differ by a constant term that is independant of the potential.

The above derivation relies on the fact that the potential energy is invariant under a translation of the whole system. This assumption obviously holds for a hard disks mixture. It is also valid for the Einstein molecule because the external field applied by the springs is tethered to the carrier. Hence, when the whole system is translated (including the carrier), the lattice sites are translated by the same amount and the potential energy is unchanged.

This means that the free energy term $\frac{1}{N} \ln \frac{\Lambda_T^2}{V}$ coming from fixing the position of the carrier exactly cancels with the term coming from freeing the position of particle 1 in the real solid. Hence, it is enough to add the free-energy difference between the ideal Einstein molecule with one particle fixed and the solid with one particle fixed directly to the free-energy of the ideal Einstein molecule (with no particle fixed) to get the free-energy of the solid (with no particle fixed).

$$
F_\text{sol} = F_\text{ein-mol-id} + \Delta F_1^* + \Delta F_2^*
$$

with $\Delta F_1^*$ the free-energy difference between the ideal Einstein molecule with one particle fixed and the interacting Einstein molecule with one particle fixed, and $\Delta F_2^*$ the free energy difference between the interacting Einstein molecule with one particle fixed and the solid with one particle fixed.

### $\Delta F_1^*$

We want to compute the free-energy difference between the ideal Einstein molecule with one particle fixed and the interacting Einstein molecule with one particle fixed.

The energy of the interacting Einstein molecule with one particle fixed is given by $U_\text{ein-sol}^* = U_\text{sol}^* + U_\text{ein-mol-id}^*$.

$$
\Delta F_1^* = F_\text{ein-sol}^* - F_\text{ein-mol-id}^*\\
= -k_BT \ln\left(\frac{\int{d\vec{r_1}\ldots d\vec{r_N} e^{-\beta U_\text{ein-sol}^*}}}{\int{d\vec{r_1}\ldots d\vec{r_N} e^{-\beta U_\text{ein-mol-id}^*}}}\right)\\
= -k_BT \ln\left(\frac{\int{d\vec{r_1}\ldots d\vec{r_N} e^{-\beta (U_\text{ein-sol}^* - U_\text{ein-mol-id}^*)} e^{-\beta U_\text{ein-mol-id}^*}}}{\int{d\vec{r_1}\ldots d\vec{r_N} e^{-\beta U_\text{ein-mol-id}^*}}}\right)\\
= -k_BT \ln(\langle e^{-\beta (U_\text{ein-sol}^* - U_\text{ein-mol-id}^*)} \rangle_{\text{ein-mol-id}^*} )\\
= -k_BT \ln(\langle e^{-\beta U_\text{sol}^*} \rangle_{\text{ein-mol-id}^*} )
$$

The average in the logarithm is computed using a NVT simulation of the ideal Einstein molecule with one particle fixed.

The exponential term takes only value 0 (when at least one overlap is present) or 1 (when the configuration has no overlap). Thus, the average term takes values in $[0, 1]$ and $\Delta F_1^*$ is positive.
This makes sense since turning on hard core repulsion on top of the harmonic springs reduces the entropy per particle, leading to a higher free-energy.

### $\Delta F_2^*$

The free-energy difference between the interacting Einstein molecule with one particle fixed and the solid of interest with one particle fixed is computed by Hamiltonian integration.

We consider a system interacting with the potential energy

$$
U_\lambda^* = \lambda U_\text{sol}^* + (1-\lambda) U_\text{ein-sol}^* = U_\text{sol}^* + (1 - \lambda) U_\text{ein-mol-id}^*
$$

that interpolates between the two systems as $\lambda$ varies from 0 to 1.
The derivative of the free-energy of this system reads

$$
\frac{\partial F_\lambda^*}{\partial \lambda} = -k_BT \frac{\partial Z_\lambda^*}{\partial \lambda} \frac{1}{Z_\lambda^*}
$$

and

$$
\frac{\partial Z_\lambda^*}{\partial \lambda} = -\beta \frac{1}{\Lambda_T^{2(N-1)}} \frac{1}{N_1!} \frac{1}{N_2!} \int{ d\vec{r_2} \ldots d\vec{r_N}(- U_\text{ein-mol-id}^*) e^{-\beta U_\lambda^*}}
$$

Hence

$$
\frac{\partial F_\lambda^*}{\partial \lambda} = - \langle U_\text{ein-mol-id}^* \rangle_{(1-\lambda)\gamma_\text{max}}
$$

The average value of the energy of the ideal Einstein molecule with one particle fixed and spring stiffness $\gamma_\text{max}$ is computed over a composite system with one particle fixed and parameter $\lambda$ (*ie.* spring stiffness $(1 - \lambda) \gamma_\text{max}$ + hard repulsion ).

The free-energy difference between the interacting molecule with one particle fixed and the solid with one particle fixed is obtained as

$$
\Delta F_2^* = - \int_0^1{d\lambda \langle U_\text{ein-mol-id}^* \rangle_{(1-\lambda)\gamma_\text{max}}}
$$

Change variable $\lambda$ for $\gamma = (1 - \lambda)\gamma_\text{max}$, the effective spring constant of the composite system.

$$
\Delta F_2^* = - \int_0^{\gamma_\text{max}} {\frac{d\gamma}{\gamma_\text{max}} \langle U_\text{ein-mol-id}^* \rangle_\gamma} \\
= - \int_0^{\gamma_\text{max}}{d\gamma \left\langle \sum_{i=2}^N (r_i - r_{i_0})^2 \right\rangle_\gamma}\\
= - \int_0^{\gamma_\text{max}}{\frac{d\gamma }{\gamma} \left\langle \gamma \sum_{i=2}^N (r_i - r_{i_0})^2 \right\rangle_\gamma}\\
= - \int_0^{\gamma_\text{max}}{ \frac{d\gamma}{\gamma} \left\langle U^* \right\rangle_\gamma}\\
$$

Since the hard core energy is always zero, the integrand can be obtained from the average energy of a composite system with spring constant $\gamma$. The average energy can be computed with NVT MC simulations.
Changes of variable can be performed to smooth the integrand further if necessary.
>[name=Etienne] I will first compute the integrand as is, and see how it looks before tweaking it further.

### Putting everything together

The calculation of the solid Helmoltz-free-energy by the Frenkel-Ladd method with Einstein molecule approach boils down to the following formula

$$
\frac{\beta F_\text{sol}}{N} = \frac{\beta F_\text{ein-mol-id}}{N} + \frac{\beta \Delta F_1^*}{N} + \frac{\beta \Delta F_2^*}{N}\\
= \frac{1}{N}\ln\left(\frac{\Lambda_T^2 N_1}{V}\right) + \left(1 - \frac{1}{N}\right)\ln\left(\frac{\Lambda_T^2\beta\gamma}{\pi}\right) -\frac{1}{N} \ln(\langle e^{-\beta U_\text{sol}^*} \rangle_{\text{ein-mol-id}^*} ) - \frac{\beta}{N} \int_0^{\gamma_\text{max}}{ \frac{d\gamma}{\gamma} \left\langle U^* \right\rangle_\gamma}
$$

### Discussion

Defects typically don't appear in simulated systems prepared in a perfect crystalline state. In particular, vacancy-interstitial defects pairs that could form and then diffuse are usually suppressed by the spring energy during Hamiltonian integration (even for low string stiffness). Therefore, the Frenkel-Ladd method allows us to compute the free-energy of perfect crystals.

Defects can be taken into account by inputing them in the starting structure. Tighting particles to the neirest lattice site instead of one specific site allows defects to diffuse and many configurations to be sampled in one simulation. A proper study of the defects impact on phases stability is a work in itself since, for instance, the equilibrium concentration of defects varies with pressure.

### Implementation of $\Delta F_1^*$ calculation
[10/04/2021]

$$
\Delta F_1^*= -k_BT \ln(\langle e^{-\beta U_\text{sol}^*} \rangle_{\text{ein-mol-id}^*} )
$$

To compute $\Delta F_1^*$, we need to sample equilibrium configuration of the Einstein molecule with one particle fixed, and compute the average of $\exp(-\beta U_\text{sol}^*)$ over those configurations.

The exponential term takes only 2 values: 1 if the configuration of the Einstein molecule results in a hard disk configuration with no overlap, and 0 otherwise.

The equilibrium configurations of the Einstein molecule with one particle fixed can be sampled by NVT MC simulation. Many simplifications arise in the code since particles do no interact with each other.

The probability of finding a particle at a position $\vec{r_i}$ is proportionnal to $\exp(-\beta \gamma (\vec{r_i} - \vec{r_{i_0}})^2)$. Hence, for a given spring stiffness $\gamma$, a new equilibrium configuration of the Einstein molecule can be generated by simply drawing new positions for all particles (except the carrier) from a Gaussian distribution with standard deviation $(2 \beta \gamma)^{-1/2}$.

I first implemented the MC version by modifying my previous MC code, and then implemented the Gaussian sampling.

#### Validation

For validation purposes, I consider an ideal Einstein molecule in a S1 configuration with 100 particles of each type.

The main check is the equipartition theorem, that predicts an average energy of $199 k_BT$ since the carrier is fixed.

For the MC sampling, the system is equilibrated for 1000 sweeps. Then, $10^4$ configurations are dumped with a spacing of 10 sweeps.
The measured average energy is $198.8 \pm 0.2 k_BT$, which is consistent with the equipartition theorem. The standard deviation is $13.9 k_BT$. Moreover, the Flyvbjerg-Peterson analysis displayed below indicates that the sampled configurations are uncorrelated because no increase in the error occurs when blocking steps are performed.


For comparison, I sampled $10^4$ configurations using Gaussian random variables, and obtained consistent results. The average energy is $198.9 \pm 0.2 k_BT$ and the standard deviation is $14.1 k_BT$. With this procedure, the sampled configurations are uncorrelated by construction, as confirmed by FP analysis.


#### Performances

The Gaussian sampling method is about 10 times faster than the MC simulation (order of minutes vs less than a second).
More importantly, the direct sampling produces uncorrelated configurations, so there is no need to worry about having a long enough equilibration or dumping configurations too frequently.

Therefore, I will use direct sampling in the following.

#### Choosing the spring stiffness

When the springs are too loose, all sampled configurations contain overlap. Therefore, $\langle \exp(-\beta U_\text{sol}) \rangle = 0$ and the free energy change from the ideal Einstein molecule with one particle fixed and the hard disk solid $\Delta F_1^* = -k_BT \ln(\langle e^{-\beta U_\text{sol}^*} \rangle_{\text{ein-mol-id}^*} )$ is infinite.
When the springs are super-tight, no overlap can occur and $\Delta F_1^*$ drops to 0, but the calculation of $\Delta F_2^*$ might become less accurate because the integral support becomes larger.

At this point, it is hard for me to grasp the balance: maybe taking a super high spring stiffness does not impact much the accuracy of the latter calculations.
In any case, it is interesting to have an idea of how the average term varies as the spring stiffness is changed.

I used a S1 solid as a test case, with a unit cell length of $1.2 \sigma_{LL}$. This corresponds to a fairly low packing fraction of 0.66.


The transition between an always overlapping system and a never overlapping system spans about $200\,(\beta \sigma^2)^{-1}$. When the number of particles in the solid increases, the transition is shifted to larger stiffnesses, but does not seem to become sharper.

The spring stiffness is directly related to the standard deviation of the displacement of the particles around their lattice site $\text{std} = (2 \beta \gamma)^{-1/2}$.
The two vertical dashed lines highlight the spring stiffness that correspond displacement of half the closest distance between large disks, or half of the closest distance between a large and a small particle in the perfect solid.
They do not provide a satisfactory choice for the spring stiffness.

We propose to choose the initial spring stiffness such that $\langle \exp(-\beta U_\text{sol})\rangle_\gamma = 0.9$. Considering the speed of the calculations for each point, this value can easily be determined with a bisection for instance.

#### Small side puzzle

[23/04/2021]

Can we understand the scaling of $\Delta F_1^*$ with the number of particles ?

Let $\mathcal{P}_o^{(N)}$ be the probability of drawing a configuration with at least one overlap in a Einstein molecule of $N$ particles.
Let $\mathcal{P}_{n-o}^{(N)}$ be the probability of drawing a configuration with no overlap in a Einstein molecule of $N$ particles.

Then
$$
\langle e^{-\beta U_\text{sol}} \rangle_\gamma^{(N)} = 0 \times \mathcal{P}_o^{(N)} + 1 \times \mathcal{P}_{n-o}^{(N)}
$$
is the the probability that a randomly chosen configuration is non-overlapping.

Let's then consider a system of $4N$ particles.

$$
\langle e^{-\beta U_\text{sol}} \rangle_\gamma^{(4N)} = 0 \times \mathcal{P}_o^{(4N)} + 1 \times \mathcal{P}_{n-o}^{(4N)}
$$

Since the particles are tighted to their lattice sites, the box with $4N$ particles can be seen a 4 sub-boxes of $N$ particles.

The probability of sampling a non-overlapping configuration is the probability of sampling non-overlapping configurations in each of the 4 sub-boxes plus a correction term.
$$
\langle e^{-\beta U_\text{sol}} \rangle_\gamma^{(4N)} = \left(\mathcal{P}_\text{n-o}^{(N)}\right)^4 + \mathcal{P}_\text{bound}
$$

$\mathcal{P}_\text{bound}$ is the correction term accounting for the possible overlaps generated at the boundaries between boxes, as well as the possibility for the number of particles to fluctuate inside each box.
Since moves far away from the lattice sites are exponentially rare, most of the contribution for this correction come from the particles at the boundary of the boxes.
This correction should become small for high spring stiffness and large systems.

Therefore we expect $\langle e^{-\beta U_\text{sol}} \rangle_\gamma^{(4N)} \sim \left(\langle e^{-\beta U_\text{sol}} \rangle_\gamma^{(N)}\right)^4$ at large $\gamma$.

The next figures shows this scaling in simulations of Einstein molecules with $200$ and $800$ particles.
The scaling is actually surprinsingly good all the way down to the lowest spring stiffnesses.

### Implementation of $\Delta F_2^*$ calculation
[29/04/2021]

$$
\Delta F_2^* = - \int_0^1{d\lambda \langle U_\text{ein-mol-id}^* \rangle_\lambda}\\
= - \frac{1}{\gamma_\text{max}} \int_0^{\gamma_\text{max}}{d\gamma \langle U_\text{ein-mol-id}^* \rangle_\lambda}
$$

We perform NVT MC simulations of the interacting Einstein molecule and compute the average value of the springs energy.

As before, for testing purposes, I used a S1 solid with a cell edge length of $1.2 \sigma_{LL}$.
The following figure displays the average spring energy per particle in solids of 200 and 800 particles.


For large spring stiffnesses, the average energy per particle goes to $(N-1) / N$ as expected from the equipartition theorem, since the carrier of the lattice is fixed.
Finite size effects seem to be tiny.

The integrand is computed from the average energy.

$\Delta F_2^*$ is obtained by integration of the simulated points using Simpson's rule.

### First full free energy calculation

Using the test data above, we can compute the Helmoltz free energy of the S1 solid at a packing fraction of $\phi = 0.661 \sigma_{LL}^{-2}$, for two different system sizes.

#### $N = 200$, $\gamma_\text{max} = 450 k_BT/\sigma_{LL}^2$

$$
\frac{\beta F_\text{id}}{N} = \frac{1}{N}\ln\left(\frac{\Lambda_T^2}{V}\right) + \frac{1}{N}\ln(N_1) + \left(1 - \frac{1}{N}\right)\ln\left(\frac{\Lambda_T^2\beta\gamma}{\pi}\right) \\
=\frac{1}{200}\ln\left(\frac{1}{144}\right) + \frac{1}{200}\ln(100) + \left(1 - \frac{1}{200}\right)\ln\left(\frac{450}{\pi}\right) \\
= 4.93787
$$

$$
\frac{\beta \Delta F_1^*}{N} = \frac{0.055618}{200} = 0.00027809
$$

$$
\frac{\beta \Delta F_2^*}{N} = -2.653425
$$

$$
\frac{\beta F}{N} = \frac{\beta F_\text{id}}{N} + \frac{\beta \Delta F_1^*}{N} + \frac{\beta \Delta F_2^*}{N} = 2.28468
$$

#### $N = 800$, $\gamma_\text{max} = 550 k_BT/\sigma_{LL}^2$

$$
\frac{\beta F_\text{id}}{N} = \frac{1}{N}\ln\left(\frac{\Lambda_T^2}{V}\right) + \frac{1}{N}\ln(N_1) + \left(1 - \frac{1}{N}\right)\ln\left(\frac{\Lambda_T^2\beta\gamma}{\pi}\right) \\
=\frac{1}{800}\ln\left(\frac{1}{576}\right) + \frac{1}{800}\ln(400) + \left(1 - \frac{1}{800}\right)\ln\left(\frac{550}{\pi}\right) \\
=5.158276
$$

$$
\frac{\beta \Delta F_1^*}{N} = \frac{0.042281}{800} = 5.285125 \times 10^{-5}
$$

$$
\frac{\beta \Delta F_2^*}{N} = -2.861572
$$

$$
\frac{\beta F}{N} = \frac{\beta F_\text{id}}{N} + \frac{\beta \Delta F_1^*}{N} + \frac{\beta \Delta F_2^*}{N} = 2.296757
$$

The difference is comparable to that found by Vega *et al.* for the free energy of FCC $0.02 k_B T / N$ between systems of 256 and 1372 particles, the larger system having the larger free energy.

Taking a maximum spring stiffness of $\gamma_\text{max} = 600 k_BT/\sigma_{LL}^2$ for both systems gives
$$
\frac{\beta F}{N} = 2.280960 \qquad \text{for 200 particles}\\
\frac{\beta F}{N} = 2.292628 \qquad \text{for 800 particles}\\
$$

The free energy should be independant of the spring stiffness at which we start the integration.

At the considered density, the S1 solid melts when the springs are turned off. This is probably not the best place to test the algorithm.

Now that I have all the parts ready (but possibly wrong), I will wrap them in a single script to be run efficiently.
I will:
* Try again at higher density
* Use larger systems
* Use monodisperse hard disks for validation

### Scripting
[07/05/2021]

I wrote a script that takes as input an initial configuration, and computes the Helmoltz free energy of the corresponding structure by Frenkel-Ladd method.

The final value depends on the number of integrand samples for $\Delta F_2^*$. For a S1 solid with 200 particles, we get:

| Integrand samples | $\frac{\beta F}{N}$ | Computing time (s) | Integrand plot |
| -------- | -------- | -------- | --------- |
| 30 | 3.258640154019 | 21 | |
| 100 | 3.328710634037 | 42 ||
| 500 | 3.404298636859 | 161 | |

This variation is due to a poor sampling of the low spring stiffness part of the integrand, where it actually varies a lot.

>[name=frank] If I'm not mistaken, the integrand itself does not diverge at gamma = 0, right? Since the energy is gamma * (sum of square displacements), dividing out the gamma does not lead to undefined behavior. It is, however, true that the behavior at the start is very strong, and not very flat. Could you just put the first datapoint at 0? (and measure the sum of squared displacements without multiplying and dividing by gamma?) That might help flatten things out. Otherwise (or perhaps additionally), either of the options you list below sound like good options. For fitting, it might be good to see if you can determine what the rough scaling is near zero of the integrand (from the data).

Possible solutions:
* Non uniform sampling. Maybe the logarithmic change of variable could help.
* Fit and integrate the continuous function.

#### Improved sampling
[12/05/2021]

Following suggestion of Vega et al., I tried a logarithmic sampling for the spring stiffness. Steps are performed in $\ln(\gamma + c)$, with $c = e^{3.5}$. This results in a finer sampling at low $\gamma$.
On top of this, I also tried performing the full variable change in the integral, ie. changing the integrand and the integral bounds.

$$
\Delta F_2^* = - \int_{0}^{\gamma_\text{max}}{d\gamma \left\langle \sum_{i=1}^N (r_i - r_{i_0})^2 \right\rangle_\gamma} \\
= - \int_{\ln(c)}^{\ln(\gamma_\text{max} + c)}{d[\ln(\gamma + c)] (\gamma + c) \left\langle \sum_{i=1}^N (r_i - r_{i_0})^2 \right\rangle_\gamma}
$$

This results in equally spaced points on a curve that should be smoother, and hence integrated with better accuracy. The next figure displays the free energy of a S1 solid of 200 particles at a density of $1.7 /\sigma_{LL}^2$. The errorbars depict the standard deviation estimated from 10 realisations at each point.


The next table shows the integrand sampling directly for the same system.

| Linear, 30 pts | Linear, 100 pts | Log, 30 pts | Log + variable change, 30 pts | Log, 100 pts | Log + variable change, 100 pts |
| -------- | -------- | -------- | -------- | -------- | -------- |
|  |  | | |  |  |

The logarithmic sampling yields significantly better results than the linear sampling. Performing the full change of variable before integration brings no measurable gain in standard deviation.

The next figure displays the impact of the constant $c$ in the logarithmic sampling $\gamma \rightarrow \ln(\gamma + c)$.

The following table depicts the points distribution associated to each of the samplings.

| | $c = e^{3.5}$ | $c = e^{3}$ | $c = e^{2}$ | $c = e^{1}$ |
| -------- | -------- | -------- | -------- | -------- |
| Log sampling | | | ||
| Log sampling + var. change| | || |

It seems that a smaller constant performs better, however, we still need to confirm that the average free energy is converging to the right one when $c$ is decreased.
The best value of the constant might actually change with density and the structure under consideration.

For the moment, I will use logarithmic sampling for the integrand with $c = e^1$ and full variable change.

Next step:
* Try estimate the error on $\Delta F_2^*$ without running the whole thing again. The uncertainty on $\Delta F_1^*$ is easy to get, and $F_{id}$ is exact, so that's really the only missing part.

#### Error estimate on $\Delta F_2^*$
[08/06/2021]

$F = F_{id} + \Delta F_1^* + \Delta F_2^*$. The uncertainty on $\Delta F_1^*$ is directly obtained from the simulation and $F_{id}$ is exact. Therefore, the uncertainty on $\Delta F_2^*$ is the only missing part to get the uncertainty of the whole free energy.

$\Delta F_2^*$ is obtained from integration of a set of points. Running simulations to get those points is the most time-consuming part of the Frenkel-Ladd procedure. Thus, estimating the uncertainty on the computed free energy by repeating the whole procedure many times at the same point and taking the standard deviation is not very efficient.

I tried two methods to estimate the error on $\Delta F_2^*$:
* Generate fake data points and integrate those many times
* Analytical propagation of errors

##### Fake data points

Each point comes with its uncertainty from the simulation. To get more data points without actually running more simulations, I draw a new set of point from Gaussian distributions centred at the simulated points, with a width equal to the uncertainty at that point. This provides fake new realisations of the simulated data points, from which $\Delta F_2^*$ can be computed.
Finally, the standard deviation of the integrals on the fake points is taken as the uncertainty on $\Delta F_2^*$.

To assess the reliability of the method, I computed 100 times the free energy of a S1 solid of 200 particles at density $1.7\sigma_{LL}^{-2}$, using 30 points logarithmically sampled for the calculation of $\Delta F_2^*$, and compared the standard deviation of this *real* realisations with the error estimated as described above.
The average measured free energy per particle is $3.3999 k_BT$ with a measured standard deviation of $0.0067k_BT$ on the one hand. The estimated error on the other hand is $0.0069 k_BT$. The agreement is perfect after rounding to significant digits.

The simulated points can be seen as obtained after one step of a random walk around the ideal curve, with a gaussian step of width given by the measured error on the points. Creating the fake data sets from the simulated points actually samples the second step of this random walk. For most reasonable functions, I expect that this will lead to a slight overestimation of the error, which is safe for errorbars.

##### Analytical propagation of errors

The Simpson rule is used to integrate the data points.
It has an analytical expression that can be used to propagate the uncertainty from the points to the integral directly with the standard formula
$$
\delta f(x_1, \ldots, x_n) = \sqrt{\sum_{i=1}^n{\left(\frac{\partial f}{\partial x_i} \delta x_i\right)^2}}.
$$

The general Simpson formula for $N$ unequally spaced points found on Wikipedia can be reorganised to factorise the contribution of each data point $(x_i, f_i)$. The intervals are denoted $h_i = x_{i+1} - x_i$. Boundary points, odd and even bulk points carry different contributions.
$$
I(\{x_i\}, \{f_i\}) = f_0 \frac{h_0 + h_1}{6} \left(2 - \frac{h_1}{h_0}\right) + f_N \frac{h_{N-2} - h_{N-1}}{6} \left(2 - \frac{h_{N-2}}{h_{N-1}}\right) +\\
\sum_{k = 0}^{N / 2 - 1} f_{2k + 1}\frac{(h_{2k} + h_{2k + 1})^3}{6 h_{2k} h_{2k - 1}} + \\
\sum_{k = 1}^{N / 2 - 1} f_{2k} \left[ \frac{h_{2k} + h_{2k + 1}}{6} \left(2 - \frac{h_{2k + 1}}{h_{2k}}\right) + \frac{h_{2k - 2} + h_{2k - 1}}{6} \left(2 - \frac{h_{2k - 2}}{h_{2k - 1}}\right) \right]
$$

Simpson rule requires an even number of intervals. An odd number of intervals can be handled in different ways by the `simpson` method of scipy. I use the default `even = "avg"` option. The integral is computed twice. The first one uses Simpson integration on the $N-1$ first intervals and trapezoidal rule on the last intervals, and the second one uses Simpson integration on the $N-1$ last inteval and uses trapezoidal rule on the first interval. These two integrals are then averaged to get the final result. This method is more precise because, on top of the average, points carry an index of different parity in each of the integrals.

The trapezoidal rule on one interval reads
$$
I_t = (x_{i+1} - x_i) \frac{f_{i+1} + f_i}{2}
$$
from which the contribution to the error can be easily obtained.

I implemented the error propagation with the `avg` scheme. It yields an uncertainty of $0.0058 k_BT$ for the same system as above. This is smaller than the measured error of $0.0067 k_BT$. I suspect there are approximations behind the formula of error propagation that start to play when applied to a large number of variables.

I will use the fake dataset method to estimate the error on $\Delta F_2^*$. Although it is more expensive than the analytical method, it yields excellent results, and calculations take a few miliseconds: optimisation is not worth the trouble :-).

#### Validation on the hexagonal solid
[17/06/2021]

It is surprisingly hard to find reference values for the free energy of the hexagonal packing of hard disks.

* W. G. Hoover and F. H. Ree, J. Chem. Phys. **49** 8, 3609-3617 (1968)
The authors try to assess the performance of the single cell occupancy method close to and below melting for hard spheres and hard disks. In particular, they want to measure the *communal entropy* which is the entropy lost by restricting particles to individual cells. To do so, they measure the equation of state for the single cell occupancy system and integrate it to get the free energy which is then compared to that of the dense fluid obtained by a Padé approximant. In table IV, they report the entropy of hard disks as a function of density relative to close packing. I assume this means $\rho / \rho_0$ with $\rho_0 \sigma^2 = 1/\sqrt{3}$ the density of the hexagonal packing. To further the matter, the reported results are extrapolated to the thermodynamic limit by a procedure that is presented in the paragraph below table IV. The highest reported relative density is $0.8$, at which the hexagonal solid melts in my simulations. So far, I couldn't get more of this paper.

I could not find other papers reporting free energies directly.
Another strategy: use a high density equation of state and integrate it to check free energy *difference* between 2 densities. I could also directly measure it in EDMD simulations.

* B. J. Alder, W. G. Hoover and D. A. Young, J. Chem. Phys. **49** 8, 3688-3696 (1968)
Provide a high density equation of state for hard disks.

The next figure displays the free energy of hard disks for various densities computed with the Frenkel-Ladd method. The x axis is the actual density $\rho = N / V$. The points are distributed between *relative* densities $0.8$ to $0.98$.
The orange curve corresponds to the free energy obtained with cell theory approximation.
The green and red points are obtained by computing the free energy difference between first (or last) point in the dataset and the point of interest by thermodynamic integration of the equation of state. This difference is then added to the free energy of the reference point computed with the Frenkel-Ladd method. In this representation, the points are almost indistinguishable.
The equation of state is calculated from a solid of 1080 disks (same size as the solid used for FL calculation) and sampled with EDMD simulations, with a density step of $0.012\sigma^{-2}$.
The purple point corresponds to the free energy value reported in Schilling and Schmid, **131**, 231102 (2009) for a solid of 100 hard disks at pressure $10 \;k_BT\sigma^{-2}$. The corresponding density is obtained from the simulated equation of state. The resulting value is compatible with our Frenkel-Ladd calculations within the errorbars (that come mostly from the equation of state). The reported system is much smaller than the one we use, so finite size effects are also expected.


[24/06/2021]

The next figure displays the free energy obtained from thermodynamic integration minus the free energy calculated with Frenkel-Ladd method. The errorbars on the thermodynamic integration points correspond to 1 standard deviation of the results, obtained by generating 50 fake equations of state from the errors of the simulated one.
These points are all compatible with Frenkel-Ladd calculation within 2 standard deviations.


This result validates the FL calculation for free energy differences down to $0.01\; k_B T /N$.
There is still the possiblity for a systematic error in the absolute value, although the good agreement with cell approximation at large density is encouraging.

Next steps

* Compute equation of state to transform to Gibbs free energies. => EDMD ? MC ?

## Solids equations of state

[26/08/2021]

Frenkel-Ladd method becomes less reliable close to melting because at low spring stiffnesses, the systems might become diffusive, leading to a divergence at small $\gamma$ in the calculation of $\Delta F_2^*$.
Therefore, the standard procedure to compute solid free energies is to compute the absolute free energy at a reference density by Frenkel-Ladd method, and then integrate the equation of state to compute the free energy at all solid densities.
This usually works better, because it is possible to superheat (or "overexpand" in our case) the solids in simulations, owing to the periodic boundary conditions that force homogeneous melting rather than melting from a boundary. This allows to compute the free energy down to the melting point and in part of the metastable region at even lower densities.

### Hexagonal solid

The simulated system is initialised as a hexagonal crystal of $1840$ disks in a rectangular box. A rectangular unit cell with aspect ratio $\sqrt3$ is repeated $40$ times in the $x$ direction, and $23$ times in the $y$ direction to keep the simulation box as close to a square as possible, as depicted in the following figure.


The equation of state is obtained via EDMD simulations at fixed density. The following figure displays the obtained EOS. Red points are detected as not equilibrated by the linear fit method (see Miscelaneous). They coincide with what appears to be a Wan der Walls loop indicating a first order transition.
In fact, as recently settled, this loop is a finite size effect. The hard disks freezing transition is more complicated : first order transition from fluid to hexatic, followed by a continuous transition from hexatic to solid. See for instance Bernard and Krauth, Phys. Rev. Lett. **107**, 155704 (2011).


### S1 solid

The equation of state is obtained from EDMD simulations of a square solid of 1024 small and 1024 large disks in a S1 structure.


Again, points in what appears as a transition region are not equilibrated. Simulation snapshots exhibit coexistences in this regime. The loop is larger (in pressure) than in the monodisperse case. I don't know the nature of the transition, but this gets me curious :).

Next step: Integrate the equation of state to obtain the Gibbs free-energy of the solids as function of p.

### Sigma solid

[26/10/2021]

The sigma phase is the first approximant to the ideal square-triangles QC.

It has a composition $x_S = 1/3$ with 12 particles in the unit cell (8 large, 4 small).

#### Structure

The following figure depicts the unit cell of the ideal structure. Dotted line highlight the underlying square-triangle tiling.


Let $a$ be the edge length of the tile. Then, the edge length of the square unit cell is $L = a (1 + \sqrt{3})$
The coordinates of each particle, in units of the unit cell length are given by
$$
(x_0, y_0) = \left(0 ,\frac{1}{2(1 + \sqrt{3})}\right)\\
(x_1, y_1) = \left(0 ,\frac{1 + 2 \sqrt{3}}{2(1 + \sqrt{3})}\right)\\
(x_2, y_2) = \left(\frac{\sqrt{3}}{2(1 + \sqrt{3})}, 0\right)\\
(x_3, y_3) = \left(\frac{2 + \sqrt{3}}{2(1 + \sqrt{3})}, 0\right)\\
(x_4, y_4) = \left(\frac{1}{2} ,\frac{\sqrt{3}}{2(1 + \sqrt{3})}\right)\\
(x_5, y_5) = \left(\frac{1}{2} ,\frac{2 + \sqrt{3}}{2(1 + \sqrt{3})}\right)\\
(x_6, y_6) = \left(\frac{1}{2(1 + \sqrt{3})}, \frac{1}{2}\right)\\
(x_7, y_7) = \left(\frac{1 + 2 \sqrt{3}}{2(1 + \sqrt{3})}, \frac{1}{2}\right)\\
(x_8, y_8) = \left(\frac{1}{4} ,\frac{1}{4} \right)\\
(x_9, y_9) = \left(\frac{1}{4} ,\frac{3}{4} \right)\\
(x_{10}, y_{10}) = \left(\frac{3}{4} ,\frac{3}{4} \right)\\
(x_{11}, y_{11}) = \left(\frac{3}{4} ,\frac{1}{4} \right)\\
$$

We repeat the unit cell 13 times in each direction to obtain an initial configuration with 2028 particles.

|||
|:-:|:-:|
|Initial configuration with 2028 particles.|Bond network between large particles.|

The maximum density is $\rho_\text{max} = 12 / (1 + \sqrt{3}) \approx 1.607695 \;\sigma_{LL}^{-2}$.

#### Equation of state

The equation of state is measured from EDMD simulations starting from the ideal solid configuration, and running for $100\;\sqrt{\beta m \sigma_{LL}^2}$.


### QC approximants

[15/11/2021]

The quasicrystal can only be called such in an infinite system. Our simulation boxes are periodic, so all we can investigate are larger and larger approximants, to approach the QC properties.

#### Stampfli-3

We use Stampfli inflation procedure to generate the initial configuration.
At each step, the tiling is inflated by a factor $2 + \sqrt{3}$ and each vertex is replaced by a "dodecagonal wheel" cluster. Two orientations of the dodecagonal wheel are possible. Stampfli's procedure states that the orientation should be taken at random for each vertex. Here, for simplicity, I only used one orientation.

>[Etienne] what is the implication of not randomising the DDwheel orientation ?

Three inflation steps with a simple square as seed result in the following tiling. It contains $4262$ particles with a small disk number fraction of $1351 / 4262 \approx 0.31698732989$ very close to the QC composition of $x_{QC} = \sqrt{3} / (2 + 2\sqrt{3}) \approx 0.3169872981$.

||  |
|:-:|:-:|
|Stampfli approximant after 3 inflation steps from a square seed.|Bond network between large disks in the Stampfli approximant.|

Starting from this initial configuration, the equation of state of the structure is measured pressures in EDMD simulations up to $\tau = 100 \; \sqrt{\beta m \sigma_{LL}^2}$ for various densities.


The transition is notably smoother than in EOS of the other solids.

## Solids Gibbs free-energy

[05/10/2021]

The equation of state of the solid phases is measured in EDMD simulations which provide a precise measurement of $p(\rho)$. The simulation points are fitted or interpolated to obtain a continuous function which can be numerically inverted to $\rho(p)$.

With the equation of state and the knowledge of the Helmoltz free-energy at a reference density $F(\rho_0)$, computed with the Frenkel-Ladd method, we compute the Gibbs free-energy of the solids phases as

$$
\frac{\beta G(p)}{N} = \beta \int_{\rho_0}^{\rho(p)} \frac{p(\rho')}{\rho'^2} d\rho' + \frac{\beta F(\rho_0)}{N} + \frac{\beta p}{\rho(p)}.
$$

### Hexagonal solid

The equation of state is measured from EDMD simulations of a hexagonal solid of 1840 particles. The EOS measurement is very accurate with an error on the mean pressure of the order of $0.1 \;k_BT \sigma^{-2}$ in the vicinity of the melting transition. Hence, interpolation with cubic splines seems to be more appropriate than fit with a high degree polynomial to smooth the integrand computed from the equation of state, as illustrated with the two following figures.

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

The Gibbs free-energy computed with the above formula is depicted in the following figure. The reference point is $\rho_0 = 1.108512$ with $\beta F(\rho_0)/N = 6.4507 \pm 0.0031$.


The Gibbs free-energy of the hexagonal solid of small disks is readily obtained by observing that the reduced free energy should be the same for both solids at the same reduced pressure.
Thus, we only need to rescale the $x$ axis by $(\sigma_{LL}/\sigma_{SS})^2 = 1 / q^2$.

> [25/10/2021]
> Actually, this is not correct. The thermal wavelength needs to be rescaled as well. From the Frenkel-Ladd expression of the free-energy, it appears that in addition to the $x$-axis rescaling, a term $\ln(q^2)$ needs to be substracted.

[12/10/2021]

Comparison with the Gibbs free-energy computed with the cell approximation shows good agreement and provides an estimate of the magnitude of the error made by the cell approximation.

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

### S1 solid

[06/10/2021]

The equation of state is measured with a good accuracy so we use cubic splines interpolation rather than fit to smooth the curve.

The next figure displays the interpolated integrand computed from the equation of state.


At a reference density of $1.9\;\sigma_{LL}^{-2}$, the Helmoltz free-energy obtained from the Frenkel-Ladd method is $5.150 \pm 0.001 \; k_BT/N$. A system of 2048 particles is used for the Frenkel-Ladd calculation.
Using this reference value, we obtain the following Gibbs free energy for the S1 solid.


[12/10/2021]

Again, the free-energy computed with cell approximation is in good agreement with the one obtain with Frenkel-Ladd method and thermodynamic integration.

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

### Sigma solid

[28/10/2021]

At a reference density of $1.50 \sigma_{LL}^{-2}$, the Frenkel-Ladd calculation gives a Helmoltz free-energy of $4.9619 \pm 0.0009 \; k_BT/N$.

The Gibbs free-energy as a function of pressure is computed by thermodynamic integration of the equation of state.

|||
|:-:|:-:|
|Integrand computed from the simulated equation of state.| Gibbs free-energy.|

### Stampfli-3

[15/11/2021]

At a reference density of $1.53 \sigma_{LL}^{-2}$, the Frenkel-Ladd calculation gives a Helmoltz free-energy of $6.279 \pm 0.002 \; k_BT/N$.

The Gibbs free-energy as a function of pressure is computed by thermodynamic integration of the equation of state.

||
|:-:|:-:|
|Integrand computed from the simulated equation of state.| Gibbs free-energy.|

The Gibbs free-energy thus computed does not incorporate the tiling configurational entropy term relevant for random tiling quasicrystals.
From P. A. Kalugin, J. Phys. A: Math. Gen. **11**, 3599-3614 (1994), the entropy per vertex of the random tiling of squares and triangles is

$$
s_v = \ln\left(\frac{108}{\cos^2\gamma}\right) + \sqrt{3} \cos\left(\frac{\gamma}{3}\right) \ln\left(\frac{2 \cos(\gamma/3) - \sqrt{3}}{2 \cos(\gamma/3) + \sqrt{3}}\right) \\
-\sin\left(\frac{\gamma}{3}\right) \ln\left(\frac{2 - \cos(2\gamma/3) + 3\sin(\gamma/3)}{2 - \cos(2\gamma/3) - 3\sin(\gamma/3)}\right)
$$

with $\gamma$ given by

$$
\sin\left(\frac{\gamma}{3}\right) = \frac{\sqrt{2 \alpha_t - 1}}{\sqrt{3} + (2 - \sqrt{3})\alpha_t}
$$

and $\alpha_t$ is the fraction of area covered by triangles in the tiling.

At the QC composition, squares and triangles cover the same area, so $\alpha_t = 0.5$. Hence, the entropy per tiling vertex reduces to

$$
s_v = \ln(108) + \sqrt{3} \ln\left(\frac{2 - \sqrt{3}}{2 + \sqrt{3}}\right)
$$

The entropy per particle in our binary disk systems is obtained from the entropy per tiling vertex as

$$
\frac{\beta S}{N} = s_v (1 - x_{QC}) \approx 0.08199773528456454
$$

The Gibbs free-energy of the random tiling QC is then obtain by substracting the term $TS$ from the free-energy computed by Frenkel-Ladd method and thermodynamic integration.

Next step:
* Build finite pressure phase diagram

## Finite pressure phase diagram

### First dirty results
[12/10/2021]

The following figure displays the location of the phase boundaries obtained by common tangent construction on the free-energies computed as described above. Only fluid, S1 and hexagonal phases are taken into account.


There are several reasons for not being happy:

* At some reduced pressures (eg: 30), transition points are missed. The following figure shows that this is due to the fluid range being too short for a suitable common tangent to be found. At this pressure, S1 lies still above the fluid curve, while no tangent can be found between fluid and HexL.


The following figure shows the extent of the pressure-composition plane covered by the available functions for the Gibbs free-energy of the candidate structures. The ranges are limited by simulated data availability, as well as phase transitions observed in simulations.


The overlap between fluid and solid regions seems quite good.

* The fluid-solid transition lines are wiggly. This is probably caused by the wiggly derivatives of the cubic splines.

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

* The reduced solidification pressure for HexL is about 11.5, while it should be about 9.2. The reduced solidification pressure for HexS however is in line with the expected value.

This suggests that the cubic spline interpolation might not be a good idea after all. I will try replacing the interpolation by a fit. I could also try improving the fluid free-energy by sampling more finely the composition axis.

### Smooth derivatives

[25/10/2021]

The wiggly derivative can be fixed by using a spline interpolation that allows the value of the interpolating spline to deviate slightly from the interpolated points, or replace interpolation by fit altogether.

I am now fitting a polynomial of degree $\text{deg} = \min(5, \text{Number of simulated points} - 1)$ to get the fluid Gibbs free-energy at constant pressure. The derivatives are now smooth.

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

In addition, I introduced cutoffs to avoid fitting the fluid equation of state in the high pressure, poorly-equilibrated regime.
After calculating a first phase diagram using rough cutoffs, I set them to values larger, but close to the predicted solidification line. This smooths the fluid-solid transition line further for $x_S < 0.5$.

With these improvements, the phase diagram up to $x_S = 0.5$ looks much better.
On the next figure, the black dots correspond to phase changes in the common tangent constructions.
The colored lines highlight the ranges over which the Gibbs free-energies of the various phases can be computed.


In addition to the smoother the solidification line, the predicted freezing pressure for Hex$_L$ is about $9.1 \;k_BT/N$ as expected.

### Corrected Hex$_S$ free-energy scaling

I made a mistake when rescaling Hex$_L$'s free-energy to obtain Hex$_S$'s one.

While it is correct that the two hexagonal solids have the same reduced Gibbs free-energy at the same reduced pressure, I forgot to take into account that the thermal wavelength had to be rescaled as well (see first term in the Frenkel-Ladd expression for the free-energy, thanks Frank for spotting this !). Therefore, in addition to the rescaling of the pressure axis, one needs to substract $\ln(q^2)$ from the free-energy of Hex$_L$ to obtain that of Hex$_S$.

With this correction, the freezing pressure at $x_S = 1$ is correctly predicted around $9.1 \times 0.46^2 \approx 43\; k_BT/N$.

### Fluid issues for $x_S > 0.5$

For $x_S > 0.5$, the phase diagram looks wrong.


The fluid appears ridiculously stable and the solidification lines look like they are never going to meet.

The next figure displays the infinite pressure phase diagram at $\Delta = 0.0709$. The red vertical line highlights the size ratio $q = 0.46$.


From this, we expect T1 and H3 phases to also be stable in our finite pressure phase diagram. This will certainly restrict the fluid stability region.

At $q = 0.46$ and $\Delta = 0.0709$, in the limit of infinite pressure, the two phases have the following unit cells.

|||
|:-:|:-:|
|T1|H3|

T1 unfortunately does not fit in an orthogonal box. Adding support for non orthogonal boxes in simulation codes used by the Frenkel-Ladd method would take time, for a result of moderate interest. Instead, we choose to focus on the first half of the phase diagram, and explore the properties of the random tiling QC more thoroughly.

### Adding sigma phase

[28/10/2021]

The sigma phase is stable nowhere, the phase diagram is unchanged.


The following figure displays the free-energy difference between the sigma solid and the stable coexistence of Hex$_L$ and S1.


### Adding QC

[15/11/2021]

Without its configurational entropy, the RTQC is not stable. The following figure displays the phase diagram with the quasicrystal approximated by 3 times Stampfli inflated square.


The following figure displays the Gibbs free-energy difference between the Stampfli approximant and the Hex$_L$ + S1 coexistence at the same composition, with and without the tiling entropy contribution.


Next steps:
* Compute the free energy difference of various RTQC
* Compute the free energy of whatever NAHD system we like :)

## Perfect vs random QC

We want to measure the free-energy difference, if any, between the perfect dodecagonal QC and its randomised counterparts.
To this end, we need to generate perfect QC approximants of larger and larger sizes. Random QC are obtained by running zipper moves on the perfect approximants.

### Approximants of the perfect QC
[19/04/2022]

Approximants are generated with Schlottmann inflation rule, as revisited by Sadoc (see Sadoc Struct. Chem. 2016 and Impéror Soft Matter 2021).

|$n = 1$|$n = 2$|$n = 3$|$n = 4$|
|:-:|:-:|:-:|:-:|
|||||
|$N_L = 15$, $N_S = 7$|$N_L = 209$, $N_S = 97$|$N_L = 2911$, $N_S = 1351$|$N_L = 40545$, $N_S = 18817$|

### Free energy of the perfect approximants
[25/04/2022]

To compute the Helmoltz free-energy of the perfect quasicrystal, we compute the free-energy of larger and larger approximants with Frenkel-Ladd method, and then extrapolate to infinite system sizes.

To prevent the system from becoming diffusive when springs are turned off in the Hamiltonian integration, we need to work at "high enough" density. Looking at the finite pressure phase diagram above, we expect the system to be a dense solid for pressures above $30 k_B T \sigma_{LL}^{-2}$. Judging from the equation of state in "Solids equations of state > QC approximants", this corresponds to densities above $\sim 1.45 \sigma_{LL}^{-2}$. I choose $\rho = 1.5 \sigma_{LL}^{-2}$ to be on the safe side (and because it is a nice round number).

At this density, the measured Helmoltz free-energies are:

| Inflation step $n$ | 1 | 2 | 3 | 4 |
|:-:|:-:|:-:|:-:|:-:|
| $\frac{\beta F}{N}(\rho = 1.5)$ | $5.71 \pm 0.06$ | $5.487 \pm 0.007$ | $5.502 \pm 0.002$ | $5.5033 \pm 0.0004$ |

The first approximant is too small to be useful in the estimation of the finite size scaling. The free-energies of the larger approximants however seem to follow a linear finite size scaling as expected.


#### Equilibration checks

The calculation of $\Delta F_2^*$ in Frenkel-Ladd procedure involves measurements of the average square displacement $d2$ of particles for various spring stiffness. Measurements are taken every 10 Monte Carlo cycles during $N_{c}$ cycles, after $N_c$ cycles of equilibration. Here, we compare results obtained for $N_c = 10^3$ and $N_c = 10^4$.

The following figures display the measured average square displacement as a function of spring stiffness, for each of the approximants. Points in red are flagged as non-equilibrated by my automatic equilibration check (see Miscelaneous note on the HackMD).

|Inflation step| $n = 2$ | $n = 3$ | $n = 4$ |
|:-:|:-:|:-:|:-:|
|$N_c = 10^3$||||
|$N_c = 10^4$||||

The following figures show the square displacement as a function of MC time at $\gamma = 0$ for the 3 approximants.

|Inflation step| $n = 2$ | $n = 3$ | $n = 4$ |
|:-:|:-:|:-:|:-:|
|$N_c = 10^3$||||
|$N_c = 10^4$||||

At low spring stiffness, the systems are not properly equilibrated. Hence, $d2$ is underestimated for these points. This results in an overestimated free-energy since $\Delta F_2^*$ is the negative integral of $d2$.
The overall effect seems to be of the order of a few $10^{-3} k_BT$ on the final free energy per particle.

|$N_c$| $n = 1$ | $n = 2$ | $n = 3$ |
|:-:|:-:|:-:|:-:|
|$10^3$| $5.478 \pm 0.006 k_BT$ | $5.502 \pm 0.002 k_BT$ | $5.5045 \pm 0.0006 k_BT$ |
|$10^4$| $5.480 \pm 0.002 k_BT$ | $5.5019 \pm 0.0006 k_BT$ | $5.5034 \pm 0.0002 k_BT$ |


How far we are from equilibrium is unclear. No particles are seen to escape their local cage in the snapshots.

[27/04/2022]

To get a better impression of the equilibration time for d2, I ran long simulations for the 3 approximants. Following figures display the measured d2/N as a function of MC cycles for small spring stiffness.

|$\beta\gamma\sigma_{LL}^2 = 0$|$\beta\gamma\sigma_{LL}^2 = 1$|$\beta\gamma\sigma_{LL}^2 = 3$|
|:-:|:-:|:-:|
||||

It seems that the larger systems are more "wiggly" than the small ones, which is not surprising because the influence of the single fixed particle is smaller.
As the systems grow fluctuations tend to decrease, but remain rather important. Large fluctuations seem to correspond to large scale collective motions rather than individual rare events. In particular, no particles are observed to escape their local cage.
As expected, larger systems take longer to equilibrate. It seems that the largest approximant reaches some kind of equilibrium after a few $10^5$ MC cycles, while the other equilibrate faster.
Very small (but non-zero) spring stiffness help a lot to kill the large fluctuations and dramatically decrease equilibration time.

[01/05/2022]

I modified the script to allow for automatic increase of the simulation time, until the simulations are classified as equilibrated by the automatic checks.

The following figures summarises the free energies calculated with various parameters with the Frenkel-Ladd script.


The following figure displays the automatically-adjusted duration of the equilibrated run as a function of spring stiffness for the approximant $n = 4$, with a base $N_c = 10^3$.


The equilibration times are in line with the previous tests with longer simulations. Equilibration time seems to grow exponentially at low spring stiffness.
Walltime duration of this free-energy calculation was 17 hours. Should the computation time become an issue, we could use unrestricted EDMD simulation at $\beta \gamma \sigma_{LL}^2 = 0$ and measure d2 by substracting the position of the carrier from the rest of the lattice. This should equilibrate much faster.

### Perfect QC free-energy
[06/05/2022]

The following figure shows the free-energy of the perfect approximants measured on 30 independant runs for each size. The black points and errorbars correspond to average and standard deviation across the runs.

For these runs, we use the automatic equilibration scheme described above.


It is interesting to note that the errorbar estimate obtained by propagating errors in the Frenkel-Ladd calculation are essentially the same as the measured standard deviation on the 30 samples.

[09/05/2022]

In order to understand better the origin of the statistical variability, I computed for each spring stiffness, $\Delta d2$ the difference between the largest and smallest $d2$ measured in at least 25 samples.

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

Dashed lines are linear fit to the data points. The slope is about $-0.82$ in the three cases.
At large spring stiffness, $\Delta d2$ "bounces". I have no clue of what could cause this.
The points at $\gamma = 0$ for $n = 3$ and $4$ lie outside of the linear trend, indicating that the spread is significantly larger there. Small spring stiffness dominate the variability in the calculation of $\Delta F_2^*$, so there is no need to use a longer base length for the simulation. $N_c = 10^3$ is enough to measure $d2$ with high accuracy at large $\gamma$.

To reduce errorbars, I ran many independent Frenkel-Ladd calculations. The following figures display the distribution of the resulting values and the mean values with their errorbars.

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

Within the errorbars, points no longer follow a straight scaling as a function of $1 / N$.
It seems that all 3 distributions are slightly skewed to the left (their skewness is indeed slightly negative).

### Random QC

Random QC are generated from the perfect configurations by applying 1000 zipper moves. Free-energies are then computed with Frenkel-Ladd procedure, and same parameters as for the perfect case.

The following figure displays the free-energies of both perfect and random QC at different sizes. Frenkel-Ladd calculation is performed 30 times on the perfect configuration at each size. For the biggest approximants, floating point errors accumulate in the zipper move randomisation, leading to crashes or particle loss in a fraction of cases. For the following figure, I managed to get 10 samples of the first approximant, 9 of the second, and 8 of the largest one.


Within the current errorbars, free-energies of perfect and random QC are indistinguishable. This bolsters the idea that the quasicrystals we observed are indeed well described by the random tiling picture: all configurations of the random tiling are almost degenerate, and act as a single macroscopic phase with finite configurational entropy.
However, the perfect QC seems to be consistently on the low end of the spectrum. Maybe some additionnal numerical efforts could provide an estimate of the free-energy difference between perfect QC and the bulk of the random-tiling ensemble. Such a tiny difference would probably not give any realistic advantage to the perfect QC over the other configurations, but could open the question of the extra vibrational modes, maybe linked to self-similarity symmetry, that provide the tiny extra entropy.

The errorbars could be reduced by simply running more independant runs.

For the random QC, so far, every configuration is used once. I could also run multiple calculations on each configuration to lower the statistical error on the average free-energy of that configuration.
At this point, doesn't make more sens to try tweaking the Frenkel-Ladd parameters to run longer, more precise calculation ? How far can we meaningfully go in this precision quest ? What are the possible sources of systematic errors that could blurr everything anyway ?

[14/05/2022]

As done for the perfect QC, I let the cluster accumulate lots of statistics to improve errorbars on 5 different random tilings at each size.


The free energy of the random tilings appears about $10^{-3} k_BT$ larger than that of the perfect QC. This entropy difference is of the order of the difference between HCP and FCC in 3D HS : it is small, but measurable !

The free-energies of the random tilings are indistinguishable within the current errorbars of the order of $10^{-4} k_BT$.

This suggests that the perfect tiling is indeed special in the family. Could special vibration modes, maybe linked to the self-similarity symmetry, explain the extra vibrational entropy ?

The degenracy of the random tilings bolsters the idea that the random tiling is indeed a proper phase, with non-zero configurational entropy. This configurational entropy stabilises it over the perfect configuration in probably most soft matter QCs.

Paradox ? Assuming that all random configurations have their own, different *real* free-energies. If we were able to tell them apart based on their free-energy, we would consider them as separate phases. Then, the RTQC does not exist anymore, and the stable phase becomes the perfect QC. Can we expect all random configurations to have the same free-energy in the thermodynamic limit ? This feels wrong since we can make arbitrarily small phason perturbations to the perfect tiling. I guess this points to the role of defects in QC (including phasonic defects to the perfect configurations). I should also get better understanding of what are "equivalent" QC structures in the thermodynamic limit.

### Sigma seed
[16/05/2022]

In order to have QC of intermediate sizes, Frank suggested trying to inflate a sigma seed.
This turns out more difficult than expected.

### QC with Hex symmetry
[14/06/2022]

The larger entropy of the DDQC could be due to:
* the different local environments that can exist in the random configurations but not in the inflated tiling (eg : 4444)
* special excitations associated to the self-similarity property.

In order to assess these contributions, we compare the free-energy of the perfect dodecagonal QC and perfect hexagonal QC. Both have the same local environments **(I think, but need a proper check)**, but different symmetries.


The hexagonal QC appears to have slightly less entropy than the dodecagonal one. Self-similarity, or 12-fold symmetry, seems to be favourable for vibrational entropy !

### Finite size corrections

[22/09/2022]

I have tried applying the finite size correction of Frankel and Smit (adding a $\ln(N)/2N$ to the free energy). I don't think that the result is better and propose using the uncorrected version.

|||
|:-:|:-:|
|No correction|Correction|

### Litterature

[10/06/2022]
* Reinhardt, A., Romano, F. & Doye, J. P. K. Computing Phase Diagrams for a Quasicrystal-Forming Patchy-Particle System. Phys. Rev. Lett. **110**, 255503 (2013).
Phase diagram calculation for a system of 2D pentavalent patchy particles.Arguably the first precise numerical calculation of a QC.
The system forms a QC12 in a region of the phase diagram (competes with Sigma, Z, approximant and fluid).
Authors perform free energy calculations for the QC with direct coexistence with the fluid. By substracting the enthalpic part, they can obtain the entropy term. By subsequently substracting the entropy of the approximant, they are able to estimate the tiling entropy as ~0.113 kT/N, which is actually close to the exact value of the maximally random tiling !

[11/06/2022]
* Kiselev, A., Engel, M. & Trebin, H.-R. Confirmation of the Random Tiling Hypothesis for a Decagonal Quasicrystal. Phys. Rev. Lett. **109**, 225502 (2012).
A remarkable paper ! The reference list is a golden mine.
Authors investigate the stability of a RTQC10 in Lennard-Jones-Gauss system. They compute separately the energetic, phononic (vibrational entropy) and phasonic contributions to the free energy, as a function of temperature, and extract the phason elastic constants as a function of temperature.
They show that at low temperatures, one of the phason elastic constant is negative, hence the RTQC cannot be stable and an approximant is formed instead. At higher temperatures, both phason elastic constants become positive, and a crystal-quasicrystal transition is observed.
In this regime, the free-energy is quadratic in the phason strain, as predicted by the random tiling hypothesis.
Frenkel-Ladd calculation is used to get the phonon free-energy. If I understand correctly, the phason contribution is obtained from accelerated simulations where the flips are generated following an Ising-like Hamiltonian ansatz that assumes uncorrelated flips (confirmed a posteriori).
Authors argue that the phonon and energy contributions prevent the RTQC to be the ground state at low temperatures. Our model with hard particles doesn’t have this feature of course, so could we argue that it corresponds to a purely entropic QC ground state ?

* Haji-Akbari, A. et al. Disordered, quasicrystalline and crystalline phases of densely packed tetrahedra. Nature **462**, 773–777 (2009).
Authors report formation of a DDQC in 3D MC simulations of hard tetrahedra system.
The QC is a stack of square triangle random tilings. Since interactions are hard, entropy is the only driving force for quasicrystallisation.

[03/08/2022]

* Henley, C. L. Clusters, phason elasticity and entropic stabilization: a theoretical perspective. Philos Mag **86**, 1123-1129 (2006)
A personal perspective by Henley on some fundamental QC questions. He is mostly concerned with atomic QC, but the questions he raises are equally relevant for soft QC. The author shows that in atomic QC, interaction via Friedel oscillations is essentially the same thing as Hume-Rothery stabilisation. He discusses the cluster models as purely descriptive models. The clusters are not physically consistent entities in the vast majority of cases. Moreover, he states very clearly the issue of "ideal vs random QC":
> The two fundamental competing scenarios of the stabilization of quasicrystals are not exactly ‘entropy’ versus ‘energy’, as we often loosely say. Rather, the question is whether the model is in the qualitative class that contains quasiperiodic ideal tilings that have purely Bragg peaks, or the class that contains the maximally random tiling in which long range order is an emergent phenomenon. This distinction is a rigorous one from the viewpoint of statistical mechanics, because these two states are separated by a phase transition, but that is of no help in distinguishing them experimentally.

He also makes the distinction between tiling and vibrational entropy, suggesting that the later might have been overlooked in the stability picture of QC:
> The evidence is abundant that many quasicrytals are high-temperature phases, ergo stabilized by entropy – but, in many cases, not the tiling configurational entropy, for that is too small. So, the larger entropy of vibrations or of chemical disorder is the only candidate to affect the phase diagram. In the past, we brushed aside such entropy contributions, by claiming they must have a similar value in the quasicrystal phase as they do in the approximant phases. The latter phases were assumed to be the immediate competitors of the quasicrystal phase in the phase diagram, so the difference in vibrational or substitutional free energies would largely cancel. But perhaps we need to examine more carefully just how the quasicrystal ordering may affect these entropic terms.

Finally, Henley mentions the stoichiometry freedom that arises in random tilings with more than 3 types of tiles. To him, this is a relatively unphysical feature:
>In the simplest cases, e.g. the square-triangle tiling with dodecagonal symmetry or the rhombohedral tiling with ico symmetry, there are two kinds of tiles, and the (irrational) ratio of their numbers is fixed by the symmetry. The decoration is deterministic, so there is a unique stoichiometry. However, tilings related to real quasicrystals often have more than two tiles, e.g. the ‘canonical cell tiling’ with ico symmetry or the ‘hexagon-boat-star’ tiling with decagonal symmetry. Then within the sum rules fixed by symmetry, one kind of tile can be traded for a combination of the others; sometimes the decoration is such that this trade leaves the atom count unchanged, so our model is line compound in this case too. In the case that the atom content changed, though, the alloy composition would indeed be variable; but as Elser once suggested, such a decoration is undesirable since phason dynamics would be coupled to, and slowed by, mass diffusion, contrary to the observed relaxation in good quasicrystals.