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

w3 = _cp(numpy.einsum('kc,jkbc->bj', fov, t2) * 2 + fov.T)
    w3 -= numpy.einsum('kc,kjbc->bj', fov, t2)
    w3 += reduce(numpy.dot, (t1.T, fov, t1.T))
    w4 = fov.copy()

    eris_ovvv = _cp(eris.ovvv)
    eris_ovvv = lib.unpack_tril(eris_ovvv.reshape(nov,-1))
    eris_ovvv = eris_ovvv.reshape(nocc,nvir,nvir,nvir)

    wovvv = numpy.empty((nocc,nvir,nvir,nvir))
    t2tmp = numpy.empty((nocc,nvir,nocc,nvir))
    for i in range(nocc):
        wovvv[i] = eris_ovvv[i].transpose(1,0,2) * 2
        t2tmp[i] = t2[i].transpose(2,0,1)
    #:wovvv += numpy.einsum('jabd,kjdc->kabc', eris_ovvv, t2) * -1.5
    tmp = lib.dot(t2tmp.reshape(nov,-1), wovvv.reshape(-1,nvir**2),
                  -1.5/2).reshape(-1,nvir,nvir,nvir)
    g2ovvv = tmp
    for i in range(nocc):
        wovvv[i] -= eris_ovvv[i].transpose(1,2,0)
        wovvv[i] += tmp[i].transpose(1,2,0)
        g2ovvv[i] = eris_ovvv[i]*2
        g2ovvv[i] -= eris_ovvv[i].transpose(1,2,0)
    tmp = t2tmp = None

    w1 += numpy.einsum('jcba,jc->ba', eris_ovvv, t1*2)
    w1 -= numpy.einsum('jabc,jc->ba', eris_ovvv, t1)
    #:w3 += numpy.einsum('kdcb,kjdc->bj', eris_ovvv, theta)
    theta = numpy.empty(t2.shape)
    for i in range(nocc):
        theta[i] = t2[i] * 2
        theta[i] -= t2[i].transpose(0,2,1)

eris_oovv = eris_ovov = eris_oOvV = taubuf = None

    tmp, tau = tau, None
    t2ibja, eris_oVOv = eris_oVOv, None
    for i in range(nocc):
        t2ibja[i] = t2[i].transpose(2,0,1)
    #: t2new += numpy.einsum('ibkc,kcja->ijab', woVoV, t2ibja)
    lib.dot(woVoV.reshape(nov,-1), t2ibja.reshape(-1,nov), 1, tmp.reshape(nov,-1))
    for i in range(nocc):
        t2new[i] += tmp[i].transpose(1,2,0)

    #: t2new[j0:j1] += numpy.einsum('iakc,kcjb->ijab', woVoV[j0:j1], t2iajb)
    t2iajb = t2ibja
    for i in range(nocc):
        t2iajb[i] = t2[i].transpose(1,0,2)
    lib.dot(woVoV.reshape(nov,-1), t2iajb.reshape(-1,nov), 1, tmp.reshape(nov,-1))
    for i in range(nocc):
        t2new[i] += tmp[i].transpose(1,0,2)
    t2ibja = t2iajb = woVoV = tmp = None
    time1 = log.timer_debug1('contract loop', *time0)

    woooo = None
    ft_ij = foo + numpy.einsum('ja,ia->ij', .5*t1, fov)
    ft_ab = fvv - numpy.einsum('ia,ib->ab', .5*t1, fov)
    #: t2new += numpy.einsum('ijac,bc->ijab', t2, ft_ab)
    #: t2new -= numpy.einsum('ki,kjab->ijab', ft_ij, t2)
    lib.dot(t2.reshape(-1,nvir), ft_ab.T, 1, t2new.reshape(-1,nvir), 1)
    lib.dot(ft_ij.T, t2.reshape(nocc,-1),-1, t2new.reshape(nocc,-1), 1)

    mo_e = fock.diagonal()
    eia = mo_e[:nocc,None] - mo_e[None,nocc:]
    t1new += numpy.einsum('ib,ab->ia', t1, fvv)

def make_rho(ao_l, ao_h, dm_lh, dm_hl):
                ngrids = ao_l[0].shape[0]
                rho = numpy.empty((4,ngrids))
                c0 = lib.dot(ao_l[0], dm_lh)
                rho[0] = dot_bra(ao_h[0], c0)
                for i in range(1, 4):
                    rho[i] = dot_bra(ao_h[i], c0)
                c0 = lib.dot(ao_h[0], dm_hl)
                for i in range(1, 4):
                    rho[i] += dot_bra(ao_l[i], c0)
                return rho * 2  # *2 for dm_lh+dm_hl.T

