From 69ccb0ddd4e041ee78994f90c98a2d4ec610e841 Mon Sep 17 00:00:00 2001
From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com>
Date: Wed, 7 Feb 2024 08:09:04 +0000
Subject: [PATCH 01/11] Minor changes to default values and documentation

---
 src/qibochem/driver/molecule.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/qibochem/driver/molecule.py b/src/qibochem/driver/molecule.py
index a5995b4..f17930a 100644
--- a/src/qibochem/driver/molecule.py
+++ b/src/qibochem/driver/molecule.py
@@ -72,7 +72,7 @@ def __init__(self, geometry=None, charge=0, multiplicity=1, basis=None, xyz_file
 
         # For HF embedding
         self.active = active  #: Iterable of molecular orbitals included in the active space
-        self.frozen = None
+        self.frozen = None  #: Iterable representing the occupied molecular orbitals removed from the simulation
 
         self.inactive_energy = None
         self.embed_oei = None
@@ -341,11 +341,11 @@ def hamiltonian(
         if ham_type is None:
             ham_type = "sym"
         if oei is None:
-            oei = self.oei
+            oei = self.oei if self.embed_oei is None else self.embed_oei
         if tei is None:
-            tei = self.tei
+            tei = self.tei if self.embed_tei is None else self.embed_tei
         if constant is None:
-            constant = 0.0
+            constant = 0.0 if self.inactive_energy is None else self.inactive_energy
         if ferm_qubit_map is None:
             ferm_qubit_map = "jw"
 

From 97965512b5f3119d677634d0a8e33d4afc89e6b2 Mon Sep 17 00:00:00 2001
From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com>
Date: Wed, 7 Feb 2024 08:34:48 +0000
Subject: [PATCH 02/11] Streamline current HF embedding testing

---
 tests/test_molecule.py | 556 ++++++++++++++++++++---------------------
 1 file changed, 276 insertions(+), 280 deletions(-)

diff --git a/tests/test_molecule.py b/tests/test_molecule.py
index e3038b1..d587585 100644
--- a/tests/test_molecule.py
+++ b/tests/test_molecule.py
@@ -1,280 +1,276 @@
-"""
-Test Molecule class functions
-"""
-
-from pathlib import Path
-
-import numpy as np
-import openfermion
-import pytest
-from qibo import gates, models
-from qibo.hamiltonians import SymbolicHamiltonian
-from qibo.symbols import X, Y, Z
-
-from qibochem.driver import hamiltonian
-from qibochem.driver.molecule import Molecule
-from qibochem.measurement.expectation import expectation
-
-
-def test_run_pyscf():
-    """PySCF driver"""
-    # Hardcoded benchmark results
-    # Change to run PySCF directly during a test?
-    h2_ref_energy = -1.117349035
-    h2_ref_hcore = np.array([[-1.14765024, -1.00692423], [-1.00692423, -1.14765024]])
-
-    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
-    h2.run_pyscf()
-    assert h2.e_hf == pytest.approx(h2_ref_energy)
-    assert np.allclose(h2.hcore, h2_ref_hcore)
-
-
-def test_run_pyscf_molecule_xyz():
-    """Pyscf driver with xyz file"""
-    file_path = Path("tests", "data", "lih.xyz")
-    if not file_path.is_file():
-        with open(file_path, "w") as file_handler:
-            file_handler.write("2\n0 1\nLi 0.0 0.0 0.0\nH 0.0 0.0 1.2\n")
-    lih_ref_energy = -7.83561582555692
-    lih = Molecule(xyz_file=file_path)
-    lih.run_pyscf()
-
-    assert lih.e_hf == pytest.approx(lih_ref_energy)
-
-
-def test_run_pyscf_molecule_xyz_charged():
-    file_path = Path("tests", "data", "h2.xyz")
-    if not file_path.is_file():
-        with open(file_path, "w") as file_handler:
-            file_handler.write("2\n \nH 0.0 0.0 0.0\nH 0.0 0.0 0.7\n")
-    h2_ref_energy = -1.117349035
-    h2 = Molecule(xyz_file=file_path)
-    h2.run_pyscf()
-
-    assert h2.e_hf == pytest.approx(h2_ref_energy)
-
-
-def test_molecule_custom_basis():
-    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.2))], 0, 1, "6-31g")
-    mol.run_pyscf()
-    assert np.isclose(mol.e_hf, -7.94129296352493)
-
-
-def test_hf_embedding_1():
-    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.2))])
-    mol.run_pyscf()
-    ref_oei = mol.oei
-    ref_tei = mol.tei
-    mol.hf_embedding()
-    embed_oei = mol.embed_oei
-    embed_tei = mol.embed_tei
-    assert np.allclose(embed_oei, ref_oei)
-    assert np.allclose(embed_tei, ref_tei)
-
-
-def test_hf_embedding_2():
-    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.2))])
-    mol.run_pyscf()
-    mol.frozen = [0]
-    mol.active = [1, 2]
-    mol.hf_embedding()
-    assert mol.n_active_orbs == 4
-    assert mol.n_active_e == 2
-
-
-def test_fermionic_hamiltonian():
-    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
-    h2.run_pyscf()
-    h2_ferm_ham = hamiltonian.fermionic_hamiltonian(h2.oei, h2.tei, h2.e_nuc)
-
-    # a^\dagger_0 a_0
-    ref_one_body_tensor = np.array(
-        [
-            [-1.27785301, 0.0, 0.0, 0.0],
-            [0.0, -1.27785301, 0.0, 0.0],
-            [0.0, 0.0, -0.4482997, 0.0],
-            [0.0, 0.0, 0.0, -0.4482997],
-        ]
-    )
-
-    assert np.isclose(h2_ferm_ham[()], 0.7559674441714287)
-    assert np.allclose(h2_ferm_ham.one_body_tensor, ref_one_body_tensor)
-
-
-def test_fermionic_hamiltonian_2():
-    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7414))])
-    h2.run_pyscf()
-
-    h2_ferm_ham_1 = h2.hamiltonian("f", h2.oei, h2.tei)
-    h2_qub_ham_jw_1 = h2.hamiltonian("q", h2.oei, h2.tei, ferm_qubit_map="jw")
-    h2_qub_ham_bk_1 = h2.hamiltonian("q", h2.oei, h2.tei, ferm_qubit_map="bk")
-    # h2_mol_ham has format of InteractionOperator
-    h2_mol_ham = hamiltonian.fermionic_hamiltonian(h2.oei, h2.tei, h2.e_nuc)
-    h2_ferm_ham_2 = openfermion.transforms.get_fermion_operator(h2_mol_ham)
-    h2_qub_ham_jw_2 = openfermion.jordan_wigner(h2_mol_ham)
-    h2_qub_ham_bk_2 = openfermion.bravyi_kitaev(h2_mol_ham)
-
-    assert h2_ferm_ham_2.isclose(h2_ferm_ham_1)
-    assert h2_qub_ham_jw_2.isclose(h2_qub_ham_jw_1)
-    assert h2_qub_ham_bk_2.isclose(h2_qub_ham_bk_1)
-
-
-def test_parse_pauli_string_1():
-    pauli_string = ((0, "X"), (1, "Y"))
-    qibo_pauli_string = hamiltonian.parse_pauli_string(pauli_string, 0.5)
-    ref_pauli_string = "0.5*X0*Y1"
-    assert str(qibo_pauli_string) == ref_pauli_string
-
-
-def test_parse_pauli_string_2():
-    qibo_pauli_string = hamiltonian.parse_pauli_string(None, 0.1)
-    assert str(qibo_pauli_string) == "0.1"
-
-
-def test_qubit_hamiltonian():
-    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
-    h2.run_pyscf()
-    h2_ferm_ham = hamiltonian.fermionic_hamiltonian(h2.oei, h2.tei, h2.e_nuc)
-    h2_qubit_ham_jw = hamiltonian.qubit_hamiltonian(h2_ferm_ham, "jw")
-    h2_qubit_ham_bk = hamiltonian.qubit_hamiltonian(h2_ferm_ham, "bk")
-    ref_h2_qubit_ham_jw = {
-        (): -0.04207897647782238,
-        ((0, "Z"),): 0.17771287465139918,
-        ((1, "Z"),): 0.1777128746513992,
-        ((2, "Z"),): -0.24274280513140478,
-        ((3, "Z"),): -0.24274280513140478,
-        ((0, "Z"), (1, "Z")): 0.17059738328801052,
-        ((0, "Z"), (2, "Z")): 0.12293305056183809,
-        ((0, "Z"), (3, "Z")): 0.16768319457718972,
-        ((1, "Z"), (2, "Z")): 0.16768319457718972,
-        ((1, "Z"), (3, "Z")): 0.12293305056183809,
-        ((2, "Z"), (3, "Z")): 0.17627640804319608,
-        ((0, "X"), (1, "X"), (2, "Y"), (3, "Y")): -0.04475014401535165,
-        ((0, "X"), (1, "Y"), (2, "Y"), (3, "X")): 0.04475014401535165,
-        ((0, "Y"), (1, "X"), (2, "X"), (3, "Y")): 0.04475014401535165,
-        ((0, "Y"), (1, "Y"), (2, "X"), (3, "X")): -0.04475014401535165,
-    }
-    ref_h2_qubit_ham_bk = {
-        ((0, "Z"),): 0.17771287465139923,
-        (): -0.04207897647782244,
-        ((0, "Z"), (1, "Z")): 0.17771287465139918,
-        ((2, "Z"),): -0.24274280513140484,
-        ((1, "Z"), (2, "Z"), (3, "Z")): -0.24274280513140484,
-        ((0, "Y"), (1, "Z"), (2, "Y")): 0.04475014401535165,
-        ((0, "X"), (1, "Z"), (2, "X")): 0.04475014401535165,
-        ((0, "X"), (1, "Z"), (2, "X"), (3, "Z")): 0.04475014401535165,
-        ((0, "Y"), (1, "Z"), (2, "Y"), (3, "Z")): 0.04475014401535165,
-        ((1, "Z"),): 0.17059738328801052,
-        ((0, "Z"), (2, "Z")): 0.12293305056183809,
-        ((0, "Z"), (1, "Z"), (2, "Z")): 0.16768319457718972,
-        ((0, "Z"), (1, "Z"), (2, "Z"), (3, "Z")): 0.16768319457718972,
-        ((0, "Z"), (2, "Z"), (3, "Z")): 0.12293305056183809,
-        ((1, "Z"), (3, "Z")): 0.17627640804319608,
-    }
-
-    jw_array = np.array([terms for terms in h2_qubit_ham_jw.terms.values()])
-    bk_array = np.array([terms for terms in h2_qubit_ham_bk.terms.values()])
-
-    ref_jw_array = np.array([terms for terms in ref_h2_qubit_ham_jw.values()])
-    ref_bk_array = np.array([terms for terms in ref_h2_qubit_ham_bk.values()])
-
-    assert np.allclose(jw_array, ref_jw_array)
-    assert np.allclose(bk_array, ref_bk_array)
-
-    # incorrect mapping circuit
-    with pytest.raises(KeyError):
-        hamiltonian.qubit_hamiltonian(h2_ferm_ham, "incorrect")
-
-
-def test_symbolic_hamiltonian():
-    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
-    h2.run_pyscf()
-    h2_ferm_ham = hamiltonian.fermionic_hamiltonian(h2.oei, h2.tei, h2.e_nuc)
-    h2_qubit_ham = hamiltonian.qubit_hamiltonian(h2_ferm_ham, "jw")
-    h2_sym_ham = hamiltonian.symbolic_hamiltonian(h2_qubit_ham)
-    ref_ham = (
-        -0.0420789764778224
-        - 0.0447501440153516 * X(0) * X(1) * Y(2) * Y(3)
-        + 0.0447501440153516 * X(0) * Y(1) * Y(2) * X(3)
-        + 0.0447501440153516 * Y(0) * X(1) * X(2) * Y(3)
-        - 0.0447501440153516 * Y(0) * Y(1) * X(2) * X(3)
-        + 0.177712874651399 * Z(0)
-        + 0.170597383288011 * Z(0) * Z(1)
-        + 0.122933050561838 * Z(0) * Z(2)
-        + 0.16768319457719 * Z(0) * Z(3)
-        + 0.177712874651399 * Z(1)
-        + 0.16768319457719 * Z(1) * Z(2)
-        + 0.122933050561838 * Z(1) * Z(3)
-        - 0.242742805131405 * Z(2)
-        + 0.176276408043196 * Z(2) * Z(3)
-        - 0.242742805131405 * Z(3)
-    )
-    ref_sym_ham = SymbolicHamiltonian(ref_ham)
-
-    assert np.allclose(h2_sym_ham.matrix, ref_sym_ham.matrix)
-
-
-def test_hamiltonian_input_error():
-    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
-    h2.e_nuc = 0.0
-    h2.oei = np.random.rand(4, 4)
-    h2.tei = np.random.rand(4, 4, 4, 4)
-    with pytest.raises(NameError):
-        h2.hamiltonian("ihpc")
-
-
-# Commenting out since not actively supporting PSI4 at the moment
-# @pytest.mark.skip(reason="Psi4 doesn't offer pip install, so needs to be installed through conda or manually.")
-# def test_run_psi4():
-#     """PSI4 driver"""
-#     # Hardcoded benchmark results
-#     h2_ref_energy = -1.117349035
-#     h2_ref_hcore = np.array([[-1.14765024, -1.00692423], [-1.00692423, -1.14765024]])
-#
-#     h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
-#     h2.run_psi4()
-#
-#     assert h2.e_hf == pytest.approx(h2_ref_energy)
-#     assert np.allclose(h2.hcore, h2_ref_hcore)
-
-
-def test_expectation_value():
-    """Tests generation of molecular Hamiltonian and its expectation value using a JW-HF circuit"""
-    # Hardcoded benchmark results
-    h2_ref_energy = -1.117349035
-
-    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
-    h2.run_pyscf()
-
-    # JW-HF circuit
-    circuit = models.Circuit(h2.nso)
-    circuit.add(gates.X(_i) for _i in range(h2.nelec))
-    # Molecular Hamiltonian and the HF expectation value
-    hamiltonian = h2.hamiltonian()
-    hf_energy = expectation(circuit, hamiltonian)
-
-    # assert h2.e_hf == pytest.approx(hf_energy)
-    assert h2_ref_energy == pytest.approx(hf_energy)
-
-
-def test_eigenvalues():
-    dummy = Molecule()
-    # FermionOperator test:
-    ferm_ham = sum(
-        openfermion.FermionOperator(f"{_i}^ {_i}") + openfermion.FermionOperator(f"{_i} {_i}^") for _i in range(4)
-    )
-    ham_matrix = openfermion.get_sparse_operator(ferm_ham).toarray()
-    assert np.allclose(dummy.eigenvalues(ferm_ham), sorted(np.linalg.eigvals(ham_matrix))[:6])
-
-    # QubitOperator test:
-    qubit_ham = openfermion.QubitOperator("Z0 Z1")
-    ham_matrix = openfermion.get_sparse_operator(qubit_ham).toarray()
-    assert np.allclose(dummy.eigenvalues(qubit_ham), sorted(np.kron(np.array([1.0, -1.0]), np.array([1.0, -1.0])))[:2])
-
-    # SymbolicHamiltonian test:
-    sym_ham = SymbolicHamiltonian(Z(0) * Z(1))
-    assert np.allclose(dummy.eigenvalues(sym_ham), sorted(np.kron(np.array([1.0, -1.0]), np.array([1.0, -1.0]))))
-
-    # Unknown Hamiltonian type
-    with pytest.raises(TypeError):
-        dummy.eigenvalues(0.0)
+"""
+Test Molecule class functions
+"""
+
+from pathlib import Path
+
+import numpy as np
+import openfermion
+import pytest
+from qibo import gates, models
+from qibo.hamiltonians import SymbolicHamiltonian
+from qibo.symbols import X, Y, Z
+
+from qibochem.driver import hamiltonian
+from qibochem.driver.molecule import Molecule
+from qibochem.measurement.expectation import expectation
+
+
+def test_run_pyscf():
+    """PySCF driver"""
+    # Hardcoded benchmark results
+    # Change to run PySCF directly during a test?
+    h2_ref_energy = -1.117349035
+    h2_ref_hcore = np.array([[-1.14765024, -1.00692423], [-1.00692423, -1.14765024]])
+
+    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    h2.run_pyscf()
+    assert h2.e_hf == pytest.approx(h2_ref_energy)
+    assert np.allclose(h2.hcore, h2_ref_hcore)
+
+
+def test_run_pyscf_molecule_xyz():
+    """Pyscf driver with xyz file"""
+    file_path = Path("tests", "data", "lih.xyz")
+    if not file_path.is_file():
+        with open(file_path, "w") as file_handler:
+            file_handler.write("2\n0 1\nLi 0.0 0.0 0.0\nH 0.0 0.0 1.2\n")
+    lih_ref_energy = -7.83561582555692
+    lih = Molecule(xyz_file=file_path)
+    lih.run_pyscf()
+
+    assert lih.e_hf == pytest.approx(lih_ref_energy)
+
+
+def test_run_pyscf_molecule_xyz_charged():
+    file_path = Path("tests", "data", "h2.xyz")
+    if not file_path.is_file():
+        with open(file_path, "w") as file_handler:
+            file_handler.write("2\n \nH 0.0 0.0 0.0\nH 0.0 0.0 0.7\n")
+    h2_ref_energy = -1.117349035
+    h2 = Molecule(xyz_file=file_path)
+    h2.run_pyscf()
+
+    assert h2.e_hf == pytest.approx(h2_ref_energy)
+
+
+def test_molecule_custom_basis():
+    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.2))], 0, 1, "6-31g")
+    mol.run_pyscf()
+    assert np.isclose(mol.e_hf, -7.94129296352493)
+
+
+def test_hf_embedding():
+    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.2))])
+    mol.run_pyscf()
+    ref_oei = mol.oei
+    ref_tei = mol.tei
+    mol.hf_embedding()
+    embed_oei = mol.embed_oei
+    embed_tei = mol.embed_tei
+    assert np.allclose(embed_oei, ref_oei)
+    assert np.allclose(embed_tei, ref_tei)
+
+    mol.frozen = [0]
+    mol.active = [1, 2]
+    mol.hf_embedding()
+    assert mol.n_active_orbs == 4
+    assert mol.n_active_e == 2
+
+
+def test_fermionic_hamiltonian():
+    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    h2.run_pyscf()
+    h2_ferm_ham = hamiltonian.fermionic_hamiltonian(h2.oei, h2.tei, h2.e_nuc)
+
+    # a^\dagger_0 a_0
+    ref_one_body_tensor = np.array(
+        [
+            [-1.27785301, 0.0, 0.0, 0.0],
+            [0.0, -1.27785301, 0.0, 0.0],
+            [0.0, 0.0, -0.4482997, 0.0],
+            [0.0, 0.0, 0.0, -0.4482997],
+        ]
+    )
+
+    assert np.isclose(h2_ferm_ham[()], 0.7559674441714287)
+    assert np.allclose(h2_ferm_ham.one_body_tensor, ref_one_body_tensor)
+
+
+def test_fermionic_hamiltonian_2():
+    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7414))])
+    h2.run_pyscf()
+
+    h2_ferm_ham_1 = h2.hamiltonian("f", h2.oei, h2.tei)
+    h2_qub_ham_jw_1 = h2.hamiltonian("q", h2.oei, h2.tei, ferm_qubit_map="jw")
+    h2_qub_ham_bk_1 = h2.hamiltonian("q", h2.oei, h2.tei, ferm_qubit_map="bk")
+    # h2_mol_ham has format of InteractionOperator
+    h2_mol_ham = hamiltonian.fermionic_hamiltonian(h2.oei, h2.tei, h2.e_nuc)
+    h2_ferm_ham_2 = openfermion.transforms.get_fermion_operator(h2_mol_ham)
+    h2_qub_ham_jw_2 = openfermion.jordan_wigner(h2_mol_ham)
+    h2_qub_ham_bk_2 = openfermion.bravyi_kitaev(h2_mol_ham)
+
+    assert h2_ferm_ham_2.isclose(h2_ferm_ham_1)
+    assert h2_qub_ham_jw_2.isclose(h2_qub_ham_jw_1)
+    assert h2_qub_ham_bk_2.isclose(h2_qub_ham_bk_1)
+
+
+def test_parse_pauli_string_1():
+    pauli_string = ((0, "X"), (1, "Y"))
+    qibo_pauli_string = hamiltonian.parse_pauli_string(pauli_string, 0.5)
+    ref_pauli_string = "0.5*X0*Y1"
+    assert str(qibo_pauli_string) == ref_pauli_string
+
+
+def test_parse_pauli_string_2():
+    qibo_pauli_string = hamiltonian.parse_pauli_string(None, 0.1)
+    assert str(qibo_pauli_string) == "0.1"
+
+
+def test_qubit_hamiltonian():
+    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    h2.run_pyscf()
+    h2_ferm_ham = hamiltonian.fermionic_hamiltonian(h2.oei, h2.tei, h2.e_nuc)
+    h2_qubit_ham_jw = hamiltonian.qubit_hamiltonian(h2_ferm_ham, "jw")
+    h2_qubit_ham_bk = hamiltonian.qubit_hamiltonian(h2_ferm_ham, "bk")
+    ref_h2_qubit_ham_jw = {
+        (): -0.04207897647782238,
+        ((0, "Z"),): 0.17771287465139918,
+        ((1, "Z"),): 0.1777128746513992,
+        ((2, "Z"),): -0.24274280513140478,
+        ((3, "Z"),): -0.24274280513140478,
+        ((0, "Z"), (1, "Z")): 0.17059738328801052,
+        ((0, "Z"), (2, "Z")): 0.12293305056183809,
+        ((0, "Z"), (3, "Z")): 0.16768319457718972,
+        ((1, "Z"), (2, "Z")): 0.16768319457718972,
+        ((1, "Z"), (3, "Z")): 0.12293305056183809,
+        ((2, "Z"), (3, "Z")): 0.17627640804319608,
+        ((0, "X"), (1, "X"), (2, "Y"), (3, "Y")): -0.04475014401535165,
+        ((0, "X"), (1, "Y"), (2, "Y"), (3, "X")): 0.04475014401535165,
+        ((0, "Y"), (1, "X"), (2, "X"), (3, "Y")): 0.04475014401535165,
+        ((0, "Y"), (1, "Y"), (2, "X"), (3, "X")): -0.04475014401535165,
+    }
+    ref_h2_qubit_ham_bk = {
+        ((0, "Z"),): 0.17771287465139923,
+        (): -0.04207897647782244,
+        ((0, "Z"), (1, "Z")): 0.17771287465139918,
+        ((2, "Z"),): -0.24274280513140484,
+        ((1, "Z"), (2, "Z"), (3, "Z")): -0.24274280513140484,
+        ((0, "Y"), (1, "Z"), (2, "Y")): 0.04475014401535165,
+        ((0, "X"), (1, "Z"), (2, "X")): 0.04475014401535165,
+        ((0, "X"), (1, "Z"), (2, "X"), (3, "Z")): 0.04475014401535165,
+        ((0, "Y"), (1, "Z"), (2, "Y"), (3, "Z")): 0.04475014401535165,
+        ((1, "Z"),): 0.17059738328801052,
+        ((0, "Z"), (2, "Z")): 0.12293305056183809,
+        ((0, "Z"), (1, "Z"), (2, "Z")): 0.16768319457718972,
+        ((0, "Z"), (1, "Z"), (2, "Z"), (3, "Z")): 0.16768319457718972,
+        ((0, "Z"), (2, "Z"), (3, "Z")): 0.12293305056183809,
+        ((1, "Z"), (3, "Z")): 0.17627640804319608,
+    }
+
+    jw_array = np.array([terms for terms in h2_qubit_ham_jw.terms.values()])
+    bk_array = np.array([terms for terms in h2_qubit_ham_bk.terms.values()])
+
+    ref_jw_array = np.array([terms for terms in ref_h2_qubit_ham_jw.values()])
+    ref_bk_array = np.array([terms for terms in ref_h2_qubit_ham_bk.values()])
+
+    assert np.allclose(jw_array, ref_jw_array)
+    assert np.allclose(bk_array, ref_bk_array)
+
+    # incorrect mapping circuit
+    with pytest.raises(KeyError):
+        hamiltonian.qubit_hamiltonian(h2_ferm_ham, "incorrect")
+
+
+def test_symbolic_hamiltonian():
+    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    h2.run_pyscf()
+    h2_ferm_ham = hamiltonian.fermionic_hamiltonian(h2.oei, h2.tei, h2.e_nuc)
+    h2_qubit_ham = hamiltonian.qubit_hamiltonian(h2_ferm_ham, "jw")
+    h2_sym_ham = hamiltonian.symbolic_hamiltonian(h2_qubit_ham)
+    ref_ham = (
+        -0.0420789764778224
+        - 0.0447501440153516 * X(0) * X(1) * Y(2) * Y(3)
+        + 0.0447501440153516 * X(0) * Y(1) * Y(2) * X(3)
+        + 0.0447501440153516 * Y(0) * X(1) * X(2) * Y(3)
+        - 0.0447501440153516 * Y(0) * Y(1) * X(2) * X(3)
+        + 0.177712874651399 * Z(0)
+        + 0.170597383288011 * Z(0) * Z(1)
+        + 0.122933050561838 * Z(0) * Z(2)
+        + 0.16768319457719 * Z(0) * Z(3)
+        + 0.177712874651399 * Z(1)
+        + 0.16768319457719 * Z(1) * Z(2)
+        + 0.122933050561838 * Z(1) * Z(3)
+        - 0.242742805131405 * Z(2)
+        + 0.176276408043196 * Z(2) * Z(3)
+        - 0.242742805131405 * Z(3)
+    )
+    ref_sym_ham = SymbolicHamiltonian(ref_ham)
+
+    assert np.allclose(h2_sym_ham.matrix, ref_sym_ham.matrix)
+
+
+def test_hamiltonian_input_error():
+    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    h2.e_nuc = 0.0
+    h2.oei = np.random.rand(4, 4)
+    h2.tei = np.random.rand(4, 4, 4, 4)
+    with pytest.raises(NameError):
+        h2.hamiltonian("ihpc")
+
+
+# Commenting out since not actively supporting PSI4 at the moment
+# @pytest.mark.skip(reason="Psi4 doesn't offer pip install, so needs to be installed through conda or manually.")
+# def test_run_psi4():
+#     """PSI4 driver"""
+#     # Hardcoded benchmark results
+#     h2_ref_energy = -1.117349035
+#     h2_ref_hcore = np.array([[-1.14765024, -1.00692423], [-1.00692423, -1.14765024]])
+#
+#     h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+#     h2.run_psi4()
+#
+#     assert h2.e_hf == pytest.approx(h2_ref_energy)
+#     assert np.allclose(h2.hcore, h2_ref_hcore)
+
+
+def test_expectation_value():
+    """Tests generation of molecular Hamiltonian and its expectation value using a JW-HF circuit"""
+    # Hardcoded benchmark results
+    h2_ref_energy = -1.117349035
+
+    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    h2.run_pyscf()
+
+    # JW-HF circuit
+    circuit = models.Circuit(h2.nso)
+    circuit.add(gates.X(_i) for _i in range(h2.nelec))
+    # Molecular Hamiltonian and the HF expectation value
+    hamiltonian = h2.hamiltonian()
+    hf_energy = expectation(circuit, hamiltonian)
+
+    # assert h2.e_hf == pytest.approx(hf_energy)
+    assert h2_ref_energy == pytest.approx(hf_energy)
+
+
+def test_eigenvalues():
+    dummy = Molecule()
+    # FermionOperator test:
+    ferm_ham = sum(
+        openfermion.FermionOperator(f"{_i}^ {_i}") + openfermion.FermionOperator(f"{_i} {_i}^") for _i in range(4)
+    )
+    ham_matrix = openfermion.get_sparse_operator(ferm_ham).toarray()
+    assert np.allclose(dummy.eigenvalues(ferm_ham), sorted(np.linalg.eigvals(ham_matrix))[:6])
+
+    # QubitOperator test:
+    qubit_ham = openfermion.QubitOperator("Z0 Z1")
+    ham_matrix = openfermion.get_sparse_operator(qubit_ham).toarray()
+    assert np.allclose(dummy.eigenvalues(qubit_ham), sorted(np.kron(np.array([1.0, -1.0]), np.array([1.0, -1.0])))[:2])
+
+    # SymbolicHamiltonian test:
+    sym_ham = SymbolicHamiltonian(Z(0) * Z(1))
+    assert np.allclose(dummy.eigenvalues(sym_ham), sorted(np.kron(np.array([1.0, -1.0]), np.array([1.0, -1.0]))))
+
+    # Unknown Hamiltonian type
+    with pytest.raises(TypeError):
+        dummy.eigenvalues(0.0)

