How to use the pyscf.lib.direct_sum function in pyscf

l2new += m3 * .5
    l1new += numpy.einsum('ijab,jb->ia', m3, t1) * 2
    l1new -= numpy.einsum('ijba,jb->ia', m3, t1)
    saved_woooo = m3 = None
    #time1 = log.timer_debug1('lambda pass [%d:%d]'%(p0, p1), *time1)

    eia = lib.direct_sum('i-a->ia', mo_e_o, mo_e_v)
    l1new /= eia

#    l2new = l2new + l2new.transpose(1,0,3,2)
#    l2new /= lib.direct_sum('ia+jb->ijab', eia, eia)
#    l2new += l2
    for i in range(nocc):
        if i > 0:
            l2new[i,:i] += l2new[:i,i].transpose(0,2,1)
            l2new[i,:i] /= lib.direct_sum('a,jb->jab', eia[i], eia[:i])
            l2new[:i,i] = l2new[i,:i].transpose(0,2,1)
        l2new[i,i] = l2new[i,i] + l2new[i,i].T
        l2new[i,i] /= lib.direct_sum('a,b->ab', eia[i], eia[i])

    time0 = log.timer_debug1('update l1 l2', *time0)
    return l1new, l2new

imds.l1a_t += numpy.einsum('ijkabc,jbkc->ia', rw.conj(),
                               numpy.asarray(eris.OVOV)) / eia * .25
    imds.l1b_t += numpy.einsum('ijkabc,iakc->jb', rw.conj(),
                               numpy.asarray(eris.ovOV)) / eIA * .5
    wvd = r4(w * 2 + v) / d3
    l2_t  = numpy.einsum('ijkabc,iaec->jkbe', wvd, numpy.asarray(eris.get_ovVV()).conj())
    l2_t -= numpy.einsum('ijkabc,iakm->jmbc', wvd, numpy.asarray(eris.ovOO).conj())
    l2_t = l2_t + l2_t.transpose(1,0,3,2)
    l2_t += numpy.einsum('ijkabc,ai->jkbc', rw, fvo)
    imds.l2bb_t += l2_t.conj() / lib.direct_sum('ia+jb->ijab', eIA, eIA) * .5
    l2_t  = numpy.einsum('ijkabc,kceb->ijae', wvd, numpy.asarray(eris.get_OVVV()).conj())
    l2_t += numpy.einsum('ijkabc,kcea->ijeb', wvd, numpy.asarray(eris.get_OVvv()).conj())
    l2_t -= numpy.einsum('ijkabc,kcjm->imab', wvd, numpy.asarray(eris.OVOO).conj())
    l2_t -= numpy.einsum('ijkabc,kcim->mjab', wvd, numpy.asarray(eris.OVoo).conj())
    l2_t += numpy.einsum('ijkabc,bj->ikac', rw, fVO)
    imds.l2ab_t += l2_t.conj() / lib.direct_sum('ia+jb->ijab', eia, eIA) * .5

    return imds

B  -= numpy.einsum('mRnQ,RiQj->minj', tmp, RRiQj)
    tmp = numpy.einsum('mPnk,PQ->mQnk', RmPnk, fPQ)
    B  -= numpy.einsum('mQnk,Qikj->minj', tmp, RQikj)
    tmp = numpy.einsum('mknc,kl->mlnc', Rmknc, fPQ[:nocc,:nocc])
    B  += numpy.einsum('mlnc,licj->minj', tmp, Rlicj)
    tmp = numpy.einsum('mpna,pq->mqna', Rmpna, fPQ[:nmo,:nmo])
    B  -= numpy.einsum('mqna,qiaj->minj', tmp, Rqiaj)
    tmp = numpy.einsum('mknc,kP->mPnc', Rmknc, fPQ[:nocc])
    tmp1= numpy.einsum('mPnc,Picj->minj', tmp, RPicj)
    tmp = numpy.einsum('mcnb,cp->mpnb', Rmcnb, fPQ[nmo:,:nmo])
    tmp1+= numpy.einsum('mpnb,pibj->minj', tmp, Rpibj)
    B  -= tmp1 + tmp1.transpose(1,0,3,2)
    B   = B + B.transpose(2,3,0,1)
    B  += tau

    e_mn = lib.direct_sum('i+j->ij', mo_energy[:nocc], mo_energy[:nocc])
    tmp = numpy.einsum('mknl,kilj->minj', tminj, B)
    emp2_f12+= numpy.einsum('minj,minj', tmp, tminj) * 2
    emp2_f12-= numpy.einsum('minj,nimj', tmp, tminj)
    tmp = numpy.einsum('mknl,kilj->minj', tminj, X)
    emp2_f12-= numpy.einsum('mn,minj,minj', e_mn, tmp, tminj) * 2
    emp2_f12+= numpy.einsum('mn,minj,nimj', e_mn, tmp, tminj)
    return emp2_f12

