# Problem with stress of spherical surfaces

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

---

## Calculating the stress (the old way)

### Theory
To obtain the stress, we need to decompose the pressure into the normal and tangential pressure, i.e.
$$
\bar{P}(r) = \frac{1}{3} P_N(r) + \frac{2}{3} P_T(R),
$$
where $\bar{P}$ is the total pressure, $P_N$ is the normal pressure, and $P_T$ is the tangential pressure.
The total pressure (in a specific radial bin) is obtained via
$$
\beta \bar{P}(r) = \rho(r) \left( 1 - \frac{\sum_\text{col} \mathbf{v}_{ij} \cdot \mathbf{r}_{ij} }{ 3 N(r) \Delta t} \right),
$$
where $\rho(r)$ is the number density of the bin, $N(r)$ is the number of particles in the bin, $\Delta t$ is the time interval, $\mathbf{v}_{ij} = \mathbf{v}_{j}-\mathbf{v}_{i}$ (before collosion), and $\mathbf{r}_{ij} = \mathbf{r}_{j}-\mathbf{r}_{i}$.
We can decompose the latter into $\mathbf{v}_{ij} = \mathbf{v}_{ij}^\perp + \mathbf{v}_{ij}^\parallel$ and $\mathbf{r}_{ij} = \mathbf{r}_{ij}^\perp + \mathbf{r}_{ij}^\parallel$.
Using this, we can then write
$$
\beta P_N(r) = \rho(r) \left( 1 - \frac{\sum_\text{col} v_{ij}^\perp r_{ij}^\perp }{ N(r) \Delta t} \right),
$$
and
$$
\beta P_T(r) = \rho(r) \left( 1 - \frac{\sum_\text{col} \mathbf{v}_{ij}^\parallel \cdot \mathbf{r}_{ij}^\parallel }{ 2 N(r) \Delta t} \right) = \rho(r) \left( 1 - \frac{\sum_\text{col} ( \mathbf{v}_{ij} \cdot \mathbf{r}_{ij} - v_{ij}^\perp r_{ij}^\perp ) }{ 2 N(r) \Delta t} \right).
$$

### Example pressure profiles
Using this we can obtain the profiles, like the two examples shown below. What is strange though, is that, for small $r$, $P_N >> P_T$. Even though in bulk they both should be equal to $\bar{P}$. Especially for the profile on the right, we see that $P_N$ is significantly larger than $P_T$ even for $r/\sigma\approx10$.

| $\rho_f\sigma^3=0.960$ | $\rho_f\sigma^3=0.947$ |
| -------- | -------- |
|  |  |

### Finding out why it is wrong
So to check this behavior, I looked at a system which is fully fluid and calculated the radial profiles from the center of the box. In this situation, $P_N(r)$ and $P_T(r)$ should both be equal to $\bar{P}(r)$ for all $r$. This means that $\frac{v_{ij}^\perp r_{ij}^\perp}{\mathbf{v}_{ij} \cdot \mathbf{r}_{ij}}$ should be equal to $1/3$. However, as shown below, we see that for small $r$, this is not the case. Looking at $r_{ij}^\perp$ and $v_{ij}^\perp/|v_{ij}|$ separately (note that $|r_{ij}|=\sigma$), we also see that even though they both correctly go to $1/2$ for large $r$, for smalll $r$ they are larger than $1/2$. So apparently the normal direction of $\mathbf{v}_{ij}$ and $\mathbf{r}_{ij}$ is favored over the tangential direction for small $r$.

To check that this is a real feature and not just a mistake in my code, I ran the same code, but always took the $x$-direction as the normal direction. This resulted in the correct behavior for all $r$, see below. (Just ignore the very first bin as the statistics of that one are just very low).

Next, I set out to find the scaling with which the normal direction is facored over the tangential direction for $\mathbf{v}_{ij}$ and $\mathbf{r}_{ij}$. The loglog-plot underneath shows that all follow $1/r^2$ for large $r$.

