# Convolutional Bond Orientational Order Parameter (cBOP)

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

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

Alptug, Laura

---

## To do
- [x] replace deltas into isotropic ($\beta = 0$) Fisher- Bingham (FB) distributions
- [x] variances are proportional to solid angles
- [x] follow the usual local BOP procedure
- [ ] calculate Dellago and Boattini averages
- [ ] Simulate various systems, take snapshots at equilibrium
- [ ] calculate local and averaged (convolutional) BOPs
- [ ] Plot distributions
- [ ] plot q4-q6 plane
- [ ] Plot Jensen–Shannon divergence for all
- [ ] Celebrate if our JS-div is lower :smile:
- [ ] if successful, repeat for w's

---
---
## BOP extension
We start by the complex conjugate (cc) of BOPs whose reason will be obvious later on. In addition, rotational inv. version is also invariant under cc:
$$
\begin{aligned}
q_{lm}(i) ={}& \frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b(i)} Y^{m*}_{l}\left( \theta_{ij},\phi_{ij} \right)\\
={}&
\frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b(i)}
\int_0^{2\pi}\mathrm{d}\phi\int_0^{\pi}\mathrm{d}\theta
Y^{m*}_{l}\left( \theta,\phi \right) \delta(\theta-\theta_{ij})\delta(\phi-\phi_{ij})\\
={}&
\frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b(i)}
\int_0^{2\pi}\mathrm{d}\phi\int_0^{\pi}\mathrm{d}\theta \sin(\theta)
Y^{m*}_{l}\left( \theta,\phi \right) \frac{1}{\sin(\theta)}\delta(\theta-\theta_{ij})\delta(\phi-\phi_{ij})\\
={}&
\frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b(i)}
\int_{S^2}\mathrm{d}\Omega
Y^{m*}_{l}\left( \theta,\phi \right) \frac{1}{\sin(\theta)}\delta(\theta-\theta_{ij})\delta(\phi-\phi_{ij}).
\end{aligned}
$$

Here, we see that BOPs assume (projected) point particles, whose pdf is given by:

$$
\begin{aligned}
\rho_i(\theta,\phi) ={}&
\frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b(i)}
\frac{1}{\sin(\theta)}\delta(\theta-\theta_{ij})\delta(\phi-\phi_{ij})
\\={}&
\sum_{l=0}^\infty\sum_{m=-l}^{l}q_{lm}(i)Y^{m}_{l}\left( \theta,\phi \right)
\end{aligned}
$$

To make BOPs robust under fluctuations, we will replace our point particles with "fuzzy" particles. To make the particles fuzzy, we need to convolve the density function with a kernel. Assuming fluctuations in the system follow a Gaussian distribution, we choose the kernel to be an isotropic Gaussian on $S^2$. This distribution is isotropic Fisher-Bingham five-parameter (FB5) distribution (also known as the Kent distribution.)

$$
g(\hat{\mathbf{r}},\kappa_{ij},\hat{\mathbf{r}}_{ij}) = \frac{1}{C(\kappa_{ij})}e^{\kappa_{ij} \hat{\mathbf{r}}\cdot \hat{\mathbf{r}}_{ij}}
$$
where
$$
C(\kappa_{ij}) = 2\pi\frac{\Gamma(1/2)}{\Gamma(1)\sqrt{\kappa_{ij}/2}}I_{1/2}(\kappa_{ij}) = \frac{4 \pi }{\kappa_{ij}}\sinh(\kappa_{ij}).
$$
Making use of SANN definition of NNs, we choose $\kappa_{ij}$ such that
:::info
I have a justification for this choice. Take a particle and 2 NNs,
displace the particles by the same amount
project them down to a sphere centered at our particle
the closer NN will experience more angular change
therfore the variance of the closer should be larger under thermal fluctuations!
:::

$$
\kappa^{-1}_{ij} \propto 2\pi\left(1-\frac{r_{ij}}{R_i}\right)
$$
where
$$
R_i = \frac{1}{N_b(i)-2}\sum_{j\in\mathcal{N}_b(i)}r_{ij}
$$
Notice that $R_i \rightarrow \lambda R_i$ as $r_{ij} \rightarrow \lambda r_{ij}$, therefore, $\kappa_{ij}$ stays invariant under rescaling. This is our main difference from SOAPs. We are not breaking scale invariance and we do not need a radial function.

As a result, we define the new NN density function $\tilde\rho_i(\theta,\phi)$ as

