Repository: quantumlib/OpenFermion
Pull request #89 (2017-10-28 16:34:31)

Data structure for quadratic Hamiltonians
Subclass of PolynomialTensor specifically for quadratic Hamiltonians (e.g. mean field models), both particle-number-conserving and not. I considered the use of the Majorana basis but decided that I really only need to use it internally. So this class uses the usual fermionic ladder operators for representation.

This is pretty basic so far, but I plan to add routines for conversion to and from FermionOperator, as well as methods to obtain circuits for ground state preparation. Just thought I should pause and ask for feedback before continuing.

googlebot: So there's good news and bad news.

:thumbsup: The good news is that everyone that needs to sign a CLA (the pull request submitter and all commit authors) have done so.  Everything is all good there.

:confused: The bad news is that it appears that one or more commits were authored by someone other than the pull request submitter.  We need to confirm that they're okay with their commits being contributed to this project.  Please have them confirm that here in the pull request.

*Note to project maintainer: This is a terminal state, meaning the `cla/google` commit status will not change from this State. It's up to you to confirm consent of the commit author(s) and merge this pull request when appropriate.*

<!-- need_author_consent -->
babbush: Okay, I'll merge this one then.

Patch:
diff --git a/src/openfermion/ops/__init__.py b/src/openfermion/ops/__init__.py
index 6b1b1ee..5500e88 100644
--- a/src/openfermion/ops/__init__.py
+++ b/src/openfermion/ops/__init__.py
@@ -20,3 +20,4 @@ from ._polynomial_tensor import (PolynomialTensor,
                                  two_body_basis_change)
 from ._interaction_operator import InteractionOperator
 from ._interaction_rdm import InteractionRDM
