Repository: quantumlib/OpenFermion
Pull request #599 (2020-05-26 18:22:34)

Construct reduced Hamiltonian
The reduced Hamiltonian has only two-body fermionic correlators.  In the
N-particle space this is an isospectral operator to the original Hamiltonian.


Patch:
diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index bf863ef..6141e91 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -38,6 +38,8 @@ from ._operator_utils import (chemist_ordered, count_qubits,
 from ._qubit_tapering_from_stabilizer import (reduce_number_of_terms,
                                               taper_off_qubits)
 
+from ._reduced_hamiltonian import make_reduced_hamiltonian
+
 from ._rdm_mapping_functions import (kronecker_delta,
                                      map_two_pdm_to_two_hole_dm,
                                      map_two_pdm_to_one_pdm,
diff --git a/src/openfermion/utils/_reduced_hamiltonian.py b/src/openfermion/utils/_reduced_hamiltonian.py
new file mode 100644
index 0000000..fb3e27c
--- /dev/null
+++ b/src/openfermion/utils/_reduced_hamiltonian.py
@@ -0,0 +1,46 @@
+from itertools import product
+import numpy as np
+from openfermion.ops import InteractionOperator
+
+
+def make_reduced_hamiltonian(molecular_hamiltonian: InteractionOperator,
+                             n_electrons: int) -> InteractionOperator:
+    r"""
+    Construct the reduced Hamiltonian.
+
+    This Hamiltonian is equivalent to the electronic structure Hamiltonian
+    but contains only two-body terms.  To do this, the operator now depends
+    on the number of particles being simulated.  We use the RDM sum rule to
+    lift the 1-body terms to the two-body space.
+
+    Derivation:
+        use the fact that i^l = (1/(n -1)) sum_{jk}\delta_{jk}i^ j^ k l
+                          i^l = (-1/(n -1)) sum_{jk}\delta_{jk}j^ i^ k l
+                          i^l = (-1/(n -1)) sum_{jk}\delta_{jk}i^ j^ l k
+                          i^l = (1/(n -1)) sum_{jk}\delta_{jk}j^ i^ l k
+
+        Rewrite each one-body term as an even weighting of all four 2-RDM
+        elements with delta functions. Then rearrange terms so that each ijkl
+        term gets a sum of permuted one-body terms multiplied by delta functions
+        One should notice that this results in the same formula if one was to
+        apply the wedge product!
+
+    Args:
+        molecular_hamiltonian: operator to write reduced hamiltonian for
+        n_electrons: number of electrons in the system
+    Returns:
+        InteractionOperator with a zero one-body component.
+    """
+    constant = molecular_hamiltonian.constant
+    h1 = molecular_hamiltonian.one_body_tensor
+    h2 = molecular_hamiltonian.two_body_tensor
+
+    delta = np.eye(h1.shape[0])
+    k2 = np.zeros_like(h2)
+    normalization = 1 / (4 * (n_electrons - 1))
+    for i, j, k, l in product(range(h1.shape[0]), repeat=4):
+        k2[i, j, k, l] = normalization * (
+            h1[i, l] * delta[j, k] + h1[j, k] * delta[i, l] -
+            h1[i, k] * delta[j, l] - h1[j, l] * delta[i, k]) + h2[i, j, k, l]
+
+    return InteractionOperator(constant, np.zeros_like(h1), k2)
\ No newline at end of file


Test patch:
diff --git a/src/openfermion/utils/_reduced_hamiltonian_test.py b/src/openfermion/utils/_reduced_hamiltonian_test.py
new file mode 100644
index 0000000..1cf802b
--- /dev/null
+++ b/src/openfermion/utils/_reduced_hamiltonian_test.py
@@ -0,0 +1,54 @@
+import os
+import numpy as np
+import openfermion as of
+from openfermion.config import DATA_DIRECTORY
+from openfermion.hamiltonians import MolecularData
+from openfermion.ops import InteractionOperator
+from openfermion.utils._reduced_hamiltonian import make_reduced_hamiltonian
+
+
+def test_mrd_return_type():
+    filename = os.path.join(DATA_DIRECTORY, "H2_sto-3g_singlet_0.7414.hdf5")
+    molecule = MolecularData(filename=filename)
+    reduced_ham = make_reduced_hamiltonian(molecule.get_molecular_hamiltonian(),
+                                           molecule.n_electrons)
+
+    assert isinstance(reduced_ham, InteractionOperator)
+
+
+def test_constant_one_body():
+    filename = os.path.join(DATA_DIRECTORY, "H2_sto-3g_singlet_0.7414.hdf5")
+    molecule = MolecularData(filename=filename)
+    reduced_ham = make_reduced_hamiltonian(molecule.get_molecular_hamiltonian(),
+                                           molecule.n_electrons)
+
+    assert np.isclose(reduced_ham.constant, molecule.nuclear_repulsion)
+    assert np.allclose(reduced_ham.one_body_tensor, 0)
+
+
+def test_fci_energy():
+    filename = os.path.join(DATA_DIRECTORY, "H2_sto-3g_singlet_0.7414.hdf5")
+    molecule = MolecularData(filename=filename)
+    reduced_ham = make_reduced_hamiltonian(molecule.get_molecular_hamiltonian(),
+                                           molecule.n_electrons)
+    np_ham = of.get_number_preserving_sparse_operator(
+        of.get_fermion_operator(reduced_ham),
+        molecule.n_qubits,
+        num_electrons=molecule.n_electrons,
+        spin_preserving=True)
+
+    w, _ = np.linalg.eigh(np_ham.toarray())
+    assert np.isclose(molecule.fci_energy, w[0])
+
+    filename = os.path.join(DATA_DIRECTORY, "H1-Li1_sto-3g_singlet_1.45.hdf5")
+    molecule = MolecularData(filename=filename)
+    reduced_ham = make_reduced_hamiltonian(molecule.get_molecular_hamiltonian(),
+                                           molecule.n_electrons)
+    np_ham = of.get_number_preserving_sparse_operator(
+        of.get_fermion_operator(reduced_ham),
+        molecule.n_qubits,
+        num_electrons=molecule.n_electrons,
+        spin_preserving=True)
+
+    w, _ = np.linalg.eigh(np_ham.toarray())
+    assert np.isclose(molecule.fci_energy, w[0])