SOAP polynomial RBF error

[pujaltes]

Hello,

Thanks for providing this amazing resource!

I believe I found an error in the construction of the polynomial RBF for the soap descriptor when trying to replicate figure 5b (below) in the original Dscribe paper:

The following code should provide the exact same result as in the paper:

from dscribe.descriptors import SOAP
import matplotlib.pyplot as plt

soap = SOAP(
    species=["H"],
    r_cut=3,
    n_max=4,
    l_max=2,
    rbf="polynomial",
)
rx, gss = soap.get_basis_poly(3, 4)
plt.plot(rx, gss.T)

However, this instead results in the following RBFs which are not orthonormal:

The error seems to originate from the get_basis_poly method, specifically in the construction of the Lowdwin orthogonalization matrix S:

dscribe/dscribe/descriptors/soap.py

Lines 1076 to 1085 in 8379685

	S = np.zeros((n_max, n_max), dtype=np.float64)
	for i in range(1, n_max + 1):
	for j in range(1, n_max + 1):
	S[i - 1, j - 1] = (2 * (r_cut) ** (7 + i + j)) / (
	(5 + i + j) * (6 + i + j) * (7 + i + j)
	)
	# Get the beta factors that orthonormalize the set with Löwdin
	# orthonormalization
	betas = sqrtm(np.linalg.inv(S))

The calculation used for S in the for loop does not follow the original equations from Bartók

$$ \begin{aligned} &g_n(r)=\sum_{\alpha=1}^{n_{max }} W_{n \alpha} \phi_\alpha(r) \\ &\phi_\alpha(r)=\frac{(r_{cut}-r)^{\alpha + 2}}{N_\alpha} \\ &N_\alpha=\sqrt{\int_0^{r_{\text {cut }}}\left(r_{\text {cut }}-r\right)^{2(\alpha+2)} d r}=\sqrt{\frac{r_{\text {cut }}^{2 \alpha+5}}{2 \alpha+5}} \\ &S_{\alpha \beta}=\int_0^{r_{\text {cut }}} \phi_\alpha(r) \phi_\beta(r) d r=\frac{\sqrt{(5+2 \alpha)(5+2 \beta)}}{5+\alpha+\beta}, \\ &\mathbf{W}=\mathbf{S}^{-1 / 2} . \end{aligned} $$

Modifying the above section of code to the following produces the correct/expected result:

S = np.zeros((n_max, n_max), dtype=np.float64)
norm_factors = np.zeros(n_max, dtype=np.float64)
for i in range(1, n_max + 1):
    norm_factors[i - 1] = np.sqrt(r_cut ** (2 * i + 5) / (2 * i + 5))
    for j in range(1, n_max + 1):
        S[i - 1, j - 1] = np.sqrt((5 + 2 * j) * (5 + 2 * i)) / (5 + i + j)

# Get the beta factors that orthonormalize the set with Löwdin
# orthonormalization
betas = sqrtm(np.linalg.inv(S)) / norm_factors[None, :]

I would be happy to correct this and open a pull request if the error is confirmed.

[lauri-codes]

Hi @pujaltes!

Thank you for your note and keen observations on the code.

I had a quick look on the subject, and found a probable explanation. The equation by Bartok et. al. in the original paper has a small problem in that the integrals over the radial component do not contain the spherical coordinate Jacobian (r^2). This is mentioned in the errata that the authors have posted, and you can find it here. After you include it, you should find that our implementation is correct (when checking orthogonality, you also have to include the r^2 in the integration).

I'm not sure about the figures in our paper: I will have to take a deeper look at whether we have plotted them with the incorrect data.

[lauri-codes]

Here is an example that demonstrates that the radial basis functions are orthonormal (within numerical accuracy) when integrated under spherical coordinates:

import numpy as np
import matplotlib.pyplot as plt
from dscribe.descriptors import SOAP

n_max = 4
soap = SOAP(
    species=["H"],
    r_cut=3,
    n_max=n_max,
    l_max=2,
    rbf="polynomial",
)
r, g = soap.get_basis_poly(3, n_max)

integrals = np.empty((n_max, n_max))
for i in range(n_max):
    for j in range(n_max):
        # note the r^2 coming from the Jacobian of the spherical coordinates
        integrals[i, j] = np.trapz(g[i, :]*g[j, :]*r**2, r) 
print(integrals)

This means that the plot in our original paper is done using the incorrect basis set, but our current implementation should be correct.

[pujaltes]

@lauri-codes thank you very much for taking a look at this and for your detailed response! This makes a lot of sense, thanks for clearing it up.