Spin Spin Coupling can result in "cannot reshape array" if not all atoms are included in the coupling

[oliver-s-lee]

The exception is raised in rhf.py here:

def make_pso(sscobj, mol, mo1, mo_coeff, mo_occ, nuc_pair=None):
    if nuc_pair is None: nuc_pair = sscobj.nuc_pair
    para = []
    nocc = numpy.count_nonzero(mo_occ> 0)
    nvir = numpy.count_nonzero(mo_occ==0)
    atm1lst = sorted(set([i for i,j in nuc_pair]))
    atm2lst = sorted(set([j for i,j in nuc_pair]))
    atm1dic = dict([(ia,k) for k,ia in enumerate(atm1lst)])
    atm2dic = dict([(ia,k) for k,ia in enumerate(atm2lst)])
    mo1 = mo1.reshape(len(atm1lst),3,nvir,nocc) # <---------------- Exception here
    h1 = make_h1_pso(mol, mo_coeff, mo_occ, atm1lst)
    h1 = numpy.asarray(h1).reshape(len(atm1lst),3,nvir,nocc)
    for i,j in nuc_pair:
        # PSO = -Tr(Im[h1_ov], Im[mo1_vo]) + cc = 2 * Tr(Im[h1_vo], Im[mo1_vo])
        e = numpy.einsum('xij,yij->xy', h1[atm1dic[i]], mo1[atm2dic[j]])
        para.append(e*4)  # *4 for +c.c. and double occupnacy
    return numpy.asarray(para) * nist.ALPHA**4

if atm1lst and atm2lst are of different lengths. This is fairly easy to achieve if we decide we don't want to include all atoms in our coupling calculation (because we want to exclude far-away pairs, for example).

To illustrate:

from pyscf import gto, scf, dft
from pyscf.prop import ssc
# Pyridine
mol = gto.M(atom='''
            C         -2.41089        0.65767       -0.00003
            H         -1.27347       -2.54099       -0.00007
            C         -0.03591        0.62028       -0.00005
            C         -1.25643       -1.45879       -0.00006
            C         -2.45359       -0.73924       -0.00005
            H         -3.40367       -1.25766       -0.00005
            N         -1.21296        1.30258       -0.00003
            H         -3.32963        1.22882       -0.00002
            C         -0.03721       -0.77729       -0.00006
            H          0.89607       -1.32536       -0.00007
            H          0.90034        1.16223       -0.00005''',
            basis='def2-svp')
mf = dft.RKS(mol).set(xc='b3lyp').run()

coupling = ssc.RKS(mf)
coupling.with_fcsd = True
coupling.with_fc = True
# Pairs within 3.5 A of each other
coupling.nuc_pair = [(1, 0), (2, 0), (2, 1), (3, 0), (3, 1), (3, 2), (4, 0), (4, 1), (4, 2), (4, 3), (5, 0), (5, 1), (5, 3), (5, 4), (7, 0), (7, 2), (7, 3), (7, 4), (7, 5), (8, 0), (8, 1), (8, 2), (8, 3), (8, 4), (8, 5), (9, 1), (9, 2), (9, 3), (9, 4), (9, 8), (10, 0), (10, 2), (10, 3), (10, 8), (10, 9)]
coupling.kernel()
# ValueError: cannot reshape array of size 44352 into shape (9,3,88,21)

I'm struggling to fully follow the logic of this module, but it looks like mo1 = mo1.reshape(len(atm1lst),3,nvir,nocc) is using the wrong atom list? Changing to mo1 = mo1.reshape(len(atm2lst),3,nvir,nocc) seems to resolve the problem.

[oliver-s-lee]

Possible fix in #15