From 96d5334825f4dc6f0333503898b2a10e2b46a62a Mon Sep 17 00:00:00 2001
From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com>
Date: Wed, 7 Feb 2024 08:58:31 +0000
Subject: [PATCH 03/11] Add code for active argument (only)

---
 src/qibochem/driver/molecule.py | 77 ++++++++++++++++++++++++++++-----
 1 file changed, 65 insertions(+), 12 deletions(-)

diff --git a/src/qibochem/driver/molecule.py b/src/qibochem/driver/molecule.py
index f17930a..4079e30 100644
--- a/src/qibochem/driver/molecule.py
+++ b/src/qibochem/driver/molecule.py
@@ -264,34 +264,87 @@ def _inactive_fock_matrix(self, frozen):
                     inactive_fock[_p][_q] += 2 * self.tei[_orb][_p][_q][_orb] - self.tei[_orb][_p][_orb][_q]
         return inactive_fock
 
-    def hf_embedding(self, active=None, frozen=None):
+    def _active_space(self, active, frozen):
         """
-        Turns on HF embedding for a given active/frozen space, and fills in the class attributes:
-            ``inactive_energy``, ``embed_oei``, and ``embed_tei``.
+        Helper function to check the input for active/frozen space and define the default values
+        for them where necessary
 
         Args:
             active: Iterable representing the active-space for quantum simulation
             frozen: Iterable representing the occupied orbitals to be removed from the simulation
+
+        Returns:
+            _active, _frozen: Iterables representing the active/frozen space
         """
