add bispectrum coefficients to matminer

[computron]

Bispectrum coefficients are widely used in the community, including in the recent GAP and SNAP potentials. IMO they are one of the most important site descriptors to have in matminer as a reference for other site descriptors.

The formalism of these coefficients is here:

1. Bartók, A. P., Payne, M. C., Kondor, R. & Csányi, G. Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons. Phys. Rev. Lett. 104, 1–4 (2010).

Implementation of these coefficients is already in the LAMMPS code in C++:

https://github.com/lammps/lammps/blob/654fec164d8f95507326086c58f0571087518ef9/src/SNAP/sna.cpp

with further notes from Aidan Thompson, LAMMPs developer

You should ignore derivatives/forces, ignore OMP variants, and assume diagonalstyle == 3

The key functions are, in order

compute_ui() // hyperspherical harmonic contributions summed over neighbors j of atom i
compute_zi() // pairwise products of Ui's
compute_bi() // bispectrum components = triplewise products of Ui's
copy_bi2bvec() // copy bispectrum components to 1D list of floats

We should either have a way to call the LAMMPs code to get these coefficients (will be a little complex for users to compile and also to implement the call), or reimplement them in Python

[athomps]

I recommend the term bispectrum component, slightly less generic than coefficient. More information on the LAMMPS algorithm is provided here:

Thompson, Swiler, Trott, Foiles, Tucker, J Comp Phys, 285, 316 (2015).

[WardLT]

The implementation for bispectrum components will likely be similar to the BondOrientationalOrder class. The main hurdles to implementing it will be developing functions for the hyperspherical harmonics and coupling constants (H). I've already done some research on these:

• Hyperspherical harmonics can be expressed as products of Legendre polynomials, which are implemented in scipy.
• The coupling constants (equation A.2 of the paper @athomps mentioned) are defined as products of Clebsch-Gordan Coefficients, which are available in sympy. Sympy is already a dependency of matminer and can be expensive, so you'll want to cache the results.

I'm not 100% on either of these suggestions but, hopefully, they give a good starting point. I hope to put some work in on this in the next few weeks, so let me know if anyone is interested in starting on this before me.

[athomps]

The BondOrientationalOrder class is indeed very elegant and it would make sense to start by doing something similar for bispectrum of hyperspherical harmonics on 3-sphere in R^4 (HSH).

Unfortunately, I could not quite make sense of the expression for HSH as products of Legendre functions. The expression in R^n due to Atsushi specializes correctly for n=3 (spherical harmonics) but when I try it for n=4 (HSH) I get this expression Here is the latex source: $Y_{l_1l_2l_3}(\theta_1, \theta_2, \theta3) = Y{l_2}^{-l_1}(\theta_2, -\theta_1) \sqrt{\frac{(l_3+1)}{\sin \theta_3}\frac{(l_3+1_2+1)!}{(l_3-l2)!}} P{l_3+\frac{1}{2}}^{-l_2-\frac{1}{2}}(\cos\theta_3)$ The problem is that the indices to the Legendre function are half integer, which does not make sense to me. Maybe I am misunderstanding something in Atsushi's notation. I looked in the very extensive handbook by Varshalovich, Moskalev, Khersonskii (VMK), but did not see any formula resembling Atsushi's.

The only practical way I know to compute a particular HSH is to generate the full set up to some order using the recursion method implemented in LAMMPS, based on formulae given in VMK Section 4.8.2

[ardunn]

superseded by #635