How to rewrite this Cirq problem into Qiskit?

Question

I have a following example problem in Cirq, representing a very simple operator given by

$$ \hat{H} = a^\dagger_0a_2 + a_0a^\dagger_1 + 0.9\, a^\dagger_0a_1a^\dagger_2a^\dagger_3 + a^\dagger_0a_1a_2a_3 $$

and an example ansatz, obtaining the expectation value in the end.

import cirq
import numpy as np
import openfermion as of
import sympy

H = of.FermionOperator('0^ 2', 1) + of.FermionOperator('0 1^', 1) + of.FermionOperator('0^ 1 2^ 3^', 0.9) \
    + of.FermionOperator('0^ 1 2 3', -0.1)

H = of.FermionOperator('0^ 1', 1) + of.FermionOperator('0 2^', 1) + of.FermionOperator('0^ 2 1^ 3^', 0.9) \
    + of.FermionOperator('0^ 2 1 3', -0.1)

Hmat = of.get_sparse_operator(H, n_qubits=4).A

print(H)
print(sorted(np.linalg.eig(Hmat)[0]))

circuit = cirq.Circuit()
qubits = cirq.LineQubit.range(4)

alpha = sympy.Symbol('alpha')
beta = sympy.Symbol('beta')
gamma = sympy.Symbol('gamma')
delta = sympy.Symbol('delta')

circuit.append(cirq.ry(alpha).on(qubits[0]))
circuit.append(cirq.ry(beta).on(qubits[1]))
circuit.append(cirq.ry(gamma).on(qubits[2]))
circuit.append(cirq.ry(delta).on(qubits[3]))
circuit.append(cirq.X(qubits[2]).controlled_by(qubits[0]))

full_circ = cirq.resolve_parameters(circuit, {'alpha': 1, 'beta': 1, 'gamma': 1, 'delta': 1})
sv = full_circ.final_state_vector()

print(sv)
print(sv.conj().T @ Hmat @ sv)

This Cirq code is giving the following output:

1.0 [0 2^] +
1 [0^ 1] +
-0.1 [0^ 2 1 3] +
0.9 [0^ 2 1^ 3^]
[(-0.547722557505166+5.551115123125783e-17j), (-7.994927024450326e-17-0.5477225575051663j), 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j, (4.85722573273506e-17+0.5477225575051662j), (0.5477225575051663+0j)]
[0.5931328 +0.j 0.32402992+0.j 0.32402992+0.j 0.17701835+0.j
 0.32402992+0.j 0.17701835+0.j 0.17701835+0.j 0.09670557+0.j
 0.17701835+0.j 0.09670557+0.j 0.32402992+0.j 0.17701835+0.j
 0.09670557+0.j 0.05283049+0.j 0.17701835+0.j 0.09670557+0.j]
(0.11040049624731614+0j)

I'd like to be able to rewrite this example into an equivalent code in Qiskit, where, as I'm aware is a different spin-orbital mapping, so that I switched 1st and 2nd qubits in ansatz. The code looks like this:

from qiskit import QuantumCircuit
from qiskit.circuit import ParameterVector
from qiskit.quantum_info import Statevector
from qiskit_nature.operators.second_quantization import FermionicOp
import numpy as np

H = FermionicOp([('+_0 -_2', 1), ('-_0 +_1', 1), ('+_0 -_1 +_2 +_3', 0.9), ('+_0 -_1 -_2 -_3', -0.1)], register_length=4)


print(H)
print(sorted(np.linalg.eig(H.to_matrix().A)[0]))

ansatz = QuantumCircuit(4)
p = ParameterVector('theta', length=4)
for i, q in enumerate(ansatz.qubits):
    ansatz.ry(p[i], q)
ansatz.cx(0, 1)

params = [1]*4

vecA = Statevector(ansatz.bind_parameters(params)).data

print(vecA)
print((vecA.conj().T @ H.to_matrix() @ vecA).real)

But the output is different:

Fermionic Operator
register length=4, number terms=4
  (1+0j) * ( +_0 -_2 )
+ (1+0j) * ( -_0 +_1 )
+ (0.9+0j) * ( +_0 -_1 +_2 +_3 )
+ (-0.1+0j) * ( +_0 -_1 -_2 -_3 )
[(-0.547722557505166+5.551115123125783e-17j), (-7.997636868137208e-17-0.5477225575051663j), 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j, (4.85722573273506e-17+0.5477225575051662j), (0.5477225575051663+0j)]
[0.5931328 +0.j 0.17701835+0.j 0.32402992+0.j 0.32402992+0.j
 0.32402992+0.j 0.09670557+0.j 0.17701835+0.j 0.17701835+0.j
 0.32402992+0.j 0.09670557+0.j 0.17701835+0.j 0.17701835+0.j
 0.17701835+0.j 0.05283049+0.j 0.09670557+0.j 0.09670557+0.j]
-0.02401568533842704

Why are the expectation values different? Is the operator rewritten in an incorrect way or does the ansatz need different treatment, when rewritten to Qiskit?