-        # Default arguments for active and frozen
+        n_orbs = self.norb
+        n_occ_orbs = self.nalpha
+
+        _active, _frozen = None, None
         if active is None:
+            if frozen is None:
+                # Default active: Full set of orbitals, frozen: empty list
+                _active = list(range(n_orbs))
+                _frozen = []
+        else:
+            # Check that active argument is valid
+            assert max(active) < n_orbs and min(active) >= 0, "Active space must be between 0 and the number of MOs"
+            _active = active
+            # active, frozen arguments both given
+            if frozen is not None:
+                # Check that active/frozen arguments don't overlap
+                assert not (set(active) & set(frozen)), "Active and frozen space cannot overlap"
+                # Frozen space must be occupied orbitals
+                assert max(frozen) + 1 < n_occ_orbs and min(frozen) >= 0, "Frozen orbital must be occupied orbitals"
+                # All occupied orbitals have to be in active or frozen
+                assert all(
+                    _occ in set(active + frozen) for _occ in range(n_occ_orbs)
+                ), "All occupied orbitals have to be in either the active or frozen space"
+                # Hopefully no more problems with the input
+                _frozen = frozen
+            # frozen argument not given
+            else:
+                # Default frozen: All occupied orbitals not in active
+                _frozen = [_i for _i in range(n_occ_orbs) if _i not in _active]
+
+        """
             if self.active is None:
-                active = list(range(self.norb))
+                _active = list(range(self.norb))
             else:
-                active = self.active
+                _active = self.active
         if frozen is None:
             if self.frozen is None:
-                frozen = [_i for _i in range(self.nalpha) if _i not in active]
+                _frozen = [_i for _i in range(self.nalpha) if _i not in _active]
             else:
-                frozen = self.frozen
+                _frozen = self.frozen
 
         # Check that arguments are valid
-        assert max(active) < self.norb and min(active) >= 0, "Active space must be between 0 " "and the number of MOs"
-        if frozen:
-            assert not (set(active) & set(frozen)), "Active and frozen space cannot overlap"
-            assert max(frozen) + 1 < self.nelec // 2 and min(frozen) >= 0, (
+        assert max(_active) < self.norb and min(_active) >= 0, "Active space must be between 0 " "and the number of MOs"
+        if _frozen:
+            assert not (set(_active) & set(_frozen)), "Active and frozen space cannot overlap"
+            assert max(_frozen) + 1 < self.nelec // 2 and min(_frozen) >= 0, (
                 "Frozen orbitals must" " be occupied orbitals"
             )
+        """
+        return _active, _frozen
+
+    def hf_embedding(self, active=None, frozen=None):
+        """
+        Turns on HF embedding for a given active/frozen space, and fills in the class attributes:
+            ``inactive_energy``, ``embed_oei``, and ``embed_tei``.
+
+        Args:
+            active: Iterable representing the active-space for quantum simulation
+            frozen: Iterable representing the occupied orbitals to be removed from the simulation
+        """
+        # Default arguments for active and frozen if no arguments given
+        if active is None and frozen is None:
+            _active, _frozen = self._active_space(self.active, self.frozen)
+        else:
+            # active/frozen arguments given, process them using _active_space similarly
+            _active, _frozen = self._active_space(active, frozen)
+        # Update the class attributes with the checked arguments
+        self.active = _active
+        self.frozen = _frozen
 
         # Build the inactive Fock matrix first
         inactive_fock = self._inactive_fock_matrix(frozen)

