diff --git a/compile_with_cma.sh b/compile_with_cma.sh
new file mode 100755
index 0000000..1609da2
--- /dev/null
+++ b/compile_with_cma.sh
@@ -0,0 +1,46 @@
+#!/bin/bash
+
+# SLURM directives
+#SBATCH --job-name=VQE12qxDBQA
+#SBATCH --output=boost_%A_%a.log
+#SBATCH --array=0-29  # Adjusted dynamically based on the number of tasks.
+
+OPTIMIZATION_METHOD="cma"
+OPTIMIZATION_CONFIG="{ \"maxiter\": 50}"
+
+# Define the epoch values
+EPOCHS=(100 200 300)
+
+# Dynamically calculate the folders
+FOLDERS=($(ls -d ./results/scaling_data/test_nqubits_12/*/))
+NUM_FOLDERS=${#FOLDERS[@]}
+NUM_EPOCHS=${#EPOCHS[@]}
+
+if [ $NUM_FOLDERS -eq 0 ]; then
+    echo "No folders found in ./results/scaling_data/test_nqubits_12. Exiting."
+    exit 1
+fi
+
+# Total number of tasks
+TOTAL_TASKS=$((NUM_FOLDERS * NUM_EPOCHS))
+
+# Adjust SLURM array size dynamically
+if [ "$SLURM_ARRAY_TASK_ID" -ge "$TOTAL_TASKS" ]; then
+    echo "Task ID exceeds available tasks. Exiting."
+    exit 1
+fi
+
+# Determine the folder and epoch for this task
+FOLDER_INDEX=$((SLURM_ARRAY_TASK_ID / NUM_EPOCHS))
+EPOCH_INDEX=$((SLURM_ARRAY_TASK_ID % NUM_EPOCHS))
+
+FOLDER=${FOLDERS[$FOLDER_INDEX]}
+EPOCH=${EPOCHS[$EPOCH_INDEX]}
+
+echo "Processing folder: $FOLDER with epoch: $EPOCH"
+
+# Run the Python script for this folder and epoch
+python3 load_vqe_and_rotate.py      --backend numpy \
+                                    --path "$FOLDER" \
+                                    --epoch $EPOCH --steps 3 --optimization_method $OPTIMIZATION_METHOD \
+                                    --optimization_config "$OPTIMIZATION_CONFIG"
diff --git a/load_vqe_and_rotate.py b/load_vqe_and_rotate.py
new file mode 100644
index 0000000..793eb71
--- /dev/null
+++ b/load_vqe_and_rotate.py
@@ -0,0 +1,271 @@
+import argparse
+import json
+import logging
+import pathlib
+import time
+from copy import deepcopy
+
+import numpy as np
+import qibo
+
+from qibo import hamiltonians
+from qibo.backends import construct_backend
+from qibo.quantum_info.metrics import fidelity
+
+from boostvqe import ansatze
+from boostvqe.models.dbi import double_bracket_evolution_oracles
+from boostvqe.models.dbi.double_bracket_evolution_oracles import (
+    FrameShiftedEvolutionOracle,
+    IsingNNEvolutionOracle,
+    MagneticFieldEvolutionOracle,
+    XXZ_EvolutionOracle,
+)
+from boostvqe.models.dbi.group_commutator_iteration_transpiler import (
+    DoubleBracketRotationType,
+    GroupCommutatorIterationWithEvolutionOracles,
+)
+from boostvqe.utils import (
+    OPTIMIZATION_FILE,
+    PARAMS_FILE,
+    optimize_D,
+    select_recursion_step_gd_circuit,
+)
+
+logging.basicConfig(level=logging.INFO)
+qibo.set_backend("numpy")
+
+
+def dump_config(config: dict, path):
+    config["path"] = config["path"]
+    config["db_rotation"] = config["db_rotation"].name
+    (path / "config.json").write_text(json.dumps(config, indent=4))
+
+
+def main(args):
+    """Loading trained VQE, hyperoptimize DBI and apply the DBQA."""
+    # load the VQE results path
+    path = pathlib.Path(args.path)
+    dump_path = (
+        path
+        / f"{args.db_rotation.name}_{args.optimization_method}_{args.epoch}e_{args.steps}s"
+    )
+    dump_path.mkdir(parents=True, exist_ok=True)
+
+    config = json.loads((path / OPTIMIZATION_FILE).read_text())
+    dump_config(deepcopy(vars(args)), path=dump_path)
+
+    if args.optimization_config is None:
+        opt_options = {}
+    else:
+        opt_options = json.loads(args.optimization_config)
+
+    try:
+        params = np.load(path / f"parameters/params_ite{args.epoch}.npy")
+    except FileNotFoundError:
+        params = np.array(
+            np.load(path / PARAMS_FILE, allow_pickle=True).tolist()[0][args.epoch]
+        )
+
+    nqubits = config["nqubits"]
+    nlayers = config["nlayers"]
+    vqe_backend = construct_backend(backend=args.backend, platform=args.platform)
+    # TODO: remove delta hardcoded
+    print
+    hamiltonian = getattr(hamiltonians, config["hamiltonian"])(
+        nqubits=nqubits, delta=0.5, backend=vqe_backend
+    )
+
+    # construct circuit from parsed ansatz name
+    circ = getattr(ansatze, args.circuit_ansatz)(config["nqubits"], config["nlayers"])
+
+    vqe = ansatze.VQE(
+        circuit=circ,
+        hamiltonian=hamiltonian,
+    )
+    vqe.circuit.set_parameters(params)
+
+    base_oracle = XXZ_EvolutionOracle.from_nqubits(
+        nqubits=nqubits, delta=0.5, steps=args.steps, order=args.order
+    )
+    oracle = FrameShiftedEvolutionOracle.from_evolution_oracle(
+        before_circuit=vqe.circuit.invert(),
+        after_circuit=vqe.circuit,
+        base_evolution_oracle=base_oracle,
+    )
+
+    gci = GroupCommutatorIterationWithEvolutionOracles(
+        oracle,
+        args.db_rotation,
+    )
+
+    # TODO: remove hardcoded magnetic field
+    eo_d_type = getattr(double_bracket_evolution_oracles, args.eo_d)
+
+    print(
+        f"The gci mode is {gci.double_bracket_rotation_type} rotation with {eo_d_type.__name__} as the oracle.\n"
+    )
+    metadata = {}
+
+    for gci_step_nmb in range(args.steps):
+        logging.info(
+            "\n################################################################################\n"
+            + f"Optimizing GCI step {gci_step_nmb+1} with optimizer {args.optimization_method}"
+            + "\n################################################################################\n"
+        )
+        it = time.time()
+        if args.optimization_method == "sgd":
+            params = (
+                [4 - np.sin(x / 3) for x in range(nqubits)]
+                if eo_d_type == MagneticFieldEvolutionOracle
+                else [4 - np.sin(x / 3) for x in range(nqubits)] + nqubits * [1]
+            )
+            mode, best_s, best_b, eo_d = select_recursion_step_gd_circuit(
+                gci,
+                mode=args.db_rotation,
+                eo_d_type=eo_d_type,
+                params=params,
+                step_grid=np.linspace(1e-5, 2e-2, 30),
+                lr_range=(1e-3, 1),
+                nmb_gd_epochs=opt_options["gd_epochs"],
+                threshold=1e-4,
+                max_eval_gd=30,
+            )
+
+            opt_dict = {"sgd_extras": "To be defined"}
+
+        else:
+            if gci_step_nmb == 0:
+                p0 = [0.01]
+                if eo_d_type == MagneticFieldEvolutionOracle:
+                    p0.extend([4 - np.sin(x / 3) for x in range(nqubits)])
+                elif eo_d_type == IsingNNEvolutionOracle:
+                    p0.extend(
+                        [4 - np.sin(x / 3) for x in range(nqubits)] + nqubits * [1]
+                    )
+
+            else:
+                p0 = [best_s]
+                p0.extend(best_b)
+            optimized_params, opt_dict = optimize_D(
+                params=p0,
+                gci=gci,
+                eo_d_type=eo_d_type,
+                mode=args.db_rotation,
+                method=args.optimization_method,
+                **opt_options,
+            )
+            best_s = optimized_params[0]
+            best_b = optimized_params[1:]
+            eo_d = eo_d_type.load(best_b)
+
+        step_data = dict(
+            best_s=best_s,
+            eo_d_name=eo_d.__class__.__name__,
+            eo_d_params=eo_d.params,
+        )
+        logging.info(f"Total optimization time required: {time.time() - it} seconds")
+        gci.mode_double_bracket_rotation = args.db_rotation
+
+        gci(best_s, eo_d, args.db_rotation)
+
+        this_report = report(vqe, hamiltonian, gci, best_s, eo_d, args.db_rotation)
+        print_report(this_report)
+        metadata[gci_step_nmb + 1] = this_report | step_data | opt_dict
+
+    (dump_path / "boosting_data.json").write_text(json.dumps(metadata, indent=4))
+
+
+def report(vqe, hamiltonian, gci, step, eo_d, mode):
+    energies = hamiltonian.eigenvalues()
+    ground_state_energy = float(energies[0])
+    vqe_energy = float(hamiltonian.expectation(vqe.circuit().state()))
+    gci_loss = float(gci.loss(step, eo_d, mode))
+    gap = float(energies[1] - energies[0])
+
+    return (
+        dict(
+            nqubits=hamiltonian.nqubits,
+            gci_loss=float(gci_loss),
+            vqe_energy=float(vqe_energy),
+            target_energy=ground_state_energy,
+            diff_vqe_target=vqe_energy - ground_state_energy,
+            diff_gci_target=gci_loss - ground_state_energy,
+            gap=gap,
+            diff_vqe_target_perc=abs(vqe_energy - ground_state_energy)
+            / abs(ground_state_energy)
+            * 100,
+            diff_gci_target_perc=abs(gci_loss - ground_state_energy)
+            / abs(ground_state_energy)
+            * 100,
+            fidelity_witness_vqe=1 - (vqe_energy - ground_state_energy) / gap,
+            fidelity_witness_gci=1 - (gci_loss - ground_state_energy) / gap,
+            fidelity_vqe=fidelity(vqe.circuit().state(), hamiltonian.ground_state()),
+            fidelity_gci=fidelity(
+                gci.get_composed_circuit()().state(), hamiltonian.ground_state()
+            ),
+        )
+        | gci.get_gate_count_dict()
+    )
+
+
+def print_report(report: dict):
+    print(
+        f"\
+    The target energy is {report['target_energy']}\n\
+    The VQE energy is {report['vqe_energy']} \n\
+    The DBQA energy is {report['gci_loss']}. \n\
+    The difference is for VQE is {report['diff_vqe_target']} \n\
+    and for the DBQA {report['diff_gci_target']} \n\
+    which can be compared to the spectral gap {report['gap']}.\n\
+    The relative difference is \n\
+        - for VQE {report['diff_vqe_target_perc']}% \n\
+        - for DBQA {report['diff_gci_target_perc']}%.\n\
+    The energetic fidelity witness of the ground state is: \n\
+        - for the VQE  {report['fidelity_witness_vqe']} \n\
+        - for DBQA {report['fidelity_witness_gci']}\n\
+    The true fidelity is \n\
+        - for the VQE  {report['fidelity_vqe']}\n\
+        - for DBQA {report['fidelity_gci']}\n\
+                    "
+    )
+    print(
+        f"The boosting circuit used {report['nmb_cnot']} CNOT gates coming from compiled XXZ evolution and {report['nmb_cz']} CZ gates from VQE.\n\
+For {report['nqubits']} qubits this gives n_CNOT/n_qubits = {report['nmb_cnot_relative']} and n_CZ/n_qubits = {report['nmb_cz_relative']}"
+    )
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Boosting VQE with DBI.")
+    parser.add_argument("--backend", default="numpy", type=str, help="Qibo backend")
+    parser.add_argument("--platform", default=None, type=str, help="Qibo backend platform")
+    parser.add_argument("--path", type=str, help="Output folder")
+    parser.add_argument(
+        "--epoch", default=-1, type=int, help="VQE epoch where DBI will be applied."
+    )
+    parser.add_argument(
+        "--eo_d",
+        default="IsingNNEvolutionOracle",
+        help="Evolution Oracle D operator. Can be either MagneticFieldEvolutionOracle or IsingNNEvolutionOracle.",
+    )
+    parser.add_argument("--steps", default=3, type=int, help="DBI steps")
+    parser.add_argument("--order", default=2, type=int, help="Suzuki-Trotter order")
+    parser.add_argument(
+        "--db_rotation",
+        type=lambda arg: DoubleBracketRotationType[arg],
+        choices=DoubleBracketRotationType,
+        default="group_commutator_third_order_reduced",
+        help="DB rotation type.",
+    )
+    parser.add_argument(
+        "--optimization_method", default="sgd", type=str, help="Optimization method"
+    )
+    parser.add_argument(
+        "--optimization_config",
+        type=str,
+        help="Options to customize the optimizer training.",
+    )
+    parser.add_argument(
+        "--circuit_ansatz", default="hw_preserving", type=str, help="Circuit ansatz"
+    )
+    args = parser.parse_args()
+    main(args)
diff --git a/results/scaling_4_12.zip b/results/scaling_4_12.zip
new file mode 100644
index 0000000..543a422
Binary files /dev/null and b/results/scaling_4_12.zip differ
diff --git a/run_vqes.sh b/run_vqes.sh
new file mode 100644
index 0000000..0e47106
--- /dev/null
+++ b/run_vqes.sh
@@ -0,0 +1,34 @@
+#!/bin/bash
+#SBATCH --output=output_%A_%a.out      # Output file for each array task
+#SBATCH --error=error_%A_%a.err        # Error file for each array task
+#SBATCH --array=0-59                   # Total tasks: 6 nqubits * 10 seeds = 60
+
+# Define variables for the parameters
+CIRCUIT_ANSAZ="hw_preserving"
+NLAYERS=5
+BACKEND="qiboml"
+PLATFORM="tensorflow"
+OPTIMIZER="sgd"
+OPTIMIZER_OPTIONS="{ \"optimizer\": \"Adam\", \"learning_rate\": 0.05, \"nmessage\": 1, \"nepochs\": 400 }"
+
+# Define the nqubits and seed values
+NQUBITS_VALUES=(2 4 6 8 10 12)
+SEED_VALUES=(0 1 2 3 4 5 6 7 8 9)
+
+# Calculate indices for nqubits and seed based on SLURM_ARRAY_TASK_ID
+NQUBITS_INDEX=$((SLURM_ARRAY_TASK_ID / 10))
+SEED_INDEX=$((SLURM_ARRAY_TASK_ID % 10))
+
+# Set nqubits and seed based on the calculated indices
+NQUBITS=${NQUBITS_VALUES[NQUBITS_INDEX]}
+SEED=${SEED_VALUES[SEED_INDEX]}
+
+# Set a specific job name for each task with the current nqubits
+#SBATCH --job-name=vqe_${NQUBITS}
+
+# Run the Python script with the specified arguments, including nqubits and seed
+python train_vqes.py    --output_folder "results/test_nqubits_${NQUBITS}" --circuit_ansatz "$CIRCUIT_ANSAZ" \
+                        --nqubits "$NQUBITS" --nlayers "$NLAYERS" --backend $BACKEND \
+                        --platform $PLATFORM --optimizer $OPTIMIZER \
+                        --optimizer_options "$OPTIMIZER_OPTIONS" \
+                        --seed $SEED
diff --git a/src/boostvqe/ansatze.py b/src/boostvqe/ansatze.py
index 2867ad0..69fe558 100644
--- a/src/boostvqe/ansatze.py
+++ b/src/boostvqe/ansatze.py
@@ -2,6 +2,7 @@
 from qibo import gates, get_backend
 from qibo.backends import construct_backend
 from qibo.models import Circuit
+from qibo.config import raise_error
 
 from boostvqe.training_utils import vqe_loss