diff --git a/doc/source/tutorials/ansatz.rst b/doc/source/tutorials/ansatz.rst
index 9d8cc3d..3dc4d4f 100644
--- a/doc/source/tutorials/ansatz.rst
+++ b/doc/source/tutorials/ansatz.rst
@@ -88,6 +88,10 @@ The following example demonstrates how the energy of the H2 molecule is affected
     |   0.2 |   0.673325849299  |
     -----------------------------
 
+..
+    Placeholder for basis rotation ansatz documentation
+
+
 .. _UCC Ansatz:
 
 Unitary Coupled Cluster Ansatz
diff --git a/examples/br_example.py b/examples/br_example.py
index 36ba772..c0813b5 100644
--- a/examples/br_example.py
+++ b/examples/br_example.py
@@ -4,10 +4,11 @@
 """
 
 import numpy as np
+from qibo import Circuit, gates, models
 from qibo.optimizers import optimize
 from scipy.optimize import minimize
 
-from qibochem.ansatz.basis_rotation import br_circuit
+from qibochem.ansatz import basis_rotation
 from qibochem.ansatz.hf_reference import hf_circuit
 from qibochem.driver.molecule import Molecule
 from qibochem.measurement.expectation import expectation
@@ -19,7 +20,6 @@
 except ModuleNotFoundError:
     mol.run_psi4()
 
-
 # Diagonalize H_core to use as the guess wave function
 
 # First, symmetric orthogonalization of overlap matrix S using np:
@@ -47,35 +47,32 @@
 )
 
 
-def electronic_energy(parameters):
-    """
-    Loss function (Electronic energy) for the basis rotation ansatz
-    """
-    circuit = hf_circuit(mol.nso, mol.nelec)
-    circuit += br_circuit(mol.nso, parameters, mol.nelec)
-
-    return expectation(circuit, hamiltonian)
+def basis_rotation_circuit(mol, parameters=0.0):
 
+    nqubits = mol.nso
+    occ = range(0, mol.nelec)
+    vir = range(mol.nelec, mol.nso)
 
-# Build  a basis rotation_circuit
-params = np.random.rand(mol.nelec * (mol.nso - mol.nelec))  # n_occ * n_virt
+    U, theta = basis_rotation.unitary(occ, vir, parameters=parameters)
 
-best, params, extra = optimize(electronic_energy, params)
+    gate_angles, final_U = basis_rotation.givens_qr_decompose(U)
+    gate_layout = basis_rotation.basis_rotation_layout(nqubits)
+    gate_list, ordered_angles = basis_rotation.basis_rotation_gates(gate_layout, gate_angles, theta)
 
-print("\nResults using Qibo optimize:")
-print(f" HF energy: {mol.e_hf:.8f}")
-print(f"VQE energy: {best:.8f} (Basis rotation ansatz)")
-# print()
-# print("Optimized parameters:", params)
+    circuit = Circuit(nqubits)
+    for _i in range(mol.nelec):
+        circuit.add(gates.X(_i))
+    circuit.add(gate_list)
 
+    return circuit, gate_angles
 
-params = np.random.rand(mol.nelec * (mol.nso - mol.nelec))  # n_occ * n_virt
 
-result = minimize(electronic_energy, params)
-best, params = result.fun, result.x
+br_circuit, qubit_parameters = basis_rotation_circuit(mol, parameters=0.1)
+vqe = models.VQE(br_circuit, hamiltonian)
+vqe_result = vqe.minimize(qubit_parameters)
 
-print("\nResults using scipy.optimize:")
 print(f" HF energy: {mol.e_hf:.8f}")
-print(f"VQE energy: {best:.8f} (Basis rotation ansatz)")
-# print()
-# print("Optimized parameters:", params)
+print(f"VQE energy: {vqe_result[0]:.8f} (Basis rotation ansatz)")
+print()
+print("Optimized qubit parameters:\n", vqe_result[1])
+print("Optimizer message:\n", vqe_result[2])
diff --git a/src/qibochem/ansatz/__init__.py b/src/qibochem/ansatz/__init__.py
index 8a70e56..a94b14c 100644
--- a/src/qibochem/ansatz/__init__.py
+++ b/src/qibochem/ansatz/__init__.py
@@ -1,9 +1,8 @@
 from qibochem.ansatz.basis_rotation import (
-    br_circuit,
-    givens_rotation_gate,
-    givens_rotation_parameters,
-    swap_matrices,
-    unitary_rot_matrix,
+    basis_rotation_gates,
+    basis_rotation_layout,
+    givens_qr_decompose,
+    unitary,
 )
 from qibochem.ansatz.hardware_efficient import hea
 from qibochem.ansatz.hf_reference import bk_matrix, bk_matrix_power2, hf_circuit
diff --git a/src/qibochem/ansatz/basis_rotation.py b/src/qibochem/ansatz/basis_rotation.py
index c5a540a..e79205d 100644
--- a/src/qibochem/ansatz/basis_rotation.py
+++ b/src/qibochem/ansatz/basis_rotation.py
@@ -3,145 +3,280 @@
 """
 
 import numpy as np
