# Pressure inside a solid nucleus (with Carlos)

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

---

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

Marjolein, Laura, Frank

---

## Theoretical background

### Young-Laplace equation (+ extension)
Let's assume that we have a solid nucleus in contact with the liquid in a system in equilibrium at constant $N$, $V$, $T$. We can write the Helmholtz free energy for this system as
$$ F = N_s \mu_s - P_s V_s + N_l \mu_l - P_l V_l + N_i \mu_i + \gamma A,
$$
where $N_s$, $N_l$, and $N_i$ are the number of particles of respectively the solid, liquid, and interface, $\mu_s$, $\mu_l$, and $\mu_i$ are the chemical potentials, $P_s$ and $P_l$ are the pressures of the solid and liquid, $\gamma$ is the interfacial free energy, and $A$ is the surface area of the solid-liquid interface. In equilibrium, the chemical potential must be homogeneous, i.e. $\mu_s=\mu_l=\mu_i\equiv\mu$. Consequently, we can simplify $F$ to
$$ F = N\mu - P_s V_s - P_l V_l + \gamma A.
$$
Then, assuming that the solid nucleus is spherical with a radius $R$, we obtain
$$ F = N\mu - P_s \frac{4\pi}{3}R^3 - P_l \left(V - \frac{4\pi}{3}R^3 \right) + \gamma 4\pi R^2.
$$

To minimize $F$ we now have two options:
1. Take $\rho_s$ as constant and minimize w.r.t. $R$.
2. Take $N_s$ as constant and minimize w.r.t. both $\rho_s$ and $R$.

Note that by keeping $\rho_s$ constant, $N_s$ is no longer a free variable as it is fully determined by $\rho_s$ and $R$. Hence, $R$ is then the only variable you need to minize with respect to in that case.

#### Option 1: Keep $\rho_s$ constant
When we keep $\rho_s$ constant, minimizing $F$ w.r.t. $R$ leads to
$$ P_s - P_l = \frac{2\gamma}{R} + \left(\frac{\partial \gamma}{\partial R}\right)_{\rho_s}.
$$
However, there is some arbitrariness in defining the radius $R$ of the solid nucleus, and thus in defining $\gamma$. One choice is to take the Gibbs dividing surface for which $R$ is chosen such that $\gamma$ is a minimum w.r.t. changes in the location of the dividing surface. In this convention
$$\left.\left(\frac{\partial \gamma}{\partial R}\right)_{\rho_s} \right|_{R=R^*}=0 ,$$
such that you get the Young-Laplace equation
$$ \Delta P \equiv P_s - P_l = \frac{2\gamma}{R^*}.
$$
This equation implies that for $\gamma>0$, we must have that $P_s>P_l$.

#### Option 2: Keep $N_s$ constant
When we keep $N_s$ constant, then minimizing w.r.t. both $R$ and $\rho_s$ leads to
$$ P_s - P_l = \frac{2}{R} \left[ \gamma + \frac{R}{2}\left(\frac{\partial \gamma}{\partial R}\right)_{\rho_s} -\frac{3}{2}\rho_s \left(\frac{\partial \gamma}{\partial \rho_s}\right)_{R}\right],
$$
where we used that $\delta V_s = 4\pi R^2\delta R = - \frac{N_s}{\rho_s^2} \delta\rho_s$ .
Then using the Gibbs dividing surface again, we find that
$$ \Delta P \equiv P_s - P_l = \frac{2}{R^*} \left[ \gamma -\frac{3}{2}\rho_s \left(\frac{\partial \gamma}{\partial \rho_s}\right)_{R}\right] = \frac{2f}{R^*},
$$
where we call $f=\gamma -\frac{3}{2}\rho_s \left(\frac{\partial \gamma}{\partial \rho_s}\right)_{R}$ the surface stress.
Thus, if $\left(\frac{\partial \gamma}{\partial \rho_s}\right)_{R}$ is positive and "large enough", then we can have that $P_s<P_l$.

---

## Configurations from Carlos

