### Bond Orientational Order Parameters at 3D

The class `static.boo.boo_3d()` calculates the bond orientational order (BOO) parameters in three-dimensions. All calculations start from measuring the $l$-fold symmetry BOO as a $2l+1$ complex vector $q_{lm}$:

$$
q_{{lm}}(i) = \frac{1}{N_i} \sum_i^{N_i} Y_{{lm}} (\theta(\vec r_{ij}), \phi(\vec r_{ij})) 
\tag{1}
$$

where $Y_{lm}$ are the spherical harmonics and $N_i$ the number of bonds of particles $i$. This is the original one with all ($N_i$) neighboring particles treated the same. However, in some cases, we need to treat these neighboring particles differently based on their properties with respect to the center particle $i$. Then we can obtain the weighted $2l+1$ complex vector $q_{lm}$ as

$$
q_{{lm}}(i) = \sum_i^{N_i} \frac{A_j}{\sum A_j} Y_{{lm}} (\theta(\vec r_{ij}), \phi(\vec r_{ij})) 
\tag{2}
$$

where $A_j$ is the property of the $j-$th neighbor. Therefore, the normalized weights $\frac{A_j}{\sum A_j}$ are used in the calculation. We can treat $A_j=1$ for the case shown in equation (1). The weights are provided in the same way as for the neighbors.

Furthermore, dealing with $q_{lm}$ by coarse-graining is sometimes helpful by

$$
Q_{lm}(i)=\frac {1}{N_i+1} \left(q_{lm}(i) + \sum_{j=1} ^{N_i} q_{lm}(j) \right) 
\tag{3}
$$

Therefore, we can define the local complex vectors in a flexible way for different purposes. In principal, we name below the local quantities with small letters (associate with $q_{lm}$) and coarse-grained ones with large letters (associate with $Q_{lm}$). In the calculations one can set `coarse_graining=True` to calculate for the latter case.

By either using the complex vetcor $q_{lm}$ or $Q_{lm}$, various BOO parameters can be calculated. For example, 
+ the rotational invariants q_l (local, $q_l$) and Q_l (coarse-grained, $Q_l$) by `boo_3d.ql_Ql()`, 

$$
q_l=\sqrt{\frac{4\pi}{2l+1} \sum_{m=-l}^{l} \left| q_{lm} \right|^2}
\tag{4}
$$

+ local coherence of the complex vectors sij for $q_{lm}$ and Sij for $Q_{lm}$ by `boo_3d.sij_ql_Ql()`,

$$
s(i,j)=\frac {\sum_{m=-l}^l q_{lm}(i) q_{lm}(j)}{\sqrt{\sum_{m=-l}^l \left| q_{lm}(i) \right|^2}
\sqrt{\sum_{m=-l}^l \left| q_{lm}(j) \right|^2}}
\tag{5}
$$

+  w_l/w_cap_l (local, $w_l$, $\hat w_l$) and W_l/W_cap_l (coarse-grained, $W_l$, 
$\hat W_l$) by `boo_3d.w_W_cap()` 

$$
w_{l} = \sum_{m_{1} + m_{2} + m_{3} = 0}^{}
\left(
\begin{matrix}
l & l & l \\
m_{1} & m_{2} & m_{3} \
\end{matrix}
\right)
q_{lm_{1}}q_{lm_{2}}q_{lm_{3}} 
$$

$$
\hat{w_l}=w_l \left(\sum_{m=-l}^l \left| q_{lm} \right|^2 \right) ^{-\frac {3}{2}} 
\tag{7}
$$

+  spatial correlation of the complex vectors by `boo_3d.spatial_corr()`

$$
G_l(r)=\frac{4\pi}{2l+1} \frac{\sum_{ij} \sum_{m=-l}^l Q_{lm}(i) Q_{lm}^*(j) \delta(\vec r_{ij} - \vec r)} {\sum_{ij} \delta(\vec r_{ij} - \vec r)} 
\tag{8}
$$

+  time correlation of the complex vectors by `boo_3d.time_corr()`