def _fftn_blas(f, mesh):
    Gx = np.fft.fftfreq(mesh[0])
    Gy = np.fft.fftfreq(mesh[1])
    Gz = np.fft.fftfreq(mesh[2])
    expRGx = np.exp(np.einsum('x,k->xk', -2j*np.pi*np.arange(mesh[0]), Gx))
    expRGy = np.exp(np.einsum('x,k->xk', -2j*np.pi*np.arange(mesh[1]), Gy))
    expRGz = np.exp(np.einsum('x,k->xk', -2j*np.pi*np.arange(mesh[2]), Gz))
    out = np.empty(f.shape, dtype=np.complex128)
    buf = np.empty(mesh, dtype=np.complex128)
    for i, fi in enumerate(f):
        buf[:] = fi.reshape(mesh)
        g = lib.dot(buf.reshape(mesh[0],-1).T, expRGx, c=out[i].reshape(-1,mesh[0]))
        g = lib.dot(g.reshape(mesh[1],-1).T, expRGy, c=buf.reshape(-1,mesh[1]))
        g = lib.dot(g.reshape(mesh[2],-1).T, expRGz, c=out[i].reshape(-1,mesh[2]))
    return out.reshape(-1, *mesh)

def add_xc_(v, ao_l, ao_h, idx_l, idx_h, wv):
            for k in range(nkpts):
                aow = numpy.einsum('npi,np->pi', ao_l[k][:4], wv[:4])
                v1  = lib.dot(aow.conj().T, ao_h[k][0])
                aow = numpy.einsum('npi,np->pi', ao_h[k][1:4], wv[1:4])
                v1 += lib.dot(ao_l[k][0].conj().T, aow)
                aow = numpy.einsum('pi,p->pi', ao_h[k][1], wv[4], out=aow)
                v1 += lib.dot(ao_l[k][1].conj().T, aow)
                aow = numpy.einsum('pi,p->pi', ao_h[k][2], wv[4], out=aow)
                v1 += lib.dot(ao_l[k][2].conj().T, aow)
                aow = numpy.einsum('pi,p->pi', ao_h[k][3], wv[4], out=aow)
                v1 += lib.dot(ao_l[k][3].conj().T, aow)
                v[k,idx_l[:,None],idx_h] += v1

saved.wooov = wooov
    woooo = wooov = None

    w3 += numpy.einsum('bc,jc->bj', w1, t1)
    w3 -= numpy.einsum('kj,kb->bj', w2, t1)

    eris_ooov = _cp(eris.ooov)
    g2ooov = eris_ooov * 2
    g2ooov -= eris_ooov.transpose(2,0,1,3)
    #:tmp = numpy.einsum('kjla,jb->kabl', g2ooov, t1)
    #:wovvv = numpy.einsum('kabl,lc->kabc', tmp, t1)
    #:wovvv += numpy.einsum('kjla,jlbc->kabc', g2ooov, t2)
    tmp = lib.dot(g2ooov.reshape(nocc,-1).T, t1).reshape(-1,nocc,nvir,nvir).transpose(0,2,3,1)
    lib.dot(_cp(tmp.reshape(-1,nocc)), t1, 1, wovvv.reshape(-1,nvir), 1)
    tmp = None
    lib.dot(_cp(g2ooov.transpose(0,2,1,3).reshape(nocc**2,-1)).T,
            t2.reshape(nocc**2,-1), 1, wovvv.reshape(nov,-1), 1)
    g2ooov = eris_ooov = vOVov = vOvOv = None

    saved.wovvv = wovvv
    saved.w1 = w1
    saved.w2 = w2
    saved.w3 = w3
    saved.w4 = w4
    return saved

lib.dot(kLI[p0:p1].T, pqkR.T,  1, j3cI[k][naux:], 1)

            for k, ji in enumerate(adapted_ji_idx):
                if is_zero(kpt) and gamma_point(adapted_kptjs[k]):
                    v = fuse(j3cR[k])
                else:
                    v = fuse(j3cR[k] + j3cI[k] * 1j)
                if j2ctag == 'CD':
                    v = scipy.linalg.solve_triangular(j2c, v, lower=True, overwrite_b=True)
                    feri['j3c/%d/%d'%(ji,istep)] = v
                else:
                    feri['j3c/%d/%d'%(ji,istep)] = lib.dot(j2c, v)

                # low-dimension systems
                if j2c_negative is not None:
                    feri['j3c-/%d/%d'%(ji,istep)] = lib.dot(j2c_negative, v)

def trans(aoi, aoj, fac=1):
        if id(aoi) == id(aoj) and iden_coeffs(mo_coeffs[0], mo_coeffs[1]):
            moi = moj = numpy.asarray(lib.dot(mo_coeffs[0].T,aoi.T), order='C')
        else:
            moi = numpy.asarray(lib.dot(mo_coeffs[0].T, aoi.T), order='C')
            moj = numpy.asarray(lib.dot(mo_coeffs[1].T, aoj.T), order='C')
        mo_pairs_G = numpy.empty((nmoi,nmoj,ngrids), dtype=numpy.complex128)
        for i in range(nmoi):
            mo_pairs_G[i] = tools.fft(fac * moi[i].conj() * moj, mydf.mesh)
        mo_pairs_G = mo_pairs_G.reshape(-1,ngrids).T
        return mo_pairs_G

