Concept check - gradients of SOAP

[g7dhaliwal]

Hi,

I have a fundamental question regarding gradients of SOAP descriptor. I used GAP to obtain SOAP descriptors for graphene dataset (as provided at libatoms.org) using following code. Corresponding to a central atom 'i', sum of gradients of all its neighbouring atoms 'j' comes out to be zero. Is it expected that gradients of central atom 'i', summed over all its neighbouring atoms 'j', to be zero?

The code below ouputs the maximum value of gradients along x,y,z, which comes out to be zero. Since SOAP for a central atom 'i' is calculated by summing over contributions from all the neighbouring 'j' atoms, I am under the impression that gradient should be non-zero.

Am I wrong in computing gradients for central atoms 'i' by summing over all the neighbouring atoms 'j'?
I will appreciate of your help.

Code

desc = descriptors.Descriptor("soap cutoff=4.2 l_max=4 n_max=4 atom_sigma=0.5 ")
trainInputConfig.set_cutoff(desc.cutoff()+wCut)
trainInputConfig.calc_connect()
d = desc.calc(trainInputConfig, grad = True)
gradIndex = d['grad_index_0based']
gradSOAP = d['grad']

gradMatrix = np.zeros((numAtoms,numFeatures,3))
for i in range(gradIndex.shape[0]):
    iAtom = gradIndex[i,0]
    gradMatrix[iAtom] = gradMatrix[iAtom]+gradSOAP[i].T
print(np.max(gradMatrix,axis=(0,1)))

[g7dhaliwal]

Hi, A follow-up to the above query. I tested the above code for few other datasets (downloaded from libatoms.org) and the sum of gradients of SOAP descriptors is zero for them as well. Another thing I noticed that few of the j-indices are repeating for a particular central atom 'i'.

I will greatly appreciate if you can help me with understanding and rectifying this issue.

[gabor1]

hi!

I think I can explain this.

Suppose you have a periodic structure (all structures in QUIP are defined as periodic) with N atoms, and you ask for the SOAP descriptors for that atoms structures, you will get N vectors back, one for each atom i considered as the “central” atom for the purposes of computing the SOAP descriptor of its environment. Let’s suppose that there are M components of these SOAP vector, indexed by m

The gradients are returned as follows. Consider an atom i, we want to know the derivative of its m-th SOAP component with respect to the atomic positions of atom j. That’s what grad_descriptor_out gives you. So why are there multiple j-s? because they correspond to various PERIODIC IMAGES of j. QUIP does not use the minimum image convention, so all periodic images within cutoff are accounted for.

What you’ve done in your script is chosen a given j, and added up all the contributions to different i-s (i.e. different centers). It’s sort of a strange thing to do, not corresponding to any physical thing that I can think of. I can’t immediately see why it would be zero, but it doesn’t especially bother me.

Gabor

On 10 May 2019, at 18:57, g7dhaliwal notifications@github.com wrote:

Hi,

A follow-up to the above query. I tested the above code for few other datasets (downloaded from libatoms.org) and the sum of gradients of SOAP descriptors is zero for them as well. Another thing I noticed that few of the j-indices are repeating for a particular central atom 'i'.

I am certain about missing some fundamental concept here and will greatly appreciate about your guidance.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

-- Gábor

Gábor Csányi Professor of Molecular Modelling Engineering Laboratory Pembroke College University of Cambridge

Pembroke College supports CARA. A Lifeline to Academics at Risk. http://www.cara.ngo/

[g7dhaliwal]

Hi Gabor,

Thanks a lot for the explanation regarding periodic images.