+from ._quadratic_hamiltonian import QuadraticHamiltonian
diff --git a/src/openfermion/ops/_quadratic_hamiltonian.py b/src/openfermion/ops/_quadratic_hamiltonian.py
new file mode 100644
index 0000000..285e9ab
--- /dev/null
+++ b/src/openfermion/ops/_quadratic_hamiltonian.py
@@ -0,0 +1,104 @@
+#   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.
+
+"""Class and functions to store and manipulate Hamiltonians that are quadratic
+in the fermionic ladder operators."""
+from __future__ import absolute_import
+
+import numpy
+
+from openfermion.config import EQ_TOLERANCE
+from openfermion.ops import PolynomialTensor
+
+
+class QuadraticHamiltonianError(Exception):
+    pass
+
+
+class QuadraticHamiltonian(PolynomialTensor):
+    """Class for storing Hamiltonians that are quadratic in the fermionic
+    ladder operators. The operators stored in this class take the form::
+
+        \sum_{p, q} (M_{pq} - \mu \delta_{pq}) a^\dagger_p a_q
+        + 1 / 2 \sum_{p, q} (A_{pq} a^\dagger_p a^\dagger_q + h.c.)
+        + constant
+
+    where
+        M is a Hermitian n_qubits x n_qubits matrix.
+        A is an antisymmetric n_qubits x n_qubits matrix.
+        \mu is a float representing the chemical potential term.
+        \delta_{pq} is the Kronecker delta symbol.
+
+    We separate the chemical potential \mu from M so that we can use it
+    to adjust the expectation value of the total number of particles.
+
+    Attributes:
+        n_qubits: An int giving the number of qubits.
+        constant: A constant term in the operator given as a float.
+    """
+
+    def __init__(self, constant, hermitian_part,
+                 antisymmetric_part=None, chemical_potential=None):
+        """
+        Initialize the QuadraticHamiltonian class.
+
+        Args:
+            constant(float): A constant term in the operator.
+            hermitian_part(ndarray): The matrix M, which represents the
+                coefficients of the particle-number-conserving terms.
+                This is an n_qubits x n_qubits numpy array of complex numbers.
+            antisymmetric_part(ndarray): The matrix A, which represents the
+                coefficients of the non-particle-number-conserving terms.
+                This is an n_qubits x n_qubits numpy array of complex numbers.
+            chemical_potential(float): The chemical potential \mu.
+        """
+        n_qubits = hermitian_part.shape[0]
+
+        # Initialize combined Hermitian part
+        if not chemical_potential:
+            combined_hermitian_part = hermitian_part
+        else:
+            combined_hermitian_part = (
+                    hermitian_part -
+                    chemical_potential * numpy.eye(n_qubits))
+
+        # Initialize antisymmetric part
+        if antisymmetric_part is None:
+            antisymmetric_part = numpy.zeros((n_qubits, n_qubits), complex)
+
+        super(QuadraticHamiltonian, self).__init__(
+                constant,
+                {(1, 0): combined_hermitian_part,
+                 (1, 1): .5 * antisymmetric_part,
+                 (0, 0): -.5 * antisymmetric_part.conj()})
+        self.chemical_potential = chemical_potential
+
+    def combined_hermitian_part(self):
+        """Return the Hermitian part including the chemical potential."""
+        return self.n_body_tensors[1, 0]
+
+    def hermitian_part(self):
+        """Return the Hermitian part not including the chemical potential."""
+        hermitian_part = self.combined_hermitian_part()
+        if self.chemical_potential:
+            hermitian_part += (self.chemical_potential *
+                               numpy.eye(self.n_qubits))
+        return hermitian_part
+
+    def antisymmetric_part(self):
+        """Return the antisymmetric part."""
+        return 2. * self.n_body_tensors[1, 1]
+
+    def conserves_particle_number(self):
+        """Return whether this Hamiltonian conserves particle number."""
+        discrepancy = numpy.max(numpy.abs(self.antisymmetric_part()))
+        return discrepancy < EQ_TOLERANCE
diff --git a/src/openfermion/utils/_operator_utils.py b/src/openfermion/utils/_operator_utils.py
index 1af9fdf..beeff34 100644
--- a/src/openfermion/utils/_operator_utils.py
+++ b/src/openfermion/utils/_operator_utils.py
@@ -35,8 +35,7 @@ def count_qubits(operator):
     """Compute the minimum number of qubits on which operator acts.
 
     Args:
-        operator: QubitOperator, InteractionOperator, FermionOperator,
-            PolynomialTensor, or InteractionRDM.
+        operator: FermionOperator, QubitOperator, or PolynomialTensor.
 
     Returns:
         num_qubits (int): The minimum number of qubits on which operator acts.
@@ -62,10 +61,8 @@ def count_qubits(operator):
                     num_qubits = term[-1][0] + 1
         return num_qubits
 
-    # Handle InteractionOperator, InteractionRDM, PolynomialTensor.
-    elif isinstance(operator, (InteractionOperator,
-                               InteractionRDM,
-                               PolynomialTensor)):
+    # Handle PolynomialTensor
+    elif isinstance(operator, PolynomialTensor):
         return operator.n_qubits
 
     # Raise for other classes.


Test patch:
diff --git a/src/openfermion/ops/_interaction_operator_test.py b/src/openfermion/ops/_interaction_operator_test.py
index 15f8025..3e3016a 100644
--- a/src/openfermion/ops/_interaction_operator_test.py
+++ b/src/openfermion/ops/_interaction_operator_test.py
@@ -10,7 +10,7 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 
-"""Tests for interaction_operators.py."""
+"""Tests for interaction_operator.py."""
 from __future__ import absolute_import
 
 import numpy
@@ -19,7 +19,7 @@ import unittest
 from openfermion.ops import InteractionOperator
 
 
-class InteractionOperatorsTest(unittest.TestCase):
+class InteractionOperatorTest(unittest.TestCase):
 
     def setUp(self):
         self.n_qubits = 5
diff --git a/src/openfermion/ops/_quadratic_hamiltonian_test.py b/src/openfermion/ops/_quadratic_hamiltonian_test.py
new file mode 100644
index 0000000..e852031
--- /dev/null
+++ b/src/openfermion/ops/_quadratic_hamiltonian_test.py
@@ -0,0 +1,84 @@
+#   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 quadratic_hamiltonian.py."""
+from __future__ import absolute_import
+
+import numpy
+import unittest
+
+from openfermion.ops import QuadraticHamiltonian
+
+
+class QuadraticHamiltoniansTest(unittest.TestCase):
+
+    def setUp(self):
+        self.n_qubits = 5
+        self.constant = 1.
+        self.chemical_potential = 2.
+
+        # Obtain random Hermitian and antisymmetric matrices
+        rand_mat = numpy.random.randn(self.n_qubits, self.n_qubits)
+        self.hermitian_mat = rand_mat + rand_mat.T.conj()
+        rand_mat = numpy.random.randn(self.n_qubits, self.n_qubits)
+        self.antisymmetric_mat = rand_mat - rand_mat.T
+
+        self.combined_hermitian = (
+                self.hermitian_mat -
+                self.chemical_potential * numpy.eye(self.n_qubits))
+
+        # Initialize a particle-number-conserving Hamiltonian
+        self.quad_ham_pc = QuadraticHamiltonian(
+                self.constant, self.hermitian_mat)
+
+        # Initialize a non-particle-number-conserving Hamiltonian
+        self.quad_ham_npc = QuadraticHamiltonian(
+                self.constant, self.hermitian_mat, self.antisymmetric_mat,
+                self.chemical_potential)
+
+    def test_combined_hermitian_part(self):
+        """Test getting the combined Hermitian part."""
+        combined_hermitian_part = self.quad_ham_pc.combined_hermitian_part()
+        for i in numpy.ndindex(combined_hermitian_part.shape):
+            self.assertAlmostEqual(self.hermitian_mat[i],
+                                   combined_hermitian_part[i])
+
+        combined_hermitian_part = self.quad_ham_npc.combined_hermitian_part()
+        for i in numpy.ndindex(combined_hermitian_part.shape):
+            self.assertAlmostEqual(self.combined_hermitian[i],
+                                   combined_hermitian_part[i])
+
+    def test_hermitian_part(self):
+        """Test getting the Hermitian part."""
+        hermitian_part = self.quad_ham_pc.hermitian_part()
+        for i in numpy.ndindex(hermitian_part.shape):
+            self.assertAlmostEqual(self.hermitian_mat[i], hermitian_part[i])
+
+        hermitian_part = self.quad_ham_npc.hermitian_part()
+        for i in numpy.ndindex(hermitian_part.shape):
+            self.assertAlmostEqual(self.hermitian_mat[i], hermitian_part[i])
+
+    def test_antisymmetric_part(self):
+        """Test getting the antisymmetric part."""
+        antisymmetric_part = self.quad_ham_pc.antisymmetric_part()
+        for i in numpy.ndindex(antisymmetric_part.shape):
+            self.assertAlmostEqual(0., antisymmetric_part[i])
+
+        antisymmetric_part = self.quad_ham_npc.antisymmetric_part()
+        for i in numpy.ndindex(antisymmetric_part.shape):
+            self.assertAlmostEqual(self.antisymmetric_mat[i],
+                                   antisymmetric_part[i])
+
+    def test_conserves_particle_number(self):
+        """Test checking whether Hamiltonian conserves particle number."""
+        self.assertTrue(self.quad_ham_pc.conserves_particle_number())
+        self.assertFalse(self.quad_ham_npc.conserves_particle_number())