From 801699936ad09515dfefc1a590cfd18b57d3545e Mon Sep 17 00:00:00 2001
From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com>
Date: Wed, 7 Feb 2024 09:29:16 +0000
Subject: [PATCH 04/11] Finish function to define active space rigorously

which hopefully works...
---
 src/qibochem/driver/molecule.py | 47 +++++++++++++--------------------
 1 file changed, 19 insertions(+), 28 deletions(-)

diff --git a/src/qibochem/driver/molecule.py b/src/qibochem/driver/molecule.py
index 4079e30..63b16b3 100644
--- a/src/qibochem/driver/molecule.py
+++ b/src/qibochem/driver/molecule.py
@@ -281,50 +281,41 @@ def _active_space(self, active, frozen):
 
         _active, _frozen = None, None
         if active is None:
+            # No arguments given
             if frozen is None:
                 # Default active: Full set of orbitals, frozen: empty list
                 _active = list(range(n_orbs))
                 _frozen = []
+            # Only frozen argument given
+            else:
+                if frozen:
+                    # Non-empty frozen space must be occupied orbitals
+                    assert max(frozen) + 1 < n_occ_orbs and min(frozen) >= 0, "Frozen orbital must be occupied orbitals"
+                _frozen = frozen
+                # Default active: All orbitals not in frozen
+                _active = [_i for _i in range(n_orbs) if _i not in _frozen]
+        # active argument given
         else:
             # Check that active argument is valid
             assert max(active) < n_orbs and min(active) >= 0, "Active space must be between 0 and the number of MOs"
             _active = active
