[BUG] Unexpected overlap integral with xtb's hydrogen basis set

[RaphaelRobidas]

Describe the bug
I have been playing with xtb's basis set and ran some simple calculations. However, the overlap integrals calculated with gbasis do not match those of dxtb along with "manual" calculations.

To Reproduce
I focused on the simplest case of hydrogen, which has two s shells in xtb's basis set:

BASIS "ao basis" SPHERICAL PRINT
#
H    s
      7.611997207059968      0.18536103626168826
      1.392901705880201      0.23771678223106302
      0.3869633462504825      0.1863220559732388
      0.12842965592697322      0.044589693726254674
H    s
      10.256286070314902      -1.3186544678254288
      0.622796532587539      1.6038777007695835
      0.23910076678539138      0.6013230101772801
      7.611997207059968      -0.9809043198633103
      1.392901705880201      -1.2579635035340504
      0.3869633462504825      -0.9859898999050776
      0.12842965592697322      -0.23596233641959458
END

According to dxtb, the overlap matrix of a system containing only one hydrogen is:

tensor([[1.0000e+00, 8.5040e-10],
        [8.5040e-10, 1.0000e+00]], dtype=torch.float64)

So both shells appear to have null net overlap, within numerical error. Now doing the same with gbasis and custom code:

import numpy as np

from gbasis.parsers import parse_nwchem, make_contractions
from gbasis.integrals.overlap import overlap_integral

bs = parse_nwchem("H_xtb.bs")
basis = make_contractions(
    bs, ["H"], np.array([[0.0, 0.0, 0.0]]), coord_types="spherical"
)

overlap = overlap_integral(basis)
print(f"gbasis overlap: {overlap[0,1]}")


def gaussian_prod(coeff1, exp1, R1, coeff2, exp2, R2):
    p = exp1 + exp2
    q = R1 * exp1 + R2 * exp2

    new_coeff = coeff1 * coeff2 * np.exp(-(exp1 * exp2 * np.linalg.norm(R1 - R2)) / p)
    return new_coeff, p, q / p


S = 0
for c1, exp1 in zip(basis[0]._coeffs, basis[0]._exps):
    for c2, exp2 in zip(basis[1]._coeffs, basis[1]._exps):
        new_coeff, p, q = gaussian_prod(
            c1[0],  # The coefficients are 1-element arrays
            exp1,
            basis[0]._coord,
            c2[0],
            exp2,
            basis[1]._coord,
        )
        S += new_coeff * (np.pi / p) ** (3 / 2)
print(f"Manual overlap: {S}")

with H_xtb.bs being the basis set specifications above. This code outputs the following:

gbasis overlap: -0.6910726848588858
Manual overlap: 8.504019710642297e-10

Expected behaviour
I would expect the overlap integral from gbasis to be close to zero as well. Maybe I am using gbasis wrong somehow, but this seems very odd.

System Information:

• Python version: 3.10.12
• NumPy version: 1.26.4
• SciPy version: 1.11.4
• gbasis version: 3e063ad

[leila-pujal]

Hi @RaphaelRobidas, thanks for opening the issue and your detailed explanation. The problem here I think is that, for your basis set H_xtb.bs, the normalization for each primitive in the contraction is already included in the primitive's coefficients. I could get your overlap result back with Gbasis when I change the return value here to be:

return np.ones((1,self.exps.shape[0]))

In our docs, we state Gbasis works with normalized primitives but it looks we may be able to add an option so the user can specify if primitives have already normalized coefficients. Maybe @PaulWAyers has more thoughts on this. @RaphaelRobidas let me know if you have any questions about this.

[RaphaelRobidas]

Hello @leila-pujal, thanks for looking into this! This makes a lot of sense, I hadn't realized that the coefficients were already normalized, my bad! In any case, knowing this, this issue is resolved for me.

[PaulWAyers]

@leila-pujal it would be good to add the feature to support unnormalized primitives. We, like most software, had opted to normalize the primitives as when the primitives are unnormalized one often ends up with contraction coefficients that span many orders of magnitude, with the consequent potential for numerical ill-conditioning. (I'm not sure how big an issue this is in practice; we didn't check, but it seemed to be a potential issue.)

It's sensible to add this as a feature, though I would argue that gbasis itself shouldn't really change, just that we should be able to pre-process unnormalized primitives to normalized primitives and then undo the result at the end. So basically we'd add a utility module to help this.

We should probably make sure this is documented more clearly, too, though the fact this issue will be searchable in the repository will help on that.

[leila-pujal]

@PaulWAyers, I think that is a good plan. I believe that option could easily be added to contractions.py module. @RaphaelRobidas, glad this solved your problem! Thanks for closing the issue. Please, feel free to open new issues if you have any other doubt that comes up from using Gbasis.

[PaulWAyers]

@leila-pujal can you make a new issue/feature-request for contractions.py for this?

[leila-pujal]

Yes, no problem, I will open an issue and add some info on how it should be implemented.