Added functionality to compute one-norm of qubit/majorana Hamiltonian from molecular integrals

src/openfermion/__init__.py

@@ -34,7 +34,7 @@
     make_reduced_hamiltonian,
 )
 
-from openfermion.functionals import contextuality
+from openfermion.functionals import (contextuality, get_one_norm)
 
 from openfermion.hamiltonians import (
     FermiHubbardModel,

src/openfermion/functionals/__init__.py

@@ -11,3 +11,5 @@
 #   limitations under the License.
 
 from .contextuality import is_contextual
+
+from .get_one_norm import get_one_norm

src/openfermion/functionals/get_one_norm.py

@@ -0,0 +1,130 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""
+Function to calculate the 1-Norm of a molecular Hamiltonian after
+fermion-to-qubit transformation from a MolecularData class.
+See https://arxiv.org/abs/2103.14753
+"""
+
+import numpy as np
+
+from openfermion import MolecularData
+
+
+def get_one_norm(mol_or_int, return_constant=True):
+    r"""
+
+    Returns the 1-Norm of the Hamiltonian described in
+    https://arxiv.org/abs/2103.14753 after a fermion-to-qubit
+    transformation given nuclear constant, one-body (2D np.array)
+    and two-body (4D np.array) integrals.
+
+    Parameters
+    ----------
+
+    mol_or_int : Tuple of (constant, one_body_integrals, two_body_integrals)
+    constant : Nuclear repulsion or adjustment to constant shift in Hamiltonian
+                from integrating out core orbitals
+    one_body_integrals : An array of the one-electron integrals having
+                shape of (n_orb, n_orb).
+    two_body_integrals : An array of the two-electron integrals having
+                shape of (n_orb, n_orb, n_orb, n_orb).
+
+    -----OR----
+
+    mol_or_int : MolecularData class object
+
+    -----------
+
+    return_constant (optional) : If False, do not return the constant term in
+                in the (majorana/qubit) Hamiltonian
+
+    Returns
+    -------
+    one_norm : 1-Norm of the Qubit Hamiltonian
+    """
+    if isinstance(mol_or_int, MolecularData):
+        return _get_one_norm_mol(mol_or_int, return_constant)
+    else:
+        if return_constant:
+            constant, one_body_integrals, two_body_integrals = mol_or_int
+            return _get_one_norm(constant, one_body_integrals,
+                                 two_body_integrals)
+        else:
+            _, one_body_integrals, two_body_integrals = mol_or_int
+            return _get_one_norm_woconst(one_body_integrals, two_body_integrals)
+
+
+def _get_one_norm_mol(molecule, return_constant):
+    if return_constant:
+        return _get_one_norm(molecule.nuclear_repulsion,
+                             molecule.one_body_integrals,
+                             molecule.two_body_integrals)
+    else:
+        return _get_one_norm_woconst(molecule.one_body_integrals,
+                                     molecule.two_body_integrals)
+
+
+def _get_one_norm(constant, one_body_integrals, two_body_integrals):
+    n_orb = one_body_integrals.shape[0]
+
+    htilde = constant
+    for p in range(n_orb):
+        htilde += one_body_integrals[p, p]
+        for q in range(n_orb):
+            htilde += ((1 / 2 * two_body_integrals[p, q, q, p]) -
+                       (1 / 4 * two_body_integrals[p, q, p, q]))
+
+    htildepq = np.zeros(one_body_integrals.shape)
+    for p in range(n_orb):
+        for q in range(n_orb):
+            htildepq[p, q] = one_body_integrals[p, q]
+            for r in range(n_orb):
+                htildepq[p, q] += ((two_body_integrals[p, r, r, q]) -
+                                   (1 / 2 * two_body_integrals[p, r, q, r]))
+
+    one_norm = abs(htilde) + np.sum(np.absolute(htildepq))
+
+    for p in range(n_orb):
+        for q in range(n_orb):
+            for r in range(n_orb):
+                for s in range(n_orb):
+                    if p > q and r > s:
+                        one_norm += 1 / 2 * abs(two_body_integrals[p, q, r, s] -
+                                                two_body_integrals[p, q, s, r])
+                    one_norm += 1 / 4 * abs(two_body_integrals[p, q, r, s])
+
+    return one_norm
+
+
+def _get_one_norm_woconst(one_body_integrals, two_body_integrals):
+    n_orb = one_body_integrals.shape[0]
+
+    htildepq = np.zeros(one_body_integrals.shape)
+    for p in range(n_orb):
+        for q in range(n_orb):
+            htildepq[p, q] = one_body_integrals[p, q]
+            for r in range(n_orb):
+                htildepq[p, q] += ((two_body_integrals[p, r, r, q]) -
+                                   (1 / 2 * two_body_integrals[p, r, q, r]))
+
+    one_norm = np.sum(np.absolute(htildepq))
+
+    for p in range(n_orb):
+        for q in range(n_orb):
+            for r in range(n_orb):
+                for s in range(n_orb):
+                    if p > q and r > s:
+                        one_norm += 1 / 2 * abs(two_body_integrals[p, q, r, s] -
+                                                two_body_integrals[p, q, s, r])
+                    one_norm += 1 / 4 * abs(two_body_integrals[p, q, r, s])
+    return one_norm

src/openfermion/functionals/get_one_norm_test.py

@@ -0,0 +1,44 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""Tests for get_one_norm."""
+
+import os
+import pytest
+from openfermion import get_one_norm, MolecularData, jordan_wigner
+from openfermion.config import DATA_DIRECTORY
+
+filename = os.path.join(DATA_DIRECTORY, "H1-Li1_sto-3g_singlet_1.45.hdf5")
+molecule = MolecularData(filename=filename)
+molecular_hamiltonian = molecule.get_molecular_hamiltonian()
+qubit_hamiltonian = jordan_wigner(molecular_hamiltonian)
+
+
+def test_one_norm_from_molecule():
+    assert qubit_hamiltonian.induced_norm() == pytest.approx(
+        get_one_norm(molecule))
+
+
+def test_one_norm_from_ints():
+    ints = (molecule.nuclear_repulsion, molecule.one_body_integrals,
+            molecule.two_body_integrals)
+    assert qubit_hamiltonian.induced_norm() == pytest.approx(get_one_norm(ints))
+
+
+def test_one_norm_woconst():
+    one_norm_woconst = (qubit_hamiltonian.induced_norm() -
+                        abs(qubit_hamiltonian.constant))
+    ints = (molecule.nuclear_repulsion, molecule.one_body_integrals,
+            molecule.two_body_integrals)
+    assert one_norm_woconst == pytest.approx(
+        get_one_norm(molecule, return_constant=False))
+    assert one_norm_woconst == pytest.approx(
+        get_one_norm(ints, return_constant=False))