-            # active, frozen arguments both given
-            if frozen is not None:
+            # frozen argument not given
+            if frozen is None:
+                # Default frozen: All occupied orbitals not in active
+                _frozen = [_i for _i in range(n_occ_orbs) if _i not in _active]
+            # active, frozen arguments both given:
+            else:
                 # Check that active/frozen arguments don't overlap
                 assert not (set(active) & set(frozen)), "Active and frozen space cannot overlap"
-                # Frozen space must be occupied orbitals
-                assert max(frozen) + 1 < n_occ_orbs and min(frozen) >= 0, "Frozen orbital must be occupied orbitals"
+                if frozen:
+                    # Non-empty frozen space must be occupied orbitals
+                    assert max(frozen) + 1 < n_occ_orbs and min(frozen) >= 0, "Frozen orbital must be occupied orbitals"
                 # All occupied orbitals have to be in active or frozen
                 assert all(
                     _occ in set(active + frozen) for _occ in range(n_occ_orbs)
                 ), "All occupied orbitals have to be in either the active or frozen space"
                 # Hopefully no more problems with the input
                 _frozen = frozen
-            # frozen argument not given
-            else:
-                # Default frozen: All occupied orbitals not in active
-                _frozen = [_i for _i in range(n_occ_orbs) if _i not in _active]
-
-        """
-            if self.active is None:
-                _active = list(range(self.norb))
-            else:
-                _active = self.active
-        if frozen is None:
-            if self.frozen is None:
-                _frozen = [_i for _i in range(self.nalpha) if _i not in _active]
-            else:
-                _frozen = self.frozen
-
-        # Check that arguments are valid
-        assert max(_active) < self.norb and min(_active) >= 0, "Active space must be between 0 " "and the number of MOs"
-        if _frozen:
-            assert not (set(_active) & set(_frozen)), "Active and frozen space cannot overlap"
-            assert max(_frozen) + 1 < self.nelec // 2 and min(_frozen) >= 0, (
-                "Frozen orbitals must" " be occupied orbitals"
-            )
-        """
         return _active, _frozen
 
     def hf_embedding(self, active=None, frozen=None):