def _gamma2_intermediates(mycc, t1, t2, l1, l2, eris=None):
    dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov = \
            gccsd_rdm._gamma2_intermediates(mycc, t1, t2, l1, l2)
    if eris is None: eris = mycc.ao2mo()

    nocc, nvir = t1.shape
    bcei = numpy.asarray(eris.ovvv).conj().transpose(3,2,1,0)
    majk = numpy.asarray(eris.ooov).conj().transpose(2,3,0,1)
    bcjk = numpy.asarray(eris.oovv).conj().transpose(2,3,0,1)

    mo_e = eris.mo_energy
    eia = mo_e[:nocc,None] - mo_e[nocc:]
    d3 = lib.direct_sum('ia+jb+kc->ijkabc', eia, eia, eia)

    t3c =(numpy.einsum('jkae,bcei->ijkabc', t2, bcei)
        - numpy.einsum('imbc,majk->ijkabc', t2, majk))
    t3c = t3c - t3c.transpose(0,1,2,4,3,5) - t3c.transpose(0,1,2,5,4,3)
    t3c = t3c - t3c.transpose(1,0,2,3,4,5) - t3c.transpose(2,1,0,3,4,5)
    t3c /= d3

    t3d = numpy.einsum('ia,bcjk->ijkabc', t1, bcjk)
    t3d += numpy.einsum('ai,jkbc->ijkabc', eris.fock[nocc:,:nocc], t2)
    t3d = t3d - t3d.transpose(0,1,2,4,3,5) - t3d.transpose(0,1,2,5,4,3)
    t3d = t3d - t3d.transpose(1,0,2,3,4,5) - t3d.transpose(2,1,0,3,4,5)
    t3d /= d3

    goovv  = numpy.einsum('kc,ijkabc->ijab', t1.conj(), t3c).conj() * (1./4)
    dovov += goovv.transpose(0,2,1,3) - goovv.transpose(0,3,1,2)

Hr2bb = Hr2bb + Hr2bb.transpose(0,1,3,2)
    Hr2bb = Hr2bb + Hr2bb.transpose(1,0,2,3)
    Hr2aa *= .5
    Hr2bb *= .5
    Hr2baaa = lib.direct_sum('Ijb+a->Ijba', Ijb, Fva)
    Hr2aaba = lib.direct_sum('ijb+A->ijAb', ijb, Fvb)
    Hr2aaba+= Fva.reshape(1,1,1,-1)
    Hr2abbb = lib.direct_sum('iJB+A->iJBA', iJB, Fvb)
    Hr2bbab = lib.direct_sum('IJB+a->IJaB', IJB, Fva)
    Hr2bbab+= Fvb.reshape(1,1,1,-1)
    Hr2baaa = Hr2baaa + Hr2baaa.transpose(0,1,3,2)
    Hr2abbb = Hr2abbb + Hr2abbb.transpose(0,1,3,2)
    Hr2baaa+= lib.direct_sum('-I+jab->Ijab', Fob, jab)
    Hr2baaa-= Foa.reshape(1,-1,1,1)
    tmpaaba = lib.direct_sum('-i+jAb->ijAb', Foa, jAb)
    Hr2abbb+= lib.direct_sum('-i+JAB->iJAB', Foa, JAB)
    Hr2abbb-= Fob.reshape(1,-1,1,1)
    tmpbbab = lib.direct_sum('-I+JaB->IJaB', Fob, JaB)
    Hr2aaba+= tmpaaba + tmpaaba.transpose(1,0,2,3)
    Hr2bbab+= tmpbbab + tmpbbab.transpose(1,0,2,3)
    tmpaaba = tmpbbab = None
    Hr2aa += Wovaa.reshape(1,nocca,1,nvira)
    Hr2aa += Wovaa.reshape(nocca,1,1,nvira)
    Hr2aa += Wovaa.reshape(nocca,1,nvira,1)
    Hr2aa += Wovaa.reshape(1,nocca,nvira,1)
    Hr2ab += Wovbb.reshape(1,noccb,1,nvirb)
    Hr2ab += Wovab.reshape(nocca,1,1,nvirb)
    Hr2ab += Wovaa.reshape(nocca,1,nvira,1)
    Hr2ab += Wovba.reshape(1,noccb,nvira,1)
    Hr2bb += Wovbb.reshape(1,noccb,1,nvirb)
    Hr2bb += Wovbb.reshape(noccb,1,1,nvirb)
    Hr2bb += Wovbb.reshape(noccb,1,nvirb,1)

