Repository: PennyLaneAI/pennylane
Pull request #2042 (2021-12-17 14:51:08)

Add function to taper Hartree Fock state for qubit tapering
**Context:**
This PR adds function to obtain the tapered Hartree Fock state corresponding to the tapered Hamiltonian.

**Description of the Change:**
1. Adds `taper_hartree_fock` function, which tapers HF state based on symmetries generated from the molecular hamiltonian.
2. Modifies `_generate_qubit_operator`  function, to return single creation/anhilation operator under JW transform.  

**Benefits:**
**Possible Drawbacks:**
**Related GitHub Issues:**


github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#2042](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (297bf7a) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/a3f2611760b270c01b863e94f6d7eb706f6afcb4?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (a3f2611) will **decrease** coverage by `2.05%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@            Coverage Diff             @@
##           master    #2042      +/-   ##
==========================================
- Coverage   99.16%   97.10%   -2.06%     
==========================================
  Files         222      225       +3     
  Lines       16942    17339     +397     
==========================================
+ Hits        16800    16837      +37     
- Misses        142      502     +360     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/hf/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2hmL19faW5pdF9fLnB5) | `100.00% <ø> (ø)` | |
| [pennylane/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL19faW5pdF9fLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/hf/hamiltonian.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2hmL2hhbWlsdG9uaWFuLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/hf/tapering.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2hmL3RhcGVyaW5nLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/optimize/qng.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wdGltaXplL3FuZy5weQ==) | `100.00% <100.00%> (ø)` | |
| [pennylane/templates/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RlbXBsYXRlcy9fX2luaXRfXy5weQ==) | `100.00% <100.00%> (ø)` | |
| [pennylane/templates/subroutines/qpe.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RlbXBsYXRlcy9zdWJyb3V0aW5lcy9xcGUucHk=) | `100.00% <100.00%> (ø)` | |
| [pennylane/templates/tensornetworks/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RlbXBsYXRlcy90ZW5zb3JuZXR3b3Jrcy9fX2luaXRfXy5weQ==) | `100.00% <100.00%> (ø)` | |
| [pennylane/templates/tensornetworks/mps.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RlbXBsYXRlcy90ZW5zb3JuZXR3b3Jrcy9tcHMucHk=) | `100.00% <100.00%> (ø)` | |
| [pennylane/transforms/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvX19pbml0X18ucHk=) | `100.00% <100.00%> (ø)` | |
| ... and [32 more](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [a3f2611...297bf7a](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2042?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).


Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 07da939ec..54c9b5d1f 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -151,25 +151,39 @@
   [(#1974)](https://github.com/PennyLaneAI/pennylane/pull/1974)
   [(#2041)](https://github.com/PennyLaneAI/pennylane/pull/2041)
 
-  With this functionality, a molecular Hamiltonian can be transformed to a new Hamiltonian that acts
-  on a reduced number of qubits.
+  With this functionality, a molecular Hamiltonian and the corresponding Hartree-Fock (HF) state can be transformed to a new Hamiltonian and HF state that acts on a reduced number of qubits, respectively.
 
   ```python
-  symbols = ["H", "H"]
-  geometry = np.array([[0.0, 0.0, -0.69440367], [0.0, 0.0, 0.69440367]])
-  mol = qml.hf.Molecule(symbols, geometry)
-  H = qml.hf.generate_hamiltonian(mol)(geometry)
-  generators, paulix_ops = qml.hf.generate_symmetries(H, len(H.wires))
-  paulix_sector = qml.hf.optimal_sector(H, generators, mol.n_electrons)
-  H_tapered = qml.hf.transform_hamiltonian(H, generators, paulix_ops, paulix_sector)
-  ```
+  from pennylane import hf
+  from pennylane import numpy as np
 
+  symbols = ["He", "H"]
+  geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.4588684632]])
+  mol = hf.Molecule(symbols, geometry, charge=1)
+  H = hf.generate_hamiltonian(mol)(geometry)
+  n_qubits, n_elec = len(H.wires), mol.n_electrons
+
+  generators, paulix_ops = hf.generate_symmetries(H, n_qubits)
+  opt_sector = hf.optimal_sector(H, generators, n_elec)
+  H_tapered = hf.transform_hamiltonian(H, generators, paulix_ops, opt_sector)
+  hf_tapered = hf.transform_hf(generators, paulix_ops, paulix_sector,
+                                    n_elec, n_qubits)
+  ```
   ```pycon
   >>> print(H_tapered)
-    ((-0.321034397355719+0j)) [I0]
-  + ((0.1809270275619743+0j)) [X0]
-  + ((0.7959678503870796+0j)) [Z0]
+    ((-1.7997297644914574+0j)) [I0]
+  + ((-0.10492941956079854+0j)) [X0]
+  + ((0.10492941956079856+0j)) [X1]
+  + ((0.5675134088336165+0j)) [Z1]
+  + ((0.5675134088336168+0j)) [Z0]
+  + ((-0.14563730440190722+0j)) [Y0 Y1]
+  + ((-0.10492941933657857+0j)) [X0 Z1]
+  + ((0.09337410512815508+0j)) [Z0 Z1]
+  + ((0.10492941933657857+0j)) [Z0 X1]
+  >>> print(hf_tapered)
+  tensor([1, 1], requires_grad=True)
   ```
+  
 
 * Added the adjoint method for the metric tensor.
   [(#1992)](https://github.com/PennyLaneAI/pennylane/pull/1992)
@@ -537,6 +551,6 @@
 
 This release contains contributions from (in alphabetical order):
 
-Juan Miguel Arrazola, Ali Asadi, Esther Cruz, Christian Gogolin, Christina Lee, Olivia Di Matteo, Diego Guala,
-Anthony Hayes, Edward Jiang, Josh Izaac, Ankit Khandelwal, Korbinian Kottmann, Jay Soni, Antal Száva,
+Juan Miguel Arrazola, Ali Asadi, Utkarsh Azad, Esther Cruz, Christian Gogolin Christina Lee, Olivia Di Matteo, Diego Guala,
+Anthony Hayes, Josh Izaac, Soran Jahangiri, Edward Jiang, Ankit Khandelwal, Korbinian Kottmann, Jay Soni, Antal Száva,
 David Wierichs, Shaoming Zhang
diff --git a/pennylane/grouping/utils.py b/pennylane/grouping/utils.py
index 8fbc52472..b7d8ccdf3 100644
--- a/pennylane/grouping/utils.py
+++ b/pennylane/grouping/utils.py
@@ -445,7 +445,7 @@ def string_to_pauli_word(pauli_string, wire_map=None):
     # Special case: all-identity Pauli
     if pauli_string == "I" * len(wire_map):
         first_wire = list(wire_map.keys())[0]
-        return Identity(wire_map[first_wire])
+        return Identity(first_wire)
 
     pauli_word = None
 
diff --git a/pennylane/hf/__init__.py b/pennylane/hf/__init__.py
index 4d6410ce2..5bf4a451c 100644
--- a/pennylane/hf/__init__.py
+++ b/pennylane/hf/__init__.py
@@ -49,6 +49,7 @@ from .tapering import (
     generate_paulis,
     generate_symmetries,
     get_generators,
-    transform_hamiltonian,
     optimal_sector,
+    transform_hamiltonian,
+    transform_hf,
 )
diff --git a/pennylane/hf/hamiltonian.py b/pennylane/hf/hamiltonian.py
index 8cf60b5e8..160a7651d 100644
--- a/pennylane/hf/hamiltonian.py
+++ b/pennylane/hf/hamiltonian.py
@@ -308,6 +308,9 @@ def _generate_qubit_operator(op):
     >>> q
     ([(0.5+0j), (-0.5+0j)], [[], [(0, 'Z')]]) # corresponds to :math:`\frac{1}{2}(I_0 - Z_0)`
     """
+    if len(op) == 1:
+        op = [((op[0], 1),)]
+
     if len(op) == 2:
         op = [((op[0], 1), (op[1], 0))]
 
diff --git a/pennylane/hf/tapering.py b/pennylane/hf/tapering.py
index 12ab39da6..b1ef7c949 100644
--- a/pennylane/hf/tapering.py
+++ b/pennylane/hf/tapering.py
@@ -21,29 +21,38 @@ import numpy
 import autograd.numpy as anp
 import pennylane as qml
 from pennylane import numpy as np
+from pennylane.hf.hamiltonian import _generate_qubit_operator
 
 
-def _binary_matrix(terms, num_qubits):
+def _binary_matrix(terms, num_qubits, wire_map=None):
     r"""Get a binary matrix representation of the Hamiltonian where each row corresponds to a
     Pauli term, which is represented by a concatenation of Z and X vectors.
 
     Args:
         terms (Iterable[Observable]): operators defining the Hamiltonian
         num_qubits (int): number of wires required to define the Hamiltonian
+        wire_map (dict): dictionary containing all wire labels used in the Pauli words as keys, and
+            unique integer labels as their values
+
     Returns:
         array[int]: binary matrix representation of the Hamiltonian of shape
         :math:`len(terms) \times 2*num_qubits`
 
     **Example**
 
-    >>> terms = [qml.PauliZ(wires=[0]) @ qml.PauliX(wires=[1]),
-    ...          qml.PauliZ(wires=[0]) @ qml.PauliY(wires=[2]),
-    ...          qml.PauliX(wires=[0]) @ qml.PauliY(wires=[3])]
-    >>> _binary_matrix(terms, 4)
+    >>> wire_map = {'a':0, 'b':1, 'c':2, 'd':3}
+    >>> terms = [qml.PauliZ(wires=['a']) @ qml.PauliX(wires=['b']),
+    ...          qml.PauliZ(wires=['a']) @ qml.PauliY(wires=['c']),
+    ...          qml.PauliX(wires=['a']) @ qml.PauliY(wires=['d'])]
+    >>> _binary_matrix(terms, 4, wire_map=wire_map)
     array([[1, 0, 0, 0, 0, 1, 0, 0],
            [1, 0, 1, 0, 0, 0, 1, 0],
            [0, 0, 0, 1, 1, 0, 0, 1]])
     """
+    if wire_map is None:
+        all_wires = qml.wires.Wires.all_wires([term.wires for term in terms], sort=True)
+        wire_map = {i: c for c, i in enumerate(all_wires)}
+
     binary_matrix = np.zeros((len(terms), 2 * num_qubits), dtype=int)
     for idx, term in enumerate(terms):
         ops, wires = term.name, term.wires
@@ -51,9 +60,9 @@ def _binary_matrix(terms, num_qubits):
             ops = [ops]
         for op, wire in zip(ops, wires):
             if op in ["PauliX", "PauliY"]:
-                binary_matrix[idx][wire + num_qubits] = 1
+                binary_matrix[idx][wire_map[wire] + num_qubits] = 1
             if op in ["PauliZ", "PauliY"]:
-                binary_matrix[idx][wire] = 1
+                binary_matrix[idx][wire_map[wire]] = 1
 
     return binary_matrix
 
@@ -217,7 +226,8 @@ def generate_paulis(generators, num_qubits):
     for row in range(bmat.shape[0]):
         bmatrow = bmat[row]
         bmatrest = np.delete(bmat, row, axis=0)
-        for col in range(bmat.shape[1] // 2):
+        # reversing the order to priortize removing higher index wires first
+        for col in range(bmat.shape[1] // 2)[::-1]:
             # Anti-commutes with the (row) and commutes with all other symmetries.
             if bmatrow[col] and np.array_equal(
                 bmatrest[:, col], np.zeros(bmat.shape[0] - 1, dtype=int)
@@ -343,6 +353,7 @@ def _simplify(h, cutoff=1.0e-12):
     """
     s = []
     wiremap = dict(zip(h.wires, range(len(h.wires) + 1)))
+
     for term in h.terms[1]:
         term = qml.operation.Tensor(term).prune()
         s.append(qml.grouping.pauli_word_to_string(term, wire_map=wiremap))
@@ -452,13 +463,16 @@ def transform_hamiltonian(h, generators, paulix_ops, paulix_sector):
                 val[i] *= paulix_sector[idx]
 
     o = []
-    wires_tap = [h.wires[i] for i in range(len(h.wires)) if i not in paulix_wires]
+    wires_tap = [i for i in h.wires if i not in paulix_wires]
     wiremap_tap = dict(zip(wires_tap, range(len(wires_tap) + 1)))
+
     for i in range(len(h.terms[0])):
         s = qml.grouping.pauli_word_to_string(h.terms[1][i], wire_map=wiremap)
         wires = [x for x in h.wires if x not in paulix_wires]
         o.append(
-            qml.grouping.string_to_pauli_word("".join([s[i] for i in wires]), wire_map=wiremap_tap)
+            qml.grouping.string_to_pauli_word(
+                "".join([s[wiremap[i]] for i in wires]), wire_map=wiremap_tap
+            )
         )
 
     c = anp.multiply(val, h.terms[0])
@@ -526,3 +540,76 @@ def optimal_sector(qubit_op, generators, active_electrons):
         perm.append(coeff)
 
     return perm
+
+
+def transform_hf(generators, paulix_ops, paulix_sector, num_electrons, num_wires):
+    r"""Transform a Hartree-Fock state with a Clifford operator and taper qubits.
+
+    The fermionic operators defining the molecule's Hartree-Fock (HF) state are first mapped onto a qubit operator
+    using the Jordan-Wigner encoding. This operator is then transformed using the Clifford operators :math:`U`
+    obtained from the :math:`\mathbb{Z}_2` symmetries of the molecular Hamiltonian resulting in a qubit operator
+    that acts non-trivially only on a subset of qubits. A new, tapered HF state is built on this reduced subset
+    of qubits by placing the qubits which are acted on by a Pauli-X or Pauli-Y operators in state :math:`|1\rangle`
+    and leaving the rest in state :math:`|0\rangle`.
+
+    Args:
+        generators (list[Hamiltonian]): list of generators of symmetries, taus, for the Hamiltonian
+        paulix_ops (list[Operation]):  list of single-qubit Pauli-X operators
+        paulix_sector (list[int]): list of eigenvalues of Pauli-X operators
+        num_electrons (int): number of active electrons in the system for generating the Hartree-Fock bitstring
+        num_wires (int): number of wires in the system for generating the Hartree-Fock bitstring
+
+    Returns:
+        array(int): tapered Hartree-Fock state
+
+    **Example**
+
+    >>> from pennylane import hf
+    >>> symbols = ['He', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.4588684632]])
+    >>> mol = hf.Molecule(symbols, geometry, charge=1)
+    >>> H = hf.generate_hamiltonian(mol)(geometry)
+    >>> n_qubits, n_elec = len(H.wires), mol.n_electrons
+    >>> generators, paulix_ops = hf.generate_symmetries(H, n_qubits)
+    >>> paulix_sector = hf.optimal_sector(H, generators, n_elec)
+    >>> hf.transform_hf(generators, paulix_ops, paulix_sector,
+    ...                 n_elec, n_qubits)
+    tensor([1, 1], requires_grad=True)
+    """
+
+    pauli_map = {"I": qml.Identity, "X": qml.PauliX, "Y": qml.PauliY, "Z": qml.PauliZ}
+
+    # build the untapered Hartree Fock state
+    hf = np.where(np.arange(num_wires) < num_electrons, 1, 0)
+
+    # convert the HF state to a corresponding HF observable under the JW transform
+    fermop_terms = []
+    for idx, bit in enumerate(hf):
+        if bit:
+            op_coeffs, op_str = _generate_qubit_operator([idx])
+            op_terms = []
+            for term in op_str:
+                op_term = pauli_map[term[0][1]](term[0][0])
+                for tm in term[1:]:
+                    op_term @= pauli_map[tm[1]](tm[0])
+                op_terms.append(op_term)
+            op_term = qml.Hamiltonian(np.array(op_coeffs), op_terms)
+        else:
+            op_term = qml.Hamiltonian([1], [qml.Identity(idx)])
+        fermop_terms.append(op_term)
+
+    ferm_op = functools.reduce(lambda i, j: _observable_mult(i, j), fermop_terms)
+
+    # taper the HF observable using the symmetries obtained from the molecular hamiltonian
+    fermop_taper = transform_hamiltonian(ferm_op, generators, paulix_ops, paulix_sector)
+    fermop_mat = _binary_matrix(fermop_taper.ops, len(fermop_taper.wires))
+
+    # iterate over the terms in tapered HF observable and build the tapered HF state
+    tapered_hartree_fock = []
+    for col in fermop_mat.T[fermop_mat.shape[1] // 2 :]:
+        if 1 in col:
+            tapered_hartree_fock.append(1)
+        else:
+            tapered_hartree_fock.append(0)
+
+    return np.array(tapered_hartree_fock).astype(int)


Test patch:
diff --git a/tests/hf/test_tapering.py b/tests/hf/test_tapering.py
index aa88fbb87..9609d7740 100644
--- a/tests/hf/test_tapering.py
+++ b/tests/hf/test_tapering.py
@@ -14,6 +14,8 @@
 """
 Unit tests for functions needed for qubit tapering.
 """
+import functools
+import scipy
 import pytest
 import pennylane as qml
 from pennylane import numpy as np
@@ -27,8 +29,9 @@ from pennylane.hf.tapering import (
     generate_paulis,
     generate_symmetries,
     get_generators,
-    transform_hamiltonian,
     optimal_sector,
+    transform_hamiltonian,
+    transform_hf,
 )
 
 
@@ -86,6 +89,17 @@ from pennylane.hf.tapering import (
                 ]
             ),
         ),
+        (
+            [
+                qml.PauliZ(wires=["a"]) @ qml.PauliX(wires=["b"]),
+                qml.PauliZ(wires=["a"]) @ qml.PauliY(wires=["c"]),
+                qml.PauliX(wires=["a"]) @ qml.PauliY(wires=["d"]),
+            ],
+            4,
+            np.array(
+                [[1, 0, 0, 0, 0, 1, 0, 0], [1, 0, 1, 0, 0, 0, 1, 0], [0, 0, 0, 1, 1, 0, 0, 1]]
+            ),
+        ),
     ],
 )
 def test_binary_matrix(terms, num_qubits, result):
@@ -537,3 +551,143 @@ def test_exceptions_optimal_sector(symbols, geometry, generators, num_electrons,
 
     with pytest.raises(ValueError, match=msg_match):
         optimal_sector(hamiltonian, generators, num_electrons)
+
+
+@pytest.mark.parametrize(
+    ("generators", "paulix_ops", "paulix_sector", "num_electrons", "num_wires", "result"),
+    [
+        (
+            [
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(1)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(2)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(3)]),
+            ],
+            [qml.PauliX(wires=[1]), qml.PauliX(wires=[2]), qml.PauliX(wires=[3])],
+            (1, -1, -1),
+            2,
+            4,
+            np.array([1]),
+        ),
+        (
+            [
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(2) @ qml.PauliZ(4)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(1) @ qml.PauliZ(3) @ qml.PauliZ(5)]),
+            ],
+            [qml.PauliX(wires=[2]), qml.PauliX(wires=[3])],
+            (-1, -1),
+            2,
+            6,
+            np.array([1, 1, 0, 0]),
+        ),
+        (
+            [
+                qml.Hamiltonian([1.0], [qml.PauliZ(6) @ qml.PauliZ(7)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(8) @ qml.PauliZ(9)]),
+                qml.Hamiltonian(
+                    [1.0],
+                    [
+                        qml.PauliZ(wires=[0])
+                        @ qml.PauliZ(wires=[2])
+                        @ qml.PauliZ(wires=[4])
+                        @ qml.PauliZ(wires=[6])
+                        @ qml.PauliZ(wires=[8])
+                        @ qml.PauliZ(wires=[10])
+                    ],
+                ),
+                qml.Hamiltonian(
+                    [1.0],
+                    [
+                        qml.PauliZ(wires=[1])
+                        @ qml.PauliZ(wires=[3])
+                        @ qml.PauliZ(wires=[5])
+                        @ qml.PauliZ(wires=[6])
+                        @ qml.PauliZ(wires=[8])
+                        @ qml.PauliZ(wires=[11])
+                    ],
+                ),
+            ],
+            [
+                qml.PauliX(wires=[7]),
+                qml.PauliX(wires=[9]),
+                qml.PauliX(wires=[0]),
+                qml.PauliX(wires=[1]),
+            ],
+            (1, 1, 1, 1),
+            4,
+            12,
+            np.array([1, 1, 0, 0, 0, 0, 0, 0]),
+        ),
+    ],
+)
+def test_transform_hf(generators, paulix_ops, paulix_sector, num_electrons, num_wires, result):
+    r"""Test that transform_hf returns the correct result."""
+
+    tapered_hf_state = transform_hf(
+        generators,
+        paulix_ops,
+        paulix_sector,
+        num_electrons,
+        num_wires,
+    )
+    assert np.all(tapered_hf_state == result)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "charge"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.40104295]], requires_grad=True),
+            0,
+        ),
+        (
+            ["He", "H"],
+            np.array(
+                [[0.0, 0.0, 0.0], [0.0, 0.0, 1.4588684632]],
+                requires_grad=True,
+            ),
+            1,
+        ),
+        (
+            ["H", "H", "H"],
+            np.array(
+                [[-0.84586466, 0.0, 0.0], [0.84586466, 0.0, 0.0], [0.0, 1.46508057, 0.0]],
+                requires_grad=True,
+            ),
+            1,
+        ),
+    ],
+)
+def test_hf_energy(symbols, geometry, charge):
+    r"""Test that HF energy obtained from the tapered Hamiltonian and tapered Hartree Fock state is consistent."""
+    mol = qml.hf.Molecule(symbols, geometry, charge)
+    hamiltonian = qml.hf.generate_hamiltonian(mol)(geometry)
+    hf_state = np.where(np.arange(len(hamiltonian.wires)) < mol.n_electrons, 1, 0)
+    generators, paulix_ops = generate_symmetries(hamiltonian, len(hamiltonian.wires))
+    paulix_sector = optimal_sector(hamiltonian, generators, mol.n_electrons)
+
+    hamiltonian_tapered = transform_hamiltonian(hamiltonian, generators, paulix_ops, paulix_sector)
+    hf_state_tapered = transform_hf(
+        generators, paulix_ops, paulix_sector, mol.n_electrons, len(hamiltonian.wires)
+    )
+
+    # calculate the HF energy <\psi_{HF}| H |\psi_{HF}> for tapered and untapered Hamiltonian
+    o = np.array([1, 0])
+    l = np.array([0, 1])
+    state = functools.reduce(lambda i, j: np.kron(i, j), [l if s else o for s in hf_state])
+    state_tapered = functools.reduce(
+        lambda i, j: np.kron(i, j), [l if s else o for s in hf_state_tapered]
+    )
+
+    energy = (
+        scipy.sparse.coo_matrix(state)
+        @ qml.utils.sparse_hamiltonian(hamiltonian)
+        @ scipy.sparse.coo_matrix(state).T
+    ).toarray()
+    energy_tapered = (
+        scipy.sparse.coo_matrix(state_tapered)
+        @ qml.utils.sparse_hamiltonian(hamiltonian_tapered)
+        @ scipy.sparse.coo_matrix(state_tapered).T
+    ).toarray()
+
+    assert np.isclose(energy, energy_tapered)