Yet, even though we know that they go as $1/r^2$ for large $r$, it is difficult to find the full behavior. The main problem for this is that the trend does not diverge in the limit of $r\to 0$. Hence, next we focus just on $v_{ij}^\perp r_{ij}^\perp$, and instead of dividing by $\mathbf{v}_{ij} \cdot \mathbf{r}_{ij}$, we now divide by $\mathbf{v}_{ij}^\parallel \cdot \mathbf{r}_{ij}^\parallel$. This ensures that the trend diverges for small $r$. As shown underneath (blue line), we see that $\frac{v_{ij}^\perp r_{ij}^\perp}{\mathbf{v}_{ij}^\parallel \cdot \mathbf{r}_{ij}^\parallel}$ can be described by $\frac{1}{4\pi r^2}$. And doing some algebra we can then derive that $\frac{v_{ij}^\perp r_{ij}^\perp}{\mathbf{v}_{ij} \cdot \mathbf{r}_{ij}}$ is described by $\frac{1+2\pi r^2}{1+6\pi r^2}$.

**Note 1:** For $r<\sigma$ the blue line differs quite a bit from the black dashed line. However, I think that the reason for this is just the weirdness of binning for small $r$. What I mean with this is that, when two particles collide, I assign half of the momentum transfer to the bin of particle $i$ and half to the bin of particl $j$. For $r<\sigma$ this way of assigning to bins can become a bit weird.
**Note 2:** I have to admit that I just guessed the $1/4\pi r^2$ behavior. I still can't explain why the blue line has to follow this trend.


Underneath I show, for the same nuclei as shown at the start, the corrected pressure profiles. These profiles are obtained by correcting the momentum transfer in the normal direction with the term $\frac{1}{3} \frac{1+6\pi r^2}{1+2\pi r^2}$. We see that in the bulk crystal and bulk fluid, $P_N=P_T=\bar{P}$.

| $\rho_f\sigma^3=0.960$ | $\rho_f\sigma^3=0.947$ |
| -------- | -------- |
|  |  |

With these corrected pressure profiles, I actually don't think that we need to normalize the stress profile before integrating it.
Integrating the left (non normalized) profile over $r/\sigma\in[2.4,13.6]$ we obtain the stress $\beta f\sigma^2=-0.70$. And integrating the right profile over $r/\sigma\in[10.0,31.1]$ we obtain the stress $\beta f\sigma^2=-0.69$.

**Update 21 Dec**
The problem is caused by the fact that the vector $\mathbf{r}_{ij}$ is not infinitesimal. To illustrate this, see the figures below. If you take $\mathbf{r}_{ij}$ to be infinitesimal and assign the radial part of $\mathbf{r}_{ij}$ to the bin corresponding to the center of collision, then you get the figure on the left. However, if you take $\mathbf{r}_{ij}$ to have a real physical length and you assign the radial part of $\mathbf{r}_{ij}$ to the bins corresponding to the first and second particle, then you get the figure on the right.
Note: these figures where made by drawing $10^9$ random pairs of particles. The yellow line is just a sanity check and indicates the part of $\mathbf{r}_{ij}$ in the x-direction.

| Infinitesimal $\mathbf{r}_{ij}$ | Real $\mathbf{r}_{ij}$ |
| -------- | -------- |
|  |  |

The trend in the right figure corresponds with the things seen before (see figure below).


So a bit ignorant you could decide to solve our problem by assigning the impulse transfers to the bin of the center of collission instead of the bin over particles $i$ and $j$. However, I don't think that this would go well, as you need the density of each bin to calculate the pressures in the end and depending on the binsize, this can fluctuate quite a bit for the first bins. So the better thing to do it to find the analytical expression for the bias towards normal directions.