-from openfermion.linalg.givens_rotations import givens_decomposition
+
+# from openfermion.linalg.givens_rotations import givens_decomposition
 from qibo import gates, models
 from scipy.linalg import expm
 
+from qibochem.ansatz import ucc
+from qibochem.driver import hamiltonian
+
 # Helper functions
 
 
-def unitary_rot_matrix(parameters, occ_orbitals, virt_orbitals):
+def unitary(occ_orbitals, virt_orbitals, parameters=None):
     r"""
-    Returns the unitary rotation matrix U = exp(\kappa) mixing the occupied and virtual orbitals,
-         where \kappa is an anti-Hermitian matrix
-
+    Returns the unitary rotation matrix :math: `U = \exp(\kappa)` mixing the occupied and virtual orbitals.
+    Orbitals are arranged in alternating spins, e.g. for 4 occupied orbitals [0,1,2,3], the spins are arranged as [0a, 0b, 1a, 1b]
+    Dimension for the array of rotation parameters is len(occ_orbitals_a)*len(virt_orbitals_a) + len(occ_orbitals_b)*len(virt_orbitals_b)
+    The current implementation of this function only accommodates systems with all electrons paired, number of alpha and beta spin electrons are equal
     Args:
-        parameters: List/array of rotation parameters. dim = len(occ_orbitals)*len(virt_orbitals)
         occ_orbitals: Iterable of occupied orbitals
         virt_orbitals: Iterable of virtual orbitals
+        parameters: List/array of rotation parameters
+            dimension = len(occ_orbitals)*len(virt_orbitals)
+    Returns:
+        exp(k): Unitary matrix of Givens rotations, obtained by matrix exponential of skew-symmetric
+                kappa matrix
     """
-    n_orbitals = len(occ_orbitals) + len(virt_orbitals)
 
+    # conserve_spin has to be true for SCF/basis_rotation cases, else expm(k) is not unitary
+    ov_pairs = ucc.generate_excitations(1, occ_orbitals, virt_orbitals, conserve_spin=True)
+    # print('ov_pairs presort', ov_pairs)
+    ov_pairs = ucc.sort_excitations(ov_pairs)
+    # print('ov_pairs sorted ', ov_pairs)
+    n_theta = len(ov_pairs)
+    if parameters is None:
+        print("basis rotation: parameters not specified")
+        print("basis rotation: using default occ-virt rotation parameter value = 0.0")
+        parameters = np.zeros(n_theta)
+    elif isinstance(parameters, float):
+        print("basis rotation: using uniform value of", parameters, "for each parameter value")
+        parameters = np.full(n_theta, parameters)
+    else:
+        if len(parameters) != n_theta:
+            raise IndexError("parameter array specified has bad size or type")
+        print("basis rotation: loading parameters from input")
+
+    n_orbitals = len(occ_orbitals) + len(virt_orbitals)
     kappa = np.zeros((n_orbitals, n_orbitals))
