# How to perform below operation in Qiskit?

I want to implement the below equation in Qiskit.

$$(A \otimes B).\rho.(B^\dagger \otimes A^\dagger)$$ where $$\rho$$ is a density matrix and $$A$$ and $$B$$ are CNOT gates. $$A=\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \end{bmatrix}$$ $$B=\begin{bmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix}.$$ The final $$pf$$ I got using the codes below is an operator, not a density matrix, as $$dot()$$ returns an operator. Hence, I can not use $$\texttt{state_fidelity()}$$ on it. The dimension of $$\rho$$ is 16. I can not use $$evolve()$$ since it is not any normal operator usage. Please let me know if there is any other method I can implement the above operation. I will be grateful.

p1=[[1,0,0,0],[0,0,0,1],[0,0,1,0],[0,1,0,0]]
p2=[[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]]
pa1=qi.Operator(p1)
pa2=qi.Operator(p2)
pa3=(pa2.transpose()).tensor(pa1.transpose())
pf=((pa1.tensor(pa2)).dot(rho)).dot(pa3)


where qi is :

import qiskit.quantum_info as qi


The constructor of Qiskit's DensityMatrix class accepts any object with a to_matrix() method. And since Operator class has a to_matrix(), you can simply convert pf to a DensityMatrix that can be used with state_fidelity():
F = state_fidelity(qi.DensityMatrix(pf), rho2)