Non zero gradients when all neighbors are periodic images

[Luthaf]

Running the code below prints False, but I would expect it to print True.

import ase
from ase import build
​
import numpy as np
​
from rascal.representations import SphericalExpansion
​
se = SphericalExpansion(
    interaction_cutoff=4.5,
    cutoff_smooth_width=0.5,
    max_radial=6,
    max_angular=6,
    gaussian_sigma_type="Constant",
    gaussian_sigma_constant=0.3,
    radial_basis="GTO",
    compute_gradients=True,
)
​
frame = ase.build.bulk("H", "fcc", 2.2)
​
representation = se.transform(frame)
gradients = representation.get_features_gradient(se)
​
print(np.all(gradients == 0.0))

This construct a structure containing a single atom, and then try to compute the spherical expansion using a large cutoff, meaning the atom have a lot of neighbors, all being periodic images of itself.

---

I'm wondering whether this is related to https://github.com/cosmo-epfl/librascal/issues/324#issuecomment-824934846: there are also duplicated entries inside the gradients (the gradient shape is (141, 3, 294) while I would expect it to be (1, 3, 294) since there is one environment with only one neighbor).

---

Changing the cutoff in related tests to 5.5 (from 2.5) make the tests fail, so I'm relatively sure this is a bug.

https://github.com/cosmo-epfl/librascal/blob/43a410c6d7cea3d17852e9c8c6fe332ec46f13df/tests/test_calculator.hh#L428

[max-veit]

Related: #115 #179 #213

Might be a good idea to do a bisection to see when this was introduced, because I'm pretty sure these tests were working with a larger cutoff when I first added them.

[agoscinski]

Hello all,

I think the output needs to be partially multiplied with the direction vectors so we get zero gradients for this frame. Firstly we need to consider the sum of the gradient, because we also output the self contribution (this has been already discussed with @Luthaf last week but it was not enough to get zeros).

Secondly I think the direction vector is not included in the computation for the spherical expansion gradient. I have printed two pair gradients of opposing ghost neighbours (central atom is at position (0.5, 0.5, 0.5) and we have a cubic cell [1, 1, 1])

Opposing neighbours have same l=1 components (n_max=1)

center atom tag 5  position -0.5  0.5  0.5
xyz sph exp gradients for pair atom tag (0, 5)
     l=0                     l=1
     m=0         m=1         m=0         m=-1
[ 0.00823064  0.          0.         -0.01274583]
[ 0.          0.00064061 -0.          0.        ]
[ 0.         -0.          0.00064061  0.        ]

center atom tag 22 position  1.5  0.5  0.5
xyz sph exp gradients for pair atom tag (0, 22)
     l=0                     l=1
     m=0         m=1         m=0         m=-1
[-0.00823064  0.         -0.         -0.01274583]
[ 0.          0.00064061 -0.          0.        ]
[ 0.         -0.          0.00064061 -0.        ]

So we need to multiply l mod 2 = 1 components with the direction vector and then it works out

import ase
from ase import build
import numpy as np
from rascal.representations import SphericalExpansion

hypers = dict(
    interaction_cutoff=4.5,
    cutoff_smooth_width=0.5,
    max_radial=6,
    max_angular=6,
    gaussian_sigma_type="Constant",
    gaussian_sigma_constant=0.3,
    radial_basis="GTO",
    compute_gradients=True
)
se = SphericalExpansion(**hypers)

frame = ase.build.bulk("H", "fcc", 2.2)
#frame = ase.Atoms('H', positions=[(0.5,0.5,0.5)], cell=[1,1,1])
frame.pbc = True

representation = se.transform(frame)
gradients = representation.get_features_gradient(se)
dir_vec = representation.get_direction_vectors()

# (i, anlm) -> (i, a, n, l, m)
def relayout(x, hypers):
    nmax = hypers['max_radial']
    lmax = hypers['max_angular']
    n = len(x)
    nel = x.shape[1]//(nmax*(lmax+1)**2)
    rx = x.reshape((n, nel, nmax, (lmax+1)**2))
    nx = np.zeros((n, nel, nmax, lmax+1, 2*lmax+1))   # lazy layout for easier manipulation
    for l in range(lmax+1):
        nx[:,:,:,l,:2*l+1] = rx[:,:,:,l**2:(l+1)**2]
    return nx

gradients = relayout(gradients, hypers)

