def METHOD_NAME(mo_coeff):
    nao, nmo = mo_coeff.shape
    mo_a = mo_coeff[:nao//2]
    mo_b = mo_coeff[nao//2:]
    orbspin = getattr(mo_coeff, 'orbspin', None)
    if orbspin is None:
        eri  = with_df.ao2mo(mo_a, kpt, compact=False)
        eri += with_df.ao2mo(mo_b, kpt, compact=False)
        eri1 = with_df.ao2mo((mo_a,mo_a,mo_b,mo_b), kpt, compact=False)
        eri += eri1
        eri += eri1.T
        eri = eri.reshape([nmo]*4)
    else:
        # If GHF orbitals have orbspin labels, alpha and beta orbitals
        # occupy different columns. Here merging them into one set of
        # orbitals then zero out spin forbidden MO integrals
        mo = mo_a + mo_b
        eri  = with_df.ao2mo(mo, kpt, compact=False).reshape([nmo]*4)
        sym_forbid = (orbspin[:,None] != orbspin)
        eri[sym_forbid,:,:] = 0
        eri[:,:,sym_forbid] = 0
    return eri