diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index 9c5610e30..335db0506 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -43,6 +43,9 @@
                                       jordan_wigner_dual_basis_hamiltonian,
                                       wigner_seitz_length_scale)
 
+from ._rdm_equality_constraints import (one_body_fermion_constraints,
+                                        two_body_fermion_constraints)
+
 from ._sparse_tools import (expectation,
                             expectation_computational_basis_state,
                             get_density_matrix,
@@ -51,6 +54,7 @@
                             is_hermitian,
                             jordan_wigner_sparse,
                             jw_hartree_fock_state,
+                            jw_number_restrict_operator,
                             qubit_operator_sparse,
                             sparse_eigenspectrum)
 
diff --git a/src/openfermion/utils/_rdm_equality_constraints.py b/src/openfermion/utils/_rdm_equality_constraints.py
new file mode 100644
index 000000000..de5ad21e2
--- /dev/null
+++ b/src/openfermion/utils/_rdm_equality_constraints.py
@@ -0,0 +1,115 @@
+#   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.
+
+"""Constraints on fermionic reduced density matrices"""
+from openfermion.ops import FermionOperator
+
+
+def one_body_fermion_constraints(n_orbitals, n_fermions):
+    """Generates one-body positivity constraints on fermionic RDMs.
+
+        The specific constraints implemented are known positivity constraints
+        on the one-fermion reduced density matrices. Constraints are generated
+        in the form of FermionOperators whose expectation value is known to be
+        zero for any N-Representable state. Generators are used for efficiency.
+
+    Args:
+        n_orbitals(int): number of spin-orbitals on which operators act.
+        n_fermions(int): number of fermions in the system.
+
+    Yields:
+        Constraint is a FermionOperator with zero expectation value.
+    """
+    # One-RDM trace condition.
+    constraint_operator = FermionOperator()
+    for i in range(n_orbitals):
+        constraint_operator += FermionOperator(((i, 1), (i, 0)))
+    if len(constraint_operator.terms):
+        constraint_operator -= FermionOperator((), n_fermions)
+        yield constraint_operator
+
+    # One-RDM Hermiticity condition.
+    for i in range(n_orbitals):
+        for j in range(i + 1, n_orbitals):
+            constraint_operator = FermionOperator(((i, 1), (j, 0)))
+            constraint_operator -= FermionOperator(((j, 1), (i, 0)))
+            if len(constraint_operator.terms):
+                yield constraint_operator
+
+
+def two_body_fermion_constraints(n_orbitals, n_fermions):
+    """Generates two-body positivity constraints on fermionic RDMs.
+
+        The specific constraints implemented are known positivity constraints
+        on the two-fermion reduced density matrices. Constraints are generated
+        in the form of FermionOperators whose expectation value is known to be
+        zero for any N-Representable state. Generators are used for efficiency.
+
+    Args:
+        n_orbitals(int): number of spin-orbitals on which operators act.
+        n_fermions(int): number of fermions in the system.
+
+    Yields:
+        Constraint is a FermionOperator with zero expectation value.
+    """
+    # Two-body trace condition.
+    constraint_operator = FermionOperator()
+    for i in range(n_orbitals):
+        for j in range(n_orbitals):
+            constraint_operator += FermionOperator(
+                ((i, 1), (j, 1), (j, 0), (i, 0)))
+    if len(constraint_operator.terms):
+        constraint_operator -= FermionOperator(
+            (), n_fermions * (n_fermions - 1))
+        yield constraint_operator
+
+    # Two-body Hermiticity condition.
+    for ij in range(n_orbitals ** 2):
+        i, j = (ij // n_orbitals), (ij % n_orbitals)
+        for kl in range(ij + 1, n_orbitals ** 2):
+            k, l = (kl // n_orbitals), (kl % n_orbitals)
+            constraint_operator = FermionOperator(
+                ((i, 1), (j, 1), (l, 0), (k, 0)))
+            constraint_operator -= FermionOperator(
+                ((k, 1), (l, 1), (j, 0), (i, 0)))
+            if len(constraint_operator.terms):
+                yield constraint_operator
+
+    # Contraction to One-RDM from Two-RDM.
+    for i in range(n_orbitals):
+        for j in range(n_orbitals):
+            constraint_operator = FermionOperator()
+            for p in range(n_orbitals):
+                constraint_operator += FermionOperator(
+                    ((i, 1), (p, 1), (p, 0), (j, 0)))
+            constraint_operator += FermionOperator(
+                ((i, 1), (j, 0)), -(n_fermions - 1))
+            if len(constraint_operator.terms):
+                yield constraint_operator
+
+    # Linear relations between two-particle matrices.
+    for ij in range(n_orbitals ** 2):
+        i, j = (ij // n_orbitals), (ij % n_orbitals)
+        for kl in range(ij, n_orbitals ** 2):
+            k, l = (kl // n_orbitals), (kl % n_orbitals)
+
+            # G-matrix condition.
+            constraint_operator = FermionOperator(
+                ((i, 1), (k, 0)), 1.0 * (j == l))
+            constraint_operator += FermionOperator(
+                ((i, 1), (l, 1), (k, 0), (j, 0)))
+            constraint_operator += FermionOperator(
+                ((i, 1), (l, 1), (j, 0), (k, 0)))
+            constraint_operator -= FermionOperator(
+                ((i, 1), (k, 0)), 1.0 * (j == l))
+            if len(constraint_operator.terms):
+                yield constraint_operator
diff --git a/src/openfermion/utils/_rdm_equality_constraints_test.py b/src/openfermion/utils/_rdm_equality_constraints_test.py
new file mode 100644
index 000000000..d3f4c8642
--- /dev/null
+++ b/src/openfermion/utils/_rdm_equality_constraints_test.py
@@ -0,0 +1,76 @@
+#   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 _rdm_equality_constraints.py"""
+import unittest
+import os
+
+from openfermion.config import THIS_DIRECTORY
+from openfermion.transforms import get_interaction_operator
+from openfermion.utils import (MolecularData,
+                               one_body_fermion_constraints,
+                               two_body_fermion_constraints)
+
+
+class FermionConstraintsTest(unittest.TestCase):
+
+    def setUp(self):
+
+        # Setup.
+        n_atoms = 2
+        geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))]
+        basis = 'sto-3g'
+        multiplicity = 1
+        filename = os.path.join(THIS_DIRECTORY, 'data',
+                                'H2_sto-3g_singlet_0.7414')
+        molecule = MolecularData(
+            geometry, basis, multiplicity, filename=filename)
+        molecule.load()
+        self.n_fermions = molecule.n_electrons
+        self.n_orbitals = molecule.n_qubits
+
+        # Get molecular Hamiltonian.
+        self.molecular_hamiltonian = molecule.get_molecular_hamiltonian()
+        self.fci_rdm = molecule.get_molecular_rdm(use_fci=1)
+
+    def test_one_body_constraints(self):
+        for constraint in one_body_fermion_constraints(
+                self.n_orbitals, self.n_fermions):
+            interaction_operator = get_interaction_operator(
+                constraint, self.n_orbitals)
+            constraint_value = self.fci_rdm.expectation(interaction_operator)
+            self.assertAlmostEqual(constraint_value, 0.)
+            for term, coefficient in constraint.terms.items():
+                if len(term) == 2:
+                    self.assertTrue(term[0][1])
+                    self.assertFalse(term[1][1])
+                else:
+                    self.assertEqual(term, ())
+
+    def test_two_body_constraints(self):
+        for constraint in two_body_fermion_constraints(
+                self.n_orbitals, self.n_fermions):
+            interaction_operator = get_interaction_operator(
+                constraint, self.n_orbitals)
+            constraint_value = self.fci_rdm.expectation(interaction_operator)
+            self.assertAlmostEqual(constraint_value, 0.)
+            for term, coefficient in constraint.terms.items():
+                if len(term) == 2:
+                    self.assertTrue(term[0][1])
+                    self.assertFalse(term[1][1])
+                elif len(term) == 4:
+                    self.assertTrue(term[0][1])
+                    self.assertTrue(term[1][1])
+                    self.assertFalse(term[2][1])
+                    self.assertFalse(term[3][1])
+                else:
+                    self.assertEqual(term, ())
diff --git a/src/openfermion/utils/_sparse_tools.py b/src/openfermion/utils/_sparse_tools.py
index c7f78e7ff..b50927b47 100644
--- a/src/openfermion/utils/_sparse_tools.py
+++ b/src/openfermion/utils/_sparse_tools.py
@@ -598,7 +598,7 @@ def expectation_three_body_db_operator_computational_basis_state(
                                  dual_basis_action[4][0] % 2 ==
                                  orbital3 % 2)):
                             expectation_value += numpy.exp(-1j * (
-                                k1ad + k2bf + k3ce)) 
+                                k1ad + k2bf + k3ce))
 
                         # Handle -\delta_{ad} \delta_{be} \delta_{cf} after FT.
                         # The Fourier transform is spin-conserving.