From 0ed41320e40dfb6c98c06c32272cf8b580dd696e Mon Sep 17 00:00:00 2001
From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com>
Date: Wed, 7 Feb 2024 09:32:09 +0000
Subject: [PATCH 05/11] Add tests for new _active_space function

---
 tests/test_molecule.py | 30 ++++++++++++++++++++++++++++++
 1 file changed, 30 insertions(+)

diff --git a/tests/test_molecule.py b/tests/test_molecule.py
index d587585..c2c4768 100644
--- a/tests/test_molecule.py
+++ b/tests/test_molecule.py
@@ -60,6 +60,36 @@ def test_molecule_custom_basis():
     assert np.isclose(mol.e_hf, -7.94129296352493)
 
 
+def test_define_active_space():
+    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.2))])
+    mol.nalpha = 2
+    mol.norb = 6
+    # Default arguments: Nothing given
+    assert mol._active_space(None, None) == (list(range(mol.norb)), [])
+    # Default frozen argument if active given
+    assert mol._active_space([1, 2, 5], None) == ([1, 2, 5], [0])
+    # Default active argument if frozen given
+    assert mol._active_space(None, [0]) == (list(range(1, 6)), [0])
+    # active, frozen arguments both given
+    assert mol._active_space([0, 1, 2, 3], []) == (list(range(4)), [])
+
+
+def test_define_active_space_assertions():
+    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.2))])
+    mol.nalpha = 2
+    mol.norb = 6
+
+    # Invalid active argument
+    with pytest.raises(AssertionError):
+        _ = mol._active_space([10], None)
+    # Invalid frozen argument
+    with pytest.raises(AssertionError):
+        _ = mol._active_space(None, [100])
+    # active/frozen spaces overlap
+    with pytest.raises(AssertionError):
+        _ = mol._active_space([0, 1], [0])
+
+
 def test_hf_embedding():
     mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.2))])
     mol.run_pyscf()