l1new -= numpy.einsum('jkia,kj->ia', ooov, tmp)
    tmp = fov - einsum('kbja,jb->ka', ovov, t1)
    l1new -= numpy.einsum('ik,ka->ia', mij, tmp)
    tmp = fov - einsum('ibjc,jb->ic', ovov, t1)
    l1new -= numpy.einsum('ca,ic->ia', mba, tmp)

    tmp = einsum('kacd,jkbd->jacb', ovvv, t2)
    l1new -= einsum('ijcb,jacb->ia', l2, tmp)
    tmp = einsum('lkic,jlbc->jkib', ooov, t2)
    l1new += einsum('kjab,jkib->ia', l2, tmp)

    mo_e = eris.fock.diagonal()
    eia = lib.direct_sum('i-j->ij', mo_e[:nocc], mo_e[nocc:])
    l1new /= eia
    l1new += l1
    l2new /= lib.direct_sum('ia+jb->ijab', eia, eia)
    l2new += l2

    time0 = log.timer_debug1('update l1 l2', *time0)
    return l1new, l2new

# bbb
    d3 = lib.direct_sum('ia+jb+kc->ijkabc', eIA, eIA, eIA)
    w = numpy.einsum('ijae,kceb->ijkabc', t2bb, numpy.asarray(eris.get_OVVV()).conj())
    w-= numpy.einsum('imab,kcjm->ijkabc', t2bb, numpy.asarray(eris.OVOO.conj()))
    v = numpy.einsum('jbkc,ia->ijkabc', numpy.asarray(eris.OVOV).conj(), t1b)
    v+= numpy.einsum('jkbc,ai->ijkabc', t2bb, fVO) * .5

    rw = p6(r6(w)) / d3
    wvd = p6(w + v) / d3
    gOO = numpy.einsum('iklabc,jklabc->ij', wvd.conj(), rw) * .125
    gVV = numpy.einsum('ijkacd,ijkbcd->ab', wvd, rw.conj()) * .125
    gVO = numpy.einsum('jkbc,ijkabc->ai', t2bb.conj(), rw) * .125

    # baa
    d3 = lib.direct_sum('ia+jb+kc->ijkabc', eIA, eia, eia)
    w  = numpy.einsum('jIeA,kceb->IjkAbc', t2ab, numpy.asarray(eris.get_ovvv()).conj()) * 2
    w += numpy.einsum('jIbE,kcEA->IjkAbc', t2ab, numpy.asarray(eris.get_ovVV()).conj()) * 2
    w += numpy.einsum('jkbe,IAec->IjkAbc', t2aa, numpy.asarray(eris.get_OVvv()).conj())
    w -= numpy.einsum('mIbA,kcjm->IjkAbc', t2ab, numpy.asarray(eris.ovoo).conj()) * 2
    w -= numpy.einsum('jMbA,kcIM->IjkAbc', t2ab, numpy.asarray(eris.ovOO).conj()) * 2
    w -= numpy.einsum('jmbc,IAkm->IjkAbc', t2aa, numpy.asarray(eris.OVoo).conj())
    v  = numpy.einsum('jbkc,IA->IjkAbc', numpy.asarray(eris.ovov).conj(), t1b)
    v += numpy.einsum('kcIA,jb->IjkAbc', numpy.asarray(eris.ovOV).conj(), t1a)
    v += numpy.einsum('kcIA,jb->IjkAbc', numpy.asarray(eris.ovOV).conj(), t1a)
    v += numpy.einsum('jkbc,AI->IjkAbc', t2aa, fVO) * .5
    v += numpy.einsum('kIcA,bj->IjkAbc', t2ab, fvo) * 2

    rw = r4(w) / d3
    wvd = (w + v) / d3
    goo += numpy.einsum('kilabc,kjlabc->ij', wvd.conj(), rw) * .25
    goo += numpy.einsum('kliabc,kljabc->ij', wvd.conj(), rw) * .25

