Calculating reduced density matrices from VQE-UCCSD: Error operator not hermitian.

[abhishekkhedkar09]

Environment

• Qiskit Nature version:
• Python version:
• Operating system:

What is happening?

I get the following error, please find below my input file.

If I check for `es_problem.second_q_ops()[1]['RDM(1, 0)'].is_hermitian() for example, it does return False.

the corresponding JW mapped qubit_op is as well not unitary - checked.

How can one compute 1RDM's from UCCSD simulation?

Traceback (most recent call last):
  File "/home/abhishek/envs/gpuqiskit/lib/python3.10/site-packages/qiskit_algorithms/observables_evaluator.py", line 72, in estimate_observables
    expectation_values = estimator_job.result().values
  File "/home/abhishek/envs/gpuqiskit/lib/python3.10/site-packages/qiskit/primitives/primitive_job.py", line 51, in result
    return self._future.result()
  File "/usr/lib/python3.10/concurrent/futures/_base.py", line 458, in result
    return self.__get_result()
  File "/usr/lib/python3.10/concurrent/futures/_base.py", line 403, in __get_result
    raise self._exception
  File "/usr/lib/python3.10/concurrent/futures/thread.py", line 58, in run
    result = self.fn(*self.args, **self.kwargs)
  File "/home/abhishek/envs/gpuqiskit/lib/python3.10/site-packages/qiskit_aer/primitives/estimator.py", line 172, in _call
    return self._compute_with_approximation(
  File "/home/abhishek/envs/gpuqiskit/lib/python3.10/site-packages/qiskit_aer/primitives/estimator.py", line 455, in _compute_with_approximation
    circuit.save_expectation_value(observable, self._layouts[i])
  File "/home/abhishek/envs/gpuqiskit/lib/python3.10/site-packages/qiskit_aer/library/save_instructions/save_expectation_value.py", line 208, in save_expectation_value
    instr = SaveExpectationValue(
  File "/home/abhishek/envs/gpuqiskit/lib/python3.10/site-packages/qiskit_aer/library/save_instructions/save_expectation_value.py", line 67, in __init__
    raise ValueError("Input operator is not Hermitian.")
ValueError: Input operator is not Hermitian.

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/abhishek/sim/UCCSD/densities/H2/input.py", line 103, in <module>
    uccsd_result = vqe_solver.compute_minimum_eigenvalue(qubit_op, aux_qubit_op)
  File "/home/abhishek/envs/gpuqiskit/lib/python3.10/site-packages/qiskit_algorithms/minimum_eigensolvers/vqe.py", line 214, in compute_minimum_eigenvalue
    aux_operators_evaluated = estimate_observables(
  File "/home/abhishek/envs/gpuqiskit/lib/python3.10/site-packages/qiskit_algorithms/observables_evaluator.py", line 74, in estimate_observables
    raise AlgorithmError("The primitive job failed!") from exc
qiskit_algorithms.exceptions.AlgorithmError: 'The primitive job failed!'

How can we reproduce the issue?

import qiskit_nature
from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.formats.molecule_info import MoleculeInfo
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.drivers.electronic_structure_driver import MethodType
from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit_nature.second_q.algorithms import GroundStateEigensolver
from qiskit_nature.second_q.circuit.library import UCCSD, HartreeFock
from qiskit_nature.second_q.properties import ElectronicDensity
from qiskit_algorithms.optimizers import SLSQP
from qiskit_aer.primitives import Estimator
from qiskit_algorithms.minimum_eigensolvers import VQE
import json
from time import time
from pyscf.scf import UHF
from pyscf.fci import FCI
from pyscf.gto import mole, loads
import numpy as np
BOHR = DistanceUnit.BOHR
qiskit_nature.settings.use_pauli_sum_op = False

mol = pyscf.M(atom="H 0 0 0; H 0 0 0.74")
mol.build()
mf = UHF(mol)
mf.kernel()
start = time()
fci = FCI(mf)
fci.kernel()
end = time()
print(mf.e_tot, fci.e_tot, mf.converged, fci.converged)
print("time taken for pyscf fci calculation ", end - start, "seconds")

molecule = MoleculeInfo(
    symbols=mol.elements,
    coords=mol.atom_coords(),
    multiplicity=mol.spin + 1,
    charge=mol.charge,
    units=BOHR,
)

driver = PySCFDriver.from_molecule(molecule, basis="sto3g", method=MethodType.UHF)
es_problem = driver.run()
density = ElectronicDensity.identity(es_problem.num_spatial_orbitals)
#density = ElectronicDensity.from_orbital_occupation(
#    es_problem.orbital_occupations, es_problem.orbital_occupations_b
#)
es_problem.properties.electronic_density = density

mapper = JordanWignerMapper()
qubit_op = mapper.map(es_problem.second_q_ops()[0])
aux_qubit_op = mapper.map(es_problem.second_q_ops()[1])

# Set up the mapper and solver
initial_state = HartreeFock(
    es_problem.num_spatial_orbitals,
    es_problem.num_particles,
    mapper,
)
ansatz = UCCSD(
    es_problem.num_spatial_orbitals,
    es_problem.num_particles,
    mapper,
    initial_state=initial_state,
)
optimizer = SLSQP(maxiter=100, disp=True)
estimator = Estimator(
    approximation=True,
    run_options={"shots": None},
    backend_options={
        "method": "statevector",
        # "device": "GPU",
        "max_parallel_threads": 24,
    },
)

# Create the VQE solver
vqe_solver = VQE(
    estimator,
    ansatz,
    optimizer,
    initial_point=[0] * ansatz.num_parameters,
)

uccsd_result = vqe_solver.compute_minimum_eigenvalue(qubit_op, aux_qubit_op)
result = es_problem.interpret(uccsd_result)

# Extract the 1-RDM for alpha and beta spins
rdm1_a = result.electronic_density.alpha
rdm1_b = result.electronic_density.beta

print("1-RDM for alpha spin:")
print(rdm1_a)
print("\n1-RDM for beta spin:")
print(rdm1_b)

What should happen?

It should be possible to compute the RDMs computed by using each fermionic Hamiltonian term, transforming them and computing the elements one-by-one...

Any suggestions?

No response