From 2ef5fd48e9677b32663015d1723ba193895c6cc1 Mon Sep 17 00:00:00 2001
From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com>
Date: Wed, 7 Feb 2024 09:49:23 +0000
Subject: [PATCH 06/11] Bugfix: Correct hf_embedding function

Oops, forgot to do so and test locally before pushing
---
 src/qibochem/driver/molecule.py | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/src/qibochem/driver/molecule.py b/src/qibochem/driver/molecule.py
index 63b16b3..a9615d5 100644
--- a/src/qibochem/driver/molecule.py
+++ b/src/qibochem/driver/molecule.py
@@ -338,22 +338,20 @@ def hf_embedding(self, active=None, frozen=None):
         self.frozen = _frozen
 
         # Build the inactive Fock matrix first
-        inactive_fock = self._inactive_fock_matrix(frozen)
+        inactive_fock = self._inactive_fock_matrix(self.frozen)
 
         # Calculate the inactive Fock energy
         # Only want frozen part of original OEI and inactive Fock matrix
-        _oei = self.oei[np.ix_(frozen, frozen)]
-        _inactive_fock = inactive_fock[np.ix_(frozen, frozen)]
+        _oei = self.oei[np.ix_(self.frozen, self.frozen)]
+        _inactive_fock = inactive_fock[np.ix_(self.frozen, self.frozen)]
         self.inactive_energy = np.einsum("ii->", _oei + _inactive_fock)
 
         # Keep only the active part
-        self.embed_oei = inactive_fock[np.ix_(active, active)]
-        self.embed_tei = self.tei[np.ix_(active, active, active, active)]
+        self.embed_oei = inactive_fock[np.ix_(self.active, self.active)]
+        self.embed_tei = self.tei[np.ix_(self.active, self.active, self.active, self.active)]
 
-        # Update class attributes
-        self.active = active
-        self.frozen = frozen
-        self.n_active_orbs = 2 * len(active)
+        # Update other class attributes
+        self.n_active_orbs = 2 * len(self.active)
         self.n_active_e = self.nelec - 2 * len(self.frozen)
 
     def hamiltonian(

From 5f4a8f4d3101748664bb2f4c61764e907bdfa31f Mon Sep 17 00:00:00 2001
From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com>
Date: Thu, 8 Feb 2024 07:52:04 +0000
Subject: [PATCH 07/11] Simplify hf_embedding test

Change LiH to H2, check removal of all virtual orbitals
---
 tests/test_molecule.py | 23 ++++++++++-------------
 1 file changed, 10 insertions(+), 13 deletions(-)

diff --git a/tests/test_molecule.py b/tests/test_molecule.py
index c2c4768..58a9b6b 100644
--- a/tests/test_molecule.py
+++ b/tests/test_molecule.py
@@ -91,21 +91,18 @@ def test_define_active_space_assertions():
 
 
 def test_hf_embedding():
-    mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.2))])
+    mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
     mol.run_pyscf()
-    ref_oei = mol.oei
-    ref_tei = mol.tei
-    mol.hf_embedding()
-    embed_oei = mol.embed_oei
-    embed_tei = mol.embed_tei
-    assert np.allclose(embed_oei, ref_oei)
-    assert np.allclose(embed_tei, ref_tei)
-
-    mol.frozen = [0]
-    mol.active = [1, 2]
-    mol.hf_embedding()
-    assert mol.n_active_orbs == 4
+    # Remove all virtual orbitals from the active space
+    mol.hf_embedding(active=[0])
+    # Check that the class attributes have been updated correctly
+    assert mol.frozen == []
+    assert mol.n_active_orbs == 2
     assert mol.n_active_e == 2
+    # OEI/TEI (in MO basis) for the occupied orbitals should remain unchanged
+    dim = mol.n_active_orbs // 2
+    assert np.allclose(mol.embed_oei, mol.oei[:dim, :dim])
+    assert np.allclose(mol.embed_tei, mol.tei[:dim, :dim, :dim, :dim])
 
 
 def test_fermionic_hamiltonian():

From 4f1464ae35ab4bbb9d80eb09c657e155b0fa0a58 Mon Sep 17 00:00:00 2001
From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com>
Date: Thu, 8 Feb 2024 09:32:49 +0000
Subject: [PATCH 08/11] Add tutorial for HF embedding into docs (draft)

---
 doc/source/tutorials/ansatz.rst   |  1 +
 doc/source/tutorials/molecule.rst | 34 ++++++++++++++++++++++++++++++-
 2 files changed, 34 insertions(+), 1 deletion(-)

diff --git a/doc/source/tutorials/ansatz.rst b/doc/source/tutorials/ansatz.rst
index ae3ff31..3352179 100644
--- a/doc/source/tutorials/ansatz.rst
+++ b/doc/source/tutorials/ansatz.rst
@@ -88,6 +88,7 @@ The following example demonstrates how the energy of the H2 molecule is affected
     |   0.2 |   0.673325849299  |
     -----------------------------
 
+.. _UCC Ansatz:
 
 Unitary Coupled Cluster Ansatz
 ------------------------------
diff --git a/doc/source/tutorials/molecule.rst b/doc/source/tutorials/molecule.rst
index 198c21f..4df364f 100644
--- a/doc/source/tutorials/molecule.rst
+++ b/doc/source/tutorials/molecule.rst
@@ -46,9 +46,41 @@ Qibochem offers the functionality to interface with `PySCF <https://pyscf.org/>`
     # Running PySCF
     h2.run_pyscf()
 
-The default level of theory is HF/STO-3G, and upon executing the PySCF driver for a given molecule, several molecular quantities are calculated and stored in the Molecule class.
+The default level of theory is HF/STO-3G, and upon executing the PySCF driver for a given molecule, several molecular quantities are calculated and stored in the ``Molecule`` class.
 These include the:
 
 * converged Hartree-Fock energy
 * optimized molecular orbital (MO) coefficients
 * one- and two-electron integrals
