Python _ccsd.unpack_tril函数代码示例

    def add_wvvVV_(self, t1, t2, eris, t2new_tril, max_memory=2000):
        time0 = time.clock(), time.time()
        nocc, nvir = t1.shape
        #: tau = t2 + numpy.einsum('ia,jb->ijab', t1, t1)
        #: t2new += numpy.einsum('ijcd,acdb->ijab', tau, vvvv)
        tau = numpy.empty((nocc*(nocc+1)//2,nvir,nvir))
        p0 = 0
        for i in range(nocc):
            tau[p0:p0+i+1] = numpy.einsum('a,jb->jab', t1[i], t1[:i+1])
            tau[p0:p0+i+1] += t2[i,:i+1]
            p0 += i + 1
        time0 = logger.timer_debug1(self, 'vvvv-tau', *time0)

        p0 = 0
        outbuf = numpy.empty((nvir,nvir,nvir))
        for a in range(nvir):
            buf = _ccsd.unpack_tril(eris.vvvv[p0:p0+a+1], out=outbuf[:a+1])
            #: t2new_tril[i,:i+1, a] += numpy.einsum('xcd,cdb->xb', tau[:,:a+1], buf)
            lib.numpy_helper._dgemm('N', 'N', nocc*(nocc+1)//2, nvir, (a+1)*nvir,
                                    tau.reshape(-1,nvir*nvir), buf.reshape(-1,nvir),
                                    t2new_tril.reshape(-1,nvir*nvir), 1, 1,
                                    0, 0, a*nvir)

            #: t2new_tril[i,:i+1,:a] += numpy.einsum('xd,abd->xab', tau[:,a], buf[:a])
            if a > 0:
                lib.numpy_helper._dgemm('N', 'T', nocc*(nocc+1)//2, a*nvir, nvir,
                                        tau.reshape(-1,nvir*nvir), buf.reshape(-1,nvir),
                                        t2new_tril.reshape(-1,nvir*nvir), 1, 1,
                                        a*nvir, 0, 0)
            p0 += a+1
            time0 = logger.timer_debug1(self, 'vvvv %d'%a, *time0)
        return t2new_tril

def make_intermediates(mycc, t1, t2, eris):
    saved = ccsd_lambda.make_intermediates(mycc, t1, t2, eris)

    nocc, nvir = t1.shape
    eris_ovvv = _cp(eris.ovvv)
    eris_ovvv = _ccsd.unpack_tril(eris_ovvv.reshape(nocc*nvir,-1))
    eris_ovvv = eris_ovvv.reshape(nocc,nvir,nvir,nvir)
    mo_e = mycc._scf.mo_energy
    eia = lib.direct_sum('i-a->ia',mo_e[:nocc], mo_e[nocc:])
    d3 = lib.direct_sum('ia,jb,kc->ijkabc', eia, eia, eia)
    eris_ovoo = eris.ovoo
    w =(numpy.einsum('iabf,kjcf->ijkabc', eris_ovvv, t2)
      - numpy.einsum('iajm,mkbc->ijkabc', eris_ovoo, t2)) / d3
    v = numpy.einsum('iajb,kc->ijkabc', eris.ovov, t1) / d3 * .5
    w = ccsd_t.p6_(w)
    v = ccsd_t.p6_(v)
    rwv = ccsd_t.r6_(w*2+v)
    jov = numpy.einsum('jbkc,ijkabc->ia', eris.ovov, ccsd_t.r6_(w))
    joovv = numpy.einsum('iabf,ijkabc->kjcf', eris_ovvv, rwv)
    joovv-= numpy.einsum('iajm,ijkabc->mkbc', eris.ovoo, rwv)
    joovv = joovv + joovv.transpose(1,0,3,2)

    saved.jov = jov
    saved.joovv = joovv
    return saved

def kernel(cc, eris, t1=None, t2=None, max_memory=2000, verbose=logger.INFO):
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(cc.stdout, verbose)

    if t1 is None: t1 = cc.t1
    if t2 is None: t2 = cc.t2

    nocc, nvir = t1.shape
    nmo = nocc + nvir
    mo_e = cc._scf.mo_energy
    eia = mo_e[:nocc,None] - mo_e[nocc:]
    eiajb = eia.reshape(-1,1) + eia.reshape(-1)
    d3 = eia.reshape(-1,1) + eiajb.reshape(-1)
    d3 = d3.reshape(nocc,nvir,nocc,nvir,nocc,nvir).transpose(0,2,4,1,3,5)

    eris_ovvv = _ccsd.unpack_tril(eris.ovvv.reshape(nocc*nvir,-1))
    eris_ovvv = eris_ovvv.reshape(nocc,nvir,nvir,nvir)
    eris_ovoo = eris.ovoo
    eris_ovov = eris.ovov
    w =(numpy.einsum('iabf,kjcf->ijkabc', eris_ovvv, t2)
      - numpy.einsum('iajm,mkbc->ijkabc', eris_ovoo, t2))
    #r = r6_(w)
    #v = numpy.einsum('iajb,kc->ijkabc', eris_ovov, t1)
    #wvd = p6_(w + .5 * v) / d3
    #et = numpy.einsum('ijkabc,ijkabc', wvd, r) * 2

#    wvd = (w + .5 * v) / d3
#    et =(numpy.einsum('ijkabc,ijkabc', wvd, r)
#       + numpy.einsum('ikjacb,ijkabc', wvd, r)
#       + numpy.einsum('jikbac,ijkabc', wvd, r)
#       + numpy.einsum('jkibca,ijkabc', wvd, r)
#       + numpy.einsum('kijcab,ijkabc', wvd, r)
#       + numpy.einsum('kjicba,ijkabc', wvd, r)) * 2

    r = r6_(p6_(w))
    v = numpy.einsum('iajb,kc->ijkabc', eris_ovov, t1)
    wvd = p6_(w + .5 * v) / d3
    wt = numpy.zeros((nvir,nvir,nvir))
    for i in range(nvir):
        for j in range(i+1):
            for k in range(j+1):
                wt[i,j,k] = 1
                if i == j or j == k:
                    wt[i,j,k] = 1./2
                    if i == k:
                        wt[i,j,k] *= 1./3
    et = numpy.einsum('ijkabc,ijkabc,abc', wvd, r, wt) * 2
    return et

def IX_intermediates(mycc, t1, t2, l1, l2, d1=None, d2=None, eris=None):
    if d1 is None:
        doo, dov, dvo, dvv = ccsd_t_rdm.gamma1_intermediates(mycc, t1, t2, l1, l2, eris)
    else:
        doo, dov, dvo, dvv = d1
    if d2 is None:
# Note gamma2 are in Chemist's notation
        d2 = ccsd_t_rdm.gamma2_intermediates(mycc, t1, t2, l1, l2, eris)
    else:
        dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov = d2
    if dovvv is None:
        dovvv = dvvov.transpose(2,3,0,1)
    elif dvvov is None:
        dvvov = dovvv.transpose(2,3,0,1)

    nocc, nvir = t1.shape
    eris_ovvv = _ccsd.unpack_tril(_cp(eris.ovvv).reshape(nocc*nvir,-1))
    eris_ovvv = eris_ovvv.reshape(nocc,nvir,nvir,nvir)
    eris_vvvv = pyscf.ao2mo.restore(1, _cp(eris.vvvv), nvir)
    dvvvv = pyscf.ao2mo.restore(1, _cp(dvvvv), nvir)

# Note Ioo is not hermitian
    Ioo  =(numpy.einsum('jakb,iakb->ij', dovov, eris.ovov)
         + numpy.einsum('kbja,iakb->ij', dovov, eris.ovov))
    Ioo +=(numpy.einsum('jabk,iakb->ij', dovvo, eris.ovov)
         + numpy.einsum('kbaj,iakb->ij', dovvo, eris.ovov)
         + numpy.einsum('jkab,ikab->ij', doovv, eris.oovv)
         + numpy.einsum('kjba,ikab->ij', doovv, eris.oovv))
    Ioo +=(numpy.einsum('jmlk,imlk->ij', doooo, eris.oooo) * 2
         + numpy.einsum('mjkl,imlk->ij', doooo, eris.oooo) * 2)
    Ioo +=(numpy.einsum('jlka,ilka->ij', dooov, eris.ooov)
         + numpy.einsum('klja,klia->ij', dooov, eris.ooov))
    Ioo += numpy.einsum('abjc,icab->ij', dvvov, eris_ovvv)
    Ioo += numpy.einsum('ljka,lika->ij', dooov, eris.ooov)
    Ioo *= -1

# Note Ivv is not hermitian
    Ivv  =(numpy.einsum('ibjc,iajc->ab', dovov, eris.ovov)
         + numpy.einsum('jcib,iajc->ab', dovov, eris.ovov))
    Ivv +=(numpy.einsum('jcbi,iajc->ab', dovvo, eris.ovov)
         + numpy.einsum('ibcj,iajc->ab', dovvo, eris.ovov)
         + numpy.einsum('jibc,jiac->ab', doovv, eris.oovv)
         + numpy.einsum('ijcb,jiac->ab', doovv, eris.oovv))
    Ivv +=(numpy.einsum('bced,aced->ab', dvvvv, eris_vvvv) * 2
         + numpy.einsum('cbde,aced->ab', dvvvv, eris_vvvv) * 2)
    Ivv +=(numpy.einsum('dbic,icda->ab', dvvov, eris_ovvv)
         + numpy.einsum('dcib,iadc->ab', dvvov, eris_ovvv))
    Ivv += numpy.einsum('bcid,idac->ab', dvvov, eris_ovvv)
    Ivv += numpy.einsum('jikb,jika->ab', dooov, eris.ooov)
    Ivv *= -1

    Ivo  =(numpy.einsum('kajb,kijb->ai', dovov, eris.ooov)
         + numpy.einsum('kbja,jikb->ai', dovov, eris.ooov))
    Ivo +=(numpy.einsum('acbd,icbd->ai', dvvvv, eris_ovvv) * 2
         + numpy.einsum('cadb,icbd->ai', dvvvv, eris_ovvv) * 2)
    Ivo +=(numpy.einsum('jbak,jbik->ai', dovvo, eris.ovoo)
         + numpy.einsum('kabj,jbik->ai', dovvo, eris.ovoo)
         + numpy.einsum('jkab,jkib->ai', doovv, eris.ooov)
         + numpy.einsum('kjba,jkib->ai', doovv, eris.ooov))
    Ivo +=(numpy.einsum('dajc,idjc->ai', dvvov, eris.ovov)
         + numpy.einsum('dcja,jidc->ai', dvvov, eris.oovv))
    Ivo += numpy.einsum('abjc,ibjc->ai', dvvov, eris.ovov)
    Ivo += numpy.einsum('jlka,jlki->ai', dooov, eris.oooo)
    Ivo *= -1

    Xvo  =(numpy.einsum('kj,kjia->ai', doo, eris.ooov) * 2
         + numpy.einsum('kj,kjia->ai', doo, eris.ooov) * 2
         - numpy.einsum('kj,kija->ai', doo, eris.ooov)
         - numpy.einsum('kj,ijka->ai', doo, eris.ooov))
    Xvo +=(numpy.einsum('cb,iacb->ai', dvv, eris_ovvv) * 2
         + numpy.einsum('cb,iacb->ai', dvv, eris_ovvv) * 2
         - numpy.einsum('cb,icab->ai', dvv, eris_ovvv)
         - numpy.einsum('cb,ibca->ai', dvv, eris_ovvv))
    Xvo +=(numpy.einsum('icjb,jbac->ai', dovov, eris_ovvv)
         + numpy.einsum('jcib,jcab->ai', dovov, eris_ovvv))
    Xvo +=(numpy.einsum('iklj,ljka->ai', doooo, eris.ooov) * 2
         + numpy.einsum('kijl,ljka->ai', doooo, eris.ooov) * 2)
    Xvo +=(numpy.einsum('ibcj,jcab->ai', dovvo, eris_ovvv)
         + numpy.einsum('jcbi,jcab->ai', dovvo, eris_ovvv)
         + numpy.einsum('ijcb,jacb->ai', doovv, eris_ovvv)
         + numpy.einsum('jibc,jacb->ai', doovv, eris_ovvv))
    Xvo +=(numpy.einsum('ijkb,jakb->ai', dooov, eris.ovov)
         + numpy.einsum('kjib,kjab->ai', dooov, eris.oovv))
    Xvo += numpy.einsum('dbic,dbac->ai', dvvov, eris_vvvv)
    Xvo += numpy.einsum('jikb,jakb->ai', dooov, eris.ovov)
    Xvo += Ivo
    return Ioo, Ivv, Ivo, Xvo

def update_amps(cc, t1, t2, eris):
    time0 = time.clock(), time.time()
    log = logger.Logger(cc.stdout, cc.verbose)
    nocc, nvir = t1.shape
    nov = nocc*nvir
    fock = eris.fock
    t1new = numpy.zeros_like(t1)
    t2new = numpy.zeros_like(t2)

#** make_inter_F
    fov = fock[:nocc,nocc:].copy()

    foo = fock[:nocc,:nocc].copy()
    foo[range(nocc),range(nocc)] = 0
    foo += .5 * numpy.einsum('ia,ja->ij', fock[:nocc,nocc:], t1)

    fvv = fock[nocc:,nocc:].copy()
    fvv[range(nvir),range(nvir)] = 0
    fvv -= .5 * numpy.einsum('ia,ib->ab', t1, fock[:nocc,nocc:])

    #: woooo = numpy.einsum('la,ikja->ikjl', t1, eris.ooov)
    eris_ooov = _cp(eris.ooov)
    foo += numpy.einsum('kc,jikc->ij', 2*t1, eris_ooov)
    foo += numpy.einsum('kc,jkic->ij',  -t1, eris_ooov)
    woooo = lib.dot(eris_ooov.reshape(-1,nvir), t1.T).reshape((nocc,)*4)
    woooo = lib.transpose_sum(woooo.reshape(nocc*nocc,-1), inplace=True)
    woooo += _cp(eris.oooo).reshape(nocc**2,-1)
    woooo = _cp(woooo.reshape(nocc,nocc,nocc,nocc).transpose(0,2,1,3))
    time1 = log.timer_debug1('woooo', *time0)

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

    fvv += numpy.einsum('kc,kcba->ab', 2*t1, eris_ovvv)
    fvv += numpy.einsum('kc,kbca->ab',  -t1, eris_ovvv)

    #: woVoV = numpy.einsum('ka,ijkb->ijba', t1, eris.ooov)
    #: woVoV -= numpy.einsum('jc,icab->ijab', t1, eris_ovvv)
    woVoV = lib.dot(_cp(eris_ooov.transpose(0,1,3,2).reshape(-1,nocc)), t1)
    woVoV = woVoV.reshape(nocc,nocc,nvir,nvir)

#: tau = t2 + numpy.einsum('ia,jb->ijab', t1, t1)
#: tmp = numpy.einsum('ijcd,kcdb->kijb', tau, eris.ovvv)
#: t2new += numpy.einsum('ka,kijb->jiba', -t1, tmp)
    tau = make_tau(t2, t1, t1)
    tmp = numpy.empty((nocc,nocc,nocc,nvir))
    for k in range(nocc):
        tmp[k] = lib.dot(tau.reshape(-1,nvir**2),
                         eris_ovvv[k].reshape(-1,nvir)).reshape(nocc,nocc,nvir).transpose(1,0,2)
        lib.dot(t1, eris_ovvv[k].reshape(nvir,-1), -1, woVoV[k].reshape(nocc,-1), 1)
    lib.dot(tmp.reshape(nocc,-1).T, t1, -1, t2new.reshape(-1,nvir), 1)
    tmp = None

#: wOVov += numpy.einsum('iabc,jc->ijab', eris.ovvv, t1)
#: wOVov -= numpy.einsum('jbik,ka->jiba', eris.ovoo, t1)
#: t2new += woVoV.transpose()
    #: wOVov = -numpy.einsum('jbik,ka->ijba', eris.ovoo, t1)
    wOVov, tau = tau, None
    lib.dot(_cp(_cp(eris.ooov).transpose(0,2,3,1).reshape(-1,nocc)), t1,
            -1, wOVov.reshape(-1,nvir))
    #: wOVov += numpy.einsum('iabc,jc->jiab', eris_ovvv, t1)
    lib.dot(t1, eris_ovvv.reshape(-1,nvir).T, 1, wOVov.reshape(nocc,-1), 1)
    for i in range(nocc):
        t2new[i] += wOVov[i].transpose(0,2,1)

#: theta = t2.transpose(0,1,3,2) * 2 - t2
#: t1new += numpy.einsum('ijcb,jcba->ia', theta, eris.ovvv)
    theta = numpy.empty((nocc,nocc,nvir,nvir))
    for i in range(nocc):
        theta[i] = t2[i].transpose(0,2,1) * 2
        theta[i] -= t2[i]
        lib.dot(_cp(theta[i].transpose(0,2,1).reshape(nocc,-1)),
                eris_ovvv[i].reshape(-1,nvir), 1, t1new, 1)
    eris_ovvv = None

    eris_ovov = _cp(eris.ovov)

    for i in range(nocc):
        t2new[i] += eris_ovov[i].transpose(1,0,2) * .5

    fov += numpy.einsum('kc,iakc->ia', t1, eris_ovov) * 2
    fov -= numpy.einsum('kc,icka->ia', t1, eris_ovov)

#: theta = t2.transpose(1,0,2,3) * 2 - t2
#: t1new += numpy.einsum('jb,ijab->ia', fov, theta)
#: t1new -= numpy.einsum('ikjb,kjab->ia', eris.ooov, theta)
    t1new += numpy.einsum('jb,jiab->ia', fov, theta)
    #: t1new -= numpy.einsum('kijb,kjab->ia', eris.ooov, theta)
    lib.dot(_cp(eris_ooov.transpose(1,0,2,3).reshape(nocc,-1)),
            theta.reshape(-1,nvir), -1, t1new, 1)
    eris_ooov = None

#: wOVov += eris.ovov.transpose(0,1,3,2)
#: theta = t2.transpose(1,0,2,3) * 2 - t2
#: tau = theta - numpy.einsum('ic,kb->ikcb', t1, t1*2)
#: wOVov += .5 * numpy.einsum('jakc,ikcb->jiba', eris.ovov, tau)
#: wOVov -= .5 * numpy.einsum('jcka,ikcb->jiba', eris.ovov, t2)
#: t2new += numpy.einsum('ikca,kjbc->ijba', theta, wOVov)
    wOVov = _cp(wOVov.transpose(0,3,1,2))
#.........这里部分代码省略.........

def make_intermediates(mycc, t1, t2, eris, max_memory=2000):
    log = logger.Logger(mycc.stdout, mycc.verbose)
    nocc, nvir = t1.shape
    nov = nocc * nvir
    foo = eris.fock[:nocc,:nocc]
    fov = eris.fock[:nocc,nocc:]
    fvv = eris.fock[nocc:,nocc:]

    class _Saved(object):
        def __init__(self):
            self._tmpfile = tempfile.NamedTemporaryFile()
            self.ftmp = h5py.File(self._tmpfile.name)
        def __del__(self):
            if hasattr(self, 'ftmp'):
                self.ftmp.close()
                self._tmpfile = None
    saved = _Saved()
    saved.woooo = saved.ftmp.create_dataset('woooo', (nocc,nocc,nocc,nocc), 'f8')
    saved.wooov = saved.ftmp.create_dataset('wooov', (nocc,nocc,nocc,nvir), 'f8')
    saved.wOVov = saved.ftmp.create_dataset('wOVov', (nocc,nvir,nocc,nvir), 'f8')
    saved.wOvOv = saved.ftmp.create_dataset('wOvOv', (nocc,nvir,nocc,nvir), 'f8')
    saved.wovvv = saved.ftmp.create_dataset('wovvv', (nocc,nvir,nvir,nvir), 'f8')

# As we don't have l2 in memory, hold tau temporarily in memory
    w1 = fvv - numpy.einsum('ja,jb->ba', fov, t1)
    w2 = foo + numpy.einsum('ib,jb->ij', fov, t1)
    w3 = 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()

    _tmpfile = tempfile.NamedTemporaryFile()
    fswap = h5py.File(_tmpfile.name)

    time1 = time.clock(), time.time()
    max_memory = max_memory - lib.current_memory()[0]
    unit = max(nocc*nvir**2*4 + nvir**3*2,
               nvir**3*3 + nocc*nvir**2,
               nocc*nvir**2*6 + nocc**2*nvir + nocc**3 + nocc**2*nvir)
    blksize = max(ccsd.BLKMIN, int(max_memory*.95e6/8/unit))
    log.debug1('ccsd lambda make_intermediates: block size = %d, nocc = %d in %d blocks',
               blksize, nocc, int((nocc+blksize-1)//blksize))
    for istep, (p0, p1) in enumerate(prange(0, nocc, blksize)):
        eris_ovvv = _cp(eris.ovvv[p0:p1])
        eris_ovvv = _ccsd.unpack_tril(eris_ovvv.reshape((p1-p0)*nvir,-1))
        eris_ovvv = eris_ovvv.reshape(p1-p0,nvir,nvir,nvir)
        w1 += numpy.einsum('jcba,jc->ba', eris_ovvv, t1[p0:p1]*2)
        w1 -= numpy.einsum('jabc,jc->ba', eris_ovvv, t1[p0:p1])
        #:w3 += numpy.einsum('kdcb,kjdc->bj', eris_ovvv, theta)
        for i in range(p1-p0):
            theta = t2[p0+i] * 2
            theta -= t2[p0+i].transpose(0,2,1)
            w3 += lib.dot(eris_ovvv[i].reshape(-1,nvir).T,
                          _cp(theta.reshape(nocc,-1)).T)
        theta = None
        #:wOVov = numpy.einsum('jbcd,kd->jbkc', eris_ovvv, t1)
        #:wOvOv = numpy.einsum('jdcb,kd->jbkc', eris_ovvv, -t1)
        wOVov = lib.dot(eris_ovvv.reshape(-1,nvir),
                        t1.T).reshape(-1,nvir,nvir,nocc).transpose(0,1,3,2).copy()
        g2ovvv = _cp(eris_ovvv.transpose(0,2,3,1))
        wOvOv = lib.dot(g2ovvv.reshape(-1,nvir),
                        -t1.T).reshape(-1,nvir,nvir,nocc).transpose(0,1,3,2).copy()
        for i in range(p1-p0):
            g2ovvv[i] *= 2
            g2ovvv[i] -= eris_ovvv[i].transpose(1,0,2)
        wooov = numpy.empty((p1-p0,nocc,nocc,nvir))
        woooo = numpy.empty((p1-p0,nocc,nocc,nocc))
        eris_ovov = _cp(_cp(eris.ovov[p0:p1]).transpose(0,2,1,3))
        for j0, j1 in prange(0, nocc, blksize):
            tau = _ccsd.make_tau(t2[j0:j1], t1[j0:j1], t1)
            #:wooov[:,j0:j1] = numpy.einsum('icbd,jkbd->ijkc', g2ovvv, tau)
            #:woooo[:,:,j0:j1] = numpy.einsum('icjd,klcd->ijkl', eris_ovov, tau)
            tmp = lib.dot(g2ovvv.reshape(-1,nvir**2), tau.reshape(-1,nvir**2).T)
            wooov[:,j0:j1] = tmp.reshape(-1,nvir,j1-j0,nocc).transpose(0,2,3,1)
            woooo[:,:,j0:j1] = lib.dot(eris_ovov.reshape(-1,nvir**2),
                                       tau.reshape(-1,nvir**2).T).reshape(-1,nocc,j1-j0,nocc)
        eris_ovov = eris_ovvv = g2ovvv = tau = tmp = None
#==== mem usage nocc*nvir**2*2 + nocc**2*nvir + nocc**3 + nvir**3*2 + nocc*nvir**2*2

        eris_ooov = _cp(eris.ooov[p0:p1])
        w2[p0:p1] += numpy.einsum('ijkb,kb->ij', eris_ooov, t1) * 2
        w2 -= numpy.einsum('kjib,kb->ij', eris_ooov, t1[p0:p1])
        #:w3 -= numpy.einsum('kjlc,klbc->bj', eris_ooov, theta)
        for i in range(p1-p0):
            theta = t2[p0+i].transpose(0,2,1) * 2
            theta -= t2[p0+i]
            w3 -= lib.dot(theta.reshape(-1,nvir).T, eris_ooov[i].reshape(nocc,-1).T)
        theta = None
        #:woooo += numpy.einsum('ikjc,lc->ijkl', eris_ooov, t1)
        #:wOvOv += numpy.einsum('jklb,lc->jbkc', eris_ooov, t1)
        woooo += lib.dot(eris_ooov.reshape(-1,nvir),
                         t1.T).reshape((-1,nocc,nocc,nocc)).transpose(0,2,1,3)
        for i in range(p1-p0):
            lib.dot(_cp(eris_ooov[i].transpose(2,0,1)).reshape(-1,nocc),
                    t1, 1, wOvOv[i].reshape(-1,nvir), 1)
            wooov[i] += eris_ooov[i].transpose(1,0,2)*2
            wooov[i] -= eris_ooov[i]

        eris_ovoo = _cp(eris.ovoo[p0:p1])
        #:woooo += numpy.einsum('icjl,kc->ijkl', eris_ovoo, t1)
#.........这里部分代码省略.........

def update_amps(mycc, t1, t2, l1, l2, eris=None, saved=None, max_memory=2000):
    if saved is None:
        saved = make_intermediates(mycc, t1, t2, eris)
    time1 = time0 = time.clock(), time.time()
    log = logger.Logger(mycc.stdout, mycc.verbose)
    nocc, nvir = t1.shape
    nov = nocc * nvir
    foo = eris.fock[:nocc,:nocc]
    fov = eris.fock[:nocc,nocc:]
    fvv = eris.fock[:nocc,:nocc]

    #:mba = numpy.einsum('klca,klcb->ba', l2, t2*2-t2.transpose(0,1,3,2))
    #:mij = numpy.einsum('ikcd,jkcd->ij', l2, t2*2-t2.transpose(0,1,3,2))
    #:theta = t2*2 - t2.transpose(0,1,3,2)
    theta = _ccsd.make_0132(t2, t2, 2, -1)
    mba = lib.dot(theta.reshape(-1,nvir).T, l2.reshape(-1,nvir))
    mij = lib.dot(l2.reshape(nocc,-1), theta.reshape(nocc,-1).T)
    theta = None
    mba1 = numpy.einsum('jc,jb->bc', l1, t1) + mba
    mij1 = numpy.einsum('kb,jb->kj', l1, t1) + mij
    mia1 =(t1 + numpy.einsum('kc,jkbc->jb', l1, t2) * 2
         - numpy.einsum('kc,jkcb->jb', l1, t2)
         - reduce(numpy.dot, (t1, l1.T, t1))
         - numpy.einsum('bd,jd->jb', mba, t1)
         - numpy.einsum('lj,lb->jb', mij, t1))

    tmp = mycc.add_wvvVV(numpy.zeros_like(l1), l2, eris, max_memory)
    l2new = numpy.empty((nocc,nocc,nvir,nvir))
    ij = 0
    for i in range(nocc):
        for j in range(i):
            tmp1 = tmp[ij] * .5  # *.5 because of l2+l2.transpose(1,0,3,2) later
            l2new[i,j] = tmp1
            l2new[j,i] = tmp1.T
            ij += 1
        l2new[i,i] = tmp[ij] * .5
        ij += 1
    l1new =(numpy.einsum('ijab,jb->ia', l2new, t1) * 4
          - numpy.einsum('jiab,jb->ia', l2new, t1) * 2)
    tmp = tmp1 = None

    l1new += fov
    l1new += numpy.einsum('ib,ba->ia', l1, saved.w1)
    l1new -= numpy.einsum('ja,ij->ia', l1, saved.w2)
    l1new -= numpy.einsum('ik,ka->ia', mij, saved.w4)
    l1new -= numpy.einsum('ca,ic->ia', mba, saved.w4)
    l1new += numpy.einsum('ijab,bj->ia', l2, saved.w3) * 2
    l1new -= numpy.einsum('ijba,bj->ia', l2, saved.w3)

    l2new += numpy.einsum('ia,jb->ijab', l1, saved.w4)
    #:l2new += numpy.einsum('jibc,ca->jiba', l2, saved.w1)
    #:l2new -= numpy.einsum('kiba,jk->jiba', l2, saved.w2)
    lib.dot(l2.reshape(-1,nvir), saved.w1, 1, l2new.reshape(-1,nvir), 1)
    lib.dot(saved.w2, l2.reshape(nocc,-1),-1, l2new.reshape(nocc,-1), 1)

    eris_ooov = _cp(eris.ooov)
    l1new -= numpy.einsum('jkia,kj->ia', eris_ooov, mij1) * 2
    l1new += numpy.einsum('ikja,kj->ia', eris_ooov, mij1)
    #:l2new -= numpy.einsum('ka,kijb->jiba', l1, eris_ooov)
    lib.dot(_cp(eris_ooov.transpose(0,2,1,3).reshape(nocc,-1)).T,
            l1, -1, l2new.reshape(-1,nvir), 1)
    eris_ooov = None

    tau = _ccsd.make_tau(t2, t1, t1)
    #:l2tau = numpy.einsum('ijcd,klcd->ijkl', l2, tau)
    l2tau = lib.dot(l2.reshape(nocc**2,-1),
                    tau.reshape(nocc**2,-1).T).reshape((nocc,)*4)
    tau = None
    #:l2t1 = numpy.einsum('jidc,kc->ijkd', l2, t1)
    l2t1 = lib.dot(l2.reshape(-1,nvir), t1.T).reshape(nocc,nocc,nvir,nocc)
    l2t1 = _cp(l2t1.transpose(1,0,3,2))

    max_memory = max_memory - lib.current_memory()[0]
    unit = max(nvir**3*2+nocc*nvir**2, nocc*nvir**2*5)
    blksize = min(nocc, max(ccsd.BLKMIN, int(max_memory*.95e6/8/unit)))
    log.debug1('block size = %d, nocc = %d is divided into %d blocks',
               blksize, nocc, int((nocc+blksize-1)/blksize))
    for p0, p1 in prange(0, nocc, blksize):
        eris_ovvv = _cp(eris.ovvv[p0:p1])
        eris_ovvv = _ccsd.unpack_tril(eris_ovvv.reshape((p1-p0)*nvir,-1))
        eris_ovvv = eris_ovvv.reshape(p1-p0,nvir,nvir,nvir)

        l1new[p0:p1] += numpy.einsum('iabc,bc->ia', eris_ovvv, mba1) * 2
        l1new[p0:p1] -= numpy.einsum('ibca,bc->ia', eris_ovvv, mba1)
        #:l2new[p0:p1] += numpy.einsum('ic,jbac->jiba', l1, eris_ovvv)
        tmp = lib.dot(l1, eris_ovvv.reshape(-1,nvir).T)
        l2new[p0:p1] += tmp.reshape(nocc,-1,nvir,nvir).transpose(1,0,2,3)
        tmp = None
        m4buf = numpy.empty((blksize,nocc,nvir,nvir))
        eris_ovvv = _cp(eris_ovvv.transpose(0,2,1,3).reshape(-1,nvir**2))
        for j0, j1 in prange(0, nocc, blksize):
            #:m4 = numpy.einsum('ijkd,kadb->ijab', l2t1[j0:j1,:,p0:p1], eris_ovvv)
            m4 = m4buf[:j1-j0]
            lib.dot(_cp(l2t1[j0:j1,:,p0:p1].reshape((j1-j0)*nocc,-1)),
                    eris_ovvv, 1, m4.reshape(-1,nvir**2))
            l2new[j0:j1] -= m4
            l1new[j0:j1] -= numpy.einsum('ijab,jb->ia', m4, t1) * 2
            l1new -= numpy.einsum('ijab,ia->jb', m4, t1[j0:j1]) * 2
            l1new += numpy.einsum('jiab,jb->ia', m4, t1[j0:j1])
            l1new[j0:j1] += numpy.einsum('jiab,ia->jb', m4, t1)
#.........这里部分代码省略.........

def _rdm2_mo2ao(mycc, d2, dm1, mo_coeff):
    log = logger.Logger(mycc.stdout, mycc.verbose)
    time1 = time.clock(), time.time()
    dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov = d2
    nocc, nvir = dovov.shape[:2]
    nov = nocc * nvir
    nao, nmo = mo_coeff.shape
    nao_pair = nao * (nao+1) // 2
    nvir_pair = nvir * (nvir+1) //2
    mo_coeff = numpy.asarray(mo_coeff, order='F')
    def _trans(vin, i0, icount, j0, jcount, out=None):
        nrow = vin.shape[0]
        fdrv = getattr(_ccsd.libcc, 'AO2MOnr_e2_drv')
        pao_loc = ctypes.POINTER(ctypes.c_void_p)()
        ftrans = ctypes.c_void_p(_ctypes.dlsym(_ccsd.libcc._handle, 'AO2MOtranse2_nr_s1'))
        fmmm = ctypes.c_void_p(_ctypes.dlsym(_ccsd.libcc._handle, 'CCmmm_transpose_sum'))
        fdrv(ftrans, fmmm,
             out.ctypes.data_as(ctypes.c_void_p),
             vin.ctypes.data_as(ctypes.c_void_p),
             mo_coeff.ctypes.data_as(ctypes.c_void_p),
             ctypes.c_int(nrow), ctypes.c_int(nao),
             ctypes.c_int(i0), ctypes.c_int(icount),
             ctypes.c_int(j0), ctypes.c_int(jcount),
             pao_loc, ctypes.c_int(0))
        return out

    blksize = 4
    bufop = numpy.empty((nocc*nmo,nao_pair))
    bufvv = numpy.empty((nvir_pair,nao_pair))
    pool = numpy.empty((blksize,nmo,nmo,nmo))
    for p0, p1 in prange(0, nocc, blksize):
        buf1 = pool[:p1-p0]
        buf1[:,:nocc,:nocc,:nocc] = doooo[p0:p1]
        buf1[:,:nocc,:nocc,nocc:] = dooov[p0:p1]
        buf1[:,:nocc,nocc:,:nocc] = 0
        buf1[:,:nocc,nocc:,nocc:] = doovv[p0:p1]
        buf1[:,nocc:,:nocc,:nocc] = 0
        buf1[:,nocc:,:nocc,nocc:] = dovov[p0:p1]
        buf1[:,nocc:,nocc:,:nocc] = dovvo[p0:p1]
        buf1[:,nocc:,nocc:,nocc:] = dovvv[p0:p1]
        for i in range(p1-p0):
            buf1[i,p0+i,:,:] += dm1
            buf1[i,:,:,p0+i] -= dm1 * .5
        _trans(buf1.reshape(-1,nmo**2), 0, nmo, 0, nmo, bufop[p0*nmo:p1*nmo])

    pool = pool.ravel()[:blksize*nvir**3].reshape((blksize*nvir,nvir,nvir))
    for p0, p1 in prange(0, nvir_pair, blksize*nvir):
        buf1 = _ccsd.unpack_tril(_cp(dvvvv[p0:p1]), out=pool[:p1-p0])
        _trans(buf1, nocc, nvir, nocc, nvir, bufvv[p0:p1])
    pool = buf1 = None

    # trans ij of d_ijkl
    bufaa = numpy.empty((nao_pair,nao_pair))
    tmpbuf1 = numpy.empty((BLKSIZE,nmo,nmo))
    tmpbuf2 = numpy.empty((BLKSIZE,nvir_pair))
    tmpbuf3 = numpy.empty((BLKSIZE,nvir,nvir))
    tmpbuf1[:,nocc:,:nocc] = 0
    for i0, i1 in prange(0, nao_pair, BLKSIZE):
        buf1 = tmpbuf1[:i1-i0].reshape(i1-i0,-1)
        for j0, j1 in prange(0, nocc*nmo, BLKSIZE):
            buf1[:,j0:j1] = bufop[j0:j1,i0:i1].T
        buf1 = tmpbuf1[:i1-i0]
        tmp  = tmpbuf2[:i1-i0]
        for j0, j1 in prange(0, nvir_pair, BLKSIZE):
            tmp[:,j0:j1] = bufvv[j0:j1,i0:i1].T
        buf1[:,nocc:,nocc:] = _ccsd.unpack_tril(tmp, out=tmpbuf3[:i1-i0])
        _trans(buf1, 0, nmo, 0, nmo, bufaa[i0:i1])
# dm2 + dm2.transpose(2,3,0,1)
        for j0, j1, in prange(0, i0, BLKSIZE):
            bufaa[i0:i1,j0:j1] += bufaa[j0:j1,i0:i1].T
            bufaa[j0:j1,i0:i1] = bufaa[i0:i1,j0:j1].T
        bufaa[i0:i1,i0:i1] = bufaa[i0:i1,i0:i1] + bufaa[i0:i1,i0:i1].T
    tmpbuf1 = tmpbuf2 = tmpbuf3 = bufop = bufvv = pool = None

    idx = numpy.arange(nao)
    idx = idx*(idx+1)//2 + idx
    bufaa[:,idx] *= .5

    return bufaa

def _rdm2_mo2ao(mycc, d2, dm1, mo_coeff, fsave=None, max_memory=2000):
    log = logger.Logger(mycc.stdout, mycc.verbose)
    if fsave is None:
        _dm2file = tempfile.NamedTemporaryFile()
        fsave = h5py.File(_dm2file.name, 'w')
    else:
        _dm2file = None
    time1 = time.clock(), time.time()
    dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov = d2
    nocc, nvir = dovov.shape[:2]
    nov = nocc * nvir
    nao, nmo = mo_coeff.shape
    nao_pair = nao * (nao+1) // 2
    nvir_pair = nvir * (nvir+1) //2
    mo_coeff = numpy.asarray(mo_coeff, order='F')
    def _trans(vin, i0, icount, j0, jcount, out=None):
        nrow = vin.shape[0]
        if out is None:
            out = numpy.empty((nrow,nao_pair))
        fdrv = getattr(_ccsd.libcc, 'AO2MOnr_e2_drv')
        pao_loc = ctypes.POINTER(ctypes.c_void_p)()
        ftrans = ctypes.c_void_p(_ctypes.dlsym(_ccsd.libcc._handle, 'AO2MOtranse2_nr_s1'))
        fmmm = ctypes.c_void_p(_ctypes.dlsym(_ccsd.libcc._handle, 'CCmmm_transpose_sum'))
        fdrv(ftrans, fmmm,
             out.ctypes.data_as(ctypes.c_void_p),
             vin.ctypes.data_as(ctypes.c_void_p),
             mo_coeff.ctypes.data_as(ctypes.c_void_p),
             ctypes.c_int(nrow), ctypes.c_int(nao),
             ctypes.c_int(i0), ctypes.c_int(icount),
             ctypes.c_int(j0), ctypes.c_int(jcount),
             pao_loc, ctypes.c_int(0))
        return out

# transform dm2_ij to get lower triangular (dm2+dm2.transpose(0,1,3,2))
    _tmpfile = tempfile.NamedTemporaryFile()
    fswap = h5py.File(_tmpfile.name)
    max_memory1 = max_memory - lib.current_memory()[0]
    blksize = max(1, int(max_memory1*1e6/8/(nmo*nao_pair+nmo**3+nvir**3)))
    iobuflen = int(256e6/8/(blksize*nmo))
    log.debug1('_rdm2_mo2ao pass 1: blksize = %d, iobuflen = %d', blksize, iobuflen)
    fswap.create_group('o')  # for h5py old version
    pool1 = numpy.empty((blksize,nmo,nmo,nmo))
    pool2 = numpy.empty((blksize,nmo,nao_pair))
    bufd_ovvv = numpy.empty((blksize,nvir,nvir,nvir))
    for istep, (p0, p1) in enumerate(prange(0, nocc, blksize)):
        buf1 = pool1[:p1-p0]
        buf1[:,:nocc,:nocc,:nocc] = doooo[p0:p1]
        buf1[:,:nocc,:nocc,nocc:] = dooov[p0:p1]
        buf1[:,:nocc,nocc:,:nocc] = 0
        buf1[:,:nocc,nocc:,nocc:] = doovv[p0:p1]
        buf1[:,nocc:,:nocc,:nocc] = 0
        buf1[:,nocc:,:nocc,nocc:] = dovov[p0:p1]
        buf1[:,nocc:,nocc:,:nocc] = dovvo[p0:p1]
        d_ovvv = bufd_ovvv[:p1-p0]
        ao2mo.outcore._load_from_h5g(dovvv, p0*nvir, p1*nvir,
                                     d_ovvv.reshape(-1,nvir**2))
        buf1[:,nocc:,nocc:,nocc:] = d_ovvv
        for i in range(p0, p1):
            buf1[i-p0,i,:,:] += dm1
            buf1[i-p0,:,:,i] -= dm1 * .5
        buf2 = pool2[:p1-p0].reshape(-1,nao_pair)
        _trans(buf1.reshape(-1,nmo**2), 0, nmo, 0, nmo, buf2)
        ao2mo.outcore._transpose_to_h5g(fswap, 'o/%d'%istep, buf2, iobuflen)
    pool1 = pool2 = bufd_ovvv = None
    time1 = log.timer_debug1('_rdm2_mo2ao pass 1', *time1)

    fswap.create_group('v')  # for h5py old version
    pool1 = numpy.empty((blksize*nvir,nao_pair))
    pool2 = numpy.empty((blksize*nvir,nvir,nvir))
    for istep, (p0, p1) in enumerate(prange(0, nvir_pair, blksize*nvir)):
        buf1 = _cp(dvvvv[p0:p1])
        buf2 = _ccsd.unpack_tril(buf1, out=pool2[:p1-p0])
        buf1 = _trans(buf2, nocc, nvir, nocc, nvir, out=pool1[:p1-p0])
        ao2mo.outcore._transpose_to_h5g(fswap, 'v/%d'%istep, buf1, iobuflen)
    pool1 = pool2 = None
    time1 = log.timer_debug1('_rdm2_mo2ao pass 2', *time1)

# transform dm2_kl then dm2 + dm2.transpose(2,3,0,1)
    max_memory1 = max_memory - lib.current_memory()[0]
    blksize = max(nao, int(max_memory1*1e6/8/(nao_pair+nmo**2)))
    iobuflen = int(256e6/8/blksize)
    log.debug1('_rdm2_mo2ao pass 3: blksize = %d, iobuflen = %d', blksize, iobuflen)
    gsave = fsave.create_group('dm2')
    for istep, (p0, p1) in enumerate(prange(0, nao_pair, blksize)):
        gsave.create_dataset(str(istep), (nao_pair,p1-p0), 'f8')
    diagidx = numpy.arange(nao)
    diagidx = diagidx*(diagidx+1)//2 + diagidx
    pool1 = numpy.empty((blksize,nmo,nmo))
    pool2 = numpy.empty((blksize,nvir_pair))
    pool3 = numpy.empty((blksize,nvir,nvir))
    pool4 = numpy.empty((blksize,nao_pair))
    for istep, (p0, p1) in enumerate(prange(0, nao_pair, blksize)):
        buf1 = pool1[:p1-p0]
        ao2mo.outcore._load_from_h5g(fswap['o'], p0, p1,
                                     buf1[:,:nocc].reshape(p1-p0,-1))
        buf2 = ao2mo.outcore._load_from_h5g(fswap['v'], p0, p1, pool2[:p1-p0])
        buf3 = _ccsd.unpack_tril(buf2, out=pool3[:p1-p0])
        buf1[:,nocc:,nocc:] = buf3
        buf1[:,nocc:,:nocc] = 0
        buf2 = _trans(buf1, 0, nmo, 0, nmo, out=pool4[:p1-p0])
#.........这里部分代码省略.........

#.........这里部分代码省略.........
        eris_ooov = None
        d_oovv = _ccsd.precontract(d_oovv.reshape(-1,nvir,nvir)).reshape(p1-p0,nocc,-1)

        d_ovvv = numpy.empty((p1-p0,nvir,nvir,nvir))
        ao2mo.outcore._load_from_h5g(dovvv, p0*nvir, p1*nvir,
                                     d_ovvv.reshape(-1,nvir**2))
        #:Ivo += numpy.einsum('jadc,jidc->ai', d_ovvv, eris_oovv)
        for i in range(p1-p0):
            Ivo += lib.dot(d_ovvv[i].reshape(nvir,-1), eris_oovv[i].reshape(nocc,-1).T)
        eris_oovv = None

        # tril part of (d_ovvv + d_ovvv.transpose(0,1,3,2))
        c_ovvv = _ccsd.precontract(d_ovvv.reshape(-1,nvir,nvir))
        ao2mo.outcore._transpose_to_h5g(fswap, 'c_vvov/%d'%istep, c_ovvv, iobuflen)
        c_ovvv = c_ovvv.reshape(-1,nvir,nvir_pair)
        eris_ovx = _cp(eris.ovvv[p0:p1])
        ao2mo.outcore._transpose_to_h5g(fswap, 'e_vvov/%d'%istep,
                                        eris_ovx.reshape(-1,nvir_pair), iobuflen)
        #:Xvo += numpy.einsum('jibc,jabc->ai', d_oovv, eris_ovvv)
        #:Ivv += numpy.einsum('ibdc,iadc->ab', d_ovvv, eris_ovvv)
        for i in range(p1-p0):
            lib.dot(eris_ovx[i].reshape(nvir,-1),
                    d_oovv[i].reshape(nocc,-1).T, 1, Xvo, 1)
            lib.dot(eris_ovx[i].reshape(nvir,-1),
                    c_ovvv[i].reshape(nvir,-1).T, 1, Ivv, 1)
        c_ovvv = d_oovv = None

        eris_ovvo = numpy.empty((p1-p0,nvir,nvir,nocc))
        for i in range(p1-p0):
            d_ovvv[i] = _ccsd.sum021(d_ovvv[i])
            eris_ovvo[i] = eris_ovov[i].transpose(0,2,1)
        #:Ivo += numpy.einsum('abjc,ibjc->ai', d_ovvv, eris_ovov)
        Ivo += lib.dot(d_ovvv.reshape(-1,nvir).T, eris_ovvo.reshape(-1,nocc))
        eris_ovvo = eris_ovov = None

        eris_ovvv = _ccsd.unpack_tril(eris_ovx.reshape(-1,nvir_pair))
        eris_ovx = None
        eris_ovvv = eris_ovvv.reshape(p1-p0,nvir,nvir,nvir)
        #:Ivv += numpy.einsum('icdb,icda->ab', d_ovvv, eris_ovvv)
        #:Xvo += numpy.einsum('jibc,jabc->ai', d_oovv, eris_ovvv)
        Ivv += lib.dot(eris_ovvv.reshape(-1,nvir).T, d_ovvv.reshape(-1,nvir))
        Xvo[:,p0:p1] +=(numpy.einsum('cb,iacb->ai', dvv, eris_ovvv) * 4
                      - numpy.einsum('cb,icba->ai', dvv+dvv.T, eris_ovvv))

        d_ovvo = _cp(fswap['dovvo'][p0:p1])
        #:Xvo += numpy.einsum('jbic,jbca->ai', d_ovov, eris_ovvv)
        lib.dot(eris_ovvv.reshape(-1,nvir).T, d_ovvo.reshape(-1,nocc), 1, Xvo, 1)

        d_ovvv = d_ovvo = eris_ovvv = None

    max_memory1 = max_memory - lib.current_memory()[0]
    unit = nocc*nvir**2 + nvir**3*2.5
    blksize = max(ccsd.BLKMIN, int(max_memory1*1e6/8/unit))
    log.debug1('IX_intermediates pass 2: block size = %d, nocc = %d in %d blocks',
               blksize, nocc, int((nocc+blksize-1)/blksize))
    for p0, p1 in prange(0, nvir, blksize):
        off0 = p0*(p0+1)//2
        off1 = p1*(p1+1)//2
        d_vvvv = _cp(dvvvv[off0:off1]) * 4
        for i in range(p0, p1):
            d_vvvv[i*(i+1)//2+i-off0] *= .5
        d_vvvv = _ccsd.unpack_tril(d_vvvv)
        eris_vvvv = _ccsd.unpack_tril(_cp(eris.vvvv[off0:off1]))
        #:Ivv += numpy.einsum('decb,deca->ab', d_vvvv, eris_vvvv) * 2
        #:Xvo += numpy.einsum('dbic,dbca->ai', d_vvov, eris_vvvv)
        lib.dot(eris_vvvv.reshape(-1,nvir).T, d_vvvv.reshape(-1,nvir), 2, Ivv, 1)
        #:d_vvvv = _cp(d_vvvv + d_vvvv.transpose(0,1,3,2))
        d_vvov = numpy.empty((off1-off0,nocc,nvir))
        ao2mo.outcore._load_from_h5g(fswap['c_vvov'], off0, off1, d_vvov.reshape(-1,nov))
        d_vvvo = _cp(d_vvov.transpose(0,2,1))
        lib.dot(eris_vvvv.reshape(-1,nvir).T, d_vvvo.reshape(-1,nocc), 1, Xvo, 1)
        d_vvov = eris_vvvv = None

        eris_vvov = numpy.empty((off1-off0,nocc,nvir))
        ao2mo.outcore._load_from_h5g(fswap['e_vvov'], off0, off1,
                                     eris_vvov.reshape(-1,nov))
        eris_vvvo = _cp(eris_vvov.transpose(0,2,1))
        #:Ioo += numpy.einsum('abjc,abci->ij', d_vvov, eris_vvvo)
        #:Ivo += numpy.einsum('dbca,dbci->ai', d_vvvv, eris_vvvo) * 2
        lib.dot(d_vvvv.reshape(-1,nvir).T, eris_vvvo.reshape(-1,nocc), 2, Ivo, 1)
        lib.dot(eris_vvvo.reshape(-1,nocc).T, d_vvvo.reshape(-1,nocc), 1, Ioo, 1)
        eris_vvov = eris_vovv = d_vvvv = None

    del(fswap['e_vvov'])
    del(fswap['c_vvov'])
    del(fswap['dovvo'])
    fswap.close()
    _tmpfile = None

    if d2 is None:
        for key in fd2intermediate.keys():
            del(fd2intermediate[key])
        fd2intermediate.close()
        _d2tmpfile = None

    Ioo *= -1
    Ivv *= -1
    Ivo *= -1
    Xvo += Ivo
    return Ioo, Ivv, Ivo, Xvo

def make_intermediates(mycc, t1, t2, eris):
    log = logger.Logger(mycc.stdout, mycc.verbose)
    nocc, nvir = t1.shape
    nov = nocc * nvir
    foo = eris.fock[:nocc,:nocc]
    fov = eris.fock[:nocc,nocc:]
    fvv = eris.fock[nocc:,nocc:]

    class _Saved(object):
        pass
    saved = _Saved()

# As we don't have l2 in memory, hold tau temporarily in memory
    w1 = fvv - numpy.einsum('ja,jb->ba', fov, t1)
    w2 = foo + numpy.einsum('ib,jb->ij', fov, t1)
    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 = _ccsd.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)
        lib.dot(eris_ovvv[i].reshape(-1,nvir).T,
                _cp(theta[i].reshape(nocc,-1)).T, 1, w3, 1)

    theta = _cp(theta.transpose(0,2,1,3))
    #:vkbca = numpy.einsum('jdca,kbjd->kbca', g2ovvv, theta)
    vkbca = lib.dot(_cp(theta.reshape(nov,-1)),
                    g2ovvv.reshape(-1,nvir*nvir)).reshape(-1,nvir,nvir,nvir)
    for i in range(nocc):
        wovvv[i] += vkbca[i].transpose(2,0,1)
        wovvv[i] -= vkbca[i].transpose(2,1,0) * .5
    vkabc = None

    #:wOVov = numpy.einsum('jbcd,kd->jbkc', eris_ovvv, t1)
    #:wOvOv = numpy.einsum('jdcb,kd->jbkc', eris_ovvv, -t1)
    wOVov = lib.dot(eris_ovvv.reshape(-1,nvir),
                    t1.T).reshape(-1,nvir,nvir,nocc).transpose(0,1,3,2).copy()
    for i in range(nocc):
        g2ovvv[i] = eris_ovvv[i].transpose(1,2,0) * 2
    wOvOv = lib.dot(g2ovvv.reshape(-1,nvir),
                    -t1.T, .5).reshape(-1,nvir,nvir,nocc).transpose(0,1,3,2).copy()
    for i in range(nocc):
        g2ovvv[i] -= eris_ovvv[i].transpose(1,0,2)
    eris_ovov = _cp(_cp(eris.ovov).transpose(0,2,1,3))
    tau = _ccsd.make_tau(t2, t1, t1)
    #:wooov[:,j0:j1] = numpy.einsum('icbd,jkbd->ijkc', g2ovvv, tau)
    #:woooo[:,:,j0:j1] = numpy.einsum('icjd,klcd->ijkl', eris_ovov, tau)
    tmp = lib.dot(g2ovvv.reshape(-1,nvir**2), tau.reshape(-1,nvir**2).T)
    wooov = _cp(tmp.reshape(-1,nvir,nocc,nocc).transpose(0,2,3,1))
    woooo = lib.dot(eris_ovov.reshape(-1,nvir**2),
                    tau.reshape(-1,nvir**2).T).reshape(-1,nocc,nocc,nocc)
    eris_ovov = eris_ovvv = g2ovvv = tau = tmp = None

    eris_ooov = _cp(eris.ooov)
    eris_ovoo = _cp(eris.ovoo)
    #:woooo += numpy.einsum('icjl,kc->ijkl', eris_ovoo, t1)
    #:wOVov += numpy.einsum('jblk,lc->jbkc', eris_ovoo, -t1)
    for i in range(nocc):
        woooo[i] += lib.dot(t1, eris_ovoo[i].reshape(nvir,-1)).reshape((nocc,)*3).transpose(1,0,2)
    lib.dot(eris_ovoo.reshape(-1,nocc), t1, -1, wOVov.reshape(-1,nvir), 1)
    #:wooov -= numpy.einsum('ibjl,lkcb->ijkc', eris_ovoo*1.5, t2)
    t2tmp = numpy.empty((nocc,nvir,nocc,nvir))
    for i in range(nocc):
        t2tmp[i] = t2[i].transpose(2,0,1)
    tmp_ooov = _cp(-eris_ooov.transpose(2,0,1,3)).reshape(-1,nov)
    lib.dot(tmp_ooov, t2tmp.reshape(nov,-1), 1.5, wooov.reshape(-1,nov), 1)
    t2tmp = None

    g2ooov, tmp_ooov = tmp_ooov.reshape(nocc,nocc,nocc,nvir), None
    g2ooov += eris_ooov * 2
    #:vikjc = numpy.einsum('iklb,jlcb->ikjc', g2ooov, theta)
    vikjc = lib.dot(g2ooov.reshape(-1,nov), theta.reshape(-1,nov).T)
    vikjc = vikjc.reshape(nocc,nocc,nocc,nvir)
    wooov += vikjc.transpose(0,2,1,3)
#.........这里部分代码省略.........

def update_amps(mycc, t1, t2, l1, l2, eris=None, saved=None):
    if saved is None:
        saved = make_intermediates(mycc, t1, t2, eris)
    time1 = time0 = time.clock(), time.time()
    log = logger.Logger(mycc.stdout, mycc.verbose)
    nocc, nvir = t1.shape
    nov = nocc * nvir
    foo = eris.fock[:nocc,:nocc]
    fov = eris.fock[:nocc,nocc:]
    fvv = eris.fock[:nocc,:nocc]

    #:mba = numpy.einsum('klca,klcb->ba', l2, t2*2-t2.transpose(0,1,3,2))
    #:mij = numpy.einsum('ikcd,jkcd->ij', l2, t2*2-t2.transpose(0,1,3,2))
    #:theta = t2*2 - t2.transpose(0,1,3,2)
    theta = _ccsd.make_0132(t2, t2, 2, -1)
    mba = lib.dot(theta.reshape(-1,nvir).T, l2.reshape(-1,nvir))
    mij = lib.dot(l2.reshape(nocc,-1), theta.reshape(nocc,-1).T)
    theta = None
    mba1 = numpy.einsum('jc,jb->bc', l1, t1) + mba
    mij1 = numpy.einsum('kb,jb->kj', l1, t1) + mij
    mia1 =(t1 + numpy.einsum('kc,jkbc->jb', l1, t2) * 2
         - numpy.einsum('kc,jkcb->jb', l1, t2)
         - reduce(numpy.dot, (t1, l1.T, t1))
         - numpy.einsum('bd,jd->jb', mba, t1)
         - numpy.einsum('lj,lb->jb', mij, t1))

    tmp = mycc.add_wvvVV(numpy.zeros_like(l1), l2, eris)
    l2new = numpy.empty((nocc,nocc,nvir,nvir))
    ij = 0
    for i in range(nocc):
        for j in range(i):
            tmp1 = tmp[ij] * .5  # *.5 because of l2+l2.transpose(1,0,3,2) later
            l2new[i,j] = tmp1
            l2new[j,i] = tmp1.T
            ij += 1
        l2new[i,i] = tmp[ij] * .5
        ij += 1
    l1new =(numpy.einsum('ijab,jb->ia', l2new, t1) * 4
          - numpy.einsum('jiab,jb->ia', l2new, t1) * 2)
    tmp = tmp1 = None

    l1new += eris.fock[:nocc,nocc:]
    l1new += numpy.einsum('ib,ba->ia', l1, saved.w1)
    l1new -= numpy.einsum('ja,ij->ia', l1, saved.w2)
    l1new -= numpy.einsum('ik,ka->ia', mij, saved.w4)
    l1new -= numpy.einsum('ca,ic->ia', mba, saved.w4)
    l1new += numpy.einsum('ijab,bj->ia', l2, saved.w3) * 2
    l1new -= numpy.einsum('ijba,bj->ia', l2, saved.w3)

    l2new += numpy.einsum('ia,jb->ijab', l1, saved.w4)
    #:l2new += numpy.einsum('jibc,ca->jiba', l2, saved.w1)
    #:l2new -= numpy.einsum('kiba,jk->jiba', l2, saved.w2)
    lib.dot(l2.reshape(-1,nvir), saved.w1, 1, l2new.reshape(-1,nvir), 1)
    lib.dot(saved.w2, l2.reshape(nocc,-1),-1, l2new.reshape(nocc,-1), 1)

    eris_ooov = _cp(eris.ooov)
    l1new -= numpy.einsum('jkia,kj->ia', eris_ooov, mij1) * 2
    l1new += numpy.einsum('ikja,kj->ia', eris_ooov, mij1)
    #:l2new -= numpy.einsum('ka,kijb->jiba', l1, eris_ooov)
    lib.dot(_cp(eris_ooov.transpose(0,2,1,3).reshape(nocc,-1)).T,
            l1, -1, l2new.reshape(-1,nvir), 1)
    eris_ooov = None

    tau = _ccsd.make_tau(t2, t1, t1)
    #:l2tau = numpy.einsum('ijcd,klcd->ijkl', l2, tau)
    l2tau = lib.dot(l2.reshape(nocc**2,-1),
                    tau.reshape(nocc**2,-1).T).reshape((nocc,)*4)
    tau = None
    l2t1 = numpy.einsum('ijcd,kc->ijkd', l2, t1)

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

    l1new += numpy.einsum('iabc,bc->ia', eris_ovvv, mba1) * 2
    l1new -= numpy.einsum('ibca,bc->ia', eris_ovvv, mba1)
    #:l2new += numpy.einsum('ic,jbac->jiba', l1, eris_ovvv)
    tmp = lib.dot(l1, eris_ovvv.reshape(-1,nvir).T).reshape(nocc,-1,nvir,nvir)
    for i in range(nocc):
        l2new[i] += tmp[i].transpose(0,2,1)
    #:m4 = numpy.einsum('ijkd,kadb->ijab', l2t1, eris_ovvv)
    m4 = tmp
    lib.dot(_cp(l2t1.reshape(nocc*nocc,-1)),
            _cp(eris_ovvv.transpose(0,2,1,3).reshape(-1,nvir**2)),
            1, m4.reshape(nocc*nocc,-1))
    l2new -= m4
    l1new -= numpy.einsum('ijab,jb->ia', m4, t1) * 2
    l1new -= numpy.einsum('ijab,ia->jb', m4, t1) * 2
    l1new += numpy.einsum('jiab,jb->ia', m4, t1)
    l1new += numpy.einsum('jiab,ia->jb', m4, t1)
    eris_ovvv = tmp = None

    eris_ovov = _cp(eris.ovov)
    l1new += numpy.einsum('jb,iajb->ia', l1, eris_ovov) * 2
    #:l2new -= numpy.einsum('jbic,ca->jiba', eris_ovov, mba1)
    #:l2new -= numpy.einsum('kajb,ik->ijab', eris_ovov, mij1)
    tmp = lib.dot(eris_ovov.reshape(-1,nvir), mba1).reshape(nocc,nvir,nocc,nvir)
    lib.dot(mij1, eris_ovov.reshape(nocc,-1), 1, tmp.reshape(nocc,-1), 1)
    tmp_oovv = numpy.empty((nocc,nocc,nvir,nvir))
    for i in range(nocc):
#.........这里部分代码省略.........

def update_amps(cc, t1, t2, eris, max_memory=2000):
    time0 = time.clock(), time.time()
    log = logger.Logger(cc.stdout, cc.verbose)
    nocc, nvir = t1.shape
    nov = nocc*nvir
    fock = eris.fock
    t1new = numpy.zeros_like(t1)
    t2new = numpy.zeros_like(t2)

#** make_inter_F
    fov = fock[:nocc,nocc:].copy()

    foo = fock[:nocc,:nocc].copy()
    foo[range(nocc),range(nocc)] = 0
    foo += .5 * numpy.einsum('ia,ja->ij', fock[:nocc,nocc:], t1)

    fvv = fock[nocc:,nocc:].copy()
    fvv[range(nvir),range(nvir)] = 0
    fvv -= .5 * numpy.einsum('ia,ib->ab', t1, fock[:nocc,nocc:])

    #: woooo = numpy.einsum('la,ikja->ikjl', t1, eris.ooov)
    eris_ooov = _cp(eris.ooov)
    foo += numpy.einsum('kc,jikc->ij', 2*t1, eris_ooov)
    foo += numpy.einsum('kc,jkic->ij',  -t1, eris_ooov)
    woooo = lib.dot(eris_ooov.reshape(-1,nvir), t1.T).reshape((nocc,)*4)
    woooo = lib.transpose_sum(woooo.reshape(nocc*nocc,-1), inplace=True)
    woooo += _cp(eris.oooo).reshape(nocc**2,-1)
    woooo = _cp(woooo.reshape(nocc,nocc,nocc,nocc).transpose(0,2,1,3))
    eris_ooov = None
    time1 = log.timer_debug1('woooo', *time0)

    unit = _memory_usage_inloop(nocc, nvir)*1e6/8
    max_memory = max_memory - lib.current_memory()[0]
    blksize = max(BLKMIN, int(max_memory*.95e6/8/unit))
    log.debug1('block size = %d, nocc = %d is divided into %d blocks',
               blksize, nocc, int((nocc+blksize-1)//blksize))

    for p0, p1 in prange(0, nocc, blksize):
# ==== read eris.ovvv ====
        eris_ovvv = _cp(eris.ovvv[p0:p1])
        eris_ovvv = _ccsd.unpack_tril(eris_ovvv.reshape((p1-p0)*nvir,-1))
        eris_ovvv = eris_ovvv.reshape(p1-p0,nvir,nvir,nvir)

        fvv += numpy.einsum('kc,kcba->ab', 2*t1[p0:p1], eris_ovvv)
        fvv += numpy.einsum('kc,kbca->ab',  -t1[p0:p1], eris_ovvv)

    #: tau = t2 + numpy.einsum('ia,jb->ijab', t1, t1)
    #: tmp = numpy.einsum('ijcd,kdcb->ijbk', tau, eris.ovvv)
    #: t2new += numpy.einsum('ka,ijbk->ijba', -t1, tmp)
        #: eris_vvov = eris_ovvv.transpose(1,2,0,3).copy()
        eris_vvov = _cp(eris_ovvv.transpose(2,1,0,3).reshape(nvir*nvir,-1))
        tmp = numpy.empty((nocc,nocc,p1-p0,nvir))
        taubuf = numpy.empty((blksize,nocc,nvir,nvir))
        for j0, j1 in prange(0, nocc, blksize):
            tau = make_tau(t2[j0:j1], t1[j0:j1], t1, 1, out=taubuf[:j1-j0])
            lib.dot(tau.reshape(-1,nvir*nvir), eris_vvov, 1,
                    tmp[j0:j1].reshape((j1-j0)*nocc,-1), 0)
        tmp = _cp(tmp.transpose(0,1,3,2).reshape(-1,p1-p0))
        lib.dot(tmp, t1[p0:p1], -1, t2new.reshape(-1,nvir), 1)
        tau = tmp = eris_vvov = None
        #==== mem usage blksize*(nvir**3*2+nvir*nocc**2*2)

    #: wOVov += numpy.einsum('iabc,jc->ijab', eris.ovvv, t1)
    #: wOVov -= numpy.einsum('jbik,ka->jiba', eris.ovoo, t1)
    #: t2new += woVoV.transpose()
        #: wOVov = -numpy.einsum('jbik,ka->ijba', eris.ovoo[p0:p1], t1)
        tmp = _cp(eris.ovoo[p0:p1].transpose(2,0,1,3))
        wOVov = lib.dot(tmp.reshape(-1,nocc), t1, -1)
        tmp = None
        wOVov = wOVov.reshape(nocc,p1-p0,nvir,nvir)
        #: wOVov += numpy.einsum('iabc,jc->jiab', eris_ovvv, t1)
        lib.dot(t1, eris_ovvv.reshape(-1,nvir).T, 1, wOVov.reshape(nocc,-1), 1)
        t2new[p0:p1] += wOVov.transpose(1,0,2,3)

        eris_ooov = _cp(eris.ooov[p0:p1])
        #: woVoV = numpy.einsum('ka,ijkb->ijba', t1, eris.ooov[p0:p1])
        #: woVoV -= numpy.einsum('jc,icab->ijab', t1, eris_ovvv)
        woVoV = lib.dot(_cp(eris_ooov.transpose(0,1,3,2).reshape(-1,nocc)), t1)
        woVoV = woVoV.reshape(p1-p0,nocc,nvir,nvir)
        for i in range(eris_ovvv.shape[0]):
            lib.dot(t1, eris_ovvv[i].reshape(nvir,-1), -1,
                    woVoV[i].reshape(nocc,-1), 1)

    #: theta = t2.transpose(0,1,3,2) * 2 - t2
    #: t1new += numpy.einsum('ijcb,jcba->ia', theta, eris.ovvv)
        theta = numpy.empty((p1-p0,nocc,nvir,nvir))
        for i in range(p1-p0):
            theta[i] = t2[p0+i].transpose(0,2,1) * 2
            theta[i] -= t2[p0+i]
            lib.dot(_cp(theta[i].transpose(0,2,1).reshape(nocc,-1)),
                    eris_ovvv[i].reshape(-1,nvir), 1, t1new, 1)
        eris_ovvv = None
        time2 = log.timer_debug1('ovvv [%d:%d]'%(p0, p1), *time1)
        #==== mem usage blksize*(nvir**3+nocc*nvir**2*4)

# ==== read eris.ovov ====
        eris_ovov = _cp(eris.ovov[p0:p1])
        #==== mem usage blksize*(nocc*nvir**2*4)

        for i in range(p1-p0):
#.........这里部分代码省略.........

    def __init__(self, cc, mo_coeff=None, method='incore'):
        cput0 = (time.clock(), time.time())
        moidx = numpy.ones(cc.mo_energy.size, dtype=numpy.bool)
        if isinstance(cc.frozen, (int, numpy.integer)):
            moidx[:cc.frozen] = False
        elif len(cc.frozen) > 0:
            moidx[numpy.asarray(cc.frozen)] = False
        if mo_coeff is None:
            self.mo_coeff = mo_coeff = cc.mo_coeff[:,moidx]
            self.fock = numpy.diag(cc.mo_energy[moidx])
        else:  # If mo_coeff is not canonical orbital
            self.mo_coeff = mo_coeff = mo_coeff[:,moidx]
            dm = cc._scf.make_rdm1(cc.mo_coeff, cc.mo_occ)
            fockao = cc._scf.get_hcore() + cc._scf.get_veff(cc.mol, dm)
            self.fock = reduce(numpy.dot, (mo_coeff.T, fockao, mo_coeff))

        nocc = cc.nocc()
        nmo = cc.nmo()
        nvir = nmo - nocc
        mem_incore, mem_outcore, mem_basic = _mem_usage(nocc, nvir)
        mem_now = pyscf.lib.current_memory()[0]

        log = logger.Logger(cc.stdout, cc.verbose)
        if (method == 'incore' and cc._scf._eri is not None and
            (mem_incore+mem_now < cc.max_memory) or cc.mol.incore_anyway):
            eri1 = pyscf.ao2mo.incore.full(cc._scf._eri, mo_coeff)
            #:eri1 = pyscf.ao2mo.restore(1, eri1, nmo)
            #:self.oooo = eri1[:nocc,:nocc,:nocc,:nocc].copy()
            #:self.ooov = eri1[:nocc,:nocc,:nocc,nocc:].copy()
            #:self.ovoo = eri1[:nocc,nocc:,:nocc,:nocc].copy()
            #:self.oovv = eri1[:nocc,:nocc,nocc:,nocc:].copy()
            #:self.ovov = eri1[:nocc,nocc:,:nocc,nocc:].copy()
            #:ovvv = eri1[:nocc,nocc:,nocc:,nocc:].copy()
            #:self.ovvv = numpy.empty((nocc,nvir,nvir*(nvir+1)//2))
            #:for i in range(nocc):
            #:    for j in range(nvir):
            #:        self.ovvv[i,j] = lib.pack_tril(ovvv[i,j])
            #:self.vvvv = pyscf.ao2mo.restore(4, eri1[nocc:,nocc:,nocc:,nocc:], nvir)
            nvir_pair = nvir * (nvir+1) // 2
            self.oooo = numpy.empty((nocc,nocc,nocc,nocc))
            self.ooov = numpy.empty((nocc,nocc,nocc,nvir))
            self.ovoo = numpy.empty((nocc,nvir,nocc,nocc))
            self.oovv = numpy.empty((nocc,nocc,nvir,nvir))
            self.ovov = numpy.empty((nocc,nvir,nocc,nvir))
            self.ovvv = numpy.empty((nocc,nvir,nvir_pair))
            self.vvvv = numpy.empty((nvir_pair,nvir_pair))
            ij = 0
            outbuf = numpy.empty((nmo,nmo,nmo))
            for i in range(nocc):
                buf = _ccsd.unpack_tril(eri1[ij:ij+i+1], out=outbuf[:i+1])
                for j in range(i+1):
                    self.oooo[i,j] = self.oooo[j,i] = buf[j,:nocc,:nocc]
                    self.ooov[i,j] = self.ooov[j,i] = buf[j,:nocc,nocc:]
                    self.oovv[i,j] = self.oovv[j,i] = buf[j,nocc:,nocc:]
                ij += i + 1
            ij1 = 0
            for i in range(nocc,nmo):
                buf = _ccsd.unpack_tril(eri1[ij:ij+i+1], out=outbuf[:i+1])
                self.ovoo[:,i-nocc] = buf[:nocc,:nocc,:nocc]
                self.ovov[:,i-nocc] = buf[:nocc,:nocc,nocc:]
                for j in range(nocc):
                    self.ovvv[j,i-nocc] = lib.pack_tril(_cp(buf[j,nocc:,nocc:]))
                for j in range(nocc, i+1):
                    self.vvvv[ij1] = lib.pack_tril(_cp(buf[j,nocc:,nocc:]))
                    ij1 += 1
                ij += i + 1
        else:
            cput1 = time.clock(), time.time()
            _tmpfile1 = tempfile.NamedTemporaryFile()
            _tmpfile2 = tempfile.NamedTemporaryFile()
            self.feri1 = h5py.File(_tmpfile1.name)
            orbo = mo_coeff[:,:nocc]
            orbv = mo_coeff[:,nocc:]
            nvpair = nvir * (nvir+1) // 2
            self.oooo = self.feri1.create_dataset('oooo', (nocc,nocc,nocc,nocc), 'f8')
            self.ooov = self.feri1.create_dataset('ooov', (nocc,nocc,nocc,nvir), 'f8')
            self.ovoo = self.feri1.create_dataset('ovoo', (nocc,nvir,nocc,nocc), 'f8')
            self.oovv = self.feri1.create_dataset('oovv', (nocc,nocc,nvir,nvir), 'f8')
            self.ovov = self.feri1.create_dataset('ovov', (nocc,nvir,nocc,nvir), 'f8')
            self.ovvv = self.feri1.create_dataset('ovvv', (nocc,nvir,nvpair), 'f8')

            self.feri2 = h5py.File(_tmpfile2.name, 'w')
            pyscf.ao2mo.full(cc.mol, orbv, self.feri2, verbose=log)
            self.vvvv = self.feri2['eri_mo']
            cput1 = log.timer_debug1('transforming vvvv', *cput1)

            tmpfile3 = tempfile.NamedTemporaryFile()
            with h5py.File(tmpfile3.name, 'w') as feri:
                pyscf.ao2mo.general(cc.mol, (orbo,mo_coeff,mo_coeff,mo_coeff),
                                    feri, verbose=log)
                cput1 = log.timer_debug1('transforming oppp', *cput1)
                eri1 = feri['eri_mo']
                outbuf = numpy.empty((nmo,nmo,nmo))
                for i in range(nocc):
                    buf = _ccsd.unpack_tril(_cp(eri1[i*nmo:(i+1)*nmo]), out=outbuf)
                    self.oooo[i] = buf[:nocc,:nocc,:nocc]
                    self.ooov[i] = buf[:nocc,:nocc,nocc:]
                    self.ovoo[i] = buf[nocc:,:nocc,:nocc]
                    self.oovv[i] = buf[:nocc,nocc:,nocc:]
                    self.ovov[i] = buf[nocc:,:nocc,nocc:]
#.........这里部分代码省略.........

#.........这里部分代码省略.........

    blksize = 4
    d_ovov = numpy.empty((nocc,nvir,nocc,nvir))
    for p0, p1 in prange(0, nocc, blksize):
        d_ovov[p0:p1] = _cp(dovov[p0:p1])
        d_ovvo = _cp(dovvo[p0:p1])
        for i in range(p0,p1):
            d_ovov[i] += d_ovvo[i-p0].transpose(0,2,1)
    d_ovvo = None
    #:d_ovov = d_ovov + d_ovov.transpose(2,3,0,1)
    lib.transpose_sum(d_ovov.reshape(nov,nov), inplace=True)
    #:Ivo += numpy.einsum('jbka,jbki->ai', d_ovov, eris.ovoo)
    Ivo += lib.dot(d_ovov.reshape(-1,nvir).T, _cp(eris.ovoo).reshape(-1,nocc))
    #:Ioo += numpy.einsum('jakb,iakb->ij', d_ovov, eris.ovov)
    #:Ivv += numpy.einsum('jcib,jcia->ab', d_ovov, eris.ovov)
    Ioo += lib.dot(eris_ovov.reshape(nocc,-1), d_ovov.reshape(nocc,-1).T)
    Ivv += lib.dot(eris_ovov.reshape(-1,nvir).T, d_ovov.reshape(-1,nvir))

    nvir_pair = nvir * (nvir+1) // 2
    bufe_ovvv = numpy.empty((blksize,nvir,nvir,nvir))
    bufc_ovvv = numpy.empty((blksize,nvir,nvir_pair))
    bufc_ovvv.data = bufe_ovvv.data
    c_vvvo = numpy.empty((nvir_pair,nvir,nocc))
    for p0, p1 in prange(0, nocc, blksize):
        d_ovvv = numpy.empty((p1-p0,nvir,nvir,nvir))
        #:Ivo += numpy.einsum('jadc,jidc->ai', d_ovvv, eris_oovv)
        for i in range(p1-p0):
            lib.dot(dovvv[p0+i].reshape(nvir,-1),
                    eris_oovv[p0+i].reshape(nocc,-1).T, 1, Ivo, 1)

        c_ovvv = bufc_ovvv[:p1-p0]
        # tril part of (d_ovvv + d_ovvv.transpose(0,1,3,2))
        _ccsd.precontract(dovvv[p0:p1].reshape(-1,nvir,nvir), out=c_ovvv)
        for i0, i1, in prange(0, nvir_pair, BLKSIZE):
            for j0, j1 in prange(0, nvir, BLKSIZE//(p1-p0)+1):
                c_vvvo[i0:i1,j0:j1,p0:p1] = c_ovvv[:,j0:j1,i0:i1].transpose(2,1,0)
        eris_ovx = _cp(eris.ovvv[p0:p1])
        #:Xvo += numpy.einsum('jibc,jabc->ai', d_oovv, eris_ovvv)
        #:Ivv += numpy.einsum('ibdc,iadc->ab', d_ovvv, eris_ovvv)
        for i in range(p1-p0):
            lib.dot(eris_ovx[i].reshape(nvir,-1),
                    d_oovv[p0+i].reshape(nocc,-1).T, 1, Xvo, 1)
            lib.dot(eris_ovx[i].reshape(nvir,-1),
                    c_ovvv[i].reshape(nvir,-1).T, 1, Ivv, 1)

        eris_ovvv = bufe_ovvv[:p1-p0]
        _ccsd.unpack_tril(eris_ovx.reshape(-1,nvir_pair),
                          out=eris_ovvv.reshape(-1,nvir**2))
        eris_ovx = None
        #:Xvo += numpy.einsum('icjb,acjb->ai', d_ovov, eris_vvov)
        d_ovvo = _cp(d_ovov[p0:p1].transpose(0,1,3,2))
        lib.dot(eris_ovvv.reshape(-1,nvir).T, d_ovvo.reshape(-1,nocc), 1, Xvo, 1)

        e_ovvo, d_ovvo = d_ovvo, None
        for i in range(p1-p0):
            d_ovvv[i] = _ccsd.sum021(dovvv[p0+i])
            e_ovvo[i] = eris_ovov[p0+i].transpose(0,2,1)
        #:Ivo += numpy.einsum('jcab,jcib->ai', d_ovvv, eris_ovov)
        #:Ivv += numpy.einsum('icdb,icda->ab', d_ovvv, eris_ovvv)
        lib.dot(d_ovvv.reshape(-1,nvir).T,
                e_ovvo[:p1-p0].reshape(-1,nocc), 1, Ivo, 1)
        lib.dot(eris_ovvv.reshape(-1,nvir).T, d_ovvv.reshape(-1,nvir), 1, Ivv, 1)

        Xvo[:,p0:p1] +=(numpy.einsum('cb,iacb->ai', dvv, eris_ovvv) * 4
                      - numpy.einsum('cb,icba->ai', dvv+dvv.T, eris_ovvv))
    d_oovv = d_ovvv = bufc_ovvv = bufe_ovvv = None
    eris_ovov = eris_ovvv = eris_oovv = e_ovvo = None

    eris_ovvv = _cp(eris.ovvv)
    bufe_vvvo = numpy.empty((blksize*nvir,nvir,nocc))
    bufe_vvvv = numpy.empty((blksize*nvir,nvir,nvir))
    bufd_vvvv = numpy.empty((blksize*nvir,nvir,nvir))
    for p0, p1 in prange(0, nvir, blksize):
        off0 = p0*(p0+1)//2
        off1 = p1*(p1+1)//2
        d_vvvv = _cp(dvvvv[off0:off1]) * 4
        for i in range(p0, p1):
            d_vvvv[i*(i+1)//2+i-off0] *= .5
        d_vvvv = _ccsd.unpack_tril(d_vvvv, out=bufd_vvvv[:off1-off0])
        eris_vvvv = _ccsd.unpack_tril(eris.vvvv[off0:off1], out=bufe_vvvv[:off1-off0])
        #:Ivv += numpy.einsum('decb,deca->ab', d_vvvv, eris_vvvv) * 2
        #:Xvo += numpy.einsum('icdb,acdb->ai', d_ovvv, eris_vvvv)
        lib.dot(eris_vvvv.reshape(-1,nvir).T, d_vvvv.reshape(-1,nvir), 2, Ivv, 1)
        d_vvvo = _cp(c_vvvo[off0:off1])
        lib.dot(eris_vvvv.reshape(-1,nvir).T, d_vvvo.reshape(-1,nocc), 1, Xvo, 1)

        #:Ioo += numpy.einsum('abjc,abci->ij', d_vvov, eris_vvvo)
        #:Ivo += numpy.einsum('dbca,dbci->ai', d_vvvv, eris_vvvo) * 2
        eris_vvvo = bufe_vvvo[:off1-off0]
        for i0, i1 in prange(off0, off1, BLKSIZE):
            for j0, j1, in prange(0, nvir, BLKSIZE//nocc+1):
                eris_vvvo[i0-off0:i1-off0,j0:j1,:] = eris_ovvv[:,j0:j1,i0:i1].transpose(2,1,0)
        lib.dot(eris_vvvo.reshape(-1,nocc).T, d_vvvo.reshape(-1,nocc), 1, Ioo, 1)
        lib.dot(d_vvvv.reshape(-1,nvir).T, eris_vvvo.reshape(-1,nocc), 2, Ivo, 1)

    Ioo *= -1
    Ivv *= -1
    Ivo *= -1
    Xvo += Ivo
    return Ioo, Ivv, Ivo, Xvo