Qiskit noisy FakeManila simulator + ZNE

Question

I am trying to run a simple VQE calculation on H2 using FakeManila noise model and ZNE for error mitigation; see code below. The problem is that I can never reach not even close the ground electronic state energy (~-1.8551550878043606 Eh) at this geometry (as obtained at FCI/sto-3g level and also with ibm_qasm_simulator). I was wondering whether someboby more experience could help me in spotting any possible problem.

import numpy as np

from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.drivers import PySCFDriver

import qiskit_nature
from qiskit_nature.second_q.circuit.library import HartreeFock, UCCSD
from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit.algorithms.optimizers import COBYLA

from qiskit.providers.fake_provider import FakeManila
from qiskit_aer.noise import NoiseModel

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    Estimator,
    Options,
    Session
)

# Avoid using the deprecated `PauliSumOp` object
qiskit_nature.settings.use_pauli_sum_op = False

# using qiskit runtime service
service = QiskitRuntimeService()

# run on simulator
backend = service.backend("ibmq_qasm_simulator")

# Import a noise model
fake_backend = FakeManila()
noise_model = NoiseModel.from_backend(fake_backend)

options = Options()
# simulator options
options.simulator = {
    "noise_model": noise_model,
    "seed_simulator": 42
}
# error supression options
options.optimization_level = 3
# error mitigation options
options.resilience_level = 2  # ZNE
# zne options
# options.resilience.noise_amplifier = 'CxAmplifier'
# options.resilience.noise_factors = tuple(range(1, 9, 2))
# options.resilience.extrapolator = 'LinearExtrapolator'
# execution options
options.execution.shots = 6000

# Use estimator to get the expected values
estimator = Estimator(backend=backend, options=options)

# Calculate qubit hamiltonian
driver = PySCFDriver(
    atom="H 0 0 0; H 0 0 0.737166",
    basis="sto3g",
    charge=0,
    spin=0,
    unit=DistanceUnit.ANGSTROM,
)

problem = driver.run()
print(f"Reference energy: {problem.reference_energy}")

nuc_rep_energy = problem.nuclear_repulsion_energy
print(f"Nuclear repulsion energy: {nuc_rep_energy}")

hamiltonian = problem.hamiltonian
second_q_op = hamiltonian.second_q_op()

mapper = JordanWignerMapper()
qubit_op = mapper.map(second_q_op)

# Set up the variational form/ansatz
n_active_electrons = (1, 1)  # => (n_alpha, n_beta)
n_active_spatial_orbitals = 2

reference_state = HartreeFock(
    n_active_spatial_orbitals,
    n_active_electrons,
    mapper,
)

# print(reference_state.draw())

ansatz = UCCSD(
    n_active_spatial_orbitals,
    n_active_electrons,
    mapper,
    initial_state=reference_state
)

# print(ansatz.decompose().draw())


def cost_func(params):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters

    Returns:
        float: Energy estimate
    """
    job = estimator.run(ansatz, qubit_op, parameter_values=params)
    result = job.result()
    energy = result.values[0]
    print("=== COST FUNCTION SUMMARY ===")
    print(f">>> Job ID: {job.job_id()}")
    print(f">>> Job Status: {job.status()}")
    print("=============================")
    print(f">>> Job Input: {job.inputs}")
    print("=============================")
    print(f">>> Backend: {job.backend()}")
    print(f">>> {result}")
    print("=============================")
    print(f">>> Expectation value (Hartree): {energy}")
    print(f">>> Total ground state energy (Hartree): {energy+nuc_rep_energy}")
    print("=============================\n")
    return energy


with Session(service=service, backend=backend) as session:

    initial_theta = np.array([1.57079357, 1.57087253, 1.45852109])
    # cost_func(initial_theta)

    optimizer = COBYLA()
    res = optimizer.minimize(
        cost_func,
        x0=initial_theta
    )
    print(res)

    session.close()

```