def MakeNuclearSqDistanceIntegrals(self, Smh=True):
    'calculate <mu|(r-rA)^2|nu> integrals in main basis, for each individual atomic core.\n      Returns nAo x nAo x nAtoms array.'
    Args = []
    if Smh:
        Args.append('--orb-trafo=Smh')
    Outputs = []
    Outputs.append(('--save-rsqN', 'RSQ_N'))
    RsqN = self._InvokeBfint(Args, Outputs)[0]
    nOrb = int(((RsqN.shape[0] ** 0.5) + 0.5))
    assert ((nOrb ** 2) == RsqN.shape[0])
    assert (RsqN.shape[1] == len(self.Atoms))
    return RsqN.reshape(nOrb, nOrb, RsqN.shape[1])
calculate <mu|(r-rA)^2|nu> integrals in main basis, for each individual atomic core.
Returns nAo x nAo x nAtoms array.
def MakeNuclearSqDistanceIntegrals(self, Smh=True):
    'calculate <mu|(r-rA)^2|nu> integrals in main basis, for each individual atomic core.\n      Returns nAo x nAo x nAtoms array.'
    Args = []
    if Smh:
        Args.append('--orb-trafo=Smh')
    Outputs = []
    Outputs.append(('--save-rsqN', 'RSQ_N'))
    RsqN = self._InvokeBfint(Args, Outputs)[0]
    nOrb = int(((RsqN.shape[0] ** 0.5) + 0.5))
    assert ((nOrb ** 2) == RsqN.shape[0])
    assert (RsqN.shape[1] == len(self.Atoms))
    return RsqN.reshape(nOrb, nOrb, RsqN.shape[1])