catch bug due to numpy multicasting, write test for it

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

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)]))