In the code, I am actually summing up contributions from all the j- atoms for a particular central atom. I am assuming that each row of d['grad_index_0based] has form as ( i, j). Below is the code with comments

#initialize tensor with all the zeros - to store gradient values
gradMatrix = np.zeros((N, N,M,3))
for i in range(gradIndex.shape[0]):
    iAtom = gradIndex[i,0] # Index for i atom - central atom
    jAtom = gradIndex[i,1] # Index for j atom 
    # the line below will ensure sum over all the periodic images of a particular 'j'
    gradMatrix[iAtom, jAtom, :, :] = gradMatrix[iAtom,jAtom,:,:]+gradSOAP[i].T

#Sum over all the neighbouring atoms j
gradMatrix = np.sum(gradMatrix,axis=1)

Gurjot

[gabor1]

Ok, so it's the gradient of a particular soap element wrt other neighboring atoms, summed over all atoms. Why is this an interesting quantity?

-- Gábor

On 10 May 2019, at 21:35, g7dhaliwal notifications@github.com wrote:

Hi Gabor,

Thanks a lot for the explanation regarding multiple images.

In the code, I am actually summing up contributions from all the j- atoms for a particular central atom. I am assuming that each row of d['grad_index_0based] has form as ( i, j). Below is the code with comments

initialize tensor with all the zeros - to store gradient values

gradMatrix = np.zeros((N, N,M,3)) for i in range(gradIndex.shape[0]): iAtom = gradIndex[i,0] # Index for i atom - central atom jAtom = gradIndex[i,1] # Index for j atom

the line below will ensure sum over all the periodic images of a particular 'j'

gradMatrix[iAtom, jAtom, :, :] = gradMatrix[iAtom,jAtom,:,:]+gradSOAP[i].T

Sum over all the neighbouring atoms j

gradMatrix = np.sum(gradMatrix,axis=1)

Gurjot

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

[gabor1]

I guess this can't generally be true. If we imagine just one neighbor existing, then that would imply that the gradient of soap wrt that neighbor is zero, which in general is not true for l=0 components.

-- Gábor

On 10 May 2019, at 21:35, g7dhaliwal notifications@github.com wrote:

Hi Gabor,

Thanks a lot for the explanation regarding multiple images.

In the code, I am actually summing up contributions from all the j- atoms for a particular central atom. I am assuming that each row of d['grad_index_0based] has form as ( i, j). Below is the code with comments

initialize tensor with all the zeros - to store gradient values

gradMatrix = np.zeros((N, N,M,3)) for i in range(gradIndex.shape[0]): iAtom = gradIndex[i,0] # Index for i atom - central atom jAtom = gradIndex[i,1] # Index for j atom

the line below will ensure sum over all the periodic images of a particular 'j'

gradMatrix[iAtom, jAtom, :, :] = gradMatrix[iAtom,jAtom,:,:]+gradSOAP[i].T

Sum over all the neighbouring atoms j

gradMatrix = np.sum(gradMatrix,axis=1)

Gurjot

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

[g7dhaliwal]

Please correct me but I think this quantity is equal to total gradient of soap element (for a central atom i) with respect to x, y, z direction.

[gabor1]

Total gradient withy respect to what?

You are adding up quantities that are derivarives with respect to different things (when you sum over j). That's not the gradient of anything with respect to anything I don't think.

-- Gábor

On 10 May 2019, at 21:50, g7dhaliwal notifications@github.com wrote:

Please correct me but I think this quantity is equal to total gradient of soap element (for a central atom i) with respect to x, y, z direction.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

[g7dhaliwal]

ok. I am seeking gradient of a soap element for atom 'i' w.r.t 'x', 'y', 'z'. Going with GAP output, my first thought was to sum it up along 'j' atoms, which is terrible idea.

With the current implementation of GAP, is it possible for me to obtain total gradient of soap element w.r.t 'x', 'y' 'z' directions?

[gabor1]

X,y,z of which atom? i? Ie the same around which that soap is calculated? That would be obtained by summing up the [i, i] elements of your matrix.

Why are you interested in this vector?

-- Gábor

On 10 May 2019, at 22:08, g7dhaliwal notifications@github.com wrote:

ok. I am seeking gradient of a soap element for atom 'i' w.r.t 'x', 'y', 'z'. Going with GAP output, my first thought was to sum it up along 'j' atoms, which is terrible idea.

With the current implementation of GAP, is it possible for me to obtain gradient of soap element w.r.t 'x', 'y' 'z' directions?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

[gabor1]

Do you agree that to get the forces on a given atom i, you should be doing the derivative of anything and everything with respect to positions of atom i, and not with respect to any other atomic position ?

For a linear model, the total energy is the sum of a vector of coefficients c, dotted into the soap vector, and this is evaluated for every atom. So the kinds of terms you need are the derivative of p(i) with respect to i (where I denote with p(.) the soap vector for a given environment), and the derivatives of p(j) with respect to positions of i, for neighbours j of i.

Also, the above derivatives aren't just added together, they also need to be multiplied by the linear model coefficients before adding them... where is this in your code?

-- Gábor

On 11 May 2019, at 01:34, g7dhaliwal notifications@github.com wrote:

For a linear model as discussed in this paper https://doi.org/10.1016/j.jcp.2014.12.018, I was getting the sum as posted in the above code

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

[g7dhaliwal]

Yes, I agree. In my local copy of the code, I have been considering linear coefficients. My mistake was summing up all the i atoms.

Thanks a lot for your valuable comments and constructive feedback. I appreciate greatly of your time.