two.py

def twoint_kernel(
        RA, RB, RC, RD, lA, lB, lC, lD,
        aA, aB, aC, aD, cA, cB, cC, cD
    ):
    """ Two-electron integral kernel

        Calculated with primitive gaussian type orbitals.

        This function implements equation 13 from:
        DOI: 10.1002/ejtc.36

        Arguments:
        mu --
        m --
        nu --
        n --
        sigma --
        s --
        lamda --
        l --

        Returns:
        two-electron integral
    """

    # coordinates and distances
    dRAB = RA - RB
    dRCD = RC - RD

    # intermediate variables and definitions
    S1 = aA + aB
    S2 = aC + aD
    S4 = S1 + S2
    S = S1 * S2 / S4

    P = (aA*RA + aB*RB)/S1
    Q = (aC*RC + aD*RD)/S2
    dPQ = P - Q
    t = S * dPQ.dot(dPQ)
    LL = 2.0*(S*PI_inv)**(1.0/2.0)

    if lA == 0 and lB == 0 and lC == 0 and lD == 0:
        # primitive <ss|ss> integral
        g0000 = integrated_incomplete_gamma(0, t)

        sab = ss_primitive_overlap(aA, aB, cA, cB, dRAB.dot(dRAB))
        scd = ss_primitive_overlap(aC, aD, cC, cD, dRCD.dot(dRCD))
        return LL * sab * scd * g0000

    raise NotImplementedError("Angular momenta beyond L=0 is not implemented")
    return 0.0