Using librascal for thousands of structures of 250 atoms each

[wilsonbill]

I generated 3000 structures with 250 atoms in each structure with corresponding energy and force data to use for training a kernel ridge regression model. I selected 1000 "sparse points" from the training set to create the kernel matrix using SOAP features. By "sparse points" I mean feature vectors of sites in structures. How many "sparse points" or "inducing points" are used to create the GAP_2017_6_17_60_4_3_56_165 potential mentioned in [1]? [1]Deringer, V. L., Bernstein, N., Bartók, A. P., Cliffe, M. J., Kerber, R. N., Marbella, L. E., ... & Csányi, G. (2018). Realistic atomistic structure of amorphous silicon from machine-learning-driven molecular dynamics. The journal of physical chemistry letters, 9(11), 2879-2885.

Does the following code look right for generating the kernel matrix and the kernel gradient matrix? I'm only using a single sparseFrame and 2 or 3 training frames since I run out of memory otherwise.

hypers = dict(soap_type="PowerSpectrum", interaction_cutoff=iCut, max_radial=radial, max_angular=angular, gaussian_sigma_constant=gSigma, gaussian_sigma_type="Constant", cutoff_smooth_width=sCut, normalize=False, radial_basis="GTO", compute_gradients=True, expansion_by_species_method='structure wise', ) soap = SphericalInvariants(hypers) managSparse= soap.transform(sparseFrames) soap1 = SphericalInvariants(hypers) manag_train= soap1.transform(trainFrames)

from rascal.models.sparse_points import SparsePoints spPts=SparsePoints(soap) spPts.extend(managSparse, [[siteIdx]]) kernel = Kernel(soap, name='GAP', zeta=zeta, target_type='Structure', kernel_type='Sparse') X_sparse=spPts KNM = kernel(managers_train, X_sparse) KNM_down = kernel(managers_train, X_sparse, grad=(True, False))

[agoscinski]

Assuming you mean this model https://github.com/libAtoms/silicon-testing-framework/blob/master/models/GAP/gp_iter6_sparse9k.xml it is 9000 sparse points as far as I understand quip model files.

I am not sure why there is such a massive RAM usage for just 3 train frames and 1 sparse frame, even with 250 atoms. How many atomic species does your system have? For such large structures it could be beneficial for RAM usage to use expansion_by_species_method='environment wise',, but I don't think that this is the problem here. The RAM usage will highly depend on your max_radial, max_angular and with gradients also on the cutoff.

It would be helpful if your code snippet contains the actual parameters values or how they were computed to get a feeling for your memory consumption. Maybe even a snippet of your dataset (for example the 4 frames you used) so we can analyze the memory consumption, since I don't see anything obviously wrong in your code.

[felixmusil]

Hello, the feature matrix for your whole training set is quite heavy (there are as many gradient features as there are atoms times 3). My advice it to compute the sparse points (without gradients) and then compute the representation (with gradients) and the kernel matrix one structure at a time. Then your memory footprint will be mostly tied to the size of the kernel matrix and not to the feature matrix. This is implemented in the feat/gaptools branch if you want a snippet.

[wilsonbill]

files.zip

Thank you for the responses. I found feat/gaptools. Can you suggest a particular python file that computes the kernel one structure at a time? If there is a workflow that has been used with such heavy loads that would be most helpful. The following computes the inducing points without gradients and the other structures with gradients which is what I understand Felix to be saying.

zeta= 1
iterat1=0
iterat=6
perIter1=3

radial=16
angular=16
interaction_cutoff=5.0;
gaussian_sigma_constant=0.4;
cutoff_smooth_width=0.5;

import numpy as np
import createFeat3
from rascal.models import Kernel, train_gap_model

import pymatgen.core.structure as pm

exCoord=np.load('coord29.npy')
exBox=np.load('box29.npy')
natom= int(exCoord.shape[1]/3)

j=28
siteIndex=236
if True:
        h2oCell= exBox[j].reshape([3,3])
        h2oLatt= pm.Lattice(h2oCell)
        h2oFrac= h2oLatt.get_fractional_coords(exCoord[j].reshape((3,natom)).T)
        #h2oFrac= h2oLatt.get_fractional_coords(exCoord[j].reshape((natom,3)))
        #hFrac = the fractional coordinants of h atoms
        #hFrac= h2oFrac[typat==1]
        rStruct= pm.Structure(pm.Lattice(h2oCell),['Si']*natom,h2oFrac, to_unit_cell=True)

from rascal.utils import from_dict, to_dict, CURFilter, dump_obj, load_obj, get_score, print_score
iCut=interaction_cutoff
gSigma=gaussian_sigma_constant
sCut= cutoff_smooth_width

hypers = dict(soap_type="PowerSpectrum",
                interaction_cutoff=iCut, 
                max_radial=radial, 
                max_angular=angular, 
                gaussian_sigma_constant=gSigma,
                gaussian_sigma_type="Constant",
                cutoff_smooth_width=sCut,
                normalize=False,
                radial_basis="GTO",
                compute_gradients=False,
                expansion_by_species_method='structure wise',
                )
from rascal.representations import SphericalInvariants