$$
\begin{aligned}
\tilde\rho_i(\theta,\phi) ={}&
\frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b(i)}g(\hat{\mathbf{r}},\kappa_{ij},\hat{\mathbf{r}}_{ij})
\\={}&
\sum_{l=0}^\infty\sum_{m=-l}^{l}\tilde q_{lm}(i)Y^{m}_{l}\left( \theta,\phi \right) ,
\end{aligned}
$$
and in the second line, we have defined cBOPs implicitly. We invert that expression to obtain cBOPs explicitly,

$$
\tilde q_{lm}(i) =
\frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b(i)}
\int_{S^2}\mathrm{d}\Omega
Y^{m*}_{l}\left( \theta,\phi \right) g(\hat{\mathbf{r}},\kappa_{ij},\hat{\mathbf{r}}_{ij}).
$$

This integral is quite non-trivial but luckily, someone else has already calculated the general formula: https://arxiv.org/abs/1501.04395

### Alem et al. simplification
:::warning
This part is a complete copy of Alem et al.
:::
In this section, we derive an analytic expression for the computation of
spherical harmonic coefficients of the FB5 distribution. For convenience, we
determine the spherical harmonic coefficients of the FB5 distribution with
mean~(centre) $\hat{\mu}^0$ located on $z$-axis and major axis
$\hat{\eta}_1^0$ and minor axis $\hat{\eta}_2^0$ aligned along $x$-axis and
$y$-axis respectively, that is,

$$
\hat{\mu}^0=[0~0~1]^T\quad
\hat{\eta}_1^0=[1~0~0]^T\quad \hat{\eta}_2^0=[0~1~0]^T.
$$
With the centre, major and minor axes as given in \eqref{Eq:paramters_standard},
the FB5 distribution, referred to as the \emph{standard} Fisher-Bingham~(FB)
distribution, has the pdf given by
$$
f(\hat{x};\kappa,\beta) := \frac{1} { C(\kappa , \beta)} e^{\kappa\cos\theta +
\beta\sin^2\theta\cos 2\phi},
$$
which is related to the FB5 distribution pdf
$g(\hat{x};\kappa,\beta,\hat{\mu},\mathbf A)$ given in
\eqref{eq:kent_rotated} through the rotation operator
$\mathcal{D}(\varphi,\vartheta,\omega)$ as

$$
g(\hat{x};\kappa,\beta,\hat{\mu},\mathbf A)=\left(\mathcal{D}(\varphi,\vartheta,\omega)f\right) (\hat{x
}; \kappa,\beta) = f(\mathbf{R}^{-1}\hat{x}; \kappa,\beta),
$$
where the rotation matrix is related to $\hat{\mu},~\hat{\eta}_1,~\hat{\eta}_2$~(parameters of FB5 distribution) as
$$
\mathbf{R}=[\hat{\eta}_1 , \hat{\eta}_2, \hat{\mu}],
$$
and the Euler angles $(\varphi,\vartheta,\omega)$ are related to the rotation matrix $\mathbf{R}$ through \eqref{eq:rot_matrix}.

$$
f_{\ell}^{m} = \int f(\hat{x
}; \kappa,\beta) \overline{Y_\ell^m(\hat{x})}
ds(\hat{x}).
$$

$$
\begin{split}
f_{\ell}^{m} &=
\frac{\pi i^{-m}}{C(\kappa,\beta)}\sqrt{\frac{2\ell+1}{2\kappa}}
\sum_{n=0}^{\infty} (2n+1)I_{n+1/2}(\kappa)
\sum_{u=-n}^{n} \left(d_{u,0}^{n}({\pi}/{2})\right)^2
\sum_{t=0}^{\infty}\frac{(\beta/2)^{2t+m/2}}{ \Gamma(t+1)\Gamma(t+m/2+1)}
~\times\\
& \hspace{70mm}\sum_{u'=-\ell}^{\ell} d_{u',0}^{\ell}({\pi}/{2})
d_{u',m}^{\ell}({\pi}/{2})\, G(4t+m+1,u+u'),
\end{split}
$$
where
$$
\begin{aligned}
G(p,q) & := \int_{0}^{\pi} (\sin\theta)^p e^{i q \theta} d\theta, \nonumber \\
& = \frac{\pi e^{iq\pi/2}\Gamma(p+2)}
{2^{p}(p+1)\Gamma(\frac{p+q+2}{2})\Gamma(\frac{p-q+2}{2})}\,.
\end{aligned}
$$
:::warning
end of copy, all the expressions are for $m \geq 0$, negative $m$s are given by $f^{-m}_l = (-1)^m f^{m*}_l$
:::

Now I take the limit $\beta \rightarrow 0$:

$$
\begin{split}
f_{\ell}^{m} &=
\frac{\pi i^{-m}}{C(\kappa)}\sqrt{\frac{2\ell+1}{2\kappa}}
\sum_{n=0}^{\infty} (2n+1)I_{n+1/2}(\kappa)
\sum_{u=-n}^{n} \left(d_{u,0}^{n}({\pi}/{2})\right)^2
\sum_{t=0}^{\infty}
\frac{\delta_{t,-m/4}}{ \Gamma(t+1)\Gamma(t+m/2+1)}
~\times\\
& \hspace{70mm}\sum_{u'=-\ell}^{\ell} d_{u',0}^{\ell}({\pi}/{2})
d_{u',m}^{\ell}({\pi}/{2})\, G(4t+m+1,u+u'),
\end{split}
$$

for $m \geq 0$ and $t \geq 0$, $\delta_{t,-m/4}\rightarrow \delta_{t,0}\delta_{m,0}$, Then, $f_{\ell}^{m}=\delta_{m,0}f_{\ell}^{0}$:

$$
\begin{split}
f_{\ell}^{0} &=
\frac{\pi}{C(\kappa)}\sqrt{\frac{2\ell+1}{2\kappa}}
\sum_{n=0}^{\infty} (2n+1)I_{n+1/2}(\kappa)
\sum_{u=-n}^{n} \left(d_{u,0}^{n}({\pi}/{2})\right)^2
%
~\times\\
& \hspace{70mm}\sum_{u'=-\ell}^{\ell} d_{u',0}^{\ell}({\pi}/{2})
d_{u',m}^{\ell}({\pi}/{2})\, G(1,u+u'),
\end{split}
$$
where
$$
\begin{aligned}
G(1,q) & := \int_{0}^{\pi} \sin\theta e^{i q \theta} d\theta, \nonumber \\
& = \frac{\pi e^{iq\pi/2}}
{2\Gamma(\frac{3+q}{2})\Gamma(\frac{3-q}{2})}\\
& =\frac{2 e^{iq\pi/2}}
{1-q^2}\cos\left(\frac{\pi q}{2}\right)\\
& = \frac{1+e^{iq\pi}}{1-q^2}
\end{aligned}
$$
Then,

$$
\tilde{q}_{lm}(i) =\frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b(i)} D^l_{m,0}(\varphi,\vartheta,\omega)f_l^0(\kappa_{ij})
$$
for
$$
\mathbf{R}(\theta,\phi) =
\begin{pmatrix}
-\sin\phi & \cos\phi\cos\theta & \cos\phi\sin\theta\\
\cos\phi & \sin\phi\cos\theta & \sin\phi\sin\theta\\
0 & -\sin\theta & \cos\theta
\end{pmatrix}
$$
$\vartheta=\theta$, $\varphi = \phi$, $\omega = -\pi/2$. Then,
$$
\begin{aligned}
\tilde{q}_{lm}(i) ={}&\frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b(i)} D^l_{m,0}(\phi,\theta,-\pi/2)f_l^0(\kappa_{ij})\\
={}&\frac{1}{N_b(i)} \sum_{j\in\mathcal{N}_b(i)} e^{-im\phi_{ij}} d^l_{m,0}(\theta_{ij})f_l^0(\kappa_{ij})\\
={}&\frac{1}{N_b(i)}
\sqrt{\frac{(l-m)!}{(l+m)!}}
\sum_{j\in\mathcal{N}_b(i)} e^{-im\phi_{ij}} P^l_m(\cos\theta_{ij})f_l^0(\kappa_{ij})\\
={}&\frac{1}{N_b(i)}
\sqrt{
\frac{4\pi}{2l+1}
}
\sum_{j\in\mathcal{N}_b(i)} Y_l^{m*}(\theta_{ij},\phi_{ij})f_l^0(\kappa_{ij})
\end{aligned}
$$

:::danger
This is computational hell, let's take a shortcut for now
:::

## Rot-inv cBOPs

We define rotationally invariant cBOPs in the usual way:

$$
\tilde{q}_l(i) := \sqrt{\frac{4\pi}{2l+1}\sum_{m=-l}^l |\tilde q_{lm}(i)|^2}
$$

then,

