[BUG] Fails evaluating density of open shell systems

[marco-2023]

Describe the bug

The function evaluate_density from module gbasis.evals.density fails to evaluate $\rho(r)$ in a grid of points for open-shell systems.

To Reproduce

# Load open shell system from fchk
#--------------------------------------------------
from iodata import load_one

# fchk of oxygen atom fetched from :
# https://github.com/theochem/chemtools/blob/master/chemtools/data/atom_08_O_N08_M3_ub3lyp_ccpvtz_g09.fchk
mol = load_one("atom_08_O_N08_M3_ub3lyp_ccpvtz_g09.fchk")

# Importing grid dependencies, and creating grid
#----------------------------------------------------------------------
# Make Becke-Lebedev molecular grid (using preset grid)
from grid.becke import BeckeWeights
from grid.molgrid import MolGrid
from grid.onedgrid import GaussChebyshev
from grid.rtransform import BeckeRTransform

oned = GaussChebyshev(100)
rgrid = BeckeRTransform(1e-4, 1.5).transform_1d_grid(oned)
grid = MolGrid.from_preset(mol.atnums, mol.atcoords, rgrid, "coarse", BeckeWeights())

# Evaluate the electron density for the grid points
#----------------------------------------------------------------------

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

one_rdm = np.dot(mol.mo.coeffs * mol.mo.occs, mol.mo.coeffs.T)
basis = from_iodata(mol)
density = evaluate_density(one_rdm, basis[0], grid.points, coord_type=basis[1])

Expected behaviour

density should be a numpy.ndarray with the values of the electron density for each point of the grid.

Screenshots

System Information:

• OS: "openSUSE Leap 15.5"
• Python version: 3.11.4
• NumPy version: 1.25.2
• SciPy version: 1.11.1

Additional context

[leila-pujal]

Hi @marco-2023, it is great you are using Gbasis. I saw your issue this morning and I tested your code because so far I haven't had problems with Gbasis but it is true it has not been in active development so it could very well be that a bug like this could have gone undetected after the merging of some PR. Anyway, I tried your code and as far as I tested it worked fine for me. Below I added the sniped I used and the outputs so you can check if this is the same thing you tried with your code. I initially copy the code you added to the issue and the only thing missing was importing Numpy. I am checking the package versions you are using and I have older versions installed. This could be a version issue where something allowed in older version of Numpy for instance is not allowed in newer ones and it is causing the output of Nans. If you could try to downgrade to my versions and see if the issue is fix that would be great and it would mean it is most likely a version problem. Let me know what you think.

• Python version: 3.7.12
• NumPy version: 1.21.6
• SciPy version: 1.7.3

# Load open shell system from fchk
#--------------------------------------------------
import numpy as np  # Added this in your original code in the issue
from iodata import load_one

# fchk of oxygen atom fetched from :
# https://github.com/theochem/chemtools/blob/master/chemtools/data/atom_08_O_N08_M3_ub3lyp_ccpvtz_g09.fchk
mol = load_one("atom_08_O_N08_M3_ub3lyp_ccpvtz_g09.fchk")

# Importing grid dependencies, and creating grid
#----------------------------------------------------------------------
# Make Becke-Lebedev molecular grid (using preset grid)
from grid.becke import BeckeWeights
from grid.molgrid import MolGrid
from grid.onedgrid import GaussChebyshev
from grid.rtransform import BeckeRTransform

oned = GaussChebyshev(100)
rgrid = BeckeRTransform(1e-4, 1.5).transform_1d_grid(oned)
grid = MolGrid.from_preset(mol.atnums, mol.atcoords, rgrid, "coarse", BeckeWeights())

# Evaluate the electron density for the grid points
#----------------------------------------------------------------------

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

one_rdm = np.dot(mol.mo.coeffs * mol.mo.occs, mol.mo.coeffs.T)
#one_rdm = np.dot(mol.mo.coeffsa * mol.mo.occsa, mol.mo.coeffsa.T)
basis = from_iodata(mol)
density = evaluate_density(one_rdm, basis[0], grid.points, coord_type=basis[1])
print('density', density)
print('number electrons', grid.integrate(density))

Output

(5324,)
density [299.08773677 299.08773677 299.08773741 ...   0.           0.
   0.        ]
number electrons 8.000000000788521

[PaulWAyers]

Could it be related to issue #129 ???

[leila-pujal]

@PaulWAyers, @marco-2023, I just took another look and yes @PaulWAyers, it looks it is because of issue #129. I test Marco's code with both, old and new versions of Scipy and the normalization is messed up as it was pointed out in #129:

gbasis/contractions.py

gbasis/gbasis/contractions.py

Lines 409 to 418 in 2eb9e49

	"""
	exponents = self.exps[np.newaxis, :]
	angmom_components_cart = self.angmom_components_cart[:, :, np.newaxis]
	return (
	(2 * exponents / np.pi) ** (3 / 4)
	* ((4 * exponents) ** (self.angmom / 2))
	/ np.sqrt(np.prod(factorial2(2 * angmom_components_cart - 1), axis=1))
	)

return is [[inf inf inf inf inf inf inf]] for the new scipy because factorial2 returning 0 is in the denominator. It returns [[981.89875275 236.62701252 77.87758373 30.1343725 12.90544497 5.90358204 2.80266636]] for the old Scipy.
Then in gbasis/evals/_deriv.py norm is this array of inf and generates Nan in the eval_deriv_contractions when it computes the tensordot. This is not the case with the values array obtained with the old scipy.

gbasis/gbasis/evals/_deriv.py

Lines 140 to 142 in 2eb9e49

	norm = norm.T[:, :, np.newaxis]
	return np.tensordot(prim_coeffs, norm * zeroth_part * deriv_part, (0, 0))

This problem is avoidable by downgrading the versions for now but I guess either we should update/check gbasis for newer versions of its dependencies or update the README file to constraint the available versions to use.

[PaulWAyers]

I'm going to close this issue referencing the underlying issue, #129