# Nucleation of repulsive hard-core Yukawa particles

###### tags: `PhD projects` `nucleation`

---

### Owners (the only one with the permission to edit the main text, others can comment)

Marjolein, Laura

---

## To do

- [ ] Write paper on barriers.
- [ ] Write paper on the study of the fluid before nucleation of the FCC crystal.
- [ ] Maybe use the inherent structure to say something more about the fluid before nucleation of the BCC crystal?

---

## Background

### Repulsive hard-core Yukawa

The interaction potential between colloids of diameter $\sigma$ and charge $Ze$ in a solvent with dielectric constant $\epsilon_s$ is given by
$$
\beta\phi(r) = \begin{cases}
\beta\epsilon \; \frac{e^{-\kappa\sigma (r/\sigma -1)} }{r/\sigma} & \quad \text{for } r\geq\sigma,\\
\infty & \quad \text{for } r<\sigma,
\end{cases}
$$
where $\kappa$ is the inverse Debye screening length and $\epsilon$ is the value of the potential at contact given by
$$
\beta\epsilon = \frac{Z^2\lambda_B/\sigma}{(1+\kappa\sigma/2)^2},
$$
with $\lambda_B=e^2/\epsilon_s k_BT$ the Bjerrum length, $k_B$ the Boltzmann constant, $T$ the temperature, and $\beta=1/k_BT$.

The hard-core Yukawa particles form an FCC crystal for strong screening and a BCC crystal for weak screening. For example, see the phase diagram for $\beta\epsilon=20$, taken from **[Hynninen & Dijkstra (2003)]**.

### Crystal nucleation
The nucleation barrier according to CNT (classical nucleation theory) is given by
$$ \Delta G(R) = 4\pi\gamma R^2 - \frac{4}{3}\pi |\Delta\mu | \rho_s R^3, $$
where $\gamma$ is the fluid-solid free-energy surface density, $|\Delta\mu |$ is the difference in chemical potential between the fluid and solid phases, and $\rho_s$ is the density of the solid phase. The height of the nucleation barrier is given by
$$ \Delta G^* = \frac{16\pi \gamma^3}{3|\Delta\mu|^2 \rho_s^2}
$$

#### Umbrella sampling
To compute the nucleation barriers, we use umbrella sampling. We perform the umbrella sampling using the biasing potential
$$
\beta U_\text{bias}\left(n(\mathbf{r}^N) \right) = \frac{\lambda}{2} \left( n(\mathbf{r}^N) -n_c \right)^2,
$$
where $\lambda$ is a coupling parameter, $n(\mathbf{r}^N)$ is the size of the largest cluster present in the configuration $\mathbf{r}^N$, and $n_c$ is a target cluster size.

#### Classifying a particle as solid-like or fluid-like
To determine whether a particle is fluid-like or solid-like, we make use of the 6-fold bond orientational order parameters
$$
q_{6m}(i) = \frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b} Y_{6m}\left( \theta_{ij},\phi_{ij} \right),
$$
where $N_b(i)$ is the number of neighbors of particle $i$, $\,\mathcal{N}_b$ is the set of neighbors of $i$, $\,Y_{6m}\left( \theta,\phi \right)$ are the spherical harmonics with $m\in[-6,6]$, and $\theta_{ij}$ and $\phi_{ij}$ are the polar and azimuthal angles of $\mathbf{r}_{ij}=\mathbf{r}(j)-\mathbf{r}(i)$.
The number of connections $\xi(i)$ of particle $i$ can then be determined via
$$
\xi(i) = \sum_{j\in\mathcal{N}_b} H\left(d_6(i,j) -d_c \right),
$$
where $H$ is the Heaviside step function, $d_c$ is the dot-product cutoff and $d_6(i,j)$ is the dot product given by
$$
d_6(i,j) = \frac{ \sum_{m=-6}^{m=6} q_{6m}(i) q_{6m}^*(j) }{ \left(\sum_{m=-6}^{m=6} |q_{6m}(i)|^2 \right)^{1/2} \left(\sum_{m=-6}^{m=6} |q_{6m}(j)|^2 \right)^{1/2} },
$$
with $^*$ indicating the complex conjugate. The neighbors of particle $i$ are defined as all particles $j$ with $|\mathbf{r}_{ij}|<r_c$ and a cluster contains all solid-like particles that have a solid-like neighbor in the same cluster.

A particle is labeled as solid-like when $\xi\geq \xi_c$, with $\xi_c$ a cutoff value.

#### Fitting the barrier to CNT
We fit our obtained nucleation barriers to CNT by assuming that the measured radius, $R(\xi_c)$, differs from $R_\text{CNT}$ with a constant shift $\alpha(\xi_c)$
$$
R(\xi_c) = R_\text{CNT} + \alpha(\xi_c).
$$
Then, using that the cluster size $n(\xi_c)$ can be related to $R(\xi_c)$ by
$$
n(\xi_c) = \frac{4\pi R(\xi_c)^3 \rho_s}{3},
$$
we obtain the equation
$$
\beta\Delta G(n,\xi_c) = 4\pi\beta\gamma \left[ \left( \frac{3n(\xi_c)}{4\pi \rho_s} \right)^{1/3} - \alpha(\xi_c) \right]^2 - \frac{4\pi}{3}\beta|\Delta\mu | \rho_s \left[ \left( \frac{3n(\xi_c)}{4\pi \rho_s} \right)^{1/3} - \alpha(\xi_c) \right]^3 + \delta,
$$
where $\delta$ is an added constant shift. This shift is added for the reason that the nucleation barrier obtained via umbrella sampling is only expected to match CNT near the top of the barrier.

#### Nucleation rate from the barrier
From the nucleation barrier, the nucleation rate $k$ can then be estimated using
$$
k \approx \rho f_{n^*} \sqrt{\frac{| \beta \Delta G''(n^*) |}{2\pi}} e^{-\beta\Delta G(n^*)},
$$
with $n^*$ the number of particles in the critical nucleus, $\rho$ the density of the supersaturated fluid, $f_{n^*}$ the rate at which particles are attached to the critical nucleus, and $\Delta G''(n)$ the second the derivative of the free-energy barrier.

The attachment rate $f_{n^*}$ can be determined from the mean squared deviation (MSD) of the cluster size at the top of the barrier
$$
f_{n^*} = \frac{1}{2} \frac{\langle \Delta n^2(t)\rangle}{t},
$$
where $\Delta n^2(t) = \left( n(t) - n^* \right)^2$.

### Nucleation rate using brute force
The nucleation rate can also be calculated by calculating the average time $\langle t \rangle$ it takes a system of volume $V$ to form a critical nucleus
$$
k = \frac{1}{V \langle t \rangle}.
$$
To find $\langle t \rangle$, one simply has to run numerous brute force simulations (MD or kMC) and wait for them to nucleate.

---

## Nucleation barriers

### Nucleation of the FCC crystal of "too soft" particles
These are some of the results I got for my Master's thesis.
Underneath, a table containing the calculate freezing and melting packing fraction of various hard-core Yukawa systems. The mapping to hard spheres was performed using the freezing packing fraction (with $\eta^{HS}_F=0.494$).


Underneath, a table giving the pressures at which the nucleation barriers where computed and some variables of the nucleation barriers gotten from fitting to CNT. The umbrella sampling was performed using $0.01 \leq \lambda \leq 0.06$ and particles were classified as solid or fluid using the cutoff values $d_c=0.7$ and $\xi_c=6$.


Underneath, the nucleation barriers for the various systems obtained via umbrella sampling ($\xi_c=6$). Although the barriers are clearly effected by the finite size of the simulation box ($N=2048$), one can clearly see that the barriers are mostly dependent on $1/\kappa\sigma$ and only weakly dependent on $\beta\epsilon$. Furthermore, as the barriers of even the "least soft" particles (those with $1/\kappa\sigma=0.05$) are approximately 1.5 times higher than the hard-sphere barrier, we conclude that the studied hard-core Yukawa particles are "too soft" to be considered as hard spheres.


I only calculated the attachment rate $f_{n^*}/D_0$, with $D_0=\Delta r^2_{max}/6$ the free diffusion coefficient, for the barriers with $1/\kappa\sigma=0.05$, as the other barriers are clearly influenced by the finite size of the system. The table underneath gives $f_{n^*}/D_0$ and the resulting nucleation rate $k\sigma^5/D_0$.

---

### Nucleation of the FCC crystal of "less soft" particles
To try to find when soft particles become "too soft" to be considered as hard spheres, we calculate the nucleation barriers (and rates) of hard-core Yukawa particles with $1/\kappa\sigma=0.04,0.03,0.02,$ and $0.01$. To give an idea of the "softness" of these particles, underneath you can compare the soft repulsive tail of these particles with that of WCA particles (for which the nucleation barrier/rate maps well to hard spheres).


#### Freezing and melting packing fractions & mapping to hard spheres
Underneath, the freezing and melting packing fractions of various hard-core Yukawa systems. The uncertainty of the packing fraction is 0.001 and the effective hard-sphere melting packing fraction is calculated using the freezing packing fraction $\eta^{HS}_F=0.492$.

Underneath, values of the various hard-core Yukawa system mapped to the hard-sphere system at $\beta P\sigma^3=17.0$.

#### Umbrella sampling with $N=2048$
Underneath, the resulting nucleation barriers of the various hard-core Yukawa systems obtained via umbrella sampling with $N=2048$ particles ($\xi_c=6$, $\lambda=0.020$). We see that the barriers with $1/\kappa\sigma=0.01$ and $\beta\epsilon=8$ or $39$ resemble the hard-sphere barrier. So these could be considered as "hard" spheres. However, that the other barriers with $1/\kappa\sigma=0.01$ are so much higher is strange. Futhermore, the barrier of $\beta\epsilon=39$ and $1/\kappa\sigma=0.03$ is abnormally high.
To check these strange behaviors, I performed the umbrella sampling again with other variable, i.e. different $\lambda$, different number of MC cycles between the application of the biasing potential, $r_c$ rounded on two decimals into of one, but the strange behavior of thse barriers remained. Hence, I believe that the problem is the finite size of the system.


#### Umbrella sampling with $N=10976$
Nucleation barriers obtained via umbrella sampling of systems with $N=10976$ particles ($\xi_c=6$, $\lambda=0.020$). The barriers are not that different from the barriers obtained with $N=2048$.


The barrier of $\beta\epsilon=39$ and $1/\kappa\sigma=0.03$ is still abnormally high. I discovered that this was caused by an equation of state that was slightly off.
From the **old** equation of state: $\eta_F=0.357$ and $\beta P \sigma^3 (\eta^*)=12.4$.
From the **new** equation of state: $\eta_F=0.361$ and $\beta P \sigma^3 (\eta^*)=13.2$.

To be sure that the other equation of state are correct, I redid all the simulations for the equations of state (making sure to let them initialize longer and run them for longer). Furthermore, I improved the fitting of the equations state. The new freezing and melting packing fractions are shown underneath, as well as the new barriers and their information. These barriers were obtained by averaging over four different umbrella simulations each starting with a different configuration.




It is important to note that, even though the difference between the "old" and "new" freezing packing fractions is small, i.e. for most state points less than 0.001, the difference between "old" and "new" pressures can be quite significant. This demonstrates that this method of mapping to hard spheres results in a relatively significant uncertainty in the pressure and consequenly a significant uncertainty in (the height of) the nucleation barrier.

As for the "too soft" spheres, the barriers of these "less soft" spheres are mostly dependent on $1/\kappa\sigma$. The ones with $1/\kappa\sigma=0.01$ are equal to or even lower than the hard-sphere barrier.

In the above table, the interfacial free-energy density $\beta\gamma\sigma^2$ is obtained by fitting the nucleation barrier to CNT.

I calculated the attachment rate $f_{n^*}/6D_l$ for some of the barriers. $f_{n^*}/6D_l$ and the resulting nucleation rate $k\sigma^5/6D_l$ are given in the table underneath.


#### Nucleation barriers at fixed supersaturation
From the results of the previous section, we see that the barriers that map well to the barrier of hard spheres (those with $1/\kappa\sigma=0.01$) all have a supersaturation almost equal to the supersaturation of the hard-spheres system. Furthermore, when we increase the softness, i.e. increase $1/\kappa\sigma$, we see that the supersaturation decreases. Thus, maybe comparing systems at equal supersaturation will be a better method for comparing to the hard-spheres system?

To investigate whether this is the case, we calculated some more nucleation barriers for $\beta|\Delta\mu|\approx0.536$ (corresponding to hard spheres at $\beta P\sigma^3=17.0$) and $\beta|\Delta\mu|\approx0.585$ (corresponding to hard spheres at $\beta P\sigma^3=17.5$). The resulting barriers with their information are given below, with a) $\beta|\Delta\mu|\approx0.536$ (first block of rows in the table) and b) $\beta|\Delta\mu|\approx0.585$ (second block of rows in the table).

Again, we see that $\beta\epsilon$ does not have a large influence on the nucleation, whereas increasing $1/\kappa\sigma$ now decreases the height of the nucleation barrier.



---

### Nucleation of the FCC crystal but with WCA cutoff?
We maybe have this feeling that the reason why these nearly-hard hard-core Yukawa particles do not map well to hard spheres is because of still relative long-ranged tail of the interaction potential compared to that of the WCA potential. So for the WCA potential the interaction is cutoff at $r/\sigma=2^{1/6}\approx1.122$. Underneath we compare the WCA potential with $\beta\epsilon=40$ to that of hard-core Yukawa particles with $\beta\epsilon=39$ for different values of $1/\kappa\sigma$. Note that for the bottom row we scaled $r$ with the effective hard-sphere diameter obtained from mapping the freezing packing fractions (see next section).

| "normal" hard-core Yukawa| hard-core Yukawa truncated <br> at $r/\sigma=2^{1/6}$ and shifted |
| -------- | -------- |
|  |  |
|  |  |

