ncrubin/n2_correct_active_space_vs_UHF

## n2_correct_active_space_vs_UHF
from pyscf import gto, scf, mcscf, ao2mo, fci
import pyscf
from pyscf.fci.cistring import make_strings
import numpy as np
import openfermion as of
from openfermion.chem.molecular_data import spinorb_from_spatial
import fqe
from fqe.openfermion_utils import integrals_to_fqe_restricted
from functools import reduce
import matplotlib.pyplot as plt


def get_molecule(bd):
    mol = gto.M()
    mol.atom = "N 0 0 0; N 0 0 {}".format(bd)
    mol.basis = 'sto3g'
    mol.spin = 0
    mol.build()
    return mol


def pyscf_to_fqe_wf(pyscf_cimat, pyscf_mf=None, norbs=None, nelec=None):
    if pyscf_mf is None:
        assert norbs is not None
        assert nelec is not None
    else:
        mol = pyscf_mf.mol
        nelec = mol.nelec
        norbs = pyscf_mf.mo_coeff.shape[1]

    norb_list = tuple(list(range(norbs)))
    n_alpha_strings = [x for x in make_strings(norb_list, nelec[0])]
    n_beta_strings = [x for x in make_strings(norb_list, nelec[1])]

    fqe_wf_ci = fqe.Wavefunction([[sum(nelec), nelec[0] - nelec[1], norbs]])
    fqe_data_ci = fqe_wf_ci.sector((sum(nelec), nelec[0] - nelec[1]))
    fqe_graph_ci = fqe_data_ci.get_fcigraph()
    fqe_orderd_coeff = np.zeros((fqe_graph_ci.lena(), fqe_graph_ci.lenb()), dtype=np.complex128)
    for paidx, pyscf_alpha_idx in enumerate(n_alpha_strings):
        for pbidx, pyscf_beta_idx in enumerate(n_beta_strings):
            # if np.abs(civec[paidx, pbidx]) > 1.0E-3:
            #     print(np.binary_repr(pyscf_alpha_idx, width=10), np.binary_repr(pyscf_beta_idx, width=10), civec[paidx, pbidx])
            fqe_orderd_coeff[fqe_graph_ci.index_alpha(
                pyscf_alpha_idx), fqe_graph_ci.index_beta(pyscf_beta_idx)] = \
                pyscf_cimat[paidx, pbidx]

    fqe_data_ci.coeff = fqe_orderd_coeff
    return fqe_wf_ci


