From c57ad618186f1a1c647effba62e489c08ea32bb7 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Thu, 29 Feb 2024 11:08:55 +0100
Subject: [PATCH 01/29] Update to 1.9.4 (#330)

* fixing syntax issue in post_init of dataclass (#327)

* fixing syntax issue in post_init of dataclass
* phoenics dropped due to maintenance resources
* more convenient randomization initialization for OO, avoiding numpy warnings

* Update version.py (#329)
---
 requirements.txt                              |   1 -
 src/tequila/apps/adapt/adapt.py               |   2 +-
 src/tequila/optimizers/__init__.py            |  15 +-
 src/tequila/optimizers/optimizer_gpyopt.py    |   2 -
 src/tequila/optimizers/optimizer_phoenics.py  | 348 ------------------
 .../quantumchemistry/chemistry_tools.py       |   4 +-
 .../quantumchemistry/madness_interface.py     |   2 +-
 .../quantumchemistry/orbital_optimizer.py     |  10 +-
 src/tequila/quantumchemistry/qc_base.py       |   2 +-
 src/tequila/version.py                        |   2 +-
 tests/test_noise_opt.py                       |  14 -
 tests/test_phoenics.py                        |  68 ----
 12 files changed, 13 insertions(+), 457 deletions(-)
 delete mode 100644 src/tequila/optimizers/optimizer_phoenics.py
 delete mode 100644 tests/test_phoenics.py

diff --git a/requirements.txt b/requirements.txt
index 6f9ae975..0dd2e329 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -20,7 +20,6 @@ qulacs # default simulator (best integration), remove if the installation gives
 #qibo <= 0.1.1 # can not be installed in the same environment as gpyopt
 
 #optional optimizers
-#phoenics # version on PyPi  isc urrently broken, we recommend to install from source (AAG github)
 #gpyopt # not in combination with qibo as quantum backend
 
 #optional third party libraries
diff --git a/src/tequila/apps/adapt/adapt.py b/src/tequila/apps/adapt/adapt.py
index 6b8dd9e9..4b9a5549 100644
--- a/src/tequila/apps/adapt/adapt.py
+++ b/src/tequila/apps/adapt/adapt.py
@@ -21,7 +21,7 @@ class AdaptParameters:
     degeneracy_threshold: float = 5.e-4
     silent: bool = False
     
-    def __post__init__(self):
+    def __post_init__(self):
         # avoid stacking of same operator-types in a row
         if "method_options" in self.optimizer_args:
             if "gtol" in self.optimizer_args["method_options"]:
diff --git a/src/tequila/optimizers/__init__.py b/src/tequila/optimizers/__init__.py
index 0022f7de..032606c1 100644
--- a/src/tequila/optimizers/__init__.py
+++ b/src/tequila/optimizers/__init__.py
@@ -16,7 +16,7 @@ class _Optimizers:
     methods: list = None
 
 
-SUPPORTED_OPTIMIZERS = ['scipy', 'phoenics', 'gpyopt', 'gd']
+SUPPORTED_OPTIMIZERS = ['scipy', 'gpyopt', 'gd']
 INSTALLED_OPTIMIZERS = {}
 INSTALLED_OPTIMIZERS['scipy'] = _Optimizers(cls=OptimizerSciPy,
                                             minimize=minimize_scipy,
@@ -37,19 +37,6 @@ class _Optimizers:
 except ImportError:
     has_gpyopt = False
 
-has_phoenics = False
-try:
-    from tequila.optimizers.optimizer_phoenics import OptimizerPhoenics
-    from tequila.optimizers.optimizer_phoenics import minimize as minimize_phoenics
-
-    INSTALLED_OPTIMIZERS['phoenics'] = _Optimizers(cls=OptimizerPhoenics,
-                                                   minimize=minimize_phoenics,
-                                                   methods=OptimizerPhoenics.available_methods())
-    has_phoenics = True
-except ImportError:
-    has_phoenics = False
-
-
 def show_available_optimizers(module=None):
     """
     Returns
diff --git a/src/tequila/optimizers/optimizer_gpyopt.py b/src/tequila/optimizers/optimizer_gpyopt.py
index 50834e07..cfde41b9 100644
--- a/src/tequila/optimizers/optimizer_gpyopt.py
+++ b/src/tequila/optimizers/optimizer_gpyopt.py
@@ -4,8 +4,6 @@
 import numbers
 from tequila.objective.objective import Variable
 import warnings
-
-warnings.simplefilter("ignore")
 import GPyOpt
 from GPyOpt.methods import BayesianOptimization
 import numpy as np
diff --git a/src/tequila/optimizers/optimizer_phoenics.py b/src/tequila/optimizers/optimizer_phoenics.py
deleted file mode 100644
index e2aa77f5..00000000
--- a/src/tequila/optimizers/optimizer_phoenics.py
+++ /dev/null
@@ -1,348 +0,0 @@
-from tequila.objective.objective import Objective
-from tequila.optimizers.optimizer_base import Optimizer, OptimizerResults, dataclass
-import typing
-import numbers
-from tequila.objective.objective import Variable
-import copy
-import warnings
-import pickle
-import time
-from tequila import TequilaException
-warnings.simplefilter("ignore")
-with warnings.catch_warnings():
-    warnings.filterwarnings("ignore", category=DeprecationWarning)
-    warnings.filterwarnings("ignore")
-import phoenics
-
-import numpy as np
-from numpy import pi as pi
-from tequila.simulators.simulator_api import compile_objective
-import os
-
-#numpy, tf, etc can get real, real, real, noisy here. We suppress it.
-warnings.filterwarnings('ignore', category=DeprecationWarning)
-os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'
-warnings.filterwarnings('ignore', category=FutureWarning)
-
-@dataclass
-class PhoenicsResults(OptimizerResults):
-
-    observations: list = None
-    phoenics_instance: phoenics.Phoenics = None
-
-class OptimizerPhoenics(Optimizer):
-    """
-    wrapper to allow optimization of objectives with Phoenics, a bayesian optimizer.
-    See: https://github.com/aspuru-guzik-group/phoenics
-    """
-    @classmethod
-    def available_methods(cls):
-        return ["phoenics"]
-
-    def __init__(self, maxiter, backend=None, save_history=True, minimize=True,
-                 samples=None, silent=None, noise=None, device=None):
-        self._minimize = minimize
-
-        super().__init__(backend=backend, maxiter=maxiter, samples=samples,
-                         noise=noise,device=device,
-                         save_history=save_history, silent=silent)
-
-    def _process_for_sim(self, recommendation, passive_angles):
-        """
-        convert from the phoenics suggestion format to a version recognizable by objectives.
-        Parameters
-        ----------
-        recommendation: dict:
-            the a phoenics suggestion.
-        passive_angles: dict:
-            passive angles not optimized over.
-
-        Returns
-        -------
-        dict:
-            dict of Bariable, float pairs.
-        """
-        rec = copy.deepcopy(recommendation)
-        for part in rec:
-            for k, v in part.items():
-                part[k] = v.item()
-            if passive_angles is not None:
-                for k, v in passive_angles.items():
-                    part[k] = v
-        return rec
-
-    def _process_for_phoenics(self, pset, result, passive_angles=None):
-        """
-        Convert results of a call to an objective into a form interpretable by phoenics.
-        Parameters
-        ----------
-        pset: dict:
-            the parameters evaluated, as a dictionary
-        result:
-            the result of calling some objective, using pset as parameters.
-        passive_angles: dict, optional:
-            passive_angles, not optimized over.
-
-        Returns
-        -------
-        dict:
-            the a dictionary, formatted as phoenics prefers it, for use as an 'observation'.
-        """
-        new = copy.deepcopy(pset)
-        for k, v in new.items():
-            new[k] = np.array([v], dtype=np.float32)
-        if passive_angles is not None:
-            for k in passive_angles.keys():
-                del new[k]
-        new['Energy'] = result
-
-        return new
-
-    def _make_phoenics_object(self, objective, passive_angles=None, conf=None, *args, **kwargs):
-        """
-        instantiate phoenics, to perform optimization.
-
-        Parameters
-        ----------
-        objective: Objective:
-            the objective to optimize over.
-        passive_angles: dict, optional:
-            a dictionary of angles not to optimize over.
-        conf: optional:
-            a user built configuration object or file, from which to initialize a phoenics object.
-            For advanced users only.
-        args
-        kwargs
-
-        Returns
-        -------
-        phoenics.Phoenics
-            a phoenics object configured to optimize an objective.
-        """
-        if conf is not None:
-            if hasattr(conf, 'readlines'):
-                bird = phoenics.Phoenics(config_file=conf)
-            else:
-                bird = phoenics.Phoenics(config_dict=conf)
-
-            return bird
-        op = objective.extract_variables()
-        if passive_angles is not None:
-            for i, thing in enumerate(op):
-                if thing in passive_angles.keys():
-                    op.remove(thing)
-
-        config = {"general": {"auto_desc_gen": "False", "batches": 5, "boosted": "False", "parallel": "False", "scratch_dir":os.getcwd()}}
-        config['parameters'] = [
-            {'name': k, 'periodic': 'True', 'type': 'continuous', 'size': 1, 'low': 0, 'high': 2 * pi} for k in op]
-        if self._minimize is True:
-            config['objectives'] = [{"name": "Energy", "goal": "minimize"}]
-        else:
-            config['objectives'] = [{"name": "Energy", "goal": "maximize"}]
-
-        for k,v in kwargs.items():
-            if hasattr(k, "lower") and k.lower() in config["general"]:
-                config["general"][k.lower()] = v
-
-        bird = phoenics.Phoenics(config_dict=config)
-        return bird
-
-    def __call__(self, objective: Objective,
-                 maxiter=None,
-                 variables: typing.List[Variable] = None,
-                 initial_values: typing.Dict[Variable, numbers.Real] = None,
-                 previous=None,
-                 phoenics_config=None,
-                 file_name=None,
-                 *args,
-                 **kwargs):
-        """
-        Perform optimization with phoenics.
-
-        Parameters
-        ----------
-        objective: Objective
-            the objective to optimize.
-        maxiter: int:
-            (Default value = None)
-            if not None, overwrite the init maxiter with new number.
-        variables: list:
-            (Default value = None)
-            which variables to optimize over. If None: all of the variables in objective are used.
-        initial_values: dict:
-            (Default value = None)
-            an initial point to begin optimization from. Random, if None.
-        previous:
-            previous observations, formatted for phoenics, to use in optimization. For use by advanced users.
-        phoenics_config:
-            a config for a phoenics object.
-        file_name:
-            a file
-        args
-        kwargs
-
-        Returns
-        -------
-        PhoenicsResults:
-            the results of optimization by phoenics.
-
-        """
-
-        objective = objective.contract()
-        active_angles, passive_angles, variables = self.initialize_variables(objective,
-                                                               initial_values=initial_values,
-                                                               variables=variables)
-
-        if maxiter is None:
-            maxiter = 10
-
-        obs = []
-        bird = self._make_phoenics_object(objective, passive_angles, phoenics_config, *args, **kwargs)
-        if previous is not None:
-            if type(previous) is str:
-                try:
-                    obs = pickle.load(open(previous, 'rb'))
-                except:
-                    print(
-                        'failed to load previous observations, which are meant to be a pickle file. Starting fresh.')
-            elif type(previous) is list:
-                if all([type(k) == dict for k in previous]):
-                    obs = previous
-                else:
-                    print('previous observations were not in the correct format (list of dicts). Starting fresh.')
-
-
-
-        if not (type(file_name) == str or file_name == None):
-            raise TequilaException('file_name must be a string, or None. Recieved {}'.format(type(file_name)))
-
-        best = None
-        best_angles = None
-
-        # avoid multiple compilations
-        compiled_objective = compile_objective(objective=objective, backend=self.backend,
-                                               device=self.device,
-                                               samples=self.samples, noise=self.noise)
-
-        if not self.silent:
-            print('phoenics has recieved')
-            print("objective: \n")
-            print(objective)
-            print("noise model : {}".format(self.noise))
-            print("samples     : {}".format(self.samples))
-            print("maxiter     : {}".format(maxiter))
-            print("variables   : {}".format(objective.extract_variables()))
-            print("passive var : {}".format(passive_angles))
-            print('now lets begin')
-        for i in range(0, maxiter):
-            with warnings.catch_warnings():
-                np.testing.suppress_warnings()
-                warnings.simplefilter("ignore")
-                warnings.filterwarnings("ignore", category=FutureWarning)
-
-            precs = bird.recommend(observations=obs)
-
-            runs = []
-            recs = self._process_for_sim(precs, passive_angles=passive_angles)
-
-            start = time.time()
-            for j, rec in enumerate(recs):
-                En = compiled_objective(variables=rec, samples=self.samples, noise=self.noise)
-                runs.append((rec, En))
-                if not self.silent:
-                    if self.print_level > 2:
-                        print("energy = {:+2.8f} , angles=".format(En), rec)
-                    else:
-                        print("energy = {:+2.8f}".format(En))
-            stop = time.time()
-            if not self.silent:
-                print("Quantum Objective evaluations: {}s Wall-Time".format(stop-start))
-
-            for run in runs:
-                angles = run[0]
-                E = run[1]
-                if best is None:
-                    best = E
-                    best_angles = angles
-                else:
-                    if self._minimize:
-                        if E < best:
-                            best = E
-                            best_angles = angles
-                    else:
-                        if E > best:
-                            best = E
-                            best_angles = angles
-
-                if self.save_history:
-                    self.history.energies.append(E)
-                    self.history.angles.append(angles)
-                obs.append(self._process_for_phoenics(angles, E, passive_angles=passive_angles))
-
-        if file_name is not None:
-            with open(file_name, 'wb') as file:
-                pickle.dump(obs, file)
-
-        if not self.silent:
-            print("best energy after {} iterations : {:+2.8f}".format(self.maxiter, best))
-        return PhoenicsResults(energy=best, variables=best_angles, history=self.history, observations=obs,phoenics_instance=bird)
-
-
-def minimize(objective: Objective,
-             maxiter: int = None,
-             samples: int = None,
-             variables: typing.List = None,
-             initial_values: typing.Dict = None,
-             backend: str = None,
-             noise=None,
-             device: str = None,
-             previous: typing.Union[str, list] = None,
-             phoenics_config: typing.Union[str, typing.Dict] = None,
-             file_name: str = None,
-             silent: bool = False,
-             *args,
-             **kwargs) -> PhoenicsResults:
-    """
-
-    Parameters
-    ----------
-    objective: Objective:
-        The tequila objective to optimize
-    initial_values: typing.Dict[typing.Hashable, numbers.Real], optional:
-        Initial values as dictionary of Hashable types (variable keys) and floating point numbers.
-        If given None they will be randomized.
-    variables: typing.List[typing.Hashable], optional:
-         List of Variables to optimize
-    samples: int, optional:
-         samples/shots to take in every run of the quantum circuits (None activates full wavefunction simulation)
-    maxiter: int:
-         how many iterations of phoenics to run.
-         Note that this is NOT identical to the number of times the circuit will run.
-    backend: str, optional:
-         Simulator backend, will be automatically chosen if set to None
-    noise: NoiseModel, optional:
-         a noise model to apply to the circuits of Objective.
-    device: optional:
-        the device from which to (potentially, simulatedly) sample all quantum circuits employed in optimization.
-    previous: optional:
-        Previous phoenics observations. If string, the name of a file from which to load them. Else, a list.
-    phoenics_config: optional:
-        a pre-made phoenics configuration. if str, the name of a file from which to load it; Else, a dictionary.
-        Individual keywords of the 'general' sections can also be passed down as kwargs
-    file_name: str, optional:
-        where to save output to, if save_to_file is True.
-    kwargs: dict:
-        Send down more keywords for single replacements in the phoenics config 'general' section, like e.g. batches=5,
-        boosted=True etc
-    Returns
-    -------
-    PhoenicsResults:
-        the result of an optimization by phoenics.
-    """
-
-    optimizer = OptimizerPhoenics(samples=samples, backend=backend,
-                                  noise=noise,device=device,
-                                  maxiter=maxiter, silent=silent)
-    return optimizer(objective=objective, initial_values=initial_values, variables=variables, previous=previous,
-                     maxiter=maxiter,
-                     phoenics_config=phoenics_config, file_name=file_name, *args, **kwargs)
diff --git a/src/tequila/quantumchemistry/chemistry_tools.py b/src/tequila/quantumchemistry/chemistry_tools.py
index 4d9b2e79..ff1f3ff1 100644
--- a/src/tequila/quantumchemistry/chemistry_tools.py
+++ b/src/tequila/quantumchemistry/chemistry_tools.py
@@ -569,7 +569,7 @@ def _verify_ordering_of(self, trials=100):
         return True
 
     def __init__(self, elems: numpy.ndarray = None, active_indices: list = None, ordering: str = None,
-                 size_full: int = None):
+                 size_full: int = None, verify=False):
         """
         Parameters
         ----------
@@ -611,7 +611,7 @@ def __init__(self, elems: numpy.ndarray = None, active_indices: list = None, ord
         if self.order == 4:
             if ordering is None:
                 ordering = self.identify_ordering()
-            else:
+            elif verify:
                 try: # some RDMs are really sloppy (depends on backend)
                     auto_ordering=self.identify_ordering()
                     if auto_ordering is not ordering:
diff --git a/src/tequila/quantumchemistry/madness_interface.py b/src/tequila/quantumchemistry/madness_interface.py
index e56d58e3..9273046c 100644
--- a/src/tequila/quantumchemistry/madness_interface.py
+++ b/src/tequila/quantumchemistry/madness_interface.py
@@ -122,7 +122,7 @@ def __init__(self, parameters: ParametersQC,
             h = "failed"
             g = "failed"
 
-        if "failed" in h or "failed" in g:
+        if (isinstance(h, str) and "failed" in h) or (isinstance(g, str) and "failed" in g):
             status = "found {}_htensor.npy={}\n".format(name, "failed" not in h)
             status += "found {}_gtensor.npy={}\n".format(name, "failed" not in g)
             try:
diff --git a/src/tequila/quantumchemistry/orbital_optimizer.py b/src/tequila/quantumchemistry/orbital_optimizer.py
index 06b36553..a7bca8d7 100644
--- a/src/tequila/quantumchemistry/orbital_optimizer.py
+++ b/src/tequila/quantumchemistry/orbital_optimizer.py
@@ -118,13 +118,15 @@ def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=N
         print(wrapper)
     if initial_guess is not None:
         if hasattr(initial_guess, "lower"):
-            if "random" in initial_guess.lower():
-                scale = 0.1
+            if "random" or "near_zero" in initial_guess.lower():
+                scale = 1.e-3
+                if "random" in initial_guess.lower():
+                    scale = 1.0
                 loc = 0.0
                 if "scale" in kwargs:
-                    scale = kwargs["scale"]
+                    scale = float(initial_guess.split("scale")[1].split("_")[0].split("=")[1])
                 if "loc" in kwargs:
-                    loc = kwargs["loc"]
+                    loc = float(initial_guess.split("loc")[1].split("_")[0].split("=")[1])
                 initial_guess = numpy.eye(no) + numpy.random.normal(scale=scale, loc=loc, size=no * no).reshape(no, no)
             else:
                 raise Exception("Unknown initial_guess={}".format(initial_guess.lower()))
diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index f2da7b28..9689b042 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -1939,7 +1939,7 @@ def _reset_rdm(rdm):
             self._rdm2 = _assemble_rdm2_spinful(evals_2) if get_rdm2 else self._rdm2
 
         if get_rdm2:
-            rdm2 = NBodyTensor(elems=self.rdm2, ordering="dirac")
+            rdm2 = NBodyTensor(elems=self.rdm2, ordering="dirac", verify=False)
             rdm2.reorder(to=ordering)
             rdm2 = rdm2.elems
             self._rdm2 = rdm2
diff --git a/src/tequila/version.py b/src/tequila/version.py
index f4531a81..a54a7564 100644
--- a/src/tequila/version.py
+++ b/src/tequila/version.py
@@ -1,2 +1,2 @@
-__version__ = "1.9.3"
+__version__ = "1.9.4"
 __author__ = "Tequila Developers "
diff --git a/tests/test_noise_opt.py b/tests/test_noise_opt.py
index 043159ed..a48b9c66 100644
--- a/tests/test_noise_opt.py
+++ b/tests/test_noise_opt.py
@@ -59,20 +59,6 @@ def test_bit_flip_scipy_hessian(p, method):
     result = tq.optimizer_scipy.minimize(objective=O, samples=1, backend=simulator, method=method, noise=NM, tol=1.e-4,
                                          silent=False)
 
-
-@pytest.mark.skipif(len(samplers) == 0, reason="Missing necessary backends")
-@pytest.mark.skipif(not tq.optimizers.has_phoenics, reason="Missing phoenics installation")
-@pytest.mark.parametrize("p", numpy.random.uniform(0.1, .4, 1))
-def test_bit_flip_phoenics(p):
-    simulator = numpy.random.choice(samplers)
-    qubit = 0
-    H = paulis.Qm(qubit)
-    U = gates.Rx(target=qubit, angle=tq.Variable('a'))
-    O = ExpectationValue(U=U, H=H)
-    NM = BitFlip(p, 1)
-    result = tq.optimizers.optimizer_phoenics.minimize(objective=O, maxiter=3, samples=1, backend=simulator, noise=NM)
-
-
 @pytest.mark.skipif(len(samplers) == 0, reason="Missing necessary backends")
 @pytest.mark.skipif(not tq.optimizers.has_gpyopt, reason="Missing gpyopt installation")
 @pytest.mark.parametrize("p", numpy.random.uniform(0.1, .4, 1))
diff --git a/tests/test_phoenics.py b/tests/test_phoenics.py
deleted file mode 100644
index 0cea65eb..00000000
--- a/tests/test_phoenics.py
+++ /dev/null
@@ -1,68 +0,0 @@
-import pytest, numpy
-import tequila as tq
-import multiprocessing as mp
-
-# Get QC backends for parametrized testing
-import select_backends
-simulators = select_backends.get()
-samplers = select_backends.get(sampler=True)
-
-has_phoenics = 'phoenics' in tq.INSTALLED_OPTIMIZERS
-
-@pytest.mark.dependencies
-def test_dependencies():
-    assert 'phoenics' in tq.INSTALLED_OPTIMIZERS
-
-
-@pytest.mark.skipif(condition=not has_phoenics, reason="you don't have phoenics")
-@pytest.mark.parametrize("simulator", simulators)
-def test_execution(simulator):
-    U = tq.gates.Rz(angle="a", target=0) \
-        + tq.gates.X(target=2) \
-        + tq.gates.Ry(angle="b", target=1, control=2) \
-        + tq.gates.Trotterized(angles=["c", "d"],
-                               generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2) \
-        + tq.gates.Trotterized(angles=[1.0, 2.0],
-                               generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2) \
-        + tq.gates.ExpPauli(angle="a", paulistring="X(0)Y(1)Z(2)")
-
-    H = 1.0 * tq.paulis.X(0) + 2.0 * tq.paulis.Y(1) + 3.0 * tq.paulis.Z(2)
-    O = tq.ExpectationValue(U=U, H=H)
-    result = tq.minimize(method="phoenics", objective=O, maxiter=1, backend=simulator)
-
-
-@pytest.mark.skipif(condition=not has_phoenics, reason="you don't have phoenics")
-@pytest.mark.parametrize("simulator", samplers)
-def test_execution_shot(simulator):
-    U = tq.gates.Rz(angle="a", target=0) \
-        + tq.gates.X(target=2) \
-        + tq.gates.Ry(angle="b", target=1, control=2) \
-        + tq.gates.Trotterized(angles=["c", "d"],
-                               generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2) \
-        + tq.gates.Trotterized(angles=[1.0, 2.0],
-                               generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2) \
-        + tq.gates.ExpPauli(angle="a", paulistring="X(0)Y(1)Z(2)")
-    H = 1.0 * tq.paulis.X(0) + 2.0 * tq.paulis.Y(1) + 3.0 * tq.paulis.Z(2)
-    O = tq.ExpectationValue(U=U, H=H)
-    mi = 2
-    result = tq.minimize(method="phoenics", objective=O, maxiter=mi, backend=simulator)
-
-
-@pytest.mark.skipif(condition=not has_phoenics, reason="you don't have phoenics")
-@pytest.mark.parametrize("simulator", simulators)
-def test_one_qubit_wfn(simulator):
-    U = tq.gates.Trotterized(angles=["a"], steps=1, generators=[tq.paulis.Y(0)])
-    H = tq.paulis.X(0)
-    O = tq.ExpectationValue(U=U, H=H)
-    result = tq.minimize(method="phoenics", objective=O, maxiter=8, backend=simulator)
-    assert (numpy.isclose(result.energy, -1.0, atol=1.e-2))
-
-
-@pytest.mark.skipif(condition=not has_phoenics, reason="you don't have phoenics")
-@pytest.mark.parametrize("simulator", samplers)
-def test_one_qubit_shot(simulator):
-    U = tq.gates.Trotterized(angles=["a"], steps=1, generators=[tq.paulis.Y(0)])
-    H = tq.paulis.X(0)
-    O = tq.ExpectationValue(U=U, H=H)
-    result = tq.minimize(method="phoenics", objective=O, maxiter=3, backend=simulator, samples=10000)
-    assert (numpy.isclose(result.energy, -1.0, atol=1.e-1))

From a752e8c2a37618a5105baa391053cf8ba4a7d93c Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 18 Mar 2024 17:32:02 +0100
Subject: [PATCH 02/29] Update qasm.py (#334)

Fixes issue #332
---
 src/tequila/circuit/qasm.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/tequila/circuit/qasm.py b/src/tequila/circuit/qasm.py
index 65878f76..7eb567ee 100644
--- a/src/tequila/circuit/qasm.py
+++ b/src/tequila/circuit/qasm.py
@@ -362,7 +362,7 @@ def parse_command(command: str, custom_gates_map: Dict[str, QCircuit], qregister
                   control=get_qregister(args[0], qregisters),
                   target=get_qregister(args[1], qregisters))
     if name in ("s", "t", "sdg", "tdg"):
-        g = gates.Phase(pi / (2 if name.startswith("s") else 4),
+        g = gates.Phase(angle=pi / (2 if name.startswith("s") else 4),
                      control=None,
                      target=get_qregister(args[0], qregisters))
         if name.find("dg") != -1:

From 3ded2f39227e162942b0596ae2736f2c2adf23be Mon Sep 17 00:00:00 2001
From: davibinco <59845652+davibinco@users.noreply.github.com>
Date: Wed, 20 Mar 2024 14:04:44 +0100
Subject: [PATCH 03/29]  added methods to create annihilation, creation, sz,
 sp, sm and s2 operators in qubit representation (#336)

* added methods to create annihilation, creation, sz, sp, sm and s2 operators in qubit representation
* Added documentation and make number operator method
---
 src/tequila/quantumchemistry/qc_base.py | 59 +++++++++++++++++++++++++
 1 file changed, 59 insertions(+)

diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index 9689b042..669a3b37 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -645,6 +645,65 @@ def n_electrons(self) -> int:
         """
         return 2 * len(self.integral_manager.active_reference_orbitals)
 
+    def make_annihilation_op(self, orbital, coefficient=1.0):
+        """
+        Compute annihilation operator on orbital=orbital in qubit representation
+        """
+        assert orbital<=self.n_orbitals*2
+        aop = openfermion.ops.FermionOperator(f'{orbital}', coefficient)
+        return self.transformation(aop)
+
+    def make_creation_op(self, orbital, coefficient=1.0):
+        """
+        Compute creation operator on orbital=orbital in qubit representation
+        """
+        assert orbital<=self.n_orbitals*2
+        cop = openfermion.ops.FermionOperator(f'{orbital}^', coefficient)
+        return self.transformation(cop)
+
+    def make_number_op(self, orbital):
+        """
+        Compute number operator on orbital=orbital in qubit representation
+        """
+        num_op = self.make_creation_op(orbital) * self.make_annihilation_op(orbital)
+        return num_op
+    
+    def make_sz_op(self):
+        """
+        Compute the spin_z operator of the molecule in qubit representation
+        """
+        sz = QubitHamiltonian()
+        for i in range(0, self.n_orbitals * 2, 2):
+            one = 0.5 * self.make_creation_op(i) * self.make_annihilation_op(i)
+            two = 0.5 * self.make_creation_op(i+1) * self.make_annihilation_op(i+1)
+            sz += (one - two)
+        return sz
+
+    def make_sp_op(self):
+        """
+        Compute the spin+ operator of the molecule in qubit representation
+        """
+        sp = QubitHamiltonian()
+        for i in range(self.n_orbitals):
+            sp += self.make_creation_op(i*2) * self.make_annihilation_op(i*2 + 1)
+        return sp
+
+    def make_sm_op(self):
+        """
+        Compute the spin- operator of the molecule in qubit representation
+        """
+        sm = QubitHamiltonian()
+        for i in range(self.n_orbitals):
+            sm += self.make_creation_op(i*2 + 1) * self.make_annihilation_op(i*2)
+        return sm
+
+    def make_s2_op(self):
+        """
+        Compute the spin^2 operator of the molecule in qubit representation
+        """
+        s2_op = self.make_sm_op() * self.make_sp_op() + self.make_sz_op() * (self.make_sz_op() + 1)
+        return s2_op
+
     def make_hamiltonian(self, *args, **kwargs) -> QubitHamiltonian:
         """
         Parameters

From 423c8e7bf147565c5a1bb457165a1572f322a1e1 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Wed, 20 Mar 2024 22:41:21 +0100
Subject: [PATCH 04/29] Update qasm.py (#335)

* Update qasm.py
* Update ci_chemistry_psi4.yml (psi4 --> conda-forge channel)
* Update README.md
---
 .github/workflows/ci_chemistry_psi4.yml |  4 +--
 README.md                               |  4 +--
 src/tequila/circuit/qasm.py             | 34 +++++++++++++++++--------
 3 files changed, 28 insertions(+), 14 deletions(-)

diff --git a/.github/workflows/ci_chemistry_psi4.yml b/.github/workflows/ci_chemistry_psi4.yml
index ef9c1e3a..98481972 100644
--- a/.github/workflows/ci_chemistry_psi4.yml
+++ b/.github/workflows/ci_chemistry_psi4.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.8]
+        python-version: [3.9]
 
     steps:
     - uses: actions/checkout@v2
@@ -50,7 +50,7 @@ jobs:
         source $HOME/.bashrc
         source $CONDABASE/bin/activate
         conda activate test_psi4
-        conda install psi4 python=3.8 -c psi4 
+        conda install psi4 python=3.9 -c conda-forge
         python -m pip install --upgrade pip
         python -m pip install -r requirements.txt
         python -m pip install -e .
diff --git a/README.md b/README.md
index 3cec3281..b415c3f3 100644
--- a/README.md
+++ b/README.md
@@ -14,7 +14,7 @@ Tequila can execute the underlying quantum expectation values on state of the ar
 - [talks and slides](https://kottmanj.github.io/talks_and_material/)  
 
 # Installation
-Recommended Python version is 3.8-3.9.   
+Recommended Python version is 3.9 - 3.10.   
 Tequila supports linux, osx and windows. However, not all optional dependencies are supported on windows.  
 
 ## Install from PyPi
@@ -267,7 +267,7 @@ Currently supported
 ### [Psi4](https://github.com/psi4/psi4).
 In a conda environment this can be installed with
 ```bash
-conda install psi4 -c psi4
+conda install psi4 -c conda-forge
 ```
 Here is a small [tutorial](https://nbviewer.org/github/tequilahub/tequila-tutorials/blob/main/chemistry/ChemistryModule.ipynb) that illustrates the usage.
 
diff --git a/src/tequila/circuit/qasm.py b/src/tequila/circuit/qasm.py
index 7eb567ee..16a35b57 100644
--- a/src/tequila/circuit/qasm.py
+++ b/src/tequila/circuit/qasm.py
@@ -313,19 +313,33 @@ def parse_command(command: str, custom_gates_map: Dict[str, QCircuit], qregister
         return apply_custom_gate(custom_circuit=custom_circuit, qregisters_values=qregisters_values)
 
     if name in ("x", "y", "z", "h", "cx", "cy", "cz", "ch"):
-        return QCircuit.wrap_gate(gates.impl.QGateImpl(name=(name[1] if name[0] == 'c' else name).upper(),
-                                            control=get_qregister(args[0], qregisters) if name[0] == 'c' else None,
-                                            target=get_qregister(args[1 if name[0] == 'c' else 0], qregisters)))
+        target = get_qregister(args[0], qregisters)
+        control = None
+        if name[0].lower() == 'c':
+            control = get_qregister(args[0], qregisters)
+            target = get_qregister(args[1], qregisters)
+            name = name[1]
+        G = getattr(gates, name.upper())
+        return G(control=control, target=target)
+
     if name in ("ccx", "ccy", "ccz"):
-        return QCircuit.wrap_gate(gates.impl.QGateImpl(name=name.upper()[2],
-                                            control=[get_qregister(args[0], qregisters), get_qregister(args[1], qregisters)],
-                                            target=get_qregister(args[2], qregisters)))
+        G = getattr(gates, name[2].upper())
+        control = [get_qregister(args[0], qregisters), get_qregister(args[1], qregisters)]
+        target = get_qregister(args[2], qregisters)
+        return G(control=control, target=target)
+
     if name.startswith("rx(") or name.startswith("ry(") or name.startswith("rz(") or \
         name.startswith("crx(") or name.startswith("cry(") or name.startswith("crz("):
-        return QCircuit.wrap_gate(gates.impl.RotationGateImpl(axis=name[2 if name[0] == 'c' else 1],
-                                                   angle=get_angle(name)[0],
-                                                   control=get_qregister(args[0], qregisters) if name[0] == 'c' else None,
-                                                   target=get_qregister(args[1 if name[0] == 'c' else 0], qregisters)))
+        angle = get_angle(name)[0]
+        i = name.find('(')
+        name = name[0:i]
+        name = name.upper()
+        name = [x for x in name]
+        name[-1] = name[-1].lower()
+        name = "".join(name)
+        G = getattr(gates, name)
+        return G(angle=angle,control=get_qregister(args[0], qregisters) if name[0] == 'C' else None,target=get_qregister(args[1 if name[0] == 'C' else 0], qregisters))
+            
     if name.startswith("U("):
         angles = get_angle(name)
         return gates.U(theta=angles[0], phi=angles[1], lambd=angles[2],

From a3d532b7f041d4df932f6104f5e1622cbea6f05e Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 25 Mar 2024 10:56:11 +0100
Subject: [PATCH 05/29] Update README.md

---
 README.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/README.md b/README.md
index 3cec3281..f51e8b75 100644
--- a/README.md
+++ b/README.md
@@ -240,6 +240,10 @@ A.G. Cadavid, I. Montalban, A. Dalal, E. Solano, N.N. Hegade
 Efficient DCQO Algorithm within the Impulse Regime for Portfolio Optimization  
 [arxiv:2308.15475](https://arxiv.org/abs/2308.15475)  
 
+P.W.K. Jensen, E.R. Kjellgren, P. Reinholdt, K.M. Ziems, S. Coriani, J. Kongsted, S. Sauer  
+Quantum Equation of Motion with Orbital Optimization for Computing Molecular Properties in Near-Term Quantum Computing  
+[arxiv:2312.12386](https://arxiv.org/abs/2312.12386)  
+
 Let us know, if you want your research project and/or tutorial to be included in this list!
 
 # Dependencies

From a43ce8b6eff5f54e0480ebcd1239b8c6b154755e Mon Sep 17 00:00:00 2001
From: davibinco <59845652+davibinco@users.noreply.github.com>
Date: Tue, 26 Mar 2024 14:50:23 +0100
Subject: [PATCH 06/29] fix doc in fermionic operators for molecule (#337)

* fix doc in fermionic operators for molecule, now spin-orbital is specified
* add doc for spin-orbital order
---
 src/tequila/quantumchemistry/qc_base.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index 669a3b37..9f089cc4 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -647,7 +647,8 @@ def n_electrons(self) -> int:
 
     def make_annihilation_op(self, orbital, coefficient=1.0):
         """
-        Compute annihilation operator on orbital=orbital in qubit representation
+        Compute annihilation operator on spin-orbital in qubit representation
+        Spin-orbital order is always (up,down,up,down,...)
         """
         assert orbital<=self.n_orbitals*2
         aop = openfermion.ops.FermionOperator(f'{orbital}', coefficient)
@@ -655,7 +656,8 @@ def make_annihilation_op(self, orbital, coefficient=1.0):
 
     def make_creation_op(self, orbital, coefficient=1.0):
         """
-        Compute creation operator on orbital=orbital in qubit representation
+        Compute creation operator on spin-orbital in qubit representation
+        Spin-orbital order is always (up,down,up,down,...)
         """
         assert orbital<=self.n_orbitals*2
         cop = openfermion.ops.FermionOperator(f'{orbital}^', coefficient)
@@ -663,7 +665,8 @@ def make_creation_op(self, orbital, coefficient=1.0):
 
     def make_number_op(self, orbital):
         """
-        Compute number operator on orbital=orbital in qubit representation
+        Compute number operator on spin-orbital in qubit representation
+        Spin-orbital order is always (up,down,up,down,...)
         """
         num_op = self.make_creation_op(orbital) * self.make_annihilation_op(orbital)
         return num_op

From a547f3c046c417f05687f30a7b04ec28362cd356 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 1 Apr 2024 23:45:59 +0200
Subject: [PATCH 07/29] Update README.md

---
 README.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/README.md b/README.md
index f51e8b75..9c350865 100644
--- a/README.md
+++ b/README.md
@@ -240,6 +240,10 @@ A.G. Cadavid, I. Montalban, A. Dalal, E. Solano, N.N. Hegade
 Efficient DCQO Algorithm within the Impulse Regime for Portfolio Optimization  
 [arxiv:2308.15475](https://arxiv.org/abs/2308.15475)  
 
+A. Anand, K. Brown  
+Hamiltonians, groups, graphs and ansätze  
+[arxiv:2312.17146](https://arxiv.org/abs/2312.17146)  
+
 P.W.K. Jensen, E.R. Kjellgren, P. Reinholdt, K.M. Ziems, S. Coriani, J. Kongsted, S. Sauer  
 Quantum Equation of Motion with Orbital Optimization for Computing Molecular Properties in Near-Term Quantum Computing  
 [arxiv:2312.12386](https://arxiv.org/abs/2312.12386)  

From f4df0fd9b2e2a4d2987192c662e6e3b5b7cd8174 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Fri, 5 Apr 2024 09:29:09 +0200
Subject: [PATCH 08/29] orbital names are better tracked through
 transformations (#338)

* orbital names are better tracked through transformations
* keep basis info when re-initializing with mol.from_tequila
* orbital-optimization gives back molecule in the original basis
---
 .../quantumchemistry/chemistry_tools.py       | 28 ++++++----
 .../quantumchemistry/orbital_optimizer.py     | 14 +++--
 .../quantumchemistry/pyscf_interface.py       |  1 +
 src/tequila/quantumchemistry/qc_base.py       | 51 +++++++++++++------
 4 files changed, 65 insertions(+), 29 deletions(-)

diff --git a/src/tequila/quantumchemistry/chemistry_tools.py b/src/tequila/quantumchemistry/chemistry_tools.py
index ff1f3ff1..8cebef81 100644
--- a/src/tequila/quantumchemistry/chemistry_tools.py
+++ b/src/tequila/quantumchemistry/chemistry_tools.py
@@ -804,20 +804,18 @@ class IntegralManager:
     _one_body_integrals: numpy.ndarray = None
     _two_body_integrals: NBodyTensor = None
     _constant_term: float = None
-    _basis_type: str = "unknown"
     _basis_name: str = "unknown"
     _orbital_type: str = "unknown" # e.g. "HF", "PNO", "native"
     _orbital_coefficients: numpy.ndarray = None
     _active_space: ActiveSpaceData = None
     _orbitals: typing.List[OrbitalData] = None
 
-    def __init__(self, one_body_integrals, two_body_integrals, basis_type="custom",
+    def __init__(self, one_body_integrals, two_body_integrals,
                  basis_name="unknown", orbital_type="unknown",
                  constant_term=0.0, orbital_coefficients=None, active_space=None, overlap_integrals=None, orbitals=None, *args, **kwargs):
         self._one_body_integrals = one_body_integrals
         self._two_body_integrals = two_body_integrals
         self._constant_term = constant_term
-        self._basis_type = basis_type
         self._basis_name = basis_name
         self._orbital_type = orbital_type
 
@@ -956,9 +954,16 @@ def transform_to_native_orbitals(self):
         """
         c = self.get_orthonormalized_orbital_coefficients()
         self.orbital_coefficients=c
-        self._orbital_type="orthonormalized-{}-basis".format(self._orbital_type)
+        self._orbital_type="orthonormalized-{}-basis".format(self._basis_name)
 
-    def transform_orbitals(self, U):
+    def is_unitary(self, U):
+        if len(U.shape) != 2: return False
+        if U.shape[0] != U.shape[1]: return False
+        test = (U.conj().T).dot(U) - numpy.eye(U.shape[0])
+        if not numpy.isclose(numpy.linalg.norm(test), 0.0): return False
+        return True
+
+    def transform_orbitals(self, U, name=None):
         """
         Transform orbitals
         Parameters
@@ -969,10 +974,12 @@ def transform_orbitals(self, U):
         -------
         updates the structure with new orbitals: c = cU
         """
-        c = self.orbital_coefficients
-        c = numpy.einsum("ix, xj -> ij", c, U, optimize="greedy")
-        self.orbital_coefficients = c
-        self._orbital_type += "-transformed"
+        assert self.is_unitary(U)
+        self.orbital_coefficients = numpy.einsum("ix, xj -> ij", self.orbital_coefficients, U, optimize="greedy")
+        if name is None:
+            self._orbital_type += "-transformed"
+        else:
+            self._orbital_type = name
 
     def get_integrals(self, orbital_coefficients=None, ordering="openfermion", ignore_active_space=False, *args, **kwargs):
         """
@@ -1001,7 +1008,9 @@ def get_integrals(self, orbital_coefficients=None, ordering="openfermion", ignor
             active_integrals = get_active_space_integrals(one_body_integrals=h, two_body_integrals=g,
                                                           occupied_indices=self._active_space.frozen_reference_orbitals,
                                                           active_indices=self._active_space.active_orbitals)
+
             c = active_integrals[0] + c
+
             h = active_integrals[1]
             g = NBodyTensor(elems=active_integrals[2], ordering="openfermion")
         g.reorder(to=ordering)
@@ -1070,7 +1079,6 @@ def __str__(self):
         return result
 
     def print_basis_info(self, *args, **kwargs) -> None:
-        print("{:15} : {}".format("basis_type", self._basis_type), *args, **kwargs)
         print("{:15} : {}".format("basis_name", self._basis_name), *args, **kwargs)
         print("{:15} : {}".format("orbital_type", self._orbital_type), *args, **kwargs)
         print("{:15} : {}".format("orthogonal", self.basis_is_orthogonal()), *args, **kwargs)
diff --git a/src/tequila/quantumchemistry/orbital_optimizer.py b/src/tequila/quantumchemistry/orbital_optimizer.py
index a7bca8d7..5972134a 100644
--- a/src/tequila/quantumchemistry/orbital_optimizer.py
+++ b/src/tequila/quantumchemistry/orbital_optimizer.py
@@ -37,7 +37,7 @@ def __call__(self, local_data, *args, **kwargs):
         self.iterations += 1
 
 def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=None, silent=False,
-                      vqe_solver_arguments=None, initial_guess=None, return_mcscf=False, use_hcb=False, molecule_factory=None, molecule_arguments=None, *args, **kwargs):
+                      vqe_solver_arguments=None, initial_guess=None, return_mcscf=False, use_hcb=False, molecule_factory=None, molecule_arguments=None, restrict_to_active_space=True, *args, **kwargs):
     """
 
     Parameters
@@ -78,7 +78,12 @@ def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=N
     if pyscf_arguments is None:
         pyscf_arguments = {"max_cycle_macro": 10, "max_cycle_micro": 3}
     no = molecule.n_orbitals
-    pyscf_molecule = QuantumChemistryPySCF.from_tequila(molecule=molecule, transformation=molecule.transformation)
+
+    if not isinstance(molecule, QuantumChemistryPySCF):
+        pyscf_molecule = QuantumChemistryPySCF.from_tequila(molecule=molecule, transformation=molecule.transformation)
+    else:
+        pyscf_molecule = molecule
+
     mf = pyscf_molecule._get_hf()
     result=OptimizeOrbitalsResult()
     mc = mcscf.CASSCF(mf, pyscf_molecule.n_orbitals, pyscf_molecule.n_electrons)
@@ -140,10 +145,11 @@ def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=N
         mc.kernel()
     # make new molecule
 
-    transformed_molecule = pyscf_molecule.transform_orbitals(orbital_coefficients=mc.mo_coeff)
+    mo_coeff = mc.mo_coeff
+    transformed_molecule = pyscf_molecule.transform_orbitals(orbital_coefficients=mo_coeff, name="optimized")
     result.molecule=transformed_molecule
     result.old_molecule=molecule
-    result.mo_coeff=mc.mo_coeff
+    result.mo_coeff=mo_coeff
     result.energy=mc.e_tot
     
     if return_mcscf:
diff --git a/src/tequila/quantumchemistry/pyscf_interface.py b/src/tequila/quantumchemistry/pyscf_interface.py
index bed3cc5d..69b681f9 100644
--- a/src/tequila/quantumchemistry/pyscf_interface.py
+++ b/src/tequila/quantumchemistry/pyscf_interface.py
@@ -75,6 +75,7 @@ def __init__(self, parameters: ParametersQC,
             kwargs["two_body_integrals"] = g_ao
             kwargs["one_body_integrals"] = h_ao
             kwargs["orbital_coefficients"] = mo_coeff
+            kwargs["orbital_type"] = "hf"
 
             if "nuclear_repulsion" not in kwargs:
                 kwargs["nuclear_repulsion"] = mol.energy_nuc()
diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index 9f089cc4..d8a6882b 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -94,7 +94,7 @@ def __init__(self, parameters: ParametersQC,
         else:
             self.integral_manager = self.initialize_integral_manager(active_orbitals=active_orbitals,
                                                                      reference_orbitals=reference_orbitals,
-                                                                     orbitals=orbitals, frozen_orbitals=frozen_orbitals, orbital_type=orbital_type, *args,
+                                                                     orbitals=orbitals, frozen_orbitals=frozen_orbitals, orbital_type=orbital_type, basis_name=self.parameters.basis_set, *args,
                                                                      **kwargs)
         
         if orbital_type is not None and orbital_type.lower() == "native":
@@ -109,15 +109,23 @@ def __init__(self, parameters: ParametersQC,
 
     @classmethod
     def from_tequila(cls, molecule, transformation=None, *args, **kwargs):
-        c, h1, h2 = molecule.get_integrals()
+        c = molecule.integral_manager.constant_term
+        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 transformation is None:
             transformation = molecule.transformation
+        parameters = molecule.parameters
         return cls(nuclear_repulsion=c,
                    one_body_integrals=h1,
                    two_body_integrals=h2,
-                   n_electrons=molecule.n_electrons,
+                   overlap_integrals = S,
+                   orbital_coefficients = molecule.integral_manager.orbital_coefficients,
+                   active_orbitals= active_orbitals,
                    transformation=transformation,
-                   parameters=molecule.parameters, *args, **kwargs)
+                   orbital_type=molecule.integral_manager._orbital_type,
+                   parameters=parameters, *args, **kwargs)
 
     def supports_ucc(self):
         """
@@ -543,11 +551,13 @@ def initialize_integral_manager(self, *args, **kwargs):
 
         return manager
 
-    def transform_orbitals(self, orbital_coefficients, *args, **kwargs):
+    def transform_orbitals(self, orbital_coefficients, ignore_active_space=False, name=None, *args, **kwargs):
         """
         Parameters
         ----------
-        orbital_coefficients: second index is new orbital indes, first is old orbital index (summed over)
+        orbital_coefficients: second index is new orbital indes, first is old orbital index (summed over), indices are assumed to be defined on the active space
+        ignore_active_space: if true orbital_coefficients are not assumed to be given in the active space
+        name: str, name the new orbitals
         args
         kwargs
 
@@ -556,9 +566,20 @@ def transform_orbitals(self, orbital_coefficients, *args, **kwargs):
         New molecule with transformed orbitals
         """
 
+        U = numpy.eye(self.integral_manager.orbital_coefficients.shape[0])
+        # mo_coeff by default only acts on the active space
+        active_indices = [o.idx_total for o in self.integral_manager.active_orbitals]
+
+        if ignore_active_space:
+            U = orbital_coefficients
+        else:
+            for kk,k in enumerate(active_indices):
+                for ll,l in enumerate(active_indices):
+                    U[k][l] = orbital_coefficients[kk][ll]
+
         # can not be an instance of a specific backend (otherwise we get inconsistencies with classical methods in the backend)
         integral_manager = copy.deepcopy(self.integral_manager)
-        integral_manager.transform_orbitals(U=orbital_coefficients)
+        integral_manager.transform_orbitals(U=U, name=name)
         result = QuantumChemistryBase(parameters=self.parameters, integral_manager=integral_manager, transformation=self.transformation)
         return result
     
@@ -583,7 +604,7 @@ def use_native_orbitals(self, inplace=False):
         else:
             integral_manager = copy.deepcopy(self.integral_manager)
             integral_manager.transform_to_native_orbitals()
-            result = QuantumChemistryBase(parameters=self.parameters, integral_manager=integral_manager, orbital_type="native", transformation=self.transformation)
+            result = QuantumChemistryBase(parameters=self.parameters, integral_manager=integral_manager, transformation=self.transformation)
             return result
 
 
@@ -867,13 +888,13 @@ def hcb_to_me(self, U=None, condensed=False):
         """
         if U is None:
             U = QCircuit()
-
-        # consistency
-        consistency = [x < self.n_orbitals for x in U.qubits]
-        if not all(consistency):
-            warnings.warn(
-                "hcb_to_me: given circuit is not defined on the first {} qubits. Is this a HCB circuit?".format(
-                    self.n_orbitals))
+        else:
+            ups = [self.transformation.up(i.idx) for i in self.orbitals]
+            consistency = [x in ups for x in U.qubits]
+            if not all(consistency):
+                warnings.warn(
+                    "hcb_to_me: given circuit is not defined on all first {} qubits. Is this a HCB circuit?".format(
+                        self.n_orbitals))
 
         # map to alpha qubits
         if condensed:

From ebd38cbed628306d2555916d43e6f689ab87517b Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Fri, 5 Apr 2024 10:43:04 +0200
Subject: [PATCH 09/29] more basis info (#339)

---
 src/tequila/quantumchemistry/__init__.py        | 2 ++
 src/tequila/quantumchemistry/chemistry_tools.py | 9 ++++++---
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/tequila/quantumchemistry/__init__.py b/src/tequila/quantumchemistry/__init__.py
index f566522f..e04398c7 100644
--- a/src/tequila/quantumchemistry/__init__.py
+++ b/src/tequila/quantumchemistry/__init__.py
@@ -95,6 +95,8 @@ def Molecule(geometry: str = None,
     if backend is None:
         if basis_set is None or basis_set.lower() in ["madness", "mra", "pno"]:
             backend = "madness"
+            basis_set = "mra"
+            parameters.basis_set = basis_set
             if orbital_type is not None and orbital_type.lower() not in ["pno", "mra-pno"]:
                 warnings.warn("only PNOs supported as orbital_type without basis set. Setting to pno - You gave={}".format(orbital_type), TequilaWarning)
             orbital_type = "pno"
diff --git a/src/tequila/quantumchemistry/chemistry_tools.py b/src/tequila/quantumchemistry/chemistry_tools.py
index 8cebef81..82b72c0d 100644
--- a/src/tequila/quantumchemistry/chemistry_tools.py
+++ b/src/tequila/quantumchemistry/chemistry_tools.py
@@ -1078,13 +1078,16 @@ def __str__(self):
             result += str(x) + "\n"
         return result
 
-    def print_basis_info(self, *args, **kwargs) -> None:
+    def print_basis_info(self, print_coefficients=True, *args, **kwargs) -> None:
         print("{:15} : {}".format("basis_name", self._basis_name), *args, **kwargs)
         print("{:15} : {}".format("orbital_type", self._orbital_type), *args, **kwargs)
-        print("{:15} : {}".format("orthogonal", self.basis_is_orthogonal()), *args, **kwargs)
-        print("{:15} : {}".format("functions", self.one_body_integrals.shape[0]), *args, **kwargs)
+        print("{:15} : {}".format("orthogonal basis", self.basis_is_orthogonal()), *args, **kwargs)
+        print("{:15} : {}".format("basis functions", self.one_body_integrals.shape[0]), *args, **kwargs)
+        print("{:15} : {}".format("active orbitals", [o.idx_total for o in self.active_orbitals]), *args, **kwargs)
         print("{:15} : {}".format("reference", [x.idx_total for x in self.reference_orbitals]), *args, **kwargs)
 
+        if not print_coefficients: return
+
         print("Current Orbitals", *args, **kwargs)
         for i,x in enumerate(self.orbitals):
             print(x, *args, **kwargs)

From 21c443abc379a3c7a4bff77db3711419e420f5e7 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Fri, 5 Apr 2024 10:57:02 +0200
Subject: [PATCH 10/29] Update version.py

---
 src/tequila/version.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/tequila/version.py b/src/tequila/version.py
index a54a7564..cb9458d9 100644
--- a/src/tequila/version.py
+++ b/src/tequila/version.py
@@ -1,2 +1,2 @@
-__version__ = "1.9.4"
+__version__ = "1.9.6-dev"
 __author__ = "Tequila Developers "

From 086ff99b9afd9f809f6d04788bb00f37b56a3d3b Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Fri, 5 Apr 2024 10:57:20 +0200
Subject: [PATCH 11/29] Update version.py

---
 src/tequila/version.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/tequila/version.py b/src/tequila/version.py
index cb9458d9..09ed92e6 100644
--- a/src/tequila/version.py
+++ b/src/tequila/version.py
@@ -1,2 +1,2 @@
-__version__ = "1.9.6-dev"
+__version__ = "1.9.5"
 __author__ = "Tequila Developers "

From 2577c7692ded5a1a7e4144a9c57c650054a79e19 Mon Sep 17 00:00:00 2001
From: JdelArco98 <javierdelarcosantos@gmail.com>
Date: Thu, 18 Apr 2024 11:50:23 +0200
Subject: [PATCH 12/29] Fixing Pyscf_mol.from_tq() (#342)

* Fixing Pyscf_mol.from_tq()
---
 README.md                               | 8 ++++++++
 src/tequila/quantumchemistry/qc_base.py | 3 ++-
 src/tequila/version.py                  | 2 +-
 3 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index b415c3f3..37a0275e 100644
--- a/README.md
+++ b/README.md
@@ -240,6 +240,14 @@ A.G. Cadavid, I. Montalban, A. Dalal, E. Solano, N.N. Hegade
 Efficient DCQO Algorithm within the Impulse Regime for Portfolio Optimization  
 [arxiv:2308.15475](https://arxiv.org/abs/2308.15475)  
 
+A. Anand, K. Brown  
+Hamiltonians, groups, graphs and ansätze  
+[arxiv:2312.17146](https://arxiv.org/abs/2312.17146)  
+
+P.W.K. Jensen, E.R. Kjellgren, P. Reinholdt, K.M. Ziems, S. Coriani, J. Kongsted, S. Sauer  
+Quantum Equation of Motion with Orbital Optimization for Computing Molecular Properties in Near-Term Quantum Computing  
+[arxiv:2312.12386](https://arxiv.org/abs/2312.12386)  
+
 Let us know, if you want your research project and/or tutorial to be included in this list!
 
 # Dependencies
diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index d8a6882b..93d51c0e 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -125,7 +125,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):
         """
diff --git a/src/tequila/version.py b/src/tequila/version.py
index cb9458d9..862db8be 100644
--- a/src/tequila/version.py
+++ b/src/tequila/version.py
@@ -1,2 +1,2 @@
-__version__ = "1.9.6-dev"
+__version__ = "1.9.6dev"
 __author__ = "Tequila Developers "

From a1651f6803d2e739d3cf8a778a4d03c52b1f802d Mon Sep 17 00:00:00 2001
From: "Phillip W. K. Jensen" <46899159+pwkj@users.noreply.github.com>
Date: Wed, 24 Apr 2024 09:57:10 +0200
Subject: [PATCH 13/29] RDMs transformation (#341)

* rdms transformation
* rdm_trafo default from False to None
---
 src/tequila/quantumchemistry/qc_base.py | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index 93d51c0e..305dc242 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -1722,7 +1722,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):
         """
         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.
@@ -1750,6 +1751,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
         -------
@@ -1994,8 +1998,13 @@ 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:
+            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

From c22991c8ff449e23a3bc5097f1f45b6ac8b3fb82 Mon Sep 17 00:00:00 2001
From: JdelArco98 <francisco.del.arco.santos@uni-a.de>
Date: Fri, 26 Apr 2024 20:23:48 +0200
Subject: [PATCH 14/29] fix U.export_to() personalized (#344)

---
 src/tequila/circuit/qpic.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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)))

From 6595f6926af0d7b0c3c312bd92e9099f8881db6d Mon Sep 17 00:00:00 2001
From: Deadly Artist <44524817+DeadlyArtist@users.noreply.github.com>
Date: Fri, 26 Apr 2024 20:25:06 +0200
Subject: [PATCH 15/29] Add Givens Methods (#343)

* add automatic construction of orbital rotation circuits
---
 src/tequila/quantumchemistry/qc_base.py | 194 +++++++++++++++++++++++-
 src/tequila/tools/random_generators.py  |  17 +++
 tests/test_chemistry.py                 |  54 +++++++
 3 files changed, 262 insertions(+), 3 deletions(-)

diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index 305dc242..560c6b91 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
@@ -2137,6 +2136,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()
@@ -2157,3 +2208,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/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/tests/test_chemistry.py b/tests/test_chemistry.py
index 86b2dc33..c4369442 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
@@ -725,3 +727,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)

From a98c38ed2b8638510f7e0f0c3737cb230cf02267 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 3 Jun 2024 11:37:48 +0200
Subject: [PATCH 16/29] updated github runners to 3.9 3.10

---
 .github/workflows/ci_backends.yml          |  2 +-
 .github/workflows/ci_basic.yml             |  2 +-
 .github/workflows/ci_basic_autograd.yml    |  8 ++++----
 .github/workflows/ci_chemistry_madness.yml |  2 +-
 .github/workflows/ci_chemistry_pyscf.yml   |  2 +-
 .github/workflows/ci_conda_madness.yml     |  2 +-
 .github/workflows/ci_ml.yml                |  2 +-
 .github/workflows/ci_optimizers.yml        |  2 +-
 .github/workflows/ci_pyquil.yml            |  2 +-
 .github/workflows/pypi.yml                 |  2 +-
 .github/workflows/release.yml              | 22 ----------------------
 11 files changed, 13 insertions(+), 35 deletions(-)
 delete mode 100644 .github/workflows/release.yml

diff --git a/.github/workflows/ci_backends.yml b/.github/workflows/ci_backends.yml
index c9c5852c..9ae5d9ab 100644
--- a/.github/workflows/ci_backends.yml
+++ b/.github/workflows/ci_backends.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7, 3.8]
+        python-version: ['3.9']
 
     steps:
     - uses: actions/checkout@v2
diff --git a/.github/workflows/ci_basic.yml b/.github/workflows/ci_basic.yml
index d57d144b..80551785 100644
--- a/.github/workflows/ci_basic.yml
+++ b/.github/workflows/ci_basic.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7, 3.8]
+        python-version: ['3.9', '3.10']
 
     steps:
     - uses: actions/checkout@v2
diff --git a/.github/workflows/ci_basic_autograd.yml b/.github/workflows/ci_basic_autograd.yml
index bfcb7dfc..82c64f29 100644
--- a/.github/workflows/ci_basic_autograd.yml
+++ b/.github/workflows/ci_basic_autograd.yml
@@ -1,7 +1,7 @@
 # This workflow will install Python dependencies, run tests and lint with a variety of Python versions
 # For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions
 
-name: Tequila-Test-Basic-Autograd
+name: Tequila-Test-Basic-JAX
 
 on:
   push:
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7, 3.8, 3.9]
+        python-version: ['3.9', '3.10']
 
     steps:
     - uses: actions/checkout@v2
@@ -28,8 +28,8 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         pip install -r requirements.txt
-        pip uninstall -y jax jaxlib
-        pip install autograd
+        pip uninstall autograd
+        pip install -y jax jaxlib
     - name: Lint with flake8
       run: |
         pip install flake8
diff --git a/.github/workflows/ci_chemistry_madness.yml b/.github/workflows/ci_chemistry_madness.yml
index 2dd7e736..7c903227 100644
--- a/.github/workflows/ci_chemistry_madness.yml
+++ b/.github/workflows/ci_chemistry_madness.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7]
+        python-version: ['3.10']
 
     steps:
     - uses: actions/checkout@v2
diff --git a/.github/workflows/ci_chemistry_pyscf.yml b/.github/workflows/ci_chemistry_pyscf.yml
index 64d75b44..f68e6c64 100644
--- a/.github/workflows/ci_chemistry_pyscf.yml
+++ b/.github/workflows/ci_chemistry_pyscf.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7]
+        python-version: [3.9]
 
     steps:
     - uses: actions/checkout@v2
diff --git a/.github/workflows/ci_conda_madness.yml b/.github/workflows/ci_conda_madness.yml
index fc70be6d..ed042d4a 100644
--- a/.github/workflows/ci_conda_madness.yml
+++ b/.github/workflows/ci_conda_madness.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7]
+        python-version: [3.9]
 
     steps:
     - uses: actions/checkout@v2
diff --git a/.github/workflows/ci_ml.yml b/.github/workflows/ci_ml.yml
index aeb00d54..e649b377 100644
--- a/.github/workflows/ci_ml.yml
+++ b/.github/workflows/ci_ml.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7, 3.8]
+        python-version: [3.9]
 
     steps:
     - uses: actions/checkout@v2
diff --git a/.github/workflows/ci_optimizers.yml b/.github/workflows/ci_optimizers.yml
index 0d544f07..58f11e07 100644
--- a/.github/workflows/ci_optimizers.yml
+++ b/.github/workflows/ci_optimizers.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7, 3.8]
+        python-version: [3.9]
 
     steps:
     - uses: actions/checkout@v2
diff --git a/.github/workflows/ci_pyquil.yml b/.github/workflows/ci_pyquil.yml
index 92e39757..926ec8dc 100644
--- a/.github/workflows/ci_pyquil.yml
+++ b/.github/workflows/ci_pyquil.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7]
+        python-version: [3.9]
 
     steps:
     - uses: actions/checkout@v2
diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml
index 4f3f6f3e..09e44946 100644
--- a/.github/workflows/pypi.yml
+++ b/.github/workflows/pypi.yml
@@ -30,7 +30,7 @@ jobs:
         # remove qulacs from dependencies (issues with windows and mac)
         # users need to install themselves if they want it
        
-        cat requirements.txt | sed "s|qulacs|#qulacs|g" > tmp.txt
+        #cat requirements.txt | sed "s|qulacs|#qulacs|g" > tmp.txt
         rm requirements.txt
         mv tmp.txt requirements.txt
         python setup.py sdist bdist_wheel
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
deleted file mode 100644
index 1c80674b..00000000
--- a/.github/workflows/release.yml
+++ /dev/null
@@ -1,22 +0,0 @@
-name: 'Tweet when released'
-on:
-  release:
-    types: [released]
-
-jobs:
-  tweet:
-    runs-on: ubuntu-latest
-    steps:
-      - name: Tweet
-        id: tweet
-        uses: doomspec/auto-tweet-v2@v0.1.0
-        env:
-          CONSUMER_API_KEY: ${{ secrets.TWITTER_CONSUMER_API_KEY }}
-          CONSUMER_API_SECRET_KEY: ${{ secrets.TWITTER_CONSUMER_API_SECRET_KEY }}
-          ACCESS_TOKEN: ${{ secrets.TWITTER_ACCESS_TOKEN }}
-          ACCESS_TOKEN_SECRET: ${{ secrets.TWITTER_ACCESS_TOKEN_SECRET }}
-        with:
-          text: |
-            New version released: ${{ github.event.release.name }}
-            ${{ github.event.release.html_url }}
-      - run: echo ${{ steps.tweet.outputs.response }}

From 478fffb3db56e6fa7266b78cafa4aa6dfe3df776 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 3 Jun 2024 11:39:16 +0200
Subject: [PATCH 17/29] slight update in readme

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 37a0275e..f7640ac7 100644
--- a/README.md
+++ b/README.md
@@ -14,7 +14,7 @@ Tequila can execute the underlying quantum expectation values on state of the ar
 - [talks and slides](https://kottmanj.github.io/talks_and_material/)  
 
 # Installation
-Recommended Python version is 3.9 - 3.10.   
+Recommended Python version is 3.9 (3.10).   
 Tequila supports linux, osx and windows. However, not all optional dependencies are supported on windows.  
 
 ## Install from PyPi

From b8fd4df6da6fa59aa9fb5550c146299be31331da Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 3 Jun 2024 11:44:30 +0200
Subject: [PATCH 18/29] Update ci_optimizers.yml

---
 .github/workflows/ci_optimizers.yml | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/.github/workflows/ci_optimizers.yml b/.github/workflows/ci_optimizers.yml
index 58f11e07..a36b9d9e 100644
--- a/.github/workflows/ci_optimizers.yml
+++ b/.github/workflows/ci_optimizers.yml
@@ -23,12 +23,7 @@ jobs:
     - name: Set up Python ${{ matrix.python-version }}
       uses: actions/setup-python@v1
       with:
-        python-version: ${{ matrix.python-version }}
-    - name: Install quantum backends
-      run: |
-              pip install "cirq" "qiskit>=0.30" "qulacs"
-              pip uninstall pyquil -y
-              
+        python-version: ${{ matrix.python-version }}              
     - name: Install and test GPyOpt interface (no qibo)
       run: |  
               python -m pip install --upgrade pip 

From 6fab0bca0bfcba7063830162a5704b8b0007dce2 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 3 Jun 2024 11:46:03 +0200
Subject: [PATCH 19/29] Update ci_basic_autograd.yml

---
 .github/workflows/ci_basic_autograd.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/ci_basic_autograd.yml b/.github/workflows/ci_basic_autograd.yml
index 82c64f29..84a2d2e0 100644
--- a/.github/workflows/ci_basic_autograd.yml
+++ b/.github/workflows/ci_basic_autograd.yml
@@ -28,7 +28,7 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         pip install -r requirements.txt
-        pip uninstall autograd
+        pip uninstall autograd -y
         pip install -y jax jaxlib
     - name: Lint with flake8
       run: |

From 89b3c89d03a0b474268bc875cc558c29ae990abb Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 3 Jun 2024 11:49:14 +0200
Subject: [PATCH 20/29] Update ci_optimizers.yml

---
 .github/workflows/ci_optimizers.yml | 36 +----------------------------
 1 file changed, 1 insertion(+), 35 deletions(-)

diff --git a/.github/workflows/ci_optimizers.yml b/.github/workflows/ci_optimizers.yml
index a36b9d9e..9cc938c2 100644
--- a/.github/workflows/ci_optimizers.yml
+++ b/.github/workflows/ci_optimizers.yml
@@ -1,35 +1 @@
-# This workflow will install Python dependencies, run tests and lint with a variety of Python versions
-# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions
-
-name: Tequila-Test-Optimizers
-
-on:
-  push:
-    branches: [ master, devel ]
-  pull_request:
-    branches: [ master, devel ]
-
-jobs:
-
-  build:
-
-    runs-on: ubuntu-latest
-    strategy:
-      matrix:
-        python-version: [3.9]
-
-    steps:
-    - uses: actions/checkout@v2
-    - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v1
-      with:
-        python-version: ${{ matrix.python-version }}              
-    - name: Install and test GPyOpt interface (no qibo)
-      run: |  
-              python -m pip install --upgrade pip 
-              pip install --upgrade -r requirements.txt
-              pip install --upgrade -r requirements_gpyopt.txt
-              pip install -e .
-
-              pytest tests/test_gpyopt.py --slow
-
+# safe to delete

From 7e54e4ccfa47baf620984e43acb2dbd23b004d78 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 3 Jun 2024 11:50:01 +0200
Subject: [PATCH 21/29] Update requirements_gpyopt.txt

---
 requirements_gpyopt.txt | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/requirements_gpyopt.txt b/requirements_gpyopt.txt
index ac7fc7b0..e201bab6 100644
--- a/requirements_gpyopt.txt
+++ b/requirements_gpyopt.txt
@@ -1,3 +1,5 @@
+### GPyOpt no longer maintained ... still works with python 3.7
+### take care
 ### requirements for gpyopt optimizer
 cycler>=0.10.0
 decorator>=4.0.10

From 3350c8b7636b0aeea23dfcad40232e47888579ea Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 3 Jun 2024 11:51:04 +0200
Subject: [PATCH 22/29] Update ci_basic_autograd.yml

---
 .github/workflows/ci_basic_autograd.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/ci_basic_autograd.yml b/.github/workflows/ci_basic_autograd.yml
index 84a2d2e0..25799934 100644
--- a/.github/workflows/ci_basic_autograd.yml
+++ b/.github/workflows/ci_basic_autograd.yml
@@ -29,7 +29,7 @@ jobs:
         python -m pip install --upgrade pip
         pip install -r requirements.txt
         pip uninstall autograd -y
-        pip install -y jax jaxlib
+        pip install jax jaxlib
     - name: Lint with flake8
       run: |
         pip install flake8

From 887b828a66ccb3bac767d12e3326b9f25f797ce2 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Mon, 3 Jun 2024 12:11:19 +0200
Subject: [PATCH 23/29] dropping GPyOpt support -- no maintenance for 4 years

---
 .github/workflows/ci_optimizers.yml |  1 -
 README.md                           |  8 +-------
 requirements_gpyopt.txt             | 19 -------------------
 3 files changed, 1 insertion(+), 27 deletions(-)
 delete mode 100644 .github/workflows/ci_optimizers.yml
 delete mode 100644 requirements_gpyopt.txt

diff --git a/.github/workflows/ci_optimizers.yml b/.github/workflows/ci_optimizers.yml
deleted file mode 100644
index 9cc938c2..00000000
--- a/.github/workflows/ci_optimizers.yml
+++ /dev/null
@@ -1 +0,0 @@
-# safe to delete
diff --git a/README.md b/README.md
index f7640ac7..63c528f7 100644
--- a/README.md
+++ b/README.md
@@ -251,11 +251,7 @@ Quantum Equation of Motion with Orbital Optimization for Computing Molecular Pro
 Let us know, if you want your research project and/or tutorial to be included in this list!
 
 # Dependencies
-Support for additional optimizers or quantum backends can be activated by intalling them in your environment.
-Tequila will then detect them automatically.
-Currently those are: [Phoenics](https://github.com/aspuru-guzik-group/phoenics)
- and [GPyOpt](https://sheffieldml.github.io/GPyOpt/).
-Quantum backends are treated in the same way.
+Support for specific backends (quantum simulators, optimizers, quantum chemistry) can be activated by intalling them in your environment.
 
 ## Quantum Backends
 Currently supported
@@ -389,5 +385,3 @@ Tequila runs on Mac OSX.
 You might get in trouble with installing qulacs since it currently does not work with Apple's clang compiler.
 You need to install latest GNU compile (at least gcc-7 and g++7) and set them as default before installing qulacs over pip.
 
-## Qibo and GPyOpt
-Currently you can't use Qibo and GPyOpt within the same environment.
diff --git a/requirements_gpyopt.txt b/requirements_gpyopt.txt
deleted file mode 100644
index e201bab6..00000000
--- a/requirements_gpyopt.txt
+++ /dev/null
@@ -1,19 +0,0 @@
-### GPyOpt no longer maintained ... still works with python 3.7
-### take care
-### requirements for gpyopt optimizer
-cycler>=0.10.0
-decorator>=4.0.10
-numpy>=1.11.2
-six>=1.10.0
-python-dateutil>=2.6.0
-paramz>=0.7.0
-GPy>=1.8
-matplotlib>=1.5.3
-pyparsing>=2.1.10
-pytz>=2016.7
-scipy>=0.18.1
-mock>=2.0.0
-PyDOE >= 0.3.0
-sobol_seq >=0.1
-emcee>=2.2.1
-gpyopt

From 2d68f8470decf478f6f649efe0a66e8a791c029f Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Fri, 7 Jun 2024 00:34:35 +0200
Subject: [PATCH 24/29] control not passed form QubitExcitation to
 QubitExcitationImpl (#346)

---
 src/tequila/circuit/gates.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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))
 
 
 """

From 61eaee1bac40fcf7b114c102d4bb738925b2c9a9 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Tue, 18 Jun 2024 14:07:36 +0200
Subject: [PATCH 25/29] fixing error in the initialization of binary
 hamiltonians closes #394 (#351)

---
 src/tequila/grouping/binary_rep.py      |  2 +-
 src/tequila/quantumchemistry/qc_base.py | 12 +++++++++---
 2 files changed, 10 insertions(+), 4 deletions(-)

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/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index 560c6b91..5f174bb4 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -1487,7 +1487,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)
@@ -1722,7 +1723,7 @@ def rdm2(self):
 
     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,
-                     rdm_trafo: QubitHamiltonian = None):
+                     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.
@@ -2000,7 +2001,12 @@ def _build_2bdy_operators_hcb() -> list:
         
         # Compute expected values
         if rdm_trafo is None:
-            evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)
+            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)

From 48a00706895a193353a8e9209d22c7d78c606f5e Mon Sep 17 00:00:00 2001
From: JdelArco98 <francisco.del.arco.santos@uni-a.de>
Date: Thu, 20 Jun 2024 00:46:34 +0200
Subject: [PATCH 26/29] Fermionic Excitation Implementation  (#350)

* Fermionic Excitation Implementation
* made indices in Fermionic excitations more consistent
* Update chemistry_tools.py
---
 .../quantumchemistry/chemistry_tools.py       | 128 +++++++++++++++++-
 src/tequila/quantumchemistry/qc_base.py       |   8 +-
 tests/test_chemistry.py                       |   6 +-
 3 files changed, 130 insertions(+), 12 deletions(-)

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/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index 5f174bb4..f0f896e3 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -441,10 +441,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))
 
diff --git a/tests/test_chemistry.py b/tests/test_chemistry.py
index c4369442..41825cad 100644
--- a/tests/test_chemistry.py
+++ b/tests/test_chemistry.py
@@ -173,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()
@@ -370,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 = {}

From f5ee3648924f91451e6376b0748dabb333630af5 Mon Sep 17 00:00:00 2001
From: JdelArco98 <francisco.del.arco.santos@uni-a.de>
Date: Wed, 26 Jun 2024 11:13:41 +0200
Subject: [PATCH 27/29] Update qc_base.from tequila (#353)

---
 src/tequila/quantumchemistry/qc_base.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index f0f896e3..2064c04e 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -112,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

From 477af27d7bc4aca45207d185a9d8b740b1283b08 Mon Sep 17 00:00:00 2001
From: "J. S. Kottmann" <jakob.kottmann@gmail.com>
Date: Fri, 28 Jun 2024 14:31:48 +0200
Subject: [PATCH 28/29] Update simulator_qiskit.py (#354)

fixing u1, cu1 deprecation in qiskit
---
 src/tequila/simulators/simulator_qiskit.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/tequila/simulators/simulator_qiskit.py b/src/tequila/simulators/simulator_qiskit.py
index ded4e784..6438a4ca 100755
--- a/src/tequila/simulators/simulator_qiskit.py
+++ b/src/tequila/simulators/simulator_qiskit.py
@@ -137,7 +137,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 +169,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),
         }
 

From 981c0a7c72978331c3195521da4a37a664430e74 Mon Sep 17 00:00:00 2001
From: Benjamin Meyer <46791542+benni-tec@users.noreply.github.com>
Date: Fri, 28 Jun 2024 15:19:27 +0200
Subject: [PATCH 29/29] quantumchemistry: added encoding transformations (#345)

* quantum chemistry: added encoding transformations
* Update CI
---
 src/tequila/quantumchemistry/encodings.py     | 188 +++++++++++++-----
 .../quantumchemistry/madness_interface.py     |   5 +-
 src/tequila/quantumchemistry/qc_base.py       |   8 +-
 tests/test_chemistry.py                       |   8 +-
 tests/test_chemistry_madness.py               |   2 +-
 5 files changed, 160 insertions(+), 51 deletions(-)

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 2064c04e..7364a58e 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -1247,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:
diff --git a/tests/test_chemistry.py b/tests/test_chemistry.py
index 41825cad..d24c6bfa 100644
--- a/tests/test_chemistry.py
+++ b/tests/test_chemistry.py
@@ -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
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)