def compute_rdms(self, U: QCircuit=None, variables: Variables=None, spin_free: bool=True, get_rdm1: bool=True, get_rdm2: bool=True, ordering='dirac'):
    "\n        Computes the one- and two-particle reduced density matrices (rdm1 and rdm2) given\n        a unitary U. This method uses the standard ordering in physics as denoted below.\n        Note, that the representation of the density matrices depends on the qubit transformation\n        used. The Jordan-Wigner encoding corresponds to 'classical' second quantized density\n        matrices in the occupation picture.\n\n        We only consider real orbitals and thus real-valued RDMs.\n        The matrices are set as private members _rdm1, _rdm2 and can be accessed via the properties rdm1, rdm2.\n\n        .. math :\n            \\text{rdm1: } \\gamma^p_q = \\langle \\psi | a^p a_q | \\psi \\rangle\n                                     = \\langle U 0 | a^p a_q | U 0 \\rangle\n            \\text{rdm2: } \\gamma^{pq}_{rs} = \\langle \\psi | a^p a^q a_s a_r | \\psi \\rangle\n                                             = \\langle U 0 | a^p a^q a_s a_r | U 0 \\rangle\n\n        Parameters\n        ----------\n        U :\n            Quantum Circuit to achieve the desired state \\psi = U |0\\rangle, non-optional\n        variables :\n            If U is parametrized, then need to hand over a set of fixed variables\n        spin_free :\n            Set whether matrices should be spin-free (summation over spin) or defined by spin-orbitals\n        get_rdm1, get_rdm2 :\n            Set whether either one or both rdm1, rdm2 should be computed. If both are needed at some point,\n            it is recommended to compute them at once.\n\n        Returns\n        -------\n        "
    if (U is None):
        raise TequilaException('Need to specify a Quantum Circuit.')
    if (type(self.transformation).__name__ == 'BravyiKitaevFast'):
        raise TequilaException('The Bravyi-Kitaev-Superfast transformation does not support general FermionOperators yet.')
    n_SOs = (2 * self.n_orbitals)
    n_MOs = self.n_orbitals
    if (U is None):
        raise TequilaException('Need to specify a Quantum Circuit.')

    def _get_of_op(operator_tuple):
        ' Returns operator given by a operator tuple as OpenFermion - Fermion operator '
        op = openfermion.FermionOperator(operator_tuple)
        return op

    def _get_qop_hermitian(of_operator) -> QubitHamiltonian:
        ' Returns Hermitian part of Fermion operator as QubitHamiltonian '
        qop = self.transformation(of_operator)
        (real, imag) = qop.split(hermitian=True)
        if real:
            return real
        elif (not real):
            raise TequilaException('Qubit Hamiltonian does not have a Hermitian part. Operator ={}'.format(of_operator))

    def _build_1bdy_operators_spinful() -> list:
        ' Returns spinful one-body operators as a symmetry-reduced list of QubitHamiltonians '
        ops = []
        for p in range(n_SOs):
            for q in range((p + 1)):
                op_tuple = ((p, 1), (q, 0))
                op = _get_of_op(op_tuple)
                ops += [op]
        return ops

    def _build_2bdy_operators_spinful() -> list:
        ' Returns spinful two-body operators as a symmetry-reduced list of QubitHamiltonians '
        ops = []
        for p in range(n_SOs):
            for q in range(p):
                for r in range(n_SOs):
                    for s in range(r):
                        if (((p * n_SOs) + q) >= ((r * n_SOs) + s)):
                            op_tuple = ((p, 1), (q, 1), (s, 0), (r, 0))
                            op = _get_of_op(op_tuple)
                            ops += [op]
        return ops

    def _build_1bdy_operators_spinfree() -> list:
        ' Returns spinfree one-body operators as a symmetry-reduced list of QubitHamiltonians '
        ops = []
        for p in range(n_MOs):
            for q in range((p + 1)):
                op_tuple = (((2 * p), 1), ((2 * q), 0))
                op = _get_of_op(op_tuple)
                op_tuple = ((((2 * p) + 1), 1), (((2 * q) + 1), 0))
                op += _get_of_op(op_tuple)
                ops += [op]
        return ops

    def _build_2bdy_operators_spinfree() -> list:
        ' Returns spinfree two-body operators as a symmetry-reduced list of QubitHamiltonians '
        ops = []
        for (p, q, r, s) in product(range(n_MOs), repeat=4):
            if ((((p * n_MOs) + q) >= ((r * n_MOs) + s)) and ((p >= q) or (r >= s))):
                op_tuple = ((((2 * p), 1), ((2 * q), 1), ((2 * s), 0), ((2 * r), 0)) if ((p != q) and (r != s)) else '0.0 []')
                op = _get_of_op(op_tuple)
                op_tuple = ((((2 * p), 1), (((2 * q) + 1), 1), (((2 * s) + 1), 0), ((2 * r), 0)) if (((2 * p) != ((2 * q) + 1)) and ((2 * r) != ((2 * s) + 1))) else '0.0 []')
                op += _get_of_op(op_tuple)
                op_tuple = (((((2 * p) + 1), 1), ((2 * q), 1), ((2 * s), 0), (((2 * r) + 1), 0)) if ((((2 * p) + 1) != (2 * q)) and (((2 * r) + 1) != (2 * s))) else '0.0 []')
                op += _get_of_op(op_tuple)
                op_tuple = (((((2 * p) + 1), 1), (((2 * q) + 1), 1), (((2 * s) + 1), 0), (((2 * r) + 1), 0)) if ((p != q) and (r != s)) else '0.0 []')
                op += _get_of_op(op_tuple)
                ops += [op]
        return ops

    def _assemble_rdm1(evals) -> numpy.ndarray:
        '\n            Returns spin-ful or spin-free one-particle RDM built by symmetry conditions\n            Same symmetry with or without spin, so we can use the same function\n            '
        N = (n_MOs if spin_free else n_SOs)
        rdm1 = numpy.zeros([N, N])
        ctr: int = 0
        for p in range(N):
            for q in range((p + 1)):
                rdm1[(p, q)] = evals[ctr]
                rdm1[(q, p)] = rdm1[(p, q)]
                ctr += 1
        return rdm1

    def _assemble_rdm2_spinful(evals) -> numpy.ndarray:
        ' Returns spin-ful two-particle RDM built by symmetry conditions '
        ctr: int = 0
        rdm2 = numpy.zeros([n_SOs, n_SOs, n_SOs, n_SOs])
        for p in range(n_SOs):
            for q in range(p):
                for r in range(n_SOs):
                    for s in range(r):
                        if (((p * n_SOs) + q) >= ((r * n_SOs) + s)):
                            rdm2[(p, q, r, s)] = evals[ctr]
                            rdm2[(r, s, p, q)] = rdm2[(p, q, r, s)]
                            ctr += 1
        for p in range(n_SOs):
            for q in range(p):
                for r in range(n_SOs):
                    for s in range(r):
                        rdm2[(p, q, s, r)] = ((- 1) * rdm2[(p, q, r, s)])
                        rdm2[(q, p, r, s)] = ((- 1) * rdm2[(p, q, r, s)])
                        rdm2[(q, p, s, r)] = rdm2[(p, q, r, s)]
        return rdm2

    def _assemble_rdm2_spinfree(evals) -> numpy.ndarray:
        ' Returns spin-free two-particle RDM built by symmetry conditions '
        ctr: int = 0
        rdm2 = numpy.zeros([n_MOs, n_MOs, n_MOs, n_MOs])
        for (p, q, r, s) in product(range(n_MOs), repeat=4):
            if ((((p * n_MOs) + q) >= ((r * n_MOs) + s)) and ((p >= q) or (r >= s))):
                rdm2[(p, q, r, s)] = evals[ctr]
                rdm2[(r, s, p, q)] = rdm2[(p, q, r, s)]
                ctr += 1
        for (p, q, r, s) in product(range(n_MOs), repeat=4):
            if ((p >= q) or (r >= s)):
                rdm2[(q, p, s, r)] = rdm2[(p, q, r, s)]
        return rdm2
    qops = []
    if spin_free:
        qops += (_build_1bdy_operators_spinfree() if get_rdm1 else [])
        qops += (_build_2bdy_operators_spinfree() if get_rdm2 else [])
    else:
        qops += (_build_1bdy_operators_spinful() if get_rdm1 else [])
        qops += (_build_2bdy_operators_spinful() if get_rdm2 else [])
    qops = [_get_qop_hermitian(op) for op in qops]
    evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)

    def _reset_rdm(rdm):
        if (rdm is not None):
            if (spin_free and (rdm.shape[0] != n_MOs)):
                return None
            if ((not spin_free) and (rdm.shape[0] != n_SOs)):
                return None
        return rdm
    self._rdm1 = _reset_rdm(self._rdm1)
    self._rdm2 = _reset_rdm(self._rdm2)
    if get_rdm1:
        len_1 = (((n_MOs * (n_MOs + 1)) // 2) if spin_free else ((n_SOs * (n_SOs + 1)) // 2))
    else:
        len_1 = 0
    (evals_1, evals_2) = (evals[:len_1], evals[len_1:])
    self._rdm1 = (_assemble_rdm1(evals_1) if get_rdm1 else self._rdm1)
    if spin_free:
        self._rdm2 = (_assemble_rdm2_spinfree(evals_2) if get_rdm2 else self._rdm2)
    else:
        self._rdm2 = (_assemble_rdm2_spinful(evals_2) if get_rdm2 else self._rdm2)
    if get_rdm2:
        rdm2 = NBodyTensor(elems=self.rdm2, ordering='dirac')
        rdm2.reorder(to=ordering)
        rdm2 = rdm2.elems
        self._rdm2 = rdm2
    if get_rdm1:
        if get_rdm2:
            return (self.rdm1, self.rdm2)
        else:
            return self.rdm1
    elif get_rdm2:
        return self.rdm2
    else:
        warnings.warn('compute_rdms called with instruction to not compute?', TequilaWarning)
Computes the one- and two-particle reduced density matrices (rdm1 and rdm2) given
a unitary U. This method uses the standard ordering in physics as denoted below.
Note, that the representation of the density matrices depends on the qubit transformation
used. The Jordan-Wigner encoding corresponds to 'classical' second quantized density
matrices in the occupation picture.

We only consider real orbitals and thus real-valued RDMs.
The matrices are set as private members _rdm1, _rdm2 and can be accessed via the properties rdm1, rdm2.

.. math :
    \text{rdm1: } \gamma^p_q = \langle \psi | a^p a_q | \psi \rangle
                             = \langle U 0 | a^p a_q | U 0 \rangle
    \text{rdm2: } \gamma^{pq}_{rs} = \langle \psi | a^p a^q a_s a_r | \psi \rangle
                                     = \langle U 0 | a^p a^q a_s a_r | U 0 \rangle

Parameters
----------
U :
    Quantum Circuit to achieve the desired state \psi = U |0\rangle, non-optional
variables :
    If U is parametrized, then need to hand over a set of fixed variables
spin_free :
    Set whether matrices should be spin-free (summation over spin) or defined by spin-orbitals
get_rdm1, get_rdm2 :
    Set whether either one or both rdm1, rdm2 should be computed. If both are needed at some point,
    it is recommended to compute them at once.

Returns
-------
def compute_rdms(self, U: QCircuit=None, variables: Variables=None, spin_free: bool=True, get_rdm1: bool=True, get_rdm2: bool=True, ordering='dirac'):
    "\n        Computes the one- and two-particle reduced density matrices (rdm1 and rdm2) given\n        a unitary U. This method uses the standard ordering in physics as denoted below.\n        Note, that the representation of the density matrices depends on the qubit transformation\n        used. The Jordan-Wigner encoding corresponds to 'classical' second quantized density\n        matrices in the occupation picture.\n\n        We only consider real orbitals and thus real-valued RDMs.\n        The matrices are set as private members _rdm1, _rdm2 and can be accessed via the properties rdm1, rdm2.\n\n        .. math :\n            \\text{rdm1: } \\gamma^p_q = \\langle \\psi | a^p a_q | \\psi \\rangle\n                                     = \\langle U 0 | a^p a_q | U 0 \\rangle\n            \\text{rdm2: } \\gamma^{pq}_{rs} = \\langle \\psi | a^p a^q a_s a_r | \\psi \\rangle\n                                             = \\langle U 0 | a^p a^q a_s a_r | U 0 \\rangle\n\n        Parameters\n        ----------\n        U :\n            Quantum Circuit to achieve the desired state \\psi = U |0\\rangle, non-optional\n        variables :\n            If U is parametrized, then need to hand over a set of fixed variables\n        spin_free :\n            Set whether matrices should be spin-free (summation over spin) or defined by spin-orbitals\n        get_rdm1, get_rdm2 :\n            Set whether either one or both rdm1, rdm2 should be computed. If both are needed at some point,\n            it is recommended to compute them at once.\n\n        Returns\n        -------\n        "
    if (U is None):
        raise TequilaException('Need to specify a Quantum Circuit.')
    if (type(self.transformation).__name__ == 'BravyiKitaevFast'):
        raise TequilaException('The Bravyi-Kitaev-Superfast transformation does not support general FermionOperators yet.')
    n_SOs = (2 * self.n_orbitals)
    n_MOs = self.n_orbitals
    if (U is None):
        raise TequilaException('Need to specify a Quantum Circuit.')

    def _get_of_op(operator_tuple):
        ' Returns operator given by a operator tuple as OpenFermion - Fermion operator '
        op = openfermion.FermionOperator(operator_tuple)
        return op

    def _get_qop_hermitian(of_operator) -> QubitHamiltonian:
        ' Returns Hermitian part of Fermion operator as QubitHamiltonian '
        qop = self.transformation(of_operator)
        (real, imag) = qop.split(hermitian=True)
        if real:
            return real
        elif (not real):
            raise TequilaException('Qubit Hamiltonian does not have a Hermitian part. Operator ={}'.format(of_operator))

    def _build_1bdy_operators_spinful() -> list:
        ' Returns spinful one-body operators as a symmetry-reduced list of QubitHamiltonians '
        ops = []
        for p in range(n_SOs):
            for q in range((p + 1)):
                op_tuple = ((p, 1), (q, 0))
                op = _get_of_op(op_tuple)
                ops += [op]
        return ops

    def _build_2bdy_operators_spinful() -> list:
        ' Returns spinful two-body operators as a symmetry-reduced list of QubitHamiltonians '
        ops = []
        for p in range(n_SOs):
            for q in range(p):
                for r in range(n_SOs):
                    for s in range(r):
                        if (((p * n_SOs) + q) >= ((r * n_SOs) + s)):
                            op_tuple = ((p, 1), (q, 1), (s, 0), (r, 0))
                            op = _get_of_op(op_tuple)
                            ops += [op]
        return ops

    def _build_1bdy_operators_spinfree() -> list:
        ' Returns spinfree one-body operators as a symmetry-reduced list of QubitHamiltonians '
        ops = []
        for p in range(n_MOs):
            for q in range((p + 1)):
                op_tuple = (((2 * p), 1), ((2 * q), 0))
                op = _get_of_op(op_tuple)
                op_tuple = ((((2 * p) + 1), 1), (((2 * q) + 1), 0))
                op += _get_of_op(op_tuple)
                ops += [op]
        return ops

    def _build_2bdy_operators_spinfree() -> list:
        ' Returns spinfree two-body operators as a symmetry-reduced list of QubitHamiltonians '
        ops = []
        for (p, q, r, s) in product(range(n_MOs), repeat=4):
            if ((((p * n_MOs) + q) >= ((r * n_MOs) + s)) and ((p >= q) or (r >= s))):
                op_tuple = ((((2 * p), 1), ((2 * q), 1), ((2 * s), 0), ((2 * r), 0)) if ((p != q) and (r != s)) else '0.0 []')
                op = _get_of_op(op_tuple)
                op_tuple = ((((2 * p), 1), (((2 * q) + 1), 1), (((2 * s) + 1), 0), ((2 * r), 0)) if (((2 * p) != ((2 * q) + 1)) and ((2 * r) != ((2 * s) + 1))) else '0.0 []')
                op += _get_of_op(op_tuple)
                op_tuple = (((((2 * p) + 1), 1), ((2 * q), 1), ((2 * s), 0), (((2 * r) + 1), 0)) if ((((2 * p) + 1) != (2 * q)) and (((2 * r) + 1) != (2 * s))) else '0.0 []')
                op += _get_of_op(op_tuple)
                op_tuple = (((((2 * p) + 1), 1), (((2 * q) + 1), 1), (((2 * s) + 1), 0), (((2 * r) + 1), 0)) if ((p != q) and (r != s)) else '0.0 []')
                op += _get_of_op(op_tuple)
                ops += [op]
        return ops

    def _assemble_rdm1(evals) -> numpy.ndarray:
        '\n            Returns spin-ful or spin-free one-particle RDM built by symmetry conditions\n            Same symmetry with or without spin, so we can use the same function\n            '
        N = (n_MOs if spin_free else n_SOs)
        rdm1 = numpy.zeros([N, N])
        ctr: int = 0
        for p in range(N):
            for q in range((p + 1)):
                rdm1[(p, q)] = evals[ctr]
                rdm1[(q, p)] = rdm1[(p, q)]
                ctr += 1
        return rdm1

    def _assemble_rdm2_spinful(evals) -> numpy.ndarray:
        ' Returns spin-ful two-particle RDM built by symmetry conditions '
        ctr: int = 0
        rdm2 = numpy.zeros([n_SOs, n_SOs, n_SOs, n_SOs])
        for p in range(n_SOs):
            for q in range(p):
                for r in range(n_SOs):
                    for s in range(r):
                        if (((p * n_SOs) + q) >= ((r * n_SOs) + s)):
                            rdm2[(p, q, r, s)] = evals[ctr]
                            rdm2[(r, s, p, q)] = rdm2[(p, q, r, s)]
                            ctr += 1
        for p in range(n_SOs):
            for q in range(p):
                for r in range(n_SOs):
                    for s in range(r):
                        rdm2[(p, q, s, r)] = ((- 1) * rdm2[(p, q, r, s)])
                        rdm2[(q, p, r, s)] = ((- 1) * rdm2[(p, q, r, s)])
                        rdm2[(q, p, s, r)] = rdm2[(p, q, r, s)]
        return rdm2

    def _assemble_rdm2_spinfree(evals) -> numpy.ndarray:
        ' Returns spin-free two-particle RDM built by symmetry conditions '
        ctr: int = 0
        rdm2 = numpy.zeros([n_MOs, n_MOs, n_MOs, n_MOs])
        for (p, q, r, s) in product(range(n_MOs), repeat=4):
            if ((((p * n_MOs) + q) >= ((r * n_MOs) + s)) and ((p >= q) or (r >= s))):
                rdm2[(p, q, r, s)] = evals[ctr]
                rdm2[(r, s, p, q)] = rdm2[(p, q, r, s)]
                ctr += 1
        for (p, q, r, s) in product(range(n_MOs), repeat=4):
            if ((p >= q) or (r >= s)):
                rdm2[(q, p, s, r)] = rdm2[(p, q, r, s)]
        return rdm2
    qops = []
    if spin_free:
        qops += (_build_1bdy_operators_spinfree() if get_rdm1 else [])
        qops += (_build_2bdy_operators_spinfree() if get_rdm2 else [])
    else:
        qops += (_build_1bdy_operators_spinful() if get_rdm1 else [])
        qops += (_build_2bdy_operators_spinful() if get_rdm2 else [])
    qops = [_get_qop_hermitian(op) for op in qops]
    evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)

    def _reset_rdm(rdm):
        if (rdm is not None):
            if (spin_free and (rdm.shape[0] != n_MOs)):
                return None
            if ((not spin_free) and (rdm.shape[0] != n_SOs)):
                return None
        return rdm
    self._rdm1 = _reset_rdm(self._rdm1)
    self._rdm2 = _reset_rdm(self._rdm2)
    if get_rdm1:
        len_1 = (((n_MOs * (n_MOs + 1)) // 2) if spin_free else ((n_SOs * (n_SOs + 1)) // 2))
    else:
        len_1 = 0
    (evals_1, evals_2) = (evals[:len_1], evals[len_1:])
    self._rdm1 = (_assemble_rdm1(evals_1) if get_rdm1 else self._rdm1)
    if spin_free:
        self._rdm2 = (_assemble_rdm2_spinfree(evals_2) if get_rdm2 else self._rdm2)
    else:
        self._rdm2 = (_assemble_rdm2_spinful(evals_2) if get_rdm2 else self._rdm2)
    if get_rdm2:
        rdm2 = NBodyTensor(elems=self.rdm2, ordering='dirac')
        rdm2.reorder(to=ordering)
        rdm2 = rdm2.elems
        self._rdm2 = rdm2
    if get_rdm1:
        if get_rdm2:
            return (self.rdm1, self.rdm2)
        else:
            return self.rdm1
    elif get_rdm2:
        return self.rdm2
    else:
        warnings.warn('compute_rdms called with instruction to not compute?', TequilaWarning)