def test_h4_rdm(self):
     mol = gto.Mole()
     mol.verbose = 0
     mol.atom = [
         ['H', ( 1.,-1.    , 0.   )],
         ['H', ( 0.,-1.    ,-1.   )],
         ['H', ( 1.,-0.5   , 0.   )],
         ['H', ( 0.,-1.    , 1.   )],
     ]
     mol.charge = 2
     mol.spin = 2
     mol.basis = '6-31g'
     mol.build()
     mf = scf.UHF(mol).set(init_guess='1e').run(conv_tol=1e-14)
     ehf0 = mf.e_tot - mol.energy_nuc()
     mycc = uccsd.UCCSD(mf).run()
     mycc.solve_lambda()
     eri_aa = ao2mo.kernel(mf._eri, mf.mo_coeff[0])
     eri_bb = ao2mo.kernel(mf._eri, mf.mo_coeff[1])
     eri_ab = ao2mo.kernel(mf._eri, [mf.mo_coeff[0], mf.mo_coeff[0],
                                     mf.mo_coeff[1], mf.mo_coeff[1]])
     h1a = reduce(numpy.dot, (mf.mo_coeff[0].T, mf.get_hcore(), mf.mo_coeff[0]))
     h1b = reduce(numpy.dot, (mf.mo_coeff[1].T, mf.get_hcore(), mf.mo_coeff[1]))
     efci, fcivec = direct_uhf.kernel((h1a,h1b), (eri_aa,eri_ab,eri_bb),
                                      h1a.shape[0], mol.nelec)
     dm1ref, dm2ref = direct_uhf.make_rdm12s(fcivec, h1a.shape[0], mol.nelec)
     t1, t2 = mycc.t1, mycc.t2
     l1, l2 = mycc.l1, mycc.l2
     rdm1 = mycc.make_rdm1(t1, t2, l1, l2)
     rdm2 = mycc.make_rdm2(t1, t2, l1, l2)
     self.assertAlmostEqual(abs(dm1ref[0] - rdm1[0]).max(), 0, 6)
     self.assertAlmostEqual(abs(dm1ref[1] - rdm1[1]).max(), 0, 6)
     self.assertAlmostEqual(abs(dm2ref[0] - rdm2[0]).max(), 0, 6)
     self.assertAlmostEqual(abs(dm2ref[1] - rdm2[1]).max(), 0, 6)
     self.assertAlmostEqual(abs(dm2ref[2] - rdm2[2]).max(), 0, 6)

    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 _make_eris_incore(mp, mo_coeff=None, ao2mofn=None, verbose=None):
    eris = _PhysicistsERIs(mp, mo_coeff)

    nocc = mp.nocc
    nao, nmo = eris.mo_coeff.shape
    nvir = nmo - nocc
    orbspin = eris.orbspin

    if callable(ao2mofn):
        orbo = eris.mo_coeff[:,:nocc]
        orbv = eris.mo_coeff[:,nocc:]
        orbo = lib.tag_array(orbo, orbspin=orbspin)
        eri = ao2mofn((orbo,orbv,orbo,orbv)).reshape(nocc,nvir,nocc,nvir)
    else:
        orboa = eris.mo_coeff[:nao//2,:nocc]
        orbob = eris.mo_coeff[nao//2:,:nocc]
        orbva = eris.mo_coeff[:nao//2,nocc:]
        orbvb = eris.mo_coeff[nao//2:,nocc:]
        if orbspin is None:
            eri  = ao2mo.kernel(mp._scf._eri, (orboa,orbva,orboa,orbva))
            eri += ao2mo.kernel(mp._scf._eri, (orbob,orbvb,orbob,orbvb))
            eri1 = ao2mo.kernel(mp._scf._eri, (orboa,orbva,orbob,orbvb))
            eri += eri1
            eri += eri1.T
            eri = eri.reshape(nocc,nvir,nocc,nvir)
        else:
            co = orboa + orbob
            cv = orbva + orbvb
            eri = ao2mo.kernel(mp._scf._eri, (co,cv,co,cv)).reshape(nocc,nvir,nocc,nvir)
            sym_forbid = (orbspin[:nocc,None] != orbspin[nocc:])
            eri[sym_forbid,:,:] = 0
            eri[:,:,sym_forbid] = 0

    eris.oovv = eri.transpose(0,2,1,3) - eri.transpose(0,2,3,1)
    return eris

    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 test_ao2mo_with_mol_cart(self):
        pmol = mol.copy()
        pmol.cart = True
        nao = pmol.nao_nr()
        numpy.random.seed(1)
        mo = numpy.random.random((nao,4))
        eri = ao2mo.kernel(pmol, mo)
        self.assertAlmostEqual(lib.finger(eri), -977.99841341828437, 9)

        eri = ao2mo.kernel(mol, mo, intor='int2e_cart')
        self.assertAlmostEqual(lib.finger(eri), -977.99841341828437, 9)

    def test_r_outcore_eri(self):
        n2c = mol.nao_2c()
        numpy.random.seed(1)
        mo = numpy.random.random((n2c,n2c)) + numpy.random.random((n2c,n2c))*1j
        eriref = trans(eri0, [mo]*4)
        ftmp = tempfile.NamedTemporaryFile()

        ao2mo.kernel(mol, mo, ftmp.name, intor='int2e_spinor', max_memory=10, ioblk_size=5)
        with ao2mo.load(ftmp) as eri1:
            self.assertAlmostEqual(lib.finger(eri1), -550.72966498073129-1149.3561026721848j, 8)
            self.assertAlmostEqual(abs(eri1-eriref.reshape(n2c**2,n2c**2)).max(), 0, 9)

        eri1 = ao2mo.kernel(mol, (mo[:,:2], mo[:,:4], mo[:,:2], mo[:,:4]),
                            ftmp.name, intor='int2e_spinor')
        with ao2mo.load(ftmp) as eri1:
            self.assertAlmostEqual(abs(eri1-eriref[:2,:4,:2,:4].reshape(8,8)).max(), 0, 9)

        ftmp = lib.H5TmpFile()
        ao2mo.kernel(mol, (mo[:,:2], mo[:,:4], mo[:,:2], mo[:,:4]),
                     ftmp, intor='int2e_spinor', aosym='s1')
        with ao2mo.load(ftmp) as eri1:
            self.assertAlmostEqual(abs(eri1-eriref[:2,:4,:2,:4].reshape(8,8)).max(), 0, 9)

        eri1 = ao2mo.kernel(mol, (mo[:,:2], mo[:,:4], mo[:,:4], mo[:,:2]),
                            intor='int2e_spinor', aosym='s2ij')
        self.assertAlmostEqual(abs(eri1-eriref[:2,:4,:4,:2].reshape(8,8)).max(), 0, 9)

        eri1 = ao2mo.kernel(mol, (mo[:,:2], mo[:,:4], mo[:,:2], mo[:,:4]),
                            intor='int2e_spinor', aosym='s2kl')
        self.assertAlmostEqual(abs(eri1-eriref[:2,:4,:2,:4].reshape(8,8)).max(), 0, 9)

        eri1 = ao2mo.kernel(mol, mo[:,:0], intor='int2e_spinor')
        self.assertTrue(eri1.size == 0)

    def test_rdm_h4(self):
        mol = gto.Mole()
        mol.verbose = 7
        mol.output = '/dev/null'
        mol.atom = [
            ['O', ( 0., 0.    , 0.   )],
            ['H', ( 0., -0.757, 0.587)],
            ['H', ( 0., 0.757 , 0.587)],]
        mol.spin = 2
        mol.basis = 'sto-3g'
        mol.build()
        mf = scf.UHF(mol).run(conv_tol=1e-14)
        myci = ucisd.UCISD(mf)
        ecisd, civec = myci.kernel()
        self.assertAlmostEqual(ecisd, -0.033689623198003449, 8)

        nmoa = nmob = nmo = mf.mo_coeff[1].shape[1]
        nocc = (6,4)
        ci0 = myci.to_fcivec(civec, nmo, nocc)
        ref1, ref2 = fci.direct_uhf.make_rdm12s(ci0, nmo, nocc)
        rdm1 = myci.make_rdm1(civec)
        rdm2 = myci.make_rdm2(civec)
        self.assertAlmostEqual(abs(rdm1[0]-ref1[0]).max(), 0, 6)
        self.assertAlmostEqual(abs(rdm1[1]-ref1[1]).max(), 0, 6)
        self.assertAlmostEqual(abs(rdm2[0]-ref2[0]).max(), 0, 6)
        self.assertAlmostEqual(abs(rdm2[1]-ref2[1]).max(), 0, 6)
        self.assertAlmostEqual(abs(rdm2[2]-ref2[2]).max(), 0, 6)

        dm1a = numpy.einsum('ijkk->ji', rdm2[0]) / (mol.nelectron-1)
        dm1a+= numpy.einsum('ijkk->ji', rdm2[1]) / (mol.nelectron-1)
        self.assertAlmostEqual(abs(rdm1[0] - dm1a).max(), 0, 9)
        dm1b = numpy.einsum('kkij->ji', rdm2[2]) / (mol.nelectron-1)
        dm1b+= numpy.einsum('kkij->ji', rdm2[1]) / (mol.nelectron-1)
        self.assertAlmostEqual(abs(rdm1[1] - dm1b).max(), 0, 9)

        eri_aa = ao2mo.kernel(mf._eri, mf.mo_coeff[0], compact=False).reshape([nmoa]*4)
        eri_bb = ao2mo.kernel(mf._eri, mf.mo_coeff[1], compact=False).reshape([nmob]*4)
        eri_ab = ao2mo.kernel(mf._eri, [mf.mo_coeff[0], mf.mo_coeff[0],
                                        mf.mo_coeff[1], mf.mo_coeff[1]], compact=False)
        eri_ab = eri_ab.reshape(nmoa,nmoa,nmob,nmob)
        h1a = reduce(numpy.dot, (mf.mo_coeff[0].T, mf.get_hcore(), mf.mo_coeff[0]))
        h1b = reduce(numpy.dot, (mf.mo_coeff[1].T, mf.get_hcore(), mf.mo_coeff[1]))
        e2 = (numpy.einsum('ij,ji', h1a, rdm1[0]) +
              numpy.einsum('ij,ji', h1b, rdm1[1]) +
              numpy.einsum('ijkl,ijkl', eri_aa, rdm2[0]) * .5 +
              numpy.einsum('ijkl,ijkl', eri_ab, rdm2[1])      +
              numpy.einsum('ijkl,ijkl', eri_bb, rdm2[2]) * .5)
        e2 += mol.energy_nuc()
        self.assertAlmostEqual(myci.e_tot, e2, 9)

    def test_rdm1(self):
        mol = gto.Mole()
        mol.verbose = 0
        mol.atom = [
            ['O', ( 0., 0.    , 0.   )],
            ['H', ( 0., -0.757, 0.587)],
            ['H', ( 0., 0.757 , 0.587)],]
        mol.basis = '631g'
        mol.build()
        mf = scf.RHF(mol).run(conv_tol=1e-14)
        myci = ci.CISD(mf)
        myci.frozen = 1
        myci.kernel()

        nmo = myci.nmo
        nocc = myci.nocc
        d1 = cisd._gamma1_intermediates(myci, myci.ci, nmo, nocc)
        myci.max_memory = 0
        d2 = cisd._gamma2_intermediates(myci, myci.ci, nmo, nocc, True)
        dm2 = cisd.ccsd_rdm._make_rdm2(myci, d1, d2, with_dm1=True, with_frozen=True)
        dm1 = numpy.einsum('ijkk->ij', dm2)/(mol.nelectron-1)
        h1 = reduce(numpy.dot, (mf.mo_coeff.T, mf.get_hcore(), mf.mo_coeff))
        eri = ao2mo.restore(1, ao2mo.kernel(mf._eri, mf.mo_coeff), nmo+1)
        e1 = numpy.einsum('ij,ji', h1, dm1)
        e1+= numpy.einsum('ijkl,ijkl', eri, dm2) * .5
        e1+= mol.energy_nuc()
        self.assertAlmostEqual(e1, myci.e_tot, 7)

def fci(mol, printroots=4, visualize=False, preservedict=True):
    """Where is my docstring?"""
    if not preservedict:
        __hamdict = dict()
    print("PYCI")
    print("method: FCI")
    print("Number of eigenvalues: ", printroots)
    print("")
    Na, Nb = mol.nelec
    E_nuc = mol.energy_nuc()
    nao = mol.nao_nr()
    myhf = scf.RHF(mol)
    E_hf = myhf.kernel()
    mo_coefficients = myhf.mo_coeff
    h1e = reduce(np.dot, (mo_coefficients.T, myhf.get_hcore(), mo_coefficients))
    print("transforming eris")
    eri = ao2mo.kernel(mol, mo_coefficients)
    print("generating all determinants")
    fulldetlist_sets = gen_dets_sets(nao, Na, Nb)
    ndets = len(fulldetlist_sets)
    print("constructing full Hamiltonian")
    full_hamiltonian = construct_hamiltonian(ndets, fulldetlist_sets, h1e, eri)
    eig_vals, eig_vecs = sp.sparse.linalg.eigsh(full_hamiltonian, k=2 * printroots)
    eig_vals_sorted = np.sort(eig_vals)[:printroots]
    print("")
    print("first {:} pyci eigvals".format(printroots))
    for i in eig_vals_sorted + E_nuc:
        print(i)
    # if visualize:
    #     print(len(fulldetlist_sets),len(eig_vecs[:,np.argsort(eig_vals)[0]]))
    #     newdet = [i for i in zip(list(fulldetlist_sets),eig_vecs[:,np.argsort(eig_vals)[0]])]
    #     visualize_sets(newdet,nao,Na,Nb,"FCI")
    print("Completed FCI!")

    return eig_vals_sorted + E_nuc

    def test_rdm_h4(self):
        mol = gto.Mole()
        mol.verbose = 7
        mol.output = '/dev/null'
        mol.atom = [
            ['O', ( 0., 0.    , 0.   )],
            ['H', ( 0., -0.757, 0.587)],
            ['H', ( 0., 0.757 , 0.587)],]
        mol.spin = 2
        mol.basis = 'sto-3g'
        mol.build()
        mf = scf.RHF(mol).run(conv_tol=1e-14)
        myci = ci.GCISD(mf)
        eris = myci.ao2mo()
        ecisd, civec = myci.kernel(eris=eris)
        self.assertAlmostEqual(ecisd, -0.035165114624046617, 8)

        nmo = eris.mo_coeff.shape[1]
        rdm1 = myci.make_rdm1(civec, nmo, mol.nelectron)
        rdm2 = myci.make_rdm2(civec, nmo, mol.nelectron)

        mo = eris.mo_coeff[:7] + eris.mo_coeff[7:]
        eri = ao2mo.kernel(mf._eri, mo, compact=False).reshape([nmo]*4)
        eri[eris.orbspin[:,None]!=eris.orbspin,:,:] = 0
        eri[:,:,eris.orbspin[:,None]!=eris.orbspin] = 0
        h1e = reduce(numpy.dot, (mo.T, mf.get_hcore(), mo))
        h1e[eris.orbspin[:,None]!=eris.orbspin] = 0
        e2 = (numpy.einsum('ij,ji', h1e, rdm1) +
              numpy.einsum('ijkl,ijkl', eri, rdm2) * .5)
        e2 += mol.energy_nuc()
        self.assertAlmostEqual(myci.e_tot, e2, 9)

        dm1 = numpy.einsum('ijkk->ji', rdm2)/(mol.nelectron-1)
        self.assertAlmostEqual(abs(rdm1 - dm1).max(), 0, 9)

    def test_general(self):
        nao = mol.nao
        mo = numpy.random.random((nao,4))

        eri = mol.intor('int2e', aosym='s8')
        eri = ao2mo.kernel(eri, [mo]*4, compact=False)
        self.assertEqual(eri.shape, (16,16))

        h5file = lib.H5TmpFile()
        ao2mo.kernel(mol, [mo]*4, h5file, intor='int2e', dataname='eri')
        self.assertEqual(h5file['eri'].shape, (10,10))

        ftmp = tempfile.NamedTemporaryFile()
        ao2mo.kernel(mol, [mo]*4, ftmp, intor='int2e', dataname='eri')
        with ao2mo.load(ftmp.name, 'eri') as eri:
            self.assertEqual(eri.shape, (10,10))

    def test_rdm(self):
        mol = gto.Mole()
        mol.verbose = 5
        mol.output = '/dev/null'
        mol.atom = [
            ['O', ( 0., 0.    , 0.   )],
            ['H', ( 0., -0.757, 0.587)],
            ['H', ( 0., 0.757 , 0.587)],]
        mol.basis = {'H': 'sto-3g',
                     'O': 'sto-3g',}
        mol.build()
        mf = scf.RHF(mol).run()
        myci = ci.CISD(mf)
        ecisd, civec = myci.kernel()
        self.assertAlmostEqual(ecisd, -0.048878084082066106, 12)

        nmo = mf.mo_coeff.shape[1]
        nocc = mol.nelectron//2
        rdm1 = myci.make_rdm1(civec)
        rdm2 = myci.make_rdm2(civec)
        self.assertAlmostEqual(finger(rdm1), 2.2685199485281689, 9)
        self.assertAlmostEqual(finger(rdm2), 6.9406226727865876, 9)
        self.assertAlmostEqual(abs(rdm2-rdm2.transpose(2,3,0,1)).sum(), 0, 9)
        self.assertAlmostEqual(abs(rdm2-rdm2.transpose(1,0,3,2)).sum(), 0, 9)
        h1e = reduce(numpy.dot, (mf.mo_coeff.T, mf.get_hcore(), mf.mo_coeff))
        h2e = ao2mo.kernel(mf._eri, mf.mo_coeff)
        h2e = ao2mo.restore(1, h2e, nmo)
        e2 = (numpy.dot(h1e.flatten(),rdm1.flatten()) +
              numpy.dot(h2e.transpose(0,2,1,3).flatten(),rdm2.flatten()) * .5)
        self.assertAlmostEqual(ecisd + mf.e_tot - mol.energy_nuc(), e2, 9)
        dm1 = numpy.einsum('ikjk->ij', rdm2)/(mol.nelectron-1)
        self.assertAlmostEqual(abs(rdm1 - dm1).sum(), 0, 9)

    def test_contract(self):
        mol = gto.M()
        mol.nelectron = 6
        nocc, nvir = mol.nelectron//2, 4
        nmo = nocc + nvir
        nmo_pair = nmo*(nmo+1)//2
        mf = scf.RHF(mol)
        numpy.random.seed(12)
        mf._eri = numpy.random.random(nmo_pair*(nmo_pair+1)//2)
        mf.mo_coeff = numpy.random.random((nmo,nmo))
        mf.mo_energy = numpy.arange(0., nmo)
        mf.mo_occ = numpy.zeros(nmo)
        mf.mo_occ[:nocc] = 2
        dm = mf.make_rdm1()
        vhf = mf.get_veff(mol, dm)
        h1 = numpy.random.random((nmo,nmo)) * .1
        h1 = h1 + h1.T
        mf.get_hcore = lambda *args: h1

        myci = ci.CISD(mf)
        eris = myci.ao2mo(mf.mo_coeff)
        eris.ehf = (h1*dm).sum() + (vhf*dm).sum()*.5

        c2 = numpy.random.random((nocc,nocc,nvir,nvir)) * .1
        c2 = c2 + c2.transpose(1,0,3,2)
        c1 = numpy.random.random(nocc*nvir+1) * .1
        c0, c1 = c1[0], c1[1:].reshape(nocc,nvir)
        civec = myci.amplitudes_to_cisdvec(c0, c1, c2)
        hcivec = ci.cisd.contract(myci, civec, eris)
        self.assertAlmostEqual(lib.finger(hcivec), 2059.5730673341673, 9)
        e2 = ci.cisd.dot(civec, hcivec+eris.ehf*civec, nmo, nocc)
        self.assertAlmostEqual(e2, 7226.7494656749295, 9)

        rdm2 = myci.make_rdm2(civec, nmo, nocc)
        self.assertAlmostEqual(lib.finger(rdm2), 2.0492023431953221, 9)

        def fcicontract(h1, h2, norb, nelec, ci0):
            g2e = fci.direct_spin1.absorb_h1e(h1, h2, norb, nelec, .5)
            ci1 = fci.direct_spin1.contract_2e(g2e, ci0, norb, nelec)
            return ci1
        ci0 = myci.to_fcivec(civec, nmo, mol.nelec)
        self.assertAlmostEqual(abs(civec-myci.from_fcivec(ci0, nmo, nocc*2)).max(), 0, 9)
        h2e = ao2mo.kernel(mf._eri, mf.mo_coeff)
        h1e = reduce(numpy.dot, (mf.mo_coeff.T, h1, mf.mo_coeff))
        ci1 = fcicontract(h1e, h2e, nmo, mol.nelec, ci0)
        ci2 = myci.to_fcivec(hcivec, nmo, mol.nelec)
        e1 = numpy.dot(ci1.ravel(), ci0.ravel())
        e2 = ci.cisd.dot(civec, hcivec+eris.ehf*civec, nmo, nocc)
        self.assertAlmostEqual(e1, e2, 9)

        dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov = \
                ci.cisd._gamma2_intermediates(myci, civec, nmo, nocc)
        self.assertAlmostEqual(lib.finger(numpy.array(dovov)), 0.02868859991188923, 9)
        self.assertAlmostEqual(lib.finger(numpy.array(dvvvv)),-0.05524957311823144, 9)
        self.assertAlmostEqual(lib.finger(numpy.array(doooo)), 0.01014399192065793, 9)
        self.assertAlmostEqual(lib.finger(numpy.array(doovv)), 0.02761239887072825, 9)
        self.assertAlmostEqual(lib.finger(numpy.array(dovvo)), 0.09971200182238759, 9)
        self.assertAlmostEqual(lib.finger(numpy.array(dovvv)), 0.12777531252787638, 9)
        self.assertAlmostEqual(lib.finger(numpy.array(dooov)), 0.18667669732858014, 9)

def grad_elec(mc, mf_grad):
    mf = mf_grad.base
    mol = mf_grad.mol
    mo_energy = mc.mo_energy
    mo_coeff = mc.mo_coeff

    ncore = mc.ncore
    ncas = mc.ncas
    nocc = ncore + ncas
    nelecas = mc.nelecas
    nao, nmo = mo_coeff.shape

    hcore_deriv = mf_grad.hcore_generator(mol)
    s1 = mf_grad.get_ovlp(mol)

    casdm1, casdm2 = mc.fcisolver.make_rdm12(mc.ci, ncas, nelecas)
    dm1 = numpy.zeros((nmo,nmo))
    dm1[numpy.diag_indices(ncore)] = 2
    dm1[ncore:nocc,ncore:nocc] = casdm1

    dm2 = numpy.zeros((nmo,nmo,nmo,nmo))
    for i in range(ncore):
        for j in range(ncore):
            dm2[i,i,j,j] += 4
            dm2[i,j,j,i] -= 2
        dm2[i,i,ncore:nocc,ncore:nocc] = casdm1 * 2
        dm2[ncore:nocc,ncore:nocc,i,i] = casdm1 * 2
        dm2[i,ncore:nocc,ncore:nocc,i] =-casdm1
        dm2[ncore:nocc,i,i,ncore:nocc] =-casdm1
    dm2[ncore:nocc,ncore:nocc,ncore:nocc,ncore:nocc] = casdm2

    h1 = reduce(numpy.dot, (mo_coeff.T, mc._scf.get_hcore(), mo_coeff))
    h2 = ao2mo.kernel(mf._eri, mo_coeff, compact=False).reshape([nmo]*4)

# Generalized Fock, according to generalized Brillouin theorm
# Adv. Chem. Phys., 69, 63
    gfock = numpy.dot(h1, dm1)
    gfock+= numpy.einsum('iqrs,qjsr->ij', h2, dm2)
    gfock = (gfock + gfock.T) * .5
    dme0 = reduce(numpy.dot, (mo_coeff[:,:nocc], gfock[:nocc,:nocc], mo_coeff[:,:nocc].T))

    dm1 = reduce(numpy.dot, (mo_coeff, dm1, mo_coeff.T))
    dm2 = lib.einsum('ijkl,pi,qj,rk,sl->pqrs', dm2,
                     mo_coeff, mo_coeff, mo_coeff, mo_coeff)

    eri_deriv1 = mol.intor('int2e_ip1', comp=3).reshape(3,nao,nao,nao,nao)

    atmlst = range(mol.natm)
    aoslices = mol.aoslice_by_atom()
    de = numpy.zeros((len(atmlst),3))

    for k, ia in enumerate(atmlst):
        shl0, shl1, p0, p1 = aoslices[ia]
        h1ao = hcore_deriv(ia)
        de[k] += numpy.einsum('xij,ij->x', h1ao, dm1)
        de[k] -= numpy.einsum('xij,ij->x', s1[:,p0:p1], dme0[p0:p1]) * 2
        de[k] -= numpy.einsum('xijkl,ijkl->x', eri_deriv1[:,p0:p1], dm2[p0:p1]) * 2

    return de

def cisd(mol, printroots=4, visualize=False, preservedict=True):
    """Where is my docstring?"""
    if not preservedict:
        __hamdict = dict()
    print("PYCI")
    print("method: CISD")
    print("Number of eigenvalues: ", printroots)
    print("")
    Na, Nb = mol.nelec
    E_nuc = mol.energy_nuc()
    # get number of spatial atomic orbitals which we assume is the
    # same as the number of spatial molecular orbitals which is true
    # when using RHF, as long as there are no linear dependencies
    nao = mol.nao_nr()
    myhf = scf.RHF(mol)
    E_hf = myhf.kernel()
    mo_coefficients = myhf.mo_coeff
    # create matrix of 1 electron integrals
    h1e = reduce(np.dot, (mo_coefficients.T, myhf.get_hcore(), mo_coefficients))
    print("transforming eris")
    # get 2 electron integrals
    eri = ao2mo.kernel(mol, mo_coefficients)
    # create dictionary of hamiltonian elements
    hamdict = dict()
    # create Hartree-Fock determinant bit representation
    hfdet = (frozenset(range(Na)), frozenset(range(Nb)))
     # initalize set of target determinants
    targetdetset = set()
    # initalize set of core determinants
    coreset = {hfdet}
    print("\nHartree-Fock Energy: ", E_hf)
    # for each determinant in the core set of determinants...
    for idet in coreset:
        # ...generate all single and double excitations and add them
        # to the target set
        targetdetset |= set(gen_singles_doubles(idet, nao))
    # create a set from the union of the target and core determinants
    targetdetset |= coreset
    # update the dictionary of hamiltonian elements with those of the
    # target determinants
    hamdict.update(populatehamdict(targetdetset, targetdetset, hamdict, h1e, eri))
    # construct the hamiltonian in the target space
    targetham = getsmallham(list(targetdetset), hamdict)
    # diagonalize the hamiltonian yielding the eigenvalues and
    # eigenvectors
    eig_vals, eig_vecs = sp.sparse.linalg.eigsh(targetham, k=2 * printroots)
    eig_vals_sorted = np.sort(eig_vals)[:printroots] # sort the eigenvalues
    print("")
    print("first {:} pyci eigvals".format(printroots))
    for i in eig_vals_sorted + E_nuc:
        print(i)
    print(np.shape(eig_vecs))
    # if visualize:
    #     newdet = [i for i in zip(list(targetdetset),eig_vecs[:,np.argsort(eig_vals)[0]])]
    #     visualize_sets(newdet,nao,Na,Nb,"CISD")
    print("Completed CISD!")

    return eig_vals_sorted + E_nuc

def create_PYSCF_fcidump():
    myhf = scf.RHF(mol)
    E = myhf.kernel()
    c = myhf.mo_coeff
    h1e = reduce(np.dot, (c.T, myhf.get_hcore(), c))
    eri = ao2mo.kernel(mol, c)
    pt.fcidump.from_integrals('fcidump.txt', h1e, eri, c.shape[1],mol.nelectron, ms=0)
    cisolver = fci.FCI(mol, myhf.mo_coeff)
    print('E(HF) = %.12f, E(FCI) = %.12f' % (E,(cisolver.kernel()[0] + mol.energy_nuc())))

    def test_h4(self):
        mol = gto.Mole()
        mol.verbose = 7
        mol.output = '/dev/null'
        mol.atom = [
            ['H', ( 1.,-1.    , 0.   )],
            ['H', ( 0.,-1.    ,-1.   )],
            ['H', ( 1.,-0.5   , 0.   )],
            ['H', ( 0.,-1.    , 1.   )],
        ]
        mol.charge = 2
        mol.spin = 2
        mol.basis = '3-21g'
        mol.build()
        mf = scf.RHF(mol).run(conv_tol=1e-14)
        myci = ci.CISD(mf)
        myci.kernel()
        self.assertAlmostEqual(myci.e_tot, -0.50569591904536926, 8)

        mf = scf.UHF(mol).run(conv_tol=1e-14)
        myci = ci.CISD(mf)
        ecisd = myci.kernel()[0]
        self.assertAlmostEqual(ecisd, -0.03598992307519909, 8)
        self.assertAlmostEqual(myci.e_tot, -0.50569591904536926, 8)

        eris = myci.ao2mo()
        ecisd = myci.kernel(eris=eris)[0]
        eri_aa = ao2mo.kernel(mf._eri, mf.mo_coeff[0])
        eri_bb = ao2mo.kernel(mf._eri, mf.mo_coeff[1])
        eri_ab = ao2mo.kernel(mf._eri, [mf.mo_coeff[0], mf.mo_coeff[0],
                                        mf.mo_coeff[1], mf.mo_coeff[1]])
        h1a = reduce(numpy.dot, (mf.mo_coeff[0].T, mf.get_hcore(), mf.mo_coeff[0]))
        h1b = reduce(numpy.dot, (mf.mo_coeff[1].T, mf.get_hcore(), mf.mo_coeff[1]))
        efci, fcivec = fci.direct_uhf.kernel((h1a,h1b), (eri_aa,eri_ab,eri_bb),
                                             h1a.shape[0], mol.nelec)
        self.assertAlmostEqual(mf.e_tot+ecisd, efci+mol.energy_nuc(), 9)
        dm1ref, dm2ref = fci.direct_uhf.make_rdm12s(fcivec, h1a.shape[0], mol.nelec)
        rdm1 = myci.make_rdm1(myci.ci, myci.get_nmo(), myci.get_nocc())
        rdm2 = myci.make_rdm2(myci.ci, myci.get_nmo(), myci.get_nocc())
        self.assertAlmostEqual(abs(dm1ref[0] - rdm1[0]).max(), 0, 5)
        self.assertAlmostEqual(abs(dm1ref[1] - rdm1[1]).max(), 0, 5)
        self.assertAlmostEqual(abs(dm2ref[0] - rdm2[0]).max(), 0, 5)
        self.assertAlmostEqual(abs(dm2ref[1] - rdm2[1]).max(), 0, 5)
        self.assertAlmostEqual(abs(dm2ref[2] - rdm2[2]).max(), 0, 5)

def _make_eris_incore(mycc, mo_coeff=None, ao2mofn=None):
    cput0 = (time.clock(), time.time())
    eris = _PhysicistsERIs()
    eris._common_init_(mycc, mo_coeff)
    nocc = eris.nocc
    nao, nmo = eris.mo_coeff.shape

    if callable(ao2mofn):
        eri = ao2mofn(eris.mo_coeff).reshape([nmo]*4)
    else:
        assert(eris.mo_coeff.dtype == np.double)
        mo_a = eris.mo_coeff[:nao//2]
        mo_b = eris.mo_coeff[nao//2:]
        orbspin = eris.orbspin
        if orbspin is None:
            eri  = ao2mo.kernel(mycc._scf._eri, mo_a)
            eri += ao2mo.kernel(mycc._scf._eri, mo_b)
            eri1 = ao2mo.kernel(mycc._scf._eri, (mo_a,mo_a,mo_b,mo_b))
            eri += eri1
            eri += eri1.T
        else:
            mo = mo_a + mo_b
            eri = ao2mo.kernel(mycc._scf._eri, mo)
            if eri.size == nmo**4:  # if mycc._scf._eri is a complex array
                sym_forbid = (orbspin[:,None] != orbspin).ravel()
            else:  # 4-fold symmetry
                sym_forbid = (orbspin[:,None] != orbspin)[np.tril_indices(nmo)]
            eri[sym_forbid,:] = 0
            eri[:,sym_forbid] = 0

        if eri.dtype == np.double:
            eri = ao2mo.restore(1, eri, nmo)

    eri = eri.reshape(nmo,nmo,nmo,nmo)
    eri = eri.transpose(0,2,1,3) - eri.transpose(0,2,3,1)

    eris.oooo = eri[:nocc,:nocc,:nocc,:nocc].copy()
    eris.ooov = eri[:nocc,:nocc,:nocc,nocc:].copy()
    eris.oovv = eri[:nocc,:nocc,nocc:,nocc:].copy()
    eris.ovov = eri[:nocc,nocc:,:nocc,nocc:].copy()
    eris.ovvo = eri[:nocc,nocc:,nocc:,:nocc].copy()
    eris.ovvv = eri[:nocc,nocc:,nocc:,nocc:].copy()
    eris.vvvv = eri[nocc:,nocc:,nocc:,nocc:].copy()
    return eris

    def test_ump2_dm(self):
        pt = mp.MP2(mf)
        emp2, t2 = pt.kernel()
        dm1 = pt.make_rdm1()
        dm2 = pt.make_rdm2()
        gpt = mp.GMP2(mf).run()
        dm1ref = gpt.make_rdm1()
        dm2ref = gpt.make_rdm2()
        ia = gpt._scf.mo_coeff.orbspin == 0
        ib = gpt._scf.mo_coeff.orbspin == 1
        mo_a, mo_b = mf.mo_coeff
        nmoa = mo_a.shape[1]
        nmob = mo_b.shape[1]
        nocca, noccb = mol.nelec

        self.assertTrue(numpy.allclose(dm1[0], dm1ref[ia][:,ia]))
        self.assertTrue(numpy.allclose(dm1[1], dm1ref[ib][:,ib]))
        self.assertTrue(numpy.allclose(dm2[0], dm2ref[ia][:,ia][:,:,ia][:,:,:,ia]))
        self.assertTrue(numpy.allclose(dm2[2], dm2ref[ib][:,ib][:,:,ib][:,:,:,ib]))
        self.assertTrue(numpy.allclose(dm2[1], dm2ref[ia][:,ia][:,:,ib][:,:,:,ib]))

        hcore = mf.get_hcore()
        eriaa = ao2mo.kernel(mf._eri, mo_a, compact=False).reshape([nmoa]*4)
        eribb = ao2mo.kernel(mf._eri, mo_b, compact=False).reshape([nmob]*4)
        eriab = ao2mo.kernel(mf._eri, (mo_a,mo_a,mo_b,mo_b), compact=False)
        eriab = eriab.reshape([nmoa,nmoa,nmob,nmob])
        h1a = reduce(numpy.dot, (mo_a.T.conj(), hcore, mo_a))
        h1b = reduce(numpy.dot, (mo_b.T.conj(), hcore, mo_b))
        e1 = numpy.einsum('ij,ji', h1a, dm1[0])
        e1+= numpy.einsum('ij,ji', h1b, dm1[1])
        e1+= numpy.einsum('ijkl,ijkl', eriaa, dm2[0]) * .5
        e1+= numpy.einsum('ijkl,ijkl', eriab, dm2[1])
        e1+= numpy.einsum('ijkl,ijkl', eribb, dm2[2]) * .5
        e1+= mol.energy_nuc()
        self.assertAlmostEqual(e1, pt.e_tot, 9)

        vhf = mf.get_veff(mol, mf.make_rdm1())
        h1a = reduce(numpy.dot, (mo_a.T, hcore+vhf[0], mo_a))
        h1b = reduce(numpy.dot, (mo_b.T, hcore+vhf[1], mo_b))
        dm1[0][numpy.diag_indices(nocca)] -= 1
        dm1[1][numpy.diag_indices(noccb)] -= 1
        e = numpy.einsum('pq,qp', h1a, dm1[0])
        e+= numpy.einsum('pq,qp', h1b, dm1[1])
        self.assertAlmostEqual(e, -emp2, 9)

    def test_rdm_real(self):
        mol = gto.M()
        mol.verbose = 0
        nocc = 6
        nvir = 10
        mf = scf.GHF(mol)
        nmo = nocc + nvir
        npair = nmo*(nmo//2+1)//4
        numpy.random.seed(12)
        mf._eri = numpy.random.random(npair*(npair+1)//2)*.3
        hcore = numpy.random.random((nmo,nmo)) * .5
        hcore = hcore + hcore.T + numpy.diag(range(nmo))*2
        mf.get_hcore = lambda *args: hcore
        mf.get_ovlp = lambda *args: numpy.eye(nmo)
        mf.mo_coeff = numpy.eye(nmo)
        mf.mo_occ = numpy.zeros(nmo)
        mf.mo_occ[:nocc] = 1
        dm1 = mf.make_rdm1()
        mf.e_tot = mf.energy_elec()[0]
        myci = gcisd.GCISD(mf).run()
        dm1 = myci.make_rdm1()
        dm2 = myci.make_rdm2()

        nao = nmo // 2
        mo_a = mf.mo_coeff[:nao]
        mo_b = mf.mo_coeff[nao:]
        eri  = ao2mo.kernel(mf._eri, mo_a)
        eri += ao2mo.kernel(mf._eri, mo_b)
        eri1 = ao2mo.kernel(mf._eri, (mo_a,mo_a,mo_b,mo_b))
        eri += eri1
        eri += eri1.T
        eri = ao2mo.restore(1, eri, nmo)
        h1 = reduce(numpy.dot, (mf.mo_coeff.T.conj(), hcore, mf.mo_coeff))
        e1 = numpy.einsum('ij,ji', h1, dm1)
        e1+= numpy.einsum('ijkl,ijkl', eri, dm2) * .5
        self.assertAlmostEqual(e1, myci.e_tot, 7)

        self.assertAlmostEqual(abs(dm2-dm2.transpose(1,0,3,2).conj()).max(), 0, 9)
        self.assertAlmostEqual(abs(dm2-dm2.transpose(2,3,0,1)       ).max(), 0, 9)
        self.assertAlmostEqual(abs(dm2+dm2.transpose(2,1,0,3)       ).max(), 0, 9)
        self.assertAlmostEqual(abs(dm2+dm2.transpose(0,3,2,1)       ).max(), 0, 9)