+
+
+Embedding the quantum electronic structure calculation
+------------------------------------------------------
+
+In quantum chemistry, most *ab initio* calculations start with a Hartree-Fock (HF) calculation.
+The obtained HF wave function is then used as a starting point to apply post-HF methods to improve upon the treatment of electron correlation.
+An example of a post-HF method that can be run on a quantum computer is the :ref:`Unitary Coupled Cluster method <UCC Ansatz>`.
+Unfortunately, the current level of quantum hardware are unable to run these methods for molecules that are too large.
+
+One possible approach to reduce the hardware required is to embed the quantum simulation into a classically computed environment.
+(see `Rossmanek et al. <https://doi.org/10.1063/5.0029536/>`_)
+Essentially, the occupancy of the core *1s* orbitals for the heavy atoms in a molecule is effectively constant; the same might hold for higher virtual orbitals.
+As such, these orbitals can be removed from the simulation of the molecule without there being a significant change in the molecular properties of interest.
+The remaining orbitals left in the simulation are then known as the active space of the molecule.
+
+An example of how to apply this using Qibochem for the LiH molecule is given below.
+
+.. code-block:: python
+
+    from qibochem.driver.molecule import Molecule
+
+    # Inline definition of H2
+    h2 = Molecule([('H', (0.0, 0.0, 0.0)), ('H', (0.0, 0.0, 0.74804))])
+    # Running PySCF
+    h2.run_pyscf()
+    frozen = [0] # The electrons on the 1s orbital of Li in LiH are frozen - removed from the quantum simulation
+    h2.hf_embedding(frozen=frozen)
+
+The code block above will re-calculate the one- and two-electron integrals for the given active space, and store the result back in the ``Molecule`` class.
+
+.. Note: Not sure if this is the best place to describe HF embedding? But where else would be good?

From f45545e72338e82364846bb965e15dbd2c20d996 Mon Sep 17 00:00:00 2001
From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com>
Date: Thu, 8 Feb 2024 09:42:44 +0000
Subject: [PATCH 09/11] Edit Hamiltonian tutorial for HF embedding

---
 doc/source/tutorials/hamiltonian.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/doc/source/tutorials/hamiltonian.rst b/doc/source/tutorials/hamiltonian.rst
index 3546cc7..99a6970 100644
--- a/doc/source/tutorials/hamiltonian.rst
+++ b/doc/source/tutorials/hamiltonian.rst
@@ -74,6 +74,7 @@ Lastly, to carry out quantum simulations of the molecular electronic structure u
     -0.107280411608667 - 0.0454153248117008*X0*X1*Y2*Y3 + 0.0454153248117008*X0*Y1*Y2*X3 + 0.0454153248117008*Y0*X1*X2*Y3 - 0.0454153248117008*Y0*Y1*X2*X3 + 0.170182611817142*Z0 + 0.16830546187935*Z0*Z1 + 0.120163790529127*Z0*Z2 + 0.165579115340828*Z0*Z3 + 0.170182611817142*Z1 + 0.165579115340828*Z1*Z2 + 0.120163790529127*Z1*Z3 - 0.219750654392482*Z2 + 0.174043557614182*Z2*Z3 - 0.219750654392482*Z3
 
 By default, the molecular Hamiltonian is returned as a ``SymbolicHamiltonian``, `i.e.` if no arguments are given in :code:`mol.hamiltonian()`.
+In addition, if HF embedding has been applied, the embedded values of the one-/two- electron integrals will be used to construct the molecular Hamiltonian as well.
 
 Otherwise, using the ``"f"``/``"ferm"`` and ``"q"``/``"qubit"`` arguments will return the molecular Hamiltonian as an OpenFermion ``FermionOperator`` and ``QubitOperator`` respectively.
 Additional information about the data structure of these two classes can be found `here <https://quantumai.google/openfermion/tutorials/intro_to_openfermion>`_.

From b2831bb1f0a89ecff6a795e36c2b72c8335a9cd8 Mon Sep 17 00:00:00 2001
From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com>
Date: Thu, 8 Feb 2024 09:56:01 +0000
Subject: [PATCH 10/11] Add tests for code coverage

---
 tests/test_molecule.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/tests/test_molecule.py b/tests/test_molecule.py
index 58a9b6b..198b45b 100644
--- a/tests/test_molecule.py
+++ b/tests/test_molecule.py
@@ -72,6 +72,8 @@ def test_define_active_space():
     assert mol._active_space(None, [0]) == (list(range(1, 6)), [0])
     # active, frozen arguments both given
     assert mol._active_space([0, 1, 2, 3], []) == (list(range(4)), [])
+    # active, frozen arguments both given
+    assert mol._active_space([1, 2, 3], [0]) == (list(range(1, 4)), [0])
 
 
 def test_define_active_space_assertions():
@@ -91,10 +93,10 @@ def test_define_active_space_assertions():
 
 
 def test_hf_embedding():
-    mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))], active=[0])
     mol.run_pyscf()
     # Remove all virtual orbitals from the active space
-    mol.hf_embedding(active=[0])
+    mol.hf_embedding()
     # Check that the class attributes have been updated correctly
     assert mol.frozen == []
     assert mol.n_active_orbs == 2

From 665138697401b9279ed8d98017be8c0a38415c29 Mon Sep 17 00:00:00 2001
From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com>
Date: Wed, 14 Feb 2024 02:34:19 +0000
Subject: [PATCH 11/11] Update Molecule API documentation

---
 src/qibochem/driver/molecule.py | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/qibochem/driver/molecule.py b/src/qibochem/driver/molecule.py
index a9615d5..9cba67a 100644
--- a/src/qibochem/driver/molecule.py
+++ b/src/qibochem/driver/molecule.py
@@ -363,16 +363,23 @@ def hamiltonian(
         ferm_qubit_map=None,
     ):
         """
-        Builds a molecular Hamiltonian using the one-/two- electron integrals
+        Builds a molecular Hamiltonian using the one-/two- electron integrals. If HF embedding has been applied,
+        (i.e. the ``embed_oei``, ``embed_tei``, and ``inactive_energy`` attributes are all not ``None``), the
+        corresponding values for the molecular integrals will be used instead.
 
         Args:
             ham_type: Format of molecular Hamiltonian returned. The available options are:
                 ``("f", "ferm")``: OpenFermion ``FermionOperator``,
                 ``("q", "qubit")``: OpenFermion ``QubitOperator``, or
                 ``("s", "sym")``: Qibo ``SymbolicHamiltonian`` (default)
-            oei: 1-electron integrals. Default: ``self.oei`` (MO basis)
-            tei: 2-electron integrals in 2ndQ notation. Default: ``self.tei`` (MO basis)
-            constant: For inactive Fock energy if embedding used. Default: 0.0
+            oei: 1-electron integrals (in the MO basis). The default value is the ``oei`` class attribute , unless
+                the ``embed_oei`` attribute exists and is not ``None``, then ``embed_oei`` is used.
+            tei: 2-electron integrals in the second-quantization notation (and MO basis). The default value is the
+                ``tei`` class attribute , unless the ``embed_tei`` attribute exists and is not ``None``, then ``embed_tei``
+                is used.
+            constant: Constant value to be added to the electronic energy. Mainly used for adding the inactive Fock
+                energy if HF embedding was applied. Default: 0.0, unless the ``inactive_energy`` class attribute exists
+                and is not ``None``, then ``inactive_energy`` is used.
             ferm_qubit_map: Which fermion to qubit transformation to use.
                 Must be either ``jw`` (default) or ``bk``