#### Mapping with $\eta_F$
I redid all simulations for the equation of states, Einstein integration, and thermodynamic integration for $1/\kappa\sigma=0.02$ and $1/\kappa\sigma=0.04$. Underneath are the results for the freezing and melting packing fraction, as well as the effective melting packing fraction $\eta^\text{eff}_M = (\eta^\text{HS}_F / \eta_F) \, \eta_M$ and effective hard-sphere diameter $\sigma_\text{eff} /\sigma = (\eta^\text{HS}_F / \eta_F)^{1/3}$. The last three columns give, for state points used for the new umbrella sampling, the packing fraction of the supersaturated fluid $\eta^*=1.088\eta_F$ and corresponding pressure and supersaturation. To compare, the numbers for the hard-spheres system, WCA, and hard-core Yukawa system with "normal" cutoff are given as well. Here with "normal" cutoff I mean when the error/shift of the interaction potential is less then $10^{-5}k_BT$. The numbers for the WCA system I got from [Laura's paper]( https://doi.org/10.1063/1.3572059).


Underneath are the barriers for these systems with a WCA cutoff for the potential. To compare, the original barriers (with normal cutoff) are shown in yellow and red, and the corresponding new barriers (with WCA cutoff) are show in dark yellow and dark red, respectively.
We can see that using the WCA cutoff brings to barrier of $1/\kappa\sigma=0.04$ down from a height of $25.7k_BT$ to $20.2k_BT$, which is significantly closer to the hard-spheres barrier. In contrast, however, the barrier of $1/\kappa\sigma=0.02$ increases from $20.7k_BT$ to $22.5k_BT$ when using the WCA cutoff.


#### Other mappings?
Underneath I compare the supersaturation as a function of two ways of mapping to the hard-sphere system. The first row shows $\beta|\Delta\mu|$ as a function of $\eta/\eta_F-1$, and the second row as a function of $(\eta-\eta_F)/(\eta_M-\eta_F)$.
Note that the first maps the freezing packing fractions on top of each other, which does not necessarily means that the melting packing fractions are mapped on top of each other. To illustrate this: for hard spheres $\eta_M/\eta_F-1=$ 0.105, whereas for hard-core Yukawa particles it is 0.101, 0.093, 0.87, and 0.78 for $1/\kappa\sigma=$ 0.01, 0.02, 0.03 and 0.04, respectively.
The second method maps the entire freezing to melting region on top of each other. *Note that, for the figures of the second row, I fixed the range of the horizontal axis such that it corresponds to the range of the hard-spheres data of the figures of the first row.*

Comparing the two rows, we see that mapping the entire freezing to melting region generally lowers the supersaturion of the hard-core Yukawa systems compared to mapping just the freezing packing fractions. This means that mapping the entire freezing to melting region is an even worse method than mapping just the freezing packing fractions, as it would result in even higher nucleation barriers for the hard-core Yukawa systems.

Lastly, I wanted to mention that we are NOT looking for a mapping that produces the same supersaturations as for the hard-spheres system. In a previous section, I showed the nucleation barriers for fixed supersaturation and their heights decreased with increasing $1/\kappa\sigma$.

| normal cutoff | compare with WCA cutoff |
| -------- | -------- |
|  |  |
|  |  |

What if we look at the supersaturation as a function of the pressure? Underneath the first row shows the supersaturation as a function of $\beta P \sigma_\text{eff}^3$, where $\sigma_\text{eff}^3$ was determined by mapping the freezing packing fractions. The second row shows the supersaturation as a function of $P/P_F$, where $P_F$ is the pressure at freezing.
*Note that I again fixed the range of the horizontal axis such that it corresponds to the range of the hard-spheres data of the figures above.*

| normal cutoff | compare with WCA cutoff |
| -------- | -------- |
|  |  |
|  |  |

And what about the equation of state itself? Underneath the scaled equations of state and the scaled derivatives of the equiations of state.

| Equation of state | Derivative|
| ----------- | -------- |
|  |  |

#### Structure of the fluid
We wanted to find out if there is any structural difference between the fluid of the hard-spheres system, WCA system ($\beta\epsilon=40$), and different hard-core Yukawa systems (both with the "normal" cutoff and WCA cutoff). To this end, I calculated the $g(r)$ and $S(q)$ of the fluid at the freezing and melting packing fraction of the different systems. Underneath are many firgures of the $g(r)$ and $S(q)$ which try to highlight the similarities and differences. The LEFT column always shows the results at the freezing packing fraction and the RIGHT column at the melting packing fraction. Note that I used $\eta_F$ and $\eta_M$ reported in the table above, but that I rounded them on the third decimal as the error in the reported $\eta_F$ and $\eta_M$ is approximately 0.001. Furthermore, note that I did my best to have the same axes range for the left and right figure on a single row.
Lastly, note that when I use the length scale $\sigma_\text{eff}$, I mean the effective hard-sphere diameter obtained from mapping the freezing packing fractions: $\sigma_\text{eff} /\sigma = (\eta^\text{HS}_F / \eta_F)^{1/3}$. These are also reported in the table above.

Let us look at the **radial distribution function** first. For all systems I calculated the $g(r)$ using a binsize of $0.02\sigma$.
Underneath I compare the resulting $g(r)$ of the hard-spheres system, WCA system, and different hard-core Yukawa systems with $\beta\epsilon=39$ (and normal cutoff). Going from top to bottom, I show the $g(r)$ as a function of $r/\sigma$, the $g(r)$ as a function of $r/\sigma_\text{eff}$, the first peak as a function of $r/\sigma_\text{eff}$, and the second peak as a function of $r/\sigma_\text{eff}$.
In general, I would say that scaling with $\sigma_\text{eff}$ works really well. The most significant differences between the "practically hard" systems (i.e. HS, WCA, and $1/\kappa\sigma=0.01$) and "slightly softer" systems (i.e. $1/\kappa\sigma>0.01$) can be seen in the broadening (and subsequent lowering) of the first peak both at $\eta_F$ and at $\eta_M$. Notice, however, that for the second peak, there are only sigficant differences at $\eta_M$ and not at $\eta_F$.
*Note: to get a smooth curve for the first peak of the $g(r)$ for the figures underneath, I used a binsize of $0.005\sigma$ for $r\leq1.08\sigma_\text{eff}$ to calculate the $g(r)$, and a binsize of $0.02\sigma$ for $r>1.08\sigma_\text{eff}$.*

| at $\eta_F$ | at $\eta_M$ |
| -------- | -------- |
|  |  |
|  |  |
|  |  |
|  |  |

Now, lets compare with the $g(r)$ of the hard-core Yukawa systems with a WCA cutoff. Both for $1/\kappa\sigma=0.02$ and $1/\kappa\sigma=0.04$, I would say that using the WCA cutoff brings the $g(r)$ closer to that of HS or WCA. With this I mean that the broadening of the first peak is less than for the system with a normal cutoff. The difference in the second peak at $\eta_M$ is less pronounced as well.

| at $\eta_F$ | at $\eta_M$ |
| -------- | -------- |
|  |  |
|  |  |
|  |  |
|  |  |

Next, we turn our attention to the **structure factor**. I calculated $S(q)$ by integrating the $g(r)$. Again let's start with comparing the hard-core Yukawa systems with the normal cutoff to HS and WCA. From top to bottom: $S(q)$ as a function of $q\sigma$, $S(q)$ as a function of $q\sigma_\text{eff}/\pi$, the first peak as a function of $q\sigma_\text{eff}/\pi$, the 2nd-5th peak as a function of $q\sigma_\text{eff}/\pi$, and the 5th-8th peak as a function of $q\sigma_\text{eff}/\pi$.
Again, we see that scaling with $\sigma_\text{eff}$ works really well. The most significant differences between the "practically hard" systems (i.e. HS, WCA, and $1/\kappa\sigma=0.01$) and "slightly softer" systems (i.e. $1/\kappa\sigma>0.01$) can be seen in the lowering of higher order peaks both at $\eta_F$ and at $\eta_M$. Notice, however, that for the first peak, there are only sigficant differences at $\eta_M$ and not at $\eta_F$.

| at $\eta_F$ | at $\eta_M$ |
| -------- | -------- |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |

Now, lets compare to the hard-core Yukawa systems with a WCA cutoff. Again, we see that using the WCA cutoff brings the $S(q)$ closer to that of HS or WCA.

| at $\eta_F$ | at $\eta_M$ |
| -------- | -------- |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |

To further show that the effective hard-spheres diameters obtained from the freezing packing fractions dead on map the structure, I also determined $\sigma_\text{eff}$ by mapping the $g(r)$ and $S(q)$ on top of the one of hard spheres via a mean-squared error reduction. The table underneath shows the resulting $\sigma_\text{eff}$ for the three different methods.


----

### Nucleation of the BCC crystal
Next, I looked at the nucleation of the BCC phase of the hard-core Yukawa particles. Picking the state points at which to study nucleation is quite tricky, as the BCC phase is enclosed by the fluid and FCC phase in the phase diagrams. Consequently, for many combinations of $\beta\epsilon$ and $1/\kappa\sigma$, the BCC phase transitions into the FCC phase before a high enough supersaturation can be reached.
In the end, I picked different combinations of $\beta\epsilon$ and $1/\kappa\sigma$ and calculated the nucleation barriers for three different values of the supersaturation. Underneath, the resulting nucleation barriers and a table containing all the information of the barriers.


Interestingly, the nucleation barriers seem to fully dependent on the supersaturation and not on $\beta\epsilon$ or $1/\kappa\sigma$. This is different to what we saw for the nucleation barriers of the FCC phase, for which the nucleation barrier (at fixed $\beta|\Delta\mu|$) decreased with increasing $1/\kappa\sigma$.
Note that for equal supersaturation, the value of $\beta\gamma/\rho_s^{2/3}$ is the approximately same as well. This is no surprise, as these barriers also have similar heights and the height is given by $\Delta G^* =\frac{16\pi\gamma^3}{3|\Delta\mu|^2\rho_s^2}$ .

#### Nucleation barriers for more intermediate values of the supersaturation
I also computed the nucleation barriers for more intermediate values of the supersaturation. I only did this for $\beta\epsilon=81$ and $1/\kappa\sigma=0.40$, as the above results show that for the nucleation barriers of the BCC phase are fully dependent on the supersaturation and not on $\beta\epsilon$ or $1/\kappa\sigma$.

Underneath are the resulting nucleation barriers, the barrier height as a function of the supersaturation, and a table with all information.




#### Nucleation barriers in $NVT$ instaed of $NPT$ esemble
As the packing fraction of the metastable fluid and BCC crystal differ only in the third (or even fourth) decimal, I wanted to check whether performing the simulations in the $NVT$ ensemble instead of the $NPT$ ensemble made any difference. Underneath, the resulting nucleation barriers for two state points:
- $\beta\epsilon=20$, $1/\kappa\sigma=0.30$, and $\beta|\Delta\mu|=0.30$
- $\beta\epsilon=81$, $1/\kappa\sigma=0.40$, and $\beta|\Delta\mu|=0.35$

The barriers of the $NVT$ and $NPT$ ensemble do not significantly differ from each other within the error of the (umbrella sampling) method.

*Note: the NPT barriers are averaged (over two independent runs for each umbrella window), whereas the NVT barriers are non-averaged.*


#### Nucleation barriers using different order parameters
Lastly, I wanted to check if the height of the nucleation barrier of the BCC phase is actually independent of the choice of order parameter (as already shown for the FCC phase). For all simulations in this chapter I used the dot-product of $q_6$ with cutoff value $d_c=0.7$ to determine if a connection is solid-like and then define a particles as solid when it has at least $\xi_c=6$ solid-like connection. I thus want to vary the order parameter by varying $d_c$ and $\xi_c$.

To get a feeling for $d_6(i,j)$ and $\xi(i)$ in the metastable fluid and BCC phase, underneath I give the histograms for a fluid and BCC crystal both with $\beta\epsilon=81$ and $1/\kappa\sigma=0.40$ and at $\eta=0.1342$. For the histogram of $d_6(i,j)$, I zoomed in on the interesting region. $\xi(i)$ is determined using $d_c=0.7$.

 

I computed the nucleation barriers for the values $d_c=[0.675,0.700,0.725,0.750]$ (using $\xi_c=6$) and $\xi_c=[5,6,7,8,9]$ (using $d_c=0.700$) for two state points:
- $\beta\epsilon=20$, $1/\kappa\sigma=0.30$, and $\beta|\Delta\mu|=0.30$
- $\beta\epsilon=81$, $1/\kappa\sigma=0.40$, and $\beta|\Delta\mu|=0.35$

The resulting nucleation barriers are shown underneath. We see that, within the error of the method, the height of the barrier is independent of the specific choice of order parameter.
*Note that these barriers, except the one of $d_c=0.700$ and $\xi_c=6$, are not averaged; hence, they are not the prettiest barriers on this page.*

| | As a function of $d_c$ <br> (with $\xi_c=6$) | As a function of $\xi_c$ <br> (with $d_c=0.700$) |
| -------- | -------- | -------- |
| $\beta\epsilon=81$ <br> $1/\kappa\sigma=0.40$ <br> $\beta$ l$\Delta\mu$l$=0.35$ |  |  |
| $\beta\epsilon=20$ <br> $1/\kappa\sigma=0.30$ <br> $\beta$ l$\Delta\mu$l$=0.30$ |  |  |

With these barriers for different order parameters, we can also show how powerful the fitting equation is. For example, underneath is a table with some information on the fits of the system with $\beta\epsilon=20$, $1/\kappa\sigma=0.30$, and $\beta|\Delta\mu|=0.300$ shown above. Here $\alpha$ is the difference between the measured radius $R(d_c,\xi_c)$ and the radius according to CNT $R_{CNT}$. We see that, as the nucleus size $n^*$ decreases, $\alpha$ decreases as well. As a result all fits find $R_{CNT}/\sigma\approx3.0$. Additionally, this fitting equation makes sure the fitted interfacial free-energy surface density $\beta\gamma\sigma^2$ does not depend on the order parameter. For all barriers, we find $\beta\gamma\sigma^2\approx0.33$.
*Note that the barriers that we fitted are not averaged; hence, they are not the greatest barriers. For averaged barriers, $R_{CNT}/\sigma$ and $\beta\gamma\sigma^2$ would probably agree even more for the different order paramters.*


##### Completely different order parameter

Additionally, I also wanted to check if using a completely different order parameter would influence the height of the barrier. Thus, instead of using the $q_6$ dot-product, I wanted to see what using the $q_8$ dot-product would do.

First, to get a feeling $d_8(i,j)$ in the metastable fluid and BCC phase, underneath I give the histograms for the same configurations as used for the $d_6(i,j)$ histogram shown above.
Again, the histrogram is zoomed in on the interesting region. $\xi(i)$ is determined using $d_c=0.7$.


Thus, looking at this histogram, we see that a cutoff somewhere between 0.5 and 0.6 would probably work. Testing on a configuration with a nucleus it in, I found that $d_c=0.55$ (and $\xi_c=6$) gives approximately the same classification as the original $q_6$ dot-product with $d_c=0.7$ and $\xi_c=6$. Thus, I used $d_c=0.55$ to compute the histrogram for $\xi$, which is shown below.


From these histograms I picked three values for $d_c$ (with $\xi_c=6$) to compute the nucleation barriers with: $d_c=[0.525,0.550,0.575]$. And four values for $\xi_c$ (with $d_c=0.55$): $\xi_c=[5,6,7,8]$. I computed the nucleation barriers for the state point:
- $\beta\epsilon=20$, $1/\kappa\sigma=0.30$, and $\beta|\Delta\mu|=0.30$

The resulting nucleation barriers together with a table containing the fitting parameters are shown underneath. Again, we that changing the order parameter (wihtin reasonable bounds) does not influence the height of the nucleation barrier.
*Note that these barriers are not averaged; hence, they are not the prettiest barriers on this page.*

| As a function of $d_c$ <br> (with $\xi_c=6$) | As a function of $\xi_c$ <br> (with $d_c=0.550$) |
| -------- | -------- |
|  |  |


---
---

## Structure analysis

To study the structure, I make use of bond-orientational order parameters (BOPs). I tried using the local BOPs, different variants of averaged BOPs, and combinations of both local and averaged. For now I found that just using the averaged BOPs introduced by Lechner and Dellago works best.
$$ q_{lm}(i) = \frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b(i)} Y_l^m(\mathbf{r}_{ij}),
$$
where $\mathcal{N}_b(i)$ is the set of the $N_b(i)$ nearest neighbors of particle $i$ determined using the SANN algorithm.
$$ \bar{q}_{lm}(i) = \frac{1}{N_b(i)+1} \sum_{k\in\{i,\mathcal{N}_b(i)\}} q_{lm}(k)
$$
$$ \bar{q}_l(i) = \sqrt{\frac{4\pi}{2l+1} \sum_{m=-l}^l |\bar{q}_{lm}(i)|^2}
$$