-    # Get all possible (occ, virt) pairs
-    orbital_pairs = [(_o, _v) for _o in occ_orbitals for _v in virt_orbitals]
-    # occ/occ and virt/virt orbital mixing doesn't change the energy, so can ignore
-    for _i, (_occ, _virt) in enumerate(orbital_pairs):
+
+    for _i, (_occ, _virt) in enumerate(ov_pairs):
         kappa[_occ, _virt] = parameters[_i]
         kappa[_virt, _occ] = -parameters[_i]
-    return expm(kappa)
 
+    return expm(kappa), parameters
 
-def swap_matrices(permutations, n_qubits):
-    r"""
-    Build matrices that permute (swap) the columns of a given matrix
 
+# clements qr routines
+def givens_qr_decompose(U):
+    r"""
+    Clements scheme QR decompose a unitary matrix U using Givens rotations
+    see arxiv:1603.08788
     Args:
-        permutations [List(tuple, ), ]: List of lists of 2-tuples permuting two consecutive numbers
-            e.g. [[(0, 1), (2, 3)], [(1, 2)], ...]
-        n_qubits: Dimension of matrix (\exp(\kappa) matrix) to be operated on i.e. # of qubits
-
+        U: unitary rotation matrix
     Returns:
-        List of matrices that carry out \exp(swap) for each pair in permutations
+        z_angles: array of rotation angles
+        U: final unitary after QR decomposition, should be an identity
     """
-    exp_swaps = []
-    _temp = np.array([[-1.0, 1.0], [1.0, -1.0]])
-    for pairs in permutations:
-        gen = np.zeros((n_qubits, n_qubits))
-        for pair in pairs:
-            # Add zeros to pad out the (2, 2) matrix to (n_qubits, n_qubits) with 0.
-            gen += np.pad(_temp, ((pair[0], n_qubits - pair[1] - 1),), "constant")
-        # Take matrix exponential, remove all entries close to zero, and convert to real matrix
-        exp_mat = expm(-0.5j * np.pi * gen)
-        exp_mat[np.abs(exp_mat) < 1e-10] = 0.0
-        exp_mat = np.real_if_close(exp_mat)
-        # Add exp_mat to exp_swaps once cleaned up
-        exp_swaps.append(exp_mat)
-    return exp_swaps
-
-
-def givens_rotation_parameters(n_qubits, exp_k, n_occ):
-    r"""
-    Get the parameters of the Givens rotation matrix obtained after decomposing \exp(\kappa)
 
-    Args:
-        n_qubits: Number of qubits (i.e. spin-orbitals)
-        exp_k: Unitary rotation matrix
-        n_occ: Number of occupied orbitals
-    """
-    # List of 2-tuples to link all qubits via swapping
-    swap_pairs = [
-        [(_i, _i + 1) for _i in range(_d % 2, n_qubits - 1, 2)]
-        # Only 2 possible ways of swapping: [(0, 1), ... ] or [(1, 2), ...]
-        for _d in range(2)
-    ]
-    # Get their corresponding matrices
-    swap_unitaries = swap_matrices(swap_pairs, n_qubits)
-
-    # Create a copy of exp_k for rotating
-    exp_k_shift = np.copy(exp_k)
-    # Apply the column-swapping transformations
-    for u_rot in swap_unitaries:
-        exp_k_shift = u_rot @ exp_k_shift
-    # Only want the n_occ by n_qubits (n_orb) sub-matrix
-    qmatrix = exp_k_shift.T[:n_occ, :]
-
-    # Apply Givens decomposition of the submatrix to get a list of individual Givens rotations
-    # (orb1, orb2, theta, phi), where orb1 and orb2 are rotated by theta(real) and phi(imag)
-    qdecomp, _, _ = givens_decomposition(qmatrix)
-    # Lastly, reverse the list of Givens rotations (because U = G_k ... G_2 G_1)
-    return list(reversed(qdecomp))
-
-
-def givens_rotation_gate(n_qubits, orb1, orb2, theta):
-    """
-    Circuit corresponding to a Givens rotation between two qubits (spin-orbitals)
+    def move_step(U, row, col, row_change, col_change):
+        """
+        internal function to move a step in Givens QR decomposition of
+        unitary rotation matrix
+        """
+        row += row_change
+        col += col_change
+        return U, row, col
 
