How to calculate the Hartree-Fock energy in Qiskit?

I am wondering how to calculate the Hartree-Fock energy efficiently using Qiskit. For one can get the Hartree-Fock state efficiently in Qiskit, however, it seems that it is not so obvious to get the energy efficiently using the state.

• What do you mean? Aug 28 '21 at 16:17
• Please clarify what you want to do. Do you want to calculate an estimate ground state using HF Or do you want to calculate the ground state of some molecule? If yes to the latter, do you want an exact estimation or VQE? Aug 29 '21 at 0:46
• I want to calculate the energy using the Hartree-Fock method. Aug 29 '21 at 11:58
• Confirm if this is what you want, in this order: (1) load the data of some specific molecule (2) produce a circuit that generates the Hartree-Fock ground state (3) read the energy of associated to that state. Aug 29 '21 at 12:27
• My question is about the third step that how to read the energy associated to the state. Sep 1 '21 at 13:07

As you are asking specifically for the evaluation of the energy only, I will be brief. I will assume that you have a init_state (a quantum circuit) that produces the the Hartree-Fock wavefunction or any other wavefunction you like to test. I could not find a Qiskit function that provides the energy expectation value of a given wavefunction, given some molecular Hamiltonian. You could look at the source code of Qiskit VQE, but here I propose a quick solution using qiskit.algorithms.VQE directly :

from qiskit.algorithms import VQE
algorithm = VQE(ansatz,quantum_instance)
result = algorithm.compute_minimum_eigenvalue(qubit_op).eigenvalue
print(result)

where quantum_instance is the backend, ansatz is your Hartree-Fock circuit and qubit_op is of class OperatorBase and is basically the operators after mapping. This answer is based on , for more details on qubit_op see section 3.

For ansatz you will need a parametrized circuit, but you can just add a dummy parameter to your init_state. Something like:

from qiskit.circuit import Parameter, QuantumCircuit, QuantumRegister
theta = Parameter('a')
n = qubit_op.num_qubits
qc = QuantumCircuit(qubit_op.num_qubits)
qc.rz(theta*0,0)
ansatz = qc
ansatz.compose(init_state, front=True, inplace=True)

This answers is inspired from  IBM Quantum Challenge 2021, problem 5.

Minimum working example

Here I leave an example for the $$\mathrm{H}_2$$ molecule, qubitized using Jordan–Wigner mapping:

from qiskit import Aer
from qiskit_nature.drivers import PySCFDriver
from qiskit_nature.problems.second_quantization.electronic import ElectronicStructureProblem
from qiskit_nature.mappers.second_quantization import JordanWignerMapper
from qiskit_nature.converters.second_quantization.qubit_converter import QubitConverter
from qiskit_nature.circuit.library import HartreeFock
from qiskit.circuit import Parameter, QuantumCircuit, QuantumRegister
from qiskit.algorithms import VQE

#Choose a backend
backend = Aer.get_backend('statevector_simulator')

molecule = "H .0 .0 .0; H .0 .0 0.739"
driver = PySCFDriver(atom=molecule)
qmolecule = driver.run()

#Build the electronic structure problem
problem = ElectronicStructureProblem(driver)

# Generate the second-quantized operators
second_q_ops = problem.second_q_ops()

# Hamiltonian
main_op = second_q_ops

#Choose some mapping and conver it to qubits
mapper = JordanWignerMapper()
converter = QubitConverter(mapper=mapper)

# The fermionic operators are mapped to qubit operators
num_particles = (problem.molecule_data_transformed.num_alpha,
problem.molecule_data_transformed.num_beta)
qubit_op = converter.convert(main_op, num_particles=num_particles)
#Create the Hartree-Fock circuit

num_particles = (problem.molecule_data_transformed.num_alpha,
problem.molecule_data_transformed.num_beta)
num_spin_orbitals = 2 * problem.molecule_data_transformed.num_molecular_orbitals
init_state = HartreeFock(num_spin_orbitals, num_particles, converter)

#Create dummy parametrized circuit
theta = Parameter('a')
n = qubit_op.num_qubits
qc = QuantumCircuit(qubit_op.num_qubits)
qc.rz(theta*0,0)
ansatz = qc
ansatz.compose(init_state, front=True, inplace=True)

#Pass it through VQE
algorithm = VQE(ansatz,quantum_instance=backend)
result = algorithm.compute_minimum_eigenvalue(qubit_op).eigenvalue
print(result)

Which results in

(-1.832882751436994+0j)

which is in hartrees.

Warning: Qiskit Nature has been updated since  and some of the code presented here may not run depending on the version (but it works now online on the IBM Quantum Lab).