IndexError using QubitConverter to get a qubit operator

cpossel commented 3 years ago

Information

Qiskit Terra version: 0.17.4
Python version: 3.8.10
Operating system: Windows 10

What is the current behavior?

Script terminates with IndexError when converting a ElectronicStructureProblem to a qubit operator using convert method from a QubitConverter object. Here is the whole traceback:

Traceback (most recent call last):
  File "c:\Users\poc\Documents\qiskit_scripts\test_with_error.py", line 48, in <module>
    qubit_op = converter.convert(main_op, num_particles=num_particles)
  File "C:\Users\poc\Anaconda3\envs\qiskit-0-27-psi4-jupyter\lib\site-packages\qiskit_nature\converters\second_quantization\qubit_converter.py", line 169, in convert
    reduced_op = self._two_qubit_reduce(qubit_op, num_particles)
  File "C:\Users\poc\Anaconda3\envs\qiskit-0-27-psi4-jupyter\lib\site-packages\qiskit_nature\converters\second_quantization\qubit_converter.py", line 272, in _two_qubit_reduce
    reduced_op = cast(PauliSumOp, two_q_reducer.convert(qubit_op))
  File "C:\Users\poc\Anaconda3\envs\qiskit-0-27-psi4-jupyter\lib\site-packages\qiskit\opflow\converters\two_qubit_reduction.py", line 95, in convert
    return z2_symmetries.taper(operator)
  File "C:\Users\poc\Anaconda3\envs\qiskit-0-27-psi4-jupyter\lib\site-packages\qiskit\opflow\primitive_ops\tapered_pauli_sum_op.py", line 360, in taper
    tapered_ops = self._taper(operator, self._tapering_values)
  File "C:\Users\poc\Anaconda3\envs\qiskit-0-27-psi4-jupyter\lib\site-packages\qiskit\opflow\primitive_ops\tapered_pauli_sum_op.py", line 377, in _taper
    spo = SparsePauliOp.from_list(pauli_list).simplify(atol=0.0)
  File "C:\Users\poc\Anaconda3\envs\qiskit-0-27-psi4-jupyter\lib\site-packages\qiskit\quantum_info\operators\symplectic\sparse_pauli_op.py", line 378, in from_list
    num_qubits = len(PauliTable._from_label(obj[0][0]))
  File "C:\Users\poc\Anaconda3\envs\qiskit-0-27-psi4-jupyter\lib\site-packages\qiskit\quantum_info\operators\symplectic\pauli_table.py", line 928, in _from_label
    if label[0] == '+':
IndexError: string index out of range

Steps to reproduce the problem

Run the example file in the README.md from https://github.com/Qiskit/qiskit-nature, but use a single helium atom (He) instead of H2.

In my case I slightly adjusted the example (switch to PSI4), so here's the script I used:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import psi4
from qiskit.algorithms.optimizers import L_BFGS_B, COBYLA
from qiskit.utils import QuantumInstance
from qiskit_nature.converters.second_quantization import QubitConverter
from qiskit_nature.drivers import Molecule, PSI4Driver
from qiskit_nature.drivers.fermionic_driver import HFMethodType
from qiskit_nature.mappers.second_quantization import ParityMapper
from qiskit_nature.problems.second_quantization.electronic import ElectronicStructureProblem

# Use PSI4, a classical computational chemistry software
# package, to compute the one-body and two-body integrals in
# electronic-orbital basis, necessary to form the Fermionic operator
# molecule = Molecule(  # H2 --> works perfectly fine (E=-1.137306035696Hartree)
#     geometry=[["H", [0.0, 0.0, 0.0]], ["H", [0.0, 0.0, 0.735]]], charge=0, multiplicity=1,
# )
molecule = Molecule(geometry=[["He", [0.0, 0.0, 0.0]]], charge=0, multiplicity=1,)  # He --> error
if molecule.multiplicity == 1:
    hf_method = HFMethodType.RHF
else:
    hf_method = HFMethodType.UHF
driver = PSI4Driver(
    molecule=molecule,
    basis="sto-3g",  # basis sets (from fast to accurate): "sto-3g", "3-21g", "6-31g", ...
    hf_method=hf_method,
)

problem = ElectronicStructureProblem(driver)

# generate the second-quantized operators
second_q_ops = problem.second_q_ops()
main_op = second_q_ops[0]

num_particles = (
    problem.molecule_data_transformed.num_alpha,
    problem.molecule_data_transformed.num_beta,
)

num_spin_orbitals = 2 * problem.molecule_data.num_molecular_orbitals

# setup the mapper and qubit converter
mapper = ParityMapper()
converter = QubitConverter(mapper=mapper, two_qubit_reduction=True)

# map to qubit operators
qubit_op = converter.convert(main_op, num_particles=num_particles)

# setup the initial state for the ansatz
from qiskit_nature.circuit.library import HartreeFock

# Ansatzes(?)/initial states (from fast to accurate): HartreeFock , ...(?), UCC, ...(?)
init_state = HartreeFock(num_spin_orbitals, num_particles, converter)

# setup the ansatz for VQE
from qiskit.circuit.library import TwoLocal

ansatz = TwoLocal(num_spin_orbitals, ["ry", "rz"], "cz")

# add the initial state
ansatz.compose(init_state, front=True, inplace=True)

# setup the classical optimizer for VQE
optimizer = L_BFGS_B(iprint=0)
# optimizer = COBYLA(disp=True)

from qiskit import Aer

sim_backend = Aer.get_backend("statevector_simulator")
backend = sim_backend

quantum_instance = QuantumInstance(backend=backend, optimization_level=0,)

# setup and run VQE
from qiskit.algorithms import VQE

algorithm = VQE(ansatz, optimizer=optimizer, quantum_instance=quantum_instance)

result = algorithm.compute_minimum_eigenvalue(qubit_op)
print(result.eigenvalue.real)

electronic_structure_result = problem.interpret(result)
print(electronic_structure_result)

The commented molecule H2 works perfectly fine, but some other molecules (He, H, H-, Ne) throw the mentioned error.

What is the expected behavior?

Convert the ElectronicStructureProblem to a qubit operator without error. In case of the example script, finally print the ground state energy.

Suggested solutions

Actually I have no idea, where exactly the problem comes from. It might even be a bug in earlier parts that propagate till pauli_table.py, and thus might belong to the qiskit-nature package. Though I don't have enough insight to figure that out. Any help and suggestions would be appreciated!

woodsp-ibm commented 3 years ago

I moved it to Nature - you may be right that the fix may come down to logic that is in Terra, but I think Nature is the 'right' place to have this for now as things are investigated

cpossel commented 3 years ago

Okay. Then the whole version info might be useful instead of only the Terra version:

mrossinek commented 3 years ago

The fact that this is an IndexError when accessing label[0] seems to suggest that the length of label is 0 which would mean that the number of qubits is 0. This sounds correct when simulating a single He atom in combination with two_qubit_reduction=True. We can add a check to the QubitConverter to skip this reduction when there are 2 or fewer qubits to begin with.

However, I don't think simulating e.g. Ne should cause this. Could you double check that please?

cpossel commented 3 years ago

I rechecked it and you're completely right. Ne works fine. In the other cases (H, H-, He) the error still persists. I tested to switch to two_qubit_reduction=False and then also H, H-, and He worked fine in agreement with your explanation. I think the check you proposed (or even only a warning message) would be helpful for all those who start with such small examples and try out the possible options without (yet) understanding all of them in detail.

qiskit-community / qiskit-nature