# Polydisperse coexistence (Lotte)

###### tags: `Masters projects`

###### Notes on Markdown (also see https://hackmd.io/SUXGdcptTaqJIjD3a13j2A )
<!--
:::success
making a green note
:::

:::info
:::

:::danger
:::

:::warning
:::
`literals`

  print('Hello, world!')
  # See if you can set arbitrary line number or 
  # change syntax highlighting for different language
  # in this code block.

--->

## See results at bottom of the page

## Overview

1) For each trial fluid = $\{p, \rho_f\}$, find $P_f$ and $\delta \mu(\sigma)$
a) $P_f$ measured
b) $\delta \mu(\sigma)$ $3^{rd}$ order polyniomal with coefficients fitted

:::info
Actually, pressure measured in semi-grand ensemble.
:::

2) Find the crystal = $\{\rho_\chi, \frac{a}{c}, \frac{a}{b}\}$ that matches the fluid in the semi-grand ensemble with enforced $\delta \mu(\sigma)$
a) $P_f = P_\chi$
b) $\mathbf{P_{xx}} = \mathbf{P_{yy}} = \mathbf{P_{zz}}$ and $\mathbf{P_{xy}} = \mathbf{P_{xz}} = \mathbf{P_{yz}} = 0$

:::info
Use symmetry argument of the crystal to first set a/b = 1. When coexistence phase point is found, repeat some measurements for a/b > 1 too show Pxx is not Pyy in that case.
:::

3) Direct coexistence to determine the proper fluid/crystal pair

## To do
setting up
- [ ] install openGl on MacOS, make new executable file
- [ ] familiarize with all the flags options for running
- [ ] try velgrind for debugging
- [ ] godbolt.org for optimization explorer
- [ ] get used to "debugging" function in VS code

organisational
- [ ] make plan how to learn using Mathematica

read
- [ ] reread paper on optimization EDMD methods
- [ ] check out Franks PhD boekje
- [ ] quasicrystal paper marjolein dijkstra

understand in more detail/do derivation
- [ ] Why mass of dif sized particles does not matter
- [ ] influence of interface, thermodynamic limit

## Keep in mind
- free energy of the different layers very close together, might end up with a mix of layers, like random HCP

## Notes on meeting 18/08/2025

#### Lotte's EDMD Specifics

- cell list, no neighbour list
- binary tree, no bucketing
- max. hard sphere size, polydisperse, but fixed distribution in time
- cubic lattice
- two types of events: collisions or cell change
- schedule only first event of each particle
- check Frank's code for inspo

#### Laves Unit Cell Code Specifics
- set number of unit cells in each direction
- make sure multiple rotations are possible (but limited number of these)
- tune a/b/c distances
- no polydisperisity (yet) just simple size ratio around 0.75

<!-- #### EDMC-poly check out

- [ ] kickback (where you need BOP6 etc.)
- [x] linked list details, cellist is a pointer to a particle, p->next points to anohter part in that cell,
- [ ] why p->vx,vy,vz set randomly in randomparticles if you later just make it randommovement
- [x] neighbourlist, xn
- [x] list counters en merge counters (//The first is the average number of events in the first that gets moved into the event tree after each block. //The second is the average length of the overflow list.)
- [ ] final pressure measurement = zero
- [x] because you put event back into the eventpool, you can first removeevent() and only then excecute if you like
- [x] still have cell lists, but no cell crossings. Just when you have to define a new neighbourhood, you first reassign the cell the particle is in. Can't be sure that the first collision will be with one of these neighbours. But if a far away particle comes at you very quickly, it will schedule its own neighbourhood events until it meets you and then things will be sorted again.
- [ ] I let it exit:
if (tmin < 0) printf("[EDMC-poly] WARNING negative time: type: %d tmin: %lf time: %lf p1: %d p2:%d\n", type, tmin , time, p1->number, partner->number); exit(0);
if (partner->number == p1->number && type == 0) printf("[EDMC-poly] WARNING p:%d colliding with itself, tmin: %lf type: %d\n", p1->number, tmin, type); exit(0);
- [ ] const double mR = 1; // the growing/shrinking velocity is also changed upon collision, where does this factor come from? Also mass/"radius mass" ratio does not matter?
- [ ] for fun: check out tailored version of findradial -->

## Theory

#### Virial expression for HS

$$
\beta P_{k l} / \rho=\delta_{k,l}-\frac{\sum_c m v_{i j}^k r_{i j}^l}{N\left(t_b-t_a\right)},
$$

(From Laura and Franks paper:)
where $P_{k l}$ indicates the $k l$-component of the pressure tensor (with $k, l \in\{x, y, z\}), \rho=N / V$ is the number density, and the sum is taken over all collisions in a time interval $\left[t_a ; t_b\right]$. Additionally, $r_{i j}^k$ and $v_{i j}^k$ are the $k$-components of the relative distance and velocity vectors between the colliding particles $i$ and $j$, respectively.

And then, pressure scalar is:

$$
P = \frac{1}{3}\cdot\mathrm{Tr}(\mathbf{P})
$$

#### Surface forces
Derivation of why chemical potential' can be linked to ensemble average of -outwards force on particle surface: see supplementary methods.

(From Laura and Franks paper:)
We can write the ensemble average:
$$
\left\langle f_{R_i}\right\rangle=\frac{-1}{t_b-t_a} \sum_k m_i\left|\left(\delta \mathbf{v}_i\right)_k\right|
$$
where $\left(\delta \mathbf{v}_i\right)_k$ is the change in velocity of particle $i$ during collision $k$, and the sum is taken over all collisions in the time interval $\left[t_a, t_b\right]$.

#### The mixing contribution (on hold)

The (Gaussian) distribution of the polydisperse system is as follows:

$$
\mathcal{P}(\sigma)=\frac{1}{p \bar{\sigma} \sqrt{2 \pi}} \exp \left(-\frac{1}{2}\left(\frac{\sigma-\bar{\sigma}}{p \bar{\sigma}}\right)^2\right)
$$

The mixing contribution of the chemical potential as a function of particle size is $\beta \mu_{\mathrm{mix}}(\sigma) = -\mathrm{log}(\mathcal{P}(\sigma)\lambda)$.

Take $\lambda = \bar{\sigma}$ (as written in Supplementary Methods), then,

\begin{align}
\mathrm{log}(\mathcal{P}(\sigma)\lambda)
&=\mathrm{log}\left(\frac{\lambda}{p \bar{\sigma} \sqrt{2 \pi}} \exp \left(-\frac{1}{2}\left(\frac{\sigma-\bar{\sigma}}{p \bar{\sigma}}\right)^2\right)\right) \\
&= \mathrm{log}\left(\frac{\bar{\sigma}}{p \bar{\sigma} \sqrt{2 \pi}}\right) -\frac{1}{2}\left(\frac{\sigma-\bar{\sigma}}{p \bar{\sigma}}\right)^2\\
&= -\frac{1}{2}\mathrm{log}\left({p^2 2 \pi}\right) -\frac{1}{2}\left(\frac{\sigma/\bar{\sigma}-1}{p }\right)^2
\end{align}

$\beta \mu_{\mathrm{mix}}(\sigma) = \frac{1}{2}\mathrm{log}\left({p^2 2 \pi}\right) + \frac{1}{2}\left(\frac{\sigma/\bar{\sigma}-1}{p }\right)^2$

This is used to find the coefficients of the potential $V(R)$ that produces a Gaussian polydispersity fluid. In combination with the surface force measurements.

#### Relative chemical potential

Relative to what?
That is different for each plot, no one size where it is relative to. I know $\delta\mu_{\mathrm{conf}}(0) = 0$, because we get that from integration and do not set a constant value. And depending on the polydispersity, the $\mu_{\mathrm{mix}}$ function differs.

Look at it this way, if there are no surface forces, mix is all you have. So the mix will be all that is responsible for the gaussian distribution. High mu, means that the resevoir you bring your system into contact with, has a high mu for that size. So there will be flow from high mu to low mu. So you will get a lot of particles of the highest mu.

#### Neighbour list

I am used to cell lists, which are still used in the EDMC-poly code. But additionally, there is a neighbour list.

Particles i and j are neighbours when:

$$
\left|\mathbf{r}_{\text {neigh }, i}-\mathbf{r}_{\text {neigh }, j}\right|<(1+\alpha) \sigma_{i j},
$$

with $\alpha > 1$. This way, you don't have to worry about particles i and j meeting if they are not in each others neighbourhood. Until one of them has moved by $\alpha R_{i,j},$ then recalculate the neighbourhood.

How does this work with polydispersity/continuously fluctuating size? (EDMC-poly code: rm = 2 * maxparticleradius * shellsize;) Just take the max diameter for determining neighbours (pfew).

## Results

### Monodisperse pressure measurements
Simulation = Lotte's EDMD, monodisperse FCC crystal at range of packing fractions: ([0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7])

From pressure **tensor** measurements, isotropy can be checked. Indeed off-diagonals are zero and diagonal elements are equal to one another.

- [ ] **Figure in thesis.** Isotropy check. Now, check if Pxx, Pyy and Pzz are equal within error bar. And check if off-diagonal elements are indeed 0 within error bar.
- [ ] Additional check, deformed FCC crystal should give non-isotropic pressure tensor

The pressure scalar can be compared to theory value for monodisperse fluid phase (see fig). Good agreement

- [ ] **Figure in thesis.** Carnahan Starling comparison, with error bars. Note: simulation starts with a crystal so phase transition is pushed towards lower packing fraction.
- [ ] redo measurements with fixed bug

### Surface forces (polydisperse)

In paper: 0.94, 0.95, 0.96, 0.97, 0.98 in unit of average diameter, $\bar{\sigma}$, at p = 0.06

In my code: units of max diameter, $\sigma$, so the measured force and particle sizes are in terms of $\beta f_{\sigma_i}\sigma$ and $\sigma$ respectively. Unit conversion done in python script, not C.

Also note the more complex relation between packing fraction and number density in terms of $\rho\bar{\sigma}^3$. For each set or loaded packing fraction and polydispersity, the number density in units $\rho\bar{\sigma}^3$ are written in a general file, which is read by python script.

Simulation: p = 0.06, total simulation time = 400. (takes around 15(?) minutes), measures surface force per every $\Delta t = 1$. Averages taken from t = 30. (Equilibrium time estimated from pressure measurements over simulation time.)

Data in folder ./data0209-2

Surface force measurements from my code is nicely described by the fitted line from Frank and Laura's paper.

Performed the same quadratic 2D fit, like described in Frank and Laura's paper equation (6), on my data myself. This results in very similar lines, while the fitting coefficients themselves are rather different. (This means it just found a different set of coefficients describing very similar behaviour. Not a problem.)


- [x] extra densities to improve fit
- [ ] make bins to get a line with error area (like in paper)
- [ ] redo measurements with fixed bug

### Laves phases
#### $\mathbf{{MgCu}_2}$

First, from primitive vectors on [here](https://http://web.archive.org/web/20111228194711/http://cst-www.nrl.navy.mil/lattice/struk/c15.html).
primitive unit cell has 2 Mg and 4 Cu.


Then, mask out a cubic unit cell of size (a, a, a) in the middle of the large piece of crystal, to find coordinates. This gives one cubic unit cell, which has 8 Mg and 16 Cu. (Coordinates listed in CubicUnitCellCu.c)


These unit cells can be stacked to get full crystal.

:::success
**Cubic MgCu**: can be used for rectangular box simulations
:::

In order to get the layered orientation from [here](https://www.nature.com/articles/nmat1841/figures/1), I made a large piece of crystal from cubic unit cells, rotated the coordinates and masked again. This way, I obtained a new rectangular unit cell which I refer to as the "Layered" unit cell, which has 12 Mg and 24 Cu.


(Note: the nature image is not a unit cell, instead slightly more than one unit cell.)

:::success
**Layered MgCu**: can be used for rectangular box simulations
:::

#### $\mathbf{{MgZn}_2}$

Again, started from [primitive vectors](http://web.archive.org/web/20111223180234/http://cst-www.nrl.navy.mil/lattice/struk/c14.html), which have 4 Mg and 8 Zn.

Complication: free parameters z2 and x1 need to be chosen. Taken from Frank's notebook.

Large piece of primitive crystal masked (a, sqrt(3)b, c) to find the rectangular unit cell again, with 8 Mg and 16 Zn.

Note that this is actually already the same as the layered unit cell [nature ](https://www.nature.com/articles/nmat1841/figures/1) can be found in the crystal in this orientation.


:::success
**Rectangular MgZn2 = Layered MgZn2** can be used for rectangular box simulations
:::

#### $\mathbf{{MgNi}_2}$

Same procedure as with MgNi2. Rectangular unit cell = layered unit cell and has 16 Mg and 32 Ni.


:::success
**Rectangular MgNi2 = Layered MgNi2** can be used for rectangular box simulations
:::

### Fluid in semi-grand canonical

Using the EDMC_poly.c code, simulated fluids of a specific density and $V(R)$. Those $V(R)$'s are based on surface forces from the canonical ensemble simulations of Gaussian size distributed system. So the resulting polydispersity in the semi-grand ensemble should match that number.

So, I simulated, EDMC_poly code, maxtime = 1000, 16384 particles, with the coefficients given for density = 0.9582 and p = [0.04, 0.06, 0.08, 0.10]. All these sizes written in mov.osph document every $\Delta t = 1$. (Made sure that the measurements were uncorrelated by checking the spearman coefficient.) Note that indeed mean size was 1 in reduced units. All the sizes binned into a histogram (sqrt(#measurements) number of bins). From these points, Gaussian curve was fitted:


Indeed, the fitted Gaussian gave a peak value of 1 $\bar{\sigma}$. The measured polydispersity indeed corresponds to the set polydispersity.


For each of these trial fluids (density, delta_mu), I will measure the pressure. This is used later.

### Finding the happy box

#### Intermediate: canonical ensemble

I first checked the stability and pressure tensor isotropy of the MgCu2 phase in the canonical regime (my own EDMD code).

+ FCC fluid = isotropic
+ FCC crystal = isotropic
+ Laves is isotropic for number density 1.8*
+ Laves is **not** isotropic for high number density 1.92*



The instability/anisotropy at highest density shown above probably because this is very high pressure. At the highest pressure, Laves phase is **not ** actually stable. Rather, FCCs and FCCl take over. Therefore, the Laves crystal deforms as it wants to change phase.

When changing the seed, the direction of the anisotropy changed too, or this splitting of {Pii} and {Pjj,Pkk} did not occur at all.

Logically, I will henceforth stay at reasonable pressures/densities.

In any case we may conclude that for this bidisperse MgCu2 system, the happy box is cubic. (a = b = c)

But, for each trial fluid with set [density, polydispersity], there is a potential $V(R)$, which will also be the potential for the simulated crystal.

To find the happy box, we have to see when Pxx = Pyy = Pzz and the off diagonal elements are zero in the crystal, at a certain $V(R)$.

So we may find a different a/c ratio for each trial fluid.

:::warning
\* Be sure to use the right number density! In my own code, I set the maximum diameter to one, instead of the average diameter. The reported number densities above are therefore different. CubicUnitCellCu_new.c makes sure the average diameter is 1 and sets the right number density in units of $\bar{\sigma}$. So this code should be used to make init files for the EDMC_poly code.
:::

#### Add polydispersity

:::info
From here number density is reported in the right units, namely the average diameter $\bar{\sigma}$. See warning above.
:::

For the box a = b = c, we see isotropic pressure tensor of a stable Laves crystal at number density 1.2. (At polydispersities 0.08 and 0.14, not 0.04.) As expected, this happy box for MgCu2 is cubic.

The figures below are for $V(R)$ with fluid at number density 1.2 (which is really high). At some point, for lower number densities, $V(R)$ changes such that binary structure is no longer kept.

We can track the sizes of the particles as the simulation runs, so we can see how the mean size goes up a little. This is not problematic, it can very well be that a crystal is in stable coexistence with a fluid of a different composition, inclusing mean size.


We start with a strictly bidisperse initial configuration (so only 2 sizes exact). At p = 0.14, all the simulation does is spread out these peaks a little, but the binary character remains.


At a lower polydispersity, the $V(R)$ does not allow the Laves phase to remain stable. What you see instead is a single peak in the size distribution. Conclusion is that at this number density, a p = 0.06 fluid can not coexist with a Laves phase.


:::warning
Note that the sizes measurements are now correlated. Do different runs and take an average to have any real statistical strength.
:::

#### Weird at low density

A weird phenomenon appears at very low number density (0.6). Something about the low density makes the size distribution curves not gaussian but sort of truncated on the high diameter side. The fitted polydispersity also does not correspond to the set polydispersity either. This may be due to a concavity in the $V(R)$, which could be introducing a fault in the EDMC code. In any case, we will not use these low pressures, so I stopped investigating it for now, but see below for some comparisons.

Normal, fine for number density 1.0 at p = 0.10 and 0.14:

Normal, fine for number density 0.6 at p = 0.10 and 0.14:

:::warning
Probably a good idea to check the polydispersity in the used trial fluid simulations.
:::

### Intermediate "happy box" results a/c

To find the happy box we need Pxx = Pyy = Pzz and the bidisperse character to be preserved. We adjust a/c ratio to see how pressure is influenced.

Below, shown at $\rho$ = 1.0:
a = 1 and c = 1 to 1.08 preserves bidisperse character


As you start to deform the box in the z direction you see how Pzz (green) start to deviate from Pxx and Pyy (red and blue)



a = 1 and c = 1.09 and 1.1 loses bidisperse character


This means, if you push c/a too hard, you lose the bidisperse character: melt the Laves crystal. (I do not recall looking at these movies.) With the melting, the balance between P's is restored.

The off diagonal elements don't really show different behaviour. All still very small values, not like one value is on average much higher.

### Visualise the intersection method

Just to show the process, I have prettier figures now.

For one set of {$p, \rho_f$}, we measure trial crystals {$\rho_\chi, c$}. (Where I just put a = 1 so "value of c" is really c/a.)

**For each $\rho_\chi$.** Find intersect to get true c value and P value.

Find intersection with trial fluid pressure to get the final value of $\rho_\chi$.

:::warning
If c value differs for each density, you might want to also do the interpolation the other way around, first density, then c value, to get the exact right value for c. But probably the c value will be 1 for now anyways.
:::

### Results

Fluid simulation:
- 10000 particles
- 1000 simtime
- polydispersity = 0.10, 0.12, 0.14
- density = 0.92 to 1.08
- write verbose every $\Delta t$ = 1
- write movie every $\Delta t$ = 100
total = 27

:::danger
At too high fluiddensities, polydispersity is not properly retrieved. This is perhaps due to fcc-crystallization of the fluid? So, measuring coexistence at too high fluiddensities (see Frank and Laura paper), is not reasonable. See the section below for some pictures.
:::

Crystal simulation:
- 3000 particles
- 1000 simtime
- polydispersity = 0.10, 0.12, 0.14
- density = 0.84 to 1.1
- fluiddensity = 0.92 to 1.08
- c/a = 0.98, 0.99, 1.0, 1.01, 1.02
- write verbose every $\Delta t$ = 1
- write movie every $\Delta t$ = 10
total = 1890

In the first version of the code "EDMC_poly_fluid.c" and "EDMC_crystal", there was a bug: t<0 was checked and exited instead of t<-eps. Which meant 57 simulations crashed midway trhough the run. These simulations are stored in "crashed" folder and those phase points were rerun.

Additional crystal simulation:
- 3000 particles
- 1000 simtime
- polydispersity = **0.06, 0.08**
- density = 0.84 to 1.1
- fluiddensity = **0.98** to 1.08
- c/a = 0.98, 0.99, 1.0, 1.01, 1.02
- write verbose every $\Delta t$ = 1
- write movie every $\Delta t$ = 10
total = 840

Additional fluid simulation:
- 10000 particles
- 1000 simtime
- polydispersity = **0.06, 0.08**
- density = **0.98** to 1.08
- write verbose every $\Delta t$ = 1
- write movie every $\Delta t$ = 100
total = 18

#### Polydispersity fail for impossible fluid densities

Consider the freezing line of fluid-FCC:


In my sweep I take fluidensities from 0.98 to 1.08. The lower bound is because otherwise, bidispersity of the crystal is not conserved. The upper bound is based on the freezing of the fluid.

But as seen in the figure from Antoine's paper, for lower polydispersities 1.08 is already in the fcc-fluid regime of the phase diagram. So those points are not likely to give any good fluid-Laves coexistences.

**Most extreme example: p = 0.06, fluiddensity = 1.08:**

Of course, this is the most extreme example. But, at p = 0.14 and fluiddensity = 1.08, the same phenomenon may influence the measurements more subtly.

One place where this shows up, is the size distribution. The measured polydispersity no longer matches the set polydispersity if there is a phase transition to FCC going on. This makes sense and is not problematic.

**Most extreme example: p = 0.06, fluiddensity = 1.08:**

Also, this effect of not quite getting the right polydispersity is seen in p = 0.14 and fluiddensity = 1.08:

**More subtle: p = 0.14, fluiddensity = 1.08:**

In contrast, at actual fluid densities (like fluiddensity = 0.92), the polydispersity that is outputted corresponds very well to the one that was set.


(Okay, still not amazing, but maybe I should write more mov. snapshots, because this is based on only 9 measurements of 10.000 particles, because I did not want to make too much useless data.)

Additionally, I look at the spearman correlation between the size of a particle at different times. But I did not analyse this quantitatively yet.

For those fluid phase points that are perhaps wanting to form an FCC, you probably can't trust the pressure measurement either. So I have to watch out with what I call "possible coexistence". For example, for p = 0.08 these points found were with fluiddensity = 1.04, 1.06 and 1.08. But that is very unlikely to be right, since at those densities we expect a crystal.

Solution: Once I find my phase point that shows stable direct coexistence, I will rerun the fluid simulation to check if the polydispersity is proper at that fluiddensity.

### Direct Coexistence simulations

- CubicCu4x4x6 = 2304 particles
- Also, depending on both densities, roughly that amount also in fluid phase
- N = approx. 4000-5000
- 10000 simtime
- phase points (poly, fluiddens, dens, c = 1):
:::spoiler
"0.08 1.0400 1.051791 1.000011"
"0.08 1.0600 1.067873 0.999989"
"0.08 1.0800 1.083883 0.999971"
"0.10 1.0000 0.990915 0.999927"
"0.10 1.0200 1.004218 1.000039"
"0.10 1.0400 1.017049 0.999991"
"0.10 1.0600 1.030453 0.999988"
"0.10 1.0800 1.044624 1.000035"
"0.12 0.9800 0.944322 1.000031"
"0.12 1.0000 0.953331 1.000002"
"0.12 1.0200 0.962378 0.999992"
"0.12 1.0400 0.972416 0.999995"
"0.12 1.0600 0.984174 1.000042"
"0.12 1.0800 1.000016 0.999979"
"0.14 0.9600 0.896771 0.999975"
"0.14 0.9800 0.902952 1.000021"
"0.14 1.0000 0.911315 0.999987"
"0.14 1.0200 0.921468 1.000002"
"0.14 1.0400 0.933553 0.999995"
"0.14 1.0600 0.949971 1.000012"
"0.14 1.0800 0.975107 1.000011"
:::
- write verbose every $\Delta t$ = 1
- write movie every $\Delta t$ = 100
total = 21

#### Remaining composition split (test points)

While all direct coex simulations are running on the cluster I already have some shorter runs to work with on my laptop. I will show the difference between a stable-ish coexistence point and a non-stable coexistence point.

Stable:
poly = 0.14
fluiddensity = 1.00
crystal density = 0.9113
N = 4832

Non-stable:
poly = 0.14
fluiddensity = 0.98
crystal density = 0.9030
N = 4804

Both simtime = 2000
Both start looking something like this:

After t = 2000, unstable looks like:

While stable looks like:

More importantly, looking at the size distribution of both phases (rough estimate). Just take zcoord < zsizebox/2 and zcoord > zsizebox/2 for now.

Unstable shows little difference between the two sides:

While stable shows the expected behaviour of coexisting phases:

This is an example of the difference between meeting the stability requirement and not meeting it. The stable one keeps the Laves phase and the unstable one melts into a fluid.

Of course when I find the exact right phase point, I will make this composition comparison less rough.

### Direct coex (first round) results by eye and rough estimate like above

"0.08 1.0400 1.051791 1.000011" = fluid
"0.08 1.0600 1.067873 0.999989" = fluid
"0.08 1.0800 1.083883 0.999971" = crystal
"0.10 1.0000 0.990915 0.999927" = fluid 
"0.10 1.0200 1.004218 1.000039" = fluid
"0.10 1.0400 1.017049 0.999991" = crystal
"0.10 1.0600 1.030453 0.999988" = crystal
"0.10 1.0800 1.044624 1.000035" = crystal
"0.12 0.9800 0.944322 1.000031" = fluid
"0.12 1.0000 0.953331 1.000002" = fluid
"0.12 1.0200 0.962378 0.999992" = little fluid, lots of crystal
"0.12 1.0400 0.972416 0.999995" = crystal
"0.12 1.0600 0.984174 1.000042" = crystal
"0.12 1.0800 1.000016 0.999979" = crystal
"0.14 0.9600 0.896771 0.999975" = fluid
"0.14 0.9800 0.902952 1.000021" = fluid
"0.14 1.0000 0.911315 0.999987" = both, almost 50/50
"0.14 1.0200 0.921468 1.000002" = little fluid, lots of crystal
"0.14 1.0400 0.933553 0.999995" = little fluid, lots of crystal
"0.14 1.0600 0.949971 1.000012" = crystal 
"0.14 1.0800 0.975107 1.000011" = crystal 

### Pressure measurements

If the fluid density is too low, the crystal melts to become a liquid. From the pressure measurements, you see a good Pxx = Pyy = Pzz in that case that is **lower** than Pf. This is because $\rho_\chi < \rho_f$ and melting of the crystal lowers the global density and thus the (isotropic) pressure.

If the fluid density is too high, the fluid will crystalize. Because $\rho_f > \rho_\chi$, the particles in their crystal form will take up more space than in their liquid form. Due to the box size restriction in the x and y direction, the crystal will have to compress in the z-direction, leading to Pzz > Pxx,yy.

The right phase point should give Pzz = Pxx,yy = Pf.

p = 0.14, fd = 0.98 (too low)

p = 0.14, fd = 1.0 (sorta dead right)

p = 0.14, fd = 1.02 (too high)

#### Direct coex (second round)
Tested points:
:::spoiler
"0.08 1.0620 1.069456 1.000032"
"0.08 1.0640 1.071112 1.000002"
"0.08 1.0660 1.072813 1.000004"
"0.08 1.0680 1.074398 0.999997"
"0.08 1.0700 1.075943 1.000000"
"0.08 1.0720 1.077429 1.000015"
"0.08 1.0740 1.079207 0.999991"
"0.08 1.0760 1.080790 1.000002"
"0.08 1.0780 1.082383 0.999983"
"0.10 1.0220 1.005786 0.999972"
"0.10 1.0240 1.006723 1.000031"
"0.10 1.0260 1.008125 0.999997"
"0.10 1.0280 1.009417 1.000025"
"0.10 1.0300 1.010662 0.999998"
"0.10 1.0320 1.012038 0.999996"
"0.10 1.0340 1.013415 0.999993"
"0.10 1.0360 1.014695 1.000005"
"0.10 1.0380 1.016175 1.000063"
"0.12 1.0020 0.954042 1.000020"
"0.12 1.0040 0.955222 1.000004"
"0.12 1.0060 0.956313 0.999999"
"0.12 1.0080 0.956904 1.000017"
"0.12 1.0100 0.957678 0.999973"
"0.12 1.0120 0.958896 1.000015"
"0.12 1.0140 0.959664 0.999972"
"0.12 1.0160 0.960735 0.999956"
"0.12 1.0180 0.961555 0.999989"
"0.14 0.9880 0.905994 1.000006"
"0.14 0.9900 0.906871 1.000028"
"0.14 0.9920 0.907886 1.000017"
"0.14 0.9940 0.908682 0.999979"
"0.14 0.9960 0.909371 1.000000"
"0.14 0.9980 0.910438 0.999992"
"0.14 1.0000 0.911300 0.999987"
"0.14 1.0020 0.912161 0.999944"
"0.14 1.0040 0.913197 1.000010"
"0.14 1.0060 0.913751 0.999986"
"0.14 1.0080 0.915116 1.000010"
:::

Best found coex (match Pf and 50/50 of the box):

p = 0.08
fd = 1.064
P = 22.2
slowly crystalizing over time
metastable wrt FCC

p = 0.10
fd = 1.032
P = 19.5
interesting movie, very much melting and crystallizing at boundary
prob also metastable

p = 0.12
fd = 1.0120
P = 18.7
also fluctuates, but stays pretty good half/half!

p = 0.14
fd = 0.998
P = 18.2
good half/half

So around these points, I run a third round. Those coex, I will run longer, with different seeds (x5) to get better statistics.

:::info
From now on, I decided to only run c/a = 1, because it was quite obvious and it means 1/5 the runs.)
:::

#### Third Round of finding coex

5 runs
100.000 time eachs
11 fluiddensities for each of the 4 polydispersities

:::spoiler
"0.08 1.0626 1.070155 1.000000" full fluid
"0.08 1.0628 1.070129 1.000000"
"0.08 1.0630 1.070566 1.000000"
"0.08 1.0632 1.070727 1.000000"
"0.08 1.0634 1.070891 1.000000"
"0.08 1.0636 1.070930 1.000000"
"0.08 1.0638 1.070947 1.000000"
"0.08 1.0640 1.070950 1.000000"
"0.08 1.0642 1.071331 1.000000"
"0.08 1.0644 1.071674 1.000000"
"0.08 1.0646 1.071890 1.000000" full crystal, but which crystal?
"0.10 1.0310 1.011448 1.000000"
"0.10 1.0312 1.011471 1.000000"
"0.10 1.0314 1.011597 1.000000"
"0.10 1.0316 1.011648 1.000000"
"0.10 1.0318 1.011994 1.000000"
"0.10 1.0320 1.012119 1.000000"
"0.10 1.0322 1.012100 1.000000"

run 1: ends in crystal with droplet of fluid inside
run 2, 3, 4: ends in fluid

"0.10 1.0324 1.012275 1.000000"
"0.10 1.0326 1.012613 1.000000"
"0.10 1.0328 1.012747 1.000000"
"0.10 1.0330 1.012654 1.000000"
"0.12 1.0110 0.958584 1.000000"
"0.12 1.0112 0.958656 1.000000"
"0.12 1.0114 0.958260 1.000000"
"0.12 1.0116 0.958517 1.000000"
"0.12 1.0118 0.959060 1.000000"
"0.12 1.0120 0.959027 1.000000"
"0.12 1.0122 0.958758 1.000000"
"0.12 1.0124 0.959240 1.000000"
"0.12 1.0126 0.959246 1.000000"
"0.12 1.0128 0.959366 1.000000"
"0.12 1.0130 0.959347 1.000000"
"0.14 0.9970 0.910855 1.000000" more fluid but stable ceox
"0.14 0.9972 0.910005 1.000000"
"0.14 0.9974 0.910260 1.000000"
"0.14 0.9976 0.909980 1.000000"
"0.14 0.9978 0.910379 1.000000"
"0.14 0.9980 0.910512 1.000000"
"0.14 0.9982 0.910141 1.000000"
"0.14 0.9984 0.910233 1.000000"
"0.14 0.9986 0.910737 1.000000"
"0.14 0.9988 0.911062 1.000000"
"0.14 0.9990 0.910949 1.000000" stable coexistence
:::
### Checking trial fluids

Like mentioned in section "Polydispersity fail for impossible fluid densities", some phase points lead to irregular size distributions in the trial fluid simulations. Possible cause is for example FCC crystallization or something happening at very low pressures which I did not investigate further.

To see if this issue is relevant, I check the polydispersity (Gaussian fits) of the fluid simulations at the most promising coexistence fluid densities.

fit_polydispersity => 0.08051068289704127
fit_polydispersity => 0.09960695146091296
fit_polydispersity => 0.11973497264072354
fit_polydispersity => 0.13942283575723538

The Spearman correlation between a particle's size over time is also very weak, indicating a well-functioning fluid indeed.

### Results third direct coex run
I study the direct coexistences by:
- looking at the mov.
- looking at the Pii's over time
- looking at the Pzz-Pfluid intersect

Because I now have 5 runs, I can add error bars on my pressure measurements.

#### For p = 0.14

The coexistence is very stable. In the movie you clearly see that both the fluid and the crystal remain.

The pressure intersect shows where the true coexistence point is.

At this point, the pressure over time looks like this, for the 5 different runs:

Intersection pressure is at around 18.24 which is **lower** than FCC.

#### For p = 0.12
Similarly for p = 0.12. The volume fluid/crystal varies, but coexistence remains in all 5 runs. The pressure plots over time show how all runs behave similarly.


Intersection pressure is at around 18.62 which is **a little higher** than FCC.

#### For p = 0.10

This is a little weirder. First of all, the pressure curve looks funny.



This indicates that runs like run0 became full fluid, while runs like run1 became almost full crystal.

I say almost because there is a droplet of fluid in the crystal. (See movie snapshot below.)


#### For p = 0.08

Fully weird.


Interpretation of pressure graphs below (verified with mov)
run0: ??
run1&2: becomes a full crystal
run3&4: full fluid

I am going to make a BOPS analysis now, to see if FCC is intervening in run0. Perhaps we are looking at an FCC/Laves coexistence or something along those lines.
-> skipping this scavenger hunt for now, because it will not be stable anyways


screenshot of "??" simulation:

:::info
Note that in p = 0.08, melting is indicated by an overall pressure rise and crystallisation is indicated by an overall pressure drop. In the case of p = 0.12, melting is indicated by a pressure drop. This is not a problem but interesting to relate to density difference etc. we find later on.
:::

### BOPs
for fcc:
cutoff 1.4 gives **coordination number 12** and bops:
0.000000
0.000001
0.000000
**0.190941 = literature value**
0.000000
**0.574524 = literature value**
0.000000
0.403915
0.000000
0.012857
0.000000
0.600083

for CubicCu:
two types, both different BOP fingerprint
:::success
BOPs of small and large particles are quite different. Difference between average BOPs of Laves and FCC is more subtle.
:::

### Higher polydispersities

We may conclude from the p = 0.12 measurement, that this Laves coexistence is metastable. For p = 0.14 however, pressure indicates that Laves crystal is stable over FCC. We therefore drop all remaining questions about the p = 0.08 and 0.10 results. We want to measure even more extreme polydispersities: p = 0.16, 0.18 and 0.20. For this, we have to go back a few steps, to find the right externally enforces chemical potential function that brings about this potential.

:::info
Remember to increase the max number of particles per cell if you increase polydispersity.
:::

### Comparison my fits antoines fits

coefficients in $V(R)$ with my fits, for rho 0.99:
-0 24.463029758464 -1 -107.559821747802 -2 64.069041887501 -3 -69.988729136357
versus antoines:
-0 24.463029758464 -1 -108.507725266054 -2 66.025917328823 -3 -71.321590150951

### test points for 16 to 20, round 1

:::spoiler
"0.16 0.9800 0.858512 1.000000"
slowly, but goes to full fluid
"0.16 0.9820 0.859556 1.000000"

"0.16 0.9840 0.859401 1.000000"

"0.16 0.9860 0.861199 1.000000"

"0.16 0.9880 0.861596 1.000000"

"0.16 0.9900 0.862400 1.000000"
stays almost 50/50
"0.16 0.9920 0.863695 1.000000"

"0.16 0.9940 0.864492 1.000000"

"0.16 0.9960 0.865421 1.000000"

"0.16 0.9980 0.866424 1.000000"

"0.16 1.0000 0.867511 1.000000"
does not go full crystal (yet), but less fluid than crystal

from pressures: around 0.9865
so between 0.9860 and 0.9880 (but on the lower side of that) so next round: 0.9862 until 0.9874 with steps of 0.0002

"0.18 0.9600 0.805511 1.000000"
fluid
"0.18 0.9620 0.806464 1.000000"

"0.18 0.9640 0.806730 1.000000"

"0.18 0.9660 0.806421 1.000000"

"0.18 0.9680 0.807788 1.000000"

"0.18 0.9700 0.808936 1.000000"

"0.18 0.9720 0.809661 1.000000"

"0.18 0.9740 0.810309 1.000000"
looks very 50/50

"0.18 0.9760 0.810928 1.000000"

"0.18 0.9780 0.811139 1.000000"

"0.18 0.9800 0.813019 1.000000"

"0.18 0.9820 0.814184 1.000000"

"0.18 0.9840 0.815011 1.000000"

"0.18 0.9860 0.815478 1.000000"

"0.18 0.9880 0.817082 1.000000"

"0.18 0.9900 0.818721 1.000000"
blijf wel wat, maar niet heel veel fluid

pressure tells us it is around 0.9725
so in between 0.9720 and 0.9740 (on the low side)
second run: 0.9722 until 0.9734 with steps of 0.0002

"0.20 0.9400 0.753166 1.000000"
full fluid
"0.20 0.9420 0.753801 1.000000"

"0.20 0.9440 0.753945 1.000000"

"0.20 0.9460 0.753635 1.000000"

"0.20 0.9480 0.754564 1.000000"

"0.20 0.9500 0.754746 1.000000"

"0.20 0.9520 0.755830 1.000000"

"0.20 0.9540 0.756890 1.000000"

"0.20 0.9560 0.757527 1.000000"
pretty much 50/50
"0.20 0.9580 0.758753 1.000000"

"0.20 0.9600 0.758801 1.000000"

"0.20 0.9620 0.760852 1.000000"

"0.20 0.9640 0.761546 1.000000"

"0.20 0.9660 0.761479 1.000000"

"0.20 0.9680 0.763465 1.000000"

"0.20 0.9700 0.763737 1.000000"
more crystal than fluid

pressure tells us: around 0.9577
between 0.9560 and 0.9580
second run: 0.9568 until 0.9578 with steps of 0.0002

:::

### Nmaxneighbrous

test max neigbours (now made it 120 for safety, maybe let it print how many it is usually) = i checked it out now, and it seems to be fine to let it be 100, it never really comes up to be higher than 60 in a run of 100 tau for p = 0.20 at density = 0.99

### MgCu2 final results

After I found the "single phase - coexistence" intersect with satisfactory precision, I simulated a final trial crystal and fluid to just get all the features of the coexistence with more accuracy: like the packing fraction. To this end I printed more snapshots and measured the mean packing fraction after equilibration, etc.


## onto MgZn and MgNi

For these crystals, unit cell measurements are not 1: 1 :1.

Importantly, I turned the crystal generated from my Cubic.c code 90 degrees in the x,z plane. Such that the growth direction in an elongated box (enlongated in the z direction) would not be able to make a mixed laves phase. This means the dimers are nicely shown in the x,y plane: the interface in the direct coexistence simulation.

note that this also means that i vary a now, not c. So i get values of a_set, being slightly lower than 1.

### Bidisperse

To see if the unit cell is also deformed in the absense of polydispersity, I check out the bidisperse system with my EDMD code. I meaure the pressure tensor and see at what value of a/c (keeping b=c), the pressure tensor is isotropic. Results show that the (cubic) MgCu is not deformed, whereas the other 2 are.


Compare this to the polydisperse coexistence points:


### Spontaneous (first thoughts)

We want to check if we can use our relative chemical potential method to spontaneously form the crystals, with brute force simulation.

To start, I am running a long simulation at the exact coexistence phase point for MgCu that I found, only simulating the crystal. By this I mean, I set the density of the simulation to be that of the crystal in coexistence, but, the initial configuration is a fluid. The enforced relative chemical potential is set by the virtual fluiddensity given to the simulation, which I set to be the fluiddensity of the fluid in the true coexitence. This way we will see if MgCu pops up. The simulation is with 5000 particles.

4 nov: checked the evolution of the size distribution of the simulations: for no polydispersity the coex sim showed two peaks yet.

This was perhaps a silly attempt. The fluid stays at a much lower packing fraction (and pressure) than the coexistence pressure. In coexistence of course, the global number density/packing fraction is higher than this single phase. Right now, it seems like we find in coexistence with the fluid of {p, fd} to be another fluid with slightly higher mean size and lower number density. Surely, the absolute chemical potentials will not be equal, right? But question is if the simulation managed to get the deltamu(sigma) right.

For the ultimately found coexistence points $\rho_{\chi, coex}$, $\delta\mu_{coex}$ for each polydispersity, long run. gives:

$\eta_f < \eta_{spont} < \eta_\chi$
$P_{spont} < P_f = P_\chi$
$p_f < p_{spont}$

:::info
There is another fluid, with number density equal to coexistence crystal that:

has equal relative chem. pot.

but probably higher μN
and lower PV

#1 ? lower or higher F = (PV + μN)
#2 ? dynamically accessible
:::

If you enforce a dmu(sigma) corresponding to a higher fluiddensity, you basically push them towards larger sizes.

But if you make the actual density lower, you get this metastable thing, where a fluid just appears with a slightly higher mean size/polydispersity. With a lower density.

But, if you do coexistence simulation with a too high density, the system fully crystallizes. So that is hopeful that with starting with a higher fluiddensity this might appear.

Not as straightforward as it sounds though:
If the fluid partially crystallizes, the remaining fluid will have a higher number density, because the crystal has a lower number density (with higher mean particle size). This remaining fluid with a higher number density, is still under the same contraint of delta mu(sigma). Therefore, the mean particle size will become a little smaller. Note that in the nucleated crystal, the mean size is actually a little higher. So the chances of the crystal growing more are getting more and more slim. To be fair, the same is true if you start and keep a global size distribution, and even with the swap moves the Laves phases should spontaneously form.

### AB13

From Rinske's work we know: some crystal planes are more useful than others for the interface.

We prefer this orientation:

over this orientation:

Checked in the bidisperse case: AB13 is isotropic at a=b=c.

Check this once for the polydisperse case, but then switch to a=b=c style method.

### more phase diagrams

Aside from overall mean size, I extracted some more phase diagram data from the coexistence CRYSTAL and FLUID.

Namely, mean size of the large and the small particles. From fitting a smallGaussian + 0.5 * largeGaussian.

The ratio of these two:

And the polydispersity of the crystal: note, this is just the standard deviation/mean size. Not from a gaussian fit.

You can do this as a fraction of the fluid mean size, or of the cyrstal mean size:


Here I used fits of gaussians, and getting the 2 means of the double peaked size distribution. But, the two peaks do not have to be two gaussian per se: just look at the FCC fractionation. I switched to skewed Gaussian in the piece of plot_laves_per_type.py code for example. Next cycle, use these fits and get the mean from a skewed gaussian.

### Extreme standard deviation of p = 0.12 for Ni and Zn

This is because there are some simulations that fully turn into a fluid and some that stay coex.



These points not taken into account for the fits.

Note that I "zoomed in" further in the end, so these points should not matter anyways.

### Spontaneous (plan)

Going to write a (dmu, Ntot, p, T) ensemble code, where we use cluster volume moves. See paper "A cluster algorithm for Monte Carlo simulation at constant pressure".

What we'll do is, we'll start with a fluid of a pressure higher than the coexistence pressure in a cubic box. The corresponding number density sets the relative chemical potential such that the fluid will remain gaussian.

Crystallization will start, and since the crystal number density is smaller, we'll see that at constant pressure, the volume of the box increases. The fluid density will be kept the same fluid density as at the start.

Note that to get the pressure that corresponds to a certain fluiddensity and polydispersity (thus relative chemical potential function), we first have to measure the pressure in an NVT ensemble.

:::info
In order to be able to do this, you should be able to measure in the NVT ensemble at this density long enough such that you can measure the pressure. At very high densities, you'll get crystallization in the NVT as well, so this becomes more difficult.
:::

### More results: new phase diagrams

Additional phases:

We will add $\mathrm{AlB_2}$ to the phase diagram and go to even higher polydispersities to see if AB13 will win it.

It is interesting to note that AB13 formed spontaneously in the paper of Bommeneni at high polydispersities, but is not stable with respect to the Laves phases according to our results.

The question is, what difference in simulation leads to this? Is it the fact that they don't replenish the fluid to be have a Gaussian distribution (the fact that they have a different chemical potential function), or that fact that they work at constant packing fraction?

Is it possible that AB13 does not have a coexistence with a fluid, but is nonetheless more stable than the Laves phase at some pressures.

We know that AB13 is not stable if you only have swap moves. This seems reasonable, since it is difficult select such a small size ratio bidisperse mixture out of a Gaussian fluid distribution. But why does it form in the case where there is a general potential?

### Clustered volume moves (MC)

#### The reason you can not just do a cutoff that defines a cluster, but you can use a cutoff in the determination of the clusters as a first step

The reason the "cutoff" is now okay: because due to the probabilistic nature, you **can** still make the same cluster after a move, always. The fact that you can now make additional clusters is not a problem, and is taken into account in the acceptance probability.

#### Replacing sigma with sigma_ij

In the paper "A cluster algorithm for Monte Carlo simulation at constant pressure". They described the method for a monodisperse system. Any mention of $\sigma$ I can just replace with $\sigma_{ij}$ for my polydisperse case.

This includes the $\sigma$ in the set bond cutoff distance. If you would not do this, $b(r_{ij})$ can get negative etc. But it does mean that larger particles have a little larger cutoff shell, so they are possible to cluster with particles a little further away.

Discussed this with Laura and Frank, not an issue.

#### Volume or log volume space

We can either do a random walk through volume or through log volume space. Going for logV now, but should check whether this is actually best.

Additionally, the cutoff parameter should be optimized.

:::warning
I copied the value of the cutoff parameter wrongly from the paper. They find optimal value 0.025, I have been using 0.001! Fix this or at least look what the efficiency boost is.
:::

#### Changing $\Delta V$ during simulation
We aim to keep acceptance rate of volume moves around 35%, by changing the $\Delta V$ during the simulation (not too often though). This breaks detailed balance, but is negligible if you don't do it too often. (Says so in book of Frenkel page 56.)

#### Acceptance rule

For random walk in V space:
$A(V'|V;\chi_c) = \mathrm{min}\bigg[1, (V'/V)^{N_c} \,\mathrm{exp}\big[-\beta p (V'-V)\big]\frac{\omega_c(\chi_c|\mathbf{X'})}{\omega_c(\chi_c|\mathbf{X})}\bigg]$

For random walk in logV space:
$A(V'|V;\chi_c) = \mathrm{min}\bigg[1, (V'/V)^{N_c+1} \,\mathrm{exp}\big[-\beta p (V'-V)\big]\frac{\omega_c(\chi_c|\mathbf{X'})}{\omega_c(\chi_c|\mathbf{X})}\bigg]$

#### Using neighbourlist
In what circumstances can you use the neighbourlist to define which particles are close enough such that you might want to bond them into a cluster?

The only issue occurs if both particles are the maximum size and they are both at the edge of their neighbour list. Or, that they get to be that after the move. To solve this, we make the neighbourhood lists based on a hypothetical max particle size, that is maxparticlessize + cutoff.

For checking a new configuration for overlaps/omega factor, it is impossible to simply use the old neighbourlist. Calculating the new neighbourlist however, is not very efficient, chances are that you will have to reset it all anyways, if you reject the volume move. Instead, we **do** change the cell list, but only recalculate the neighbour list if we accept the move.

#### Check efficiency boost due to cluster algorithm

:::warning
You should check whether the cluster algorithm actually helps boost the simulation time / volume correlation time.
:::

#### Omega term extreme values

I noticed the omega term gets to be really big. I checked however if the values are still more often close to one, and this is the case. I made a little histogram of a collection of omega values and saw that the frequency of values omega (for omega > 1) and 1/omega (for omega < 1) nicely decayed with highest point at 1.

Explained why sometimes omega_term can get really big:

In that case omega_new >> omega_old, so perhaps omega_old is really small, because b(r12) was close to 1, but 1 and 2 were not bonded. In the new configuration, they could have moves further apart.

Explained why sometimes omega_term can get really small:

In the new configurations, particles can get to be really close together, while they were not that close in the old configuration. In that case b(r12') is really close to one, making the omega_term small.

:::warning
Is it sometimes too big for a double precision? Does that even matter then? The acceptance probability will be essentially either 1 or 0 in that case.
:::

### Semi-grand ensemble free energy

You can just view this as a Helmholtz free energy without interaction U, but them with some added U due to the external field. So free energy = $\sum_i U(\sigma_i) - TS$.

Should also do: legendre transform. Do you need a lagrange multiplier?

### Spontaneous formation results
Some results of the NPT formation simulations.

crystallizes in one of the two runs in ... simulation time

0.16:
dens = 1.04 NOT in 10^5
dens = 1.06 in 10^4

0.18:
dens = 1.04 in 10^5
dens = 1.06 in 10^4

0.20:
dens = 1.02 NOT in 10^5
dens = 1.04 in 10^4
dens = 1.06 in 10^4

Now I'll do runs of 10^6, this time 5 of each. Only for the top 2 non-crystallized densities. These pressure are between 20.8 and 25.0. I do this to see if the system crystallizes at pressures closer to the coexistence pressure (around 18).

Checking if the non-crystalizes run variants of a crystalized phase point are acting as a fluid: for example 18% at 1.04 crystalizes in run 0 but not in run 1. Based on the movement of the particles through the box, it does seem like it is just still a fluid at this time.

*Note that even though i might write NOT crystallizes, this doesnt mean that the size distribution isn't acting a little strange

After these lower pressure runs: I'll do some more runs to get many seeds different crystallization. Question: which phases are formed (implement phase identification for polydisperse systems) and perhaps who knows see if change in (local) size distribution is a precursor for nucleation.

### Size differences between small particles
We are testing trial crystals that we know are formed in some bidisperse systems. In our case, we have a polydisperse system and the size ratio's adjust themselves to the best possible value. In binary crystals, particles of the same size do not always have equal environments. This means they might prefer slightly different sizes. Clearly, this preference does not prevent them from forming a binary crystal, they are fine having the same size, but it does not have to be the optimum if they are allowed to cCheck higher polydispersitieshange.

This is the clearest in the case of ico-AB13. The 13 small particles that form the icosahedron, consist of 12 shell particles and 1 centre particles. The centre particle, turns out, has a slightly smaller mean size.


The size ratio of the two small particle types is around 0.904. Frank: "For reference, to get an icosahedron to be close-packed with all 12 shell particles touching both the center and their 5 shell-neighbors, the size ratio is sqrt(2+\phi)-1 ~= 0.9021 (with \phi the golden ratio). So pretty close to the ratio you are seeing here!"

We of course want to see if there is a similar effect in the Laves phases.

##### Laves:
The different types in the Laves phases are a little harder to distinguish.

##### MgCu
Should not have difference, I don't see them either.

##### MgZn and MgNi
On first view, we can identify two types by eye:

MgZn: 

MgNi: 

At first I thought dark blue and red would be the same, but as it will turn out, no!

I checked the environments by counting the neighbours are certain distances. Indeed the coloured particles are different in immediate surroundings. But, the distances of neighbours are not the only thing that can be different about an environment. The angles between neighbours could be different.

First, we might have to consider only the neighbour shell.
:::warning
Improve this analysis.
:::

Firstly however, we can just look for the differences we can spot for each different particle in a unit cell. Two or more might overlap, but at least differences should be clear.

The difference seems most prominent for MgZn.

However, we encounter something annoying when simulating the coexistence cyrstal for a long time:

##### MgZn

In MgZn polydispersity 12%, there is diffusion in the crystal. This prevents easy access to the mean sizes per site.

Normally going from frame 0 to last frame looks like this:

For p = 0.12 (and p = 0.18) it looks like this:

Due to this, the particles that start out small, eventually show a slight double peak in size distribution, because they take the large particles' places.

The same thing (but less extreme) goes for p = 18%. This messes up the green/red size ratio

If this is not an issue, you clearly see the different types of small particles have slightly different mean sizes. (So this is for p = 0.14, 0.16 and 0.20)


This is a very zoomed in picture of the Gaussians for each small particle in the unit cell (total of 16). Clearly, 4 are slightly bigger. The 4 purple lines are the green coloured particles. So indeed the different environments result in different mean size. That being said, the sizeratio here is 0.996-0.997.

-------------------------------
##### MgNi
In MgNi, things are steady across the board. So this changing of sites is not a problem.


Identifying the ones on the left side (persistent over all polydispersities), they are the dark blue particles in this unit cell:


The size ratio is around 0.998.

For one measurement (=polydispersity 0.20), all three types seem to be split.


:::warning
This is mostly reversed engineered. You should neatly look into the voronoi cell of all particles and see if this classification makes any sense from principles.
:::

### Going to 30% polydispersity

Try to shoot for measuring surface forces around 0.94-0.97. Hopefully there is not too much wiggling if you run them long enough in that case

For broad distributions of many particles, size values become negative. What we do instead is simulate smaller systems for longer at these polydispersities, such that we have less values with higher statistical value to fit from. This is a fluid we are talking about, so finite size effects should not be problematic for this.

### Spontaneous changes:
Made some changes in the NpT code:
- removed the temporal overlap check every so many steps: instead trust on the negative time check
- in determining overlap in the new configuration, do not allow for a small error in the overlap
- fixed the neighbour checking for after omega term to be cell list checking, only small differences in omega term, still expect previous results to be fine

For now:
- optimize the delta cutoff value
- just see whenever the max volume step is largest after optimization of that max volume step parameter: that is the cluster cutoff value that is best

:::warning
Do this more rigorously
:::

### AlB2
Tested for fd = 0.9 and 1.0
16-20 %
each their own crystal density

for 16%: no stable
for 18%: only when fd = 1.0 and d < 0.86
for 20%: only when fd = 1.0

value for a is around 1.01, so go to array of 1.00, 1.01, 1.02

For 18%: from coexistence simulations it shows that becomes fluid at and below fd = 1.02 in direct coex. simulations. At this density the pressure is already 1.02, so I am not going to try to find the right points for 16 and 18%.

The main issue for this is that I don't want to go near the glass regime when I am measuring my fluid pressures. Also, it will not be the most stable phase here anyways.

At 20%, I see it is around fd = 1.02 stable. This does indicate that at higher polydispersities this phase might become more stable and therefore I will pick this back up when I have my higher polydispersity surface force results.

We then know we should consider 20% - 30%.
For a value 1.0 1.01 1.02
Start with fluid densities around 1, check the crystal densities that match this pressure for each polydispersity.
At these crystal densities, vary the fluid density from perhaps 0.94 to 1.02 and check their coexistences (short).

### Using skewed gaussian as a fit on the particle sizes
Helps to seperate the two peaks better:


But, note the mean of a skewed gaussian.


where


For the laves phases, this is much nicer fit and I think more accurate estimation of the mean size of the small and large particles.

However, the mixed skewed Gaussian is not working great on the AB13 crystal. Luckily, the seperation between large and small is really big here, so I can just split them, fit two seperate skewed gaussians and take the mean based on a cutoff classification. Then also check that there is indeed a 1:13 ratio.

### Phase identification: Frank Kasper bonds

The three laves phases are part of a class called Frank Kasper phases. These phases all excist of a mixture of 4 FK polyhedra, that are the Voronoi cells of the particles. These are called Z12, Z14, Z15 and Z16 and are composed of 12 pentagonal and 0, 2, 3 and 4 hexagonal faces respectively.

Two particles have a Frank Kasper bond between them if they share 6 neighbours. These six neighbours dictate a hexagonal face on both those particles' Voronoi cell.

In the Laves phases, the small particles have a Z12 environment, while the Voronoi cell of a large particle is Z16. The large particles hence have 4 FK bonds, while the small particles don't have any. Other Frank Kasper phases can easily have much different type cells for each particles.

Each Laves phase will have a different way these FK bonds are structured. For example, I think in MgCu2, they form a cubic diamond network like carbon does. But there is also a hexagonal diamond variant.

**Plan of action:** From the perfect Laves phases, draw the FK bonds and notice the structure they form.
From the cyrstallized configurations, determine the FK bonds and draw them. This will help you visualise which Laves phase has formed.
Also, check to see if there are any Z14 or Z15 polyhedra present (you should see this by checking if any atom has 2 or 3 FK bonds). This would mean there is a whole new non-Laves FK phase appearing!
You will likely see a combination of Laves phases, will it be clear which is which from these bonds visualisation?
Ultimately, a quantization would be great.

:::info
Also look into the (un)supervised machine learning option.
:::

### Cuboctahedron

How would I find out if my ico-NaZn13 has become a cub-NaZn13? With cuboctahedron of 12 small particles around another small particles in the large particle holes.


Just checking one in the pure crystal phase simulations:

And also checked one in the coexistence simulation:

### NPT results

polydispersity = 0.16, density = 1.04, pressure = 23.86
4/5 v polydispersity = 0.16, density = 1.05, pressure = 25.10

polydispersity = 0.18, density = 1.02, pressure = 23.11
polydispersity = 0.18, density = 1.03, pressure = 24.35

polydispersity = 0.20, density = 1.02, pressure = 24.96
5/5 v polydispersity = 0.20, density = 1.03, pressure = 26.33

So need a pressure of 25 of higher. So that is coexistence pressure +7.

It appears that for the intermediately slow (p = 0.16 with d = 1.05 and p = 0.20 for d = 1.03) crystallisation, some nice coherent crystals are formed!


The max. 1 million time simulations of 1.04, 1.02 and 1.02 respectively are not showing any crystallisation yet.

It does indeed seem to be the case that the quick high pressure crystallisation there is more fluid phase left. Could be because simulations are just shorter, or more rapid crystallization prevents a more durable crystal phase.

Here, it seems like the grain boundary consist of extra small particles. Perhaps I even see this in the extra wide tail of the size distribution. (see below)


Drawing FK bonds:

For MgCu, this gives:

https://www.atomic-scale-physics.de/lattice/struk/a4.html

For MgZn, this gives: hexagonal diamond structure
https://www.atomic-scale-physics.de/lattice/struk/hexdia.html

For MgNi, it is a blend (not one in the link below i think) https://www.atomic-scale-physics.de/lattice/struk/carbon.html#sp3

### next points:
Running bigger systems:

polydispersity = 0.16, density = 1.0500, pressure = 25.10
polydispersity = 0.16, density = 1.0510, pressure = 25.23
polydispersity = 0.16, density = 1.0520, pressure = 25.36

polydispersity = 0.18, density = 1.0360, pressure = 25.13
polydispersity = 0.18, density = 1.0370, pressure = 25.27
polydispersity = 0.18, density = 1.0380, pressure = 25.39

polydispersity = 0.20, density = 1.0220, pressure = 25.22
polydispersity = 0.20, density = 1.0230, pressure = 25.36
polydispersity = 0.20, density = 1.0240, pressure = 25.51

I now also go to lower polydispersities, non-laves FK phases are probably more present there.

"0.08000000 1.09100000 25.17593131"
"0.08000000 1.09200000 25.30297595"
"0.08000000 1.09300000 25.41812434"
"0.10000000 1.08400000 25.15703918"
"0.10000000 1.08500000 25.26573731"
"0.10000000 1.08600000 25.39130171"
"0.12000000 1.06900000 25.14445181"
"0.12000000 1.07000000 25.27719500"
"0.12000000 1.07100000 25.41926011"
"0.14000000 1.06300000 25.12699892"
"0.14000000 1.06400000 25.26397400"
"0.14000000 1.06500000 25.38764398"

0.08 1.0930 25.41 (2) => FCC

Actually the box volume does not look very suddenly changed. I did not write a lot of snapshots per time, so can't yet look at it very nicely. Makes sense though, you can see how for FCC at 0.08 both number densities are very similar.

0.12 1.069 (1) => start of laves
0.12 1.070 (2) => Laves
0.12 1.071 (1) => Laves
0.12 1.071 (2) => start of Laves

0.14 1.064 (0) => laves
0.14 1.065 (2) => start of laves

... but more now

### Other Frank Kasper phases

In order to quickly spot any non-Laves FK phase, I plot a histogram of the number of FK bonds per particle. Note: systems with the same histogram fingerprint as the Laves phases are still possibly not a Laves phase.

Fluid:

Laves:

Interesting; there seems to be a dip at the CN = 13, even in a fluid. Mainly for the higher polydispersities. Maybe even in the fluid this forbidden coordination number is not often formed. This would mean there are some of the 4 actual polyhedra in the "fluid" already.

I checked the mean of all runs & densities (in red):
p = 0.20 

p = 0.08


We also see a FCC crystallize. There are clearly no FK bonds there.

### Michael's snapshot

Number density of the snapshot: 1.118832857

The distribution is not double peaked:


While in the paper it reads:


Where you do see a shoulder in there. Probably they used more than 1 snapshot for their distrubution, so I will just assume it is lack of precision.

If I place the snapshot in my EDMC simulation at constant volume, at polydispersity = 0.10 and number density 1.118832857 (of the snapshot), I find that it is stable.

The new size distribution is:

Note: for the comparison below, I used a chemical potential function that corresponds to a gaussian fluid of p = 0.10 and a fluid density equal to the crystal density. This means we are not looking at a coexistence point, only trying to make a fair estimation of which phase leads to a lower pressure under the same circumstances.

Putting the phase in my ensemble, I also read off the pressure of the snapshot: 30.04.

To compare, I ran a Laves phase simulation at the same pressure and number density. The pressure of MgCu2 at 0.10 and 1.118832857 gives P = 30.22. So the pressure does seem to be a little lower in the complex phase.

This difference in pressure is comparable to the difference in pressure we find between Laves phases freezing lines. So actually, checking the pressure of MgNi would be helpful here. (A little more involved because you need to optimize for the box measures. Measured it to be: 0.99697566) Interestingly, MgZn at this phase point results in a crystal pressure of 30.23.

There is a clear anisotropy of the pressure tensor in this simulation.


If I place the snapshot in an EDMD code, such that the particles no longer change size, pressure tensor is also anisotropic in a similar direction:

So indeed, the snapshot without any size changes also is strained.

Note that the pressure is lower in the EDMD because the packing fraction is significantly in the EDMC variant (probably, due to the large particles appearing).

Laves phases at this phase point:
0.6390462775
0.6387345625
Complex phase at this phase point in my field:
0.6273016
Complex phase initially:

Comparing the size distributions:

Pressure in the EDMD: 22.74814168516667
packing fraction in the EDMD: 0.61

Note: the size distribution of EDMD is the same of course as the initial snapshot (it did not look like this at first, but that was a binning difference).

The g(r) of the snapshot looks like this:


smoothed over, to find minimum, looks like this:


The histogram of number of FK bonds then looks like this:

Remember my crystallized Laves looked something like this:

And a fluid like this:
p = 0.20 
or
p = 0.08

Weirdly enough, particles with 1 FK bond is not actually expected in FK crystals. It does not correspond to one of the four FK polyhedra. So there being relatively many particles with 1 FK bonds is strange. Of course, they could be form the fluid. This is why I check:

In order to check which of these particles are in the crystal and which are in the fluid, I seperate them based on the mobility (for now just take the stdev of the positions x + y + z. Could only do this this simple becuase particles are not placed back in box.), like this:


That way, I can check the histograms of the fluid and crystal regions seperately.

Assumption herein is that in both our ensembles, the same particles are crystal and fluid.


or in one figure:


### New surface force measurements

p = 0.22
p = 0.24
p = 0.26
p = 0.28
p = 0.30

Even though these **should** all be of equal statistics, the high polydispersity, high density simulations are not getting more sharp.

I think this is due to glassy dynamics that prevent improvement of statistics of particle surroundings. (Frank agrees.)

For now, I will just try to make the best of it, fitting from the lower denisities for each polydispersity. I can check if the fit is good enough by just simulating some fluids and seeing if the size distribution is alright.

If that does not work, I could make a EDMD with swap moves, that should improve the statistics.

<!-- Using only the first 12, 10, 8, 6 and 6 densities for p = 0.22, 0.24, 0.26, 0.28, 0.30 respectively gives these fits:

p = 0.28

Clearly, the fit is not very nice for the higher densities. But also, can't be too sure about that. Lets compare it to the fit we get if we use all the data we have. Just also use the low statistic stuff, I get slightly different fits:
Now only have to check if the lower densities are not messed up by it. Seen from the other side:

In a way, I think this is better. I'll go some middle ground I think. -->

So my approach is: get as many of the densities involved, while keeping the fit of the lowest density reasonable. This I estimate to be optimal for 12,12,12,10,10

In the EDMC code I use for the fluid, I now set:
double maxparticleradius = 1.0;
double minparticleradius = 0.01;
Let's see how that goes.

The resulting size distributions of fluids of densities: 0.93, ..., 1.0 are shown below, with a Gaussian fit.


The places that were left empty are because of simulation crash. The crash was due to a "Local minimum found outside of considered particle size range", according to the error message. I did not look into the calculation much further, because from the dmu curve you can probably see what happens.

The breaking down of the code happens at the point where a maximum in the relative chemical potential function appears: there is an increasing part in the function of dmu. As written in the SM of Antoines paper, the code assumes the relative chemical potential function to be a monotonically decreasing function. So hence the breakdown. (I checked this with sympy function, not just looked at it by eye.)


However, importantly, you can see the size distribution already is not right before the point the dmu function shows a maximum. So we are going to have to fix the fit in the first place, not just allow for maxima in the V(R) in the EDMC code.

In any case, I made the mistake of thinking getting the fits at high as possible density would be the best approach, and extrapolating to lower densities would then be easy.

But this turns out not to be the case, partly due to the bad statistics.

Since I know I will probably be looking at lower fluid densities the higher I go in polydispersity, I can actually even say I care more about that regime. Luckily, this lower density regime is also the one that can more quickly give me nice surface force measurements.

So I started up some more surface force measurements simulations. Hopefully, this gives me better fits for lower density region, such that my fluids are of the right size distribution.

They already seem to be much quicker, so this might work out in time!

Note that it is still possible that I find the true chemical potential to have a maximum, and then we need to adjust the EDMC code still.

<!-- see if there is a chance they are just not yet equilibrated

Are they maybe crystallizing? Look at the movies. What I notice about simulation in CHECKFLUID of poly 0.30, density 0.92 (whose size distribution is way off): the big ones are barely moving. Actually that is not that weird, but does it mean that it is a jammed system? No, because the little ones move almost on a whole other timescale. Doesn't mean large ones are actually stationary, right? In our simulation, large ones do not have less possibility to grow than small ones, right?
-->

Result: going to lower densities for determining surface forces works well. Leave this for what it is.

### surface to surface distance ($s(r)$)
I wrote an additional analysis code, that makes neighbours based on surface-to-surface distance. From those neighbours, again you can determine the FK bonds in the system.

For a cutoff distance for the surface-to-surface distance, I do something similar to the minimum in $g(r)$. But normalisation for spherical shell volumes is more complicated (see below). This needs improvement.

The $s(r)$ neighbours work nicely for Laves:

You can see that the code picks out more FK = 4, which is the Laves crystal. So probably more of the Laves crystal is nicely detected.

It does not work as well for the Complex snapshot though:
The "complex cyrstal" consists of both FK = 4 and FK = 1, so it is not like there is more crystal detected. It just seems like there are less neighbours.

But, maybe I can improve the cutoff determination, currently not correcting for expected number of particles in a shell for a gas. This is harder to do because the normalisation depends on the two radii, as it is dependent on the centre-centre distance.

Fix this my normalising every "hit" right away, when you still know both particle radii.

So fix this normalisation and see if that fixed the FK bonds drawing in the complex snapshot without ruining the win you get in the Laves phase

Fixing the normalisation was nice in the sense that the cutoff indeed shifted a little. But, the complex phase histogram looked pretty much the same. With $s(r)$ came more FK = 0, less FK = 1, 2, 3, 4.

### also drawing particles with CN = 12

### effect of NVT simulation on the snapshot

Notice how the snapshot might be clearing up under the influence of the NVT simulation:

If there we actually more "quasicrystal" columns, you would also see an increase in #FK = 1. But 2 and 3 are really just the fluid. So it might be that the NVT is just making it such that there is laves growing out of the fluid region.

There are even some regions that seem like they have elongated the columns in the NVT ensemble. So perhaps i will turn on a quite long simulation just in the NVT ensemble if NPT fails.


#### No crystallization for p = 0.10

Funnily enough, no spontaneous crystallisation shows for 10% polydispersity. Even when I look at higher pressures.

Normally, all fluids crystallize around the value of P = 25. But even if I raise the pressure to P = 25.6, nothing happens for the p = 0.10 case.

#### Complex phase in NPT

I expect there to not be much crystallisation. The box will not likely shrink because it is filled with crystal already. So probably, it works as well in the NPT as in the NVT ensemble. But we can try it anyways.

Complex phase in NVT ensemble melted to fluid for number densities of 1.025 and lower.

From the NVT measurements, we know the packing fraction of the snapshot in our ensemble (original snapshot was 0.61).

DENSITY: 1.0500, PRESSURE: 21.694253493333
avg packing fraction 0.5980866233333333
DENSITY: 1.0750, PRESSURE: 24.505105380000
avg packing fraction 0.6107198133333334
DENSITY: 1.1000, PRESSURE: 27.569915486667
avg packing fraction 0.62081905
DENSITY: 1.1188, PRESSURE: 30.052352797949
avg packing fraction 0.6276589428205128

Only the 1.1188 seems to not fully have equilibrated, the other do seem like they are at a stable packing fraction.

This is nice, because it means that are measurements are around the packing fraction that they used in their ensemble. Granted, the size distribution is still different.

### the new fits, for the higher polydispersities


As you can see, the new fits work better than the ones I got from the higher densities. They are not perfect though.

For p = 0.22, 0.24 and 0.26, the shape seems all good for now, but you see that the fitted polydispersity is a little low. Option 1) this is because of the imperfect dmu fit. Option 2) this is a fault in the measurement and I should equilibrate longer etc. (unlikely).

In any case, I will go forward, and check the found fluid densities of my points to see if they are in the reasonable range. I don't want to be testing my coexistence with a slightly off polydispersity.

For p = 0.28 and p = 0.30, the shape itself also is distorted.

The reason being maybe that at p = 0.30 with the high densities the glassy dynamics start to show up again. The fit for p = 0.30 was made based on surface force measurements between density of 0.85 and 0.90. So clearly, the extrapolation is not very nice, on the high density side at least.

So again, I expect the coexistence points to be at densities where this dmu function works, but this is important to check.

### intermediate results in high polydispersity regime

- For next pure phases simulation: shift the crystal density range a little. (See especially high fluiddensity high p simulations.)

- AB13 from 28% and higher needs higher maxsize (now 1.5).

- NOTE: that for MgZn (and MgNi), the value of c is substantially different if you fit it from 3 points around zero wrt 5 points around zero. (0.99, 1.00, 1.01) or (0.98, 0.99, 1.00, 1.01, 1.02)
- So in final points, probably better do do all 5 here

Point of interest: AlB2 seems to have the unit cell measurements more strongly depend on the crystal density. Also, the a/c ratio goes from > 1 for p <= 0.22 to a/c < 1.0 for p >= 0.24. This is funny and nice to analyse further along other things that change like size ratio

## spontaneous formation

p = 0.10 makes sense that it crystallizes a little later becuase both fcc and laves are not so stable there

maybe that is why the snapshot phase appears here, because laves is not that stable

these are both formed from fluid:


this second one is form ./FORMATION/NPT/DATA/COMPLEX/POLY0.10/DENS1.1000/RUN0/FK1or4.sqs, almost full box length channel!

i'll look for more in the many other runs i have not checked out yet

### New direct coexistence data
In the first round of direct coexistence: some melting was a little slow, leading to skew in pressure, that pressure was not equilibrated out yet. But thats okay, because those points are not around the true coex value anyways. I set the eq steps to 5000

Found a region for stable AB13!

The fluid densities of the coexistence are not quite in the fully safe zone concerning their gaussian shape and polydispersity. For now, I decided to continue with the process because AB13 and Laves seem to be on the okay side of things. However, AlB2 is not going to be accurate.

### snapshot in NpT

works as well as in NVT, nothing too exciting

Please note: COMPLEX2 is actually from initial complex snapshot, COMPLEX (now POLY0.10ALSO) was a little mistake where is was intialised from a fluid.

I could not change the naming on Hera. So I removed all the misnamed files there. On PC, it is correctly named. (So also on the Odin backup.) If used to be that I accidentally initiated "COMPLEX" in NPT as a fluid, making it just another version of POLY0.10 with different seeding. Useful runs, though. That is why I kept it. But changing names with MobaXterm was disasterous.

Running Michaels snapshot at NPT did not make the phase cleaner, at least no more than in the NVT (makes sense).

## oS276
{12.6, 5.2, 4}

Putting Frank's snapshot in FK analysis: checks out

Symmetry in thermal crystal = perfect crystal, no visible difference.

note how tiling in the x,y plane has a neglicable effect on the free energy in the thermodynamic limit, due to the z direction degree of freedom

<!-- stable at 0.8, for p = 0.12 p = 0.14
stable at 0.9, for p = 0.14, almost for 0.12

broad scoop of fd's, might affect the pressure in the crystal

for p = 0.12:
0.8: P = 16.0
0.9: P = 17.0 (almost stable)

for p = 0.14:
0.8: P = 17.5
0.9: P = 18.2

p = 0.12 fd = 1.0: Pfluid = 18.5
p = 0.14 fd = 1.0: Pfluid = 18.2

p = 0.10 around fd = 1.02 and d = 0.7
p = 0.12 around fd = 1.01
p = 0.14 around fd = 1 and d = 0.9
p = 0.16 around fd = 0.985
p = 0.18 around fd = 0.97 -->

its melting at a too low fd, makes sense:
low fd: less intense potential: large particles favoured less: smaller particles on average: melts at lower crystal densities

! p = 0.14/0.16 fd = 1.10 is glassy, not a great size distribution
So you can't do that high fluid densities
Too high pressures anyways

#### Separating types for oS276

Because they are so clearly different, it is not strange they are separatable in size:


### COEXIstence of oS276

p = 0.08:
1.05; fully melt
1.07; crystal grows

p = 0.10:
1.03; fully melt
1.05; fully crystal
1.07; fully crystal

p = 0.12:
1.00; fully melt
1.05; fully crystal

p = 0.14, fd = 1.00: fully melts
p = 0.14, fd = 1.05: crystal grows

p = 0.16, fd = 1.00: pretty much 50/50
p = 0.16, fd = 1.05: crystal grows

p = 0.14: fd = 1.02, P = 20.4, crystal density somewhere between 0.92 and 0.95
p = 0.16: fd = 1.01, P = 20.6, crystal density somewhere between 0.88 and 0.92

fluid of p = 0.08 and p = 0.10 at even uptil fd = 1.07 look greatly gaussian

p = 0.08 and p = 0.10 are not stable at fd = 1.00

p = 0.08 at fd = 1.03, not stable
p = 0.10 at fd = 1.03, stable: crystal density 1.011, P = 19.4

p = 0.08 at fd = 1.05, stable: crystal density 1.058, P = 20.7
p = 0.10 at fd = 1.05, stable: crystal density 1.026, P = 21.4

p = 0.08 at fd = 1.07, stable: crystal density 1.075, P = 22.8
p = 0.10 at fd = 1.07, stable: crystal density 1.040, P = 23.5

Predicted fd values:

p = 0.08: 1.06
p = 0.10: 1.035
p = 0.12: 1.02
p = 0.14: 1.02
p = 0.16: 1.01

with pressures around 20

<!-- Probably, 0.10 is not going to work out. For the same reason it is difficult to get FCC or Laves. That's what makes it difficult to compare though, because it was created at 0.10. Perhaps for that exact reason? Laves is destabilized/hardly accessable and thus a metastable phase can appear. -->

write about different types

#### Hidden problem: determination of a/c

MgCu and AB13: fully symmetric
MgZn, MgNi and AlB2: vary a/c
oS276: vary both a/c and b/c

I have accidentally been using fits of order 3, on my interpolation of the a/c values for MgZn, MgNi and AlB2. By checking the pure crystal phase, I noticed that the pressure thereof is not sufficiently isotropic. Quite a relevant difference: 18.38 versus 18.53.

The estimate is likely still in the right ballpark. Therefore I can probably use the same pure phases, but just run a new final direct coexistence measurement.

With linear fit, the difference is still not zero, but 0.04 instead of 0.2. So I like it better.

#### To do's phase diagram points:

MgCu:
Only high polydispersities: final accuracy at multiple runs
- [ ] wait: coex, multiple runs

MgZn and MgNi:
- [ ] wait: the right coexistence points at 0.0005 accuracy for high polydispersity
- [ ] determine range for 0.0001 accuracy thereof
- [ ] run pure crystal and fluid for the high polydispersities (0.0001 accuracy)
- [ ] do multiple runs for high polydispersities
- [ ] redo the multiple runs for low polydispersities, now with the linear fit (just hope the range is okay)
- [ ] find final stable coexistence, check that crystal for isotropy!!

<!-- Only have to run fitting_c.py code to find new coex. points. Run these direct coexistence simulations.
- For low polydispersity regime: 5 a/c points
- For high polydispersity regime: 3 a/c points
- Do this for a single run, not yet multiple runs
- Check to see if indeed, the difference between xx and yy is now smaller. -->

AB13:
All polydispersities: more accurate, multiple runs
- [ ] wait: coex, multiple runs for high polydispersity
Note: did not go to higher accuracy for AB13 at low polydispersity due to metastability

AlB2:
- [ ] wait: for direct coexistences, determine pressures for in phase diagram
- [ ] Analyse unstrained unit cell measurements for AlB2 and the change in sign of (a/c-1). Compare with size ratio changing. Check if this is still the case with these new fits.

oS276: all polydispersities
- [ ] wait: for direct coex
- [ ] more accurate pure phase points
- [ ] coex higher accuracy
Note: accuracy of 0.0005 is good enough, since it is metastable
- [ ] find final stable coexistence, check that crystal for isotropy

General, after:
- Run crystal points to find all the accurate things for in the phase diagrams
- Maybe try to get laves for 0.08, for a metastable point
- Figure out if 10% is really that impossible
- check what the number density difference was between fluid and crystal for the p = 0.10 case. For hypothesis of why oS276 is formed over Laves at 10% polydispersity. (Note: we actually find more Laves at 10% than oS276. Perhaps the reason should more so be found in )
- actually could we try harder to find meta stable laves in 10%?

### This week: Fixing the Laves phases QC, and writing

- process last bits of phase diagramm incl. higher p and oS276
- running the final pure phases
- while waiting, read QC stuff
- QC meeting (Monday)
- clean up hackMD, perhaps start a new file for cleaner results

### NEXT WEEK:

- once all phase diagram is finished: show the worst gaussian shapes that are part of a stable coexistence
- once you find proper area for oS276: check if longer runs of crystal are actually fully stable
- fix mobaxterm
- clear some data off PC
- clear off some data from hera
- is it worth it to delete all the .c files from the data?

#### Later:

NPT:
- fix whatever bug is caused in the higher polydisperse NPT code
- check to see if NVT is done and what the size distribution looks like, at these higher polydispersities
- turn on NpT at high polydispersities, see if AB13 is formed spontaneously

? - In paper it says: "Only coordination numbers (CN) 12, 14, 15, and 16 occur in oS276, which means it is a FK phase [19]." So does that mean they are never seeing #FK = 1 or just that the other side of that particle does not look isocahedral.

Formation optimization:
- optimize for the del parameter
- analysis of the speed up due to volume cluster moves (CPU or simulation time)

Find coexistence point with fluid in the glassy regime? With the EDMC code, we have the extra freedom, so this will not be quite visible. But in the EDMD we clearly see glassy dynamics from a certain density onwards.

overlap in spontaneous??

Read up on alternative for FK bonds: machine learning like in laura's paper.

Note: Calling it chemical potential still kind of bugs me a little bit. because in chemistry, we like to say that if things can react into one another, they have to have equal chemical potential. Maybe what we call the enforced chemical potential is more like a $\mu_0$. Where then the configurational/density dependent part should be added so that all species indeed have an equal chemical potential. Measure relative chemical potential in the crystal. (might have a problem with the mixing contribution, different anyways in a crystal? indistuinguishability?)

check out NaZn-cuboctahedral, should we add it?

check out Wahnstro ̈m Lennard-Jones mixture

check out all the references about alternatives in the semigrand ensemble, both SG simulations and free energy calculations & why do they struggle at high polydispersity

compare chemical potential function with nikki

see if bidispersity is a precursor for crystallization

do the ensemble potential derivation, lagrange multiplier

for my sanity, check the lowest density FLUID simulations, whether their polydispersity is nicely Gaussian

Write_movie is triggered at very end = bad for statistics

look for supplementary information in brute force paper, look for the Frank Kasper phase

check if $V(R)$ fit was at approx. the right fluid densities

at some point you are going to need estimation of the error in phase diagram

Fix:
- do #runs, not on which nodes