$$
\begin{aligned}
\sum_{m=-l}^l |\tilde q_{lm}(i)|^2 ={}&
\frac{1}{N_b(i)^2} \sum_{j\in\mathcal{N}_b(i)}\sum_{k\in\mathcal{N}_b(i)}
\int_{S^2}\mathrm{d}\Omega\int_{S^2}\mathrm{d}\Omega^\prime
g(\hat{\mathbf{x}},\kappa_{ij},\hat{\mathbf{r}}_{ij})g(\hat{\mathbf{y}},\kappa_{ik},\hat{\mathbf{r}}_{ik})
\sum_{m=-l}^l
Y^{m}_{l}\left( \hat{\mathbf{x}}\right)
Y^{m*}_{l}\left( \hat{\mathbf{y}} \right) \\
={}&
\frac{2l+1}{4\pi}
\frac{1}{N_b(i)^2} \sum_{j\in\mathcal{N}_b(i)}\sum_{k\in\mathcal{N}_b(i)}
\int_{S^2}\mathrm{d}\Omega\int_{S^2}\mathrm{d}\Omega^\prime
g(\hat{\mathbf{x}},\kappa_{ij},\hat{\mathbf{r}}_{ij})g(\hat{\mathbf{y}},\kappa_{ik},\hat{\mathbf{r}}_{ik})P_l(\hat{\mathbf{x}}\cdot\hat{\mathbf{y}})
\end{aligned}
$$
Maybe this solves
$$
\begin{aligned}
\sum_{m=-l}^l |\tilde q_{lm}(i)|^2
={}&
\frac{1}{N_b(i)^2} \sum_{j\in\mathcal{N}_b(i)}\sum_{k\in\mathcal{N}_b(i)}
P_l(\hat{\mathbf{r}}_{ij}\cdot\hat{\mathbf{r}}_{ik})f_l^0(\kappa_{ij})f_l^{0*}(\kappa_{ik})
\end{aligned}
$$

Define

$$
\begin{split}
h_\ell(\kappa) &:= \frac{\sqrt{2\pi\kappa}}{4 \sinh(\kappa)}
\sum_{n=0}^{\infty} (2n+1)I_{n+1/2}(\kappa)
\sum_{u=-n}^{n} \left(d_{u,0}^{n}({\pi}/{2})\right)^2
%
~\times\\
& \hspace{70mm}\sum_{u'=-\ell}^{\ell} d_{u',0}^{\ell}({\pi}/{2})
d_{u',0}^{\ell}({\pi}/{2})\, G(1,u+u'),
\end{split}
$$
then
$$
f_\ell^0(\kappa) = \sqrt{\frac{2\ell+1}{4\pi}} h_\ell(\kappa)
$$.

Finally, we have
$$
\tilde{q}_l(i) = \frac{1}{N_b(i)}\sqrt{
\sum_{j\in\mathcal{N}_b(i)}\sum_{k\in\mathcal{N}_b(i)}
P_l(\hat{\mathbf{r}}_{ij}\cdot\hat{\mathbf{r}}_{ik})h_l(\kappa_{ij})h_l^{*}(\kappa_{ik})
}
$$

Now, I will make a bold statement. Define
$$
F_\ell^n :=
(2n+1)
\sum_{u=-n}^{n} \left(d_{u,0}^{n}({\pi}/{2})\right)^2
\sum_{u'=-\ell}^{\ell}
\left( d_{u',0}^{\ell}({\pi}/{2})\right)^2
G(1,u+u'),
$$
After some Mathematica investigation, I conjecture
$$
F_\ell^n = 2 \delta_{\ell,n}.
$$
Therefore,
$$
\begin{aligned}
h_\ell(\kappa) ={}& \frac{\sqrt{2\pi\kappa}}{2 \sinh(\kappa)}
\sum_{n=0}^{\infty}I_{n+1/2}(\kappa)\\
={}&\frac{\sqrt{2\pi\kappa}}{4 \sinh(\kappa)}e^\kappa\mathrm{Erf}(\sqrt{2\kappa})
\end{aligned}
$$

## Final Results

$$
\begin{aligned}
\tilde{q}_{lm}(i) ={}&
\frac{1}{H(i)}
\sum_{j\in\mathcal{N}_b(i)} Y_l^{m*}(\theta_{ij},\phi_{ij})h(\kappa_{ij})
\end{aligned}
$$

$$
\tilde{q}_l(i) = \frac{1}{H(i)}\sqrt{
\sum_{j\in\mathcal{N}_b(i)}\sum_{k\in\mathcal{N}_b(i)}
P_l(\hat{\mathbf{r}}_{ij}\cdot\hat{\mathbf{r}}_{ik})h(\kappa_{ij})h(\kappa_{ik})
}
$$

where
$$
h(\kappa)=\frac{\sqrt{2\pi\kappa}}{2 -2e^{-2\kappa}}\mathrm{Erf}(\sqrt{2\kappa})
$$

$$
H(i) = \sum_{j\in\mathcal{N}_b(i)} h(\kappa_{ij})
$$

$$
\kappa^{-1}_{ij} \propto 2\pi\left(1-\frac{r_{ij}}{R_i}\right)
$$

$$
R_i = \frac{1}{N_b(i)-2}\sum_{j\in\mathcal{N}_b(i)}r_{ij}
$$