Primitive integrals/contraction coefficients

[maxscheurer]

Hi, I'm currently experimenting with some integrals for implementing GOSTSHYP in pyscf. For that, I need to build a custom "fakemol" and basis set, so I've been using pyscf's fakemol_for_charges to do the setup. The reason for this is that these special basis functions are not normalized.

Along the way, I've been wondering which integrals libcint actually computes under the hood, so I've tried to "manually" compute the overlap of a primitive s-function with itself, with exponent of one, contraction coefficient of one, centered at the origin.

I was wondering why the result from libcint is not identical to what I get with sympy? What am I forgetting/missing here?

import numpy as np
from pyscf import gto
from sympy import symbols, exp, integrate, oo, N


def fakemol_with_exponents(coords, exponents):
    nbas = coords.shape[0]

    fakeatm = np.zeros((nbas, gto.mole.ATM_SLOTS), dtype=np.int32)
    fakebas = np.zeros((nbas, gto.mole.BAS_SLOTS), dtype=np.int32)
    fakeenv = [0] * gto.mole.PTR_ENV_START
    ptr = gto.mole.PTR_ENV_START
    fakeatm[:, gto.mole.PTR_COORD] = np.arange(ptr, ptr + nbas * 3, 3)
    fakeenv.append(coords.ravel())
    ptr += nbas * 3
    fakebas[:, gto.mole.ATOM_OF] = np.arange(nbas)
    fakebas[:, gto.mole.NPRIM_OF] = 1
    fakebas[:, gto.mole.NCTR_OF] = 1
    # approximate point charge with gaussian distribution exp(-expnt*r^2)
    fakebas[:, gto.mole.PTR_EXP] = ptr + np.arange(nbas) * 2
    fakebas[:, gto.mole.PTR_COEFF] = ptr + np.arange(nbas) * 2 + 1
    coeff = np.ones_like(exponents)  # / (2 * np.sqrt(np.pi) * gaussian_int(2, expnt))
    fakeenv.append(np.vstack((exponents, coeff)).T.ravel())

    fakemol = gto.Mole()
    fakemol._atm = fakeatm
    fakemol._bas = fakebas
    fakemol._env = np.hstack(fakeenv)
    fakemol._built = True
    return fakemol


gmol2 = fakemol_with_exponents(
    np.array(
        [
            [0, 0, 0],
        ]
    ),
    np.array([1.0]),
)
print(gmol2.intor("int1e_ovlp")[0, 0])

x, y, z = symbols("x, y, z")
ref = integrate(
    exp(-1.0 * (x ** 2 + y ** 2 + z ** 2)) ** 2, (x, -oo, oo), (y, -oo, oo), (z, -oo, oo)
)
print(N(ref))

Output:

0.15666426716443752
1.96870124321530

Thanks for helping with this (probably stupid) question 😁

[sunqm]

1. Overlap is defined as the product of two Gaussians. You should compare to

ref = integrate(
    exp(-2.0 * (x ** 2 + y ** 2 + z ** 2)) ** 2, (x, -oo, oo), (y, -oo, oo), (z, -oo, oo)
)

2. for s-type orbitals, int1e_ovlp includes the normalization factor for angular part 1/4pi

[maxscheurer]

Thanks, the 1/(4pi) did the trick 👍 Regarding 1., I already had a product of two Gaussians, note the ** 2 after the exp function. I'll leave this issue open for follow-up questions during my implementation if that's okay.