---
# System prepended metadata

title: Binary coexistence
tags: [PhD projects]

---

# Binary coexistence

###### tags: `PhD projects`
## To do meeting 16/07
- Change fitting in mathematica. Use tangent to determine point and then get average of points before that. 
- Simulate FCC small particle crystal. 
- Change explanation in thesis on why or not minimum pressure
- Check whether the update p function has no influence
- Check free energy calculations FCC
## Current phase diagram
![Screenshot from 2025-07-04 15-48-19](https://hackmd.io/_uploads/HyxaiUHSxe.png)


![Screenshot from 2025-05-26 09-42-06](https://hackmd.io/_uploads/HJ1Do5WMeg.png)

## Motivation
We want to show that the long box method can be used to obtain the phase coexsistence of binary hard-spheres mixtures. This project is a continuation of Alessandro's work. 

## Theory
We use the long box method for this paper. This means that we put crystals at different densities and a fluid next to each other in a box (along the $z$-direction), measure the $zz$-pressure and see when it matches the pressure of the bulk crystal. 

For binary systems, there is the extra complication that apart from the crystal density, we don't know the composition of the fluid. Initially we therefore want a global composition that is inside the coexistence regime. After equilibration, we then measure the composition in the fluid. 

We start by donsidering the phase diagram of binary hard spheres with a size ratio of $\sigma_B/\sigma_A = 0.58$ (see, M. D. Eldridge, P. A. Madden, and D. Frenkel, Entropy-driven formation of a superlattice in a hard-sphere
binary mixture, Nature 365, 35 (1993)).

![Screenshot from 2024-12-11 14-12-36](https://hackmd.io/_uploads/rkzegMvNJx.png)


Good explantion of Gibbs phase rule: 
https://www.hyper-tvt.ethz.ch/fundamentals-thermo-gibbs.html#:~:text=DEGREES%20OF%20FREEDOM%20%3D%20NUMBER%20OF,F%20%3D%20C%20%2D%20P%20%2B%202

$$F = C-P+2$$
Number of constraints is $P(C-1)+2P$ (where this last term comes from P and T). Then number of degrees of fredom is $C(P-1) + 2(P-1)$.

So F is number of independent degrees of freedom. If F=1, but we have $P$ and $T$ fixed, we would have to be exactly at the correct point to a have a coexistence. Every deviation of P would have to lead to a change in T, but that is not possible, because T is fixed.

So for monodisperse, there are two independent variables, namely pressure and temperature. For two phases, $F = 1-2+2=1$ so e.g. $T(P)$. That means in NVT, two-phase coexistence is stable, but tripple point isn't (only in NVS ensemble). In NPT, P and T are both fixed, so a small disterbance of T cannot be met by a small disterbance by P, so two-phase coexistence is not stable.

For binary, there are three independent variables, $P$, $T$ and $\chi$. In NPT. So for two phases, $F=2-2+2=2$, So e.g. $P(T)$ and $\chi(T)$, but $P$ and $\chi$ do not depend on each other. So we can make a phase diagram of $P$ against $T$ or $P$ against $\chi$.

For common tangent it is useful if two constraints are left to be determined. Namely a line is associated with two degrees of freedom. In NVT, a binary system would have a common tangent plane since still the chemical potential of the two species and the pressure have to be fixed. In NPT there is one degree of freedom less to be fixed, because temperature and pressure are already set.

Note that given our NVT, or NPT or NVS ensemble a specific coexistence will form with determined densities and pressure and temperature. The only question is whether that is stable or not, i.e. whether small deviations will lead to a non-stable phase. 

And note that common tangent works in general because free energies are extensive.

In general constraints are chemical potential, temeprature and pressure (i.e all internal variables). Degrees of freedom are density, composition and entropy density. However, if you fix the intensive variable, you immediatly fix the corresponding degree of freedom. I.e. fixing $P$ fixes the densities of the two phases by the equation of state. Fixing $T$ fixes the entropy density. That also means that composition can only be a degree of freedom as long as we fix $N$, in the $\mu_\alpha\mu_\beta VT$ ensemble composition is no longer a degree of freedom. 

### Maximum packing fraction
Source: https://www.researchgate.net/profile/Daan-Frenkel/publication/27693183_The_stability_of_the_AB13_crystal_in_a_binary_hard_sphere_system/links/0deec51780b7470a42000000/The-stability-of-the-AB13-crystal-in-a-binary-hard-sphere-system.pdf
The maximum or close-packed packaging fractions for these two crystal structures are: for $\alpha = 0.58$, $\eta_{cp}(AB2) = 0.776$ and $\eta_{cp}(AB13) = 0.713$, and for $\alpha = 0.61$, $\eta_{cp}(AB2) = 0.750$ and
$\eta_{cp}(AB13) 0.685$
## Methodology
### EOS bulk crystals
We need to measure the EOS for the bulk crystals. In order to do so, for the AB2 crystal, we need to measure the ratio between the lattice constants in yz direction and the lattic constant 
#### Hexspacing
*Found on cluster. Codes in /home/users/5670772/COEX_longbox/CODES/CODE_hexspacing/ and Data in /home/users/5670772/COEX_longbox/DATA/AB2_alpha_0.58_hexspacing/ Analysis is currently found on my computer /home/rinskealkemade/Documents/LongBox_binary/AB2_alpha_0.58_hexspacing*
![Screenshot from 2024-12-18 13-48-45](https://hackmd.io/_uploads/B1Ja4HlByg.png)
In the $zy$-plane we have to lattice constants equal to 
$L_a$ in the $z$-direction and $\sqrt{3}/2 L_a$ in the $y$-direction. In the $x-$direction the lattice paramter is equal to $L_a\cdot R$, where we have to find the ratio $R$.
Note that the volume of one unit cell is thus equal to $L_a^3\sqrt{3}/2R$. 

##### Results
*Found on /home/users/5670772/COEX_longbox/DATA/AB2_alpha_0.58_hexspacing/*
To find the hexagonal spacing we change the ratio between the lattice constant in the $x$ and $z$ direction (see below, $L_a$ is lattice constant in the $z$-direction, while $L_a\cdot \text{ratio}$ is the spacing in the $x$-direction).
![Screenshot from 2025-01-07 09-53-04](https://hackmd.io/_uploads/SJZFoD9LJx.png)

| Ratio pressures | Lattice constant |
| -------- | -------- |
| ![Screenshot from 2025-01-07 09-51-21](https://hackmd.io/_uploads/HJLzsvq8Jl.png) | ![Screenshot from 2025-01-07 09-53-46](https://hackmd.io/_uploads/SkPiivcIyx.png)|

Where the fit is equal to 
![Screenshot from 2025-01-07 09-55-28](https://hackmd.io/_uploads/rJQM3w9UJe.png)

This leads to isotropic  rpessure (*Se /home/rinskealkemade/Documents/LongBox_binary/DATA/Analysis_EOS/Analysis_EOS.nb*)


| $P_{xx}/P_{zz}$ | $P_{xx}/P_{yy}$  | 
| -------- | -------- |
| ![Screenshot from 2025-07-17 11-44-59](https://hackmd.io/_uploads/S1WIUHILle.png)| ![Screenshot from 2025-07-17 11-44-53](https://hackmd.io/_uploads/H1v8UB88xx.png) |

#### AB13
We start the 13 small particles in the AB13 crystal at a semi-random distance from the big particle (the particles will not lead to overlap, but it is also not the equilibrated distance). Laura however sad that these particles will fast enough equilibrate to the correct distance. 
#### Results EOS
*Code found on ODIN in /home/users/5670772/COEX_longbox/CODES/CODE_EOS/ and data in /home/users/5670772/COEX_longbox/DATA/EOS_\* for the different crystals.*
:::danger
Put the notebook still on ODIN, currently found on my computer.
::::
We plot the EOS, and compare it with literature data  (all obtained from figure 2 from the paper https://arxiv.org/pdf/0708.0497. For FCC, we follow the Alder data form 1964, for AB2 we follow the N = 648 form the paper itself and for AB13 we follow Eldridge from 1990).
![Screenshot from 2025-01-07 12-06-19](https://hackmd.io/_uploads/H1Qact5L1x.png)
Note that our systems consisted of (FCC) 2916, (AB13) 3024, (AB2) 3000 particles. Simulations where run for $t/\tau = 10000$
::: info
We see a slight disagreement with the data for AB2 and AB13. Since FCC agress completely with the literature data and with the Speedy equation, I conclude that the code works correctly. It could be finite size effects, or maybe a wrong hexspacing. Look into this more.  
:::

::: danger
In the end match the EOS number of unit cells with the number of unit cells we put in the $x$- and $y$- direction.
:::

### Obtaining the initial box
We have two different parameters that we set beforehand, namely the global packing fraction $\eta_G$ and the global composition $\chi_G$ that should lie in the coexistence regime. 

Then we set the number of crystal particles $N_X$ and choose a range of crystal packing fractions $\eta_X$. Furthermore we set the fluid volume $V_F$ to twice $V_X$. Furthermore, we make the crystal box as square as possible, such that $L_x\approx L_y\approx L_z$

Given these four parameters, all other quantities follow. Namely we know that

- $V_X = \frac{\pi}{6\eta_X} N_X (1-\chi_X + \alpha^3\chi_X)$, with $\chi_X$ the crystal composition that is known.
- Furthermore, we know that it must be true that ![Screenshot from 2024-12-11 14-27-17](https://hackmd.io/_uploads/HJBB7GPE1l.png) such that
 ![Screenshot from 2024-12-11 14-27-46](https://hackmd.io/_uploads/B1GvQGwV1e.png)
 
 
Given these parameters, we initialize a fluid with the correct composition, number of particles, density, and such that $L_x$ and $L_y$ are equal to the $L_x$ and $L_y$ of the crystal. In order to later 'glue' the fluid to the crystal, we decrease $L_z$ by $2\sigma_A$ and $L_{x/y}$ by $0.05\sigma_A$. When the fluid is equilibrated, 'glue' the fluid to the desired crystal and let the box equilibrate.

## Performing simulations
*Code can be found on ODIN /home/users/5670772/COEX_longbox/CODES/CODE_longboxsim/*

We perform our simulation using the EDMD code of Frank for binary crystals (which we slightly altered such that it could measure pressure tensors as well as the pressure). During the simulation we measure both the pressure and the composition.

We found out that the simulation must run for quite long to obtain trustworthy composition and pressure profiles. This also implies that doing post analysis is very costly. Instead, we do the analysis during the simulation (where we equilibrate the system for $1/4^{th}$ of the maximum time and then start measuring). Simultaneously we also write out the individual not-shifted and shifted composition profiles and pressure every $1000$ steps which takes up way less diskspace, but still gives us an idea of how the composition changes over time and whether the system is equilibrated.

For now (04/03/2025) we let the FCC simulation run for $t/\tau =10^6$, while we let AB2 and AB13 run for $t/\tau =5\cdot10^6$

::: warning
Very stupidly in the code I have named the shifted composition file *composition.dat* while, the not-shifted *composition_shift.dat*. I did not want to change is during the analysis, but this should be changed afterwards. 
:::
### Measuring Composition
We measure the composition during the simulation every $t/\tau = 1$ as a function of the $z$-coordinate. The binning is spaced such that it size of one bin is equal to the lattice constant of the crystal. This means that probably there does not fit a whole number of bins in the system. To solve this we throw out the remaining smaller bin. 

Since we want the somposition profiles to overlap, we do a fourier transform of the composition profile given by 
$$
z = \sum_i^{N_b} \chi(z_i)e^{-2\pi i \cdot z_i/L_z},
$$
where $z_i$ is the $z$-coordinate of bin $i$, $\chi(z_i)$ is the corresponding composition and $L_z$ is the boxlenght in the $z$-direction.

From this $z$ we can compute
$$
\phi = \arctan\left(\frac{\text{Im}(z)}{\text{Re}(z)}\right)
$$ 
The phase associated with this profile is then given by 
$$
\text{shift} = \begin{cases}
\frac{\phi}{2\pi}L_z &\text{ if } \text{Im}(z)> 0 \text{ and } \text{Re}(z)>0\\
\frac{\pi+ \phi}{2\pi}L_z &\text{ if } \text{Re}(z)< 0\\
\frac{\phi + 3\pi/2}{2\pi}L_z &\text{ if } \text{Im}(z)< 0 \text{ and } \text{Re}(z)>0\\
\end{cases}
$$
The amount of bins that the system than has to shift is given by 
$$
\text{bin}_s = \text{floor}(\text{shift}/dR_b);
$$
with $dR_b$ the binsize. The new bin $\text{bin}_n$ is then given by 
$$
\text{bin}_n = \text{bin}_o - \text{bin}_s,
$$
with $\text{bin}_o$ the old bin. Of course we have periodic boundaries, such that if $\text{bin}_n < 0$, we set $\text{bin}_n += N_b$.

## Pressure minimum 
We saw that, contrary to the single component system, the crossing between the pressure in the $z$-direction and the bulk crystal EOS does not have to be in the minimum of the long box pressure. We explained this in the following way: 
     
![image](https://hackmd.io/_uploads/BkItgHEike.png =250x)

Say we have a binary crystal where for the statepoint indicated by the red point want to find the coexistence densities and compositions (here we use eta1 and eta2 to indicate the global point we are considering, but one could evenly well use a global composition and global density representation). If we know the phasediagram, we can already see that the system will phase seperate into two phases determined by the blue and the yellow dot which are connected by a tieline. 
     
To determine the two phases ourselves, we do longbox simulations. Here we keep the red dot the same (i.e. the global statepoint stayes the same), but we change the density of the crystal. Therefore, our measurements of the coexistences will result in points connected by different (metastable) tielines that all cross the red dot. 
![image](https://hackmd.io/_uploads/rydTgrEiye.png =250x)

Of course, I did not take into acount the effect of cyrstal deformations here. The yellow dots will thus in fact not lie exactly where they lie now because the free-energy landscape changed. However, let's assume for now that that influence is smaller than the effect of the different densities the crystal has and we can still use this cartoon.
 
If we then look at the pressure as a function of crystal density we find something like the figure below, where the lowest point in the pressure does not have to correspond to the crystal we are looking at, since the crystal associated with that point will only coexists with a fluid at a differen composition than the fluid associated with the phase point we are considering. 
![image](https://hackmd.io/_uploads/r1vx-rEoJx.png =250x)

 
We could also perform the longbox simulation in a different way. Namely, changing the global composition and density such that we always end up with the same fluid coexisting with different crystals. This would have to be via a process of trial and error, because we do not know where thes points lie, but it could in principle be done. Then our metastable tielines would look like
![image](https://hackmd.io/_uploads/BJJbbBVoyx.png =250x)
Vice versa we could also keep crystal density fixed, but change fluid density. 
In both cases we do not expect a minimum in pressure as a function of the variable we change. The reason is that in all cases we do not move along a tieline. The tieline is purely dictated by equal pressure, chemical potential and temperature. As a result pressure minimization has nothing to do with it. Only if we could walk along the tieline and choose a range of wrong fluids/crystals, would we see the wrong minimum. The problem is that that is not possible: We have no power to hold the wrong fluid, it can simply change it density/composition by growing or shrinking the crystal. For the crystal this is different since we can hold its density by the boundary conditions. However if we would move along a tieline that does not coincide with the line the crystal compoisition line (which would imply crystal composition is the same as fluid composition), we cannot move along a tieline without changing the composition of the crystal. i.e. the monodisperse case is very special. 
## Observations
During the simulations we made some observations and alterations that we will list below.

- For the FCC crystal we saw that the crystal can deform in the following way (especialy for higher packing fractions)
 ![image](https://hackmd.io/_uploads/BkX7lZqF1g.png)  
    This means that one has to look by eye at the snapshots and see whether this is not happening before including that datapoint. 
	
-	For the higher FCC coexistences, we see the following happing 
 	![Screenshot from 2025-03-07 16-43-48](https://hackmd.io/_uploads/SytINqdoJg.png)
	Should we also try to turn the crystal for the FCC crystal?
	
	The point $\chi_G =0.5$, $\eta_G=0.62$ and $\eta_X=0.62$ is a good example which can be found in *DATA/WRONG_POINTS_FOR_PAPER/FCC_MIXING*
	
	This is probably due to the fact that the box is very long. Therefore we do use a slightly less elongated box. Old datapoints can be found on ODIN */home/users/5670772/COEX_longbox/DATA/OLD/FCC_FLUID_LONGERBOX/*
	
	
- ::: warning
  For the FCC crystal we sometimes see the following happing 
	![Screenshot from 2025-03-17 11-11-18](https://hackmd.io/_uploads/BJqUSuBnkl.png)
	Leading to
	![Screenshot from 2025-03-17 11-11-33](https://hackmd.io/_uploads/HyDDBOr2yl.png)

	The point $\chi_G =0.2$, $\eta_G=0.57$ and $\eta_X=0.575$ is a good example. **This seems to be most pronounced for higher global densities.**
	What do we do with this?  
	
	For the lower $\eta_X$ and $\chi_G =0.2$, $\eta_G=0.58$ it is even worse...
	:::
- For the FCC crystal we also saw that the obtained coexistence compositions are slightly below the original phase diagram. This could be finite size effects. Therefore we increased the unit cells of the box in the $x$- and $y$-directions for the smalles composition of $0.1$. 

	:::success
	This can be explained with inaccuracies in the original paper. When we do our own free energy calculations (see below), our data points fit.
	:::
	- With more datapoints the data eems to distort faster.
	For $\eta_G = 0.54, \chi_G =0.1$ (snapshots are from $\eta_X =0.62$)

	|     | 6x6 |8 x 8 unit cells |
	| --- | -------- | -------- |
	| Pressure    | ![Screenshot from 2025-03-07 16-50-43](https://hackmd.io/_uploads/SkbJLqdiJg.png)|  ![Screenshot from 2025-03-07 16-49-16](https://hackmd.io/_uploads/Bykpr9_o1e.png) |
	| Snapshot    | ![Screenshot from 2025-03-07 16-55-50](https://hackmd.io/_uploads/B1rfPq_ske.png)|![Screenshot from 2025-03-07 16-55-19](https://hackmd.io/_uploads/r1Ugwcdo1x.png)|
	
	So probably more unit cells in the $x$- and $y$-direction would be a good idea
	
	The point $\chi_G =0.1$, $\eta_G=0.54$ and $\eta_X=0.62$ is a good example of deforming which can be found in *DATA/WRONG_POINTS_FOR_PAPER/FCC_Deforming* 

  
- In the case of the AB2 crystal there can occur a weird deformation where an hexagonal layer of the big particles is replaced by smaller particles. This occurs at statepoints quite close to the packing fraction where we expect the coexistence (see picture below).

    ![Screenshot from 2025-02-25 11-11-52](https://hackmd.io/_uploads/rk4dPMi9ke.png)
    As a results the compositoin profiles becomes almost two step (see below)
    ![Screenshot from 2025-02-25 11-09-13](https://hackmd.io/_uploads/SybkwMi51g.png)
    To solve this we rotated the crystal such that the hexagonal planes no longer lay along the $xz$-plane, but instead in the $xy$-plane. Another advantage of this rotation is that the hexagonal spacing now lies in the $z$-directions, meaning that it is no longer dictated by the box. This makes the simulation less sensitive to a slightly wrong hexspacing. 
    
    ::: warning
    These simulationes are stored in *DATA/WRONG_POINTS_FOR_PAPER/rotated_AB2_FLUID_alpha_0.580_WRONG_HEXSPACING* such that we could put it in the SI of the paper. However in the end, the hexspacing turned out to be wrong, so please **do not use them in the paper**. I let them there because it is a good example of what can go wrong!	
    :::
	
- In general we saw that both AB2 and AB13 need to equilibrate for a very long time. 

- In the case of AB13 the shifting did not work properply because the values are so close together. To solve this problem, I substracted the composition of the crystal of during the shift. This solved the problem. I do the same for the AB2 crystal. (In *DATA/WRONG_POINTS_FOR_PAPER/AB2_FLUID_SHIFT* in the notebook in that folder *compare_shift_substract_comp_X.nb* we compare the shift profiles between not substracting 2/3 (composition.dat and composition_test.dat respectively). We indeed see that the composition profiels improve slightly, although it is not that much)

- For the AB13 crystal I include slightly more particles (max 6000 instead of 4000) because one unit cell contains so many particles that we need more particles to get a reasonable size crystal. 
- We equilibrate the initial fluid for $t/\tau = 10^4$. This is probably on the short side (we saw while testing the Mansoori fluid equations that for higher packing fractions the fluid can have to equilibrate longer). We still thihnk it is okay to use this fluid though because after we glue the fluid together with the crystal we let it again equilibrate for a very long time. Moreover, we also check whether the final pressure and composition profiles are equilibrated. 
- Initially I set the number of unit cells in $x$- and $y$-direction to be equal to 6 which made the box very long. However, since the finite size effects are determined by the smalles direction, I increased the number of unit cells to 8 to make the box slightly less elongated.
::: info
All the datapoints below as well as the mathematica notebook can be foun in the folder /COEX_longbox/DATA/WRONG_POINTS_FOR_PAPER/rotated_AB2_FLUID_alpha_0.580/. 

NOTE THAT THIS DATA IS OBTAINED WITH HEXSPACING WITH TOO LITTLE DECIMALS. HOWEVER, BECAUSE THE HEXSPACING IS IN THE Z-DIRECTION THE CRYSTAL WILL BE ABLE TO ADJUST THE SPACING TO THE CORRECT SPACING, SO THIS IS NOT BAD.  
- For the AB2 crystal we see sudden jumps in the pressure (a lowering). For e.g. $\eta_G = 0.58$, $\chi_G=0.8$ and $\eta_X=0.64$ we see this
 ![Screenshot from 2025-05-07 16-30-34](https://hackmd.io/_uploads/B1_2CyKggg.png)
 I checked and this jump is in all pressure components simultaneously. I also counted the crystal layers. We start with 20, then grow quickly to 21 
 ![Screenshot from 2025-05-07 16-32-08](https://hackmd.io/_uploads/SJbWJxYxgx.png)
 Then we stay there for a very long time, although sometimes a layer tries to grow (see below)
 ![Screenshot from 2025-05-07 16-23-48](https://hackmd.io/_uploads/HyVMygtxlx.png)
 After which after $3.5\cdot 10^6$ timesteps suddenly we grow the extra layer
 ![Screenshot from 2025-05-07 16-37-38](https://hackmd.io/_uploads/rylrxxtxxe.png)
 In *Molecular simulation of homogeneous crystal nucleation of AB2 solid phase from a binary hard sphere mixture - Praveen Kumar Bommineni; Sudeep N. Punnathanam* they say that AB2 has a lower interfacial tension than FCC. They say: 
 "**The nucleation of the AB2 phase has both a lower free-energy barrier and higher particle attachment rates when compared with the FCC phase. The favoring of the nucleation of the AB2 phase is attributed to the absence of difference in its composition with that of the fluid phase.**" 
 But if this is the case, why is the interface for AB2 way sharper than for FCC...
 
 

| AB2 $\eta_G =0.595$                                                                | AB2 $\eta_G =0.58$                                                                | FCC $\eta_G =0.53$, $\chi_G=0.1$                                                  |
| --------------------------------------------------------------------------------- | --------------------------------------------------------------------------------- | --------------------------------------------------------------------------------- |
| ![Screenshot from 2025-05-13 11-41-52](https://hackmd.io/_uploads/ByrxV9g-gg.png)| ![Screenshot from 2025-05-13 09-46-05](https://hackmd.io/_uploads/HynBFdxWll.png) | ![Screenshot from 2025-05-13 09-45-55](https://hackmd.io/_uploads/HJmIt_eWgx.png) |
|![Screenshot from 2025-05-13 11-40-06](https://hackmd.io/_uploads/B17pmcgZge.png)|![Screenshot from 2025-05-13 11-42-41](https://hackmd.io/_uploads/Hy77N5eZll.png)|![Screenshot from 2025-05-13 11-43-19](https://hackmd.io/_uploads/H1arE9gZee.png)|


- For the AB13 we see the same thing happen, although the crystal jumps less. Below we plot the measured pressure for a system at $\chi_G=0.95$ and $\eta_G=0.55$ together with the number of layers the crystal has
  ![EOS_layers](https://hackmd.io/_uploads/ByXXqx9gel.png)
On top of that there are also a lot of defects (see composition profile for the same statepoint and $\eta_X=0.555$)
![Screenshot from 2025-05-08 11-46-22](https://hackmd.io/_uploads/HJCY6g9gll.png)

**I put this in the paper. Notebooks to obtain the figures in the paper can be found in /home/rinskealkemade/Documents/LongBox_binary/DATA/WRONG_POINTS_FOR_PAPER/Analysis_composition_\*_0.58_sturdy_interface_chi_G_\*_eta_G_\*.nb. The corresponding figures for the good interfaces can be found in the folders of AB13 and AB2 (called OUTPUT_snapshot\*)**
:::
- For the FCC crystals close to/in the AB2 regime we see spontaneous formation of AB2. This is for $\chi_G = 0.6$, $\eta_G=0.62, $\eta_X = 0.65$.


| All particles| Big particles| Small particles |
| -------- | -------- | -------- |
| ![Screenshot from 2025-07-02 10-01-22](https://hackmd.io/_uploads/H1Er_DMBll.png) | ![Screenshot from 2025-07-02 10-01-18](https://hackmd.io/_uploads/BkDEuwMrlx.png)     | ![Screenshot from 2025-07-02 10-01-14](https://hackmd.io/_uploads/HyxLOPMrgg.png)|
- For FCC at $\eta_G =0.58$, $\chi_G=0.4$ and $\eta_X = 0.69$, we see this very cool kind of deformation of an AB crystal: (Note that this is far from the 'good' crystal, so probably this makes the deformed crystal a little bit more happy)


| Column 1 | Column 2 |
| -------- | -------- |
| ![Screenshot from 2025-07-08 11-24-51](https://hackmd.io/_uploads/SJNPEP9Hel.png)| ![Screenshot from 2025-07-08 11-24-40](https://hackmd.io/_uploads/H1qPNw9Sge.png)| 

*In /home/rinskealkemade/Documents/LongBox_binary/DATA/WRONG_POINTS_FOR_PAPER/FCC_planer_defect this data point can also be found, since we want to put the defect in the paper*
- The small FCC crystal can have weird composition profiles (see e.g. $\
eta_G = 0.52$ and $\chi_G=0.998$. In hink these are finite size effects that can be attributed to the small number of big particles. )
![Screenshot from 2025-08-18 14-19-36](https://hackmd.io/_uploads/B1oFq9etge.png)

### Turned AB13 EOS
*Found on ODIN /home/users/5670772/COEX_longbox/TEST/DATA_CHECK_AB13/*
After we decided we wanted to turn AB13, I made a new code that generated the crystal. 

To test the code, we measured the EOS again, see below, where dots are measured in our simulation and the line comes from the earlier mentioned paper. We did not wait for the highest packing fractions to be finished because the data clearly agrees.
![Screenshot from 2025-05-15 16-33-10](https://hackmd.io/_uploads/SJWSi_QZge.png)




### Results turned crystal
*Found on ODIN /home/users/5670772/COEX_longbox/DATA/EXPLORE_20250602/*
::: danger
Note that the AB2 simulations were performed with the wrong hexspacing, i.e. too little decimals. I.e. don't use this data for real
:::
#### AB2
| $\chi_G$ | $\eta_G$ | Composition                                                                       | Pressure | Remarks    |
| -------- | -------- | --------------------------------------------------------------------------------- | -------- | --- |
| 0.8     | 0.58    |![Screenshot from 2025-05-26 10-00-41](https://hackmd.io/_uploads/By23Ji-Gel.png)|![Screenshot from 2025-05-26 10-02-50](https://hackmd.io/_uploads/rJjVeoWzlg.png)|Equilibates quite fast, run for $t/\tau= 4\cdot 10^6$|
|     | 0.585    |![Screenshot from 2025-05-26 10-03-50](https://hackmd.io/_uploads/ryY_loZfgl.png)|![Screenshot from 2025-05-26 10-04-37](https://hackmd.io/_uploads/Hywjxs-fll.png)|Equilibates slower, the initial pressure is too high|
|    | 0.59    |![Screenshot from 2025-05-26 10-06-22](https://hackmd.io/_uploads/ByTbWo-Gxe.png)|![Screenshot from 2025-05-26 10-05-28](https://hackmd.io/_uploads/r12JbsWMlg.png)||
||0.595|![Screenshot from 2025-05-26 10-07-04](https://hackmd.io/_uploads/SJIVWoWGee.png)|![Screenshot from 2025-05-26 10-07-24](https://hackmd.io/_uploads/r1iS-j-Gge.png)|

So we see that AB2, the higher global densities needs to grow a lot of crystal. It would maybe make more sense to increase the composition a bit
#### AB13 
| $\chi_G$ | $\eta_G$ | Composition                                                                       | Pressure | Remarks    |
| -------- | -------- | --------------------------------------------------------------------------------- | -------- | --- |
| 0.95     | 0.52     | ||Too little crystal, pressure has to increase a lot and crystal melts|
||0.53     | ||Too little crystal, pressure has to increase a lot and crystal melts|
|      | 0.54     | ![Screenshot from 2025-05-26 09-45-37](https://hackmd.io/_uploads/HkXNn9-zxx.png) | ![Screenshot from 2025-05-26 09-47-33](https://hackmd.io/_uploads/BJIj29Wzel.png)|  Not too many defects, run for $t/\tau~ 6\cdot10^6$. However, we are still really at the low side of the phase boundary. Maybe overall density should be slightly higher  |
|    | 0.55           | ![Screenshot from 2025-05-26 09-48-20](https://hackmd.io/_uploads/Bk6R3qbfge.png)|  ![Screenshot from 2025-05-26 09-48-53](https://hackmd.io/_uploads/rJ7x69WGxl.png)| Initial pressure is very far from final pressure, crystal has too grow too much   |
||0.56|![Screenshot from 2025-05-26 09-50-25](https://hackmd.io/_uploads/HkbL69-fge.png)|![Screenshot from 2025-05-26 09-50-54](https://hackmd.io/_uploads/H1buT5bfxx.png)|Same story, too high overall density. Leads to fast equilibration, but bad fluid statistics|
|0.96|0.52|||Too little crystal, pressure has to increase a lot and crystal melts|
||0.53|||Already a lot better, but the crystal still has to decrease. |
||0.54|||This actually seemse really good, we essentially start at the correct pressure|
||0.55|||This one starts at a too high global density, pressure has to decrease|
|0.965|0.52|||This one starts at a too low global density, pressure has to increase|
||0.53|||This one starts at a too low global density, pressure has to increase|
||0.54|![Screenshot from 2025-05-26 09-54-14](https://hackmd.io/_uploads/rJvNR5bMeg.png)|![Screenshot from 2025-05-26 09-54-54](https://hackmd.io/_uploads/B1fw05ZMxe.png)|We start almost exactly right pressure wise|
|0.55|||Small particles in the fluid crystallize too, i.e. we are in the FCC-AB13 coexistence|

So for $\chi_G = 0.95$ the sweetspot is around $\eta_G=0.54$. Below it has to grow the crystal a lot and above it needs to decrease it a lot. 


#### Remark
For $\eta_\chi=0.965$ and $\eta_G=0.55$ we see this. Here we see that the lower points are associated with a crystallized fluid, i.e. the small particles cyrstallize. 
![Screenshot from 2025-06-02 16-56-17](https://hackmd.io/_uploads/HJ5hsNoMxx.png)

## Final datapoints
For every datapoint where we see a correct composition, we write down what the minimum and maximum $\eta_X$ are.
:::warning
Do we really want to fit a low order polynomial through the data. It feels like sometimes the composition might not be a simple quadratic fit, but something with mroe behaviour.  

This is e.g. $\eta_G=0.61$, $\chi_G=0.6$ for FCC.
![Screenshot from 2025-09-22 12-04-29](https://hackmd.io/_uploads/SJ3NJsCile.png)
I can also not trivailly say why it should be such a simple function in general.
:::

### FCC
| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           | Finished analysis    | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------ |
| $0.005$     | $0.52$         | $8$           | $0.525$            | $0.615$            | $0.525$       | $0.57$        | :white_check_mark: |    :ballot_box_with_check: |        |
*******************************************
| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           | Finished analysis  | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------ |
| $0.01$      | $0.51$         | $8$           | $0.525$            | $0.615$            | $0.525$       | $0.57$        | :white_check_mark: |      :ballot_box_with_check:|        |

*******************************************
| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           |  Finished analysis    | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------ |
| $0.03$      | $0.51$         | $8$           | $0.525$            | $0.615$            | $0.525$       | $0.56$        | :white_check_mark: |    :ballot_box_with_check:  |        ||

*******************************************

| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           | Finished analysis | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------ |
| $0.05$      | $0.51$         | $8$           | $0.525$            | $0.625$            | $0.525$       | $0.56$        | :white_check_mark: |      :ballot_box_with_check:| Simulations with $\eta> 0.56$ did not finish       |
|             | $0.52$         | $8$           | $0.525$            | $0.625$            |          $0.525$       | $0.57$             | :white_check_mark: | :ballot_box_with_check:    |Datapoint included that did not run completly        |

*******************************************
| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           |     Finished analysis| Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------ |
| $0.075$     | $0.52$         | $8$           | $0.525$            | $0.615$            | $0.525$       | $0.57$        | :white_check_mark: | :ballot_box_with_check:      |        |

*******************************************
| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           | Finished analysis    | Reason                                                                                                        |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------------------------------------------------------------------------------------------------------------- |
| $0.1$       | $0.52$         | $8$           | $0.545$            | $0.635$            | $0.545$       | $0.590$       | :white_check_mark: |:ballot_box_with_check:       | For higher $\eta$, crystal deformed  and very little crystal                                                  |
|             | $0.53$         | $8$           | $0.545$            | $0.635$            | $0.545$       | $0.61$        | :white_check_mark: | :ballot_box_with_check:      | Initial fluid at $\eta=0.545$ is partly crystalized, but it melts again, so system is equilibrated in the end |
|             | $0.54$         | $8$           | $0.545$            | $0.635$            | $0.545$       | $0.615$       | :white_check_mark: | :ballot_box_with_check:      | Same                                                                                                          |

*******************************************

| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           | Finished Anlaysis    | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------ |
| $0.15$      | $0.53$         | $8$           | $0.55$             | $0.63$             | $0.55$        | $0.610$       | :white_check_mark: | :ballot_box_with_check:      |        |
|             | $0.54$         | $8$           | $0.55$             | $0.63$             | $0.55$        | $0.615$       | :white_check_mark: | :ballot_box_with_check:      |        |
|             | $0.55$         | $8$           | $0.55$             | $0.63$             | $0.56$        | $0.625$       | :white_check_mark: | :ballot_box_with_check:      |  $\eta=0.565$ ran twice, so the pressure tensor has twice as much inputs. The average is based on the last simulation      |
*******************************************

| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           | Finished Analysis    | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------ |
| $0.2$       | $0.55$         | $8$           | $0.57$             | $0.67$             | $0.57$        | $0.635$       | :white_check_mark: |  :ballot_box_with_check:     | Same   |
|             | $0.56$         | $8$           | $0.57$             | $0.67$             | $0.57$        | $0.645$       | :white_check_mark: |  :ballot_box_with_check:     | Same   |
|             | $0.595$         | $8$           | $0.57$             | $0.67$             | $0.57$        | $0.65$        | :white_check_mark: |:ballot_box_with_check:       | Lower $\eta$ have a lot of defects   |
*******************************************

| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           | Finished Analysis    | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------ |
| $0.3$       | $0.56$         | $6$           | $0.57$             | $0.71$             | $0.59$        | $0.655$       | :white_check_mark: | :ballot_box_with_check:      | Same   |
|             | $0.57$         | $8$           | $0.57$             | $0.71$             | $0.59$        | $0.665$       | :white_check_mark: | :ballot_box_with_check:  | Same   |
|             | $0.58$         | $6$           | $0.57$             | $0.71$             | $0.61$        | $0.670$       | :white_check_mark: |  :ballot_box_with_check:      | Lower $\eta$ have a lot of defects |


*******************************************
| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           | Finished Analysis       | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | ----------------------- | ------ |
| $0.35$      | $0.58$         | $8$           | $0.60$             | $0.70$             | $0.60$        | $0.675$       | :white_check_mark: | :ballot_box_with_check: |        |
                
*******************************************
| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           | Finished Analysis    | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------ |
| $0.4$       | $0.58$         | $8$           | $0.620$            | $0.7$              | $0.625$       | $0.68$        | :white_check_mark: | :ballot_box_with_check:    |        |
|             | $0.59$         | $8$           | $0.620$            | $0.7$              | $0.625$       | $0.68$        | :white_check_mark: | :ballot_box_with_check:    |        |
|             | $0.60$         | $8$           | $0.620$            | $0.7$              | $0.635$       | $0.68$        | :white_check_mark: |:ballot_box_with_check:|        |
|             | $0.61$         | $8$           | $0.625$            | $0.7$              | $0.645$       | $0.68$        | :white_check_mark: |  :ballot_box_with_check:   |        |


*******************************************
| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished           | Finished Analysis    | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------ |
| $0.5$       | $0.60$         | $8$           | $0.65$             | $0.71$             | $0.650$       | $0.70$        | :white_check_mark: | :ballot_box_with_check:    |        |
|             | $0.61$         | $8$           | $0.65$             | $0.71$             | $0.655$       | $0.70$        | :white_check_mark: | :ballot_box_with_check:    |     A lot of simulations finished early, but we rean them extremely long, so it is okay   |


*******************************************
| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used | Finished  | Finished Analysis    | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------- | ------------------ | --- | ------ |
|  $0.55$  |     $0.61$         | $8$           | $0.65$             | $0.71$             |              $0.66$       | $0.70$                 |:white_check_mark: | :ballot_box_with_check:     |   Datapoint included that did not run completly     |
*******************************************
| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used | $\eta_F$ used            | Finished | Finished Analysis    | Reason |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------- | ------------------------ | -------- | --- | ------ |
| $0.60$      | $0.61$         | $8$           | $0.65$             | $0.71$             | $0.67$       | $0.70$      |         :white_check_mark: | :ballot_box_with_check:    |      |


### AB2

*******************************************
:::info
Since we only checked the axis ratio from $\eta=0.62$, in the analysis we only include values of $\eta>0.625$
:::
| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | Finished           | Finished analysis | Remark    |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | ------------------ | ----------------- | --- |
| $0.750$     | $0.560$        | $10$          | $0.610$            | $0.660$            | :white_check_mark: | :ballot_box_with_check: ||
|             | $0.570$        | $10$          | $0.610$            | $0.660$            | :white_check_mark: | :ballot_box_with_check:   | Datapoint included that did not run completly    |
|             | $0.580$        | $10$          | $0.610$            | $0.660$            | :white_check_mark: | :ballot_box_with_check:  |     |
|             | $0.590$        | $10$          | $0.6275$           | $0.660$            | :white_check_mark: |:ballot_box_with_check:               |     |
| $0.800$     | $0.575$        | $10$          | $0.630$            | $0.665$            | :white_check_mark: |:ballot_box_with_check:       |     |
|             | $0.580$        | $10$          | $0.630$            | $0.665$            | :white_check_mark: |  :ballot_box_with_check:   |     |
$0.810$       | $0.58$         | $10$           | $0.630$             | $0.665$             |  :white_check_mark:|:ballot_box_with_check: | Datapoint included that did not run completly|
|       | $0.585$         | $10$           | $0.630$             | $0.665$             |    :white_check_mark:|:ballot_box_with_check:    | |
| $0.825$       | $0.58$         | $10$           | $0.640$             | $0.665$             |:white_check_mark: | :ballot_box_with_check:  | 
|        | $0.585$         | $10$           | $0.640$             | $0.665$             | :white_check_mark: | :ballot_box_with_check:  |$\eta=0.66$ is maesured twice. However, the average only goes over last one. so maybe we should remove the pressure tensor here?| 
|        | $0.59$         | $10$           | $0.640$             | $0.675$             |:white_check_mark:    | :ballot_box_with_check:  | |
|       | $0.595$         | $10$           | $0.645$             | $0.675$             |  :white_check_mark:    | :ballot_box_with_check: |
| $0.850$       | $0.585$         | $10$           | $0.655$             | $0.680$             |  :white_check_mark:|:ballot_box_with_check:| 
|       | $0.590$         | $10$           | $0.655$             | $0.680$             | :white_check_mark:|:ballot_box_with_check: |

### AB13
| Composition | Global density | Unit cells xy |$\eta_L$ simulated | $\eta_H$ simulated | Finished    | Finished Analysis |    Reason  |
| ----------- | -------------- |-------------- | ------------------ | ------------------ | --- | ------- | --- |
| $0.9400$     | $0.530$        | $10$               | $0.565$            |      $0.590$ |   :white_check_mark:    |:ballot_box_with_check: |Datapoint included that did not run completly
|     | $0.540$        | $10$               | $0.565$            |      $0.590$ |   :white_check_mark:     |:ballot_box_with_check:|Datapoint included that did not run completly|
| $0.9500$     | $0.530$        | $10$               | $0.565$            |      $0.590$ |   :white_check_mark:     |:ballot_box_with_check:|Datapoint included that did not run completly
|             | $0.540$        | $10$               | $0.565$            |     $0.590$ |:white_check_mark:|:ballot_box_with_check:|Datapoint included that did not run completly|
| $0.9525$      | $0.540$        | $10$               | $0.570$            |      $0.610$ |   ||
| $0.9550$      | $0.540$        | $10$               | $0.570$            |      $0.610$ |      ||
| $0.9575$      | $0.540$        | $10$               | $0.570$            |      $0.610$ |  :white_check_mark:   |:ballot_box_with_check:|
| $0.9600$      | $0.540$        | $10$               | $0.575$            |      $0.615$ |  :white_check_mark:    |:ballot_box_with_check:|
| $0.9625$      | $0.540$        | $10$               | $0.575$            |      $0.615$ |   :white_check_mark:   |:ballot_box_with_check:|
| $0.965$      | $0.540$        | $10$               | $0.575$            |      $0.5975$ |   :white_check_mark: ||



### FCC small

| Composition | Global density | Unit cells xy | $\eta_L$ simulated | $\eta_H$ simulated | $\eta_L$ used  | $\eta_F$ used  | Finished    | Finished final analysis |
| ----------- | -------------- | ------------- | ------------------ | ------------------ | -------------- | -------------- | --- | ---------------------------------------------------------------------------------------- |
| $0.975$       | $0.53$        | $8$           | $0.54$            |  $0.6100$           |$0.545$|$0.5625$| :white_check_mark:  | :ballot_box_with_check:    |
| $0.98$       | $0.53$        | $8$           | $0.5575$            |  $0.6100$           |$0.540$|$0.5650$|:white_check_mark:   |   :ballot_box_with_check:  |
| $0.985$       | $0.53$        | $8$           | $0.5575$            |  $0.6100$           |$0.540$|$0.5650$|  :white_check_mark:   | :ballot_box_with_check:    |
| $0.990$       | $0.53$        | $8$           | $0.5575$            |  $0.6100$           |$0.540$|$0.5650$|:white_check_mark:   | :ballot_box_with_check:    |
| $0.995$       | $0.53$        | $8$           | $0.5575$            |  $0.6100$           |$0.540$|$0.5650$|  :white_check_mark:  | :ballot_box_with_check:     |
| $0.996$       | $0.53$        | $8$           | $0.5575$            |  $0.6100$           |$0.540$|$0.5650$| :white_check_mark:  |  :ballot_box_with_check:  |
| $0.997$       |  $0.53$        | $8$           | $0.5575$            |  $0.6100$           |$0.540$|$0.5650$| :white_check_mark:    |   :ballot_box_with_check:|
| $0.998$       | $0.53$        | $8$           | $0.5575$            |  $0.6100$           |$0.540$|$0.5650$|  :white_check_mark:   |   :ballot_box_with_check:|


::: danger
Think about different fitting strategies, tangens might not be the best one
:::

#### REDO


| X   | $\chi_G$ | $\eta_G$ | $\eta_X$     |Date     |
| --- | -------- | -------- | --- | --- |
| AB2 | $0.75$     |  $0.58$    | $0.6125$     |**2025-06-10**  Started X    |
| AB2 | $0.85$     |  $0.59$    | $0.6675$     |**2025-06-10** Started X    |
| AB2 | $0.80$     |  $0.575$    | $0.6375$     |**2025-06-16**  Started X |
| AB2 | $0.75$     |  $0.58$    | $0.6350$     |**2025-06-16** Started X |  
| AB2 | $0.825$     |  $0.585$    | $0.6450$     |**2025-06-16** Started X |
| AB2 | $0.825$     |  $0.585$    | $0.64$     |**2025-06-10**  Started X  |
| AB2 | $0.750$     |  $0.590$    | $0.6275$     |**2025-06-16** Started X |
| AB13 | $0.955$ |  $0.540$    | $0.6025$    |**2025-06-18**  Started X  |
| AB13 | $0.945$ |  $0.540$    | $0.585$    |**2025-06-30**  Started X  |
| AB2 | $0.75$ |  $0.56$    | $0.62$    |**2025-07-01** Started X  |
| AB2 | $0.81$ |  $0.585$    | $0.655$    |**2025-07-10** Started X  |
| AB2 | $0.825$ |  $0.59$    | $0.6525$    |**2025-07-10** Started X|

## Check phase diagram with theory
*Found on my computer/home/rinskealkemade/Documents/LongBox_binary/TEST/EOS_Binary_Fluid Mansoori_Fluid.nb*

To check our data points, we also check our results with theory. For this we use the existing free energies for both the FCC crystal and the binary fluid. 

### Fluid
Based on *Equilibrium Thermodynamic Properties of the Mixture of Hard Spheres by G. A. Mansoori, N. F. Carnahan, K. E. Starling, and T. W. Leland Jr.*.

We know that the mansoor compressibility is given by ($\phi$ is packing fraction) (see also \url{https://www.researchgate.net/profile/Daan-Frenkel/publication/27693183_The_stability_of_the_AB13_crystal_in_a_binary_hard_sphere_system/links/0deec51780b7470a42000000/The-stability-of-the-AB13-crystal-in-a-binary-hard-sphere-system.pdf})
$$
Z^\text{Mans}(\phi) = \frac{\beta P}{\rho} = \frac{1+\phi+\phi^2-3\phi(y_1+y_2\phi)-y_3\phi^3}{(1-\phi)^3}
$$
where 
![Screenshot from 2025-02-26 15-45-01](https://hackmd.io/_uploads/Byjgtsh91g.png)

For later we also write $\eta$ in terms of $\rho$, using that $\chi = N_B/N$ and $1-\chi = N_A/N$ and $\sigma_B = \alpha\sigma_A$, we find
$$
\begin{align*}
\phi &= \phi_A +\phi_B \\
&= \frac{\pi}{6}(\sigma_A^3 N_A/V+ \sigma_B^3 N_B/V)\\
&= \frac{\pi}{6}((1-\chi)+\alpha^3\chi)\rho\sigma_A^3 
\end{align*}
$$

and thus
$$
\begin{align*}
\rho &= \frac{\phi}{\frac{\pi}{6}((1-\chi)+\alpha^3\chi)\sigma_A^3 }
\end{align*}
$$
From this we can find the Mansoori pressure in terms of the packing fraction
$$
\frac{\beta P}{\rho} =\beta P\sigma_A^3\cdot \frac{\pi}{6\phi}((1-\chi)+\alpha^3\chi)=Z^\text{Mans}(\phi)
$$
$$
\beta P^\text{Mans}(\phi)\sigma_A^3 = \frac{6}{\pi}\phi \frac{Z^\text{Mans}(\phi)}{(1-\chi)+\alpha^3\chi}
$$
We checked that this pressure corresponds to the CS pressure and it does 
![Screenshot from 2025-02-26 15-59-20](https://hackmd.io/_uploads/BkxP2sn51e.png)
#### Check with EOS fluid measurements
*Found on ODIN /home/users/5670772/COEX_longbox/TEST/EOS_Binary_Fluid*

We also measure the EOS with actual simulations for different $\chi_G$ and $\eta_G$. These measurements agree quite well, altough there turns out to be a systematic difference. Moreover, for $\eta = 0.50$ we see that the system starts crystallizing. 

| Theory with experiments | Difference theory and experiments|
| -------- | -------- |
|![Screenshot from 2025-03-07 15-03-59](https://hackmd.io/_uploads/Hymyau_s1x.png) | ![Screenshot from 2025-03-07 15-05-13](https://hackmd.io/_uploads/BJr46dOoyl.png)|

Here different colors indicate different $\chi_G$ (from 0 to 1 and red to blue), lines indicate theory, and points indicate measurements from simulations.


### Helmholtz Free energy
To obtain the free energy of the crystal we use the fact that the excess free energy is given by 
$$
\frac{\beta F^{ex}}{N} = \int_0^\rho\left(\frac{\beta P}{\rho'}-1\right)\frac{d\rho'}{\rho'} = \int_0^\phi\left(Z^{Mans}(\phi)-1\right)\frac{d\phi'}{\phi'}
$$
Note that the $-1$ is there in the integral because we consider the term w.r.t an ideal gas which has $Z=1$, if we would not substract that we would not have to add the ideal term later.

::: danger
If we compute this from our obtained $Z^{Mans}(\phi)$ we obtain the same expression as the origional Mansoori paper (where $\xi =\phi$)
![Screenshot from 2025-02-26 16-08-53](https://hackmd.io/_uploads/r1O9Rsh9kl.png)
However, Daan Frenkel has an expression that misses this first term.
![Screenshot from 2025-03-03 14-16-01](https://hackmd.io/_uploads/ByQsjX7syg.png)
Why is this different? For now I take the first expression, since that is also what I obtain when I do the calculations. This also makes sense becuase at $\phi= 0$ the excess free energy should go to zero, which it doesn't for the DF term.
:::
And since the free energy of a mixture is given by 
$$
\beta F^{Id}/N = \log\left(\frac{6\phi}{\pi((1-\chi)+\alpha^3\chi)}\right)-1+ \chi\log\chi +(1-\chi)\log(1-\chi) 
$$
such that our total free energy is given by the sum of the two $F= F^{Id} + F^{Ex}$.

#### Difference Mansoori free energy and free energy from measurements
From the obtained EOS for binary particles, we can also compute the free energy. To do so, we fit the EOS and integrate that as a function of the packing fraction. (where we force the lines to go through zero and we force the linear component to equal the ideal gas law).

Below we show the fitted functions together with the Mansoori theory. On the right we plot again the difference between the two, but now for the fitted lines.

| Theory with experiments | Difference theory and experiments|
| -------- | -------- |
|![Screenshot from 2025-03-07 15-12-42](https://hackmd.io/_uploads/B1T11tujJg.png)(colored lines are Mansoorifor different $\chi$, black dashed is us) |![Screenshot from 2025-03-07 15-17-52](https://hackmd.io/_uploads/r19Qetds1x.png)|

Here different colors indicate different $\chi_G$ (from 0 to 1 and red to blue), lines indicate theory, and points indicate measurements from simulations.
##### Virial coefficient

Paper on second virial coefficient binary mixture https://journals.aps.org/pre/pdf/10.1103/PhysRevE.80.051122

It says 
$$
B_n(s) = \sum_{i = 0}^2{n \choose i} B_n^i(s)(1-\chi)^{n-i}\chi^i
$$
where $n$ is the virial coefficient we are considering, and $s=\sigma_1/\sigma_2$.
Then $B_2^0 = \frac{2\pi}{3}$ $B_2^1 = \frac{2\pi}{3}\left(\frac{1+s}{2}\right)^3$ adn  $B_2^2 = \frac{2\pi}{3}s^3$. From the paper, we also include $B_3$.

Given the expression for the pressure, we then want to fit the integrand, i.e. fit
$$
\text{integrand}(\phi, \chi) = \frac{\beta P \sigma^3/\rho - 1}{\phi}
$$
We perform the same transformations on the measured data and then fit the integrand (see below for the fit and the datapoints). 
![Screenshot from 2025-03-11 12-02-49](https://hackmd.io/_uploads/rkhP_cao1x.png)

We then integrate these fitted functions to obtain the free energy.We compute the difference vbetween the measured and the analytical free energy. Note that this is a very small difference
| Difference free eenergy |Total free energy|
| -------- | -------- |
|![Screenshot from 2025-03-11 12-03-24](https://hackmd.io/_uploads/rk9t_qTokx.png)|![Screenshot from 2025-03-11 12-03-43](https://hackmd.io/_uploads/H1fod5asye.png)|

### Gibbs free energy
Then we want to obtain the Gibbs free energy. We know
$$
\begin{align*}
\beta G/N &= \beta F/N + \beta P V/N\\
& = \beta F/N + Z^\text{Mans}\\
\end{align*}
$$
The above equation we have as a function of $\phi$. However, we want to write it down in terms of a certain pressure, because in the Gibbs ensemble, we control the pressure.  Let's say we set $\beta P \sigma_A^3 = A$. Then we want to solve for a spectra of $\chi$, for which $\phi$
$$
\beta P \sigma_A^3 [\phi, \chi]= A
$$


In the end we want to solve for a specific pressure for which $\chi_F$ the tangent of the gibbs free energy is equal to the gibbs free energy of the crystal at $\chi_X=0$ (for FCC), so when 
$$
G(\chi_F) -G'(\chi_F)\cdot\chi_F = G_X 
$$

We can shift $G$ by a random function $a\chi$, since this term cancels between $G(\chi_F)$ and $G'(\chi_F)$.
### Crystal
To obtain the crystal Gibbs free energy, we use the Frenkel and Ladd method. We know that the pressure of the solid is given by 
$$
P\beta\sigma^3 =\frac{6\phi}{\pi} \left(\frac{3}{1-\frac{6\phi}{\sqrt{2}\pi}}- \frac{0.5921(\frac{6\phi}{\sqrt{2}\pi}-0.7072)}{\frac{6\phi}{\sqrt{2}\pi}-0.601}\right),
$$
Then we know that 
$$
\frac{\beta F^{ex}(\phi)}{N} = \frac{\beta F^{ex}(\phi_\text{ref})}{N}+\int_{\phi_\text{ref}}^\phi \frac{Z(\phi')}{\phi'}d\phi'
$$
Note that here we do not have the $-1$ term because now the excess energy free energy at the reference point is explicitely included. From the primitive function of P (which we call $P^F$) we can then comput the excess free energy of the FCC crystal. We know that (https://pubs.aip.org/aip/jcp/article/112/12/5339/473564/Finite-size-corrections-to-the-free-energies-of)
$$
\beta F^{ex}(\rho)/N= 5.91889 \text{ with }\rho = 1.04086 \text{ or } \phi = 0.544993 
$$

Then we find a Helmholtz free energy of (using that $\beta F^{id}/N= \log(\rho)-1$)

$$
\begin{align*}
\frac{\beta F}{N} &= \frac{\beta F^{ex}(\phi)}{N} + \frac{\beta F^{id}(\phi)}{N} \\\\\
&=(5.91889 + P^F(\phi)-P^F(\phi_\text{ref})) + \log(6 \cdot 0.544993 /\pi)-1 
\end{align*}
$$
Andt thus a Gibbs free energy of (for a certain pressure)
$$
\begin{align*}
\frac{\beta G}{N} = \frac{\beta F}{N} + \frac{6\beta P\sigma^3}{\phi(\beta P\sigma^3)}
\end{align*}
$$
We in the end obtain the following phase diagram (black is read from the paper, red indicates our values). As we see especially for lower $\phi$ our measurements deviate from the paper values. We do so however that our crossingpoint with the $y$-axis lies closer toknown coexistence pressure for single component hard spheres of $\beta P\sigma^3=11.54$.
![Screenshot from 2025-03-07 16-16-25](https://hackmd.io/_uploads/S1NyCY_oyx.png)


::: danger
Also check this phase diagram since they might have better points 
Molecular simulation of homogeneous crystal nucleation of
AB2 solid phase from a binary hard sphere mixture 
Praveen Kumar Bommineni; Sudeep N. Punnathanam
Sounded nice, but they have the wrong size ratio
:::

Discuss:
Our free energy calulcations (black) are slightly different from our measurements yellow. However, what we also wee is that our line fits better with the pure transition from Fluid to FCC-A (red point).
![Screenshot from 2025-06-23 13-23-5](https://hackmd.io/_uploads/H1jSt28Veg.png)
This is the difference
![Screenshot from 2025-06-23 15-01-54](https://hackmd.io/_uploads/r18DlAU4le.png)
****
****
## To check
- [ ] Check whether hexspacing is correct by checking whether the pressure is really isotropic. 
- [ ] Check whether pressure tensor trace gives pressure.
- [ ] Check whether composition profiles are really computed correctly. Mainly because AB13 is slightly blow the line.
## To do
- [ ] Put EOS notebook on cluster
- [ ] Look which AB13 simulations give coexistences and stop the others (and see whether they could do with less simulation time?)
- [ ] Look at increased number cells of FCC and see whether it works
- [ ] Increase number of unit cells for AB2, if we see that this turning of the crystal worked and see whether we can do with less simulation time.
- [ ] Also increase number of unit cells for AB13 and FCC (for the statepoints that we did not yet look at)
- [ ] EOS with same number of unit cells along the short axis?

## To keep in mind
 - https://link.springer.com/article/10.1140/epjp/i2013-13010-8 says: A second problem, that is more relevant for Monte Carlo than for Molecular Dynamics simulations, is that two-phase systems equilibrate slowly. In particular, coexistence between two phases requires that the pressures are equal. In Molecular Dynamics simulations, the rate of pressure equilibration is determined by the speed of sound —typically such equilibration is rapid. In contrast, in Monte Carlo simulations, pressure equilibrates through diffusion, and this is slower.
 The reason that we cannot use and $NPT$ ensemble simulation is because
 ![Screenshot from 2025-05-08 13-00-57](https://hackmd.io/_uploads/BJwlJz5gxx.png)

## Questions
- So surface tension has influence on the potential. Surface stress is about the pressure. Surface tension is always positive, however, surface stress for a specific plane is not. As a result the pressure inside a nucleus can be positive. However, this does not mean that the interface becomes very big, because still the surface tension is positive. 
- Math of why interace in the z-direction has no influence
## Literature
- Binary coexistence: https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-023-42713-5/MediaObjects/41467_2023_42713_MOESM1_ESM.pdf


## To discuss

