def writeIntegralFile(shciobj, h1eff, eri_cas, ncas, nelec, ecore=0):
    if isinstance(nelec, (int, numpy.integer)):
        if shciobj.spin is None:
            nelecb = nelec // 2
        else:
            nelecb = (nelec - shciobj.spin) // 2
        neleca = nelec - nelecb
    else :
        neleca, nelecb = nelec

    if shciobj.groupname is not None and shciobj.orbsym is not []:
# First removing the symmetry forbidden integrals. This has been done using
# the pyscf internal irrep-IDs (stored in shciobj.orbsym)
        orbsym = numpy.asarray(shciobj.orbsym) % 10
        pair_irrep = (orbsym.reshape(-1,1) ^ orbsym)[numpy.tril_indices(ncas)]
        sym_forbid = pair_irrep.reshape(-1,1) != pair_irrep.ravel()
        eri_cas = ao2mo.restore(4, eri_cas, ncas)
        eri_cas[sym_forbid] = 0
        eri_cas = ao2mo.restore(8, eri_cas, ncas)
        # Convert the pyscf internal irrep-ID to molpro irrep-ID
        orbsym = numpy.asarray(symmetry.convert_orbsym(shciobj.groupname, orbsym))
    else:
        orbsym = []
        eri_cas = ao2mo.restore(8, eri_cas, ncas)

    if not os.path.exists(shciobj.runtimedir):
        os.makedirs(shciobj.runtimedir)

    # The name of the FCIDUMP file, default is "FCIDUMP".
    integralFile = os.path.join(shciobj.runtimedir, shciobj.integralfile)
    tools.fcidump.from_integrals(integralFile, h1eff, eri_cas, ncas,
                                 neleca+nelecb, ecore, ms=abs(neleca-nelecb),
                                 orbsym=orbsym)
    return integralFile

def contract_2e(eri, fcivec, nsite, nelec, nphonon):
    neleca, nelecb = _unpack_nelec(nelec)
    link_indexa = cistring.gen_linkstr_index(range(nsite), neleca)
    link_indexb = cistring.gen_linkstr_index(range(nsite), nelecb)
    cishape = make_shape(nsite, nelec, nphonon)

    ci0 = fcivec.reshape(cishape)
    t1a = numpy.zeros((nsite,nsite)+cishape)
    t1b = numpy.zeros((nsite,nsite)+cishape)
    for str0, tab in enumerate(link_indexa):
        for a, i, str1, sign in tab:
            t1a[a,i,str1] += sign * ci0[str0]
    for str0, tab in enumerate(link_indexb):
        for a, i, str1, sign in tab:
            t1b[a,i,:,str1] += sign * ci0[:,str0]

    g2e_aa = ao2mo.restore(1, eri[0], nsite)
    g2e_ab = ao2mo.restore(1, eri[1], nsite)
    g2e_bb = ao2mo.restore(1, eri[2], nsite)
    t2a = numpy.dot(g2e_aa.reshape(nsite**2,-1), t1a.reshape(nsite**2,-1))
    t2a+= numpy.dot(g2e_ab.reshape(nsite**2,-1), t1b.reshape(nsite**2,-1))
    t2b = numpy.dot(g2e_ab.reshape(nsite**2,-1).T, t1a.reshape(nsite**2,-1))
    t2b+= numpy.dot(g2e_bb.reshape(nsite**2,-1), t1b.reshape(nsite**2,-1))

    t2a = t2a.reshape((nsite,nsite)+cishape)
    t2b = t2b.reshape((nsite,nsite)+cishape)
    fcinew = numpy.zeros(cishape)
    for str0, tab in enumerate(link_indexa):
        for a, i, str1, sign in tab:
            fcinew[str1] += sign * t2a[a,i,str0]
    for str0, tab in enumerate(link_indexb):
        for a, i, str1, sign in tab:
            fcinew[:,str1] += sign * t2b[a,i,:,str0]
    return fcinew.reshape(fcivec.shape)

