Repository: tequilahub/tequila (language: python)

Issue:
Pr master sync



Discussion / hints:

Fix patch:
diff --git a/src/tequila/circuit/gates.py b/src/tequila/circuit/gates.py
index 3662e607..f39a9cb0 100644
--- a/src/tequila/circuit/gates.py
+++ b/src/tequila/circuit/gates.py
@@ -940,7 +940,7 @@ def QubitExcitation(angle: typing.Union[numbers.Real, Variable, typing.Hashable]
     except:
         raise Exception("QubitExcitation: Needs an even number of targets")
 
-    return QCircuit.wrap_gate(QubitExcitationImpl(angle=angle, target=target, assume_real=assume_real, compile_options=compile_options))
+    return QCircuit.wrap_gate(QubitExcitationImpl(angle=angle, target=target, assume_real=assume_real, compile_options=compile_options, control=control))
 
 
 """
diff --git a/src/tequila/circuit/qpic.py b/src/tequila/circuit/qpic.py
index 18faaffa..be2e394e 100644
--- a/src/tequila/circuit/qpic.py
+++ b/src/tequila/circuit/qpic.py
@@ -224,7 +224,7 @@ def export_to(circuit,
             'always_use_generators': True,
             'group_together': "BARRIER"
         }
-    elif not hasattr("style", "items"):
+    elif not hasattr(style, "items"):
         raise Exception(
             "style needs to be `tequila`, or `standard` or `generators` or a dictionary, you gave: {}".format(
                 str(style)))
diff --git a/src/tequila/grouping/binary_rep.py b/src/tequila/grouping/binary_rep.py
index 6aea305d..ac15efd0 100644
--- a/src/tequila/grouping/binary_rep.py
+++ b/src/tequila/grouping/binary_rep.py
@@ -26,7 +26,7 @@ def init_from_qubit_hamiltonian(cls, hamiltonian: QubitHamiltonian, n_qubits=Non
                 del Hof.terms[()]
                 hamiltonian = QubitHamiltonian.from_openfermion(Hof)
         if n_qubits is None:
-            n_qubits = hamiltonian.n_qubits
+            n_qubits = max(hamiltonian.qubits)+1
         binary_terms = [
             BinaryPauliString(
                 p.binary(n_qubits).binary,
diff --git a/src/tequila/quantumchemistry/chemistry_tools.py b/src/tequila/quantumchemistry/chemistry_tools.py
index 82b72c0d..faa36843 100644
--- a/src/tequila/quantumchemistry/chemistry_tools.py
+++ b/src/tequila/quantumchemistry/chemistry_tools.py
@@ -2,12 +2,12 @@
 import typing
 import warnings
 from dataclasses import dataclass
-
+from copy import deepcopy
+from numbers import Real
 import numpy
 
-from tequila import BitString, QCircuit, TequilaException
+from tequila import BitString, QCircuit, TequilaException,Variable,compile_circuit
 from tequila.circuit import gates
-
 try:
     from openfermion.ops.representations import get_active_space_integrals  # needs openfermion 1.3
 except ImportError as E:
@@ -50,16 +50,132 @@ def __init__(self, generator, p0, transformation, indices=None, *args, **kwargs)
         self._name = "FermionicExcitation"
         self.transformation = transformation
         self.indices = indices
-
+        if not hasattr(indices[0],"__len__"):
+            self.indices = [(indices[2 * i], indices[2 * i+1]) for i in range(len(indices) // 2)]
+        self.sign = self.format_excitation_variables(self.indices)
+        self.indices = self.format_excitation_indices(self.indices)
     def compile(self, *args, **kwargs):
         if self.is_convertable_to_qubit_excitation():
             target = []
             for x in self.indices:
                 for y in x:
                     target.append(y)
-            return gates.QubitExcitation(target=target, angle=-self.parameter, control=self.control)
+            return gates.QubitExcitation(target=target, angle=self.parameter, control=self.control)
+        else:
+            if self.transformation.lower().strip("_") == "jordanwigner":
+                return self.fermionic_excitation(angle=self.sign*self.parameter, indices=self.indices, control=self.control,opt=False)
+            else:
+                return gates.Trotterized(generator=self.generator, control=self.control, angle=self.parameter, steps=1)
+    def format_excitation_indices(self, idx):
+        """
+        Consistent formatting of excitation indices
+        idx = [(p0,q0),(p1,q1),...,(pn,qn)]
+        sorted as: p0<p1<pn and pi<qi
+        :param idx: list of index tuples describing a single(!) fermionic excitation
+        :return: list of index tuples
+        """
+
+        idx = [tuple(sorted(x)) for x in idx]
+        idx = sorted(idx, key=lambda x: x[0])
+        return list(idx)
+    def format_excitation_variables(self, idx):
+        """
+        Consistent formatting of excitation variable
+        idx = [(p0,q0),(p1,q1),...,(pn,qn)]
+        sorted as: pi<qi and p0 < p1 < p2
+        :param idx: list of index tuples describing a single(!) fermionic excitation
+        :return: sign of the variable with re-ordered indices
+        """
+        sig = 1
+        for pair in idx:
+            if pair[1]>pair[0]:
+                sig *= -1
+        for pair in range(len(idx)-1):
+            if idx[pair+1][0]>idx[pair][0]:
+                sig *= -1
+        return sig
+    def cCRy(self, target: int, dcontrol: typing.Union[list, int], control: typing.Union[list, int],
+             angle: typing.Union[Real, Variable, typing.Hashable], case: int = 1) -> QCircuit:
+        '''
+        Compilation of CRy as on https://doi.org/10.1103/PhysRevA.102.062612
+        If not control passed, Ry returned
+        Parameters
+        ----------
+        case: if 1 employs eq. 12 from the paper, if 0 eq. 13
+        '''
+        if control is not None and not len(control):
+            control = None
+        if isinstance(dcontrol, int):
+            dcontrol = [dcontrol]
+        if not len(dcontrol):
+            return compile_circuit(gates.Ry(angle=angle, target=target, control=control))
         else:
-            return gates.Trotterized(generator=self.generator, control=self.control, angle=self.parameter, steps=1)
+            if isinstance(angle, str):
+                angle = Variable(angle)
+            U = QCircuit()
+            aux = dcontrol[0]
+            ctr = deepcopy(dcontrol)
+            ctr.pop(0)
+            if case:
+                U += self.cCRy(target=target, dcontrol=ctr, angle=angle / 2, case=1, control=control) + gates.H(
+                    aux) + gates.CNOT(target, aux)
+                U += self.cCRy(target=target, dcontrol=ctr, angle=-angle / 2, case=0, control=control) + gates.CNOT(
+                    target, aux) + gates.H(aux)
+            else:
+                U += gates.H(aux) + gates.CNOT(target, aux) + self.cCRy(target=target, dcontrol=ctr, angle=-angle / 2,
+                                                                        case=0, control=control)
+                U += gates.CNOT(target, aux) + gates.H(aux) + self.cCRy(target=target, dcontrol=ctr, angle=angle / 2,
+                                                                        case=1, control=control)
+            return U
+
+    def fermionic_excitation(self, angle: typing.Union[Real, Variable, typing.Hashable], indices: typing.List,
+                             control: typing.Union[int, typing.List] = None, opt: bool = True) -> QCircuit:
+        '''
+            Excitation [(i,j),(k,l)],... compiled following https://doi.org/10.1103/PhysRevA.102.062612
+            opt: whether to optimized CNOT H CNOT --> Rz Rz CNOT Rz
+        '''
+        lto = []
+        lfrom = []
+        if isinstance(indices,tuple) and not hasattr(indices[0],"__len__"):
+            indices = [(indices[2 * i], indices[2 * i + 1]) for i in range(len(indices) // 2)]
+        for pair in indices:
+            lfrom.append(pair[0])
+            lto.append(pair[1])
+        Upair = QCircuit()
+        if isinstance(angle, str) or isinstance(angle, tuple):
+            angle = Variable(angle)
+        for i in range(len(lfrom) - 1):
+            Upair += gates.CNOT(lfrom[i + 1], lfrom[i])
+            Upair += gates.CNOT(lto[i + 1], lto[i])
+            Upair += gates.X(lto[i]) + gates.X(lfrom[i])
+        Upair += gates.CNOT(lto[-1], lfrom[-1])
+        crt = lfrom[::-1] + lto
+        Uladder = QCircuit()
+        pairs = lfrom + lto
+        pairs.sort()
+        orbs = []
+        for o in range(len(pairs) // 2):
+            orbs += [*range(pairs[2 * o] + 1, pairs[2 * o + 1])]
+        if len(orbs):
+            for o in range(len(orbs) - 1):
+                Uladder += gates.CNOT(orbs[o], orbs[o + 1])
+            Uladder += compile_circuit(gates.CZ(orbs[-1], lto[-1]))
+        crt.pop(-1)
+        if control is not None and (isinstance(control, int) or len(control) == 1):
+            if isinstance(control, int):
+                crt.append(control)
+            else:
+                crt = crt + control
+            control = []
+        Ur = self.cCRy(target=lto[-1], dcontrol=crt, angle=angle, control=control)
+        Upair2 = Upair.dagger()
+        if opt:
+            Ur.gates.pop(-1)
+            Ur.gates.pop(-1)
+            Upair2.gates.pop(0)
+            Ur += gates.Rz(numpy.pi / 2, target=lto[-1]) + gates.Rz(-numpy.pi / 2, target=lfrom[-1])
+            Ur += gates.CNOT(lto[-1], lfrom[-1]) + gates.Rz(numpy.pi / 2, target=lfrom[-1]) + gates.H(lfrom[-1])
+        return Upair + Uladder + Ur + Uladder.dagger() + Upair2
 
     def __str(self):
         if self.indices is not None:
diff --git a/src/tequila/quantumchemistry/encodings.py b/src/tequila/quantumchemistry/encodings.py
index 12d5ed46..d31ae761 100644
--- a/src/tequila/quantumchemistry/encodings.py
+++ b/src/tequila/quantumchemistry/encodings.py
@@ -2,16 +2,21 @@
 Collections of Fermion-to-Qubit encodings known to tequila
 Most are Interfaces to OpenFermion
 """
+import abc
+
+from tequila import TequilaException
 from tequila.circuit.circuit import QCircuit
-from tequila.circuit.gates import X
+from tequila.circuit.gates import X, CNOT
 from tequila.hamiltonian.qubit_hamiltonian import QubitHamiltonian
 import openfermion
+import numpy
+
 
 def known_encodings():
     # convenience for testing and I/O
-    encodings= {
-        "JordanWigner":JordanWigner,
-        "BravyiKitaev":BravyiKitaev,
+    encodings = {
+        "JordanWigner": JordanWigner,
+        "BravyiKitaev": BravyiKitaev,
         "BravyiKitaevFast": BravyiKitaevFast,
         "BravyiKitaevTree": BravyiKitaevTree,
         "TaperedBravyiKitaev": TaperedBravyKitaev
@@ -22,14 +27,14 @@ def known_encodings():
                  "ReorderedBravyiKitaev": lambda **kwargs: BravyiKitaev(up_then_down=True, **kwargs),
                  "ReorderedBravyiKitaevTree": lambda **kwargs: BravyiKitaevTree(up_then_down=True, **kwargs),
                  }
-    return {k.replace("_","").replace("-","").upper():v for k,v in encodings.items()}
+    return {k.replace("_", "").replace("-", "").upper(): v for k, v in encodings.items()}
 
-class EncodingBase:
 
+class EncodingBase(metaclass=abc.ABCMeta):
     # true if the encoding is fully integrated
     # false: can only do special things (like building the Hamiltionian)
     # but is not consistent with UCC gate generation
-    _ucc_support=False
+    _ucc_support = False
 
     @property
     def supports_ucc(self):
@@ -37,20 +42,20 @@ def supports_ucc(self):
 
     @property
     def name(self):
-        prefix=""
+        prefix = ""
         if self.up_then_down:
-            prefix="Reordered"
+            prefix = "Reordered"
         if hasattr(self, "_name"):
-            return prefix+self._name
+            return prefix + self._name
         else:
-            return prefix+type(self).__name__
+            return prefix + type(self).__name__
 
     def __init__(self, n_electrons, n_orbitals, up_then_down=False, *args, **kwargs):
         self.n_electrons = n_electrons
         self.n_orbitals = n_orbitals
         self.up_then_down = up_then_down
 
-    def __call__(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> QubitHamiltonian:
+    def __call__(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> QubitHamiltonian:
         """
         :param fermion_operator:
             an openfermion FermionOperator
@@ -58,7 +63,8 @@ def __call__(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs
             The openfermion QubitOperator of this class ecoding
         """
         if self.up_then_down:
-            op = openfermion.reorder(operator=fermion_operator, order_function=openfermion.up_then_down, num_modes=2*self.n_orbitals)
+            op = openfermion.reorder(operator=fermion_operator, order_function=openfermion.up_then_down,
+                                     num_modes=2 * self.n_orbitals)
         else:
             op = fermion_operator
 
@@ -73,18 +79,18 @@ def up(self, i):
         if self.up_then_down:
             return i
         else:
-            return 2*i
+            return 2 * i
 
     def down(self, i):
         if self.up_then_down:
-            return i+self.n_orbitals
+            return i + self.n_orbitals
         else:
-            return 2*i+1
+            return 2 * i + 1
 
-    def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
+    def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
         raise Exception("{}::do_transform: called base class".format(type(self).__name__))
 
-    def map_state(self, state:list, *args, **kwargs) -> list:
+    def map_state(self, state: list, *args, **kwargs) -> list:
         """
         Expects a state in spin-orbital ordering
         Returns the corresponding qubit state in the class encoding
@@ -112,7 +118,7 @@ def map_state(self, state:list, *args, **kwargs) -> list:
         # default is a lazy workaround, but it workds
         n_qubits = 2 * self.n_orbitals
 
-        spin_orbitals = sorted([i for i,x in enumerate(state) if int(x)==1])
+        spin_orbitals = sorted([i for i, x in enumerate(state) if int(x) == 1])
 
         string = "1.0 ["
         for i in spin_orbitals:
@@ -128,26 +134,82 @@ def map_state(self, state:list, *args, **kwargs) -> list:
         key = list(wfn.keys())[0].array
         return key
 
-    def hcb_to_me(self, *args, **kwargs):
-        return None
+    @abc.abstractmethod
+    def me_to_jw(self) -> QCircuit:
+        """
+        This method needs to be implemented to enable default conversions via Jordan-Wigner
+        """
+        pass
+
+    # independent conversion methods, these are used for default conversions
+    # arXiv:1808.10402 IV. B. 2, Eq. 57
+    # original: https://doi.org/10.1063/1.4768229
+    def _jw_to_bk(self) -> QCircuit:
+        U = QCircuit()  # Constructs empty circuit
+
+        flipper = False
+        for i in range(self.n_orbitals * 2):
+            # even qubits only hold their own value
+            if i % 2 == 0:
+                continue
+
+            # sum always includes the last qubit
+            U += CNOT(control=i - 1, target=i)
+
+            # every second odd qubit ties together with the last odd qubit
+            if flipper:
+                U += CNOT(control=i - 2, target=i)
+
+            flipper = not flipper
+
+            # we have now created the 4x4 blocks on the diagonal of this operators matrix
+
+            # every power of 2 connects to the last power of 2
+            # this corresponds to the last row in the recursive definitions being all 1s
+            x = numpy.log2(i + 1)
+            if x.is_integer() and x >= 3:
+                x = int(x)
+                U += CNOT(control=2 ** (x - 1) - 1, target=i)
+
+        return U
+
+    def _hcb_to_jw(self):
+        U = QCircuit()
+        for i in range(self.n_orbitals):
+            U += X(target=self.down(i), control=self.up(i))
+        return U
+
+    # Convenience Methods
+    def jw_to_me(self) -> QCircuit:
+        return self.me_to_jw().dagger()
+
+    def me_to_bk(self) -> QCircuit:
+        return self.me_to_jw() + self._jw_to_bk()
+
+    def bk_to_me(self) -> QCircuit:
+        return self.me_to_bk().dagger()
+
+    def hcb_to_me(self) -> QCircuit:
+        return self._hcb_to_jw() + self.jw_to_me()
 
     def __str__(self):
         return type(self).__name__
 
+
 class JordanWigner(EncodingBase):
     """
     OpenFermion::jordan_wigner
     """
-    _ucc_support=True
+    _ucc_support = True
 
-    def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
+    def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
         return openfermion.jordan_wigner(fermion_operator, *args, **kwargs)
 
-    def map_state(self, state:list, *args, **kwargs):
-        state = state + [0]*(self.n_orbitals-len(state))
-        result = [0]*len(state)
+    def map_state(self, state: list, *args, **kwargs):
+        state = state + [0] * (self.n_orbitals - len(state))
+        result = [0] * len(state)
         if self.up_then_down:
-            return [state[2*i] for i in range(self.n_orbitals)] + [state[2*i+1] for i in range(self.n_orbitals)]
+            return [state[2 * i] for i in range(self.n_orbitals)] + [state[2 * i + 1] for i in range(self.n_orbitals)]
         else:
             return state
 
@@ -157,15 +219,34 @@ def hcb_to_me(self, *args, **kwargs):
             U += X(target=self.down(i), control=self.up(i))
         return U
 
+    def me_to_jw(self) -> QCircuit:
+        return QCircuit()
+
+    def jw_to_me(self) -> QCircuit:
+        return QCircuit()
+
+
 class BravyiKitaev(EncodingBase):
     """
     Uses OpenFermion::bravyi_kitaev
     """
 
-    _ucc_support=True
+    _ucc_support = True
+
+    def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
+        return openfermion.bravyi_kitaev(fermion_operator, n_qubits=self.n_orbitals * 2)
+
+    def me_to_jw(self) -> QCircuit:
+        return self._jw_to_bk().dagger()
 
-    def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
-        return openfermion.bravyi_kitaev(fermion_operator, n_qubits=self.n_orbitals*2)
+    def jw_to_me(self) -> QCircuit:
+        return self._jw_to_bk()
+
+    def bk_to_me(self) -> QCircuit:
+        return QCircuit()
+
+    def me_to_bk(self) -> QCircuit:
+        return QCircuit()
 
 
 class BravyiKitaevTree(EncodingBase):
@@ -173,28 +254,37 @@ class BravyiKitaevTree(EncodingBase):
     Uses OpenFermion::bravyi_kitaev_tree
     """
 
-    _ucc_support=True
+    _ucc_support = True
+
+    def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
+        return openfermion.bravyi_kitaev_tree(fermion_operator, n_qubits=self.n_orbitals * 2)
+
+    def me_to_jw(self) -> QCircuit:
+        raise TequilaException("{}::me_to_jw: unimplemented".format(type(self).__name__))
 
-    def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
-        return openfermion.bravyi_kitaev_tree(fermion_operator, n_qubits=self.n_orbitals*2)
 
 class BravyiKitaevFast(EncodingBase):
     """
     Uses OpenFermion::bravyi_kitaev_tree
     """
 
-    _ucc_support=False
+    _ucc_support = False
 
-    def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
+    def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
         n_qubits = openfermion.count_qubits(fermion_operator)
-        if n_qubits != self.n_orbitals*2:
-            raise Exception("BravyiKitaevFast transformation currently only possible for full Hamiltonians (no UCC generators).\nfermion_operator was {}".format(fermion_operator))
+        if n_qubits != self.n_orbitals * 2:
+            raise Exception(
+                "BravyiKitaevFast transformation currently only possible for full Hamiltonians (no UCC generators).\nfermion_operator was {}".format(
+                    fermion_operator))
         op = openfermion.get_interaction_operator(fermion_operator)
         return openfermion.bravyi_kitaev_fast(op)
 
-class TaperedBravyKitaev(EncodingBase):
+    def me_to_jw(self) -> QCircuit:
+        raise TequilaException("{}::me_to_jw: unimplemented".format(type(self).__name__))
 
-    _ucc_support=False
+
+class TaperedBravyKitaev(EncodingBase):
+    _ucc_support = False
 
     """
     Uses OpenFermion::symmetry_conserving_bravyi_kitaev (tapered bravyi_kitaev_tree arxiv:1701.07072)
@@ -202,6 +292,7 @@ class TaperedBravyKitaev(EncodingBase):
     See OpenFermion Documentation for more
     Does not work for UCC generators yet
     """
+
     def __init__(self, n_electrons, n_orbitals, active_fermions=None, active_orbitals=None, *args, **kwargs):
         if active_fermions is None:
             self.active_fermions = n_electrons
@@ -209,7 +300,7 @@ def __init__(self, n_electrons, n_orbitals, active_fermions=None, active_orbital
             self.active_fermions = active_fermions
 
         if active_orbitals is None:
-            self.active_orbitals = n_orbitals*2 # in openfermion those are spin-orbitals
+            self.active_orbitals = n_orbitals * 2  # in openfermion those are spin-orbitals
         else:
             self.active_orbitals = active_orbitals
 
@@ -217,17 +308,20 @@ def __init__(self, n_electrons, n_orbitals, active_fermions=None, active_orbital
             raise Exception("Don't pass up_then_down argument to {}, it can't be changed".format(type(self).__name__))
         super().__init__(n_orbitals=n_orbitals, n_electrons=n_electrons, up_then_down=False, *args, **kwargs)
 
-    def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
-        if openfermion.count_qubits(fermion_operator) != self.n_orbitals*2:
+    def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator:
+        if openfermion.count_qubits(fermion_operator) != self.n_orbitals * 2:
             raise Exception("TaperedBravyiKitaev not ready for UCC generators yet")
-        return openfermion.symmetry_conserving_bravyi_kitaev(fermion_operator, active_orbitals=self.active_orbitals, active_fermions=self.active_fermions)
+        return openfermion.symmetry_conserving_bravyi_kitaev(fermion_operator, active_orbitals=self.active_orbitals,
+                                                             active_fermions=self.active_fermions)
 
-    def map_state(self, state:list, *args, **kwargs):
-        non_tapered_trafo = BravyiKitaevTree(up_then_down=True, n_electrons=self.n_electrons, n_orbitals=self.n_orbitals)
+    def map_state(self, state: list, *args, **kwargs):
+        non_tapered_trafo = BravyiKitaevTree(up_then_down=True, n_electrons=self.n_electrons,
+                                             n_orbitals=self.n_orbitals)
         key = non_tapered_trafo.map_state(state=state, *args, **kwargs)
-        n_qubits = self.n_orbitals*2
+        n_qubits = self.n_orbitals * 2
         active_qubits = [i for i in range(n_qubits) if i not in [n_qubits - 1, n_qubits // 2 - 1]]
         key = [key[i] for i in active_qubits]
         return key
 
-
+    def me_to_jw(self) -> QCircuit:
+        raise TequilaException("{}::me_to_jw: unimplemented".format(type(self).__name__))
diff --git a/src/tequila/quantumchemistry/madness_interface.py b/src/tequila/quantumchemistry/madness_interface.py
index 9273046c..7b291760 100644
--- a/src/tequila/quantumchemistry/madness_interface.py
+++ b/src/tequila/quantumchemistry/madness_interface.py
@@ -420,7 +420,10 @@ def make_upccgsd_ansatz(self, name="UpCCGSD", label=None, direct_compiling=None,
 
         """
         # check if the used qubit encoding has a hcb transformation
-        have_hcb_trafo = self.transformation.hcb_to_me() is not None
+        try:
+            have_hcb_trafo = self.transformation.hcb_to_me() is not None
+        except:
+            have_hcb_trafo = False
         name = name.upper()
 
         # Default Method
diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index d8a6882b..7364a58e 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -17,7 +17,7 @@
 
 import typing, numpy, numbers
 from itertools import product
-
+import tequila.grouping.fermionic_functions as ff
 
 
 try:
@@ -32,8 +32,7 @@
     except Exception as E:
         raise Exception("{}\nIssue with Tequila Chemistry: Please update openfermion".format(str(E)))
 import warnings
-
-
+OPTIMIZED_ORDERING = "Optimized"
 class QuantumChemistryBase:
     """
     Base Class for tequila chemistry functionality
@@ -113,7 +112,11 @@ def from_tequila(cls, molecule, transformation=None, *args, **kwargs):
         h1 = molecule.integral_manager.one_body_integrals
         h2 = molecule.integral_manager.two_body_integrals
         S = molecule.integral_manager.overlap_integrals
-        active_orbitals = [o.idx_total for o in molecule.integral_manager.active_orbitals]
+        if "active_orbitals" not in kwargs:
+            active_orbitals = [o.idx_total for o in molecule.integral_manager.active_orbitals]
+        else:
+            active_orbitals = kwargs["active_orbitals"]
+            kwargs.pop("active_orbitals")
         if transformation is None:
             transformation = molecule.transformation
         parameters = molecule.parameters
@@ -125,7 +128,8 @@ def from_tequila(cls, molecule, transformation=None, *args, **kwargs):
                    active_orbitals= active_orbitals,
                    transformation=transformation,
                    orbital_type=molecule.integral_manager._orbital_type,
-                   parameters=parameters, *args, **kwargs)
+                   parameters=parameters,
+                   reference_orbitals= molecule.integral_manager.active_space.reference_orbitals,*args, **kwargs)
 
     def supports_ucc(self):
         """
@@ -441,10 +445,14 @@ def make_excitation_gate(self, indices, angle, control=None, assume_real=True, *
 
         generator = self.make_excitation_generator(indices=indices, remove_constant_term=control is None)
         p0 = self.make_excitation_generator(indices=indices, form="P0", remove_constant_term=control is None)
-
+        if self.transformation.up_then_down:
+            idx = []
+            for pair in indices:
+                idx.append((pair[0]//2+(pair[0]%2)*self.n_orbitals,pair[1]//2+(pair[1]%2)*self.n_orbitals))
+        else:idx = indices
         return QCircuit.wrap_gate(
             FermionicGateImpl(angle=angle, generator=generator, p0=p0,
-                              transformation=type(self.transformation).__name__.lower(), indices=indices,
+                              transformation=type(self.transformation).__name__.lower(), indices=idx,
                               assume_real=assume_real,
                               control=control, **kwargs))
 
@@ -1239,7 +1247,13 @@ def make_upccgsd_ansatz(self,
         indices = self.make_upccgsd_indices(key=name)
 
         # check if the used qubit encoding has a hcb transformation
-        have_hcb_trafo = self.transformation.hcb_to_me() is not None
+        have_hcb_trafo = True
+        try:
+            if self.transformation.hcb_to_me() is None:
+                have_hcb_trafo = False
+        except:
+            have_hcb_trafo = False
+        
 
         # consistency checks for optimization
         if have_hcb_trafo and hcb_optimization is None and include_reference:
@@ -1487,7 +1501,8 @@ def make_uccsd_ansatz(self, trotter_steps: int = 1,
         factor = 1.0 / trotter_steps
         for step in range(trotter_steps):
             for idx, angle in indices.items():
-                UCCSD += self.make_excitation_gate(indices=idx, angle=factor * angle)
+                converted = [(idx[2 * i], idx[2 * i + 1]) for i in range(len(idx) // 2)]
+                UCCSD += self.make_excitation_gate(indices=converted, angle=factor * angle)
         if hasattr(initial_amplitudes,
                    "lower") and initial_amplitudes.lower() == "mp2" and parametrized and add_singles:
             # mp2 has no singles, need to initialize them here (if not parametrized initializling as 0.0 makes no sense though)
@@ -1721,7 +1736,8 @@ def rdm2(self):
             return None
 
     def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_free: bool = True,
-                     get_rdm1: bool = True, get_rdm2: bool = True, ordering="dirac", use_hcb: bool = False):
+                     get_rdm1: bool = True, get_rdm2: bool = True, ordering="dirac", use_hcb: bool = False,
+                     rdm_trafo: QubitHamiltonian = None, decompose=None):
         """
         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.
@@ -1749,6 +1765,9 @@ def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_fre
         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.
+        rdm_trafo :
+            The rdm operators can be transformed, e.g., a^dagger_i a_j -> U^dagger a^dagger_i a_j U,
+            where U represents the transformation. The default is set to None, implying that U equas the identity.
 
         Returns
         -------
@@ -1993,8 +2012,18 @@ def _build_2bdy_operators_hcb() -> list:
         # Transform operator lists to QubitHamiltonians
         if (not use_hcb):
             qops = [_get_qop_hermitian(op) for op in qops]
+        
         # Compute expected values
-        evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)
+        if rdm_trafo is None:
+            if decompose is not None:
+                print("MANIPULATED")
+                X = decompose(H=qops, U=U)
+                evals = simulate(X, variables=variables)
+            else:
+                evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)
+        else:
+            qops = [rdm_trafo.dagger()*qops[i]*rdm_trafo for i in range(len(qops))]
+            evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)
 
         # Assemble density matrices
         # If self._rdm1, self._rdm2 exist, reset them if they are of the other spin-type
@@ -2127,6 +2156,58 @@ def perturbative_f12_correction(self, rdm1: numpy.ndarray = None, rdm2: numpy.nd
                                                    n_ri=n_ri, external_info=external_info, **kwargs)
         return correction.compute()
 
+    def n_rotation(self, i, phi):
+        '''
+        Creates a quantum circuit that applies a phase rotation based on phi to both components (up and down) of a given qubit.
+
+        Parameters:
+        - i (int): The index of the qubit to which the rotation will be applied.
+        - phi (float): The rotation angle. The actual rotation applied will be multiplied with -2 for both components.
+
+        Returns:
+        - QCircuit: A quantum circuit object containing the sequence of rotations applied to the up and down components of the specified qubit.
+        '''
+
+        # Generate number operators for the up and down components of the qubit.
+        n_up = self.make_number_op(2*i)
+        n_down = self.make_number_op(2*i+1)
+
+        # Start a new circuit and apply rotations to each component.
+        circuit = gates.GeneralizedRotation(generator = n_up, angle=-2*phi)
+        circuit += gates.GeneralizedRotation(generator = n_down, angle=-2*phi)
+        return circuit
+        
+    def get_givens_circuit(self, unitary, tol = 1e-12, ordering = OPTIMIZED_ORDERING):
+        '''
+        Constructs a quantum circuit from a given real unitary matrix using Givens rotations.
+        
+        This method decomposes a unitary matrix into a series of Givens and Rz (phase) rotations,
+        then constructs and returns a quantum circuit that implements this sequence of rotations.
+
+        Parameters:
+        - unitary (numpy.array): A real unitary matrix representing the transformation to implement.
+        - tol (float): A tolerance threshold below which matrix elements are considered zero.
+        - ordering (list of tuples or 'Optimized'): Custom ordering of indices for Givens rotations or 'Optimized' to generate them automatically.
+
+        Returns:
+        - QCircuit: A quantum circuit implementing the series of rotations decomposed from the unitary.
+        '''
+        # Decompose the unitary matrix into Givens and phase (Rz) rotations.
+        theta_list, phi_list = get_givens_decomposition(unitary, tol, ordering)
+
+        # Initialize an empty quantum circuit.
+        circuit = QCircuit()
+
+        # Add all Rz (phase) rotations to the circuit.
+        for phi in phi_list:
+            circuit += self.n_rotation(phi[1], phi[0])
+
+        # Add all Givens rotations to the circuit.
+        for theta in reversed(theta_list):
+            circuit += self.UR(theta[1], theta[2], theta[0]*2)
+
+        return circuit
+
 
     def print_basis_info(self):
         return self.integral_manager.print_basis_info()
@@ -2147,3 +2228,140 @@ def __str__(self) -> str:
         result += "\nmore information with: self.print_basis_info()\n"
 
         return result
+
+def givens_matrix(n, p, q, theta):
+    '''
+    Construct a complex Givens rotation matrix of dimension n by theta between rows/columns p and q.
+    '''
+    '''
+    Generates a Givens rotation matrix of size n x n to rotate by angle theta in the (p, q) plane. This matrix can be complex
+
+    Parameters:
+    - n (int): The size of the Givens rotation matrix.
+    - p (int): The first index for the rotation plane.
+    - q (int): The second index for the rotation plane.
+    - theta (float): The rotation angle.
+
+    Returns:
+    - numpy.array: The Givens rotation matrix.
+    '''
+    matrix = numpy.eye(n)  # Matrix to hold complex numbers
+    cos_theta = numpy.cos(theta)
+    sin_theta = numpy.sin(theta)
+
+    # Directly assign cosine and sine without complex phase adjustment
+    matrix[p, p] = cos_theta
+    matrix[q, q] = cos_theta
+    matrix[p, q] = sin_theta
+    matrix[q, p] = -sin_theta
+
+    return matrix
+
+def get_givens_decomposition(unitary, tol = 1e-12, ordering = OPTIMIZED_ORDERING, return_diagonal = False):
+    '''
+    Decomposes a real unitary matrix into Givens rotations (theta) and Rz rotations (phi).
+
+    Parameters:
+    - unitary (numpy.array): A real unitary matrix to decompose. It cannot be complex.
+    - tol (float): Tolerance for considering matrix elements as zero. Elements with absolute value less than tol are treated as zero.
+    - ordering (list of tuples or 'Optimized'): Custom ordering of indices for Givens rotations or 'Optimized' to generate them automatically.
+    - return_diagonal (bool): If True, the function also returns the diagonal matrix as part of the output.
+
+    Returns:
+    - list: A list of tuples, each representing a Givens rotation. Each tuple contains the rotation angle theta and indices (i,j) of the rotation.
+    - list: A list of tuples, each representing an Rz rotation. Each tuple contains the rotation angle phi and the index (i) of the rotation.
+    - numpy.array (optional): The diagonal matrix after applying all Givens rotations, returned if return_diagonal is True.
+    '''
+    U = unitary # no need to copy as we don't modify the original
+    U[abs(U) < tol] = 0 # Zeroing out the small elements as per the tolerance level.
+    n = U.shape[0]
+
+    # Determine optimized ordering if specified.
+    if ordering == OPTIMIZED_ORDERING:
+        ordering = ff.depth_eff_order_mf(n)
+
+    theta_list = []
+    phi_list = []
+
+    def calcTheta(U, c, r):
+        '''Calculate and apply the Givens rotation for a specific matrix element.'''
+        t = numpy.arctan2(-U[r,c], U[r-1,c])
+        theta_list.append((t, r, r-1))
+        g = givens_matrix(n,r,r-1,t)
+        U = numpy.dot(g, U)
+
+        return U
+
+    # Apply and store Givens rotations as per the given or computed ordering.
+    if ordering is None:
+        for c in range(n):
+            for r in range(n-1, c, -1):
+                U = calcTheta(U, c, r)
+    else:
+        for r, c in ordering:
+            U = calcTheta(U, c, r)
+    
+    # Calculating the Rz rotations based on the phases of the diagonal elements.
+    # For real elements this means a 180 degree shift, i.e. a sign change.
+    for i in range(n):
+        ph = numpy.angle(U[i,i])
+        phi_list.append((ph, i))
+
+    # Filtering out rotations without significance.
+    theta_list_new = []
+    for i, theta in enumerate(theta_list):
+        if abs(theta[0] % (2*numpy.pi)) > tol:
+            theta_list_new.append(theta)
+    
+    phi_list_new = []
+    for i, phi in enumerate(phi_list):
+        if abs(phi[0]) > tol:
+            phi_list_new.append(phi)
+        
+    if return_diagonal:
+        # Optionally return the resulting diagonal
+        return theta_list_new, phi_list_new, U
+    else:
+        return theta_list_new, phi_list_new
+    
+def reconstruct_matrix_from_givens(n, theta_list, phi_list, to_real_if_possible = True, tol = 1e-12):
+    '''
+    Reconstructs a matrix from given Givens rotations and Rz diagonal rotations.
+    This function is effectively an inverse of get_givens_decomposition, and therefore only works with data in the same format as its output.
+
+    Parameters:
+    - n (int): The size of the unitary matrix to be reconstructed.
+    - theta_list (list of tuples): Each tuple contains (angle, i, j) representing a Givens rotation of `angle` radians, applied to rows/columns `i` and `j`.
+    - phi_list (list of tuples): Each tuple contains (angle, i), representing an Rz rotation by `angle` radians applied to the `i`th diagonal element.
+    - to_real_if_possible (bool): If True, converts the matrix to real if its imaginary part is effectively zero.
+    - tol (float): The tolerance whether to swap a complex rotation for a sign change.
+
+    Returns:
+    - numpy.ndarray: The reconstructed complex or real matrix, depending on the `to_real_if_possible` flag and matrix composition.
+    '''
+    # Start with an identity matrix
+    reconstructed = numpy.eye(n, dtype=complex)
+    
+    # Apply Rz rotations for diagonal elements
+    for phi in phi_list:
+        angle, i = phi
+        # Directly apply a sign flip if the rotation angle is π
+        if numpy.isclose(angle, numpy.pi, atol=tol):
+            reconstructed[i, i] *= -1
+        else:
+            reconstructed[i, i] *= numpy.exp(1j * angle)
+        
+    # Apply Givens rotations in reverse order
+    for theta in reversed(theta_list):
+        angle, i, j = theta
+        g = givens_matrix(n, i, j, angle)
+        reconstructed = numpy.dot(g.conj().T, reconstructed) # Transpose of Givens matrix applied to the left
+
+    # Convert matrix to real if its imaginary part is negligible unless disabled via to_real_if_possible
+    if to_real_if_possible:
+        # Directly apply a sign flip if the rotation angle is π
+        if numpy.all(reconstructed.imag == 0):
+            # Convert to real by taking the real part
+            reconstructed = reconstructed.real
+
+    return reconstructed
diff --git a/src/tequila/simulators/simulator_api.py b/src/tequila/simulators/simulator_api.py
index 63506ca6..9db3991a 100755
--- a/src/tequila/simulators/simulator_api.py
+++ b/src/tequila/simulators/simulator_api.py
@@ -36,9 +36,12 @@
     HAS_QISKIT = True
     INSTALLED_SIMULATORS["qiskit"] = BackendTypes(BackendCircuitQiskit, BackendExpectationValueQiskit)
     INSTALLED_SAMPLERS["qiskit"] = BackendTypes(BackendCircuitQiskit, BackendExpectationValueQiskit)
-    INSTALLED_NOISE_SAMPLERS["qiskit"] = BackendTypes(BackendCircuitQiskit, BackendExpectationValueQiskit)
+    from tequila.simulators.simulator_qiskit import HAS_NOISE as HAS_QISKIT_NOISE
+    if HAS_QISKIT_NOISE:
+        INSTALLED_NOISE_SAMPLERS["qiskit"] = BackendTypes(BackendCircuitQiskit, BackendExpectationValueQiskit)
 except ImportError:
     HAS_QISKIT = False
+    HAS_QISKIT_NOISE = False
 
 HAS_QIBO = True
 try:
@@ -133,6 +136,8 @@ def show_available_simulators():
                                                              str(k in INSTALLED_SAMPLERS),
                                                              str(k in INSTALLED_NOISE_SAMPLERS),
                                                              str(k in INSTALLED_BACKENDS)))
+    if HAS_QISKIT and not HAS_QISKIT_NOISE:
+        print("missing qiskit_aer: no noisy simulation")
 
 
 def pick_backend(backend: str = None, samples: int = None, noise: NoiseModel = None, device=None,
diff --git a/src/tequila/simulators/simulator_qiskit.py b/src/tequila/simulators/simulator_qiskit.py
index ded4e784..9d2130aa 100755
--- a/src/tequila/simulators/simulator_qiskit.py
+++ b/src/tequila/simulators/simulator_qiskit.py
@@ -3,12 +3,20 @@
 from tequila import TequilaException, TequilaWarning
 from tequila import BitString, BitNumbering, BitStringLSB
 from tequila.utils.keymap import KeyMapRegisterToSubregister
-import qiskit, numpy, warnings
-import qiskit.providers.aer.noise as qiskitnoise
 from tequila.utils import to_float
-import qiskit.test.mock.backends
-from qiskit.providers.ibmq import IBMQBackend
+import qiskit, numpy, warnings
 
+HAS_NOISE=True
+try:
+    from qiskit_aer import noise as qiskitnoise
+except:
+    HAS_NOISE = False
+
+HAS_IBMQ=True
+try:
+    from qiskit.providers.ibmq import IBMQBackend
+except:
+    HAS_IBMQ=False
 
 def get_bit_flip(p):
     """
@@ -137,7 +145,7 @@ class BackendCircuitQiskit(BackendCircuit):
     }
 
     numbering = BitNumbering.LSB
-
+            
     def __init__(self, abstract_circuit: QCircuit, variables, qubit_map=None, noise=None,
                  device=None, *args, **kwargs):
         """
@@ -169,7 +177,7 @@ def __init__(self, abstract_circuit: QCircuit, variables, qubit_map=None, noise=
             'Rx': (lambda c: c.rx, lambda c: c.mcrx),
             'Ry': (lambda c: c.ry, lambda c: c.mcry),
             'Rz': (lambda c: c.rz, lambda c: c.mcrz),
-            'Phase': (lambda c: c.u1, lambda c: c.cu1),
+            'Phase': (lambda c: c.p, lambda c: c.cp),
             'SWAP': (lambda c: c.swap, lambda c: c.cswap),
         }
 
diff --git a/src/tequila/tools/random_generators.py b/src/tequila/tools/random_generators.py
index 07f3a7cd..434022fc 100644
--- a/src/tequila/tools/random_generators.py
+++ b/src/tequila/tools/random_generators.py
@@ -2,6 +2,7 @@
 from tequila.circuit import gates
 from tequila.circuit.circuit import QCircuit
 from tequila.hamiltonian.qubit_hamiltonian import QubitHamiltonian
+from scipy.stats import unitary_group, ortho_group
 
 def make_random_circuit(n_qubits: int, rotation_gates: list=['rx', 'ry', 'rz'], n_rotations: int=None,
                         enable_controls: bool=None) -> QCircuit:
@@ -75,3 +76,19 @@ def make_random_hamiltonian(n_qubits: int , paulis: list=['X','Y','Z'], n_ps: in
     
     H = QubitHamiltonian(ham)
     return H
+
+def generate_random_unitary(size, complex = False):
+    '''
+    Generates a random unitary (or furthermore orthogonal if complex is False) matrix of a specified size.
+
+    Parameters:
+    - size (int): The size of the unitary matrix to be generated.
+    - complex (bool, optional): Whether the unitary should be complex.
+
+    Returns:
+    - numpy.ndarray: A randomly generated unitary matrix.
+    '''
+    if complex:
+        return unitary_group.rvs(size)
+    else:
+        return ortho_group.rvs(size)
\ No newline at end of file
diff --git a/src/tequila/version.py b/src/tequila/version.py
index 09ed92e6..88cc42f3 100644
--- a/src/tequila/version.py
+++ b/src/tequila/version.py
@@ -1,2 +1,2 @@
-__version__ = "1.9.5"
+__version__ = "1.9.6"
 __author__ = "Tequila Developers "


Test patch:
diff --git a/tests/test_chemistry.py b/tests/test_chemistry.py
index 86b2dc33..d24c6bfa 100644
--- a/tests/test_chemistry.py
+++ b/tests/test_chemistry.py
@@ -8,6 +8,8 @@
 from tequila.objective import ExpectationValue
 from tequila.quantumchemistry.encodings import known_encodings
 from tequila.simulators.simulator_api import simulate
+import tequila.quantumchemistry.qc_base as qcb
+import tequila.tools.random_generators as rg
 
 HAS_PYSCF = "pyscf" in qc.INSTALLED_QCHEMISTRY_BACKENDS
 HAS_PSI4 = "psi4" in qc.INSTALLED_QCHEMISTRY_BACKENDS
@@ -171,8 +173,6 @@ def test_ucc_singles_psi4():
 def do_test_ucc(qc_interface, parameters, result, trafo, backend="qulacs"):
     # check examples for comments
     psi4_interface = qc_interface(parameters=parameters, transformation=trafo)
-
-    hqc = psi4_interface.make_hamiltonian()
     amplitudes = psi4_interface.compute_ccsd_amplitudes()
     U = psi4_interface.make_uccsd_ansatz(trotter_steps=1, initial_amplitudes=amplitudes, include_reference_ansatz=True)
     variables = amplitudes.make_parameter_dictionary()
@@ -368,9 +368,9 @@ def test_hamiltonian_reduction(backend):
 
 @pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4/pyscf not found")
 @pytest.mark.parametrize("assume_real", [True, False])
-@pytest.mark.parametrize("trafo", ["jordan_wigner", "bravyi_kitaev", "tapered_bravyi_kitaev"])
+@pytest.mark.parametrize("trafo", ["jordan_wigner", "bravyi_kitaev", "reordered_jordan_wigner"])
 def test_fermionic_gates(assume_real, trafo):
-    mol = tq.chemistry.Molecule(geometry="H 0.0 0.0 0.7\nLi 0.0 0.0 0.0", basis_set="sto-3g")
+    mol = tq.chemistry.Molecule(geometry="H 0.0 0.0 0.7\nLi 0.0 0.0 0.0", basis_set="sto-3g",transformation=trafo)
     U1 = mol.prepare_reference()
     U2 = mol.prepare_reference()
     variable_count = {}
@@ -587,8 +587,14 @@ def test_spa_consistency(geometry, name, optimize, transformation):
     mol = tq.Molecule(geometry=geometry, basis_set="sto-3g",
                       transformation=transformation).use_native_orbitals()
 
-    if mol.transformation.hcb_to_me() is None:
+    # test compares HCB and non-HCB SPA implementations
+    # conversion needs to be implemented for comparisson
+    # if not, HCB is not used in  circuit optimization anyways
+    try:
+        mol.transformation.hcb_to_me()
+    except:
         return
+        
     # doesn't need to make physical sense for the test
     edges = []
     i = 0
@@ -725,3 +731,55 @@ def test_orbital_optimization_hcb(geometry):
     assert numpy.isclose(opt1.energy,opt2.energy,atol=1.e-5)
     assert time1 < time2
     assert (numpy.isclose(opt1.mo_coeff,opt2.mo_coeff, atol=1.e-5)).all()
+
+@pytest.mark.parametrize("transformation", ["JordanWigner", "ReorderedJordanWigner", "BravyiKitaev", "BravyiKitaevTree"])
+@pytest.mark.parametrize("size", [2, 8])
+def test_givens_on_molecule(size, transformation):
+    # dummy one-electron integrals
+    h = numpy.ones(shape=[size,size])
+    # dummy two-electron integrals
+    g = numpy.ones(shape=[size, size, size, size])
+
+    U = rg.generate_random_unitary(size)
+
+    # transformed integrals
+    th = (U.T.dot(h)).dot(U)
+    tg = numpy.einsum("ijkx, xl -> ijkl", g, U, optimize='greedy')
+    tg = numpy.einsum("ijxl, xk -> ijkl", tg, U, optimize='greedy')
+    tg = numpy.einsum("ixkl, xj -> ijkl", tg, U, optimize='greedy')
+    tg = numpy.einsum("xjkl, xi -> ijkl", tg, U, optimize='greedy')
+
+    # original molecule/H
+    mol = tq.Molecule(geometry="He 0.0 0.0 0.0", nuclear_repulsion=0.0, one_body_integrals=h, two_body_integrals=g, basis_set="dummy", transformation=transformation)
+    H = mol.make_hamiltonian()
+    # transformed molecule/H
+    tmol = tq.Molecule(geometry="He 0.0 0.0 0.0", nuclear_repulsion=0.0, one_body_integrals=th, two_body_integrals=tg,basis_set="dummy", transformation=transformation)
+    tH = tmol.make_hamiltonian()
+
+    # transformation in qubit space (this corresponds to the U above)
+    UR = mol.get_givens_circuit(U) # Works!
+
+    # test circuit
+    circuit = rg.make_random_circuit(size)
+
+    # create expectation values and see if they are the same
+    E1 = tq.ExpectationValue(U=circuit, H=tH)
+    E2 = tq.ExpectationValue(U=circuit + UR, H=H)
+
+    result1 = tq.simulate(E1)
+    result2 = tq.simulate(E2)
+    
+    assert numpy.isclose(result1, result2)
+
+@pytest.mark.parametrize("size", [2, 8])
+def test_givens_decomposition(size):
+    # generate random unitary
+    unitary = rg.generate_random_unitary(size)
+
+    # decompose givens
+    theta_list, phi_list = qcb.get_givens_decomposition(unitary)
+
+    # reconstruct original unitary from givens
+    reconstructed_matrix = qcb.reconstruct_matrix_from_givens(unitary.shape[0], theta_list, phi_list)
+    
+    assert numpy.allclose(unitary, reconstructed_matrix)
diff --git a/tests/test_chemistry_madness.py b/tests/test_chemistry_madness.py
index ae22ea7e..b398e26c 100644
--- a/tests/test_chemistry_madness.py
+++ b/tests/test_chemistry_madness.py
@@ -135,7 +135,7 @@ def test_madness_upccgsd(trafo):
     E = tq.ExpectationValue(H=H, U=U)
     assert (len(E.extract_variables()) == 2)
     variables = result.variables
-    if "bravyi" in trafo.lower():
+    if "bravyikitaevtree" in trafo.lower():
         # signs of angles change in BK compared to JW-like HCB
         variables = {k: -v for k, v in variables.items()}
     energy = tq.simulate(E, variables)