soap = SphericalInvariants(**hypers)
from pymatgen.io.ase import AseAtomsAdaptor as AAA
frames=[AAA.get_atoms(rStruct)]
managSparse= soap.transform(frames)
from rascal.models.sparse_points import SparsePoints

spPts=SparsePoints(soap)
spPts.extend(managSparse, [[siteIndex]])

structs=[]
for j in range(0,4):
        h2oCell= exBox[j].reshape([3,3])
        h2oLatt= pm.Lattice(h2oCell)
        h2oFrac= h2oLatt.get_fractional_coords(exCoord[j].reshape((3,natom)).T)
        #h2oFrac= h2oLatt.get_fractional_coords(exCoord[j].reshape((natom,3)))
        #hFrac = the fractional coordinants of h atoms
        #hFrac= h2oFrac[typat==1]
        rStruct= pm.Structure(pm.Lattice(h2oCell),['Si']*natom,h2oFrac, to_unit_cell=True)
        structs.append(rStruct)

hypers = dict(soap_type="PowerSpectrum",
                interaction_cutoff=iCut, 
                max_radial=radial, 
                max_angular=angular, 
                gaussian_sigma_constant=gSigma,
                gaussian_sigma_type="Constant",
                cutoff_smooth_width=sCut,
                normalize=False,
                radial_basis="GTO",
                compute_gradients=True,
                expansion_by_species_method='structure wise',
                )

frames= [AAA.get_atoms(struct) for struct in structs]
soap = SphericalInvariants(**hypers)
managers_train= soap.transform(frames)
from rascal.models import Kernel, train_gap_model
kernel = Kernel(soap, name='GAP', zeta=zeta, target_type='Structure', kernel_type='Sparse')

KNM = kernel(managers_train, spPts)
KNM_down = kernel(managers_train, spPts, grad=(True, False))

[felixmusil]

Sorry for the delay. You can have a look at this example notebook librascal/examples/zundel_i-PI.ipynb which shows how to fit a potential. Internally, it computes the kernel matrix progressively so you should not have problems with the memory anymore.

[agoscinski]

What @felixmusil mentioned is the function compute_KNM the part of

kernel = Kernel(soap, name='GAP', zeta=zeta, target_type='Structure', kernel_type='Sparse')

KNM = compute_KNM(tqdm(frames_train, leave=True, desc="Computing kernel matrix"), X_sparse, kernel, soap)

model = train_gap_model(kernel, frames_train, KNM, X_sparse, y_train, energy_baseline, 
                        grad_train=-f_train, lambdas=[1e-12, 1e-12], jitter=1e-13)

It could be that you still run out of memory, since your parameters and structures are large. In that case try to reduce your max_radial parameter to maybe 10 (the model you mention uses max_angular=12 max_radial=10). The parameters you use should be no problem for a HPC with 100GB RAM, but for a personal computer it might be at the limit.

[ceriottm]

you can use

hypers = get_optimal_radial_basis_hypers(hypers, frames,
expanded_max_radial=20)

on a smaller set of frames to compute an optimal radial basis following https://aip.scitation.org/doi/abs/10.1063/5.0057229

that should give you good results with a smaller nmax

On Mon, 1 Nov 2021 at 21:57, agoscinski @.***> wrote:

What @felixmusil https://github.com/felixmusil mentioned is the function compute_KNM the part of

kernel = Kernel(soap, name='GAP', zeta=zeta, target_type='Structure', kernel_type='Sparse')

KNM = compute_KNM(tqdm(frames_train, leave=True, desc="Computing kernel matrix"), X_sparse, kernel, soap)

model = train_gap_model(kernel, frames_train, KNM, X_sparse, y_train, energy_baseline, grad_train=-f_train, lambdas=[1e-12, 1e-12], jitter=1e-13)

It could be that you still run out of memory, since your parameters and structures are large. In that case try to reduce your max_radial parameter to maybe 10 (the model you mention uses max_angular=12 max_radial=10). The parameters you use should be no problem for a HPC with 100GB RAM, but for a personal computer it might be at the limit.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/cosmo-epfl/librascal/issues/390#issuecomment-956583500, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAIREZY3MTVYV3T7LR4JSUTUJ35KTANCNFSM5G57CI7Q .

[wilsonbill]

Thank you both again for the very helpful replies. The file librascal/examples/zundel_i-PI.ipynb looks like it would be exactly what I need. However I can't find functions like sparsify_environments() or calculate_representation() even though I am looking at branch feat/gaptools.

I see the following line in gp_iter6_sparse9k.xml which has many of the SOAP parameters. \soap l_max=12 n_max=10 atom_sigma=0.5 zeta=4 cutoff=5.0 cutoff_transition_width=1.0 central_weight=1.0 n_sparse=9000...

[wilsonbill]

Thank you Michele for the pointer to the function for choosing parameters.

I found sparsify_environments in librascal/bindings/rascal/models/gaptools.py. I was using the github search which doesn't seem to seach nonmaster branches.

Creating the kernel matrix is embarassing parallelizable so I've made some code using singularity, mpi and PBS batch scripts to run it on an HPC. I will see if changing l_max=12 n_max=10 reduces the load so that an HPC isn't necessary.