def writeIntegralFile(DMRGCI, h1eff, eri_cas, ncas, nelec, ecore=0):
    if isinstance(nelec, (int, numpy.integer)):
        neleca = nelec//2 + nelec%2
        nelecb = nelec - neleca
    else :
        neleca, nelecb = nelec

    # The name of the FCIDUMP file, default is "FCIDUMP".
    integralFile = os.path.join(DMRGCI.runtimeDir, DMRGCI.integralFile)
    if DMRGCI.groupname is not None and DMRGCI.orbsym is not []:
# First removing the symmetry forbidden integrals. This has been done using
# the pyscf internal irrep-IDs (stored in DMRGCI.orbsym)
        orbsym = numpy.asarray(DMRGCI.orbsym) % 10
        pair_irrep = (orbsym.reshape(-1,1) ^ orbsym)[numpy.tril_indices(ncas)]
        sym_forbid = pair_irrep.reshape(-1,1) != pair_irrep.ravel()
        eri_cas = ao2mo.restore(4, eri_cas, ncas)
        eri_cas[sym_forbid] = 0
        eri_cas = ao2mo.restore(8, eri_cas, ncas)
       #orbsym = numpy.asarray(dmrg_sym.convert_orbsym(DMRGCI.groupname, DMRGCI.orbsym))
       #eri_cas = pyscf.ao2mo.restore(8, eri_cas, ncas)
# Then convert the pyscf internal irrep-ID to molpro irrep-ID
        orbsym = numpy.asarray(dmrg_sym.convert_orbsym(DMRGCI.groupname, orbsym))
    else:
        orbsym = []
        eri_cas = ao2mo.restore(8, eri_cas, ncas)
    if not os.path.exists(DMRGCI.scratchDirectory):
        os.makedirs(DMRGCI.scratchDirectory)
    if not os.path.exists(DMRGCI.runtimeDir):
        os.makedirs(DMRGCI.runtimeDir)

    tools.fcidump.from_integrals(integralFile, h1eff, eri_cas, ncas,
                                 neleca+nelecb, ecore, ms=abs(neleca-nelecb),
                                 orbsym=orbsym)
    return integralFile

    def test_nr_transe1incore(self):
        eri_ao = _vhf.int2e_sph(mol._atm, mol._bas, mol._env)
        eriref = ao2mo.restore(1, eri_ao, nao)
        eriref = numpy.einsum("ijpl,pk->ijkl", eriref, mo)
        eriref = numpy.einsum("ijkp,pl->ijkl", eriref, mo)
        eriref = ao2mo.restore(4, eriref, nao)

        eri_ao = ao2mo.restore(8, eri_ao, nao)
        ftrans1 = f1pointer("AO2MOtranse1_incore_s8")
        fmmm = f1pointer("AO2MOmmm_nr_s2_s2")
        eri1 = numpy.empty((naopair, naopair))
        libao2mo1.AO2MOnr_e1incore_drv(
            ftrans1,
            fmmm,
            eri1.ctypes.data_as(ctypes.c_void_p),
            eri_ao.ctypes.data_as(ctypes.c_void_p),
            mo.ctypes.data_as(ctypes.c_void_p),
            ctypes.c_int(0),
            ctypes.c_int(naopair),
            ctypes.c_int(nao),
            ctypes.c_int(0),
            ctypes.c_int(nao),
            ctypes.c_int(0),
            ctypes.c_int(nao),
        )
        self.assertTrue(numpy.allclose(eri1, eriref))

    def test_mp2_contract_eri_dm(self):
        nocc = mol.nelectron//2
        nmo = mf.mo_energy.size
        nvir = nmo - nocc

        pt = mp.mp2.MP2(mf)
        emp2, t2 = pt.kernel()
        eri = ao2mo.restore(1, ao2mo.kernel(mf._eri, mf.mo_coeff), nmo)
        hcore = mf.get_hcore()
        rdm1 = pt.make_rdm1()
        rdm2 = pt.make_rdm2()
        h1 = reduce(numpy.dot, (mf.mo_coeff.T, hcore, mf.mo_coeff))
        e1 = numpy.einsum('ij,ji', h1, rdm1)
        e1+= numpy.einsum('ijkl,ijkl', eri, rdm2) * .5
        e1+= mol.energy_nuc()
        self.assertAlmostEqual(e1, pt.e_tot, 9)

        pt.frozen = 2
        pt.max_memory = 1
        emp2, t2 = pt.kernel(with_t2=False)
        eri = ao2mo.restore(1, ao2mo.kernel(mf._eri, mf.mo_coeff), nmo)
        hcore = mf.get_hcore()
        rdm1 = pt.make_rdm1()
        rdm2 = pt.make_rdm2()
        h1 = reduce(numpy.dot, (mf.mo_coeff.T, hcore, mf.mo_coeff))
        e1 = numpy.einsum('ij,ji', h1, rdm1)
        e1+= numpy.einsum('ijkl,ijkl', eri, rdm2) * .5
        e1+= mol.energy_nuc()
        self.assertAlmostEqual(e1, pt.e_tot, 9)

    def __init__(self, myci, mo_coeff, method='incore'):
        mol = myci.mol
        mf = myci._scf
        nocc = myci.nocc
        nmo = myci.nmo
        nvir = nmo - nocc
        if mo_coeff is None:
            self.mo_coeff = mo_coeff = myci.mo_coeff
        if (method == 'incore' and mf._eri is not None):
            eri = ao2mo.kernel(mf._eri, mo_coeff, verbose=myci.verbose)
        else:
            eri = ao2mo.kernel(mol, mo_coeff, verbose=myci.verbose)
        eri = ao2mo.restore(1, eri, nmo)
        eri = eri.reshape(nmo,nmo,nmo,nmo)

        self.oooo = eri[:nocc,:nocc,:nocc,:nocc]
        self.vvoo = eri[nocc:,nocc:,:nocc,:nocc]
        self.vooo = eri[nocc:,:nocc,:nocc,:nocc]
        self.voov = eri[nocc:,:nocc,:nocc,nocc:]
        self.vovv = lib.pack_tril(eri[nocc:,:nocc,nocc:,nocc:].reshape(-1,nvir,nvir))
        self.vvvv = ao2mo.restore(4, eri[nocc:,nocc:,nocc:,nocc:].copy(), nvir)

        dm = mf.make_rdm1()
        vhf = mf.get_veff(mol, dm)
        h1 = mf.get_hcore(mol)
        self.fock = reduce(numpy.dot, (mo_coeff.T, h1 + vhf, mo_coeff))