$$
C_l(t)=\frac{4\pi}{2l+1} \frac {\left< \sum_i^N \sum_{m=-l}^l Q_{lm}(i, t) Q_{lm}^*(i, 0) \right>} {\left< \sum_i^N \sum_{m=-l}^l \left| Q_{lm}(i, 0) \right|^2 \right>} \tag{9}
$$

In equations (4)-(9), one can replace $q_{lm}$ by $Q_{lm}$ to calculate the corresponding coarse-grained values, or vice versa.

One example of using $s(i,j)$ with $q_{lm}$ for $l=6$ in computer simulations is to detect the bcc crystal nuclei. The bond between the center particle $i$ and a neighbor $j$ is considered crystalline if $s(i,j)>c$ (such as c=0.7). A particle is crystalline if the number of crystalline bond is larger than a threshold (usually half of the coordination number).

In these calculations, one important input is the neighbor list of each particle. Currently, the module accepts data generated by the `neighbors` module. One can also reweight BOO (see equation (2)) by using Voronoi polyhedron face area. Other weighting methods are acceptable as long as the correct format to account for the center-neighbor property is given, similar to the neighboring particles.

Note that this module accounts for both orthogonal and triclinic cells.


##### 1. `boo_3d` class

**Input Arguments**

- `snapshots` (`reader.reader_utils.Snapshots`): snapshot object of input trajectory (returned by reader.dump_reader.DumpReader)
- `l` (`int`): degree of spherical harmonics
- `neighborfile` (`str`): file name of particle neighbors (see module `neighbors`)
- `weightsfile` (`str`): file name of particle-neighbor weights (see module `neighbors`), default `None`.
  one typical example is Voronoi face area of the polyhedron; this file should be consistent with `neighborfile`. If provided, `weightsfile` will be read to calculate the weighted BOO, otherwise calculate the original BOO.
- `ppp` (`npt.NDArray`): the periodic boundary conditions, setting `1` for yes and `0` for no, default `np.array([1,1,1])`,
- `Nmax` (`int`): maximum number for neighbors, default 30

**Return**:
- None

**Example**
```python
from PyMatterSim.reader.dump_reader import DumpReader
from PyMatterSim.reader.reader_utils import DumpFileType
from PyMatterSim.neighbors.voropp_neighbors import cal_voro
from PyMatterSim.static.boo import boo_3d

filename = 'dump.atom'
readdump = DumpReader(filename, ndim=3, filetype=DumpFileType.LAMMPS, moltypes=None)
readdump.read_onefile()

cal_voro(readdump.snapshots, outputfile='test')
boo = boo_3d(readdump.snapshots,
             l=6,
             neighborfile='test.neighbor.dat'
             weightsfile='test.facearea.dat')
```

##### 2. `qlm_Qlm()`
`qlm_Qlm()` method gives $q_{lm}$ and $Q_{lm}$ values of different particles in complex vectors as numpy array. The results in
different snapshots are stored in an array for each of them, separately. Both of them are returned as ($q_{lm}$, $Q_{lm}$). Note that this method will be called automatically when initiallizing the `boo_3d()` class.

**Input Arguments**
- None

**Return**
- `smallqlm` (`npt.NDArray`): $q_{lm}$ in vector complex number, in numpy array with shape `[nsnapshots, nparticle, 2l+1]`
- `largeQlm` (`npt.NDArray`): $Q_{lm}$ in vector complex number, in numpy array with shape `[nsnapshots, nparticle, 2l+1]`


**Example**
```python
qlm, Qlm = boo.qlm_Qlm()
```

##### 3. `ql_Ql()`
`ql_Ql()` method calculates $q_l$ (local) or $Q_l$ (coarse-grained).

**Input Arguments**
- `coarse_graining` (`bool`): whether use coarse-grained $Q_{lm}$ or local $q_{lm}$, default `False`
- `outputfile` (`str`): file name for $q_l$ or $Q_l$ results, default `None`. To reduce storage size and ensure loading speed, save npy file as default with extension ".npy". If the file extension is ".dat" or ".txt", also saved a text file.

**Return**
- calculated $q_l$ or $Q_l$ in `npt.NDArray` with shape `[nsnapshots, nparticle]`

