Repository: quantumlib/OpenFermion
Pull request #704 (2021-01-19 21:00:35)

Doci hamiltonian

google-cla[bot]: We found a Contributor License Agreement for you (the sender of this pull request), but were unable to find agreements for all the commit author(s) or Co-authors.  If you authored these, maybe you used a different email address in the git commits than was used to sign the CLA ([login here](https://cla.developers.google.com/) to double check)?  If these were authored by someone else, then they will need to sign a CLA as well, and confirm that they're okay with these being contributed to Google.
In order to pass this check, please resolve this problem and then comment `@googlebot I fixed it.`. If the bot doesn't comment, it means it doesn't think anything has changed.

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F704) for more info**.

<!-- need_author_cla -->
ncrubin: @obriente I merged master and fixed some commit issues but it seems like this branch somehow has your very long commit history in it.  Is it possible to compress like you did with the other PRs?
google-cla[bot]: We found a Contributor License Agreement for you (the sender of this pull request), but were unable to find agreements for all the commit author(s) or Co-authors.  If you authored these, maybe you used a different email address in the git commits than was used to sign the CLA ([login here](https://cla.developers.google.com/) to double check)?  If these were authored by someone else, then they will need to sign a CLA as well, and confirm that they're okay with these being contributed to Google.
In order to pass this check, please resolve this problem and then comment `@googlebot I fixed it.`. If the bot doesn't comment, it means it doesn't think anything has changed.

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F704) for more info**.

<!-- need_author_cla -->
obriente: Hi Nick, taking a look through the commit history I think these are all relevant to this PR. (I think I've fixed the issues on my personal OF fork now.) Probably easiest to just leave this as is?
ncrubin: @googlebot I fixed it.
google-cla[bot]: We found a Contributor License Agreement for you (the sender of this pull request), but were unable to find agreements for all the commit author(s) or Co-authors.  If you authored these, maybe you used a different email address in the git commits than was used to sign the CLA ([login here](https://cla.developers.google.com/) to double check)?  If these were authored by someone else, then they will need to sign a CLA as well, and confirm that they're okay with these being contributed to Google.
In order to pass this check, please resolve this problem and then comment `@googlebot I fixed it.`. If the bot doesn't comment, it means it doesn't think anything has changed.

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F704) for more info**.

<!-- need_author_cla -->
ncrubin: @cvjjm Can you consent?
cvjjm: I have the OK from Legal, but it was decided that we need to sign a company wide CLA and for that we need a caretaker and documentation.... The wheels are turning but it will take a few more days.
googlebot: We found a Contributor License Agreement for you (the sender of this pull request), but were unable to find agreements for all the commit author(s) or Co-authors.  If you authored these, maybe you used a different email address in the git commits than was used to sign the CLA ([login here](https://cla.developers.google.com/) to double check)?  If these were authored by someone else, then they will need to sign a CLA as well, and confirm that they're okay with these being contributed to Google.
In order to pass this check, please resolve this problem and then comment `@googlebot I fixed it.`. If the bot doesn't comment, it means it doesn't think anything has changed.

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F704) for more info**.

<!-- need_author_cla -->
cvjjm: Changes look good to me. Unfortunately this PR is now also stuck because of the CLA. I already requested a CLA via the link https://cla.developers.google.com/ last Friday, but we haven't received anything to sign since. Have you had this happen previously? Unfortunately there is no contact on that page where I could ask what the status is.
googlebot: All (the pull request submitter and all commit authors) CLAs are signed, **but** one or more commits were authored or co-authored by someone other than the pull request submitter.

We need to confirm that all authors are ok with their commits being contributed to this project.  Please have them confirm that by leaving a comment that contains only `@googlebot I consent.` in this pull request.

*Note to project maintainer:* There may be cases where the author cannot leave a comment, or the comment is not properly detected as consent.  In those cases, you can manually confirm consent of the commit author(s), and set the `cla` label to `yes` (if enabled on your project).

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F704) for more info**.

<!-- need_author_consent -->
ncrubin: @cvjjm It looks like you can now consent.  Covestro is showing up in the CLA list.
cvjjm: @googlebot I consent.
google-cla[bot]: CLAs look good, thanks!

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F704) for more info**.

<!-- ok -->
cvjjm: Looks like we are finally ready to go!
google-cla[bot]: CLAs look good, thanks!

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F704) for more info**.

<!-- ok -->
ncrubin: @cvjjm Oh no! tests failed when I updated the branch.  This is also happening for on of my PRs where I changed 1 line.  I will get to the bottom of this.  Something else is broken...not this PR.
ncrubin: In fact I think it is related to your _issmall PR.  Was that passing tests?
ncrubin: @cvjjm okay it is not your PR.  The "bug" was introduced by the the 1.6.1 release of scipy (which was released 2 days ago). There is a change that fixes a `dtype` propagation in initializing sparse matrices. The tests that are failing all use sparse_tools.py methods that are initializing a sparse matrix with `csc_matrix`.   Scipy's PR# 13403 fixes a bug with sparse matrix initialization.  When tuples are passed to csc_matrix (which is a subclass of `_cs_matrix`) they eventually call coo_matrix for initialization.  coo_matrix WAS using `np.array(data)` to cast the `data` for the `coo_matrix`.  np.array is very clever sometimes.  If you provide `data` as a list of `float`s and `np.complex128`s it is smart enough to cast the resulting array as the higher dtype--`np.complex128`.  This is bad because it was ignoring the dtype that was passed to the `coo_matrix` constructor.  Thus you are at the whim of however smart np.array is at determining dtypes.  For us this was a blessing.  We were getting defaulted into the right 'np.complex128' without knowing it.  We weren't specifying the dtype so it defaulted to `float` and then at some point was changed to `np.complex128` when `np.array`'s cleverness kicked in. When the new scipy release got cut this bug was fixed and all of a sudden our tests that tried to initialize complex sparse matrices started throwing errors about casting to floats.  We were trying to pass data with complex numbers to `coo_matrix` but without specifying the dtype.   Thus scipy caught our bug of not correctly specifying the datatype.    

More info can be found here (https://github.com/scipy/scipy/pull/13403).

I will open a PR to fix the sparse routines and then circle back and merge this PR.
google-cla[bot]: CLAs look good, thanks!

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F704) for more info**.

<!-- ok -->
google-cla[bot]: CLAs look good, thanks!

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F704) for more info**.

<!-- ok -->

Patch:
diff --git a/src/openfermion/__init__.py b/src/openfermion/__init__.py
index 952902c..1529b5b 100644
--- a/src/openfermion/__init__.py
+++ b/src/openfermion/__init__.py
@@ -186,6 +186,7 @@ from openfermion.ops import (
     InteractionRDMError,
     QuadraticHamiltonian,
     QuadraticHamiltonianError,
+    DOCIHamiltonian,
     antisymmetric_canonical_form,
 )
 
diff --git a/src/openfermion/chem/molecular_data.py b/src/openfermion/chem/molecular_data.py
index e51321e..c0f0d1f 100644
--- a/src/openfermion/chem/molecular_data.py
+++ b/src/openfermion/chem/molecular_data.py
@@ -871,30 +871,9 @@ class MolecularData(object):
 
         # Get integrals.
         one_body_integrals, two_body_integrals = self.get_integrals()
-
-        # Determine core constant
-        core_constant = 0.0
-        for i in occupied_indices:
-            core_constant += 2 * one_body_integrals[i, i]
-            for j in occupied_indices:
-                core_constant += (2 * two_body_integrals[i, j, j, i] -
-                                  two_body_integrals[i, j, i, j])
-
-        # Modified one electron integrals
-        one_body_integrals_new = numpy.copy(one_body_integrals)
-        for u in active_indices:
-            for v in active_indices:
-                for i in occupied_indices:
-                    one_body_integrals_new[u, v] += (
-                        2 * two_body_integrals[i, u, v, i] -
-                        two_body_integrals[i, u, i, v])
-
-        # Restrict integral ranges and change M appropriately
-        return (core_constant,
-                one_body_integrals_new[numpy.ix_(active_indices,
-                                                 active_indices)],
-                two_body_integrals[numpy.ix_(active_indices, active_indices,
-                                             active_indices, active_indices)])
+        return reps.get_active_space_integrals(one_body_integrals,
+                                               two_body_integrals,
+                                               occupied_indices, active_indices)
 
     def get_molecular_hamiltonian(self,
                                   occupied_indices=None,
@@ -953,17 +932,15 @@ class MolecularData(object):
                 raise MissingCalculationError(
                     'Missing FCI RDM in {}'.format(self.filename) +
                     'Run FCI calculation before loading FCI RDMs.')
-            else:
-                one_rdm = self.fci_one_rdm
-                two_rdm = self.fci_two_rdm
+            one_rdm = self.fci_one_rdm
+            two_rdm = self.fci_two_rdm
         else:
             if self.cisd_energy is None:
                 raise MissingCalculationError(
                     'Missing CISD RDM in {}'.format(self.filename) +
                     'Run CISD calculation before loading CISD RDMs.')
-            else:
-                one_rdm = self.cisd_one_rdm
-                two_rdm = self.cisd_two_rdm
+            one_rdm = self.cisd_one_rdm
+            two_rdm = self.cisd_two_rdm
 
         # Truncate.
         one_rdm[numpy.absolute(one_rdm) < EQ_TOLERANCE] = 0.
diff --git a/src/openfermion/ops/__init__.py b/src/openfermion/ops/__init__.py
index 59250cc..a844b63 100644
--- a/src/openfermion/ops/__init__.py
+++ b/src/openfermion/ops/__init__.py
@@ -43,5 +43,6 @@ from .representations import (
     InteractionRDMError,
     QuadraticHamiltonian,
     QuadraticHamiltonianError,
+    DOCIHamiltonian,
     antisymmetric_canonical_form,
 )
diff --git a/src/openfermion/ops/representations/__init__.py b/src/openfermion/ops/representations/__init__.py
index 21785ea..37ef64d 100644
--- a/src/openfermion/ops/representations/__init__.py
+++ b/src/openfermion/ops/representations/__init__.py
@@ -9,6 +9,8 @@ from .diagonal_coulomb_hamiltonian import DiagonalCoulombHamiltonian
 from .interaction_operator import (
     InteractionOperator,
     InteractionOperatorError,
+    get_tensors_from_integrals,
+    get_active_space_integrals,
 )
 
 from .interaction_rdm import (
@@ -21,3 +23,4 @@ from .quadratic_hamiltonian import (
     QuadraticHamiltonianError,
     antisymmetric_canonical_form,
 )
+from .doci_hamiltonian import DOCIHamiltonian
diff --git a/src/openfermion/ops/representations/doci_hamiltonian.py b/src/openfermion/ops/representations/doci_hamiltonian.py
new file mode 100644
index 0000000..12fcd03
--- /dev/null
+++ b/src/openfermion/ops/representations/doci_hamiltonian.py
@@ -0,0 +1,456 @@
+#   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 interaction operators."""
+import numpy
+
+from openfermion.ops import QubitOperator
+from openfermion.ops.representations import (PolynomialTensor,
+                                             get_tensors_from_integrals)
+
+COEFFICIENT_TYPES = (int, float, complex)
+
+
+class DOCIHamiltonian(PolynomialTensor):
+    r"""Class for storing DOCI hamiltonians which are defined to be
+    restrictions of fermionic operators to doubly occupied configurations.
+    As such they are by nature hard-core boson Hamiltonians, but the
+    hard-core Boson algebra is identical to the Pauli algebra, which is
+    why it is convenient to represent DOCI hamiltonians as QubitOperators
+
+    Note that the operators stored in this class take the form:
+
+        .. math::
+
+            constant + \sum_{p} h^{(r1)}_{p, p}/2 (1 - \sigma^Z_p) +
+            \sum_{p \neq q} h^{(r1)}_{p, q}/4
+                * (\sigma^X_p \sigma^X_q + \sigma^Y_p \sigma^Y_q) +
+            \sum_{p \neq q} h^{(r2)}_{p, q}/4 (1 - \sigma^Z_p -
+            \sigma^Z_q + \sigma^Z_p \sigma^Z_q)
+            =
+            constant + \sum_{p} h_{p, p} N_p +
+            \sum_{p \neq q} w_{p, p} N_p N_q +
+            \sum_{p \neq q} v_{p, p} P_p^\dagger P_q,
+
+    where
+
+        .. math::
+
+            N_p = (1 - \sigma^Z_p)/2,
+            P_p = a_{i,\beta} a_{i,\alpha},
+            h_p = h^{(r1)}_{p, p} = \langle p|h|p \rangle =
+                2 I^{(1)}_{p, p} + I^{(2)}_{p, p, p, p},
+            w_{p, q} = h^{(r2)}_{p, q} = 2 \langle pq|v|pq \rangle -
+                                         \langle pq|v|qp \rangle =
+                                        2 I^{(2)}_{p, q, q, p} -
+                                        I^{(2)}_{p, q, p, q},
+            v_{p, q} = h^{(r1)}_{p, q} = \langle pp|v|qq \rangle =
+                I^{(2)}_{p, p, q, q},
+
+    with (:math:`I^{(1)}_{p, q}`) and (:math:`I^{(2)}_{p, q, r, s}`) are the one
+    and two body electron integrals and (:math:`h`) and (:math:`v`) are the
+    coefficients of the corresponding InteractionOperator
+
+        .. math::
+
+            constant + \sum_{p, q} h_{p, q} a^\dagger_p a_q +
+            \sum_{p, q, r, s} h_{p, q, r, s} a^\dagger_p a^\dagger_q a_r a_s.
+
+
+    Attributes:
+        constant: The constant offset.
+        hr1: The coefficients of (:math:`h^{r1}_{p, q}`).
+            This is an n_qubits x n_qubits numpy array of floats.
+        hr2: The coefficients of (:math:`h^{r2}_{p, q}`).
+            This is an n_qubits x n_qubits numpy array of floats.
+    """
+
+    def __init__(self, constant, hc, hr1, hr2):
+        r"""
+        Initialize the DOCIHamiltonian class.
+
+        Args:
+            constant: A constant term in the operator given as a
+                float. For instance, the nuclear repulsion energy.
+            hc: the coefficients of (:math:`h^{(c)}_{p}`)
+            hr1: The coefficients of (:math:`h^{(r1)}_{p, q}`).
+               This is an n_qubits x n_qubits numpy array of floats.
+            hr2: The coefficients of (:math:`h^{(r2)}_{p, q}`).
+                This is an n_qubits x n_qubits array of floats.
+        """
+        super(DOCIHamiltonian, self).__init__(None)
+
+        self._n_qubits = hc.shape[0]
+        self._constant = constant
+        self._hr1 = hr1
+        self._hr2 = hr2
+        self._hc = hc
+
+    @property
+    def qubit_operator(self):
+        """Return the QubitOperator representation of this DOCI Hamiltonian"""
+        return self.identity_part + self.z_part + self.xy_part
+
+    def xx_term(self, p, q):
+        """Returns the XX term on a single pair of qubits as a QubitOperator
+        Arguments:
+            p, q [int] -- qubit indices
+        Returns:
+            [QubitOperator] -- XX term on the chosen qubits.
+        """
+        return QubitOperator("X" + str(p) + " X" + str(q), self.hr1[p, q] / 2)
+
+    def yy_term(self, p, q):
+        """Returns the YY term on a single pair of qubits as a QubitOperator
+        Arguments:
+            p, q [int] -- qubit indices
+        Returns:
+            [QubitOperator] -- YY term on the chosen qubits.
+        """
+        return QubitOperator("Y" + str(p) + " Y" + str(q), self.hr1[p, q] / 2)
+
+    def z_term(self, p):
+        """Returns the Z term on a single qubit as a QubitOperator
+        Arguments:
+            p [int] -- qubit index
+        Returns:
+            [QubitOperator] -- Z term on the chosen qubit.
+        """
+        return QubitOperator("Z" + str(p),
+                             self.hc[p] / 2 + sum(self.hr2[:, p]) / 2)
+
+    def zz_term(self, p, q):
+        """Returns the ZZ term on a single pair of qubits as a QubitOperator
+        Arguments:
+            p, q [int] -- qubit indices
+        Returns:
+            [QubitOperator] -- ZZ term on the chosen qubits.
+        """
+        return QubitOperator("Z" + str(p) + " Z" + str(q), self.hr2[p, q] / 2)
+
+    @property
+    def identity_part(self):
+        """Returns identity term of this operator (i.e. trace-ful term)
+        in QubitOperator form.
+        """
+        return QubitOperator(
+            (), self.constant + numpy.sum(self.hc) / 2 +
+            numpy.sum(self.hr2) / 4 + numpy.sum(numpy.diag(self.hr2)) / 4)
+
+    @property
+    def xx_part(self):
+        """Returns the XX part of the QubitOperator representation of this
+        DOCIHamiltonian
+        """
+        return sum([
+            self.xx_term(p, q)
+            for p in range(self.n_qubits)
+            for q in range(p + 1, self.n_qubits)
+        ])
+
+    @property
+    def yy_part(self):
+        """Returns the YY part of the QubitOperator representation of this
+        DOCIHamiltonian
+        """
+        return sum([
+            self.yy_term(p, q)
+            for p in range(self.n_qubits)
+            for q in range(p + 1, self.n_qubits)
+        ])
+
+    @property
+    def xy_part(self):
+        """Returns the XX+YY part of the QubitOperator representation of this
+        DOCIHamiltonian
+        """
+        return self.xx_part + self.yy_part
+
+    @property
+    def zz_part(self):
+        """Returns the ZZ part of the QubitOperator representation of this
+        DOCIHamiltonian
+        """
+        return sum([
+            self.zz_term(p, q)
+            for p in range(self.n_qubits)
+            for q in range(p + 1, self.n_qubits)
+        ])
+
+    @property
+    def z_part(self):
+        """Return the Z and ZZ part of the QubitOperator representation of this
+        DOCI Hamiltonian"""
+        return self.zz_part + sum(
+            [self.z_term(p) for p in range(self.n_qubits)])
+
+    @property
+    def hc(self):
+        return self._hc
+
+    @hc.setter
+    def hc(self, value):
+        self._hc = value
+
+    @property
+    def hr1(self):
+        """The value of hr1."""
+        return self._hr1
+
+    @hr1.setter
+    def hr1(self, value):
+        self._hr1 = value
+
+    @property
+    def hr2(self):
+        """The value of hr2."""
+        return self._hr2
+
+    @hr2.setter
+    def hr2(self, value):
+        self._hr2 = value
+
+    # Override base class
+    @property
+    def constant(self):
+        return self._constant
+
+    @constant.setter
+    def constant(self, value):
+        self._constant = value
+
+    # Override base class to generate on the fly
+    @property
+    def n_body_tensors(self):
+        one_body_coefficients, two_body_coefficients =\
+            get_tensors_from_doci(self.hc, self.hr1, self.hr2)
+        n_body_tensors = {
+            (): self.constant,
+            (1, 0): one_body_coefficients,
+            (1, 1, 0, 0): two_body_coefficients
+        }
+        return n_body_tensors
+
+    @n_body_tensors.setter
+    def n_body_tensors(self, value):
+        raise TypeError('Raw edits of the n_body_tensors of a DOCIHamiltonian '
+                        'is not allowed. Either adjust the hc/hr1/hr2 terms '
+                        'or cast to another PolynomialTensor class.')
+
+    def _get_onebody_term(self, key, index):
+        if key[0] != 1 or key[1] != 0:
+            raise IndexError('DOCIHamiltonian class only contains '
+                             'one-body terms in the (1, 0) sector.')
+        if index[0] != index[1]:
+            raise IndexError('DOCIHamiltonian class only contains '
+                             'diagonal one-body electron integrals.')
+        return self.hc[index[0] // 2] / 2
+
+    def _get_twobody_term(self, key, index):
+        if key[0] != 1 or key[1] != 1 or key[2] != 0 or key[3] != 0:
+            raise IndexError('DOCIHamiltonian class only contains '
+                             'two-body terms in the (1, 1, 0, 0) sector.')
+        if index[0] == index[3] and index[1] == index[2]:
+            return self.hr2[index[0] // 2, index[1] // 2] / 2
+        if index[0] // 2 == index[1] // 2 and index[2] // 2 == index[3] // 2:
+            return self.hr1[index[0] // 2, index[2] // 2] / 2
+        raise IndexError('DOCIHamiltonian class only contains '
+                         'two-electron integrals corresponding '
+                         'to a double excitation.')
+
+    # Override base class
+    def __getitem__(self, args):
+        """Look up matrix element.
+
+        Args:
+            args: Tuples indicating which coefficient to get. For instance,
+                `my_tensor[(6, 1), (8, 1), (2, 0)]`
+                returns
+                `my_tensor.n_body_tensors[1, 1, 0][6, 8, 2]`
+        """
+        if len(args) == 0:
+            return self.constant
+        index = tuple([operator[0] for operator in args])
+        key = tuple([operator[1] for operator in args])
+        if len(index) == 2:
+            return self._get_onebody_term(key, index)
+        if len(index) == 4:
+            return self._get_twobody_term(key, index)
+        raise IndexError('DOCIHamiltonian class only contains '
+                         'one and two-electron and constant terms.')
+
+    # Override root class
+    def __setitem__(self, args, value):
+        # This is not well-defined for a DOCIHamiltonian as we want to keep
+        # certain tensor elements within the class the same. Better to
+        # make the user update the hr1/hr2/hc terms --- if they really
+        # want to play around with the n_body_tensors here they should
+        # be casting this to a raw PolynomialTensor.
+        raise TypeError('Raw edits of the n_body_tensors of a DOCIHamiltonian '
+                        'is not allowed. Either adjust the hc/hr1/hr2 terms '
+                        'or cast to another PolynomialTensor class.')
+
+    # Override root class
+    def __iadd__(self, addend):
+        if isinstance(addend, COEFFICIENT_TYPES):
+            self.constant += addend
+            return self
+        if not issubclass(type(addend), DOCIHamiltonian):
+            raise TypeError('Invalid type.')
+        if self.n_qubits != addend.n_qubits:
+            raise TypeError('Invalid tensor shape.')
+        self.hc += addend.hc
+        self.hr1 += addend.hr1
+        self.hr2 += addend.hr2
+        self.constant += addend.constant
+        return self
+
+    def __isub__(self, subtrahend):
+        if isinstance(subtrahend, COEFFICIENT_TYPES):
+            self.constant -= subtrahend
+            return self
+        if not issubclass(type(subtrahend), DOCIHamiltonian):
+            raise TypeError('Invalid type.')
+        if self.n_qubits != subtrahend.n_qubits:
+            raise TypeError('Invalid tensor shape.')
+        self.hc -= subtrahend.hc
+        self.hr1 -= subtrahend.hr1
+        self.hr2 -= subtrahend.hr2
+        self.constant -= subtrahend.constant
+        return self
+
+    def __imul__(self, multiplier):
+        if not isinstance(multiplier, COEFFICIENT_TYPES):
+            raise TypeError('Invalid type.')
+        self.hc *= multiplier
+        self.hr1 *= multiplier
+        self.hr2 *= multiplier
+        self.constant *= multiplier
+        return self
+
+    def __itruediv__(self, dividend):
+        if not isinstance(dividend, COEFFICIENT_TYPES):
+            raise TypeError('Invalid type.')
+        self.hc /= dividend
+        self.hr1 /= dividend
+        self.hr2 /= dividend
+        self.constant /= dividend
+        return self
+
+    @classmethod
+    def from_integrals(cls, constant, one_body_integrals, two_body_integrals):
+        hc, hr1, hr2 = get_doci_from_integrals(one_body_integrals,
+                                               two_body_integrals)
+        return cls(constant, hc, hr1, hr2)
+
+    @classmethod
+    def zero(cls, n_qubits):
+        return cls(0, numpy.zeros((n_qubits,), dtype=numpy.complex128),
+                   numpy.zeros((n_qubits,) * 2, dtype=numpy.complex128),
+                   numpy.zeros((n_qubits,) * 2, dtype=numpy.complex128))
+
+
+def get_tensors_from_doci(hc, hr1, hr2):
+    '''Makes the one and two-body tensors from the DOCI Hamiltonian matrices
+
+    Args:
+        hc [numpy array]: The single-particle DOCI terms in matrix form
+        hr1 [numpy array]: The off-diagonal DOCI Hamiltonian terms in matrix
+            form
+        hr2 [numpy array]: The diagonal DOCI Hamiltonian terms in matrix form
+
+    Returns:
+        one_body_coefficients [numpy array]: The corresponding one-body
+            tensor for the electronic structure Hamiltonian
+        two_body_coefficients [numpy array]: The corresponding two body
+            tensor for the electronic structure Hamiltonian
+    '''
+    one_body_integrals, two_body_integrals =\
+        get_projected_integrals_from_doci(hc, hr1, hr2)
+    one_body_coefficients, two_body_coefficients = get_tensors_from_integrals(
+        one_body_integrals, two_body_integrals)
+    return one_body_coefficients, two_body_coefficients
+
+
+def get_projected_integrals_from_doci(hc, hr1, hr2):
+    '''Generates a set of atomic integrals corresponding to this Hamiltonian
+
+    Makes the one and two-body integrals from the DOCI projectionof the hr1,
+    hr2, and hc matrices. This is technically not well-defined, as hr2 is
+    not generated in a one-to-one fashion. In particular, here we assume
+    hr2 corresponds entirely to a pqqp term, while this could be equally well
+    added to the pqpq term. This implies that calling
+
+    get_doci_from_integrals(*get_projected_integrals_from_doci(hc, hr1, hr2))
+
+    should return the same hc, hr1, and hr2, but there is no such guarantee
+    for
+
+    get_projected_integrals_from_doci(*get_doci_from_integrals(
+        one_body_integrals, two_body_integrals))
+
+    Args:
+        hc [numpy array]: The single-particle DOCI terms in matrix form
+        hr1 [numpy array]: The off-diagonal DOCI Hamiltonian terms in matrix
+            form
+        hr2 [numpy array]: The diagonal DOCI Hamiltonian terms in matrix form
+
+    Returns:
+        projected_onebody_integrals [numpy array]: The corresponding one-body
+            integrals for the electronic structure Hamiltonian
+        projected_twobody_integrals [numpy array]: The corresponding two body
+            integrals for the electronic structure Hamiltonian
+    '''
+    n_qubits = hr1.shape[0]
+    projected_onebody_integrals = numpy.zeros((n_qubits, n_qubits))
+    projected_twobody_integrals = numpy.zeros(
+        (n_qubits, n_qubits, n_qubits, n_qubits))
+    for p in range(n_qubits):
+        projected_onebody_integrals[p, p] = hc[p] / 2
+        projected_twobody_integrals[p, p, p, p] = hr2[p, p]
+        for q in range(n_qubits):
+            if p == q:
+                continue
+            projected_twobody_integrals[p, q, q, p] = hr2[p, q] / 2
+            projected_twobody_integrals[p, p, q, q] = hr1[p, q]
+
+    return projected_onebody_integrals, projected_twobody_integrals
+
+
+def get_doci_from_integrals(one_body_integrals, two_body_integrals):
+    r"""Construct a DOCI Hamiltonian from electron integrals
+
+    Args:
+        one_body_integrals [numpy array]: one-electron integrals
+        two_body_integrals [numpy array]: two-electron integrals
+
+    Returns:
+        hc [numpy array]: The single-particle DOCI terms in matrix form
+        hr1 [numpy array]: The off-diagonal DOCI Hamiltonian terms in matrix
+            form
+        hr2 [numpy array]: The diagonal DOCI Hamiltonian terms in matrix form
+    """
+
+    n_qubits = one_body_integrals.shape[0]
+    hc = numpy.zeros(n_qubits)
+    hr1 = numpy.zeros((n_qubits, n_qubits))
+    hr2 = numpy.zeros((n_qubits, n_qubits))
+
+    for p in range(n_qubits):
+        hc[p] = 2 * one_body_integrals[p, p]
+        for q in range(n_qubits):
+            hr2[p, q] = (2 * two_body_integrals[p, q, q, p] -
+                         two_body_integrals[p, q, p, q])
+            if p == q:
+                continue
+            hr1[p, q] = two_body_integrals[p, p, q, q]
+
+    return hc, hr1, hr2
diff --git a/src/openfermion/ops/representations/interaction_operator.py b/src/openfermion/ops/representations/interaction_operator.py
index c1ccbe8..d10f2d8 100644
--- a/src/openfermion/ops/representations/interaction_operator.py
+++ b/src/openfermion/ops/representations/interaction_operator.py
@@ -15,6 +15,7 @@ import itertools
 import numpy
 
 from openfermion.ops.representations.polynomial_tensor import PolynomialTensor
+from openfermion.config import EQ_TOLERANCE
 
 
 class InteractionOperatorError(Exception):
@@ -150,3 +151,117 @@ def _symmetric_two_body_terms(quad, complex_valued):
         yield q, r, s, p
         yield s, p, q, r
         yield r, q, p, s
+
+
+def get_tensors_from_integrals(one_body_integrals, two_body_integrals):
+    '''Converts one and two-body integrals into tensor form
+
+    Arguments:
+        one_body_integrals [numpy array] -- the one-body integrals
+            of the given Hamiltonian
+        two_body_integrals [numpy array] -- the two-body integrals
+            of the given Hamiltonian
+    '''
+
+    n_qubits = 2 * one_body_integrals.shape[0]
+
+    # Initialize Hamiltonian coefficients.
+    one_body_coefficients = numpy.zeros((n_qubits, n_qubits))
+    two_body_coefficients = numpy.zeros(
+        (n_qubits, n_qubits, n_qubits, n_qubits))
+    # Loop through integrals.
+    for p in range(n_qubits // 2):
+        for q in range(n_qubits // 2):
+
+            # Populate 1-body coefficients. Require p and q have same spin.
+            one_body_coefficients[2 * p, 2 * q] = one_body_integrals[p, q]
+            one_body_coefficients[2 * p + 1, 2 * q +
+                                  1] = one_body_integrals[p, q]
+            # Continue looping to prepare 2-body coefficients.
+            for r in range(n_qubits // 2):
+                for s in range(n_qubits // 2):
+
+                    # Mixed spin
+                    two_body_coefficients[2 * p, 2 * q + 1, 2 * r + 1, 2 *
+                                          s] = (two_body_integrals[p, q, r, s] /
+                                                2.)
+                    two_body_coefficients[2 * p + 1, 2 * q, 2 * r, 2 * s +
+                                          1] = (two_body_integrals[p, q, r, s] /
+                                                2.)
+
+                    # Same spin
+                    two_body_coefficients[2 * p, 2 * q, 2 * r, 2 *
+                                          s] = (two_body_integrals[p, q, r, s] /
+                                                2.)
+                    two_body_coefficients[2 * p + 1, 2 * q + 1, 2 * r +
+                                          1, 2 * s +
+                                          1] = (two_body_integrals[p, q, r, s] /
+                                                2.)
+
+    # Truncate.
+    one_body_coefficients[
+        numpy.absolute(one_body_coefficients) < EQ_TOLERANCE] = 0.
+    two_body_coefficients[
+        numpy.absolute(two_body_coefficients) < EQ_TOLERANCE] = 0.
+
+    return one_body_coefficients, two_body_coefficients
+
+
+def get_active_space_integrals(one_body_integrals,
+                               two_body_integrals,
+                               occupied_indices=None,
+                               active_indices=None):
+    """Restricts a molecule at a spatial orbital level to an active space
+
+    This active space may be defined by a list of active indices and
+        doubly occupied indices. Note that one_body_integrals and
+        two_body_integrals must be defined
+        n an orthonormal basis set.
+
+    Args:
+        one_body_integrals: One-body integrals of the target Hamiltonian
+        two_body_integrals: Two-body integrals of the target Hamiltonian
+        occupied_indices: A list of spatial orbital indices
+            indicating which orbitals should be considered doubly occupied.
+        active_indices: A list of spatial orbital indices indicating
+            which orbitals should be considered active.
+
+    Returns:
+        tuple: Tuple with the following entries:
+
+        **core_constant**: Adjustment to constant shift in Hamiltonian
+        from integrating out core orbitals
+
+        **one_body_integrals_new**: one-electron integrals over active
+        space.
+
+        **two_body_integrals_new**: two-electron integrals over active
+        space.
+    """
+    # Fix data type for a few edge cases
+    occupied_indices = [] if occupied_indices is None else occupied_indices
+    if (len(active_indices) < 1):
+        raise ValueError('Some active indices required for reduction.')
+
+    # Determine core constant
+    core_constant = 0.0
+    for i in occupied_indices:
+        core_constant += 2 * one_body_integrals[i, i]
+        for j in occupied_indices:
+            core_constant += (2 * two_body_integrals[i, j, j, i] -
+                              two_body_integrals[i, j, i, j])
+
+    # Modified one electron integrals
+    one_body_integrals_new = numpy.copy(one_body_integrals)
+    for u in active_indices:
+        for v in active_indices:
+            for i in occupied_indices:
+                one_body_integrals_new[u, v] += (
+                    2 * two_body_integrals[i, u, v, i] -
+                    two_body_integrals[i, u, i, v])
+
+    # Restrict integral ranges and change M appropriately
+    return (core_constant,
+            one_body_integrals_new[numpy.ix_(active_indices, active_indices)],
+            two_body_integrals[numpy.ix_(active_indices, active_indices,
+                                         active_indices, active_indices)])
diff --git a/src/openfermion/ops/representations/polynomial_tensor.py b/src/openfermion/ops/representations/polynomial_tensor.py
index f92791b..7a340b0 100644
--- a/src/openfermion/ops/representations/polynomial_tensor.py
+++ b/src/openfermion/ops/representations/polynomial_tensor.py
@@ -127,17 +127,21 @@ class PolynomialTensor(object):
         """Initialize the PolynomialTensor class.
 
         Args:
-            n_body_tensors(dict): A dictionary storing the tensors describing
-                n-body interactions.
+            n_body_tensors(dict or None): A dictionary storing the tensors
+                describing n-body interactions. If None, n_body_tensors are
+                assumed to be generated on-the-fly by other data (for
+                subclassing purposes).
         """
-        self.n_body_tensors = n_body_tensors
-
-        # Set n_qubits
-        key_iterator = iter(n_body_tensors.keys())
-        key = next(key_iterator)
-        if key == ():
+        self._n_body_tensors = n_body_tensors
+        if n_body_tensors is None:
+            self._n_qubits = None
+        else:
+            # Set n_qubits
+            key_iterator = iter(n_body_tensors.keys())
             key = next(key_iterator)
-        self.n_qubits = n_body_tensors[key].shape[0]
+            if key == ():
+                key = next(key_iterator)
+            self._n_qubits = n_body_tensors[key].shape[0]
 
     @property
     def constant(self):
@@ -149,6 +153,18 @@ class PolynomialTensor(object):
         """Set the value of the constant term."""
         self.n_body_tensors[()] = value
 
+    @property
+    def n_body_tensors(self):
+        return self._n_body_tensors
+
+    @n_body_tensors.setter
+    def n_body_tensors(self, value):
+        self._n_body_tensors = value
+
+    @property
+    def n_qubits(self):
+        return self._n_qubits
+
     def __getitem__(self, args):
         """Look up matrix element.
 


Test patch:
diff --git a/src/openfermion/ops/representations/doci_hamiltonian_test.py b/src/openfermion/ops/representations/doci_hamiltonian_test.py
new file mode 100644
index 0000000..5c89a58
--- /dev/null
+++ b/src/openfermion/ops/representations/doci_hamiltonian_test.py
@@ -0,0 +1,219 @@
+#   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 doci_hamiltonian.py."""
+import os
+import unittest
+
+import numpy
+
+from openfermion.config import EQ_TOLERANCE
+from openfermion.chem.molecular_data import MolecularData
+from openfermion.config import DATA_DIRECTORY
+from openfermion.transforms import jordan_wigner
+from openfermion.linalg import get_sparse_operator
+from openfermion.ops.representations.doci_hamiltonian import (
+    DOCIHamiltonian, get_tensors_from_doci, get_projected_integrals_from_doci,
+    get_doci_from_integrals)
+
+
+class IntegralTransformsTest(unittest.TestCase):
+
+    def setUp(self):
+        self.geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))]
+        self.basis = 'sto-3g'
+        self.multiplicity = 1
+        self.filename = os.path.join(DATA_DIRECTORY, 'H2_sto-3g_singlet_0.7414')
+        self.molecule = MolecularData(self.geometry,
+                                      self.basis,
+                                      self.multiplicity,
+                                      filename=self.filename)
+        self.molecule.load()
+
+    def test_integrals_self_inverse(self):
+        hc, hr1, hr2 = get_doci_from_integrals(self.molecule.one_body_integrals,
+                                               self.molecule.two_body_integrals)
+        proj_one_body, proj_two_body = get_projected_integrals_from_doci(
+            hc, hr1, hr2)
+        hc_test, hr1_test, hr2_test = get_doci_from_integrals(
+            proj_one_body, proj_two_body)
+        self.assertTrue(numpy.allclose(hc, hc_test))
+        self.assertTrue(numpy.allclose(hr1, hr1_test))
+        print(hr2)
+        print(hr2_test)
+        self.assertTrue(numpy.allclose(hr2, hr2_test))
+
+    def test_integrals_to_doci(self):
+        one_body_integrals = self.molecule.one_body_integrals
+        two_body_integrals = self.molecule.two_body_integrals
+        hc, hr1, hr2 = get_doci_from_integrals(one_body_integrals,
+                                               two_body_integrals)
+        self.assertEqual(hc.shape[0], 2)
+        self.assertEqual(hr1.shape[0], 2)
+        self.assertEqual(hr2.shape[0], 2)
+
+        for p in range(2):
+            self.assertEqual(
+                hc[p] + hr2[p, p],
+                2 * one_body_integrals[p, p] + two_body_integrals[p, p, p, p])
+            for q in range(2):
+                if p != q:
+                    self.assertEqual(hr1[p, q], two_body_integrals[p, p, q, q])
+                    self.assertEqual(
+                        hr2[p, q], 2 * two_body_integrals[p, q, q, p] -
+                        two_body_integrals[p, q, p, q])
+
+
+class DOCIHamiltonianTest(unittest.TestCase):
+
+    def setUp(self):
+        self.geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))]
+        self.basis = 'sto-3g'
+        self.multiplicity = 1
+        self.filename = os.path.join(DATA_DIRECTORY, 'H2_sto-3g_singlet_0.7414')
+        self.molecule = MolecularData(self.geometry,
+                                      self.basis,
+                                      self.multiplicity,
+                                      filename=self.filename)
+        self.molecule.load()
+
+    def test_n_body_tensor_errors(self):
+        doci_hamiltonian = DOCIHamiltonian.zero(n_qubits=2)
+        with self.assertRaises(TypeError):
+            doci_hamiltonian.n_body_tensors = 0
+        with self.assertRaises(IndexError):
+            _ = doci_hamiltonian[((0, 0), (0, 0))]
+        with self.assertRaises(IndexError):
+            _ = doci_hamiltonian[((0, 0), (0, 0), (0, 0), (0, 0))]
+        with self.assertRaises(IndexError):
+            _ = doci_hamiltonian[((1, 1), (0, 0))]
+        with self.assertRaises(IndexError):
+            _ = doci_hamiltonian[((0, 1), (2, 1), (3, 0), (8, 0))]
+
+    def test_errors_operations(self):
+        doci_hamiltonian = DOCIHamiltonian.zero(n_qubits=2)
+        doci_hamiltonian2 = DOCIHamiltonian.zero(n_qubits=3)
+        with self.assertRaises(TypeError):
+            doci_hamiltonian += 'a'
+        with self.assertRaises(TypeError):
+            doci_hamiltonian -= 'a'
+        with self.assertRaises(TypeError):
+            doci_hamiltonian *= 'a'
+        with self.assertRaises(TypeError):
+            doci_hamiltonian /= 'a'
+        with self.assertRaises(TypeError):
+            doci_hamiltonian += doci_hamiltonian2
+        with self.assertRaises(TypeError):
+            doci_hamiltonian -= doci_hamiltonian2
+
+    def test_adding_constants(self):
+        doci_hamiltonian = DOCIHamiltonian.zero(n_qubits=2)
+        doci_hamiltonian += 2
+        self.assertAlmostEqual(doci_hamiltonian.constant, 2)
+        doci_hamiltonian -= 3
+        self.assertAlmostEqual(doci_hamiltonian.constant, -1)
+
+    def test_basic_operations(self):
+        doci_hamiltonian1 = DOCIHamiltonian.zero(n_qubits=2)
+        doci_hamiltonian2 = DOCIHamiltonian.from_integrals(
+            constant=self.molecule.nuclear_repulsion,
+            one_body_integrals=self.molecule.one_body_integrals,
+            two_body_integrals=self.molecule.two_body_integrals)
+        self.assertTrue(doci_hamiltonian2 == doci_hamiltonian1 +
+                        doci_hamiltonian2)
+        self.assertTrue(doci_hamiltonian1 -
+                        doci_hamiltonian2 == doci_hamiltonian2 / -1)
+        self.assertTrue(doci_hamiltonian2 * 0 == doci_hamiltonian1)
+
+    def test_error(self):
+        doci_hamiltonian = DOCIHamiltonian.from_integrals(
+            constant=self.molecule.nuclear_repulsion,
+            one_body_integrals=self.molecule.one_body_integrals,
+            two_body_integrals=self.molecule.two_body_integrals)
+        with self.assertRaises(TypeError):
+            doci_hamiltonian[((1, 0), (0, 1))] = 1
+        with self.assertRaises(IndexError):
+            _ = doci_hamiltonian[((1, 0),)]
+        with self.assertRaises(IndexError):
+            _ = doci_hamiltonian[((1, 1), (0, 0))]
+        with self.assertRaises(IndexError):
+            _ = doci_hamiltonian[((0, 1), (1, 1), (0, 0), (2, 0))]
+
+    def test_getting_setting_constant(self):
+        doci_hamiltonian = DOCIHamiltonian.zero(n_qubits=2)
+        doci_hamiltonian.constant = 1
+        self.assertEqual(doci_hamiltonian[()], 1)
+
+    def test_getting_setting_1body(self):
+        doci_hamiltonian = DOCIHamiltonian.zero(n_qubits=2)
+        doci_hamiltonian.hc[0] = 2
+        doci_hamiltonian.hc[1] = 4
+        self.assertEqual(doci_hamiltonian[((0, 1), (0, 0))], 1)
+        self.assertEqual(doci_hamiltonian[((1, 1), (1, 0))], 1)
+        self.assertEqual(doci_hamiltonian[((2, 1), (2, 0))], 2)
+        self.assertEqual(doci_hamiltonian[((3, 1), (3, 0))], 2)
+
+    def test_getting_setting_hr2(self):
+        doci_hamiltonian = DOCIHamiltonian.zero(n_qubits=2)
+        doci_hamiltonian.hr2[0, 0] = 2
+        doci_hamiltonian.hr2[1, 1] = 4
+        self.assertEqual(doci_hamiltonian[((0, 1), (1, 1), (1, 0), (0, 0))], 1)
+        self.assertEqual(doci_hamiltonian[((1, 1), (0, 1), (0, 0), (1, 0))], 1)
+        self.assertEqual(doci_hamiltonian[((2, 1), (3, 1), (3, 0), (2, 0))], 2)
+        self.assertEqual(doci_hamiltonian[((3, 1), (2, 1), (2, 0), (3, 0))], 2)
+
+        doci_hamiltonian.hr2[0, 1] = 2
+        self.assertEqual(doci_hamiltonian[((0, 1), (2, 1), (2, 0), (0, 0))], 1)
+        self.assertEqual(doci_hamiltonian[((0, 1), (3, 1), (3, 0), (0, 0))], 1)
+        self.assertEqual(doci_hamiltonian[((1, 1), (2, 1), (2, 0), (1, 0))], 1)
+        self.assertEqual(doci_hamiltonian[((1, 1), (3, 1), (3, 0), (1, 0))], 1)
+
+    def test_getting_setting_hr1(self):
+        doci_hamiltonian = DOCIHamiltonian.zero(n_qubits=2)
+        doci_hamiltonian.hr1[0, 1] = 2
+        self.assertEqual(doci_hamiltonian[(0, 1), (1, 1), (3, 0), (2, 0)], 1)
+        self.assertEqual(doci_hamiltonian[(1, 1), (0, 1), (2, 0), (3, 0)], 1)
+
+    def test_from_integrals_to_qubit(self):
+        hamiltonian = jordan_wigner(self.molecule.get_molecular_hamiltonian())
+        doci_hamiltonian = DOCIHamiltonian.from_integrals(
+            constant=self.molecule.nuclear_repulsion,
+            one_body_integrals=self.molecule.one_body_integrals,
+            two_body_integrals=self.molecule.two_body_integrals).qubit_operator
+
+        hamiltonian_matrix = get_sparse_operator(hamiltonian).toarray()
+        doci_hamiltonian_matrix = get_sparse_operator(
+            doci_hamiltonian).toarray()
+        diagonal = numpy.real(numpy.diag(hamiltonian_matrix))
+        doci_diagonal = numpy.real(numpy.diag(doci_hamiltonian_matrix))
+        position_of_doci_diag_in_h = [0] * len(doci_diagonal)
+        for idx, doci_eigval in enumerate(doci_diagonal):
+            closest_in_diagonal = None
+            for idx2, eig in enumerate(diagonal):
+                if closest_in_diagonal is None or abs(eig - doci_eigval) < abs(
+                        closest_in_diagonal - doci_eigval):
+                    closest_in_diagonal = eig
+                    position_of_doci_diag_in_h[idx] = idx2
+            assert abs(closest_in_diagonal - doci_eigval) < EQ_TOLERANCE, (
+                "Value " + str(doci_eigval) + " of the DOCI Hamiltonian " +
+                "diagonal did not appear in the diagonal of the full " +
+                "Hamiltonian. The closest value was " +
+                str(closest_in_diagonal))
+
+        sub_matrix = hamiltonian_matrix[numpy.ix_(position_of_doci_diag_in_h,
+                                                  position_of_doci_diag_in_h)]
+        assert numpy.allclose(doci_hamiltonian_matrix, sub_matrix), (
+            "The coupling between the DOCI states in the DOCI Hamiltonian " +
+            "should be identical to that between these states in the full " +
+            "Hamiltonian bur the DOCI hamiltonian matrix\n" +
+            str(doci_hamiltonian_matrix) +
+            "\ndoes not match the corresponding sub-matrix of the full " +
+            "Hamiltonian\n" + str(sub_matrix))
diff --git a/src/openfermion/ops/representations/interaction_operator_test.py b/src/openfermion/ops/representations/interaction_operator_test.py
index 4fd9ed9..d664533 100644
--- a/src/openfermion/ops/representations/interaction_operator_test.py
+++ b/src/openfermion/ops/representations/interaction_operator_test.py
@@ -17,6 +17,14 @@ import itertools
 import numpy
 
 from openfermion.ops.representations import InteractionOperator
+from openfermion.ops.representations import get_active_space_integrals
+
+
+class ActiveSpaceIntegralsTest(unittest.TestCase):
+
+    def test_raises_error(self):
+        with self.assertRaises(ValueError):
+            _ = get_active_space_integrals(None, None, [], [])
 
 
 class InteractionOperatorTest(unittest.TestCase):
diff --git a/src/openfermion/ops/representations/polynomial_tensor_test.py b/src/openfermion/ops/representations/polynomial_tensor_test.py
index dc4f2d0..c608309 100644
--- a/src/openfermion/ops/representations/polynomial_tensor_test.py
+++ b/src/openfermion/ops/representations/polynomial_tensor_test.py
@@ -176,6 +176,11 @@ class PolynomialTensorTest(unittest.TestCase):
             two_body_spinful
         })
 
+    def test_set_n_body_tensors(self):
+        pt_temp = copy.deepcopy(self.polynomial_tensor_a)
+        pt_temp.n_body_tensors = {(): 3}
+        self.assertEqual(pt_temp.constant, 3)
+
     def test_setitem_0body(self):
         self.polynomial_tensor_a[()] = 1
         self.assertEqual(self.polynomial_tensor_a.n_body_tensors[()], 1)