I got two configurations from Carlos (IV and V). Both were obtained using MD simulations (GROMACS) of the $NVT$ ensemble of pseudo-hard spheres (PHS):
$$ u_\text{PHS}(r) = \begin{cases}
50 \left( \frac{50}{49}\right)^{49} \epsilon \;\left( \left(\frac{\sigma}{r}\right)^{50} - \left(\frac{\sigma}{r}\right)^{49} \right) + \epsilon &\quad\text{ for } r< \left( \frac{50}{49}\right)\sigma,\\
0 &\quad\text{ else.}
\end{cases}
$$
They use the mislabeling criterion to find a threshold value for $\bar{q}_6$ to label the particles as liquid or solid. They find the threshold $\bar{q}_{6,t}=0.372$, where neighboring particles are determined as all particles within a radial cutoff of $1.33\sigma$.

Underneath some information on the two configurations. For more information see **[Montero de Hijes et al., JCPC 124 (2020)](https://doi.org/10.1021/acs.jpcc.0c00816)**.

| label | $V/\sigma^3$ | $N$ | $\rho\sigma^3$ | $\eta$ | $\langle N_\text{sol}\rangle$ | $\langle N_\text{sol}\rangle/N$ | $\langle \beta P\sigma^3\rangle$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| IV | 49599.9 | 48207 | 0.97192 | 0.50889 | 5604 | 0.116 | 12.739 |
| V | 49599.9 | 48357 | 0.97494 | 0.51048 | 8602 | 0.178 | 12.579 |

In **[Montero de Hijes et al., JCP 153 (2020)](https://doi.org/10.1063/5.0032602)**, they determine the pressure in- and outside the solid nucleus for configurations IV, VII, and VIII. Underneath the information that they gathered on configuration IV.

| label | $\rho_\text{liq}\sigma^3$ | $\beta P_\text{liq}\sigma^3$ | $\rho_\text{sol}\sigma^3$ | $\beta P_\text{sol}\sigma^3$ | $\rho_\text{sol}^\mu\sigma^3$ | $\beta P_\text{sol}^\mu\sigma^3$ | $\beta\Delta P^\mu\sigma^3$ | $2\beta\gamma_s\sigma^3/R_s$ | $R_s/\sigma$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| IV | 0.9619 | 12.7437 | 1.0613 | 12.6046 | 1.067 | 12.8627 | 0.1190 | 0.1164 | 10.791 |

Here $\rho_\text{liq(sol)}\sigma^3$ and $\beta P_\text{liq(sol)}\sigma^3$ are the measured density and pressure of the liquid (solid) phase. Note that indeed $\beta P_\text{sol}\sigma^3 < \beta P_\text{liq}\sigma^3$, which would mean a negative surface tension. To solve this problem, they say that they need to the chemical potential $\mu$ corresponding to $\beta P_\text{liq}\sigma^3$ of a bulk liquid, and then take $\beta P_\text{sol}^\mu\sigma^3$ corresponding to a bulk solid at this same $\mu$. This pressure they call the thermodynamical pressure. Then from the Young-Laplace equation you get
$$
\Delta P = P_\text{sol}^\mu -P_\text{liq} = \frac{2\gamma_s}{R_s},
$$
where $R_s$ is the radius of the solid nucleus corresponding to the surface tension $\gamma_s$.

---

## Equilibrating the configurations for true HS

### Getting rid of overlaps
As the configurations are of PHS, there were many overlaps in them. So, in order to get rid of these overlaps, I shrank the particles using a shrinking factor 1.0001 times the minimum separating distance in the configuration. Then I quickly grew the particles to their original size using Frank's EDMD grow code.

Underneath are the resulting configurations before and after getting rid of the overlaps. Whether a particle is solid or fluid is determined using the $q_6$ dot product, with $d_c=0.7$, $\xi_c=9$, and $r_\text{cut}=1.45\sigma$. This is slightly different from Carlos' mislabeling, but it is similar enough. Note that I use $r_\text{cut}=1.45\sigma$ instead of $r_\text{cut}=1.33\sigma$, as the former better captures the first minimum of the $g(r)$.

| **IV original** ($N_\text{sol}=5693$) | **IV no overlap** ($N_\text{sol}=5706$) |
| -------- | -------- |
|  |  |
| **V original** ($N_\text{sol}=8473$) | **V no overlap** ($N_\text{sol}=8468$) |
|  |  |

### Checking the stability of the cluster

Since Carlos used PHS and not true HS, we need to check if the configurations are stable, or if we need to slightly adjust the packing fraction. To this end, I used both MC and EDMD simulations of the $NVT$ ensemble initialized at a range of different packing fractions and measured the cluster size as a function of time. Note that I started out with MC instead of EDMD, as I already had an MC code lying around for this and my EDMD needed some more debugging.

#### MC
Underneath, you can see for cluster IV (left) and V (right) the cluster size over time for some of the packing fractions. The horizontal dashed line in all figures indicates the cluster size of the starting configuration, i.e. 5706 for IV and 8468 for V, and the solid black line indicates the simulation at the original packing fraction. To give an idea for the duration of these simulations, the long-time diffusion in the fluid of these systems is approximately $1.1\cdot10^4$ MC steps. So the simulations ran for approximately 250 $\tau_d$.

| **IV:** $\eta\in[0.5088,0.5080]$ | **V:** $\eta\in[0.5104,0.5096]$ |
| -------- | -------- |
|  |  |
| **IV:** $\eta\in[0.5078,0.5090]$ | **V:** $\eta\in[0.5094,0.5086]$ |
|  |  |

#### EDMD
Similar as for MD, we show for EDMD the cluster size for different packing fractions. The time is given in terms of long-time diffusion time $\tau_d$. Note that we ran these simulations for almost 20 times longer than the MC simulations. The reason for this is partly because EDMD is a lot faster than MC, but mostly because I did not have good feeling yet for the MSD in the EDMD simulations when I started these simulations. One more side note: I was a bit lazy and used $\tau_d$ determined for $\eta=0.50889$ and $\eta=0.51048$ to scale the time axis of the runs at other packing fractions. But, since the packing fractions are so close, this shouldn't result in a big shift.

For the IV clusters, the simulations at $\eta=0.5070$ and $\eta=0.5071$ actaully fully melt. For the V cluster, the the simulation at the original packing fraction and those at $\eta=0.5103$ and $\eta=0.5102$ change from a spherical cluster to a system spanning cylinder at $t\approx3200\tau_d$.

| **IV:** $\eta\in[0.5088,0.5080]$ | **V:** $\eta\in[0.5104,0.5096]$ |
| -------- | -------- |
|  |  |
| **IV:** $\eta\in[0.5078,0.5090]$ | **V:** $\eta\in[0.5094,0.5086]$ |
|  |  |

#### Average cluster size
Underneath, I show the average cluster size as a function of packing fraction for the IV and V cluster in the MC and EDMD simulations. The dashed lines indicate linear fits through the data point (omitting the data points of clusters that melt or become system spanning).

Thus, if we do not care that much about the slight growth of the cluster, then we can use the original configurations, i.e those at $\eta=0.50889$ and $\eta=0.51048$. However, for V we already saw with the EDMD that this cluster then (eventually) becomes a system spanning cylinder. So, if we want the cluster to stay approximately the same size as in the original configuration, then we need to use the configurations at approximately $\eta=0.5074$ and $\eta=0.5089$ for the IV and V cluster, respectively.

| IV | V |
| -------- | -------- |
|  |  |

## Local density profile
For both clusters, I calculated the local density as a function of the radial distance w.r.t. the center-of-mass of the solid cluster. For this I used the EDMD simulation at $\eta=0.5074$ for IV and at $\eta=0.5089$ for V. Underneath, the resulting local density. The dots are the data points obtained using a bin size of $0.5\sigma$, the solid line is a tanh fit through the data points, and the horizontal dashed lines indicate the obtained bulk densities for the solid and fluid phase. The table below the figure gives these bulk densities together with the corresponding pressures (obtained from the Carnahan-Starling and Speedy's equations). Note that both $\rho_\text{sol}\sigma^3$ and $\rho_\text{liq}\sigma^3$ are greater than the HS freezing and melting densities, i.e. $\rho_F\sigma^3=0.940$ and $\rho_M\sigma^3=1.037$, and that the pressure of the liquid phase is higher than the pressure of the solid phase.

| IV ($\eta=0.5074$) | V ($\eta=0.5089$) |
| -------- | -------- |
|  |  |

| Cluster | $\rho_\text{sol}\sigma^3$ | $\rho_\text{liq}\sigma^3$ | $\beta P_\text{sol}\sigma^3$ | $\beta P_\text{liq}\sigma^3$ | $\beta \mu_\text{sol}$ | $\beta \mu_\text{liq}$ |
| ------- | -------- | ------- | ------ | -------- |-------- |-------- |
| IV | 1.058 | 0.958 | 12.43 | 12.58 | 17.14| 16.89 |
| V | 1.056 | 0.956 | 12.33 | 12.46 | 17.02 | 16.81 |

> [name=Frank] I added chemical potentials (based on CS and Speedy) because I was curious. Might want to double-check and fix any errors due to using rounded densities.

Some notes:
1. I calculated the solid bulk density counting all particles within a radial distance of $7.5\sigma$ for IV and $8.5\sigma$ for V, and the fluid bulk density counting all particles outside a radial distance of $15\sigma$ for IV and $16\sigma$ for V. This ensures that the solid bulk density is not affected by the fluctuations in the local density of the nucleus caused by the crystal structure.
2. In order to calculate the local density for $r>L/2$, where $L$ is the box size, I just corrected the volume of the spherical shell. This allows you to calculate the local density up to $r=1/2\sqrt{3}L$. In the above situation the box had sides $L=36.7778\sigma$ for IV and $L=36.7797\sigma$ for V.

### Local density for different global densities
Underneath are more results on the local density for different global densities. We clearly see from the density profiles that the cluster size increases for increasing global density.

| IV | V |
| -------- | -------- |
|  |  |

Furthermore, we plot the liquid and solid "bulk" densities as a function of the global packing fraction. We see that for increasing global density, the bulk densities of the liquid and solid decrease, albeit slowly. They are both still significantly higher than the freezing and melting denisities, i.e. $\rho_F\sigma^3=0.940$ and $\rho_M\sigma^3=1.037$. But at least the two different cluster configurations follow the same trend.

| liquid bulk density | solid bulk density |
| -------- | -------- |
|  |  |

### Cluster size fluctuations and center-of-mass position
Lastly, if we want to compute the pressure inside the solid nucleus, then it is useful to get a better idea of the fluctuations in the cluster size and the drift of the center-of-mass of the cluster. Underneath, I show the cluster size and center-of-mass location during the simulation of IV at $\eta=0.5074$ and V at $\eta=0.5989$. For IV, the cluster size fluctuates around an average of 5790 particles (indicated by the horizontal dashed line) with a standard deviation of $\sim$ 440 particles. Making the approximation that the radius is given by $R=\left( \frac{3 N}{4\pi\rho_\text{sol}} \right)^{1/3}$, we get that the cluster has a radius of $13.2\pm4 \;\sigma$. So these fluctuations will probably not form a problem for measuring the pressure inside the cluster. The shifting center-of-mass of the the cluster on the other hand is more problematic and can not be ignored.

| IV: Cluster size | IV: Center-of-mass location |
| -------- | -------- |
|  |  |
| **V: Cluster size** | **V: Center-of-mass location** |
|  |  |

---

## Local pressure profile

I use EDMD simulations to measure the local pressure profiles. Within these simulations you can get the pressure from the momentum transfer during collisions. The global pressure can be calculated with
$$ \beta P/\rho = 1 +
\frac{\sum_\text{collisions} \mathbf{p}_{ij} \cdot \mathbf{r}_\text{ij}}{ 3N \Delta t},
$$
where $\rho=N/V$ is the global density, $\sum_\text{collisions}$ indicates the sum over all collisions during a time interval $\Delta t$, $\mathbf{r}_{ij}=\mathbf{r}_{j}-\mathbf{r}_{i}$ is the center-to-center distance at collision, and $\mathbf{p}_{ij}$ is the momentum of transfer of a collision between particle $i$ and $j$. Note that for a monodisperse HS system with diameter $\sigma$ and mass $m$, we simply have $\mathbf{p}_{ij} \cdot \mathbf{r}_\text{ij} = -m \mathbf{v}_{ij} \cdot \mathbf{r}_\text{ij}$, where $\mathbf{v}_{ij}=\mathbf{v}_{j}-\mathbf{v}_{i}$ (before collision).

To meassure the local pressure profile, I divide the system into spherical shells around the center-of-mass of the nucleus of thickness $\Delta_\text{bin}$. Instead of the total impulse transfer, one then needs to keep track of the sum of the impulse transfer inside each shell, as well as the local density inside each shell. Note though that collisions can happen at the boundary between two shells. In order to handle this correctly, one should add $\mathbf{p}_{ij} \cdot \mathbf{r}_\text{ij}/2$ to the shell in which $\mathbf{r}_i$ lies and $\mathbf{p}_{ij} \cdot \mathbf{r}_\text{ij}/2$ to the shell in which $\mathbf{r}_j$ lies.

I measured the momentum transfer during a simulation of $10^4\tau$, which corresponds to approximately $10^3\tau_d$, for a radial bin width of $\Delta_\text{bin}=0.5\sigma$. Since the center of mass of the nucleus slowly shifts during the simulation, I update it each $0.5\tau$. During this update, I also measure the number of particles in each shell. At the end of the simulation I calculate the local pressure profile using the average number of particles per shell and the total momenta transfer per shell. The resulting density and pressure profiles for IV at $\eta=0.5074$ and V at $\eta=0.5089$ are shown below. Notice that we obtained the profiles for $r<L/2$, where $L$ is the box size.

| IV: density profile | V: density profile |
| -------- | -------- |
|  |  |
| **IV: pressure profile** | **V: pressure profile** |
|  |  |

The horizontal dashed lines indicate the obtained densities and pressures in the "bulk" phases. The table below gives the values of these bulk properties. Additionally, the table gives the globally averaged pressure $\langle\beta P\sigma^3\rangle$, as well as the "bulk" pressures calculated from the "bulk" densities using Speedy's and the Carnahan-Starling equation. Note that $\langle\beta P\sigma^3\rangle = \beta P_\text{liq}\sigma^3$, and that $\beta P_\text{sol}^\text{Speedy}\sigma^3$ and $\beta P_\text{liq}^\text{CS}\sigma^3$ agree well with the measured pressures. Furthermore, note that the pressure inside the solid nucleus is indeed smaller than the pressure of the surrounding fluid.

| Cluster | $\rho_\text{sol}\sigma^3$ | $\rho_\text{liq}\sigma^3$ | $\beta P_\text{sol}\sigma^3$ | $\beta P_\text{liq}\sigma^3$ | $\langle\beta P\sigma^3\rangle$ | $\beta P_\text{sol}^\text{Speedy}\sigma^3$ | $\beta P_\text{liq}^\text{CS}\sigma^3$ |
| ------- | -------- | ------- | ------ |------ | ------ | ------ | ------ |
| IV | 1.0583 | 0.9582 | 12.483 | 12.617 | 12.617 | 12.483 | 12.604 |
| V | 1.0555 | 0.9558 | 12.366 | 12.474 | 12.475 | 12.315 | 12.462 |

Some notes:
- For IV, the bulk properties of the solid phase where obtained within a sphere of radius $7\sigma$, and the bulk properties of the fluid phase within a spherical shell of $14\sigma < r < 18\sigma$.
- For V, the bulk properties of the solid phase where obtained within a sphere of radius $9\sigma$, and the bulk properties of the fluid phase within a spherical shell of $16\sigma < r < 18\sigma$.

---

## Critical cluster?
To check if the clusters are indeed critical clusters in the $NPT$ ensemble, I started 30 MC simulations in the $NPT$ ensemble. To make the volume changes more efficient, I use the algorithm introduced in **[Almarza, JCP 130 (2009)](https://doi.org/10.1063/1.3133328)**.

At first, I started these simulations before I did the (long) simulations for obtaining the pressure profiles reported above. Hence, I ran these simulations at $\beta P\sigma^3=12.593$ for IV and $\beta P\sigma^3=12.477$ for V. For IV approximately 2/3 of the simulations grow and 1/3 melt, and for V approximately 1/3 grow and 2/3 melt. According to Willem the transition for such large clusters is very sharp, so these ratios are close enough to 50/50 grow/melt that we can indeed state that they are critical clusters.

However, just to be sure, I also ran simulations using the "new" pressures reported in the previous section, i.e. $\beta P\sigma^3=12.617$ for IV and $\beta P\sigma^3=12.475$ for V. Note that, for IV, this is a slightly higher pressure than the one I used for the first set of runs. Consequently, we should observe an increase in the number of runs that grow (just as we need). For V, the pressure is slightly lower than the pressure used for the first set of runs, which should increase the number of runs that melt (just as we need). For V, we now observe a perfect 50/50 grow/melt. For IV, approximately 5/6 grow and 1/6 melt.

The table underneath summarizes the results.

| Cluster | $\beta P\sigma^3$ | grow | melt |
| ---- | ---- | ---- | ---- |
| IV | 12.593 | 7/30 | 23/30 |
| IV | 12.617 | 26/30 | 4/30 |
| V | 12.477 | 21/30 | 9/30 |
| V | 12.475 | 15/30 | 15/30 |

---

## Measuring the surface stress

Following Mullins, J. Chem. Phys. **81**, 1436 (1984), the pressure difference between the inside and outside of a spherical crystal cluster is determined by the surface stress $f$, rather than the surface tension $\gamma$. In particular:
$$P_X - P_f = 2f/R,$$
with
$$f = \gamma + A \left(\frac{d \gamma}{d A}\right)_{N_X}$$
the surface stress. The derivative on the right hand side is taken at a fixed number of unit cells in the crystal cluster, i.e. the cluster is compressed without changing the number of particles in the cluster. (Technically, some change in defect concentration as a result of compression could change $N_X$ in the treatment by Mullins, here I'm simplifying to ignore that.)

To measure the surface stress, we perform EDMD simulations of hard spheres in a long box containing a direct coexistence. The long axis of the box and the normal of the interfaces are along the $z$-axis. In the initial configuration, we vary the density of the initial crystal in order to tune the lattice spacing in the $xy$-plane. We then look for the spot where the crystal phase is unstressed, i.e. where the pressure tensor inside the crystal phase is expected to be a constant times the identity matrix. This happens when the pressure in the $z$-direction (blue line below) coincides with the bulk unstressed crystal pressure at the same density (dashed line below, using Speedy). For this coexistence, the only contribution to the global pressure tensor that is anisotropic comes from the interface.

We can then directly calculate $f = L_z (P^z - P^{xy}) / 2$.

| Hexagonal plane | Square plane |
| -------- | -------- |
|  |  |
| $\beta f \sigma^2 \simeq -1.0$| $\beta f \sigma^2 \simeq -0.15$|
| $\beta \Delta P \sigma^3 (R = 10 \sigma) = -0.2$ | $\beta \Delta P \sigma^3 (R = 10 \sigma) = -0.03$

The obtained values are in decent agreement with the measurements by Davidchack and Laird (J. Chem. Phys. **108**, 9452 (1998)) who used a similar approach. (They obtained -0.7 and -0.17, respectively, with significant error bars. Note that these calculations may be sensitive to small details (e.g. Speedy accuracy for the chosen system size)).

Also of interest is the observation of Cacciuto, Auer, and Frenkel (PRL **93**, 166105 (2004)) who looked at nuclei in system of hard spheres with short-ranged depletion attractions. They highlighted the role of surface stress in determining the pressure difference (noting that it can change the sign) and point out the relationship between the surface tension, radius of the surface of tension, and the pressure difference between the fluid and crystal phase at the chemical potential of the fluid. For this, they cite Mullins (see above).

---
### Quick summary of Mullins

$$\Omega = \Omega^X + \Omega^F + \Omega^I$$

Consider two variations. First, vary the number of particles in the cluster without straining/deforming it. To first order in the strain, this leads to:
$$P^X(\mu) - P^F(\mu) = 2 \gamma / R,$$
with $\gamma$ measured at the surface of tension (i.e. such that its derivative with respect to the dividing surface position is 0). Note that the crystal pressure is measured

The second variation is a change in volume of the crystal without changing the number of unit cells. This leads to:
$$P^X - P^F = 2f / R$$

A rederivation of this is also given in this paper by Cacciuto and Frenkel:
https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.93.166105