Numpy operations on sparsed derivatives about dscribe

[lauri-codes]

Hi @VondrakMar!

We use the sparse library for the multi-dimensional sparse arrays. You can take a look at their API to see what kind of operations are supported: I do see at least einsum and moveaxis there, but I have not tried myself.

Another option is to convert to other formats for the operations. Converting to numpy arrays may only be possible for parts of the data, but if you can work with 2D slices, then you can also use scipy.sparse to perform certain tasks.

There are some additional details in our docs.

[VondrakMar]

Hey @lauri-codes,
really sorry for a late reply.

I rewrote my code with a sparse library. However, are derivatives transformed to the sparse output in one go AFTER all of them are computed?

I mean, if my derivatives of 1 system are filling my memory even before finishing their computation, will they fill the memory the same way even with sparse = True?

(I went very briefly through the code, and I think this is what is happening? But I am really not sure...)

Thank you again for your answer,
Best,
Martin

[lauri-codes]

Hi @VondrakMar:

Good question: The implementation is such that at most we keep the dense form for one system at a time per process. So when using multiple systems, each system gets translated to the sparse form when stored to the final result before moving on to the next system.

It would be possible to use the sparse form throughout the process, but this would be quite tedious to implement, as the sparse arrays are not a simple drop-in replacement on the C++ side. I guess you are hitting the limits of your system memory: how big systems are you dealing with and how many n_jobs are you using?

[VondrakMar]

Hey @lauri-codes
I have ~200 atoms. I am using 144 GB or 240 GB of RAM (both are too small).
I am calculating only 1 system, so n_job = 1.

However, I am pretty sure that it should be possible to fit sparsed derivatives into the memory.

Now I think that derivatives_single is what I should use.

Just calculate neighbor list with cutoff extended by get_cutoff_padding and then use this method for calculating only non-zero entries...

However, thank you for your answer,
Best,
Martin

[VondrakMar]

Hey @lauri-codes,
I am not sure, if I should open a new issue, or just write it here, but...

I did change my code to use derivatives_single, however for some reason that I don't understand, attach = True is not working as it should. I am attaching a code.

import numpy as np
from ase.build import bulk,make_supercell
from dscribe.descriptors import SOAP


a2 = bulk('Cu', 'fcc', a=3.6, orthorhombic=True)
newMol = make_supercell(a2,[[5,0,0],[0,5,0],[0,0,5]])
soap = SOAP(
    n_max = 4,
    l_max = 3,
    r_cut = 2,
    sigma = 1,
    periodic=True,
    species=["Cu"]
derW,descW = soap.derivatives(newMol,attach=True)
der,desc = soap.derivatives_single(newMol,centers=[5],indices=[4],attach=True)
print(der[0][0][0]/derW[5][4][0]) # 1s
der,desc = soap.derivatives_single(newMol,centers=[5],indices=[5],attach=True)
print(der[0][0][0]/derW[5][5][0]) # not 1s
derW,descW = soap.derivatives(newMol,attach=False)
der,desc = soap.derivatives_single(newMol,centers=[5],indices=[5],attach=True)
print(der[0][0][0]/derW[5][5][0]) # 1s again

derW,descW = soap.derivatives(newMol,attach=True)
der,desc = soap.derivatives_single(newMol,centers=[0],indices=[0],attach=True)
print(der[0][0][0]/derW[0][0][0]) # 1s again

If indices and centers are the same (and are not 0s), it is calculated as if attach = False. Do you have a idea how it can be repair, or what I am doing wrong? Also, I think this is happening only for periodic systems.
Thank you,
Best,
Martin