SOAP Theory

Authors

Felix Musil, Max Veit, Michael Willatt, Klim Goldshtein

Last updated

15 Aug 2020

The mathematics behind the analytical radial integrals, and Cartesian gradients, of the SOAP descriptor 8 is described below.

From the density to the density expansion

Density expansion:

\[ \begin{align}\begin{aligned}\newcommand{\bvec}[1]{\mathbf{#1}} \newcommand{\CHF}[3]{{}_{1}F_{1}\left(#1 ,#2, #3 \right)} \newcommand{\CHLF}[2]{{}_{0}F_{1}\left(#1,#2\right)} \newcommand{\GA}[1]{\Gamma\left(#1\right)} \newcommand{\rij}{r_{ij}} \newcommand{\phiij}{\phi_{ij}} \newcommand{\uij}{u_{ij}} \newcommand{\thetaij}{\theta_{ij}} \newcommand{\br}{\mathbf{r}} \newcommand{\bhr}{\hat{\mathbf{r}}} \newcommand{\bm}[1]{\mathbf{#1}} \newcommand{\bu}{\mathbf{u}} \newcommand{\Rhat}{{\hat{R}}} \newcommand{\bk}{\mathbf{k}} \newcommand{\bK}{\mathbf{K}} \newcommand{\bt}{\mathbf{t}} \newcommand{\by}{\mathbf{y}} \newcommand{\bw}{\mathbf{w}} \newcommand{\bz}{\mathbf{z}} \newcommand{\bx}{\mathbf{x}} \newcommand{\bPsi}{\mathbf{\Psi}} \newcommand{\bPhi}{\mathbf{\Phi}} \newcommand{\dbr}{\textrm{d}\br} \newcommand{\dint}{\textrm{d}} \newcommand{\drhat}{\textrm{d}\hat{R}} \newcommand{\cald}{\mathcal{D}} \newcommand{\CX}{\mathcal{X}} \newcommand{\calk}{\mathcal{K}} \newcommand{\CA}{\mathcal{A}} \newcommand{\CB}{\mathcal{B}} \newcommand{\loss}{\ell} \newcommand{\diag}{\operatorname{diag}} % Stuff that was in the physics package (but mathjax can't use, obvs) \newcommand{\braket}[2]{\left\langle #1 \middle| #2 \right \rangle} \newcommand{\dd}{\mathrm{d}\;} \newcommand{\norm}[1]{\left\| #1 \right\|} \newcommand{\dv}[2]{\frac{\mathrm{d}\, #1}{\mathrm{d}\, #2}} \newcommand{\grad}{\nabla} \newcommand{\iket}[2]{\ket{#1}_{#2}} \newcommand{\ibraket}[3]{\braket{#1}{#2}_{#3}} \newcommand{\CN}{\mathcal{N}} \newcommand{\gaus}[2]{\CN_{#1}\left(#2\right)}\\\rho^{\alpha}_i(\mathbf{r})=\sum_{j\in \alpha} \mathcal{N}_{\sigma}\left(\mathbf{r}-\mathbf{r}_{ij}\right) = \sum_{nlm} \braket{\alpha n \ell m}{\mathcal{X}_i} B_{nlm}(\mathbf{r}),\end{aligned}\end{align} \]

where \(\alpha\) is a tag refers to the species of the considered atoms, \(\sum_{j \in \alpha}\) defines a sum over the neighbours of atom \(i\) of species \(\alpha\), \(\mathcal{N}_{\sigma}\) a Gaussian centered around \(0\) and variance \(\sigma^2\), \(B_{nlm}(\mathbf{r}) = R_n(r) Y_{\ell}^m(\hat{\mathbf{r}})\) defines a complete orthonormal basis set, \(Y_{\ell}^m\) a spherical harmonic (SH), \(r=\norm{\mathbf{r}}\) and \(\hat{\mathbf{r}}=\mathbf{r}/r\).

The following derivation aims at deriving expressions for the coefficients of the density expansion and their derivative with respect to atomic coordinate for several basis sets.

Density coefficients: angular integration

Spherical harmonics are the only orthonormal basis set of \(S^2\) so this part should apply except for real space basis.

We use the orthonormality of the basis set to compute the expressiont for the density coefficients and express the resulting integral over \(\rm I\!R^3\) in spherical coordinates using \(\norm{\mathbf{r}-\mathbf{r}_{ij}}^2=\mathbf{r}^2+\mathbf{r}_{ij}^2-2\norm{\mathbf{r}}\norm{\mathbf{r}_{ij}}\cos{\theta}\):

\[\begin{split}\begin{split} \braket{\alpha n \ell m}{\mathcal{X}_i}=& \sum_{j \in \alpha} c^{ij}_{n \ell m} = \sum_{j \in \alpha} \int_{\rm I\!R^3} \exp\left[-a\left(\mathbf{r}-\mathbf{r}_{ij}\right)^2\right]B_{nlm}(\mathbf{r})\\ =&\sum_{j \in \alpha} \int_{0}^{\infty}r^2 \exp\left[-a\left(r^2+r_{ij}^2\right)\right] g_n(r) \int_{-1}^{1}\mathrm{d}\left(\cos{\theta}\right) \\ & \int_0^{2\pi}\mathrm{d}\phi \exp\left[2arr_{ij}\cos{\theta}\right]Y_{\ell}^{m}\left(\hat{\mathrm{R}}\hat{\bm{q}}\right),\\ \end{split}\end{split}\]

where \(\hat{\mathrm{R}} = \hat{\mathrm{R}}_{ZYZ}\left(\alpha_{ij},\beta{ij},0\right)\) is the ZYZ-Euler matrix that rotate \(\hat{e}_z\) onto \(\bm{q}_{ij}\), \(a=\frac{1}{2\sigma^2}\). Note that global normalization constants are omitted because of a normalization at the end.

We use the following convention for the SH

\[Y_{\ell}^{m}\left(\hat{\mathbf{r}}\right)=Y_{\ell}^{m}\left(\theta,\phi\right)=A_{\ell}^{m}e^{im\phi}P^{m}_{\ell}\left(\cos\theta\right),\]

where \(A_{\ell}^{m} =\sqrt{\frac{(\ell-m)!(2l+1)}{4\pi(\ell+m)!}}\) and \(\hat{q}\) is the direction vector defined by the angle \(\theta\) and \(\phi\). Note that the phase factor \((-1)^\ell\) is included in the definition of the Associated Legendre Polynomials (ALPs). The set of spherical harmonics is calculated in the | librascal/src/math/spherical_harmonics.cc file, which is provided with explanations. The total number of SH in the set is \((l_{max} +1)^2\).

The integration over the angular part yields

\[\begin{split}\begin{split} \int_{-1}^{1}\mathrm{d}\left(\cos{\theta}\right) \exp\left[arr_{ij}\cos{\theta}\right] \int_0^{2\pi}\mathrm{d}\phi Y_{\ell}^{m}\left(\hat{\mathrm{R}}\left(\alpha_{ij},\beta{ij},0\right)\hat{\bm{q}}\right) =& 4\pi Y_\ell^m \left(\beta_{ij},\alpha_{ij}\right) \mathsf{i}_{\ell}\left(arr_{ij}\right), \\ =& Y_\ell^m \left(\hat{\mathbf{r}}_{ij}\right) \mathsf{i}_{\ell}\left(arr_{ij}\right), \end{split}\end{split}\]

where \(\mathsf{i}_{\ell}\) stands for the modified spherical Bessel function of the first kind. The intermediate steps are detailed in the following paragraphs.

Integration over \(\phi\)

The integration over the polar angle cancels out all orders of \(m\) from the SH

\[\int_0^{2\pi}\mathrm{d}\phi Y_{\ell}^{m}\left(\theta,\phi\right) = \sqrt{\pi\left(2l+1\right)}\mathrm{P}_{\ell}^{m}\left(\cos{\theta}\right) \delta_{m0},\]

since

\[\begin{split} \int_0^{2\pi}\mathrm{d}\phi \exp\left[im\phi\right] = 2\pi \delta_{0m}. \end{split}\]

Nevertheless, the rotation of the spherical harmonic breaks down into a linear combination of spherical harmonics. The coefficents are the entries of the Wigner D-matrix constructed from the Euler angles of the rotation matrix \(\hat{R}\).

\[\begin{split}\begin{split} Y_{\ell}^{m}\left(\hat{\mathrm{R}}\hat{\mathbf{r}}\right) = & \sum_{m'=-\ell}^{\ell} \mathrm{D}_{mm'}^\ell\left(\hat{R}\,\right) Y_{\ell}^{m}\left(\hat{\mathbf{r}}\right), \\ \mathrm{D}_{m0}^\ell\left(\alpha,\beta,\gamma\right) =& \sqrt{\frac{4\pi}{2l+1}} Y_l^m\left(\beta,\alpha\right). \\ \end{split}\end{split}\]

Thus, the polar integral over the rotated SH simplifies into

\[\begin{split}\begin{split} \int_0^{2\pi}\mathrm{d}\phi\, Y_{\ell}^{m}\left(\hat{\mathrm{R}}\hat{\mathbf{r}}\right) =& \sum_{m'=-\ell}^{\ell} \mathrm{D}_{mm'}^\ell\left(\alpha_{ij},\beta_{ij},0\right) \sqrt{\pi(2l+1)} P_{\ell}^{m'}\left(\cos{\theta}\right) \delta_{m'0} \\ =& 2\pi Y_\ell^m \left(\beta_{ij},\alpha_{ij}\right) P_{\ell}^{0}\left(\cos{\theta}\right). \end{split} \label{eq:int-phi}\end{split}\]

Integration over \(\theta\)

The modified spherical Bessel function of the first kind (MSBF) admit the following integral representation

\[\mathsf{i}_n\left(z\right)= \frac{1}{2}\int_{-1}^{1}\mathrm{d}x \exp\left(zx\right)\mathrm{P}_{n}^{0}\left(x\right),\]

which can be shown using the reference relations 1 2 3:

\[\begin{split}\begin{aligned} \mathsf{j}_n\left(z\right) =& \frac{(-i)^n}{2}\int_{-1}^{1}\mathrm{d}x \exp\left[izx\right]P_{n}^{0}\left(x\right), \label{eq:bessel-1}\\ \mathsf{i}_n\left(z\right)=& (-i)^{n} \mathsf{j}_n\left(iz\right), \label{eq:bessel-3}\\ \mathsf{i}_n\left(z\right)=& (-1)^{n} \mathsf{i}_n\left(-z\right), \label{eq:bessel-4}\end{aligned}\end{split}\]

1

http://dlmf.nist.gov/10.54.E2

2

http://dlmf.nist.gov/10.47.E12

3

http://dlmf.nist.gov/10.47.E16

\(j_n\) is the spherical Bessel function of the first kind. The integral over the polar angle is then given by

\[\begin{split} \int_{-1}^{1}\mathrm{d}\left(\cos{\theta}\right) \exp\left[2arr_{ij}\cos{\theta}\right]P_{\ell}^{0}\left(\cos{\theta}\right) =& 2 \mathsf{i}_{\ell}(2arr_{ij}). \end{split}\]

Density coefficients: Radial integration

Summing up the results from the previous section:

\[c^{ij}_{n\ell m} = 4\pi Y_{\ell}^m(\hat{\mathbf{r}}_{ij}) \exp\left[-ar^2_{ij}\right] \underbrace{\int_0^\infty \dd{r} r^2 R_n(r) e^{-ar^2} \mathsf{i}_{\ell}\left(2a r r_{ij}\right)}_{=\text{I}_{n\ell}^{ij}} ,\]

we identify \(\text{I}_{n\ell}^{ij}\) as the last term to simplify for particular choices of radial basis functions.

GTO like radial basis

The Gaussian Type Orbital radial basis is defined

\[R^{GTO}_{n}(r) = \mathcal{N}_n\ r^{n} \exp[-br^2],\]

where \(b=\frac{1}{2\sigma_n^2}\), \(\sigma_n = (r_\text{cut}-\delta r_\text{cut}) \max(\sqrt{n},1)/n_\text{max}\) and the normalization factor is given by

\[\mathcal{N}_n^2 = \frac{2(1)}{\sigma_n^{2n + 3}\Gamma(n + 3/2)}.\]

The overlap between GTO radial basis is:

\[\int_0^\infty R^{GTO}_{n}(r) R^{GTO}_{n^\prime}(r) \dd{r}= 2 \left(\frac{1}{2 \sigma_{n}^2}+\frac{1}{2 \sigma_{n^\prime}^2} \right)^{-\frac{1}{2} (3+n+n^\prime)} \Gamma(\frac{3+n+n^\prime}{2})\]

This equals what we use in the implementation

\[\int_0^\infty R^{GTO}_{n}(r) R^{GTO}_{n^\prime}(r) \dd{r}= N_n N_{n^\prime} \left(\frac{1}{2 \sigma_{n}^2}+\frac{1}{2 \sigma_{n^\prime}^2} \right)^{-\frac{1}{2} (3+n+n^\prime)} \Gamma(\frac{3+n+n^\prime}{2})\]

The radial integral becomes

\[I^{ij\,\text{GTO}}_{nl}= \mathcal{N}_n \frac{\sqrt{\pi}}{4} \frac{\GA{\frac{n+\ell+k+3}{2}}}{\GA{\ell+\frac{3}{2}}}a^\ell \rij^\ell(a+b)^{-\frac{n+k+\ell+3}{2}} \CHF{\frac{n+\ell+k+3}{2}}{\ell+\frac{3}{2}}{\frac{a^2 \rij^2}{a+b}}, \label{eq:rad-int-gto-1}\]

which yields the following expression for the neighbour contribution

\[\begin{split}c^{ij\,\text{GTO}}_{n\ell m}=& (\pi)^{\frac{3}{2}} \mathcal{N}_n \frac{\GA{\frac{n+\ell+3}{2}}}{\GA{\ell+\frac{3}{2}}} (a+b)^{-\frac{n+\ell+3}{2}} \\ & Y_{\ell}^m(\bhr_{ij}) \exp\left[-ar^2_{ij}\right] (a\rij)^\ell \CHF{\frac{n+\ell+3}{2}}{\ell+\frac{3}{2}}{\frac{a^2 \rij^2}{a+b}}. \label{eq:density-gto}\end{split}\]

where \(\Gamma\) is the Gamma function, and \({}_1F_1\) is the confluent hypergeometric function of the first kind.

The neighbour contribution is calculated in

file librascal/src/representations/ representation_manager_spherical_expansion.hh,

function compute_neighbour_contribution, line 338.

The steps of the derivation are detailed in the next paragraph.

Analytic radial integral

We write an integral representation of the confluent hypergeometric function \(\CHF{a}{b}{z}\) (CHF) in terms of MSBF:

\[\CHF{a}{\ell+\frac{3}{2}}{x} = \frac{2x^{-\frac{\ell}{2}}}{\sqrt{\pi}}\frac{\GA{\ell+\frac{3}{2}}}{\GA{a}}\int_0^\infty e^{-t} t^{a-1-\frac{\ell}{2}} \mathsf{i}_{\ell}(2\sqrt{xt})\dd{t}, \label{eq:chf-int}\]

using these relations 4 5 6

\[\begin{split}\begin{align} \CHF{a}{b}{z} = & \frac{1}{\GA{a}} \int_0^\infty e^{-t}t^{a-1}\CHLF{b}{zt}\dd{t},\\ I_l(z) =& \frac{(\frac{z}{2})^{\ell}}{\GA{l+1}} \CHLF{\ell+1}{\frac{z^2}{4}},\\ \mathsf{i}_{\ell}(z) =& \sqrt{\frac{\pi}{2z}}I_{\ell+1/2}(z),\\ \mathsf{i}_{\ell}(z) =& \sqrt{\frac{\pi}{4}}\frac{(\frac{z}{2})^{\ell}}{\GA{\ell+\frac{3}{2}}} \CHLF{\ell+\frac{3}{2}}{\frac{z^2}{4}},\\ \CHLF{\ell+\frac{3}{2}}{xt}=& \sqrt{\frac{4}{\pi}}\GA{\ell+\frac{3}{2}} x^{-\frac{\ell}{2}}t^{-\frac{\ell}{2}}\mathsf{i}_{\ell}(2\sqrt{xt}), \end{align}\end{split}\]

4

http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F1/07/01/01/0002/ http://dlmf.nist.gov/16.5.E3

5

https://en.wikipedia.org/wiki/Generalized_hypergeometric_function#The_series_0F1

6(1,2)

http://mathworld.wolfram.com/ModifiedSphericalBesselFunctionoftheFirstKind.html

where \(I_\ell\) is the modified Bessel function and \(\CHLF{b}{z}\) is the limit conflent hypergeometric function.

The module for calculating \(\CHF{..}{..}{..}\) is located in librascal/src/math/hyp1f1.hh.

The radial integral with GTO radial basis function is:

\[I^{ij\,\text{GTO}}_{nl}=\int_0^\infty \dd{r} r^{2+k} g^{\text{GTO}}_n(r) e^{-\frac{r^2}{2\sigma^2}} \mathsf{i}_{\ell}\left(r r_{ij} / \sigma^2\right) = \mathcal{N}_n \int_0^\infty \mathrm{d}r r^{2+k+n} e^{-r^2(a+b)} \mathsf{i}_{\ell}\left(2a r r_{ij}\right), \label{eq:rad-int-gto-0}\]

with \(k\) an additional power of \(r\) that will be non zero for the derivative. We partially identify the terms between [eq:chf-int] and [eq:rad-int-gto-0]:

\[\begin{split}\begin{aligned} t =& r^2(a+b),\\ \dd{t} =& 2 r \dd{r} (a+b),\\ x = & \frac{a^2 r_{ij}^2}{a+b},\end{aligned}\end{split}\]

to change the integrand of the radial integral

\[I^{ij\,\text{GTO}}_{nl}= \mathcal{N}_n \int_0^\infty \frac{\dd{t}}{2(a+b)} (a+b)^{-\frac{n+k+1}{2}} t^{\frac{n+k+1}{2}} e^{-t} \mathsf{i}_{\ell}\left(2\sqrt{xt}\right), \label{eq:rad-int-gto-01}\]

and identify the last term

\[\begin{aligned} a =& \frac{n+\ell+k+3}{2}.\end{aligned}\]

Numerical Integration of the Radial Integral

The numerical integration does not rely on a specific form of the radial basis

\[\text{I}_{n\ell}^{ij} = \sum_{k=1}^{K} \omega_k r_k^2 R_n(r_k) e^{-ar_k^2} \mathsf{i}_{\ell}\left(2a r_k r_{ij}\right),\]

where the \(\omega_k\) are the quadrature weights evaluated at the quadrature nodes \(r_k\). Depending on the quadrature rule, the following shifting formula is useful,

\[\int_a^b f(x)\,\dd{x} \approx \frac{b-a}{2} \sum_{i=1}^n w_i f\left(\frac{b-a}{2}x_i + \frac{a+b}{2}\right).\]

Discrete Variable Representation

In the special case of the the DVR radial basis 7 with Gauss-Legendre quadrature rule, the radial integral simplifies into:

\[\text{I}_{n\ell}^{ij} = \frac{r_c}{2} \sqrt{\omega_n} x_n^2 e^{-ax_n^2} \mathsf{i}_{\ell}\left(2a x_n r_{ij}\right),\]

where \(x_n=\frac{r_c}{2}r_n+\frac{r_c}{2}\).

7

Light, J. C., & Carrington, T. (2007). Discrete-Variable Representations and their Utilization (pp. 263–310). John Wiley & Sons, Ltd. https://doi.org/10.1002/9780470141731.ch4

Gradient of the density coefficients with respect to the Cartesian coordinates

The density coefficients can be split into two parts: one that depends on the choice of radial basis function (\(\text{I}_{n\ell}^{ij}\)) and the rest:

\[c^{ij}_{n\ell m} = Y_{\ell}^m(\hat{\mathbf{r}}_{ij}) \exp\left[-ar^2_{ij}\right] \text{I}_{n\ell}^{ij} = D^{ij}_{\ell m} C^{ij} \text{I}_{n\ell}^{ij},\]

where \(C^{ij}\) is the Gaussian exponential factor and \(\bar{D}^{ij}_{\ell m} = \bar{Y}_{\ell,m}(\hat{r}_{ij})\) is the spherical harmonic, see eq. [eq:real-spherical-harmonics]. Note the constant factors are omitted.

The following derivations end up with this formula that does not depend on the radial basis:

\[\begin{split}\begin{aligned} \grad_i\,c^{ij}_{\alpha n \ell m} =& 2a c^{ij}_{\alpha n \ell m} \mathbf{r}_{ij}\nonumber\\ &{} + C \bar{D}^{ij}_{\ell m} \cdot \grad_i \text{I}_{n\ell}^{ij}\nonumber\\ &{} + N_{n \ell}A_{n\ell} B_\ell C \cdot \grad_i\,\bar{D}^{ij}_{\ell,m},\end{aligned}\end{split}\]

where \(\grad_i\bar{D}^{ij}_{\ell,m} = \grad_i \bar{Y}_{\ell,m}(\hat{r}_{ij})\) is defined in [eq:dbx0,eq:dbx1,eq:dbx2,eq:dby0,eq:y1,eq:dby2,eq:dbz0,eq:dbz1,eq:dbz2].

Terms common to the different radial basis

Gaussian

\[\begin{gathered} \dv{C^{ij}}{r_{ij}} = -2ar_{ij}C^{ij}\end{gathered}\]

Length

So for the radial terms, we just use the derivatives of the radius \(r_{ij}\) wrt the Cartesian coordinates:

\[\begin{split}\begin{gathered} \dv{ r_{ij}}{ \{x_i, y_i, z_i\}} = -\frac{\{x_{ij}, y_{ij}, z_{ij}\}}{r_{ij}}\\ \grad_i\,r_{ij} = \frac{-\mathbf{r}_{ij}}{r_{ij}}\\ \text{where }\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i\end{gathered}\end{split}\]

Spherical Harmonics

The derivative of the spherical harmonic can be expressed in a few different ways. The versions below are in terms of the original harmonic with possibly different \(m\) values. The \(z\) component is:

\[\begin{split}\begin{aligned} \frac{\partial D_{\ell m}}{\partial z_i} &= \frac{-\sqrt{1-u^2}}{2r}\big(e^{i\phi}\sqrt{(\ell+m)(\ell-m+1)}Y_l^{m-1}(\hat{r})\nonumber\\ &\qquad\qquad - e^{-i\phi}\sqrt{(\ell-m)(\ell+m+1)}Y_l^{m+1}(\hat{r})\big)\nonumber\\ &= \frac{-\sin{\theta}}{2r_{ij}}(\cos(m\phi) + i\sin(m\phi)) \\ &\qquad\qquad \left(\sqrt{(\ell+m)(\ell - m + 1)}\sqrt{\frac{2\ell+1}{4\pi}\frac{(\ell-m+1)!}{(\ell+m-1)!}} P_l^{m-1}(\cos{\theta})\right.\nonumber\\ &\qquad\qquad\qquad \left. {} - \sqrt{(\ell-m)(\ell + m + 1)}\sqrt{\frac{2\ell+1}{4\pi}\frac{(\ell-m-1)!}{(\ell+m+1)!}} P_l^{m+1}(\cos{\theta})\right) \end{aligned}\end{split}\]

But remember, we’re actually using the real spherical harmonics:

[eq:real-spherical-harmonics]

\[\begin{split}\begin{aligned} \left.\begin{aligned} \bar{Y}_{\ell m}(\hat{r}_{ij}) &= \cos(m\phi) \bar{P}_\ell^m(\cos{\theta})\\ \bar{Y}_{\ell,-m}(\hat{r}_{ij}) &= \sin(m\phi) \bar{P}_\ell^m(\cos{\theta}) \end{aligned}\right\}&\text{ for }m > 0\\ \bar{Y}_{\ell,0}(\hat{r}_{ij}) = \frac{1}{\sqrt{2}} \bar{P}_\ell^0(\cos{\theta})&\end{aligned}\end{split}\]

where

\[\bar{P}_\ell^m(\cos{\theta}) = \sqrt{\frac{2\ell + 1}{2\pi}\frac{(\ell - m)!}{(\ell + m)!}}P_\ell^m(\cos{\theta}).\]

So we can write

\[\begin{split}\begin{aligned} \frac{\partial \bar{D}_{\ell m}}{\partial z_i} &= \frac{-\sin\theta}{2r_{ij}}\cos(m\phi)\left(\sqrt{(\ell + m)(\ell - m + 1)}\bar{P}_\ell^{m-1}(\cos\theta) - \sqrt{(\ell - m)(\ell + m + 1)}\bar{P}_\ell^{m+1}(\cos\theta)\right) \label{eq:dbz0}\\ \frac{\partial \bar{D}_{\ell,-m}}{\partial z_i} &= \frac{-\sin\theta}{2r_{ij}}\sin(m\phi)\left(\sqrt{(\ell + m)(\ell - m + 1)}\bar{P}_\ell^{m-1}(\cos\theta) - \sqrt{(\ell - m)(\ell + m + 1)}\bar{P}_\ell^{m+1}(\cos\theta)\right)\label{eq:dbz1}\\ \frac{\partial \bar{D}_{\ell,0}}{\partial z_i} &= \frac{\sin\theta}{r_{ij}} \sqrt{\frac{\ell(\ell + 1)}{2}}\bar{P}_\ell^{1}(\cos\theta))\label{eq:dbz2}\end{aligned}\end{split}\]

(the last one comes from the identity \(\sqrt{\frac{(\ell+m)!}{(\ell-m)!}}P_\ell^{-m} = (-1)^m \sqrt{\frac{(\ell - m)!}{(\ell + m)!}}P_l^m(\cos\theta)\) with \(m=1\)).

The \(x\) component is:

\[\begin{split}\begin{aligned} \frac{\partial \bar{D}_{\ell m}}{\partial x_i} &= \frac{-m\sin\phi}{\sqrt{x_{ij}^2 + y_{ij}^2}} \bar{D}_{\ell,-m} + \frac{\cos\phi \cos\theta}{2r_{ij}}\cos(m\phi)\left( \sqrt{(\ell + m)(\ell - m + 1)}\bar{P}_\ell^{m-1}(\cos\theta)\right.\nonumber\\ &\qquad\qquad\qquad\left. {} - \sqrt{(\ell - m)(\ell + m + 1)}\bar{P}_\ell^{m+1}(\cos\theta)\right)\label{eq:dbx0}\\ \frac{\partial \bar{D}_{\ell,-m}}{\partial x_i} &= \frac{m\sin\phi}{\sqrt{x_{ij}^2 + y_{ij}^2}} \bar{D}_{\ell,m} + \frac{\cos\phi \cos\theta}{2r_{ij}}\sin(m\phi)\left( \sqrt{(\ell + m)(\ell - m + 1)}\bar{P}_\ell^{m-1}(\cos\theta)\right.\nonumber\\ &\qquad\qquad\qquad\left. {} - \sqrt{(\ell - m)(\ell + m + 1)}\bar{P}_\ell^{m+1}(\cos\theta)\right)\label{eq:dbx1}\\ \frac{\partial \bar{D}_{\ell,0}}{\partial x_i} &= \frac{-\cos\phi \cos\theta}{r_{ij}}\sqrt{\frac{\ell(\ell+1)}{2}}\bar{P}_\ell^1(\cos\theta)\label{eq:dbx2}\end{aligned}\end{split}\]

and for the \(y\) component, similarly:

\[\begin{split}\begin{aligned} \frac{\partial \bar{D}_{\ell m}}{\partial y_i} &= \frac{m\cos\phi}{\sqrt{x_{ij}^2 + y_{ij}^2}} \bar{D}_{\ell,-m} + \frac{\sin\phi \cos\theta}{2r_{ij}}\cos(m\phi)\left( \sqrt{(\ell + m)(\ell - m + 1)}\bar{P}_\ell^{m-1}(\cos\theta)\right.\nonumber\\ &\qquad\qquad\qquad\left. {} - \sqrt{(\ell - m)(\ell + m + 1)}\bar{P}_\ell^{m+1}(\cos\theta)\right)\label{eq:dby0}\\ \frac{\partial \bar{D}_{\ell,-m}}{\partial y_i} &= \frac{-m\cos\phi}{\sqrt{x_{ij}^2 + y_{ij}^2}} \bar{D}_{\ell,m} + \frac{\sin\phi \cos\theta}{2r_{ij}}\sin(m\phi)\left( \sqrt{(\ell + m)(\ell - m + 1)}\bar{P}_\ell^{m-1}(\cos\theta)\right.\nonumber\\ &\qquad\qquad\qquad\left. {} - \sqrt{(\ell - m)(\ell + m + 1)}\bar{P}_\ell^{m+1}(\cos\theta)\right)\label{eq:dby1}\\ \frac{\partial \bar{D}_{\ell,0}}{\partial y_i} &= \frac{-\sin\phi \cos\theta}{r_{ij}}\sqrt{\frac{\ell(\ell+1)}{2}}\bar{P}_\ell^1(\cos\theta)\label{eq:dby2}\end{aligned}\end{split}\]

The formulæ above have a singularity at the poles for \(m \neq 0\), so use the following identity:

\[\begin{split}\begin{gathered} \frac{m}{\sqrt{x_{ij}^2 + y_{ij}^2}} \begin{pmatrix}\bar{Y}_{\ell, -m}(\hat{r}_{ij})\\ \bar{Y}_{\ell, m}(\hat{r}_{ij})\end{pmatrix} = \frac{-1}{2z_{ij}}\begin{pmatrix}\sin(m\phi)\\\cos(m\phi)\end{pmatrix} \left(\sqrt{(\ell+m)(\ell - m + 1)}\bar{P}_\ell^{m-1}(\cos\theta) \right.\\ \left. {} + \sqrt{(\ell - m)(\ell + m + 1)}\bar{P}_\ell^{m+1}(\cos\theta)\right)\end{gathered}\end{split}\]

to shift the singularity to the equator (\(z=0\)). In the code derivatives of spherical harmonics is computed in the feat/soap_gradients branch, librascal/src/math/spherical_harmonics.hh

GTO like radial basis

We rewrite [eq:rad-int-gto-1]

\[I^{ij\,\text{GTO}}_{nl} = N_{n\ell} \cdot A_{n\ell} \cdot B_\ell ,\]

where \(B_{\ell} = r_{ij}^{\ell}\), \(A_{n\ell} = \CHF{\frac{n + \ell + 3}{2}}{\ell+\frac{3}{2}}{\frac{a^2 r_{ij}^2}{a+b}}\), \(N_{n \ell} = \frac{\mathcal{N}_n}{4} a^\ell\left(a+b\right)^{-\frac{n + \ell + 3}{2}} \frac{\Gamma\left(\frac{n + \ell + 3}{2}\right)}{\Gamma\left(\frac{3}{2} + \ell\right)}\), \(\mathcal{N}_n = \sqrt{\frac{2}{\sigma_n^{2n+3}\Gamma\left(n + \frac{3}{2}\right)}}\). Note that some constant multiplying factors of \(\pi\) have been omitted.

\(B_{\ell}\)

\[\dv{B_\ell}{r_{ij}} = \frac{\ell}{r_{ij}} B_\ell\]

CHF

for the hypergeometric term:

\[\dv{A_{n \ell}}{r_{ij}} = \frac{\frac{n + \ell + 3}{2}}{\left(\ell + \frac{3}{2}\right)} \frac{2a^2 r_{ij}}{a+b} \CHF{\frac{n + \ell + 5}{2}}{\ell+\frac{5}{2}}{\frac{a^2 r_{ij}^2}{a+b}}\]

which is not proportional to \(A_{n \ell}\), or even to \(A_{n+1,\ell + 1}\) – so just recompute it explicitly.

GTO formula for practical computation

Finally, putting the radial and angular components together, we get:

\[\begin{split}\begin{aligned} \grad_i\,c^{ij}_{\alpha n \ell m} &= c^{ij}_{\alpha n \ell m}\left(-\frac{\ell}{r_{ij}^2} + 2a\right)\br_{ij}\nonumber\\ &{} + N_{n \ell}B_\ell C \bar{D}_{\ell m} \cdot \frac{\frac{n + \ell + 3}{2}}{\left(\ell + \frac{3}{2}\right)} \frac{2a^2}{a+b} \CHF{\frac{n + \ell + 5}{2}}{\ell+\frac{5}{2}}{\frac{a^2 r_{ij}^2}{a+b}} \bvec{r}_{ij}\nonumber\\ &{} + N_{n \ell}A_{n\ell} B_\ell C \cdot \nabla_i\,\bar{D}_{\ell,m} \end{aligned}\end{split}\]

where the gradient of the spherical harmonic has already been computed separately using the equations above.

Gradient of the coefficients is calculated in feat/soap_gradients branch,

file librascal/src/representations/representation_manager_spherical_expansion.hh,

function compute_neighbour_derivative, line 420.

Numerical Integration

Using the recurrence relation of the MSBF 6:

\[\dv{\mathsf{i}_{\ell}(x)}{x} = \frac{1}{2\ell+1}[\ell\mathsf{i}_{\ell-1}(x)+(\ell+1)\mathsf{i}_{\ell+1}(x)],\]

the gradient of the radial integral becomes:

\[\grad_i \text{I}_{n\ell}^{ij} = -\frac{2a}{2\ell+1}\sum_{k=1}^{K} \omega_k r_k^3 R_n(r_k) e^{-ar_k^2} [\ell\mathsf{i}_{\ell-1}(2a r_k r_{ij})+(\ell+1)\mathsf{i}_{\ell+1}(2a r_k r_{ij})] \hat{\mathbf{r}}_{ij}.\]

In the case of the DVR radial basis:

\[\text{I}_{n\ell}^{ij} = -\frac{2a\sqrt{\omega_n}}{2\ell+1}\frac{r_c}{2} x_n^3 e^{-ax_n^2} [\ell\mathsf{i}_{\ell-1}(2a x_n r_{ij})+(\ell+1)\mathsf{i}_{\ell+1}(2a x_n r_{ij})] \hat{\mathbf{r}}_{ij},\]

where \(x_n=\frac{r_c}{2}r_n+\frac{r_c}{2}\).

Beyond SOAP

SOAP can be seen as one of the simplest members of a hierarchy of “density correlation features”, that are obtained by appropriately symmetrizing tensor products of the neighbor density. The formalism, first introduced in Ref. 9, leads to formulas to evaluate discretized versions of these features as a combination of the density expansion coefficients, and includes also equivariant features, that transform as spherical harmonics under rotation 10. A comparatively simple expression to compute these higher-order features iteratively, based on an angular momentum combination relation has been discussed in Ref. 11, as part of the N-body iterative contraction of equivariants (NICE) framework. Even though these higher-order features are not the main focus of librascal, you can find some utilities that compute them starting from the expansion coefficients. These are part of the of the bindings/rascal/utils/cg_utils.py, and are demonstrated and tersely documented in the example notebooks examples/equivariant_demo.ipynb and examples/nice_demo.ipynb.

8

A. P. Bartók, R. Kondor, and G. Csányi (2013) On representing chemical environments. Physical Review B 87(18), 184115. https://doi.org/10.1103/PhysRevB.87.184115

9

M. J. Willatt, F. Musil, and M. Ceriotti (2019) Atom-density representations for machine learning. Journal of Chemical Physics 150(15), 154110. https://doi.org/10.1063/1.5090481

10

A. Grisafi, D. M. Wilkins, G. Csányi, and M. Ceriotti (2018) Symmetry-Adapted Machine Learning for Tensorial Properties of Atomistic Systems. Physical Review Letters 120(3), 036002. https://doi.org/10.1103/PhysRevLett.120.036002

11

J. Nigam, S. Pozdnyakov, and M. Ceriotti (2020) Recursive evaluation and iterative contraction of N -body equivariant features. J. Chem. Phys. 153(12), 121101. https://doi.org/10.1063/5.0021116