-    Args:
-        n_qubits: Number of qubits in the quantum circuit
-        orb1, orb2: Orbitals used for the Givens rotation
-        theta: Rotation parameter
+    def row_op(U, row, col):
+        """
+        internal function to zero out a row using Givens rotation
+        with angle z
+        """
+        Uc = U.copy()
+        srow = row - 1
+        scol = col
+        z = np.arctan2(-Uc[row][col], Uc[srow][scol])
+        temp_srow = Uc[srow, :]
+        temp_row = Uc[row, :]
+        new_srow = np.cos(z) * temp_srow - np.sin(z) * temp_row
+        new_row = np.sin(z) * temp_srow + np.cos(z) * temp_row
+        Uc[srow, :] = new_srow
+        Uc[row, :] = new_row
+        return z, Uc
+
+    def col_op(U, row, col):
+        """
+        internal function to zero out a column using Givens rotation
+        with angle z
+        """
+        Uc = U.copy()
+        srow = row
+        scol = col + 1
+        z = np.arctan2(-Uc[row][col], Uc[srow][scol])
+        temp_scol = Uc[:, scol]
+        temp_col = Uc[:, col]
+        new_scol = np.cos(z) * temp_scol - np.sin(z) * temp_col
+        new_col = np.sin(z) * temp_scol + np.cos(z) * temp_col
+        Uc[:, scol] = new_scol
+        Uc[:, col] = new_col
+        return z, Uc
+
+    N = len(U[0])
+
+    z_array = []
+
+    # start QR from bottom left element
+    row = N - 1
+    col = 0
+    z, U = col_op(U, row, col)
+    z_array.append(z)
+    # traverse the U matrix in diagonal-zig-zag manner until the main
+    # diagonal is reached
+    # if move = up, do a row op
+    # if move = diagonal-down-right, do a row op
+    # if move = right, do a column op
+    # if move = diagonal-up-left, do a column op
+    while row != 1:
+        U, row, col = move_step(U, row, col, -1, 0)
+        z, U = row_op(U, row, col)
+        z_array.append(z)
+        while row < N - 1:
+            U, row, col = move_step(U, row, col, 1, 1)
+            z, U = row_op(U, row, col)
+            z_array.append(z)
+        if col != N - 2:
+            U, row, col = move_step(U, row, col, 0, 1)
+            z, U = col_op(U, row, col)
+            z_array.append(z)
+        # else:
+        #    break
+        while col > 0:
+            U, row, col = move_step(U, row, col, -1, -1)
+            z, U = col_op(U, row, col)
+            z_array.append(z)
+
+    return z_array, U
 
+
+def basis_rotation_layout(N):
+    r"""
+    generates the layout of the basis rotation circuit for Clements scheme QR decomposition
+    Args:
+        N: number of qubits/modes
     Returns:
-        circuit: Qibo Circuit object representing a Givens rotation between orb1 and orb2
+        A:  NxN matrix, with -1 being null, 0 is the control and integers
+            1 or greater being the index for angles in clements QR decomposition of
+            the unitary matrix representing the unitary transforms that
+            rotate the basis
     """
