diff --git a/kwant/system.py b/kwant/system.py
index ed15e86829d4c067e34792d5663745b0793cbb9b..e1ca468e4d09d6ecf9d2840a41b2783c64b91e3a 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -125,12 +125,12 @@ class System(object):
                     n_j = to_coord.get(j)
                     if n_j is None:
                         continue
-                    try:
-                        h_sub[to_off[n_j] : to_off[n_j + 1],
-                              from_off[n_i] : from_off[n_i + 1]] = \
-                              ham(j, i) if j != i else diag[i]
-                    except ValueError:
+                    h = ham(j, i) if j != i else diag[i]
+                    shape = (1, 1) if np.isscalar(h) else h.shape
+                    if shape != (to_norb[n_j], from_norb[n_i]):
                         raise ValueError(msg.format(i, j))
+                    h_sub[to_off[n_j] : to_off[n_j + 1],
+                            from_off[n_i] : from_off[n_i + 1]] = h
             return h_sub
 
         gr = self.graph
diff --git a/kwant/tests/test_system.py b/kwant/tests/test_system.py
index 2411b0d771f48ffd3917d3bc25ef5903214e52a9..893b4637cf519f1a7e96ca1297b6fb213bbf2f75 100644
--- a/kwant/tests/test_system.py
+++ b/kwant/tests/test_system.py
@@ -54,6 +54,10 @@ def test_hamiltonian_submatrix():
     sys2 = sys.finalized()
     assert_raises(ValueError, sys2.hamiltonian_submatrix)
     assert_raises(ValueError, sys2.hamiltonian_submatrix, None, None, True)
+    sys[(0,), (2,)] = 1
+    sys2 = sys.finalized()
+    assert_raises(ValueError, sys2.hamiltonian_submatrix)
+    assert_raises(ValueError, sys2.hamiltonian_submatrix, None, None, True)
 
 def test_energies():
     sys = kwant.Builder(kwant.TranslationalSymmetry([(-1, 0)]))