The definition of a 2-RDM is $D_{ijkl}=\langle\psi|a_i^{\dagger}a_j^{\dagger}a_k a_l|\psi\rangle,i,,j,k,l\in[0,N-1]$, where $N$ is the number of qubits/ orbitals in the quantum system.I am trying to use the following code to get the 2-RDM for a given state and the eigenvalue for the 2-RDM,

from qiskit import *
import numpy as np

#Operator Imports
from qiskit.opflow import Z, X, I

#Circuit imports
from qiskit_nature.circuit.library import HartreeFock
from qiskit import Aer
from qiskit_nature.mappers.second_quantization import JordanWignerMapper
from qiskit_nature.converters.second_quantization import QubitConverter
from qiskit.algorithms import VQE, NumPyMinimumEigensolver
import matplotlib.pyplot as plt
import matplotlib
from qiskit.tools.visualization import circuit_drawer
from qiskit_nature.operators.second_quantization import FermionicOp
from qiskit.opflow import Z, X, I, StateFn, CircuitStateFn, SummedOp
from qiskit.opflow.converters import CircuitSampler
from qiskit.utils import QuantumInstance
from qiskit.opflow import expectations
import itertools

mapper = JordanWignerMapper()
converter = QubitConverter(mapper=mapper, two_qubit_reduction=False)
num_particles = (2, 2)
num_spin_orbitals = 6
init_state = HartreeFock(num_spin_orbitals, num_particles, converter)

init_state1 = QuantumCircuit(num_spin_orbitals)


def get_rdm(init_state):
    two_rdm = np.zeros((num_spin_orbitals,) * 4)
    for i, j, k, l in itertools.product(range(num_spin_orbitals), repeat=4):
        # if i != j and k != l:
        s = "+_{i} +_{j} -_{k} -_{l}".format(i=str(i),j=str(j),k=str(k),l=str(l))
        fermi_term = FermionicOp(s, register_length=num_spin_orbitals)
        qubit_term = converter.convert(fermi_term, num_particles=num_particles)

        # Evaluate the Hamiltonian term w.r.t. the given state
        # temp = (~init_state @ qubit_term @ init_state).eval()
        temp = ~StateFn(qubit_term) @ CircuitStateFn(primitive=init_state, coeff=1.)
        temp = int(temp.eval())

        two_rdm[i,j,k,l] = temp

    return two_rdm

# rdm1 = rdm(init_state)
# print('--------------')
two_rdm = get_rdm(init_state1)

# Transpose of the 2-RDM
two_rdm = np.transpose(two_rdm, (0, 1, 3, 2))
two_rdm = two_rdm.reshape((num_spin_orbitals ** 2, num_spin_orbitals ** 2))
# SVD of the 2-RDM
a, b = np.linalg.eig(two_rdm)

u, s, vh = np.linalg.svd(two_rdm, full_matrices=True)

For the Hartree Fock state, the eigenvalues of the corresponding 2-RDM output six non-zero elements, which all are $2$. In my opinion, the eigenvalues of the 2-RDM stand for the meaning that probability of two electron residing in the corresponding two orbitals. Therefore, the eigenvalue should be from zero to one.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.