diff --git a/src/openfermion/utils/rdm_mapping_functions.py b/src/openfermion/utils/rdm_mapping_functions.py
index 1ac4173c5..f7b26886a 100644
--- a/src/openfermion/utils/rdm_mapping_functions.py
+++ b/src/openfermion/utils/rdm_mapping_functions.py
@@ -202,4 +202,6 @@ def map_particle_hole_dm_to_one_pdm(phdm, num_particles, num_basis_functions):
     Returns:
         opdm (numpy.ndarray): the 1-RDM transformed from a 1-RDM.
     """
+    if num_particles > num_basis_functions:
+        raise ValueError('Number of particles cannot exceed number of basis functions.')
     return numpy.einsum('prrq', phdm) / (num_basis_functions - num_particles + 1)
diff --git a/src/openfermion/utils/rdm_mapping_functions_test.py b/src/openfermion/utils/rdm_mapping_functions_test.py
index 5b9216818..0ee0e9630 100644
--- a/src/openfermion/utils/rdm_mapping_functions_test.py
+++ b/src/openfermion/utils/rdm_mapping_functions_test.py
@@ -171,6 +171,14 @@ def test_phdm_conversions_h2_631g(self):
         test_tpdm = map_particle_hole_dm_to_two_pdm(true_phdm, molecule.fci_one_rdm)
         assert numpy.allclose(test_tpdm, molecule.fci_two_rdm)
 
+    def test_map_particle_hole_dm_to_one_pdm_edge_cases(self):
+        # Trigger division by zero
+        phdm = numpy.zeros((2, 2, 2, 2))
+        num_particles = 3
+        num_basis_functions = 2
+        with self.assertRaises(ValueError):
+            map_particle_hole_dm_to_one_pdm(phdm, num_particles, num_basis_functions)
+
     def test_phdm_conversions_h2_sto3g(self):
         filename = "H2_sto-3g_singlet_1.4.hdf5"
         molecule = MolecularData(filename=os.path.join(DATA_DIRECTORY, filename))