### Comparing PCA and UML
For all the particles in the system I calculate first eight BOPs ($\bar{q}_l$ with $l\in[1,8]$). To reduce the dimensionality, I tried using both principal component analysis (PCA) and unsupervised machine learning (UML). For the UML I used the network described by Emanuele in his JCP **151** (2019) article. I find that PCA and UML do extremely similar things.
For example, underneath you can see the PCA (top row) and UML (bottom row) decomposition performed on a handful of snapshots containing a metastable fluid of hard spheres ($\beta P\sigma^3=17.0$). To highlight the similarites between the two methods, I colored the data points using $\bar{q}_6$ (first column), $\bar{q}_7$ (second column), and the probability $P_{GMM}$ to belong to one of two Gaussian distributions (third column).

 |  |  |
--- | --- | --- |
 |  |  |

As the methods are this similar, I decided to focus on PCA, since PCA is "simpler" and provides valuable information on the importance of each BOP.

### Brute force nucleation events
To study if there are any precursors for nucleation in the metastable fluid, I performed MC simulations of the NVT-ensemble and waited for them to spontaneously nucleate. I did this for hard spheres, nearly-hard hard-core Yukawa particles (crystalizing into FCC), and soft hard-core Yukawa particles (crystalizing into BCC). For all systems, I performed the simulations at a packing fraction for which the nucleation barrier had a height of approximately 15-16$k_BT$, e.g. I used $\eta^*=0.5385$ for hard spheres (corresponding to $\beta P\sigma^3=17.5$).

Next, I "trained" the PCA models for each system using just the snapshots before nucleation, so configurations containing mostly fluid particles (and maybe some small solid nuclei).

Furthermore, I determined local packing fractions using the volume of the Voronoi cells obtained from voro++.

#### Similarities of PCA for the different systems
The first thing I noticed is that the PCA is very similar for the different systems. For example, underneath you can see the contribution of each BOP to the frist and second principal component, as well as the average value of each BOP in the metastable fluid. The "HC Yukawa fluid (FCC)" depicted in these figures has $\beta\epsilon=81, 1/\kappa\sigma=0.01$ and the "HC Yukawa fluid (BCC)" has $\beta\epsilon=81, 1/\kappa\sigma=0.40$.

It seems that the fluid does not "know" if it will nucleate into FCC or BCC.

|  |  |
| -------- | -------- |
| Contribution to PC-1 | Contribution to PC-2 |

|  |
| -------- |
| Mean of each BOP |

#### Hard spheres
To perform the brute force nucleation, I started 10 MC simulations of the NVT ensemble from the same starting configuration. All simulations were performed at $\eta^*=0.5385$ (corresponding to $\beta P\sigma^3=17.5$) and the nucleation was tracked using $r_c=1.40\sigma$, $d_c=0.7$, and $\xi_c=6$. The critical nucleus for this state point contains $n^*=75$ particles. The maximum distance a particle could be displaced in any direction during a single MC step was $\Delta r=0.05\sigma$, resulting in a long-time diffusion time of $\tau=\sigma^2/6D_l = 44.9(1)\cdot10^3$ MC cycles (in the metastable fluid).

To study the local environment of the region in which the nucleus is born, I pick a coordinate that captures the center of the nucleus at the start of nucleation well, and then for each snapshot select all the particles inside a sphere of radius $1.5r_c$ around this coordinate. On average 40 particles are selected. The local properties (such as packing fraction, BOPs, and PC-1) of the nucleation region are then calculated by averaging their values over the selected particles.
*Note that I just picked the radius to be $1.5r_c$ such that it would automatically scale with the system, and 40 particles seemed like a reasonable number to me.*

To illustrate this, underneath you can see four snapshots during a nucleation event *(note: event 3)*. The particles classified as fluid are depicted at a quarter of their actual size, and the particle inside the selected region are colored red. The time $t_0$ is an approximate guess for the start of the nucleation event.

