Bond order parameter - lan496の日記

結晶やガラスの局所的な環境を測る秩序変数のひとつにBond order parameterがある。

journals.aps.org

原子 $i$ の配位数を $N_{b}$ として、隣接する原子 $j (j=1, \cdots ,N_{b})$ との距離ベクトルを $\mathbf{r}_{ij}$ とする。

このとき次の量を考える。

$\overline{Q}_{lm} = \frac{1}{N_{b}} \sum_{j=1}^{N_{b}} Y_{lm}(\hat{ \mathbf{r} }_{ij} )$

ここで $Y_{lm}$ は球面調和関数。

この量自体はSO(3)の作用に関して不変でないが、次のように組み合わせると不変になる。

$Q_{l} := \sqrt{ \frac{ 4\pi }{ 2l + 1} \sum_{ m=-l }^{ l } | \overline{Q}_{lm} |^{2} }$

この $Q_{l}$ はBond(-orientational) order parameter(BOP)と呼ばれていて、例えばpymatgen.analysis.local_envに実装されている。

しかし、BOPを組み合わせて得られる次の不変量(論文ではthird-order invariant と呼ばれていた)の手頃な実装はなさそうだった。

$$ W_{l} = \sum_{ m_{1} + m_{2} + m_{3} } \begin{pmatrix} l & l & l \\ m_{1} & m_{2} & m_{3} \end{pmatrix} \overline{Q}_{ lm_{1} } \overline{Q}_{ lm_{2} } \overline{Q}_{ lm_{3} } $$ $\hat{ W_{l} } = \frac{ W_{l} }{ \left( \sum_{ m=-l }{ l } | \overline{Q}_{lm} |^{2} } \right)^{3/2} }$

ということで以下のとおり実装した。

折角書いたのでpymatgenにプルリクを投げようと思ったが、該当する秩序変数を求めるクラスが魔境と化していて何もわからないので止めておいた。

gist.github.com