Which Hamiltonian object to use #102
Description
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