|  |  |
| -------- | -------- |
| Before nucleation <br> ( $(t-t_0)/\tau = -0.13$ ) | Nuclues of $n=32$ <br> ( $(t-t_0)/\tau = 0.53$ ) |
|  |  |
| Nuclues of $n=95$ <br> ( $(t-t_0)/\tau = 1.07$ ) | Nuclues of $n=142$ <br> ( $(t-t_0)/\tau = 1.43$ ) |

Thus, by focussing only on the particles inside the selected region, I was able to track the local properties as a function of the "time" for each nucleation event. Underneath, you can see the average local packing fraction, BOPs (and $w_l$), and first and second principal components for two events. The other events show similar behavior.

Some important things to notice:
* The packing fraction and structural ordering (i.e. $\bar{q}_6$, $\bar{q}_8$, PC-1, and PC-2) increase simultaneously.
* There is no increase in packing fraction or structure before nucleation starts.
* Once nucleation starts, there is a very steep increase in the first principal component (PC-1).
* Nucleation happens vary rapidly, i.e. inside 1 unit of long-time diffusion time the nucleus has grown from 0 to past its critical size.
* The first principal component (PC-1) is a great order parameter for nucleation.
* Increase in $\bar{q}_4$, combined with a negative $\bar{w}_6$ and $\bar{w}_4$ suggests that indeed an FCC crystal is formed.

| Event 3 | Event 7 |
| -------- | -------- |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |

To better see that the local packing fraction and first principal component increase simultaneously, underneath you see $\eta_{local}$ and PC-1 plotted together. The range of the y-axes is choosen such that PC-1=0.0 (which is the mean value in the fluid) and $\eta_{local}=\eta^*$ coincide.
Additionally, the last row shows the Pearson correlation between the inverse of the Voronoi volume and PC-1 for the particles in the selected region as a function of time. We see that in the fluid the two are uncorrelated.

| Event 3 | Event 7 |
| -------- | -------- |
|  |  |
|  |  |

To further show that there is no structural precursor for nucleation present in the fluid, I computed the autocorrelation functions (ACF) for the first and second principal component (using only configurations before nucleation starts). Underneath is the resulting ACF. Note that the ACF is normalized such that at time 0 it is equal to 1.
We clearly see that the correlation quickly decays. Much too fast for there to be any talk of a structural precursor.

|  |  |
| -------- | -------- |
| Autocorrelation function of PC-1 and PC-2 | Same with logarithmic y-axis |

Lastly, I also tried to find a precursor by looking at all of the TCC clusters and comparing them with the first principal component. Underneath are some examples of TCC clusters in the selected region for event 3, the results of other clusters are similar. The second column shows the population, which is the number of particles that are part of at least one cluster divided by the number of particles in the region. The third column shows the average number of clusters that a particle is involved in. The gray, dashed horizontal lines indicate the average value of the TCC clusters in the metastable fluid. The black, dashed horizontal lines indicate Pc-1=0.0 (which is the average value of PC-1 in the fluid).
*Note that the range of the y-axis for the number of clusters per particle is chosen such that the average value in the fluid (gray dashed line) coincides with PC-1=0.0 (black dashed line). There are, however, three exceptions for which this was not possible: 7A, 7K, and 8B.*

Again, we see that there is no perceivable change before nucleation starts. (Also for the clusters not shown here). For many clusters, the increase/decrease of the population and number of clusters per particle starts the moment nucleation starts, meaning that they increase/decrease simultaneously with PC-1. The average number of 6A clusters per particle is a great example of this. Some larger clusters, however, have a small delay of ~0.25$\tau$. You can best see this delay in the number of clusters per particle of the BCC_9, 11F, and FCC clusters.

Some other specific notes:
* The population and number of clusters per particle increase for cluster types that are also found in a bulk FCC phase. These are: 6A (sp4c), BCC_9, and FCC.
* There is also an increase in the population and number of clusters per particle for cluster types that are not really found in a bulk FCC phase. These are: 8A, 9K, 11F, 12E, HCP.
However, when the nucleation is tracked for longer, these clusters will at some point start to decrease again. Note that all these cluster types are found in both bulk HCP and bulk BCC phases, except 8A which is only found in bulk BCC.
* Note that 8K is also a cluster that is only found in bulk BCC, but during nucleation this cluster remains at the same low values as in the fluid.
* Then there are the cluster types that are not found in a bulk FCC phase that show a decrease in population and number of clusters per particle. These are: 7A (sp5c), 7K, 8B, 9B, 10B, 11C, 11E, 12B, 12D.

| **6A (sp4c)** |  |  |
| -------- | -------- | -------- |
| **7A (sp5c)** |  |  |
| **7K** |  |  |
| **8A** |  |  |
| **8B** |  |  |
| **BCC_9**|  |  |
| **11F** |  |  |
| **FCC** |  | ) |
| **HCP** |  |  |

To better show that there is no precursor in the fluid that gets "averaged out" when I take the average over the ~40 particles in the studied region, I also visualized the region and colored the particles according to their local properties. For example, underneath you can see the studied region of the event above from $(t-t_0)/\tau=-0.38$ to $0.69$.
*Note: for PC-1 the color indicates the deviation from the mean (0.0) with standard deviation $\sigma=0.05$.*

| $\eta_{local}$ | 6A | FCC |
| -------- | -------- | -------- |
|  |  |  |


|  |  |  |
| -------- | -------- | -------- |
| **PC-1** | **8A** | **11F** |

Additionally, we could also just plot the local properties of a handful of these particles over time. Underneath you can see PC-1 and the number of 6A, 8A, and FCC clusters for 10 particles inside the nucleus. Each of the differently colored thin lines corresponds to one of the 10 particles, and the thick dark blue line indicates the mean of these 10 particles.
Again, we see that there is no change in the fluid before nucleation starts and that everything starts to change as soon as nucleation starts.

|  |  |
| -------- | -------- |
| PC-1| 6A |
|  |  |
| 8A | FCC |

#### Hard-core Yukawa (FCC)

Next, I studied multiple nucleation events of a handful of nearly-hard hard-core Yukawa systems. Similar to hard spheres, I performed the brute force nucleation by starting 10 MC simulations of the NVT ensemble from the same starting configuration. The packing fraction ($\eta^*$) used for each hard-core Yukawa system are given in the table underneath together with the corresponding pressure, critical nucleus size, and barrier height. For comparison, the top row gives the data for the hard-sphere system discussed above.

To track the nucleation events, I used $d_c=0.7$ and $\xi_c=6$ (same as for hard spheres above), and the radial cutoff $r_c$ for each system is given in the table as well. Lastly, the table provides the maximum displacement of a particle in any direction during a single MC step $\Delta r/\sigma$ and the resulting long-time diffusion time $\tau=\sigma^2/6D_l$ (in the metastable fluid).


Again, I track the local properties as a function of the "time" for each nucleation event by focussing only on the particles inside the selected region. Underneath, you can see the average local packing fraction, BOPs (and $w_l$), and first and second principal components for an event of $\beta\epsilon=81$ and $1/\kappa\sigma=0.01$ (left) and an event of $\beta\epsilon=8$ and $1/\kappa\sigma=0.04$ (right). The last row gives the autocorrelation function for the first and second principal component. The other events show similar behavior.

As for hard spheres, I find that the packing fraction and structural ordering increase simultaneously when nucleation starts, that there is no precursor for nucleation in the metastable fluid, and that nucleation happens extremely fast.

| $\beta\epsilon=81$ and $1/\kappa\sigma=0.01$ (event 6) | $\beta\epsilon=8$ and $1/\kappa\sigma=0.04$ (event 9) |
| -------- | -------- |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
| |  |
|  |  |

I also did the full TCC analysis for the nucleation event on the left. The results are very similar to the one in the hard spheres section. So also for the hard-core Yukawa system the decrease/increase of the population and number of clusters per particle starts when nucleation starts, meaning that they increase/decrease simultaneously with PC-1. There is perceivable change in the TCC clusters before nucleation starts.

One thing I didn't get to see in the TCC analysis of the hard spheres event, was the decrease of clusters of type 8A, 9K, 11F, 12E, and HCP. These clusters first increase when nucleation starts, but are expected to decrease again when the nucleus has crossed a certain size. I didn't see this for the hard spheres event, as the simulations are stopped once the nucleus reaches a size of approximately 200 particles. However, for this nucleation event of hard-core Yukawa particles I did observe the expected decrease in these clusters. See the figures underneath of $\beta\epsilon=81$ and $1/\kappa\sigma=0.01$ (event 6).
*Note that the range of the y-axis for the number of clusters per particle is chosen such that the average value in the fluid (gray dashed line) coincides with PC-1=0.0 (black dashed line).*

| **6A** |  |  |
| -------- | -------- | -------- |
| **8A** |  |  |
| **9K** |  |  |
| **11F** |  |  |
| **12E** |  |  |
| **FCC** |  |  |
| **HCP** |  |  |
| **8K** |  | |

I also visualized the region and colored the particles according to their local properties for the nucleation event above. Underneath you can see the event from $(t-t_0)/\tau=-0.29$ to $0.67$. Again we see that there is no precursor in the fluid that gets "averaged out" when I take the average over the ~40 particles in the studied region.
*Note: for PC-1 the color indicates the deviation from the mean (0.0) with standard deviation $\sigma=0.05$.*

| $\eta_{local}$ | 6A | FCC |
| -------- | -------- | -------- |
|  |  |  |


|  |  |  |
| -------- | -------- | -------- |
| **PC-1** | **8A** | **11F** |

#### Hard-core Yukawa (BCC)

Lastly, I studied multiple nucleation events of three soft hard-core Yukawa systems that form BCC crystals. Again, I performed the brute force nucleation by starting 10 MC simulations of the NVT ensemble from the same starting configuration. The table underneath gives all the necessary information on the systems.


As before, I track the local properties as a function of the "time" for each nucleation event by focussing only on the particles inside the selected region. Underneath, you can see the average local BOPs (and $w_l$), and first and second principal components for an event of $\beta\epsilon=81$ and $1/\kappa\sigma=0.40$ (left) and an event of $\beta\epsilon=39$ and $1/\kappa\sigma=0.35$ (right). The last row gives the autocorrelation function for the first and second principal component. The other events show similar behavior.
*Note: there is no figure tracking the packing fraction as the difference in packing fraction between the fluid and solid phase is much smaller than the fluctuations in the local packing fraction.*

Similar to the nucleation of the FCC crystal, I find that generally nucleation happens extremely fast and that the first principal component (PC-1) rapidly increases when nucleation starts. Thus, PC-1 is a great order parameter for nucleation of the BCC phase as well. Furthermore, almost no increase in $\bar{q}_4$, combined with a positive $\bar{w}_6$ suggests that indeed a BCC crystal is formed.

| $\beta\epsilon=81$ and $1/\kappa\sigma=0.40$ (event 5) | $\beta\epsilon=39$ and $1/\kappa\sigma=0.35$ (event 3) |
| -------- | -------- |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |

I also did a full TCC analysis on the nucleation event above on the left to see if TCC picked up on any precursors. Underneath are some examples TCC clusters in the selected region of this event.
As for hard spheres and nearly-hard hard-core Yukawa particles, all clusters show no perceivable change in the population and average number of cluster per particle before nucleation starts. The decrease/increase of the population and number of clusters per particle starts when nucleation starts, meaning that they increase/decrease simultaneously with PC-1.
*Note that the range of the y-axis for the number of clusters per particle is chosen such that the average value in the fluid (gray dashed line) coincides with PC-1=0.0 (black dashed line). There are, however, three exceptions for which this was not possible: 7A, 7K, and 8B.*

Many cluster types behave the same as for the nucleation of the FCC crystal. These are: 6A, 7A, 8B, 9B, 10B, 11C, 11E, 12B, 12D. Some notes on the clusters that behave differently (or slightly differntly):
* 8A, 9K, 11F, and 12E, which increased and then decreased for the FCC nucleation, only increase for the BCC nucleation. This is to be expected, since these can all be found in a bulk BCC phase.
* HCP and FCC clusters are sometimes found in the BCC nucleus.
* 7K remains at a population of ~1.0 (for FCC nucleation it decreases), but the number of clusters per particles slightly decreases.
* The number of 8K clusters increases during nucleation (for FCC nucleation it remains at a low population).
* The population of BCC_9 behaves the same as for FCC nucleation, but the number of clusters per particle increases to a lower average, i.e. ~6-7 instead of ~12-15.

| **6A (sp4c)** |  |  |
| -------- | -------- | -------- |
| **7A (sp5c)** |  |  |
| **7K** |  |  |
| **8A** |  |  |
| **8B** |  |  |
| **8K** |  |  |
| **9K** |  |  |
| **BCC_9**|  |  |
| **11F** |  |  |
| **12E** |  |  |
| **FCC** |  |  |
| **HCP** |  |  |

#### Failed nucleation

With these simulations I also have quite some data on "failed" nucleation event, i.e. nuclei that melted into fluid again. I have not extensively studied these failed nuclei, but there are some things that I have noticed by just briefly looking at a small collection of them.
* There are quite some nuclei that just do not managed to cross the nucleation barrier. These nuclei look the same as nuclei that do manage to cause crystal nucleation. With this I mean that they have approximately (nice) spherical shapes and have "normal" values for the local properties (packing fraction, BOPs, PCA, and TCC). So I guess that they really just melt again because they did not manage to cross the nucleation barrier.
* Some failed nuclei reach sizes greater than the critical nucleus size before melting again; however, in all the case that I have observed this, the nucleus reaches this high number of particles due to strings/arms of solid-like particle that are attached to a central nucleus, which is often smaller than the critical nucleus. In a couple of cases the nucleus is just all "stringy" and ugly. In other some other cases, two (or more) small nuclei try to merge or connect to each other, which the simulation then counts as one big nucleus. The "big nucleus" then melts, when the smaller nuclei disconnect and separately melt back into fluid. Underneath are two examples of such nuclei.

 
**Left:** Example of two small, but connected nuclei, that together are counted as one big nucleus $>n^*$. In this case $n=117$. *(HS: run 2)*
**Right:** Example of a small nucleus with "arms" of solid-like partices attached to it. In this case $n=170$. *(HC Yukawa $\beta\epsilon=8$, $1/\kappa\sigma=0.01$: run 6)*

### Correlation with TCC

Not only can we track the average local properties of a region during a nucleation event, it can also be useful to just look at the snapshots colored with the local properties. For example, underneath you can see a snapshot *(HS: event 3, $(t-t_0)/\tau = 1.01$ )* where each particle is colored with its value for

* **Top left**: the first principal component. Here blue/red indicates two standard deviations below/above the mean value. *(Note: mean=0.0 and $\sigma=0.05$).*
* **Top right**: the local packing fraction determined using the volume of the Voronoi cell. For this system $\eta^*=0.5385$ and $\eta_s=0.598$.
* **Bottem left**: the number of 11F clusters a particle is involved in. Gray particles are not part of a 11F cluster.
* **Bottom right**: the number of 9B clusters a particle is involved in. Gray particles are not part of a 9B cluster.

|  |  |
| -------- | -------- |
| PC-1:  | $\eta_\text{local}$:  |
|  |  |
| 11F:  | 9B:  |

Comparing the top left figure (PC-1) with the two bottom figure, we observe some correlation between the regions with high/low PC-1 and the regions with 11F and 9B particles. In particular, we see that regions with high PC-1 contain many 11F clusters but few or even none 9B clusters, whereas regions with low PC-1 contain many 9B clusters but few or even none 11F clusters.

Just by eye, we can observe similar correlations for many of the TCC clusters. To determine the actual correlation between the first principal component and the number of TCC clusters, I calculated the Spearman correlation between the two.

Underneath you can see the Spearman correlation between the first principal component (PC-1) and the number of TCC clusters of a certain type a particle is involved in ($n_{TCC}$) for a metastable fluid of hard spheres, nearly-hard hard-core Yukawa particles (forms FCC), and soft hard-core Yukawa particles (forms BCC). The correlations were calculated using 10 configurations (all before nucleation, so containing mostly fluid particles and maybe some small solid nuclei), and the errorbars using the difference between the first 5 and second 5 configurations. The blue bars indicate clusters that consist of one or more square pyramidal subclusters, while the green bars indicate clusters that consist of one or moretetrahedral subclusters. Gray clusters contain neither, or both.

*Note that I am not entirely sure yet if TCC determines the small clusters sp3b, sp3c, and sp4b correctly. So that could explain why sp3b is the only green bar with possitive correlation and sp4b the only blue bar with negative correlation.*

The hard spheres and nearly-hard hard-core Yukawa particles have very similar correlations. Generally speaking, we see that the blue bars (square pyramidal) have a positive correlation, whereas the green bars (tetrahedral) have a negative correlation.
Futhermore, the correlations for the soft hard-core Yukawa particles show some differences with the correlations for hard spheres; however, the overall trend is similar.

|  |  |
| -------- | -------- |
| Hard spheres | Hard-core Yukawa ($\beta\epsilon=81$, $1/\kappa\sigma=0.01$) (->FCC) |
| |  |
| Hard-core Yukawa ($\beta\epsilon=8$, $1/\kappa\sigma=0.04$) (->FCC) | Hard-core Yukawa ($\beta\epsilon=81$, $1/\kappa\sigma=0.40$) (->BCC) |

