def _sort_t2_vooo(mycc, orbsym, t1, t2, eris):
    ovoo = numpy.asarray(eris.ovoo)
    nocc, nvir = t1.shape
    if mycc.mol.symmetry:
        orbsym = numpy.asarray(orbsym, dtype=numpy.int32)
        o_sorted = _irrep_argsort(orbsym[:nocc])
        v_sorted = _irrep_argsort(orbsym[nocc:])
        mo_energy = eris.fock.diagonal()
        mo_energy = numpy.hstack((mo_energy[:nocc][o_sorted],
                                  mo_energy[nocc:][v_sorted]))
        t1T = numpy.asarray(t1.T[v_sorted][:,o_sorted], order='C')

        t2T = lib.transpose(t2.reshape(nocc**2,-1))
        _ccsd.libcc.CCsd_t_sort_t2(t2.ctypes.data_as(ctypes.c_void_p),
                                   t2T.ctypes.data_as(ctypes.c_void_p),
                                   orbsym.ctypes.data_as(ctypes.c_void_p),
                                   ctypes.c_int(nocc), ctypes.c_int(nvir))
        t2T = t2.reshape(nvir,nvir,nocc,nocc)
        vooo = numpy.empty((nvir,nocc,nocc,nocc))
        _ccsd.libcc.CCsd_t_sort_vooo(vooo.ctypes.data_as(ctypes.c_void_p),
                                     ovoo.ctypes.data_as(ctypes.c_void_p),
                                     orbsym.ctypes.data_as(ctypes.c_void_p),
                                     ctypes.c_int(nocc), ctypes.c_int(nvir))
        ovoo = None
    else:
        t1T = t1.T.copy()
        t2T = lib.transpose(t2.reshape(nocc**2,-1))
        t2T = lib.transpose(t2T.reshape(-1,nocc,nocc), axes=(0,2,1), out=t2)
        t2T = t2T.reshape(nvir,nvir,nocc,nocc)
        vooo = ovoo.transpose(1,0,2,3).copy()
        mo_energy = numpy.asarray(eris.fock.diagonal(), order='C')
    return mo_energy, t1T, t2T, vooo

def cholesky_eri(mol, auxbasis='weigend+etb', aosym='s1', verbose=0):
    t0 = (time.clock(), time.time())
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mol.stdout, verbose)
    auxmol = incore.format_aux_basis(mol, auxbasis)

    j2c = incore.fill_2c2e(mol, auxmol)
    log.debug('size of aux basis %d', j2c.shape[0])
    t1 = log.timer('2c2e', *t0)
    low = scipy.linalg.cholesky(j2c, lower=True)
    j2c = None
    t1 = log.timer('Cholesky 2c2e', *t1)

    j3c_ll = aux_e2(mol, auxmol, intor='cint3c2e_spinor', aosym=aosym)
    j3c_ss = aux_e2(mol, auxmol, intor='cint3c2e_spsp1_spinor', aosym=aosym)
    t1 = log.timer('3c2e', *t1)
    cderi_ll = scipy.linalg.solve_triangular(low, j3c_ll.T, lower=True,
                                             overwrite_b=True)
    cderi_ss = scipy.linalg.solve_triangular(low, j3c_ss.T, lower=True,
                                             overwrite_b=True)
    # solve_triangular return cderi in Fortran order
    cderi = (lib.transpose(cderi_ll.T), lib.transpose(cderi_ss.T))
    log.timer('cholesky_eri', *t0)
    return cderi

def eval_ao_kpts(cell, coords, kpts=None, deriv=0, relativity=0,
                 shl_slice=None, non0tab=None, out=None, verbose=None, **kwargs):
    '''
    Returns:
        ao_kpts: (nkpts, ngs, nao) ndarray
            AO values at each k-point
    '''
    if kpts is None:
        if 'kpt' in kwargs:
            sys.stderr.write('WARN: _KNumInt.eval_ao function finds keyword '
                             'argument "kpt" and converts it to "kpts"\n')
            kpts = kpt
        else:
            kpts = numpy.zeros((1,3))
    kpts = numpy.reshape(kpts, (-1,3))
    nkpts = len(kpts)
    ngrids = len(coords)

    if non0tab is None:
        non0tab = numpy.ones(((ngrids+BLKSIZE-1)//BLKSIZE, cell.nbas),
                             dtype=numpy.int8)

    ao_loc = cell.ao_loc_nr()
    nao = ao_loc[-1]
    comp = (deriv+1)*(deriv+2)*(deriv+3)//6
    ao_kpts = [numpy.zeros((ngrids,nao,comp), dtype=numpy.complex128, order='F')
               for k in range(nkpts)]
    out_ptrs = (ctypes.c_void_p*nkpts)(
            *[x.ctypes.data_as(ctypes.c_void_p) for x in ao_kpts])
    coords = numpy.asarray(coords, order='C')
    Ls = cell.get_lattice_Ls()
    expLk = numpy.exp(1j * numpy.asarray(numpy.dot(Ls, kpts.T), order='C'))

    drv = getattr(libpbc, 'PBCval_sph_deriv%d' % deriv)
    drv(ctypes.c_int(ngrids), ctypes.c_int(BLKSIZE),
        Ls.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(len(Ls)),
        expLk.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(nkpts),
        (ctypes.c_int*2)(0, cell.nbas),
        ao_loc.ctypes.data_as(ctypes.c_void_p),
        out_ptrs, coords.ctypes.data_as(ctypes.c_void_p),
        non0tab.ctypes.data_as(ctypes.c_void_p),
        cell._atm.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(cell.natm),
        cell._bas.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(cell.nbas),
        cell._env.ctypes.data_as(ctypes.c_void_p))

    for k, kpt in enumerate(kpts):
        mat = ao_kpts[k].transpose(2,0,1)
        if comp == 1:
            aos = lib.transpose(mat[0].T)
        else:
            aos = numpy.empty((comp,ngrids,nao), dtype=numpy.complex128)
            for i in range(comp):
                lib.transpose(mat[i].T, out=aos[i])

        if abs(kpt).sum() < 1e-9:  # gamma point
            aos = aos.real.copy()

        ao_kpts[k] = aos
    return ao_kpts

 def test_transpose(self):
     a = numpy.random.random((400,900))
     self.assertAlmostEqual(abs(a.T - lib.transpose(a)).max(), 0, 12)
     b = a[:400,:400]
     c = numpy.copy(b)
     self.assertAlmostEqual(abs(b.T - lib.transpose(c,inplace=True)).max(), 0, 12)
     a = a.reshape(40,10,-1)
     self.assertAlmostEqual(abs(a.transpose(0,2,1) - lib.transpose(a,(0,2,1))).max(), 0, 12)

def contract_2e(eri, civec_strs, norb, nelec, link_index=None):
    ci_coeff, nelec, ci_strs = _unpack(civec_strs, nelec)
    if link_index is None:
        link_index = _all_linkstr_index(ci_strs, norb, nelec)
    cd_indexa, dd_indexa, cd_indexb, dd_indexb = link_index
    na, nlinka = cd_indexa.shape[:2]
    nb, nlinkb = cd_indexb.shape[:2]

    eri = ao2mo.restore(1, eri, norb)
    eri1 = eri.transpose(0,2,1,3) - eri.transpose(0,2,3,1)
    idx,idy = numpy.tril_indices(norb, -1)
    idx = idx * norb + idy
    eri1 = lib.take_2d(eri1.reshape(norb**2,-1), idx, idx) * 2
    fcivec = ci_coeff.reshape(na,nb)
    # (bb|bb)
    if nelec[1] > 1:
        mb, mlinkb = dd_indexb.shape[:2]
        fcivecT = lib.transpose(fcivec)
        ci1T = numpy.zeros((nb,na))
        libfci.SCIcontract_2e_aaaa(eri1.ctypes.data_as(ctypes.c_void_p),
                                   fcivecT.ctypes.data_as(ctypes.c_void_p),
                                   ci1T.ctypes.data_as(ctypes.c_void_p),
                                   ctypes.c_int(norb),
                                   ctypes.c_int(nb), ctypes.c_int(na),
                                   ctypes.c_int(mb), ctypes.c_int(mlinkb),
                                   dd_indexb.ctypes.data_as(ctypes.c_void_p))
        ci1 = lib.transpose(ci1T, out=fcivecT)
    else:
        ci1 = numpy.zeros_like(fcivec)
    # (aa|aa)
    if nelec[0] > 1:
        ma, mlinka = dd_indexa.shape[:2]
        libfci.SCIcontract_2e_aaaa(eri1.ctypes.data_as(ctypes.c_void_p),
                                   fcivec.ctypes.data_as(ctypes.c_void_p),
                                   ci1.ctypes.data_as(ctypes.c_void_p),
                                   ctypes.c_int(norb),
                                   ctypes.c_int(na), ctypes.c_int(nb),
                                   ctypes.c_int(ma), ctypes.c_int(mlinka),
                                   dd_indexa.ctypes.data_as(ctypes.c_void_p))

    h_ps = numpy.einsum('pqqs->ps', eri)
    eri1 = eri * 2
    for k in range(norb):
        eri1[:,:,k,k] += h_ps/nelec[0]
        eri1[k,k,:,:] += h_ps/nelec[1]
    eri1 = ao2mo.restore(4, eri1, norb)
    # (bb|aa)
    libfci.SCIcontract_2e_bbaa(eri1.ctypes.data_as(ctypes.c_void_p),
                               fcivec.ctypes.data_as(ctypes.c_void_p),
                               ci1.ctypes.data_as(ctypes.c_void_p),
                               ctypes.c_int(norb),
                               ctypes.c_int(na), ctypes.c_int(nb),
                               ctypes.c_int(nlinka), ctypes.c_int(nlinkb),
                               cd_indexa.ctypes.data_as(ctypes.c_void_p),
                               cd_indexb.ctypes.data_as(ctypes.c_void_p))

    return _as_SCIvector(ci1.reshape(ci_coeff.shape), ci_strs)

    def make_kpt(kpt):  # kpt = kptj - kpti
        # search for all possible ki and kj that has ki-kj+kpt=0
        kk_match = numpy.einsum('ijx->ij', abs(kk_table + kpt)) < 1e-9
        kpti_idx, kptj_idx = numpy.where(kk_todo & kk_match)
        nkptj = len(kptj_idx)
        log.debug1('kpt = %s', kpt)
        log.debug2('kpti_idx = %s', kpti_idx)
        log.debug2('kptj_idx = %s', kptj_idx)
        kk_todo[kpti_idx,kptj_idx] = False
        if swap_2e and not is_zero(kpt):
            kk_todo[kptj_idx,kpti_idx] = False

        max_memory1 = max_memory * (nkptj+1)/(nkptj+5)
        blksize = max(int(max_memory1*4e6/(nkptj+5)/16/nao**2), 16)
        bufR = numpy.empty((blksize*nao**2))
        bufI = numpy.empty((blksize*nao**2))
        # Use DF object to mimic KRHF/KUHF object in function get_coulG
        mydf.exxdiv = exxdiv
        vkcoulG = mydf.weighted_coulG(kpt, True, mydf.gs)
        kptjs = kpts[kptj_idx]
        # <r|-G+k_rs|s> = conj(<s|G-k_rs|r>) = conj(<s|G+k_sr|r>)
        for k, pqkR, pqkI, p0, p1 \
                in mydf.ft_loop(mydf.gs, kpt, kptjs, max_memory=max_memory1):
            ki = kpti_idx[k]
            kj = kptj_idx[k]
            coulG = numpy.sqrt(vkcoulG[p0:p1])

# case 1: k_pq = (pi|iq)
#:v4 = numpy.einsum('ijL,lkL->ijkl', pqk, pqk.conj())
#:vk += numpy.einsum('ijkl,jk->il', v4, dm)
            pqkR *= coulG
            pqkI *= coulG
            pLqR = lib.transpose(pqkR.reshape(nao,nao,-1), axes=(0,2,1), out=bufR)
            pLqI = lib.transpose(pqkI.reshape(nao,nao,-1), axes=(0,2,1), out=bufI)
            iLkR = numpy.empty((nao*(p1-p0),nao))
            iLkI = numpy.empty((nao*(p1-p0),nao))
            for i in range(nset):
                iLkR, iLkI = zdotNN(pLqR.reshape(-1,nao), pLqI.reshape(-1,nao),
                                    dmsR[i,kj], dmsI[i,kj], 1, iLkR, iLkI)
                zdotNC(iLkR.reshape(nao,-1), iLkI.reshape(nao,-1),
                       pLqR.reshape(nao,-1).T, pLqI.reshape(nao,-1).T,
                       1, vkR[i,ki], vkI[i,ki], 1)

# case 2: k_pq = (iq|pi)
#:v4 = numpy.einsum('iLj,lLk->ijkl', pqk, pqk.conj())
#:vk += numpy.einsum('ijkl,li->kj', v4, dm)
            if swap_2e and not is_zero(kpt):
                iLkR = iLkR.reshape(nao,-1)
                iLkI = iLkI.reshape(nao,-1)
                for i in range(nset):
                    iLkR, iLkI = zdotNN(dmsR[i,ki], dmsI[i,ki], pLqR.reshape(nao,-1),
                                        pLqI.reshape(nao,-1), 1, iLkR, iLkI)
                    zdotCN(pLqR.reshape(-1,nao).T, pLqI.reshape(-1,nao).T,
                           iLkR.reshape(-1,nao), iLkI.reshape(-1,nao),
                           1, vkR[i,kj], vkI[i,kj], 1)
            pqkR = pqkI = coulG = pLqR = pLqI = iLkR = iLkI = None

 def loop(self):
     Lpq = None
     coulG = self.weighted_coulG()
     for pqkR, pqkI, p0, p1 in self.pw_loop(aosym='s2', blksize=self.blockdim):
         vG = numpy.sqrt(coulG[p0:p1])
         pqkR *= vG
         pqkI *= vG
         Lpq = lib.transpose(pqkR, out=Lpq)
         yield Lpq
         Lpq = lib.transpose(pqkI, out=Lpq)
         yield Lpq

def _sort_t2_vooo_(mycc, orbsym, t1, t2, eris):
    assert(t2.flags.c_contiguous)
    vooo = numpy.asarray(eris.ovoo).transpose(1,0,3,2).conj().copy()
    nocc, nvir = t1.shape
    if mycc.mol.symmetry:
        orbsym = numpy.asarray(orbsym, dtype=numpy.int32)
        o_sorted = _irrep_argsort(orbsym[:nocc])
        v_sorted = _irrep_argsort(orbsym[nocc:])
        mo_energy = eris.mo_energy
        mo_energy = numpy.hstack((mo_energy[:nocc][o_sorted],
                                  mo_energy[nocc:][v_sorted]))
        t1T = numpy.asarray(t1.T[v_sorted][:,o_sorted], order='C')
        fvo = eris.fock[nocc:,:nocc]
        fvo = numpy.asarray(fvo[v_sorted][:,o_sorted], order='C')

        o_sym = orbsym[o_sorted]
        oo_sym = (o_sym[:,None] ^ o_sym).ravel()
        oo_sorted = _irrep_argsort(oo_sym)
        #:vooo = eris.ovoo.transpose(1,0,2,3)
        #:vooo = vooo[v_sorted][:,o_sorted][:,:,o_sorted][:,:,:,o_sorted]
        #:vooo = vooo.reshape(nvir,-1,nocc)[:,oo_sorted]
        oo_idx = numpy.arange(nocc**2).reshape(nocc,nocc)[o_sorted][:,o_sorted]
        oo_idx = oo_idx.ravel()[oo_sorted]
        oo_idx = (oo_idx[:,None]*nocc+o_sorted).ravel()
        vooo = lib.take_2d(vooo.reshape(nvir,-1), v_sorted, oo_idx)
        vooo = vooo.reshape(nvir,nocc,nocc,nocc)

        #:t2T = t2.transpose(2,3,1,0)
        #:t2T = ref_t2T[v_sorted][:,v_sorted][:,:,o_sorted][:,:,:,o_sorted]
        #:t2T = ref_t2T.reshape(nvir,nvir,-1)[:,:,oo_sorted]
        t2T = lib.transpose(t2.reshape(nocc**2,-1))
        oo_idx = numpy.arange(nocc**2).reshape(nocc,nocc).T[o_sorted][:,o_sorted]
        oo_idx = oo_idx.ravel()[oo_sorted]
        vv_idx = (v_sorted[:,None]*nvir+v_sorted).ravel()
        t2T = lib.take_2d(t2T.reshape(nvir**2,nocc**2), vv_idx, oo_idx, out=t2)
        t2T = t2T.reshape(nvir,nvir,nocc,nocc)
        def restore_t2_inplace(t2T):
            tmp = numpy.zeros((nvir**2,nocc**2), dtype=t2T.dtype)
            lib.takebak_2d(tmp, t2T.reshape(nvir**2,nocc**2), vv_idx, oo_idx)
            t2 = lib.transpose(tmp.reshape(nvir**2,nocc**2), out=t2T)
            return t2.reshape(nocc,nocc,nvir,nvir)
    else:
        fvo = eris.fock[nocc:,:nocc].copy()
        t1T = t1.T.copy()
        t2T = lib.transpose(t2.reshape(nocc**2,nvir**2))
        t2T = lib.transpose(t2T.reshape(nvir**2,nocc,nocc), axes=(0,2,1), out=t2)
        mo_energy = numpy.asarray(eris.mo_energy, order='C')
        def restore_t2_inplace(t2T):
            tmp = lib.transpose(t2T.reshape(nvir**2,nocc,nocc), axes=(0,2,1))
            t2 = lib.transpose(tmp.reshape(nvir**2,nocc**2), out=t2T)
            return t2.reshape(nocc,nocc,nvir,nvir)
    t2T = t2T.reshape(nvir,nvir,nocc,nocc)
    return mo_energy, t1T, t2T, vooo, fvo, restore_t2_inplace

 def get_Lpq(shls_slice, col0, col1, buf):
     # Be cautious here, _ri.nr_auxe2 assumes buf in F-order
     Lpq = _ri.nr_auxe2(intor, atm, bas, env, shls_slice, ao_loc,
                        's2ij', 1, cintopt, buf).T
     if mydf.charge_constraint:
         Lpq = numpy.ndarray(shape=(naux+1,col1-col0), buffer=buf)
         Lpq[naux,:] = ovlp[col0:col1]
         Lpq1 = solve(Lpq)
         assert(Lpq1.flags.f_contiguous)
         lib.transpose(Lpq1.T, out=Lpq)
         return Lpq[:naux]
     else:
         return solve(Lpq)

def trans_e1_incore(eri_ao, mo, ncore, ncas):
    nmo = mo.shape[1]
    nocc = ncore + ncas
    eri1 = ao2mo.incore.half_e1(eri_ao, (mo,mo[:,:nocc]), compact=False)
    eri1 = eri1.reshape(nmo,nocc,-1)

    klppshape = (0, nmo, 0, nmo)
    klpashape = (0, nmo, ncore, nocc)
    aapp = numpy.empty((ncas,ncas,nmo,nmo))
    for i in range(ncas):
        _ao2mo.nr_e2(eri1[ncore+i,ncore:nocc], mo, klppshape,
                      aosym='s4', mosym='s1', out=aapp[i])
    ppaa = lib.transpose(aapp.reshape(ncas*ncas,-1)).reshape(nmo,nmo,ncas,ncas)
    aapp = None

    papa = numpy.empty((nmo,ncas,nmo,ncas))
    for i in range(nmo):
        _ao2mo.nr_e2(eri1[i,ncore:nocc], mo, klpashape,
                      aosym='s4', mosym='s1', out=papa[i])

    pp = numpy.empty((nmo,nmo))
    j_cp = numpy.zeros((ncore,nmo))
    k_pc = numpy.zeros((nmo,ncore))
    for i in range(ncore):
        _ao2mo.nr_e2(eri1[i,i:i+1], mo, klppshape, aosym='s4', mosym='s1', out=pp)
        j_cp[i] = pp.diagonal()
    j_pc = j_cp.T.copy()

    pp = numpy.empty((ncore,ncore))
    for i in range(nmo):
        klshape = (i, i+1, 0, ncore)
        _ao2mo.nr_e2(eri1[i,:ncore], mo, klshape, aosym='s4', mosym='s1', out=pp)
        k_pc[i] = pp.diagonal()
    return j_pc, k_pc, ppaa, papa

    def __getitem__(self, s):
        dat = self.dat
        if isinstance(dat, h5py.Group):
            v = numpy.hstack([dat[str(i)][s] for i in range(len(dat))])
        else: # For mpi4pyscf, pyscf-1.5.1 or older
            v = numpy.asarray(dat[s])

        nao = int(numpy.sqrt(v.shape[-1]))
        v1 = lib.transpose(v.reshape(-1,nao,nao), axes=(0,2,1)).conj()
        return v1.reshape(v.shape)

    def loop(self, blksize=None):
        cell = self.cell
        if cell.dimension == 2 and cell.low_dim_ft_type != 'inf_vacuum':
            raise RuntimeError('ERIs of PBC-2D systems are not positive '
                               'definite. Current API only supports postive '
                               'definite ERIs.')

        if blksize is None:
            blksize = self.blockdim
        # coulG of 1D and 2D has negative elements.
        coulG = self.weighted_coulG()
        Lpq = None
        for pqkR, pqkI, p0, p1 in self.pw_loop(aosym='s2', blksize=blksize):
            vG = numpy.sqrt(coulG[p0:p1])
            pqkR *= vG
            pqkI *= vG
            Lpq = lib.transpose(pqkR, out=Lpq)
            yield Lpq
            Lpq = lib.transpose(pqkI, out=Lpq)
            yield Lpq

 def store(b, label):
     if b.ndim == 3 and b.flags.f_contiguous:
         b = lib.transpose(b.T, axes=(0,2,1)).reshape(naux,-1)
     else:
         b = b.reshape((-1,naux)).T
     if tag == 'cd':
         cderi = scipy.linalg.solve_triangular(low, b, lower=True,
                                               overwrite_b=True)
     else:
         cderi = lib.dot(low.T, b)
     feri[label] = cderi

    def gen_int3c(job_id, ish0, ish1):
        dataname = 'j3c-chunks/%d' % job_id
        i0 = ao_loc[ish0]
        i1 = ao_loc[ish1]
        dii = i1*(i1+1)//2 - i0*(i0+1)//2
        dij = (i1 - i0) * nao
        if j_only:
            buflen = max(8, int(max_memory*1e6/16/(nkptij*dii+dii)))
        else:
            buflen = max(8, int(max_memory*1e6/16/(nkptij*dij+dij)))
        auxranges = balance_segs(aux_loc[1:]-aux_loc[:-1], buflen)
        buflen = max([x[2] for x in auxranges])
        buf = numpy.empty(nkptij*dij*buflen, dtype=dtype)
        buf1 = numpy.empty(dij*buflen, dtype=dtype)

        naux = aux_loc[-1]
        for kpt_id, kptij in enumerate(kptij_lst):
            key = '%s/%d' % (dataname, kpt_id)
            if aosym_s2[kpt_id]:
                shape = (naux, dii)
            else:
                shape = (naux, dij)
            if gamma_point(kptij):
                fswap.create_dataset(key, shape, 'f8')
            else:
                fswap.create_dataset(key, shape, 'c16')

        naux0 = 0
        for istep, auxrange in enumerate(auxranges):
            log.alldebug2("aux_e1 job_id %d step %d", job_id, istep)
            sh0, sh1, nrow = auxrange
            sub_slice = (ish0, ish1, 0, cell.nbas, sh0, sh1)
            if j_only:
                mat = numpy.ndarray((nkptij,dii,nrow), dtype=dtype, buffer=buf)
            else:
                mat = numpy.ndarray((nkptij,dij,nrow), dtype=dtype, buffer=buf)
            mat = int3c(sub_slice, mat)

            for k, kptij in enumerate(kptij_lst):
                h5dat = fswap['%s/%d'%(dataname,k)]
                v = lib.transpose(mat[k], out=buf1)
                if not j_only and aosym_s2[k]:
                    idy = idxb[i0*(i0+1)//2:i1*(i1+1)//2] - i0 * nao
                    out = numpy.ndarray((nrow,dii), dtype=v.dtype, buffer=mat[k])
                    v = numpy.take(v, idy, axis=1, out=out)
                if gamma_point(kptij):
                    h5dat[naux0:naux0+nrow] = v.real
                else:
                    h5dat[naux0:naux0+nrow] = v
            naux0 += nrow

def make_rdm2s(civec_strs, norb, nelec, link_index=None, **kwargs):
    r'''Spin separated 2-particle density matrices.
    The return values include three density matrices:
    (alpha,alpha,alpha,alpha), (alpha,alpha,beta,beta), (beta,beta,beta,beta)

    2pdm[p,q,r,s] = :math:`\langle p^\dagger r^\dagger s q\rangle`
    '''
    ci_coeff, nelec, ci_strs = _unpack(civec_strs, nelec)
    if link_index is None:
        cd_indexa = cre_des_linkstr(ci_strs[0], norb, nelec[0])
        dd_indexa = des_des_linkstr(ci_strs[0], norb, nelec[0])
        cd_indexb = cre_des_linkstr(ci_strs[1], norb, nelec[1])
        dd_indexb = des_des_linkstr(ci_strs[1], norb, nelec[1])
    else:
        cd_indexa, dd_indexa, cd_indexb, dd_indexb = link_index
    na, nlinka = cd_indexa.shape[:2]
    nb, nlinkb = cd_indexb.shape[:2]

    fcivec = ci_coeff.reshape(na,nb)
    # (bb|aa) and (aa|bb)
    dm2ab = rdm.make_rdm12_spin1('FCItdm12kern_ab', fcivec, fcivec,
                                 norb, nelec, (cd_indexa,cd_indexb), 0)[1]
    # (aa|aa)
    dm2aa = numpy.zeros([norb]*4)
    if nelec[0] > 1:
        ma, mlinka = dd_indexa.shape[:2]
        libfci.SCIrdm2_aaaa(libfci.SCIrdm2kern_aaaa,
                            dm2aa.ctypes.data_as(ctypes.c_void_p),
                            fcivec.ctypes.data_as(ctypes.c_void_p),
                            fcivec.ctypes.data_as(ctypes.c_void_p),
                            ctypes.c_int(norb),
                            ctypes.c_int(na), ctypes.c_int(nb),
                            ctypes.c_int(ma), ctypes.c_int(mlinka),
                            dd_indexa.ctypes.data_as(ctypes.c_void_p))
    # (bb|bb)
    dm2bb = numpy.zeros([norb]*4)
    if nelec[1] > 1:
        mb, mlinkb = dd_indexb.shape[:2]
        fcivecT = lib.transpose(fcivec)
        libfci.SCIrdm2_aaaa(libfci.SCIrdm2kern_aaaa,
                            dm2bb.ctypes.data_as(ctypes.c_void_p),
                            fcivecT.ctypes.data_as(ctypes.c_void_p),
                            fcivecT.ctypes.data_as(ctypes.c_void_p),
                            ctypes.c_int(norb),
                            ctypes.c_int(nb), ctypes.c_int(na),
                            ctypes.c_int(mb), ctypes.c_int(mlinkb),
                            dd_indexb.ctypes.data_as(ctypes.c_void_p))
    return dm2aa, dm2ab, dm2bb

def _trans(mo, ncore, ncas, fload, cvcv=None, ao_loc=None):
    nao, nmo = mo.shape
    nocc = ncore + ncas
    nvir = nmo - nocc
    nav = nmo - ncore

    if cvcv is None:
        cvcv = numpy.zeros((ncore*nvir,ncore*nvir))
    pacv = numpy.empty((nmo,ncas,ncore*nvir))
    aapp = numpy.empty((ncas,ncas,nmo*nmo))
    papa = numpy.empty((nmo,ncas,nmo*ncas))
    vcv = numpy.empty((nav,ncore*nvir))
    apa = numpy.empty((ncas,nmo*ncas))
    vpa = numpy.empty((nav,nmo*ncas))
    app = numpy.empty((ncas,nmo*nmo))
    for i in range(ncore):
        buf = fload(i, i+1)
        klshape = (0, ncore, nocc, nmo)
        _ao2mo.nr_e2(buf, mo, klshape,
                      aosym='s4', mosym='s1', out=vcv, ao_loc=ao_loc)
        cvcv[i*nvir:(i+1)*nvir] = vcv[ncas:]
        pacv[i] = vcv[:ncas]

        klshape = (0, nmo, ncore, nocc)
        _ao2mo.nr_e2(buf[:ncas], mo, klshape,
                      aosym='s4', mosym='s1', out=apa, ao_loc=ao_loc)
        papa[i] = apa
    for i in range(ncas):
        buf = fload(ncore+i, ncore+i+1)
        klshape = (0, ncore, nocc, nmo)
        _ao2mo.nr_e2(buf, mo, klshape,
                      aosym='s4', mosym='s1', out=vcv, ao_loc=ao_loc)
        pacv[ncore:,i] = vcv

        klshape = (0, nmo, ncore, nocc)
        _ao2mo.nr_e2(buf, mo, klshape,
                      aosym='s4', mosym='s1', out=vpa, ao_loc=ao_loc)
        papa[ncore:,i] = vpa

        klshape = (0, nmo, 0, nmo)
        _ao2mo.nr_e2(buf[:ncas], mo, klshape,
                      aosym='s4', mosym='s1', out=app, ao_loc=ao_loc)
        aapp[i] = app
    #lib.transpose(aapp.reshape(ncas**2, -1), inplace=True)
    ppaa = lib.transpose(aapp.reshape(ncas**2,-1))
    return (ppaa.reshape(nmo,nmo,ncas,ncas), papa.reshape(nmo,ncas,nmo,ncas),
            pacv.reshape(nmo,ncas,ncore,nvir), cvcv)

    def compress(Lpq):
        ncol = Lpq.shape[-1]
        auxchgs = numpy.zeros((naochg,ncol), dtype=Lpq.dtype)
        for i in range(auxcell.nbas):
            l  = auxcell.bas_angular(i)
            ia = auxcell.bas_atom(i)
            p0 = modchg_offset[ia,l]
            if p0 >= 0:
                nc = auxcell.bas_nctr(i)
                i0,i1 = aux_loc[i:i+2]
                lchg = numpy.einsum('imn->mn', Lpq[i0:i1].reshape(nc,-1,ncol))
                auxchgs[p0:p0+len(lchg)] -= lchg

        if Lpq.flags.f_contiguous:
            Lpq = lib.transpose(Lpq.T)
        Lpq = numpy.vstack((Lpq, auxchgs))
        return Lpq

def cholesky_eri(mol, auxbasis='weigend+etb', auxmol=None,
                 int3c='int3c2e', aosym='s2ij', int2c='int2c2e', comp=1,
                 verbose=0, fauxe2=aux_e2):
    '''
    Returns:
        2D array of (naux,nao*(nao+1)/2) in C-contiguous
    '''
    assert(comp == 1)
    t0 = (time.clock(), time.time())
    log = logger.new_logger(mol, verbose)
    if auxmol is None:
        auxmol = addons.make_auxmol(mol, auxbasis)

    j2c = auxmol.intor(int2c, hermi=1)
    naux = j2c.shape[0]
    log.debug('size of aux basis %d', naux)
    t1 = log.timer('2c2e', *t0)

    j3c = fauxe2(mol, auxmol, intor=int3c, aosym=aosym).reshape(-1,naux)
    t1 = log.timer('3c2e', *t1)

    try:
        low = scipy.linalg.cholesky(j2c, lower=True)
        j2c = None
        t1 = log.timer('Cholesky 2c2e', *t1)
        cderi = scipy.linalg.solve_triangular(low, j3c.T, lower=True,
                                              overwrite_b=True)
    except scipy.linalg.LinAlgError:
        w, v = scipy.linalg.eigh(j2c)
        idx = w > LINEAR_DEP_THR
        v = (v[:,idx] / numpy.sqrt(w[idx]))
        cderi = lib.dot(v.T, j3c.T)

    j3c = None
    if cderi.flags.f_contiguous:
        cderi = lib.transpose(cderi.T)
    log.timer('cholesky_eri', *t0)
    return cderi

    def __enter__(self):
        self.feri = h5py.File(self.cderi, 'r')
        kpti_kptj = numpy.asarray(self.kpti_kptj)
        kptij_lst = self.feri['%s-kptij'%self.label].value
        dk = numpy.einsum('kij->k', abs(kptij_lst-kpti_kptj))
        k_id = numpy.where(dk < 1e-6)[0]
        if len(k_id) > 0:
            dat = self.feri['%s/%d' % (self.label,k_id[0])]
        else:
            # swap ki,kj due to the hermiticity
            kptji = kpti_kptj[[1,0]]
            dk = numpy.einsum('kij->k', abs(kptij_lst-kptji))
            k_id = numpy.where(dk < 1e-6)[0]
            if len(k_id) == 0:
                raise RuntimeError('%s for kpts %s is not initialized.\n'
                                   'Reset attribute .kpts then call '
                                   '.build() to initialize %s.'
                                   % (self.label, kpti_kptj, self.label))
# FIXME: memory usage
            dat = self.feri['%s/%d' % (self.label, k_id[0])].value
            nao = int(numpy.sqrt(dat.shape[1]))
            dat = lib.transpose(dat.reshape(-1,nao,nao), axes=(0,2,1)).conj()
            dat = dat.reshape(-1,nao**2)
        return dat

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):
#.........这里部分代码省略.........