In **[Nakamura et al. (2015)](https://www.sciencedirect.com/science/article/pii/S0010465514004007#br000075)** they address this problem and give the solution (Eqs. 18-20), which they actually obtained from **[Nakamura et al. (2011)](https://pubs.aip.org/aip/jcp/article/135/9/094106/983507)**. The sad news is that you cannot correct for this bias in the post processing. You really need to do it for each individual collision as the correction depends also on the specific $\mathbf{v}_{ij}$ and $\mathbf{r}_{ij}$ of the collision.

---

## Calculating the stress (the correct way)

### Theory / methods
Based on the derivation in **[Nakamura et al. (2015)](https://www.sciencedirect.com/science/article/pii/S0010465514004007#br000075)** and **[Nakamura et al. (2011)](https://pubs.aip.org/aip/jcp/article/135/9/094106/983507)**, I changed my code to correct for the "normal bias" of each collision. Here is a brief summary of the theory and method.

For a collision between particle $i$ and $j$, we have have the properties
$\mathbf{r}_{ij}=\mathbf{r}_j-\mathbf{r}_i$ and $\mathbf{v}_{ij}=\mathbf{v}_j-\mathbf{v}_i$. Note that, since we consider monodisperse hard spheres, $|\mathbf{r}_{ij}|/\sigma=1$. Then, instead of assigning the momentum transfer (i.e. $-\mathbf{v}_{ij}\cdot\mathbf{r}_{ij}$) to the bin of particle $i$ and $j$, we use the Irving-Kirkwood contour. This contour is a straight path from $\mathbf{r}_i$ to $\mathbf{r}_j$ given by
$$
\mathbf{\ell}(\lambda) = \mathbf{r}_i + \lambda \mathbf{r}_{ij},
$$
where $0\leq\lambda\leq1$. Note that for all calculation we take the reference frame of the nucleus. That means that $\mathbf{r}_i$ points from the center-of-mass of the nucleus to particle $i$.

In order to get a better grip on the large effect of decomposing at different points along the collision contour, I show three cases underneath. In all three cases, the vector $\mathbf{r}_{ij}$ is decomposed into its normal (red) and tangential direction (blue) parts. For the left case, this is done at $\lambda=0$. For the middle case, this is done at $\lambda=1/2$. And for the right case, this is done at $\lambda=1$. The large differences in the ratio of the normal and tangential component are clearly visible.


We thus assign a fraction of the momentum transfer to each bin along the contour according to the specific path $\mathbf{\ell}_a$ to $\mathbf{\ell}_b$ inside each bin. Here $\mathbf{\ell}_a\equiv\mathbf{\ell}(\lambda_a)$ and $\mathbf{\ell}_b\equiv\mathbf{\ell}(\lambda_b)$ indicate either an intersection with a bin boundary or a terminal point of the contour (i.e. $\lambda_a=0$ or $\lambda_b=1$). The contribution to the momentum transfer in the normal and tangential direction of a specific bin are then given by
$$
\begin{align}
\Pi_{ij}^N &= -\mathbf{v}_{ij}\cdot(\mathbf{\ell}_b-\mathbf{\ell}_a) - \left[ G_{ij}(\mathbf{v}_{ij},\mathbf{\ell}_b) - G_{ij}(\mathbf{v}_{ij},\mathbf{\ell}_a) \right] , \tag{1}\\
\Pi_{ij}^T &= \frac{1}{2} \left[ G_{ij}(\mathbf{v}_{ij},\mathbf{\ell}_b) - G_{ij}(\mathbf{v}_{ij},\mathbf{\ell}_a) \right] , \tag{2} \\
G_{ij}(\mathbf{v}_{ij},\mathbf{\ell}_a) &= \frac{-(\mathbf{v}_{ij}\cdot\mathbf{r}_{ij}) |\mathbf{r}_{ij}\times\mathbf{r}_i| }{|\mathbf{r}_{ij}|^2} \arctan\left(\frac{\mathbf{r}_{ij}\cdot\mathbf{\ell}_a}{|\mathbf{r}_{ij}\times\mathbf{r}_i|}\right) \\
&+ \frac{-(\mathbf{v}_{ij}\times\mathbf{r}_{ij}) \cdot (\mathbf{r}_{ij}\times\mathbf{r}_i ) }{|\mathbf{r}_{ij}|^2} \;\ln|\mathbf{\ell}_a|, \tag{3}
\end{align}
$$
These equations are taken from Eq.(18)-(20) in **[Nakamura et al. (2015)](https://www.sciencedirect.com/science/article/pii/S0010465514004007#br000075)**. Note that in this paper Eq. (20) has a typo of a cross-product instead of a dot-product in the second term. I will not derive these equations here. For a derivation, see Appendix B of **[Nakamura et al. (2011)](https://pubs.aip.org/aip/jcp/article/135/9/094106/983507)**. Note that both papers use the force $\mathbf{F}$ which I have replaced with $-\mathbf{v}$, and that their definition of particles $i$ and $j$ is the other way around.

We can write these expressions in a simpler way by introducing the variables
$$
\begin{align}
\mathbf{\omega} &= \mathbf{r}_{ij}\times\mathbf{r}_i, \\
\alpha(\lambda) &= \mathbf{r}_{i}\cdot\mathbf{\ell}(\lambda).
\end{align}
$$
Note that in **[Nakamura et al. (2011)](https://pubs.aip.org/aip/jcp/article/135/9/094106/983507)** the variable $\alpha$ is called $\sigma$, but I wanted to avoid confusion with the particle diameter $\sigma$ here. In my code it is also called sigma though.
Equations (1)-(3) can then be rewritten as
$$
\begin{align}
\Pi_{ij}^N &= -\mathbf{v}_{ij}\cdot\mathbf{r}_{ij}\frac{\alpha_b-\alpha_a}{|\mathbf{r}_{ij}|^2} - \left[ G_{ij}(\mathbf{v}_{ij}) \right]^b_a , \tag{4}\\
\Pi_{ij}^T &= \frac{1}{2} \left[ G_{ij}(\mathbf{v}_{ij}) \right]^b_a , \tag{5} \\
\left[ G_{ij}(\mathbf{v}_{ij}) \right]^b_a &= \frac{-(\mathbf{v}_{ij}\cdot\mathbf{r}_{ij}) |\mathbf{\omega}| }{|\mathbf{r}_{ij}|^2} \left[\arctan\frac{\alpha_b}{|\mathbf{\omega}|} -\arctan\frac{\alpha_a}{|\mathbf{\omega}|} \right] + \frac{-(\mathbf{v}_{ij}\times\mathbf{r}_{ij}) \cdot \mathbf{\omega} }{|\mathbf{r}_{ij}|^2} \; \ln\frac{|\mathbf{\ell}_b|}{|\mathbf{\ell}_a|}, \tag{6}
\end{align}
$$
where $\alpha_a\equiv\alpha(\lambda_a)$ and $\alpha_b\equiv\alpha(\lambda_b)$.

The convenient thing about the above notation is that the crossings with the boundaries of the bin, i.e. $\alpha_a$ and $\alpha_b$, can be very easily obtained by solving
$$
\alpha_a = \pm \sqrt{|\mathbf{r}_{ij}|^2 R_a^2 - |\mathbf{\omega}|^2}, \tag{7}
$$
where $R_a$ indicates the radius of the boundary of the bin that is crossed.
To obtain the relation of Eq. (7), we use the vector identity $|\mathbf{a}\times\mathbf{b}|=a^2 b^2 - (\mathbf{a}\cdot\mathbf{b})^2$ such that we can write
$$
|\mathbf{\omega}|^2 = |\mathbf{r}_{ij}|^2 |\mathbf{\ell}(\lambda)|^2 - \alpha(\lambda)^2.
$$
Solving this quadratic equation for $\alpha$ gives you Eq. (7).
The sign ($+$ or $-$) of Eq. (7) depends on the type of collision path. There are 3 types or cases for the collision:


- **Type 1**: $\lambda_0<0$
or equivalently $\mathbf{r}_{ij}\cdot\mathbf{r}_i>0$
Thus: all $\alpha$ positive
- **Type 2**: $\lambda_0>1$
or equivalently $\mathbf{r}_{ij}\cdot\mathbf{r}_i < -|\mathbf{r}_{ij}|^2$
Thus: all $\alpha$ negative
- **Type 3**: $0\leq\lambda_0\leq1$
or equivalently $- |\mathbf{r}_{ij}|^2 \leq \mathbf{r}_{ij}\cdot\mathbf{r}_i \leq 0$
Thus: $\alpha$ positive for $\lambda<\lambda_0$ and $\alpha$ negative for $\lambda>\lambda_0$ and

Thus, by making a list of all crossings $\alpha(\lambda)$ for a collision, we can easily calculate the normal and tangential contribution to each bin along the path between particles $i$ and $j$ with Eqs. (4)-(6). **[Nakamura et al. (2011)](https://pubs.aip.org/aip/jcp/article/135/9/094106/983507)** even has a psuedo code for finding these $\alpha$'s in Appendix D. But becareful because it contains some mistakes in the for-loop boundaries.

### Checks
To check that everything works (as it is a bit more sensitive to bugs then the original method), I first checked my code on a simulation containing just fluid ($\eta=0.45$). Below I compare the results for the old method and new method. The old method assigns half of the impulse transfer to bin of particle $i$ and half to bin of particle $j$ and has no post correction.
In the first row, I show the fraction of the impulse tranfer in the normal direction in comparision to the total impulse transfer. In an isotropic fluid, this should be 1/3 regardless of the radial distance. We clearly see that the old method has this bias for the normal direction for small $r$, but that the new method correctly follows 1/3 for the entire range of $r$.
In the second row, I show the total, normal, and tangential pressure as a function of $r$. Again, we see that the old method has a bias towards the normal direction for small $r$ (note that the total pressure is fine though!), and that the new method correctly calculates the different pressures.

| Old method | New method |
| -------- | -------- |
|  |  |
|  |  |

To also see the difference on an actual system with a crystal nucleus, underneath I show for one of the small systems the difference between the new and the old method. Interestingly, the normal pressure seems to be less sensitive to the density fluctuations due to the crystal lattice. Unfortunatetly, this makes chosing the integration range for obtaining the surface stress a bit more difficult.
If I integrate the profile below over $r/\sigma\in[3.5,13.0]$, I obtain the stress $\beta f\sigma^2=-0.70$.

| Old method | New method |
| -------- | -------- |
|  |  |

The above system contains only a small nucleus of around 4200 particles. Some examples of the profiles for larger nuclei are shown below. Integrating the left profile over $r/\sigma\in[5.5,16.5]$ results in a stress of $\beta f\sigma^2=-0.71$. Integrating the right profile over $r/\sigma\in[7.2,19.9]$ results in a stress of $\beta f\sigma^2=-0.69$.

| $N_x\approx10.8\cdot10^3$ particles | $N_x\approx17.5\cdot10^3$ particles |
| -------- | -------- |
|  |  |

### New results
I reran all simulations for all nuclei. Below are the new results. Both column show the same data points, the only difference is in the model (black dashed line). On the right, the model is trained on the new data (only Frank's nuclei). On the left, the model is trained on the old data of the simulation in which we did not decompose the pressure (also only Frank's nuclei). The differences in the models are very small and (almost) invisible. Only the stress shows a visible difference in the model.

In general, the results agree extremely well with the model. And, even though there is a bit of a spread in the newly obtained surface stresses, they generally agree reasonably well with the theoretical trend. Keep in mind that the surface stresses are obtained by integrating over a stress profile. Hence, the error in the surface stresses are easily around $0.01-0.02 \; k_BT\sigma^{-2}$.

From the model, we obtain that $\beta\gamma_0\sigma^3 = 0.5488$.
(For the model trained on the old data this was $\beta\gamma_0\sigma^3 = 0.5483$.)

| Model on old data | Model on this data |
| -------- | -------- |
|  |  |

---