SINGROUP / dscribe

DScribe is a python package for creating machine learning descriptors for atomistic systems.
https://singroup.github.io/dscribe/
Apache License 2.0
398 stars 89 forks source link

'nan' value occur in acsf descriptors #70

Closed zjujdj closed 2 years ago

zjujdj commented 3 years ago

hello, i calculated the acsf descriptors for the atoms in a protein using following codes(dscribe==1.1.0), and i was confused by 'nan' value could occur for some atoms, does anyone know why? thanks.

# codes

import torch
from rdkit import Chem
import numpy as np
import pickle
from ase import Atoms
with open(r'C:\Users\15766\Desktop\1A2K_complex.3450_19.71_0', 'rb') as f:
    ligand, pocket = pickle.load(f)
atom_ls = []
atom_ls.extend([atom.GetSymbol() for atom in ligand.GetAtoms()])
atom_ls.extend([atom.GetSymbol() for atom in pocket.GetAtoms()])
atom_positions = np.concatenate([ligand.GetConformer().GetPositions(), pocket.GetConformer().GetPositions()], axis=0)
mol = Atoms(symbols=atom_ls, positions=atom_positions)
from dscribe.descriptors import ACSF
acsf = ACSF(
    species=['C', 'O', 'N', 'S', 'P', 'F', 'Cl', 'Br', 'I'],
    rcut=12.0,
    g2_params=[[4.0, 3.17]],
    g4_params=[[0.1, 3.14, 1]],
)
res = acsf.create(mol)
res_th = torch.tensor(res, dtype=torch.float)
print(torch.any(torch.isnan(res_th)))
print(np.where(torch.isnan(res_th)))
res_th[1690]

# outputs

tensor(True)
(array([1690], dtype=int64), array([19], dtype=int64))
Out[92]: 
tensor([2.7536e+01, 1.7760e+00, 5.9737e+00, 7.7800e-01, 4.8037e+00, 4.7535e-07,
        0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        3.9842e+00,        nan, 1.3123e-02, 8.6276e-02, 2.6853e-02, 8.8496e-04,
        0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 0.0000e+00, 0.0000e+00])

# original file 1A2K_complex.3450_19.71_0

lauri-codes commented 3 years ago

Hi @zjujdj!

I think this is a case of some kind of a floating point overflow. Maybe at some point your G4 terms get a really large value for an exponent (since you have a very large system and large cutoff), which causes numerical problems. On my machine simply slightly decreasing rcut got rid of the problem. I can try and see if this particular case is fixable by increasing the precision on some of our internal variables.

zjujdj commented 3 years ago

Hi @lauri-codes : many thanks for your timely reply, i tried a smller rcut value of 8A, it did slove many problematic cases, but not all. my system are some big ppi complexes, and i guess 8A is a reasonable cutoff to determine atom-atom interaction. Anyway, thanks for your suggestions.

aneepi commented 3 years ago

Hello!

I have also encountered similar problems. I tried to describe local environments of atoms in complex nanoparticles. If I use g4 and g5 parameters of ACSF, I get infs and nans. By changing cutoff some infs change to nans or vice versa. This happens with just cutoff of 5 Å, which I think is relatively short distance, because most chemical bonds are within 2 to 3Å. Interestingly infs and nans do not occur on version 0.2.2 but at least 0.4.0 and 1.1.0 are both producing infs and nans regularly if systems are large.

zjujdj commented 3 years ago

@aneepi hi Thanks for your information! maybe i can try the former version of this package. here i tired other similar package such as torchani, i never get such error, but torchani is very memory-consuming when deal with large system. Overall, dscribe is better.

aneepi commented 3 years ago

Hello @lauri-codes

I would also want to add previous comment that I sent. On version 0.2.2 the corresponding g4 and g5 descriptions are very small. They are in my case between 1e-3 and 1e-6. For me this feels like some weird memory handling issue in C or Fortran bindings in newer versions. Especially because g1, g2 and g3 did not have this issue. Furthermore, g4 and g5 work well with small systems on newer versions too. It seems that big systems break something.

Nice if my comment helped you @zjujdj : ) Seems like keeping different older versions in virtual environments was useful.

lauri-codes commented 2 years ago

Hi @aneepi, @zjujdj,

Please try out the new version v1.2.0: it uses larger precision for ACSF internally and should get rid of the nan values in most cases.

zjujdj commented 2 years ago

@lauri-codes , Hi, thanks very much! i'll try the new version later!

lauri-codes commented 2 years ago

Ok, great! Just re-open this issue if you still experience any problems.