Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Which Hamiltonian object to use #102

Open
ncrubin opened this issue Jul 19, 2021 · 3 comments
Open

Which Hamiltonian object to use #102

ncrubin opened this issue Jul 19, 2021 · 3 comments

Comments

@ncrubin
Copy link
Collaborator

ncrubin commented Jul 19, 2021

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
@shiozaki
Copy link
Collaborator

Let me cc @klgunst and @awhite862

@ncrubin
Copy link
Collaborator Author

ncrubin commented Jul 19, 2021

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

@shiozaki
Copy link
Collaborator

shiozaki commented Jul 19, 2021

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants