diff --git a/kwant/operator.pyx b/kwant/operator.pyx
index 1dd99ebbb547d754c8ff324f621f3add9fe8bbd0..2b9da33985d5c930455e5a37514b6d802d4211f1 100644
--- a/kwant/operator.pyx
+++ b/kwant/operator.pyx
@@ -628,13 +628,17 @@ cdef class Density(_LocalOperator):
         cdef int unique_onsite = not callable(self.onsite)
         # prepare onsite matrices
         cdef complex[:, :] _tmp_mat
-        cdef complex *M_a
+        cdef complex *M_a = NULL
         cdef BlockSparseMatrix M_a_blocks
+
         if unique_onsite:
             _tmp_mat = self.onsite
             M_a = <complex*> &_tmp_mat[0, 0]
+        elif self._bound_onsite:
+            M_a_blocks = self._bound_onsite
         else:
-            M_a_blocks = self._bound_onsite or self._eval_onsites(args)
+            M_a_blocks = self._eval_onsites(args)
+
         # loop-local variables
         cdef gint a, a_s, a_norbs
         cdef gint i, j, w
@@ -727,15 +731,23 @@ cdef class Current(_LocalOperator):
         # prepare onsite matrices and hamiltonians
         cdef int unique_onsite = not callable(self.onsite)
         cdef complex[:, :] _tmp_mat
-        cdef complex *M_a,
-        cdef complex *H_ab
+        cdef complex *M_a = NULL
+        cdef complex *H_ab = NULL
         cdef BlockSparseMatrix M_a_blocks, H_ab_blocks
+
         if unique_onsite:
             _tmp_mat = self.onsite
             M_a = <complex*> &_tmp_mat[0, 0]
+        elif self._bound_onsite:
+            M_a_blocks = self._bound_onsite
+        else:
+            M_a_blocks = self._eval_onsites(args)
+
+        if self._bound_hamiltonian:
+            H_ab_blocks = self._bound_hamiltonian
         else:
-            M_a_blocks = self._bound_onsite or self._eval_onsites(args)
-        H_ab_blocks = self._bound_hamiltonian or self._eval_hamiltonian(args)
+            H_ab_blocks = self._eval_hamiltonian(args)
+
         # main loop
         cdef gint a, a_s, a_norbs, b, b_s, b_norbs
         cdef gint i, j, k, w
@@ -837,15 +849,23 @@ cdef class Source(_LocalOperator):
         # prepare onsite matrices and hamiltonians
         cdef int unique_onsite = not callable(self.onsite)
         cdef complex[:, :] _tmp_mat
-        cdef complex *M_a,
-        cdef complex *H_aa
+        cdef complex *M_a = NULL
+        cdef complex *H_aa = NULL
         cdef BlockSparseMatrix M_a_blocks, H_aa_blocks
+
         if unique_onsite:
             _tmp_mat = self.onsite
             M_a = <complex*> &_tmp_mat[0, 0]
+        elif self._bound_onsite:
+            M_a_blocks = self._bound_onsite
         else:
-            M_a_blocks = self._bound_onsite or self._eval_onsites(args)
-        H_aa_blocks = self._bound_hamiltonian or self._eval_hamiltonian(args)
+            M_a_blocks = self._eval_onsites(args)
+
+        if self._bound_hamiltonian:
+            H_aa_blocks = self._bound_hamiltonian
+        else:
+            H_aa_blocks = self._eval_hamiltonian(args)
+
         # main loop
         cdef gint a, a_s, a_norbs
         cdef gint i, j, k, w