# Why VQE gives an incorrect ground state energy for a chain of 4 hydrogen atoms?

I am a bit hesitant to ask this very specific question, as I feel other people need not benefit from it. But since I have struggled for a while, and I think I should get some help.

So I am using VQE in qiskit to calculate the ground sate energy of a chain of hydrogen atoms, but it appears that the result is in-consistent with the result from exact diagonalization. The code works well for other molecules like H2, LiH, so this is confusing. I guess the question boils down to how to set the threshold for the VQE. I have attached the code below, and many thanks for the help!

from qiskit import BasicAer
import logging
from qiskit.chemistry import set_qiskit_chemistry_logging
set_qiskit_chemistry_logging(logging.ERROR)

# chemistry related modules
from qiskit.chemistry import FermionicOperator
from qiskit.chemistry.drivers import PySCFDriver, UnitsType

from qiskit.aqua.algorithms import VQE, NumPyEigensolver
import numpy as np
from qiskit.chemistry.components.variational_forms import UCCSD
from qiskit.chemistry.components.initial_states import HartreeFock
from qiskit.aqua.components.optimizers import L_BFGS_B
from qiskit.aqua.operators import Z2Symmetries

def get_qubit_op(atom,basis,map_type ):
driver = PySCFDriver(atom=atom, unit=UnitsType.ANGSTROM,
charge=0, spin=0, basis=basis)
molecule = driver.run()
num_particles = molecule.num_alpha + molecule.num_beta
num_spin_orbitals = molecule.num_orbitals * 2
ferOp = FermionicOperator(h1=molecule.one_body_integrals, h2=molecule.two_body_integrals)

qubitOp = ferOp.mapping(map_type=map_type, threshold=0.00000001)
qubitOp = Z2Symmetries.two_qubit_reduction(qubitOp, num_particles)

return qubitOp, num_particles, num_spin_orbitals

import timeit

start = timeit.default_timer()

atom = 'H .0 .0 .0; H .0 .0 1.5 ; H .0 .0 3.0 ; H .0 .0 4.5 '

basis='sto3g'

map_type = 'parity'

qubitOp, num_particles, num_spin_orbitals = get_qubit_op(atom,basis,map_type )

print('Ground state energy is' , NumPyEigensolver( qubitOp ).run().eigenvalues )

init_state = HartreeFock( num_spin_orbitals , num_particles , map_type )

# set the backend for the quantum computation=
backend = BasicAer.get_backend('statevector_simulator')

# setup a classical optimizer for VQE
optimizer = L_BFGS_B()
print( init_state.bitstr )

var_form_vqe = UCCSD(
num_orbitals=num_spin_orbitals,
num_particles=num_particles,
initial_state=init_state,
qubit_mapping=map_type
)

algorithm_vqe = VQE(qubitOp, var_form_vqe, optimizer  )

result_vqe = algorithm_vqe.run(backend)

print( 'eigenvalue = ' , result_vqe['eigenvalue' ] )

stop = timeit.default_timer()

print('The run time of this part: ', stop - start)


The output is below, and as you can see, they differ quite significantly.

Ground state energy is [-3.52488449+5.88070795e-18j]
[False False  True False False  True]
eigenvalue =  (-3.523526951494827+0j)
The run time of this part:  57.303660957000034

• I can't see this setting here, but I would begin with increasing the trotter order (i.e. having more runs of the ansatz circuit). But may only on a simulator, the circuit will be too big for existing ibmq chips. Oct 6 '20 at 22:31
• Hi @mavzolej, thanks for the comments. Actually this is one of my puzzles for such simulation, is this code run on an actual quantum machine, or is it just a simulator? Here the quantum circuit only involve 6 qubits, so in either case, it should be within the limit. What other specific setting of the VQE would you like to know? What do you mean by "having more runs of the ansatz circuit"? Thanks!
– fagd
Oct 6 '20 at 22:41
• backend = BasicAer.get_backend('statevector_simulator') - in this line you choose the classical simulator to be your backend. So, all the observables are evaluated exactly. The next step could be the QASM backend which would calculate observables using sampling from the exact probability distribution. And then you move to the real hardware. Oct 6 '20 at 23:06
• I would suggest that you to read about UCC here and also google smth on trotterization, so that you understood what's meant by "trotter order" in Ibid. Oct 6 '20 at 23:08
• Thanks for both the comments. Now I understand what you mean by trotterization order. It seems to corresponds to the "reps" parameter in the UCCSD, and I will try it out. Thanks for reminding me that.
– fagd
Oct 6 '20 at 23:27

I think this has to do with the fact that the var_form that you picked, UCCSD, is not expressible enough. What I mean by that is UCCSD doesn't have the sufficient complexity to generate the ground state wavefunction.

Maybe you should try UCCSDT or some other var_form that have greater expressibility.

• Hi @KAJ226, thanks for the answer. It is a bit surprised that UCCSD is not general enough, given that H_4 is a pretty small system. For "UCCSDT" you mentioned, I could not find it in the documentation of qiskit, could you please elaborate a bit more? Thanks!
– fagd
Oct 6 '20 at 22:44
• I don't think UCCSDT is available in Qiskit. But I think "Single" and "Double" excitation might be enough here, actually...As suggested by the comment to your question, increasing the trotter order will help. As UCCSD operator is actually an exponential of a matrix with terms that do not commute. Therefore it has to be approximated. Qiskit takes the approximation order to be 1 if you don't specify it specifically. Here is the link to the source code: qiskit.org/documentation/stubs/… Oct 6 '20 at 23:18
• I guess you mean the parameter "reps" in UCCSD. Thanks for pointing that out, as I almost forgot in the trotterization, there is a factor of n in the exponential. I will try to make the parameter larger.
– fagd
Oct 6 '20 at 23:26
• I mean the "num_time_slice". Although I don't know why they have the parameter "reps" for UCCSD to be honest . That is usually for those heuristic ansatz. For UCCSD , which is a chemically motivated ansatz, it doesn't make sense to repeat it multiple time in that sense, in my opinion. I didn't get any better result though. It is about the same as before. Maybe you did. Oct 6 '20 at 23:36
• No, I tried reps = 2,3, the result is still roughly -3.5236... If you try different "num_time_slice" and still get the same result, then it is really confusing.
– fagd
Oct 6 '20 at 23:42