def _gamma2_intermediates(cc, t1, t2, l1, l2):
    d2 = ccsd_rdm._gamma2_intermediates(cc, t1, t2, l1, l2)
    nocc, nvir = t1.shape
    if cc.frozen is None or cc.frozen is 0:
        dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov = d2
        dvvov = dovvv.transpose(2,3,0,1)
        dvvvv = ao2mo.restore(1, d2[1], nvir)
        return dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov
    nocc0 = numpy.count_nonzero(cc.mo_occ>0)
    nvir0 = cc.mo_occ.size - nocc0
    OA, VA, OF, VF = index_frozen_active(cc)
    dovov = numpy.zeros((nocc0,nvir0,nocc0,nvir0))
    dvvvv = numpy.zeros((nvir0,nvir0,nvir0,nvir0))
    doooo = numpy.zeros((nocc0,nocc0,nocc0,nocc0))
    doovv = numpy.zeros((nocc0,nocc0,nvir0,nvir0))
    dovvo = numpy.zeros((nocc0,nvir0,nvir0,nocc0))
    dovvv = numpy.zeros((nocc0,nvir0,nvir0,nvir0))
    dooov = numpy.zeros((nocc0,nocc0,nocc0,nvir0))
    dovov[OA[:,None,None,None],VA[:,None,None],OA[:,None],VA] = d2[0]
    dvvvv[VA[:,None,None,None],VA[:,None,None],VA[:,None],VA] = ao2mo.restore(1, d2[1], nvir)
    doooo[OA[:,None,None,None],OA[:,None,None],OA[:,None],OA] = d2[2]
    doovv[OA[:,None,None,None],OA[:,None,None],VA[:,None],VA] = d2[3]
    dovvo[OA[:,None,None,None],VA[:,None,None],VA[:,None],OA] = d2[4]
    dovvv[OA[:,None,None,None],VA[:,None,None],VA[:,None],VA] = d2[6]
    dooov[OA[:,None,None,None],OA[:,None,None],OA[:,None],VA] = d2[7]
    dvvov = dovvv.transpose(2,3,0,1)
    return dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov

def mp2(mol, h1e, eri, mo, nelec, with_1pdm, with_e2frag):
    eri1 = ao2mo.restore(8, eri, mo.shape[1])
    hf_energy, mo_energy, mo_coeff, mo_occ = simple_hf(h1e, eri1, mo, nelec)
    mf = scf.RHF(mol)
    mf.get_hcore = lambda mol: h1e
    mf._eri = eri
    mymp2 = mp.MP2(mf)
    mymp2.nocc = nelec // 2
    mymp2.nmo = len(mo_energy)
    emp2, t2 = mymp2.kernel(mo_energy, mo_coeff)

    if with_1pdm:
        rdm1 = mymp2.make_rdm1(t2)
        for i in range(nelec // 2):
            rdm1[i, i] += 2
        rdm1 = reduce(numpy.dot, (mo_coeff, rdm1, mo_coeff.T))
    else:
        rdm1 = None

    if with_e2frag:
        norb = mo_coeff.shape[0]
        eri1 = part_eri_hermi(eri, norb, with_e2frag)
        eri1 = ao2mo.incore.full(eri1, mo_coeff)
        eri1 = ao2mo.restore(1, eri1, norb)
        rdm2 = mymp2.make_rdm2(t2)
        for i in range(nelec // 2):
            for j in range(nelec // 2):
                rdm2[i, i, j, j] += 4
                rdm2[i, j, i, j] += -2
        e2frag = 0.5 * numpy.dot(rdm2.reshape(-1), eri1.reshape(-1))
    else:
        e2frag = None
    return hf_energy, emp2 + hf_energy, e2frag, rdm1

    def test_nroutcore_eri(self):
        ftmp = tempfile.NamedTemporaryFile(dir=lib.param.TMPDIR)
        erifile = ftmp.name
        eri_ao = ao2mo.restore(1, mol.intor('int2e', aosym='s2kl'), nao)
        eriref = numpy.einsum('pjkl,pi->ijkl', eri_ao, mo)
        eriref = numpy.einsum('ipkl,pj->ijkl', eriref, mo)
        eriref = numpy.einsum('ijpl,pk->ijkl', eriref, mo)
        eriref = numpy.einsum('ijkp,pl->ijkl', eriref, mo)

        mos = (mo[:,:4], mo[:,:3], mo[:,:3], mo[:,:2])
        ao2mo.outcore.general(mol, mos, erifile, dataname='eri_mo',
                              intor='int2e', aosym=1)
        with ao2mo.load(erifile) as eri1:
            eri1 = numpy.asarray(eri1).reshape(4,3,3,2)
        self.assertTrue(numpy.allclose(eri1, eriref[:4,:3,:3,:2]))

        ao2mo.outcore.full(mol, mo, erifile, dataname='eri_mo',
                              intor='int2e_sph', aosym='s2ij', comp=1)
        with ao2mo.load(erifile, 'eri_mo') as eri:
            eri1 = ao2mo.restore(1, numpy.array(eri), nao)
        eri1 = eri1.reshape(nao,nao,nao,nao)
        self.assertTrue(numpy.allclose(eri1, eriref))

        mos = (mo[:,:3], mo[:,:3], mo[:,:3], mo[:,:2])
        ao2mo.outcore.general(mol, mos, erifile, dataname='eri_mo',
                              intor='int2e_sph', aosym='s4', comp=1,
                              compact=False)
        with ao2mo.load(erifile, 'eri_mo') as eri1:
            eri1 = numpy.asarray(eri1).reshape(3,3,3,2)
        self.assertTrue(numpy.allclose(eri1, eriref[:3,:3,:3,:2]))

        ao2mo.outcore.full(mol, mo[:,:0], erifile,
                           intor='int2e', aosym='1', comp=1)
        with ao2mo.load(erifile, 'eri_mo') as eri:
            self.assertTrue(eri.size == 0)

def part_eri_hermi(eri, norb, nimp):
    eri1 = ao2mo.restore(4, eri, norb)
    for i in range(eri1.shape[0]):
        tmp = lib.unpack_tril(eri1[i])
        tmp[nimp:] = 0
        eri1[i] = lib.pack_tril(tmp + tmp.T)
    eri1 = lib.transpose_sum(eri1, inplace=True)
    return ao2mo.restore(8, eri1, norb) * 0.25

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 test_ccsd(self):
        mol = gto.M()
        mf = scf.RHF(mol)
        mcc = cc.CCSD(mf)
        numpy.random.seed(12)
        mcc.nocc = nocc = 5
        mcc.nmo = nmo = 12
        nvir = nmo - nocc
        eri0 = numpy.random.random((nmo,nmo,nmo,nmo))
        eri0 = ao2mo.restore(1, ao2mo.restore(8, eri0, nmo), nmo)
        fock0 = numpy.random.random((nmo,nmo))
        fock0 = fock0 + fock0.T + numpy.diag(range(nmo))*2
        t1 = numpy.random.random((nocc,nvir))
        t2 = numpy.random.random((nocc,nocc,nvir,nvir))
        t2 = t2 + t2.transpose(1,0,3,2)
        l1 = numpy.random.random((nocc,nvir))
        l2 = numpy.random.random((nocc,nocc,nvir,nvir))
        l2 = l2 + l2.transpose(1,0,3,2)

        eris = cc.ccsd._ChemistsERIs()
        eris.oooo = eri0[:nocc,:nocc,:nocc,:nocc].copy()
        eris.ovoo = eri0[:nocc,nocc:,:nocc,:nocc].copy()
        eris.oovv = eri0[:nocc,:nocc,nocc:,nocc:].copy()
        eris.ovvo = eri0[:nocc,nocc:,nocc:,:nocc].copy()
        idx = numpy.tril_indices(nvir)
        eris.ovvv = eri0[:nocc,nocc:,nocc:,nocc:][:,:,idx[0],idx[1]].copy()
        eris.vvvv = ao2mo.restore(4,eri0[nocc:,nocc:,nocc:,nocc:],nvir)
        eris.fock = fock0
        eris.mo_energy = fock0.diagonal()

        saved = ccsd_lambda.make_intermediates(mcc, t1, t2, eris)
        l1new, l2new = ccsd_lambda.update_lambda(mcc, t1, t2, l1, l2, eris, saved)
        self.assertAlmostEqual(abs(l1new).sum(), 38172.7896467303, 8)
        self.assertAlmostEqual(numpy.dot(l1new.flatten(), numpy.arange(35)), 739312.005491083, 8)
        self.assertAlmostEqual(numpy.dot(l1new.flatten(), numpy.sin(numpy.arange(35))), 7019.50937051188, 8)
        self.assertAlmostEqual(numpy.dot(numpy.sin(l1new.flatten()), numpy.arange(35)), 69.6652346635955, 8)

        self.assertAlmostEqual(abs(l2new).sum(), 72035.4931071527, 8)
        self.assertAlmostEqual(abs(l2new-l2new.transpose(1,0,3,2)).sum(), 0, 9)
        self.assertAlmostEqual(numpy.dot(l2new.flatten(), numpy.arange(35**2)), 48427109.5409886, 7)
        self.assertAlmostEqual(numpy.dot(l2new.flatten(), numpy.sin(numpy.arange(35**2))), 137.758016736487, 8)
        self.assertAlmostEqual(numpy.dot(numpy.sin(l2new.flatten()), numpy.arange(35**2)), 507.656936701192, 8)

        mcc.max_memory = 0
        saved = ccsd_lambda.make_intermediates(mcc, t1, t2, eris)
        l1new, l2new = ccsd_lambda.update_lambda(mcc, t1, t2, l1, l2, eris, saved)
        self.assertAlmostEqual(abs(l1new).sum(), 38172.7896467303, 8)
        self.assertAlmostEqual(numpy.dot(l1new.flatten(), numpy.arange(35)), 739312.005491083, 8)
        self.assertAlmostEqual(numpy.dot(l1new.flatten(), numpy.sin(numpy.arange(35))), 7019.50937051188, 8)
        self.assertAlmostEqual(numpy.dot(numpy.sin(l1new.flatten()), numpy.arange(35)), 69.6652346635955, 8)

        self.assertAlmostEqual(abs(l2new).sum(), 72035.4931071527, 8)
        self.assertAlmostEqual(abs(l2new-l2new.transpose(1,0,3,2)).sum(), 0, 9)
        self.assertAlmostEqual(numpy.dot(l2new.flatten(), numpy.arange(35**2)), 48427109.5409886, 7)
        self.assertAlmostEqual(numpy.dot(l2new.flatten(), numpy.sin(numpy.arange(35**2))), 137.758016736487, 8)
        self.assertAlmostEqual(numpy.dot(numpy.sin(l2new.flatten()), numpy.arange(35**2)), 507.656936701192, 8)

def gen_hop_uhf_external(mf):
    mol = mf.mol
    mo_a, mo_b = mf.mo_coeff
    mo_ea, mo_eb = mf.mo_energy
    mo_occa, mo_occb = mf.mo_occ
    nmo = mo_a.shape[1]
    nocca = numpy.count_nonzero(mo_occa)
    noccb = numpy.count_nonzero(mo_occb)
    nvira = nmo - nocca
    nvirb = nmo - noccb

    eri_aa = ao2mo.restore(1, ao2mo.full(mol, mo_a), nmo)
    eri_ab = ao2mo.restore(1, ao2mo.general(mol, [mo_a,mo_a,mo_b,mo_b]), nmo)
    eri_bb = ao2mo.restore(1, ao2mo.full(mol, mo_b), nmo)
    # alpha -> alpha
    haa =-numpy.einsum('abji->iajb', eri_aa[nocca:,nocca:,:nocca,:nocca])
    haa+= numpy.einsum('ajbi->iajb', eri_aa[nocca:,:nocca,nocca:,:nocca])
    for a in range(nvira):
        for i in range(nocca):
            haa[i,a,i,a] += mo_ea[nocca+a] - mo_ea[i]
    # beta -> beta
    hbb =-numpy.einsum('abji->iajb', eri_bb[noccb:,noccb:,:noccb,:noccb])
    hbb+= numpy.einsum('ajbi->iajb', eri_bb[noccb:,:noccb,noccb:,:noccb])
    for a in range(nvirb):
        for i in range(noccb):
            hbb[i,a,i,a] += mo_eb[noccb+a] - mo_eb[i]

    nova = nocca * nvira
    novb = noccb * nvirb
    h1 = numpy.zeros((nova+novb,nova+novb))
    h1[:nova,:nova] = haa.transpose(1,0,3,2).reshape(nova,nova)
    h1[nova:,nova:] = hbb.transpose(1,0,3,2).reshape(novb,novb)
    def hop1(x):
        return h1.dot(x)

    h11 =-numpy.einsum('abji->iajb', eri_ab[nocca:,nocca:,:noccb,:noccb])
    for a in range(nvira):
        for i in range(noccb):
            h11[i,a,i,a] += mo_ea[nocca+a] - mo_eb[i]
    h22 =-numpy.einsum('jiab->iajb', eri_ab[:nocca,:nocca,noccb:,noccb:])
    for a in range(nvirb):
        for i in range(nocca):
            h22[i,a,i,a] += mo_eb[noccb+a] - mo_ea[i]
    h12 =-numpy.einsum('ajbi->iajb', eri_ab[nocca:,:nocca,noccb:,:noccb])
    h21 =-numpy.einsum('biaj->iajb', eri_ab[nocca:,:nocca,noccb:,:noccb])

    n1 = noccb * nvira
    n2 = nocca * nvirb
    h2 = numpy.empty((n1+n2,n1+n2))
    h2[:n1,:n1] = h11.transpose(1,0,3,2).reshape(n1,n1)
    h2[n1:,n1:] = h22.transpose(1,0,3,2).reshape(n2,n2)
    h2[:n1,n1:] = h12.transpose(1,0,3,2).reshape(n1,n2)
    h2[n1:,:n1] = h21.transpose(1,0,3,2).reshape(n2,n1)
    def hop2(x):
        return h2.dot(x)
    return hop1, hop2

    def test_kernel(self):
        h1 = h1e + h1e.T
        eri = .5* ao2mo.restore(1, ao2mo.restore(8, h2e, norb), norb)
        h = fci.direct_spin1.pspace(h1, eri, norb, nelec, np=5000)[1]
        eref, c0 = numpy.linalg.eigh(h)

        sol = fci.direct_nosym.FCI()
        e, c1 = sol.kernel(h1, eri, norb, nelec, max_space=40)
        self.assertAlmostEqual(eref[0], e, 9)
        self.assertAlmostEqual(abs(c0[:,0].dot(c1.ravel())), 1, 9)

def _get_VVVV(eris):
    if eris.VVVV is None and getattr(eris, 'VVL', None) is not None:  # DF eris
        VVL = np.asarray(eris.VVL)
        nvir = int(np.sqrt(eris.VVL.shape[0]*2))
        return ao2mo.restore(1, lib.dot(VVL, VVL.T), nvir)
    elif len(eris.VVVV.shape) == 2:
        nvir = int(np.sqrt(eris.VVVV.shape[0]*2))
        return ao2mo.restore(1, np.asarray(eris.VVVV), nvir)
    else:
        return eris.VVVV

def absorb_h1e(h1e, eri, norb, nelec, fac=1):
    '''Modify 2e Hamiltonian to include 1e Hamiltonian contribution.
    '''
    if not isinstance(nelec, (int, numpy.number)):
        nelec = sum(nelec)
    h2e = ao2mo.restore(1, eri.copy(), norb)
    f1e = h1e - numpy.einsum('jiik->jk', h2e) * .5
    f1e = f1e * (1./(nelec+1e-100))
    for k in range(norb):
        h2e[k,k,:,:] += f1e
        h2e[:,:,k,k] += f1e
    return ao2mo.restore(4, h2e, norb) * fac

def contract_2e(eri, civec_strs, norb, nelec, link_index=None, orbsym=None):
    ci_coeff, nelec, ci_strs = selected_ci._unpack(civec_strs, nelec)
    if link_index is None:
        link_index = selected_ci._all_linkstr_index(ci_strs, norb, nelec)
    cd_indexa, dd_indexa, cd_indexb, dd_indexb = link_index
    na, nlinka = nb, nlinkb = cd_indexa.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
    lib.transpose_sum(eri1, inplace=True)
    eri1 *= .5
    eri1, dd_indexa, dimirrep = selected_ci_symm.reorder4irrep(eri1, norb, dd_indexa, orbsym, -1)
    fcivec = ci_coeff.reshape(na,nb)
    ci1 = numpy.zeros_like(fcivec)
    # (aa|aa)
    if nelec[0] > 1:
        ma, mlinka = mb, mlinkb = dd_indexa.shape[:2]
        libfci.SCIcontract_2e_aaaa_symm(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),
                                        dimirrep.ctypes.data_as(ctypes.c_void_p),
                                        ctypes.c_int(len(dimirrep)))

    h_ps = numpy.einsum('pqqs->ps', eri) * (.5/nelec[0])
    eri1 = eri.copy()
    for k in range(norb):
        eri1[:,:,k,k] += h_ps
        eri1[k,k,:,:] += h_ps
    eri1 = ao2mo.restore(4, eri1, norb)
    lib.transpose_sum(eri1, inplace=True)
    eri1 *= .5
    eri1, cd_indexa, dimirrep = selected_ci_symm.reorder4irrep(eri1, norb, cd_indexa, orbsym)
    # (bb|aa)
    libfci.SCIcontract_2e_bbaa_symm(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_indexa.ctypes.data_as(ctypes.c_void_p),
                                    dimirrep.ctypes.data_as(ctypes.c_void_p),
                                    ctypes.c_int(len(dimirrep)))

    lib.transpose_sum(ci1, inplace=True)
    return selected_ci._as_SCIvector(ci1.reshape(ci_coeff.shape), ci_strs)

    def test_large_ci(self):
        norb = 6
        nelec = (3,3)
        numpy.random.seed(10)
        h1e = numpy.random.random((norb,norb))
        h1e = h1e + h1e.T
        g2e = numpy.random.random((norb,norb,norb,norb))
        eri = .5* ao2mo.restore(1, ao2mo.restore(8, g2e, norb), norb)

        sol = fci.direct_spin1.FCI(mol)
        e, c1 = sol.kernel(h1e, eri, norb, nelec)
        self.assertAlmostEqual(e, -4.742664117996546, 9)
        val = sol.large_ci(c1, norb, nelec)
        self.assertAlmostEqual(val[0][0], c1[0,0], 9)

    def test_restore4(self):
        self.assertTrue(numpy.allclose(a4, ao2mo.restore(4, a1, n)))
        self.assertTrue(numpy.allclose(a4, ao2mo.restore('4', a4, n)))
        self.assertTrue(numpy.allclose(a4, ao2mo.restore('s4', a8, n)))
        self.assertTrue(numpy.allclose(a4, ao2mo.restore('s4', a2ij, n)))
        self.assertTrue(numpy.allclose(a4, ao2mo.restore('s4', a2kl, n)))

        self.assertTrue(numpy.allclose(b4, ao2mo.restore(4, b1, n)))
        self.assertTrue(numpy.allclose(b4, ao2mo.restore('4', b4, n)))
        self.assertTrue(numpy.allclose(b4, ao2mo.restore('4', b2ij, n)))
        self.assertTrue(numpy.allclose(b4, ao2mo.restore('4', b2kl, n)))

    def test_uccsd_rdm2_mo2ao(self):
        mol = gto.Mole()
        mol.verbose = 0
        mol.atom = [
            [8 , (0. , 0.     , 0.)],
            [1 , (0. , -0.757 , 0.587)],
            [1 , (0. , 0.757  , 0.587)]]
        mol.spin = 2
        mol.basis = '631g'
        mol.build(0, 0)
        mf = scf.UHF(mol)
        mf.conv_tol_grad = 1e-8
        mf.kernel()
        mycc = cc.UCCSD(mf)
        mycc.diis_start_cycle = 1
        mycc.conv_tol = 1e-10
        eris = mycc.ao2mo()
        ecc, t1, t2 = mycc.kernel(eris=eris)
        l1, l2 = mycc.solve_lambda(eris=eris)
        fdm2 = lib.H5TmpFile()
        d2 = cc.uccsd_rdm._gamma2_outcore(mycc, t1, t2, l1, l2, fdm2, True)

        nao = mycc.mo_coeff[0].shape[0]
        ref = cc.uccsd_rdm._make_rdm2(mycc, None, d2, with_dm1=False)
        aa = lib.einsum('ijkl,pi,qj,rk,sl->pqrs', ref[0], mycc.mo_coeff[0],
                        mycc.mo_coeff[0], mycc.mo_coeff[0], mycc.mo_coeff[0])
        ab = lib.einsum('ijkl,pi,qj,rk,sl->pqrs', ref[1], mycc.mo_coeff[0],
                        mycc.mo_coeff[0], mycc.mo_coeff[1], mycc.mo_coeff[1])
        bb = lib.einsum('ijkl,pi,qj,rk,sl->pqrs', ref[2], mycc.mo_coeff[1],
                        mycc.mo_coeff[1], mycc.mo_coeff[1], mycc.mo_coeff[1])
        aa = aa + aa.transpose(0,1,3,2)
        aa = aa + aa.transpose(1,0,2,3)
        aa = ao2mo.restore(4, aa, nao) * .5
        bb = bb + bb.transpose(0,1,3,2)
        bb = bb + bb.transpose(1,0,2,3)
        bb = ao2mo.restore(4, bb, nao) * .5
        ab = ab + ab.transpose(0,1,3,2)
        ab = ab + ab.transpose(1,0,2,3)
        ab = ao2mo.restore(4, ab, nao) * .5
        ref = (aa+ab, bb+ab.T)
        rdm2 = uccsd_grad._rdm2_mo2ao(mycc, d2, mycc.mo_coeff)
        self.assertAlmostEqual(abs(ref[0]-rdm2[0]).max(), 0, 10)
        self.assertAlmostEqual(abs(ref[1]-rdm2[1]).max(), 0, 10)
        uccsd_grad._rdm2_mo2ao(mycc, d2, mycc.mo_coeff, fdm2)
        self.assertAlmostEqual(abs(ref[0]-fdm2['dm2aa+ab'].value).max(), 0, 10)
        self.assertAlmostEqual(abs(ref[1]-fdm2['dm2bb+ab'].value).max(), 0, 10)
        self.assertAlmostEqual(lib.finger(rdm2[0]), -1.6247203743431637, 7)
        self.assertAlmostEqual(lib.finger(rdm2[1]), -0.44062825991527471, 7)