**Example**
```python
Ql = boo.ql_Ql(coarse_graining=True, outputfile='./results/Ql.npy')
```

##### 4. `sij_ql_Ql()`
`sij_ql_Ql()` method calculates the orientation correlation of $q_{lm}$ or $Q_{lm}$, named as $s_{ij}$ or $S_{ij}$. 

**Input Arguments**
- `coarse_graining` (`bool`): whether use coarse-grained $Q_{lm}$ or local $q_{lm}$, default `False`
- `c` (`float`): cutoff defining bond property, such as solid or not, default 0.7
- `outputqlQl` (`str`): npy file name of $q_l$ or $Q_l$ from particle-level sij, default `None`; The header is "id sum_sij num_neighbors", available from logger
- `outputsij` (`str`): npy file name for original $s_{ij}$ of $q_l$ or $Q_l$, default `None`; The header is "id CN sij", available from logger

**Return**
- aggregated $s_{ij}$ or $S_{ij}$ to particle-level in `npt.NDArray` with shape `[nsnapshots*nparticle, 3]`. This result is also stored in `outputqlQl` file for each particle.
- calculated $s_{ij}$ or $S_{ij}$ in `npt.NDArray` with shape `[nsnapshots*nparticle, 2+max_neighbors]`, where `max_neighbors` is the maximum number of neighbors and the first two columns are particle index and coordination number. The returned result is also stored in `outputsij` file for each pair

In `outputqlQl` file, the first column is particle index, the second is the sum of $s_{ij}$, and third one is number of neighbors.

`outputsij` file stores the $s_{ij}$ value of each particle with each of its neighbors. Results in different snapshots are written in sequence, with one header "id cn sij".

**Example**
```python
sij_Ql = boo.sij_ql_Ql(coarse_graining=True, 
                       outputqlQl='sum_sij.npy',
                       outputsij='sij_Ql.npy')
```

##### 5. `w_W_cap()`
`w_W_cap()` calculates the Wigner 3-j symbol boo based on qlm or Qlm.

**Input Arguments**
- `coarse_graining` (`bool`): whether use coarse-grained $Q_{lm}$ or local $q_{lm}$, default `False`
- `outputw` (`str`): txt file name for w (original) based on $q_{lm}$ or $Q_{lm}$, default `None`
- `outputwcap` (`str`): file name for wcap (normalized) based on qlm or Qlm, default None. To reduce storage size and ensure loading speed, save npy file as default with extension ".npy". If the file extension is ".dat" or ".txt", also saved a text file.

**Return**
- `w_W` (`np.adarray`), calculated $w$ or $W$ with shape `[nsnapshots, nparticle]`
- `w_W_cap` (`np.adarray`), calculated $w_{cap}$ or $W_{cap}$ with shape `[nsnapshots, nparticle]`

**Example**
```python
W, W_cap = boo.w_W_cap(coarse_graining=True, 
                       outputw='W.dat',
                       outputwcap='W_cap.npy')
```

##### 6. `spatial_corr()`
`spatial_corr()` method calculates spatial correlation function of $q_{lm}$ or $Q_{lm}$

**Input Arguments**
- `coarse_graining` (`bool`): whether use coarse-grained $Q_{lm}$ or local $q_{lm}$, default `False`
- `rdelta` (`float`): bin size in calculating g(r) and Gl(r), default 0.01
- `outputfile` (`str`): csv file name for $g_l$, default `None`
  
**Example**
```python
Gl_Q = boo.spatial_corr(coarse_graining=True, outputfile='Gl_Q.csv')
```

##### 7. `time_corr()`
`time_corr()` method calculates time correlation of $q_{lm}$ or $Q_{lm}$

**Input Arguments**
- `coarse_graining` (`bool`): whether use coarse-grained $Q_{lm}$ or local $q_{lm}$, default `False`
- `dt` (`float`): timestep used in user simulations, default 0.002
- `outputfile` (`str`): csv file name for time correlation results, default `None`

**Example**
```python
Gl_time = boo.time_corr(coarse_graining=True, outputfile='Gl_time.csv')
```