When running the VQE simulation for $H_2$ over a range of distances as detailed in the Simulating Molecules using VQE section of the Qiskit Textbook, the VQE energy deviation from the exact energy increases with distance as shown in the figure below: H_2 Energy(Hartree ) vs Distance(Angstrom)

How is this relation($Error \propto distance$) explained and how can it's effects be minimized?

My VQE code is as follows:

distances = np.arange(0.2, 5, 0.5) # for a finer exact_energies plot run with step of 0.1
exact_energies = []
vqe_energies = []

optimizer = SPSA(maxiter=1000)

for dist in distances:
    molecule = "H .0 .0 -" + str(dist) + "; H .0 .0 " + str(dist)
    driver = PySCFDriver(atom = molecule, unit=UnitsType.ANGSTROM, charge=0, spin=0, basis='sto3g')
    qmolecule = driver.run()
    num_particles = qmolecule.num_alpha + qmolecule.num_beta
    qubitOp = FermionicOperator(h1=qmolecule.one_body_integrals, h2=qmolecule.two_body_integrals).mapping(map_type='parity')
    qubitOp = Z2Symmetries.two_qubit_reduction(qubitOp, num_particles)
    result = NumPyEigensolver(qubitOp).run()
    var_form = EfficientSU2(qubitOp.num_qubits, entanglement="linear")
    vqe = VQE(qubitOp, var_form, optimizer)
    vqe_result = np.real(vqe.run(backend)['eigenvalue'])

I expected the error to be very low as maxiter = 1000. Also note the backend is a statevector_simulator.

  • $\begingroup$ Also I found it strange that the exact energy plot for $H_2$ does not have the concave shape as other curves such as $LiH$ as can be seen in the aforementioned section in the Qiskit Textbook. $\endgroup$
    – MShakeG
    Commented Sep 28, 2020 at 19:01

2 Answers 2


The reason why you don't have the concave curve is because you are not taken into account of the nuclear repulsion (energy shift) term.

Note that the electronic Hamiltonian (using Born-Oppenheimer approximation) is:

$$ H = -\sum_i \dfrac{\nabla^2_{r_i}}{2} - \sum \dfrac{Z_i}{|R_i - r_j|} + \sum \dfrac{1}{|r_i - r_j|} + \sum \dfrac{Z_iZ_j}{|R_i - R_j|} $$

Since the nuclei are fixed, the last term is a constant. This term is the nuclear repulsion. You must add it in at the end to get the electronic ground state energy.

As to why you get better result at some distance compare to other is because of the structure of your var_form. You used the same var_form for every geometric configuration here. But in reality, a var_form might be good at one geometry configuration and bad at a different geometry configuration. This is because you don't expect that single circuit structure can prepare an arbitrary quantum state... At some distance, the state might be more entangled, and the var_form you picked might not be suffice to generate that state.

  • $\begingroup$ thank you the detailed answer, I tried using the get_qubit_op function defined in the aforementioned chapter of the textbook for $H_2$ in order to obtain the shift, however unsuccessfully. I could only get it to work for $LiH$ and $BeH_2$. $\endgroup$
    – MShakeG
    Commented Sep 29, 2020 at 15:28
  • $\begingroup$ Rather than using defining qubitOp through FermionicOperator, you can use the Hamiltonian class that is built in Qiskit. You can extract the nuclear repulsion energy from there. Here is the link: qiskit.org/documentation/stubs/… $\endgroup$
    – KAJ226
    Commented Sep 29, 2020 at 16:56

I managed to resolve the issue; I'm not sure if this applies to the other VQE optimizers(e.g. SLSQP) but by placing optimizer = SPSA(maxiter=1000) inside the loop in order to re-initialize for each iteration of distance the issue(i.e. $Error∝distance$) is resolved.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.