-    # Define 2x2 sqrt(iSwap) matrix
-    iswap_mat = np.array([[1.0, 1.0j], [1.0j, 1.0]]) / np.sqrt(2.0)
-    # Build and add the gates to circuit
-    circuit = models.Circuit(n_qubits)
-    circuit.add(gates.GeneralizedfSim(orb1, orb2, iswap_mat, 0.0, trainable=False))
-    circuit.add(gates.RZ(orb1, -theta))
-    circuit.add(gates.RZ(orb2, theta + np.pi))
-    circuit.add(gates.GeneralizedfSim(orb1, orb2, iswap_mat, 0.0, trainable=False))
-    circuit.add(gates.RZ(orb1, np.pi, trainable=False))
-    return circuit
-
-
-def br_circuit(n_qubits, parameters, n_occ):
-    """
-    Google's basis rotation circuit, applied between the occupied/virtual orbitals. Forms the exp(kappa) matrix, decomposes
-    it into Givens rotations, and sets the circuit parameters based on the Givens rotation decomposition. Note: Supposed
-    to be used with the JW fermion-to-qubit mapping
 
-    Args:
-        n_qubits: Number of qubits in the quantum circuit
-        parameters: Rotation parameters for exp(kappa); Must have (n_occ * n_virt) parameters
-        n_occ: Number of occupied orbitals
+    def move_step(A, row, col, row_change, col_change):
+        """
+        internal function to move a step in gate layout of br circuit
+        """
+        row += row_change
+        col += col_change
+        return A, row, col, updown
 
-    Returns:
-        Qibo ``Circuit``: Circuit corresponding to the basis rotation ansatz between the occupied and virtual orbitals
+    def check_top(A, row, col, row_change, col_change):
+        """
+        check if we are at the top of layout matrix,
+        return false if we are at the top, i.e. row == 1
+        (row 0 is not assigned any operation; just control)
+        """
+        move_step(A, row, col, row_change, col_change)
+        return row > 1
+
+    def check_bottom(A, row, col, row_change, col_change):
+        """
+        check if we are at the bottom of layout matrix
+        return false if we are at the bottom, i.e. row == N-1
+        """
+        N = len(A)
+        move_step(A, row, col, row_change, col_change)
+        return row < N - 1
+
+    def jump_right(A, row, col, updown):
+        """ """
+        N = len(A)
+        row = N - 2 - col
+        col = N - 1
+        updown = -1
+        return A, row, col, updown
+
+    def jump_left(A, row, col, updown):
+        N = len(A)
+        row = N + 1 - col
+        col = 0
+        updown = 1
+        return A, row, col, updown
+
+    def assign_element(A, row, col, k):
+        A[row - 1][col] = 0
+        A[row][col] = k
+        return A, row, col
+
+    nangles = (N - 1) * N // 2  # half-triangle
+    A = np.full([N, N], -1, dtype=int)
+
+    # first step
+    row = 1
+    col = 0
+    k = 1
+    A, row, col = assign_element(A, row, col, k)
+    k += 1
+    updown = 1
+
+    while k <= nangles:
+
+        if updown == 1:
+            if check_top(A, row, col, -1, 1):
+                A, row, col, updown = move_step(A, row, col, -1, 1)
+                A, row, col = assign_element(A, row, col, k)
+            else:  # jump
+                # print('jump_right')
+                A, row, col, updown = jump_right(A, row, col, updown)
+                A, row, col = assign_element(A, row, col, k)
+        elif updown == -1:
+            if check_bottom(A, row, col, 1, -1):
+                A, row, col, updown = move_step(A, row, col, 1, -1)
+                A, row, col = assign_element(A, row, col, k)
+            else:  # jump
+                # print('jump left')
+                A, row, col, updown = jump_left(A, row, col, updown)
+                A, row, col = assign_element(A, row, col, k)
+        # else:
+        #    raise ValueError("Bad direction")
+
+        # print(row, col)
+        k += 1
+
+    return A
+
+
+def basis_rotation_gates(A, z_array, parameters):
+    r"""
+    places the basis rotation gates on circuit in the order of Clements scheme QR decomposition
+    Args:
+        A:
+            NxN matrix, with -1 being null, 0 is the control and integers
+            1 or greater being the index for angles in clements QR decomposition of
+            the unitary matrix representing the unitary transforms that
+            rotate the basis
+        z_array:
+            array of givens rotation angles in order of traversal from
+            QR decomposition
+        parameters:
+            array of parameters in order of traversal from QR decomposition
+    Outputs:
+        gate_list:
+            list of gates which implement the basis rotation using Clements scheme QR decomposition
+        ordered_angles:
+            list of angles ordered by sequence of singles excitation gates added to circuit
     """
-    assert len(parameters) == (n_occ * (n_qubits - n_occ)), "Need len(parameters) == (n_occ * n_virt)"
-    # Unitary rotation matrix \exp(\kappa)
-    exp_k = unitary_rot_matrix(parameters, range(n_occ), range(n_occ, n_qubits))
-    # List of Givens rotation parameters
-    g_rotation_parameters = givens_rotation_parameters(n_qubits, exp_k, n_occ)
-
-    # Start building the circuit
-    circuit = models.Circuit(n_qubits)
-    # Add the Givens rotation circuits
-    for g_rot_parameter in g_rotation_parameters:
-        for orb1, orb2, theta, _phi in g_rot_parameter:
-            assert np.allclose(_phi, 0.0), "Unitary rotation is not real!"
-            circuit += givens_rotation_gate(n_qubits, orb1, orb2, theta)
-    return circuit
+    N = len(A[0])
+    gate_list = []
+    ordered_angles = []
+    #
+    for j in range(N):
+        for i in range(N):
+            if A[i][j] == 0:
+                # print("CRY", i, i + 1, A[i+1][j]-1, z_array[A[i + 1][j] - 1])
+                gate_list.append(gates.CNOT(i + 1, i))
+                gate_list.append(gates.CRY(i, i + 1, z_array[A[i + 1][j] - 1]))
+                ordered_angles.append(z_array[A[i + 1][j] - 1])
+                gate_list.append(gates.CNOT(i + 1, i))
+
+    return gate_list, ordered_angles
diff --git a/src/qibochem/ansatz/ucc.py b/src/qibochem/ansatz/ucc.py
index 9888be1..df38bbc 100644
--- a/src/qibochem/ansatz/ucc.py
+++ b/src/qibochem/ansatz/ucc.py
@@ -10,13 +10,13 @@
 
 
 def expi_pauli(n_qubits, theta, pauli_string):
-    """
+    r"""
     Build circuit representing exp(i*theta*pauli_string)
 
     Args:
         n_qubits: No. of qubits in the quantum circuit
         theta: parameter
-        pauli_string: OpenFermion QubitOperator object, e.g. X_0 Y_1 Z_2 X_4
+        pauli_string: OpenFermion QubitOperator object, e.g. :math: `X_0 Y_1 Z_2 X_4`
 
     Returns:
         circuit: Qibo Circuit object representing exp(i*theta*pauli_string)
@@ -119,9 +119,9 @@ def ucc_circuit(n_qubits, excitation, theta=0.0, trotter_steps=1, ferm_qubit_map
 
 
 def mp2_amplitude(excitation, orbital_energies, tei):
-    """
+    r"""
     Calculate the MP2 guess amplitude for a single UCC circuit: 0.0 for a single excitation.
-        for a double excitation (In SO basis): t_{ij}^{ab} = (g_{ijab} - g_{ijba}) / (e_i + e_j - e_a - e_b)
+        for a double excitation (In SO basis): :math:`t_{ij}^{ab} = (g_{ijab} - g_{ijba}) / (e_i + e_j - e_a - e_b)`
 
     Args:
         excitation: Iterable of spin-orbitals representing a excitation. Must have either 2 or 4 elements exactly,
diff --git a/tests/test_basis_rotation.py b/tests/test_basis_rotation.py
index c046e53..8f4dc8a 100644
--- a/tests/test_basis_rotation.py
+++ b/tests/test_basis_rotation.py
@@ -4,98 +4,142 @@
 
 import numpy as np
 import pytest
-from qibo import Circuit, gates
-from qibo.optimizers import optimize
+from qibo import Circuit, gates, models
 
 from qibochem.ansatz import basis_rotation
 from qibochem.driver.molecule import Molecule
 from qibochem.measurement.expectation import expectation
 
+# from qibo.optimizers import optimize
 
-def test_givens_rotation_gate():
-    n_qubits = 2
-    orb1 = 0
-    orb2 = 1
-    theta = -0.1
-    circuit = basis_rotation.givens_rotation_gate(n_qubits, orb1, orb2, theta)
-    ref_u = np.array(
+
+def test_unitary():
+    """
+    Test for basis_rotation.unitary()
+    """
+    N = 6
+    occ = range(0, 2)
+    vir = range(2, 6)
+
+    preset_params = [-0.1, -0.2, -0.3, -0.4]
+    toomany_params = [0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.11, 0.12]
+
+    U1, theta1 = basis_rotation.unitary(occ, vir)
+    U2, theta2 = basis_rotation.unitary(occ, vir, parameters=0.1)
+    U3, theta3 = basis_rotation.unitary(occ, vir, parameters=preset_params)
+
+    ref_U2 = np.array(
         [
-            [
-                -1.0,
-                0.0,
-                0.0,
-                0.0,
-            ],
-            [0.0, 0.99500417, 0.09983342, 0.0],
-            [0.0, -0.09983342, 0.99500417, 0.0],
-            [0.0, 0.0, 0.0, -1.0],
+            [0.99001666, 0.0, 0.099667, 0.0, 0.099667, 0.0],
+            [0.0, 0.99001666, 0.0, 0.099667, 0.0, 0.099667],
+            [-0.099667, 0.0, 0.99500833, 0.0, -0.00499167, 0.0],
+            [0.0, -0.099667, 0.0, 0.99500833, 0.0, -0.00499167],
+            [-0.099667, 0.0, -0.00499167, 0.0, 0.99500833, 0.0],
+            [0.0, -0.099667, 0.0, -0.00499167, 0.0, 0.99500833],
         ]
     )
 
-    assert np.allclose(circuit.unitary(), ref_u)
-
-
-def test_givens_rotation_parameters():
-    u = np.array(
-        [[0.99001666, 0.099667, 0.099667], [-0.099667, 0.99500833, -0.00499167], [-0.099667, -0.00499167, 0.99500833]]
+    ref_U3 = np.array(
+        [
+            [0.95041528, 0.0, -0.09834165, 0.0, -0.29502494, 0.0],
+            [0.0, 0.9016556, 0.0, -0.19339968, 0.0, -0.38679937],
+            [0.09834165, 0.0, 0.99504153, 0.0, -0.01487542, 0.0],
+            [0.0, 0.19339968, 0.0, 0.98033112, 0.0, -0.03933776],
+            [0.29502494, 0.0, -0.01487542, 0.0, 0.95537375, 0.0],
+            [0.0, 0.38679937, 0.0, -0.03933776, 0.0, 0.92132448],
+        ]
     )
-    n_occ = 1
-    n_qubits = 3
-    params = basis_rotation.givens_rotation_parameters(n_qubits, u, n_occ)
-    ref_params = [((0, 1, -1.4709635780470989, 0.0),), ((1, 2, 1.4704623293305714, 0.0),)]
 
-    assert np.allclose(params, ref_params)
+    identity = np.eye(6)
+
+    assert np.allclose(U1 @ U1.T, identity)
+    assert np.allclose(U2 @ U2.T, identity)
+    assert np.allclose(U3 @ U3.T, identity)
+    assert np.allclose(U1, identity)
+    assert np.allclose(U2, ref_U2)
+    assert np.allclose(U3, ref_U3)
 
+    with pytest.raises(IndexError):
+        Ux, thetax = basis_rotation.unitary(occ, vir, parameters=toomany_params)
 
-def test_swap_matrices():
-    """Test for swap_matrices"""
-    permutations = [[(0, 1), (2, 3)]]
-    n_qubits = 4
-    smat = basis_rotation.swap_matrices(permutations, n_qubits)
-    ref_smat = np.array([[0.0, 1.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0]])
 
-    assert np.allclose(smat, ref_smat)
+def test_givens_qr_decompose():
+    """
+    Test for basis_rotation.givens_qr_decompose()
+    """
+    N = 6
+    occ = range(0, 2)
+    vir = range(2, 6)
 
+    U2, theta2 = basis_rotation.unitary(occ, vir, parameters=0.1)
+    z_angles, final_U2 = basis_rotation.givens_qr_decompose(U2)
 
-def test_unitary_rot_matrix():
-    """Test for unitary rotation matrix"""
-    occ = [0]
-    vir = [1, 2]
-    parameters = np.zeros(len(occ) * len(vir))
-    parameters += 0.1
-    u = basis_rotation.unitary_rot_matrix(parameters, occ, vir)
-    ref_u = np.array(
-        [[0.99001666, 0.099667, 0.099667], [-0.099667, 0.99500833, -0.00499167], [-0.099667, -0.00499167, 0.99500833]]
+    ref_z = np.array(
+        [
+            -3.141592653589793,
+            -1.5707963267948966,
+            -2.356194490192345,
+            0.0,
+            -1.5707963267948966,
+            -1.5207546393123066,
+            -1.5707963267948966,
+            -1.5707963267948954,
+            -3.000171297352484,
+            -2.356194490192345,
+            0.0,
+            -1.5707963267948966,
+            -1.5707963267948954,
+            -0.09995829685982476,
+            -1.5207546393123068,
+        ]
     )
 
-    assert np.allclose(u, ref_u)
+    assert np.allclose(z_angles, ref_z)
+    assert np.allclose(np.eye(6), final_U2)
+
+
+def test_basis_rotation_layout():
+
+    N = 10
+    ref_A = np.array(
+        [
+            [0, -1, 0, -1, 0, -1, 0, -1, 0, -1],
+            [1, 0, 6, 0, 15, 0, 28, 0, 45, 0],
+            [0, 5, 0, 14, 0, 27, 0, 44, 0, 29],
+            [4, 0, 13, 0, 26, 0, 43, 0, 30, 0],
+            [0, 12, 0, 25, 0, 42, 0, 31, 0, 16],
+            [11, 0, 24, 0, 41, 0, 32, 0, 17, 0],
+            [0, 23, 0, 40, 0, 33, 0, 18, 0, 7],
+            [22, 0, 39, 0, 34, 0, 19, 0, 8, 0],
+            [0, 38, 0, 35, 0, 20, 0, 9, 0, 2],
+            [37, -1, 36, -1, 21, -1, 10, -1, 3, -1],
+        ]
+    )
 
+    A = basis_rotation.basis_rotation_layout(N)
+    assert np.allclose(A, ref_A)
 
-def test_br_ansatz():
-    """Test of basis rotation ansatz against hardcoded HF energies"""
-    h2_ref_energy = -1.117349035
 
-    mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
-    mol.run_pyscf(max_scf_cycles=1)
-    # Use an un-converged wave function to build the Hamiltonian
-    mol_sym_ham = mol.hamiltonian()
+def test_basis_rotation_gates():
 
-    # Define quantum circuit
-    circuit = Circuit(mol.nso)
-    circuit.add(gates.X(_i) for _i in range(mol.nelec))
+    mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.9)), ("H", (0.0, 0.0, 1.8))], 1, 1)
+    mol.run_pyscf(max_scf_cycles=100)
+    ham = mol.hamiltonian("sym", mol.oei, mol.tei, 0.0, "jw")
 
-    # Add basis rotation ansatz
-    # Initialize all at zero
-    parameters = np.zeros(mol.nelec * (mol.nso - mol.nelec))  # n_occ * n_virt
-    circuit += basis_rotation.br_circuit(mol.nso, parameters, mol.nelec)
+    nqubits = mol.nso
+    occ = range(0, mol.nelec)
+    vir = range(mol.nelec, mol.nso)
 
-    def electronic_energy(parameters):
-        """
-        Loss function (Electronic energy) for the basis rotation ansatz
-        """
-        circuit.set_parameters(parameters)
-        return expectation(circuit, mol_sym_ham)
+    U, theta = basis_rotation.unitary(occ, vir, parameters=0.1)
+    gate_angles, final_U = basis_rotation.givens_qr_decompose(U)
+    gate_layout = basis_rotation.basis_rotation_layout(nqubits)
+    gate_list, qubit_parameters = basis_rotation.basis_rotation_gates(gate_layout, gate_angles, theta)
 
-    hf_energy, parameters, _extra = optimize(electronic_energy, parameters)
+    circuit = Circuit(nqubits)
+    for _i in range(mol.nelec):
+        circuit.add(gates.X(_i))
+    circuit.add(gate_list)
 
-    assert hf_energy == pytest.approx(h2_ref_energy)
+    vqe = models.VQE(circuit, ham)
+    res = vqe.minimize(qubit_parameters)
+    assert np.isclose(res[0], mol.e_hf, atol=0.1)