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