diff --git a/src/openfermion/linalg/sparse_tools.py b/src/openfermion/linalg/sparse_tools.py
index 28e5731a2..44aa51a5e 100644
--- a/src/openfermion/linalg/sparse_tools.py
+++ b/src/openfermion/linalg/sparse_tools.py
@@ -700,7 +700,7 @@ def expectation_computational_basis_state(operator, computational_basis_state):
                   If operator is a FermionOperator, it must be normal-ordered.
         computational_basis_state (scipy.sparse vector / list): normalized
             computational basis state (if scipy.sparse vector), or list of
-            occupied orbitals.
+            zeros and ones for occupied and unoccupied orbitals, respectively.
 
     Returns:
         A real float giving expectation value.
@@ -714,6 +714,9 @@ def expectation_computational_basis_state(operator, computational_basis_state):
     if not isinstance(operator, FermionOperator):
         raise TypeError('operator must be a FermionOperator.')
 
+    if not operator.is_normal_ordered():
+        raise ValueError('operator must be a normal ordered.')
+
     occupied_orbitals = computational_basis_state
 
     if not isinstance(occupied_orbitals, list):
@@ -730,7 +733,8 @@ def expectation_computational_basis_state(operator, computational_basis_state):
             expectation_value += operator.terms.get(((i, 1), (i, 0)), 0.0)
 
             for j in range(i + 1, len(occupied_orbitals)):
-                expectation_value -= operator.terms.get(((j, 1), (i, 1), (j, 0), (i, 0)), 0.0)
+                if occupied_orbitals[j]:
+                    expectation_value -= operator.terms.get(((j, 1), (i, 1), (j, 0), (i, 0)), 0.0)
 
     return expectation_value
 
diff --git a/src/openfermion/linalg/sparse_tools_test.py b/src/openfermion/linalg/sparse_tools_test.py
index 30a555e20..11b3f5ed8 100644
--- a/src/openfermion/linalg/sparse_tools_test.py
+++ b/src/openfermion/linalg/sparse_tools_test.py
@@ -727,6 +727,17 @@ def test_expectation_qubit_operator_not_implemented(self):
                 QubitOperator(), csc_matrix(([1], ([6], [0])), shape=(16, 1))
             )
 
+    def test_expectation_bad_operator_order(self):
+        operator = (
+            FermionOperator('2^ 2', 1.9)
+            + FermionOperator('2^ 1')
+            + FermionOperator('2^ 1 2 1^', -1.7)
+        )
+        state = [0, 1, 1]
+
+        with self.assertRaises(ValueError):
+            expectation_computational_basis_state(operator, state)
+
 
 class ExpectationDualBasisOperatorWithPlaneWaveBasisState(unittest.TestCase):
     def setUp(self):