OpenFermion-FQE icon indicating copy to clipboard operation
OpenFermion-FQE copied to clipboard

Which Hamiltonian object to use

Open ncrubin opened this issue 4 years ago • 3 comments

If I want to build a Hamiltonian from UHF I would do the following in pyscf + OpenFermion + FQE. This uses the build_hamiltonian function to translate the FermionOperator into an FQE Hamiltonian. A more direct way would be to populate a FQE Hamiltonian directly (a SSO Hamiltonian object?).

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

where spinorb_from_spatial_blocks builds the spin-orbital Hamiltonian from h1e_a, h1e_b, g2eaa, g2e_bb, etc...

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

ncrubin avatar Jul 19 '21 15:07 ncrubin

Let me cc @klgunst and @awhite862

shiozaki avatar Jul 19 '21 15:07 shiozaki

To be more specific, which Hamiltonian object in FQE should I try to populate with the h1a, h1b, teiaa, teiab, teibb tensors?

ncrubin avatar Jul 19 '21 16:07 ncrubin

It is a SSO hamiltonian. The dimensions of the arrays have to be (norb*2, norb*2) and (norb*2, norb*2, norb*2, norb*2). I believe that h1a and h1b should be (:norb, :norb) and (norb:, norb:) part of the one body, teiaa, teiab, teibb are (:norb, :norb, :norb, :norb), [(:norb, norb:, :norb, norb:) and (norb:, :norb, norb:, :norb)], and (norb:, norb:, norb:, norb:) respectively.

shiozaki avatar Jul 19 '21 19:07 shiozaki