large_num (float):
            Number to replace the padded elements.
    '''
    def get_idx(x0,x1,kx,n0_p):
        return np.logical_and(n0_p[kx] >= x0, n0_p[kx] < x1)

    assert(all([len(x) == 5 for x in [pindices,qindices]]))
    p0,p1,kp,mo_e_p,nonzero_p = pindices
    q0,q1,kq,mo_e_q,nonzero_q = qindices
    fac_p, fac_q = fac

    epq = large_num * np.ones((p1-p0,q1-q0), dtype=mo_e_p[0].dtype)
    idxp = get_idx(p0,p1,kp,nonzero_p)
    idxq = get_idx(q0,q1,kq,nonzero_q)
    n0_ovp_pq = np.ix_(nonzero_p[kp][idxp], nonzero_q[kq][idxq])
    epq[n0_ovp_pq] = lib.direct_sum('p,q->pq', fac_p*mo_e_p[kp][p0:p1],
                                               fac_q*mo_e_q[kq][q0:q1])[n0_ovp_pq]
    return epq

u1a -= numpy.einsum('iajk,kj->ia', ovoo, moo1)
    u1a -= numpy.einsum('iaJK,KJ->ia', ovOO, mOO1)
    u1b -= numpy.einsum('iajk,kj->ia', OVOO, mOO1)
    u1b -= numpy.einsum('IAjk,kj->IA', OVoo, moo1)

    tmp  = fova - numpy.einsum('kbja,jb->ka', ovov, t1a)
    tmp += numpy.einsum('kaJB,JB->ka', ovOV, t1b)
    u1a -= lib.einsum('ik,ka->ia', moo, tmp)
    u1a -= lib.einsum('ca,ic->ia', mvv, tmp)
    tmp  = fovb - numpy.einsum('kbja,jb->ka', OVOV, t1b)
    tmp += numpy.einsum('jbKA,jb->KA', ovOV, t1a)
    u1b -= lib.einsum('ik,ka->ia', mOO, tmp)
    u1b -= lib.einsum('ca,ic->ia', mVV, tmp)

    eia = lib.direct_sum('i-j->ij', mo_ea_o, mo_ea_v)
    eIA = lib.direct_sum('i-j->ij', mo_eb_o, mo_eb_v)
    u1a /= eia
    u1b /= eIA

    u2aa /= lib.direct_sum('ia+jb->ijab', eia, eia)
    u2ab /= lib.direct_sum('ia+jb->ijab', eia, eIA)
    u2bb /= lib.direct_sum('ia+jb->ijab', eIA, eIA)

    time0 = log.timer_debug1('update l1 l2', *time0)
    return (u1a,u1b), (u2aa,u2ab,u2bb)

iJb =-np.einsum('ibJE,iJbE->iJb', eris_ovOV, t2ab)
    IjB =-np.einsum('jeIB,jIeB->IjB', eris_ovOV, t2ab)
    jab = np.einsum('kajb,jkab->jab',      ovov, t2aa)
    JAB = np.einsum('kajb,jkab->jab',      OVOV, t2bb)
    jAb =-np.einsum('jbKA,jKbA->jAb', eris_ovOV, t2ab)
    JaB =-np.einsum('kaJB,kJaB->JaB', eris_ovOV, t2ab)
    jaB =-np.einsum('jaKB,jKaB->jaB', eris_ovOV, t2ab)
    JAb =-np.einsum('kbJA,kJbA->JAb', eris_ovOV, t2ab)
    eris_ovov = eris_ovOV = eris_OVOV = ovov = OVOV = None
    Hr2aa = lib.direct_sum('ijb+a->ijba', ijb, Fva)
    Hr2bb = lib.direct_sum('ijb+a->ijba', IJB, Fvb)
    Hr2ab = lib.direct_sum('iJb+A->iJbA', iJb, Fvb)
    Hr2ab+= lib.direct_sum('iJB+a->iJaB', iJB, Fva)
    Hr2aa+= lib.direct_sum('-i+jab->ijab', Foa, jab)
    Hr2bb+= lib.direct_sum('-i+jab->ijab', Fob, JAB)
    Hr2ab+= lib.direct_sum('-i+JaB->iJaB', Foa, JaB)
    Hr2ab+= lib.direct_sum('-I+jaB->jIaB', Fob, jaB)
    Hr2aa = Hr2aa + Hr2aa.transpose(0,1,3,2)
    Hr2aa = Hr2aa + Hr2aa.transpose(1,0,2,3)
    Hr2bb = Hr2bb + Hr2bb.transpose(0,1,3,2)
    Hr2bb = Hr2bb + Hr2bb.transpose(1,0,2,3)
    Hr2aa *= .5
    Hr2bb *= .5
    Hr2baaa = lib.direct_sum('Ijb+a->Ijba', Ijb, Fva)
    Hr2aaba = lib.direct_sum('ijb+A->ijAb', ijb, Fvb)
    Hr2aaba+= Fva.reshape(1,1,1,-1)
    Hr2abbb = lib.direct_sum('iJB+A->iJBA', iJB, Fvb)
    Hr2bbab = lib.direct_sum('IJB+a->IJaB', IJB, Fva)
    Hr2bbab+= Fvb.reshape(1,1,1,-1)
    Hr2baaa = Hr2baaa + Hr2baaa.transpose(0,1,3,2)
    Hr2abbb = Hr2abbb + Hr2abbb.transpose(0,1,3,2)
    Hr2baaa+= lib.direct_sum('-I+jab->Ijab', Fob, jab)