#:goovv += numpy.einsum('iald,jlbd->ijab', mOVov*2+mOvOv, t2) * .5
    t2a, tau2 = tau2.reshape(nocc,nvir,nocc,nvir), None
    for i in range(nocc):
        t2a[i] = t2[i].transpose(1,0,2)
    tmp = mOVov*2
    tmp += mOvOv
    lib.dot(tmp.reshape(-1,nov), t2a.reshape(nov,-1), .5, goovv.reshape(nov,-1), 1)
    t2a = tmp = None
    for i in range(nocc):
        goovv[i] = goovv[i] * 2 - goovv[i].transpose(2,1,0)
    dovov = goovv
    goooo = goovv = None

    #:gvovv += numpy.einsum('aick,kb->aibc', pvOVo, t1)
    mOVov = lib.transpose(mOVov.reshape(nov,-1))
    gvovv = lib.dot(mOVov.reshape(nocc,-1).T, t1).reshape(nvir,nocc,nvir,nvir)
    mOVov = None
    tmp = numpy.einsum('ja,jb->ab', l1, t1)
    #:gvovv += numpy.einsum('ab,ic->aibc', tmp, t1)
    #:gvovv += numpy.einsum('ba,ic->aibc', mvv, t1*.5)
    for i in range(nvir):
        gvovv[i] += numpy.einsum('b,ic->icb', tmp[i], t1)
        gvovv[i] += numpy.einsum('b,ic->icb', mvv[:,i]*.5, t1)
        gvovv[i] = gvovv[i].transpose(0,2,1)

    #:gvovv += numpy.einsum('ja,jibc->aibc', l1, t2)
    #:gvovv += numpy.einsum('jibc,ja->aibc', l2, t1)
    #:gvovv -= numpy.einsum('aibk,kc->aibc', pvOvO, t1)
    mOvOv = lib.transpose(mOvOv.reshape(nov,-1))
    lib.dot(mOvOv.reshape(nocc,-1).T, t1, -1, gvovv.reshape(-1,nvir), 1)
    mOvOv = None
    lib.dot(l1.T, t2.reshape(nocc,-1), 1, gvovv.reshape(nvir,-1), 1)

#:mOVov += numpy.einsum('ikac,jkbc->jbia', l2, theta)
    l2a = numpy.empty((nocc,nvir,nocc,nvir))
    t2a = numpy.empty((nocc,nvir,nocc,nvir))
    for i in range(nocc):
        l2a[i] = l2[i].transpose(2,0,1)
        t2a[i] = t2[i].transpose(2,0,1)
    mOvOv = lib.dot(t2a.reshape(-1,nov), l2a.reshape(-1,nov).T).reshape(nocc,nvir,nocc,nvir)
    for i in range(nocc):
        t2a[i] = t2[i].transpose(1,0,2)
    mOVov = lib.dot(t2a.reshape(-1,nov), l2a.reshape(-1,nov).T, -1).reshape(nocc,nvir,nocc,nvir)
    theta = t2a
    for i in range(nocc):
        l2a[i] = l2[i].transpose(1,0,2)
        theta[i] *= 2
        theta[i] -= t2[i].transpose(2,0,1)
    lib.dot(theta.reshape(-1,nov), l2a.reshape(nov,-1).T, 1, mOVov.reshape(nov,-1), 1)
    theta = l2a = t2a = None
    moo =(numpy.einsum('jdld->jl', mOvOv) * 2
        + numpy.einsum('jdld->jl', mOVov))
    mvv =(numpy.einsum('lbld->bd', mOvOv) * 2
        + numpy.einsum('lbld->bd', mOVov))
    mia =(numpy.einsum('kc,ikac->ia', l1, t2) * 2
        - numpy.einsum('kc,ikca->ia', l1, t2))
    mab = numpy.einsum('kc,kb->cb', l1, t1)
    mij = numpy.einsum('kc,jc->jk', l1, t1) + moo*.5

    gooov = numpy.zeros((nocc,nocc,nocc,nvir))
    tau = _ccsd.make_tau(t2, t1, t1)
    #:goooo = numpy.einsum('ijab,klab->klij', l2, tau)*.5
    goooo = lib.dot(tau.reshape(-1,nvir**2), l2.reshape(-1,nvir**2).T, .5)
    goooo = goooo.reshape(-1,nocc,nocc,nocc)
    doooo = _cp(make_theta(goooo).transpose(0,2,1,3))