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

Question

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

Answer 1

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.