Underneath is also the average population of TCC clusters in the different metastable fluids. The populations are calculated using the same 10 configurations as for the correlation, and the errorbars using the difference between the first 5 and second 5 configurations.

To make it complete, underneath is a comparision between the average number of clusters per particle in the different metastable fluids.
*Note: the bars of the 7T_s, 7T_z, and 6Z clusters lies outside the range of this plot.*


---

## Comparing metastable fluids: FCC vs BCC

We wanted to get a better idea of if there are differences in the structural fluctuations in the metastable fluid of a system that forms the FCC phase and a system that forms the BCC phase. To this end, I ran MC simulations in the NVT ensemble of a hard-core Yukawa fluid with $\beta\epsilon=81$ and $1/\kappa\sigma=0.01$, which I will refer to as the "FCC fluid", and a system with $\beta\epsilon=81$ and $1/\kappa\sigma=0.40$, which I will refer to as the "BCC fluid". For both systems I performed the simulations at a range of different packing fractions (and thus at a range of different super saturations). Underneath are the packing fractions ($\eta^*$) for which I performed the simulations together with the corresponding pressure and supersaturation.

| FCC fluid <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.01$| BCC fluid <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.40$ |
| -------- | -------- |
|  |  |

*Note: although the supersaturations differ with almost a factor of two, the last row of both tables corresponds to a nuclation barrier with a height of approximately $16k_BT$. (More specific, $16.3k_BT$ for the FCC fluid and $15.3k_BT$ for the BCC fluid)*

For each packing fraction I ran 5 independent (short) simulations, and trained a PCA model on the BOPs (Lechner & Dellago) of the final configuration of each simulation. Note that I trained a PCA model on the combined data of the FCC and BCC fluid. I did this because I wanted to compare fluctuations of the first principal component (PC-1) of both fluids, and using separate PCA models for the FCC and BCC fluids would result in (slight) differences in the PCA decomposition, which would influence this comparison.

Underneath you can see the scatterplot of PC-1 against PC-2. The red arrows indicate how each BOP contributes to the PCs, and the histogram along PC-1 and PC-2 are given in gray.
64.0% of info is explained by PC-1 and 9.96% by PC-2.


### Averages and standard deviations
I first started out with just looking at the averages and standard deviations of the different order parameters in the different fluids. I calculated these using the last 10 snapshots of each simulation, such that I had 50 snapshots for each packing fraction, with $N=10976$ and $N=11664$ for the FCC and BCC fluid respectively. Underneath you can see the average BOPs (first row), standard deviation of the BOPs (second row), average PC-1 (third row), and standard deviation of PC-1 (fourth row) as a function of the packing fraction for both the FCC and BCC fluid.

We see that the average and standard deviation of the BOPs slightly increase with increasing packing fraction. The figures of PC-1 are a bit noisier as PC-1 is more sensitive to small solid-like clusters in the fluid, but we still see the same general increase with increasing packing fraction.

| FCC fluid <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.01$| BCC fluid <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.40$ |
| -------- | -------- |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |

### PDF of $\bar{w}_6$ and $\bar{w}_4$ for high $\bar{q}_6$
This is something that John Russo uses in some of his nucleation papers to stress the point that it is clear from the metastable fluid that FCC is formed and not BCC or HCP. For example, for the nucleation of hard spheres **[[Russo & Tanaka, Sci Rep 2 (2012)]](https://doi.org/10.1038/srep00505)**, or for the nucleation of the Gaussian core model **[[Russo & Tanaka, Soft Matter 8 (2012)]](https://doi.org/10.1039/C2SM07007C)**.

What he does is plot the probability distribut function (PDF) of $\bar{w}_6$ and $\bar{w}_4$ of the metastable fluid for particles with $\bar{q}_6>\bar{q}_{6,t}$, where $\bar{q}_{6,t}$ is some threshold value. The idea behind this is that particles with $\bar{w}_6>0$ are more BCC-like, whereas particles with $\bar{w}_6<0$ are FCC/HCP-like. Next, $\bar{w}_4$ can be use to distinguish between FCC and HCP. Particles with $\bar{w}_4>0$ are more HCP-like, whereas particles with $\bar{w}_4<0$ are FCC-like.

Underneath you can see the resulting PDFs for the different systems studied for their spontaneous nucleation. We use a range of different thresholds: $\bar{q}_{6,t}=[0.25,0.26,0.27,0.28,0.29,0.30,0.31,0.32]$. The darkness of the lines indicates the increase in $\bar{q}_{6,t}$, as well as the black arrow.

We clearly see for the systems that should nucleate the FCC phase, i.e. top three rows of figures, that the metastable fluid tents more towards negative $\bar{w}_6$ (and thus FCC/HCP) for increasing $\bar{q}_{6,t}$ and more toward negative $\bar{w}_4$ (and thus FCC). However, when we look at the system that should nucleat the BCC phase, then we see something similar. Although the peak of the PDF of $\bar{w}_6$ lies at positive $\bar{w}_6$ for low $\bar{q}_{6,t}$, the peak shifts to negative $\bar{w}_6$ for higher $\bar{q}_{6,t}$. This would then indicate that the particles with high $\bar{6}$ tent more towards the FCC/HCP-like instead of the BCC-like...

| **HS** $\eta^*=0.5385$ | **HS** $\eta^*=0.5385$ |
| ------ | ------ |
|  |  |
| **Yukawa** $(\beta\epsilon=81,1/\kappa\sigma=0.01)$ </br> $\eta^*=0.4681$ | **Yukawa** $(\beta\epsilon=81,1/\kappa\sigma=0.01)$ </br> $\eta^*=0.4681$ |
|  |  |
| **Yukawa** $(\beta\epsilon=8,1/\kappa\sigma=0.04)$ </br> $\eta^*=0.4400$ | **Yukawa** $(\beta\epsilon=8,1/\kappa\sigma=0.04)$ </br> $\eta^*=0.4400$ |
|  |  |
| **Yukawa** $(\beta\epsilon=81,1/\kappa\sigma=0.40)$ </br> $\eta^*=0.1305$ | **Yukawa** $(\beta\epsilon=81,1/\kappa\sigma=0.40)$ </br> $\eta^*=0.1305$ |
|  |  |

For comparison, underneath you can see, for the different systems, the fraction of particles that has a negative $\bar{w}_6$ as a function of the $\bar{q}_6$ threshold. For the fluids that should nucleate the FCC phase, we see that the fraction of negative $\bar{w}_6$ starts out just bigger than 0.5, and in all cases only starts to increase from $\bar{q}_{6,t}\gtrsim0.15$. The fluid that should nucleate the BCC phase on the other hand, starts out with a fraction of negative $\bar{w}_6$ smaller than 0.5. Notice that it dips around $\bar{q}_{6,t}\approx0.20$, before increasing to values larger than 0.5.

| **HS** $\eta^*=0.5385$ | **Yukawa** $(\beta\epsilon=81,1/\kappa\sigma=0.01)$ </br> $\eta^*=0.4681$ |
| -------- | -------- |
|  |  |
| **Yukawa** $(\beta\epsilon=8,1/\kappa\sigma=0.04)$ </br> $\eta^*=0.4400$ | **Yukawa** $(\beta\epsilon=81,1/\kappa\sigma=0.40)$ </br> $\eta^*=0.1305$ |
|  |  |

Could this have something to do with the fact that the BCC phase has lower $\bar{q}_6$ than the FCC phase?


### Radial correlation

Next, I calculated the radial correlation of PC-1 to get an idea of the lenght scale of the structural fluctuations. For calculating the radical correlation I used

$$
S(\Delta r) = \frac{\langle x_i x_j \rangle_{\Delta r}}{\langle x_i^2 \rangle},
$$

where $x$ is the variable for which you want to calculate the correlalation (PC-1 in this case), $\langle x_i x_j \rangle_{\Delta r}$ is the average over all particle pairs $i,j$ seperated by $\Delta r$, and $\langle x_i^2 \rangle$ is the average over all particles and serves as a normalization factor.

| FCC fluid <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.01$| BCC fluid <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.40$ |
| -------- | -------- |
|  |  |
|  |  |

To better compare the results of the FCC and the BCC fluid, underneath I plot the results for the highest and lowest packing fraction studied and scale the x-axis with the density.

 

Integrating the first peak of the $g(r)$ shows that there are approximately 12.5 neighbors in the first shell for the FCC fluid, and 13.5 for the BCC fluid.

### First part of nucleation barrier

First part of the nucleation barrier for hard spheres ($\eta=0.538$), hard-core Yukawa particles with $\beta\epsilon=81$ and $1/\kappa\sigma=0.01$ ($\eta=0.468$), and hard-core Yukawa particles with $\beta\epsilon=81$ and $1/\kappa\sigma=0.40$ ($\eta=0.134$).

---

## Comparing simulation methods

### Cluster size distribution in metastable fluid

To make sure that the simulation method, "normal" MC (move acceptance of $\sim30\%$) vs. kinetic MC (move acceptance of $\sim80-90\%$), does not influence the metastable fluid, I looked at the distribution of nucleus sizes for a handful of systems. For this I simply used the data of the brute force nucleation events, as I wrote away the sizes of all clusters in the system frequently.

For each system, I tally the cluster sizes of each nucleation event up to $t=t_0$, where $t_0$ is the start of nucleation as used in the analysis of the nucleation event. Joining the data of the different simulations, I then calculate the size distribution $P(n)$ by dividing by the total number of snapshots and the number of particles in the system (same as you would do for the start of the nucleation barrier).

Underneath are the results for five different systems. The second column shows $-\log P(n)$ for the "normal" MC (blue) and kinetic MC (dashed yellow). For all systems, we see that MC and kMC are practically identical. To better demonstrate this, the third column shows the difference $|\log P_\text{MC}(n) - \log P_\text{kMC}(n) |$.

Note that for "large" $n$, i.e. $n\gtrsim 30$, the difference increases as these clusters sizes are observed less frequently. To give an idea: in $10^5$ configurations, only $\sim200$ clusters of size 50 are found, whereas $\sim3\cdot10^4$ clusters of size 10 are found.

| System | Cluster size distribution | Difference in distribution |
| -------- | -------- | -------- |
| **HS** <br> $\beta P\sigma^3=17.5$ <br> $\beta$ l$\Delta\mu$l$\approx0.58$ |  |  |
| **Yukawa (FCC)** <br> $\beta\epsilon=8$ <br> $1/\kappa\sigma=0.01$ <br> $\beta$ l$\Delta\mu$l$\approx0.58$ |  |  |
| **Yukawa (FCC)** <br> $\beta\epsilon=81$ <br> $1/\kappa\sigma=0.01$ <br> $\beta$ l$\Delta\mu$l$\approx0.58$ |  |  |
| **Yukawa (FCC)** <br> $\beta\epsilon=8$ <br> $1/\kappa\sigma=0.04$ <br> $\beta$ l$\Delta\mu$l$\approx0.54$ |  |  |
| **Yukawa (BCC)** <br> $\beta\epsilon=81$ <br> $1/\kappa\sigma=0.40$ <br> $\beta$ l$\Delta\mu$l$=0.35$ |  |  |

### Nucleation rate

To further check the (potential) influence of the simulation method, I calculated the nucleation rate for a handful of systems using four different methods:
- Using the nucleation barrier gotten from umbrella sampling and calculating the attachment rate on top of the barrier using "normal" MC (move acceptance of $\sim30\%$).
- Using the nucleation barrier gotten from umbrella sampling and calculating the attachment rate on top of the barrier using kinetic MC (move acceptance of $\sim80-90\%$).
- Doing around 30 brute force simulations using "normal" MC.
- Doing around 30 brute force simulations using kinetic MC.

The resulting nucleation rates are given in the tables underneath. The first table gives the nucleation rate in terms of the short-time diffusion coefficient $D_0=\Delta r_{max}^2/6$, where $\Delta r_{max}$ is the maximum distance a particle can be moved in the x-, y-, and z-direction during a single move in the (k)MC simulation, and the second table gives the nucleation rate in terms of the long-time diffusion coefficient $D_l$, which I calculated from the MSD in the brute force simulations.
Note: the error in the nucleation rate gotten from the nucleation barrier is approximately a factor 3. This comes from an approximate error of $1k_BT$ in the barrier height and an error of 10% in the attachment rate.


For the last row (BCC), I also calculated the nucleation rate using brute force MD simulations: $k\sigma^5/D_l=3.7(9)\cdot10^{-6}$.

Note that for the hard-core Yukawa system with $\beta\epsilon=81$ and $1/\kappa\sigma=0.01$, which forms FCC, the nucleation rates found using brute force have a questionmark behind them. The reason for this is that I did not run these simulations for long enough, such that quite some of these simulations did not nucleate.

To give some insigth into the statistics used for getting the brute force nucleation rates, underneath a table with the predicted and average nucleation times for the specific systems treated above (plus two additional systems). The second table contains some info on the number of simulations and the max run time of the simulations.

| System | Predicted $t_{nuc}$ <br> (in $\sigma^3/6D_l$) | Average $t_{nuc}$ MC <br> (in $\sigma^3/6D_l$) | Average $t_{nuc}$ KMC <br> (in $\sigma^3/6D_l$) |
| -------- | -------- | --- | -------- |
| HS | 37.5 | 46.1 | 40.4 |
| Yukawa (FCC) <br> $\beta\epsilon=8$, $1/\kappa\sigma=0.01$ | 12.5 | 22.6 | 35.8 |
| Yukawa (FCC) <br> $\beta\epsilon=20$, $1/\kappa\sigma=0.01$ | - | 69 | - |
| Yukawa (FCC) <br> $\beta\epsilon=39$, $1/\kappa\sigma=0.01$ | - | 91 | - |
| Yukawa (FCC) <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.01$ | 40.7 | 151 | 89.5 |
| Yukawa (BCC) <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.40$ | 10.1 | 24.5 | 29.7 |

| System | # sim | max run time <br> in $\sigma^3/6D_l$ | max run time in <br> $t^{predict}_{nuc}$ ($3\cdot t^{predict}_{nuc}$) | # didn't <br> nucleate |
| -------- | ---- | ---- | --- | ----- |
| HS MC | 10 <br> 20 | 105 <br> 209 | 2.8 (0.9) <br> 5.6 (1.7) | 1 <br> 1 |
| HS KMC | 30 | 127 | 3.4 (1.1) | 2 |
| Yukawa (FCC) MC <br> $\beta\epsilon=8$, $1/\kappa\sigma=0.01$ | 10 <br> 20 | 87 <br> 292 | 7.0 (2.3) <br> 23.3 (7.8) | 2 <br> 0 |
| Yukawa (FCC) KMC <br> $\beta\epsilon=8$, $1/\kappa\sigma=0.01$ | 10 <br> 20 | 80 <br> 265 | 6.4 (2.3) <br>21.2 (7.1) | 0 <br> 0 |
| Yukawa (FCC) MC <br> $\beta\epsilon=20$, $1/\kappa\sigma=0.01$ | 10 <br> 20 | 95 <br> 315 | 1.4* <br> 4.6* | 3 <br> 0 |
| Yukawa (FCC) MC <br> $\beta\epsilon=39$, $1/\kappa\sigma=0.01$ | 10 <br> 20 | 102 <br> 340 | 1.1* <br> 3.7* | 7 <br> 10 |
| Yukawa (FCC) MC <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.01$ | 50 | 337 | 8.3 (2.8) | 12 |
| Yukawa (FCC) KMC <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.01$ | 50 | 299 | 7.3 (2.4) | 15 |
| Yukawa (BCC) MC <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.40$ | 30 | 196 | 19.4 (6.5) | 0 |
| Yukawa (BCC) KMC <br> $\beta\epsilon=81$, $1/\kappa\sigma=0.40$ | 10 <br> 20 | 125 <br> 251 | 12.4 (4.1) <br> 24.7 (8.2) | 0 <br> 0 |

Where * indicates that the average nucleation time was used instead of the predicted, as I do not have a predicted rate for these two state points.

Looking at the table above, we see that many of the "harder" potentials, i.e. $\beta\epsilon=39,81$ and $1/\kappa\sigma=0.01$, did not nucleate. The reason for this is just to do with the statistics of rare events. The nucleation times follow an exponential distribution
$$
P(t) = k e^{-kt},
$$
where $k$ is the nucleation rate. The mean of this distribution is the inverse rate $k^{-1}$.

Using this distribution, I drew 100 series of 50 samples with $k=1$. Underneath you can see the total histogram of the 100x50 samples, and the histrogram of samples per serie with $t>2.0$ and $t>3.0$. As you can see, many samples have $t>2.0$. One of the series even has 13 samples with $t>2.0$, and three series have 12 samples with $t>0$. So, assuming that they only ran for the worst case of approximately 2 nucleation times, it is possible that 12/50 MC and 15/50 KMC simulations of hard-core Yukawa particles with $\beta\epsilon=81$ and $1/\kappa\sigma=0.01$ did not nucleate within that time.

| Histogram of the total 100x50 samples | Histrogram of samples per serie with $t>2.0$ and $t>3.0$ |
| -------- | -------- |
|  |  |

---