def spinorb_from_spatial_blocks(h1a, h1b, eriaa, eribb, eriab,
                                EQ_TOLERANCE=1.0E-12):
    n_qubits = 2 * h1a.shape[0]

    # Initialize Hamiltonian coefficients.
    one_body_coefficients = np.zeros((n_qubits, n_qubits))
    two_body_coefficients = np.zeros(
        (n_qubits, n_qubits, n_qubits, n_qubits))
    # Loop through integrals.
    for p in range(n_qubits // 2):
        for q in range(n_qubits // 2):

            # Populate 1-body coefficients. Require p and q have same spin.
            one_body_coefficients[2 * p, 2 * q] = h1a[p, q]
            one_body_coefficients[2 * p + 1, 2 * q + 1] = h1b[p, q]
            # Continue looping to prepare 2-body coefficients.
            for r in range(n_qubits // 2):
                for s in range(n_qubits // 2):
                    # Mixed spin
                    two_body_coefficients[2 * p, 2 * q + 1, 2 * r + 1, 2 * s] = (eriab[p, q, r, s])
                    two_body_coefficients[2 * p + 1, 2 * q, 2 * r, 2 * s + 1] = (eriab[q, p, s, r])

                    # Same spin
                    two_body_coefficients[2 * p, 2 * q, 2 * r, 2 * s] = (eriaa[p, q, r, s])
                    two_body_coefficients[2 * p + 1, 2 * q + 1, 2 * r + 1, 2 * s + 1] = (eribb[p, q, r, s])

    # Truncate.
    one_body_coefficients[
        np.absolute(one_body_coefficients) < EQ_TOLERANCE] = 0.
    two_body_coefficients[
        np.absolute(two_body_coefficients) < EQ_TOLERANCE] = 0.

    return one_body_coefficients, two_body_coefficients


def active_space_uhf(mol, mf, run_stability=False, previous_rho=None):
    cas = mcscf.CASCI(mf, ncas=8, nelecas=sum(mol.nelec) - 4)
    h1, ecore = cas.get_h1eff()
    eri = cas.get_h2cas()
    h1 = np.ascontiguousarray(h1)
    eri = np.ascontiguousarray(eri)
    eri = ao2mo.restore('s8', eri, 8)
    ecore = float(ecore)

    active_mol = gto.M()
    active_mol.nelectron = sum(cas.nelecas) - 4

    mf_uhf = scf.UHF(active_mol)
    mf_uhf.nelec = (mol.nelec[0] - 2, mol.nelec[1] - 2)
    mf_uhf.get_hcore = lambda *args: np.asarray(h1)
    mf_uhf.get_ovlp = lambda *args: np.eye(8)
    mf_uhf.energy_nuc = lambda *args: ecore
    mf_uhf._eri = eri  # ao2mo.restore('8', np.zeros((8, 8, 8, 8)), 8)
    mf_uhf.init_guess = '1e'

    if previous_rho is None:
        dma = np.eye(8)
        dmb = np.eye(8)
        dmb[-1, -1] = 0
    else:
        dma, dmb = previous_rho
    mf_uhf.max_cycle = 200
    mf_uhf.verbose = 0
    mf_uhf.kernel((dma, dmb))

    if run_stability:
        mo1 = mf_uhf.stability()[0]
        dm1 = mf_uhf.make_rdm1(mo1, mf_uhf.mo_occ)
        mf_uhf = mf_uhf.run(dm1)
        mf_uhf.stability()
    return mf_uhf


def active_space_rhf(mol, mf, run_stability=False):

    cas = mcscf.CASCI(mf, ncas=8, nelecas=sum(mol.nelec) - 4)
    h1, ecore = cas.get_h1eff()
    eri = cas.get_h2cas()
    h1 = np.ascontiguousarray(h1)
    eri = np.ascontiguousarray(eri)
    eri = ao2mo.restore('s8', eri, 8)
    ecore = float(ecore)

    active_mol = gto.M()
    active_mol.nelectron = sum(cas.nelecas) - 4

    mf_rohf = scf.RHF(active_mol)
    mf_rohf.nelec = (mol.nelec[0] - 2, mol.nelec[1] - 2)
    mf_rohf.get_hcore = lambda *args: np.asarray(h1)
    mf_rohf.get_ovlp = lambda *args: np.eye(8)
    mf_rohf.energy_nuc = lambda *args: ecore
    mf_rohf._eri = eri  # ao2mo.restore('8', np.zeros((8, 8, 8, 8)), 8)
    # mf_rohf.init_guess = '1e'
    # mf_rohf.kernel()

    norb = h1.shape[1]
    nalpha, nbeta = mol.nelec[0] - 2, mol.nelec[1] - 2
    alpha_diag = [1] * nalpha + [0] * (norb - nalpha)
    beta_diag = [1] * nbeta + [0] * (norb - nbeta)
    mf_rohf.mo_coeff = np.eye(h1.shape[0])
    mf_rohf.mo_occ = np.array(alpha_diag) + np.array(beta_diag)
    w, v = np.linalg.eigh(mf_rohf.get_fock())
    mf_rohf.mo_energy = w
    mf_rohf.e_tot = mf_rohf.energy_tot()

    if run_stability:
        mo1 = mf_rohf.stability()[0]
        dm1 = mf_rohf.make_rdm1(mo1, mf_rohf.mo_occ)
        mf_rohf = mf_rohf.run(dm1)
        mf_rohf.stability()
    return mf_rohf

def get_rhf_active_space_hamiltonian_pyscf(mf):
    assert isinstance(mf, pyscf.scf.rhf.RHF)
    cas = mcscf.CASCI(mf, ncas=8, nelecas=sum(mf.nelec) - 4)
    h1, ecore = cas.get_h1eff()
    eri = cas.get_h2cas()
    h1 = np.ascontiguousarray(h1)
    eri = np.ascontiguousarray(eri)
    eri = ao2mo.restore('s1', eri, 8)
    ecore = float(ecore)

    # See PQRS convention in OpenFermion.hamiltonians._molecular_data
    # h[p,q,r,s] = (ps|qr)
    tei = np.asarray(
        eri.transpose(0, 2, 3, 1), order='C')
    ele_ham = integrals_to_fqe_restricted(h1, tei)
    soei, stei = spinorb_from_spatial(h1, tei)
    astei = np.einsum('ijkl', stei) - np.einsum('ijlk', stei)
    active_space_ham = of.InteractionOperator(ecore, soei, 0.25 * astei)
    return active_space_ham, ele_ham


def get_uhf_hamiltonian(mf):
    """Return OpenFermion InteractionOperator and FQE Hamiltonian"""
    assert isinstance(mf, pyscf.scf.uhf.UHF)
    # EXTRACT Hamiltonian for UHF
    norb = mf.mo_energy[0].size
    mo_a = mf.mo_coeff[0]
    mo_b = mf.mo_coeff[1]
    h1e_a = reduce(np.dot, (mo_a.T, mf.get_hcore(), mo_a))
    h1e_b = reduce(np.dot, (mo_b.T, mf.get_hcore(), mo_b))
    g2e_aa = ao2mo.incore.general(mf._eri, (mo_a,) * 4, compact=False)
    g2e_aa = g2e_aa.reshape(norb, norb, norb, norb)
    g2e_ab = ao2mo.incore.general(mf._eri, (mo_a, mo_a, mo_b, mo_b), compact=False)
    g2e_ab = g2e_ab.reshape(norb, norb, norb, norb)
    g2e_bb = ao2mo.incore.general(mf._eri, (mo_b,) * 4, compact=False)
    g2e_bb = g2e_bb.reshape(norb, norb, norb, norb)
    # h1e = (h1e_a, h1e_b)
    # eri = (g2e_aa, g2e_ab, g2e_bb)
    # print(h1e[0].shape)
    # print(eri[0].shape, eri[1].shape, eri[2].shape)

    # See PQRS convention in OpenFermion.hamiltonians._molecular_data
    # h[p,q,r,s] = (ps|qr)
    g2e_aa_of = np.asarray(1. * g2e_aa.transpose(0, 2, 3, 1), order='C')
    g2e_bb_of = np.asarray(1. * g2e_bb.transpose(0, 2, 3, 1), order='C')
    g2e_ab_of = np.asarray(1. * g2e_ab.transpose(0, 2, 3, 1), order='C')

    soei, stei = spinorb_from_spatial_blocks(h1e_a, h1e_b, g2e_aa_of, g2e_bb_of,
                                             g2e_ab_of)
    astei = np.einsum('ijkl', stei) - np.einsum('ijlk', stei)
    io_uhf_ham = of.InteractionOperator(0, soei, 0.25 * astei)
    fqe_uhf_ham = fqe.build_hamiltonian(of.get_fermion_operator(io_uhf_ham),
                                        norb=norb, conserve_number=True)
    return io_uhf_ham, fqe_uhf_ham


def main():

    bond_distances = np.linspace(3, 0.85, 20)
    rhf_energy = []
    rhf_cas_energy = []
    rhf_fci_cas_ovlp = []
    uhf_energy = []
    uhf_cas_energy = []
    uhf_fci_cas_ovlp = []
    previous_rho = None
    previous_rho_uhf = None
    for bd in bond_distances:
        mol = get_molecule(bd)
        mf = scf.RHF(mol)

        if previous_rho is None:
            mf.init_guess_by_atom()
            mf.run()
            mo1 = mf.stability()[0]
            dm1 = mf.make_rdm1(mo1, mf.mo_occ)
            mf = mf.run(dm1)
            mf.stability()
            # mf.analyze(verbose=5)
            previous_rho = mf.make_rdm1()

        else:
            mf.run(dm1=previous_rho)
            mo1 = mf.stability()[0]
            dm1 = mf.make_rdm1(mo1, mf.mo_occ)
            mf = mf.run(dm1)
            mf.stability()
            previous_rho = mf.make_rdm1()


        uhf_mf_active = active_space_uhf(mol, mf, run_stability=True, previous_rho=previous_rho_uhf)
        previous_rho_uhf = uhf_mf_active.make_rdm1()
        rhf_mf_active = active_space_rhf(mol, mf, run_stability=False)
        rhf_energy.append(rhf_mf_active.e_tot)
        uhf_energy.append(uhf_mf_active.e_tot)

        print("CALCULATION RESULTS")
        print()
        print("UHF-active Space E {: 5.15f}".format(
            uhf_mf_active.e_tot - uhf_mf_active.energy_nuc()))
        print("Core energy {: 5.15f}".format(uhf_mf_active.energy_nuc()))
        print("UHF Total {: 5.15f}".format(uhf_mf_active.e_tot))
        print()

        print("RHF-active Apace E {: 5.15f}".format(
            rhf_mf_active.e_tot - rhf_mf_active.energy_nuc()))
        print("Core energy {: 5.15f}".format(rhf_mf_active.energy_nuc()))
        print("RHF-active Total {: 5.15f}".format(rhf_mf_active.e_tot))
        print()

        # uhf_active_io_ham, uhf_active_fqe_ham = get_uhf_hamiltonian(
        #     uhf_mf_active)
        # rohf_active_io_ham, rohf_active_fqe_ham = \
        #     get_rhf_active_space_hamiltonian_pyscf(rhf_mf_active)

        uhf_active_cisolver = fci.FCI(uhf_mf_active)
        # uhf_active_cisolver = fci.addons.fix_spin(uhf_active_cisolver, shift=1,
        #                                           ss=0)
        rhf_active_cisolver = fci.FCI(rhf_mf_active)
        # rhf_active_cisolver = fci.addons.fix_spin(rhf_active_cisolver,
        #                                            shift=1, ss=0)
        uhf_fci_e, uhf_fci_civec = uhf_active_cisolver.kernel(
            nelec=uhf_mf_active.nelec)
        rohf_fci_e, rohf_fci_civec = rhf_active_cisolver.kernel(
            nelec=rhf_mf_active.nelec)

        uhf_fci_fqe_wf = pyscf_to_fqe_wf(uhf_fci_civec, norbs=8,
                                         nelec=uhf_mf_active.nelec)
        rhf_fci_fqe_wf = pyscf_to_fqe_wf(rohf_fci_civec, norbs=8,
                                          nelec=rhf_mf_active.nelec)
        rhf_fci_fqe_wf.print_wfn()
        uhf_cas_energy.append(uhf_fci_e)
        rhf_cas_energy.append(rohf_fci_e)

        print("\n\nFCI-HF overlaps\n")
        fqe_uhf_wf = fqe.Wavefunction([[sum(uhf_mf_active.nelec),
                                        uhf_mf_active.nelec[0] -
                                        uhf_mf_active.nelec[1],
                                        uhf_mf_active.mo_coeff[0].shape[1]]])
        fqe_rhf_wf = fqe.Wavefunction([[sum(rhf_mf_active.nelec),
                                         rhf_mf_active.nelec[0] -
                                         rhf_mf_active.nelec[1],
                                         rhf_mf_active.mo_coeff.shape[1]]])
        fqe_uhf_wf.set_wfn(strategy='hartree-fock')
        fqe_rhf_wf.set_wfn(strategy='hartree-fock')
        print("UHF-FCI Ovlp\t{: 5.15f}".format(
            fqe.vdot(fqe_uhf_wf, uhf_fci_fqe_wf)))
        print("RHF-FCI Ovlp\t{: 5.15f}".format(
            fqe.vdot(fqe_rhf_wf, rhf_fci_fqe_wf)))
        uhf_overlap = fqe.vdot(fqe_uhf_wf, uhf_fci_fqe_wf)
        rohf_overlap = fqe.vdot(fqe_rhf_wf, rhf_fci_fqe_wf)

        uhf_fci_cas_ovlp.append(np.abs(uhf_overlap))
        rhf_fci_cas_ovlp.append(np.abs(rohf_overlap))

    plt.plot(bond_distances, rhf_energy,'C0o--', label='RHF')
    plt.plot(bond_distances, uhf_energy,'C1o--', label='UHF-incas')
    plt.plot(bond_distances, rhf_cas_energy,'C0s-', label='CASCI-RHF')
    plt.plot(bond_distances, uhf_cas_energy,'C1s-', label='CASCI-UHF')
    plt.legend()
    plt.savefig("n2_rhf_uhf_cas.png", format="PNG", dpi=300)
    plt.show()


    plt.semilogy(bond_distances, rhf_fci_cas_ovlp, 'C0o--', label='RHF')
    plt.semilogy(bond_distances, uhf_fci_cas_ovlp, 'C1o--', label='UHF-incas')
    plt.legend()
    plt.savefig("n2_rhf_uhf_cas_ovlp.png", format="PNG", dpi=300)
    plt.show()


if __name__ == "__main__":
    main()
	from pyscf import gto, scf, mcscf, ao2mo, fci
	import pyscf
	from pyscf.fci.cistring import make_strings
	import numpy as np
	import openfermion as of
	from openfermion.chem.molecular_data import spinorb_from_spatial
	import fqe
	from fqe.openfermion_utils import integrals_to_fqe_restricted
	from functools import reduce
	import matplotlib.pyplot as plt


	def get_molecule(bd):
	mol = gto.M()
	mol.atom = "N 0 0 0; N 0 0 {}".format(bd)
	mol.basis = 'sto3g'
	mol.spin = 0
	mol.build()
	return mol


	def pyscf_to_fqe_wf(pyscf_cimat, pyscf_mf=None, norbs=None, nelec=None):
	if pyscf_mf is None:
	assert norbs is not None
	assert nelec is not None
	else:
	mol = pyscf_mf.mol
	nelec = mol.nelec
	norbs = pyscf_mf.mo_coeff.shape[1]

	norb_list = tuple(list(range(norbs)))
	n_alpha_strings = [x for x in make_strings(norb_list, nelec[0])]
	n_beta_strings = [x for x in make_strings(norb_list, nelec[1])]

	fqe_wf_ci = fqe.Wavefunction([[sum(nelec), nelec[0] - nelec[1], norbs]])
	fqe_data_ci = fqe_wf_ci.sector((sum(nelec), nelec[0] - nelec[1]))
	fqe_graph_ci = fqe_data_ci.get_fcigraph()
	fqe_orderd_coeff = np.zeros((fqe_graph_ci.lena(), fqe_graph_ci.lenb()), dtype=np.complex128)
	for paidx, pyscf_alpha_idx in enumerate(n_alpha_strings):
	for pbidx, pyscf_beta_idx in enumerate(n_beta_strings):
	# if np.abs(civec[paidx, pbidx]) > 1.0E-3:
	# print(np.binary_repr(pyscf_alpha_idx, width=10), np.binary_repr(pyscf_beta_idx, width=10), civec[paidx, pbidx])
	fqe_orderd_coeff[fqe_graph_ci.index_alpha(
	pyscf_alpha_idx), fqe_graph_ci.index_beta(pyscf_beta_idx)] = \
	pyscf_cimat[paidx, pbidx]

	fqe_data_ci.coeff = fqe_orderd_coeff
	return fqe_wf_ci


	def spinorb_from_spatial_blocks(h1a, h1b, eriaa, eribb, eriab,
	EQ_TOLERANCE=1.0E-12):
	n_qubits = 2 * h1a.shape[0]

	# Initialize Hamiltonian coefficients.
	one_body_coefficients = np.zeros((n_qubits, n_qubits))
	two_body_coefficients = np.zeros(
	(n_qubits, n_qubits, n_qubits, n_qubits))
	# Loop through integrals.
	for p in range(n_qubits // 2):
	for q in range(n_qubits // 2):

	# Populate 1-body coefficients. Require p and q have same spin.
	one_body_coefficients[2 * p, 2 * q] = h1a[p, q]
	one_body_coefficients[2 * p + 1, 2 * q + 1] = h1b[p, q]
	# Continue looping to prepare 2-body coefficients.
	for r in range(n_qubits // 2):
	for s in range(n_qubits // 2):
	# Mixed spin
	two_body_coefficients[2 * p, 2 * q + 1, 2 * r + 1, 2 * s] = (eriab[p, q, r, s])
	two_body_coefficients[2 * p + 1, 2 * q, 2 * r, 2 * s + 1] = (eriab[q, p, s, r])

	# Same spin
	two_body_coefficients[2 * p, 2 * q, 2 * r, 2 * s] = (eriaa[p, q, r, s])
	two_body_coefficients[2 * p + 1, 2 * q + 1, 2 * r + 1, 2 * s + 1] = (eribb[p, q, r, s])

	# Truncate.
	one_body_coefficients[
	np.absolute(one_body_coefficients) < EQ_TOLERANCE] = 0.
	two_body_coefficients[
	np.absolute(two_body_coefficients) < EQ_TOLERANCE] = 0.

	return one_body_coefficients, two_body_coefficients


	def active_space_uhf(mol, mf, run_stability=False, previous_rho=None):
	cas = mcscf.CASCI(mf, ncas=8, nelecas=sum(mol.nelec) - 4)
	h1, ecore = cas.get_h1eff()
	eri = cas.get_h2cas()
	h1 = np.ascontiguousarray(h1)
	eri = np.ascontiguousarray(eri)
	eri = ao2mo.restore('s8', eri, 8)
	ecore = float(ecore)

	active_mol = gto.M()
	active_mol.nelectron = sum(cas.nelecas) - 4

	mf_uhf = scf.UHF(active_mol)
	mf_uhf.nelec = (mol.nelec[0] - 2, mol.nelec[1] - 2)
	mf_uhf.get_hcore = lambda *args: np.asarray(h1)
	mf_uhf.get_ovlp = lambda *args: np.eye(8)
	mf_uhf.energy_nuc = lambda *args: ecore
	mf_uhf._eri = eri # ao2mo.restore('8', np.zeros((8, 8, 8, 8)), 8)
	mf_uhf.init_guess = '1e'

	if previous_rho is None:
	dma = np.eye(8)
	dmb = np.eye(8)
	dmb[-1, -1] = 0
	else:
	dma, dmb = previous_rho
	mf_uhf.max_cycle = 200
	mf_uhf.verbose = 0
	mf_uhf.kernel((dma, dmb))

	if run_stability:
	mo1 = mf_uhf.stability()[0]
	dm1 = mf_uhf.make_rdm1(mo1, mf_uhf.mo_occ)
	mf_uhf = mf_uhf.run(dm1)
	mf_uhf.stability()
	return mf_uhf


	def active_space_rhf(mol, mf, run_stability=False):

	cas = mcscf.CASCI(mf, ncas=8, nelecas=sum(mol.nelec) - 4)
	h1, ecore = cas.get_h1eff()
	eri = cas.get_h2cas()
	h1 = np.ascontiguousarray(h1)
	eri = np.ascontiguousarray(eri)
	eri = ao2mo.restore('s8', eri, 8)
	ecore = float(ecore)

	active_mol = gto.M()
	active_mol.nelectron = sum(cas.nelecas) - 4

	mf_rohf = scf.RHF(active_mol)
	mf_rohf.nelec = (mol.nelec[0] - 2, mol.nelec[1] - 2)
	mf_rohf.get_hcore = lambda *args: np.asarray(h1)
	mf_rohf.get_ovlp = lambda *args: np.eye(8)
	mf_rohf.energy_nuc = lambda *args: ecore
	mf_rohf._eri = eri # ao2mo.restore('8', np.zeros((8, 8, 8, 8)), 8)
	# mf_rohf.init_guess = '1e'
	# mf_rohf.kernel()

	norb = h1.shape[1]
	nalpha, nbeta = mol.nelec[0] - 2, mol.nelec[1] - 2
	alpha_diag = [1] * nalpha + [0] * (norb - nalpha)
	beta_diag = [1] * nbeta + [0] * (norb - nbeta)
	mf_rohf.mo_coeff = np.eye(h1.shape[0])
	mf_rohf.mo_occ = np.array(alpha_diag) + np.array(beta_diag)
	w, v = np.linalg.eigh(mf_rohf.get_fock())
	mf_rohf.mo_energy = w
	mf_rohf.e_tot = mf_rohf.energy_tot()

	if run_stability:
	mo1 = mf_rohf.stability()[0]
	dm1 = mf_rohf.make_rdm1(mo1, mf_rohf.mo_occ)
	mf_rohf = mf_rohf.run(dm1)
	mf_rohf.stability()
	return mf_rohf

	def get_rhf_active_space_hamiltonian_pyscf(mf):
	assert isinstance(mf, pyscf.scf.rhf.RHF)
	cas = mcscf.CASCI(mf, ncas=8, nelecas=sum(mf.nelec) - 4)
	h1, ecore = cas.get_h1eff()
	eri = cas.get_h2cas()
	h1 = np.ascontiguousarray(h1)
	eri = np.ascontiguousarray(eri)
	eri = ao2mo.restore('s1', eri, 8)
	ecore = float(ecore)

	# See PQRS convention in OpenFermion.hamiltonians._molecular_data
	# h[p,q,r,s] = (ps\|qr)
	tei = np.asarray(
	eri.transpose(0, 2, 3, 1), order='C')
	ele_ham = integrals_to_fqe_restricted(h1, tei)
	soei, stei = spinorb_from_spatial(h1, tei)
	astei = np.einsum('ijkl', stei) - np.einsum('ijlk', stei)
	active_space_ham = of.InteractionOperator(ecore, soei, 0.25 * astei)
	return active_space_ham, ele_ham


	def get_uhf_hamiltonian(mf):
	"""Return OpenFermion InteractionOperator and FQE Hamiltonian"""
	assert isinstance(mf, pyscf.scf.uhf.UHF)
	# EXTRACT Hamiltonian for UHF
	norb = mf.mo_energy[0].size
	mo_a = mf.mo_coeff[0]
	mo_b = mf.mo_coeff[1]
	h1e_a = reduce(np.dot, (mo_a.T, mf.get_hcore(), mo_a))
	h1e_b = reduce(np.dot, (mo_b.T, mf.get_hcore(), mo_b))
	g2e_aa = ao2mo.incore.general(mf._eri, (mo_a,) * 4, compact=False)
	g2e_aa = g2e_aa.reshape(norb, norb, norb, norb)
	g2e_ab = ao2mo.incore.general(mf._eri, (mo_a, mo_a, mo_b, mo_b), compact=False)
	g2e_ab = g2e_ab.reshape(norb, norb, norb, norb)
	g2e_bb = ao2mo.incore.general(mf._eri, (mo_b,) * 4, compact=False)
	g2e_bb = g2e_bb.reshape(norb, norb, norb, norb)
	# h1e = (h1e_a, h1e_b)
	# eri = (g2e_aa, g2e_ab, g2e_bb)
	# print(h1e[0].shape)
	# print(eri[0].shape, eri[1].shape, eri[2].shape)

	# See PQRS convention in OpenFermion.hamiltonians._molecular_data
	# h[p,q,r,s] = (ps\|qr)
	g2e_aa_of = np.asarray(1. * g2e_aa.transpose(0, 2, 3, 1), order='C')
	g2e_bb_of = np.asarray(1. * g2e_bb.transpose(0, 2, 3, 1), order='C')
	g2e_ab_of = np.asarray(1. * g2e_ab.transpose(0, 2, 3, 1), order='C')

	soei, stei = spinorb_from_spatial_blocks(h1e_a, h1e_b, g2e_aa_of, g2e_bb_of,
	g2e_ab_of)
	astei = np.einsum('ijkl', stei) - np.einsum('ijlk', stei)
	io_uhf_ham = of.InteractionOperator(0, soei, 0.25 * astei)
	fqe_uhf_ham = fqe.build_hamiltonian(of.get_fermion_operator(io_uhf_ham),
	norb=norb, conserve_number=True)
	return io_uhf_ham, fqe_uhf_ham


	def main():

	bond_distances = np.linspace(3, 0.85, 20)
	rhf_energy = []
	rhf_cas_energy = []
	rhf_fci_cas_ovlp = []
	uhf_energy = []
	uhf_cas_energy = []
	uhf_fci_cas_ovlp = []
	previous_rho = None
	previous_rho_uhf = None
	for bd in bond_distances:
	mol = get_molecule(bd)
	mf = scf.RHF(mol)

	if previous_rho is None:
	mf.init_guess_by_atom()
	mf.run()
	mo1 = mf.stability()[0]
	dm1 = mf.make_rdm1(mo1, mf.mo_occ)
	mf = mf.run(dm1)
	mf.stability()
	# mf.analyze(verbose=5)
	previous_rho = mf.make_rdm1()

	else:
	mf.run(dm1=previous_rho)
	mo1 = mf.stability()[0]
	dm1 = mf.make_rdm1(mo1, mf.mo_occ)
	mf = mf.run(dm1)
	mf.stability()
	previous_rho = mf.make_rdm1()


	uhf_mf_active = active_space_uhf(mol, mf, run_stability=True, previous_rho=previous_rho_uhf)
	previous_rho_uhf = uhf_mf_active.make_rdm1()
	rhf_mf_active = active_space_rhf(mol, mf, run_stability=False)
	rhf_energy.append(rhf_mf_active.e_tot)
	uhf_energy.append(uhf_mf_active.e_tot)

	print("CALCULATION RESULTS")
	print()
	print("UHF-active Space E {: 5.15f}".format(
	uhf_mf_active.e_tot - uhf_mf_active.energy_nuc()))
	print("Core energy {: 5.15f}".format(uhf_mf_active.energy_nuc()))
	print("UHF Total {: 5.15f}".format(uhf_mf_active.e_tot))
	print()

	print("RHF-active Apace E {: 5.15f}".format(
	rhf_mf_active.e_tot - rhf_mf_active.energy_nuc()))
	print("Core energy {: 5.15f}".format(rhf_mf_active.energy_nuc()))
	print("RHF-active Total {: 5.15f}".format(rhf_mf_active.e_tot))
	print()

	# uhf_active_io_ham, uhf_active_fqe_ham = get_uhf_hamiltonian(
	# uhf_mf_active)
	# rohf_active_io_ham, rohf_active_fqe_ham = \
	# get_rhf_active_space_hamiltonian_pyscf(rhf_mf_active)

	uhf_active_cisolver = fci.FCI(uhf_mf_active)
	# uhf_active_cisolver = fci.addons.fix_spin(uhf_active_cisolver, shift=1,
	# ss=0)
	rhf_active_cisolver = fci.FCI(rhf_mf_active)
	# rhf_active_cisolver = fci.addons.fix_spin(rhf_active_cisolver,
	# shift=1, ss=0)
	uhf_fci_e, uhf_fci_civec = uhf_active_cisolver.kernel(
	nelec=uhf_mf_active.nelec)
	rohf_fci_e, rohf_fci_civec = rhf_active_cisolver.kernel(
	nelec=rhf_mf_active.nelec)

	uhf_fci_fqe_wf = pyscf_to_fqe_wf(uhf_fci_civec, norbs=8,
	nelec=uhf_mf_active.nelec)
	rhf_fci_fqe_wf = pyscf_to_fqe_wf(rohf_fci_civec, norbs=8,
	nelec=rhf_mf_active.nelec)
	rhf_fci_fqe_wf.print_wfn()
	uhf_cas_energy.append(uhf_fci_e)
	rhf_cas_energy.append(rohf_fci_e)

	print("\n\nFCI-HF overlaps\n")
	fqe_uhf_wf = fqe.Wavefunction([[sum(uhf_mf_active.nelec),
	uhf_mf_active.nelec[0] -
	uhf_mf_active.nelec[1],
	uhf_mf_active.mo_coeff[0].shape[1]]])
	fqe_rhf_wf = fqe.Wavefunction([[sum(rhf_mf_active.nelec),
	rhf_mf_active.nelec[0] -
	rhf_mf_active.nelec[1],
	rhf_mf_active.mo_coeff.shape[1]]])
	fqe_uhf_wf.set_wfn(strategy='hartree-fock')
	fqe_rhf_wf.set_wfn(strategy='hartree-fock')
	print("UHF-FCI Ovlp\t{: 5.15f}".format(
	fqe.vdot(fqe_uhf_wf, uhf_fci_fqe_wf)))
	print("RHF-FCI Ovlp\t{: 5.15f}".format(
	fqe.vdot(fqe_rhf_wf, rhf_fci_fqe_wf)))
	uhf_overlap = fqe.vdot(fqe_uhf_wf, uhf_fci_fqe_wf)
	rohf_overlap = fqe.vdot(fqe_rhf_wf, rhf_fci_fqe_wf)

	uhf_fci_cas_ovlp.append(np.abs(uhf_overlap))
	rhf_fci_cas_ovlp.append(np.abs(rohf_overlap))

	plt.plot(bond_distances, rhf_energy,'C0o--', label='RHF')
	plt.plot(bond_distances, uhf_energy,'C1o--', label='UHF-incas')
	plt.plot(bond_distances, rhf_cas_energy,'C0s-', label='CASCI-RHF')
	plt.plot(bond_distances, uhf_cas_energy,'C1s-', label='CASCI-UHF')
	plt.legend()
	plt.savefig("n2_rhf_uhf_cas.png", format="PNG", dpi=300)
	plt.show()


	plt.semilogy(bond_distances, rhf_fci_cas_ovlp, 'C0o--', label='RHF')
	plt.semilogy(bond_distances, uhf_fci_cas_ovlp, 'C1o--', label='UHF-incas')
	plt.legend()
	plt.savefig("n2_rhf_uhf_cas_ovlp.png", format="PNG", dpi=300)
	plt.show()


	if __name__ == "__main__":
	main()