diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index e6b6ad98..d4262ba0 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -20,6 +20,6 @@ repos:
       - id: isort
         args: ["--profile", "black"]
   - repo: https://github.com/asottile/pyupgrade
-    rev: v3.19.0
+    rev: v3.19.1
     hooks:
       - id: pyupgrade
diff --git a/README.md b/README.md
index c90b320e..d0ad654f 100644
--- a/README.md
+++ b/README.md
@@ -53,7 +53,11 @@ If you use the Qibochem plugin please refer to the documentation for citation in
 
 ## Contact
 
-For questions, comments and suggestions please contact us at [https://matrix.to/#/#qibo:matrix.org](url)
+To get in touch with the community and the developers, consider joining the Qibo workspace on Matrix:
+
+[![Matrix](https://img.shields.io/matrix/qibo%3Amatrix.org?logo=matrix)](https://matrix.to/#/#qibo:matrix.org)
+
+If you have a question about the project, contact us at [📫](mailto:qiboteam@qibo.science).
 
 ## Contributing
 
diff --git a/doc/source/api-reference/ansatz.rst b/doc/source/api-reference/ansatz.rst
index b123103a..c61ffa40 100644
--- a/doc/source/api-reference/ansatz.rst
+++ b/doc/source/api-reference/ansatz.rst
@@ -30,3 +30,15 @@ Basis rotation
 --------------
 
 .. autofunction:: qibochem.ansatz.basis_rotation.basis_rotation_gates
+
+Givens Excitation
+-----------------
+
+.. autofunction:: qibochem.ansatz.givens_excitation.givens_excitation_circuit
+
+.. autofunction:: qibochem.ansatz.givens_excitation.givens_excitation_ansatz
+
+Symmetry-Preserving
+-------------------
+
+.. autofunction:: qibochem.ansatz.symmetry.symm_preserving_circuit
diff --git a/doc/source/appendix/citing-qibochem.rst b/doc/source/appendix/citing-qibochem.rst
index b45c97a6..1a94e4a0 100644
--- a/doc/source/appendix/citing-qibochem.rst
+++ b/doc/source/appendix/citing-qibochem.rst
@@ -1,22 +1,20 @@
-Publications
-============
-
-If Qibochem has been significant in your research, and you would like to acknowledge
-the project in your academic publication, we suggest citing the following documents:
-
-Software References in Zenodo
------------------------------
-
-Wong Zicheng, Adrian Mak, Tan Le, Stefano Carrazza, Alessandro Candido, & Ye Jun. (2024). qiboteam/qibochem: Qibochem 0.0.1 (v0.0.1). Zenodo. https://doi.org/10.5281/zenodo.10473173
-
-Authorship Guideline
---------------------
-
-In order to appear as an author of a Qibochem publication (paper, proceedings, etc)
-each author must fullfil the following requirements:
-
-* Participate to the official meetings.
-
-* Contribute to the code with documented commits.
-
+Publications
+============
+
+If Qibochem has been significant in your research, and you would like to acknowledge the project in your academic publication, we suggest citing the following reference:
+
+Software References in Zenodo
+-----------------------------
+
+Wong Zi Cheng, Adrian Matthew Mak, Tan Le, Stefano Carrazza, Alessandro Candido & Ye Jun. (2024). qiboteam/qibochem: Qibochem 0.0.1 (v0.0.1). Zenodo. https://doi.org/10.5281/zenodo.10473173
+
+Authorship Guideline
+--------------------
+
+In order to appear as an author of a Qibochem publication (paper, proceedings, etc), one must fullfil the following requirements:
+
+* Participate to the official meetings.
+
+* Contribute to the code with documented commits.
+
 * Contribute to the manuscript elaboration.
diff --git a/doc/source/tutorials/adaptive.rst b/doc/source/tutorials/adaptive.rst
new file mode 100644
index 00000000..f3189f27
--- /dev/null
+++ b/doc/source/tutorials/adaptive.rst
@@ -0,0 +1,366 @@
+Adaptive methods
+================
+
+Despite VQE being cheaper than QPE, circuit depth is still a big problem for today's quantum hardware.
+
+.. code-block:: python
+
+    from qibochem.driver import Molecule
+    from qibochem.ansatz import ucc_ansatz
+
+    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))])
+    mol.run_pyscf()
+    mol.hf_embedding(active=[0, 1, 2, 5])
+    circuit = ucc_ansatz(mol)
+    print(circuit.summary())
+
+Output:
+
+.. code-block:: output
+
+    Circuit depth = 1874
+    Total number of gates = 3300
+    Number of qubits = 8
+    Most common gates:
+    cx: 1312
+    h: 1216
+    sdg: 304
+    s: 304
+    rz: 160
+    x: 4
+
+As shown in the above code block, the full UCCSD circuit for a simplified LiH/STO-3G system has a circuit depth of 1874 (!), with more than 1000 CNOT gates required!
+Hence, there is still a need to further reduce and simplify the circuit ansatzes used for running a VQE simulation.
+
+Other than designing shorter and more efficient circuit ansatzes, one alternative approach is through the use of energy gradients - for instance, through the Parameter Shift Rule on hardware - to filter and reduce the number of fermionic excitations in a circuit ansatz. [#f1]_ [#f2]_
+This is known as an adaptive method, in the sense that the quantum gates used to construct the circuit ansatz, as well as its actual structure and arrangement is not fixed, and varies depending on the molecular system under study.
+
+For example, in a H2/STO-3G system mapped with the Jordan-Wigner transformation, there are three possible spin-allowed fermionic excitations:
+two single excitations (``[0, 2]``, ``[1, 3]``) and one double excitation (``[0, 1, 2, 3]``).
+The full UCCSD circuit for this system has been shown in an earlier :ref:`example <UCC Ansatz>`, and it requires 64 CNOT gates for this simple molecular system.
+
+Let's look at the gradients of each of the individual fermionic excitations:
+
+.. code-block:: python
+
+    from qibo.derivative import parameter_shift
+
+    from qibochem.driver import Molecule
+    from qibochem.ansatz import hf_circuit, ucc_circuit
+
+    mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    mol.run_pyscf()
+    hamiltonian = mol.hamiltonian()
+
+    excitations = [[0, 1, 2, 3], [0, 2], [1, 3]]
+    circuit = hf_circuit(4, 2)
+    for excitation in excitations:
+        _circuit = circuit.copy()
+        _circuit += ucc_circuit(4, excitation)
+        n_parameters = len(_circuit.get_parameters())
+        gradients = [round(parameter_shift(_circuit, hamiltonian, index), 5) for index in range(n_parameters)]
+        print(f"Energy gradients for {excitation}: {gradients}")
+
+Output:
+
+.. code-block:: output
+
+    Energy gradients for [0, 1, 2, 3]: [0.179, -0.179, -0.179, -0.179, 0.179, 0.179, 0.179, -0.179]
+    Energy gradients for [0, 2]: [0.0, 0.0]
+    Energy gradients for [1, 3]: [0.0, 0.0]
+
+Only the doubles excitation has a non-zero energy gradient, which follows Brillouin's theorem.
+Running the VQE with only the doubles excitation then gives:
+
+.. code-block:: python
+
+    import numpy as np
+
+    from qibo.models import VQE
+
+    from qibochem.driver import Molecule
+    from qibochem.ansatz import hf_circuit, ucc_circuit
+
+    mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    mol.run_pyscf()
+    hamiltonian = mol.hamiltonian()
+
+    circuit = hf_circuit(4, 2)
+    circuit += ucc_circuit(4, [0, 1, 2, 3])
+
+    vqe = VQE(circuit, hamiltonian)
+
+    # Optimize starting from a random guess for the variational parameters
+    initial_parameters = np.random.uniform(0, 2*np.pi, len(circuit.get_parameters()))
+    best, params, extra = vqe.minimize(initial_parameters, method='BFGS', compile=False)
+
+    # Exact result
+    print(f"Exact result: {mol.eigenvalues(hamiltonian)[0]:.7f}")
+    print(f"  VQE result: {best:.7f}")
+
+Output:
+
+.. code-block:: output
+
+    Exact result: -1.1361895
+      VQE result: -1.1361895
+
+We managed to find the exact result by applying only the doubles excitation!
+
+Next, let's look at the potential savings for the simplified LiH/STO-3G system.
+To reduce the circuit depth further, we will use the more modern ansatz, the Givens excitation circuit from Arrazola et al., [#f1]_ instead of the UCC ansatz.
+
+As was done in the above example, we will start with a HF circuit, then find the gradients for each circuit ansatz corresponding to a fermionic excitation.
+After that, the excitation with the largest absolute value of the gradient will be added to the initial circuit, followed by a VQE simulation.
+
+.. code-block:: python
+
+    from qibo.derivative import parameter_shift
+    from qibo.models import VQE
+
+    from qibochem.driver import Molecule
+    from qibochem.ansatz import hf_circuit, givens_excitation_circuit, generate_excitations, sort_excitations
+
+    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))])
+    mol.run_pyscf()
+    mol.hf_embedding(active=[0, 1, 2, 5])
+    hamiltonian = mol.hamiltonian()
+
+    n_qubits = mol.n_active_orbs
+    n_elec = mol.n_active_e
+
+    circuit = hf_circuit(n_qubits, n_elec)
+
+    excitations = sort_excitations(generate_excitations(2, list(range(n_elec)), list(range(n_elec, n_qubits))))
+    excitations += sort_excitations(generate_excitations(1, list(range(n_elec)), list(range(n_elec, n_qubits))))
+
+    excitation_gradients = {}
+    for excitation in excitations:
+        _circuit = circuit.copy()
+        _circuit += givens_excitation_circuit(n_qubits, excitation)
+        n_parameters = len(_circuit.get_parameters())
+        gradient = [round(parameter_shift(_circuit, hamiltonian, index), 5) for index in range(n_parameters)]
+        print(f"Energy gradients for {excitation}: {gradient}")
+        excitation_gradients[tuple(excitation)] = gradient[0] # Gradient magnitude is equal throughout
+
+    # Find the excitation corresponding to the largest gradient, and add it to the circuit
+    max_grad = max(excitation_gradients, key=lambda x: abs(excitation_gradientis.get(x)))
+    print(f"\nExcitation with the largest gradient: {max_grad}; Gradient = {excitation_gradients[max_grad]}")
+    circuit += givens_excitation_circuit(n_qubits, max_grad)
+
+    # Run VQE with the updated circuit
+    vqe = VQE(circuit, hamiltonian)
+
+    circuit_parameters = [param for _tuple in circuit.get_parameters() for param in _tuple]
+    best, params, extra = vqe.minimize(circuit_parameters, method='BFGS', compile=False)
+
+    print(f" HF energy: {mol.e_hf:.7f}")
+    print(f"VQE result: {best:.7f}")
+
+Output:
+
+.. code-block:: output
+
+    Energy gradients for [0, 1, 4, 5]: [0.02132, -0.02132, 0.02132, -0.02132, -0.02132, 0.02132, -0.02132, 0.02132]
+    Energy gradients for [0, 1, 6, 7]: [0.00569, -0.00569, 0.00569, -0.00569, -0.00569, 0.00569, -0.00569, 0.00569]
+    Energy gradients for [2, 3, 4, 5]: [0.01136, -0.01136, 0.01136, -0.01136, -0.01136, 0.01136, -0.01136, 0.01136]
+    Energy gradients for [2, 3, 6, 7]: [0.12225, -0.12225, 0.12225, -0.12225, -0.12225, 0.12225, -0.12225, 0.12225]
+    Energy gradients for [0, 1, 4, 7]: [0.00016, -0.00016, 0.00016, -0.00016, -0.00016, 0.00016, -0.00016, 0.00016]
+    Energy gradients for [0, 1, 5, 6]: [-0.00016, 0.00016, -0.00016, 0.00016, 0.00016, -0.00016, 0.00016, -0.00016]
+    Energy gradients for [2, 3, 4, 7]: [-0.03254, 0.03254, -0.03254, 0.03254, 0.03254, -0.03254, 0.03254, -0.03254]
+    Energy gradients for [2, 3, 5, 6]: [0.03254, -0.03254, 0.03254, -0.03254, -0.03254, 0.03254, -0.03254, 0.03254]
+    Energy gradients for [0, 3, 4, 5]: [0.00029, -0.00029, 0.00029, -0.00029, -0.00029, 0.00029, -0.00029, 0.00029]
+    Energy gradients for [1, 2, 4, 5]: [-0.00029, 0.00029, -0.00029, 0.00029, 0.00029, -0.00029, 0.00029, -0.00029]
+    Energy gradients for [0, 3, 6, 7]: [0.00108, -0.00108, 0.00108, -0.00108, -0.00108, 0.00108, -0.00108, 0.00108]
+    Energy gradients for [1, 2, 6, 7]: [-0.00108, 0.00108, -0.00108, 0.00108, 0.00108, -0.00108, 0.00108, -0.00108]
+    Energy gradients for [0, 2, 4, 6]: [0.00299, -0.00299, 0.00299, -0.00299, -0.00299, 0.00299, -0.00299, 0.00299]
+    Energy gradients for [1, 3, 5, 7]: [0.00299, -0.00299, 0.00299, -0.00299, -0.00299, 0.00299, -0.00299, 0.00299]
+    Energy gradients for [0, 3, 4, 7]: [-0.00236, 0.00236, -0.00236, 0.00236, 0.00236, -0.00236, 0.00236, -0.00236]
+    Energy gradients for [0, 3, 5, 6]: [-0.00063, 0.00063, -0.00063, 0.00063, 0.00063, -0.00063, 0.00063, -0.00063]
+    Energy gradients for [1, 2, 4, 7]: [-0.00063, 0.00063, -0.00063, 0.00063, 0.00063, -0.00063, 0.00063, -0.00063]
+    Energy gradients for [1, 2, 5, 6]: [-0.00236, 0.00236, -0.00236, 0.00236, 0.00236, -0.00236, 0.00236, -0.00236]
+    Energy gradients for [0, 4]: [0.0, -0.0]
+    Energy gradients for [1, 5]: [-0.0, 0.0]
+    Energy gradients for [0, 6]: [0.0, -0.0]
+    Energy gradients for [1, 7]: [-0.0, 0.0]
+    Energy gradients for [2, 4]: [-0.0, 0.0]
+    Energy gradients for [3, 5]: [0.0, -0.0]
+    Energy gradients for [2, 6]: [-0.0, 0.0]
+    Energy gradients for [3, 7]: [0.0, -0.0]
+
+    Excitation with the largest gradient: (2, 3, 6, 7); Gradient = 0.12225
+     HF energy: -7.8605387
+    VQE result: -7.8732886
+
+
+After adding the circuit ansatz corresponding to one double excitation and running VQE,
+the resultant energy was found to be about 0.01 Hartrees lower compared to the bare HF ansatz.
+Clearly, there is still room for further improvement in the obtained energies.
+
+To do so, we further extend the circuit ansatz by adding in more excitations iteratively,
+with the excitation with the largest (in magnitute) energy gradients added on at each step.
+This can be carried out until the difference between each iteration is small (<0.001 Hartrees), or until there are no more remaining excitations to be added.
+
+.. warning::
+
+      The code block below might take a few minutes to run!
+
+.. code-block:: python
+
+    import numpy as np
+
+    from qibo.models import VQE
+    from qibo.derivative import parameter_shift
+
+    from qibochem.driver import Molecule
+    from qibochem.ansatz import hf_circuit, givens_excitation_circuit, generate_excitations, sort_excitations
+
+    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))])
+    mol.run_pyscf()
+    mol.hf_embedding(active=[0, 1, 2, 5])
+    hamiltonian = mol.hamiltonian()
+
+    exact_result = mol.eigenvalues(hamiltonian)[0]
+
+    n_qubits = mol.n_active_orbs
+    n_elec = mol.n_active_e
+
+    circuit = hf_circuit(n_qubits, n_elec)
+
+    excitations = sort_excitations(generate_excitations(2, list(range(n_elec)), list(range(n_elec, n_qubits))))
+    excitations += sort_excitations(generate_excitations(1, list(range(n_elec)), list(range(n_elec, n_qubits))))
+
+    excitation_gradient = {tuple(excitation):0.0 for excitation in excitations}
+    count = 0
+    current_energy = mol.e_hf
+    n_fixed_params = 0
+    # Iterating through all excitations; loop breaks if difference in VQE result between excitations is v small
+    while excitations:
+        print(f"Iteration {count+1}:")
+        for excitation in excitations:
+            _circuit = circuit.copy()
+            _circuit += givens_excitation_circuit(n_qubits, excitation)
+            n_parameters = len(_circuit.get_parameters())
+            gradient = [round(parameter_shift(_circuit, hamiltonian, index), 5) for index in range(n_parameters)][n_fixed_params:]
+            # print(f"Energy gradient for {excitation}: {gradient}")
+            excitation_gradient[tuple(excitation)] = gradient[0] # Gradient magnitude is equal throughout
+
+        # Find the excitation corresponding to the largest gradient, and add it to the circuit
+        max_grad = max(excitation_gradient, key=lambda x: abs(excitation_gradient.get(x)))
+        print(f"Excitation with the largest gradient: {max_grad}; Gradient = {excitation_gradient[max_grad]}")
+        circuit += givens_excitation_circuit(n_qubits, max_grad)
+        # Remove max_grad from excitations and excitation_data
+        excitations.pop(excitations.index(list(max_grad)))
+        del excitation_gradient[max_grad]
+
+        # Run VQE with the updated circuit
+        vqe = VQE(circuit, hamiltonian)
+
+        circuit_parameters = [param for _tuple in circuit.get_parameters() for param in _tuple]
+        best, params, extra = vqe.minimize(circuit_parameters, method='BFGS', compile=False)
+
+        n_fixed_params = len(params)
+
+        print(f"\nExact result: {exact_result:.7f}")
+        print(f"  VQE result: {best:.7f}")
+        energy_diff = best - current_energy
+        print(f"Difference to previous result: {energy_diff:.7f}")
+
+        if abs(energy_diff) < 1e-3:
+            print("\nEnergy has converged; exiting while loop")
+            break
+        print()
+        # Update circuit parameters and current energy
+        circuit.set_parameters(params)
+        current_energy = best
+        count += 1
+
+
+    print("\nFinal circuit:")
+    print(circuit.draw())
+    print("\nCircuit statistics:")
+    print(circuit.summary())
+
+Output:
+
+.. code-block:: output
+
+    Iteration 1:
+    Excitation with the largest gradient: (2, 3, 6, 7); Gradient = 0.12225
+
+    Exact result: -7.8770974
+      VQE result: -7.8732886
+    Difference to previous result: -0.0127499
+
+    Iteration 2:
+    Excitation with the largest gradient: (2, 3, 4, 7); Gradient = -0.03485
+
+    Exact result: -7.8770974
+      VQE result: -7.8748417
+    Difference to previous result: -0.0015531
+
+    Iteration 3:
+    Excitation with the largest gradient: (2, 3, 5, 6); Gradient = 0.03364
+
+    Exact result: -7.8770974
+      VQE result: -7.8762910
+    Difference to previous result: -0.0014493
+
+    Iteration 4:
+    Excitation with the largest gradient: (0, 1, 4, 5); Gradient = 0.02124
+
+    Exact result: -7.8770974
+      VQE result: -7.8763762
+    Difference to previous result: -0.0000853
+
+    Energy has converged; exiting while loop
+
+    Final circuit:
+    q0:     ─X──────────────────────────────────────────────────────────────────── ...
+    q1:     ─X──────────────────────────────────────────────────────────────────── ...
+    q2:     ─X───o─H─o───RY─o─────RY───X─RY─────o─RY─o─X─H─o─────o─H─o───RY─o───── ...
+    q3:     ─X───|───X───RY─|───X─RY─X─|─RY─X───|─RY─X─|───|─────|───X───RY─|───X─ ...
+    q4:     ─────|──────────|───|────|─|────|───|──────|───|───o─X─────o────|───|─ ...
+    q5:     ─────|──────────|───|────|─|────|───|──────|───|───|───────|────|───|─ ...
+    q6:     ───o─X─────o────|───|────o─o────|───|──────o───X─o─|───────|────|───|─ ...
+    q7:     ───X───H───X────X─H─o───────────o─H─X────────H───X─X───H───X────X─H─o─ ...
+
+    q0: ... ────────────────────────────────────────────────────────────────────── ...
+    q1: ... ────────────────────────────────────────────────────────────────────── ...
+    q2: ... RY───X─RY─────o─RY─o─X─H─o─────o─H─o───RY─o─────RY───X─RY─────o─RY─o─X ...
+    q3: ... RY─X─|─RY─X───|─RY─X─|───|─────|───X───RY─|───X─RY─X─|─RY─X───|─RY─X─| ...
+    q4: ... ───o─o────|───|──────o───X─o───|──────────|───|────|─|────|───|──────| ...
+    q5: ... ──────────|───|────────────|─o─X─────o────|───|────o─o────|───|──────o ...
+    q6: ... ──────────|───|────────────|─X───H───X────X─H─o───────────o─H─X─────── ...
+    q7: ... ──────────o─H─X────────H───X────────────────────────────────────────── ...
+
+    q0: ... ─────────o─H─o───RY─o─────RY───X─RY─────o─RY─o─X─H─o───
+    q1: ... ─────────|───X───RY─|───X─RY─X─|─RY─X───|─RY─X─|───|───
+    q2: ... ─H─o─────|──────────|───|────|─|────|───|──────|───|───
+    q3: ... ───|─────|──────────|───|────|─|────|───|──────|───|───
+    q4: ... ───|───o─X─────o────|───|────o─o────|───|──────o───X─o─
+    q5: ... ───X─o─X───H───X────X─H─o───────────o─H─X────────H───X─
+    q6: ... ─H───X─────────────────────────────────────────────────
+    q7: ... ───────────────────────────────────────────────────────
+
+    Circuit statistics:
+    Circuit depth = 78
+    Total number of gates = 116
+    Number of qubits = 8
+    Most common gates:
+    cx: 56
+    ry: 32
+    h: 24
+    x: 4
+
+Recall that the full UCCSD circuit for our system had a circuit depth of 1874, with more than 1000 CNOT gates required.
+In contrast, the use of a simpler circuit ansatz in conjunction with an adaptive approach allowed us to find a VQE energy that is within chemical accuracy,
+while using only 56 CNOT gates and with a final gate depth of only 78.
+This is a >20-fold reduction in the gate depth and number of CNOT gates used!
+
+
+  .. rubric:: References
+
+  .. [#f1] Arrazola, J. M. et al., 'Universal quantum circuits for quantum chemistry', Quantum, 6, (2022), 742.
+
+  .. [#f2] Schuld M. et al., 'Evaluating analytic gradients on quantum hardware', Phys. Rev. A, 99, (2019), 032331.
diff --git a/doc/source/tutorials/index.rst b/doc/source/tutorials/index.rst
index 4ce96f7f..3617a44c 100644
--- a/doc/source/tutorials/index.rst
+++ b/doc/source/tutorials/index.rst
@@ -12,3 +12,4 @@ This section provides examples of how to use Qibochem.
     hamiltonian
     ansatz
     measurement
+    adaptive
diff --git a/doc/source/tutorials/measurement.rst b/doc/source/tutorials/measurement.rst
index 58c18b9d..9bc817d6 100644
--- a/doc/source/tutorials/measurement.rst
+++ b/doc/source/tutorials/measurement.rst
@@ -42,9 +42,12 @@ In practice, the expectation value for each of the individual Pauli terms have t
 This process of obtaining the electronic energy (Hamiltonian expectation value) is still reasonable for a small system.
 However, the number of Pauli terms in a molecular Hamiltonian scales on the order of :math:`O(N^4)`, where N is the number of qubits.
 
+.. warning::
+
+    The code block below might take a few minutes to run!
+
 .. code-block:: python
 
-    # Warning: This code block might take a few minutes to run
     from qibochem.driver import Molecule
 
     # Build the N2 molecule and get the molecular Hamiltonian
diff --git a/poetry.lock b/poetry.lock
index ce5e14eb..8bad7840 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -1,4 +1,4 @@
-# This file is automatically @generated by Poetry 1.8.3 and should not be changed by hand.
+# This file is automatically @generated by Poetry 1.8.5 and should not be changed by hand.
 
 [[package]]
 name = "alabaster"
@@ -1091,13 +1091,13 @@ testing = ["Django", "attrs", "colorama", "docopt", "pytest (<9.0.0)"]
 
 [[package]]
 name = "jinja2"
-version = "3.1.4"
+version = "3.1.5"
 description = "A very fast and expressive template engine."
 optional = false
 python-versions = ">=3.7"
 files = [
-    {file = "jinja2-3.1.4-py3-none-any.whl", hash = "sha256:bc5dd2abb727a5319567b7a813e6a2e7318c39f4f487cfe6c89c6f9c7d25197d"},
-    {file = "jinja2-3.1.4.tar.gz", hash = "sha256:4a3aee7acbbe7303aede8e9648d13b8bf88a429282aa6122a993f0ac800cb369"},
+    {file = "jinja2-3.1.5-py3-none-any.whl", hash = "sha256:aba0f4dc9ed8013c424088f68a5c226f7d6097ed89b246d7749c2ec4175c6adb"},
+    {file = "jinja2-3.1.5.tar.gz", hash = "sha256:8fefff8dc3034e27bb80d67c671eb8a9bc424c0ef4c0826edbff304cceff43bb"},
 ]
 
 [package.dependencies]
diff --git a/src/qibochem/ansatz/__init__.py b/src/qibochem/ansatz/__init__.py
index 8bd4b795..cb7d9802 100644
--- a/src/qibochem/ansatz/__init__.py
+++ b/src/qibochem/ansatz/__init__.py
@@ -1,13 +1,11 @@
 from qibochem.ansatz.basis_rotation import basis_rotation_gates
+from qibochem.ansatz.givens_excitation import (
+    givens_excitation_ansatz,
+    givens_excitation_circuit,
+)
 from qibochem.ansatz.hardware_efficient import he_circuit
 from qibochem.ansatz.hf_reference import hf_circuit
 from qibochem.ansatz.qeb import qeb_circuit
-from qibochem.ansatz.ucc import (
-    generate_excitations,
-    mp2_amplitude,
-    sort_excitations,
-    ucc_ansatz,
-    ucc_circuit,
-)
-
-# TODO: Probably can move some of the functions, e.g. generate_excitations/sort_excitations to a new 'util.py'
+from qibochem.ansatz.symmetry import symm_preserving_circuit
+from qibochem.ansatz.ucc import ucc_ansatz, ucc_circuit
+from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations
diff --git a/src/qibochem/ansatz/givens_excitation.py b/src/qibochem/ansatz/givens_excitation.py
new file mode 100644
index 00000000..e891d15b
--- /dev/null
+++ b/src/qibochem/ansatz/givens_excitation.py
@@ -0,0 +1,155 @@
+"""
+Circuit ansatz for representing a fermionic excitation as a Givens rotation by Arrazola et al.
+Reference: https://doi.org/10.22331/q-2022-06-20-742
+"""
+
+from qibo import Circuit, gates
+
+from qibochem.ansatz.hf_reference import hf_circuit
+from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations
+
+# Helper functions
+
+
+def single_excitation_gate(sorted_orbitals, theta):
+    """
+    Decomposition of a Givens single excitation gate into single qubit rotations and CNOTs. In principle, should be
+    identical to gates.GIVENS(qubit0, qubit1, theta)
+
+    Args:
+        sorted_orbitals (list): Sorted list of orbitals involved in the excitation
+        theta (float): Rotation angle
+
+    Returns:
+        (list): List of gates representing the decomposition of the Givens' single excitation gate
+    """
+    result = []
+    result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[1]))
+    result.append(gates.RY(sorted_orbitals[0], 0.5 * theta))
+    result.append(gates.CNOT(sorted_orbitals[1], sorted_orbitals[0]))
+    result.append(gates.RY(sorted_orbitals[0], -0.5 * theta))
+    result.append(gates.CNOT(sorted_orbitals[1], sorted_orbitals[0]))
+    result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[1]))
+    return result
+
+
+def double_excitation_gate(sorted_orbitals, theta):
+    """
+    Decomposition of a Givens double excitation gate into single qubit rotations and CNOTs
+
+    Args:
+        sorted_orbitals (list): Sorted list of orbitals involved in the excitation
+        theta (float): Rotation angle
+
+    Returns:
+        (list): List of gates representing the decomposition of the Givens' double excitation gate
+    """
+    result = []
+    result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[3]))
+    result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[2]))
+    result.append(gates.H(sorted_orbitals[0]))
+    result.append(gates.H(sorted_orbitals[3]))
+    result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[1]))
+    result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[3]))
+    result.append(gates.RY(sorted_orbitals[0], -0.125 * theta))
+    result.append(gates.RY(sorted_orbitals[1], 0.125 * theta))
+    result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[3]))
+    result.append(gates.H(sorted_orbitals[3]))
+    result.append(gates.CNOT(sorted_orbitals[3], sorted_orbitals[1]))
+    result.append(gates.RY(sorted_orbitals[0], -0.125 * theta))
+    result.append(gates.RY(sorted_orbitals[1], 0.125 * theta))
+    result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[1]))
+    result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[0]))
+    result.append(gates.RY(sorted_orbitals[0], 0.125 * theta))
+    result.append(gates.RY(sorted_orbitals[1], -0.125 * theta))
+    result.append(gates.CNOT(sorted_orbitals[3], sorted_orbitals[1]))
+    result.append(gates.H(sorted_orbitals[3]))
+    result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[3]))
+    result.append(gates.RY(sorted_orbitals[0], 0.125 * theta))
+    result.append(gates.RY(sorted_orbitals[1], -0.125 * theta))
+    result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[1]))
+    result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[0]))
+    result.append(gates.H(sorted_orbitals[0]))
+    result.append(gates.H(sorted_orbitals[3]))
+    result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[2]))
+    result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[3]))
+    return result
+
+
+# Main function
+def givens_excitation_circuit(n_qubits, excitation, theta=None):
+    """
+    Circuit ansatz corresponding to the Givens rotation excitation from Arrazola et al.
+
+    Args:
+        n_qubits: Number of qubits in the circuit
+        excitation: Iterable of orbitals involved in the excitation; must have an even number of elements
+            E.g. ``[0, 1, 2, 3]`` represents the excitation of electrons in orbitals ``(0, 1)`` to ``(2, 3)``
+        theta (float): Rotation angle. Default: 0.0
+
+    Returns:
+        Qibo ``Circuit``: Circuit ansatz for a single Givens excitation
+    """
+    sorted_orbitals = sorted(excitation)
+    # Check size of orbitals input
+    assert len(sorted_orbitals) % 2 == 0, f"{excitation} must have an even number of items"
+
+    if theta is None:
+        theta = 0.0
+
+    circuit = Circuit(n_qubits)
+    if len(excitation) == 2:
+        circuit.add(single_excitation_gate(sorted_orbitals, theta))
+    elif len(excitation) == 4:
+        circuit.add(double_excitation_gate(sorted_orbitals, theta))
+    else:
+        raise NotImplementedError("Can only handle single and double excitations!")
+    return circuit
+
+
+def givens_excitation_ansatz(
+    molecule,
+    excitations=None,
+    include_hf=True,
+    use_mp2_guess=True,
+):
+    """
+    Convenience function for buildng a circuit corresponding to the Givens excitation ansatz with multiple excitations
+    for a given ``Molecule``. If no excitations are given, it defaults to including all possible spin-allowed
+    excitations, up to doubles.
+
+    Args:
+        molecule: The ``Molecule`` of interest.
+        excitations: List of excitations (e.g. ``[[0, 1, 2, 3], [0, 1, 4, 5]]``) used to build the
+            UCC circuit. Overrides the ``excitation_level`` argument
+        include_hf: Whether or not to start the circuit with a Hartree-Fock circuit. Default: ``True``
+        use_mp2_guess: Whether to use MP2 amplitudes or a numpy zero array as the initial guess parameter. Default: ``True``;
+            use the MP2 amplitudes as the default guess parameters
+
+    Returns:
+        Qibo ``Circuit``: Circuit corresponding to a Givens excitation circuit ansatz
+    """
+    # TODO: Consolidate/Meld this function with the ucc_ansatz function; both are largely identical
+
+    # Get the number of electrons and spin-orbitals from the molecule argument
+    n_elec = molecule.nelec if molecule.n_active_e is None else molecule.n_active_e
+    n_orbs = molecule.nso if molecule.n_active_orbs is None else molecule.n_active_orbs
+
+    # If no excitations given, defaults to all possible double and single excitations
+    if excitations is None:
+        excitations = []
+        for order in range(2, 0, -1):  # Reversed to get double excitations first, then singles
+            excitations += sort_excitations(generate_excitations(order, range(0, n_elec), range(n_elec, n_orbs)))
+    else:
+        # Some checks to ensure the given excitations are valid
+        assert all(len(_ex) in (2, 4) for _ex in excitations), "Only single and double excitations allowed!"
+
+    # Build the circuit
+    if include_hf:
+        circuit = hf_circuit(n_orbs, n_elec)  # Only works with (default) JW mapping
+    else:
+        circuit = Circuit(n_orbs)
+    for excitation in excitations:
+        theta = mp2_amplitude(excitation, molecule.eps, molecule.tei) if use_mp2_guess else None
+        circuit += givens_excitation_circuit(n_orbs, excitation, theta)
+    return circuit
diff --git a/src/qibochem/ansatz/symmetry.py b/src/qibochem/ansatz/symmetry.py
new file mode 100644
index 00000000..2a4992e3
--- /dev/null
+++ b/src/qibochem/ansatz/symmetry.py
@@ -0,0 +1,105 @@
+"""
+Symmetry-Preserving circuit ansatz from Gard et al. Reference: https://doi.org/10.1038/s41534-019-0240-1
+"""
+
+from math import factorial
+
+import numpy as np
+from qibo import Circuit, gates
+
+# Helper functions
+
+
+def a_gate(qubit1, qubit2, theta=None, phi=None):
+    """
+    Decomposition of the 'A' gate as defined in the paper, acting on qubit1 and qubit2. 'A' corresponds to the following
+    unitary matrix:
+
+    A(\\theta, \\phi) =
+    \\begin{pmatrix}
+        1 & 0 & 0 & 0 \\\\
+        0 & \\cos \\theta & e^{i \\phi} \\sin \\theta & 0 \\\\
+        0 & e^{-i \\phi} \\sin \\theta & -\\cos \\theta & 0 \\\\
+        0 & 0 & 0 & 1
+    \\end{pmatrix}
+
+    Args:
+        qubit1 (int): Index of the first qubit
+        qubit2 (int): Index of the second qubit
+        theta (float): First rotation angle. Default: 0.0
+        phi (float): Second rotation angle. Default: 0.0
+
+    Returns:
+        (list): List of gates representing the decomposition of the 'A' gate
+    """
+    if theta is None:
+        theta = 0.0
+    if phi is None:
+        phi = 0.0
+
+    # R(theta, phi) = R_z (phi + pi) R_y (theta + 0.5*pi)
+    r_gate = [gates.RY(qubit2, theta + 0.5 * np.pi), gates.RZ(qubit2, phi + np.pi)]
+    r_gate_dagger = [_gate.dagger() for _gate in r_gate][::-1]
+    return (
+        [gates.CNOT(qubit2, qubit1)]
+        + r_gate_dagger
+        + [gates.CNOT(qubit1, qubit2)]
+        + r_gate
+        + [gates.CNOT(qubit2, qubit1)]
+    )
+
+
+def x_gate_indices(n_qubits, n_electrons):
+    """Obtain the qubit indices for X gates to be added to the circuit"""
+    indices = [2 * _i for _i in range(0, min(n_electrons, n_qubits // 2))]
+    if n_electrons > n_qubits // 2:
+        indices += [2 * _i + 1 for _i in range(n_electrons - (n_qubits // 2))]
+    return sorted(indices)
+
+
+def a_gate_indices(n_qubits, n_electrons, x_gates):
+    """
+    Obtain the qubit indices for a single layer of the primitive pattern of 'A' gates in the circuit ansatz
+    """
+    assert len(x_gates) == n_electrons, f"n_electrons ({n_electrons}) != Number of X gates given! ({x_gates})"
+    # 2. Apply 'first layer' of gates on all adjacent pairs of qubits on which either X*I or I*X has been applied.
+    first_layer = [(_i, _i + 1) for _i in x_gates if _i + 1 < n_qubits and _i + 1 not in x_gates]
+    first_layer += [(_i - 1, _i) for _i in x_gates if _i - 1 >= 0 and _i - 1 not in x_gates]
+    # 3a. Apply 'second layer' of gates on adjacent pairs of qubits. Each pair includes 1 qubit acted on in the previous
+    # step and a qubit free of gates. Continue placing gates on adjacent qubits until all neighboring qubits are connected
+    second_layer = [(_i, _i + 1) for _i in range(max(pair[1] for pair in first_layer), n_qubits - 1)]
+    second_layer += [(_i - 1, _i) for _i in range(min(pair[0] for pair in first_layer), 0, -1)]
+    # 3b. The first and second layers define a primitive pattern:
+    primitive_pattern = first_layer + second_layer
+    # Need to add any missing connections between neighbouring qubits
+    primitive_pattern += [pair for _i in range(n_qubits - 1) if (pair := (_i, _i + 1)) not in primitive_pattern]
+    # 4. Repeat the primitive pattern until (n_qubits choose n_electrons) A gates are placed
+    n_gates_per_layer = len(primitive_pattern)
+    n_a_gates = factorial(n_qubits) // (factorial(n_qubits - n_electrons) * factorial(n_electrons))
+    assert (
+        n_a_gates % n_gates_per_layer == 0
+    ), f"n_a_gates ({n_a_gates}) is not a multiple of n_gates_per_layer ({n_gates_per_layer})!"
+    return (n_a_gates // n_gates_per_layer) * primitive_pattern
+
+
+# Main function
+def symm_preserving_circuit(n_qubits, n_electrons):
+    """
+    Symmetry-preserving circuit ansatz from Gard et al. (https://doi.org/10.1038/s41534-019-0240-1)
+
+    Args:
+        n_qubits: Number of qubits in the quantum circuit
+        n_electrons: Number of electrons in the molecular system
+
+    Returns:
+        Qibo ``Circuit``: Circuit corresponding to the symmetry-preserving ansatz
+    """
+    circuit = Circuit(n_qubits)
+    x_gates = x_gate_indices(n_qubits, n_electrons)
+    circuit.add(gates.X(_i) for _i in x_gates)
+    # Generate the qubit pair indices for adding A gates
+    a_gate_qubits = a_gate_indices(n_qubits, n_electrons, x_gates)
+    a_gates = [a_gate(qubit1, qubit2) for qubit1, qubit2 in a_gate_qubits]
+    # Each a_gate is a list of elementary gates, so a_gates is a nested list; need to unpack it
+    circuit.add(_gates for _a_gate in a_gates for _gates in _a_gate)
+    return circuit
diff --git a/src/qibochem/ansatz/ucc.py b/src/qibochem/ansatz/ucc.py
index 4fde1841..99a9b78f 100644
--- a/src/qibochem/ansatz/ucc.py
+++ b/src/qibochem/ansatz/ucc.py
@@ -7,6 +7,7 @@
 from qibo import Circuit, gates
 
 from qibochem.ansatz.hf_reference import hf_circuit
+from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations
 
 
 def expi_pauli(n_qubits, pauli_string, theta):
@@ -108,160 +109,6 @@ def ucc_circuit(n_qubits, excitation, theta=0.0, trotter_steps=1, ferm_qubit_map
     return circuit
 
 
-# Utility functions
-
-
-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): :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,
-            representing a single or double excitation respectively.
-        orbital_energies: eigenvalues of the Fock operator, i.e. orbital energies
-        tei: Two-electron integrals in MO basis and second quantization notation
-
-    Returns:
-        MP2 guess amplitude (float)
-    """
-    # Checks validity of excitation argument
-    assert len(excitation) % 2 == 0 and len(excitation) // 2 <= 2, f"{excitation} must have either 2 or 4 elements"
-    # If single excitation, can just return 0.0 directly
-    if len(excitation) == 2:
-        return 0.0
-
-    # Convert orbital indices to be in MO basis
-    mo_orbitals = [orbital // 2 for orbital in excitation]
-    # Numerator: g_ijab - g_ijba
-    g_ijab = (
-        tei[tuple(mo_orbitals)]  # Can index directly using the MO TEIs
-        if (excitation[0] + excitation[3]) % 2 == 0 and (excitation[1] + excitation[2]) % 2 == 0
-        else 0.0
-    )
-    g_ijba = (
-        tei[tuple(mo_orbitals[:2] + mo_orbitals[2:][::-1])]  # Reverse last two terms
-        if (excitation[0] + excitation[2]) % 2 == 0 and (excitation[1] + excitation[3]) % 2 == 0
-        else 0.0
-    )
-    numerator = g_ijab - g_ijba
-    # Denominator is directly from the orbital energies
-    denominator = sum(orbital_energies[mo_orbitals[:2]]) - sum(orbital_energies[mo_orbitals[2:]])
-    return numerator / denominator
-
-
-def generate_excitations(order, excite_from, excite_to, conserve_spin=True):
-    """
-    Generate all possible excitations between a list of occupied and virtual orbitals
-
-    Args:
-        order: Order of excitations, i.e. 1 == single, 2 == double
-        excite_from: Iterable of integers
-        excite_to: Iterable of integers
-        conserve_spin: ensure that the total electronic spin is conserved
-
-    Return:
-        List of lists, e.g. [[0, 1]]
-    """
-    # If order of excitation > either number of electrons/orbitals, return list of empty list
-    if order > min(len(excite_from), len(excite_to)):
-        return [[]]
-
-    # Generate all possible excitations first
-    from itertools import combinations  # pylint: disable=C0415
-
-    all_excitations = [
-        [*_from, *_to] for _from in combinations(excite_from, order) for _to in combinations(excite_to, order)
-    ]
-    # Filter out the excitations if conserve_spin set
-    if conserve_spin:
-        # Not sure if this filtering is exhaustive; might not remove some redundant excitations?
-        all_excitations = [
-            _ex
-            for _ex in all_excitations
-            if sum(_ex) % 2 == 0 and (sum(_i % 2 for _i in _ex[:order]) == sum(_i % 2 for _i in _ex[order:]))
-        ]
-    return all_excitations
-
-
-def sort_excitations(excitations):
-    """
-    Sorts a list of excitations according to some common-sense and empirical rules (see below).
-        The order of the excitations must be the same throughout.
-
-    Sorting order:
-    1. (For double excitations only) All paired excitations between the same MOs first
-    2. Pair up excitations between the same MOs, e.g. (0, 2) and (1, 3)
-    3. Then count upwards from smallest MO
-
-    Args:
-        excitations: List of iterables, each representing an excitation; e.g. [[1, 5], [0, 4]]
-
-    Returns:
-        List of excitations after sorting
-    """
-    # Check that all excitations are of the same order and <= 2
-    order = len(excitations[0]) // 2
-    if order > 2:
-        raise NotImplementedError("Can only handle single and double excitations!")
-    assert all(len(_ex) // 2 == order for _ex in excitations), "Cannot handle excitations of different orders!"
-
-    # Define variables for the while loop
-    copy_excitations = [list(_ex) for _ex in excitations]
-    result = []
-    prev = []
-
-    # No idea how I came up with this, but it seems to work for double excitations
-    def sorting_fn(x):
-        # Default sorting is OK for single excitations
-        return sum((order + 1 - _i) * abs(x[2 * _i + 1] // 2 - x[2 * _i] // 2) for _i in range(0, order))
-
-    # Make a copy of the list of excitations, and use it populate a new list iteratively
-    while copy_excitations:
-        if not prev:
-            # Take out all pair excitations first
-            pair_excitations = [
-                _ex
-                for _ex in copy_excitations
-                # Indices of the electrons/holes must be consecutive numbers
-                if sum(abs(_ex[2 * _i + 1] // 2 - _ex[2 * _i] // 2) for _i in range(0, order)) == 0
-            ]
-            while pair_excitations:
-                pair_excitations = sorted(pair_excitations)
-                ex_to_remove = pair_excitations.pop(0)
-                if ex_to_remove in copy_excitations:
-                    # 'Move' the first pair excitation from copy_excitations to result
-                    index = copy_excitations.index(ex_to_remove)
-                    result.append(copy_excitations.pop(index))
-
-            # No more pair excitations, only remaining excitations should have >=3 MOs involved
-            # Sort the remaining excitations
-            copy_excitations = sorted(copy_excitations, key=sorting_fn if order != 1 else None)
-        else:
-            # Check to see for excitations involving the same MOs as prev
-            _from = prev[:order]
-            _to = prev[order:]
-            # Get all possible excitations involving the same MOs
-            new_from = [_i + 1 if _i % 2 == 0 else _i - 1 for _i in _from]
-            new_to = [_i + 1 if _i % 2 == 0 else _i - 1 for _i in _to]
-            same_mo_ex = [sorted(list(_f) + list(_t)) for _f in (_from, new_from) for _t in (_to, new_to)]
-            # Remove the excitations with the same MOs from copy_excitations
-            while same_mo_ex:
-                same_mo_ex = sorted(same_mo_ex)
-                ex_to_remove = same_mo_ex.pop(0)
-                if ex_to_remove in copy_excitations:
-                    # 'Move' the first entry of same_mo_index from copy_excitations to result
-                    index = copy_excitations.index(ex_to_remove)
-                    result.append(copy_excitations.pop(index))
-            prev = None
-            continue
-        if copy_excitations:
-            # Remove the first entry from the sorted list of remaining excitations and add it to result
-            prev = copy_excitations.pop(0)
-            result.append(prev)
-    return result
-
-
 def ucc_ansatz(
     molecule,
     excitation_level=None,
diff --git a/src/qibochem/ansatz/util.py b/src/qibochem/ansatz/util.py
new file mode 100644
index 00000000..a3521af8
--- /dev/null
+++ b/src/qibochem/ansatz/util.py
@@ -0,0 +1,154 @@
+"""
+Utility functions that can be used by different ansatzes
+"""
+
+
+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): :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,
+            representing a single or double excitation respectively.
+        orbital_energies: eigenvalues of the Fock operator, i.e. orbital energies
+        tei: Two-electron integrals in MO basis and second quantization notation
+
+    Returns:
+        MP2 guess amplitude (float)
+    """
+    # Checks validity of excitation argument
+    assert len(excitation) % 2 == 0 and len(excitation) // 2 <= 2, f"{excitation} must have either 2 or 4 elements"
+    # If single excitation, can just return 0.0 directly
+    if len(excitation) == 2:
+        return 0.0
+
+    # Convert orbital indices to be in MO basis
+    mo_orbitals = [orbital // 2 for orbital in excitation]
+    # Numerator: g_ijab - g_ijba
+    g_ijab = (
+        tei[tuple(mo_orbitals)]  # Can index directly using the MO TEIs
+        if (excitation[0] + excitation[3]) % 2 == 0 and (excitation[1] + excitation[2]) % 2 == 0
+        else 0.0
+    )
+    g_ijba = (
+        tei[tuple(mo_orbitals[:2] + mo_orbitals[2:][::-1])]  # Reverse last two terms
+        if (excitation[0] + excitation[2]) % 2 == 0 and (excitation[1] + excitation[3]) % 2 == 0
+        else 0.0
+    )
+    numerator = g_ijab - g_ijba
+    # Denominator is directly from the orbital energies
+    denominator = sum(orbital_energies[mo_orbitals[:2]]) - sum(orbital_energies[mo_orbitals[2:]])
+    return numerator / denominator
+
+
+def generate_excitations(order, excite_from, excite_to, conserve_spin=True):
+    """
+    Generate all possible excitations between a list of occupied and virtual orbitals
+
+    Args:
+        order: Order of excitations, i.e. 1 == single, 2 == double
+        excite_from: Iterable of integers
+        excite_to: Iterable of integers
+        conserve_spin: ensure that the total electronic spin is conserved
+
+    Return:
+        List of lists, e.g. [[0, 1]]
+    """
+    # If order of excitation > either number of electrons/orbitals, return list of empty list
+    if order > min(len(excite_from), len(excite_to)):
+        return [[]]
+
+    # Generate all possible excitations first
+    from itertools import combinations  # pylint: disable=C0415
+
+    all_excitations = [
+        [*_from, *_to] for _from in combinations(excite_from, order) for _to in combinations(excite_to, order)
+    ]
+    # Filter out the excitations if conserve_spin set
+    if conserve_spin:
+        # Not sure if this filtering is exhaustive; might not remove some redundant excitations?
+        all_excitations = [
+            _ex
+            for _ex in all_excitations
+            if sum(_ex) % 2 == 0 and (sum(_i % 2 for _i in _ex[:order]) == sum(_i % 2 for _i in _ex[order:]))
+        ]
+    return all_excitations
+
+
+def sort_excitations(excitations):
+    """
+    Sorts a list of excitations according to some common-sense and empirical rules (see below).
+        The order of the excitations must be the same throughout.
+
+    Sorting order:
+    1. (For double excitations only) All paired excitations between the same MOs first
+    2. Pair up excitations between the same MOs, e.g. (0, 2) and (1, 3)
+    3. Then count upwards from smallest MO
+
+    Args:
+        excitations: List of iterables, each representing an excitation; e.g. [[1, 5], [0, 4]]
+
+    Returns:
+        List of excitations after sorting
+    """
+    # Check that all excitations are of the same order and <= 2
+    order = len(excitations[0]) // 2
+    if order > 2:
+        raise NotImplementedError("Can only handle single and double excitations!")
+    assert all(len(_ex) // 2 == order for _ex in excitations), "Cannot handle excitations of different orders!"
+
+    # Define variables for the while loop
+    copy_excitations = [list(_ex) for _ex in excitations]
+    result = []
+    prev = []
+
+    # No idea how I came up with this, but it seems to work for double excitations
+    def sorting_fn(x):
+        # Default sorting is OK for single excitations
+        return sum((order + 1 - _i) * abs(x[2 * _i + 1] // 2 - x[2 * _i] // 2) for _i in range(0, order))
+
+    # Make a copy of the list of excitations, and use it populate a new list iteratively
+    while copy_excitations:
+        if not prev:
+            # Take out all pair excitations first
+            pair_excitations = [
+                _ex
+                for _ex in copy_excitations
+                # Indices of the electrons/holes must be consecutive numbers
+                if sum(abs(_ex[2 * _i + 1] // 2 - _ex[2 * _i] // 2) for _i in range(0, order)) == 0
+            ]
+            while pair_excitations:
+                pair_excitations = sorted(pair_excitations)
+                ex_to_remove = pair_excitations.pop(0)
+                if ex_to_remove in copy_excitations:
+                    # 'Move' the first pair excitation from copy_excitations to result
+                    index = copy_excitations.index(ex_to_remove)
+                    result.append(copy_excitations.pop(index))
+
+            # No more pair excitations, only remaining excitations should have >=3 MOs involved
+            # Sort the remaining excitations
+            copy_excitations = sorted(copy_excitations, key=sorting_fn if order != 1 else None)
+        else:
+            # Check to see for excitations involving the same MOs as prev
+            _from = prev[:order]
+            _to = prev[order:]
+            # Get all possible excitations involving the same MOs
+            new_from = [_i + 1 if _i % 2 == 0 else _i - 1 for _i in _from]
+            new_to = [_i + 1 if _i % 2 == 0 else _i - 1 for _i in _to]
+            same_mo_ex = [sorted(list(_f) + list(_t)) for _f in (_from, new_from) for _t in (_to, new_to)]
+            # Remove the excitations with the same MOs from copy_excitations
+            while same_mo_ex:
+                same_mo_ex = sorted(same_mo_ex)
+                ex_to_remove = same_mo_ex.pop(0)
+                if ex_to_remove in copy_excitations:
+                    # 'Move' the first entry of same_mo_index from copy_excitations to result
+                    index = copy_excitations.index(ex_to_remove)
+                    result.append(copy_excitations.pop(index))
+            prev = None
+            continue
+        if copy_excitations:
+            # Remove the first entry from the sorted list of remaining excitations and add it to result
+            prev = copy_excitations.pop(0)
+            result.append(prev)
+    return result
diff --git a/tests/test_ansatz_util.py b/tests/test_ansatz_util.py
new file mode 100644
index 00000000..9a81548b
--- /dev/null
+++ b/tests/test_ansatz_util.py
@@ -0,0 +1,53 @@
+"""
+Tests for ansatz utility functions (util.py)
+"""
+
+import numpy as np
+import pytest
+
+from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations
+from qibochem.driver import Molecule
+
+
+@pytest.mark.parametrize(
+    "order,excite_from,excite_to,expected",
+    [
+        (1, [0, 1], [2, 3], [[0, 2], [1, 3]]),
+        (2, [2, 3], [4, 5], [[2, 3, 4, 5]]),
+        (3, [0, 1], [2, 3], [[]]),
+    ],
+)
+def test_generate_excitations(order, excite_from, excite_to, expected):
+    """Test generation of all possible excitations between two lists of orbitals"""
+    test = generate_excitations(order, excite_from, excite_to)
+    assert test == expected
+
+
+@pytest.mark.parametrize(
+    "order,excite_from,excite_to,expected",
+    [
+        (1, [0, 1], [2, 3, 4, 5], [[0, 2], [1, 3], [0, 4], [1, 5]]),
+        (2, [0, 1], [2, 3, 4, 5], [[0, 1, 2, 3], [0, 1, 4, 5], [0, 1, 2, 5], [0, 1, 3, 4]]),
+    ],
+)
+def test_sort_excitations(order, excite_from, excite_to, expected):
+    test = sort_excitations(generate_excitations(order, excite_from, excite_to))
+    assert test == expected
+
+
+def test_sort_excitations_triples():
+    with pytest.raises(NotImplementedError):
+        sort_excitations([[1, 2, 3, 4, 5, 6]])
+
+
+def test_mp2_amplitude_singles():
+    assert mp2_amplitude([0, 2], np.random.rand(4), np.random.rand(4, 4)) == 0.0
+
+
+def test_mp2_amplitude_doubles():
+    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    h2.run_pyscf()
+    l = mp2_amplitude([0, 1, 2, 3], h2.eps, h2.tei)
+    ref_l = 0.06834019757197053
+
+    assert np.isclose(l, ref_l)
diff --git a/tests/test_givens_excitation.py b/tests/test_givens_excitation.py
new file mode 100644
index 00000000..2344a202
--- /dev/null
+++ b/tests/test_givens_excitation.py
@@ -0,0 +1,189 @@
+"""
+Tests for the symmetry preserving ansatz from Gard et al. (DOI: https://doi.org/10.1038/s41534-019-0240-1)
+"""
+
+import numpy as np
+import pytest
+from qibo import Circuit, gates
+
+from qibochem.ansatz import hf_circuit
+from qibochem.ansatz.givens_excitation import (
+    double_excitation_gate,
+    givens_excitation_ansatz,
+    givens_excitation_circuit,
+    single_excitation_gate,
+)
+from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations
+from qibochem.driver import Molecule
+
+
+def test_single_excitation_gate():
+    # Hardcoded test
+    theta = 0.1
+
+    control_gates = [
+        gates.CNOT(0, 1),
+        gates.RY(0, 0.5 * theta),
+        gates.CNOT(1, 0),
+        gates.RY(0, -0.5 * theta),
+        gates.CNOT(1, 0),
+        gates.CNOT(0, 1),
+    ]
+    test_list = single_excitation_gate([0, 1, 2, 3], 0.1)
+
+    # Check gates are correct
+    assert all(
+        control.name == test.name and control.target_qubits == test.target_qubits
+        for control, test in zip(control_gates, test_list)
+    )
+
+
+def test_double_excitation_gate():
+    # Hardcoded test
+    theta = 0.0
+
+    control_gates = [
+        gates.CNOT(2, 3),
+        gates.CNOT(0, 2),
+        gates.H(0),
+        gates.H(3),
+        gates.CNOT(0, 1),
+        gates.CNOT(2, 3),
+        gates.RY(0, -theta),
+        gates.RY(1, theta),
+        gates.CNOT(0, 3),
+        gates.H(3),
+        gates.CNOT(3, 1),
+        gates.RY(0, -theta),
+        gates.RY(1, theta),
+        gates.CNOT(2, 1),
+        gates.CNOT(2, 0),
+        gates.RY(0, theta),
+        gates.RY(1, -theta),
+        gates.CNOT(3, 1),
+        gates.H(3),
+        gates.CNOT(0, 3),
+        gates.RY(0, theta),
+        gates.RY(1, -theta),
+        gates.CNOT(0, 1),
+        gates.CNOT(2, 0),
+        gates.H(0),
+        gates.H(3),
+        gates.CNOT(0, 2),
+        gates.CNOT(2, 3),
+    ]
+    test_list = double_excitation_gate([0, 1, 2, 3], 0.0)
+
+    # Check gates are correct
+    assert all(
+        control.name == test.name and control.target_qubits == test.target_qubits
+        for control, test in zip(control_gates, test_list)
+    )
+
+
+@pytest.mark.parametrize(
+    "excitation,expected",
+    [
+        ([0, 2], single_excitation_gate([0, 2], 0.0)),
+        ([0, 1, 2, 3], double_excitation_gate([0, 1, 2, 3], 0.0)),
+    ],
+)
+def test_givens_excitation_circuit(excitation, expected):
+    test_circuit = givens_excitation_circuit(4, excitation)
+    assert all(
+        control.name == test.name and control.target_qubits == test.target_qubits
+        for control, test in zip(expected, list(test_circuit.queue))
+    )
+    # Check parameters of parametrized gates are all 0.0
+    assert all(np.isclose(gate.parameters[0], 0.0) for gate in test_circuit.queue if gate.parameters)
+
+
+def test_givens_excitation_errors():
+    """Input excitations are single or double?"""
+    with pytest.raises(NotImplementedError):
+        test_circuit = givens_excitation_circuit(4, list(range(6)))
+
+
+def test_givens_excitation_ansatz_h2():
+    """Test the default arguments of ucc_ansatz using H2"""
+    mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    mol.run_pyscf()
+
+    # Build control circuit
+    control_circuit = hf_circuit(4, 2)
+    excitations = ([0, 1, 2, 3], [0, 2], [1, 3])
+    for excitation in excitations:
+        control_circuit += givens_excitation_circuit(4, excitation)
+
+    test_circuit = givens_excitation_ansatz(mol)
+
+    assert all(
+        control.name == test.name and control.target_qubits == test.target_qubits
+        # for control, test in zip(list(control_circuit.queue), list(test_circuit.queue))
+        for control, test in zip(control_circuit.queue, test_circuit.queue)
+    )
+    # Check that number of parametrised gates is the same
+    assert len(control_circuit.get_parameters()) == len(test_circuit.get_parameters())
+
+    # Then check that the circuit parameters are the MP2 guess parameters
+    # Get the MP2 amplitudes first, then expand the list based on the excitation type
+    mp2_guess_amplitudes = [mp2_amplitude([0, 1, 2, 3], mol.eps, mol.tei) for _ in range(8)]  # Doubles
+    mp2_guess_amplitudes += [0.0, 0.0, 0.0, 0.0]  # Singles
+    coeffs = np.array([-0.125, 0.125, -0.125, 0.125, 0.125, -0.125, 0.125, -0.125, 1.0, 1.0, 1.0, 1.0])
+    mp2_guess_amplitudes = coeffs * np.array(mp2_guess_amplitudes)
+    # Need to flatten the output of circuit.get_parameters() to compare it to mp2_guess_amplitudes
+    test_parameters = np.array([_x for _tuple in test_circuit.get_parameters() for _x in _tuple])
+    assert np.allclose(mp2_guess_amplitudes, test_parameters)
+
+
+def test_givens_excitation_ansatz_embedding():
+    """Test the default arguments of ucc_ansatz using LiH with HF embedding applied, but without the HF circuit"""
+    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))])
+    mol.run_pyscf()
+    mol.hf_embedding(active=[1, 2, 5])
+
+    # Generate all possible excitations
+    excitations = []
+    for order in range(2, 0, -1):
+        # 2 electrons, 6 spin-orbitals
+        excitations += sort_excitations(generate_excitations(order, range(0, 2), range(2, 6)))
+    # Build control circuit
+    control_circuit = hf_circuit(6, 0)
+    for excitation in excitations:
+        control_circuit += givens_excitation_circuit(6, excitation)
+
+    test_circuit = givens_excitation_ansatz(mol, include_hf=False, use_mp2_guess=False)
+
+    assert all(
+        control.name == test.name and control.target_qubits == test.target_qubits
+        for control, test in zip(list(control_circuit.queue), list(test_circuit.queue))
+    )
+    # Check that number of parametrised gates is the same
+    assert len(control_circuit.get_parameters()) == len(test_circuit.get_parameters())
+
+    # Check that the circuit parameters are all zeros
+    test_parameters = np.array([_x for _tuple in test_circuit.get_parameters() for _x in _tuple])
+    assert np.allclose(test_parameters, np.zeros(len(test_parameters)))
+
+
+def test_ucc_ansatz_excitations():
+    """Test the `excitations` argument of ucc_ansatz"""
+    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))])
+    mol.run_pyscf()
+    mol.hf_embedding(active=[1, 2, 5])
+
+    # Generate all possible excitations
+    excitations = [[0, 1, 2, 3], [0, 1, 4, 5]]
+    # Build control circuit
+    control_circuit = hf_circuit(6, 2)
+    for excitation in excitations:
+        control_circuit += givens_excitation_circuit(6, excitation)
+
+    test_circuit = givens_excitation_ansatz(mol, excitations=excitations)
+
+    assert all(
+        control.name == test.name and control.target_qubits == test.target_qubits
+        for control, test in zip(list(control_circuit.queue), list(test_circuit.queue))
+    )
+    # Check that number of parametrised gates is the same
+    assert len(control_circuit.get_parameters()) == len(test_circuit.get_parameters())
diff --git a/tests/test_symmetry_preserving_ansatz.py b/tests/test_symmetry_preserving_ansatz.py
new file mode 100644
index 00000000..72a7fa73
--- /dev/null
+++ b/tests/test_symmetry_preserving_ansatz.py
@@ -0,0 +1,81 @@
+"""
+Tests for the symmetry preserving ansatz from Gard et al. (DOI: https://doi.org/10.1038/s41534-019-0240-1)
+"""
+
+import numpy as np
+import pytest
+from qibo import Circuit, gates
+
+from qibochem.ansatz.symmetry import (
+    a_gate,
+    a_gate_indices,
+    symm_preserving_circuit,
+    x_gate_indices,
+)
+
+
+@pytest.mark.parametrize(
+    "theta,phi,expected",
+    [
+        (None, None, np.array([0.0, 0.0, -1.0, 00])),
+        (0.5 * np.pi, 0.0, np.array([0.0, 1.0, 0.0, 00])),
+        (0.5 * np.pi, np.pi, np.array([0.0, -1.0, 0.0, 00])),
+    ],
+)
+def test_a_gate(theta, phi, expected):
+    circuit = Circuit(2)
+    circuit.add(gates.X(0))
+    a_gates = a_gate(0, 1, theta=theta, phi=phi)
+    circuit.add(a_gates)
+
+    result = circuit(nshots=1)
+    state_ket = result.state()
+    assert np.allclose(state_ket, expected)
+
+
+@pytest.mark.parametrize(
+    "n_qubits,n_electrons,expected",
+    [
+        (4, 2, [0, 2]),
+        (4, 3, [0, 1, 2]),
+        (4, 4, [0, 1, 2, 3]),
+    ],
+)
+def test_x_gate_indices(n_qubits, n_electrons, expected):
+    test = x_gate_indices(n_qubits, n_electrons)
+    assert test == expected
+
+
+@pytest.mark.parametrize(
+    "n_qubits,n_electrons,x_gates,expected",
+    [
+        (4, 2, [0, 2], 2 * [(0, 1), (2, 3), (1, 2)]),
+        (6, 4, [0, 1, 2, 4], 3 * [(2, 3), (4, 5), (3, 4), (1, 2), (0, 1)]),
+    ],
+)
+def test_a_gate_indices(n_qubits, n_electrons, x_gates, expected):
+    test = a_gate_indices(n_qubits, n_electrons, x_gates)
+    assert test == expected
+
+
+@pytest.mark.parametrize(
+    "n_qubits,n_electrons",
+    [
+        (4, 2),
+        (6, 4),
+    ],
+)
+def test_symm_preserving_circuit(n_qubits, n_electrons):
+    control_circuit = Circuit(n_qubits)
+    x_gates = x_gate_indices(n_qubits, n_electrons)
+    control_circuit.add(gates.X(_i) for _i in x_gates)
+    a_gate_qubits = a_gate_indices(n_qubits, n_electrons, x_gates)
+    a_gates = [a_gate(qubit1, qubit2) for qubit1, qubit2 in a_gate_qubits]
+    control_circuit.add(_gates for _a_gate in a_gates for _gates in _a_gate)
+
+    test_circuit = symm_preserving_circuit(n_qubits, n_electrons)
+
+    assert all(
+        control.name == test.name and control.target_qubits == test.target_qubits
+        for control, test in zip(list(control_circuit.queue), list(test_circuit.queue))
+    )
diff --git a/tests/test_ucc.py b/tests/test_ucc.py
index 90ec930a..8695a2b8 100644
--- a/tests/test_ucc.py
+++ b/tests/test_ucc.py
@@ -1,5 +1,5 @@
 """
-Tests for the UCC ansatz and related functions
+Tests for the UCC ansatz
 """
 
 from functools import reduce
@@ -11,61 +11,11 @@
 
 from qibochem.ansatz import hf_circuit
 from qibochem.ansatz.qeb import qeb_circuit
-from qibochem.ansatz.ucc import (
-    expi_pauli,
-    generate_excitations,
-    mp2_amplitude,
-    sort_excitations,
-    ucc_ansatz,
-    ucc_circuit,
-)
+from qibochem.ansatz.ucc import expi_pauli, ucc_ansatz, ucc_circuit
+from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations
 from qibochem.driver import Molecule
 
 
-@pytest.mark.parametrize(
-    "order,excite_from,excite_to,expected",
-    [
-        (1, [0, 1], [2, 3], [[0, 2], [1, 3]]),
-        (2, [2, 3], [4, 5], [[2, 3, 4, 5]]),
-        (3, [0, 1], [2, 3], [[]]),
-    ],
-)
-def test_generate_excitations(order, excite_from, excite_to, expected):
-    """Test generation of all possible excitations between two lists of orbitals"""
-    test = generate_excitations(order, excite_from, excite_to)
-    assert test == expected
-
-
-@pytest.mark.parametrize(
-    "order,excite_from,excite_to,expected",
-    [
-        (1, [0, 1], [2, 3, 4, 5], [[0, 2], [1, 3], [0, 4], [1, 5]]),
-        (2, [0, 1], [2, 3, 4, 5], [[0, 1, 2, 3], [0, 1, 4, 5], [0, 1, 2, 5], [0, 1, 3, 4]]),
-    ],
-)
-def test_sort_excitations(order, excite_from, excite_to, expected):
-    test = sort_excitations(generate_excitations(order, excite_from, excite_to))
-    assert test == expected
-
-
-def test_sort_excitations_triples():
-    with pytest.raises(NotImplementedError):
-        sort_excitations([[1, 2, 3, 4, 5, 6]])
-
-
-def test_mp2_amplitude_singles():
-    assert mp2_amplitude([0, 2], np.random.rand(4), np.random.rand(4, 4)) == 0.0
-
-
-def test_mp2_amplitude_doubles():
-    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
-    h2.run_pyscf()
-    l = mp2_amplitude([0, 1, 2, 3], h2.eps, h2.tei)
-    ref_l = 0.06834019757197053
-
-    assert np.isclose(l, ref_l)
-
-
 @pytest.mark.parametrize(
     "pauli_string",
     [