sum_gradients_x = np.zeros((1, gradients.shape[1], gradients.shape[2], gradients.shape[3], gradients.shape[4]))
sum_gradients_y = np.zeros((1, gradients.shape[1], gradients.shape[2], gradients.shape[3], gradients.shape[4]))
sum_gradients_z = np.zeros((1, gradients.shape[1], gradients.shape[2], gradients.shape[3], gradients.shape[4]))

for i, sample in enumerate(dir_vec):
  for angular_l in range(hypers['max_angular']):
    if angular_l % 2 == 1:
      sum_gradients_x += gradients[3 * i + 0, :, :, angular_l:angular_l+1, :] * dir_vec[i, 0]
      sum_gradients_y += gradients[3 * i + 1, :, :, angular_l:angular_l+1, :] * dir_vec[i, 1]
      sum_gradients_z += gradients[3 * i + 2, :, :, angular_l:angular_l+1, :] * dir_vec[i, 2]
    else:
      sum_gradients_x += gradients[3 * i + 0, :, :, angular_l:angular_l+1, :]
      sum_gradients_y += gradients[3 * i + 1, :, :, angular_l:angular_l+1, :]
      sum_gradients_z += gradients[3 * i + 2, :, :, angular_l:angular_l+1, :]

print(np.allclose(sum_gradients_x, np.zeros(sum_gradients_x.shape)))
print(np.allclose(sum_gradients_y, np.zeros(sum_gradients_y.shape)))
print(np.allclose(sum_gradients_z, np.zeros(sum_gradients_z.shape)))

Output

True
True
True

I am not sure why librascal output it this way but it seems to be the solution to the problem. This needs to be definitely better documented

[agoscinski]

Small update, I tried if the same thing happens for a nonghost atoms (within the cell) and here one does not need to include the direction vector. So it is inconsistent behavior with respect to the ghosts which is a bit worrisome

import ase
from ase import build

import numpy as np
from rascal.representations import SphericalExpansion
se = SphericalExpansion(
    interaction_cutoff=1.1,
    cutoff_smooth_width=0.5,
    max_radial=6,
    max_angular=6,
    gaussian_sigma_type="Constant",
    gaussian_sigma_constant=0.3,
    radial_basis="GTO",
    compute_gradients=True,
)
frame = ase.Atoms('HHHHHHH', positions=[(1,1,1),(0.5,1,1),(1.5,1,1),(1,0.5,1),(1,1.5,1),(1,1,0.5),(1,1,1.5)], pbc=[False,False,False], cell=[2,2,2])
representation = se.transform(frame)
gradients = representation.get_features_gradient(se)
info = representation.get_gradients_info()
# shape (n_atoms*(n_neigh[i]+1)*xyz, n_feat)
# chose gradients corresponding to the atom in the middle of the two
# middle atom is at the beginning and has shape (6 [neighbours] + 1 [self contribution]) * 3 [xyz]  = 21
print("x close to 0", np.allclose(np.sum(gradients[0:21][0::3]),0))
print("y close to 0", np.allclose(np.sum(gradients[0:21][1::3]),0))
print("z close to 0", np.allclose(np.sum(gradients[0:21][2::3]),0))

Output

True
True
True

[agoscinski]

Hello,

@Luthaf do you still get an error when you increase the cutoff as you suggested in your first response? I don't get it and I checked the gradient provider in the tests it seems to ignore ghost atoms when computing the analytical gradients.

It seems that the fix of issue #115 is related to this, https://github.com/cosmo-epfl/librascal/pull/179/files#diff-fed2083c2fcd751ec9a9e5f9894ab837906968d2ad46c2e70849c743c4da5442R1183-R1191 if I change in the code in the diff

        // gives you always the within-cell atom tag of atom j
        auto atom_j_tag = atom_j.get_atom_tag();

to

        // gives you always the tag of the j atom, even if j is ghost
        auto atom_j_tag = neigh.get_atom_tag();

then the python script gives me zero without the contribution of direction vector, but the cpp tests fail for this periodic one-atom structure above. Since it seems to me that the the gradient provider RepresentationCalculatorGradientProvider ignores ghost contributions in the analytical gradient function grad_f https://github.com/cosmo-epfl/librascal/blob/db2e2445d34c196c94731249061740123f9fbc28/tests/test_calculator.hh#L972-L975 I think the gradient provider does something wrong, but I don't understand the function good enough to make a statement who does something wrong

[Luthaf]

@Luthaf do you still get an error when you increase the cutoff as you suggested in your first response?

No, I get a different error though (#379), but nothing that looks related to gradients.