[BUG] Negative Electron Density Values

[FarnazH]

Describe the bug I was integrating density^(1/3) when I noticed that the electron density computed by evaluate_density gives a negative, but very small value, for faraway points (e.g. -2.77765871e-128).

To Reproduce Helium FCHK file: helium.fchk.tar.gz

import numpy as np

from iodata import load_one

from grid.molgrid import MolGrid
from grid.becke import BeckeWeights
from grid.rtransform import BeckeTF
from grid.onedgrid import GaussChebyshevType2

from gbasis.wrappers import from_iodata
from gbasis.evals.density import evaluate_density

# load molecule
mol = load_one("helium.fchk")
dm = mol.one_rdms['scf']
print("molecular charge   = ", mol.charge)
print("atomic numbers     = ", mol.atnums)
print("atomic coordinates = ", mol.atcoords)

# make grid
oned = GaussChebyshevType2(150)
rgrid = BeckeTF(1.0e-6, 1.5).transform_1d_grid(oned)
grid = MolGrid.horton_molgrid(mol.atcoords, mol.atnums, rgrid, 350, BeckeWeights(order=3))

# evaluate density
basis, coord_types = from_iodata(mol)
dens = evaluate_density(dm, basis, grid.points, coord_type=coord_types)
print("min & max density = ", np.min(dens), np.max(dens))

index = np.where(dens < 0.0)
print("density     at points with negative densiry = ", dens[index])
print("coordinates of points with negative density = ", grid.points[index])

The output is:

molecular charge   =  0.0
atomic numbers     =  [2]
atomic coordinates =  [[0. 0. 0.]]
min & max density =  -2.77765870781763e-128 2.930734188084798
density     at points with negative densiry =  [-2.77765871e-128]
coordinates of points with negative density =  [[ 10.15753086 -15.57090172 -10.15753086]]

Expected behaviour The electron density is expected to be positive at every point.

System Information:

• OS: macOS
• Python version: 3.6.10
• NumPy version: 1.18.1
• SciPy version: 1.4.1

Additional context I am using the snippet below as a workaround, until this get fixed. dens[np.where((dens < 0.0) & (abs(dens) < 1.0e-15))] = 0.0

[PaulWAyers]

I think this is a round-off error in evaluating the density or a precision error in the 1-electron density matrix. It is sensible that when a quantity that is returned is expected to be positve, gbasis would test for negative elements, and then set them to zero (if they are negligible, perhaps using the criterion @FarnazH uses) and otherwise return an error.

Probably it is inevitable that some density values evaluate to something negative....but no reason to return those values and crash the codes of users.

[PaulWAyers]

A note the positive-definite kinetic energy density also shouldn't ever return a negative value. gbasis.density.evaluate_posdef_kinetic_energy_density

Probably the same approach is valid there.

[FarnazH]

@maximilianvz fixed this issue in PR #136 using the numpy.clip by providing 0.0 as the minimum value of array, i.e., dens.clip(min=0.0).

Now that I am thinking, the only caveat is that if there are negative values with a large magnitude in the array, they are also set to zero. Whereas the recommended workaround in this issue, dens[np.where((dens < 0.0) & (abs(dens) < 1.0e-15))] = 0.0, make sure that the negative values have a small magnitude, when replacing them with zero.

I am not sure whether we should worry about this, but wanted to mention it. @PaulWAyers and @leila-pujal, feel free to let me know what you think.

[PaulWAyers]

@FarnazH I think that 1e-15 may be ambitious; I would be happy with 1e-8. (I can imagine cases where the electron density was evaluated from a density matrix that was read from a file with limited precision in the exponents.)

The other caveat is what I mentioned, which is that the same problem may occur for the positive definite k.e. density.

[FarnazH]

@PaulWAyers, sure, I understand that 1.0e-15 was ambitious. It is just what I had used in my specific script as a workaround. My main question was whether clipping (which does not check the magnitude of negative values) is a good enough solution. Any thoughts?

@maximilianvz clipped the output of gbasis.density.evaluate_posdef_kinetic_energy_density too. So, that is taken care of.

[PaulWAyers]

I'd like it to be vectorized....what about finding the minimum value of the array first, and if the minimum value is > -1.0e-8 we clip, and otherwise we raise an error?